Pricing Energy Contracts under Regime Switching Time-Changed models

A main goal of the paper is to develop a methodology to price European option’s contracts on electricity and future oil prices. The approach is based on Fourier expansions and implements models that capture specific stylized features of the underlying assets such as stochastic volatility and random jumps. In particular, we consider a switching time-changed Levy process as an alternative to the BSM model and implement a pricing algorithm based on the expansion of its characteristic function as considered in [6].

We compare the prices we obtain with those obtained via a computationally costly, but accurate, Monte Carlo method and study sensitivities with respect to relevant parameters, e.g. maturity, strike price and initial price. Moreover, we contrast prices under the regime switching model to those given by the Black-Scholes equation and show that the prices agree when the switching model is reduced to the Black-Scholes model.

We rely on the Esscher transformation, see [8], to obtain an equivalent martingale measure(EMM) and in order to work on a risk-neutral setting. To calibrate parameters, we use market option prices and minimize the mean squared error. On the other hand, and in order to estimate parameters, we implement the method of moments, minimum distance estimation and maximum likelihood estimation techniques. For simplicity an Expectation Maximization algorithm (EM) is not considered.

Although most of these elements have been previously implemented, to the best of our knowledge, the combination consisting of the selected model class, the pricing methodology, the risk-neutral framework and the estimation/calibration approach have not been studied before in energy markets or elsewhere. It is worth noticing that non-switching models with Levy noises have been introduced in [3]. On the other hand, results for switching Levy models and their characteristic functions can be found in [4].

Many financial time series, including commodity futures seem to exhibit dramatic breaks in their behaviour, for example in the events of political changes or financial crises. Different intervals sharing similar characteristic can be grouped together under a single regime. Models that can capture such behaviour are regime switching Levy. Under such a model, the process switches randomly between different Levy processes according to an unobservable Markov chain. The regime switching time-changed Levy process is a pure jump process which captures two key features of the market: the existence of regimes and price jumps.

The organization of the paper is the following: Section 2 introduces the regime switching time-changed Levy model, we then derive its characteristic function under Gamma and Inverse Gaussian subordinators. In section 3 we discuss Monte Carlo and Fourier Cosine pricing methods. For Monte Carlo, we develop an algorithm to simulate trajectories of the process, as well as for pricing European call options by simulating many regime switching processes simultaneously. In section 4 we use calibration and various methods to estimate values of model parameters from option quotes and historical prices of oil and electricity commodities.

Let $(\Omega,\mathcal{A},(\mathcal{F}_{t})_{t\geq 0},P)$ be a filtered probability space verifying the usual conditions. For a stochastic process $(X_{t})_{t\geq 0}$ defined on the space filtered space above the functions $\varphi_{X_{t}}(u)$ and $\Psi_{X}(u)=\frac{1}{t}\log\varphi_{X_{t}}(u)$ define its characteristic function and characteristic exponent respectively. When the process has stationary and independent increments the later does not depend on $t$. $A^{T}$ is the transpose of matrix $A$ unless it is specified differently.
Let the process $\{S_{t}\}_{t\geq 0}$ represent the price of the underlying asset at any time $t>0$, $X_{t}=\log S_{t}$ represents the logarithm of the prices.
We define a continuous-time Markov chain $\{s_{t}\}_{0\leq t\leq T}$ driving the changes between regimes with the state space $E=\{1,2\}$. The switching times are described by a sequence of independent random variables $(\tau_{j})_{j\in\mathbb{N}}$ in a way that:

$$\lim_{t\rightarrow\tau^{-}_{k}}s_{t}\neq s_{\tau_{k}}\>\>\>\>\text{for}\>\>\>\>k=1,2.$$

The infinitesimal generator matrix of the continuous-time Markov chain is given by

(1)

$$Q=\begin{bmatrix}-\lambda_{12}&\lambda_{12}\\ \lambda_{21}&-\lambda_{21}\end{bmatrix}.$$

Hence:

	$\displaystyle P\{s_{t+h}=2/s_{t}=1\}$	$\displaystyle=$	$\displaystyle\lambda_{12}h+o(h)$
	$\displaystyle P\{s_{t+h}=1/s_{t}=2\}$	$\displaystyle=$	$\displaystyle\lambda_{21}h+o(h)$

