Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

On correlated measurement errors in the Schwartz–Smith two-factor model

On correlated measurement errors in the Schwartz–Smith two-factor model 1IntroductionStochastic processes have been commonly used for pricing of commodity derivatives for almost 50 years. The risk-neutral pricing theory for commodity derivatives was first developed in [7], which has become known as the Black–Scholes–Merton framework, where the commodity spot price is represented as a geometric Brownian motion (GBM), and for further details, see [8] and [21]. The principles of Black–Scholes–Merton’s framework laid the foundation for asset pricing theory. Since then, many models were developed by considering a number of factors as stochastic processes, which reflect the specifics of the commodity market. The mean-reverting process, or the Ornstein–Uhlenbeck (O–U) process, is often used for pricing of commodity derivatives. For example, in the two-factor oil contingent claims pricing model in [15], a mean-reverting factor and GBM were employed for modelling of the convenience yield and correlated oil spot price, respectively.In the Schwartz–Smith two-factor model [23], the sum of a short-term and a long-term factors, incorporating short-term deviations and long-term equilibrium price level, respectively, is equal to the spot price of a commodity. The short-term factor is assumed to tend towards zero, as it reflects short-term variations in prices from temporary changes in demand, supply, and other current market conditions, which will be corrected as the market responds over time. In addition, it is assumed that the dynamics of the long-term factor follows a Brownian motion with drift, which reflects expected permanent changes in the equilibrium price level, which can be explained by the advancement in technology for production, or any regulatory changes. The spot price of a commodity is then used to price futures contracts of different maturities jointly, under the risk-neutral probability measure. Studies in [5] and [12] develop models under the Schwartz–Smith model framework assuming both latent factors to follow an O–U process, with an additional constraint, which remedies the parameter identification problem.In this article, we present the model that incorporates dependence between futures contracts with different maturities. The novelty of our approach includes the introduction of correlations between measurement errors of different futures contracts, as well as allowing for serial correlation in each marginal measurement error. The correlations of the measurement errors along with other unknown model parameters will be jointly estimated with state variables using the Kalman filter. For illustration, we use the daily prices of European Union Allowance (EUA) futures contracts from January 2017 to April 2021, which were obtained using the Macquarie University access to Refinitiv Datascope.The European Union Emission Trading System (EU-ETS) was launched in 2005, with its aim to reduce greenhouse gas emissions from a variety of different sectors, such as agriculture, aviation, energy, and manufacturing industries across registered European nations. The implementation of the system puts obligations on those sectors to surrender one unit of EUA in order to emit one tonne of CO2{{\rm{CO}}}_{2}or equivalent gases. The history of the EUA market is relatively short, compared with other classic commodity markets such as crude oil, metal, and gas. We choose the selected period to study the recent dynamics of the EUA futures market. The selected time period covers the second half of Phase III and the beginning of Phase IV of the EU-ETS. The EU-ETS initiation and Phase II data were used in the following studies, see [4,25]. Also, the study by [13] used intra-phase and inter-phase futures data and accommodated specifics of each type contract by continuous-time diffusion models with jumps.The remaining sections are organised as follows. Section 2 reviews previous studies on modifications of the Schwartz–Smith two-factor model used for pricing of commodity derivatives, and studies on different approaches to pricing of EUA derivatives. Section 3 presents the main model that deals with both serial correlations and inter-correlations in measurement errors of logarithms of futures prices or their logarithmic returns. In Section 4, the results of simulation study are summarised, where we validate our approach to estimation of the parameters and state variables in case when both inter-correlations and serial correlations between measurement errors of different contracts are present. In Section 5, we present the results of the calibration of the two proposed models relative to the extended Schwartz–Smith model using historical daily EUA futures prices. Section 6 concludes with overall discussion of results of this study. For reference, the detailed setup of the Kalman filtering procedure in the Schwartz–Smith two-factor modelling framework is presented in Appendix A.2Literature reviewSince 2000, many researchers worked on adjustments of the original Schwartz–Smith two-factor model. A multi-commodity model was proposed in [11], with dynamics of state variables following the Ornstein–Uhlenbeck process, and the work in [6] extended the study for pricing of crude oil futures contracts, modifying the two-factor model in [23] to allow the long-term factor to follow a mean-reverting process. Further, the parameter identification problem for the two-factor model has been discussed in [5]. Both futures prices and analysts’ spot price forecasts were incorporated in the Kalman filter by [12], facilitating the term structure modelling in the crude oil market, and the study in [10] extended the analysis in the copper market. A novel method for testing the predictability of futures prices was proposed by [20], where the Kalman filter procedure has been modified to incorporate heteroscedasticity of prices and to estimate time-varying risk premium. For pricing of agricultural commodities, the performance of the Schwartz–Smith two-factor model has been studied in [24], using Fourier series as a seasonal component. An attempt at estimating covariances of measurement errors was also made, using a parametrised function of the time to maturity, but claimed that a substantial improvement in the model could not be seen. The study in [2] extended the two-factor modelling framework by incorporating explanatory variables/regression structure into the drift terms of the latent factors. The three-factor model was studied in [14], where the study allowed the deterministic seasonal component in the volatility of latent factors and used a function of inverse inventory as the third-state variable in the model. A step function was used in [22] as a seasonal component for the calibration of commodity spot and futures prices in a general multi-factor model, and a multi-factor model of commodity futures has been developed with stochastic seasonality in [19]. Under the same setup as in [23] and [24], instead of optimising the sample likelihood function, the study in [16] proposed a different estimation method, so-called a two-step least-square estimation method, which involves minimising the sum of squared residuals from the state equation.For pricing of EUA futures, the non-compliance event in terms of the total normalised emission was considered, along with the level of penalty in [9]. They used the digital nature of the terminal allowance price as the basis for modelling of the spot price process, and hence pricing of European options on EUA futures. The study in [26] developed a bivariate model in state-space form for parameter estimation through the Kalman filter, using December-maturity futures contracts from 2005 to 2012. In a recent study by [3], they have evaluated the term structure of EUA futures prices and compared performances of a single-factor GBM model by [1] and the original Schwartz–Smith two-factor model.3Main modelIn this section, we introduce the modifications for modelling of the logarithmic prices and logarithmic returns of futures prices, incorporating serial correlation and inter-correlation in measurement errors of different contracts, within the Schwartz–Smith two-factor modelling framework.The risk-neutral dynamics of short-term and long-term factors, notated as χt{\chi }_{t}and ξt{\xi }_{t}at time tt, are expressed as the following stochastic differential equations: (1)dχt=(−κχt−λχ)dt+σχdWtχ∗,{\rm{d}}{\chi }_{t}=\left(-\kappa {\chi }_{t}-{\lambda }_{\chi }){\rm{d}}t+{\sigma }_{\chi }{\rm{d}}{W}_{t}^{\chi \ast },(2)dξt=(μξ−λξ−γξt)dt+σξdWtξ∗,{\rm{d}}{\xi }_{t}=\left({\mu }_{\xi }-{\lambda }_{\xi }-\gamma {\xi }_{t}){\rm{d}}t+{\sigma }_{\xi }{\rm{d}}{W}_{t}^{\xi \ast },where κ,γ>0\kappa ,\gamma \gt 0are the speed of the mean-reversion for χt{\chi }_{t}and ξt{\xi }_{t}, respectively. σχ,σξ>0{\sigma }_{\chi },{\sigma }_{\xi }\gt 0are instantaneous volatilities of two latent variables, and λχ{\lambda }_{\chi }and λξ{\lambda }_{\xi }are the risk premia adjustments for χt{\chi }_{t}and ξt{\xi }_{t}that appear after transforming the model from the real probability measure to the risk-neutral (note that, the risk-neutral process is used for deriving the futures price). Wtχ∗{W}_{t}^{\chi \ast }and Wtξ∗{W}_{t}^{\xi \ast }are correlated standard Brownian processes, and E[dWtχ∗dWtξ∗]=ρχξdt{\mathbb{E}}{[}{\rm{d}}{W}_{t}^{\chi \ast }{\rm{d}}{W}_{t}^{\xi \ast }]={\rho }_{\chi \xi }{\rm{d}}t, where ρχξ{\rho }_{\chi \xi }is the correlation coefficient of two stochastic processes.By setting up the pricing model as a linear state-space model, two latent variables are expressed in the state equation, and the relationship between state variables and futures prices is expressed in the measurement equation. Then, we implement the Kalman filter to estimate values of latent variables and the marginal likelihood function, which are used for the estimation of model parameters. The readers are referred to [5] for a detailed setup under the assumption of measurement errors being independent for each contract.3.1Correlations in measurement errorsConsider the following linear state-space form for the Schwartz–Smith model: (3)xt=c+Gxt−1+wt,{{\bf{x}}}_{t}={\bf{c}}+{\bf{G}}{{\bf{x}}}_{t-1}+{{\bf{w}}}_{t},(4)yt=dt+Btxt+ϕvt−1+εt,εt∼N(0,Vε),{{\bf{y}}}_{t}={{\bf{d}}}_{t}+{{\bf{B}}}_{t}{{\bf{x}}}_{t}+{\boldsymbol{\phi }}{{\bf{v}}}_{t-1}+{{\boldsymbol{\varepsilon }}}_{t},\hspace{1em}{{\boldsymbol{\varepsilon }}}_{t}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{{\bf{V}}}_{\varepsilon }),where, for t=1,2,…,nt=1,2,\ldots ,nand for NNcontracts with different maturities T1<T2<…<TN{T}_{1}\lt {T}_{2}\hspace{0.33em}\lt \ldots \lt {T}_{N}, and (5)xt=χtξt,c=0μξ(1−e−γΔt)γ,G=e−κΔt00e−γΔt,yt=(lnFt,T1,lnFt,T2,…,lnFt,TN)′,dt=(A(T1−t)+g(T1),A(T2−t)+g(T2),…,A(TN−t)+g(TN))′,Bt′=e−κ(T1−t)…e−κ(TN−t)e−γ(T1−t)…e−γ(TN−t),ϕ=diag(ϕ11,ϕ12,…,ϕ1N),\begin{array}{rcl}{{\bf{x}}}_{t}& =& \left(\begin{array}{c}{\chi }_{t}\\ {\xi }_{t}\end{array}\right),\hspace{1em}{\bf{c}}=\left(\begin{array}{c}0\\ \frac{{\mu }_{\xi }\left(1-{e}^{-\gamma \Delta t})}{\gamma }\end{array}\right),\hspace{1em}{\bf{G}}=\left(\begin{array}{cc}{e}^{-\kappa \Delta t}& 0\\ 0& {e}^{-\gamma \Delta t}\end{array}\right),\\ {{\bf{y}}}_{t}& =& (\mathrm{ln}{F}_{t,{T}_{1}},\mathrm{ln}{F}_{t,{T}_{2}},\hspace{-0.18em}\ldots ,\mathrm{ln}{F}_{t,{T}_{N}})^{\prime} ,\\ {{\bf{d}}}_{t}& =& (A\left({T}_{1}-t)+g\left({T}_{1}),A\left({T}_{2}-t)+g\left({T}_{2}),\hspace{-0.18em}\ldots ,A\left({T}_{N}-t)+g\left({T}_{N}))^{\prime} ,\\ {{\bf{B}}}_{t}^{^{\prime} }& =& \left(\begin{array}{ccc}{e}^{-\kappa \left({T}_{1}-t)}& \ldots & {e}^{-\kappa \left({T}_{N}-t)}\\ {e}^{-\gamma \left({T}_{1}-t)}& \ldots & {e}^{-\gamma \left({T}_{N}-t)}\end{array}\right),\\ {\boldsymbol{\phi }}& =& \hspace{0.1em}\text{diag}\hspace{0.1em}\left({\phi }_{11},{\phi }_{12},\hspace{-0.18em}\ldots ,{\phi }_{1N}),\end{array}and A(T−t)A\left(T-t)is defined as follows: A(T−t)=−λχκ(1−e−κ(T−t))+μξ−λξγ(1−e−γ(T−t))+121−e−2κ(T−t)2κσχ2+21−e−(κ+γ)(T−t)κ+γσχσξρχξ+1−e−2γ(T−t)2γσξ2.\begin{array}{rcl}A\left(T-t)& =& -\frac{{\lambda }_{\chi }}{\kappa }(1-{e}^{-\kappa \left(T-t)})+\frac{{\mu }_{\xi }-{\lambda }_{\xi }}{\gamma }(1-{e}^{-\gamma \left(T-t)})\\ & & +\frac{1}{2}\left(\frac{1-{e}^{-2\kappa \left(T-t)}}{2\kappa }{\sigma }_{\chi }^{2}+2\frac{1-{e}^{-(\kappa +\gamma )\left(T-t)}}{\kappa +\gamma }{\sigma }_{\chi }{\sigma }_{\xi }{\rho }_{\chi \xi }+\frac{1-{e}^{-2\gamma \left(T-t)}}{2\gamma }{\sigma }_{\xi }^{2}\right).\end{array}Here, ϕ{\boldsymbol{\phi }}is a diagonal matrix that consists of autoregressive (AR) coefficients for each marginal measurement error of different contracts, and Δt\Delta tis the time difference in years between t−1t-1and tt. We may generalise AR process in measurement errors vt{{\bf{v}}}_{t}with order pp, and introduce additional ϕm{{\boldsymbol{\phi }}}_{m}matrices for m=1,…,pm=1,\ldots ,pin (4). For state and measurement errors, denoted as wt{{\bf{w}}}_{t}and vt{{\bf{v}}}_{t}, we assume that they are independent of each other, and (6)wt∼N(0,W),vt∼N(0,V),{{\bf{w}}}_{t}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{\bf{W}}),\hspace{1em}{{\bf{v}}}_{t}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{\bf{V}}),where (7)W=1−e−2κΔt2κσχ21−e−(κ+γ)Δtκ+γσχσξρχξ1−e−(κ+γ)Δtκ+γσχσξρχξ1−e−2γΔt2γσξ2,V=s112s12⋯s1Ns12s222⋯s2N⋮⋮⋱⋮s1Ns2N⋯sNN2.\begin{array}{rcl}{\bf{W}}& =& \left(\begin{array}{cc}\frac{1-{e}^{-2\kappa \Delta t}}{2\kappa }{\sigma }_{\chi }^{2}& \frac{1-{e}^{-\left(\kappa +\gamma )\Delta t}}{\kappa +\gamma }{\sigma }_{\chi }{\sigma }_{\xi }{\rho }_{\chi \xi }\\ \frac{1-{e}^{-\left(\kappa +\gamma )\Delta t}}{\kappa +\gamma }{\sigma }_{\chi }{\sigma }_{\xi }{\rho }_{\chi \xi }& \frac{1-{e}^{-2\gamma \Delta t}}{2\gamma }{\sigma }_{\xi }^{2}\\ \end{array}\right),\\ {\bf{V}}& =& \left(\begin{array}{cccc}{s}_{11}^{2}& {s}_{12}& \cdots \hspace{0.33em}& {s}_{1N}\\ {s}_{12}& {s}_{22}^{2}& \cdots \hspace{0.33em}& {s}_{2N}\\ \vdots & \vdots & \ddots & \vdots \\ {s}_{1N}& {s}_{2N}& \cdots \hspace{0.33em}& {s}_{NN}^{2}\end{array}\right).\end{array}We denote sjj2{s}_{jj}^{2}to be variances of measurement errors for contract jj, and sjk{s}_{jk}to be covariances of measurement errors between contracts jjand kk, where (j,k)=1,…,N\left(j,k)=1,\ldots ,N, and j≠kj\ne k.For estimation of covariances of measurement errors, we follow the estimation method in [18], where we estimate correlation coefficients of measurement errors between different contracts, and convert them back to covariances. We use the estimation approach introduced in [17] for correlation coefficients, which is often used in credit risk modelling. Let zj{{\bf{z}}}_{j}be the normalised prices of the contract jj, so that the vector of prices consists of z=(z1,z2,…,zN){\bf{z}}=\left({{\bf{z}}}_{1},{{\bf{z}}}_{2},\hspace{-0.18em}\ldots ,{{\bf{z}}}_{N})for NNcontracts. We assume that (8)zj=ρjε0+1−ρj2εj,j=1,…,N,{{\bf{z}}}_{j}={\rho }_{j}{\varepsilon }_{0}+\sqrt{1-{\rho }_{j}^{2}}{\varepsilon }_{j},\hspace{0.33em}j=1,\ldots ,N,where ε=(ε0,…,εN)∼N(0,IN+1){\boldsymbol{\varepsilon }}=\left({\varepsilon }_{0},\hspace{-0.18em}\ldots ,{\varepsilon }_{N})\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{{\bf{I}}}_{N+1}), ε0{\varepsilon }_{0}is a systematic component, and ε1,…,εN{\varepsilon }_{1},\ldots ,{\varepsilon }_{N}are idiosyncratic components. Therefore, (9)Corr[(zj,zk)]=ρjρk,(j,k)=1,…,N,j≠k.{\rm{Corr}}{[}({{\bf{z}}}_{j},{{\bf{z}}}_{k})]={\rho }_{j}{\rho }_{k},\hspace{0.33em}\left(j,k)=1,\ldots ,N,\hspace{0.33em}j\ne k.In this setup, we have the following correlation matrix structure: (10)R=1ρ1ρ2ρ1ρ3⋯ρ1ρNρ1ρ21ρ2ρ3⋯ρ2ρNρ1ρ3ρ2ρ31⋯ρ3ρN⋮⋮⋮⋱⋮ρ1ρNρ2ρN⋯ρN−1ρN1.{\bf{R}}=\left(\begin{array}{ccccc}1& {\rho }_{1}{\rho }_{2}& {\rho }_{1}{\rho }_{3}& \cdots \hspace{0.33em}& {\rho }_{1}{\rho }_{N}\\ {\rho }_{1}{\rho }_{2}& 1& {\rho }_{2}{\rho }_{3}& \cdots \hspace{0.33em}& {\rho }_{2}{\rho }_{N}\\ {\rho }_{1}{\rho }_{3}& {\rho }_{2}{\rho }_{3}& 1& \cdots \hspace{0.33em}& {\rho }_{3}{\rho }_{N}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {\rho }_{1}{\rho }_{N}& {\rho }_{2}{\rho }_{N}& \cdots \hspace{0.33em}& {\rho }_{N-1}{\rho }_{N}& 1\\ \end{array}\right).Using V=D12RD−12{\bf{V}}={{\bf{D}}}^{\tfrac{1}{2}}{\bf{R}}{{\bf{D}}}^{-\tfrac{1}{2}}, where D=diag(s112,…,sNN2){\bf{D}}=\hspace{0.1em}\text{diag}\hspace{0.1em}\left({s}_{11}^{2},\ldots ,{s}_{NN}^{2})is a diagonal matrix that consists of volatilities of measurement errors, the estimation will involve estimating both D{\bf{D}}and R{\bf{R}}. Hence, the modified covariance matrix V{\bf{V}}is applied in the Kalman filter and in the estimation procedure.3.2Modelling of the logarithmic returnsIn this section, we develop the measurement equations for the logarithms of relative returns on futures prices. Since the logarithmic returns are differences of the logarithmic prices at time ttand t−1t-1, we set up a linear state space model in the following way.For t=1,2,…,n−1t=1,2,\ldots ,n-1, state and measurement equations are written as follows: (11)xtr=ctr+Gtrxt−1r+wtr,{{\bf{x}}}_{t}^{r}={{\bf{c}}}_{t}^{r}+{{\bf{G}}}_{t}^{r}{{\bf{x}}}_{t-1}^{r}+{{\bf{w}}}_{t}^{r},(12)rt=yt−yt−1=dtr+Btrxtr+vtr,{{\bf{r}}}_{t}={{\bf{y}}}_{t}-{{\bf{y}}}_{t-1}={{\bf{d}}}_{t}^{r}+{{\bf{B}}}_{t}^{r}{{\bf{x}}}_{t}^{r}+{{\bf{v}}}_{t}^{r},where (13)xtr=xtxt−1,cr=c0,Gr=G0I0,wr=W000,dtr=dt−dt−1,Btr=Bt−Bt−1,vtr=vt−vt−1,vtr∼N(0,Vr).\begin{array}{rcl}{{\bf{x}}}_{t}^{r}& =& \left(\begin{array}{c}{{\bf{x}}}_{t}\\ {{\bf{x}}}_{t-1}\end{array}\right),\hspace{1em}{{\bf{c}}}^{r}=\left(\begin{array}{c}{\bf{c}}\\ {\bf{0}}\end{array}\right),\hspace{1em}{{\bf{G}}}^{r}=\left(\begin{array}{cc}{\bf{G}}& {\bf{0}}\\ {\bf{I}}& {\bf{0}}\end{array}\right),\hspace{1em}{{\bf{w}}}^{r}=\left(\begin{array}{cc}{\bf{W}}& {\bf{0}}\\ {\bf{0}}& {\bf{0}}\end{array}\right),\\ {{\bf{d}}}_{t}^{r}& =& {{\bf{d}}}_{t}-{{\bf{d}}}_{t-1},\hspace{1em}{{\bf{B}}}_{t}^{r}=\left(\begin{array}{c}{{\bf{B}}}_{t}\\ -{{\bf{B}}}_{t-1}\end{array}\right),\\ {{\bf{v}}}_{t}^{r}& =& {{\bf{v}}}_{t}-{{\bf{v}}}_{t-1},\hspace{1em}{{\bf{v}}}_{t}^{r}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{{\bf{V}}}^{r}).\end{array}Constant vectors, transition matrices, and measurement errors are from the original model setup shown in (5) and (7), with I{\bf{I}}being the identity matrix. Note that Vr{{\bf{V}}}^{r}is the covariance matrix of the measurement errors of logarithm of returns, instead of logarithm of prices. If vtr{{\bf{v}}}_{t}^{r}follows the AR process, then we can also set our measurement errors similar to Section 3.1. We can proceed with the standard Kalman filter and maximise the likelihood function using new notations accordingly.3.2.1Parameter estimationThe unknown parameter set ψ=(κ,σχ,λχ,γ,μξ,σξ,λξ,ρχξ,V,ϕ)\psi =(\kappa ,{\sigma }_{\chi },{\lambda }_{\chi },\gamma ,{\mu }_{\xi },{\sigma }_{\xi },{\lambda }_{\xi },{\rho }_{\chi \xi },{\bf{V}},{\boldsymbol{\phi }})is estimated by optimising the log-likelihood function of y{\bf{y}}, the joint distribution of (y1,y2,…,yn)({{\bf{y}}}_{1},{{\bf{y}}}_{2},\hspace{-0.18em}\ldots ,{{\bf{y}}}_{n}), with respect to ψ\psi . The log-likelihood function is (14)ℓ(ψ;y)=∑t=1np(yt∣It−1),\ell \left(\psi ;\hspace{0.33em}{\bf{y}})=\mathop{\sum }\limits_{t=1}^{n}p({{\bf{y}}}_{t}| {{\mathfrak{I}}}_{t-1}),where p(yt∣It−1)p({{\bf{y}}}_{t}| {{\mathfrak{I}}}_{t-1})is the logarithm of the conditional probability density of yt{{\bf{y}}}_{t}given information available until time t−1t-1. Now, assuming that the prediction error et=yt−E[yt∣It−1]{{\bf{e}}}_{t}={{\bf{y}}}_{t}-{\mathbb{E}}{[}{{\bf{y}}}_{t}| {{\mathfrak{I}}}_{t-1}]follows a multivariate normal distribution, the log-likelihood function in (14) can be re-expressed as follows: (15)ℓ(ψ;y)=−Nn2ln(2π)−12∑t=1n[ln(det(Lt∣t−1))+et′Lt∣t−1−1et],\ell \left(\psi ;\hspace{0.33em}{\bf{y}})=-\frac{Nn}{2}\mathrm{ln}\left(2\pi )-\frac{1}{2}\mathop{\sum }\limits_{t=1}^{n}{[}\mathrm{ln}\left(\det \left({{\bf{L}}}_{t| t-1}))+{{\bf{e}}}_{t}^{^{\prime} }{{\bf{L}}}_{t| t-1}^{-1}{{\bf{e}}}_{t}],where Lt∣t−1=Cov(et∣It−1){{\bf{L}}}_{t| t-1}={\rm{Cov}}({{\bf{e}}}_{t}| {{\mathfrak{I}}}_{t-1}), and det(⋅)\hspace{0.1em}\text{det}\hspace{0.1em}(\cdot )denotes the determinant of a matrix. Then, parameter estimates are obtained by maximising (15) with respect to ψ\psi jointly. To obtain the quantity for the log-likelihood function, latent state vectors and their covariance matrices at each time ttneed to be estimated through the Kalman filter. [5] detected the parameter identification problem within the log-likelihood function in the Kalman filter, and hence, the constraint κ≥γ\kappa \ge \gamma is considered in the optimisation procedure.4Simulation studyIn this section, we perform a simulation study to validate the new approach described in Section 3.1. We focus on validating serial correlation and inter-correlation assumptions in measurement error considered under the Schwartz–Smith model framework, to see how our novel approach performs at estimating parameters and state variables.The steps for this simulation study are as follows. (1)Set the true parameter set ψ=(κ,σχ,λχ,γ,μξ,σξ,λξ,ρχξ,V,ϕ)\psi =(\kappa ,{\sigma }_{\chi },{\lambda }_{\chi },\gamma ,{\mu }_{\xi },{\sigma }_{\xi },{\lambda }_{\xi },{\rho }_{\chi \xi },{\bf{V}},{\boldsymbol{\phi }}). Then, obtain the state and measurement variables xt{{\bf{x}}}_{t}and yt{{\bf{y}}}_{t}by simulating error terms wt{{\bf{w}}}_{t}and vt{{\bf{v}}}_{t}. vt{{\bf{v}}}_{t}is assumed to follow an AR(1) process.(2)Choose appropriate initial values, and determine appropriate feasible bounds for each parameter.(3)Conduct the optimisation procedure through the Kalman filter, which involves jointly estimating state variables and parameter estimates. Obtain parameter estimates, and estimated values for state variables.We simulate for n=200n=200, 500, 1,000, 2,000, 5,000, with five different futures contracts maturing in 1, 2, 3, 4, and 5 months. The parameter estimates are presented in Tables 1–3, with standard errors computed using Monte Carlo approach, to assess the accuracy of parameter estimates. The MATLAB code for the simulation study is available at https://github.com/Junee1992/EUA_Futures_Pricing/tree/main/Serial-Correlation.Table 1Estimates of parameters indicated in top row. Their true values are given in the last row. nnis used for number of time-steps indicated; (standard errors are provided and were computed using 1,000 paths)nnκ\kappa σχ{\sigma }_{\chi }λχ{\lambda }_{\chi }γ\gamma μξ{\mu }_{\xi }σξ{\sigma }_{\xi }λξ{\lambda }_{\xi }ρχξ{\rho }_{\chi \xi }2002.93160.10700.07152.84281.48800.16050.0290−0.9935-0.9935(0.0245)(0.0114)(0.0022)(0.0226)(0.0118)(0.0120)(0.0020)(0.0231)5000.89830.12040.00030.00270.00158.9136×10−68.9136\times 1{0}^{-6}0.0750−0.9254-0.9254(0.0230)(0.0062)(0.0016)(0.0188)(0.0098)(0.0069)(0.0015)(0.0232)1,0001.90420.11740.00020.25420.11910.03570.01820.9999(0.0226)(0.0028)(0.0013)(0.0159)(0.0082)(0.0034)(0.0013)(0.0218)2,0002.99850.05740.00011.02450.49640.11410.00000.9519(0.0223)(0.0040)(0.0010)(0.0133)(0.0067)(0.0043)(0.0007)(0.0157)5,0002.25230.05600.00061.10600.53780.11650.00000.9065(0.0189)(0.0015)(0.0006)(0.0114)(0.0057)(0.0015)(0.0006)(0.0105)True20.10.0110.50.10.010.8Table 2Estimates of autoregressive coefficients of the measurement errors. (Standard errors are provided and were computed using 1,000 paths)nnϕ11{\phi }_{11}ϕ12{\phi }_{12}ϕ13{\phi }_{13}ϕ14{\phi }_{14}ϕ15{\phi }_{15}2000.88900.93150.91100.91680.8774(0.0100)(0.0084)(0.0070)(0.0064)(0.0061)5000.91350.88090.91440.86970.8999(0.0011)(0.0007)(0.0007)(0.0006)(0.0006)1,0000.91760.92300.90310.91270.8889(0.0006)(0.0004)(0.0004)(0.0004)(0.0004)2,0000.90800.88950.90390.88710.8937(0.0293)(0.0003)(0.0003)(0.0003)(0.0003)5,0000.90430.90190.90080.90640.8956(0.0002)(0.0002)(0.0002)(0.0002)(0.0002)True0.90.90.90.90.9Table 3Estimates of standard deviations and factors of correlation coefficients of the measurement errors. (standard errors are provided and were computed using 1,000 paths)nns1{s}_{1}s2{s}_{2}s3{s}_{3}s4{s}_{4}s5{s}_{5}ρ1{\rho }_{1}ρ2{\rho }_{2}ρ3{\rho }_{3}ρ4{\rho }_{4}ρ5{\rho }_{5}2000.01200.01280.01210.01160.01140.85180.88440.86080.85720.8664(0.0001)(0.0001)(0.0001)(0.0001)(0.0001)(0.0021)(0.0013)(0.0012)(0.0013)(0.0012)5000.01010.00980.01000.00960.00920.75690.79510.80130.78730.7693(4.67×10−54.67\times 1{0}^{-5})(4.01×10−54.01\times 1{0}^{-5})(3.69×10−53.69\times 1{0}^{-5})(3.35×10−53.35\times 1{0}^{-5})(3.00×10−53.00\times 1{0}^{-5})(0.0048)(0.0040)(0.0040)(0.0034)(0.0028)1,0000.00990.00990.01060.01050.01000.79610.80730.82120.81570.8066(2.83×10−52.83\times 1{0}^{-5})(2.41×10−52.41\times 1{0}^{-5})(2.17×10−52.17\times 1{0}^{-5})(1.93×10−51.93\times 1{0}^{-5})(1.74×10−51.74\times 1{0}^{-5})(0.0024)(0.0017)(0.0015)(0.0014)(0.0012)2,0000.00960.00960.01020.01010.01000.79810.78350.81200.79520.8142(1.34×10−41.34\times 1{0}^{-4})(2.46×10−52.46\times 1{0}^{-5})(3.28×10−53.28\times 1{0}^{-5})(3.29×10−53.29\times 1{0}^{-5})(3.67×10−53.67\times 1{0}^{-5})(0.0180)(0.0089)(0.0095)(0.0094)(0.0098)5,0000.00900.00910.00910.00940.00940.90430.90190.90080.90640.8956(9.748×10−69.748\times 1{0}^{-6})(8.258×10−68.258\times 1{0}^{-6})(7.605×10−67.605\times 1{0}^{-6})(6.934×10−66.934\times 1{0}^{-6})(6.277×10−66.277\times 1{0}^{-6})(0.0008)(0.0007)(0.0006)(0.0006)(0.0005)True0.010.010.010.010.010.80.80.80.80.8Overall, parameter estimates are quite close to their true parameters, except for κ,λχ\kappa ,{\lambda }_{\chi }and λξ{\lambda }_{\xi }. κ,γ\kappa ,\gamma and ρχξ{\rho }_{\chi \xi }tend to fluctuate for n≤1,000n\le \hspace{0.1em}\text{1,000}\hspace{0.1em}; however, it stabilises as the sample size increases. Estimated AR coefficients, volatilities, and correlation coefficients are close to the true values, indicating that the model is able to capture dependencies between measurement errors of contracts with different maturities as well as their serial correlations. The estimation errors of first eight parameters are summarised in Figure 1.Figure 1Parameter estimation error plots versus n=500n=500, 1,000, 2,000, 5,000.The simulated and estimated state variables are shown in Figure 2, along with mean absolute errors (MAE) calculated for two latent variables in each panel, showing that estimation of state variables improves as nnincreases.Figure 2Actual and estimated χt{\chi }_{t}and ξt{\xi }_{t}versus n=500n=500, 1,000, 2,000, 5,000, and MAE for two factors, starting from top-left to right-bottom.5Application to EUA futures dataGiven the increasing importance of reducing carbon emissions and the role of the EU-ETS as the largest emissions trading scheme in the world, we have applied our model to this market. We consider the historical daily prices of EUA futures contracts traded in the Intercontinental Exchange from January 30, 2017, to April 30, 2021, obtained from Refinitiv Datascope. This period covers the second half of Phase III and the beginning of Phase IV of the EU-ETS. We consider the first seven available contracts for model calibration.Figure 3 illustrates the historical daily futures price dynamics through the trading phases of the EU ETS (Phases I–IV). The graph clearly illustrates the necessity of developing sufficiently versatile models to be able to accommodate the intricacies of EUA futures price dynamics across the different phases, including Phase IV. In addition, Figure 4 shows the term structure of EUA futures prices, which matures in December annually. The main difference between commonly traded commodities and EUA is that futures curves tend to increase smoothly for as we increase the maturity of futures contracts over time.Figure 3Path of prices of first available EUA futures contract for the period from April 22, 2005, to April 30, 2021.Figure 4Term structure of EUA futures prices for seven available contracts from January 30, 2017, to April 30, 2021.In the two works by [23] and [5], measurement errors are assumed to be independent. However, based on the belief that the price movement for each contract must be correlated with other available contracts during the same period within the same commodity market, we assume that measurement errors have the full covariance matrix. After investigating the statistics of measurement errors for EUA future prices, we found that measurement errors of contracts with different maturities are highly correlated with each other, as close to 1. The correlation matrix of measurement errors of EUA future prices using Model 1 is as follows: R=10.99990.99990.99990.99950.99990.999510.99990.99980.99940.99990.999510.99970.99930.99990.999610.99920.99980.999410.99960.998210.99941.{\boldsymbol{R}}=\left(\begin{array}{ccccccc}1& 0.9999& 0.9999& 0.9999& 0.9995& 0.9999& 0.9995\\ & 1& 0.9999& 0.9998& 0.9994& 0.9999& 0.9995\\ & & 1& 0.9997& 0.9993& 0.9999& 0.9996\\ & & & 1& 0.9992& 0.9998& 0.9994\\ & & & & 1& 0.9996& 0.9982\\ & & & & & 1& 0.9994\\ & & & & & & 1\end{array}\right).We also observe that each series of measurement errors follow AR(1) process, with all AR coefficients being highly significant (with pp-value <0.0001\lt 0.0001). Hence, the price data of EUA futures contracts are suitable to test the model that is developed in this study.For comparative analysis, we use two different models. The following two models consider estimating inter-correlation between measurement errors in different settings: Model 1– Modelling of logarithmic returns;Model 2– Modelling of logarithmic returns with serially correlated measurement errors.For goodness-of-fit assessment, we present the performance of each model using root mean-squared error (RMSE). The results are summarised in Table 4 for futures contracts Cj{C}_{j}, j=1,2,…,7j=1,2,\ldots ,7. In both models, the logarithmic returns are converted to the logarithm of prices. In temrs of RMSE solely, Model 1 performs better than Model 2. However, using model selection criteria (Akaike and Bayesian information criterion), Model 2 is preferable to Model 1.Table 4Akaike (AIC), Bayesian information criterion (BIC), and RMSE for Models 1 and 2Model 1Model 2C1{C}_{1}0.02380.0281C2{C}_{2}0.02300.0264C3{C}_{3}0.02360.0262C4{C}_{4}0.02690.0288C5{C}_{5}0.04930.0507C6{C}_{6}0.06950.0703C7{C}_{7}0.08960.0900AIC−8.9906×104-8.9906\times 1{0}^{4}−8.9985×104-8.9985\times 1{0}^{4}BIC−8.9796×104-8.9796\times 1{0}^{4}−8.9840×104-8.9840\times 1{0}^{4}Next, we provide results of out-of-sample predictions using 30/50-day out-of-sample windows in Tables 5 and 6. In each setting, we considered four different scenarios, where we repeat parameter estimation and out-of-sample prediction for every 1 day, 5 days, 10 days, and 30 or 50 days. We use the deseasonalised price data from January 30, 2017, to February 18, 2021, n=1,040n=\hspace{0.1em}\text{1,040}\hspace{0.1em}business days.Table 5RMSE for 1, 5, and 10 days ahead predictions in 30-day out-of-sample window using 30/1/2017–18/2/2021 dataModel 1Model 21 day5 days10 days30 days1 day5 days10 days30 daysC1{C}_{1}0.04470.04480.04500.04740.04400.04310.04480.0467C2{C}_{2}0.04490.04510.04530.04840.04390.04350.04520.0472C3{C}_{3}0.04480.04500.04520.04890.04390.04360.04520.0473C4{C}_{4}0.04430.04550.04470.04880.04380.04340.04500.0471C5{C}_{5}0.04490.04510.04530.05060.04440.04400.04550.0476C6{C}_{6}0.04540.04570.04590.05060.04500.04470.04620.0483C7{C}_{7}0.04610.04640.04660.05150.04580.04550.04690.0490Table 6RMSE for 1, 5, and 10 days ahead predictions in 50-day out-of-sample window using 30/1/2017–18/2/2021 dataModel 1Model 21 day5 days10 days50 days1 day5 days10 days50 daysC1{C}_{1}0.06440.07110.13520.07330.07720.06510.15690.0732C2{C}_{2}0.06700.07210.13730.07490.07970.06680.15830.0741C3{C}_{3}0.06860.07320.13940.07650.08120.06860.15990.0751C4{C}_{4}0.06890.07310.14020.07670.08150.06920.16030.0749C5{C}_{5}0.07030.07420.14200.07810.08260.07070.16150.0759C6{C}_{6}0.07110.07470.14300.07880.08330.07170.16200.0764C7{C}_{7}0.07200.07540.14410.07950.08390.07260.16260.0769From Table 5, we found that, in general, Model 2 performed better than Model 1, although differences in RMSE were minimal. We used Diebold–Mariano test for detecting a significant difference in out-of-sample predictability between the two models. For 30-day out-of-sample predictions, by using all possible pairs of series of marginal predictions under the same hh-days-ahead predictions for h=1,5,10,30h=1,5,10,30, we had pp-values ≈0.31\approx 0.31. Hence, there is no significant difference in the two models’ out-of-sample predictions.In Table 6, we found that Model 1 performed better for out-of-sample predictions for h=1,10h=1,10day, whereas Model 2 obtained the lower RMSE for h=5,50h=5,50. The Diebold–Mariano test showed no significant difference in 50 days out-of-sample predictions between the two models again. This pattern persisted in other subsets of data from Phases III and IV.6ConclusionIn this article, we have presented the modified Schwartz–Smith two-factor model, which can be used for modelling of logarithmic returns of futures contracts, incorporating serial correlation and inter-correlation of measurement errors. We compared the two different models that use the logarithm of futures prices and the logarithmic returns, applying them to EUA futures price data from Phase III and Phase IV of EU-ETS.The simulation study has illustrated that our novel approach was able to jointly estimate both parameters and state variables when serial correlations and inter-correlations were present in measurement errors. We have illustrated that the parameter and state variables estimates converged as the sample size increases. The maximum likelihood method was used for the estimation of coefficients of Gaussian AR processes used for modelling of serial correlations in measurement errors.Finally, we discussed the results of the comparative analysis of the two models in the context of EUA futures data. The goodness-of-fit and out-of-sample performances in predicting futures prices were discussed for each model. Overall, models for logarithmic returns showed a good performance in terms of RMSE for out-of-sample predictions. The full Model 2, for logarithmic returns with serial correlations, performed better than its reduced-form Model 1 for calibration of data in terms of AIC and BIC, and for predicting for 30 days out-of-sample window.The empirical results of the study emphasised the necessity of considering both serial correlations and inter-correlations of measurement errors for modelling of futures prices. For parameter estimation in Model 2, we used the two-step approach. First, we obtained the parameter estimates similar to Model 1, then we proceeded with estimating parameters of AR processes used for modelling of serial correlations. By using the real data through the incorporation of serial correlations in the measurement errors, we showed that the proposed models for logarithmic return capture the price movement more reliably. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Dependence Modeling de Gruyter

