Pricing of Commodity and Energy Derivatives for Polynomial Processes

Abstract

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.

1. 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.

2. Background on Polynomial Jump-Diffusion Processes

We first recall the definition of a polynomial jump-diffusion process with values inE⊆^Rd . Our presentation is adopted from Filipović and Larsson [8].

Throughout the paper, we let(Ω,F,P)be a probability space equipped with a filtration_(Ft)t≥0 satisfying the usual conditions (see [12]). For a subsetE⊆^Rd,d∈N, we denote_Poln(E)the set of real-valued polynomials on E with order at mostn∈_N0:=N∪{0}, whilePol(E)is the algebra of all polynomials on E. To be precise, a polynomialp∈Pol(E)is defined as the restriction_p=q|Eof a polynomialq∈_Poln(^Rd), with degreedeg(p)=min{deg(q):p=q_|E,q=Pol(^Rd)}.

Consider a stochastic process_(X(t))t≥0with state space E being a special semimartingale and having the property that

f(X(t))−f(X(0))−_∫0tGf(X(s))ds

is a local martingale for allf∈_Cb2(^Rd), the space of bounded twice continuously differentiable real-valued functions on^Rd. Here,

Gf(x)=a^(x)⊤∇f(x)+12Trσ(x)^∇2f(x)+_∫^Rd (f(x+z)−f(x)−^z⊤∇f(x))ℓ(x,dz)

for measurable functionsa:Rd→Rdandσ:Rd→S+d,S+dis the space of positive-semidefinite symmetricd×dmatrices. Furthermore,ℓ(x,·)is the Lévy jump measure onRd, with propertiesℓ(x,{0})=0and∫Rd min(|z|,|z|2)ℓ(x,dz)<∞for allx∈Rd. Additionally, we assume thatGf(x)=0forx∈Efor allf∈Pol(Rd)with the property thatf(x)=0forx∈Eand that the Lévy measure has finite moments of all orders, i.e.,

∫Rd |z|nℓ(x,dz)<∞

for allx∈^Rdandn≥2 . Filipović and Larsson [8] referred to this asG being well defined. Following Filipović and Larsson ([8], Definition 1), the E-valued jump-diffusion process_(X(t))t≥0with extended generatorGis said to be a polynomial process ifGmaps_Poln(E)into itself for eachn∈N . By Lemma 1 in Filipović and Larsson [8], a characterisation of the polynomial processes is given as

a∈Pol1(E)σ+∫Rd zz⊤ℓ(·,dz)∈Pol2(E)∫Rd zαℓ(·,dz)∈Pol|α|(E)

for all multi-indicesα=(_α1,…,_αd)∈_N0dwhere|α|:=_α1+⋯+_αd≥3. Moreover,^zα=_z1α1 ⋯_zdαd . 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_Poln(E), denoted by

_Hn,d(x)=^{(1,_v1(x),…,_vK(x))⊤}

forx∈EandK:=K(n,d):=dim_Poln(E)−1. Sometimes, we may for notational reasons write_v0(x)for the basis vector 1. In the basis_Hn,d(x), we have_vi∈_Poln(E),i=0,1,…,K.

Remark 1.

The dimension of_Poln(^R2)is easily shown to satisfydim_Poln(^R2)=n+1+dim_Poln−1(^R2). Sincedim_Pol0(^R2)=1, we find that

dim_Poln(^R2)=1+∑i=1n+1i=12(n+2)(n+1).

In general, we havedim_Poln(^Rd)=n+dn (see, e.g., [8]).

Forp∈_Poln(E), let_pn∈^RK+1be the coefficient vector such thatp(x)=_{pn⊤ Hn,d}(x). Furthermore, let_Gn,dbe the(K+1)×(K+1)-matrix representation ofGrestricted to_Poln(E), which is determined by

G_Hn,d(x)=_{Gn,d Hn,d}(x).

It follows that

Gp(x)=_{pn⊤ Gn,d Hn,d}(x).

We end this section with the main result of importance for our analysis ([8] Theorem 1):

Theorem 1.

Assume_(X(t))t≥0is an E-valued polynomial process. Then, for anyn∈_N0andp∈_Poln(E), it holds that

Ep(X(T))|_Ft=_pn⊤exp(_Gn,d(T−t))_Hn,d(X(t))

for anyT≥t, with_pn∈^RK+1being such thatp(x)=_{pn⊤ Hn,d}(x).

In this paper, we mostly deal with polynomial processes having state spaceE=^Rd. However, from time to time, we also encounter other state spaces in the discussion.

3. 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. 3.1. Commodity “Spot” Dynamics

The classical commodity spot model is given by the Schwartz model (see [13]): let

S(t)=exp(X(t))

where X follows an Ornstein–Uhlenbeck process

dX(t)=(μ−αX(t))dt+σdB(t).

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 two-factor 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

S(t)=exp(^γ⊤X(t))

whereγ∈^Rdand X follows a d-dimensional Ornstein–Uhlenbeck process

dX(t)=(μ−AX(t))dt+dL(t).

forμ∈^Rd, A ad×d-matrix and L a d-dimensional Lévy process. Ornstein–Uhlenbeck processes in^Rd 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

dS(t)=α(L(t)−S(t))dt+σS(t)dB(t),

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 letst↦Λ(t)be a positive-valued seasonality function and assumes

S(t)=Λ(t)exp(^γ⊤X(t)).

A simple rewriting gives,

S(t)=exp(λ(t)+^γ⊤X(t)).

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

S(t)=λ(t)+^γ⊤X(t)

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γ=_u1, the canonical unit vector in^Rd,_u1⊤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,

dX(t)=α(μ−X(t))dt+σX(t)(1−X(t))dB(t)

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. 3.2. Plain-Vanilla Forward Contracts

The forward priceF(t,T)at timet≥0of a plain-vanilla forward contract delivering at timeT≥tis given as

F(t,T)=E[S(T)|_Ft],

assuming that S is integrable. We suppose in this work that the risk-free interest rate is fixed to ber>0. Thus, forwards and futures prices are identical, and we do not make any distinction between them. Notice also that we assumePto 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

F(t,T)=∑n=0∞1n!E[^{(λ(T)+γ⊤X(T))n}|_Ft]

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.

3.3. 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

F(t,T)=Emax(S(T)−c,0)|_Ft

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_T1≤T≤_T2rather than at a fixed time point T. This means that the forward price becomes simply a sum (or integral) over the specified time period,

F(t,_T1,_T2)=∑T=_T1T2F(t,T)

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. 3.4. 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

C(t,τ):=^{e−r(τ−t)}E[max(F(τ,_T1,_T2)−K,0)|_Ft]

