Subordinated affine structure models for commodity future prices

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.

Keywords:

• stable distributions
• affine structure models
• latent regression models

AMS Subject Classification:

• 62G05
• 62G07
• 62G32

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 (1973) for pricing options and corporate liabilities. Among the many models developed, the widely used and referenced study on commodities is the work by Schwartz (1997). 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 1.

Table 1. Various commodity jump models

Author(s) (Year)	Study	Result
G´omez-Valle and Martnez-Rodr`ıguez (2017)	Jump size distribution	More accurate than a diffusion model
Kateregga, Mataramvura, and Taylor (2017)	Jump distribution	Better estimation method than Maximum Likelihood
Kyriakou et al. (2016)	Jumps & Heston-type stochastic vol.	Evidence of jumps and stochastic vol. in oil
Prokopczuk, Symeonidis, and Simen (2015)	Characterisation of jumps	Jump importance to vol.
Schmitz, Wang, and Kimn (2014)	Stoch. vol., price jumps, seasonality & stoch. cost of carry	A model that capture all the four features.
Maslyuka, Rotarub, and Dokumentovc (2013)	Frequency of price discontinuities	Evidence of discontinuities
Deng (2000)	Multiple jumps, regime–switching & stoch. vol.	Prices of various commodities.

The current research builds on findings from Kyriakou, Nomikos, Pouliasis, and Papapostolou (2016) 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 (1988), 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 $\delta$) 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) $F_{t,T}=S_t{e}^{(r-c-\delta )(T-t)},$(1)

where $c$ accounts for the storage costs, $r$ is the periodically continuously compounded interest rate, $S_t$ 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 (1930) and Hicks (1939). It asserts that the future prices are given by the discounted (by the risk premium) expected future spot price:

(2) $F_{t,T}={\mathrm{E}}_t[S_T]e^{-r\gamma [T-t]},$(2)

where $\gamma$ is the risk premium and ${\mathrm{E}}_t[\cdot ]=\mathrm{E}[\cdot |{\mathcal{F}}_t]$, ${\mathcal{F}}_t$ 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 (1997), Schwartz and Smith (2000) and the jump models of Kyriakou et al. (2016).

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 $\alpha$-stable process, $\alpha \in (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 (2010) and Yang, Lin, Luk, and Nahar (2014). 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. (2017). A stable or $\alpha$-stable process, $\alpha \in (0,2]$, belongs to the general class of Lévy distributions. It has a limiting distribution with a definitive exponent parameter $\alpha$ that determines its shape. The following two definitions follow from Rachev (2003).

Definition 3.1 Let $X_1,X_2,\cdots ,X_n$ be independent and identically distributed random variables and suppose a random variable $S$ defined by

(3) $\frac{1}{a_n}\left(\sum _{i=1}^n{X}_i-b_n\right)\to S,$(3)

where “→” represents weak convergence in distribution, $a_n$ is a positive constant and $b_n$ is real. Then, $S$ is a stable process. The constants $a_n$ and $b_n$ 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 $a_n$ and $b_n$ do not necessarily need to be finite provides the generalised central limit theorem.

Definition 3.2 (Generalised Central Limit Theorem) Suppose $X_1,X_2\cdots$, denotes a sequence of identically distributed independent random variables from some arbitrary distribution and let sequences $a_n\in R$ and $b_n\in R^{+}$. Then, we define a sequence

(4) $Z_n:=\frac{1}{b_n}\left(\sum _{i=1}^n{X}_i-a_n\right)$(4)

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

(5) $\mathrm{P}(Z_n<x)\to H(x),n\to \mathrm{\infty },$(5)

where H(x) is some limiting distribution.

The traditional central limit theorem assumes finite mean $a:=\mathrm{E}[X_i]$ and finite variance ${\sigma }^2:=\mathrm{V}ar[X_i]$ and defines the sequence of sums

(6) $Z_n:=\frac{1}{\sigma \sqrt{n}}\left(\sum _{i=1}^n{X}_i-na\right),$(6)

such that the distribution functions of $Z_n$ weakly converge to $h^{sG}(x)$. That is

(7) $\mathrm{P}(x_1<Z_n<x_2)\to {\int }_{x_1}^{x_2}{h}^{sG}(x)dx,n\to \mathrm{\infty }$(7)

where $h^{sG}(x)$ denotes the standard Gaussian density

(8) $h^{sG}(x)=\frac{1}{\sqrt{2\pi }}exp(-x^2/2).$(8)

Suppose the identically distributed independent random variables $X_i$ equal to a positive constant $c$ almost surely and the sequences $a_n$ and $b_n$ in (4) are defined by $a_n=(n-1)c$ and $b_n=1$, then $Z_n$ is also equal to $c$ for all $n>0$ almost surely. In this case, the random variables $X_i$ are mutually independent and as a consequence, the limiting distribution for the sums $Z_n$ 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 (2008). 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 (1994) An $\alpha$-stable distribution is a four-parameter family of distributions denoted by $S(\alpha ,\beta ,\nu ,\mu )$ where

(1) $\alpha \in (0,2]$ is the characteristic exponent responsible for the tail of the distribution.

(2) $\beta \in [-1,1]$ is responsible for skewness.

(3) $\nu >0$ is the scale parameter (sometimes referred to as variance when $\alpha =2$).

(4) $\mu \in \mathbb{R}$ is the location (sometimes referred to as mean).

Densities of $\alpha$-stable distributions do not have closed-form representations except the Gaussian ($\alpha =2$), Cauchy ($\alpha =1$) and Inverse Gaussian or Pearson ($\alpha =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 (1964), 1980, Zolotarev (1986)). The distribution functions for the different $\alpha$ values have been tabulated in Dumouchel (1971), Fama and Roll (1968) and Holt and Crow (1973).

Definition 3.4 (Gajda and Wyłoman`Ska (2012)) Let $S_t$ and $L_t$ denote an α-stable process and its respective inverse. Then for $t\in [0,T]$, we define the process $L_t$ as¹

