A Self-Exciting Modelling Framework for Forward Prices in Power Markets

We propose and investigate two model classes for forward power price dynamics, based on continuous branching processes with immigration, and on Hawkes processes with exponential kernel, respectively. The models proposed exhibit jumps clustering features. Models of this kind have been already proposed for the spot price dynamics, but the main purpose of the present work is to investigate the performances of such models in describing the forward dynamics. We adopt a Heath-Jarrow-Morton approach in order to capture the whole forward curve evolution. By examining daily data in the French power market, we perform a goodness-of-fit test and we present our conclusions about the adequacy of these models in describing the forward prices evolution.


Introduction
Energy markets, and in particular, electricity markets, exhibit very peculiar features. The historical series of both futures and spot prices include seasonality, mean-reversion, spikes and small fluctuations. One can alternatively describe the power price dynamics by modelling the spot or the forward price. In the former case, the spot price can be obtained as a limit of the forward price when the maturity is close to the current time, in the latter case it is possible to derive the forward price from the spot by computing the conditional expectation with respect to a suitable risk-neutral measure of the spot price at the maturity.
After the pioneering paper by Schwartz [37], where an Ornstein-Uhlenbeck dynamics is assumed to describe the spot price behaviour, several different approaches have been investigated in order to describe the power price evolution. A comprehensive literature review until 2008 is offered in the book by Benth et al. [5]. A similar effort has been devoted to identify reliable models for the forward price dynamics, and a huge amount of literature is available focusing on Heath-Jarrow-Morton type models as in Benth et al. [4] and in Filimonov et al. [30], in the attempt to provide a description of the whole forward curves dynamics, in analogy with forward interest rates in fixedincome markets as in Heath et al. [23]. Some of the classical models proposed include jumps and/or stochastic volatility. Benth and Paraschiv [3] propose a random field approach based on Gaussian random fields by adopting the Musiela parametrization in order to describe the forward curve dynamics. Empirical evidence suggests that in many assets prices often jumps appear in cluster, thus requiring the introduction of jump processes exhibiting a clustering or self-exciting behaviour.
Kiesel and Paraschiv [27] recently presented a systematic empirical investigation of electricity intraday prices and suggested an approach based on Hawkes processes in order to describe the power price dynamics with jump clustering features. Self-exciting features in electricity prices attracted already some attention by several authors: Herrera and Gonzalez [24] proposed a selfexcited model for electricity spot prices, while Christensen et al. [8], Clements et al. [9] pointed out that time between spikes has a significant impact on the likelihood of future occurrences, thus providing a strong support to models including self-exciting properties.
The large class of models available in the literature describing the power price dynamics is then widening in order to include models exhibiting self-exciting features.
We also mention the paper by Jiao et al. [25], where a model based on continuous branching processes with immigration for power spot prices was proposed, and the forward prices computed with respect to a suitable structure preserving equivalent martingale measure.
Eyjolfsson and Tjøshteim [14] describe a class of Hawkes processes and present an empirical investigation based on data from UK power market supporting Hawkes-type models for spot prices.
The purpose of the present paper is to investigate if self-exciting features can arise in the power forward prices evolution as well, and in order to perform this investigation we shall focus on two different model classes: the Continuous Branching Processes with Immigration (CBI henceforth) and the Hawkes processes. While CBI processes are always affine, Hawkes processes in general are not, but when the kernel describing the intensity dynamics is of exponential type they are, and this feature makes the Hawkes processes with exponential kernel appealing from the modeling point of view. By considering the two model classes mentioned before, i.e. CBI and Hawkes processes, we then want to provide the description of the full term structure of power forward prices, following a Heath-Jarrow-Morton approach.
Power is a flow commodity, this meaning that instantaneous forward contracts are not directly traded on the market, but futures (sometimes called flow forward) are. So, in order to perform any kind of inference on the model proposed, it is necessary to extract the relevant information on the forward dynamics included in the futures prices. This can be done by applying suitable optimization procedures proposed in the literature, eventually modified in order to provide the best performances in the case under examination. These procedures are far from trivial from the computational point of view and require a careful implementation of the optimization step. We deliberately chose to work on daily data, in order to show how self-exciting effects can arise not only on a small time scale, but also at a coarser level.
The paper is organized as follows: in Section 2 we introduce the processes on which our models are based and in Section 3 we present and discuss the models proposed for the forward power price dynamics. In Section 4 we discuss the dynamics of Futures contracts when the forward dynamics is assumed to be given by the models introduced. From Sections 5 to 8, we provide the theoretical background and numerical results relative to the calibration/parameters' estimation for the model proposed. In the final section we provide some concluding remarks and discuss future extensions of the present work.