for a call option with strike K and exercise timeτ, written on a forward delivering over the period[_T1,_T2], whereτ≤_T1. 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 onF(τ,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 numerical 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.

4. 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

Ep(λ(t)+^γ⊤X(t))|_Fs

for somet≥s≥0, with X being an^Rd-valued polynomial process,γ∈^Rdandλ(t)a deterministic function. We assumep∈_Poln(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_Poln(R). Secondly, we consider polynomials on real-valued affine transformations of multi-dimensional polynomial processes. This requires an understanding of how one transformsp(^γ⊤x)for a real-valued polynomial p into a polynomial on^Rd, withγ,x∈^Rd. We can assign a linear transformation between the polynomial bases in_Poln(R)and_Poln(^Rd). In the next two lemmas, we spell this out.

The first lemma is a simple consequence of the binomial formula:

Lemma 1.

Let_Mn(x)=^{(1,x,x2,…,xn)⊤}be a vector of monomials onRup to ordern∈N. Then,

_Mn(λ+x)=_Λn(λ)_Mn(x)

where_Λn∈^{R(n+1)×(n+1)}is a lower triangular matrix with elements_(Λn(λ))ij=i−1j−1^λi−jfori=1,2,…n+1andj≤i.

Proof.

Letm≤n. Then, by the binomial formula,

^(λ+x)m=∑k=0mmk^{λm−k xk}

This yields the result.  □

If we are given a basis_Hn(x):=_Hn,1(x)of polynomials in_Poln(R),

_Hn(x):=^{(_h0(x):=1,_h1(x),…,_hn(x))⊤}

where_hk(x)∈_Poln(R), then one can find an invertible matrix_Cn∈^{R(n+1)×(n+1)}such that

_Hn(x)=_{Cn Mn}(x)

In forward pricing, we consider polynomials which are shifted by the seasonal functionλ(t),p(λ(t)+x). Hence, if_pn∈^Rn+1is such thatp(x)=_{pn⊤ Hn}(x), for_Hn(x)as in Lemma above, then we find

p(λ+x)=_{pn⊤ Hn}(λ+x)=_{pn⊤ Cn Mn}(λ+x)=_{pn⊤ Cn Λn}(λ)_Mn(x)=_pn ^(λ)⊤ _Hn(x)

with

_pn ^(λ)⊤:=_{pn⊤ Cn Λn}(λ)_Cn−1.

We have the following technical result on changing from univariate to multivariate polynomial bases:

Lemma 2.

There exists an(n+1)×(K+1)-dimensional matrix_Γn,dsuch that

_Hn(^γ⊤x)=_{Γn,d Hn,d}(x)

for anyγ,x∈^Rd. Here,K+1is the dimension of_Poln(^Rd).

Proof.

From (15), we have that_Hn(^γ⊤x)=_{Cn Mn}(^γ⊤x). Moreover, by the multinomial formula, it holds for1≤k≤n,

^(γ⊤x)k=^{(_{γ1 x1}+_{γ2 x2}+⋯+_{γd xd})k}=∑_k1+⋯+_kd=kk_k1,_k2,…,_{kdΠi=1d γiki xiki} ,

where

k_k1,_k2,…,_kd=k!_k1!⋯_kd!

is the multinomial coefficient. Hence,^(γ⊤x)k∈_Polk(^Rd)and_Mn(^γ⊤x)is ann+1-dimensional vector of polynomials in_Poln(^Rd). Therefore, we can find a matrix_Γ˜n,dsuch that_Mn(^γ⊤x)=_{Γ˜n,d Hn,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)=0and 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 dimensionK+1of_Poln(^Rd).

4.1. Plain-Vanilla Forward Prices We find the following result, collecting the notation introduced in this section:

Proposition 1.

Letp∈_Poln(^Rd). ForT≥t≥0, it holds,

Ep(λ(t)+^γ⊤X(T))|_Ft=_pn ^(λ(T))⊤ _Γn,dexp(_Gn,d(T−t))_Hn,d(X(t))

where_pn(λ) is given in (16),_Γn,din Lemma 2 and_Gn,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∈^Rdandλ∈R,

p(λ+^γ⊤x)=p^(λ)⊤ _Hn(^γ⊤x)=_pn ^(λ)⊤ _{Γn,d Hn,d}(x).

Hence,

Ep(λ(T)+^γ⊤X(T))|_Ft=_pn ^(λ(T))⊤ _Γn,dE_Hn,d(X(T))|_Ft

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

F(t,T)=_p1 ^(λ(T))⊤ _Γ1,dexp(_G1,d(T−t))_H1,d(X(t))

where_p1:=^{(_p1,0,_p1,1)⊤}is the vector such thatx=_{p1,0 h0}(x)+_{p1,1 h1}(x). Of course, if we choose_H1(x)=_M1(x), then_p1=^(0,1)⊤,_C1=_I2, the2×2-identity matrix and

_Λ1(λ)=10λ1

Then,

_p1 ^(λ)⊤=(λ,1).

Furthermore,_Γ1,dis in this case the matrix mapping_M1(^γ⊤x)=^{(1,_{γ1 x1}+⋯+_{γd xd})⊤}into_H1,d(x)=^{(_v0(x),…,_vd(x))⊤}, with_vi(x)being polynomials of order 1 on^Rd. If_v0(x)=1and_vi(x)=_xifori=1,…,d, then

_Γ1,d=10⋯00_γ1⋯_γd

Thus, we have a simple relationship when the monomials are chosen as the basis of_H1,d(x). Indeed, the forward price is

F(t,T)=(λ(T),_γ1,…,_γd)exp(_G1,d(T−t))1_X1(t)··_Xd(t)

Remark 2.

Proposition 1 is a simple extension of the formulas by Kleisinger-Yu et al. [11] who focussed on the case ofp∈_Pol2(^Rd) . 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))withX∈Randp∈_Pol5(R). Here, we ignore seasonality for simplicity in the discussion. Then, usingτ:=T−t, the time to maturity, we get

f(t,τ):=F(t,t+τ)=^θ⊤exp(_G5,1τ)_H5,1(X(t))

whereθ:=_{Γ5,1⊤ p5}∈^R6. Using the monomials as basis, we have_H5,1(x)=^{(1,x,x2,x3,x4,x5)⊤}. Moreover, as X is polynomial, we have that its drifta(x):=_a0+_a1xis in_Pol1(R)and diffusionσ(x)=_σ0+_σ1x+_σ2 ^x2is in_Pol2(R). If_vk(x)=^xk,k=0,…,5, it is easy to see thatG_v0(x)=0,G_v1(x)=a(x)andG_vk(x)=a(x)k_vk−1(x)+12^σ2(x)k(k−1)_vk−1(x)fork=2,…,5. Hence,_G5,1will be a lower triangular matrix with zero in the first row and diagonal elementsk_a1+12k(k−1)_σ2fork=1,…,5. (Note that first diagonal element is zero, the second is_a1, the next is2_a1+_σ2, etc. to the last which is5_a1+10_σ2.) The distinct eigenvalues of the matrix then becomes_λ0=0and_λk=k_a1+12k(k−1)_σ2fork=1,…,5. We can find a basis of^R6of eigenvectors_w0,_w1,…,_w5, where_w0can be chosen to be the first canonical basis vector of^R6with 1 in first coordinate and zeros otherwise. From this, we find

exp(_G5,1τ)_H5,1(X(t))=_w0+∑k=15^e_λkτ _{wk⊤ H5,1}(X(t))_wk

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_{p5⊤ w0}whenever the eigenvalues_λkare negative. Hence, we obtain a flat forward curve in the long end whenk_a1+12k(k−1)_σ2<0fork=1,…,5. For example, if_σ2=0, then this is achieved when_a1<0 . On the other hand, in the Jacobi model suggested by Ware [10],_σ1=−1, which again implies that_a1<0. This is indeed the case in his model.

