An Option Valuation Formula for Stochastic Volatility Driven by GARCH Processes

Zhongmin Qian¹, Xingcheng Xu^2,3*

¹Mathematical Institute, University of Oxford, Oxford, United Kingdom.
²Shanghai Artificial Intelligence Laboratory, Shanghai, China.
³School of Mathematical Sciences, Peking University, Beijing, China.

DOI: 10.4236/jmf.2023.132015                

Abstract

We have developed a practical and elegant closed-form option pricing formula for general GARCH models using a risk-neutral argument. To estimate the parameters, we propose a procedure and utilize Monte Carlo simulation to calculate the prices. Our formula has been successfully applied to S&P 500 index options and Chinese SSE 50 ETF options, providing empirical evidence that it outperforms the Black-Scholes formula with constant volatility in both the U.S. and Chinese financial markets. While there may be other equivalent martingale measures in this setting, our formula serves as a useful reference for pricing options.

Keywords

Option Pricing, Stochastic Volatility, GARCH, Risk Premia

1. Introduction

In this paper, we present a solution to the question posed by [1] regarding option pricing based on GARCH models. We achieve this by deriving a closed formula.

The theoretical valuation formula for options was established by [2] and [3] using the geometric Brownian motion model for the return of the underlying financial instrument. The work [4] provided a straightforward discrete model approach. Black-Scholes’ formula expresses the price of an option for a financial derivative instrument in terms of the underlying financial asset price, maturity, striking price, risk-free rate, and volatility of the underlying financial instrument. The only unknown variable in the option formula is the volatility of the financial asset, which is assumed to be constant in Black-Scholes-Merton’s model. Therefore, the theoretical valuation formula demonstrates the option price in terms of volatility, while the “implied” volatility provides a prediction of the option price in practice.

Financial economists and practitioners are increasingly concerned with modeling volatility in asset returns. The assumption of constant volatility in Black-Scholes-Merton’s model is unrealistic in actual financial markets, and volatility in real markets is often random. Engle’s autoregressive conditional heteroskedasticity (ARCH) model for time-varying volatility represented a significant breakthrough in the statistical modeling of volatility and option valuation. Engle’s revolutionary approach allows for the systematic explanation of volatility movements over time. Since the seminal work of [5] , many extensions of ARCH models have been developed. The first significant generalization of ARCH models was the family of generalized autoregressive conditional heteroskedasticity (GARCH) models introduced by [6] , which allows for a flexible lag length. Another significant development was the family of exponential GARCH (EGARCH) models proposed by [7] . These models recognize asymmetrical responses to past forecast errors. Many researchers have proposed generalizations of GARCH models, and there is now a wide range of ARCH models available for modeling volatilities.

These GARCH models incorporate key features such as nonlinearity, asymmetry, and long memory properties, making them essential for modeling and forecasting volatility. For further reference, readers may refer to papers such as [8] and others. For simplicity, we refer to these ARCH models as “GARCH” in what follows.

GARCH models have been widely used in finance to measure time-varying volatility and forecast the future movement of an underlying financial asset. This is particularly important for investigating the trade-off between risk and return in financial markets. One notable model in this field is the GARCH-in-mean (GARCH-M) family, with a representative equation given by

$R_t=\mu +c{\sigma }_t^2+{\sigma }_t{z}_t,$ (1)

where $R_t$ represents the logarithm return of the asset, ${\sigma }_t^2$ is a GARCH process, and $\left\{z_t\right\}$ denotes Gaussian white noise at discrete time t. The GARCH-M model was first introduced by [9] . The general form of a GARCH-M model is often written as

$R_t={\mu }_t-\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t,$ (2)

where ${\mu }_t$ may depend linearly or nonlinearly on past returns $\left\{R_{t-1}\mathrm{,}{R}_{t-2}\mathrm{,}\cdots \right\}$ and volatility ${\sigma }_t^2$ . For further details, refer to Section 2 below.

A natural question is how to value option instruments based on GARCH models. Many researchers use the “plug-in” method to value options, which involves evaluating options using the Black-Scholes pricing formula with the forecasting values of volatility, such as the average time-varying volatility of the underlying security over the life of the option, in place of the constant volatility ${\sigma }^2$ . In practice, Monte Carlo simulations of future price paths are applied, and the expected average per-period volatility ${\hat{\sigma }}_{t,T}^2$ between the current time t and the expiration time T is estimated, then plugged into the Black-Scholes pricing formula. Recently, the option pricing formula has been used by combining it with stochastic volatility models, such as those developed in [10] . The Hull-White price is equal to the expected Black-Scholes price integrated over the average instantaneous variance during the life of the option, which calculates the price as

$C_t=\frac{1}{N}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}{C}^{BS}\left({\hat{\sigma }}_{t,T}^2\left(n\right),S_t,K,r\right),$ (3)

where $C^{BS}\left({\sigma }^2\mathrm{,}{S}_t\mathrm{,}K\mathrm{,}r\right)$ is the Black-Scholes price, and ${\hat{\sigma }}_{t,T}^2\left(n\right)$ is the average volatility per period for the nth Monte Carlo simulation of the future price paths. For further details, refer to [11] and other similar sources.

Although the performance of these approaches for option pricing formulas may be good, they are not satisfactory, and we would like to ask what the theoretical option pricing formula is for general GARCH-M models. This is one of the questions raised by [1] .

The work [12] was the first researcher to implement a local risk neutralization within the framework of GARCH models using a quadratic utility function of a representative agent’s risk preference. Subsequently, [13] [14] also adopted the risk neutralization approach for GARCH models, which we believe deserves further investigation.

In this paper, we present closed-form option pricing formulas (see Equations (14) and (21)) based on risk-neutral arguments. These formulas are both elegant and practical for applications, and lead to new insights, such as the relationship among the forecast of future volatility, the return rate, and the risk premia. Empirically, the theoretical option prices based on GARCH models can also be easily applied in real financial markets. We also discuss the estimation problems for these models, and interested readers can refer to [15] [16] [17] , and others for further details on statistical inference.

We have evaluated the performance of our option pricing formula through numerical simulations and empirical applications to the S&P 500 composite stock price index and Chinese SSE 50 ETF. Our results have been compared to those obtained via Black-Scholes’ pricing formula. However, there is still much to be explored and developed in the realm of option pricing based on GARCH-M models. While there may be other equivalent martingale measures in this setting, our formula serves as a useful reference for pricing options. In future work, we plan to investigate other equivalent martingale measures and option valuation for high-frequency econometric data following a GARCH model. For interested readers, we suggest referring to [18] [19] , and [20] for some hints on this aspect.

Our contribution extends previous work on option pricing with GARCH models, which relied on the “plug-in” method or stochastic volatility models. The contributions of our proposed formula are that it leads to new insights such as the relationship among the forecast of future volatility, the return rate, and the risk premia. We have evaluated the performance of our formula through numerical simulations and empirical applications to real financial markets. It serves as a useful reference for pricing options, but there is still much to be explored and developed in the realm of option pricing based on GARCH-M models.

The outline for the remaining sections of the paper is as follows. Section 2 presents a pricing formula for general GARCH-M models. In Section 3, we analyze the option pricing behavior using empirical data from the daily S&P 500 stock index for three GARCH-M type data-generating mechanisms. In Section 4, we evaluate the performance of our theoretical formula using real market prices of S&P 500 stock index options and Chinese SSE 50 ETF options with recent market data. We also compare these prices with those given by the Black-Scholes model.

2. Option Pricing for GARCH-M Models

In this section, we will present option pricing formulas for GARCH-M models. The first model we consider is described by the following equations

$R_t=\mu -\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t\mathrm{,}t\in {\mathbb{Z}}_{+}\mathrm{,}$ (4)

${\sigma }_t^2~\text{GARCHtypeprocess}\mathrm{,}$ (5)

where $R_t=\log \frac{S_t}{S_{t-1}}$ is the logarithm return of assets, $\mu$ is the constant return rate, and $\left\{z_t\mathrm{,}t\in {\mathbb{Z}}_{+}\right\}$ is Gaussian white noise. In this paper, we also call the conditional variance process ${\sigma }_t^2$ as “GARCH type” process. The model above is the simplest one in GARCH-M family, which is however a good representative of this family for option pricing. The important feature of this model is that the return on assets has stochastic volatility which is driven by GARCH processes. Based on risk-neutral arguments, we give the pricing formula in Section 2.1.

Apart from the stochastic volatility in the above model, the return rate can also be generalized such that it depends on the past returns and the volatility of the logarithm return of assets. The first family of generalized return models are the so-called AR(k)-GARCH-in-mean (AR-GARCH-M) models in which the return rate is given by the autoregressive process and the risk premia. That is, they are determined by the equations:

$R_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+V\left({\sigma }_t^2\right)+{\sigma }_t{z}_t,$ (6)

${\sigma }_t^2~\text{GARCHtypeprocess}\mathrm{,}$ (7)

where $V\left({\sigma }_t^2\right)$ is the risk-return relation/risk premia describing the trade-off between risk and return, for example, $V\left({\sigma }_t^2\right)=c{\sigma }_t^2$ , $c\log \left({\sigma }_t^2\right)$ , $c_1{\sigma }_t^2+c_2\sin \left(\omega +c_3{\sigma }_t^2\right)$ etc. (see e.g. [15] for a discussion). Let

${\mu }_t\triangleq {\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+V\left({\sigma }_t^2\right)+\frac{1}{2}{\sigma }_t^2.$ (8)

Then the model can be rewritten as Equations (4) and (5). A closed form of the option pricing formula can be derived too by a similar method used for the first model. More generally, the return rate may depend (linearly or non-linearly) on the past returns and the volatility in a general functional form, that is,

${\mu }_t=\mu \left({\left\{R_s\right\}}_{s<t},{\sigma }_t^2\right).$

The option pricing Formula (21) in the following part is still valid.

2.1. Option Pricing for GARCH-M Models

Let us derive an option pricing formula based on GARCH-M models, described by Equations (4) and (5), where the stochastic volatility ${\sigma }_t^2$ belongs to GARCH family. ${\sigma }_t^2$ may be an ARCH (p), GARCH (p, q), or EGARCH type volatility process. That is,

${\sigma }_t^2={\alpha }_0+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\alpha }_j{\sigma }_{t-j}^2{z}_{t-j}^2,$ (9)

