677
Views
0
CrossRef citations to date
0
Altmetric
Research Article

Subordinated affine structure models for commodity future prices

ORCID Icon, & | (Reviewing editor)
Article: 1512360 | Received 02 Apr 2018, Accepted 11 Aug 2018, Published online: 07 Sep 2018

Abstract

To date the existence of jumps in different sectors of the financial market is certain and the commodity market is no exception. While there are various models in literature on how to capture these jumps, we restrict ourselves to using subordinated Brownian motion by an α-stable process, α ∈ (0,1), as the source of randomness in the spot price model to determine commodity future prices, a concept which is not new either. However, the key feature in our pricing approach is the new simple technique derived from our novel theory for subordinated affine structure models. Different from existing filtering methods for models with latent variables, we show that the commodity future price under a one factor model with a subordinated random source driver, can be expressed in terms of the subordinator which can then be reduced to the latent regression models commonly used in population dynamics with their parameters easily estimated using the expectation maximisation method. In our case, the underlying joint probability distribution is a combination of the Gaussian and stable densities.

AMS Subject Classification:

PUBLIC INTEREST STATEMENT

This paper is entitled Subordinated Affine Structure Models for Commodity Futures Prices. The aim is to provide a mathematical background on a special family of processes that could capture unusual patterns such as extreme events in the commodities market. Such events could range from wild fires on farms to Tsunamis and as a result the commodities spot prices would be affected significantly. This realisation provides the motivation of the study in this paper. As a risk mitigation process, it is undeniable that investors in the commodities market would be interested in how this kind of risk could be captured in the future pricing model and more so, how the model would be calibrated to historical data. The current paper seeks to answer these questions and establishing a novel approach on pricing commodities futures.

1. Literature

A large volume of literature on commodities market has been published since the invention of the continuous benchmark model of Black and Scholes (Citation1973) for pricing options and corporate liabilities. Among the many models developed, the widely used and referenced study on commodities is the work by Schwartz (Citation1997). From the latter, numerous models have been developed as a result of the growing commodities market in terms of volumes traded and complexities of their contracts over the years. We give an account of the various literatures relevant to this research in Table .

Table 1. Various commodity jump models

The current research builds on findings from Kyriakou, Nomikos, Pouliasis, and Papapostolou (Citation2016) by extending the results to subordinated Brownian motion.

A selection of key commodity jump models that have developed overtime.

2. Introduction

Commodities exhibit distinctive features that a good model ought to capture. For instance to estimate the commodities market as closely as possible, one has to factor in jumps in the underlying spot price. However, models designed with a jump component are non-trivial. In this research, we derive a relatively easy estimation method for commodities prices using subordination as a proxy for introducing jumps. Other features include mean-reversion, contango, backwardation and seasonality. Commodities also experience extreme volatility and price spikes resulting in heavy-tailed distribution of the returns. Commodity markets are unique compared with other markets such as equity, bond, currency or interest rate markets in the sense that most commodities are real physical assets that are produced, transported, stored and consumed. They are not assets valued on long-lived companies like in equity markets.

As indicated in Fama and French (Citation1988), commodity pricing can be approached from two perspectives, the theory of storage which explains why high supplies and inventories running at minimum would result into contango, low futures and spot price volatilities, and in turn futures premiums being equivalent to full storage costs. On the other hand, why low supplies and enhanced production inventory levels yield to backwardation, a rise in volatilities of the spot and the nearby future prices. Another feature explained by this theory is the periodically continuously compounded convenience yield (usually denoted by δ) on inventory which is the benefit of holding a physical commodity as opposed to having a future contract of its delivery at some future time and second, the cost of storage. The future price motivated by the theory of storage is given by

(1) Ft,T=Ste(rcδ)(Tt),(1)

where c accounts for the storage costs, r is the periodically continuously compounded interest rate, St is the current spot price and T is the maturity date of the future contract.

The second perspective is the theory of expected risk premium discussed in Keynes (Citation1930) and Hicks (Citation1939). It asserts that the future prices are given by the discounted (by the risk premium) expected future spot price:

(2) Ft,T=Et[ST]erγ[Tt],(2)

where γ is the risk premium and Et[]=E[|Ft], Ft is the filtration up to time t.

A number of models based on the latter approach have been developed over the years to mimic the market as closely as possible for various commodities. This includes Schwartz’s common continuous stochastic factor models Schwartz (Citation1997), Schwartz and Smith (Citation2000) and the jump models of Kyriakou et al. (Citation2016).

The motivation and contribution of this paper are based on the existing erratic features in electricity and energy markets, where jumps are evident, resulting in skewed distributions of the spot prices. We consider a subordinated Brownian motion by an α-stable process, α(0,1), as the source of randomness in the underlying asset to model commodity future prices. The stunning feature in our pricing approach is the new simple technique derived from our novel approach for subordinated affine structure models.

We show that the affine property is attainable and applicable to generalised commodity spot models, and as illustration, we consider a stochastic differential equation with subordinated Brownian motion as the source of randomness to derive the commodity future prices. It is argued in some existing literature that the likelihood function exists in integrated form for models with singular noise meanwhile for cases of partially observed processes, a filtering technique is required, see for instance Date and Ponomareva (Citation2010) and Yang, Lin, Luk, and Nahar (Citation2014). However, the work presented in this paper provides a new approach of pricing commodity futures for models with latent variables using the maximum expectation maximisation. We show that the commodity future price under a one-factor model with a subordinated random source driver can be expressed in terms of the subordinator, which can then be reduced to the latent regression models commonly used in population dynamics with their parameters estimated easily using the expectation maximisation method. In our case, the underlying joint probability distribution is a combination of Gaussian and stable densities.

The rest of the paper is organised as follows. The following Section 3 introduces some features of stable processes essential to this work. In Section 4, we review the concept of affine models and extend the idea of obtaining Laplace transforms of random processes to subordinated processes. In Section 5, we derive our pricing formulas for commodity futures using the results derived in Section 4. In Section 6, we discuss the numerical implementation of our one-factor commodity future models. Section 7 concludes.

3. Stable processes

The discussion in this section is mainly based on Kateregga et al. (Citation2017). A stable or α-stable process, α(0,2], belongs to the general class of Lévy distributions. It has a limiting distribution with a definitive exponent parameter α that determines its shape. The following two definitions follow from Rachev (Citation2003).

Definition 3.1 Let X1,X2,,Xn be independent and identically distributed random variables and suppose a random variable S defined by

(3) 1ani=1nXibnS,(3)

where “→” represents weak convergence in distribution, an is a positive constant and bn is real. Then, S is a stable process. The constants an and bn need not to be finite.

Definition 3.1 allows the modelling of a number of natural phenomena beyond normality using stable distributions. The fact that an and bn do not necessarily need to be finite provides the generalised central limit theorem.

Definition 3.2 (Generalised Central Limit Theorem) Suppose X1,X2, denotes a sequence of identically distributed independent random variables from some arbitrary distribution and let sequences anR and bnR+. Then, we define a sequence

(4) Zn:=1bni=1nXian(4)

of sums Zn such that their distribution functions weakly converge to some limiting distribution. That is