On correlated measurement errors in the Schwartz–Smith two-factor model

Loading next page...
 
/lp/de-gruyter/on-correlated-measurement-errors-in-the-schwartz-smith-two-factor-c5N6RD8FdU
Publisher
de Gruyter
Copyright
© 2022 Jun S. Han et al., published by De Gruyter
ISSN
2300-2298
eISSN
2300-2298
DOI
10.1515/demo-2022-0106
Publisher site
See Article on Publisher Site

Abstract

1IntroductionStochastic processes have been commonly used for pricing of commodity derivatives for almost 50 years. The risk-neutral pricing theory for commodity derivatives was first developed in [7], which has become known as the Black–Scholes–Merton framework, where the commodity spot price is represented as a geometric Brownian motion (GBM), and for further details, see [8] and [21]. The principles of Black–Scholes–Merton’s framework laid the foundation for asset pricing theory. Since then, many models were developed by considering a number of factors as stochastic processes, which reflect the specifics of the commodity market. The mean-reverting process, or the Ornstein–Uhlenbeck (O–U) process, is often used for pricing of commodity derivatives. For example, in the two-factor oil contingent claims pricing model in [15], a mean-reverting factor and GBM were employed for modelling of the convenience yield and correlated oil spot price, respectively.In the Schwartz–Smith two-factor model [23], the sum of a short-term and a long-term factors, incorporating short-term deviations and long-term equilibrium price level, respectively, is equal to the spot price of a commodity. The short-term factor is assumed to tend towards zero, as it reflects short-term variations in prices from temporary changes in demand, supply, and other current market conditions, which will be corrected as the market responds over time. In addition, it is assumed that the dynamics of the long-term factor follows a Brownian motion with drift, which reflects expected permanent changes in the equilibrium price level, which can be explained by the advancement in technology for production, or any regulatory changes. The spot price of a commodity is then used to price futures contracts of different maturities jointly, under the risk-neutral probability measure. Studies in [5] and [12] develop models under the Schwartz–Smith model framework assuming both latent factors to follow an O–U process, with an additional constraint, which remedies the parameter identification problem.In this article, we present the model that incorporates dependence between futures contracts with different maturities. The novelty of our approach includes the introduction of correlations between measurement errors of different futures contracts, as well as allowing for serial correlation in each marginal measurement error. The correlations of the measurement errors along with other unknown model parameters will be jointly estimated with state variables using the Kalman filter. For illustration, we use the daily prices of European Union Allowance (EUA) futures contracts from January 2017 to April 2021, which were obtained using the Macquarie University access to Refinitiv Datascope.The European Union Emission Trading System (EU-ETS) was launched in 2005, with its aim to reduce greenhouse gas emissions from a variety of different sectors, such as agriculture, aviation, energy, and manufacturing industries across registered European nations. The implementation of the system puts obligations on those sectors to surrender one unit of EUA in order to emit one tonne of CO2{{\rm{CO}}}_{2}or equivalent gases. The history of the EUA market is relatively short, compared with other classic commodity markets such as crude oil, metal, and gas. We choose the selected period to study the recent dynamics of the EUA futures market. The selected time period covers the second half of Phase III and the beginning of Phase IV of the EU-ETS. The EU-ETS initiation and Phase II data were used in the following studies, see [4,25]. Also, the study by [13] used intra-phase and inter-phase futures data and accommodated specifics of each type contract by continuous-time diffusion models with jumps.The remaining sections are organised as follows. Section 2 reviews previous studies on modifications of the Schwartz–Smith two-factor model used for pricing of commodity derivatives, and studies on different approaches to pricing of EUA derivatives. Section 3 presents the main model that deals with both serial correlations and inter-correlations in measurement errors of logarithms of futures prices or their logarithmic returns. In Section 4, the results of simulation study are summarised, where we validate our approach to estimation of the parameters and state variables in case when both inter-correlations and serial correlations between measurement errors of different contracts are present. In Section 5, we present the results of the calibration of the two proposed models relative to the extended Schwartz–Smith model using historical daily EUA futures prices. Section 6 concludes with overall discussion of results of this study. For reference, the detailed setup of the Kalman filtering procedure in the Schwartz–Smith two-factor modelling framework is presented in Appendix A.2Literature reviewSince 2000, many researchers worked on adjustments of the original Schwartz–Smith two-factor model. A multi-commodity model was proposed in [11], with dynamics of state variables following the Ornstein–Uhlenbeck process, and the work in [6] extended the study for pricing of crude oil futures contracts, modifying the two-factor model in [23] to allow the long-term factor to follow a mean-reverting process. Further, the parameter identification problem for the two-factor model has been discussed in [5]. Both futures prices and analysts’ spot price forecasts were incorporated in the Kalman filter by [12], facilitating the term structure modelling in the crude oil market, and the study in [10] extended the analysis in the copper market. A novel method for testing the predictability of futures prices was proposed by [20], where the Kalman filter procedure has been modified to incorporate heteroscedasticity of prices and to estimate time-varying risk premium. For pricing of agricultural commodities, the performance of the Schwartz–Smith two-factor model has been studied in [24], using Fourier series as a seasonal component. An attempt at estimating covariances of measurement errors was also made, using a parametrised function of the time to maturity, but claimed that a substantial improvement in the model could not be seen. The study in [2] extended the two-factor modelling framework by incorporating explanatory variables/regression structure into the drift terms of the latent factors. The three-factor model was studied in [14], where the study allowed the deterministic seasonal component in the volatility of latent factors and used a function of inverse inventory as the third-state variable in the model. A step function was used in [22] as a seasonal component for the calibration of commodity spot and futures prices in a general multi-factor model, and a multi-factor model of commodity futures has been developed with stochastic seasonality in [19]. Under the same setup as in [23] and [24], instead of optimising the sample likelihood function, the study in [16] proposed a different estimation method, so-called a two-step least-square estimation method, which involves minimising the sum of squared residuals from the state equation.For pricing of EUA futures, the non-compliance event in terms of the total normalised emission was considered, along with the level of penalty in [9]. They used the digital nature of the terminal allowance price as the basis for modelling of the spot price process, and hence pricing of European options on EUA futures. The study in [26] developed a bivariate model in state-space form for parameter estimation through the Kalman filter, using December-maturity futures contracts from 2005 to 2012. In a recent study by [3], they have evaluated the term structure of EUA futures prices and compared performances of a single-factor GBM model by [1] and the original Schwartz–Smith two-factor model.3Main modelIn this section, we introduce the modifications for modelling of the logarithmic prices and logarithmic returns of futures prices, incorporating serial correlation and inter-correlation in measurement errors of different contracts, within the Schwartz–Smith two-factor modelling framework.The risk-neutral dynamics of short-term and long-term factors, notated as χt{\chi }_{t}and ξt{\xi }_{t}at time tt, are expressed as the following stochastic differential equations: (1)dχt=(−κχt−λχ)dt+σχdWtχ∗,{\rm{d}}{\chi }_{t}=\left(-\kappa {\chi }_{t}-{\lambda }_{\chi }){\rm{d}}t+{\sigma }_{\chi }{\rm{d}}{W}_{t}^{\chi \ast },(2)dξt=(μξ−λξ−γξt)dt+σξdWtξ∗,{\rm{d}}{\xi }_{t}=\left({\mu }_{\xi }-{\lambda }_{\xi }-\gamma {\xi }_{t}){\rm{d}}t+{\sigma }_{\xi }{\rm{d}}{W}_{t}^{\xi \ast },where κ,γ>0\kappa ,\gamma \gt 0are the speed of the mean-reversion for χt{\chi }_{t}and ξt{\xi }_{t}, respectively. σχ,σξ>0{\sigma }_{\chi },{\sigma }_{\xi }\gt 0are instantaneous volatilities of two latent variables, and λχ{\lambda }_{\chi }and λξ{\lambda }_{\xi }are the risk premia adjustments for χt{\chi }_{t}and ξt{\xi }_{t}that appear after transforming the model from the real probability measure to the risk-neutral (note that, the risk-neutral process is used for deriving the futures price). Wtχ∗{W}_{t}^{\chi \ast }and Wtξ∗{W}_{t}^{\xi \ast }are correlated standard Brownian processes, and E[dWtχ∗dWtξ∗]=ρχξdt{\mathbb{E}}{[}{\rm{d}}{W}_{t}^{\chi \ast }{\rm{d}}{W}_{t}^{\xi \ast }]={\rho }_{\chi \xi }{\rm{d}}t, where ρχξ{\rho }_{\chi \xi }is the correlation coefficient of two stochastic processes.By setting up the pricing model as a linear state-space model, two latent variables are expressed in the state equation, and the relationship between state variables and futures prices is expressed in the measurement equation. Then, we implement the Kalman filter to estimate values of latent variables and the marginal likelihood function, which are used for the estimation of model parameters. The readers are referred to [5] for a detailed setup under the assumption of measurement errors being independent for each contract.3.1Correlations in measurement errorsConsider the following linear state-space form for the Schwartz–Smith model: (3)xt=c+Gxt−1+wt,{{\bf{x}}}_{t}={\bf{c}}+{\bf{G}}{{\bf{x}}}_{t-1}+{{\bf{w}}}_{t},(4)yt=dt+Btxt+ϕvt−1+εt,εt∼N(0,Vε),{{\bf{y}}}_{t}={{\bf{d}}}_{t}+{{\bf{B}}}_{t}{{\bf{x}}}_{t}+{\boldsymbol{\phi }}{{\bf{v}}}_{t-1}+{{\boldsymbol{\varepsilon }}}_{t},\hspace{1em}{{\boldsymbol{\varepsilon }}}_{t}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{{\bf{V}}}_{\varepsilon }),where, for t=1,2,…,nt=1,2,\ldots ,nand for NNcontracts with different maturities T1<T2<…<TN{T}_{1}\lt {T}_{2}\hspace{0.33em}\lt \ldots \lt {T}_{N}, and (5)xt=χtξt,c=0μξ(1−e−γΔt)γ,G=e−κΔt00e−γΔt,yt=(lnFt,T1,lnFt,T2,…,lnFt,TN)′,dt=(A(T1−t)+g(T1),A(T2−t)+g(T2),…,A(TN−t)+g(TN))′,Bt′=e−κ(T1−t)…e−κ(TN−t)e−γ(T1−t)…e−γ(TN−t),ϕ=diag(ϕ11,ϕ12,…,ϕ1N),\begin{array}{rcl}{{\bf{x}}}_{t}& =& \left(\begin{array}{c}{\chi }_{t}\\ {\xi }_{t}\end{array}\right),\hspace{1em}{\bf{c}}=\left(\begin{array}{c}0\\ \frac{{\mu }_{\xi }\left(1-{e}^{-\gamma \Delta t})}{\gamma }\end{array}\right),\hspace{1em}{\bf{G}}=\left(\begin{array}{cc}{e}^{-\kappa \Delta t}& 0\\ 0& {e}^{-\gamma \Delta t}\end{array}\right),\\ {{\bf{y}}}_{t}& =& (\mathrm{ln}{F}_{t,{T}_{1}},\mathrm{ln}{F}_{t,{T}_{2}},\hspace{-0.18em}\ldots ,\mathrm{ln}{F}_{t,{T}_{N}})^{\prime} ,\\ {{\bf{d}}}_{t}& =& (A\left({T}_{1}-t)+g\left({T}_{1}),A\left({T}_{2}-t)+g\left({T}_{2}),\hspace{-0.18em}\ldots ,A\left({T}_{N}-t)+g\left({T}_{N}))^{\prime} ,\\ {{\bf{B}}}_{t}^{^{\prime} }& =& \left(\begin{array}{ccc}{e}^{-\kappa \left({T}_{1}-t)}& \ldots & {e}^{-\kappa \left({T}_{N}-t)}\\ {e}^{-\gamma \left({T}_{1}-t)}& \ldots & {e}^{-\gamma \left({T}_{N}-t)}\end{array}\right),\\ {\boldsymbol{\phi }}& =& \hspace{0.1em}\text{diag}\hspace{0.1em}\left({\phi }_{11},{\phi }_{12},\hspace{-0.18em}\ldots ,{\phi }_{1N}),\end{array}and A(T−t)A\left(T-t)is defined as follows: A(T−t)=−λχκ(1−e−κ(T−t))+μξ−λξγ(1−e−γ(T−t))+121−e−2κ(T−t)2κσχ2+21−e−(κ+γ)(T−t)κ+γσχσξρχξ+1−e−2γ(T−t)2γσξ2.\begin{array}{rcl}A\left(T-t)& =& -\frac{{\lambda }_{\chi }}{\kappa }(1-{e}^{-\kappa \left(T-t)})+\frac{{\mu }_{\xi }-{\lambda }_{\xi }}{\gamma }(1-{e}^{-\gamma \left(T-t)})\\ & & +\frac{1}{2}\left(\frac{1-{e}^{-2\kappa \left(T-t)}}{2\kappa }{\sigma }_{\chi }^{2}+2\frac{1-{e}^{-(\kappa +\gamma )\left(T-t)}}{\kappa +\gamma }{\sigma }_{\chi }{\sigma }_{\xi }{\rho }_{\chi \xi }+\frac{1-{e}^{-2\gamma \left(T-t)}}{2\gamma }{\sigma }_{\xi }^{2}\right).\end{array}Here, ϕ{\boldsymbol{\phi }}is a diagonal matrix that consists of autoregressive (AR) coefficients for each marginal measurement error of different contracts, and Δt\Delta tis the time difference in years between t−1t-1and tt. We may generalise AR process in measurement errors vt{{\bf{v}}}_{t}with order pp, and introduce additional ϕm{{\boldsymbol{\phi }}}_{m}matrices for m=1,…,pm=1,\ldots ,pin (4). For state and measurement errors, denoted as wt{{\bf{w}}}_{t}and vt{{\bf{v}}}_{t}, we assume that they are independent of each other, and (6)wt∼N(0,W),vt∼N(0,V),{{\bf{w}}}_{t}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{\bf{W}}),\hspace{1em}{{\bf{v}}}_{t}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{\bf{V}}),where (7)W=1−e−2κΔt2κσχ21−e−(κ+γ)Δtκ+γσχσξρχξ1−e−(κ+γ)Δtκ+γσχσξρχξ1−e−2γΔt2γσξ2,V=s112s12⋯s1Ns12s222⋯s2N⋮⋮⋱⋮s1Ns2N⋯sNN2.\begin{array}{rcl}{\bf{W}}& =& \left(\begin{array}{cc}\frac{1-{e}^{-2\kappa \Delta t}}{2\kappa }{\sigma }_{\chi }^{2}& \frac{1-{e}^{-\left(\kappa +\gamma )\Delta t}}{\kappa +\gamma }{\sigma }_{\chi }{\sigma }_{\xi }{\rho }_{\chi \xi }\\ \frac{1-{e}^{-\left(\kappa +\gamma )\Delta t}}{\kappa +\gamma }{\sigma }_{\chi }{\sigma }_{\xi }{\rho }_{\chi \xi }& \frac{1-{e}^{-2\gamma \Delta t}}{2\gamma }{\sigma }_{\xi }^{2}\\ \end{array}\right),\\ {\bf{V}}& =& \left(\begin{array}{cccc}{s}_{11}^{2}& {s}_{12}& \cdots \hspace{0.33em}& {s}_{1N}\\ {s}_{12}& {s}_{22}^{2}& \cdots \hspace{0.33em}& {s}_{2N}\\ \vdots & \vdots & \ddots & \vdots \\ {s}_{1N}& {s}_{2N}& \cdots \hspace{0.33em}& {s}_{NN}^{2}\end{array}\right).\end{array}We denote sjj2{s}_{jj}^{2}to be variances of measurement errors for contract jj, and sjk{s}_{jk}to be covariances of measurement errors between contracts jjand kk, where (j,k)=1,…,N\left(j,k)=1,\ldots ,N, and j≠kj\ne k.For estimation of covariances of measurement errors, we follow the estimation method in [18], where we estimate correlation coefficients of measurement errors between different contracts, and convert them back to covariances. We use the estimation approach introduced in [17] for correlation coefficients, which is often used in credit risk modelling. Let zj{{\bf{z}}}_{j}be the normalised prices of the contract jj, so that the vector of prices consists of z=(z1,z2,…,zN){\bf{z}}=\left({{\bf{z}}}_{1},{{\bf{z}}}_{2},\hspace{-0.18em}\ldots ,{{\bf{z}}}_{N})for NNcontracts. We assume that (8)zj=ρjε0+1−ρj2εj,j=1,…,N,{{\bf{z}}}_{j}={\rho }_{j}{\varepsilon }_{0}+\sqrt{1-{\rho }_{j}^{2}}{\varepsilon }_{j},\hspace{0.33em}j=1,\ldots ,N,where ε=(ε0,…,εN)∼N(0,IN+1){\boldsymbol{\varepsilon }}=\left({\varepsilon }_{0},\hspace{-0.18em}\ldots ,{\varepsilon }_{N})\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{{\bf{I}}}_{N+1}), ε0{\varepsilon }_{0}is a systematic component, and ε1,…,εN{\varepsilon }_{1},\ldots ,{\varepsilon }_{N}are idiosyncratic components. Therefore, (9)Corr[(zj,zk)]=ρjρk,(j,k)=1,…,N,j≠k.{\rm{Corr}}{[}({{\bf{z}}}_{j},{{\bf{z}}}_{k})]={\rho }_{j}{\rho }_{k},\hspace{0.33em}\left(j,k)=1,\ldots ,N,\hspace{0.33em}j\ne k.In this setup, we have the following correlation matrix structure: (10)R=1ρ1ρ2ρ1ρ3⋯ρ1ρNρ1ρ21ρ2ρ3⋯ρ2ρNρ1ρ3ρ2ρ31⋯ρ3ρN⋮⋮⋮⋱⋮ρ1ρNρ2ρN⋯ρN−1ρN1.{\bf{R}}=\left(\begin{array}{ccccc}1& {\rho }_{1}{\rho }_{2}& {\rho }_{1}{\rho }_{3}& \cdots \hspace{0.33em}& {\rho }_{1}{\rho }_{N}\\ {\rho }_{1}{\rho }_{2}& 1& {\rho }_{2}{\rho }_{3}& \cdots \hspace{0.33em}& {\rho }_{2}{\rho }_{N}\\ {\rho }_{1}{\rho }_{3}& {\rho }_{2}{\rho }_{3}& 1& \cdots \hspace{0.33em}& {\rho }_{3}{\rho }_{N}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {\rho }_{1}{\rho }_{N}& {\rho }_{2}{\rho }_{N}& \cdots \hspace{0.33em}& {\rho }_{N-1}{\rho }_{N}& 1\\ \end{array}\right).Using V=D12RD−12{\bf{V}}={{\bf{D}}}^{\tfrac{1}{2}}{\bf{R}}{{\bf{D}}}^{-\tfrac{1}{2}}, where D=diag(s112,…,sNN2){\bf{D}}=\hspace{0.1em}\text{diag}\hspace{0.1em}\left({s}_{11}^{2},\ldots ,{s}_{NN}^{2})is a diagonal matrix that consists of volatilities of measurement errors, the estimation will involve estimating both D{\bf{D}}and R{\bf{R}}. Hence, the modified covariance matrix V{\bf{V}}is applied in the Kalman filter and in the estimation procedure.3.2Modelling of the logarithmic returnsIn this section, we develop the measurement equations for the logarithms of relative returns on futures prices. Since the logarithmic returns are differences of the logarithmic prices at time ttand t−1t-1, we set up a linear state space model in the following way.For t=1,2,…,n−1t=1,2,\ldots ,n-1, state and measurement equations are written as follows: (11)xtr=ctr+Gtrxt−1r+wtr,{{\bf{x}}}_{t}^{r}={{\bf{c}}}_{t}^{r}+{{\bf{G}}}_{t}^{r}{{\bf{x}}}_{t-1}^{r}+{{\bf{w}}}_{t}^{r},(12)rt=yt−yt−1=dtr+Btrxtr+vtr,{{\bf{r}}}_{t}={{\bf{y}}}_{t}-{{\bf{y}}}_{t-1}={{\bf{d}}}_{t}^{r}+{{\bf{B}}}_{t}^{r}{{\bf{x}}}_{t}^{r}+{{\bf{v}}}_{t}^{r},where (13)xtr=xtxt−1,cr=c0,Gr=G0I0,wr=W000,dtr=dt−dt−1,Btr=Bt−Bt−1,vtr=vt−vt−1,vtr∼N(0,Vr).\begin{array}{rcl}{{\bf{x}}}_{t}^{r}& =& \left(\begin{array}{c}{{\bf{x}}}_{t}\\ {{\bf{x}}}_{t-1}\end{array}\right),\hspace{1em}{{\bf{c}}}^{r}=\left(\begin{array}{c}{\bf{c}}\\ {\bf{0}}\end{array}\right),\hspace{1em}{{\bf{G}}}^{r}=\left(\begin{array}{cc}{\bf{G}}& {\bf{0}}\\ {\bf{I}}& {\bf{0}}\end{array}\right),\hspace{1em}{{\bf{w}}}^{r}=\left(\begin{array}{cc}{\bf{W}}& {\bf{0}}\\ {\bf{0}}& {\bf{0}}\end{array}\right),\\ {{\bf{d}}}_{t}^{r}& =& {{\bf{d}}}_{t}-{{\bf{d}}}_{t-1},\hspace{1em}{{\bf{B}}}_{t}^{r}=\left(\begin{array}{c}{{\bf{B}}}_{t}\\ -{{\bf{B}}}_{t-1}\end{array}\right),\\ {{\bf{v}}}_{t}^{r}& =& {{\bf{v}}}_{t}-{{\bf{v}}}_{t-1},\hspace{1em}{{\bf{v}}}_{t}^{r}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,{{\bf{V}}}^{r}).\end{array}Constant vectors, transition matrices, and measurement errors are from the original model setup shown in (5) and (7), with I{\bf{I}}being the identity matrix. Note that Vr{{\bf{V}}}^{r}is the covariance matrix of the measurement errors of logarithm of returns, instead of logarithm of prices. If vtr{{\bf{v}}}_{t}^{r}follows the AR process, then we can also set our measurement errors similar to Section 3.1. We can proceed with the standard Kalman filter and maximise the likelihood function using new notations accordingly.3.2.1Parameter estimationThe unknown parameter set ψ=(κ,σχ,λχ,γ,μξ,σξ,λξ,ρχξ,V,ϕ)\psi =(\kappa ,{\sigma }_{\chi },{\lambda }_{\chi },\gamma ,{\mu }_{\xi },{\sigma }_{\xi },{\lambda }_{\xi },{\rho }_{\chi \xi },{\bf{V}},{\boldsymbol{\phi }})is estimated by optimising the log-likelihood function of y{\bf{y}}, the joint distribution of (y1,y2,…,yn)({{\bf{y}}}_{1},{{\bf{y}}}_{2},\hspace{-0.18em}\ldots ,{{\bf{y}}}_{n}), with respect to ψ\psi . The log-likelihood function is (14)ℓ(ψ;y)=∑t=1np(yt∣It−1),\ell \left(\psi ;\hspace{0.33em}{\bf{y}})=\mathop{\sum }\limits_{t=1}^{n}p({{\bf{y}}}_{t}| {{\mathfrak{I}}}_{t-1}),where p(yt∣It−1)p({{\bf{y}}}_{t}| {{\mathfrak{I}}}_{t-1})is the logarithm of the conditional probability density of yt{{\bf{y}}}_{t}given information available until time t−1t-1. Now, assuming that the prediction error et=yt−E[yt∣It−1]{{\bf{e}}}_{t}={{\bf{y}}}_{t}-{\mathbb{E}}{[}{{\bf{y}}}_{t}| {{\mathfrak{I}}}_{t-1}]follows a multivariate normal distribution, the log-likelihood function in (14) can be re-expressed as follows: (15)ℓ(ψ;y)=−Nn2ln(2π)−12∑t=1n[ln(det(Lt∣t−1))+et′Lt∣t−1−1et],\ell \left(\psi ;\hspace{0.33em}{\bf{y}})=-\frac{Nn}{2}\mathrm{ln}\left(2\pi )-\frac{1}{2}\mathop{\sum }\limits_{t=1}^{n}{[}\mathrm{ln}\left(\det \left({{\bf{L}}}_{t| t-1}))+{{\bf{e}}}_{t}^{^{\prime} }{{\bf{L}}}_{t| t-1}^{-1}{{\bf{e}}}_{t}],where Lt∣t−1=Cov(et∣It−1){{\bf{L}}}_{t| t-1}={\rm{Cov}}({{\bf{e}}}_{t}| {{\mathfrak{I}}}_{t-1}), and det(⋅)\hspace{0.1em}\text{det}\hspace{0.1em}(\cdot )denotes the determinant of a matrix. Then, parameter estimates are obtained by maximising (15) with respect to ψ\psi jointly. To obtain the quantity for the log-likelihood function, latent state vectors and their covariance matrices at each time ttneed to be estimated through the Kalman filter. [5] detected the parameter identification problem within the log-likelihood function in the Kalman filter, and hence, the constraint κ≥γ\kappa \ge \gamma is considered in the optimisation procedure.4Simulation studyIn this section, we perform a simulation study to validate the new approach described in Section 3.1. We focus on validating serial correlation and inter-correlation assumptions in measurement error considered under the Schwartz–Smith model framework, to see how our novel approach performs at estimating parameters and state variables.The steps for this simulation study are as follows. (1)Set the true parameter set ψ=(κ,σχ,λχ,γ,μξ,σξ,λξ,ρχξ,V,ϕ)\psi =(\kappa ,{\sigma }_{\chi },{\lambda }_{\chi },\gamma ,{\mu }_{\xi },{\sigma }_{\xi },{\lambda }_{\xi },{\rho }_{\chi \xi },{\bf{V}},{\boldsymbol{\phi }}). Then, obtain the state and measurement variables xt{{\bf{x}}}_{t}and yt{{\bf{y}}}_{t}by simulating error terms wt{{\bf{w}}}_{t}and vt{{\bf{v}}}_{t}. vt{{\bf{v}}}_{t}is assumed to follow an AR(1) process.(2)Choose appropriate initial values, and determine appropriate feasible bounds for each parameter.(3)Conduct the optimisation procedure through the Kalman filter, which involves jointly estimating state variables and parameter estimates. Obtain parameter estimates, and estimated values for state variables.We simulate for n=200n=200, 500, 1,000, 2,000, 5,000, with five different futures contracts maturing in 1, 2, 3, 4, and 5 months. The parameter estimates are presented in Tables 1–3, with standard errors computed using Monte Carlo approach, to assess the accuracy of parameter estimates. The MATLAB code for the simulation study is available at https://github.com/Junee1992/EUA_Futures_Pricing/tree/main/Serial-Correlation.Table 1Estimates of parameters indicated in top row. Their true values are given in the last row. nnis used for number of time-steps indicated; (standard errors are provided and were computed using 1,000 paths)nnκ\kappa σχ{\sigma }_{\chi }λχ{\lambda }_{\chi }γ\gamma μξ{\mu }_{\xi }σξ{\sigma }_{\xi }λξ{\lambda }_{\xi }ρχξ{\rho }_{\chi \xi }2002.93160.10700.07152.84281.48800.16050.0290−0.9935-0.9935(0.0245)(0.0114)(0.0022)(0.0226)(0.0118)(0.0120)(0.0020)(0.0231)5000.89830.12040.00030.00270.00158.9136×10−68.9136\times 1{0}^{-6}0.0750−0.9254-0.9254(0.0230)(0.0062)(0.0016)(0.0188)(0.0098)(0.0069)(0.0015)(0.0232)1,0001.90420.11740.00020.25420.11910.03570.01820.9999(0.0226)(0.0028)(0.0013)(0.0159)(0.0082)(0.0034)(0.0013)(0.0218)2,0002.99850.05740.00011.02450.49640.11410.00000.9519(0.0223)(0.0040)(0.0010)(0.0133)(0.0067)(0.0043)(0.0007)(0.0157)5,0002.25230.05600.00061.10600.53780.11650.00000.9065(0.0189)(0.0015)(0.0006)(0.0114)(0.0057)(0.0015)(0.0006)(0.0105)True20.10.0110.50.10.010.8Table 2Estimates of autoregressive coefficients of the measurement errors. (Standard errors are provided and were computed using 1,000 paths)nnϕ11{\phi }_{11}ϕ12{\phi }_{12}ϕ13{\phi }_{13}ϕ14{\phi }_{14}ϕ15{\phi }_{15}2000.88900.93150.91100.91680.8774(0.0100)(0.0084)(0.0070)(0.0064)(0.0061)5000.91350.88090.91440.86970.8999(0.0011)(0.0007)(0.0007)(0.0006)(0.0006)1,0000.91760.92300.90310.91270.8889(0.0006)(0.0004)(0.0004)(0.0004)(0.0004)2,0000.90800.88950.90390.88710.8937(0.0293)(0.0003)(0.0003)(0.0003)(0.0003)5,0000.90430.90190.90080.90640.8956(0.0002)(0.0002)(0.0002)(0.0002)(0.0002)True0.90.90.90.90.9Table 3Estimates of standard deviations and factors of correlation coefficients of the measurement errors. (standard errors are provided and were computed using 1,000 paths)nns1{s}_{1}s2{s}_{2}s3{s}_{3}s4{s}_{4}s5{s}_{5}ρ1{\rho }_{1}ρ2{\rho }_{2}ρ3{\rho }_{3}ρ4{\rho }_{4}ρ5{\rho }_{5}2000.01200.01280.01210.01160.01140.85180.88440.86080.85720.8664(0.0001)(0.0001)(0.0001)(0.0001)(0.0001)(0.0021)(0.0013)(0.0012)(0.0013)(0.0012)5000.01010.00980.01000.00960.00920.75690.79510.80130.78730.7693(4.67×10−54.67\times 1{0}^{-5})(4.01×10−54.01\times 1{0}^{-5})(3.69×10−53.69\times 1{0}^{-5})(3.35×10−53.35\times 1{0}^{-5})(3.00×10−53.00\times 1{0}^{-5})(0.0048)(0.0040)(0.0040)(0.0034)(0.0028)1,0000.00990.00990.01060.01050.01000.79610.80730.82120.81570.8066(2.83×10−52.83\times 1{0}^{-5})(2.41×10−52.41\times 1{0}^{-5})(2.17×10−52.17\times 1{0}^{-5})(1.93×10−51.93\times 1{0}^{-5})(1.74×10−51.74\times 1{0}^{-5})(0.0024)(0.0017)(0.0015)(0.0014)(0.0012)2,0000.00960.00960.01020.01010.01000.79810.78350.81200.79520.8142(1.34×10−41.34\times 1{0}^{-4})(2.46×10−52.46\times 1{0}^{-5})(3.28×10−53.28\times 1{0}^{-5})(3.29×10−53.29\times 1{0}^{-5})(3.67×10−53.67\times 1{0}^{-5})(0.0180)(0.0089)(0.0095)(0.0094)(0.0098)5,0000.00900.00910.00910.00940.00940.90430.90190.90080.90640.8956(9.748×10−69.748\times 1{0}^{-6})(8.258×10−68.258\times 1{0}^{-6})(7.605×10−67.605\times 1{0}^{-6})(6.934×10−66.934\times 1{0}^{-6})(6.277×10−66.277\times 1{0}^{-6})(0.0008)(0.0007)(0.0006)(0.0006)(0.0005)True0.010.010.010.010.010.80.80.80.80.8Overall, parameter estimates are quite close to their true parameters, except for κ,λχ\kappa ,{\lambda }_{\chi }and λξ{\lambda }_{\xi }. κ,γ\kappa ,\gamma and ρχξ{\rho }_{\chi \xi }tend to fluctuate for n≤1,000n\le \hspace{0.1em}\text{1,000}\hspace{0.1em}; however, it stabilises as the sample size increases. Estimated AR coefficients, volatilities, and correlation coefficients are close to the true values, indicating that the model is able to capture dependencies between measurement errors of contracts with different maturities as well as their serial correlations. The estimation errors of first eight parameters are summarised in Figure 1.Figure 1Parameter estimation error plots versus n=500n=500, 1,000, 2,000, 5,000.The simulated and estimated state variables are shown in Figure 2, along with mean absolute errors (MAE) calculated for two latent variables in each panel, showing that estimation of state variables improves as nnincreases.Figure 2Actual and estimated χt{\chi }_{t}and ξt{\xi }_{t}versus n=500n=500, 1,000, 2,000, 5,000, and MAE for two factors, starting from top-left to right-bottom.5Application to EUA futures dataGiven the increasing importance of reducing carbon emissions and the role of the EU-ETS as the largest emissions trading scheme in the world, we have applied our model to this market. We consider the historical daily prices of EUA futures contracts traded in the Intercontinental Exchange from January 30, 2017, to April 30, 2021, obtained from Refinitiv Datascope. This period covers the second half of Phase III and the beginning of Phase IV of the EU-ETS. We consider the first seven available contracts for model calibration.Figure 3 illustrates the historical daily futures price dynamics through the trading phases of the EU ETS (Phases I–IV). The graph clearly illustrates the necessity of developing sufficiently versatile models to be able to accommodate the intricacies of EUA futures price dynamics across the different phases, including Phase IV. In addition, Figure 4 shows the term structure of EUA futures prices, which matures in December annually. The main difference between commonly traded commodities and EUA is that futures curves tend to increase smoothly for as we increase the maturity of futures contracts over time.Figure 3Path of prices of first available EUA futures contract for the period from April 22, 2005, to April 30, 2021.Figure 4Term structure of EUA futures prices for seven available contracts from January 30, 2017, to April 30, 2021.In the two works by [23] and [5], measurement errors are assumed to be independent. However, based on the belief that the price movement for each contract must be correlated with other available contracts during the same period within the same commodity market, we assume that measurement errors have the full covariance matrix. After investigating the statistics of measurement errors for EUA future prices, we found that measurement errors of contracts with different maturities are highly correlated with each other, as close to 1. The correlation matrix of measurement errors of EUA future prices using Model 1 is as follows: R=10.99990.99990.99990.99950.99990.999510.99990.99980.99940.99990.999510.99970.99930.99990.999610.99920.99980.999410.99960.998210.99941.{\boldsymbol{R}}=\left(\begin{array}{ccccccc}1& 0.9999& 0.9999& 0.9999& 0.9995& 0.9999& 0.9995\\ & 1& 0.9999& 0.9998& 0.9994& 0.9999& 0.9995\\ & & 1& 0.9997& 0.9993& 0.9999& 0.9996\\ & & & 1& 0.9992& 0.9998& 0.9994\\ & & & & 1& 0.9996& 0.9982\\ & & & & & 1& 0.9994\\ & & & & & & 1\end{array}\right).We also observe that each series of measurement errors follow AR(1) process, with all AR coefficients being highly significant (with pp-value <0.0001\lt 0.0001). Hence, the price data of EUA futures contracts are suitable to test the model that is developed in this study.For comparative analysis, we use two different models. The following two models consider estimating inter-correlation between measurement errors in different settings: Model 1– Modelling of logarithmic returns;Model 2– Modelling of logarithmic returns with serially correlated measurement errors.For goodness-of-fit assessment, we present the performance of each model using root mean-squared error (RMSE). The results are summarised in Table 4 for futures contracts Cj{C}_{j}, j=1,2,…,7j=1,2,\ldots ,7. In both models, the logarithmic returns are converted to the logarithm of prices. In temrs of RMSE solely, Model 1 performs better than Model 2. However, using model selection criteria (Akaike and Bayesian information criterion), Model 2 is preferable to Model 1.Table 4Akaike (AIC), Bayesian information criterion (BIC), and RMSE for Models 1 and 2Model 1Model 2C1{C}_{1}0.02380.0281C2{C}_{2}0.02300.0264C3{C}_{3}0.02360.0262C4{C}_{4}0.02690.0288C5{C}_{5}0.04930.0507C6{C}_{6}0.06950.0703C7{C}_{7}0.08960.0900AIC−8.9906×104-8.9906\times 1{0}^{4}−8.9985×104-8.9985\times 1{0}^{4}BIC−8.9796×104-8.9796\times 1{0}^{4}−8.9840×104-8.9840\times 1{0}^{4}Next, we provide results of out-of-sample predictions using 30/50-day out-of-sample windows in Tables 5 and 6. In each setting, we considered four different scenarios, where we repeat parameter estimation and out-of-sample prediction for every 1 day, 5 days, 10 days, and 30 or 50 days. We use the deseasonalised price data from January 30, 2017, to February 18, 2021, n=1,040n=\hspace{0.1em}\text{1,040}\hspace{0.1em}business days.Table 5RMSE for 1, 5, and 10 days ahead predictions in 30-day out-of-sample window using 30/1/2017–18/2/2021 dataModel 1Model 21 day5 days10 days30 days1 day5 days10 days30 daysC1{C}_{1}0.04470.04480.04500.04740.04400.04310.04480.0467C2{C}_{2}0.04490.04510.04530.04840.04390.04350.04520.0472C3{C}_{3}0.04480.04500.04520.04890.04390.04360.04520.0473C4{C}_{4}0.04430.04550.04470.04880.04380.04340.04500.0471C5{C}_{5}0.04490.04510.04530.05060.04440.04400.04550.0476C6{C}_{6}0.04540.04570.04590.05060.04500.04470.04620.0483C7{C}_{7}0.04610.04640.04660.05150.04580.04550.04690.0490Table 6RMSE for 1, 5, and 10 days ahead predictions in 50-day out-of-sample window using 30/1/2017–18/2/2021 dataModel 1Model 21 day5 days10 days50 days1 day5 days10 days50 daysC1{C}_{1}0.06440.07110.13520.07330.07720.06510.15690.0732C2{C}_{2}0.06700.07210.13730.07490.07970.06680.15830.0741C3{C}_{3}0.06860.07320.13940.07650.08120.06860.15990.0751C4{C}_{4}0.06890.07310.14020.07670.08150.06920.16030.0749C5{C}_{5}0.07030.07420.14200.07810.08260.07070.16150.0759C6{C}_{6}0.07110.07470.14300.07880.08330.07170.16200.0764C7{C}_{7}0.07200.07540.14410.07950.08390.07260.16260.0769From Table 5, we found that, in general, Model 2 performed better than Model 1, although differences in RMSE were minimal. We used Diebold–Mariano test for detecting a significant difference in out-of-sample predictability between the two models. For 30-day out-of-sample predictions, by using all possible pairs of series of marginal predictions under the same hh-days-ahead predictions for h=1,5,10,30h=1,5,10,30, we had pp-values ≈0.31\approx 0.31. Hence, there is no significant difference in the two models’ out-of-sample predictions.In Table 6, we found that Model 1 performed better for out-of-sample predictions for h=1,10h=1,10day, whereas Model 2 obtained the lower RMSE for h=5,50h=5,50. The Diebold–Mariano test showed no significant difference in 50 days out-of-sample predictions between the two models again. This pattern persisted in other subsets of data from Phases III and IV.6ConclusionIn this article, we have presented the modified Schwartz–Smith two-factor model, which can be used for modelling of logarithmic returns of futures contracts, incorporating serial correlation and inter-correlation of measurement errors. We compared the two different models that use the logarithm of futures prices and the logarithmic returns, applying them to EUA futures price data from Phase III and Phase IV of EU-ETS.The simulation study has illustrated that our novel approach was able to jointly estimate both parameters and state variables when serial correlations and inter-correlations were present in measurement errors. We have illustrated that the parameter and state variables estimates converged as the sample size increases. The maximum likelihood method was used for the estimation of coefficients of Gaussian AR processes used for modelling of serial correlations in measurement errors.Finally, we discussed the results of the comparative analysis of the two models in the context of EUA futures data. The goodness-of-fit and out-of-sample performances in predicting futures prices were discussed for each model. Overall, models for logarithmic returns showed a good performance in terms of RMSE for out-of-sample predictions. The full Model 2, for logarithmic returns with serial correlations, performed better than its reduced-form Model 1 for calibration of data in terms of AIC and BIC, and for predicting for 30 days out-of-sample window.The empirical results of the study emphasised the necessity of considering both serial correlations and inter-correlations of measurement errors for modelling of futures prices. For parameter estimation in Model 2, we used the two-step approach. First, we obtained the parameter estimates similar to Model 1, then we proceeded with estimating parameters of AR processes used for modelling of serial correlations. By using the real data through the incorporation of serial correlations in the measurement errors, we showed that the proposed models for logarithmic return capture the price movement more reliably.

Journal

Dependence Modelingde Gruyter

Published: Jan 1, 2022

Keywords: pricing; futures; commodity; CO 2 emission allowances; Kalman filter; correlation; maximum likelihood estimation; linear state-space model; 62P05; 91B70

There are no references for this article.