Next, consider a collection of independent subordinators $L^{j}=\{L^{j}_{t}\}_{t\geq 0}$ for $j\in E$, where each subordinator $L^{j}$ is also independent of each process $X^{i}$, for $i,j\in E$. Each subordinator is characterized by two parameters $\alpha_{j},\beta_{j}>0$ which change between states on the (non-observable) Markov Chain. Each process $L^{j}$ is a pure jump process and each process $X^{j}$ has almost surely continuous paths
We define the collection of time-changed Levy processes $Y^{j}=\{Y_{t}^{j}\}_{t\geq 0}$ where

$$Y_{t}^{j}=\mu_{j}L_{t}^{j}+\sigma_{j}B_{L_{t}^{j}}.$$

with $\mu_{j}\in\mathbb{R}$ and $\sigma_{j}>0$.
There exists a natural economic interpretation of a time-changed process. Energy markets alternate at random times between calmed and frenzy periods at random times.
We now define the regime switching time-changed Levy process $Z=\{Z_{t}\}_{t\geq 0}$ as:

(2)

$$Z_{t}\equiv Y_{t}^{s_{t}}\>\>\>\>\text{where}\>\>\>\>Y^{s_{t}}_{t}=\mu_{s_{t}}L^{s_{t}}_{t}+\sigma_{s_{t}}B_{L^{s_{t}}_{t}}$$

or, in differential terms:

$$dZ_{t}\equiv dY_{t}^{s_{t}}\>\>\>\>\text{where}\>\>\>\>dY^{s_{t}}_{t}=\mu_{s_{t}}dL^{s_{t}}_{t}+\sigma_{s_{t}}dB_{L^{s_{t}}_{t}}$$

The regime switching time-changed Levy process $Z$ is assumed to be the log-price process of the underlying asset and the stochastic process of the asset price itself $\{S_{t}\}_{t\geq 0}$ is defined as:

$$S_{t}=S_{0}\exp(Z_{t}).$$

For simplicity we assume that the process will always start out at state 1 with probability 1.

Following [4], for the process $Z$ defined in equation (2) along with a Markov chain with the infinitesimal generator matrix $Q$ defined in equation (1), its characteristic function is given by:

(3)

$$\varphi_{Z}(u)=\exp(iuy_{0})\leavevmode\nobreak\ \mathbb{E}_{\mathcal{Q}}\{[1,1]\exp(t\Phi(u))[1,0]^{T}\}$$

where $y_{0}=\log S_{0}$ and $\Phi(u)$ is the matrix:

$$\Phi(u)=\left(\begin{array}[]{cc}-\lambda_{12}+\Psi_{L_{t}^{(1)}}(\mu^{1}u+\frac{1}{2}i\sigma^{2}_{1}u^{2})&\lambda_{12}\\ \lambda_{21}&-\lambda_{21}+\Psi_{L_{t}^{(2)}}(\mu_{2}u+\frac{1}{2}i\sigma^{2}_{2}u^{2})\\ \end{array}\right).$$

Notice that conditionally on $s_{t}=j$ the characteristic function of $Z^{j}$, i.e. the characteristic function of $Z$ conditionally on $s_{t}=j$, is:

$$\varphi_{Z_{t}}(u)=\varphi_{L^{j}_{t}}(\mu_{j}u+\frac{1}{2}i\sigma^{2}_{j}u^{2}).$$

To compute the exponential matrix $\Phi(u)$ we use a scaling and squaring algorithm, see [10], based on the following approximation:

$$e^{\Phi(u)}=(e^{2^{-s}\Phi(u)})^{2^{s}}\approx r_{m}(2^{-s}\Phi(u))^{2^{s}},$$

where $r_{m}(x)$ is the $[m/m]$ Pade approximant of $e^{x}$ and the nonnegative integers $m$ and $s$ are chosen in such a way as to achieve minimum error at minimal cost. A table of errors as a function of $s$ and $m$ is given in [1]. The $[k/m]$ Pade approximant for the exponential function is:

$$r_{km(x)}=p_{km}(x)/q_{km}(x)$$

where:

$$p_{km}(x)=\sum_{j=0}^{k}\frac{(k+m-j)!k!}{(k+m)!(k-j)!}\frac{x^{j}}{j!},\>\>\>q_{km}(x)=\sum_{j=0}^{m}\frac{(k+m-j)!m!}{(k+m)!(m-j)!}\frac{(-x)^{j}}{j!}.$$

The discounted price process is denoted $\tilde{S}=(\tilde{S}_{t})_{t\geq 0}$ where $\tilde{S}_{t}:=exp(-rt)S_{t}$. We have that under an EMM $\mathcal{Q}$, the discounted price process $\tilde{S}$ is a martingale under $\mathcal{Q}$ if and only if the following equation is satisfied:

(4)

$$\Psi_{Z}(-\text{i})=r.$$

See [14] for details.

We turn back to pricing, consider a European call contract with maturity at $T$ and strike price $K$. Its payoff, written in terms of the log-returns, is:

$$h(Z_{T})=(S_{0}e^{Z_{T}}-K)_{+}=K(exp(x_{0}+Z_{T})-1)_{+}$$

where $x_{0}:=\log(S_{0}/K)$.
The price of the contract at a time $t<T$ and $x=\log(\frac{S_{t}}{K})$ is denoted as $v(x,t)$ and verifies:

(7)

$$v(x,t)=e^{-r\Delta t}\leavevmode\nobreak\ \mathbb{E}_{\theta}[v(y,T)]=e^{-r\Delta t}\int_{\mathbb{R}}v(y,T)f_{Z_{T}}(y|x)dy,$$

Notice that $v(y,T)$ is the payoff at maturity time $T$ and $y:=\log(S_{T}/K)$.
The value $\Delta t=T-t$ is the time to maturity and $r$ is the risk-neutral interest rate. The function $f_{Z_{T}}(y|x)$ is the probability density function (p.d.f.) of $Z_{T}$ given the value $x=\log(S_{0}/K)$.
$\mathbb{E}_{\theta}$ is the expectation value operator with respect to an EMM $\mathcal{Q}^{\theta},\theta\in\mathbb{R}$ which is determined by an Esscher transform of the historic measure $P$. See [8] for a rationale in terms of a utility-maximization criteria.
For a stochastic process $(X_{t})_{t\geq 0}$ the latter is defined as:

(8)

$$\frac{d\mathcal{Q}^{\theta}_{t}}{dP_{t}}=\exp(\theta X_{t}-t\leavevmode\nobreak\ l_{X}(\theta)),\;0\leq t\leq T,\;\theta\in\mathbb{R}$$

where $P_{t}$ and $\mathcal{Q}^{\theta}_{t}$ are the respective restrictions of $P$ and $\mathcal{Q}^{\theta}$ to the $\sigma$-algebra $\mathcal{F}_{t}$. We define by $\varphi^{\theta}_{X_{t}}$, $\Psi^{\theta}_{X_{t}}$ and $l^{\theta}_{X}(u)$ respectively the characteristic function, characteristic exponent and moment generating function of a process $(X_{t})_{t\geq 0}$ under the probability $\mathcal{Q}^{\theta}$ obtained by an Esscher transformation as given in equation (8).
We follow a pricing approach via Fourier- Cosine Series expansion of the p.d.f. $f_{Z_{T}}$. The method has been proposed in [6].
The solution to (7) is obtained by expanding the p.d.f. $f_{Z_{T}}(./x)$ in terms of a Fourier basis under Gamma and Inverse Gaussian subordinators introduced in the previous section .
The domain of integration is truncated to a finite interval $[a,b]$ for the purposes of numerical integration, for a discussion about selecting the truncation interval and its associated error, see [11].
The Fourier-cosine expansion of $f_{Z_{T}}$ is given by:

(9)

$$f_{Z_{T}}(y|x)=\sum_{k=0}^{\infty}A_{k}(x)\cos\Big{(}k\pi\frac{y-a}{b-a}\Big{)},$$

where the first term of the summation is weighted by one-half.
The coefficients of the Fourier expansion, denoted by $A_{k}(x)$ , are approximated by:

$$A_{k}(x)=\frac{2}{b-a}\text{Re}\Big{\{}\varphi_{Z_{T}}\Big{(}\frac{k\pi}{b-a}/x\Big{)}\exp\Big{(}-\text{i}k\pi\frac{a}{b-a}\Big{)}\Big{\}}.$$