(9) $L_s:=\left\{\begin{array}{cc}inf\left\{t:S_t>s\right\}& \mathrm{i}\mathrm{f}s\in [0,S_t)\\ T& \mathrm{i}\mathrm{f}s=S_T.\end{array}\right.$(9)

$S_t$ and $L_t$ are non-decreasing and cádlág with their graphical representations given in Figure 1.

Figure 1. The top graph shows a plot of a stable process ${\mathit{S}}_{\mathit{t}}$ and the bottom graph shows its inverse process ${\mathit{L}}_{\mathit{t}}$ simulated using exponent parameter value, $\mathit{\alpha }\mathbf{=}\mathbf{0.8}$, plotted against time on the horizontal.

3.1. Density and characteristic functions

Let $(X_t,t\ge 0)$ denote a Lévy process. The characterisation of $X_t$ is usually given by the Lévy–Khintchine formula.

Definition 3.5 (Lévy–Khintchine) Applebaum (2004) Consider a Lévy process $X=(X_t{)}_{t\ge 0}$. There exists $b\in \mathbb{R}$ and $\sigma \ge 0$ such that the characteristic function of $X$ is

(10) $\mathrm{\Phi }(t):=\mathbb{E}[e^{itX}]=exp\left(itb-\frac{1}{2}{\sigma }^2{t}^2+{\int }_{\mathbb{R}-\left\{0\right\}}(e^{itx}-1-itx{\mathbb{I}}_{\left|x\right|<1})\mathrm{\Pi }(dx)\right),$(10)

where $\mathbb{I}$ is the indicator function and $\mathrm{\Pi }$ is a $\sigma$-finite measure satisfying the constraint

(11) ${\int }_{{\mathbb{R}}^d-\left\{0\right\}}min(1,|x{|}^2)\mathrm{\Pi }(dx)<\mathrm{\infty };\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{y}{\int }_{{\mathbb{R}}^d-\left\{0\right\}}\frac{|x{|}^2}{1+|x{|}^2}\mathrm{\Pi }(dx)<\mathrm{\infty }.$(11)

Definition 3.6 (The Lévy-Itô Decomposition) Applebaum (2004) If $X_t$ is a Lévy process, there exists $b\in R$, a Brownian motion $B_{\sigma }(t)$ with variance $\sigma \in R_{+}$ and an independent Poisson random measure $N$ on $R_{+}\times (R-\left\{0\right\})$ such that, for each $t\ge 0$,

(12) $X_t=bt+B_{\sigma }(t)+{\int }_{|x|<1}x\tilde{N}(t,dx)+{\int }_{|x|\ge 1}xN(t,dx),$(12)

where

(13) $b=\mathbb{E}\left[X_1-{\int }_{\left|x\right|\ge 1}xN(1,dx)\right].$(13)

To preserve the martingale property, the compound Poisson random measure is compensated as $\tilde{N}=N-t\lambda$ where $\lambda$ is a Lévy measure satisfying (11).

For process $S_t$, we require $\sigma =0$ in (10) or $B=0$ in (12) and the Lévy measure in (11) given by

(14) $\mathrm{\Pi }(x)=\frac{C}{|x{|}^{1+\alpha }}dx;C>0,$(14)

The characteristic function $\mathrm{\Phi }$ of $S_t$ is obtained using the domain of attraction of stable random variables (See Grigelionis, Kubilius, Paulauskas, Statulevicius, & Pragarauskas, 1999) and the Lévy–Khinchine representation formula (See Definition 3.5 or Applebaum (2004) for a detailed explanation), i.e.

(15) $\mathrm{\Phi }(\theta )=\mathrm{E}[exp(i\theta S)]=\left\{\begin{array}{cc}exp(-{\nu }^{\alpha }|\theta {|}^{\alpha }[1-i\beta \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\theta )tan(\frac{\pi \alpha }{2})]+i\mu \theta );& \mathrm{f}\mathrm{o}\mathrm{r}\alpha \ne 1.\\ exp(-\nu \left|\theta \right|[1+i\beta \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\theta )\frac{2}{\pi }log\left|\theta \right|]+i\mu \theta );& \mathrm{f}\mathrm{o}\mathrm{r}\alpha =1.\end{array}\right.$(15)

Using Fourier transformation, the density function of $S_t$ is given by

(16) $h_S(t,u)=\frac{1}{\pi }{\int }_0^{\mathrm{\infty }}{e}^{-iu\theta }\mathrm{\Phi }(\theta )d\theta .$(16)

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

(17) $S_u<t\iff L_t\ge u.$(17)

It follows that $F(t,u):=\mathrm{P}(S_u<t)=\mathrm{P}(L_t\ge u)={\int }_0^{\mathrm{\infty }}{h}_L(\tau ;t)\tau$ where $h_L(u,t)$ denotes the probability density function of $L_t$. Consequently

(18) $h_L(u,t)=-\frac{\mathrm{\partial }F(t,u)}{\mathrm{\partial }u}=-\frac{\mathrm{\partial }}{\mathrm{\partial }u}{\int }_{-\mathrm{\infty }}^t{h}_S(\tau ,u)d\tau .$(18)

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

(19) $h_S(t,u)=u^{-1/\alpha }h(t{u}^{-1/\alpha }),$(19)

where $h(\tau )$ is the density of a standard stable process with a Laplace transform $\tilde{h}(\tau )=exp(-{\tau }^{\alpha })$. This follows from the fact that $S_u$ has the same distribution as $u^{1/\alpha }{S}_1$. As a result, the density of the inverse stable process $L_t$ can be given in terms of the standard stable process by

(20) $h_L(u,t)=\frac{t}{\alpha }{u}^{-1-1/\alpha }h(t{u}^{-1/\alpha }).$(20)

Figure 2 shows density graphs of $S_t$ for different exponent parameter values.

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

3.2. Laplace transforms

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

(21) $\mathrm{E}[e^{-\tau X_u}]=e^{-u\varphi (\tau )},$(21)

where $\varphi$ is the Laplace exponent of $X_u$ known as the Bernstein function represented by

(22) $\varphi (\tau )=a+b\tau +{\int }_{(0,\mathrm{\infty })}(1-e^{-\tau x})\mathrm{\Pi }(dx).$(22)

where $a,b>0$ and $\mathrm{\Pi }$ is the Lévy measure on $(0,\mathrm{\infty })$ such that $\int \frac{x}{1+x}\mathrm{\Pi }(dx)<\mathrm{\infty }.$

The Laplace transform of the stable process $S_u$ is given by (see Meerschaert & Straka, 2013)

(23) ${\tilde{h}}_S(\tau ,u)={\int }_0^{\mathrm{\infty }}{e}^{-t\tau }{h}_S(t,u)dt=exp(-uC\mathrm{\Gamma }(1-\alpha ){\tau }^{\alpha })=exp(-u((\tau +\beta {)}^{\alpha }+{\beta }^{\alpha })).$(23)

where $0\le \beta \le 1$. For $C=\mathrm{\Gamma }(1-\alpha )$ (alternatively $\beta =0$), the Laplace transform simplifies to that of a standard stable process:

(24) ${\tilde{h}}_S(\tau ,u)=\mathrm{E}[e^{-\tau S_u}]=exp(-u{\tau }^{\alpha });0<\alpha <1.$(24)

The Laplace transform ${\tilde{h}}_L(u;\tau )$ of the inverse stable process $L_t$ is obtained from (18):

(25) $\begin{array}{rl}& {\tilde{h}}_L(u,\tau )=-\frac{\mathrm{\partial }}{\mathrm{\partial }u}({\tau }^{-1}exp(-u((\tau +\beta {)}^{\alpha }+{\beta }^{\alpha })))\\ & ={\tau }^{-1}((\tau +\beta {)}^{\alpha }+{\beta }^{\alpha })exp(-u((\tau +\beta {)}^{\alpha }+{\beta }^{\alpha })).\end{array}$(25)

where the Laplace transform of ${\int }_0^{t}f(y)y$ is ${\tau }^{-1}\tilde{f}(\tau )$ and $h_L(u,\tau ):=0$ for $l<0$ or $\tau <0$. Since (25) does not have the general form for a Laplace transform of a Lévy process, then $L_t$ 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 $M_u(\tau )$ and $\tilde{h}(\tau ,u)$ denote the respective moment-generating function and Laplace transform of a random variable then

(26) $M_u(\tau )=\tilde{h}(\tau ,-u)+\tilde{h}(-\tau ,u),$(26)

where

(27) $\tilde{h}(a,b)={\int }_0^{\mathrm{\infty }}{e}^{-ta}h(t,b)dt.$(27)

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

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

$M_u(\tau )=exp(u((\tau +\beta {)}^{\alpha }+{\beta }^{\alpha }))+exp(-u((-\tau +\beta {)}^{\alpha }+{\beta }^{\alpha })).$

(28) $\mathrm{E}[S_t]=M_u(0)=\alpha u{\beta }^{\alpha -1}[e^{2u{\beta }^{\alpha }}+e^{-2u{\beta }^{\alpha }}].$(28)

(29) $\mathrm{V}ar[S_t]=M_u^{{ }^{\prime \prime }}(0)-(M_u^{\prime }(0){)}^2.$(29)