Continuous Branching Processes with Immigration
We now introduce our modeling framework for the electricity price, which is based on stochastic differential equations driven by Lévy random fields. We consider a Lévy random field, which is a combination of a Gaussian random measure W and a compensated Poisson random measure N independent of W . For a background on such general stochastic equations with jumps, we refer the readers e.g. to Dawson and Li [11], Li and Ma [32] and Walsh [38].
Let us now briefly introduce all the relevant ingredients of our work and recall some preliminary results. We fix a probability space (Ω, A, P). A white noise W on R 2 + is a Gaussian random measure such that, for any Borel set A ∈ B(R 2 + ) with finite Lebesgue measure |A|, W (A) is a normal random variable of mean zero and variance |A| and if A 1 , · · · , A n are disjoint Borel sets in B(R 2 + ), then W (A 1 ), · · · , W (A n ) are mutually independent. We denote by N the Poisson random measure on R 3 + with intensity λ which is a Borel measure on R 3 + defined as the product of the Lebesgue measure on R + × R + with a Borel measure µ on R + such that Note that µ is a Lévy measure since ∞ 0 (1 ∧ z 2 )µ(dz) < +∞. Recall that for each Borel set B ∈ B(R 3 + ) with λ(B) < +∞, the random variable N (B) has the Poisson distribution with parameter λ(B). Moreover, if B i , i = 1, . . . , n are disjoint Borel sets in B(R 3 + ), then N (B 1 ), · · · , N (B n ) are mutually independent. We let N = N − λ be the compensated Poisson random measure on R 3 + associated to N .
We introduce the filtration F = (F t ) t 0 as the natural filtration generated by the Lévy random field (see Dawson and Li [11]) and satisfying the usual conditions, namely, for any Borel subset A ∈ B(R + ) and B ∈ B(R 2 + ) of finite Lebesgue measure, the processes (W ([0, t] × A), t ≥ 0) and We consider the following stochastic differential equation in the integral form. Let a, b, σ, γ ∈ R + be constant parameters. Consider the equation: (2.1) where W (ds, du) is a white noise on R 2 + with unit covariance, N (ds, du, dz) is an independent compensated Poisson random measure on R 3 + with intensity ds du µ(dz) with µ(dz) being a Lévy measure on R + and satisfying ∞ 0 (z ∧ z 2 )µ(dz) < ∞. The integrals appearing in Equation (2.1) (and in the following) are both in the sense of Walsh [38]. It follows from Dawson and Li [11,Theorem 3.1]  Our model actually belongs to the family of CBI processes. Continuous Branching Processes with Immigration (CBI) are a class of stochastic processes commonly used in modelling population dynamics as in Padoux [35]. The self-exciting features, arising from the integrals in Equation (2.1) extended on the domain [0, Y (s)) with respect to the integration variable u, describe the growth of the population due to the reproduction of the previous generations. In the present modelling framework they just describe jumps generated by previous jumps. We briefly recall the definition by Kawazu and Watanabe [26] [Def. 1.1]. A Markov process Y with state space R + is called a CBI process characterized by branching mechanism Ψ(·) and immigration rate Φ(·), if its characteristic representation is given, for p ≥ 0, by: where E y denotes the conditional expectation with respect to the initial value Y (0) = y. The function v : R + × R + → R + satisfies the following differential equation: and Ψ and Φ are functions of the variable q ≥ 0 given by with σ, γ ≥ 0, β ∈ R and π, ν being two Lévy measures such that It is proved in Dawson and Li [11,Theorem 3.1] that the process in Equation (2.1) is a CBI process with the branching mechanism Ψ given by: and the immigration rate Φ(q) = abq.
The link between CBI processes and the affine term structure models has been established by Filipović [17]. If the process Y takes values in R + he proves equivalence between the two classes. We recall that the joint Laplace transform of a CBI process Y and its integrated process, which is given in Filipovic [17,Theorem 5.3], is defined as follows: for non-negative real numbers ξ and θ, we have: where v(t, ξ, θ) is the unique solution of ∂v(t, ξ, θ) ∂t = −Ψ(v(t, ξ, θ)) + θ, v(0, ξ, θ) = ξ. (2.7)

