{'fraction_non_alphanumeric': 0.045333622889407676, 'fraction_numerical': 0.01169987512894294, 'mean_word_length': 4.0360902255639095, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A multi-factor polynomial framework for long-term electricity forwards with delivery period

We propose a multi-factor polynomial framework to model and hedge long-term electricity contracts with delivery period. This framework has various advantages: the computation of forwards and correlation between different forwards are fully explicit, and the model can be calibrated to observed electricity forward curves easily and well. Electricity markets suffer from non-storability and poor medium- to long-term liquidity. Therefore, we suggest a rolling hedge which only uses liquid forward contracts and is risk-minimizing in the sense of F\"ollmer and Schweizer. We calibrate the model to over eight years of German power calendar year forward curves and investigate the quality of the risk-minimizing hedge over various time horizons.


Introduction
Electricity differs from other energy commodities due to specific features such as limited storability, possibility of intra-day negative prices, its unique mechanism of the auction market, high liquidity of short-to medium-term trading and illiquidity of its long-term trading. Much of the academic literature is dedicated to short-to medium-term modeling of electricity spot and futures prices, as its highly frequent and huge data amount makes it ideal for empirical studies of time series analysis. However, the literature addressing the modeling of long-term electricity forwards and the corresponding hedging problems is scarce.
In this paper, we propose a mathematically tractable multi-factor polynomial diffusion framework to model long-term forwards. In this framework the computation of forwards and cross-maturity correlations are fully explicit. Fitting the model to long time series of single market electricity data works easily and well. Furthermore, we set up a rolling hedge mechanism that only uses liquid forward contracts. This allows us to address the non-storability of electricity and poor liquidity in its long-term markets. Within the setup the hedging strategy we suggest minimizes the conditional variance of the cost processes at any time, and thus is risk-minimizing in the sense of Föllmer and Schweizer. A simulation study using the estimated model shows that the risk-minimizing rolling hedge significantly reduces, yet does not fully eliminate, the variance and skew of the long-term exposures.
The proposed modeling framework has various applications in forward modeling. It can be used to smoothly extrapolate the curve to the non-liquid horizon while calibrating it to the liquid horizon; it can also be used to smooth the forward surface implied by the market once calibrated and to filter out market noise; moreover, it can be used to model the prices within the real data horizon between two quotation dates. Furthermore, the model can be extended to model multiple electricity markets and other energy markets simultaneously. It can thus serve as an alternative model for risk management purposes, and for conducting simulations. We do however not pursue such multi-market extensions in this paper.
Polynomial models have been used to solve a number of problems in finance, see [19,2,3,13,17,1,18,10,15] for references as well as [14,16] for a treatment of the underlying mathematical theory. With the exception of [20], the polynomial processes have not been used for electricity modeling. Our polynomial framework makes assumption on properties of spot and forward and not on supply-demand relation, and thus falls into the category of classical reduced-form model (see [11] for details on reduced-form model versus structural approach). It is closest to the arithmetic models of Benth et al (see [8,7,9]), and extends them by making the spot price not a linear combination but a squared combination of underlying polynomial processes. In doing so we extend the class of stochastic process on the one hand, and guarantee non-negative spot prices on the other hand.
The local risk-minimization hedging criterion of Föllmer & Schweizer 1991 is one the two main quadratic hedging approaches; see e.g. [21,25,35,34] for references of the general theory of local risk-minimization, [22] for the mean-variance hedge, and [26] for a comparison of the two approaches. In a recent paper on hedging, a locally risk-minimizing hedge was given for the arithmetic model of Benth et al. under illiquidity; see [12]. Our work differs from theirs, as we consider a rolling hedge which only uses liquid forward contracts and give explicit expression for the locally risk-minimizing hedging strategy for our modeling framework.
This paper is structured in the following way: In Section 2, we define the underlying polynomial framework, model the spot price as a quadratic function of it, and provide two main specifications. In Section 3, we briefly review the main characteristics of polynomial diffusions, with a focus on the moment formula for polynomials of degree two. In Section 4, we define electricity forwards with and without delivery period. We give pricing formulas for forwards, as well as explicit expressions for covariances and correlations between different forwards. In Section 5, to incorporate time series observations of forward prices, we specify a market price of risk function, which determines the forward price dynamics under the real-world measure P. In Section 6, we introduce a rolling hedge mechanism with liquidity constraints for hedging a long-term electricity commitment. Further, we give a rollinghedge that is locally risk-minimizing in the sense of Föllmer and Schweizer. In Section 7, we perform model estimation of a specification of the polynomial framework to a time series of real observations of power forwards using a quadratic Kalman filter. Further we simulate forward curves and investigate the quality of the risk-minimizing hedge over various time horizon.
Throughout this paper, we fix a filtered probability space (Ω, F, F t , Q), where Q is a risk-neutral probability measure used for pricing. For simplicity we assume zero interest rate and thus apply no discounting. We denote by S d the set of all symmetric d×d matrices and S d + the subset consisting of positive semidefinite matrices. We let Pol n denote the space of polynomials on R d of degree at most n.

The model
In this section we define the underlying polynomial framework. Firstly, we model the spot price S t as a quadratic function of an underlying d-dimensional state variable X t which evolves according to a polynomial diffusion. More precisely, we let where p S (x) = c + x Qx with c ∈ R + and Q ∈ S d + , κ ∈ R d×d , θ ∈ R d , W a d-dimensional Brownian motion under Q and σ : R d → R d×d is continuous. We assume that the components of the diffusion matrix a(x) := σ(x)σ(x) are polynomials of degree at most two. This ensures that X t is a polynomial diffusion, see Lemma 2.2 in [16].
We will now focus on the following two specifications.
t ) a standard two-dimensional Brownian motion. Here Y t mean-reverts at rate κ Y towards the correlated process Z t . And thus, Y t and Z t can be seen as factor processes that drive the short-end and long-end dynamics of spot prices respectively. Let α, β, c ∈ R + and let the spot price be given by This guarantees nonnegative spot price, as S t ≥ c ≥ 0. This specification is of the form Specification 2.2 (Three-factor model). We now present a specification which extends the two-factor model by modeling correlation between the underlying processes stochastically via a Jacobi process. Conditions under which the model exists and is unique are given below. Let κ Z , κ Y ∈ R, κ R , σ Z , σ Y , σ R > 0, and θ R ∈ (−1, 1). The process X t := (Z t , Y t , R t ) evolves according to the SDE t ) a standard three-dimensional Brownian motion. Let α, β, c ∈ R + and let the spot price be given by Proposition 2.3. Recall that κ R > 0, θ R ∈ (−1, 1), and assume moreover that Then for any initial condition with Z 0 ∈ R, Y 0 ∈ R and R 0 ∈ (−1, 1), there exists a unique strong solution X t = (Z t , Y t , R t ) of the SDE (5). Furthermore, this solution satisfies R t ∈ (−1, 1) for all t ≥ 0.
Proof. In the following we show the existence and uniqueness of a strong solution R t as well as its boundary non-attainment. Once this is shown, we can explicitly find X t := (Y t , Z t ) in terms of R t . Indeed, Itô's formula yields We now prove existence of a weak solution of the SDE for R t . Let ϕ(r) be a continuous function that is equal to one for r ∈ [−1, 1] and is equal to zero for |r| > 2, for example We letb(r) : andσ(r) are continuous and bounded, and hence an R-valued weak solution R t exists for the SDE dR t =b(R t )dt +σ(R t )dW (3) t ; see Theorem 4.22 of Section 5.4D in [28]. We next show that R t stays in (−1, 1) using a version of "McKean's argument". Let p(r) := 1 − r 2 and note that p(R 0 ) > 0. Further define the stopping times τ n := inf{t : p(R t ) ≤ 1 n } and τ := lim n→∞ τ n . Observe that (7)- (8) for t < τ . Consider the process Then M t is a local martingale on the stochastic interval [0, τ ). By definition, this means that for all n ∈ N, M t∧τn is a local martingale. We now show that τ = ∞ a.s. Suppose for contradiction that P(τ < ∞) > 0. Then there exists a large T < ∞ such that P(τ < T ) > 0. Note that for all t < T ∧ τ . Thus M t∧T is uniformly bounded from below, and hence a local supermartingale on the stochastic interval [0, τ ). The supermartingale convergence theorem for processes on stochastic interval [0, τ ) now gives that lim t→τ M t∧T exists in R almost surely; see e.g. the proof of Theorem 5.7 in [16]. Hence, in view of (9), − log p(R t ) is pathwise bounded above on [0, T ∧ τ ), which in turn means that τ > T a.s. This contradiction shows that R t ∈ (−1, 1) for all t ≥ 0.
t . For the existence and uniqueness of strong solutions, we note that b(.) is Lipschitz continuous, and σ(.) is Hölder continuous of order 1/2. Hence, pathwise uniqueness holds for this SDE; see Theorem 3.5(ii) in [33]. As a result, any (−1, 1)-valued solution is a strong solution by the Yamada-Watanabe theorem; see e.g. Theorem 1.7 in [33].