${\sigma }_t^2={\alpha }_0+{\displaystyle \underset{j=1}{\overset{q}{\sum }}}{\beta }_j{\sigma }_{t-j}^2+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\alpha }_j{\sigma }_{t-j}^2{z}_{t-j}^2,$ (10)

or

$\ln \left({\sigma }_t^2\right)={\alpha }_0+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\alpha }_{j}g\left(z_{t-j}\right)+{\displaystyle \underset{j=1}{\overset{q}{\sum }}}{\beta }_j\ln \left({\sigma }_{t-j}^2\right),$ (11)

$g\left(z_{t-j}\right)=\theta z_{t-j}+\lambda \left(\left|z_{t-j}\right|-\mathbb{E}\left|z_{t-j}\right|\right)\mathrm{,}$ (12)

respectively. In general, the stochastic volatility ${\sigma }_t^2$ is determined by the conditional variance of the stochastic process ${\sigma }_t{z}_t$ given the information up to time $t-1$ , that is,

${\sigma }_t^2=\text{Var}\left({\sigma }_t{z}_t\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_{t-1}\right)\mathrm{,}$ (13)

where the sigma field (information up to time t) ${\mathcal{F}}_t=\sigma \left(\left\{z_t\mathrm{,}{z}_{t-1}\mathrm{,}\cdots \right\}\right)$ .

The option pricing formula for this discrete-time stochastic volatility driven by GARCH-M models is obtained by using the idea of arbitrage arguments. The derivation and proofs are given in Appendix A.

Theorem 2.1. Let r be the one plus riskless interest rate on one period, $T-t$ be time to maturity, and K be the striking pricing. Then, for the discrete-time GARCH-M Models (4) and (5), the pricing formula for European call option is given by

$C_t=S_t{\Phi }_{T-t}\left(D_1\right)-K{r}^{-\left(T-t\right)}{\Phi }_{T-t}\left(D_2\right)\mathrm{,}$ (14)

where $S_t$ is the spot price of the underlying asset, ${\Phi }_d(\cdot )$ is d-dimensional standard normal distribution function with independent components, i.e.

${\Phi }_d\left(D\right)={\displaystyle \int }\cdots {\displaystyle {\int }_D}{\left(2\pi \right)}^{-\frac{d}{2}}{\text{e}}^{-\frac{1}{2}{\displaystyle {\sum }_{j=1}^d{u}_j^2}}\text{d}{u}_1\cdots \text{d}{u}_d.$

The integral domains D₁ and D₂ are defined as follows. Let $\rho =-\left(\mu -\log r\right)$ and

$\begin{array}{l}D=\{\left(x_1,\cdots ,x_{T-t}\right)\in {ℝ}^{T-t}:{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-\rho \left(T-t\right)\\ >-\left(\log \frac{S_t}{K{r}^{-\left(T-t\right)}}-\frac{1}{2}{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)\right)\},\end{array}$

then

$D_1=\left\{\left(u_1,\cdots ,u_{T-t}\right)\in {ℝ}^{T-t}:u_j=x_j-\frac{\rho }{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}-{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)\right\}$

and

$D_2=\left\{\left(v_1,\cdots ,v_{T-t}\right)\in {ℝ}^{T-t}:v_j=x_j-\frac{\rho }{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}\right\}.$

In the above notation, ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ , $j\ge 2$ , are obtained by substituting $\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\right)$ in the conditional heteroskedasticity ${\sigma }_{t+j}={\sigma }_{t+j}\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\mathrm{<mo>\; ;\; </mo>}{z}_t\mathrm{,}\cdots \right)$ with ${\tilde{x}}_{j-1}=\left(x_{j-1}\mathrm{,}\cdots \mathrm{,}{x}_1\right)$ , that is,

${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)={\sigma }_{t+j}\left(x_{j-1},\cdots ,x_1;z_t,z_{t-1},\cdots \right),j=2,\cdots ,T-t,$

and

${\sigma }_{t+1}\left({\tilde{x}}_0\right)\equiv {\sigma }_{t+1}\left(z_t\mathrm{,}{z}_{t-1}\mathrm{,}\cdots \right)\mathrm{.}$

We call ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right),j=1,2,\cdots$ , as the conditional volatility function evaluated at time t.

The difficulty in the closed form of the option pricing Formula (14) is the computation of the integral domains D₁ and D₂. By specifying the GARCH stochastic volatility ${\sigma }_t^2$ , we may obtain more information about the option pricing formula theoretically. Let us see some classical examples.

Example 1. Begin with the case in which the volatility is a constant, i.e. ${\sigma }_t^2\equiv {\sigma }^2$ . In this case, our formula coincides with Black-Scholes’ formula:

$C_t=S_{t}N\left(d_1\right)-K{r}^{-\left(T-t\right)}N\left(d_2\right)\mathrm{,}$ (15)

where

$d_1=\frac{\log \frac{S_t}{K{r}^{-\left(T-t\right)}}+\frac{1}{2}{\sigma }^2\left(T-t\right)}{\sigma \sqrt{T-t}}\mathrm{,}$ (16)

$d_2=\frac{\log \frac{S_t}{K{r}^{-\left(T-t\right)}}-\frac{1}{2}{\sigma }^2\left(T-t\right)}{\sigma \sqrt{T-t}},$ (17)

and $N(\cdot )$ is standard normal distribution function. As we can see that the integral domain $D_1$ is a half hyperplane, that is,

$D_1=\left\{\left(u_1,\cdots ,u_{T-t}\right)\in {ℝ}^{T-t}:\sigma {\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{u}_j>-\left[\log \frac{S_t}{K{r}^{-\left(T-t\right)}}+\frac{1}{2}{\sigma }^2\left(T-t\right)\right]\right\},$

then

${\Phi }_{T-t}\left(D_1\right)=\mathbb{P}\left(Z_{T-t}\in D_1\right)=\mathbb{P}\left(Z>-\frac{\log \frac{S_t}{K{r}^{-\left(T-t\right)}}+\frac{1}{2}{\sigma }^2\left(T-t\right)}{\sigma \sqrt{T-t}}\right)=N\left(d_1\right),$

where $Z_{T-t}~N_{T-t}\left(\mathrm{0,}{I}_{T-t}\right)$ and $Z:=\frac{1}{\sqrt{T-t}}{\displaystyle {\sum }_{j=1}^{T-t}}{Z}_{T-t}^{\left(j\right)}~N\left(0,1\right)$ . The integral domain D₂ can also be obtained similarly.

If the stochastic volatility ${\sigma }_t^2$ is modelled by a general GARCH process, we can also get some information on the integral domains.

Example 2 When the stochastic volatility is an ARCH (1) process

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2,{\alpha }_0>0,{\alpha }_1\ge 0,$

then the volatility process ${\sigma }_t^2$ has the following representation:

${\sigma }_{t+j}^2={\alpha }_0+{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\alpha }_0{\alpha }_1^i{\displaystyle \underset{\mathcal{l}=1}{\overset{i}{\prod }}}{z}_{t+j-\mathcal{l}}^2+{\alpha }_1^j{\displaystyle \underset{\mathcal{l}=1}{\overset{j-1}{\prod }}}{z}_{t+j-\mathcal{l}}^2{\sigma }_t^2{z}_t^2,\forall j\ge 1,$

where the conventions that ${\displaystyle {\sum }_{i=1}^0}:=0$ , ${\displaystyle {\prod }_{i=1}^0}:=1$ have been used. The integral domains D₁ and D₂ in Theorem 2.1 are determined for this case by

${\sigma }_{t+1}\left({\tilde{x}}_0\right)\equiv {\sigma }_{t+1}\left(z_t,z_{t-1},\cdots \right)={\left({\alpha }_0+{\alpha }_1{\sigma }_t^2{z}_t^2\right)}^{\frac{1}{2}},$

$\begin{array}{c}{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)={\sigma }_{t+j}\left(x_{j-1},\cdots ,x_1;z_t,z_{t-1},\cdots \right)\\ ={\left({\alpha }_0+{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\alpha }_0{\alpha }_1^i{\displaystyle \underset{\mathcal{l}=1}{\overset{i}{\prod }}}{x}_{j-\mathcal{l}}^2+{\alpha }_1^j{\displaystyle \underset{\mathcal{l}=1}{\overset{j-1}{\prod }}}{x}_{j-\mathcal{l}}^2{\sigma }_t^2{z}_t^2\right)}^{\frac{1}{2}}.\end{array}$

Example 3. 1) If the stochastic volatility ${\sigma }_t^2$ is the GARCH (1, 1) model, that is,

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2,{\alpha }_0>0,{\alpha }_1,{\beta }_1\ge 0.$

Then the volatility process ${\sigma }_t^2$ has the following representation:

${\sigma }_{t+j}^2={\alpha }_0+{\alpha }_0{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\displaystyle \underset{\mathcal{l}=1}{\overset{i}{\prod }}}\left({\alpha }_1{z}_{t+j-\mathcal{l}}^2+{\beta }_1\right)+{\displaystyle \underset{\mathcal{l}=1}{\overset{j-1}{\prod }}}\left({\alpha }_1{z}_{t+j-\mathcal{l}}^2+{\beta }_1\right)\left({\alpha }_1{\sigma }_t^2{z}_t^2+{\beta }_1{\sigma }_t^2\right),$

and the integral domains $D_1$ and $D_2$ in Theorem 2.1 are determined by

${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)=({\alpha }_0+{\alpha }_0{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\displaystyle \underset{\mathcal{l}=1}{\overset{i}{\prod }}}\left({\alpha }_1{x}_{j-\mathcal{l}}^2+{\beta }_1\right)+{{\displaystyle \underset{\mathcal{l}=1}{\overset{j-1}{\prod }}}\left({\alpha }_1{x}_{j-\mathcal{l}}^2+{\beta }_1\right)\left({\alpha }_1{\sigma }_t^2{z}_t^2+{\beta }_1{\sigma }_t^2\right))}^{\frac{1}{2}}.$

2) For the GARCH (p, p) models

${\sigma }_t^2={\alpha }_0+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\alpha }_j{\sigma }_{t-j}^2{z}_{t-j}^2+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\beta }_j{\sigma }_{t-j}^2,$

the stochastic volatility ${\sigma }_{t+j}^2$ is given by

