GARCH Option Pricing Models and the Variance Risk Premium

In this paper, we modify Duan’s (1995) local risk-neutral valuation relationship (mLRNVR) for the GARCH option-pricing models. In our mLRNVR, the conditional variances under two measures are designed to be different and the variance process is more persistent in the risk-neutral measure than in the physical one, so that one is able to capture the variance risk premium. Empirical estimation exercises show that the GARCH option-pricing models under our mLRNVR are able to price the SPX one-month variance swap rate, i.e., the CBOE Volatility Index (VIX) accurately. Our research suggests that one should use our mLRNVR when pricing options with GARCH models.


Introduction
In this paper, we modify the local risk-neutral valuation relationship (mLRNVR) in the GARCH option-pricing models. The GARCH option-pricing model was first introduced by Duan (1995) with a locally risk-neutral valuation relationship (LRNVR), in which the conditional variances and model parameters remained the same under the physical measure and the risk-neutral measure. Since then, Duan's LRNVR has been widely used by finance researchers and practitioners in pricing options under the GARCH framework. However, as Barone-Adesi et al. (2008) pointed out, empirical evidence showed that the restriction in Duan's LRNVR led to rather poor pricing and hedging performances. Hao and Zhang (2013) showed explicitly that the GARCH option-pricing models under Duan's LRNVR under-priced S&P 500 (SPX) one-month variance swap rate, i.e., the Chicago Board Options Exchange (CBOE) Volatility Index (VIX) by about 10%. We propose a new mLRNVR in GARCH option models in order to resolve the issue of under-pricing by the LRNVR.
Since the seminal work of Bollerslev (1986); Engle (1982), the family of GARCH volatility models has been widely used in empirical asset pricing and financial risk management partly because of the likelihood function of asset returns in the GARCH models could often be expressed in a closed-form in terms of observed data. One can estimate the model parameters by using the maximum likelihood estimation (MLE) method, which is often a challenging task for most of the stochastic volatility models. Motivated by the success of GARCH models in fitting asset returns, Duan (1995) pioneered in applying GARCH models for the SPX option pricing by proposing the LRNVR. Duan's LRNVR has been followed by many papers. Ritchken and Trevor (1999) developed an efficient lattice algorithm to price European and American options under the GARCH processes. Chernov and Ghysels (2000) estimated physical and risk-neutral probabilities by using Heston (1993) and GARCH option models. Heston and Nandi (2000) developed a closed-form option valuation formula for a spot asset whose variance follows a GARCH(p, q) process that can be correlated with the returns of the spot asset. Christoffersen and Jacobs (2004) compared a variety of GARCH option-pricing models using option prices and asset returns. Christoffersen et al. (2008) extended the Heston-Nandi model to incorporate long-run and short-run volatility components for the valuation of European options. An interesting survey on the theoretical results and empirical evidence of GARCH option-pricing could be found in Christoffersen et al. (2013b). The evaluation of European options in the GARCH models with Levy jumps was discussed in Ornthanalai (2014). Kanniainen et al. (2014) evaluated SPX options using three variant GARCH models with VIX data. Their empirical evidence showed that a joint maximum likelihood estimation using SPX returns and VIX improved the performance of pricing SPX options compared with traditional MLE using the returns data only. Wang et al. (2017) priced the VIX futures with the Heston-Nandi GARCH model. All of these papers used Duan's LRNVR in changing probability measures from a physical one to a risk-neutral one. However, there exists a problem in pricing options under the GARCH models with the LRNVR.
Empirical evidence by Chernov and Ghysels (2000); Christoffersen and Jacobs (2004) showed that the GARCH option-pricing model under the LRNVR had poor pricing and hedging performances, but they did not realize that the problem came from the LRNVR. Barone-Adesi et al. (2008) were the first to explicitly point out that the poor performance came from the restriction required by the LRNVR. They proposed a new method for pricing options based on GARCH models with filtered historical innovations, in which they did not specify directly the change of measures, and proposed approximating them by calibrating a new set of risk-neutral GARCH parameters directly on market option prices. Christoffersen et al. (2013a) developed a GARCH option model with a new pricing kernel allowing for a variance premium. They presented an explicit relationship between physical and risk-neutral conditional variances and GARCH parameters, in which risk aversion and the market price of variance risk entered into their analytical formulas. It is interesting to note that the persistence parameter remained unchanged in their change of probability measures. Hao and Zhang (2013) developed a model for the CBOE VIX by using a GARCH option-pricing model under the LRNVR. They found that the model-implied VIX was about 10% lower than the market VIX due to the lack of a variance risk premium. Given that the VIX is a proxy of the SPX one-month variance swap rate, one may conclude that GARCH option-pricing models under LRNVR under-price variance swap and hence at-the-money options by about 10%. Despite these developments of literature pointing out the problem of LRNVR, later studies, e.g., Ornthanalai (2014); Kanniainen et al. (2014); Wang et al. (2017), still used it in pricing SPX options and other derivatives. They might have realized the problem with the LRNVR, but lacked a better one to resolve the issue. Developing a new risk-neutral valuation relationship in GARCH option-pricing models becomes an urgent task. That is the focus in this paper. Parallel to our work, Huang et al. (2017Huang et al. ( 2019) used a stochastic discount factor for risk-neutralization for option pricing and VIX futures pricing, respectively. The stochastic discount factor risk-neutralization is a generalization of the LRNVR. In contrast to our modification, the persistent parameter in the models remained unchanged under the stochastic discount factor risk-neutralization.
In this paper, we propose a modified local risk-neutral valuation relationship, referred to as the mLRNVR. In our mLRNVR, conditional variances under two measures are designed to be different and the GARCH process is more persistent in the risk-neutral measure than in the physical one, so that we are able to capture the negative variance risk premium. We display the improvement of the mLRNVR compared with the LRNVR using theoretical and empirical evidences. Specifically, we find the theoretical VIX squared value as the conditional risk-neutral expectation of the arithmetic mean variance over the next 21 trading days under the mLRNVR. The GARCH implied VIX formulas are derived using the features of square-root stochastic autoregressive volatility (SR-SARV) models. We apply several calibration methods to estimate the model parameters using various sets of time series data, and compare the theoretical formula performances with the market data. Various combinations of time series of the daily closing price of the S&P 500 index and the CBOE VIX are used to find the maximum likelihood estimation of the GARCH models. The corresponding implied VIX time series are then calculated from the calibrated model. Similar to the empirical evidence in Hao and Zhang (2013) and Wang et al. (2017), when only S&P 500 returns are used for estimation, the GARCH implied VIX is consistently and significantly lower than the CBOE VIX. When the CBOE VIX is used for estimation, the implied VIX fits the statistical properties of the CBOE VIX and matches the CBOE VIX data better. The numerical results provide evidence that the GARCH option-pricing under the mLRNVR is more suitable to price variance swap. In the case of GARCH(1,1), we compare the diffusion limit of the GARCH process under the physical measure and the mLRNVR risk-neutral measure to show that the variance premium is captured in the risk-neutral dynamics.
This paper makes at least two contributions. First, by using a newly proposed mLRNVR, we have resolved the issue of under-pricing in GARCH option-pricing models under the LRNVR. Our research suggests that one should use our mLRNVR when pricing options with GARCH models. Second, we develop a model for the VIX in a GARCH framework by using the mLRNVR. Our model is able to capture negative variance risk premiums, and hence should be used in pricing VIX derivatives, including VIX futures and VIX options.
The article is structured as follows. Section 2 proposes a new risk-neutral valuation relationship for the GARCH(1,1). In Section 3, we derive a theoretical VIX formula for the GARCH(1,1) model under the mLRNVR, and extend the derivation idea to a broad class of GARCH models, which include GARCH, TGARCH, AGARCH and EGARCH models. In Section 4, we calibrate these GARCH models using various combinations of time series of the S&P 500 index and the CBOE VIX. Section 5 compares the CBOE VIX with the GARCH implied VIX obtained from the calibrated GARCH models.
In Section 6, we analyze the diffusion limit of the GARCH process under the risk-neutral measure in mLRNVR to demonstrate that the risk-neutral dynamics captures the variance risk premium. We then conclude the findings in Section 7.

GARCH Model Specification
In this paper, we consider the asset price as a discrete-time stochastic process and denote the asset price at time t as X t . It was proposed in Duan (1995) that the returns of the asset follows a conditional lognormal distribution under the physical measure P as where r is the one-period risk-free interest rate, λ 1 is the asset risk premium, and t follows a GARCH(p, q) process introduced in Bollerslev (1986) with mean zero and conditional variance h t where φ t is the information set of up to and including time t; α 0 ≥ 0, α i ≥ 0 for i = 1, 2, . . . , q and β j ≥ 0 for j = 1, 2, . . . , p. We focus on the GARCH(1,1) case, so the Equation (2) simplifies to In order to accommodate the heteroskedasticity of the asset returns process in (1), Duan (1995) introduced the LRNVR under which the expected return should be the risk-free rate and the one-step-ahead variance should be the same in both measures. Specifically, under the LRNVR, the dynamics of asset returns in the risk-neutral pricing measure Q has the form It is well documented in Bollerslev et al. (2009);Carr and Wu (2009) that variance risk premiums come from either its correlation with the return risk and return risk premium, or a separate premium of the independent variance variation. It is also established in Carr and Wu (2009) that the majority of the market variance risk premium is generated by an independent variance risk factor. Motivated by the empirical findings that the market variance risk premiums is closely related to the variance variation, we propose an alternative risk-neutral valuation relationship to the LRNVR which can capture the variance risk premium as follows.

Definition 1.
A pricing measure Q is said to satisfy the modified local risk-neutral valuation relationship (mLRNVR) if the dynamics of asset returns in the risk-neutral pricing measure Q under the mLRNVR has the following form Remark 1. In the classic continuous-time stochastic volatility models, e.g., the Heston model (Heston 1993), the instantaneous variances are the same under the physical and risk-neutral measures. However, in the discrete-time GARCH models, the conditional variances are defined in a finite-time interval which is usually one day in practice. The conditional variances should be designed to be different under the physical and risk-neutral measures, as in the mLRNVR, to incorporate the variance risk premium during the finite-time interval.
The proposed GARCH(1,1) process (5) is different to the one derived by Duan (1995) under the LRNVR. Under the mLRNVR the persistence parameter is designed to be different in the P and Q measures, whereas under the LRNVR the persistence parameter is the same in the P and Q measures. Specifically, for the dynamics of risk-neutral measure Q under the mLRNVR, the persistence parameter of conditional variance is β * 1 = β 1 − √ 2α 1 λ 2 , where λ 2 represents the variance risk premium of the asset. Note that there is a negative sign in front of the variance risk premium λ 2 which is estimated to be negative from empirical studies Bollerslev et al. (2009);Carr and Wu (2009). Hence, we expect the variance process is more persistent in the risk-neutral measure than that in the physical one. The motivation for the inclusion of the variance risk premium is discussed in Hao and Zhang (2013), where it was shown that there is no risk adjustment for the variance risk of the process in Duan (1995) from the physical measure to the risk-neutral measure under the LRNVR. It was also discussed in Barone-Adesi et al. (2008) and Christoffersen et al. (2013a) that the restriction of conditional volatility of historical and risk-neutral pricing distributions with the same model parameters leads to poor calibration results in the empirical studies (cf. Chernov and Ghysels 2000;Christoffersen et al. 2006;Hao and Zhang 2013). Therefore, it was suggested that the parameters of volatility dynamics of historical and risk-neutral pricing returns might be different in Barone-Adesi et al. (2008). We adapt the idea by modifying the persistence parameter in Q to incorporate the variance risk premium in the model. 1 The theoretical justification of the modification is further discussed in Section 6.

VIX Formulas of the GARCH Models
The CBOE introduced the VIX index in 1993. The VIX was calculated from the implied volatilities of the eight near-the-money, nearby, and second nearby S&P 100 index options based on the methodology of Whaley (Whaley 1993). The VIX was a proxy of the implied volatility of 30 calendar days at-the-money (ATM) options. In 2003, the CBOE used another theory proposed in Carr and Madan (1998) and Demeterfi et al. (1999) to design a new methodology to compute the CBOE's VIX. The new VIX is based on the prices of a portfolio of 30 calendar days out-of-the-money (OTM) S&P 500 index call and put options. The square of the new VIX represents the S&P 500 30-day variance swap rate. The old VIX has been renamed the VXO. We use VIX Mkt to denote the CBOE VIX which is computed using market option prices.
The CBOE VIX index reflects investors' expectation of the volatility of the S&P 500 in the next 30 calendar days or 21 trading days. Following Hao and Zhang (2013), we have calculated the VIX index implied by the GARCH models as the mean value of the expected variance in the n sub-periods of the next 21 trading days, that is where τ 0 represents 21 trading days and VIX Imp t stands for the model implied VIX index. In particular, we use the daily closing value data, so it implies τ 0 = n, and where the term V t = 1 252 VIX Imp t 100 2 is defined as a function of VIX Imp t to measure the expected daily variance of the S&P 500. The conditional mean of the future variance can be calculated in a broad class of GARCH models, as discussed in Hao and Zhang (2013) and Wang et al. (2017). We derive the implied VIX from the model (5) under Q by first rewriting the error terms of the process using the standard normal distribution as where t is the standard normal random variable, conditional on the information set up to and including time t − 1 under Q.

1
Our purpose is trying to propose a simple modification to the LRNVR in order to incorporate the variance risk premium under the risk-neutral measure in the GARCH models. It is noted that in Huang et al. (2019), with an additional source of uncertainty in the conditional variance under the physical measure, the Radon-Nikodym derivative was derived and the variance risk premium was incorporated in the risk-neutral measure. However, these two approaches differ, since our mLRNVR method only changes the persistence parameter without introducing an extra uncertainty.
One can rewrite the GARCH(1, 1) process (8) as a special case of the square-root stochastic autoregressive volatility (SR-SARV(1)) models introduced in Meddahi and Renault (2004) with the following form It was shown in Hao and Zhang (2013) that if the S&P 500 returns follow a SR-SARV(p) process under the risk-neutral measure, then the implied daily variance V t at time t is affine in the conditional variance h t+1 . Following similar ideas, we can obtain the long term variance as Then the conditional expectation of the variance after two periods can be obtained via the long-run variance: So the conditional expectation of the variance after n periods is given by Therefore, we can represent the expected daily variance as an affine function of h t+1 . Apart from the GARCH(1,1) model discussed above, we also consider the threshold GARCH(1,1) (TGARCH) model introduced in Glosten et al. (1993), the non-linear asymmetric GARCH(1,1) (AGARCH) model proposed in Engle and Ng (1993) and the exponential GARCH(1,1) (EGARCH) model by Nelson (Nelson 1991). The forms of the models in the physical measure P and in the risk-neutral measure Q under the mLRNVR are as follows: TGARCH(1,1) mLRNVR: EGARCH(1,1) Physical measure: mLRNVR: As shown in Hao and Zhang (2013), these widely used GARCH models are special cases of SR-SARV(p) models, and following similar derivation process as the GARCH(1,1) model, we can obtain the implied VIX formula for different GARCH models analogous to the ones obtained in Hao and Zhang (2013). For the convenience of readers, we list the implied VIX formula as follows: TGARCH(1,1) where Note that N(·) denotes the cumulative function of the normal distribution. AGARCH(1,1) where where

Data and Estimation
It was shown in Hao and Zhang (2013) that under the LRNVR, the GARCH implied VIX does not fit the market data of the CBOE VIX very well. The model was analyzed to show that the reason may be that the variance risk premium and the volatility risk price were not present in the diffusion limit of the GARCH models under the LRNVR. In the modified GARCH processes, we include the variance risk premium in the models under the mLRNVR. It is also of interest to see whether the implied VIX in the modified GARCH models fit the CBOE VIX market values better. In this section, we will investigate this question by estimating the parameters in the modified GARCH models and calculating the corresponding GARCH implied VIX time series for comparison with the CBOE VIX.
The time-series data we use for the GARCH models calibration are the closing values of the S&P 500 and the CBOE VIX ranging from 2nd January 1990 to 30th June 2017. For the daily risk-free interest rate, we use the three-month Treasury bill secondary market rate from the U.S. Federal Reserve website.
There are different methods to calibrate the models using market data. We will use the common maximum likelihood approach to estimate the parameters of the models. We can use only the S&P 500 returns data to obtain a maximum likelihood estimation of the GARCH processes under the physical measure P and fix the variance risk premium parameter λ 2 = 0, since λ 2 is not included in the GARCH models under the P measure. For the S&P 500 returns data only, the log-likelihood function ln L R for the GARCH models is given by where the conditional variance h t is updated by corresponding processes using different forms of GARCH models and X t is defined in equation (1) as the asset price. For the maximum likelihood estimation, the conditional variance for the first period is set as the variance of S&P 500 returns over the whole sample period. The stationary conditions for the GARCH processes under the physical and the risk-neutral measures are different, with the latter having stricter constraints on the parameters. Thus, we find the estimation of the parameters in the GARCH models by maximizing the corresponding log-likelihood function subject to the stationary conditions under the risk-neutral measures. We may also calibrate the GARCH models by matching the model implied VIX to the market value of the CBOE VIX, since the CBOE VIX series may contain additional information about the underlying S&P 500 returns process. To utilize both time series, we follow the assumption in Hao and Zhang (2013) that the pricing differences between the CBOE VIX and the implied VIX on a daily basis come from a normal distribution where s 2 is estimated using the sample variance of pricing differenceŝ 2 = var(VIX Mkt − VIX Imp ). Under the above assumption, the log-likelihood function corresponding to the CBOE VIX data is Alternative to the method discussed in Hao and Zhang (2013), it was suggested in Kanniainen et al. (2014) to use the following model to describe autoregressive disturbances: where µ t = VIX Mkt t − VIX Imp t and e t ∼ NID(0, σ 2 ). Based on the autoregressive disturbance process, we can estimate the variance of pricing differenceŝ 2 in terms of the autoregressive disturbance correlation ρ and the varianceσ 2 obtained from the sample. The formula can be specified asŝ 2 =σ 2 1−ρ 2 from the relation Var(µ t ) = Var(ρµ t−1 + e t ) =⇒ s 2 = ρ 2 s 2 + σ 2 .
The log-likelihood function corresponding to the CBOE VIX data is We will compare the performance of this autoregressive disturbance process to the mLRNVR process in the numerical simulation section. Apart from using the S&P 500 returns data and CBOE VIX data for calibration of the GARCH models separately, we also combine both time series to find a joint maximum likelihood estimation of the models by maximizing the joint log-likelihood function

Numerical Results
In this section, we compare the estimated parameters from different data used for calibration. 2 In particular, the output tables display the maximum likelihood estimates and the standard errors of the GARCH models. The values of the log-likelihood functions (20), (22), (25), (26) are also displayed in the tables. Although the contributions from the S&P 500 returns and the CBOE VIX as well as the joint likelihood values are reported, we maximize the function ln L R when only S&P 500 returns are used, the function ln L V when only CBOE VIX data are used and the function ln L T when both time series are used. All the output tables show the comparison results among the LRNVR process, the mLRNVR process and the LRNVR process with an autoregressive coefficient (referred to as the generalized LRNVR process in the following). We can observe that the mLRNVR process has the largest maximum likelihood values compared with the LRNVR and the generalized LRNVR processes in all GARCH models. The maximum likelihood values of the generalized LRNVR processes are generally better than those of the LRNVR processes. It is as expected since the generalized LRNVR process has an extra autoregressive coefficient ρ compared with the LRNVR process. However, the extra autoregressive coefficient ρ is not statistically significant in all GARCH models. So the generalized LRNVR process is not significantly better than the LRNVR process, and we will not discuss the numerical results of the generalized LRNVR process in details and focus on the comparison between the mLRNVR and LRNVR processes.
From the output in Table 1, we can see that the equity risk premium λ 1 increases significantly from 0.0886 (returns data used) to 0.2134 (both data used) and 0.2253 (VIX data used) in the GARCH(1,1) models when the CBOE VIX data are used for calibration. The variance risk premium λ 2 is negative and significantly different from zero as −0.3670 (both sets of data used) and −0.3514 (VIX data used). The persistence of the conditional variance β 1 increases slightly from 0.8543 (returns data used) to 2 We have also tried the out-of-sample prediction. Note that the relation between VIX 2 t and the latent process h t is linear according to Equation (13). Predicting VIX 2 t is essentially forecasting the latent process h t . This task seems to be not straightforward in the GARCH models, hence is left for further research. 0.9251 (both sets of data used) and 0.9286 (VIX data used). There is a sizable decrease of the parameter value α 1 from 0.1256 (returns data used) to 0.0474 (both sets of data used) and 0.0456 (VIX data used). Comparing the maximum likelihood result of the model under the mLRNVR and the results under the LRNVR, we can see that the maximum likelihood value increases significantly from 54,697 to 55,921 (both sets of data used) and from 33,424 to 33,662 (VIX data used).
Similar numerical results are also observed in the other types of GARCH models as displayed in Tables 2-4. Specifically, Table 2 shows that the equity risk premium λ 1 increases significantly from 0.0131 (returns data used) to 0.1160 (both sets of data used) and 0.0889 (VIX data used) in the TGARCH(1,1) model when the CBOE VIX data are used for calibration. The variance risk premium λ 2 is negative and significantly different from zero as −0.4112 (both sets of data used) and −0.3978 (VIX data used). The persistence of conditional variance, β 1 increases significantly from 0.8338 (returns data used) to 0.9561 (both sets of data used) and 0.9553 (VIX data used). There is a decrease of the parameter value α 1 from 0.0256 (returns data used) to 0.0091 (both sets of data used) and 0.0060 (VIX data used). Comparing the maximum likelihood result of the TGARCH(1,1) model under the mLRNVR and the results under the LRNVR, we can see that the maximum likelihood value increases significantly from 55,455 to 56,282 (both sets of data used) and from 33,468 to 33,795 (VIX data used). Table 3 shows the calibration results of the AGARCH(1,1) model using both returns and VIX data. If using VIX data only, it is not easy to distinguish the parameters θ and λ 1 . Therefore, the numerical results using VIX data are not displayed in the table. From Table 3 we observe that the equity risk premium λ 1 increases significantly from 0.0255 (returns data used) to 0.1158 (both sets of data used) in the AGARCH(1,1) model when the CBOE VIX data and returns are used for calibration. The variance risk premium λ 2 is negative and significantly different from zero as −0.3125 (both sets of data used). The persistence of conditional variance β 1 increases from 0.8810 (returns data used) to 0.9316 (both sets of data used). There is a big decrease in the parameter value α 1 from 0.0841 (returns data used) to 0.0380 (both sets of data used). Comparing the maximum likelihood result of the AGARCH(1,1) model under the mLRNVR and the results under the LRNVR, we can see that the maximum likelihood value increases significantly from 55,483 to 56,333 (both sets of data used). Table 4 shows that in the EGARCH(1,1) model the variance risk premium λ 2 is negative as −0.0567 (both sets of data used) and −0.0483 (VIX data used), both significantly different than zero. The persistence of conditional variance, β 1 increases slightly from 0.9792 (returns data used) to 0.9906 (both sets of data used) and 0.9895 (VIX data used). Comparing the maximum likelihood result of the EGARCH(1,1) model under the mLRNVR and the results under the LRNVR, we can see that the maximum likelihood value increases significantly from 56,399 to 57,105 (both sets of data used) and from 33,774 to 34,303 (VIX data used).

Return Only VIX Only
Return & VIX LRNVR mLRNVR LRNVR mLRNVR Duan (1995) KLY ( Table 5. The table displays the related statistics between the model implied VIX and the CBOE VIX for the GARCH models during the period from 2 January 1990 to 30 June 2017. We use the sample of 6928 daily data for the model estimation and comparison. The error is computed as the difference between CBOE VIX and the implied VIX. The mean error (ME) represents the mean daily difference between the implied VIX and the CBOE VIX. The standard error (Std.Err.) represents the standard deviation of the error. The mean absolute error (MAE) calculates the mean daily absolute difference between the implied VIX and the CBOE VIX. The mean squared error (MSE) computes the mean daily squared difference between the implied VIX and the CBOE VIX. The root mean squared error (RMSE) computes the square root of the mean squared error. The P-value is for the null hypothesis that the implied VIX and the CBOE VIX have the same mean values. Violation of one-sigma band stands for the probability that the implied VIX lies out of the one-standard-deviation band of the CBOE VIX. The correlation coefficient (Corr. Coef.) computes the linear correlation between the implied VIX and the CBOE VIX. Autocorrelation coefficients with 1, 10, and 30 days lag and higher moments of implied VIX are also reported. Note that the calibration results using returns data only are produced by using GARCH models under the local risk-neutral valuation relationship (LRNVR). The results using VIX and both sets of data are produced by using GARCH models with the modified local risk-neutral valuation relationship (mLRNVR). After obtaining the estimates of the parameters in the models, we can then calculate the conditional variance h t and compute the corresponding GARCH implied VIX. Figure 1 shows the time series of the CBOE VIX and the implied VIX of the four GARCH models estimated using returns only. Figure 2 shows the time series of the CBOE VIX and the implied VIX of the GARCH(1,1) model estimated using VIX data only. Figure 3 shows the time series of the CBOE VIX and the implied VIX of the GARCH(1,1) model estimated using both returns and VIX. Similar comparison plots are obtained for other GARCH models. Specifically, the time series of the CBOE VIX and the implied VIX of the TGARCH(1,1) model estimated with VIX data only are displayed in Figure 4. The time series of the CBOE VIX and the implied VIX of the TGARCH(1,1) and AGARCH (1,1) model estimated with both returns and VIX data are shown in Figures 5 and 6, respectively. For the EGARCH(1,1) model, Figure 7 shows the comparison between the CBOE VIX and model implied VIX with VIX data only, and Figure 8 displays the comparison between the CBOE VIX and model implied VIX with both returns and VIX data. From the list of graphs, we observe that the model implied VIX fits the CBOE VIX better under the mLRNVR compared with the LRNVR in general. In particular, the ratios of the implied VIX values to the CBOE VIX values are closer to 1 under the mLRNVR compared to those under the LRNVR, as shown in Figures 3, 5, 6 and 8. The direct comparisons of performance under LRNVR and mLRNVR for the GARCH, TGARCH, AGARCH and EGARCH models are listed in Tables 1 -4, respectively. 1990 1993 1996 1999 2002 2005 2008

Theoretical Justification
Duan studied the bivariate diffusion limit of the GARCH(1,1) model as the length of the time period tends towards zero in Duan (1996Duan ( , 1997. Applying Duan's arguments, one can show that the limiting bivariate diffusion process of the approximating GARCH(1, 1) process under the physical measure P is given by where the persistence parameter of conditional variance is defined as β * 1 = β 1 − √ 2α 1 λ 2 under the mLRNVR. The terms dW 1t and dW 2t are independent standard Brownian motions under the physical measure P. The limiting bivariate diffusion under the risk-neutral measure Q is a re-parameterization of Hull and White (1987) bivariate diffusion model as follows: d ln X t = r − 1 2 h t dt + h t dZ 1t dh t = (α 0 + (α 1 + β * 1 − 1)h t )dt + √ 2α 1 h t dZ 2t , where dZ 1t = dW 1t + λ 1 dt and dZ 2t = dW 2t + λ 2 dt are independent standard Brownian motions under the mLRNVR Q. Both equity risk premium λ 1 and variance risk premium λ 2 are present in the model under the risk-neutral measure Q. The discrete-time GARCH(1, 1) process (5) corresponds to the limiting diffusion process under the risk-neutral measure Q.

Conclusions
In this paper, we follow the GARCH option-pricing framework of Duan (1995) and propose a modified local risk-neutral valuation relationship. The new risk-neutral valuation is referred to as the mLRNVR. The advantage of the mLRNVR compared with the LRNVR commonly used in the literature (Duan (1995); Hao and Zhang (2013); Wang et al. (2017)) is that the variance risk premium is included in the risk-neutral dynamics under the mLRNVR. The absence of a variance risk premium in the risk-neutral dynamics under the LRNVR is noted in Hao and Zhang (2013), where it is shown that both empirical studies and theoretical results indicated that the GARCH models under the LRNVR did not capture the variance premium.
We then find the theoretical VIX squared value as the conditional risk-neutral expectation of the arithmetic mean variance over the next 21 trading days under the mLRNVR. Specifically, the GARCH implied VIX formulas are derived using the features of square-root stochastic autoregressive volatility (SR-SARV) models. We apply several calibration methods to estimate the model parameters using various sets of time series data, and compare the theoretical formula performances with the market data. Various combinations of the time series data of the daily closing price of the S&P 500 index and the CBOE VIX are used to find the maximum likelihood estimation of the GARCH models. The corresponding implied VIX time series are then calculated from the calibrated model. Similar to the empirical evidence in Hao and Zhang (2013) and Wang et al. (2017), when only the S&P 500 returns are used for estimation, the GARCH implied VIX is consistently and significantly lower than the CBOE VIX. When the CBOE VIX is used for estimation, the implied VIX fits the statistical properties of the CBOE VIX and matches the CBOE VIX data better. The numerical results provide evidence that the GARCH option pricing under the mLRNVR is more suitable to price volatility. In particular, the numerical results show that the EGARCH model using both returns and VIX data under the mLRNVR provides the best fit to the sample data. Therefore, we recommend that the option pricing in the GARCH framework should use the EGARCH model with returns and VIX data under the mLRNVR for the best results in the future. In the case of GARCH(1,1), we also compare the diffusion limit of the GARCH process under the physical measure and the mLRNVR risk-neutral measure to show that the variance premium is captured in the risk-neutral dynamics.