Hence, substituting (9) into equation (7) we have:

$$v(x,t)=\frac{1}{2}(b-a)e^{-r\Delta t}\sum_{k=0}^{\infty}A_{k}(x)V_{k},$$

where $V_{k}$ is the Fourier coefficients of $v(y,T)$ given by:

$$V_{k}:=\frac{2}{b-a}\int_{a}^{b}v(y,T)\cos\Big{(}k\pi\frac{y-a}{b-a}\Big{)}dy.$$

In particular, for a European call option we have:

$$V^{Call}_{k}=\frac{2K}{b-a}\int_{0}^{b}(e^{y}-1)\cos\Big{(}k\pi\frac{y-a}{b-a}\Big{)}dy.$$

For a European put option denoted by $V^{Put}_{k}$ a similar expression is found. Finally, truncating the infinite series to $N$ terms, we obtain:

$$v(x,t)\approx e^{-r\Delta t}\sum_{0\leq k<N}\text{Re}\Big{\{}\varphi_{Z_{T}}\Big{(}\frac{k\pi}{b-a}/x\Big{)}e^{-\text{i}k\pi\frac{a}{b-a}}\Big{\}}V_{k}.$$

The characteristic function of the log-return at maturity $Z_{T}$ is computed above by equation (3) under Inverse Gaussian and Gamma subordinators. As Fourier-cosine series of entire functions converges exponentially, so $N$ does not be too large to obtain good approximations. For European call options, it is found that the price is not accurate and extremely sensitive to the values of $b.$ On the other hand, it is also found that for large values of $b$, $V^{Call}_{k}$ diverges while $V^{Put}_{k}$ converges quickly and varies little with changes in the left-end of the truncation interval $a$. We therefore rely on the put-call parity which allows for the computation of the European call option using the put option.

We summarize the calculations as follows:
Algorithm

1. (1)

   Initialization: choose appropriate boundary points $a,b$, number of terms $N$ in the series expansion and contract specifications (i.e. interest rate $r$, initial stock price $S_{0}$, strike price $K$, with $x:=S_{0}/K$ and time to maturity $\Delta t$).

2. (2)

   Initialize $N\times 1$ array of payoffs $v^{Put}$ and $v^{Call}$.