(5) P(Zn < x)H(x),n,(5)
where H(x) is some limiting distribution.

The traditional central limit theorem assumes finite mean a:=E[Xi] and finite variance σ2:=Var[Xi] and defines the sequence of sums

(6) Zn:=1σni=1nXina,(6)
such that the distribution functions of Zn weakly converge to hsG(x). That is
(7) P(x1 < Zn < x2)x1x2hsG(x)dx,n(7)
where hsG(x) denotes the standard Gaussian density
(8) hsG(x)=12πexp(x2/2).(8)
Suppose the identically distributed independent random variables Xi equal to a positive constant c almost surely and the sequences an and bn in (4) are defined by an=(n1)c and bn=1, then Zn is also equal to c for all n > 0 almost surely. In this case, the random variables Xi are mutually independent and as a consequence, the limiting distribution for the sums Zn belongs to the stable family of distributions by definition. Therefore, stable distributions behave similarly to the central limit theorem for distributions with a finite second-order moment (the Gaussian), Crosby (Citation2008). This is one reason why they are regarded as stable. They are also preferred compared with all other laws such as the Normal Inverse Gaussian (NIG), Variance Gamma (VG) and other distributions from the generalised hyperbolic family because of their heavier tails.

Definition 3.3 Samoradnitsky and Taqqu (Citation1994) An α-stable distribution is a four-parameter family of distributions denoted by S(α,β,ν,μ) where

(1) α(0,2] is the characteristic exponent responsible for the tail of the distribution.

(2) β[1,1] is responsible for skewness.

(3) ν > 0 is the scale parameter (sometimes referred to as variance when α=2).

(4) μR is the location (sometimes referred to as mean).

Densities of α-stable distributions do not have closed-form representations except the Gaussian (α=2), Cauchy (α=1) and Inverse Gaussian or Pearson (α=0.5) distributions.

The analysis of stable processes is usually through their characteristic functions and Laplace or Fourier transformation. Unlike their densities, their characteristic functions always exist. Literature on their integral representations and density functions is provided in Zolotarev (Citation1964), Citation1980, Zolotarev (Citation1986)). The distribution functions for the different α values have been tabulated in Dumouchel (Citation1971), Fama and Roll (Citation1968) and Holt and Crow (Citation1973).