Hawkes Processes
A Hawkes process is a special counting process with a random intensity function. We introduce now the Hawkes processes with exponential kernel. They can be written as follows : where the last term is an Ito integral, N t is the number of jumps in the interval between 0 and t and J(dz, ds) is a Poisson random measure with intensity λ(t), satisfying the SDE: Here β > 0 is the rate of exponential decay of the influence of previous jumps on the intensity level and α the amplitude of the memory kernel, t i are the jumps times and Z i the jump sizes, which we shall assume distributed according to an exponential density with parameter δ, so that only positive jumps appear in both Equations (2.8) and (2.9), and we can writeJ(ds, dz) = J(ds, dz) − λ(s)µ(dz)ds and µ(dz) = δ exp (−δz)dz, whereJ(ds, dz) denotes the compensated version of the Poisson measure J(ds, dz). We assume the following condition holds: β − α/δ > 0, granting the non-explosiveness of the Hawkes process (see e.g. Bernis et al. [6]).
Hawkes processes with exponential kernel are the only class of Hawkes processes exhibiting both the Markov property and an affine structure (see e.g. Errais et al. [13]). The have been extensively used in order to describe the dynamics of several asset classes, including equities as in Hainaut and Moraux [20], commodities as in Eyjolfsson and Tjøsteim [14], exchange rates as in Rambaldi et al. [36] and credit risk as in Errais et al. [13].

Forward Prices Modelling
In this section we are going to introduce the two alternative models for the forward prices, that we are going to test against electricity market data. In both cases the price at time t of a forward contract with maturity T ≥ t is additive and it can be defined as follows where Λ(t) is a deterministic seasonality function that will be made precise later on, n is the number of factors used and each of the terms X i is an underlying factor, whose dynamics will be specified in the following Subsections 3.1 and 3.2.

The Forward Model based on CBI
Our first model assumes the following dynamics for the factors X i , i = 1, . . . , n: Namely, the X i 's evolve in time with respect to the historical measure P according to Equation (2.1) with immigration rate b i = 0. By recalling that the intensity of the Poisson random measure It is possible to re-write Equation (3.10) as follows: or, equivalently, as The relation between the dynamics of the forward price with respect to the historical measure P and the risk-neutral dynamics, written with respect to Q, can be easily obtained by applying the following result, proved in the paper by Jiao et al. [25,Proposition 4.1].
Proposition 3.1. Let X 1 , X 2 , · · · , X n be independent CBI processes where for each i ∈ {1, · · · , n}, X i is a CBI process under the probability measure P, with dynamics given by Eq. 3.11. Assume that the filtration F = (F t ) t≥0 is generated by the random fields W 1 , W 2 , · · · , W n and N 1 , N 2 , · · · , N n . For each i, fix η i ∈ R and ξ i ∈ R + and define Then the Doléans-Dade exponential E(U ) is a martingale under P and the probability measure Q defined by Remark 3.1. In this context, the parameters η i , ξ i can be interpreted as the Market Price of Risk associated with the diffusion/jump part of X i , i = 1, . . . , n, respectively.
Remark 3.2. In order to avoid arbitrage opportunities we shall assume that the de-seasonalized dynamics of every factor X i is a local martingale under Q and this will automatically imply that a i = 0 under Q. Since the first integral is defined with respect to the Gaussian white noise W i (ds, du) and the second integral is defined with respect to the compensated Poisson random measure N i (ds, du, dz), each process X i (t, T ) is in fact a local martingale with respect to Q.
, specifying the relations between the model parameters under the riskneutral measure Q and the historical measure P, it is clear that in the present modelling framework, for each factor X i , a mean reversion speed coefficient a i can be non-null under P and zero under Q. As far as the immigration term b i is concerned, if it vanishes under Q, it will be zero under any equivalent probability measure.
Assumption 3.1. In the estimation procedure applied to the real market data we shall assume that only one process of the type introduced in Equation (3.11) will drive the forward curve dynamics.