3. (3)

   For $k=0$ to $N-1$

   1. (a)

      Define the $k$-th element of $v^{Put}$ to be:
      $v^{Put}(k)=e^{-rT}\text{Re}\Big{\{}\varphi_{Z_{\Delta t}}(k\pi/(b-a);x)\exp(-\text{i}k\pi(a/(b-a))\Big{\}}V^{Put}_{k}$.
      

   2. (b)

      $v^{Call}(k)=v^{Put}(k)+S_{0}-Ke^{-rT}$   (put-call parity)

4. (4)

   $v^{Call}_{final}=\frac{1}{2}v(1)+\sum_{k=1}^{N-1}v(k)$  (Summation)

The algorithm is implemented in a desktop computer using MATLAB. We compare the European call payoff and running time under Monte Carlo simulation and Fourier-Cosine pricing in Table 1 for different maturities and strike prices.

In the parametric set considered (the next section discusses how to estimate the parameters), the Fourier-cosine method is on average 100 times faster than a standard Monte Carlo approach, while producing similar level of precision.
On the other hand, it is found that the difference between pricing European call options using Monte Carlo and Fourier-Cosine pricing remains constant for different strike prices. The error however grows linearly for increasing maturity times.
To compare with the price obtained via Monte Carlo we discuss the simulation procedure employed.
It is well-known that the Markov chain $\{s_{t}\}_{t\geq 0}$ spends independent and exponentially distributed random times between regime transitions. The $k-$th switching time is determined by:

$$\tau_{k}=\sum_{j=1}^{k}\Delta\tau_{j}$$

where $\Delta\tau_{j}$ are the times between regime changes. The parameters of the exponential random variables alternate between $\lambda_{12}$ and $\lambda_{21}$.
The differential equation is approximated numerically through finite differences using the Euler-Maruyama Method. Hence, if at time $t$ the process $Z$ is under regime $j\in E$, the increment of the process $Z$ during the time interval $[t,t+\Delta t)$ is given by:

$$\Delta Z^{j}_{t}=\mu^{j}\Delta L^{j}_{t}+\sigma_{j}\sqrt{\Delta L^{j}_{t}}N(0,1).$$

Notice that, for small $\Delta t$, the process remains under regime $j$ with a probability close to one.

Figure 2(a) shows a single realization of a switching time-changed Levy process (top) under an IG subordinator, as well as the underlying Markov chain (bottom). Figure 2(b) shows trajectories under a Gamma subordinator.
We devise an algorithm which computes the payoff of a European Call option by simulating $m$ independent realizations of the process $Z$ simultaneously and then estimating the price according to:

$\displaystyle v(x,T)$

$\displaystyle\simeq\exp(-rT)\frac{1}{m}\sum_{j=1}^{m}(S_{0}\exp(Z_{j,T})-K)_{+}$

where $Z_{j,T}$ is the $j$-th simulated log-price at maturity.

We can also estimate the price using confidence intervals. The confidence interval is useful because it provides a range of values that are likely to contain the population mean. Endpoints of the confidence interval are given by:

$$\bar{h}\pm z_{0.95}\frac{s}{\sqrt{m}},$$

where $\bar{h}$ is the sample mean of the simulated payoffs, $s$ is the sample standard deviation, $m$ is the sample size and $z_{0.95}$ is the $95\%$ -percentile of the normal distribution.
The confidence interval decreases as the number of simulations approaches infinity, however in the case of the Inverse Gaussian subordinator, the interval is larger at each simulation because Inverse Gaussian random variables have a larger variance than Gamma random variables when $\beta<1$.
At $m=10^{4},$ the confidence interval when the subordinator is Inverse Gaussian is $[18.7353,19.1168]$. When the subordinator is a Gamma process, the confidence interval is $[17.8768,17.9136]$.

Figure 4 indicates the behaviour of the price of a European call option for a regime switching time-changed Levy process model under Inverse Gaussian (top) and Gamma subordinators (bottom). Both payoff models are monotone in $T$ and $K$. Moreover as $T$ increases, the expected payoff increases. For $K>>S_{0}$ the probability that $S_{T}\geq K$ is very small, therefore the payoff is close to $0.$

In each figure parameters between states are held constant, except the ones indicated in the graph: $\mu_{1}=\mu_{2}=0.01,\>\>\alpha_{1}=\alpha_{2}=\beta_{1}=\beta_{2}=0.1,\>\>\sigma_{1}=\sigma_{2}=0.01,\>\lambda_{12}=\lambda_{21}=0.25,\>r=0.04,\>T=1,\>S_{0}=20\>$ and $K=1.$ Setting the parameters $\lambda_{12}=\lambda_{21}$ implies the process spends an equal amount of time in each regime, on average. The only exception made is in Figure 5, where $\lambda_{12}=1000$ and $\lambda_{21}=1/10$ so that the process would spend a majority of the time in the second regime; this makes the process an approximation to a time-changed Levy process.
For changes in $\beta_{1},\beta_{2},$ the payoff approaches an asymptote because Gamma and Inverse Gaussian random variables both have mean $\alpha/\beta,$ which approaches infinity for $\beta\rightarrow 0^{+}.$ In Figure 6 changes in the price with respect to the intensity parameters and the parameters of the underlynig subordinators are shown.

Notice that we can reduce the regime switching Levy processes to the Black-Scholes model by defining the subordinator $L_{t}=t$ and setting the parameters equal across both regimes. Setting the drift such that the discounted Black Scholes price process is a martingale $\mu_{1}=r-\frac{1}{2}\sigma^{2}_{1}$. we find that the price by a Fourier-cosine approach is consistent with the price given by the Black-Scholes formula. See Table 2.

We implement two approaches of fitting the parameters of the underlying model to financial historical data: calibration and estimation depending when option’s prices or the underlying electricity and oil prices are considered. In a calibration approach, the parameters are fitting by minimizing the quadratic error between the prices obtained numerically and option quotes. The option quotes are taken from Bloomberg’s data base for a variety of strike prices and times to maturity. The parameters are fitted using European call option quotes of West Texas Intermediate (WTI) crude oil.
In a parameter estimation, we implement the following techniques based on historical prices of oil and electricity: method of moments, minimum distance method and maximum likelihood estimation, combined with empirical estimation of the switching parameters. Specifically, we use daily historical NYMEX WTI crude oil futures (11-16-2012 to 06-05-2018) and IESO Ontario (Canada) Zone 24H electricity average spot prices (06-06-2008 to 06-05-2018). Spikes and stochastic volatility are observed in the series, Similar phenomena have been reported in other electricity markets, see [3, 7]. We assume that there are 250 trading days in a year with each trading day corresponding to $\Delta t=1/250$ of unit time.
Figures 7 and 8 plot the historic WTI oil futures and average Ontario electricity prices, as well as their log-returns. Electricity spot prices sometimes move below zero, implying a surplus of electricity produced during low demand. Because electricity produced by power suppliers must be consumed immediately, the supplier pays wholesale customers to buy the surplus energy, see [15]. All negative prices are arbitrarily modified to CAD $\$0.01$ for estimation purposes.
We also compare the empirical density function of the log-returns of each commodity to the normal distribution under the same mean and variance parameters as the historical log-return process. To this end we implement a kernel smooth technique. Hence, the estimated p.d.f. is:

$$\hat{f}(x;\theta)=\frac{1}{nh}\sum_{j=1}^{n}K\Big{(}\frac{x-z_{j}}{h}\Big{)},$$

where the function $K$ is the Gaussian kernel.
As the p.d.f. $f(x_{k};\theta)$ of the log-returns is not available in a close form we simulate the model under a parametric set $\theta$ to get the empirical p.d.f. $\hat{f}(x_{k};\theta)$ using a kernel smoothing technique, see for example [9], and continue the optimization procedure.
The value of $h,$ is chosen to equal Silverman’s quantity $h=1.06\sigma n^{-1/5}$, where $\sigma$ is the standard deviation of the log-return series, see [13]. See Figure 9 for empirical p.d.f.’s corresponding to WTI futures (left) and Ontario electricity log-return prices(right).

The kurtosis of WTI future and Ontario electricity log-return series are respectively 6.34 and 23.947, much larger than that of the normal distribution, suggesting the presence of heavier tails and extreme behavior. Significant negative skewness is also reported on both series, respectively -0.1314 and -0.1377
In our model the parameters are described by vectors $\theta^{j}=(\mu_{j},\sigma_{j},\alpha_{j},\beta_{j})$ for $j=1,2,$ while $\lambda_{12}$ and $\lambda_{21}$ reflect the parameters of the hidden Markov chain. Hence the times the chain remains in regime $j$ are independent and exponentially distributed with mean $\lambda^{j}=\frac{1}{\lambda_{jk}}$, with $k=\mod(2)+1$. We set $\Theta^{j}\subset\mathbb{R}^{4}$ to be the set of all feasible parameters for $\theta^{j}$. We assume that the two sets of parameters belong in different parameter spaces i.e. $\Theta^{1}\neq\Theta^{2}.$ The two parameters of the subordinator and the diffusion coefficient are required to be positive, therefore we add the natural constraints $\sigma_{j}>0,\alpha_{j}>0,\beta_{j}>0.$
We assume that the duration of the $j$-th observed historic regime is the most probable value i.e. it is equal to the expectation value $\lambda^{j}.$ For each commodity, the $j$-th holding-rate parameter is given by:

$$\lambda^{j}=\frac{\text{total number of days in regime}\>\>j}{\text{ number of occurrences of regime}\>\>j}$$

where it is assumed that there are 250 trading days every year.
By simple inspection of the log-return process data of oil futures we set the process to be in regime one between 11-16-2012 and 11-16-2014 as well as between 02-06-2017 and 06-05-2018; otherwise, we assume that the process is in regime two. In the case of the log-returns of electricity spot prices; we set the process be in regime two whenever the absolute value of the log-returns exceeds 3 and in regime one otherwise.
Table 3 shows the average time and daily standard deviation in the two regimes of the WTI futures and Ontario electricity series. We also include the variance of the log-returns within each regime; the different orders of magnitude between regimes justifies the use of a switching model.

By having defined the location of the regime changes and therefore estimated the values of $\lambda^{1},\lambda^{2}$, the historic log-returns are separated into two sets of data, one containing all the data points for each regime.
To calibrate the parameters within each regime we minimize the mean square error between the numeric option payoffs and European call option quotes. When the option is out of the money the option price is obtained using a Monte Carlo approach because the Fourier Cosine method exhibits significant error.
Thus, the objective function in regime $j$ is

$\displaystyle J(\theta^{j})$

$\displaystyle=\sqrt{\frac{1}{n}\sum_{T,K}(V(T,K;\theta^{j})-\hat{V}(T,K;\theta^{j}))^{2}},\quad j\in\{1,2\}.$

where, by a convenient change in notation to emphasize the dependence on the parameters we write $V(T,K;\theta^{j})$ and the option quotes $\hat{V}(T,K;\theta^{j})$, taken over a range of strike prices $K$ and maturity times $T$.
The optimal parameter is $\hat{\theta}^{j}=\operatorname*{arg\,min}_{\theta\in\Theta^{j}}J(\theta^{j})$. It is calculated using a gradient descent method.
The stopping criteria is taken to be step tolerance, taken to be equal to 1e-10. The k-th step tolerance is a lower bound on the size of the step $||\theta^{t}-\theta^{t-1}||_{2}$. The solver stops if the stopping criteria is reached, or if the maximum number of iterations (fixed to 1000 steps) is exceeded. Different initial starting points are found to give similar estimation of the parameters. Table 4 gives the estimated calibration in the case when the subordinator is a Gamma process or an Inverse Gaussian process. As expected, the volatility is higher in regime two in the case of both subordinators.

To choose the initial set of parameters we use the Method of Moments. Theoretical moments are computed from the characteristic function of the model under both subordinators considered. Matching both, empirical and theoretical moments up to order forth leads to the following non-linear system of equations, in the case of a model under an Inverse Gaussian subordinator:

	$\displaystyle{\hat{\mu}}_{1}$	$\displaystyle={\alpha}{\mu}{\Delta t}/{\beta}$
	$\displaystyle\hat{\mu}_{2}$	$\displaystyle={\alpha}{\Delta t}({\mu}^{2}+{\beta}^{2}{\sigma}^{2}+{\alpha}{\beta}{\mu}^{2}{\Delta t})/{\beta}^{3}$
	$\displaystyle\hat{\mu}_{3}$	$\displaystyle={\alpha}{\mu}{\Delta t}({3}{\mu}^{2}+{3}{\beta}^{2}{\sigma}^{2}+{3}{\alpha}{\beta}{\mu}^{2}{\Delta t}+{\alpha}^{2}{\beta}^{2}{\mu}^{2}{\Delta t}^{2}+{3}{\alpha}{\beta}^{3}{\sigma}^{2}{\Delta t})/{\beta}^{5}$
	$\displaystyle\hat{\mu}_{4}$	$\displaystyle=({\alpha}{\Delta t}({15}{\mu}^{4}+{3}{\beta}^{4}{\sigma}^{4}+{18}{\beta}^{2}{\mu}^{2}{\sigma}^{2}+{15}{\alpha}{\beta}{\mu}^{4}{\Delta t}+$
		$\displaystyle\quad\quad{6}{\alpha}^{2}{\beta}^{2}{\mu}^{4}{\Delta t}^{2}+{\alpha}^{3}{\beta}^{3}{\mu}^{4}{\Delta t}^{3}+{3}{\alpha}{\beta}^{5}{\sigma}^{4}{\Delta t}+{6}{\alpha}^{2}{\beta}^{4}{\mu}^{2}{\sigma}^{2}{\Delta t}^{2}+{18}{\alpha}{\beta}^{3}{\mu}^{2}{\sigma}^{2}{\Delta t}))/{\beta}^{7}.$

where ${\hat{\mu}}_{k}$ is the empirical k-th moment.
The system of equations is solved separately for each regime using the function solve in MATLAB based on the trust region algorithm. The results are summarized in Table \retab:method_moment IG.
A similar result is obtained for the model under a Gamma subordinator.

The method encounter difficulties to find a global minimum in the case where the empirical moments were calculated using electricity log-return prices. Changing the initial starting points resulted in varying results, which indicates the presence of local minima. Despite these shortcomings the method of moments can be used as an initial solution for a Minimum Distance approach based on the distance between the theoretical and empirical characteristics function of the log-returns, the later defined as:

$$\hat{\varphi}_{Z^{j}_{\Delta t}}(u)=\frac{1}{n}\sum_{k=1}^{n}\exp(\text{i}uz_{k})$$

for a sample $z_{1},z_{2},\ldots,z_{n}$ of $n$ log-returns of the underlying series. See [17] for details.
The objective function under regime $j$ is defined by:

$$||\varphi_{Z^{j}_{\Delta t}}(u;\theta)-\hat{\varphi}_{Z^{j}_{\Delta t}}(u)||_{2}:=\Bigg{(}\int_{-\infty}^{\infty}|\varphi_{Z^{j}_{\Delta t}}(u;\theta)-\hat{\varphi}_{Z^{j}_{\Delta t}}(u)|^{2}w(x)dx\Bigg{)}^{1/2},$$

where $w$ is the weight function $w(x)=(1/\sqrt{2\pi})\exp(-x^{2}/2)$.
Then, $\hat{\theta}$ is the minimum distance estimate of $\theta$ if

$$||\varphi_{Z^{j}_{\Delta t}}(u;\hat{\theta})-\hat{\varphi}_{Z^{j}_{\Delta t}}(u)||_{2}=\inf_{\theta\in\Theta}\{||\varphi_{Z^{j}_{\Delta t}}(u;\theta)-\hat{\varphi}_{Z^{j}_{\Delta t}}(u)||_{2}\}$$

Again, we apply the algorithm to each regime separately.
The integral is computed numerically using a global adaptive quadrature algorithm, where the interval of integration is subdivided and the integration takes place on each subdivided interval. Intervals are further subdivided if the algorithm determines that the integral is not computed to sufficient accuracy.

In table 6 estimates of the model parameters under both regimes and for the two series of underlying assets are shown.
Given a random sample $x=(x_{1},...,x_{n})$ of a random variable $X$ with an associated density function $f(x;\theta)$ of the data $x$ under the real world and unknown parameters $\theta$, maximum likelihood estimation (MLE) is a method used to estimate the vector valued parameter $\theta$ of the model by maximizing the likelihood function:

$$l(\theta;z)=\sum_{k=1}^{n}\log f(z_{k};\theta);\>\>\>\>\theta\in\Theta,$$

with respect to $\theta$. The value of $\theta$ is constrained to $\Theta\subset\mathbb{R}^{4}$, the space of all feasible values of the parameters. The maximum likelihood function $\mathcal{L}$ is primarily a function of the unknown parameters $\theta$. The maximum likelihood estimator is given by:

$$\hat{\theta}=\arg\max_{\theta\in\Theta}l(\theta;x)$$

In tables 7 and 8 estimations based on the empirical m.l.e. for both regimes and models under Gamma and Inverse Gaussian subordinators are shown.

In each case, the values of the volatility $\sigma$ are found to be higher in the second regime, hence justifying the use of a regime switching model. In nearly every method, the value of $|\mu|$ was found to be very small, which is expected as the long term deterministic contribution to the process is expected to be near zero.
In choosing constraints, we set the lower bound of $\sigma,\alpha,\beta$ to be some small number $\epsilon=10^{-6}$. We set the upper bound of $\sigma$ to 5 as the diffusion is expected to be smaller than 1 and for $\alpha,\>\>\beta,$ we set the upper bound to be 100, as the expected value of both Inverse Gaussian and Gamma random variables depends on the ratio $\alpha/\beta$ rather than any particular value for $\alpha$ and $\beta$. The drift $\mu$ is expected to be small, so in most cases, it was constrained to the set $[-1,1]$.

We price European-style options with oil and electricity prices as underlying assets under a switching Levy time-changed noise. These are realistic models that allow to incorporate stylized features in the dynamic of the prices. Our findings show that under this framework a pricing method based on a Fourier-cosine offers an efficient and accurate result when compared with a standard Monte Carlo approach.
In addition, we address the problem of parameter estimation and calibration. To this end we successfully tried different methods based on both, historic and risk-neutral measure.

This research is supported by the Natural Sciences and Engineering Research Council of Canada.