Polynomial diffusions and moment formulas
In this section we briefly review some important results regarding polynomial diffusions; see e.g. [16] for more details. We also provide a moment formula for polynomials of degree two. Consider the (extended) generator G associated to the R d -valued polynomial diffusion X t introduced in Section 2, namely for x ∈ R d and any C 2 function f . By Itô's formula, the process is a local martingale. Since the components of a(x) are polynomials of degree at most two, it follows that for any n ∈ N and any polynomial p ∈ Pol n , Gp is also polynomial of the same degree or lower degree, i.e. Gp ∈ Pol n . Fix n and let N = d+n n be the dimension of Pol n . Let H : R d → R N be a function whose components form a basis of Pol n . Then for any p ∈ Pol n where p ∈ R N is the coordinate representation of p(x), and G ∈ R N ×N the matrix representation of the generator G.
Below we give a more explicit version of the moment formula for polynomials of degree two. Note that the quantity tr(π a(x)) with π ∈ S d is quadratic in x, and thus of the form for some a 0 (π) ∈ R, a 1 (π) ∈ R d , and a 2 (π) ∈ S d that depend linearly on π.
Theorem 3.2 (Moment formula for polynomials of degree two). Let q(x) be a polynomial of the form q(x) = q 0 + q x + x Qx with q 0 ∈ R, q ∈ R d and Q ∈ S d . Further let X t satisfy (2). Then for 0 ≤ t ≤ T we have: where φ, ψ, π solve the linear ODE Proof. Define Let τ = T − t. Itô's formula along with (13) and then (14) gives: . Then with Cauchy-Schwartz inequality, for some constant C ∈ R. Together with Tonelli's theorem, this bound yields which is finite by Theorem 3.1. Hence, M (t, X t ) is a square-integrable true martingale. As a result, This is the claimed formula.