$\begin{array}{c}{\sigma }_{t+j}^2={\alpha }_0+{\alpha }_0{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}\left({\displaystyle \underset{k_1=1}{\overset{p}{\sum }}}\cdots {\displaystyle \underset{k_i=1}{\overset{p}{\sum }}}{\displaystyle \underset{\mathcal{l}=1}{\overset{i}{\prod }}}\left({\alpha }_{k_{\mathcal{l}}}{z}_{t+j-{\displaystyle {\sum }_{q=1}^{\mathcal{l}}{k}_q}}^2+{\beta }_{k_{\mathcal{l}}}\right)\right)\\ +{\displaystyle \underset{k_1=1}{\overset{p}{\sum }}}\cdots {\displaystyle \underset{k_j=1}{\overset{p}{\sum }}}{\displaystyle \underset{\mathcal{l}=1}{\overset{j}{\prod }}}\left({\alpha }_{k_{\mathcal{l}}}{z}_{t+j-{\displaystyle {\sum }_{q=1}^{\mathcal{l}}{k}_q}}^2+{\beta }_{k_{\mathcal{l}}}\right){\sigma }_{t+j-{\displaystyle {\sum }_{q=1}^j{k}_q}}^2.\end{array}$

The models include all ARCH (p) processes by taking ${\beta }_j\equiv 0$ . Substituting $\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\right)$ in the conditional heteroskedasticity

${\sigma }_{t+j}^2={\sigma }_{t+j}^2\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\mathrm{<mo>\; ;\; </mo>}{z}_t\mathrm{,}\cdots \right)$

with ${\tilde{x}}_{j-1}=\left(x_{j-1},\cdots ,x_1\right)$ , we obtain the conditional volatility function ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ , which determines the integral domains $D_i,i=1,2$ in Theorem 2.1.

Example 4. For the EGARCH (p, q) models,

$\ln \left({\sigma }_t^2\right)={\alpha }_0+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{\alpha }_j\left(\left|z_{t-j}\right|-\mathbb{E}\left|z_{t-j}\right|\right)+{\displaystyle \underset{j=1}{\overset{p}{\sum }}}{{\alpha }^{\prime }}_j{z}_{t-j}+{\displaystyle \underset{j=1}{\overset{q}{\sum }}}{\beta }_j\ln \left({\sigma }_{t-j}^2\right),$

then the representation of $\ln \left({\sigma }_{t+j}^2\right)$ is

$\begin{array}{c}\ln \left({\sigma }_{t+j}^2\right)=\left({\alpha }_0+{\displaystyle \underset{i=1}{\overset{p}{\sum }}}{\alpha }_i\left(\left|z_{t+j-i}\right|-\mathbb{E}\left|z_{t+j-i}\right|\right)+{\displaystyle \underset{i=1}{\overset{p}{\sum }}}{{\alpha }^{\prime }}_i{z}_{t+j-i}\right)\\ \times \left(1+{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{q}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_i=1}{\overset{q}{\sum }}}{\beta }_{{\mathcal{l}}_1}{\beta }_{{\mathcal{l}}_2}\cdots {\beta }_{{\mathcal{l}}_i}\right)\\ +{\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{q}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_j=1}{\overset{q}{\sum }}}{\beta }_{{\mathcal{l}}_1}{\beta }_{{\mathcal{l}}_2}\cdots {\beta }_{{\mathcal{l}}_j}\ln \left({\sigma }_{t+j-{\displaystyle {\sum }_{k=1}^j{\mathcal{l}}_k}}^2\right).\end{array}$

The conditional volatility function ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ follows.

For other GARCH type processes, e.g. IGARCH, FIGARCH, IEGARCH, FIEGARCH etc. explicit computations on the condition volatility functions can be done too. Though the conditional volatility function ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ looks very complicated, but it is easy to be implemented iteratively by computer software, e.g. Matlab, Python and so on.

Let us make some comments on the theoretical option valuation Formula (14) in Theorem 2.1. The first observation is the dependence of option price on the volatility. From our formula, we can see that the option price (at time t) depends on the future volatility of the underlying security i.e. ${\sigma }_{t+1}^2\mathrm{,}{\sigma }_{t+2}^2\mathrm{,}\cdots \mathrm{,}{\sigma }_T^2$ through the option duration. Since the volatility is varying from one period to another period, the better forecast of the future volatility is made, the more accurate the option will be valued. The GARCH type time-varying volatility is thus a method to model and forecast the dynamics of volatility. The family of GARCH models is huge, so that one has a great number of choices to select an appropriate GARCH type process for a specific underlying security to make a better forecast of the future volatility.

The second point we want to emphasize is that the return rate $\mu$ appears in the theoretical option valuation Formula (14), which is quiet different from Black-Scholes’ option pricing Formula (15). However, by an informal analysis of sensitivity on the return rate, we find that there is a small impact on valuation of options for a large range of return rate.

Figure 1 plots the relationship between daily return rate $\mu$ and option price C computed from Formula (14). We find that for reasonable return rate $\mu$ , it has little impact on the option price. Table 1 shows daily return rates and their annualized ones.

Figure 1. The effect of daily return rate $\mu$ on option price. Parameters are estimated from S&P 500 Index, see Model 1 in Section C.2. The annualized riskless interest rate is taken to 2.5%, string price $K=S_0$ (i.e. at-the-money option), and time-to-maturity $T=30$ days.

Table 1. Daily and annualized return rate $\mu$ .

The mathematical reason for appearance of return rate in our formula is that the martingale $\hat{M}$ (M with drift) in Equation (38) under risk-neutral measure $\mathbb{Q}$ is not same in law with the martingale M in Equation (35) under physical measure $\mathbb{P}$ , which is different from the “Brownian motion” case, i.e. Brownian motion with drift under measure $\mathbb{Q}$ is still Brownian motion. In econometric words, for the GARCH-M case we are considering, the return innovation distribution is changed to another distribution by changing physical measure $\mathbb{P}$ to risk-neutral measure $\mathbb{Q}$ . This is different from the constant volatility case, in which the return innovation distribution keeps unchanged under the transformation of measures.

Another point we want to explain is about the integral domains D₁ and D₂. The dimension of the integral domains changes according to time-to-maturity. When the volatility is constant, ${\Phi }_{T-t}\left(D_i\right)$ become $N\left(d_i\right)$ , $i=1,2$ , due to the Gaussianity of linear combination of Gaussian/normal random variables. The integral domains of constant volatility case are always half planes with linear partition. However, for the GARCH type volatility, the borderline is no longer linear/flat but a curve/surface. A 2-dimensional example (with GARCH (1, 1) volatiltiy) may provide a heuristic knowledge of the shape. In Figure 2, we simulated 100,000 points (all the colored points) with 2-dimensional standard normal distribution. According to the three-sigma rule of thumb, there are 99.73% of the values which lie within three standard deviations of the mean. So all the colored points cluster as shown in Figure 2. The red parts are the integral domains

Figure 2. Integral domains D₁ (the left panel) and D₂ (the right panel) in dimension 2 space.

D₁ (the left panel) and D₂ (the right panel). One can see very clearly the borderline and the shape of domains from this figure. By the way, one can actually compute out values of ${\Phi }_{T-t}\left(D_i\right)$ , which are the ratio of numbers of red points to all colored points by Manto Carlo method.

2.2. Option Pricing for AR-GARCH-M Models

Now, we consider the more general AR(k)-GARCH-M models:

$R_t={\mu }_t-\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t\mathrm{,}t\in {\mathbb{Z}}_{+}\mathrm{,}$ (18)

${\mu }_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j},$ (19)

${\sigma }_t^2~\text{GARCHtypeprocess}\mathrm{,}$ (20)

where $R_t=\log \frac{S_t}{S_{t-1}}$ is the logarithm return of assets, $\left\{z_t\mathrm{,}t\in {\mathbb{Z}}_{+}\right\}$ is the Gaussian white noise, and ${\sigma }_t^2$ belongs to GARCH family. The return rate ${\mu }_t$ can be more general as long as ${\mu }_t$ is measurable with respect to ${\mathcal{F}}_{t-1}$ . The riskless interest rate may depend on the time but should not be random.

Then we have the following European call option pricing formula, the proof of which is given in Appendix B. (For simplicity, we still assume that the riskless interest rate is a constant.)

Theorem 2.2. Let r be one plus riskless interest rate on one period, $T-t$ be time to maturity, and K be the striking pricing. Then, for the discrete-time AR(k)-GARCH-M models, the pricing formula for European call option is

$C_t=S_t{\Phi }_{T-t}\left(D_1\right)-K{r}^{-\left(T-t\right)}{\Phi }_{T-t}\left(D_2\right)\mathrm{,}$ (21)

where the integral domains D₁ and D₂ are defined as follows. Let ${\rho }_{t+j}\left({\tilde{x}}_{j-1}\right)=-\left({\mu }_{t+j}\left({\tilde{x}}_{j-1}\right)-\log r\right)$ and

$\begin{array}{l}D=\{\left(x_1,\cdots ,x_{T-t}\right)\in {ℝ}^{T-t}:{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-{\rho }_{t+j}\left({\tilde{x}}_{j-1}\right)\left(T-t\right)\\ >-\left(\log \frac{S_t}{K{r}^{-\left(T-t\right)}}-\frac{1}{2}{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)\right)\},\end{array}$

then

$D_1=\left\{\left(u_1,\cdots ,u_{T-t}\right)\in {ℝ}^{T-t}:u_j=x_j-\frac{{\rho }_{t+j}\left({\tilde{x}}_{j-1}\right)}{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}-{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)\right\}$

and

$D_2=\left\{\left(v_1,\cdots ,v_{T-t}\right)\in {ℝ}^{T-t}:v_j=x_j-\frac{{\rho }_{t+j}\left({\tilde{x}}_{j-1}\right)}{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}\right\}.$

In the above notation, ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ , $j=1,2,\cdots ,T-t$ , are conditional volatility functions evaluated at time t, and ${\mu }_{t+j}\left({\tilde{x}}_{j-1}\right)$ , $j=1,2,\cdots ,T-t$ , are conditional return rate functions evaluated at time t which are defined in the same way as conditional volatility functions.

Consider as an example, the AR (1)-GARCH (1, 1)-M model:

$R_t={\gamma }_0+{\gamma }_1{R}_{t-1}-\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2.$

Theorem 2.2 gives the closed form of the theoretical option pricing formula $C_t$ for this model, where the conditional return rate functions ${\mu }_{t+j}\left({\tilde{x}}_{j-1}\right)$ , $j=1,2,\cdots ,T-t$ , can be computed explicitly.

In fact, since ${\mu }_t={\gamma }_0+{\gamma }_1{R}_{t-1}$ ,

${\mu }_{t+j}={\displaystyle \underset{i=0}{\overset{j-1}{\sum }}}{\gamma }_0{\gamma }_1^i+{\gamma }_1^j{R}_t+{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\gamma }_1^{j-i}\left(-\frac{1}{2}{\sigma }_{t+i}^2+{\sigma }_{t+i}{z}_{t+i}\right).$ (22)

Substituting $\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\right)$ in the return rate ${\mu }_{t+j}$ with ${\tilde{x}}_{j-1}=\left(x_{j-1},\cdots ,x_1\right)$ , we may obtain the conditional return rate function

