Pricing of Commodity and Energy Derivatives for Polynomial Processes

Operating in energy and commodity markets require a management of risk using derivative products such as forward and futures, as well as options on these. Many of the popular stochastic models for spot dynamics and weather variables developed from empirical studies in commodity and energy markets belong to the class of polynomial jump diffusion processes. We derive a tailor-made framework for efficient polynomial approximation of the main derivatives encountered in commodity and energy markets, encompassing a wide range of arithmetic and geometric models. Our analysis accounts for seasonality effects, delivery periods of forwards and exotic temperature forwards where the underlying “spot” is a nonlinear function of the temperature. We also include in our derivations risk management products such as spread, Asian and quanto options.


Introduction
In commodity markets, forwards and futures are traded actively in various markets and over-the-counter as a means of hedging production, controlling price risk or for pure speculation. Other derivatives, such as plain-vanilla call and put options are also widely offered for trade at organised exchanges. Such options are typically written on forward and futures contracts, but there exist also a wide zoology of tailor-made derivatives, which have payoff structures to mitigate risk exposure towards various factors beyond price, for example temperature in a power market context. To effectively apply forward, futures and derivatives in commodity and energy markets for risk management, one must have available efficient methods for quantifying the prices of these assets. This paper provides a theoretical framework for the application of polynomials to price derivatives.
The no-arbitrage theory in financial mathematics determines prices in a complete market situation as the risk-neutral expected payoff (see, e.g., [1]). In incomplete markets, which is the typical situation in commodity and energy markets, prices can be valued by the expected payoff as well, but then under a pricing measure that must be determined (the reader is referred to the works of Geman [2] and Eydeland and Wolyniec [3] for a general commodity market analysis and that or Benth, Šaltytė Benth and Koekebakker [4] for stochastic modelling and pricing). In either case, a straightforward method for computing prices is to simulate the payoff, known as a Monte Carlo approach. In a few rare cases, even formulas are available, for example the famous Black-Scholes or Black76 formulas for call and put options or Margrabe's formula for spread options. These formulas require the very restrictive assumption of a geometric Brownian motion describing the spot dynamics. Alternatively, in a more general Markovian setting for the underlying stochastic price dynamics (for the spot price of the forward/futures prices), one can resort to the accompanying partial differential equation and numerical methods to solve this. Another approach, which requires knowledge of the characteristic function of the underlying asset dynamics, is numerical integration of Fourier-based representations of the price. For an analysis of these methods in the context of energy markets, the interested reader is referred to the work of Benth, Šaltytė Benth and Kokebakker [4].
Polynomial processes provide an attractive alternative to the pricing approaches mentioned above. A polynomial process is, roughly speaking, a process for which the conditional moments become polynomials of lower order of the process. Polynomials are extremely efficient to compute on a computer, and thus such processes provide a framework for pricing. Indeed, polynomial processes have gained, since their introduction by Cuchiero [5], a lot of attention in the research literature and in various application areas, maybe foremost within finance (see the works of Filipović and Larsson [6] and Cuchiero, Keller-Ressel and Teichmann [7] for polynomial processes in finance). Most of the models used in commodity and energy markets are within the polynomial class, and pricing of forwards, futures and derivatives are feasible using methods from polynomial process theory.
Forwards and futures are settled on commodity and energy spot prices, which can be modelled by a polynomial process or some function thereof. Such functions can be exponentials, which has a series representation in terms of monomials or, as in weather markets or power, nonlinear functions of max/min type. The latter appears for example in connection with HDD and CDD weather futures traded at the CME exchange in the US.
We analyse the problem of deriving forward and futures pricing in these different situations, where we include seasonality, multi-factor models and delivery periods into our framework. The polynomial processes are of jump-type, following the recent introduction and analysis of Filipović and Larsson [8]. When we face nonlinear functions, we rely on series expansions of the function in terms of polynomials which are related to the ratio between the distribution of the polynomial process and a target distribution. Our analysis extends in this respect the ideas and derivations by Ackerer, Filipović and Pulido [9], who studied call and put options on an asset price with stochastic volatility. Recently, Ware [10] proposed and empirically estimated a polynomial model for power prices in Alberta, Canada, while Kleisinger-Yu, Komaric, Larsson and Regez [11] studied long-term forwards in the EEX power market using some specific polynomial processes. The current paper provides a general polynomial framework for forward and futures pricing in commodity and energy markets.
We also price general options on forwards and futures. These can, in some cases, be viewed as compound options. For example, options on temperature futures are options on options on temperature. This is a new application of polynomial processes in derivatives pricing, which, as we device here, relies again on the existence of polynomials for the ratio between the polynomial process distribution and a target distribution. An example is the class of multivariate Hermite polynomials, which are closely related to the Gaussian distribution, or the Laguerre polynomials, which are associated with the Gamma distribution. We also provide series expansions for the prices of spread and quanto options, exotic derivatives which are relevant in the energy markets. It is worth noticing that the pricing formulas based on polynomial analysis gives regression-type expressions for the derivatives price in terms of the current price of the underlying (or the factors thereof). Such relations can be proven useful in statistical studies of derivatives prices.
The analysis of this paper is presented as follows. In Section 2, we give a brief account on polynomial jump-diffusions, followed up by a short survey section of different stochastic dynamical models used in pricing forwards and options in commodity and energy markets. Section 4 provides a detailed analysis of polynomial processes and pricing of commodity and energy forward contracts. In Section 5, we price options on forwards based on polynomial processes. Finally, we end the paper with conclusions and outlook.