The Forward Model Based on Hawkes Processes
As alternative to the model proposed in the previous subsection, we consider,under P, Equation (3.11) for the instantaneous forward price, where now each X i , i = 1, . . . , n satisfies a SDE of the following form: whereJ i (dz, ds) are compensated marked point process with intensity λ i (t), satisfying the SDE: We assume the jump size distributed according to an exponential density with parameter δ i for each (λ i , X i ), so we can write: (3.20) Remark 3.4. The choice of a square-root process for the diffusion part of the forward curves dynamics is motivated by the positivity requirement as well as the choice of the exponential distribution for the jumps size.
In order to make the presentation of the two model classes more homogeneous, we can introduce the Dawson-Li representation for the Hawkes-type dynamics as well and write the SDE governing the dynamics of forward prices under the historical measure P as follows: where the definition of the integrals and the notations are the same as in Subsection 2.1 and the λ i (t) evolve according to Eq. (3.19).
It is immediate to remark that the dynamics described by the two model classes look almost identical when written in the Dawson-Li representation, the main difference being the specification of the equation governing the evolution of the intensity processes. This is one of the reasons behind the choice of these two alternative models to describe the forward prices' evolution.
The dynamics just described is given with respect to the historical probability measure P. In order to obtain a description with respect to the risk-neutral measure Q we need to introduce a measure change. The following proposition provides a measure change preserving the Hawkestype dynamics. A proof can be found in Bernis et al. [7].
and define: Then the Doléans-Dade exponential E(U ) is a martingale under P and the probability measure Q defined by dQ dP Ft := E(U ) t is equivalent to P. The dynamics with respect to Q takes the following form: Remark 3.5. In this context, the parameters η i , ξ i can be interpreted as the Market Price of Risk associated with the diffusion/jump part of the i − th factor X i , respectively.
Remark 3.6. We shall assume, as for the previous model, that the de-seasonalized dynamics of X i is a local martingale under Q and this will automatically imply that the mean reversion speed c i of any X i must vanish under Q. Both the diffusion and the jump terms are in fact local martingales with respect to Q.
Remark 3.7. From the formulas in the previous lines, specifying the relations between the model parameters under the risk-neutral measure Q and the historical measure P, it is clear that in the Hawkes modeling framework, for each factor X i , a mean reversion speed coefficient c i can be nonzero under P and zero under Q. A non zero mean-reverting term can then appear in the dynamics written with respect to the historical measure P, although this term vanishes under Q.
Assumption 3.2. In the estimation procedure applied to the real market data we shall assume that only one process of the type introduced in Equation (3.11) will drive the forward curve dynamics.

The Futures Dynamics
We focus here rigorously on forward contracts delivering a quantity of energy over a finite period of time. We shall refer to them as futures, even if in the literature they are sometimes called swaps or flow forwards.
where f (t, ·) is the price at time t of the forward contract to be paid upon delivery. The value at time t of a Futures contract with delivery period [T 1 , T 2 ] is given, in our modelling framework, by (recall Equation (3.10)): (4.23) By introducing the dynamics of the factors X i into the above equation, we get the following equation describing the futures' dynamics under the risk-neutral probability Q both in the CBI framework (recall Equation (3.12)): and in the Hawkes setting (recall Equation (3.18)): Assumption 4.3. From now on, in view of our numerical analysis, we will assume that one driving factor is sufficient. Namely, we will consider the case i = 1.
In order to rule out arbitrage opportunities the prices of futures with different delivery periods must satisfy specific time-consistency relations. In particular, the value of a futures contract with delivery period [T 1 , T n ] is linked to the values of the contracts with delivery on intervals ] represents a partition of the interval [T 1 , T n ], by the following relation: (4.24) The situation is described by the following picture, describing the so called "Cascade unpacking mechanism":

Data Analysis: from Futures Prices to Forward Curves
From a theoretical point of view the contracts are settled continuously over the delivery period, as you can see from Equation (4.22), but in practice they are settled at discrete times. Assuming settlement at N points in time u 1 < u 2 < . . . < u N , with u 1 = T 1 , u N = T 2 , and ∆ i = u i+1 − u i , then the discrete version of Equation (4.22) becomes where, again, w(u, T 1 , T 2 ) = 1 T 2 −T 1 . The main goal of what follows is to provide a forward dynamics formulation starting from the futures prices that we observe in the market. What we are going to do is to build a smooth curve describing today's forward prices from quoted futures prices, according to the Heath-Jarrow-Morton framework outlined before. This is a well studied problem in literature and there are basically two approaches to do so: either fitting a parametric function to the entire yield curve by regression, or fitting all observed yields with a spline (see for example Anderson and Deacon [1] for a survey on different methods for constructing yield curves). Here we follow the second approach. Throughout the paper we will also use the notation T s i , T e i to denote the first (start) and the last (end) day of the delivery period of the i-th contract, and T s , T e in case of no ambiguity.
Our data set consists of French futures closing prices downloaded from Thomson Reuters, that span over a period of 17 years, from 2002 to 2019. These contracts are divided with respect to the duration of the delivery period into: weekly (tickers F7B1-B5), monthly (ticker F7BM), quarterly (ticker F7BQ) and yearly (ticker F7BY) contracts. For each of these we have 4 typologies of rolling contracts, namely c1, c2, c3 and c4, where c1 and c4 are the ones with the closest and the farthest delivery period, respectively (for an example to see how rolling contracts work see Section 5.2). On the market we observe the quantity F (0, T s , T e ) for every contract, for different choices of T s , T e (T e − T s = 7 days for the weekly, T e − T s = 30 days for the monthly, T e − T s = 90 days for the quarterly and T e − T s = 365 days for the yearly), where "0" is the current date, the first available being July 1, 2002, while the last available being March 15, 2019, for a total number of 4234 current dates. More precisely, for each day, representing the "0" day, we have a different number of contracts with different delivery periods, depending on the data availability of that day. We want to extract the curve f (0, u) for all the different choices of the "0" date.
Notation 5.1. When possible from now on we will write f (u) instead of f (0, u) and F (T 1 , T 2 ) instead of F (0, T 1 , T 2 ) to shorten the notation.
All the code and the computations have been implemented in MATLAB R2018a, on a CPU 2.6 GHz and 12 Gb of RAM HP Notebook with Windows 10.