The term structure of forward prices
In this section we define electricity forwards, present their pricing formulas and give expressions for covariances and correlations between different forwards. The price at time-t of an electricity forward with instantaneous delivery at time T ≥ t is given by In practice, electricity is not delivered instantaneously, but gradually over a period of time. This leads us to the following definition: the time-t price of an electricity forward with delivery period [T 1 , T 2 ), t ≤ T 1 < T 2 , is given by Note that a forward contract (financial or physical) can have settlement that takes place either before or after the delivery period. Discounting is not needed in the pricing, as the difference in cashflow can be evened out by the purchase of a bond of that time period. F (t, T 1 , T 2 , X t ) is often also referred to as swap price, as the delivery of underlying power happens over a period of time and thus the price is the averaged price over that period.
It is intuitive that a forward with delivery period is the summation of all forwards (with instantaneous delivery) that deliver at single time points within the delivery period; moreover, a forward with delivery period which collapses into one single time point should be priced the same as a forward with instantaneous delivery. The following proposition confirms this relationship between forwards with and without delivery period.
Proof. In view of (15) and (16), the first identity follows from the conditional version of Tonelli's theorem since S t is nonnegative. The second identity then follows from the fundamental theorem of calculus, using that f (t, T, X t ) is continuous in T , see Proposition 4.2 below.
The following result gives closed-form expression for the forward prices.
Proof. This follows from Theorem 3.1 and rearranging terms.
Note that G is a non-invertible matrix. Still, τ 0 e uG du is explicit; see Appendix A for the explicit computation. Specification 2.1 Recall the Specification 2.1 in Section 2. We consider the basis given by Then S t can be uniquely represented as: For any C 2 function f and x = (z, y) ∈ R 2 , the generator G is : Applying G to each element of H(X t ) gives its matrix representation, Specification 2.2 Recall Specification 2.2 in Section 2. Here the general G actually preserves a proper subspace of Pol 2 , namely the one spanned by the components of Therefore it is not necessary to include the remaining basis functions in the definition of H. Then S t can be uniquely represented as For any C 2 function f and x = (z, y, r) ∈ R 3 , the generator G is Applying G to each element of H(X t ) gives For later use, we briefly discuss the instantaneous quadratic covariation and correlations between different forwards and give explicit forms for both specifications. The instantaneous covariation between two forwards with instantaneous delivery at T 1 and T 2 is, at where We define the corresponding instantaneous correlation as with Σ(X t ) from (24). The matrices Σ(X t ) for Specification 2.1 and Specification 2.2 are given in Appendix B. Similarly, for t ≤ T 1 < T 2 and t ≤ T 3 < T 4 , the time-t instantaneous covariation of forwards with delivery periods [T 1 , and the time-t instantaneous correlation is: with Σ(X t ) from (24) and

Market price of risk specification
In order to incorporate time series observations of real-world forward curves, we must specify the forward dynamics under the real-world probability measure P. Thus, in this section, we specify a market price of risk function λ : for some γ ∈ R d and Λ ∈ S d×d , and denote the associated Radon-Nikodym density process by We choose γ and Λ such that M λ t is a true martingale. We can then define P on every finite time interval [0, T ] via its Radon-Nikodym density dP dQ | F T = M λ T . Then, by Girsanov's theorem, the P-dynamics of X t becomes . Consider now Specification 2.1. In this case M λ t is a true martingale for any choice for γ and Λ, as the following result shows.
Proposition 5.1. Let X t evolve according to (3). Then M λ t from (29) is a martingale.
Note that X t has drift and diffusion that are affine in X t ; see computations in Section 4, Appendix B and Section 7. To be explicit, let Λ = diag(λ Z , λ Y ) and γ = (γ Z , γ Y ) . Then the P-dynamics of X t is given by: This can also be written as , and σ(x) from (4).
In the case of Specification 2.2 it is a more delicate problem to determine those market price of risk parameters for which M λ t is a true martingale. Since we will not use Specification 2.2 is our empirical analysis, we do not consider this issue here.

Hedging
In this section we first describe a rolling hedge setup with constraints which addresses the illiquidity and non-storability issues when hedging a long-term electricity contract. Rolling hedges for commodities form a well-known hedging scheme; see for example [23,31]. We then briefly review the locally risk-minimizing hedge of Föllmer and Schweizer, and give a rolling hedge for our modeling framework that is risk-minimizing.

A rolling hedge setup
Suppose we have committed to deliver power from year T to year T + 1 for a large T ∈ N (e.g. T = 10 years) and our objective is to hedge this long-term electricity commitment. In our framework the time-t valuation of the commitment is Note that F t is a Q-martingale and the pricing formula (Proposition 4.2) gives explicit pricing at any t ∈ [0, T ]. In an interest rate context, the analogous hedging problem is rather easy: just buy bonds and hold them as the payout in 10 years is known in advance. For electricity the problem is more difficult for a number of reasons: • Long-term forwards are not liquidly traded (otherwise buy and hold the financial contracts as in the interest rate context); • Electricity cannot be stored without significant costs (otherwise cash and carry as for other storable commodities: simply buy the amount needed in [ T , T + 1] and hold).
• Only short-term / near-dated contracts with same delivery length is available. But its underlying commodity (electricity) is not the same as the one underlying a longterm contract because power is not storable. Some empirical evidence suggests that short-term prices carry limited information about what spot prices will be far into the future (see [24]).
One possible strategy in this case is a rolling hedge, where we take a long position in near-term contracts as a hedge, and roll the hedge going forward. The underlying assumption of this strategy is that the short-and long-term contracts are highly correlated, which is often the case in practice. Under this strategy, we first enter the first-nearby calendar forward, close it out when it matures (one day before the contract enters into the delivery period), and then enter into the new first-nearby contract, close it out on maturity, and so on, until the delivery period of the commitment starts. This mechanism is illustrated in Figure 1.
To formalize this, let us first define the price process P t containing all calendar-year forwards with a one-year delivery period (short: cal forward): . . .
Hedge here using F (t, 1, 2, X t ) Hedge here using F (t, 2, 3, X t ) Hedge here using F (t, T , T + 1, X t ) where N = T and P N t = F t . Note that each P k t is a Q-martingale by its definition (16). By Proposition 4.2, P k t can be expressed as where w 01 is defined in (28). An admissible 1 hedging strategy is an R N +1 -valued process ϕ t = (η t , ξ t ) = (η t , ξ 1 t , . . . , ξ N t ) , where η t is adapted (representing bank account) and ξ t is predictable (representing amount 1 Note that for any polynomial processes p(Xt) all moments of Pt := E Q [p(XT )|Ft] exist. Therefore, integration with respect to any moments of P is well-defined. And thus, ξ ∈ L 2 (P ), i.e. E Q T 0 ξ s d P sξs < ∞, and ϕ := (η, ξ) is admissible.
of tradable assets or hedge ratio), and satisfies The constraint (32) reflects the liquidity issue and trading rule of those markets: • only the first-nearby forwards are liquid; • a contract that has started to deliver can no longer be traded.

The value process (or the P & L) at time
The cumulative cost of the hedge up to time t is: where G t denotes the cumulative gain of the hedge up to time t: Note that the market is incomplete under the restriction (32), since there are two different Brownian motions, but only one risky asset to invest in at any given time. In an incomplete market a claim generally cannot be fully replicated at maturity by a selffinancing hedging strategy. Depending on the restriction on cash account η, one can either use a strategy that is self-financing but does not perfectly replicate the claim at maturity, or use a strategy that fully replicates the claim at maturity but needs additional investment throughout the hedge, i.e. is not self-financing. In the first case, we have residual risk and in the latter case additional cash infusion is needed. Either way, risk cannot be fully eliminated and can only be minimized. In the following we briefly review the concept of risk-minimizing strategy in the sense of Föllmer and Schweizer, and then give a rolling hedge that is locally risk-minimizing.

A locally risk minimizing hedging criterion
The risk-minimization criterion proposed and developed by Föllmer and Schweizer (see e.g. [26], [25], [35], [34], [21] for details), is to minimize the conditional variance R t (ϕ) of the cost process C t (ϕ), among all not necessarily self-financing strategies ϕ that perfectly replicate F at maturity: In our setup, (34) is equivalent to η T N = 0 and ξ N T N = 1. A strategy ϕ * is called risk-minimizing if for any ϕ that satisfies (34) we have R t (ϕ * ) ≤ R t (ϕ), Q-a.s. for every t ∈ [0, T ]; see Schweizer (page 545 in [34]). One can show that any risk-minimizing strategy is mean self-financing, i.e. C t (ϕ) is a Q-martingale. Föllmer and Schweizer showed that the existence of such a strategy ϕ is guaranteed if the price process P t is a Q-local martingale. Moreover, in the martingale case, finding such a strategy is equivalent to finding the Galtchouk-Kunita-Watanabe (GKW) decomposition of F , namely where ξ is an admissible, predictable process and L is a square-integrable Q-martingale strongly orthogonal to P with L 0 = 0. The risk-minimizing hedging strategy ϕ rm is then given by where the value process is V t (ϕ rm ) = E[ F |F t ] = F t = F 0 + t 0 ξ F s dP s + L t and the cost process is C t (ϕ rm ) = F 0 + L t . Obviously this risk-minimizing strategy satisfies V T (ϕ rm ) = F T , and the associated risk process R t (ϕ rm ) is minimal (zero) at t = T .

A risk-minimizing rolling hedge
Recall that the price process P t is a Q-martingale. Then the time-t valuation of the longterm electricity commitment F T has a GKW-decomposition as in (35). We now compute the process ξ in this decomposition. This will give us the hedging strategy. Using (35), (32) and (33), we obtain for any t ∈ [k − 1, k), k ∈ N: Therefore, the risk-minimizing hedging strategy of the tradable assets is given by where And thus, for t ∈ [k − 1, k), the cash account η rm t is then given by and the associated cost process is Remark 6.1. The risk minimizing strategy also minimizes the quadratic covariation between the claim and and the value of hedge without the cash account. Indeed, formally one has This expression is minimized by ξ k t = d P k , F t d P k ,P k t as in (37).

Empirical analysis
In this section we demonstrate the use of our polynomial framework for modeling and hedging long-term electricity forwards and analyzing their performance. Based on a time series of real observations of power forwards provided by Axpo Solutions AG, we estimate parameters of a model specification. Further we simulate forward curves and investigate the quality of risk-minimizing hedges over various time horizons.

The data
Electricity long-term contracts lack liquidity and are not available on exchange. 2 In fact, long-term forwards with delivery periods are only offered by a small group of market participants over the counter (OTC), mostly by energy producing and trading companies. The data we use are provided by Axpo Solutions AG, and come originally from Totem Markit service, which surveys prices of various electricity contracts from each member firm and in term provides market consensus prices. More concretely, the data are German calender-year baseload (Cal) forwards that are quoted monthly from January 2010 to April 2018. 3 On each quotation date, we have at most 10 quoted contracts, i.e. first to tenth nearby Cal forwards. For each quoted contract on each quotation date, we have consensus price and the price spread between the highest quoted price and the lowest quoted price. A visualization of consensus prices is given in Figure 2.  Figure (a) shows the dynamics of each nearby Cal forward contract with respect to quotation date. We see that not every contract is available on every quotation date. In Figure (b), each curve is the forward curve of a quotation date, i.e. each curve shows the prices of the first to at most tenth nearby Cal forwards of that date. For the sake of a clearer view, we take a selection of forward curves. These curves (of chosen quotation dates) are stacked and time-lagged into a day. We note two shapes of forward curves: a straight contango curve and a curve which is flat with slight backwardation at the front and contango at the back end of the curve.

Model estimation
In order to capture the dynamics of the forward curves with our model, a non-linear filter is needed for model estimation, as the forward prices are quadratic in the Gaussian underlying factor process X t . Recall that the fundamental assumption of Kalman filter is that the measurement space is linear and Gaussian in the state space. Thus, in order to work with a Kalman filter, we can either linearize the quadratic relationship between state and measurement. This leads to a so-called extended Kalman filter. Alternatively, we can augment the state to incorporate the linear and quadratic terms of X t , so that the measurements become linear in the augmented state. Inspired by the work of Monfort, Renne and Roussellet [29], we will use a timedependent version of the latter approach to estimate a discrete version of Specification 2.1 based on the data from Section 7.1. The estimation will be under P, which means that we also need to estimate the market price of risk parameters.
Note that we do not have direct access to the underlying state process X t through the available data. Indeed, at each quotation date t k , we only see the prevailing price F j k of the j-th nearby forward contract, with j = 1, . . . , 10. 4 We view F j k as a noisy observation of the model price. More precisely, we assume that where F (t k , T j , T j +1) is the model price computed using Proposition 4.2, η j k are iid standard Gaussian noise, modulated by some parameters N j k > 0. The role of N j k is to encode the trustworthiness of the price of the j-th nearby contract on quotation date t k . A large value means that the price is considered noisy and uncertain, and a small value that the price is considered accurate. The N j k are chosen based on the spreads δ j k between the highest and lowest quoted price for the j-th nearby contract on date t k . Specifically, we use where δ j denotes the time series average of the spreads δ j k for a fixed maturity j, and δ denotes the overall average of all the spreads δ j k . The use of iid noise corresponds to assuming that our model captures all systematic effects. This is a standard assumption to reduce the complexity of the estimation.
A quadratic Kalman filter for Specification 2.1 We will now overload notation in the following manner: we write X k for the state X t k at quotation date t k , and similarly for other quantities that depend on time.
Since model prices at date t k are quadratic in the state X k , we have the expression for some a j k ∈ R, B j k ∈ R 2 and C j k ∈ S 2 that can be deduced from the pricing formula in Proposition 4.2. In view of (17), and following [29], we observe that F j k is affine in the augmented state vector . Specifically, the vector of prices, F k = (F 1 k , . . . , F 10 k ) is given by where a k := (a 1 k , . . . , a 10 k ) and B k := ( B 1 k , . . . , B 10 k ) can be computed as follows: for each maturity j = 1, . . . , 10, we have (18) and G from (19). Moreover, we have defined N k := diag N 1 k , . . . , N 10 k and η k := (η 1 k , . . . , η 10 k ) . Next, the discretized (non-augmented) state dynamics is given by Here we use the market price of risk parameters Λ = diag(λ Z , λ Y ) and γ = (γ Z , γ Y ) from Section 5. The discretized dynamics of the augmented state X k is where the involved quantities are conveniently expressed using the standard vector stacking operator V ec(), Kronecker product ⊗, selection matrix H d , and duplication matrix G d . The resulting expressions are: where Σ := KK and I d is the identity matrix of size d, and Λ m is the standard commutation matrix of size m 2 × m 2 . We then let K(X k−1 ) be the Cholesky factor of Σ(X k−1 ), i.e., K(X k−1 ) K(X k−1 ) = Σ(X k−1 ). We finally define F k−1 := σ(F k−1 , F k−2 , ..., F 1 ). The filtering algorithm is then described in Algorithm 1, where we use the notation

Algorithm 1 Quadratic Kalman filtering algorithm
Anchoring: gives the prediction error. Update: Optimization with the quadratic Kalman filter for Specification 2.1 For the model estimation with the quadratic filter, we use both the Least-Squares (LS) and the Maximum Likelihood (ML) criteria. We start with LS, as it is robust and converges fast. Once a stable result is obtained, we apply ML to obtain further improvement. Moreover, we impose 1 ≥ κ Y ≥ κ Z ≥ 0 on the parameters, in line with the interpretation that Y t and Z t drive the short and the long end of the forward curve respectively and thus mean-revert at different speed. The filtered underlying process X t = (Z t , Y t ) is given in Figure 3. The estimated parameters are shown in Table 1.
In the implementation we use the R package DEoptim, which is an optimizer based on a differential evolution algorithm; see [36], [32] for details of the algorithm and https: //cran.r-project.org/web/packages/DEoptim/index.html, [4], [6], [30] [5] for use of the package. Figure 4 gives a visualization of the model estimation using Specification 2.1. We quantify the goodness of fit in terms of relative errors, both cross-sectionally at each quotation date (Figure 5a), and across time for each nearby forward contract (Figure 5b). The overall relative error, i.e. the average relative error across all contracts and quotation dates, is as low as 0.661%, indicating a very good model fit.  Table 1: Estimated parameters of Specification 2.1.
In Figure 5a we notice a single spike of the time series of averaged errors reaching almost 2% (on a quotation date in February 2016). This is due to a single dramatic price drop of a forward curve on that date that is moderately captured by our model as it is continuous and gives smooth prices.
Looking at the estimation of the time series of each nearby forward (Figure 5b), we find that the front end fit (i.e. the first nearby to the sixth nearby forward contract) works very well while the fit deteriorates for longer maturities. This occurs by construction, as the prices of contracts with very long time-to-maturity are less reliable than those on the front end of the forward curve. In the filter this is captured by the data variance N j t , which is influenced by the time series of price spread of each forward; in general N j t tends to be higher for longer time-to-maturity (i.e. larger j).
We also performed model estimation under Q. This is equivalent to assuming P = Q, meaning that the market price of risk is zero (λ(X t ) = 0). This produces different parameters than those in Table 1, but the fit remains remarkably good.

Simulation and hedging analysis
In the following, we simulate forward surfaces, run locally risk-minimizing hedging strategies on those, and analyze their performance with respect to different hedging horizons.
Simulation of forward surfaces With a given set of parameters, we generate samples of entire forward surfaces over a fixed time horizon T . This can be done efficiently by first simulating the P-dynamics of the underlying process X t = (Y t , Z t ) until year T using a simple Euler scheme (with, say, N discretization steps). We can then compute the forward price for the 1-st through L-th nearby contract at each point t ≤ T on the time grid by applying the pricing formula, Proposition 4.2. The complexity of simulating M evolutions of forward curves is of the order O((M × N ) L ). A brief pseudo code is given in Algorithm 2.
(a) averaged relative errors with respect to quotation date (b) first to third quantile of time series of relative errors with respect to rolling contract Figure 5: (color online). Relative errors of model estimation. The overall relative error (averaged over all contracts and all quotation dates) is 0.661%. In (a) the averaged relative error of forward curve on each quotation date is shown. The spike in February 2016 is caused by a large downward drop of the observed forward curve, leading the model to deviate 2% on average on that date. In (b) the distribution of relative pricing errors for each nearby contract over time is given in boxplot: each whisker gives the range from mininum value to maximum value of the time series of relative errors for that contract (outliers are removed). Each green box marks the 25th to 75th percentile of the time series. The thick black line marks the median relative error. In addition to (b), the time averaged relative errors and standard deviations for each contract are given in the table below. We see that the first to sixth nearby contracts are well estimated by the model, while the seventh to tenth nearby contracts have much larger estimation errors. This occurs by construction. The real data on the back end of the forward curve are very rare and thus have a huge price uncertainty; in particular the tenth nearby contract was only available on four quotation dates on over nine years of monthly quotation data. The uncertainty of real data is captured by the parameters N j t for each j-th nearby rolling contract in the quadratic filter.  (j∆t mod 1))G w 0,1 end for end for end for Simulation study of hedging performance We aim to evaluate hedging performance by comparing the unhedged exposures with exposures when we use the locally risk-minimizing rolling hedges from Section 6 on different hedging horizons. For this, we consider different claims F (t, T ) := F (t, T, T + 1, X t ) with T = 2, . . . , 10 years. Next, we simulate M = 5000 forward curve evolutions using the estimated parameters from Table 1. For the Euler discretization we use 120 time points per year. For the hedging we use a monthly rebalancing frequency. Finally, we compare the percentage exposure if left unhedged, i.e. A visual comparison of those exposures (hedged versus unhedged) with respect to different hedging horizons is given in Figure 6. We see that the distribution of the exposure widens with increasing hedging horizon, and that the sample standard deviation and skewness go up; see the table below Figure 6. The exposure is significantly higher if left unhedged. Moreover, in all cases, the locally risk-minimizing rolling hedge significantly reduces, but does not eliminate, the variance and skew of long-term exposures. Figure 6: (color online). Density of hedged exposure (green) versus that of unhedged exposure (red) with respect to different hedging horizons. Forwards are simulated using the P-dynamics and market price of risk. In each figure, a different forward is simulated such that the time to maturity corresponds the hedging horizon: i.e. in top left figure, we simulate a forward that matures and starts delivery in two years and compute the exposure at maturity; we then compute a risk-minimizing hedge (with two years hedging horizon), the hedged exposure, and obtain the comparison. Standard deviations and skewnesses are reported in the table below. The G-matrices arising in both specifications have a zero first column, and are therefore not invertible. This is in general the case when 1 is part of the basis H(x), as G1 = 0. Moreover, if we remove the first row and column of G, the submatrix G is invertible and upper-triangular. In the following we show a straightforward way to compute t 0 e Gs ds for such G, which helps to reduce the computational effort of evaluating the pricing formula. Proof. Let F (t) denote the claimed expression for e At . One easily checks that F (t) = AF (t) and that F (0) is the identity. This implies that F (t) = e At . The expression for t 0 e As ds is easily obtained by integrating each block of F (t).