(30) $=\alpha (\alpha -1)u{\beta }^{\alpha -2}[e^{2u{\beta }^{\alpha }}-e^{-2u{\beta }^{\alpha }}]+{\alpha }^2{u}^{\alpha }{\beta }^{2\alpha }[{\beta }^2{e}^{2u{\beta }^{\alpha }}+e^{-u{\beta }^{\alpha }}]$(30)

3.4. Subordination

Let $B_{S_t}$ denote subordinated Brownian motion where $S_t$ is an $\alpha$-stable process introduced before with $\alpha \in (0,1)$ and $B_{\cdot }$ represents standard Brownian motion with mean zero. We require $B_t$ and $S_t$ to be independent. To ensure a complete working environment, we introduce a probability space.

Let $(\mathrm{\Omega },\mathcal{F},\mathrm{P})=({\mathrm{\Omega }}^B\times {\mathrm{\Omega }}^S,{\mathcal{F}}^B\times {\mathcal{F}}^S,{\mu }^B\times {\mu }^S)$ denote a complete joint probability space endowed with a filtration $({\mathcal{F}}_t{)}_{t\ge 0}$ such that ${\mathcal{F}}_t={\mathcal{F}}_t^B\vee {\mathcal{F}}_t^S$ where ${\mathcal{F}}_t^B$ and ${\mathcal{F}}_t^S$ are filtrations generated by $B_t$ and $S_t$, respectively.

Definition 3.10 Let $B=(B_t,{\mathrm{P}}^x)=\sqrt{2}{B}_t^{\ast }$ where $(B_t^{\ast }{)}_{t\ge 0}$ is standard Brownian motion in $R$. The transition density $p(x,y;t)$ of $B$ is given by

(31) $p(x;t)=\frac{1}{2\sqrt{\pi t}}exp\left(-{\left|x\right|}^2/4t\right),t>0,x,y\in R.$(31)

The semi-group $(P_t:t\ge 0)$ of $B$ is given by

(32) $P_{t}f(x)={\mathrm{E}}_x[f(B_t)]={\int }_{R}p(x,y;t)f(y)dy,$(32)

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

(33) $\mathcal{G}f(x):=\underset{t\downarrow 0}{lim}\frac{{\mathrm{E}}_x[f(X)]-f(x)}{t}=\underset{t\downarrow 0}{lim}\frac{P_{t}f(x)-f(x)}{t}.$(33)

Definition 3.11 Suppose $Y_t:=B_{S_t},t\ge 0$ is a subordinated Brownian motion. Its’ semi-group $(Q_t:t\ge 0)$ is defined by

(34) $Q_{t}f(y,t)={\mathrm{E}}_y[f(Y_t,t)]={\mathrm{E}}_y[f(B_{S_t})]={\int }_0^{\mathrm{\infty }}{P}_{u}f(y,u)h_S(t,u)du.$(34)

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

(35) $q(y,t)={\int }_0^{\mathrm{\infty }}p(y,u)h_S(t,u)du.$(35)

Lemma 3.12 The mean and variance of $B_{S_t}$ are given by

(36) ${\mathrm{E}}_y[B_{S_t}]={\int }_0^{\mathrm{\infty }}{\mathrm{E}}_y[B_s]h_S(t,u)du=0.$(36)