${\mu }_{t+j}\left({\tilde{x}}_{j-1}\right)=\frac{{\gamma }_0\left(1-{\gamma }_1^j\right)}{1-{\gamma }_1}+{\gamma }_1^j{R}_t+{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}{\gamma }_1^{j-i}\left(-\frac{1}{2}{\sigma }_{t+i}^2\left({\tilde{x}}_{i-1}\right)+{\sigma }_{t+i}\left({\tilde{x}}_{i-1}\right)x_i\right),$ (23)

where $\left\{{\sigma }_{t+i}\left({\tilde{x}}_{i-1}\right)\right\}$ are the conditional volatility functions defined in Section 2.1.

For a general AR(k)-GARCH-M model given by

$R_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}-\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t,$

${\sigma }_t^2~\text{GARCHtypeprocesses}\mathrm{,}$

the return rate has the following representation

$\begin{array}{c}{\mu }_{t+j}={\gamma }_0+{\gamma }_0{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}\left({\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{k}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_i=1}{\overset{k}{\sum }}}{\gamma }_{{\mathcal{l}}_1}{\gamma }_{{\mathcal{l}}_2}\cdots {\gamma }_{{\mathcal{l}}_i}\right)\\ +{\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{k}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_j=1}{\overset{k}{\sum }}}{\gamma }_{{\mathcal{l}}_1}{\gamma }_{{\mathcal{l}}_2}\cdots {\gamma }_{{\mathcal{l}}_j}{R}_{t+j-{\displaystyle {\sum }_{q=1}^j{\mathcal{l}}_q}}\\ +{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}({\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{k}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_i=1}{\overset{k}{\sum }}}{\gamma }_{{\mathcal{l}}_1}{\gamma }_{{\mathcal{l}}_2}\cdots {\gamma }_{{\mathcal{l}}_i}\\ \times \left(-\frac{1}{2}{\sigma }_{t+j-{\displaystyle {\sum }_{q=1}^i{\mathcal{l}}_q}}^2+{\sigma }_{t+j-{\displaystyle {\sum }_{q=1}^i{\mathcal{l}}_q}}{z}_{t+i-{\displaystyle {\sum }_{q=1}^j{\mathcal{l}}_q}}\right)).\end{array}$ (24)

By replacing $\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\right)$ in the return rate ${\mu }_{t+j}$ with ${\tilde{x}}_{j-1}=\left(x_{j-1}\mathrm{,}\cdots \mathrm{,}{x}_1\right)$ , and replacing conditional volatility ${\sigma }_{t+j}$ with associated conditional volatility functions ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ , we may calculate the conditional return rate function ${\mu }_{t+j}\left({\tilde{x}}_{j-1}\right)$ at time t.

The risk-return relation can also be generalized such as linear or nonlinear depending on ${\sigma }_t^2$ . For example, $V\left({\sigma }_t^2\right)=c{\sigma }_t^2$ , $c\log \left({\sigma }_t^2\right)$ , $c_1{\sigma }_t^2+c_2\sin \left(\omega +c_3{\sigma }_t^2\right)$ etc. (see e.g. [15] and references therein) have been discussed. The AR(k)-GARCH-M models can be therefore generalized, such as

$R_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+V\left({\sigma }_t^2\right)+{\sigma }_t{z}_t,$ (25)

${\sigma }_t^2~\text{GARCHtypeprocess}\mathrm{,}$ (26)

where $V\left({\sigma }_t^2\right)$ is the risk premia describing the tradeoff between risk and return. We may define

${\mu }_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+V\left({\sigma }_t^2\right)+\frac{1}{2}{\sigma }_t^2,$ (27)

so that the model is still described by Equations (18)-(20).

Since ${\mu }_t$ is measurable with respect to ${\mathcal{F}}_{t-1}$ , the option pricing Formula (21) in Theorem 2.2 still works for this generalized model, with the conditional volatility functions ${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)$ evaluated at time t as in Section 2.1. The conditional return rate

$\begin{array}{c}{\mu }_{t+j}={\gamma }_0+V\left({\sigma }_{t+j}^2\right)+\frac{1}{2}{\sigma }_{t+j}^2+{\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{k}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_j=1}{\overset{k}{\sum }}}{\gamma }_{{\mathcal{l}}_1}{\gamma }_{{\mathcal{l}}_2}\cdots {\gamma }_{{\mathcal{l}}_j}{R}_{t+j-{\displaystyle {\sum }_{q=1}^j{\mathcal{l}}_q}}\\ +\left({\gamma }_0+V\left({\sigma }_{t+j}^2\right)+\frac{1}{2}{\sigma }_{t+j}^2\right){\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}\left({\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{k}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_i=1}{\overset{k}{\sum }}}{\gamma }_{{\mathcal{l}}_1}{\gamma }_{{\mathcal{l}}_2}\cdots {\gamma }_{{\mathcal{l}}_i}\right)\\ +{\displaystyle \underset{i=1}{\overset{j-1}{\sum }}}({\displaystyle \underset{{\mathcal{l}}_1=1}{\overset{k}{\sum }}}\cdots {\displaystyle \underset{{\mathcal{l}}_i=1}{\overset{k}{\sum }}}{\gamma }_{{\mathcal{l}}_1}{\gamma }_{{\mathcal{l}}_2}\cdots {\gamma }_{{\mathcal{l}}_i}\\ \times \left(-\frac{1}{2}{\sigma }_{t+j-{\displaystyle {\sum }_{q=1}^j{\mathcal{l}}_q}}^2+{\sigma }_{t+j-{\displaystyle {\sum }_{q=1}^j{\mathcal{l}}_q}}{z}_{t+i-{\displaystyle {\sum }_{q=1}^j{\mathcal{l}}_q}}\right)).\end{array}$

Like the conditional volatility ${\sigma }_{t+j}$ , the conditional return rate ${\mu }_{t+j}$ looks very complicated, it is in fact very easy to be implemented like ${\sigma }_{t+j}$ by e.g. Matlab, Python and so on. The informal analysis for the sensitivity of theoretical valuation formula with respect to return rate $\mu$ reveals little sensitivity in a large range. Thus, in real applications, we do not really need the calculation of these conditional return rate functions for the purpose of saving computation costs.

Remark 2.3. We would like to point out that our option valuation formula in Theorem 2.2 is still valid for the following model:

$R_t={\mu }_t-\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t\mathrm{,}$ (28)

${\mu }_t=\mu \left({\left\{R_s\right\}}_{s<t},{\sigma }_t^2\right),$ (29)

${\sigma }_t^2~\text{GARCHtypeprocess}\mathrm{,}$ (30)

The return rate ${\mu }_t$ may depend on the past returns and the volatility linearly or both nonlinearly.

3. Option Pricing

In this section, we analysis the option valuation formula

$C_t=S_t{\Phi }_{T-t}\left(D_1\right)-K{r}^{-\left(T-t\right)}{\Phi }_{T-t}(\; D\; 2\; )$

for European call options on the S&P 500 index. The quantities ${\Phi }_{T-t}\left(D_i\right),i=1,2$ , are computed by Monte Carlo method. Since

${\Phi }_{T-t}\left(D_i\right)={\displaystyle \int }\cdots {\displaystyle {\int }_{D_i}}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-\frac{1}{2}{\displaystyle {\sum }_{j=1}^{T-t}{y}_j^2}}\text{d}{y}_1\cdots \text{d}{y}_{T-t},$

by independently simulating a number of $\left(T-t\right)$ -dimensional standard normal random variables, one may obtain simulated values of ${\Phi }_{T-t}\left(D_i\right)$ .

In Appendix C, we provide details on our proposed procedure for estimating parameters. We utilize the three models specified in that section to investigate the pricing behavior given by our option valuation formula. The underlying asset is S&P 500 index. The present time t is the day on December 29, 2017. We consider at-the-money options for this moment, that is, the striking price $K=S\left(t\right)$ . The number of days to maturity $T-t$ varies from 1 to 60. The riskless annual continuously compounded interest rate $\tilde{r}=2.5\%$ , i.e. $r={\text{e}}^{\tilde{r}\Delta t}=1+6.85\times {10}^{-5}$ with $\Delta t=1/365$ in our formula. The results are presented in Figure 3. We also included the simulated prices calculated with the traditional Black-Scholes pricing formula (see [2] ). The roughness of pricing lines in Figure 3 is caused by the error of Monte Carlo simulation. From the figure, we can see that, for these three GARCH (1, 1)-M models, the simulated prices based on our option valuation formula are almost same and higher than the celebrated Black-Scholes prices. Since the only difference of Model 1, 2 and 3 is the return rate, it reveals to some degree that there is little sensitivity with respect to the return rate models.

Figure 3. Option prices as a function of maturity by different models.

4. Empirical Applications

In this section, we show the performance of our theoretical option valuation formula when applying it to S&P 500 index options, and SSE 50 ETF options traded on Shanghai Stock Exchange (SSE), China.

4.1. S&P 500 Index Options

Based on the in-sample analysis in the previous section, we shall compare our option valuation formula with the real financial market. We use GARCH-in-mean data-generating mechanism as Model 1 (with GARCH (1, 1) volatility) in Appendix C.2 for the underlying S&P 500 index.