Extracting Smooth Forward Curves from market data
Obtaining a smooth curve of forward prices from futures prices is a well studied problem in the literature, see for example Fleten and Lemming [18]. The initial condition for using a Heath-Jarrow-Morton approach when modelling forwards is a smooth curve describing today's forward prices, which must be extracted from the futures prices observed in the market. We will follow the approach by Benth et al. [5,Ch. 7] by imposing the following Assumption 5.4. The forward curve can be represented as the sum of two continuous functions Λ(u) and ε(u):

25)
where T s is the starting day of the settlement period for the contract with the closest delivery period and T e is the last day of the settlement period for the contract with the farthest delivery period. We interpret Λ(u) as a seasonality function and ε(u) to be an adjustment function that captures the forward curve's deviation from the seasonality.
For the specification of the seasonality function we follow Benth et al. [5], namely we define The parameter a ∈ R + is obtained by finding the minimum of the prices over all the contracts, while b is the (normalized 1 ) distance between the end of the last day of the year from the day when the minimum occurs. This procedure leads to: a = 13.600, b = 1358.038. (5.27) There are several other methods for extracting the seasonality function from the data (see for example Paraschiv [34], and Kiesel et al. [28] for an application to hourly data), but since this topics is not the main focus of our study, we prefer to stick on the well known method proposed by Benth et al. [2] and systematically described in Benth et al. [5] (Chap.7, Sect.7.2.1).
We shall see now how the adjustment function ε is obtained.

The function ε: a maximum smooth forward curve
By following the approach followed by Benth et al. [5] and we follow a maximum smoothness criterion applied to the adjustment function ε.
Remark 5.9. One may ask why the maximum smoothness criterion is applied only to the adjustment function ε and not to the entire forward function f . This ensures the presence of a seasonality pattern that, otherwise, would have possibly been smoothed out.
The properties we require for the adjustment function are that it is twice continuously differentiable and horizontal at time T e , i.e. ε (T e ) = 0. and such that the closing prices matching condition holds (this is made precise in Equation (5.34)).
We interpret the smoothest forward curve (5.25) to be the one for which ε solves the minimization problem above, with Λ chosen as in (5.26) and a, b as in (5.27).