Definition 3.4 (Gajda and Wyłoman`Ska (Citation2012)) Let St and Lt denote an α-stable process and its respective inverse. Then for t[0,T], we define the process Lt asFootnote1

(9) Ls:=inft:St>sifs[0,St)Tifs=ST.(9)

St and Lt are non-decreasing and cádlág with their graphical representations given in Figure .

Figure 1. The top graph shows a plot of a stable process St and the bottom graph shows its inverse process Lt simulated using exponent parameter value, α=0.8, plotted against time on the horizontal.

Figure 1. The top graph shows a plot of a stable process St and the bottom graph shows its inverse process Lt simulated using exponent parameter value, α=0.8, plotted against time on the horizontal.

3.1. Density and characteristic functions

Let (Xt,t0) denote a Lévy process. The characterisation of Xt is usually given by the Lévy–Khintchine formula.

Definition 3.5 (Lévy–Khintchine) Applebaum (Citation2004) Consider a Lévy process X=(Xt)t0. There exists bR and σ0 such that the characteristic function of X is

(10) Φ(t):=E[eitX]=expitb12σ2t2+R0(eitx1itxIx<1)Π(dx),(10)

where I is the indicator function and Π is a σ-finite measure satisfying the constraint

(11) Rd0min(1,|x|2)Π(dx) < ;alternativelyRd0|x|21+|x|2Π(dx) < .(11)

Definition 3.6 (The Lévy-Itô Decomposition) Applebaum (Citation2004) If Xt is a Lévy process, there exists bR, a Brownian motion Bσ(t) with variance σR+ and an independent Poisson random measure N on R+×(R0) such that, for each t0,

(12) Xt=bt+Bσ(t)+|x|<1xN˜(t,dx)+|x|1xN(t,dx),(12)

where

(13) b=EX1x1xN(1,dx).(13)

To preserve the martingale property, the compound Poisson random measure is compensated as N˜=Ntλ where λ is a Lévy measure satisfying (11).

For process St, we require σ=0 in (10) or B=0 in (12) and the Lévy measure in (11) given by

(14) Π(x)=C|x|1+αdx;C > 0,(14)

The characteristic function Φ of St is obtained using the domain of attraction of stable random variables (See Grigelionis, Kubilius, Paulauskas, Statulevicius, & Pragarauskas, Citation1999) and the Lévy–Khinchine representation formula (See Definition 3.5 or Applebaum (Citation2004) for a detailed explanation), i.e.

(15) Φ(θ)=E[exp(iθS)]=exp(να|θ|α[1iβ sign(θ)tan(πα2)]+iμθ);forα1.exp(νθ[1+iβ sign(θ)2πlogθ]+iμθ);forα=1.(15)

Using Fourier transformation, the density function of St is given by

(16) hS(t,u)=1π0eiuθΦ(θ)dθ.(16)

From Definition 3.4, it is easy to see the equivalence relation

(17) Su < tLtu.(17)

It follows that F(t,u):=P(Su < t)=P(Lt  u)=0hL(τ;t)τ where hL(u,t) denotes the probability density function of Lt. Consequently

(18) hL(u,t)=F(t,u)u=uthS(τ,u)dτ.(18)

According to Meerschaert and Straka (Citation2013), the density h(t,u) can also be given by

(19) hS(t,u)=u1/αh(tu1/α),(19)

where h(τ) is the density of a standard stable process with a Laplace transform h˜(τ)=exp(τα). This follows from the fact that Su has the same distribution as u1/αS1. As a result, the density of the inverse stable process Lt can be given in terms of the standard stable process by

(20) hL(u,t)=tαu11/αh(tu1/α).(20)

Figure shows density graphs of St for different exponent parameter values.

Figure 2. The graphs represent densities of an α-stable process for different values of the exponent parameter, α (0,2]. Observe the variation in the tail sizes and the skewness as the exponent parameter is varied.

Figure 2. The graphs represent densities of an α-stable process for different values of the exponent parameter, α ∈(0,2]. Observe the variation in the tail sizes and the skewness as the exponent parameter is varied.

3.2. Laplace transforms

Definition 3.7 Let Xu be a subordinator. The Laplace transform of Xu is defined by

(21) E[eτXu]=euϕ(τ),(21)

where ϕ is the Laplace exponent of Xu known as the Bernstein function represented by

(22) ϕ(τ)=a+bτ+(0,)(1eτx)Π(dx).(22)

where a,b > 0 and Π is the Lévy measure on (0,) such that x1+xΠ(dx)<.

The Laplace transform of the stable process Su is given by (see Meerschaert & Straka, Citation2013)

(23) h˜S(τ,u)=0etτhS(t,u)dt=exp(uCΓ(1α)τα)=exp(u((τ+β)α+βα)).(23)

where 0β1. For C=Γ(1α) (alternatively β=0), the Laplace transform simplifies to that of a standard stable process:

(24) h˜S(τ,u)=E[eτSu]=exp(uτα);0 < α < 1.(24)

The Laplace transform h˜L(u;τ) of the inverse stable process Lt is obtained from (18):

(25) h˜L(u,τ)=u(τ1exp(u((τ+β)α+βα)))=τ1((τ+β)α+βα)exp(u((τ+β)α+βα)).(25)

where the Laplace transform of 0tf(y)y is τ1f˜(τ) and hL(u,τ):=0 for l < 0 or τ < 0. Since (25) does not have the general form for a Laplace transform of a Lévy process, then Lt is not a Lévy process.

3.3. Moment-generating function

There is a relationship between a moment-generating function of a random variable and its Laplace transform.

Lemma 3.8 Let Mu(τ) and h˜(τ,u) denote the respective moment-generating function and Laplace transform of a random variable then

(26) Mu(τ)=h˜(τ,u)+h˜(τ,u),(26)

where

(27) h˜(a,b)=0etah(t,b)dt.(27)

Proof 3.9 This relationship is verified in Miller (Citation1951).

As a consequence of Lemma 3.8 and the explicit Laplace transform given by (23), we can deduce the first and second moments of St. That is

Mu(τ)=exp(u((τ+β)α+βα))+exp(u((τ+β)α+βα)).
(28) E[St]=Mu(0)=αuβα1[e2uβα+e2uβα].(28)
(29) Var[St]=Mu ′′(0)(Mu(0))2.(29)
(30) =α(α1)uβα2[e2uβαe2uβα]+α2uαβ2α[β2e2uβα+euβα](30)

3.4. Subordination

Let BSt denote subordinated Brownian motion where St is an α-stable process introduced before with α(0,1) and B represents standard Brownian motion with mean zero. We require Bt and St to be independent. To ensure a complete working environment, we introduce a probability space.

Let (Ω,F,P)=(ΩB×ΩS,FB×FS,μB×μS) denote a complete joint probability space endowed with a filtration (Ft)t  0 such that Ft=FtBFtS where FtB and FtS are filtrations generated by Bt and St, respectively.

Definition 3.10 Let B=(Bt,Px)=2Bt where (Bt)t  0 is standard Brownian motion in R. The transition density p(x,y;t) of B is given by

(31) p(x;t)=12πtexpx2/4t,t>0,x,yR.(31)

The semi-group (Pt:t0) of B is given by

(32) Ptf(x)=Ex[f(Bt)]=Rp(x,y;t)f(y)dy,(32)

where f is a non-negative Borel function on R satisfying the following generator equation:

(33) Gf(x):=limt0Ex[f(X)]f(x)t=limt0Ptf(x)f(x)t.(33)

Definition 3.11 Suppose Yt:=BSt,t0 is a subordinated Brownian motion. Its’ semi-group (Qt:t0) is defined by

(34) Qtf(y,t)=Ey[f(Yt,t)]=Ey[f(BSt)]=0Puf(y,u)hS(t,u)du.(34)

The semi-group Qt has a transition density q(y,z,t)=q(yz,t) defined by

(35) q(y,t)=0p(y,u)hS(t,u)du.(35)

Lemma 3.12 The mean and variance of BSt are given by

(36) Ey[BSt]=0Ey[Bs]hS(t,u)du=0.(36)
(37) Ey[BSt2]=0Ey[Bs2]hS(t,u)du=0uhst,udu=EySt(37)

Proof. Suppose f in Definition 3.10 and Definition 3.11 is such that f(z,t)=z. Using (32) and (34) and partitioning the time interval [0,T] such that 0 < τ1 < < τn < T, where τi are the jump times of the process BSt, we observe that (36) and (37) hold for every interval [τi,τi+1). Thus, in the limits, their sums converge, respectively, to 0 and Ex[St] on [0,T].

Lemma 3.13 Let X be a Lévy process with characteristic exponent Ψ and S an independent subordinator with Laplace exponent Φ. Then, the subordinated process XSt is a Lévy process with characteristic exponent

(38) m()=Φ(Ψ()).(38)

Proof 3.14 The proof is given in Bochner (Citation2012).

It is known that any Lévy process Xs, t < s  T with drift μ is fully determined by its characteristic function given by (see Fusai & Roncoroni, Citation2007)

(39) E[eiλXs]=eμΔ+Ψ(λ)Δ(39)

where Δ=st, μ is the drift parameter and Ψ(λ) is the characteristic exponent. A typical example is Brownian motion whose characteristic function is given by

(40) E[eiλBs]=eμΔ12σ2λ2Δ,whereΨ(λ)=12σ2λ2.(40)

The characteristic exponent of subordinated Brownian motion BSt is deduced from (23), (38), (40):

(41) m(u)=12σ2u2+βα+βα.(41)

4. Affine models

In this section, we provide an overview on affine processes. We retain some definitions and notations used in Keller-Ressel (Citation2008); Keller-Ressel, Schachermayer, and Teichmann (Citation2011) and Duffie, Filipovi´C, and Schachermayer (Citation2003):

(1) D:=R0m×Rn.

(2) U:=uCd:ReuI0,ReuJ=0,

where I:=1,,m, J:=m+1,,m+n and M:=IJ=1,,d.

(3) Ptf(x)=Ex[f(Xt)] for all xD,t0.

(4) O:=(t,u)R0×U:Psfu(0)0s[0,t].

(5) X will denote a closed state space.

Definition 4.1 (Keller-Ressel, Citation2008) A process is stochastically continuous if for any sequence tnt in R0, then the random variables Xtn converge to Xt in probability with respect to (Px)xD.

Definition 4.2 (Keller-Ressel, Citation2008) An affine process is a stochastically continuous time-homogeneous Markov process (Xt,Px) with a state space D where the characteristic function is an exponentially affine function of the state vector. In other words, there exist functions ψ0:R0×iRdC and ψ1:R0×iRdCd such that

(42) Ex[euXt]=exp(ψ0(t,u)+ψ1(t,u)x),(42)

for all xD and for all (t,u)R0×iR.

Definition 4.1 can be extended to O satisfying the following properties [Prop. 1.3, Keller-Ressel (Citation2008)]:

(1) ψ0 maps O to C where C:=uC:Reu0.

(2) ψ1 maps O to U.

(3) ψ0(0,u)=0 and ψ1(0,u)=u for all uU.

(4) ψ0 and ψ1 admit the `semi-flow property’:

ψ0(t+s,u)=ψ0(t,u)+ψ0(s,ψ1(t,u)),

ψ1(t+s,u)=ψ1(s,ψ1(t,u)), for all t,s0 with (t+s,u)O.

(5) ψ0 and ψ1 are jointly continuous on O.

(6) With the remaining arguments fixed, uIψ0(t,u) and uIψ1(t,u) are analytic functions in uI:ReuI<0;(t,u)O.

(7) Let (t,u),(t,w)O with ReuRew. Then

Reψ0(t,u)ψ1(t,Rew),

Reψ1(t,u)ψ1(t,Rew).

Lemma 4.3 An affine process (Xt)t  0 is regular if the following right derivatives exist for all uU and are continuous at u=0:

(43) F(1)u:ψ0t(t,u) t=0+,F(2)u:ψ1t(t,u)t=0+(43)

The regularity condition can be extended to O for which case the following Riccati equations hold:

(44) ψ0t(t,u)=F(1)(ψ1(t,u)),ψ0(0,u)=0,(44)
(45) ψ1t(t,u)=F(2)(ψ1(t,u)),ψ1(0,u)=u.(45)

Proof 4.4 See Keller-Ressel (Citation2008); Keller-Ressel et al. (Citation2011); Rouah and Heston (Citation2015).

Functions ψ0, ψ1 can be characterised respectively by admissible sets of parameters (a,b,c,π) and (p,q,r,μ) where π and μ are measures. The details are given in Keller-Ressel (Citation2008); Keller-Ressel et al. (Citation2011). We are interested in the affine property of the solution to the SDE:

(46) Xt=b(Xt)dt+σ(Xt)dMt,(46)

where b:XRd is continuous, σ:XRd×d is measurable such that the diffusion matrix σ(x)σT is continuous and Mt is a d-dimensional standard Lévy process such as Brownian motion.

The following theorem is one of the contributions in this paper.

Theorem 4.5 Suppose Xt is a regular affine solution to (46). Then b and σ can be expressed as:

(47) b(x)=K0+K1x1++Kdxd,KiRd(47)
(48) σ(x)σ(x)T=H0+H1x++Hdxd,HiRd×Rd,(48)

where i=0,d. Moreover, the characteristic function of Xt has a log-linear form

(49) E[eiu1Xt(1)++iudXt(d)](49)
(50) =exp(ψ0(t,u1,,ud)+ψ1(t,u1,,ud)x0(1)++ψd(t,u1,,ud)x0(d)).(50)

where uiR. The coefficients ψi are obtained by solving the system of Riccati equations:

(51) F(i)(t,ψ1,,ψd)=ψit=KiTη+12ηTHiη,i=0,,d,(51)

where ηT=(ψ1ψd).

Proof. The proof is a generalisation of the two-dimensional case in Rouah and Heston (Citation2015).  ∎

There is extensive literature on affine processes Xt where M:=Bt or M:=Bt+σ10tξss with ξt, a Poisson jump process. We are interested in the solution to

(52) dYt=b(Yt)dSt+σ(Yt)dBSt.(52)

Another contribution follows in the following theorem. We show that Yt=XSt is affine in the following theorem with d=1.

Theorem 4.6 Let (Ω,F,Px,(Ft)t0) denote a joint probability space for (St)t0, a non-decreasing affine process taking values in D and (Xt)t0,X0=x, an independent Lévy process. Define a process Yt:=XSt,Y0=y with Lévy exponent m(w) and suppose (St)t0 is regular with functional characteristics F(1)(u),F(2)(u). Then (Yt)t0 is regular affine with functional characteristics F(1)(m(w)),FX(2)(m(w)) and FY(2)(m(w))=0, u,wiR with the characteristic function given by

(53) Es[eiuYt]=exp(ψ0(t,m(w))+ψ1(t,m(w))St),(53)

for some functions ψ0 and ψ1.

Proof 4.7 The Markov property of St and the definition of its Laplace transform yields

(54) Es[ewYτFt]=Es[ewXSt]=Es[Es[exp(wXSt)σ(Ss)0st]].=Es[exp(m(w)St)].=exp(ψ0(t,m(w))+ψ1(t,m(w))St).(54)

The last equality follows from the affine property of St (see Definition 4.2).

5. Commodity future pricing

5.1. Introduction

We develop representation formulas for future prices using the concepts introduced before. The source of randomness in the models developed in this section is Brownian motion subordinated by a non-decreasing α-stable process where α(0,1). The aim is to obtain future price formulas for commodity spot price models that incorporate stochastic volatility, jumps, seasonality and mean-reversion effects.

5.2. One-factor commodity spot model

We consider a one-factor commodity spot price model given by

(55) zt=f(t)+exp(Yt),(55)

where Yt satisfies (52) and seasonality is defined according to Kyriakou et al. (Citation2016), as

(56) f(t)=δ0t+δ1sin(δ2[t+2π/264])+δ3sin(δ4[t+4π/264]),(56)

where δ0,δ1,δ2,δ3,δ4 account for deterministic regularities in the spot price dynamics.

The following theorem presents the first main contribution of this paper.

Theorem 5.1 Suppose without seasonality (i.e.f=0), the commodity spot price z given by (55) satisfies the following stochastic differential equation

(57) dzt=κ(θlnzt)ztdSt+σztdBSt,(57)

where St is an independent, non-decreasing stable process with α(0,1). Then, the future price is

(58) F(T,ST)=exp12σ2+βα+βαγ(1eκT)+14κσ212σ2+βα+βα2(1e2κT)+12σ2+βα+βαSTeκT(58)

where γ:=θσ22κ, β[0,1] denotes the skewness parameter of the subordinator St.

Proof 5.2 See Appendix A

One of the advantages of the future prices above is that we can determine the price by estimating the distribution of the latent variable ST. Moreover, this latency could be observed as jumps, volatility or extreme events such as a tsunami, fire etc. However since two-factor models have proven over the years to provide better fits, we propose one in the framework of subordination. This way we can separately represent volatility and jumps or extreme events.

5.3. The two-factor commodity spot model

In the two-factor spot model, the volatility is modelled as a stochastic process while retaining jumps in the spot model. The future prices model is given by

(59) F(t,Xt,Vt)=f(t)+E[exp(Xτ+Vτ)],(59)

We present the second main contribution of the paper in the following theorem.

Theorem 5.3 Suppose without seasonality (i.e.f=0), the commodity spot price X satisfies the set of subordinated stochastic differential equation

(60) zt=κ(θzt)St+VtBSt(1)(60)
(61) Vt=λ(εVt)t+υVtBt(2),(61)

such that BS(1),B(2)t=ρAt where At=g(t,St) some random process. The future price is:

(62) F=exp(ψ0+ψ1zt+ψ2Vt),(62)

where the details of coefficients ψ0,ψ1 and ψ2 are given in Appendix A.

6. Numerical implementation

We focus on the one-factor model to explain our approach for estimating the model parameters. The data used in this section are obtained from the US Energy Information Administration and include future prices of Crude Oil (Light-Sweet, Cushing, Oklahoma) from 30 March 1983 to 6 December 2016 (8452 observations), Reformulated Regular Gasoline (New York Harbor) from 3 December 1984 to 31 October 2006 (5492 observations), Heating Oil (New York Harbor) from 2 January 1980 to 6 December 2016 (9262 observations) and Propane (Mont Belvieu, Texas) from 17 December 1993 to 18 September 2009 (3941 observations).

The parameters in the seasonality function (56) are estimated by fitting the function to the historical spot prices. The spot prices used include Crude Oil from 2 January 1986 to 12 December 2016 (7807 observations), RBOB Regular Gasoline from 11 March 2003 to 12 December 2016 (3460 observations), No. 2 Heating Oil from 2 June 1986 to 12 December 2016 (7683 observations) and Propane from 9 July 1992 to 12 December 2016 (6133 observations).

f(t)=δ0t+δ1sin(δ2[t+2π/264])+δ3sin(δ4[t+4π/264]).

The accuracy of the fitting in Figure depends on the choice of the initial parameters δ0, δ1, δ2, δ3 and δ4.

Figure 3. Seasonality is captured by the function f(t) defined in the following table. The best fit of f(t) can be obtained by obtaining an optimal set of the δ parameters.

Figure 3. Seasonality is captured by the function f(t) defined in the following table. The best fit of f(t) can be obtained by obtaining an optimal set of the δ parameters.

6.1. Equivalent latent regression model

The de-seasonalised one-factor future price given by (58) can be written as

(63) y=a+bx+ε,(63)

where y=lnF, x=St and a and b are given by

(64) a=12σ2+βα+βαγ(1eκT)+14κσ212σ2+βα+βα2(1e2κT)(64)
(65) b=12σ2+βα+βαeκT,(65)

and ε is an independent random error distributed as N(0,Θ) with zero mean and variance Θ.

Clearly, (63) belongs to the class of latent regression models since x=St is not observable.

There is literature where this kind of problem is handled using expectation-maximisation (EM) algorithms (see Dempster, Laird, and Rubin (Citation1977)) to estimate the model parameters. The latent variable x can be considered binary where in this case the EM algorithm would give estimates for a two-component normal mixture model. On the other hand, x can be allowed to be continuous between 0 and 1 with a beta distribution as in Tarpey and Petkova (Citation2010). The EM algorithm for estimating the model parameters in this case is more involved than for the two-component mixture model and more computationally challenging, but can be done nonetheless. Basically, the latent or unobserved x variables are imputed by their conditional expectation given the outcomes y. Our approach is adaptable to the latter approach through the Dynkin–Lamperti Theorem (see Gupta and Nadarajah (Citation2004)) where the unobserved variable follows a stable distribution defined on (,) with α(0,2] and the observable variable y represents the log-returns of the future prices. The algorithm is applied to the joint likelihood of the response y. We assume the error ε is independent of the latent predictor x. The joint density for x and y is given by

(66) h(x,y)=h(y|x;a,b,Θ)S(x;α,β,ν,μ).=N(y;a+bx,Θ)S(x;α,β,ν,μ),(66)

where S(x;α,β,ν,μ) is the α stable distribution, N(y;A,B) denotes the normal distribution of a random variable y with mean A and variance B; and Θ is the variance of the outcome sample data. The marginal density of the response y is

f(y)=12πΘ1/2exp(Θ1(yβ0β1x)Θ1(yβ0β1x)/2)
(67) ×S(x;α,β,ν,μ)dx.(67)

Density (67) is an example of infinite mixture models used in ecological statistics (Fox, Negrete-Yankelevich, and Sosa (Citation2015)).

6.2. The EM algorithm

For a data set (x1,y1),,(xn,yn) in (63). The log-likelihood is derived from (66) as

(68) L(α,β,ν,μ,Θ,a,b)=n2log(2π)n2logΘi=1n(yiβ0β1xi)(yiβ0β1xi)/2+i=1nlog(S(xi;α,β,ν,μ)).(68)

Since x is not observable, the EM algorithm requires maximising the conditional expectation of the log-likelihood given the response vector y. That is

(69) E[L(α,β,ν,μ,Θ,a,b)|y],(69)

where at each iteration of the EM algorithm, the above conditional expectation is computed using the current parameter estimates. This current expectation-maximisation problem is similar to the problem handled in Tarpey and Petkova (Citation2010), which implies that a similar EM algorithm can be applied here. That is, suppose ρ=ab is a 2×p dimension coefficient matrix where each of the p columns of ρ provides the intercept and slope regression coefficients for each of the ρ response variables. Then the design matrix is denoted by X whose first column consists of ones for the intercept and the second column consists of the latent predictors xi,i=1,,n. The multivariate regression model follows:

(70) Y=Xρ+ε.(70)

Moreover, and as indicated in Tarpey and Petkova (Citation2010), the likelihood for the multivariate normal regression model can be given as

(71) L(ρ,Θ)=(2π)np/2|Θ|1/2exp(tr[Θ1(YXρ)(YXρ)]/2).(71)

The EM approach requires that we maximise the expectation of the logarithm of (71) conditional on Y with respect to ρ and Θ. This leads to the following optimal factors:

(72) ρˆ=(X˜X)1X˜Y.(72)
(73) Θˆ=YYρˆ(X˜X)ρˆ,(73)

where X˜=E[X|Y] and (X˜X)=E[XX|Y].

To implement the EM method in the R programming language, we first highlight that there are a minor differences to bear in mind before implementing the algorithm as we explain in the following.

First and foremost, the density of the predicator in our case is from a stable distribution. Recall, in general, the densities of stable processes cannot be expressed analytically, which makes it difficult to compute the log-likelihood. However, with the help of inbuilt packages in R including stabledist and StableEstim, the log likelihood can be satisfactorily estimated using Estim() to obtain the stable parameters of St, and dstable() for its corresponding stable density.

Second, from the log-likelihood expression, we notice that we only require estimates of the conditional expectations of x, x2 and logx with respect to the joint probability density given the response vector y.

On the other hand, we retain some of the steps in Tarpey and Petkova (Citation2010). The initial values for the regression parameters a and b can be obtained from fitting a two-component finite mixture model or by a preliminary search over the parameter space. Initial values for Θ can be obtained using the sample covariance matrix from the raw data.

6.3. Data for the EM algorithm

The data used is stored in a data frame with three columns containing futures log-returns, spot price log-returns and binary data of 1’s and 0’s representing whether or not a jump has occurred within a given window size (see Table ). Table shows the estimated parameters as a result. The parameters were obtained from 5000 data points of crude oil log-returns arranged as in Table . We have displayed results from only two iterations because for large data sets the code tends to be slow in addition to suffering convergence issues. However, this can be improved and using faster machines.

Table 2. Estimation of parameters in the seasonality function

Table 3. Parameters obtained from maximum likelihood method

Table 4. A snapshot of the structure of the data used

The jump occurrence due to St is determined by the method discussed in Lee and Mykland (Citation2008) (also see Maslyuka et al., Citation2013). That is the realised return at any given time is compared with a continuously estimated instantaneous volatility σti to measure local variation arising from the continuous part of the process. The volatility σti is estimated using a modified version of realised bipower variation calculated as the sum of products of consecutive absolute returns in the local window (see Barndorff-Nielsen & Shephard, Citation2004). Then, the jump detection statistic Li testing for jumps in returns occurring at a time ti within a window size K is calculated as the ratio of realised returns to estimated instantaneous volatility:

(74) LilogYti/Yti1σˆti,(74)

where Yt at t0 represents the commodity spot price and σˆti is estimated by

(75) σˆti21K2j=iK+2i1log(Ytj/Ytj1)log(Ytj1/Ytj2).(75)

Care must be taken in choosing K, it must be large enough to accurately estimate integrated volatility but small enough for the variance to be approximately constant. In other words, K should be large enough but smaller than N, the number of observations so that the effect of jumps on estimating instantaneous volatility disappears. Some authors recommend K to be computed as K=252n, where n is the daily number of observations, whereas 252 is the number of days in the (financial) year. Moreover, the window size should be such that K=O(Δtλ) with 1<λ<0.5. For high-frequency data, Lee and Mykland (Citation2008) recommend, for returns sampled at frequencies of 60, 30, 15 and 5 minutes, the corresponding values of K to be 78, 110, 156 and 270. For our case, we shall choose K=4 for crude oil future prices with returns sampled daily.

6.3.1. Detection of jumps in the data

The test statistic L follows approximately a normal distribution when the data set has no jumps and its value becomes large otherwise. According to Lee and Mykland (Citation2008), the region for L is chosen based on the distribution of its maximum. For instance, suppose a particular interval (ti1,ti] has no jumps and the distance between two consecutive observations in this interval is small (i.e. Δ0), then the maximum should converge to the Gumbel variable:

(76) maxiAˉNLicNsNξ,(76)

where ξ has a cumulative distribution function P(ξx)=exp(ex), AˉN is the set of i1,2,,N such that there is no jump in (ti1,ti] and cN, sN are defined as

(77) cN=(2logN)1/20.8logπ+log(logN)1.6(2logN)1/2.(77)
(78) sN=10.8(2logN)1/2.(78)

The test is conducted by comparing the standardised maximum of Li in (76) to the critical values from the Gumbel distribution where the null hypothesis of no jump is rejected when the jump statistic is given by

(79) Li>G1(1λ)sN+cN,(79)

where G1(1λ) is the (1λ) quantile function of the standard Gumbel distribution. Suppose λ=0.1, then we reject the null hypothesis of no jump when Li>sNη+cN where η is such that exp(eη)=1η=0.9. That is η=log(log(0.9))=2.25. Figure shows a graph of jumps detected in crude oil future prices where we have used 1’s to record a jump occurrence and 0’s for no jump.

Figure 4. Detection of jumps in crude oil future prices.

Figure 4. Detection of jumps in crude oil future prices.

Figure 5. The large spikes in the jump statistic graph reflect the extreme events in the price where as the blue strips represent the smaller jumps that might go unnoticed. The former jumps are easily detected in returns yet the latter are not that visible.

Figure 5. The large spikes in the jump statistic graph reflect the extreme events in the price where as the blue strips represent the smaller jumps that might go unnoticed. The former jumps are easily detected in returns yet the latter are not that visible.

7. Summary

We have shown that the affine property is attainable and applicable to generalised spot models. We considered a stochastic differential equation with the source of randomness as subordinated Brownian motion as a specific example to derive the futures price. Moreover, it has been argued in some existing literature that the likelihood function exists in integrated form for models with singular noise meanwhile for cases of partially observed processes, a filtering technique is required. However, the work presented in this paper provided a new approach of pricing commodity futures for models with latent variables using the maximum expectation-maximisation, without using any filtering. Our approach is easy to implement once the joint probability density is established. The numerical implementation of the two factor model is left for future work.

Correction

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Acknowledgements

This work was supported by funds from the National Research Foundation of South Africa (NRF), the African Institute for Mathematical Sciences (AIMS) and the African Collaboration for Quantitative Finance and Risk Research (ACQuFRR), which is the research section of the African Institute of Financial Markets and Risk Management (AIFMRM), which delivers postgraduate education and training in financial markets, risk management and quantitative finance at the University of Cape Town in South Africa.

All the numerical results were generated using MATLAB and the Robust Analysis software package in R accessible at www.RobustAnalysis.com.

Additional information

Funding

This work was supported by the University of Cape Town [OAPF-713];

Notes on contributors

M. Kateregga

M. Kateregga holds a Ph.D. in Mathematical Finance from the University of Cape Town. Michael is a Software Developer at Mira Networks in South Africa. His main interests include financial markets research and functional programming particularly Haskell.

Notes

1. The process L is also interpreted as the first passage time of the stable process S.

References

  • Applebaum, D. (2004). L´evy processes and stochastic calculus. Cambridge studies in advanced mathematics. Cambridge: Cambridge University Press.
  • Barndorff-Nielsen, O. E., & Shephard, N. (2004). Power and bipower variation with stochastic volatility and jumps. Journal of Financial Econometrics, 2(1), 1–37. doi:10.1093/jjfinec/nbh001
  • Black, F., & Scholes, M. (1973). The pricing of options and corporate liabilities. Journal of Political Economy, 81(3), 637–654. doi:10.1086/260062
  • Bochner, S. (2012). Harmonic analysis and the theory of probability. Dover books on mathematics. Mineola, NY: Dover Publications.
  • Crosby, J. (2008). A multi-factor jump-diffusion model for commodities. Journal of Quantitative Finance, 8(2), 181–200. doi:10.1080/14697680701253021
  • Date, P., & Ponomareva, K. (2010). Linear and nonlinear filtering in mathematical finance: A review. IMA Journal of Management Mathematics, 22(3), 195–211. doi:10.1093/imaman/dpq008
  • Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the American Statistical Association, 39, 1–38.
  • Deng, S. (2000). "Stochastic models of energy commodity prices and their applications: Mean-reversion with jumps and spikes," Working Paper PWP-073, University of California Energy Institute.
  • Duffie, D., Filipovi´C, D., & Schachermayer, W. (2003). Affine processes and applications in finance. Annals of Applied Probability, 13, 984–1053. doi:10.1214/aoap/1060202833
  • Dumouchel, W. H. (1971). Stable distributions in statistical inference. The Journal of the American Statistical Association, 78(342), 469–477.
  • Fama, E. F., & French, K. R. (1988). Permanent and temporary components of stock prices. The Journal of Political Economy, 96(2), 246–273. doi:10.1086/261535
  • Fama, E. F., & Roll, R. (1968). Some properties of symmetric stable distributions. Journal of the American Statistical Association, 63(323), 817–836.
  • Fox, G., Negrete-Yankelevich, S., & Sosa, V. (2015). Ecological statistics: Contemporary theory and application. Oxford: Oxford University Press.
  • Fusai, G., & Roncoroni, A. (2007). Implementing models in quantitative finance: Methods and cases. Springer finance. Berlin, Heidelberg: Springer.
  • G´omez-Valle, L., Habibilashkary, Z., & Martnez-Rodr`ıguez, J. (2017). A new technique to estimate the risk-neutral processes jump-diffusion commodity futures models. Journal of Computational and Applied Mathematics, 309, 435–441. doi:10.1016/j.cam.2015.12.028
  • Gajda, J., & Wyłoman`Ska, A. (2012). Geometric Brownian motion with tempered stable waiting times. Journal of Statistical Physics 1, 148(2), 296–305. doi:10.1007/s10955-012-0537-3
  • Grigelionis, B., Kubilius, J., Paulauskas, V., Statulevicius, V., & Pragarauskas, H. (1999). Probability theory and mathematical statistics: Proceedings of the seventh vilnius conference (1998), Vilnius, Lithuania, 12–18 August, 1998. Probability theory and mathematical statistics. TEV.
  • Gupta, A., & Nadarajah, S. (2004). Handbook of beta distribution and its applications. Statistics: A series of textbooks and monographs. Abingdon: Taylor & Francis.
  • Hicks, J. R. (1939). Value and Capital. Cambridge: Oxford University Press.
  • Holt, D. R., & Crow, E. L. (1973). Tables and graphs of the stable probability density function. Journal of Research of the National Bureau of Standards, Section B, 77, 143–198. doi:10.6028/jres.077B.017
  • Kateregga, M., Mataramvura, S., & Taylor, D. (2017). Parameter estimation for stable distributions with application to commodity futures log-returns. Cogent Economics & Finance, 5(1), 1318813. doi:10.1080/23322039.2017.1318813
  • Keller-Ressel, M. (2008). Affine processes theory and applications in finance. ( PhD thesis), Technischen Universita¨t Wien, Vienna, Austria.
  • Keller-Ressel, M., Schachermayer, W., & Teichmann, J. (2011). Affine processes are regular. Probability Theory and Related Fields, 151(3), 591–611. doi:10.1007/s00440-010-0309-4
  • Keynes, J. M. (1930). A treatise on money, vol II. London, United Kingdom: Macmillan.
  • Kyriakou, I., Nomikos, N. K., Pouliasis, P. K., & Papapostolou, N. C. (2016). Affine-structure models and the pricing of energy commodity derivatives. European Financial Management, 22(5), 853–881. doi:10.1111/eufm.12071
  • Lee, S., & Mykland, P. A. (2008). Jumps in financial markets: A new nonparametric test and jump clustering. Review of Financial Studies, 21(6), 2535–2563.
  • Maslyuka, S., Rotarub, K., & Dokumentovc, A. (2013). Price discontinuities in energy spot and futures prices. Monash Economics Working Papers 33-13, Monash University, Department of Economics.
  • Meerschaert, M. M., & Straka, P. (2013). Inverse stable subordinators. Mathematical Modelling of Natural Phenomena, 8(2), 1–16. doi:10.1051/mmnp/20138201
  • Miller, J. F. (1951). Moment-generating functions and laplace transforms. Journal of the Arkansas Academy of Science, 4(1), 16.
  • Prokopczuk, M., Symeonidis, L., & Simen, C. W. (2015). Do jumps matter for volatility forecasting? Evidence from energy markets. Journal of Futures Markets, 20, 1–15.
  • Rachev, S. (2003). Handbook of heavy tailed distributions in finance: Handbooks in finance. Handbooks in Finance. Amsterdam: Elsevier Science.
  • Rouah, F., & Heston, S. (2015). The heston model and its extensions in VBA. Wiley Finance. Hoboken, NJ: Wiley.
  • Samoradnitsky, G., & Taqqu, M. (1994). Stable non-gaussian random processes: Stochastic models with infinite variance. Stochastic modeling series. United Kingdom: Taylor & Francis.
  • Schmitz, A., Wang, Z., & Kimn, J.-H. (2014). A jump diffusion model for agricultural commodities with bayesian analysis. Journal of Futures Markets, 34(3), 235–260. doi:10.1002/fut.21597
  • Schwartz, E. S. (1997). The stochastic behaviour of commodity prices: Implications for valuation and hedging. The Journal of Finance, 52(3), 923–973. doi:10.1111/j.1540-6261.1997.tb02721.x
  • Schwartz, E. S., & Smith, J. E. (2000). Short-term variations and long-term dynamics in commodity prices. Management Science, 46(7), 893–911. doi:10.1287/mnsc.46.7.893.12034
  • Tarpey, T., & Petkova, E. (2010). Latent regression analysis. Statistical Modelling, 10(2), 133–158. doi:10.1177/1471082X0801000202
  • Yang, J., Lin, B., Luk, W., & Nahar, T. (2014). Particle filtering-based maximum likelihood estimation for financial parameter estimation. In 2014 24th International Conference on Field Programmable Logic and Applications (FPL), 1–4.
  • Zolotarev, V. (1986). One-dimensional stable distributions. Translations of mathematical mono-graphs. Providence, Rhode Island, United States: American Mathematical Society.
  • Zolotarev, V. M. (1964). On the representation of stable laws by integrals. Trudy Matematicheskogo Instituta Imeni V.A. Steklova, 71, 46–50.
  • Zolotarev, V. M. (1980). Statistical estimates of the parameters of stable laws. Banach Center Publications, 6(1), 359–376. doi:10.4064/-6-1-359-376

Appendix A

Proof for Theorem 5.1

By applying Itô’s formula to Yt=lnzt, it is readily seen that

(80) dYt=κ(γYt)dSt+σdBSt,(80)

where γ:=θσ22κ and the future price with maturity date T is given by (see (see Schwartz (Citation1997))

(81) F(T,zT)=E[zT]=E[eYT].(81)

Theorem 4.5 suggests an explicit representation of (81) is attainable and it can be deduced by considering first, the continuous case E[eXt]. Suppose a continuous mean-reverting model given by

(82) dXt=κ(γXt)dt+σdBt,X0=x.(82)

The corresponding affine forms of the coefficients according to Theorem 4.5 yield:

(83) K0=κγ,K1=κ,K2=0,H0=H1=0,H2=σ.(83)

Since Xt is regular affine and σ is constant for all t[0,T], then we have

(84) E[eiuXt]=exp(ψ0(t,u)+ψ1(t,u)x+ψ2(t,u)σ),uiR=C,(84)

where the functions ψ0(t,u), ψ1(t,u) and ψ2(t,u) satisfy the set of Riccati equations:

(85) ψ0t=K0Tη+12ηTH0η=κγψ1;ψ0(0,u)=0.(85)
(86) ψ1t=K1Tη+12ηTH1η=κψ1;ψ1(0,u)=iu.(86)
(87) ψ2t=K2Tη+12ηTH2η=12σψ12;ψ2(0,u)=0.(87)

where ηT=(ψ1ψ2). The solution set to the system of Riccati equations is given by

(88) ψ1(t,u)=iueκt,(88)
(89) ψ0(t,u)=iuγ(1eκt).(89)
(90) ψ2(t,u)=σu24κ(1e2κt).(90)

Using (84) where u=i, one can easily deduce E[eXτ] leading to the price of a commodity future price under a continuous model framework. Capitalising on the affine nature of Yt and Theorem 4.6, we deduce the representation for E[eYτ] as:

(91) E[eiuYt]=exp(ψ0(t,m(u))+ψ1(t,m(u))St+ψ2(t,m(u))σ),(91)

where the volatility σ is a constant and the system of Riccati equations takes the form

(92) ψ0t(t,m(u))=κγψ1;ψ0(0,m(u))=0.(92)
(93) ψ1t(t,m(u))=κψ1;ψ1(0,m(u))=m(iu).(93)
(94) ψ2t(t,m(u))=12σψ12;ψ2(0,m(u))=0.(94)

Consequently, the solution set is directly deduced from (88) (90) to obtain

(95) ψ1(t,m(u))=m(iu)eκt,(95)
(96) ψ0(t,m(u))=m(iu)γ(1eκt),(96)
(97) ψ2(t,m(u))=14κσm(iu)2(1e2κt).(97)

Setting u:=i yields

(98) E[eYt]=exp(m(1)γ(1eκt)+14κσ2m(1)2(1e2κt)+m(1)Steκt).(98)

The required result follows by substituting the Lévy exponent m(1) from (41).

Justifications for theorem 5.3

(99) ψ1(t,u1,u2)=m(iu1)eκt.(99)

(100) ϕ2(t,u1,u2)=2κm(iu1)υ2eκtj=1djm(iu1)jejκt+If(t,u1)C(u1,u2)12υ20tIf(s,u1)ds.(100)

(101) ϕ0(t,u1,u2)=θm(iu1)(1eκt)+2κλεm(iu1)υ2j=1djm(iu1)j1κ(1+j)(1eκt(1+j)).(101)

(102) +0tIf(s,u1)C(u1,u2)12υ20sIf(τ,u1)dτds,(102)

where coefficients djj=1 satisfy

(103) dj+1=i=1j1djdj11j>1ρυκdj1j>0+υ24κ21j=0(j+1κλκ).(103)

Note that to ensure that the futures price is positive real, the values of uiC,i=1,2 have to be chosen carefully which in this case it could be ui=i.

The factor C(u1,u2) is defined as

(104) C(u1,u2)=expρυm(iu1)κ+2j=1di1+jm(iu2)j+1m(iu2)+2κm(iu1)υ2j=1djm(iu1)j.(104)

Finally, the integrating factor If is such that

(105) If(t,u1)=expλtρυm(iu1)κeκt+2eκtj=1djj+1m(iu1)j+1.(105)
0tIf(τ,u1)τ=If(t,ui)(λ+ρυm(iu1)eκt2eκtj=1djj+1m(iu1)j+1
(106) expρυm(iu1)κ+2j=1djj+1m(iu1)j+1(λ+ρυm(iu1)2j=1djj+1m(iu1)j+1).(106)

We provide the proof using the following proposition and subsequent lemmas.

Proposition 8.2 The mean and variance of the model (60)–(61) are given by

(107) μ:=κ(θXt)λ(εVt),σ:=Vt0ρυVtυ1ρ2Vt,H:=σσT=VtρυVtρυVtυ2Vt.(107)

Moreover, their affine forms can be given as linear models of both X and V:

(108) μ=K0+K1Xt+K2Vt,(108)
(109) H=H0+H1Xt+H2Vt,(109)

where

(110) K0=κθλε,K1=κ0,K2=0λ.(110)
(111) H0=0000,H1=0000,H2=1ρυρυυ2.(111)

As a consequence, we deduce the following system of Riccati equations:

(112) ψ0t=K0Tη+12ηTH0η=κθψ1+λεψ2,(112)
(113) ψ1t=K1Tη+12ηTH1η=κψ1(113)
(114) ψ2t=K2Tη+12ηTH2η=λψ2+12ψ12+ρυψ1ψ2+12υ2ψ22,(114)

where ηT=(ψ1ψ2) with initial conditions ψ0(0,m(u1),m(u2))=0, ψ1(0,m(u1),m(u2))=m(iu1), ψ2(0,m(u1),m(u2))=m(iu2). The solutions take the form:

(115) ψ1(t,u1,u2)=m(iu1)eκt.(115)
(116) ψ0(t,u1,u2)=θm(iu1)(1eκt)+λε0tψ2(s,u1,u2)ds.(116)

Proof 8.3 This follows from the applications of Theorems 4.5 and 4.6.

To obtain the solution ψ2 to the Riccati, Equation (114) is not trivial. However, a similar problem has been handled in Kyriakou et al. (Citation2016), see Lemmas 8.4 and 8.6.

Lemma 8.4 Consider Proposition 8.2 and let ζ(y),yC0 be such that it satisfies

(117) dζ(y)dy=ζ(y)2+(κλκyρυκ)ζ(y)+υ24κ2,(117)

then the solution to (114) can be expressed by

(118) χ(t,u1)=2κm(iu1)υ2eκtζ(m(iu1)eκt).(118)

Moreover, the general solution to (114) takes the form

(119) ψ2=χ+1ω,(119)

where ω satisfies

(120) ωt+(λ+ρυm(iu1)eκt+υ2χ)ω=12υ2,(120)

with the general solution given by

(121) ω(t)=C12υ20tIf(s)dsIf(t),(121)

where If is an integrating factor and C is the constant of integration.

Proof 8.5 Claim (118) is verified by differentiating with respect to t and relating it to (117) and (115):

(122) χt=λχ+12υ2χ2+ρυχψ1+12ψ12.(122)

Similarly, (119) is verified by substitution into (114) and relating it to (122) resulting into

(123) 1ω2ωt=λω+12υ22χω+1ω2+ρυψ1ω,(123)

from which (120) follows. The general solution to (120) is obtained using the integrating factor

If(t,u1):=expλtρυm(iu1)κeκt2κm(iu1)eκtζ(m(iu1)eκt)dt
(124) =expλtρυm(iu1)κeκt+20m(iu1)eκtζ(y)dy.(124)

Lemma 8.6 A representation of the solution ψ2 to (114) is given by

(125) ψ2(t,u1,u2)=2κm(iu1)υ2eκtζ(m(iu1)eκt)+If(t,u1)C(u1,u2)12υ20tIf(s,u1)ds,(125)

where the constant of integration C is determined by applying ψ2(0,u1,u2)=m(iu2):

(126) C(u1,u2)=exp(ρυm(iu1)κ+20m(iu1)ζ(y)dy)m(iu2)+2κm(iu1)υ2ζ(m(iu1)).(126)

Proof 8.7 The function ζ can be expressed in the form (see Kyriakou et al. (Citation2016)):

(127) ζ(y)=j=1djyj,(127)

Functions ψ2(t,u1,u2) and ψ0(t,u1,u2) in (125) and (116), respectively, can be re-written as

(128) ψ2(t,u1,u2)=2κm(iu1)υ2eκtj=1djm(iu1)jejκt+If(t,u1)C(u1,u2)12υ20tIf(s,u1)ds.(128)
ψ0(t,u1,u2)=θm(iu1)(1eκt)
(129) +2κλεm(iu1)υ2j=1djm(iu1)j1κ(1+j)(1eκt(1+j)).(129)
(130) +0tIf(s,u1)C(u1,u2)12υ20sIf(τ,u1)dτds,(130)

where the integrating factor introduced in (124) and its integral are given by

(131) If(t,u1)=expλtρυm(iu1)κeκt+2eκtj=1djj+1m(iu1)j+1.(131)
0tIf(τ,u1)dτ=If(t,ui)(λ+ρυm(iu1)eκt2eκtj=1djj+1m(iu1)j+1
(132) expρυm(iu1)κ+2j=1djj+1m(iu1)j+1(λ+ρυm(iu1)2j=1djj+1m(iu1)j+1).(132)

Finally, the constant of integration (126) can be re-written as

(133) C(u1,u2)=expρηm(iu1)κ+2j=1di1+jm(iu2)j+1m(iu2)+2κm(iu1)η2j=1djm(iu1)j.(133)

This completes the proof.