We randomly pick two heavily traded SPX call options expired on 2018. One option is issued on October 23, 2017, exercised on February 15, 2018 with striking price $K=2650$ while another one is issued on April 26, 2017, exercised on March 15, 2018 with striking price $K=2700$ . All these historical data are provided by Bloomberg. The riskless interest rate $\tilde{r}$ over the life of the options is taken to be 2.5%, which is annual continuously compounded and approximately equals to the interest rate of US treasury bills at the same period. So the interest $r={\text{e}}^{\tilde{r}\Delta t}=1+6.85\times {10}^{-5}$ with $\Delta t=1/365$ in our theoretical valuation formula. As a remark, some informal analysis revealed little sensitivity to the choice of the riskless interest rate.

For the first SPX call option C2650, we computed the model prices from December 1, 2017 to February 15, 2018 by Model 1 data-generating mechanism in Appendix C.2 using the rolling windows method, that is, the parameters of Model 1 are estimated by the past one year daily data of the underlying index SPX before the date we are computing. We also compared the model prices calculated by the celebrated Black-Scholes option pricing formula. We showed the results in Figure 4. In this figure, the three prices—market prices, Black-Scholes model prices and GARCH-M Model 1 prices, of SPX call option C2650 are compared. For the other SPX call option C2700, we use the same approach to obtain the model prices from January 9, 2018 to March 15, 2018, which are shown in Figure 5.

The GARCH-M Model 1 prices are clearly closer to the market prices than the Black-Scholes model prices for most of the trading days. It means that the forecast of future volatility works better with our method. The performance of Model 1 is much better than the traditional method especially when the market volatility is not pleasant, e.g. February 2018. For better performance, one can try to look for a better forecast of future volatility, for example, employing FIGARCH, EGARCH process to model the volatility rather than just using GARCH (1, 1) model.

4.2. Chinese Financial Market

In this subsection, we show results by applying our option valuation formula to China financial market. Options were introduced to the Chinese financial market in February 9, 2015. Chinese SSE 50 ETF option is the first and only standardized option traded on Chinese market. The underlying asset is the SSE 50 ETF. All 50 ETF options are European options and are traded on the Shanghai Stock Exchange (SSE). For the underlying asset SSE 50 ETF, it is also traded on the Shanghai Stock Exchange and was the first ETF traded in China. This ETF tracks the SSE 50 index, which includes 50 of the most active and reputable stocks listed on the Shanghai Stock Exchange. It is one of the most heavily traded ETFs in China.

The data of the underlying asset we use are the daily close prices of SSE 50 ETF from September 1, 2016 through December 29, 2017. There are 324 prices in total. The dataset was provided by Wind Info, Inc. Figure 6 plots the price process and the associated logarithm return of SSE 50 ETF for this dataset.

Figure 4. SPX option (C2650) prices from December 1, 2017 to February 15, 2018. $\tilde{r}=2.5\%$ , $K=2650$ , issued on October 23, 2017, exercised on February 15, 2018.

Figure 5. SPX option (C2700) prices from January 9, 2018 to March 15, 2018. $\tilde{r}=2.5\%$ , $K=2700$ , issued on April 26, 2017, exercised on March 15, 2018.

Our aim in this section is to compute the market prices with our option valuation formula for Chinese financial market. We also compared the prices with the celebrated Black-Scholes option pricing formula as previous subsection.

Figure 6. Prices and log-return of SSE 50 ETF (9/1/2016-12/29/2017).

Figure 7. SSE 50 ETF options’ prices from September 1, 2017 to September 29, 2017. $\tilde{r}=3.9\%$ , $K=2.70,2.75,2.80$ , issued on August 24, 2017, exercised on October 25, 2017.

The market prices of these SSE 50 ETF options at the end of each day from September 1, 2017 to September 29, 2017 are available. There are 21 trading days in this month, and there are 9 striking prices for these options: 1 at-the-money, 4 out-of-the-money and 4 in-the-money. We use three actively traded options, whose symbols are 10000993.SH, 10000994.SH and 10000995.SH. Both started on August 24, 2017, and were exercised on October 25, 2017. The exercise prices K are RMB 2.70, 2.75 and 2.80, respectively. We calculated the model prices of these options for every trading day of the whole of September by using estimates of the parameters of GARCH-M models based on past data of underlying asset SSE 50 ETF. The riskless annual continuously compound interest rate $\tilde{r}$ is taken to be 3.9%, which is the constant interest rate of Chinese bonds at the same period.

Standing on the date September 1, 2017, on which the price of SSE 50 ETF is RMB 2.76, we want to decide the option model price on September 1 by using the past one year’s daily prices of SSE 50 ETF, i.e. the data from September 1, 2016 to September 1, 2017. By using these past one year’s data, we first obtain estimated parameters for GARCH-M models as Model 1 in Appendix C.2, then we can work out a model price with our option pricing formula. For the next trading day, we use the rolling windows method, we can compute a new option price for it. The results are shown in Figure 7.

It is apparent to see that the GARCH-M Model 1 prices are between the market prices and the Black-Scholes model prices for most of the time.

5. Conclusions

In this paper, we present option pricing formulas for general GARCH-M models based on risk-neutral arguments. These formulas are not only elegant, but also practical for real-world applications. We propose a parameter estimation procedure and use Monte Carlo simulation to evaluate the prices (see Appendix C). We demonstrate the pricing behavior of these theoretical formulas for S&P 500 index options based on three data-generating mechanisms. Empirical evidence suggests that the performance of these theoretical pricing formulas is very good compared to real market prices, and better than the celebrated Black-Scholes pricing formula, in both the U.S. stock market and the Chinese financial market.

However, our proposed formula has limitations. One challenge is the computation of the multidimensional Gaussian integral over non-standard domains D₁ and D₂, while another is the lack of consideration for implied volatility, a crucial factor in option pricing. Additionally, it may be necessary to investigate other equivalent martingale measures in this setting.

Despite these limitations, our formula serves as a valuable reference for pricing options. Further investigation into its properties and implications, as well as extending the empirical analysis to other stocks with appropriate GARCH-type stochastic volatility, would be beneficial for both research and real-world applications.

Acknowledgements

The authors thank anonymous reviewers for helpful comments.

Appendix A. Proof of Theorem 2.1

Assume that the riskless interest rate is constant, and let interest r denote one plus the riskless interest rate over one period. In risk-neutral world (the assumptions underlying the risk-neutral argument could be found in e.g. [2] [3] ), the expected rate of return on the asset would be the riskless interest rate, so we want to find a risk-neutral measure $\mathbb{Q}$ such that

${\mathbb{E}}^{\mathbb{Q}}\left[S_{t+k}\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]=r^k{S}_t\mathrm{,}\forall k\ge \mathrm{1,}t\ge 0.$ (31)

The price of an option for this underlying asset is then

$C_t=r^{-\left(T-t\right)}{\mathbb{E}}^{\mathbb{Q}}\left[C_T\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]\mathrm{.}$ (32)

For European call option, $C_T=\max \left\{0,S_T-K\right\}$ . Actually, under measure $\mathbb{Q}$ , the two processes $\left\{\frac{S_t}{r^t}\mathrm{,}t\in {\mathbb{Z}}_{+}\right\}$ and $\left\{\frac{C_t}{r^t},t=0,1,2,\cdots ,T\right\}$ are discrete-time martingales. The question left is how to construct the risk-neutral measure $\mathbb{Q}$ .

For GARCH-M models, by Equation (31) we have

${\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{{\displaystyle {\sum }_{i=t+1}^{t+k}{R}_i}}|{\mathcal{F}}_t\right]=r^k,\forall k\ge 1,t\ge 0.$ (33)

Then

${\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{{\displaystyle {\sum }_{i=t+1}^{t+k}{\sigma }_i{z}_i}+k\left(\mu -\log r\right)-\frac{1}{2}{\displaystyle {\sum }_{i=t+1}^{t+k}{\sigma }_i^2}}|{\mathcal{F}}_t\right]=1,\forall k\ge 1,t\ge 0.$ (34)

In the following, we are going to construct the risk-neutral measure $\mathbb{Q}$ satisfying Equation (34). Define measure $\mathbb{P}$ as the distribution of Gaussian white noise $\left\{z_i\mathrm{,}i\ge 0\right\}$ , and

$M_t\triangleq M_0+{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\sigma }_i{z}_i,$ (35)

$Z_t\triangleq \exp \left\{{\displaystyle \underset{i=1}{\overset{t}{\sum }}}\frac{\rho }{{\sigma }_i}{z}_i-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{t}{\sum }}}\frac{{\rho }^2}{{\sigma }_i^2}\right\},\rho =-\left(\mu -\log r\right).$ (36)

We can verify that $M_t$ and $Z_t$ are martingales under measure $\mathbb{P}$ .

Let measure $\mathbb{Q}$ be

$\mathbb{Q}\left(A\right)\triangleq {\mathbb{E}}^{\mathbb{P}}\left({\chi }_A{Z}_t\right),\forall A\in {\mathcal{F}}_t,$ (37)

with ${\chi }_A\left(\omega \right)=1$ , if $\omega \in A$ , otherwise ${\chi }_A\left(\omega \right)=0$ , and

${\hat{M}}_t\triangleq M_t-\rho t=M_0+{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\sigma }_i{z}_i+\left(\mu -\log r\right)t.$ (38)

Lemma A.1. For any $t\mathrm{,}k\in {\mathbb{Z}}_{+}$ , if any random variable X is measurable in ${\mathcal{F}}_{t+k}$ , then

${\mathbb{E}}^{\mathbb{Q}}\left[X|{\mathcal{F}}_t\right]=\frac{1}{Z_t}{\mathbb{E}}^{\mathbb{P}}\left[X{Z}_{t+k}|{\mathcal{F}}_t\right].$ (39)

Proof. For any $A\in {\mathcal{F}}_t$ , by the definition of conditional expectation, we have

$\begin{array}{c}{\mathbb{E}}^{\mathbb{Q}}\left[{\chi }_A\frac{1}{Z_t}{\mathbb{E}}^{\mathbb{P}}\left[X{Z}_{t+k}|{\mathcal{F}}_t\right]\right]={\mathbb{E}}^{\mathbb{P}}\left[{\chi }_A{\mathbb{E}}^{\mathbb{P}}\left[X{Z}_{t+k}|{\mathcal{F}}_t\right]\right]\\ ={\mathbb{E}}^{\mathbb{P}}\left[{\chi }_{A}X{Z}_{t+k}\right]\\ ={\mathbb{E}}^{\mathbb{Q}}\left[{\chi }_{A}X\right].\end{array}$