Background on Polynomial Jump-Diffusion Processes
We first recall the definition of a polynomial jump-diffusion process with values in E ⊆ R d . Our presentation is adopted from Filipović and Larsson [8].
Throughout the paper, we let (Ω, F , P) be a probability space equipped with a filtration (F t ) t≥0 satisfying the usual conditions (see [12]). For a subset E ⊆ R d , d ∈ N, we denote Pol n (E) the set of real-valued polynomials on E with order at most n ∈ N 0 := N ∪ {0}, while Pol(E) is the algebra of all polynomials on E. To be precise, a polynomial p ∈ Pol(E) is defined as the restriction p = q| E of a polynomial q ∈ Pol n (R d ), with degree deg(p) = min{deg(q) : Consider a stochastic process (X(t)) t≥0 with state space E being a special semimartingale and having the property that , the space of bounded twice continuously differentiable real-valued functions on R d . Here, with the property that f (x) = 0 for x ∈ E and that the Lévy measure has finite moments of all orders, i.e., for all x ∈ R d and n ≥ 2. Filipović and Larsson [8] referred to this as G being well defined. Following Filipović and Larsson ( [8], Definition 1), the E-valued jump-diffusion process (X(t)) t≥0 with extended generator G is said to be a polynomial process if G maps Pol n (E) into itself for each n ∈ N. By Lemma 1 in Filipović and Larsson [8], a characterisation of the polynomial processes is given as This is an "if and only if" characterisation. The attractive property of polynomial processes is that they are stable under conditional expectations of polynomials applied to the process. This property will play a key role in our analysis of derivatives pricing in energy and commodity markets. To this end, introduce a polynomial basis vector for Pol n (E), denoted by for x ∈ E and K := K(n, d) := dim Pol n (E) − 1. Sometimes, we may for notational reasons write v 0 (x) for the basis vector 1. In the basis H n,d (x), we have v i ∈ Pol n (E), i = 0, 1, . . . , K.
For p ∈ Pol n (E), let p n ∈ R K+1 be the coefficient vector such that p(x) = p n H n,d (x). Furthermore, let G n,d be the (K + 1) × (K + 1)-matrix representation of G restricted to Pol n (E), which is determined by It follows that G p(x) = p n G n,d H n,d (x).
We end this section with the main result of importance for our analysis ([8] Theorem 1): is an E-valued polynomial process. Then, for any n ∈ N 0 and p ∈ Pol n (E), it holds that for any T ≥ t, with p n ∈ R K+1 being such that p(x) = p n H n,d (x).
In this paper, we mostly deal with polynomial processes having state space E = R d . However, from time to time, we also encounter other state spaces in the discussion.

Forwards and Options on Energy and Commodities in a Polynomial Context
We review some models which are popular to apply in a commodity and energy market modelling context. The focus is on stochastic models for the spot price dynamics, including other relevant "spots" such as temperature and wind speed indices. We argue for the overwhelming evidence of polynomial processes in energy and more generally commodities markets modelling.

Commodity "Spot" Dynamics
The classical commodity spot model is given by the Schwartz model (see [13]): let where X follows an Ornstein-Uhlenbeck process Here, µ, α and σ are constants, with σ and α being positive, and B is a standard Brownian motion. The Schwartz-Smith/Gibson-Schwartz model (see [14,15]) is a twofactor extension of this which adds to X another drifted Brownian motion or Ornstein-Uhlenbeck model (see the works of Lucia and Schwartz [16] for a power market application and Prokopczuk [17] on freight futures). General multi-factor Ornstein-Uhlenbeck models which also include Lévy processes as drivers for the noise were analysed by Benth, Šaltytė-Benth and Koekebakker [4]. That is, a general spot model for commodities and energies may be where γ ∈ R d and X follows a d-dimensional Ornstein-Uhlenbeck process for µ ∈ R d , A a d × d-matrix and L a d-dimensional Lévy process. Ornstein-Uhlenbeck processes in R d belong to the class of polynomial processes. Nomikos and Soldatos [18] proposed a two-factor model of this kind for spot power prices in the NordPool electricity market, where they assume a Brownian motion and jump-driven Ornstein-Uhlenbeck process. Seasonality is added to the exponential stochastic dynamics, and the model is even extended to allow for a regime switch in the level µ of the Gaussian factor to account for the impact of reservoir filling. The Pilipović model (see [19] (Page 64)) is another class of polynomial mean-reverting two-factor process for the spot price of energies. Here, the spot price is assumed to follow the dynamics with a stochastic long-term equilibrium level L following a geometric Brownian motion and α, σ positive constants. By fixing the long-term level L to be a constant, this model coincides with what is called a GARCH diffusion by Filipović and Larsson [8] (Example 2). A typical feature of many commodity markets, in particular energy and weather markets, is seasonality. In geometric models such as (5), one typically lets t → Λ(t) be a positive-valued seasonality function and assumes A simple rewriting gives, with λ(t) := ln Λ(t). Cartea and Figueroa [20] proposed such a model for electricity spot prices in England and Wales, where the X is a univariate Ornstein-Uhlenbeck process of the form (6) with L being a Lévy prices with both Brownian motion and a compound Poisson process present. Interestingly, Mirantes, Población and Serna [21] proposed a stochastic seasonality model for the spot dynamics of Henry Hub natural gas prices. In their model, λ is assumed to follow a sum of complex Ornstein-Uhlenbeck processes in order to capture a trigonometric seasonal variation which is affected by random fluctuations. If we were to allow for complex-valued processes in our framework, this would mean that we could extend the size of the vector X by this number of complex-valued Ornstein-Uhlenbeck processes and let λ(t) = 0 in (7). So-called arithmetic models have also been proposed, taking the form Temperature (see the work of Benth and Šaltytė Benth [22] for empirical analysis and stochastic modeling) may be described by an arithmetic CARMA-process, which is given by (8) and a special case of (6) with A being a particular matrix and L is a Brownian or Lévy noise in just the last dimension and zero otherwise. For example, choosing γ = u 1 , the canonical unit vector in R d , u 1 X(t) will form a continuous autoregressive process of order d, a so-called CAR(d)-model. CAR(3)-processes have been used to model the temperature dynamics in several locations across the globe (see, e.g., the works of Härdle and Lopez-Cabrera [23] for US cities and Asian cities and Swishchuk and Cui [24] for Canadian cities). Geometric multi-factor CARMA-processes have been proposed in the context of commodity futures pricing (see the work of Paschke and Prokopczuk [25] for crude oil).
Power spot markets have a cap on the range of possible prices, which defends the introduction of a stochastic model with values in a given interval. Ware [10] proposed to model the power spot price dynamics in Alberta, Canada using a polynomial transform of the Jacobi process, with α > 0, θ ∈ [0, 1] and σ > 0. The Jacobi process takes values in the unit interval [0, 1] and is an example of a polynomial process. We remark that the Jacobi process was proposed as a stochastic volatility model by Ackerer, Filipović and Pulido [9] in financial derivatives pricing. Kleisinger-Yu et al. [11] proposed to model power spot by a quadratic polynomial of a two-factor Schwartz-Smith dynamics in a theoretical and empirical hedge study of long-term power forward contracts. They also considered a stochastic correlation between the two factors modelled as the Jacobi process.
Wind futures based on relative wind production to total capacity require a "spot" which is a process taking values in the unit interval [0, 1]. Hence, the Jacobi process is a potential candidate. However, Benth and Pircalabu [26] suggested an exponential process of the form (7) with X being the negative of a univariate Ornstein-Uhlenbeck process as in (6) driven by a subordinator Lévy process. This leads to the exponential of polynomial process with jumps, with state space [0, 1]. The Cox-Ingersoll-Ross (CIR) stochastic process was applied by Bensoussan and Brouste [27] to model 10-min wind speed data at rotor height in Wyoming, USA. The CIR stochastic dynamics is an example of a polynomial process. Benth and Rohde [28] extended the study of Bensoussan and Brouste [27] in two different directions. Considering a finite sum of squared CARMA processes, they defined a so-called CIR-CARMA model on the one hand. As an alternative to this, they defined and analysed a CARMA-process with exponential jumps. Both models fit wind speed data very well and are within the class of polynomial processes (the former in fact being a polynomial transform of a polynomial process).
There also exists some research on stochastic volatility models with polynomial structure in the context of energy markets. Kyriakou et al. [29] proposed a mean-reverting exponential model with jumps and a Heston stochastic volatility for the spot dynamics of oil prices. They calibrated their model to different refined oil futures price series in Europe and the US. Kleisinger-Yu et al. [11] modelled the correlation by a Jacobi process in a polynomial power dynamics model.
All these mentioned processes are polynomial, demonstrating from an empirical and theoretical point of view the relevance of this class of dynamical models for risk management in commodity and energy markets.

Plain-Vanilla Forward Contracts
The forward price F(t, T) at time t ≥ 0 of a plain-vanilla forward contract delivering at time T ≥ t is given as assuming that S is integrable. We suppose in this work that the risk-free interest rate is fixed to be r > 0. Thus, forwards and futures prices are identical, and we do not make any distinction between them. Notice also that we assume P to be the pricing measure, which is not necessarily equal to the market/objective probability. If the spot can be liquidly traded, the pricing measure is equal to the risk-neutral probability (i.e., the equivalent martingale measure). We notice from (10) that, when S follows a geometric model (7), we find from the power expansion of the exponential function assuming that we can interchange summation and expectation. The forward price can thus be computed by conditional expectations of polynomials of X. One can also use other polynomial expansions of the exponential function, tailor-made to the distributional properties of X. An arithmetic spot model (8) gives a very simple first-order polynomial in the conditional expectation, while Ware [10], for example, assumed the spot to be a polynomial of the Jacobi process, hence leading to a forward price of a finite sum of conditional expectations of polynomials to be computed.

Exotic Forward Contracts
Temperature forwards traded at the Chicago Mercantile Exchange (CME) are written on cooling-degree days (CDD) and heating-degree-days (HDD), which can be formulated as, respectively, call and put options on temperature spot (see Jewson and Brix [30]). The German energy exchange EEX launched in 2015 and 2016 intraday cap and floor futures, which have very similar settlement as temperature forwards (see the work of Hinderks and Wagner [31] for a discussion and analysis of these cap and floor futures). If temperature S is modelled by an arithmetic process (8), a CDD-forward price is where c is a threshold temperature, which at the CME is set to 18 • C. The max-function is obviously not polynomial, however, as presented by Ackerer et al. [9] in the context of option pricing with stochastic volatility, and discussed in a broader context in this paper below, one can represent the forward price (12) via a series of polynomials. Such series representations may be utilised in practice by a truncation of the infinite series followed by an efficient computation of the price resorting to the polynomial property of the process. We remark that, in power and weather markets, as well as in freight markets, the forwards deliver over a specified period of time T 1 ≤ T ≤ T 2 rather than at a fixed time point T. This means that the forward price becomes simply a sum (or integral) over the specified time period, For example, in the temperature market, one is summing over the daily CDD or HDD, which again is based on the average of the maximum and minimum temperature on a given day. In the power market, one is aggregating over the hourly spot prices in the delivery period. These delivery periods are typically given as weeks, months, quarters or even years.

Options in Energy and Commodities
Exchange-traded options in energy and commodity markets are typically plain-vanilla call and put options written on forward contracts. In the EEX market, e.g., call and put options are listed on forwards delivering over months and quarters and years. The price of such contracts is for a call option with strike K and exercise time τ, written on a forward delivering over the period [T 1 , T 2 ], where τ ≤ T 1 . Other markets with similar contracts include temperature derivatives at CME and options on the forward freight rate agreements (FFA) at the Baltic Exchange in London, UK.
In the oil market, e.g., the options are written on forwards with fixed-delivery time, i.e., the contract is settled on F(τ, T) with τ ≤ T. Such options are traded, e.g., at the New York Mercantile Exchange (NYMEX). NYMEX also offers trade in spread options on different blends of oil. Spreads play an important role in energy and commodity markets, for example the spark and dark spreads between, respectively, power and gas and power and coal. In the OTC-market, there exists an abundance of various options and derivatives based on such spreads but also other variations of options on several underlyings as the quanto options. Quanto options are settled on the product of two payoff functions with respect to price and a volume measure such as temperature (see, e.g., [32]).
Asian-style options on spot prices are also attractive products in energy and commodity markets. Asian options on the spot price of electricity were traded on the Nord Pool power exchange and at the Imarex shipping exchange a few decades ago (see [17,33]), and Asian-like payoff structures on spots appear in many tailor-made commodity derivatives contracts traded OTC. Fusai, Marena and Roncoroni [34] advocated Fourier-based pricing methods for discretely-monitored Asian options in corn and gas markets using a square-root (Cox-Ingersoll-Ross) stochastic process for the spot dynamics. Discrete and continuous-time Asian options in the context of crude oil were priced by an iterative numer-ical method by Kyriakou, Pouliasis and Papapostolou [35] based on the Heston stochastic volatility model and its Bates' extension where the log-price dynamics have jumps (see also [29]). The Heston along with the SABR stochastic volatility models were proposed for oil futures price dynamics by Shiraya and Takahashi [36], who derived approximation formulas by asymptotic expansions for Asian options. The models are calibrated using WTI futures American options.

Polynomial Processes and Forward Pricing
As shown in the previous section, the problem of finding the forward price in commodity and energy markets can be reduced to computing for some t ≥ s ≥ 0, with X being an R d -valued polynomial process, γ ∈ R d and λ(t) a deterministic function. We assume p ∈ Pol n (R). We see from Expression (13) that we are facing the following two problems when computing the conditional expectation. First, we shift the polynomial process by a seasonality function λ(T), which means that we need to know how the polynomials act under shifting. We express this in terms of a shift in monomials, followed by a matrix for basis change on Pol n (R). Secondly, we consider polynomials on real-valued affine transformations of multi-dimensional polynomial processes. This requires an understanding of how one transforms p(γ x) for a real-valued polynomial p into a polynomial on R d , with γ, x ∈ R d . We can assign a linear transformation between the polynomial bases in Pol n (R) and Pol n (R d ).
In the next two lemmas, we spell this out.
The first lemma is a simple consequence of the binomial formula: . , x n ) be a vector of monomials on R up to order n ∈ N. Then, where Λ n ∈ R (n+1)×(n+1) is a lower triangular matrix with elements (Λ n (λ)) ij = ( i−1 j−1 )λ i−j for i = 1, 2, . . . n + 1 and j ≤ i.
Proof. Let m ≤ n. Then, by the binomial formula, This yields the result.
If we are given a basis H n (x) := H n,1 (x) of polynomials in Pol n (R), where h k (x) ∈ Pol n (R), then one can find an invertible matrix C n ∈ R (n+1)×(n+1) such that In forward pricing, we consider polynomials which are shifted by the seasonal function λ(t), p(λ(t) + x). Hence, if p n ∈ R n+1 is such that p(x) = p n H n (x), for H n (x) as in Lemma above, then we find We have the following technical result on changing from univariate to multivariate polynomial bases: Proof. From (15), we have that H n (γ x) = C n M n (γ x). Moreover, by the multinomial formula, it holds for 1 ≤ k ≤ n, . Therefore, we can find a matrix Γ n,d such that M n (γ x) = Γ n,d H n,d (x) and the lemma follows.
Typically, a seasonality function λ(t) may be a finite series of sine and cosine functions. Since, for example, the pair of functions (sin(kt), cos(kt)) satisfies a two-dimensional first order linear system of ordinary differential equations, we see that truncated series of trigonometric functions fits into the framework of polynomial processes. This was pointed out and discussed by Filipović and Willems [37] in their study of dividend derivatives. As mentioned in Section 3, complex-valued Ornstein-Uhlenbeck processes were proposed by Mirantes, Población and Serna [21] to model stochastic seasonality. Their model can be viewed as an extension of the seasonality dynamics by Filipović and Willems [37] with additive (complex-valued) noise. More generally, one can allow for polynomial complex-valued processes of polynomial type to model stochastic seasonality. Such an approach would in our framework imply that we ignore the seasonality function λ, i.e., let λ(t) = 0 and extend the dimensionality of the vector-valued polynomial process X (as well as allowing for more generally complex-valued polynomial processes). The price we would pay for interpreting the seasonality λ as a polynomial process would be that the dimensionality d of X increases and hence the dimension K + 1 of Pol n (R d ).