By an application of Ito’s formula, one furthermore observes that the volatility off(t,τ)will satisfy the Samuelson effect, as the volatility will be determined by the diffusion term fromX(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 beyondd=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_Poln(R)to be the monomials, i.e.,_Hn(x)=_Mn(x). The basis for_Poln(^Rd)is arbitrary selected. We find:

Proposition 2.

Suppose the commodity spot dynamics is given byS(t)=exp(λ(t)+^γ⊤X(t)) as in (7). Assume thatexp(|^γ⊤X(T)|)∈^L1(P) . Then, the forward price in (11) is given by

F(t,T)=∑n=0∞1n!_{un+1⊤ Λn}(λ(T))_Γn,dexp(_Gn,d(T−t))_Hn,d(X(t)),

where_un+1is then+1-canonical unit vector in^Rn+1(i.e., the vector with 1 in coordinaten+1and zero otherwise),_Γn,din Lemma 2,_Λn(λ)in Lemma 1 and_Gn,d in (2).

Proof.

From the exponential integrability condition on^γ⊤X(T) , it follows from monotone convergence (see [38] (Theorem 2.15)) that

∑n=0∞1n!E[|λ(T)+^γ⊤X(T)^|n]=E[∑n=0∞1n!|λ(T)+^γ⊤X(T)^|n]=E[exp(|λ(T)+^γ⊤X(T)|]≤exp(|Λ(T)|)E[exp(|^γ⊤X(T)|)]<∞

Hence, in particular, we find thatE[|^γ⊤X(T)^|n]<∞for everyn∈N. As_{(^(γ⊤X(T))n)n∈N}therefore is a sequence in^L1(P)and_∑n=0∞1n!E[|^γ⊤X(T)^|n]<∞ , it follows from dominated convergence theorem (see [38] (Theorem 2.25)) that

E[exp(λ(T)+^γ⊤X(T))|_Ft]=∑n=0∞1n!E[^{(λ(T)+γ⊤X(T))n}|_Ft].

By Lemma 1, we derive

E[^{(λ(T)+γ⊤X(T))n}|_Ft]=_un+1⊤E[_Mn(λ(T)+^γ⊤X(T))|_Ft]=_{un+1⊤ Λn}(λ(T))E[_Mn(^γ⊤X(T))|_Ft].

Next, by Lemma 2 followed by the polynomial property of X yield,

E[_Mn(^γ⊤X(T))|_Ft]=_Γn,dE[_Hn,d(X(T))|_Ft]=_Γn,dexp(_Gn,d(T−t))_Hn,d(X(t))

The result follows.  □

The Proposition above allows for stating the forward price as

F(t,T)=∑n=0∞1n!_fn(t,T;X(t))

wherex↦_fn(t,T;x)∈_Poln(^Rd). 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_fn’s by regressing observed forward prices against polynomials of X. Such a study requires knowledge of the state ofX(t), which can be recovered from the spot prices. Such recovery may involve stochastic filtering ifd>1.

The exponential integrability conditionexp(|^γ⊤X(T)|)∈^L1(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

dX(t)=κ(θ−X(t))dt+2κX(t)dB(t)

for constantsκ,θ∈_R+is a polynomial process with ergodic solution being inverse Gaussian distributed with shape parameter 2 and1/θ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.

4.2. 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. Forn∈_N0, let_ξn(x)denote the nth Hermite polynomial (known as the “probabilistic” Hermite polynomial) defined as

_ξn(x)=⁽⁻¹⁾ⁿ1w(x)^dnd^xnw(x)

where

w(x)=12π^e−x2/2

is the density of the standard normal distribution function. We notice that_ξ0(x)=1. Further, define

_en(x)=_ξn(x)n!,

forn∈_N0with the usual convention that0!=1. It is known that_(en)n∈N0is an ONB of the Hilbert space_Lw2:=^L2(R,w(x)dx). Consideringf(x)=max(x−c,0), we readily find thatf∈_Lw2asf(x)=(x−c)forx>cand zero forx<cand w integrates any polynomial. Moreover, for anyf∈_Lw2, we find that

_|f|w2:=_∫R ^f2(x)w(x)dx=E[^f2(Y)]

withY∼N(0,1), a standard normal random variable. Moreover, from elementary functional analysis, we have

f(x)=∑n=0∞_∫Rf(y)_en(y)w(y)dy_en(x).

The following simple result holds:

Lemma 3.

Supposef∈_Lw2and denote by^fm(x):=_{∑n=0m fn en}(x). Then, forY∼N(0,1),_{(^fm(Y))m∈N}converges tof(Y)in^L2(P). Moreover,

E[f(Y)]=∑n=0∞_fnE[_en(Y)]

Proof.

Obviously,^fm∈_Lw2andf−^fm→0in_Lw2asm→∞. The latter means

0=limm→∞_|f−^fm|w2=limm→∞_∫R ^{(f(x)−fm(x))2}w(x)dx=limm→∞E[^{(f(Y)−fm(Y))2}]

Hence, the first claim follows. For the second claim, the Cauchy–Schwarz inequality implies

|E[f(Y)]−∑n=0m_fnE[_en(Y)]^|2=|E[f(Y)−^fm(Y)]^|2≤E[|f(Y)−^fm(Y)^|2].

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_ϕYsatisfying_ϕY(y)≤C_wa,^b2(y)fora.e.y∈R, whereC≥1is a constant and_wa,^b2is the normal density function with mean a and variance is^b2with0<b<1. Iff∈_Lw2then

E[f(Y)]=∑n=0∞_fnE[_en(Y)]

Proof.

By the Cauchy–Schwarz inequality,

|E[f(Y)]−∑n=0m_fnE[_en(Y)]^|2=|E[f(Y)−^fm(Y)]^|2≤E[|f(Y)−^fm(Y)^|2]=_∫R ^{(f(y)−fm(y))2} _ϕY(y)dy≤C_∫R ^{(f(y)−fm(y))2} _wa,^b2(y)dy=C_∫R ^{(f(y)−fm(y))2}_wa,^b2(y)w(y)w(y)dy.

Consider the positive function

u(y):=_wa,^b2(y)w(y)=exp12^y2(1−^b−2)+a^b2y−^a22^b2

As0<b<1, it holds that1−^b−2<0and thus u has a maximum value onR. It follows that

|E[f(Y)]−∑n=0m_fnE[_en(Y)]^|2≤Csupy∈Ru(y)_|f−^fm|w2.

Sincef∈_Lw2and^fmis its truncation in the basis representation, the result follows after passing to the limit.  □

We recall thatf(x)=max(x−c,0)satisfies the requirement thatf∈_Lw2. 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 representationexp(x)=_∑k=0∞1k!^xk, we may use the Hermite polynomials as series expansion for the exponential function since obviouslyexp(x)∈_Lw2.

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)dyis absolutely continuous with respect to_wa,^b2(y)dy. By assuming instead that the distribution of Y,_ΦYis dominated by that of_wa,^b2(y)dy, i.e.,Φ(dy)≤C_wa,^b2(y)dyin the sense_ΦY(U)≤C_{∫U wa,^b2}(y)dyfor every Borel setU⊂R, we find that there exists ana.e.non-negative Radon–Nikodym density_ℓY∈^L1(R,_wa,^b2(y)dy). In this case, we can define a probability density_ϕY(y):=_ℓY(y)_wa,^b2(y). However, then,_ϕY(y)≤C_wa,^b2(y),a.e.y∈Rbecause, if this is not the case there exists a measurable set U with strictly positive mass such_ϕY(y)>C_wa,^b2(y),y∈U, which implies that_ΦY(U)=_{∫U ϕY}(y)dy>C_{∫U wa,^b2}(y)dybeing 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_Hn(x)=^{(_e0(x),_e1(x),…,_en(x))⊤}for_Poln(R), where_{(en(x))n∈N0} are the normalised Hermite polynomials defined in (20). For_Poln(^Rd), we fix an arbitrary basis_Hn,d(x).

Proposition 3.

Suppose the commodity spot dynamics is given byS(t)=λ(t)+^γ⊤X(t) as in (8), where X is a d-dimensional polynomial process. Assume that the random variable^γ⊤X(T)has an_Ft-conditional probability density which is bounded by a normal density_wa,^b2 as in Lemma 4. Then, the forward price in (12) is given by

F(t,T)=∑n=0∞_{fn un+1⊤ Cn Λn}(λ(T))_{Cn−1 Γn,d}exp(_Gn,d(T−t))_Hn,d(X(t)),

where_un+1is then+1-canonical unit vector in^Rn+1(i.e., the vector with 1 in coordinaten+1and zero otherwise),_Cn is given in (15),_Γn,din Lemma 2,_Λn(λ)in Lemma 1 and_Gn,d in (2).

Proof.

We find thatf(x)=max(x−c,0)∈_Lw2, and therefore

f(x)=∑n=0∞_{fn en}(x)

From the condition on the density of^γ⊤X(T)given_Ft, 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,

F(t,T)=E[max(λ(T)+^γ⊤X(T)−c,0)|_Ft]=E[f(λ(T)+^γ⊤X(T))|_Ft]=∑n=0∞_fnE[_en(λ(T)+^γ⊤X(T))|_Ft]=∑n=0∞_{fn un+1⊤}E[_Hn(λ(T)+^γ⊤X(T))|_Ft]

By the translation of the monomial basis in Lemma 1, we find that

E[_Hn(λ(T)+^γ⊤X(T))|_Ft]=_{Cn Λn}(λ(T))_Cn−1E[_Hn(^γ⊤X(T))|_Ft]

Furthermore, invoking Lemma 2 gives

E[_Hn(^γ⊤X(T))|_Ft]=_Γn,dE[_Hn,d(X(T))|_Ft].

Finally, applying the polynomial property of X, we find

E[_Hn,d(X(T))|_Ft]=exp(_Gn,d(T−t))_Hn,d(X(t))

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 functionf∈_Lw2would do, with the only difference that the coefficients_fnappearing in the expression for F in Proposition 3 would change (as they depend explicitly on f, of course). For example, usingf(x)=exp(x), which defines a function in_Lw2, 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_fn. We recall these to be defined as

_fn:=_∫Rmax(x−c,0)_en(x)w(x)dx=E[max(Y−c,0)_en(Y)]

for Y being standard normally distributed. We can compute these coefficients once for a given functionf∈_Lw2, 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

F(t,_T1,_T2)=∑n=0∞_{fn un+1⊤ Cn}∑T=_T1T2Λn(λ(T))_{Cn−1 Γn,d}exp(_Gn,d(T−t))_Hn,d(X(t))

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 thatd=3andγ=_u1, the canonical unit vector in^R3with 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 payoffmax(^T−1 _∫0T ^γ˜⊤X˜(s)ds−c,0)at exercise time T. Here,X˜is a d-dimensional polynomial process andγ˜∈^Rd, with the spot price beingS(t)=^γ˜⊤X˜(t)(we ignore seasonality here in this short discussion). Let now X be the process in^Rd+1defined asX=(X˜,_Xd+1), where

d_Xd+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−1max(^γ⊤X(T)−Tc,0)withγ=_ud+1, the canonical unit vector in^Rd+1with one in the last coordinate and zero otherwise. In particular, assuming thatX˜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.

4.3. 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

S(t)=g(X(t);λ(t))

for some measurable functiong:^Rd×R→R, seasonality functionλand X being a d-dimensional polynomial process. Examples of relevance can be

g(x;λ):=ξ(λ+^γ⊤x)

forγ∈^Rd,λ∈Randξ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

F(t,T)=E[g(X(T);λ(T))|_Ft]

To achieve this goal, we employ a multivariate generalisation of the space_Lw2along with an integrability assumption on the conditional probability distribution ofX(T)given_Ft . 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_Lw2. To this end, letρbe a probability density function on^Rd, and ford∈N, denote by_Lρ,d2the Hilbert space of real-valued functions on^Rdfor which

|g|ρ,d2:=∫Rd g2(x)ρ(x)dx

with inner product

_{⟨g,h⟩ρ,d}:=_∫^Rd g(x)h(x)ρ(x)dx

Assume further that there exists an ONB for_Lρ,d2of polynomials, given by_{(vnd )nd∈N0d}using the multi-index notation_nd:=(_n1,…,_nd). We use the notation|_nd|=_n1+⋯+_ndfor the order of the multi-index, where it is supposed that_vnd ∈_Pol|nd|(^Rd). Furthermore,_{(vnd )|nd|≤N}is a basis of polynomials of order N, which we use as_HN,d(x). Ranking the basis functions_vnd 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^RdofX(T)givenX(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

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 thatg(·;λ(T))∈_Lρ,d2andℓ(x,·;t,T)∈_Lρ,d2for any0≤t≤T<∞andx∈R , where g and ℓ are defined, respectively, in (21) and (23). Then, we have that

^FN(t,T)→F(t,T)

(pointwise) whenN→∞ , where F is the forward price in (22) with representation

F(t,T)=∑n∈_N0dgnd (λ(T))_ℓnd (X(t);t,T)

while for anyN∈N,

^FN(t,T)=∑n∈_N0d,|_nd|≤N_gnd (λ(T))_ℓnd (X(t);t,T).

Here,

_gnd (λ(T))=_∫^Rd g(y;λ(T))_vnd (y)ρ(y)dy

and

_ℓnd (x;t,T)=_und⊤exp(_G|nd|,d(T−t))_H|nd|,d(x)

with_und ∈^{RK(|_nd|,d)+1}is such that_vnd (x)=_{und⊤ H|nd|,d}(x). We recallK(n,d)=dim_Poln(^Rd)−1 , and G given in (2).

Proof.

Notice first by the Markovian property of X thatF(t,T)=f(X(t);t,T), where

f(x;t,T):=E[g(X(T);λ(T))|X(t)=x].

We find from the assumptionsg(·;λ(T)),ℓ(x,·;t,T)∈_Lρ,d2that

f(x;t,T)=_∫^Rd g(y;λ(T))ϕ(x,dy;t,T)=_∫^Rd g(y;λ(T))ℓ(x,y;t,T)ρ(y)dy=_{⟨g(·;λ(T)),ℓ(x,·;t,T)⟩ρ,d}

Therefore, by Parseval’s identity,

f(x;t,T)=∑_nd∈_{N0dgnd ℓnd} (x;t,T)

where

_gnd (λ(T)):=_{⟨g(·,λ(T)),vnd ⟩ρ,d}

and

_ℓnd (x;t,T):=_{⟨ℓ(x,·;t,T),vnd ⟩ρ,d}

are the coefficients in the ONB representation ofg(·;λ(T))andℓ(x,·;t,T)in_Lρ,d2, respectively. Tracing back the definitions, we find

_ℓnd (x;t,T)=_{∫^Rd vnd} (y)ϕ(x,dy;t,T)=E[_vnd (X(T))|X(t)=x].

By assumption on the polynomial basis,_vnd ∈_Pol|nd|(^Rd). Thus, there is a vector with length equal to the dimension of_Pol|nd|(^Rd)such that_vnd (x)=_{und⊤ H|nd|,d}(x). We then conclude the desired form,

_ℓnd (x;t,T)=_und⊤E[_H|nd|,d(X(T))|X(t)=x]=_und⊤exp(_G|nd|,d(T−t))_H|nd|,d(x)

after appealing to the polynomial property of X. This proves the representation ofF(t,T).

Define for eachN∈Nthe approximation

^fN(x;t,T):=_{⟨g(·;λ(T)),^ℓN(x,·;t,T)⟩ρ,d}=∑_nd∈_N0d,|_nd|≤N_gnd (λ(T))_ℓnd (x;t,T)

with the notation

^ℓN(x,·;t,T):=∑_nd∈_N0d,|_nd|≤N_ℓnd (x;t,T)_vnd (·).

We observe that^ℓN(x,·;t,T)is nothing buty↦ℓ(x,y;t,T)projected down on the finite dimensional subspace of_Lρ,d2spanned by_{(vnd )nd∈N0d,|nd|≤N}. This gives us^FN(t,T).

Notice that by Parseval’s identity,

_{∞>|ℓ(x,·;t,T)|ρ,d2}=∑_nd∈_N0dℓnd2(x;t,T)

From the very definitions of f and^fN, we find by the Cauchy–Schwarz inequality and Parseval’s identity,

|f(x;τ,T)−^{fN (x;τ,T)|2}=|_{⟨g,ℓ(x,·;τ,T)−^ℓN(x,·;τ,T)⟩ρ,d} ^|2≤_{|g|ρ,d2 |ℓ(x,·;τ,T)−^ℓN(x,·;τ,T)|ρ,d2}≤_|g|ρ,d2∑_nd∈_N0d,|_nd|>N_ℓnd2(x;t,T)

In conclusion,^fN(x;t,T)→f(x;t,T)for everyx∈^RdwhenN→∞. The proof is complete.  □

We notice that the dependency on the seasonality component is merged into the coefficients_gnd (λ(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_gnd (λ(T))by numerical integration methods, for example Gaussian quadrature or Monte Carlo simulation. Indeed, we have

_gnd (λ(T))=E[g(Z;λ(T))_vnd (Z)]

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 ofX(T)conditional onX(t). Indeed, removing the polynomial property, we see that all conclusions in the theorem holds, except that

_ℓnd (x;t,T)=E[_vnd (X(T))|X(t)=x],

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ρ,d2 , 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ρ,12ifρ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(_x1)⋯w(_xd). This will provide a d-dimensional version of the space_Lw2 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 ind=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_Lw2. To this end, ford∈N, denote by_Lw,d2the Hilbert space of real-valued functions on^Rdfor which

|g|w,d2:=∫Rd g2(x)w⊗d(x)dx

with inner product

_⟨g,h⟩w,d:=_∫^Rd g(x)h(x)^w⊗d(x)dx

where we recall^w⊗d(x):=w(_x1)⋯w(_xd). An ONB for_Lw,d2is given by_{(end )nd∈N0d}using the multi-index notation_nd:=(_n1,…,_nd), and_end (x):=_en1 ⊗⋯⊗_end (x)=_en1 (_x1)⋯_end (_xd). Here, we recall_(en)n∈N0to be an ONB of_Lw2.

Let us look at some particular cases of the function g in the case of_Lw,d2, which are of relevance to commodity and energy markets. First, in an exponential spot price model, we have

g(x;λ)=exp(λ+^γ⊤x)

forγ,x∈^Rdandλ∈R. SinceR∋y↦exp(2_γiy)w(y)is integrable, we see thatg(·;λ)∈_Lw,d2. We compute the coefficients_gnd (λ(T)):

gnd (λ(T))=∫Rd exp(λ(T)+γ⊤x)end (x)w⊗d(x)dx=eλ(T) ∫Rd eγ1 x1⋯eγd xd en1 (x1)⋯end (xd)w(x1)⋯w(xd)dx1…dxd =eλ(T) ⟨eγ1·,en1 ⟩w⋯⟨eγd·,end ⟩w

In the arithmetic case, the spot takes the form

g(x,λ)=λ+^γ⊤x,

with againx,γ∈^Rdandλ∈R. SinceR∋y↦^(_γiy)2w(y)is integrable andg(·,λ)∈_Lw,d2. We find the coefficients in this arithmetic case as follows:

_gnd (λ(T))=_∫^Rd (λ(T)+^γ⊤x)_end (x)^w⊗d(x)dx=λ(T)_{∫^Rd end} (x)^w⊗d(x)dx+∑i=1d_{γi ∫^Rd en1} (_x1)w(_x1)d_x1⋯_{∫^Rd eni−1} (_xi−1)w(_xi−1)d_xi−1×_{∫R xi eni} (_xi)w(_xi)d_xi×_{∫^Rd eni+1} (_xi+1)w(_xi+1)d_xi+1⋯_{∫R end} (_xd)w(_xd)d_xd=λ(T)_{⟨e0,en1 ⟩w}⋯_{⟨e0,end ⟩w}+∑i=1d_{γi ⟨e0,en1 ⟩w}⋯_{⟨e0,eni−1 ⟩w ⟨e1,eni ⟩w ⟨e0,eni+1 ⟩w}⋯_{⟨e0,end ⟩w}.

Now, for_nj≥1it holds that_{∫R enj} (y)w(y)dy=_{⟨e0,enj ⟩w}=0. Thus, the only non-zero terms in the sum above are those where_nd=(0,…,0,1,0,…,0), where 1 is appearing in coordinate i, in which case_gnd (λ(T))=_γi. All other_ndwill give_gnd (λ(T))=0, except_nd=(0,…,0)where we get_gnd (λ(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 thatg(x;λ)=ξ(λ+^γ⊤x)forξ(z)=max(z−c,0). As the max-function grows at most linearly, we have0≤g(x;λ)≤|λ|+|^γ⊤x|, and it follows thatg(·;λ)∈_Lw,d2. We may represent the coefficients as

_gnd (λ(T))=E[max(λ(T)−c+^γ⊤Z,0)_end (Z)]

forZ∼N(0,I)with I being thed×didentity matrix. By iterated conditional expectation, conditioning on_Zi,i=1,…,d−1, we can define for R being a standard normal random variable

^gd(_z1,…,_zd−1):=E[max(λ(T)−c+_{γ1 z1}+⋯_{γd−1 zd−1}+_γdR,0)_end (R)]

and iteratively backwardsi=d−1,d−2,…,1,

^gi(_z1,…,_zi−1):=E[^gi+1(_z1,…,_zi−1,R)_eni (R)]

yielding_gnd (λ(T))=^g1.

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)|Ft 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μ∈^Rdand covariance matrixV∈^Rd×d. Here, we collapse the notation to make our discussion more transparent. Hence, we find that the likelihood ratio is

lnℓ(x,y;t,T)∼−12^{(y−μ)⊤ V−1}(y−μ)+12^y⊤y=−12^y⊤(^V−1−I)y+(^{μ⊤ V−1})y−12^{μ⊤ V−1}μ.

Here, I is thed×didentity matrix. It is evident that the functiony↦ℓ(x,y;t,T)∈_Lw,d2if and only if−(^V−1−I)−12I<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 dominating2I, 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 somed×dmatrix such thatCV^C⊤<2I. If we have available such a matrix, we can define a new stochastic processY(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, thenX(t)=^C−1Y(t), and we have

g(x;λ)=g(^C−1y;λ).

In Theorem 2, we assume thatg(^C−1·;λ(T))∈_Lw,d2. Furthermore, for the polynomial process Y we find that the likelihood ratio function is (as a matrix transformation of the Gaussian variable_X(T)|Ft ),

lnℓ(x,y;t,T)∼−12^{(y−Cμ)⊤ (CVC⊤)−1}(y−Cμ)+12^y⊤y=−12^y⊤(^(CVC⊤)−1−I)y+(^{μ⊤ C⊤ (CVC⊤)−1})y−12^{μ⊤ C⊤ (CVC⊤)−1}Cμ

Hence, we have thaty↦ℓ(x,y;t,T)∈_Lw,d2whenever C is such thatCV^C⊤<2I.

Here is an example of a re-scaling: Let_X(T)|Ft be bivariate Gaussian with covariance matrix

V=_{σ12σ1 σ2}ρ_{σ1 σ2}ρ_σ22

where_σi,_σ2are 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 ind=2. Observe that a diagonalisation of V is given by

V=10_{σ2 σ1}ρ1−^ρ2_σ1200_σ221_{σ2 σ1}ρ01−^ρ2

Let, for some positive constantc<2

C:=cmax(_σ1,_σ2)^{10_{σ2 σ1}ρ1−ρ2−1}

Then, we find

CV^C⊤=cmax^(_σ1,_σ2)2_σ1200_σ22<2I

sincec<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 scaler>0and shapem>0, given as

ξ(y)=1Γ(r)^{mryr−1 e−y/m}.

Suppose we have a target distributionϕwhich behaves as^ys−1,s>0fory∼0and^e−y/kfory∼∞, then the likelihood ratio will be^ys−1/^yr−1close to zero, andexp(−(1/k−1/m)y)fory∼∞. However, integrating the square of the likelihood function againstξ, yields finiteness whenevers>r/2andk<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:=^L2(_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ρ,22withρ(x)=w⊗ξ(_x1,_x2)=w(_x1)ξ(_x2)and the canonically generated polynomials from the respective marginal densities.

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

Consider a European option written on a plan-vanilla forward with payoffζ(F(τ,T))at timeτ≤Tfor a payoff functionζ, with the forward priceF(t,T)given as in Proposition 1; that is, for somed,n∈N, we have

F(t,T)=h^(t,T)⊤ _Hn,d(X(t))

where

h^(t,T)⊤=_pn ^(λ(T))⊤ _Γn,dexp(_Gn,d(T−t)).

In the formulation of Proposition 1,_Hn,d(x)is some basis of the nth order polynomials on^Rd. It is convenient, however, for the purpose of option pricing, to turn to the polynomial ONB_{(vnd (x))nd∈N0d}of_Lρ,d2 , as used in Section 4.3 above. We fix_Hn,d(x)to be the basis_{(vnd (x))|nd|≤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 timet≤τbe

P(t,τ,T)=E[ζ(F(τ,T))|_Ft]

where we assumeζ(F(τ,T))∈^L1(P). The following result is essentially a repetition of Theorem 2 and is therefore formulated as a corollary.

Corollary 1.

Assume for all0≤t≤Tthat^Rd∋x↦ζ(h^(t,T)⊤ _Hn,d(x))∈_Lρ,d2 with h given as in (25). Let F be given in (24) with X being a polynomial process on^Rd for which the likelihood ratio function defined in (23) satisfies^Rd∋y↦ℓ(x,y;t,T)∈_Lρ,d2 . Then, we have that

^PN(t,τ,T)→P(t,τ,T)

(pointwise) whenN→∞, whereP(t,τ,T) is the option price in (26) with representation

P(t,τ,T)=∑n∈_N0dζnd (τ,T)_ℓnd (X(t);t,τ)

while for anyN∈N,

^PN(t,τ,T)=∑n∈_N0d,|_nd|≤N_ζnd (τ,T)_ℓnd (X(t);t,τ).

Here,

_ζnd (τ,T)=_∫^Rd ζ(h^(τ,T)⊤y)_vnd (y)^w⊗d(y)dy

and

_ℓnd (x;t,τ)=_und⊤exp(_G|nd|,d(τ−t))_H|nd|,d(x)

with_und ∈^{RK(|_nd|,d)+1}is such that_vnd (x)=_{und⊤ H|nd|,d}(x). We recallK(n,d)=dim_Poln(^Rd)−1 , and G given in (2).

Proof.

The proof is identical to the argument of Theorem 2, but now usingP(t,τ,T)=f(X(t);t,τ,T)where

f(x;t,τ,T):=E[ζ(h^(τ,T)⊤X(τ))|X(t)=x]

Notice thatτnow plays the role of T in the proof of Theorem 2, and thatζis g. We also observe thatζ(F(τ,T))∈^L1(P)under the assumptions on g and X.  □

Ifζ(z)=max(z−K,0), the payoff function of a call option, we find that0≤ζ(h^(t,T)⊤ _Hn,d(x))≤K+∥h(t,T)∥∥_Hn,d(x)∥, using the notation∥·∥for the Euclidean 2-norm on, e.g.,^Rd. This shows readily thatζ(h^(t,T)⊤·)∈_Lρ,d2as long as_Lρ,d2supports polynomials of degree n. This is the case if we choose the space_Lw,d2. 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_Fi(t,T),i=1,2 which are given by (24) for two different functions_hi(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ζ(_F1(τ,_T1)−_F2(τ,_T2)), withτ≤min(_T1,_T2). Hence, we have potentially two different maturities of the forwards. Notice also that the order n and dimensionality d of_Hn,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. Ifx↦ζ((_h1 ^(t,_T1)⊤−_h2 ^(t,_T2)⊤)_Hn,d(x))∈_Lρ,d2, we can apply Corollary 1 withh(t,_T1,_T2):=_h1(t,_T1)−_h2(t,_T2). A typical example of a spread isζ(_F1−_F2)=max(_F1−_F2,0), which satisfies the regularity condition for at least_Lw,d2.

It is remarked in passing that we can also treat quanto options in a similar manner. A quanto option pays_ζ1(_F1(τ,_T1))_ζ2(_F2(τ,_T2))for two payoff functions_ζ1and_ζ2. These may be of the form of two calls, or a call and put, or two puts. Definingζ(_Hn,d(x);t,T):=_ζ1(_h1 ^(t,_T1)⊤ _Hn,d(x))_ζ2(_h2 ^(t,_T2)⊤ _Hn,d(x))and assuming^Rd∋x↦ζ(_Hn,d(x);t,T)∈_Lρ,d2puts 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

F(t,_T1,_T2)=∑T=_T1T2F(t,T)=∑T=_T1T2h^(t,T)⊤_Hn,d(X(t))

we can simply redefine the meaning ofh(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_ζnd (τ,T)and_ℓnd (x;t,τ). The latter family of coefficients,_ℓnd (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_ζnd (τ,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_ζnd (τ,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_ℓnd (x;t,τ)for rather high-dimensional matrices_G|nd|,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.

5.2. 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, ford∈N, recall the Hilbert space_Lρ,d2introduced in the previous section. Consider the “doubled” space_Lρ2,2d2which is the Hilbert space of real-valued functions on^R2dfor which

|g|ρ2,2d2:=∫R2d g2(x,y)ρ2(x,y)dxdy

with inner product

_{⟨g,h⟩ρ2,2d}:=_∫^R2d g(x,y)h(x,y)_ρ2(x,y)dxdy

where_ρ2(x,y):=^ρ⊗2(x,y):=ρ(x)ρ(y). An ONB for_Lρ2,2d2is given by_{(vnd ⊗vkd )(nd,kd)∈N02d}where we recall_{(vnd )nd∈N0d}to be an ONB of_Lρ,d2. We also recall the notationg(·;λ(T)) from (21) for the forward price

F(t,T)=E[g(X(T);λ(T))|_Ft]

in (22), where X is a polynomial process in^Rdwhich in addition is assumed to be Markovian. We further recall that we denote byϕ(x,dy;t,T)the probability distribution function ofX(T)givenX(t)=x∈^Rdfort≥T. Under suitable conditions, Theorem 2 provides an expression and approximation of F. Consider an option written on the forward with exercise timeτ≤Tand payoff functionζ(F(τ,T))for some functionζ:R→Rsuch thatζ(F(τ,T))∈^L2(P). The arbitrage-free option price at timet≤τ(for risk-free interest rate set to zero) is given byP(t,τ,T) as defined in (26).

Theorem 3.

Assumeg(·;λ(T))∈_Lρ,d2and suppose that the likelihood functionℓ(x,y;t,T) defined in (23) satisfies^Rd×^Rd∋(x,y)↦ℓ(x,y;t,T)∈_Lρ2,2d2for any0≤t≤T<∞. Ifζ:R→R is Lipschitz continuous and of linear growth, then

^PN,K(t,τ,T)→P(t,τ,T),

where

P(t,τ,T)=∑_kd∈_{N0d(ζ∘f)kd} (τ,T)_ℓkd (X(t);t,τ),

with

_(ζ∘f)kd (τ,T):=_{⟨ζ(f(·;τ,T)),vkd ⟩ρ,d}.

and, moreover,

^PN,K(t,τ,T)=∑_kd∈_N0d,|_kd|≤K_(ζ∘^fN)kd (τ,T)_ℓkd (X(t);t,τ)

Here,_ℓkd (x;t,τ)is defined in Theorem 2,

f(x;τ,T)=∑_nd∈_{N0dgnd ℓnd} (x;τ,T)

with

_gnd =_{⟨g(·;λ(T)),vnd ⟩ρ,d}

and

^fN(x;τ,T)=∑_nd∈_N0d,|_nd|≤N_{gnd ℓnd} (x;τ,T).

Proof.

By assumptionℓ(·,·;t,T)∈_Lρ2,2d2and thereforey→ℓ(x,y;t,T)∈_Lρ,d2,a.e.,x∈^Rd. Hence, sinceg(·;λ(T))∈_Lw,d2we find a polynomial expression of F as given in Theorem 2 along with an approximation^FN. From the proof of Theorem 2, we also recall the function

f(x;τ,T)=E[g(X(T);λ(T))|_Fτ]

along with its representation and series expansions found in the proof of that result.

Let us show that the mapx↦f(x;τ,T)belongs to_Lρ,d2: indeed, by definition off(x;τ,T), we find from the Cauchy–Schwarz inequality,

|f(·;τ,T)|ρ,d2=∫Rd |⟨g(·;λ(T)),ℓ(x,·;τ,T)⟩ρ,d|2ρ(x)dx≤|g|ρ,d2 ∫Rd ∫Rd ℓ2(x,y;τ,T)ρ(y)ρ(x)dydx=|g|ρ,d2 |ℓ(·,·;τ,T)|ρ2,2d2

which is finite by the assumption on the likelihood ratio function. It follows thatf(·;τ,T)∈_Lρ,d2. Moreover, from Theorem 2, we find that^fN(·;τ,T)→f(·;τ,T)in_Lρ,d2whenN→∞.

By assumption,ζhas linear growth. Hence, for some constantk>0, it follows fromf(·;τ,T)∈_Lρ,d2that

_∫^Rd ζ^(f(x;τ,T))2ρ(x)dx≤k_∫^Rd (1+^f2(x;τ,T))ρ(x)dx<∞.

In other words,x↦ζ(f(x;τ,T))∈_Lρ,d2. Hence,

c(x;t,τ,T):=E[ζ(f(X(τ);τ,T))|X(t)=x]=_∫^Rd ζ(f(y;τ,T))ϕ(x,dy;t,τ)=_{⟨ζ(f(·;τ,T)),ℓ(x,·;t,τ)⟩ρ,d}=∑_kd∈_{N0d(ζ∘f)kd} (τ,T)_ℓkd (x;t,τ)

where

_(ζ∘f)kd (τ,T):=⟨ζ_{(f(·;τ,T),vkd ⟩ρ,d}.

We findP(t,τ,T)=c(X(t);t,τ,T).

We truncate this sum at multi-indices of order up to K and consider the approximation

^cN,K(x;t,τ,T):=_{⟨ζ(^fN(·;τ,T)),^ℓK(x,·;t,τ)⟩ρ,d}.

By the triangle inequality, after subtracting and adding_{⟨ζ(f(·;τ,T)),^ℓK(x,·;t,τ)⟩ρ,d}, we find from the Cauchy–Schwarz inequality

|c(x;t,τ,T)−^cN,K(x;t,τ,T)|≤|_{⟨ζ(f(·;τ,T)),ℓ(x,·;t,τ)−^ℓK(x,·;t,τ)⟩ρ,d}|+|_{⟨ζ(f(·;τ,T))−ζ(^fN(·;τ,T)),^ℓK(x,·;t,τ)⟩ρ,d}|≤_{|ζ(f(·;τ,T))|ρ,d |ℓ(x,·;t,τ)−^ℓK(x,·;t,τ)|ρ,d}+|ζ(f(·;τ,T))−ζ(^fN(·;τ,T))_{|ρ,d |^ℓK(x,·;t,τ)|ρ,d}≤_{|ζ(f(·;τ,T))|ρ,d |ℓ(x,·;t,τ)−^ℓK(x,·;t,τ)|ρ,d}+k|f(·;τ,T)−^fN _{(·;τ,T)|ρ,d |^ℓK(x,·;t,τ)|ρ,d}

In the last inequality, we appealed to the Lipschitz-continuity ofζ(denoting the Lipschitz constantk>0). Recall by the definitions ofℓ(x,y;t,τ)and^ℓK(x,y;t,τ)and the analysis in the proof of Theorem 2 that^ℓK(x,·;t,τ)→ℓ(x,·;t,τ)in_Lρ,d2whenK→∞. Moreover,|^ℓK _{(x,·;t,τ)|ρ,d}≤_{|ℓ(x,·;t,τ)|ρ,d}(indeed, the norm of^ℓK(x,·;t,τ)converges to that ofℓ(x,·;t,τ)!) To conclude, it remains to recall that^fN(·;τ,T)converges tof(·;τ,T)in_Lρ,d2, forN→∞, as shown above.  □

To apply Theorem 3 in practice, we need to compute the coefficients_(ζ∘^fN)kd for all_kd∈_N0dsuch that|_kd|≤K. We have that^fNis again a truncated sum, but the coefficients_gnd of this is available from the computation of forward prices (or approximations thereof). This representation can be used to calculate_(ζ∘^fN)kd , 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

P(t,τ,_T1,_T2)=E[ζ(∑T=_T1T2F(τ,T))|_Ft]

To apply the arguments of Theorem 3, we must assume that_∑T=T1T2 g(·;λ(T))∈_Lρ,d2. Ifg(·;λ(T))∈_Lρ,d2for any T, then this condition holds as_Lρ,d2is a vector space. Tracing the steps in the proof of Theorem 3 results in

P(t,τ,_T1,_T2)=∑_kd∈_{N0d(ζ∘f˜)kd} (τ,_T1,_T2)_ℓkd (X(t);t,τ),

with

_(ζ∘f˜)kd (τ,T):=_{ζ(∑T=T1T2f(·;τ,T)),vkd ρ,d}.

On the other hand, as long asζis of bounded linear growth and Lipschitz continuous,^PN,K(t,τ,_T1,_T2)→P(t,τ,_T1,_T2)with

^PN,K(t,τ,_T1,_T2)=∑_kd∈_N0d,|_kd|≤K_{(ζ∘^f˜N)kd} (τ,T)_ℓkd (X(t);t,τ)

and where we find that

^f˜N(x;τ,_T1,_T2)=∑_nd∈_N0d,|_nd|≤N_gnd ∑T=_T1T2ℓnd (x;τ,T).

In fact,^f˜N(x;τ,_T1,_T2)=_∑T=T1T2 ^fN(x;τ,T), thus, to approximate the option price on a temperature HDD or CDD futures, we simply add up a finite sum of terms^fN(x;τ,T)over T when computing the coefficients_{(ζ∘^f˜N)kd} rather than using only one^fN(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_ℓkd (x;τ,T)with respect to x. This lowers the polynomial order of_ℓkd (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 aboutN=10for 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_ℓnd (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_ℓnd (x;τ,T). The convergence speed for these expressions should be further analysed. These are challenging computational problems which we leave for future studies.

6. 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 square-integrability 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.

Funding

The author acknowledges financial support from the thematic research group SPATUS funded by UiO:Energy, University of Oslo.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

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.

References

1. Björk, T. Arbitrage Theory in Continuous Time Finance, 3rd ed.; Oxford University Press: Oxford, UK, 2009.

2. Geman, H. Commodities and Commodity Derivatives; John Wiley & Sons: Chichester, UK, 2005.

3. Eydeland, A.; Wolyniec, K. Energy and Power Risk Management; Wiley-Finance; John Wiley & Sons: Hoboken, NJ, USA, 2003.

4. Benth, F.E.; Benth, J.X.; Koekebakker, S. Stochastic Modelling of Electricity and Related Markets; World Scientific: Singapore, 2008.

5. Cuchiero, C. Affine and Polynomial Processes. Ph.D. Thesis, ETH, Zürich, Switzerland, 2011.

6. Filipović, D.; Larsson, M. Polynomial diffusions and applications in finance. Financ. Stoch. 2016, 4, 931-972.

7. Cuchiero, C.; Keller-Ressel, M.; Teichmann, J. Polynomial processes and their applications to mathematical finance. Financ. Stoch. 2012, 16, 711-740.

8. Filipović, D.; Larsson, M. Polynomial jump-diffusion models. Stoch. Syst. 2020, 10, 71-97.

9. Ackerer, D.; Filipović, D.; Pulido, S. The Jacobi stochastic volatility model. Financ. Stoch. 2018, 22, 667-700.

10. Ware, T. Polynomial processes for power prices. Appl. Math. Financ. 2019, 26, 453-474.

11. Kleisinger-Yu, X.; Komaric, V.; Larsson, M.; Regez, M. A multi-factor polynomial framework for long-term electricity forwards with delivery period. SIAM J. Finan. Math. 2020, 11, 928-957.

12. Karatzas, I.; Shreve, S.E. Brownian Motion and Stochastic Calculus, 2nd ed.; Springer: New York, NY, USA, 1991.

13. Schwartz, E.S. The stochastic behaviour of commodity prices: Implications for valuation and hedging. J. Financ. 1997, 52, 923-973.

14. Gibson, R.; Schwartz, E.S. Stochastic convenience yield and the pricing of oil contingent claims. J. Financ. 1990, 45, 959-976.

15. ESchwartz, S.; Smith, J.E. Short-term variations and long-term dynamics in commodity prices. Manag. Sci. 2000, 46, 893-911.

16. Lucia, J.J.; Schwartz, E.S. Electricity prices and power derivatives: Evidence from the Nordic Power Exchange. Rev. Deriv. Res. 2001, 5, 5-50.

17. Prokopczuk, M. Pricing and hedging in the freight futures market. J. Futures Mark. 2011, 31, 440-464.

18. Nomikos, N.K.; Soldatos, O. Using affine jump diffusion models for modelling and pricing electricity derivatives. Appl. Math. Financ. 2008, 15, 41-71.

19. Pilipović, D. Energy Risk-Valuing and Managing Energy Derivatives; McGraw-Hill: New York, NY, USA, 1998.

20. Cartea, A.; Figueroa, M. Pricing in electricity markets: A mean reverting jump diffusion model with seasonality. Appl. Math. Financ. 2005, 12, 313-335.

21. Mirantes, A.G.; Población, J.; Serna, G. The stochastic seasonal behaviour of natural gas prices. Europ. Financ. Manag. 2012, 18, 410-443.

22. Benth, F.E.; Benth, J.X. Modeling and Pricing in Financial Markets for Weather Derivatives; World Scientific: Singapore, 2013.

23. Härdle, W.; Lopez-Cabrera, B. The implied market price of weather risk. Appl. Math. Financ. 2012, 18, 59-95.

24. Swishchuk, A.; Cui, K. Weather derivatives with applications to Canadian data. J. Math. Financ. 2013, 3, 81-95.

25. Paschke, R.; Prokopczuk, M. Commodity derivatives valuation with autoregressive and moving average components in the price dynamics. J. Bank. Financ. 2010, 34, 2742-2752.

26. Benth, F.E.; Pircalabu, A. A non-Gaussian Ornstein-Uhlenbeck model for pricing wind power futures. Appl. Math. Financ. 2018, 25, 36-65.

27. Bensoussan, A.; Brouste, A. Cox-Ingersoll-Ross model for wind speed modeling and forecasting. Wind Energy 2016, 19, 1355-1365.

28. Benth, F.E.; Rohde, V. On non-negative modeling with CARMA processes. J. Math. Anal. Appl. 2019, 476, 196-214.

29. Kyriakou, I.; Nomikos, N.K.; Papapostolou, N.C.; Pouliasis, P.K. Affine-structure models and the pricing of energy commodity derivatives. Eur. Financ. Manag. 2016, 22, 853-881.

30. Jewson, S.; Brix, A. Weather Derivative Valuation; Cambridge University Press: Cambridge, UK, 2005.

31. Hinderks, W.J.; Wagner, A. Pricing German energiewende products: Intraday cap/floor futures. Energy Econ. 2019, 81, 287-296.

32. Caporin, M.; Preś, J.; Torro, H. Model based Monte Carlo pricing of energy and temperature quanto options. Energy Econ. 2012, 34, 1700-1712.

33. Weron, R. Market price of risk implied by Asian-style electricity options and futures. Energy Econom. 2008, 30, 1098-1115.

34. Fusai, G.; Marena, M.; Roncoroni, A. Analytical pricing of discretely monitored Asian-style options: Theory and application to commodity markets. J. Bank. Financ. 2008, 32, 2033-2045.

35. Kyriakou, I.; Pouliasis, P.K.; Papapostolou, N.C. Jumps and stochastic volatility in crude oil prices and advances in average option pricing. Quant. Financ. 2016, 16, 1859-1873.

36. Shiraya, K.; Takahashi, A. Pricing average options on commodities. J. Futures Mark. 2011, 31, 407-439.

37. Filipović, D.; Willems, S. A Term structure model for dividends and interest rates. Math. Financ. 2020, 40, 1461-1496.

38. Folland, G.B. Analysis-Modern Techniques and Their Applications; John Wiley & Sons: Hoboken, NJ, USA, 1984.

39. Asmussen, S.; Goffard, P.-O.; Laub, P.J. Orthonormal Polynomial Expansions and Lognormal Sum Densities. In Risk and Stochastics: Ragnar Norberg; Barrieu, P., Ed.; World Scientific: Singapore, 2019; Chapter 6; pp. 127-150.

40. Eggers, H.C. From Gram-Chalier series to orthogonal polynomials. Acta Phys. Pol. 2009, 40, 1209-1215.

41. Rahman, S. Wiener-Hermite polynomial expansion for multivariate Gaussian probability measures. J. Math. Anal. Appl. 2017, 454, 303-334.

42. Szegö, G. Orthogonal Polynomials; American Mathematical Society Colloquium Publications: New York, NY, USA, 1939; Volume XXIII.

43. Ackerer, D.; Filipović, D. Option pricing with orthogonal polynomial expansions. Math. Financ. 2020, 30, 47-84.

44. Benth, F.E.; Lavagnini, S. Correlators of polynomial processes. arXiv 2020, arXiv:1906:11320.