This concluded this lemma.

Theorem A.2. Under measure $\mathbb{Q}$ , stochastic processes ${\hat{M}}_t$ and

$\exp \left\{{\hat{M}}_t-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\sigma }_i^2\right\}$

are ${\mathcal{F}}_t$ -martingales.

Proof. 1) We show that ${\hat{M}}_t$ is a $\mathbb{Q}$ -martingale below. We first prove that ${\hat{M}}_t{Z}_t$ is a $\mathbb{P}$ -martingale. Since

$\begin{array}{l}{\mathbb{E}}^{\mathbb{P}}\left[{\hat{M}}_{t+1}{Z}_{t+1}|{\mathcal{F}}_t\right]={\mathbb{E}}^{\mathbb{P}}\left[\left({\hat{M}}_t+{\sigma }_{t+1}{z}_{t+1}-\rho \right)Z_t{\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}|{\mathcal{F}}_t\right]\\ ={\hat{M}}_t{Z}_t{\mathbb{E}}^{\mathbb{P}}\left[{\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}|{\mathcal{F}}_t\right]+Z_t{\mathbb{E}}^{\mathbb{P}}\left[\left({\sigma }_{t+1}{z}_{t+1}-\rho \right){\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}|{\mathcal{F}}_t\right].\end{array}$

For the first term on the right hand side, as ${\sigma }_{t+1}^2\in {\mathcal{F}}_t$ , and $z_{t+1}$ is independent of ${\mathcal{F}}_t$ under $\mathbb{P}$ , thus

${\mathbb{E}}^{\mathbb{P}}\left[{\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}|{\mathcal{F}}_t\right]={\mathbb{E}}^{\mathbb{P}}{\left[{\text{e}}^{\frac{\rho }{y}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{y^2}}\right]}_{y={\sigma }_{t+1}}=1.$

For the second term, we have

$\begin{array}{l}{\mathbb{E}}^{\mathbb{P}}\left[\left({\sigma }_{t+1}{z}_{t+1}-\rho \right){\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}|{\mathcal{F}}_t\right]\\ ={\mathbb{E}}^{\mathbb{P}}{\left[\left(y{z}_{t+1}-\rho \right){\text{e}}^{\frac{\rho }{y}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{y^2}}\right]}_{y={\sigma }_{t+1}}\\ ={\displaystyle {\int }_{-\infty }^{\infty }}\left({\sigma }_{t+1}x-\rho \right){\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}x-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}\frac{1}{\sqrt{2\pi }}{\text{e}}^{-\frac{x^2}{2}}\text{d}x\\ ={\displaystyle {\int }_{-\infty }^{\infty }}\left({\sigma }_{t+1}x-\rho \right)\frac{1}{\sqrt{2\pi }}{\text{e}}^{-\frac{1}{2{\sigma }_{t+1}^2}{\left({\sigma }_{t+1}x-\rho \right)}^2}\text{d}x\\ ={\displaystyle {\int }_{-\infty }^{\infty }}u\frac{1}{\sqrt{2\pi {\sigma }_{t+1}^2}}{\text{e}}^{-\frac{u^2}{2{\sigma }_{t+1}^2}}\text{d}u=0.\end{array}$

Thus, we have shown that

${\mathbb{E}}^{\mathbb{P}}\left[{\hat{M}}_{t+1}{Z}_{t+1}\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]={\hat{M}}_t{Z}_t\mathrm{,}\forall t\in {\mathbb{Z}}_{+}\mathrm{.}$

For any $k\in {\mathbb{Z}}_{+}$ ,

$\begin{array}{c}{\mathbb{E}}^{\mathbb{P}}\left[{\hat{M}}_{t+k}{Z}_{t+k}\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]={\mathbb{E}}^{\mathbb{P}}\left({\mathbb{E}}^{\mathbb{P}}\left[{\hat{M}}_{t+k}{Z}_{t+k}\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_{t+k-1}\right]\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right)\\ ={\mathbb{E}}^{\mathbb{P}}\left[{\hat{M}}_{t+k-1}{Z}_{t+k-1}\mathrm{<mo>\; |\; </mo>}{\mathcal{F}}_t\right]=\cdots ={\hat{M}}_t{Z}_t\mathrm{.}\end{array}$

This means that ${\hat{M}}_t{Z}_t$ is a $\mathbb{P}$ -martingale.

To prove that ${\hat{M}}_t$ is a $\mathbb{Q}$ -martingale, we apply Lemma A.1 and the results above, then get

${\mathbb{E}}^{\mathbb{Q}}\left[{\hat{M}}_{t+k}|{\mathcal{F}}_t\right]=\frac{1}{Z_t}{\mathbb{E}}^{\mathbb{P}}\left[{\hat{M}}_{t+k}{Z}_{t+k}|{\mathcal{F}}_t\right]={\hat{M}}_t.$ (40)

We concluded the first statement of this theorem.

2) $\exp \left\{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t}{\sigma }_i^2\right\}$ is also a $\mathbb{Q}$ -martingale. In fact, by Lemma A.1,

$\begin{array}{l}{\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{{\hat{M}}_{t+1}-\frac{1}{2}{\displaystyle {\sum }_{i=1}^{t+1}{\sigma }_i^2}}|{\mathcal{F}}_t\right]=\frac{1}{Z_t}{\mathbb{E}}^{\mathbb{P}}\left[{\text{e}}^{{\hat{M}}_{t+1}-\frac{1}{2}{\displaystyle {\sum }_{i=1}^{t+1}{\sigma }_i^2}}{Z}_{t+1}|{\mathcal{F}}_t\right]\\ ={\mathbb{E}}^{\mathbb{P}}\left[e^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}+{\sigma }_{t+1}{z}_{t+1}-\rho -\frac{1}{2}{\sigma }_{t+1}^2}{\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}|{\mathcal{F}}_t\right]\\ ={\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}{\mathbb{E}}^{\mathbb{P}}{\left[{\text{e}}^{y{z}_{t+1}-\rho -\frac{1}{2}{y}^2}{\text{e}}^{\frac{\rho }{y}{z}_{t+1}-\frac{1}{2}\frac{{\rho }^2}{y^2}}\right]}_{y={\sigma }_{t+1}}\\ ={\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}{\displaystyle {\int }_{-\infty }^{+\infty }}{\text{e}}^{{\sigma }_{t+1}x-\rho -\frac{1}{2}{\sigma }_{t+1}^2}{\text{e}}^{\frac{\rho }{{\sigma }_{t+1}}x-\frac{1}{2}\frac{{\rho }^2}{{\sigma }_{t+1}^2}}\frac{1}{\sqrt{2\pi }}{\text{e}}^{-\frac{x^2}{2}}\text{d}x\\ ={\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}{\text{e}}^{-\frac{1}{2}{\sigma }_{t+1}^2}{\displaystyle {\int }_{-\infty }^{+\infty }}{\text{e}}^{{\sigma }_{t+1}x-\rho }\frac{1}{\sqrt{2\pi }}{\text{e}}^{-\frac{1}{2{\sigma }_{t+1}^2}{\left({\sigma }_{t+1}x-\rho \right)}^2}\text{d}x\\ ={\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}{\text{e}}^{-\frac{1}{2}{\sigma }_{t+1}^2}{\displaystyle {\int }_{-\infty }^{+\infty }}{\text{e}}^u\frac{1}{\sqrt{2\pi {\sigma }_{t+1}^2}}{\text{e}}^{-\frac{u^2}{2{\sigma }_{t+1}^2}}\text{d}u={\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}.\end{array}$

For any $k\in {\mathbb{Z}}_{+}$ ,

$\begin{array}{c}{\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{{\hat{M}}_{t+k}-\frac{1}{2}{\displaystyle {\sum }_{i=1}^{t+k}{\sigma }_i^2}}|{\mathcal{F}}_t\right]={\mathbb{E}}^{\mathbb{Q}}\left({\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{{\hat{M}}_{t+k}-\frac{1}{2}{\displaystyle {\sum }_{i=1}^{t+k}{\sigma }_i^2}}|{\mathcal{F}}_{t+k-1}\right]|{\mathcal{F}}_t\right)\\ =\cdots ={\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}.\end{array}$

Thus we have proved that ${\text{e}}^{{\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}}$ is a $\mathbb{Q}$ -martingale.

Therefore,

$\begin{array}{l}{\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{{\displaystyle {\sum }_{i=t+1}^{t+k}{\sigma }_i{z}_i}+k\left(\mu -\log r\right)-\frac{1}{2}{\displaystyle {\sum }_{i=t+1}^{t+k}{\sigma }_i^2}}|{\mathcal{F}}_t\right]\\ ={\mathbb{E}}^{\mathbb{Q}}\left[{\text{e}}^{\left({\hat{M}}_{t+k}-{\displaystyle {\sum }_{i=1}^{t+k}\frac{1}{2}{\sigma }_i^2}\right)-\left({\hat{M}}_t-\frac{1}{2}{\displaystyle {\sum }_{i=1}^t{\sigma }_i^2}\right)}|{\mathcal{F}}_t\right]=1.\end{array}$

That is, we have constructed a risk-neutral measure $\mathbb{Q}$ satisfying Equation (34) for these GARCH-M models. Now we present the discrete-time option pricing formula for European call. Since $\left\{\frac{C_t}{r^t},t=0,1,2,\cdots ,T\right\}$ is a martingale under risk-neutral measure $\mathbb{Q}$ , then one pricing formula is

$C_t={\mathbb{E}}^{\mathbb{Q}}\left[{\left(S_t{\text{e}}^{{\displaystyle {\sum }_{i=t+1}^T{\sigma }_i{z}_i}-\rho \left(T-t\right)-\frac{1}{2}{\displaystyle {\sum }_{i=t+1}^T{\sigma }_i^2}}-K{r}^{-\left(T-t\right)}\right)}^{+}|{\mathcal{F}}_t\right].$ (41)

The option pricing formula can also be calculated under measure $\mathbb{P}$ , that is, by Lemma A.1,