Plain-Vanilla Forward Prices
We find the following result, collecting the notation introduced in this section: where p n (λ) is given in (16), Γ n,d in Lemma 2 and G n,d is the polynomial transition matrix defined in (2).
Proof. First note that a polynomial process has finite conditional moments of all orders (see [8] (Theorem 1)). Hence, the conditional expectation is well-defined. From Lemma 2, we have for γ, x ∈ R d and λ ∈ R, Hence, The result follows from the polynomial process property of X (see Theorem 1).
If the spot dynamics follows an arithmetic model (8), then we find the forward price (10) as where Thus, we have a simple relationship when the monomials are chosen as the basis of Proposition 1 is a simple extension of the formulas by Kleisinger-Yu et al. [11] who focussed on the case of p ∈ Pol 2 (R d ). Ware [10] also discussed such formulas in his polynomial approach to power spot models.
Let us discuss the forward price dynamics in Proposition 1 from an empirical viewpoint. To put our discussion into context, we note that Ware [10] proposed, among other models, a one-factor polynomial diffusion process combined with a fifth-order polynomial as a spot dynamics, i.e., S(t) = p(X(t)) with X ∈ R and p ∈ Pol 5 (R). Here, we ignore seasonality for simplicity in the discussion. Then, using τ := T − t, the time to maturity, we get where θ := Γ 5,1 p 5 ∈ R 6 . Using the monomials as basis, we have . . , 5. Hence, G 5,1 will be a lower triangular matrix with zero in the first row and diagonal elements ka 1 + 1 2 k(k − 1)σ 2 for k = 1, . . . , 5. (Note that first diagonal element is zero, the second is a 1 , the next is 2a 1 + σ 2 , etc. to the last which is 5a 1 + 10σ 2 .) The distinct eigenvalues of the matrix then becomes λ 0 = 0 and λ k = ka 1 + 1 2 k(k − 1)σ 2 for k = 1, . . . , 5. We can find a basis of R 6 of eigenvectors w 0 , w 1 , . . . , w 5 , where w 0 can be chosen to be the first canonical basis vector of R 6 with 1 in first coordinate and zeros otherwise. From this, we find In commodity markets, one typically expects forward prices to be flat in the long end of the curve as long as spot prices are expected to possess some stationarity properties. It is evident from above that the forward prices for large τ's tend to p 5 w 0 whenever the eigenvalues λ k are negative. Hence, we obtain a flat forward curve in the long end when ka 1 + 1 2 k(k − 1)σ 2 < 0 for k = 1, . . . , 5. For example, if σ 2 = 0, then this is achieved when a 1 < 0. On the other hand, in the Jacobi model suggested by Ware [10], σ 1 = −1, which again implies that a 1 < 0. This is indeed the case in his model.
By an application of Ito's formula, one furthermore observes that the volatility of f (t, τ) will satisfy the Samuelson effect, as the volatility will be determined by the diffusion term from X(t) and scaled by exponentials e λ k τ . Whenever λ k < 0, these exponentials will tend to 1 when τ tends to zero, which gives a Samuelson effect as the forward volatility converges to the spot volatility in this case.
It is also worth noticing that the sum of exponentials in (17) gives rise to several humps in the forward term structure, that is, the curve τ → f (t, τ) may have several local maxima and minima according to the values of the eigenvalues. A hump shape behaviour of the forward curve is reasonable from an economic viewpoint, indicating differences in risk preferences of the traders along the forward curve. We refer to the work of Benth, Šaltytė Benth and Koekebakker [4] for a discussion of the various stylised facts of forward curves in power and commodity markets.
The analysis above can be generalised to arbitrary polynomials p, as well as also extending the dimension of the polynomial process beyond d = 1.
If the spot follows a geometric model (7), then, from the Taylor series expansion of the forward price in (11), one chooses the basis for Pol n (R) to be the monomials, i.e., H n (x) = M n (x). The basis for Pol n (R d ) is arbitrary selected. We find: Proposition 2. Suppose the commodity spot dynamics is given by S(t) = exp(λ(t) + γ X(t)) as in (7). Assume that exp(|γ X(T)|) ∈ L 1 (P). Then, the forward price in (11) is given by where u n+1 is the n + 1-canonical unit vector in R n+1 (i.e., the vector with 1 in coordinate n + 1 and zero otherwise), Γ n,d in Lemma 2, Λ n (λ) in Lemma 1 and G n,d in (2).
Proof. From the exponential integrability condition on γ X(T), it follows from monotone convergence (see [38] (Theorem 2.15)) that Hence, in particular, we find that E[|γ X(T)| n ] < ∞ for every n ∈ N. As ((γ X(T)) n ) n∈N therefore is a sequence in L 1 (P) and ∑ ∞ n=0 1 n! E[|γ X(T)| n ] < ∞, it follows from dominated convergence theorem (see [38] (Theorem 2.25)) that By Lemma 1, we derive Next, by Lemma 2 followed by the polynomial property of X yield, The result follows.
The Proposition above allows for stating the forward price as where x → f n (t, T; x) ∈ Pol n (R d ). In practical computations, one would of course truncate the sum. One could also make use of the (truncated) sum in a regression study, where one empirically could reveal the structure of the f n 's by regressing observed forward prices against polynomials of X. Such a study requires knowledge of the state of X(t), which can be recovered from the spot prices. Such recovery may involve stochastic filtering if d > 1.
The exponential integrability condition exp(|γ X(T)|) ∈ L 1 (P) in Proposition 2 is rather restrictive. Geometric Brownian motion, e.g., does not satisfy this condition for any γ. From Example 2 of Filipović and Larsson [8], the GARCH diffusion process for constants κ, θ ∈ R + is a polynomial process with ergodic solution being inverse Gaussian distributed with shape parameter 2 and 1/θ as scale. In the invariant case, X does not have a finite variance and hence not being exponentially integrable either. By letting θ be a geometric Brownian motion, the GARCH diffusion becomes the Pilipović model (see Pilipović [19]) briefly discussed in Section 3. This model will not in general be exponentially integrable. Ornstein-Uhlenbeck processes driven by compound Poisson processes having exponential jumps will have a gamma distributed stationary solution. This will also not be exponentially integrable, except under restrictive conditions on the parameters.