A smooth forward curve constrained by closing prices
In this subsection we present the general procedure to extract the forward dynamics in a general situation with a fixed number of contracts from the market, but we will often make references to our own case. Before presenting the algorithm we need to introduce a procedure in order to deal with overlapping periods. Let be a list of start and end dates for the settlement periods of m different futures contracts for a given day (in our case, m = 16). We need to be able to handle the problem of overlapping settlement periods to rule out arbitrage opportunities. This was a concrete issue working with our data because it happens that, in a given day, two or more contracts have delivery periods that intersect. To overcome this, we construct a new list of dates T , namely where overlapping contracts are split into sub-periods. In our case n is typically 24, T 0 denotes the starting day of the contract with the closest delivery period, while T n denotes the last day of the contract with the farthest delivery period. The procedure is illustrated in the following figure: settlement period for the first contract As we can see from Figure 2, the elements of this new list are basically the elements in T sorted in ascending order, with duplicate dates removed. The futures prices could be taken into account either by exact matching or by a constraint on the bid-ask spread prices. Dealing with closing prices, here we impose an exact matching prices on closing prices (see Equation (5.34)). From now on we denote with F C i the closing price for the future i, i ∈ {1, . . . , m}.
The adjustment functions ε is chosen in the class C, namely (with a slight abuse of notation we use ε(u; x) instead of ε(u) to stress the dependence on x) . . . a n u 4 + b n u 3 + c n u 2 + d n u + e n , u ∈ [T n−1 , T n ].
where x = [a 1 , b 1 , c 1 , d 1 , e 1 , . . . a n , b n , c n , d n , e n ] is the row vector of the coefficients of the splines that we want to find. In this way we have, roughly speaking, a spline for every settlement period. To find the unknown parameters x = [a 1 , b 1 , c 1 , d 1 , e 1 , . . . a n , b n , c n , d n , e n ] in order to fully recover the adjustment function, we need to solve the following equality constrained convex quadratic programming problem (5.29) subject to the following constraints: i) continuity of the derivatives up to second order at the knots, for j = 1, . . . , n − 1, 32) ii) flatness at the end (see Equation (5.28)) ε (T n ; x) = 0, (5.33) iii) matching of the closing prices (see Equation (4.22)), for i = 1, . . . , m, In this way the minimisation problem (5.29) has a total of 3n + m − 2 constraints (i.e., 3(n − 1) constraints from (5.30)-(5.32), one constraint from (5.33) and m constraints from (5.34)). By computing the second derivative of ε and inserting it in Equation (5.29) and integrating for every delivery period, we can rewrite the minimisation problem (5.29) as min x∈R 5n x Hx, (5.35) where for j = 1, . . . , n, and l = 1, . . . , 5.
We clearly see that the constraints (5.30)-(5.34) are linear w.r.t. x, and so they can be formulated in a matrix form as Ax = b, where A is a (3n+m−2)×5n-dimensional matrix, and b is a (3n+m− 2)-dimensional vector. Solving the problem (5.29) with the constraints (5.30)-(5.34) is equivalent to solving (5.35) with the constraints written in the form Ax = b. Let λ = [λ 1 , λ 2 , . . . , λ 3n+m−2 ] be the corresponding Lagrange multiplier vector to the constraints (5.30)-(5.34). So, we can now express (5.29) as the following unconstrained minimization problem min x∈R 5n ,λ∈R 3n+m−2 x Hx + λ (Ax − b). (5.38) Remark 5.10. The advantage of dealing with problem (5.38), instead of (5.29) with the constraints (5.30)-(5.34), is that (5.38) is a unconstrained problem that can be simply solved. Indeed the solution [x,λ] is obtained just solving the linear system The dimension of the left matrix is (8n + m − 2) × (8n + m − 2). Solving (5.39) numerically is standard, and can be done using various techniques (e.g. QR or LU factorisation) 2 . As you can see from Figure 3, the c1 contract is the closer one to the current date and its delivery period spans from 01/02/2014 to 28/02/2014. After 28/02/2014, there is a rollover from the c1 contract to the c2 contract and so on.

Numerical Results
The following Figure 4 shows the plot of the futures closing prices for c1, c2, c3 and c4 contract for the weekly contract. In the x-axis there are the different dates, while in the y-axis there is the price. As you can see the presence of seasonality is pretty strong and this could also be seen from the monthly and quarterly contracts. Since we were worried that our analysis could have been affected by the presence of the quarterly contracts being sensible to the seasonality pattern, we did the analysis in both cases, with and without the quarterly contracts. After the analysis was performed, we noticed that the results were coherent in both cases, so from now on we will focus only on contracts different from quarterly.
The algorithm presented takes approximately 40 seconds to extract the 4234 different curves, i.e. the curves f (0, u), for all the 4234 different values of "0".   (0, u), for a total number of 4234 curves. In the x-axis we have the time to maturity while in the y-axis we have the prices. Note that different colours in the curves mean different type of contracts. Note also that the further we move on the x-axis, the flatter the curves become. This is in line with the flatness constraint in Equation (5.28). The reason behind the difference in the shapes of the curves is to be investigated in the price constraint but also in the nature of the contracts, since for each day the number and the type of available contracts were different, having to deal also with overlapping settlement periods.

Jump Detection
We now want to detect the jumps. We will be only dealing with positive jumps since we have supposed that the jump size is distributed according to an exponential density, as already described in Section 3. In the next subsection we will describe an algorithm that allows to detect jumps, and, as a by-product, which also gives their size.