$C_t={\mathbb{E}}^{\mathbb{P}}\left[{\left(S_t{\text{e}}^{{\displaystyle {\sum }_{i=t+1}^T{\sigma }_i{z}_i}-\rho \left(T-t\right)-\frac{1}{2}{\displaystyle {\sum }_{i=t+1}^T{\sigma }_i^2}}-K{r}^{-\left(T-t\right)}\right)}^{+}\times {\text{e}}^{{\displaystyle {\sum }_{i=t+1}^T\frac{\rho }{{\sigma }_i}{z}_i}-\frac{1}{2}{\displaystyle {\sum }_{i=t+1}^T\frac{{\rho }^2}{{\sigma }_i^2}}}|{\mathcal{F}}_t\right].$ (42)

Since

${\sigma }_{t+j}={\sigma }_{t+j}\left(z_{t+j-1}\mathrm{,}\cdots \mathrm{,}{z}_{t+1}\mathrm{,}{z}_t\mathrm{,}\cdots \right)\mathrm{,}$ (43)

they are correlated with $\left\{z_t\mathrm{,}t\in {\mathbb{Z}}_{+}\right\}$ . The price of option $C_t$ depends on attributes of the forecasting of the path followed by $\left\{{\sigma }_{t+1}^2\mathrm{,}{\sigma }_{t+2}^2\mathrm{,}\cdots \mathrm{,}{\sigma }_T^2\right\}$ .

In the following, we give the proof of Theorem 2.1.

Proof of Theorem 2.1. Actually, from the pricing formula in Equation (42), we have the following representation:

$C_t=S_t{N}_1-K{r}^{-\left(T-t\right)}{N}_2,$ (44)

where