Exotic Forward Prices
Next, consider CDD-forwards on temperature (or a floor electricity forward) with price given in (12). We need some preparatory material on polynomial expansions of call and put payoff functions. For n ∈ N 0 , let ξ n (x) denote the nth Hermite polynomial (known as the "probabilistic" Hermite polynomial) defined as where is the density of the standard normal distribution function. We notice that ξ 0 (x) = 1. Further, define for n ∈ N 0 with the usual convention that 0! = 1. It is known that (e n ) n∈N 0 is an ONB of the Hilbert space L 2 w := L 2 (R, w(x)dx). Considering f (x) = max(x − c, 0), we readily find that f ∈ L 2 w as f (x) = (x − c) for x > c and zero for x < c and w integrates any polynomial. Moreover, for any f ∈ L 2 w , we find that with Y ∼ N (0, 1), a standard normal random variable. Moreover, from elementary functional analysis, we have The following simple result holds: Hence, the first claim follows. For the second claim, the Cauchy-Schwarz inequality implies Invoking the first claim proves the Lemma.
We extend the previous result to more general random variables Y in the next lemma: Lemma 4. Suppose Y is a random variable with probability density φ Y satisfying φ Y (y) ≤ Cw a,b 2 (y) for a.e. y ∈ R, where C ≥ 1 is a constant and w a,b 2 is the normal density function with mean a and variance is b Proof. By the Cauchy-Schwarz inequality, Consider the positive function As 0 < b < 1, it holds that 1 − b −2 < 0 and thus u has a maximum value on R. It follows that w .
Since f ∈ L 2 w and f m is its truncation in the basis representation, the result follows after passing to the limit.
We recall that f (x) = max(x − c, 0) satisfies the requirement that f ∈ L 2 w . Moreover, CARMA(p, q)-processes driven by Brownian motion are normally distributed, and hence the above result applies with φ Y being a normal distribution with variance less than 1. We also recall from Ackerer et al. [9] that the Jacobi volatility process has a distribution which is absolutely continuous with respect to the normal distribution. In addition, rather than using the Taylor series representation exp(x) = ∑ ∞ k=0 1 k! x k , we may use the Hermite polynomials as series expansion for the exponential function since obviously exp(x) ∈ L 2 w . The condition that the variance b is strictly less than one is very restrictive. However, one can overcome this by a change in the Hermite basis or by appropriate rescaling the function f . We provide a thorough discussion of this in Section 4.3, where we take a more general perspective. For the moment, we note that Ackerer et al. [9] used an affine transform of the Hermite polynomials as basis.