Description of the Algorithm
In order to detect jumps we proceed in the following way: for a fixed maturity T , we define the vertical section at maturity T , where the parameter t ranges through all the curves, i.e. t = 1, . . . , 4234. Roughly speaking, looking at Figure 5, this is nothing but the intersection between the vertical line x = T and the curves. There are several ways to detect jumps from the data, maybe the more natural one consisting in fixing a threshold Θ ∈ R + and saying that a jump occurs at timet if |Vt +1 − Vt| ≥ Θ. We follow here an iterative weighted least square approach. Define n = 4234 the total number of curves and N = {1, 2, . . . , n − 1}. The algorithm to detect jumps reads as follows: 2. Identify all the t ∈ N such that V t+1 −Vt √ Vt ≥ 3σ 1 and denote by M 1 ⊆ N this family of indices, so that m 1 = |M 1 |; 6. Stop when finding k ∈ N such that m k = |M k | = 0 (no new jumps are detected).
This procedure finds, at every iterations, new jumps. Clearly as the number of iterations increases, σ i decreases and so the jumps detected become smaller and smaller. After several tests on the data, we noticed that stopping at k ∈ N such that m k = 0 would lead to too many jumps, of which the last detected are much smaller compared to the ones discovered at the first iterations. So we chose, as a good compromise, to stop the algorithm after the first two iterations.

Jumps Analysis
We selected different values of T , namely T = 200, T = 400 and T = 700 days. This covers all the different shapes of the forward curves and so it represents a good sampling of our data. After applying the algorithm described in Subsection 6.1 we end up with the following pictures showing the size and distribution of the jumps detected at T = 200, represented by the orange vertical lines, together with the corresponding price plot, represented by the continuous blue line:  As one can see from the pictures above, at T = 200 the jumps detected are the bigger ones with respect to their amplitude, and this is not surprising looking at Figure 5, where one can clearly see that the price movements are pretty significant at T = 200. On the other hand, when T = 700, the jumps detected are quite small and this is due to the fact that at T = 700 the curves are pretty flatten, leading to prices which are close to each other.

Parameters Estimation
Before starting with the statistical tests on the two models, we still have to estimate: the size of the jumps and the parameters characterizing the drift and the volatility coefficients. We start with δ, the jumps' size. Recall that in Sections 3.2 and 3.1 we assumed for both models the jumps' size to be distributed like an exponential random variable with parameter δ > 0. Let z i be the size of the i-th jump, where i = 1, . . . , L and L is the number of jumps at the chosen maturity T ( see Table 1). Then δ can be estimated e.g. via its Maximum Likelihood Estimator: We obtain what follows (the fact that we have the smallest values ofδ at T = 200 is not surprising at all, because at the beginning jumps are bigger, as said before): T 200 400 700 δ 0.064282682 0.194300598 0.430239733 Table 2: Parameter estimation for δ.
We now need to estimate the parameter appearing in the drift coefficient of our forward dynamics.
Remark 7.11. Notice that neither in Equation (3.12) (stated in the CBI framework), nor in Equation (3.21) (given for the Hawkes case) the mean-reversion term appears, but it does if we pass under the measure P by exploiting, respectively, Propositions 3.1 and 3.2. In both cases we end up with the following dynamics under P, for a fixed T : where X denotes the factor appearing in the forward dynamics without the seasonality and with no jumps (i.e., we remove all the times at which a jump has occurred).
In order to estimate a we simply discretize the above equation, by writing a = 1 − X(t + 1, T )  As you can see from Table 3, the estimated value of a in all the three cases is really small, very close to 0.
Now it remains to estimate the volatility parameter σ, appearing in both Equations (3.12) and in Equation (3.21). By recalling the iterative algorithm presented in Subsection 6.1 and taking as σ the value of σ 2 (namely, the estimation after the second iteration), we get

Testing the Models
In this section we want to perform statistical tests concerning the intensity of the jumps. We want to check what is the best process modelling the jumps we have detected before. We test the two models based on Hawkes and branching processes, plus the Poisson, which is a toy-model: (0) Poisson process; (1) Hawkes process; (2) Self exciting branching process.
We will mainly rely on the Kolmogorov-Smirnov (KS) test, namely we will test the null hypothesis H 0 , stating that the data have the same cumulative distribution function as the one coming from one of the above models, against the alternative H 1 . We fix a significance level equal to 0.05.

Jump Intensity Estimation
Before using the KS test to check whether the jumps distribution comes from one of the three models, we need to estimate the intensity from our data. The input in all the cases will be the time occurrences of the jumps over [0, T ] (for the three different values of T ), 0 < τ 1 < τ 2 < · · · < τ N = T , where N can take the values 86, 74 and 79 depending on the chosen maturity T , as you can check from Table 1.
(0) [Poisson] The (constant) intensity, λ P > 0, is estimated as the ratio between the total number of (positive) jumps and the sum of the inter-times between two consecutive jumps.
(1) [Hawkes] Here the intensity is given by Equation (3.19), so this case will be treated in a separate subsection.
(2) [Branching] In this case the stochastic intensity λ B (t) ∝ X(t, T ) and the constant of proportionality γ is estimated, for a fixed T , as the ratio between the total number of (positive) jumps and the cumulative (de-seasonalized) forward prices.