$\begin{array}{l}{N}_1={\displaystyle \int }\cdots {\displaystyle {\int }_D}{\text{e}}^{{\displaystyle {\sum }_{j=1}^{T-t}{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j}-\rho \left(T-t\right)-\frac{1}{2}{\displaystyle {\sum }_{j=1}^{T-t}{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}\\ \times {\text{e}}^{{\displaystyle {\sum }_{j=1}^{T-t}\frac{\rho }{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}{x}_j}-\frac{1}{2}{\displaystyle {\sum }_{j=1}^{T-t}\frac{{\rho }^2}{{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-\frac{1}{2}{\displaystyle {\sum }_{j=1}^{T-t}{x}_j^2}}\text{d}{x}_1\cdots \text{d}{x}_{T-t},\end{array}$ (45)

$N_2={\displaystyle \int }\cdots {\displaystyle {\int }_D}{\text{e}}^{{\displaystyle {\sum }_{j=1}^{T-t}\frac{\rho }{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}{x}_j}-\frac{1}{2}{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}\frac{{\rho }^2}{{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-\frac{1}{2}{\displaystyle {\sum }_{j=1}^{T-t}{x}_j^2}}\text{d}{x}_1\cdots \text{d}{x}_{T-t},$ (46)

and

${\sigma }_{t+1}\left({\tilde{x}}_0\right)={\sigma }_{t+1}\left(z_t\mathrm{,}{z}_{t-1}\mathrm{,}\cdots \right)\mathrm{,}$

${\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)={\sigma }_{t+j}\left(x_{j-1},\cdots ,x_1;z_t,z_{t-1},\cdots \right),j=2,\cdots ,T-t,$

and the domain

$\begin{array}{l}D=\{\left(x_1,\cdots ,x_{T-t}\right)\in {ℝ}^{T-t}:{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-\rho \left(T-t\right)\\ >-\left(\log \frac{S_t}{K{r}^{-\left(T-t\right)}}-\frac{1}{2}{\displaystyle \underset{j=1}{\overset{T-t}{\sum }}}{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)\right)\},\end{array}$

with ${\tilde{x}}_j=\left(x_j\mathrm{,}\cdots \mathrm{,}{x}_1\right)$ and ${\tilde{x}}_0=\varnothing$ . By some calculations, we further have

$\begin{array}{c}{N}_1={\displaystyle \int }\cdots {\displaystyle {\int }_D}{\text{e}}^{{\displaystyle {\sum }_{j=1}^{T-t}\left({\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-\rho \right)}-\frac{1}{2}{\displaystyle {\sum }_{j=1}^{T-t}{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}\\ \times {\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-{\displaystyle {\sum }_{j=1}^{T-t}\frac{{\left({\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-\rho \right)}^2}{2{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}}\text{d}{x}_1\cdots \text{d}{x}_{T-t}\\ ={\displaystyle \int }\cdots {\displaystyle {\int }_D}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-{\displaystyle {\sum }_{j=1}^{T-t}\frac{{\left({\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-\rho -{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)\right)}^2}{2{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}}\text{d}{x}_1\cdots \text{d}{x}_{T-t}\end{array}$

Define

$u_j=x_j-\frac{\rho }{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)}-{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)\mathrm{,}$

we get

$N_1={\displaystyle \int }\cdots {\displaystyle {\int }_{D_1}}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-{\displaystyle {\sum }_{j=1}^{T-t}\frac{u_j^2}{2}}}\text{d}{u}_1\cdots \text{d}{u}_{T-t}={\Phi }_{T-t}\left(D_1\right),$

with the domain D₁ as stated in the theorem.

Similarly, by defining

$v_j=x_j-\frac{\rho }{{\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)},$

we have

$\begin{array}{c}{N}_2={\displaystyle \int }\cdots {\displaystyle {\int }_D}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-{\displaystyle {\sum }_{j=1}^{T-t}\frac{{({\sigma }_{t+j}\left({\tilde{x}}_{j-1}\right)x_j-\rho )}^2}{2{\sigma }_{t+j}^2\left({\tilde{x}}_{j-1}\right)}}}\text{d}{x}_1\cdots \text{d}{x}_{T-t}\\ ={\displaystyle \int }\cdots {\displaystyle {\int }_{D_2}}{\left(2\pi \right)}^{-\frac{T-t}{2}}{\text{e}}^{-{\displaystyle {\sum }_{j=1}^{T-t}\frac{v_j^2}{2}}}\text{d}{v}_1\cdots \text{d}{v}_{T-t}={\Phi }_{T-t}\left(D_2\right).\end{array}$

Thus, we concluded this theorem.

Appendix B. Proof of Theorem 2.2

As the above case in which the return rate $\mu$ is constant, we can prove

$M_t\triangleq M_0+{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\sigma }_i{z}_i,$ (47)

$Z_t\triangleq \exp \left\{{\displaystyle \underset{i=1}{\overset{t}{\sum }}}\frac{{\rho }_i}{{\sigma }_i}{z}_i-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{t}{\sum }}}\frac{{\rho }_i^2}{{\sigma }_i^2}\right\},{\rho }_i=-\left({\mu }_i-\log r_i\right),$ (48)

are martingales under measure $\mathbb{P}$ .

Let measure $\mathbb{Q}$ be

$\mathbb{Q}\left(A\right)\triangleq {\mathbb{E}}^{\mathbb{P}}\left({\chi }_A{Z}_t\right),\forall A\in {\mathcal{F}}_t,$ (49)

and

${\hat{M}}_t\triangleq M_t-{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\rho }_i=M_0+{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\sigma }_i{z}_i+{\displaystyle \underset{i=1}{\overset{t}{\sum }}}\left({\mu }_i-\log r_i\right).$ (50)

Proving line by line as the constant return rate case, we know that the stochastic processes ${\hat{M}}_t$ and

$\exp \left\{{\hat{M}}_t-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{t}{\sum }}}{\sigma }_i^2\right\}$

are ${\mathcal{F}}_t$ -martingales under measure $\mathbb{Q}$ .

For the proof of Theorem 2.2, it is almost the same line-by-line with Theorem 2.1. So we omit the proofs of these theorems here.

Appendix C. Parameter Estimation

We will describe one estimation algorithm for the models we have discussed in this paper.

C.1. Estimation Illustration

Consider the first model GARCH-M with constant return rate, for example,

$R_t=\mu -\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t,$ (51)

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2.$ (52)

There are four parameters $\mu \mathrm{,}\left({\alpha }_0\mathrm{,}{\alpha }_1\mathrm{,}{\beta }_1\right)$ to be estimated. The iterated estimation algorithm is suggested below.

Step 1: Provide a set of initial parameters ${\hat{\theta }}^{\left(0\right)}=\left({\hat{\alpha }}_0^{\left(0\right)}\mathrm{,}{\hat{\alpha }}_1^{\left(0\right)}\mathrm{,}{\hat{\beta }}_1^{\left(0\right)}\right)$ and ${\hat{\mu }}^{\left(0\right)}$ by fitting the data using a standard GARCH model with a nonzero mean, i.e.

$R_t=\mu +{\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2.$

Step 2: Compute the conditional volatility process ${\hat{\sigma }}_t^{\mathrm{2,}\left(i\right)}$ for $t=1,2,\cdots ,T$ by Equation (52) based on the parameters ${\hat{\theta }}^{\left(i\right)}=\left({\hat{\alpha }}_0^{\left(i\right)}\mathrm{,}{\hat{\alpha }}_1^{\left(i\right)}\mathrm{,}{\hat{\beta }}_1^{\left(i\right)}\right)$ obtained in the last step.

Step 3: Update ${\hat{\theta }}^{\left(i\right)}$ and ${\hat{\mu }}^{\left(i\right)}$ . That is, find ${\hat{\theta }}^{\left(i+1\right)}$ and ${\hat{\mu }}^{\left(i+1\right)}$ by the standard GARCH model with a nonzero mean

${\hat{R}}_t^{\left(i\right)}\triangleq R_t+\frac{1}{2}{\hat{\sigma }}_t^{2,\left(i\right)}=\mu +{\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2.$

Step 4: Repeat Steps 2 and 3 for a finite fixed number of iterations or until convergence.

More generally, for AR-GARCH-M models with linear risk premia, i.e.

$R_t=\mu +c{\sigma }_t^2+{\sigma }_t{z}_t,$ (53)

and

$R_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+c{\sigma }_t^2+{\sigma }_t{z}_t,$ (54)

with ${\sigma }_t^2$ , for example, a GARCH (1, 1) process, the iterated estimation algorithm is suggested below.

Step 1: Provide a set of initial parameters ${\hat{\theta }}^{\left(0\right)}=\left({\hat{\alpha }}_0^{\left(0\right)}\mathrm{,}{\hat{\alpha }}_1^{\left(0\right)}\mathrm{,}{\hat{\beta }}_1^{\left(0\right)}\right)$ by fitting the data using a standard AR(k)-GARCH model

$R_t={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+{\sigma }_t{z}_t.$

Step 2: Compute the conditional volatility process $\left\{{\hat{\sigma }}_t^{2,\left(i\right)},t=1,2,\cdots ,T\right\}$ by Equation (52) based on the parameters ${\hat{\theta }}^{\left(i\right)}=\left({\hat{\alpha }}_0^{\left(i\right)}\mathrm{,}{\hat{\alpha }}_1^{\left(i\right)}\mathrm{,}{\hat{\beta }}_1^{\left(i\right)}\right)$ obtained in the last step.

Step 3: Estimate ${\hat{\gamma }}^{\left(i\right)}=\left({\hat{\gamma }}_0^{\left(i\right)}\mathrm{,}{\hat{\gamma }}_1^{\left(i\right)}\mathrm{,}\cdots \mathrm{,}{\hat{\gamma }}_k^{\left(i\right)}\right)$ and ${\hat{c}}^{\left(i\right)}$ by the linear regression model

$\mathbb{E}\left(R_t|{\mathcal{F}}_{t-1}\right)={\gamma }_0+{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\gamma }_j{R}_{t-j}+c{\hat{\sigma }}_t^{2,\left(i\right)}.$ (55)

Step 4: Update ${\hat{\theta }}^{\left(i\right)}$ , and find ${\hat{\theta }}^{\left(i+1\right)}$ via the standard GARCH model

${\hat{R}}_t^{\left(i\right)}={\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2,$

where

${\hat{R}}_t^{\left(i\right)}\triangleq R_t-{\hat{\gamma }}_0^{\left(i\right)}-{\displaystyle \underset{j=1}{\overset{k}{\sum }}}{\hat{\gamma }}_j^{\left(i\right)}{R}_{t-j}-{\hat{c}}^{\left(i\right)}{\hat{\sigma }}_t^{2,\left(i\right)}.$

Step 5: Repeat Steps 2, 3 and 4 for a finite fixed number of iterations or until convergence.

For the estimation of the general model

$R_t=\mu \left(R\mathrm{,}{\sigma }_t^2\right)+{\sigma }_t{z}_t\mathrm{,}$ (56)

${\sigma }_t^2~\text{GARCHtypeprocess}\mathrm{,}$ (57)

one should first present a semiparametric estimation for $\mu \left(R\mathrm{,}{\sigma }_t^2\right)$ , and then use iterated estimation procedure as suggested above. As our main aim of this work is not to present the estimation problems for these models, we refer to [15] [16] and references therein for the details of estimation and convergence problems. For the exact likelihood inference of these models, see e.g. [17] , in which they suggested a Markov chain Monte Carlo algorithm to carry out the estimation probelms.

C.2. Empirical Example

The data set analyzed here consists of daily prices on the Standard and Poor’s 500 composite stock index from January 3 through December 29, 2017, for a total of $T=251$ observations. The daily prices are denoted $\left\{S_t,t=1,2,\cdots ,T\right\}$ and logarithm return as $R_t=\log \left(S_t\right)-\log \left(S_{t-1}\right)$ . The time t subscript refers to trading days.

Figure A1 plots the prices $S_t$ on the left and the logarithm return $R_t$ on the right for the whole year. Here we use three represented GARCH type models to fit these historical data.

Figure A1. Prices and log-return of S&P 500 Index during 2017.

The first model (say, Model 1) is

$R_t=\mu -\frac{1}{2}{\sigma }_t^2+{\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2.$

The iterated estimation procedure above leads to identification of the following model for $R_t=\log \left(S_t\right)-\log \left(S_{t-1}\right)$ :

$\begin{array}{l}{R}_t=6.6488\times {10}^{-4}-0.5{\sigma }_t^2+{\sigma }_t{z}_t\mathrm{,}\\ \left(2.6145\times {10}^{-4}\right)\end{array}$

$\begin{array}{l}{\sigma }_t^2=8.753\times {10}^{-7}+0.05{\sigma }_{t-1}^2{z}_{t-1}^2+0.9{\sigma }_{t-1}^2\mathrm{.}\\ \left(8.9816\times {10}^{-7}\right)\left(0.0176\right)\left(0.0165\right)\end{array}$

The numbers in the parenthesis are standard errors.

The second model (say, Model 2) is

$R_t=\mu +c{\sigma }_t^2+{\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2,$

which is specified to

$\begin{array}{l}{R}_t=5.8333\times {10}^{-4}+5.5693{\sigma }_t^2+{\sigma }_t{z}_t\mathrm{,}\\ \left(9.2205\times {10}^{-5}\right)\left(5.4418\right)\end{array}$

$\begin{array}{l}{\sigma }_t^2=8.7397\times {10}^{-7}+0.05{\sigma }_{t-1}^2{z}_{t-1}^2+0.9{\sigma }_{t-1}^2\mathrm{.}\\ \left(8.9682\times {10}^{-7}\right)\left(0.0175\right)\left(0.0162\right)\end{array}$

The third model (say, Model 3) is

$R_t={\gamma }_0+{\gamma }_1{R}_{t-1}+c{\sigma }_t^2+{\sigma }_t{z}_t,$

${\sigma }_t^2={\alpha }_0+{\alpha }_1{\sigma }_{t-1}^2{z}_{t-1}^2+{\beta }_1{\sigma }_{t-1}^2.$

The specification for this model is

$\begin{array}{l}{R}_t=8.0621\times {10}^{-4}-0.1344R_{t-1}-2.216{\sigma }_t^2+{\sigma }_t{z}_t\mathrm{,}\\ \left(0.0015\right)\left(0.0632\right)\left(31.7553\right)\end{array}$

$\begin{array}{l}{\sigma }_t^2=8.5598\times {10}^{-7}+0.05{\sigma }_{t-1}^2{z}_{t-1}^2+0.9{\sigma }_{t-1}^2.\\ \left(8.7837\times {10}^{-7}\right)\left(0.0174\right)\left(0.0158\right)\end{array}$

These results are standard in the literature. We show them as representative examples. For these three models, the parameters are robustly obtained by the iterated estimation procedure. Another fact is that the volatility processes in these three models are very stable.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Engle, R.F. (2002) New Frontiers for ARCH Models. Journal of Applied Econometrics, 17, 425-446. https://doi.org/10.1002/jae.683
[2]	Black, F., and Scholes, M. (1973) The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-659. https://doi.org/10.1086/260062
[3]	Merton, R.C. (1973) Theory of Rational Option Pricing. Bell Journal of Economics and Management Science, 4, 141-183. https://doi.org/10.2307/3003143
[4]	Cox, J., Ross, S. and Rubinstein, M. (1979) Option Pricing: A Simplified Approach. Journal of Financial Economics, 7, 229-263. https://doi.org/10.1016/0304-405X(79)90015-1
[5]	Engle, R.F. (1982) Autoregressive Conditional Heteroskedasticity with Estimates of the Variance of U.K. Inflation. Econometrica, 50, 987-1008. https://doi.org/10.2307/1912773
[6]	Bollerslev, T. (1986) Generalized Autoregressive Conditional Heteroskedasticity. Journal of Econometrics, 31, 307-327. https://doi.org/10.1016/0304-4076(86)90063-1
[7]	Nelson, D.B. (1991) Conditional Heteroskedasticity in Asset Returns: A New Approach. Econometrica, 59, 347-370. https://doi.org/10.2307/2938260
[8]	Francq, C. and Zakoian, J.-M. (2010) GARCH Models: Structure, Statistical Inference and Financial Applications. Wiley, Hoboken. https://doi.org/10.1002/9780470670057
[9]	Engle, R.F., Lilien, D.M. and Robins, R.P. (1987) Estimating the Time Varying Risk Premia in the Term Structure: The ARCH-M Model. Econometrica, 55, 391-407. https://doi.org/10.2307/1913242
[10]	Hull, J. and White, A. (1987) The Pricing of Options on Assets with Stochastic Volatilities. Journal of Finance, 42, 281-300. https://doi.org/10.1111/j.1540-6261.1987.tb02568.x
[11]	Bollerslev, T. and Mikkelsen, H. (1996) Modeling and Pricing Long Memory in Stock Market Volatility. Journal of Econometrics, 73, 151-184. https://doi.org/10.1016/0304-4076(95)01736-4
[12]	Duan, J.-C. (1995) The GARCH Option Pricing Model. Mathematical Finance, 5, 13-32. https://doi.org/10.1111/j.1467-9965.1995.tb00099.x
[13]	Zumbach, G. and Fernández, L. (2013) Fast and Realistic European ARCH Option Pricing and Hedging. Quantitative Finance, 13, 713-728. https://doi.org/10.1080/14697688.2012.750009
[14]	Zumbach, G. and Fernández, L. (2014) Option Pricing with Realistic ARCH Processes. Quantitative Finance, 14, 143-170. https://doi.org/10.1080/14697688.2013.816437
[15]	Christensen, B., Dahl, C. and Iglesias, E. (2012) Semiparametric Inference in a GARCH-in-Mean Model. Journal of Econometrics, 167, 458-472. https://doi.org/10.1016/j.jeconom.2011.09.028
[16]	Conrad, C. and Mammen, E. (2016) Asymptotics for Parametric GARCH-in-Mean Models. Journal of Econometrics, 194, 319-329. https://doi.org/10.1016/j.jeconom.2016.05.010
[17]	Fiorentini, G., Sentana, E. and Shephard, N. (2004) Likelihood-Based Estimation of Latent Generalized ARCH Structures. Econometrica, 72, 1481-1517. https://doi.org/10.1111/j.1468-0262.2004.00541.x
[18]	Aït-Sahalia Y. and Jacod, J. (2014) High-Frequency Financial Econometrics. Princeton University Press, Princeton. https://doi.org/10.23943/princeton/9780691161433.001.0001
[19]	Noureldin, D., Shephard, N. and Sheppard, K. (2012) Multivariate High-Frequency-Based Volatility (HEAVY) Models. Journal of Applied Econometrics, 27, 907-933. https://doi.org/10.1002/jae.1260
[20]	Shephard, N. and Sheppard, K. (2010) Realising the Future: Forecasting with High-Frequency-Based Volatility (HEAVY) Models. Journal of Applied Econometrics, 25, 197-231. https://doi.org/10.1002/jae.1158