Remark 3.
We notice that the condition on the probability density of Y in Lemma 4 implies that the distribution Φ Y (dy) := φ Y (y)dy is absolutely continuous with respect to w a,b 2 (y)dy. By assuming instead that the distribution of Y, Φ Y is dominated by that of w a,b 2 (y)dy, i.e., Φ(dy) ≤ Cw a,b 2 (y)dy in the sense Φ Y (U) ≤ C U w a,b 2 (y)dy for every Borel set U ⊂ R, we find that there exists an a.e. non-negative Radon-Nikodym density Y ∈ L 1 (R, w a,b 2 (y)dy). In this case, we can define a probability density φ Y (y) := Y (y)w a,b 2 (y). However, then, φ Y (y) ≤ Cw a,b 2 (y), a.e. y ∈ R because, if this is not the case there exists a measurable set U with strictly positive mass such φ Y (y) > Cw a,b 2 (y), y ∈ U, which implies that Φ Y (U) = U φ Y (y)dy > C U w a,b 2 (y)dy being a contradiction. Further notice that the constant C must be greater than or equal to 1 simply because we have distributions with total mass 1 on both sides of the bound.
In the next subsection, we take a more general perspective where the distribution of the polynomial process does not need to be bounded by a Gaussian but other suitable classes of distributions for which we can associate polynomials. At the current stage of our exposition, we focus on the Gaussian case as this is the most relevant in connection with temperature forwards, where the underlying dynamics have empirical evidence for being normally distributed (recall discussions in Section 3).
We next show a polynomial expression for the CDD-temperature forward price: to this end, choose the basis H n (x) = (e 0 (x), e 1 (x), . . . , e n (x)) for Pol n (R), where (e n (x)) n∈N 0 are the normalised Hermite polynomials defined in (20). For Pol n (R d ), we fix an arbitrary basis H n,d (x). (8), where X is a d-dimensional polynomial process. Assume that the random variable γ X(T) has an F t -conditional probability density which is bounded by a normal density w a,b 2 as in Lemma 4. Then, the forward price in (12) is given by

Proposition 3. Suppose the commodity spot dynamics is given by S(t) = λ(t) + γ X(t) as in
where u n+1 is the n + 1-canonical unit vector in R n+1 (i.e., the vector with 1 in coordinate n + 1 and zero otherwise), C n is given in (15), Γ n,d in Lemma 2, Λ n (λ) in Lemma 1 and G n,d in (2).

Proof.
We find that f (x) = max(x − c, 0) ∈ L 2 w , and therefore From the condition on the density of γ X(T) given F t , it holds from Lemma 4 that the conditional expectation is well-defined as integrability holds, and that we can commute sum and conditional expectation. That is, By the translation of the monomial basis in Lemma 1, we find that Furthermore, invoking Lemma 2 gives Finally, applying the polynomial property of X, we find The result follows.
We remark in passing that the above Proposition could also have been developed for other functions f than the one appearing for the CDD-temperature forwards. In fact, any function f ∈ L 2 w would do, with the only difference that the coefficients f n appearing in the expression for F in Proposition 3 would change (as they depend explicitly on f , of course). For example, using f (x) = exp(x), which defines a function in L 2 w , Proposition 3 provides an alternative forward price series expression to Proposition 2 for geometric spot price models.
To efficiently compute the CDD-temperature forward price by exploiting the polynomial structure of X, we truncate the infinite sum. All the matrices and vectors involved are explicitly given, except the coefficient functions f n . We recall these to be defined as for Y being standard normally distributed. We can compute these coefficients once for a given function f ∈ L 2 w , as they are independent of the polynomial process X.

Remark 4.
In the market for temperature forwards, the contracts are settled over a pre-specified period of time. In that case, a CDD-temperature forward is after appealing to the Fubini-Tonelli theorem to commute sums.
If temperature follows a CARMA-dynamics driven by a Brownian motion, then the dimension d will indicate the autoregressive order. Moreover, γ X(t) will be Gaussian and X a d-dimensional Ornstein-Uhlenbeck process, and thus the conditions in Proposition 3 hold. We recall that the temperature dynamics is conveniently modelled by a CAR(3)process (see [22]), which means that d = 3 and γ = u 1 , the canonical unit vector in R 3 with 1 in first coordinate and zero otherwise.
Interestingly, Asian options are closely related to the above CARMA-situation by the following argument: consider an Asian option with payoff max(T −1 T 0 γ X(s)ds − c, 0) at exercise time T. Here, X is a d-dimensional polynomial process and γ ∈ R d , with the spot price being S(t) = γ X(t) (we ignore seasonality here in this short discussion). Let now X be the process in R d+1 defined as X = ( X, X d+1 ), where dX d+1 (t) = γ X(t)dt.
It follows that X is a polynomial process, and we have that the Asian option payoff can be written as T −1 max(γ X(T) − Tc, 0) with γ = u d+1 , the canonical unit vector in R d+1 with one in the last coordinate and zero otherwise. In particular, assuming that X is a multivariate Gaussian Ornstein-Uhlenbeck process, we will have that X is a Gaussian Ornstein-Uhlenbeck process and we find ourselves in a situation which is closely resembling the forward price of a CDD-temperature contract analysed above.

A General Polynomial Approach to Forward Pricing
In this subsection, we take a general perspective on forward pricing, providing a unifying expression for the forward price in markets with a polynomially based "spot"process. The approach requires some additional conditions on the polynomial process, but on the other hand gives an attractive treatment of options on forwards, a topic which is analysed in Section 5.
Suppose that the "spot price" dynamics is given by for some measurable function g : R d × R → R, seasonality function λ and X being a d-dimensional polynomial process. Examples of relevance can be for γ ∈ R d , λ ∈ R and ξ being one of the following functions: ξ(x) = x (arithmetic spot model), ξ(x) = exp(x) (geometric spot model) or ξ(x) = max(x − c, 0) (spot for an exotic forward such as temperature futures). Our aim is to compute the forward price, defined as To achieve this goal, we employ a multivariate generalisation of the space L 2 w along with an integrability assumption on the conditional probability distribution of X(T) given F t . In fact, without losing any generality for practical purposes in commodity and energy markets, we assume that X is also a Markovian process. (Note that non-Markovian polynomial jump diffusion processes exist, see., e.g., [8] (Page 71).) Let us start by introducing a multi-dimensional generalisation of the space L 2 w . To this end, let ρ be a probability density function on R d , and for d ∈ N, denote by L 2 ρ,d the Hilbert space of real-valued functions on R d for which Assume further that there exists an ONB for L 2 ρ,d of polynomials, given by (v n d ) n d ∈N d 0 using the multi-index notation n d := (n 1 , . . . , n d ). We use the notation |n d | = n 1 + · · · + n d for the order of the multi-index, where it is supposed that v n d ∈ Pol |n d | (R d ). Furthermore, (v n d ) |n d |≤N is a basis of polynomials of order N, which we use as H N,d (x). Ranking the basis functions v n d according to their polynomial order is convenient and natural when doing approximations in practical applications of this theory. Next, denote by φ(x, dy; t, T) the transition probability distribution on R d of X(T) given X(t) = x. Following Filipović and Larsson [8] (Sect. 7), introduce the likelihood ratio as the function (x, y; t, T) such that φ(x, dy; t, T) = (x, y; t, T)ρ(y)dy (23) We assume that such a likelihood ratio of φ with respect to ρ exists. In the next theorem, we state a general series representation in terms of polynomials for the forward price along with a computationally convenient truncation.