The Hawkes Setting: Estimating λ
Recalling Equation (3.19), it is clear that we have to estimate three parameters: λ(0), α and β. We mainly rely on the paper by Ozaki [33] and we will find a Maximum Likelihood Estimation (MLE). The log-likelihood of a Hawkes process whose response function is of the form αe −βt , is given by log (λ(0) + αA(i)) , (8.42) where A(i) = τ j <τ i e −β(τ i −τ j ) for i ≥ 2 and A(1) = 0.
In order to estimate the parameters λ(0), α, β, we need to find the maximum of the function in (8.42), which is a real value function of three variables. The maximum was found using the command fminsearch of Matlab.
The following three tables show the parameters estimated for the jump intensity of the three different models at the maturities T = 200, T = 400 and T = 700. As far as the branching model is concerned, λ B (t) is proportional to the process X(t, T ) (for a fixed T ), and the parameter to be estimated is the γ, which is the constant ratio between the two processes.

Model
Parameters  Table 6: Parameters estimation for the three models at T = 400. Table 7: Parameters estimation for the three models at T = 700.

KS test for the models
We perform a Kolmogorov-Smirnov test in order to check which of the proposed distributions best models the jumps in our data. At the end of the subsection we will provide the p-values to conclude.

(0) [Poisson]
We check whether the jumps inter-times are drawn from an exponential distribution with parameter λ P , where λ P is the one given in Tables 5, 6 and 7. The p-value is automatically given by Matlab via the function kstest.
(1) [Hawkes] We mainly adapt the methods in Lallouache and Challet [29] to our purpose. In particular, we check if the time-deformed series of durations {θ i } i=1,...,N , defined by has an exponential distribution of parameter 1, where λ t is the intensity estimated before (the estimated parameters can be found in Tables 5, 6 and 7), and where recall that the τ i 's are the jumps arrival times. The p-value is automatically given by Matlab via the function kstest.
(2) [Branching] The procedure here is quite different from the ones adopted before and it is the object of the following subsection.

Setting the KS Test for the Branching Model
The KS test we will perform in this case was constructed based on the following classical result. where in our case λ B (s) ∝ f (s, T ), for a fixed T .
Remark 8.12. Recall that λ B (s) ∝ f (s, T ), for any fixed T . It is crucial to notice, from Proposition 8.3, that the distribution of the arrival times is independent of the factor of proportionality connecting λ B and f (·, T ).
So in this case we perform a KS test, comparing the cumulative distribution function F (t) in Equation (8.44) with the empirical one relative to the jump times. Since F (t) given in Proposition (8.3) is not a priori associated to a known distribution, we cannot use the Matlab command kstest and we have to rely on the classical theory on the KS test. The KS statistics for the test is D n = sup x∈R |S n (x) − F (x)| , (8.45) where n is the number of our data (recall Table 1), F is the cumulative distribution function in Equation (8.44) and S n (x) is the empirical cumulative distribution function of the jump arrival times. We find: In the following Figure 7 we graphically compare the two cumulative distribution functions at the three different maturities.   In order to obtain the p-value in the branching case, we apply the asymptotic results in Facchinetti [15] in the case when the dataset is greater than 35. We provide the p-values in the following  Table 9: p-value test for the three models for different T As one can see from Table 9, the hypothesis that the intensity follows a Poisson process is highly rejected, as we expected. In the branching case the hypothesis is also rejected, even if the p-value in this case was much closer to the acceptance level of 0.05. For the Hawkes case the test fails to reject the hypothesis since all the three values are above our level of acceptance.
As a conclusion we can resume the main achievement presented in the present paper. We proposed two alternative models for power forward prices evolution, based on a HJM approach, extending to forward prices dynamics two models already proposed for the spot price dynamics [14], [25] . After extracting forward curves from quoted futures prices, we proposed a parameters estimation method for both models and then we performed a test on the adequacy of the two models in describing the observed forward prices evolution. The final conclusion of our test is that the hypothesis that forward prices follow a CBI-type dynamics is rejected, while the hypothesis of a Hawkes type dynamics is not. This conclusion suggests that self-exciting effects can arise in power forward dynamics as well as in the spot dynamics, and that an approach based on Hawkes processes can capture these effects in a natural and parsimonious way.