(37) ${\mathrm{E}}_y[B_{S_t}^2]={\int }_0^{\mathrm{\infty }}{\mathrm{E}}_y[B_s^2]h_S(t,u)du={\int }_0^{\mathrm{\infty }}u{h}_s\left(t,u\right)du=E_y\left[S_t\right]$(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<{\tau }_1<\cdots <{\tau }_n<T$, where ${\tau }_i$ are the jump times of the process $B_{S_t}$, we observe that (36) and (37) hold for every interval $[{\tau }_i,{\tau }_{i+1})$. Thus, in the limits, their sums converge, respectively, to $0$ and ${\mathrm{E}}_x[S_t]$ on $[0,T]$.

Lemma 3.13 Let $X$ be a Lévy process with characteristic exponent $\mathrm{\Psi }$ and $S$ an independent subordinator with Laplace exponent $\mathrm{\Phi }$. Then, the subordinated process $X_{S_t}$ is a Lévy process with characteristic exponent

(38) $m(\cdot )=\mathrm{\Phi }(\mathrm{\Psi }(\cdot )).$(38)

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

It is known that any Lévy process $X_s$, $t<s\le T$ with drift $\mu$ is fully determined by its characteristic function given by (see Fusai & Roncoroni, 2007)

(39) $\mathrm{E}[e^{i\lambda X_s}]=e^{\mu \mathrm{\Delta }+\mathrm{\Psi }(\lambda )\mathrm{\Delta }}$(39)

where $\mathrm{\Delta }=s-t$, $\mu$ is the drift parameter and $\mathrm{\Psi }(\lambda )$ is the characteristic exponent. A typical example is Brownian motion whose characteristic function is given by

(40) $\mathrm{E}[e^{i\lambda B_s}]=e^{\mu \mathrm{\Delta }-\frac{1}{2}{\sigma }^2{\lambda }^2\mathrm{\Delta }},\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{\Psi }(\lambda )=-\frac{1}{2}{\sigma }^2{\lambda }^2.$(40)

The characteristic exponent of subordinated Brownian motion $B_{S_t}$ is deduced from (23), (38), (40):

(41) $m(u)={\left(-\frac{1}{2}{\sigma }^2{u}^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }.$(41)

4. Affine models

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

(1) $\mathcal{D}:=R_{\ge 0}^m\times R^n$.

(2) $\mathcal{U}:=\left\{u\in C^d:Re{u}_I\le 0,Re{u}_J=0\right\}$,

where $I:=\left\{1,\cdots ,m\right\}$, $J:=\left\{m+1,\cdots ,m+n\right\}$ and $M:=I\cup J=\left\{1,\cdots ,d\right\}$.

(3) $P_{t}f(x)={\mathrm{E}}^x[f(X_t)]$ for all $x\in \mathcal{D},t\ge 0$.

(4) $\mathcal{O}:=\left\{(t,u)\in R_{\ge 0}\times \mathcal{U}:P_s{f}_u(0)\ne 0\mathrm{\forall }s\in [0,t]\right\}$.

(5) $\mathcal{X}$ will denote a closed state space.

Definition 4.1 (Keller-Ressel, 2008) A process is stochastically continuous if for any sequence $t_n\to t$ in $R_{\ge 0}$, then the random variables $X_{t_n}$ converge to $X_t$ in probability with respect to $({\mathrm{P}}^x{)}_{x\in \mathcal{D}}$.

Definition 4.2 (Keller-Ressel, 2008) An affine process is a stochastically continuous time-homogeneous Markov process ($X_t,{\mathrm{P}}^x$) with a state space $\mathcal{D}$ where the characteristic function is an exponentially affine function of the state vector. In other words, there exist functions ${\psi }_0:R_{\ge 0}\times i{R}^d\to C$ and ${\psi }_1:R_{\ge 0}\times i{R}^d\to C^d$ such that

(42) ${\mathrm{E}}^x[e^{u{X}_t}]=exp({\psi }_0(t,u)+{\psi }_1(t,u)x),$(42)

for all $x\in \mathcal{D}$ and for all $(t,u)\in R_{\ge 0}\times iR$.

Definition 4.1 can be extended to $\mathcal{O}$ satisfying the following properties [Prop. 1.3, Keller-Ressel (2008)]:

(1) ${\psi }_0$ maps $\mathcal{O}$ to ${\mathrm{C}}_{-}$ where ${\mathrm{C}}_{-}:=\left\{u\in C:Reu\le 0\right\}$.

(2) ${\psi }_1$ maps $\mathcal{O}$ to $\mathcal{U}$.

(3) ${\psi }_0(0,u)=0$ and ${\psi }_1(0,u)=u$ for all $u\in \mathcal{U}$.

(4) ${\psi }_0$ and ${\psi }_1$ admit the `semi-flow property’:

● ${\psi }_0(t+s,u)={\psi }_0(t,u)+{\psi }_0(s,{\psi }_1(t,u)),$

● ${\psi }_1(t+s,u)={\psi }_1(s,{\psi }_1(t,u))$, for all $t,s\ge 0$ with $(t+s,u)\in \mathcal{O}$.

(5) ${\psi }_0$ and ${\psi }_1$ are jointly continuous on $\mathcal{O}$.

(6) With the remaining arguments fixed, $u_I\mapsto {\psi }_0(t,u)$ and $u_I\mapsto {\psi }_1(t,u)$ are analytic functions in $\left\{u_I:Re{u}_I<0;(t,u)\in \mathcal{O}\right\}$.

(7) Let $(t,u),(t,w)\in \mathcal{O}$ with $Reu\le Rew$. Then

● $Re{\psi }_0(t,u)\le {\psi }_1(t,Rew),$

● $Re{\psi }_1(t,u)\le {\psi }_1(t,Rew).$

Lemma 4.3 An affine process $(X_t{)}_{t\ge 0}$ is regular if the following right derivatives exist for all $u\in \mathcal{U}$ and are continuous at $u=0$:

(43) $F^{(1)}\left(u\right):\frac{\mathrm{\partial }{\psi }_0}{\mathrm{\partial }t}(t,u){\left|{ }_{t=0^{+}},F^{(2)}\left(u\right):\frac{\mathrm{\partial }{\psi }_1}{\mathrm{\partial }t}(t,u)\right|}_{t=0^{+}}$(43)

The regularity condition can be extended to $\mathcal{O}$ for which case the following Riccati equations hold:

(44) $\frac{\mathrm{\partial }{\psi }_0}{\mathrm{\partial }t}(t,u)=F^{(1)}({\psi }_1(t,u)),{\psi }_0(0,u)=0,$(44)

(45) $\frac{\mathrm{\partial }{\psi }_1}{\mathrm{\partial }t}(t,u)=F^{(2)}({\psi }_1(t,u)),{\psi }_1(0,u)=u.$(45)

Proof 4.4 See Keller-Ressel (2008); Keller-Ressel et al. (2011); Rouah and Heston (2015).

Functions ${\psi }_0$, ${\psi }_1$ can be characterised respectively by admissible sets of parameters $(a,b,c,\pi )$ and $(p,q,r,\mu )$ where $\pi$ and $\mu$ are measures. The details are given in Keller-Ressel (2008); Keller-Ressel et al. (2011). We are interested in the affine property of the solution to the SDE:

(46) $X_t=b(X_t)dt+\sigma (X_t)d{\mathcal{M}}_t,$(46)

where $b:\mathcal{X}\to R^d$ is continuous, $\sigma :\mathcal{X}\to R^{d\times d}$ is measurable such that the diffusion matrix $\sigma (x){\sigma }^T$ is continuous and ${\mathcal{M}}_t$ 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 $X_t$ is a regular affine solution to (46). Then $b$ and $\sigma$ can be expressed as:

(47) $b(x)=K_0+K_1{x}_1+\cdots +K_d{x}_d,K_i\in R^d$(47)

(48) $\sigma (x)\sigma (x{)}^T=H_0+H_{1}x+\cdots +H_d{x}_d,H_i\in R^d\times R^d,$(48)

where $i=0,\cdots d$. Moreover, the characteristic function of $X_t$ has a log-linear form

(49) $\mathrm{E}[e^{i{u}_1{X}_t^{(1)}+\cdots +i{u}_d{X}_t^{(d)}}]$(49)

(50) $=exp({\psi }_0(t,u_1,\cdots ,u_d)+{\psi }_1(t,u_1,\cdots ,u_d)x_0^{(1)}+\cdots +{\psi }_d(t,u_1,\cdots ,u_d)x_0^{(d)}).$(50)

where $u_i\in R$. The coefficients ${\psi }_i$ are obtained by solving the system of Riccati equations:

(51) $F^{(i)}(t,{\psi }_1,\cdots ,{\psi }_d)=\frac{\mathrm{\partial }{\psi }_i}{\mathrm{\partial }t}=K_i^T\eta +\frac{1}{2}{\eta }^T{H}_i\eta ,i=\left\{0,\cdots ,d\right\},$(51)

where ${\eta }^T=({\psi }_1\cdots {\psi }_d)$.

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

There is extensive literature on affine processes $X_t$ where $\mathcal{M}:=B_t$ or $\mathcal{M}:=B_t+{\sigma }^{-1}{\int }_0^t{\xi }_{s}s$ with ${\xi }_t$, a Poisson jump process. We are interested in the solution to

(52) $d{Y}_t=b(Y_t)d{S}_t+\sigma (Y_t)d{B}_{S_t}.$(52)

Another contribution follows in the following theorem. We show that $Y_t=X_{S_t}$ is affine in the following theorem with $d=1$.

Theorem 4.6 Let $(\mathrm{\Omega },\mathcal{F},{\mathrm{P}}^x,({\mathcal{F}}_t{)}_{t\ge 0})$ denote a joint probability space for $(S_t{)}_{t\ge 0}$, a non-decreasing affine process taking values in $\mathcal{D}$ and $(X_t{)}_{t\ge 0},X_0=x$, an independent Lévy process. Define a process $Y_t:=X_{S_t},Y_0=y$ with Lévy exponent $m(w)$ and suppose $(S_t{)}_{t\ge 0}$ is regular with functional characteristics $F^{(1)}(u),F^{(2)}(u)$. Then $(Y_t{)}_{t\ge 0}$ is regular affine with functional characteristics $F^{(1)}(m(w)),F_X^{(2)}(m(w))$ and $F_Y^{(2)}(m(w))=0$, $u,w\in i\mathrm{R}$ with the characteristic function given by

(53) ${\mathrm{E}}^s[e^{iu{Y}_t}]=exp({\psi }_0(t,m(w))+{\psi }_1(t,m(w))S_t),$(53)

for some functions ${\psi }_0$ and ${\psi }_1$.

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

(54) $\begin{array}{rl}& {\mathrm{E}}^s[e^{w{Y}_{\tau }}\left|{\mathcal{F}}_t]={\mathrm{E}}^s[e^{w{X}_{S_t}}]={\mathrm{E}}^s[{\mathrm{E}}^s[exp(w{X}_{S_t})\right|\sigma (S_s{)}_{0\le s\le t}]].\\ & ={\mathrm{E}}^s[exp(m(w)S_t)].\\ & =exp({\psi }_0(t,m(w))+{\psi }_1(t,m(w))S_t).\end{array}$(54)

The last equality follows from the affine property of $S_t$ (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 $\alpha$-stable process where $\alpha \in (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) $z_t=f(t)+exp(Y_t),$(55)

where $Y_t$ satisfies (52) and seasonality is defined according to Kyriakou et al. (2016), as

(56) $f(t)={\delta }_{0}t+{\delta }_{1}sin({\delta }_2[t+2\pi /264])+{\delta }_{3}sin({\delta }_4[t+4\pi /264]),$(56)

where ${\delta }_0,{\delta }_1,{\delta }_2,{\delta }_3,{\delta }_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) $d{z}_t=\kappa (\theta -ln{z}_t)z_{t}d{S}_t+\sigma z_{t}d{B}_{S_t},$(57)

where $S_t$ is an independent, non-decreasing stable process with $\alpha \in (0,1)$. Then, the future price is

(58) $\begin{array}{rl}& F(T,S_T)=exp\left(\left[{\left(-\frac{1}{2}{\sigma }^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }\right]\right.\gamma (1-e^{-\kappa T})+\frac{1}{4\kappa }{\sigma }^2{\left[{\left(-\frac{1}{2}{\sigma }^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }\right]}^2(1-e^{-2\kappa T})\\ & \left.+\left[{\left(-\frac{1}{2}{\sigma }^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }\right]S_T{e}^{-\kappa T}\right)\end{array}$(58)

where $\gamma :=\theta -\frac{{\sigma }^2}{2\kappa }$, $\beta \in [0,1]$ denotes the skewness parameter of the subordinator $S_t$.

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 $S_T$. 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,X_t,V_t)=f(t)+\mathrm{E}[exp(X_{\tau }+V_{\tau })],$(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) $z_t=\kappa (\theta -z_t)S_t+\sqrt{V_t}{B}_{S_t}^{(1)}$(60)

(61) $V_t=\lambda (\epsilon -V_t)t+\upsilon \sqrt{V_t}{B}_t^{(2)},$(61)

such that ${\left[B_{S_{\cdot }}^{(1)},B_{\cdot }^{(2)}\right]}_t=\rho A_t$ where $A_t=g(t,S_t)$ some random process. The future price is:

(62) $F=exp({\psi }_0+{\psi }_1{z}_t+{\psi }_2{V}_t),$(62)

where the details of coefficients ${\psi }_0,{\psi }_1$ and ${\psi }_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)={\delta }_{0}t+{\delta }_{1}sin({\delta }_2[t+2\pi /264])+{\delta }_{3}sin({\delta }_4[t+4\pi /264]).$

The accuracy of the fitting in Figure 3 depends on the choice of the initial parameters ${\delta }_0$, ${\delta }_1$, ${\delta }_2$, ${\delta }_3$ and ${\delta }_4$.

Figure 3. Seasonality is captured by the function $\mathit{f}(t)$ defined in the following table. The best fit of $\mathit{f}(t)$ can be obtained by obtaining an optimal set of the $\mathit{\delta }$ 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+\epsilon ,$(63)

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

(64) $a=\left[{\left(-\frac{1}{2}{\sigma }^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }\right]\gamma (1-e^{-\kappa T})+\frac{1}{4\kappa }{\sigma }^2{\left[{\left(-\frac{1}{2}{\sigma }^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }\right]}^2(1-e^{-2\kappa T})$(64)

(65) $b=\left[{\left(-\frac{1}{2}{\sigma }^2+\beta \right)}^{\alpha }+{\beta }^{\alpha }\right]e^{-\kappa T},$(65)

and $\epsilon$ is an independent random error distributed as $N(0,\mathrm{\Theta })$ with zero mean and variance $\mathrm{\Theta }$.

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

There is literature where this kind of problem is handled using expectation-maximisation (EM) algorithms (see Dempster, Laird, and Rubin (1977)) 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 (2010). 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 (2004)) where the unobserved variable follows a stable distribution defined on $(-\mathrm{\infty },\mathrm{\infty })$ with $\alpha \in (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 $\epsilon$ is independent of the latent predictor $x$. The joint density for $x$ and $y$ is given by

(66) $\begin{array}{rl}& h(x,y)=h(y|x;a,b,\mathrm{\Theta })S(x;\alpha ,\beta ,\nu ,\mu ).\\ & =N(y;a+bx,\mathrm{\Theta })S(x;\alpha ,\beta ,\nu ,\mu ),\end{array}$(66)

where $S(x;\alpha ,\beta ,\nu ,\mu )$ is the $\alpha$ stable distribution, $N(y;A,B)$ denotes the normal distribution of a random variable $y$ with mean $A$ and variance $B$; and $\mathrm{\Theta }$ is the variance of the outcome sample data. The marginal density of the response $y$ is

$f(y)=\int \frac{1}{\sqrt{2\pi }}{\mathrm{\Theta }}^{-1/2}exp(-{\mathrm{\Theta }}^{-1}(-y-{\beta }_0-{\beta }_{1}x{)}^{\prime }{\mathrm{\Theta }}^{-1}(y-{\beta }_0-{\beta }_{1}x)/2)$

(67) $\times S(x;\alpha ,\beta ,\nu ,\mu )dx.$(67)

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

6.2. The EM algorithm

For a data set $(x_1,y_1),\cdots ,(x_n,y_n)$ in (63). The log-likelihood is derived from (66) as

(68) $\begin{array}{rl}& L(\alpha ,\beta ,\nu ,\mu ,\mathrm{\Theta },a,b)=-\frac{n}{2}log(2\pi )-\frac{n}{2}log\left|\mathrm{\Theta }\right|\\ & -\sum _{i=1}^n(y_i-{\beta }_0-{\beta }_1{x}_i)\prime (y_i-{\beta }_0-{\beta }_1{x}_i)/2\\ & +\sum _{i=1}^{n}log(S(x_i;\alpha ,\beta ,\nu ,\mu )).\end{array}$(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) $\mathrm{E}[L(\alpha ,\beta ,\nu ,\mu ,\mathrm{\Theta },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 (2010), which implies that a similar EM algorithm can be applied here. That is, suppose $\rho =\left(\begin{array}{c}a\\ b\end{array}\right)$ is a $2\times p$ dimension coefficient matrix where each of the $p$ columns of $\rho$ provides the intercept and slope regression coefficients for each of the $\rho$ response variables. Then the design matrix is denoted by $\mathbf{X}$ whose first column consists of ones for the intercept and the second column consists of the latent predictors $x_i,i=1,\cdots ,n$. The multivariate regression model follows:

(70) $\mathbf{Y}=\mathbf{X}\rho +\epsilon .$(70)

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

(71) $L(\rho ,\mathrm{\Theta })=(2\pi {)}^{-np/2}|\mathrm{\Theta }{|}^{-1/2}exp(-tr[{\mathrm{\Theta }}^{-1}(\mathbf{Y}-\mathbf{X}\rho {)}^{\prime }(\mathbf{Y}-\mathbf{X}\rho )]/2).$(71)

The EM approach requires that we maximise the expectation of the logarithm of (71) conditional on $\mathbf{Y}$ with respect to $\rho$ and $\mathrm{\Theta }$. This leads to the following optimal factors:

(72) $\hat{\rho }=({\tilde{\mathbf{X}}}^{\prime }\mathbf{X}{)}^{-1}{\tilde{\mathbf{X}}}^{\prime }\mathbf{Y}.$(72)

(73) $\hat{\mathrm{\Theta }}={\mathbf{Y}}^{\prime }\mathbf{Y}-{\hat{\rho }}^{\prime }({\tilde{\mathbf{X}}}^{\prime }\mathbf{X})\hat{\rho },$(73)

where $\tilde{\mathbf{X}}=\mathrm{E}[\mathbf{X}|\mathbf{Y}]$ and $({\tilde{\mathbf{X}}}^{\prime }\mathbf{X})=\mathrm{E}[{\mathbf{X}}^{\prime }\mathbf{X}|\mathbf{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 $S_t$, 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$, $x^2$ 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 (2010). 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 $\mathrm{\Theta }$ 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 4). Table 2 shows the estimated parameters as a result. The parameters were obtained from 5000 data points of crude oil log-returns arranged as in Table 4. 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

Commodity	${\delta }_0$	${\delta }_1$	${\delta }_2$	${\delta }_3$	${\delta }_4$
Crude Oil	0.01097129	−0.04210297	13.29908652	18.43819425	6.28446063
Gasoline	0.001018775	−0.515953818	6.27995144	1.252457783	6.28435941
Heating Oil	0.000287309	0.159154792	12.57029401	−0.49422029	6.28182486
Propane	0.000243901	−0.252880065	12.56438012	0.071328671	6.286055367

Table 3. Parameters obtained from maximum likelihood method

Parameter Estimation
Parameter	1st Iteration	2nd Iteration
loglike	$30687$	$31133.2$
$a$	$7.74707e-05$	$1.17456e-05$
$b$	$0.00218499$	$0.00166723$
$\alpha$	$1.6605$	$1.6605$
$\beta$	$-0.0651915$	$-0.0651915$
$\nu$	$0.00576286$	$0.00576286$
$\mu$	$0.000415904$	$0.000415904$
$\mathrm{\Theta }$	$0.0106344$	$0.0110766$
$T$	$5000$ days	$5000$ days

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

Futures log-returns	Spot log-returns	Jump detection
−0.015646306	−0.020072054	0
0.0322302	0.047966668	0
0	−0.016401471	0
−0.001778683	0.001489862	0
−0.017271048	−0.019778648	1
−0.02611648	−0.022030684	1
−0.043103026	−0.033698961	0
0	−0.010387745	0
0.001808809	0.003969483	0
0.037999099	0.036829588	1
−0.030483308	−0.030057522	0
0.024152859	0.019704761	0
0.009625568	0.008368168	0

The jump occurrence due to $S_t$ is determined by the method discussed in Lee and Mykland (2008) (also see Maslyuka et al., 2013). That is the realised return at any given time is compared with a continuously estimated instantaneous volatility ${\sigma }_{t_i}$ to measure local variation arising from the continuous part of the process. The volatility ${\sigma }_{t_i}$ 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, 2004). Then, the jump detection statistic ${\mathcal{L}}_i\in$ testing for jumps in returns occurring at a time $t_i$ within a window size $K$ is calculated as the ratio of realised returns to estimated instantaneous volatility:

(74) ${\mathcal{L}}_i\equiv \frac{log{Y}_{t_i}/Y_{t_{i-1}}}{{\hat{\sigma }}_{t_i}},$(74)

where $Y_t$ at $t\ge 0$ represents the commodity spot price and ${\hat{\sigma }}_{t_i}$ is estimated by

(75) ${\hat{\sigma }}_{t_i}^2\equiv \frac{1}{K-2}\sum _{j=i-K+2}^{i-1}\left|log(Y_{t_j}/Y_{t_{j-1}})\right|\left|log(Y_{t_{j-1}}/Y_{t_{j-2}})\right|.$(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=\sqrt{252\cdot n}$, 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=\mathcal{O}(\mathrm{\Delta }{t}^{\lambda })$ with $-1<\lambda <-0.5$. For high-frequency data, Lee and Mykland (2008) 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 $\mathcal{L}$ follows approximately a normal distribution when the data set has no jumps and its value becomes large otherwise. According to Lee and Mykland (2008), the region for $\mathcal{L}$ is chosen based on the distribution of its maximum. For instance, suppose a particular interval $(t_{i-1},t_i]$ has no jumps and the distance between two consecutive observations in this interval is small (i.e. $\mathrm{\Delta }\to 0$), then the maximum should converge to the Gumbel variable:

(76) $\frac{\underset{\hspace{0.17em}}{{max}_{i\in {\bar{A}}_N}}\left|{\mathcal{L}}_i\right|-c_N}{s_N}\to \xi ,$(76)

where $\xi$ has a cumulative distribution function $\mathrm{P}(\xi \le x)=exp(-e^{-x})$, ${\bar{A}}_N$ is the set of $i\in \left\{1,2,\cdots ,N\right\}$ such that there is no jump in $(t_{i-1},t_i]$ and $c_N$, $s_N$ are defined as

(77) $c_N=\frac{{(2logN)}^{1/2}}{0.8}-\frac{log\pi +log(logN)}{1.6{(2logN)}^{1/2}}.$(77)

(78) $s_N=\frac{1}{0.8{(2logN)}^{1/2}}.$(78)

The test is conducted by comparing the standardised maximum of ${\mathcal{L}}_i$ 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) ${\mathcal{L}}_i>G^{-1}(1-\lambda )s_N+c_N,$(79)

where $G^{-1}(1-\lambda )$ is the $(1-\lambda )$ quantile function of the standard Gumbel distribution. Suppose $\lambda =0.1$, then we reject the null hypothesis of no jump when ${\mathcal{L}}_i>s_N{\eta }^{\ast }+c_N$ where ${\eta }^{\ast }$ is such that $exp(-e^{-{\eta }^{\ast }})=1-{\eta }^{\ast }=0.9$. That is ${\eta }^{\ast }=-log(-log(0.9))=2.25$. Figure 5 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 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$.

Appendix A

Proof for Theorem 5.1

By applying Itô’s formula to $Y_t=ln{z}_t$, it is readily seen that

(80) $d{Y}_t=\kappa (\gamma -Y_t)d{S}_t+\sigma d{B}_{S_t},$(80)

where $\gamma :=\theta -\frac{{\sigma }^2}{2\kappa }$ and the future price with maturity date $T$ is given by (see (see Schwartz (1997))

(81) $F(T,z_T)=\mathrm{E}[z_T]=\mathrm{E}[e^{Y_T}].$(81)

Theorem 4.5 suggests an explicit representation of (81) is attainable and it can be deduced by considering first, the continuous case $\mathrm{E}[e^{X_t}]$. Suppose a continuous mean-reverting model given by

(82) $d{X}_t=\kappa (\gamma -X_t)dt+\sigma d{B}_t,X_0=x.$(82)

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

(83) $K_0=\kappa \gamma ,K_1=-\kappa ,K_2=0,H_0=H_1=0,H_2=\sigma .$(83)

Since $X_t$ is regular affine and $\sigma$ is constant for all $t\in [0,T]$, then we have

(84) $\mathrm{E}[e^{iu{X}_t}]=exp({\psi }_0(t,u)+{\psi }_1(t,u)x+{\psi }_2(t,u)\sigma ),u\in iR=C,$(84)

where the functions ${\psi }_0(t,u)$, ${\psi }_1(t,u)$ and ${\psi }_2(t,u)$ satisfy the set of Riccati equations:

(85) $\frac{\mathrm{\partial }{\psi }_0}{\mathrm{\partial }t}=K_0^T\eta +\frac{1}{2}{\eta }^T{H}_0\eta =\kappa \gamma {\psi }_1;{\psi }_0(0,u)=0.$(85)

(86) $\frac{\mathrm{\partial }{\psi }_1}{\mathrm{\partial }t}=K_1^T\eta +\frac{1}{2}{\eta }^T{H}_1\eta =-\kappa {\psi }_1;{\psi }_1(0,u)=iu.$(86)

(87) $\frac{\mathrm{\partial }{\psi }_2}{\mathrm{\partial }t}=K_2^T\eta +\frac{1}{2}{\eta }^T{H}_2\eta =\frac{1}{2}\sigma {\psi }_1^2;{\psi }_2(0,u)=0.$(87)

where ${\eta }^T=({\psi }_1{\psi }_2)$. The solution set to the system of Riccati equations is given by

(88) ${\psi }_1(t,u)=iu{e}^{-\kappa t},$(88)

(89) ${\psi }_0(t,u)=iu\gamma (1-e^{-\kappa t}).$(89)

(90) ${\psi }_2(t,u)=\frac{\sigma u^2}{4\kappa }(1-e^{-2\kappa t}).$(90)

Using (84) where $u=-i$, one can easily deduce $\mathrm{E}[e^{X_{\tau }}]$ leading to the price of a commodity future price under a continuous model framework. Capitalising on the affine nature of $Y_t$ and Theorem 4.6, we deduce the representation for $\mathrm{E}[e^{Y_{\tau }}]$ as:

(91) $\mathrm{E}[e^{iu{Y}_t}]=exp({\psi }_0(t,m(u))+{\psi }_1(t,m(u))S_t+{\psi }_2(t,m(u))\sigma ),$(91)

where the volatility $\sigma$ is a constant and the system of Riccati equations takes the form

(92) $\frac{\mathrm{\partial }{\psi }_0}{\mathrm{\partial }t}(t,m(u))=\kappa \gamma {\psi }_1;{\psi }_0(0,m(u))=0.$(92)

(93) $\frac{\mathrm{\partial }{\psi }_1}{\mathrm{\partial }t}(t,m(u))=-\kappa {\psi }_1;{\psi }_1(0,m(u))=m(iu).$(93)

(94) $\frac{\mathrm{\partial }{\psi }_2}{\mathrm{\partial }t}(t,m(u))=\frac{1}{2}\sigma {\psi }_1^2;{\psi }_2(0,m(u))=0.$(94)

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

(95) ${\psi }_1(t,m(u))=m(iu)e^{-\kappa t},$(95)

(96) ${\psi }_0(t,m(u))=m(iu)\gamma (1-e^{-\kappa t}),$(96)

(97) ${\psi }_2(t,m(u))=\frac{1}{4\kappa }\sigma m(iu{)}^2(1-e^{-2\kappa t}).$(97)

Setting $u:=-i$ yields

(98) $\mathrm{E}[e^{Y_t}]=exp(m(1)\gamma (1-e^{-\kappa t})+\frac{1}{4\kappa }{\sigma }^{2}m(1{)}^2(1-e^{-2\kappa t})+m(1)S_t{e}^{-\kappa t}).$(98)

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

Justifications for theorem 5.3

(99) ${\psi }_1(t,u_1,u_2)=m(i{u}_1)e^{-\kappa t}.$(99)

(100) ${\varphi }_2(t,u_1,u_2)=-\frac{2\kappa m(i{u}_1)}{{\upsilon }^2}{e}^{-\kappa t}\sum _{j=1}^{\mathrm{\infty }}{d}_{j}m(i{u}_1{)}^j{e}^{-j\kappa t}+\frac{I_f(t,u_1)}{C(u_1,u_2)-\frac{1}{2}{\upsilon }^2\underset{0}{\overset{t}{\int }}{I}_f(s,u_1)ds}.$(100)

(101) $\begin{array}{rl}& {\varphi }_0(t,u_1,u_2)=\theta m(i{u}_1)(1-e^{-\kappa t})\\ & +\frac{-2\kappa \lambda \epsilon m(i{u}_1)}{{\upsilon }^2}\sum _{j=1}^{\mathrm{\infty }}{d}_{j}m(i{u}_1{)}^j\left(\frac{1}{\kappa (1+j)}(1-e^{-\kappa t(1+j)})\right).\end{array}$(101)

(102) $+{\int }_0^t\frac{I_f(s,u_1)}{C(u_1,u_2)-\frac{1}{2}{\upsilon }^2{\int }_0^s{I}_f(\tau ,u_1)d\tau }ds,$(102)

where coefficients ${\left\{d_j\right\}}_{j=1}^{\mathrm{\infty }}$ satisfy

(103) $d_{j+1}=\frac{\sum _{i=1}^{j-1}{d}_j{d}_{j-1}{1}_{j>1}-\frac{\rho \upsilon }{\kappa }{d}_j{1}_{j>0}+\frac{{\upsilon }^2}{4{\kappa }^2}{1}_{j=0}}{(j+1-\frac{\kappa -\lambda }{\kappa })}.$(103)

Note that to ensure that the futures price is positive real, the values of $u_i\in C,i=\left\{1,2\right\}$ have to be chosen carefully which in this case it could be $u_i=-i$.

The factor $C(u_1,u_2)$ is defined as

(104) $C(u_1,u_2)=\frac{exp\left(-\frac{\rho \upsilon m(i{u}_1)}{\kappa }+2\sum _{j=1}^{\mathrm{\infty }}\frac{d_i}{1+j}m{(i{u}_2)}^{j+1}\right)}{m(i{u}_2)+\frac{2\kappa m(i{u}_1)}{{\upsilon }^2}\sum _{j=1}^{\mathrm{\infty }}{d}_{j}m{(i{u}_1)}^j}.$(104)

Finally, the integrating factor $I_f$ is such that

(105) $I_f(t,u_1)=exp\left(-\lambda t-\frac{\rho \upsilon m(i{u}_1)}{\kappa }{e}^{-\kappa t}+2e^{-\kappa t}\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1}\right).$(105)

${\int }_0^t{I}_f(\tau ,u_1)\tau =\frac{I_f(t,u_i)}{(-\lambda +\rho \upsilon m(i{u}_1)e^{-\kappa t}-2e^{-\kappa t}\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1}}$

(106) $-\frac{exp\left(\frac{-\rho \upsilon m(i{u}_1)}{\kappa }+2\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1}\right)}{(-\lambda +\rho \upsilon m(i{u}_1)-2\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{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) $\mu :=\left(\begin{array}{c}\kappa (\theta -X_t)\\ \lambda (\epsilon -V_t)\end{array}\right),\sigma :=\left(\begin{array}{cc}\sqrt{V_t}& 0\\ \rho \upsilon \sqrt{V_t}& \upsilon \sqrt{1-{\rho }^2}\sqrt{V_t}\end{array}\right),H:=\sigma {\sigma }^T=\left(\begin{array}{cc}{V}_t& \rho \upsilon V_t\\ \rho \upsilon V_t& {\upsilon }^2{V}_t\end{array}\right).$(107)

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

(108) $\mu =K_0+K_1{X}_t+K_2{V}_t,$(108)

(109) $H=H_0+H_1{X}_t+H_2{V}_t,$(109)

where

(110) $K_0=\left(\begin{array}{c}\kappa \theta \\ \lambda \epsilon \end{array}\right),K_1=\left(\begin{array}{c}-\kappa \\ 0\end{array}\right),K_2=\left(\begin{array}{c}0\\ -\lambda \end{array}\right).$(110)

(111) $H_0=\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right),H_1=\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right),H_2=\left(\begin{array}{cc}1& \rho \upsilon \\ \rho \upsilon & {\upsilon }^2\end{array}\right).$(111)

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

(112) $\frac{\mathrm{\partial }{\psi }_0}{\mathrm{\partial }t}=K_0^T\eta +\frac{1}{2}{\eta }^T{H}_0\eta =\kappa \theta {\psi }_1+\lambda \epsilon {\psi }_2,$(112)

(113) $\frac{\mathrm{\partial }{\psi }_1}{\mathrm{\partial }t}=K_1^T\eta +\frac{1}{2}{\eta }^T{H}_1\eta =-\kappa {\psi }_1$(113)

(114) $\frac{\mathrm{\partial }{\psi }_2}{\mathrm{\partial }t}=K_2^T\eta +\frac{1}{2}{\eta }^T{H}_2\eta =-\lambda {\psi }_2+\frac{1}{2}{\psi }_1^2+\rho \upsilon {\psi }_1{\psi }_2+\frac{1}{2}{\upsilon }^2{\psi }_2^2,$(114)

where ${\eta }^T=({\psi }_1{\psi }_2)$ with initial conditions ${\psi }_0(0,m(u_1),m(u_2))=0$, ${\psi }_1(0,m(u_1),m(u_2))=m(i{u}_1)$, ${\psi }_2(0,m(u_1),m(u_2))=m(i{u}_2)$. The solutions take the form:

(115) ${\psi }_1(t,u_1,u_2)=m(i{u}_1)e^{-\kappa t}.$(115)

(116) ${\psi }_0(t,u_1,u_2)=\theta m(i{u}_1)(1-e^{-\kappa t})+\lambda \epsilon \underset{0}{\overset{t}{\int }}{\psi }_2(s,u_1,u_2)ds.$(116)

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

To obtain the solution ${\psi }_2$ to the Riccati, Equation (114) is not trivial. However, a similar problem has been handled in Kyriakou et al. (2016), see Lemmas 8.4 and 8.6.

Lemma 8.4 Consider Proposition 8.2 and let $\zeta (y),y\in C-\left\{0\right\}$ be such that it satisfies

(117) $\frac{d\zeta (y)}{dy}=\zeta (y{)}^2+(\frac{\kappa -\lambda }{\kappa y}-\frac{\rho \upsilon }{\kappa })\zeta (y)+\frac{{\upsilon }^2}{4{\kappa }^2},$(117)

then the solution to (114) can be expressed by

(118) $\chi (t,u_1)=-\frac{2\kappa m(i{u}_1)}{{\upsilon }^2}{e}^{-\kappa t}\zeta (m(i{u}_1)e^{-\kappa t}).$(118)

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

(119) ${\psi }_2=\chi +\frac{1}{\omega },$(119)

where $\omega$ satisfies

(120) $\frac{\mathrm{\partial }\omega }{\mathrm{\partial }t}+(-\lambda +\rho \upsilon m(i{u}_1)e^{-\kappa t}+{\upsilon }^2\chi )\omega =-\frac{1}{2}{\upsilon }^2,$(120)

with the general solution given by

(121) $\omega (t)=\frac{C-\frac{1}{2}{\upsilon }^2{\int }_0^t{I}_f(s)ds}{I_f(t)},$(121)

where $I_f$ 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) $\frac{\mathrm{\partial }\chi }{\mathrm{\partial }t}=-\lambda \chi +\frac{1}{2}{\upsilon }^2{\chi }^2+\rho \upsilon \chi {\psi }_1+\frac{1}{2}{\psi }_1^2.$(122)

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

(123) $-\frac{1}{{\omega }^2}\frac{\mathrm{\partial }\omega }{\mathrm{\partial }t}=-\frac{\lambda }{\omega }+\frac{1}{2}{\upsilon }^2\left(\frac{2\chi }{\omega }+\frac{1}{{\omega }^2}\right)+\frac{\rho \upsilon {\psi }_1}{\omega },$(123)

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

$I_f(t,u_1):=exp\left(-\lambda t-\frac{\rho \upsilon m(i{u}_1)}{\kappa }{e}^{-\kappa t}-2\kappa m(i{u}_1)\int e^{-\kappa t}\zeta (m(i{u}_1)e^{-\kappa t})dt\right)$

(124) $=exp\left(-\lambda t-\frac{\rho \upsilon m(i{u}_1)}{\kappa }{e}^{-\kappa t}+2{\int }_0^{m(i{u}_1)e^{-\kappa t}}\zeta (y)dy\right).$(124)

Lemma 8.6 A representation of the solution ${\psi }_2$ to (114) is given by

(125) ${\psi }_2(t,u_1,u_2)=-\frac{2\kappa m(i{u}_1)}{{\upsilon }^2}{e}^{-\kappa t}\zeta (m(i{u}_1)e^{-\kappa t})+\frac{I_f(t,u_1)}{C(u_1,u_2)-\frac{1}{2}{\upsilon }^2{\int }_0^t{I}_f(s,u_1)ds},$(125)

where the constant of integration $C$ is determined by applying ${\psi }_2(0,u_1,u_2)=m(i{u}_2)$:

(126) $C(u_1,u_2)=\frac{exp(-\frac{\rho \upsilon m(i{u}_1)}{\kappa }+2{\int }_0^{m(i{u}_1)}\zeta (y)dy)}{m(i{u}_2)+\frac{2\kappa m(i{u}_1)}{{\upsilon }^2}\zeta (m(i{u}_1))}.$(126)

Proof 8.7 The function $\zeta$ can be expressed in the form (see Kyriakou et al. (2016)):

(127) $\zeta (y)=\sum _{j=1}^{\mathrm{\infty }}{d}_j{y}^j,$(127)

Functions ${\psi }_2(t,u_1,u_2)$ and ${\psi }_0(t,u_1,u_2)$ in (125) and (116), respectively, can be re-written as

(128) ${\psi }_2(t,u_1,u_2)=-\frac{2\kappa m(i{u}_1)}{{\upsilon }^2}{e}^{-\kappa t}\sum _{j=1}^{\mathrm{\infty }}{d}_{j}m(i{u}_1{)}^j{e}^{-j\kappa t}+\frac{I_f(t,u_1)}{C(u_1,u_2)-\frac{1}{2}{\upsilon }^2{\int }_0^t{I}_f(s,u_1)ds}.$(128)

${\psi }_0(t,u_1,u_2)=\theta m(i{u}_1)(1-e^{-\kappa t})$

(129) $+\frac{-2\kappa \lambda \epsilon m(i{u}_1)}{{\upsilon }^2}\sum _{j=1}^{\mathrm{\infty }}{d}_{j}m(i{u}_1{)}^j\left(\frac{1}{\kappa (1+j)}(1-e^{-\kappa t(1+j)})\right).$(129)

(130) $+{\int }_0^t\frac{I_f(s,u_1)}{C(u_1,u_2)-\frac{1}{2}{\upsilon }^2{\int }_0^s{I}_f(\tau ,u_1)d\tau }ds,$(130)

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

(131) $I_f(t,u_1)=exp\left(-\lambda t-\frac{\rho \upsilon m(i{u}_1)}{\kappa }{e}^{-\kappa t}+2e^{-\kappa t}\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1}\right).$(131)

${\int }_0^t{I}_f(\tau ,u_1)d\tau =\frac{I_f(t,u_i)}{(-\lambda +\rho \upsilon m(i{u}_1)e^{-\kappa t}-2e^{-\kappa t}\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1}}$

(132) $-\frac{exp\left(\frac{-\rho \upsilon m(i{u}_1)}{\kappa }+2\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1}\right)}{(-\lambda +\rho \upsilon m(i{u}_1)-2\sum _{j=1}^{\mathrm{\infty }}\frac{d_j}{j+1}m{(i{u}_1)}^{j+1})}.$(132)

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

(133) $C(u_1,u_2)=\frac{exp\left(-\frac{\rho \eta m(i{u}_1)}{\kappa }+2\sum _{j=1}^{\mathrm{\infty }}\frac{d_i}{1+j}m{(i{u}_2)}^{j+1}\right)}{m(i{u}_2)+\frac{2\kappa m(i{u}_1)}{{\eta }^2}\sum _{j=1}^{\mathrm{\infty }}{d}_{j}m{(i{u}_1)}^j}.$(133)

This completes the proof.