Theorem 2.
Assume that g(·; λ(T)) ∈ L 2 ρ,d and (x, ·; t, T) ∈ L 2 ρ,d for any 0 ≤ t ≤ T < ∞ and x ∈ R, where g and are defined, respectively, in (21) and (23). Then, we have that Here, . We recall K(n, d) = dim Pol n (R d ) − 1, and G given in (2). We find from the assumptions g(·; λ(T)), (x, ·; t, T) ∈ L 2 ρ,d that Therefore, by Parseval's identity, are the coefficients in the ONB representation of g(·; λ(T)) and (x, ·; t, T) in L 2 ρ,d , respectively. Tracing back the definitions, we find By assumption on the polynomial basis, v n d ∈ Pol |n d | (R d ). Thus, there is a vector with length equal to the dimension of Pol . We then conclude the desired form, after appealing to the polynomial property of X. This proves the representation of F(t, T).
Define for each N ∈ N the approximation We observe that N (x, ·; t, T) is nothing but y → (x, y; t, T) projected down on the finite dimensional subspace of L 2 ρ,d spanned by (v n d ) n d ∈N d 0 ,|n d |≤N . This gives us F N (t, T). Notice that by Parseval's identity, From the very definitions of f and f N , we find by the Cauchy-Schwarz inequality and Parseval's identity, In conclusion, f N (x; t, T) → f (x; t, T) for every x ∈ R d when N → ∞. The proof is complete.
We notice that the dependency on the seasonality component is merged into the coefficients g n d (λ(T)) and is as such not material in the analysis above. We include it simply because seasonality is present in relevant models, and we prefer to have it explicit. Additionally, it highlights a difference with the other polynomial expansions which we present in this section. Furthermore, we observe that we may compute the coefficients g n d (λ(T)) by numerical integration methods, for example Gaussian quadrature or Monte Carlo simulation. Indeed, we have where Z is a d-dimensional random variable with probability density ρ.

Remark 5.
In the case X is not a polynomial process, we see by inspection of the proof of Theorem 2 that we still have an interesting representation of the forward price in terms of the polynomial moments of X(T) conditional on X(t). Indeed, removing the polynomial property, we see that all conclusions in the theorem holds, except that which will not be explicit in terms of polynomials of x. Of course, we still need the regularity assumptions of g and the likelihood ratio to hold. We can approximate the forward prices by polynomial moments of the process up to a certain order.
The main assumption in our general approach to forward pricing is the existence of a density ρ admitting a polynomial basis for L 2 ρ,d , such that there is likelihood ratio function being an element of this space. This problem is classical, and has a long history in probability and physics, where we refer to the works of Asmussen, Goffard and Laub [39] and Eggers [40] for some recent applications and studies. One thinks of ρ as the reference measure, and for a target distribution φ the goal is to have a Gram-Charlier series with efficiently computable polynomials. Following the discussion in Asmussen et al. [39], if ρ in Dimension 1 has all moments finite, there exists an orthogonal sequence of polynomials, which, moreover, defines a basis in L 2 ρ,1 if ρ has finite exponential moment. One can easily build up multivariate reference measures in general dimensions d by tensorising. For example, we may define ρ(x) := w ⊗d (x) := w(x 1 ) · · · w(x d ). This will provide a d-dimensional version of the space L 2 w based on Hermite polynomials. In Rahman [41], a general multivariate basis of Hermite polynomials are defined appealing to the Rodrigues formula, i.e., based on the derivatives of the multivariate Gaussian distribution function with mean zero and general covariance. A special case of this, choosing the covariance matrix to be the identity, leads back to the definition of ρ(x) = w ⊗d (x). Another example could be a reference measure in d = 1 defined by the gamma-distribution (see the work of Asmussen et al. [39], where Laguerre polynomials appear).
We discuss the case ρ(x) = w ⊗d (x) in some more detail. We start by introducing a multi-dimensional version of L 2 w . To this end, for d ∈ N, denote by L 2 w,d the Hilbert space of real-valued functions on R d for which where we recall w ⊗d (x) := w(x 1 ) · · · w(x d ). An ONB for L 2 w,d is given by (e n d ) n d ∈N d 0 using the multi-index notation n d := (n 1 , . . . , n d ), and e n d (x) := e n 1 ⊗ · · · ⊗ e n d (x) = e n 1 (x 1 ) · · · e n d (x d ). Here, we recall (e n ) n∈N 0 to be an ONB of L 2 w . Let us look at some particular cases of the function g in the case of L 2 w,d , which are of relevance to commodity and energy markets. First, in an exponential spot price model, we have g(x; λ) = exp(λ + γ x) for γ, x ∈ R d and λ ∈ R. Since R y → exp(2γ i y)w(y) is integrable, we see that g(·; λ) ∈ L 2 w,d . We compute the coefficients g n d (λ(T)): In the arithmetic case, the spot takes the form with again x, γ ∈ R d and λ ∈ R. Since R y → (γ i y) 2 w(y) is integrable and g(·, λ) ∈ L 2 w,d . We find the coefficients in this arithmetic case as follows: = λ(T) e 0 , e n 1 w · · · e 0 , e n d w γ i e 0 , e n 1 w · · · e 0 , e n i−1 w e 1 , e n i w e 0 , e n i+1 w · · · e 0 , e n d w . Now, for n j ≥ 1 it holds that R e n j (y)w(y)dy = e 0 , e n j w = 0. Thus, the only non-zero terms in the sum above are those where n d = (0, . . . , 0, 1, 0, . . . , 0), where 1 is appearing in coordinate i, in which case g n d (λ(T)) = γ i . All other n d will give g n d (λ(T)) = 0, except n d = (0, . . . , 0) where we get g n d (λ(T)) = λ(T). This is a very unsurprising result, of course.
Next, let us consider the case of a CDD-temperature forward or a floor electricity forward, for which we have that g(x; λ) = ξ(λ + γ x) for ξ(z) = max(z − c, 0). As the max-function grows at most linearly, we have 0 ≤ g(x; λ) ≤ |λ| + |γ x|, and it follows that g(·; λ) ∈ L 2 w,d . We may represent the coefficients as for Z ∼ N (0, I) with I being the d × d identity matrix. By iterated conditional expectation, conditioning on Z i , i = 1, . . . , d − 1, we can define for R being a standard normal random variable and iteratively backwards i = d − 1, d − 2, . . . , 1, yielding g n d (λ(T)) = g 1 . A more detailed discussion of the condition on the likelihood ratio is in place. We first recall from Ackerer et al. [9] that the Jacobi volatility model has a likelihood ratio with respect to the Gaussian density which satisfies the condition of square integrability. Hence, the Jacobi volatility model allows itself to a series expansion in terms of the Hermite polynomials, as conducted in detail by Ackerer et al. [9]. Many of the interesting polynomial models are such that X(T)| F t is Gaussian. For example, we have the two-factor models of Lucia and Schwartz or CARMA-models driven by Brownian motion. This results in a conditional distribution function φ(x, dy; t, T) being a Gaussian distribution with mean µ ∈ R d and covariance matrix V ∈ R d×d . Here, we collapse the notation to make our discussion more transparent. Hence, we find that the likelihood ratio is Here, I is the d × d identity matrix. It is evident that the function y → (x, y; t, T) ∈ L 2 w,d if and only if −(V −1 − I) − 1 2 I < 0, or, equivalently, V < 2I. This is not always true, as we can have two-factor models with independent Brownian motions having variance each strictly bigger that 2. Then, V is a diagonal matrix with variances on the diagonal which is dominating 2I, and the required integrability of the likelihood ratio fails.
In such cases, we can re-scale the polynomial process X. To this end, let C be some d × d matrix such that CVC < 2I. If we have available such a matrix, we can define a new stochastic process Y(t) := CX(t). Since any matrix transform of a polynomial process again is a polynomial process, Y is a polynomial process. If further C is invertible, then X(t) = C −1 Y(t), and we have g(x; λ) = g(C −1 y; λ).
In Theorem 2, we assume that g(C −1 ·; λ(T)) ∈ L 2 w,d . Furthermore, for the polynomial process Y we find that the likelihood ratio function is (as a matrix transformation of the Gaussian variable X(T)| F t ), Hence, we have that y → (x, y; t, T) ∈ L 2 w,d whenever C is such that CVC < 2I. Here is an example of a re-scaling: Let X(T)| F t be bivariate Gaussian with covariance matrix where σ i , σ 2 are two strictly positive constants (being the marginal standard deviations) and −1 < ρ < 1 (which is the correlation). For example, this is the situation with the two-factor model of Lucia and Schwartz, or a CARMA-model in d = 2. Observe that a diagonalisation of V is given by Let, for some positive constant c < 2 Then, we find since c < 2. In conclusion, we find a scaling C of the original polynomial process, for which the covariance matrix can be dominated by 2 times the identity. Then, the likelihood ratio has the desired integrability, but we must adjust slightly the integrability condition on g. For most interesting functions g, this is not any added restriction, as for example the cases considered above.
Rather than re-scaling, we could use the multivariate Hermite polynomials introduced by Rahman [41] for a sufficiently big covariance matrix. In the above case of re-scaling, we need to have some knowledge of the matrix C before doing the computations. However, the advantage then is that one can simply apply the standard one-dimensional Hermite polynomials as basis. An approach using the multivariate Hermite polynomials requires knowledge of a suitable covariance matrix, which essentially is such that the target density φ can be dominated by this. The multivariate Hermite polynomials can then be derived, a task that must be tailor-made to the choice of matrix.
Next, we consider a case with non-Gaussian reference probability ρ, focussing on the one-dimensional situation. We recall that factor models with Ornstein-Uhlenbeck dynamics driven by jump processes are relevant for power price and wind speed modelling. In particular, Ornstein-Uhlenbeck processes with exponential jump processes leading to invariant Γ-distributions are applied (recall discussion from [4,26,28] in Section 4, say). In addition, we have CIR-processes as a model for wind speeds as we recall from [27]. The CIR-process is skewed χ 2 -distributed at each time instant, a distribution which is closely related to the Γ-distribution. Let now ξ be the density of the Γ-distribution with scale r > 0 and shape m > 0, given as ξ(y) = 1 Γ(r)m r y r−1 e −y/m . Suppose we have a target distribution φ which behaves as y s−1 , s > 0 for y ∼ 0 and e −y/k for y ∼ ∞, then the likelihood ratio will be y s−1 /y r−1 close to zero, and exp(−(1/k − 1/m)y) for y ∼ ∞. However, integrating the square of the likelihood function against ξ, yields finiteness whenever s > r/2 and k < 2m. Such conditions were found by Asmussen et al. [39] as well. Thus, tuning the m to be sufficiently large and r to be sufficiently small, we can obtain a target Γ-distribution such that the likelihood ratio is square integrable with respect to ξ. Furthermore, as is well-known, the basis of orthogonal polynomials for L 2 ξ := L 2 (R + , ξ(y)dy) is the generalised Laguerre polynomials (see [42]). If we have a two-factor model, with one Gaussian and one exponential jump Ornstein-Uhlenbeck process, we can consider the tensorised space L 2 ρ,2 with ρ(x) = w ⊗ ξ(x 1 , x 2 ) = w(x 1 )ξ(x 2 ) and the canonically generated polynomials from the respective marginal densities.

Pricing of Options on Forwards
This section is concerned with the problem of pricing options on forwards in the framework of polynomial processes.

Options on Plain-Vanilla Forwards
Consider a European option written on a plan-vanilla forward with payoff ζ (F(τ, T)) at time τ ≤ T for a payoff function ζ, with the forward price F(t, T) given as in Proposition 1; that is, for some d, n ∈ N, we have where h(t, T) = p n (λ(T)) Γ n,d exp(G n,d (T − t)).
In the formulation of Proposition 1, H n,d (x) is some basis of the nth order polynomials on R d . It is convenient, however, for the purpose of option pricing, to turn to the polynomial ρ,d , as used in Section 4.3 above. We fix H n,d (x) to be the basis (v n d (x)) |n d |≤n from now on. Furthermore, we suppose that X is a Markovian process.
Let the price of the option (with risk-free interest rate set to zero) at time t ≤ τ be where we assume ζ(F(τ, T)) ∈ L 1 (P). The following result is essentially a repetition of Theorem 2 and is therefore formulated as a corollary.
ρ,d with h given as in (25). Let F be given in (24) with X being a polynomial process on R d for which the likelihood ratio function defined in (23) satisfies R d y → (x, y; t, T) ∈ L 2 ρ,d . Then, we have that P N (t, τ, T) → P(t, τ, T) (pointwise) when N → ∞, where P(t, τ, T) is the option price in (26) with representation while for any N ∈ N, Here, Proof. The proof is identical to the argument of Theorem 2, but now using P(t, τ, Notice that τ now plays the role of T in the proof of Theorem 2, and that ζ is g. We also observe that ζ(F(τ, T)) ∈ L 1 (P) under the assumptions on g and X.
If ζ(z) = max(z − K, 0), the payoff function of a call option, we find that 0 ≤ ζ(h(t, T) H n,d (x)) ≤ K + h(t, T) H n,d (x) , using the notation · for the Euclidean 2-norm on, e.g., R d . This shows readily that ζ(h(t, T) ·) ∈ L 2 ρ,d as long as L 2 ρ,d supports polynomials of degree n. This is the case if we choose the space L 2 w,d . Another popular class of derivatives in commodity and energy markets is spread options, which we discuss next.
Assume we have two commodity forwards, with respective forward prices F i (t, T), i = 1, 2 which are given by (24) for two different functions h i (t, T), i = 1, 2. Thus, the dynamics of both forward prices are driven by the same polynomial process X. The spread option payoff at time τ is ζ(F 1 (τ, T 1 ) − F 2 (τ, T 2 )), with τ ≤ min(T 1 , T 2 ). Hence, we have potentially two different maturities of the forwards. Notice also that the order n and dimensionality d of H n,d (x) are the same for both forwards, which is not a lack of generality as we can extend the dimensionality of both canonically, if necessary.
ρ,d , we can apply Corollary 1 with h(t, T 1 , , which satisfies the regularity condition for at least L 2 w,d . It is remarked in passing that we can also treat quanto options in a similar manner. A quanto option pays ζ 1 (F 1 (τ, T 1 ))ζ 2 (F 2 (τ, T 2 )) for two payoff functions ζ 1 and ζ 2 . These may be of the form of two calls, or a call and put, or two puts. Defining ζ(H n,d (x); t, T) := puts us again in the situation where Corollary 1 may be applied.

Remark 6.
To include forward contracts with delivery period into the pricing framework is straightforward. Since we have we can simply redefine the meaning of h(t, T) in order to apply Corollary 1.
From a computational point of view, it is important to notice that the polynomial representation of the price P in Corollary 1 is split into coefficients ζ n d (τ, T) and n d (x; t, τ). The latter family of coefficients, n d (x; t, τ), is only dependent on the underlying stochastic model and the choice of polynomial basis, and thus can be computed irrespective of the option in question. As noted by Ackerer et al. [9], these coefficients may be relatively costly to compute, but one can do this once and apply the coefficients for the numerical evaluation of different options. The option payoff is encoded in the parameters ζ n d (τ, T).
At this instance, we make some further comments on the numerical implementation and performance of polynomial pricing found in the literature. The already mentioned paper by Ackerer et al. [9] presents three numerical case studies for the Jacobi stochastic volatility model. Pricing a call option for given realistic model parameters for the stock market, they show that the price error is accurate in two decimal points for N chosen between 10 and 15, where the "exact" price is determined using Monte Carlo simulations. They also considered a forward start call option and an Asian option, where dimension is increased due to the payoff structure. In these two cases, the level N must be increased to above 15, which requires rather many coefficients to be computed. In [9], various approaches for the computation of the so-called Fourier coefficients (being ζ n d (τ, T) in our notation) are discussed, including recurrence relations based on properties of Hermite polynomials and Gaussian cubature integration. Moreover, there are references to efficient methods for the computation of matrix exponentials, which we encounter in the numerical computation of n d (x; t, τ) for rather high-dimensional matrices G |n d |,d . Further numerical studies on polynomial volatility models and option pricing can be found in the work of Ackerer and Filipović [43].
Related numerical studies based on polynomial processes are found in the works of Kleisinger-Yu et al. [11] and Benth and Lavagnini [44]. Kleisinger-Yu et al. [11] computed the quadratic risk minimising hedging strategy of long-term delivery forwards in the power markets. For polynomial models, this entails in rather explicit expressions which are efficiently computed for low-dimensional polynomial processes as stochastic models for the forward price dynamics. Benth and Lavagnini [44] aimed at the computation of correlators, which occur for example in the iterative definition of discretely-sampled path-dependent options or in a series expansion of Fourier-based pricing of options in stochastic volatility models. Polynomial processes allow for explicit matrix representations of the correlators, and numerical case studies show a good performance even for high-dimensional situations compared with Monte Carlo methods.

A General Polynomial Approach to Option Pricing
In this subsection, we price options on forwards which have price expressions developed in the context of Section 4.3. To this end, for d ∈ N, recall the Hilbert space L 2 ρ,d introduced in the previous section. Consider the "doubled" space L 2 ρ 2 ,2d which is the Hilbert space of real-valued functions on R 2d for which to be an ONB of L 2 ρ,d . We also recall the notation g(·; λ(T)) from (21) for the forward price (22), where X is a polynomial process in R d which in addition is assumed to be Markovian. We further recall that we denote by φ(x, dy; t, T) the probability distribution function of X(T) given X(t) = x ∈ R d for t ≥ T. Under suitable conditions, Theorem 2 provides an expression and approximation of F. Consider an option written on the forward with exercise time τ ≤ T and payoff function ζ(F(τ, T)) for some function ζ : R → R such that ζ(F(τ, T)) ∈ L 2 (P). The arbitrage-free option price at time t ≤ τ (for risk-free interest rate set to zero) is given by P(t, τ, T) as defined in (26). Theorem 3. Assume g(·; λ(T)) ∈ L 2 ρ,d and suppose that the likelihood function (x, y; t, T) defined in (23) satisfies R d × R d (x, y) → (x, y; t, T) ∈ L 2 ρ 2 ,2d for any 0 ≤ t ≤ T < ∞. If ζ : R → R is Lipschitz continuous and of linear growth, then and, moreover, Here, k d (x; t, τ) is defined in Theorem 2, Proof. By assumption (·, ·; t, T) ∈ L 2 ρ 2 ,2d and therefore y → (x, y; t, T) ∈ L 2 ρ,d , a.e., x ∈ R d . Hence, since g(·; λ(T)) ∈ L 2 w,d we find a polynomial expression of F as given in Theorem 2 along with an approximation F N . From the proof of Theorem 2, we also recall the function along with its representation and series expansions found in the proof of that result.
To apply Theorem 3 in practice, we need to compute the coefficients (ζ • f N ) k d for all k d ∈ N d 0 such that |k d | ≤ K. We have that f N is again a truncated sum, but the coefficients g n d of this is available from the computation of forward prices (or approximations thereof). This representation can be used to calculate (ζ • f N ) k d , which require numerical integration or possibly Monte Carlo simulation by drawing from a random variable distributed according to ρ.
Theorem 3 provides us with an approximation of call and put option prices on forwards that again have option-like structures, i.e., an approximation of compound options. As noted above, the options in the temperature market can be viewed as a class of compound options, although we recall that temperature forwards have a measurement (delivery) period which is not allowed for by a direct use of Theorem 3. However, we can easily adjust the above arguments to account for a measurement (delivery) period in the forward contract. In this case, we have that the option price in (26) is given by To apply the arguments of Theorem 3, we must assume that ∑ T 2 T=T 1 g(·; λ(T)) ∈ L 2 ρ,d . If g(·; λ(T)) ∈ L 2 ρ,d for any T, then this condition holds as L 2 ρ,d is a vector space. Tracing the steps in the proof of Theorem 3 results in P(t, τ, T 1 , On the other hand, as long as ζ is of bounded linear growth and Lipschitz continuous, P N,K (t, τ, T 1 , T 2 ) → P(t, τ, T 1 , T 2 ) with P N,K (t, τ, T 1 , and where we find that In fact, f N (x; τ, T 1 , T 2 ) = ∑ T 2 T=T 1 f N (x; τ, T), thus, to approximate the option price on a temperature HDD or CDD futures, we simply add up a finite sum of terms f N (x; τ, T) over T when computing the coefficients (ζ • f N ) k d rather than using only one f N (x; τ, T).
The Greeks, or sensitivities with respect to different parameters of the option price, are important for hedging purposes. The so-called "delta", the derivative of the option price with respect to the present value of the underlying asset, is readily computed in terms of the derivatives of the polynomials k d (x; τ, T) with respect to x. This lowers the polynomial order of k d (x; τ, T) by one, and we have available an explicit series representations and approximation under appropriate conditions. Other Greeks, for example the sensitivity with respect to the volatility, can also be represented as derivatives of these polynomials as the volatility will be inherent in the specification of X. This further emphasises the attractiveness of polynomial models and their expressions of option prices.
From a computational perspective, several challenging issues arise. Focussing on temperature futures options, we note above that CARMA-processes are suitable for modelling the temperature dynamics. This calls for higher-dimensional models, i.e., it is natural to have the dimension d to be around 3 or even higher. If we were to specify the seasonality also as a polynomial process, we would reach a much higher dimensionality of the underlying polynomial dynamics. As noted above, the numerical studies of Ackerer et al. [9] indicate that one needs the order of re-scaled Hermite polynomials to be about N = 10 for call options on a stochastic volatility model. In our situation, which is of greater dimensionality, we would expect even higher order to reach a satisfactory convergence. This at the level of approximating the forward price, where, additionally, we need to aggregate up the n d (x; τ, T) over the delivery period. Then, on the next level, we again must approximate a call option, where we also may need high-order polynomials. On the other hand, we know from the theory that the sums are converging and thus the terms must tend to zero despite the involved polynomials occurring in n d (x; τ, T). The convergence speed for these expressions should be further analysed. These are challenging computational problems which we leave for future studies.

Conclusions and Outlook
We derive polynomial series representations for forward prices and options on forwards based on polynomial models of the spot dynamics in commodity markets. Commodity markets have special features such as seasonality and delivery period, as well as exotic payment structures for forwards found in, e.g., power and temperature markets. In a review of the literature on modelling of price risk in energy and commodity, we present many different polynomial models, which we use as motivation and foundation for further derivatives pricing. We also note some empirical facts on polynomial models and issues and challenges concerning numerical applications.
When considering prices based on nonlinear payoff structures, such as those appearing in temperature futures and options on forwards, the successful derivation of a polynomial series expansion rests on the relationship between the generating probability density of a class of polynomials, the probability distribution of the polynomial process and the squareintegrability of likelihood function. It is an interesting area of further research to gain a deeper understanding of this connection. It is also interesting to analyse further numerical implementations of some of the derivatives analysed in this paper, where dimensionality becomes a challenge. Acknowledgments: Two anonymous referees are thanked for their careful reading and constructive criticism improving the presentation of the paper.

Conflicts of Interest:
The author declares no conflict of interest.