# Option Pricing with Time-Varying Volatility Risk Aversion††thanks: We are
grateful to Jens Jackwerth for valuable comments.__

Peter Reinhard Hansena and Chen Tongb,c

a _University of North Carolina & Copenhagen Business School_
b _Department of Finance, School of Economics, Xiamen University_
c _Wang Yanan Institute for Studies in Economics, Xiamen University_ Address:
University of North Carolina, Department of Economics, 107 Gardner Hall Chapel
Hill, NC 27599-3305

(_August 27, 2024_)

###### Abstract

We introduce a novel pricing kernel with time-varying _variance risk aversion_
that can explain key pricing kernel puzzles. When combined with the Heston-
Nandi GARCH model, the framework yields closed-form expressions for the VIX.
We also obtain closed-form expressions for option prices by proposing a novel
method that extrapolation from the closed-form VIX.

We estimate the model with S&P 500 returns and option prices and find a
substantial reduction in pricing errors by permitting time-variation in
volatility risk aversion. This reduction is seen for both option pricing and
VIX pricing and both in-sample and out-of-sample. The variance risk ratio
emerges as a fundamental variable in our framework and we show that it is
closely related to economic fundamentals and leading measures of sentiment and
uncertainty.

Keywords: Pricing Kernel, Volatility Risk Aversion, GARCH, Heston-Nandi,
Option Pricing, Score-driven

JEL Classification: G12, G13, C51, C52

###  1 Introduction

The pricing kernel is a fundamental concept in economics with roots in Arrow-
Debreu securities, see Arrow and Debreu, (1954). The pricing kernel is the
ratio of risk neutral state probabilities and physical state probabilities.
The pricing kernel can be estimated from option and stock prices, see e.g.
Aït-Sahalia and Lo, (2000), Jackwerth, (2000) and Rosenberg and Engle, (2002),
and this literature has established a number of stylized facts, including so-
called pricing kernel puzzles. The asset pricing model by Lucas, (1978), which
is the cornerstone of asset pricing theory, implies a pricing kernel that is
monotonically decreasing in aggregate wealth. One puzzle is that the empirical
pricing kernel is not monotonically decreasing as implied by Lucas, (1978).
The pricing kernel is typically found to be upward sloping in some regions
when plotted against the market return.

To address the non-monotonicity, Christoffersen et al., (2013) generalized the
model in Heston and Nandi, (2000) by adding a variance risk premium to the
model (in addition to the equity risk premium). Their model is tractable and
yields close-form expressions for option prices, which explains that this
model has become popular and widely adopted in this literature. The shape of
the pricing kernel in Christoffersen et al., (2013) is defined by the variance
risk premium. Specifically, a negative variance risk premium induces a
$U$-shaped pricing kernel, which is the shape they obtain in their empirical
analysis.While the framework in Christoffersen et al., (2013) can resolved one
pricing puzzle it is at odds with another stylized fact, which that the shape
of the pricing kernels is time varying.

When the pricing kernel is estimated over long sample periods, it is typically
found to have an $U$-shape, see e.g. Cuesdeanu and Jackwerth, (2018), but
during some periods it has distinctly different shapes, including hump-shaped
pricing kernels. The time variation in the pricing kernel can be expressed as
time-variation in the probability weighting function. Polkovnichenko and Zhao,
(2013) estimate the probability weighting function non-parametrically and find
evidence that it is time varying. When individuals overweight low-probability
events and underweight events with high probabilities it will result in a
probability weighting function with an inverse $S$-shaped (and a $U$-shaped
pricing kernel). This is the shape that Polkovnichenko and Zhao, (2013) obtain
empirically during most of their sample period. But for the years 2004–2007,
their estimated weighting function has a regular $S$-shaped. This implies that
tail events are underweighted, which is consistent with the negative risk
aversion function documented in Aït-Sahalia and Lo, (2000) and Jackwerth,
(2000). Similar results were obtained in Chabi-Yo and Song, (2013) and the
inverted shape is also implied by DAX 30 options during the same years, see
Grith et al., (2017). Time-variation is also documented by Kiesel and Rahe,
(2017) and by Beare and Schmidt, (2016). The latter developed a formal
statistical test for monotonicity and used it to identify days with different
shapes. Interestingly, figure 3 in Christoffersen et al., (2013) also contains
evidence of a pricing kernel whose shape varies over time. In most calendar
years the empirical shape is in line with their estimated $U$-shape. However,
an inverted $U$-shape can clearly be seen for the calendar years 2004–2007,
albeit the authors do not comment it.

In this paper, we make five main contribution. First, we propose a novel
framework for derivative pricing that is based on a pricing kernel with time-
varying risk aversion. When combined with the Heston-Nandi GARCH model, the
risk neutral conditional variance divided by the physical conditional
variance, $\eta_{t}=h_{t}^{\ast}/h_{t}$, emerges as a fundamental quantity.
This _variance risk ratio, $\eta_{t}$,_ is responsible for all time-variation
in the pricing kernel and it is functionally linked to the risk aversion
parameters, which define the curvature of the pricing kernel. We show that the
proposed framework can generate pricing kernels with the variety of shapes
found empirically. Importantly, the framework can also explain the observed
time-variation in pricing kernels.

Second, we proceed to derive closed-form expressions for the VIX and option
prices when $\eta_{t}$ is time-varying. There are important differences
between the case where the future path for $\eta_{t}$ is predetermined and the
case where it is driven by a stochastic process. In the former case, the
pricing kernel has an affine structure and close-form expressions are
available from the moment-generating function for cumulative returns. The
affine structure is lost when the future path for $\eta_{t}$ is random, and it
is not possible to derive closed-form expressions for option prices with
conventional methods. However, we derive a novel adjustment of the closed-form
expressions we have for the predetermined case., which takes advantage of the
closed-form VIX expression for both cases. Our VIX formula is valid regardless
of $\eta_{t}$ being predetermined or follow a random (autoregressive) process.
We can therefore evaluate the VIX with random $\eta_{t}$ and compare it to the
VIX with predetermined $\eta_{t}$. The two are different, sometimes by an
substantial amount, which shows the importance of accounting for random
variation in $\eta_{t}$. We can take advantage of these expressions for the
purpose of option pricing, by determining the (higher) level of volatility,
$h_{t}^{\ast}$, for which the predetermined expression yields the same VIX as
the expression with random $\eta_{t}$. We then port the adjusted volatility
level to the predetermined option pricing formula to price options when
$\eta_{t}$ is random. This analytical approximation method work remarkably
well, as we show empirically, see e.g. Figure D.1 in Appendix D.

The third contribution in this paper is a empirical implementation of the
theoretical results. This requires that we connect the autoregressive model
for $\eta_{t}$ with observed data. Because the HNG model (like other GARCH
models) is an observation-driven model it is natural to use a similar model
structure for $\eta_{t}$. For this purpose we adopt the score-driven approach
by Creal et al., (2013). Specifically, we used the first-order conditions from
the log-likelihood function to define innovations in $\eta_{t}$. This is a
natural and intuitive approach, that continuously updates the value of
$\eta_{t}$ based on empirical evidence. Intuitively this amount revising
$\eta$ in the direction that reduces pricing errors. A key advantage of
observation-driven models is that they are simpler and faster to estimate,
than the parameter-driven counterparts. This feature will be particular useful
when the model is estimates with a long sample and/or a large panel of option
prices.

The fourth main contribution is an empirical application to 30 years of daily
S&P 500 returns, the CBOE VIX, and a panel of option prices. The empirical
results shows that time variation in the pricing kernel typically reduce
pricing errors by 50%, or more, as measured by the root mean square error.
This reduction is seen for both the VIX and option prices, it is seen in-
sample as well as out-of-sample, and uniformly across a variety of
characteristic for options (moneyness, time to maturity, and contemporaneous
VIX level). The estimated time variation in the shape of the pricing kernel is
consistent with the existing literature. The pricing kernel is found to have
an inverted hump shape during the years 2004 to 2007. We also see the inverted
shape during a period around 1993, and more recently in 2017.

Fifth, we connect the variance risk ratio with key variables used in this
literature. As we have stated earlier, a variance risk ratio, $\eta_{t}$,
emerges as a fundamental variable. So it is natural to investigate if this
variable is closely related to other variables that are commonly used in this
literature. Several variables are found to be related to the pricing kernel
puzzles. For instance, investors with heterogeneous beliefs and disagreements
about the physical distribution can bring about a $U$-shaped pricing kernel.
The $U$-shape arises as a consequence of aggregation of investors who take
opposite position in the market due to different beliefs about expected
returns. Time variation in the level of disagreement would translate to time-
variation in the pricing kernel, see e.g. Shefrin, (2001, 2008), Bakshi and
Madan, (2008), Bakshi et al., (2010). A similar argument can be made sentiment
and uncertainty among investors, see e.g. Han, (2008), Polkovnichenko and
Zhao, (2013), Baker et al., (2016), and Bollerslev et al., (2018). Bali and
Zhou, (2016) argue that market uncertainty can be approximated by the variance
risk premium (VRP), see Carr and Wu, (2008). The variance risk ratio,
$\eta_{t}$, is obviously related to the VRP, which has been shown to predicts
stock returns, see Bollerslev et al., (2009). A straight forward explanation
is that preferences are state-dependent preferences, with market uncertainty
as a possible state variable, see e.g. Grith et al., (2017). We connect our
results with this literature by showing that the monthly average of the
variance risk ratio, $\eta_{t}$, is closely related to commonly used measures
of sentiment, uncertainty, disagreement, etc. as well as other economic
variables that are commonly used in this literature.

Our paper is related to a large body of literature. Our estimation method
builds on the idea of estimating parameter using information about the
physical and risk-neutral probability measures simultaneously, as advocated by
Chernov and Ghysels, (2000). Another related strand of papers are those that
sought to explain the pricing kernel puzzles by augmenting existing models
with additional state variables, see e.g. Chabi-Yo et al., (2008), Chabi-Yo,
(2012), Brown and Jackwerth, (2012) and Song and Xiu, (2016). The variance
risk ratio, $\eta_{t}$, may be interpreted as an additional state variable.
Albeit it emerges naturally from the model structure and it can inferred from
data using the observation-drive model for $\eta_{t}$. Our paper is also
related to Barone-Adesi et al., (2008), who estimate a GJR GARCH model for
returns (under the physical measure). They do not specify a pricing kernel,
but simply assume that the model under the risk neutral measure is also a GJR
GARCH model, whose time-varying parameters are calibrated by minimizing
pricing errors in daily options. We also seek to minimize pricing errors
(through maximization of the log-likelihood function) but our approach is
fundamentally different. Our starting point is a pricing kernel, which in
conjunction with model under $\mathbb{P}$, defines the model under
$\mathbb{Q}$. In the empirical comparison, Christoffersen et al., (2013)
explored the same ideas as in Barone-Adesi et al., (2008), and estimated HNG
models under both $\mathbb{P}$ and $\mathbb{Q}$. Because this is not deduced
from their pricing kernel, the referred to is an ad-hoc model. We show that
the pricing kernel with a time varying variance risk ratio, is not coherent
with (homogeneous) HNG models under both $\mathbb{P}$ $\mathbb{Q}$.

Our paper is also related to the literature on time varying risk aversion, via
the latter’s connection with the local shape of the pricing kernel. A seminal
paper with time-varying risk aversion is Campbell and Cochrane, (1999), and
building on their framework, Li, (2007) showed how time-varying risk aversion
interacts with time-varying risk premium. Moreover, González et al., (2018)
find that time-varying risk aversion is one of the main determinants of stock
market betas. A recent paper by Bekaert et al., (2020) introduced a no-
arbitrage asset pricing model with a time-varying risk aversion, for the
pricing of equities and corporate bonds where the risk aversion is driven by
the uncertainty shock.

The rest of the paper is organized as follows. We present the theoretical
model in Section 2 and derive closed-form pricing formulae for the VIX and
option prices in Section 3. In Section 4 we introduce a score-driven model for
the variance risk ratio. An empirical application with nearly 30 years of
return, VIX, and option data is presented in Section 5. In show how the
variance risk ratio relates of well known measures of sentiment, disagreement,
and uncertainty in Section 6. A summery and some concluding remarks are give
in Section 7. All proofs are presented in Appendix A, and the derivation of
terms related to the score are derived in Appendix B. A simplification of the
expression for the log-likelihood function is derived and presented in
Appendix C, and we evaluated the accuracy of our approximation method in
Appendix D.

###  2 The Model

The observed variables include the daily returns,
$R_{t}\equiv\log\left(S_{t}/S_{t-1}\right)$, where $S_{t}$ is the underlying
asset price on day $t$, and a vector of derivative prices, denoted $X_{t}$. We
will work with two filtrations,
$\mathcal{F}_{t}=\sigma(\\{R_{j},X_{j}\\},j\leq t)$ and
$\mathcal{G}_{t}=\sigma(\\{R_{j},X_{j-1}\\},j\leq t)$, where the latter arises
naturally from a factorization of the likelihood function. We obviously have
$\mathcal{G}_{t}\subset\mathcal{F}_{t}$ for all $t$. We model returns with the
classical Heston-Nandi GARCH (HNG) model, which has a convenient structure for
derivative pricing. This model is given by the two equations:

$\displaystyle R_{t+1}$ $\displaystyle=$ $\displaystyle
r+(\lambda-\tfrac{1}{2})h_{t+1}+\sqrt{h_{t+1}}z_{t+1},$ (1) $\displaystyle
h_{t+1}$ $\displaystyle=$ $\displaystyle\omega+\beta
h_{t}+\alpha(z_{t}-\gamma\sqrt{h_{t}})^{2},$ (2)

where $r$ is the risk free rate; $\lambda$ is the equity risk premium; and
$h_{t+1}=\mathrm{var}(R_{t+1}|\mathcal{F}_{t})$ is the daily conditional
variance, see Heston and Nandi, (2000). The return shock, $z_{t+1}$, is
assumed to be independent and identically distributed (iid) with a standard
normal distribution, $N(0,1)$. The HNG model captures time-variation in the
conditional variance, as is the case for ARCH and GARCH models, see Engle,
(1982) and Bollerslev, (1986). But the HNG model also allows for dependence
between returns and volatility (leverage). It is simple to verify that
$\mathrm{cov}_{t-1}(R_{t},h_{t+1})=-2\gamma\alpha h_{t}$, which shows that the
magnitude of the leverage effect is defined by $\gamma$. A leverage effect is
required to generate the empirically important smile in option prices. The
dynamic structure of the HNG model is carefully crafted to yield a closed-form
solution for option valuation, and Heston and Nandi, (2000) showed that the
continuous limit of the HNG model (as the time interval between observations
shrinks to zero) yields a variance process, $h_{t}$, that converges weakly to
a continuous-time square-root variance process, see Feller, (1951), Cox et
al., (1985), and Heston, (1993).

Option pricing with GARCH models is possible with a simple risk-neutralization
by Duan, (1995), which uses a _locally risk-neutral valuation relationship_
(LRNVR). This method is equivalent to a pricing kernel with a risk premium for
equity risk, see Huang et al., (2017). Unfortunately, the LRNVR pricing kernel
is inadequate for explaining many differences between $\mathbb{P}$ and
$\mathbb{Q}$, including the variance risk premium, see Hao and Zhang, (2013)
and Christoffersen et al., (2013). This motivated Christoffersen et al.,
(2013) to include a variance risk premium in the pricing kernel (in addition
to the equity risk premium). We review their model next and then proceed to
introduce the time-varying variance risk premium.

####  2.1 Pricing Kernel with Constant Parameters

The pricing kernel introduced in Christoffersen et al., (2013) is given by:

$\frac{M_{t+1}}{M_{t}}=\left(\frac{S_{t+1}}{S_{t}}\right)^{\phi}\exp\left[\delta+\pi
h_{t+1}+\xi(h_{t+2}-h_{t+1})\right],$

see Christoffersen et al., (2013, eq. 9). Corsi et al., (2013) showed that it
can conveniently be expressed as

$M_{t+1,t}=\frac{M_{t+1}}{\mathbb{E}_{t}M_{t+1}}=\frac{\exp[\phi R_{t+1}+\xi
h_{t+2}]}{\mathbb{E}_{t}[\exp(\phi R_{t+1}+\xi h_{t+2})]},$ (3)

where the conditional expectation is only taken with respect to the natural
filtration for returns under $\mathbb{P}$, i.e.
$\mathbb{E}_{t}(\cdot)=\mathbb{E}^{\mathbb{P}}(\cdot|R_{t},R_{t-1},\ldots)$.
The expression above shows that the pricing kernel depends on both equity risk
and variance risk, where the magnitude of variance risk is characterized by
$\xi$. The equity risk is govern by $\phi$, as well as $\xi$, because
$h_{t+2}$ depends on $R_{t+1}$. This structure implies the following dynamic
under the risk neutral measure,

$\displaystyle R_{t+1}$ $\displaystyle=$ $\displaystyle
r-\tfrac{1}{2}h_{t+1}^{*}+\sqrt{h_{t+1}^{*}}z_{t+1}^{*}$ $\displaystyle
h_{t+1}^{*}$ $\displaystyle=$
$\displaystyle\omega^{*}+\beta^{*}h_{t}^{*}+\alpha^{*}(z_{t}^{*}-\gamma^{*}\sqrt{h_{t}^{*}})^{2},$

where

$h_{t}^{*}=h_{t}\eta,\quad\omega^{*}=\omega\eta,\quad\beta^{\ast}=\beta,\quad\alpha^{*}=\alpha\eta^{2},\quad\gamma^{\ast}=\gamma-\phi,\quad\eta=(1-2\alpha\xi)^{-1}.$
(4)

The logarithm of the pricing kernel is a quadratic function of the daily
market return,

$\log\frac{M_{t}}{M_{t-1}}=-\lambda(R_{t}-r)+\xi\alpha\tfrac{(R_{t}-r)^{2}}{h_{t}}+\kappa_{0}+\kappa_{1}h_{t},$
(5)

where $\kappa_{0}=\delta+\xi\omega+\phi r$ and
$\kappa_{1}=\pi+\xi(\beta-1+\alpha(\lambda-\tfrac{1}{2}+\gamma)^{2})$, see
Christoffersen et al., (2013). This shows that $\lambda$ and $\xi$ define the
linear and quadratic term of returns, respectively. The former, $\lambda$, is
the equity risk premium, whereas $\xi$ governs the variance risk premium and
$\xi>0$ is associated with a $U$-shaped relationship between $\log
M_{t}/M_{t-1}$ and returns.

Empirically, this relationship tends to be $U$-shaped, but it is not stable
over time as can be seen in Christoffersen et al., (2013, figure 3). In Figure
1 we show parts of figures 3 and 5 in Christoffersen et al., (2013). The four
left panels are based on model-free market prices and the corresponding four
right panels are based on model prices using their estimated model.

(a) Log ratios inferred from market prices

(b) Model-based log ratios

Figure 1: Log ratios of risk-neutral one-month densities and physical one-
month histogram against log returns in monthly standard deviations. Panels (a)
and (b) are subsets of figures 3 and 5, respectively, in Christoffersen et
al., (2013),. The former is based on empirical market prices and the latter
has the corresponding model-based log ratios.

There are periods where the relationship has an inverted U-shape, such as the
years 2004-2007 (of which 2005 and 2006 are shown in Figure 1), but the model
by Christoffersen et al., (2013) does not have the required flexibility
because the quadratic coefficient is constant.

####  2.2 Dynamic Pricing Kernel by means of Variance Risk Aversion

We seek a more flexible pricing kernel that is coherent with the empirical
time-variation. The expression, (5) motivates a structure where $\xi$ is time-
varying $\xi$, because $\xi$ defines the quadratic term in (5). Next, we
introduce our model with time-varying $\xi$, while maintaining the Heston-
Nandi GARCH structure under $\mathbb{P}$.

######  Assumption 1.

Returns are given by the Heston-Nandi GARCH model, (1)-(2), $\mathbb{P}$.

For the pricing kernel we make the following assumption, where
$\mathbb{E}_{t}^{\mathbb{P}}(\cdot)=\mathbb{E}^{\mathbb{P}}(\cdot|\mathcal{G}_{t})$.

######  Assumption 2.

The pricing kernel takes the form:

$M_{t+1,t}=\frac{\exp[\phi_{t}R_{t+1}+\xi_{t}h_{t+2}]}{\mathbb{E}_{t}^{\mathbb{P}}[\exp(\phi_{t}R_{t+1}+\xi_{t}h_{t+2})]},$

where $\phi_{t}$ and $\xi_{t}$ are $\mathcal{G}_{t}$-measurable.

We can now present the dynamic properties under $\mathbb{Q}$.

######  Theorem 1.

Suppose that Assumptions 1 and 2 hold and define
$\eta_{t}\equiv(1-2\alpha\xi_{t})^{-1}$. Then:
$(i)$

$\xi_{t}=\tfrac{1}{2\alpha}\tfrac{\eta_{t}-1}{\eta_{t}},\qquad\phi_{t}=\tfrac{\eta_{t}-1}{\eta_{t}}(\gamma-\tfrac{1}{2})-\tfrac{1}{\eta_{t}}\lambda,$

where the latter is implied by no-arbitrage.
$(ii)$ The return dynamics under risk neutral probability measure $\mathbb{Q}$
is given by:

$\displaystyle R_{t+1}$ $\displaystyle=$ $\displaystyle
r-\tfrac{1}{2}h_{t+1}^{*}+\sqrt{h_{t+1}^{*}}z_{t+1}^{*},$ $\displaystyle
h_{t+1}^{*}$ $\displaystyle=$
$\displaystyle\omega_{t}^{*}+\beta_{t}^{*}h_{t}^{*}+\alpha_{t}^{*}(z_{t}^{*}-\gamma_{t-1}^{*}\sqrt{h_{t}^{*}})^{2},$

and the two risk aversion parameters both depend on the variance risk ratio,
and the relationships between parameters under $\mathbb{P}$ and $\mathbb{Q}$,
are given by $\omega_{t}^{*}=\omega\eta_{t}$,
$\beta_{t}^{*}=\beta\tfrac{\eta_{t}}{\eta_{t-1}}$,
$\alpha_{t}^{*}=\alpha\eta_{t}\eta_{t-1}$ and
$\gamma_{t}^{*}=\tfrac{1}{\eta_{t}}(\gamma+\lambda-\tfrac{1}{2})+\tfrac{1}{2}$.
$(iii)$ We have

$\eta_{t}=\frac{h_{t+1}^{\ast}}{h_{t+1}}.$ (6)

Theorem 1 uncovers several interesting features of the risk neutral
probability measure and its dynamic properties.

First, the _variance risk ratio_ , $\eta_{t}$, emerges as a fundamental
quantity.111The variance risk ratio, $\eta_{t}$, is obviously related to the
variance risk premium, which concerns future expectations of
$h_{t+1}^{\ast}-h_{t+1}$. In the present context, $\log\eta_{t}$ could serve
as a model-based measure of the variance risk premium, as an alternative to
the model-free measures by Carr and Wu, (2008) and Bollerslev et al., (2011).
In Christoffersen et al., (2013) the ratio $h_{t}^{\ast}/h_{t}$ is held
constant. So their time-variation in the variance risk premium is driven by
the concavity of the logarithm (similar to Jensen’s inequality) and defined by
the level of volatility, $h_{t}$. This variable appears in all the expressions
that relates $\mathbb{P}$ and $\mathbb{Q}$, and is straight forward to
interpret. Its functional relationship with $\xi_{t}$, shows that $\eta_{t}$
is an indirect measure of aversion to volatility risk. When $\eta$ is large,
agents require a large compensation for taking on variance risk, while a small
value of $\eta$ corresponds to a large appetite for variance risk. Note that
$\xi_{t}$ is positive if $\eta_{t}>1$, which generates a $U$-shape in (5)
whereas the inverted $U$-shape arises if $\eta_{t}<1$, which implies
$\xi_{t}<0$.222Because $\alpha>0$ in the Heston-Nandi GARCH model.

An important reason for allowing the variance risk ratio to be time-varying,
is that it permits $h_{t}^{\ast}$ to depend on information provided by
derivative prices. In conventional GARCH models, the conditional variance
depends on lagged returns, such that $h_{t+1}\in\sigma(R_{t},R_{t-1},\ldots)$.
An indirect implication of the variance risk ratio, $\eta=h_{t}^{\ast}/h_{t}$,
being constant is that $h_{t+1}^{\ast}$ is also driven exclusively by lagged
returns, and there is no way for derivative pricing to influence the dynamic
properties $h_{t+1}^{\ast}$. However, once the model permits time-variation in
$\eta$ is it possible for the information in derivative prices to influence
$h_{t+1}^{\ast}$. We take advantage of this feature in our empirical
implementation of the model.

Figure 2: Properties of cumulative returns with a low volatility risk ratio,
$\eta=0.73$, (red dashed lines) and a high volatility risk ratio, $\eta=1.69$,
(solid blue lines). Log-ratios of the $\mathbb{Q}$-to-$\mathbb{P}$ densities
for cumulative returns over one month are shown in the upper left panel, news
impact curves under $\mathbb{Q}$ are shown in the upper right panel, and
unconditional skewness and kurtosis of multi-period cumulative returns under
the risk neutral measure $\mathbb{Q}$ are shown in the two lower panels. The
values for the low and high levels of the volatility risk ratio are the $10\%$
and $90\%$ quantiles for $\eta_{t}$ in our empirical analysis in Section 5.

Second, the structure in Theorem 1 nests the Black-Scholes model as a special
case. This classical model arises when $h_{t}$ and $\eta_{t}$ are constant
with $\eta_{t}=1$, because the latter leads to the pricing kernel given by the
power utility function, see Rubinstein, (1976). When $\eta_{t}$ is time-
varying we have time variation in the risk aversion parameters, which is
similar to the habit formation models, see Campbell and Cochrane, (1999).

Third, this framework can generate a variety of shapes of the pricing kernel.
Based on the empirical estimates of the model, which we present in Section 5,
we can plot the shape of the pricing kernel for different levels of the
variance risk ratio. From the empirical distribution of $\eta_{t}$, we set
$\eta_{\mathrm{low}}=0.73$ and $\eta_{\mathrm{high}}=1.69$, which correspond
to the 10% and 90% quantiles, respectively. The upper left panel of Figure 2
displays the shape of the estimated pricing kernel for cumulated returns over
one month, when $\eta_{\mathrm{low}}=0.73$ and $\eta_{\mathrm{high}}=1.69$.
The low value of $\eta$ produce an inverted $U$-shape similar to that observed
during the years 2004 to 2007, whereas the high value of $\eta$ results in a
shape with a pronounced $U$-shape.

Fourth, the model generates time-varying leverage under $\mathbb{Q}$, which
manifests itself in a number of ways. A time-varying leverage effect will
impact the distribution of cumulative returns, and the lower panels of Figure
2 show the skewness and kurtosis of cumulative returns under $\mathbb{Q}$ as a
function of days returns are cumulated over (the $x$-axis). These results are
based on 100,000 simulations.333We simulate using the estimated “SHNG model”
based on option prices, see the second last column of Table 3. Skewness and
kurtosis are shown for the two levels of the variance risk ratio,
$\eta_{\mathrm{low}}$ and $\eta_{\mathrm{high}}$, and these result in
distinctively different levels of skewness and kurtosis under $\mathbb{Q}$.
Both skewness and kurtosis are far more pronounced for the large value of
$\eta$. Another way to illustrate the leverage effect is with the news impact
curve by Engle and Ng, (1993). The news impact curve under $\mathbb{Q}$ is
defined by
$\mathrm{NIC}(z^{\ast})=\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+1}^{\ast}|z_{t}^{\ast}=z^{\ast})-\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+1}^{\ast}|z_{t}^{\ast}=0)$,
and from Theorem 1 it follows that

$\mathrm{NIC}(z^{\ast})=\alpha\eta_{t}\eta_{t-1}z^{\ast
2}-2\alpha\gamma\eta_{t}\sqrt{h_{t}^{*}}z^{\ast},$ (7)

whose shape depends on the variance risk ratio. The upper right panel of
Figure 2 presents the news impact curve $\mathrm{NIC}(z^{\ast})$ for the two
different level of $\eta$, where we evaluate (7) with $h_{t}^{\ast}$ equal to
its unconditional mean under $\mathbb{Q}$, as estimated in Section 5, and with
$\eta_{t}$ and $\eta_{t-1}$ equal to either $\eta_{\mathrm{low}}=0.73$ (red
dashed line) or $\eta_{\mathrm{high}}=1.69$ (solid blue line). We observe the
asymmetric shape of the news impact curve, see e.g. Chen and Ghysels, (2011).
However, the level of asymmetry depends on the level of $\eta$, such that
volatility is far more responsive to negative return shocks when $\eta$ is
high, than when $\eta$ is small.

In the next section we turn to derivative pricing under the risk neutral
dynamic of Theorem 1.

###  3 Derivative Pricing with Time-varying Volatility Risk Aversion

In this section we consider cases where $\eta_{t}$ is time-varying. We derive
closed-form expression for the VIX under the assumption that $\eta_{t}$
follows an autoregressive process of order one, AR(1). We also derive closed-
form expressions for option prices when the path of $\eta_{t}$ is
predetermined until maturity date, but otherwise arbitrary. The we introduce
the adjustment, that yields approximate closed-form expressions for option
prices when $\eta_{t}$ follow an AR(1) process. [Smile-type adjustmed)

####  3.1 VIX Pricing

Once the dynamic properties of $\eta_{t}$ are know under $\mathbb{Q}$, it is
relatively straight forward to obtain closed-form expressions for VIX pricing.
This only requires expressions for $\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})$,
$k=1,\ldots,M$, because the VIX is given by
$A\times\sqrt{\frac{1}{M}\sum_{k=1}^{M}\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})}$,
where $A=100\sqrt{252}$ is the annualizing factor.

The next theorem presents VIX pricing for the case where $\eta_{t}$ is an
AR(1) process with martingale difference innovations.444Extending to the case
where $\eta_{t}$ follows an ARMA(p,q) is possible at the expense of less
transparent expressions.

######  Assumption 3.

Suppose that $\eta_{t}=(1-\theta)\zeta+\theta\eta_{t-1}+\varepsilon_{t}$ where
$|\theta|<1$ and $\varepsilon_{t}$ is independent of $R_{t}$ with
$\mathbb{E}^{\mathbb{Q}}(\varepsilon_{t}|\mathcal{G}_{t-1})=0$ and
$\mathrm{var}^{\mathbb{Q}}(\varepsilon_{t}|\mathcal{G}_{t-1})=\sigma^{2}<\infty$.

######  Theorem 2 (VIX pricing).

Suppose that Assumptions 1-3 hold. Then the VIX is given by

$\mathrm{VIX}(M,\sigma^{2},\eta_{t},h_{t+1}^{\ast})=\sqrt{a_{1}(\eta_{t},M)+a_{2}(M)\sigma^{2}+a_{3}(\eta_{t},M)h_{t+1}^{*},}$
(8)

where $a_{1}(\eta_{t},M)$, $a_{2}(M)$, and $a_{3}(\eta_{t},M)$ are given in
Appendix A.555We have suppressed the three terms dependence on $\theta$,
$\zeta$, and the parameters of the HNG model and (8) suppress an approximation
term that is negligible in practice. This term, denoted $\Delta$ in Lemma A.1,
is bounded.

Constant variance risk aversion corresponds to $\sigma^{2}=0$ (so that
$\eta_{t}=\zeta$), and the VIX formula (8) simplifies to
$\mathrm{VIX}^{2}=a_{1}(\zeta,M)+a_{3}(\zeta,M)h_{t+1}^{*}$ in this special
case. So the randomness of $\eta_{t}$, $\sigma^{2}>0$, add the term,
$a_{2}(M)\sigma^{2}$ to the expression for $\mathrm{VIX}^{2}$. Based on our
empirical estimates, this term accounts accounts for 1.3%-6.1% of
$\mathrm{VIX}^{2}$.

####  3.2 Option Pricing

In this section we derive the option pricing formula for the case where the
variance risk aversion is time-varying. It is important to distinguish between
the case where the future path for $\eta_{t}$ is predetermined and the case
where $\eta_{t}$ is not fully predictable. For an option with $M$ periods to
maturity, the relevant path of $\eta_{t}$ is
$\eta_{t+1:t+M}=(\eta_{t+1},\ldots,\eta_{t+M})$, and its expected path is
denoted
$\eta_{t+1:t+M}^{e}=(\mathbb{E}_{t}^{\mathbb{Q}}\eta_{t+1},\ldots,\mathbb{E}_{t}^{\mathbb{Q}}\eta_{t+M})$.
When $\eta_{t}$ is predetermined, we have $\eta_{t+1:t+M}=\eta_{t+1:t+M}^{e}$
and the affine structure is preserved., such that closed-form option prices
formula can be derived from the moment generating function. This formula does
not rely on a particular model for $\eta_{t+1:t+M}$, and we do not need to
specify one at this point. The VIX pricing formula revealed how uncertainty
about future path of $\eta_{t+1:t+M}$ is associates with a premium, as defined
by the term $a_{2}(M)\sigma^{2}$. With stochastic variation in $\eta_{t}$, the
pricing kernel is no longer affine. However, we can obtain an analytical
approximation, which utilizes a clever adjustment we obtain from the VIX
pricing formula.

#####  3.2.1 Case with Predetermined Variance Risk Aversion

The option pricing formula for the predetermined case can be deduced from the
moment generating function of cumulative returns. This is a relatively
straight forward extension of the result with a constant $\eta_{t}$, which is
derived in Christoffersen et al., (2013).

When $\eta_{t+1:t+M}$ is predetermined, it follows from Theorem 1 that
$\phi_{t}$, $\xi_{t}$, and the GARCH parameters under $\mathbb{Q}$,
$(\omega_{t}^{*},\alpha_{t}^{*},\beta_{t}^{*},\gamma_{t}^{*})$, are also
predetermined. This preserves the affine structure of pricing kernel and we
can obtain a closed-form pricing formula, which we denote by

$C(S_{t},M,K;0,\eta_{t+1:t+M},h_{t+1}^{\ast}).$ (9)

The zero used as the fourth argument indicates the predetermined case, where
$\sigma^{2}=\mathrm{var}_{t}^{\mathbb{Q}}[(\eta_{t+j}-\eta_{t+j}^{e})^{2}]$ is
zero. This pricing formula simplifies to that in Christoffersen et al., (2013)
in the special case where $\eta_{t+1:t+M}=(\eta,\ldots,\eta)$.

######  Theorem 3.

Suppose that Assumptions 1 and 2 hold and that $\eta_{t+1:t+M}$ is
predetermined. Then the option pricing formula is given by

$C(S_{t},M,K;0,\eta_{t+1:t+M},h_{t+1}^{\ast})=S_{t}P_{1}(t)-K\exp(-rM)P_{2}(t),$

where

$\displaystyle P_{1}(t)$ $\displaystyle=$
$\frac{1}{2}+\frac{\exp(-rM)}{\pi}\int_{0}^{\infty}{\rm
Re}\left[\frac{K^{-i\varphi}g_{t,M}^{*}(i\varphi+1)}{i\varphi
S_{t}}\right]d\varphi$, $\displaystyle P_{2}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}+\frac{1}{\pi}\int_{0}^{\infty}{\rm
Re}\left[\frac{K^{-i\varphi}g_{t,M}^{*}(i\varphi)}{i\varphi}\right]d\varphi,$

and where
$g_{t,M}^{*}(s)=\mathbb{E}_{t}^{\mathbb{Q}}(\exp(s\sum_{i=1}^{M}R_{t+i}))$ is
the moment generating function for cumulative returns under $\mathbb{Q}$. It
has the affine form,

$\displaystyle g_{t,M}^{*}(s)$
$\displaystyle=\exp\left(\text{$\mathcal{A}$}_{T}(s,M)+\mathcal{B}_{T}(s,M)h_{t+1}^{*}\right),\quad
T=t+M,$

where the expressions for $\mathcal{A}_{T}(s,m)$ and $\mathcal{B}_{T}(s,m)$
are given in the Appendix, see (A.6) and (A.7).

#####  3.2.2 Option Pricing with Stochastic Variance Risk Aversion

We now turn to the more challenging situation where the $\eta_{t+1:t+M}$ is
not fully predictable. One of the main contributions in this paper is the
option pricing formula for this case. When $\sigma>0$, the moment generating
function, $g_{t,M}^{*}(s)$, does not have an affine form. This thwarts the
standard approach to obtaining a closed-form option pricing formula, and it is
common to resort to simulation methods in this situation, which can be
computationally burdensome. We will present a simpler alternative, which
entails a simple adjustment of the existing pricing formula. We take advantage
of the fact that we have a closed-form expressions for the VIX in both the
predetermined case and the stochastic case. In the stochastic case,
$\sigma>0$, the idea is to consider the counterfactual situation, $\sigma=0$,
and determine an adjusted level of volatility, $\tilde{h}_{t+1}^{\ast}$, which
would have produced the same VIX as $h_{t+1}^{\ast}$ does with $\sigma>0$. We
can interpret $\tilde{h}_{t+1}^{\ast}$ as the $\eta$-certainty equivalent
volatility level and $\tilde{h}_{t+1}^{\ast}-h_{t+1}^{\ast}$ as a compensation
factor for uncertainty about the path of $\eta_{t}$. The central idea is to
price option with the available closed-form expressions available for the
deterministic case, but at the the adjusted volatility level,
$\tilde{h}_{t+1}^{\ast}$, rather than $h_{t+1}^{\ast}$. Because this adjusted
volatility yields the correct VIX, the hope is that this “VIX-fix” will also
price options correctly. Note that the VIX-adjustment factor depends on $M$.

######  Proposition 1 (Certainty equivalent volatility).

Suppose that Assumptions 1-3 hold. An approximate option price formula is
given by

$C(S_{t},M,K;0,\eta_{t+1:t+M}^{e},\tilde{h}_{t+1}^{\ast}),$

where $\tilde{h}_{t+1}^{\ast}=h_{t+1}^{\ast}+\psi(\eta_{t},\sigma^{2},M)$ and
the compensation factor,

$\psi(\eta_{t},\sigma^{2},M)=\frac{a_{2}(M)}{a_{3}(\eta_{t},M)}\sigma^{2},$

solves
$\mathrm{VIX}(M,\sigma^{2},\eta_{t},h_{t+1}^{\ast})=\mathrm{VIX}(M,0,\eta_{t},h_{t+1}^{\ast}+\psi(\eta_{t},\sigma^{2},M))$,
and where $a_{2}(M)$ and $a_{3}(\eta_{t},M)$ are those in Theorem 2.

Since the VIX is monotonically increasing in $h^{\ast}$ and $\sigma^{2}$, we
have
$\mathrm{VIX}(M,0,\eta_{t},h_{t+1}^{\ast})<\mathrm{VIX}(M,\sigma^{2},\eta_{t},h_{t+1}^{\ast})$,
for $\sigma^{2}>0$, and the $\eta$-certainty equivalent volatility level,
$\tilde{h}_{t+1}^{\ast}$, is therefore larger than $h_{t+1}^{\ast}$. Note that
this method does not require the use of actual VIX prices but solely rely on
the model-based expression for the VIX. The correction, $\psi$, is determined
from the model-implied expectation,
$\mathbb{E}_{t}^{\mathbb{Q}}(\sum_{k=1}^{M}h_{t+k}^{*})$, which is evaluated
at both the actual value of $\sigma^{2}$ and the counterfactual value,
$\sigma^{2}=0$.

In our empirical analysis, we find that it is important to account for the
stochastic variation in $\eta_{t}$. The magnitude of the corrections is given
in Table 1 for the same model used in Figure 2. The table presents a
decomposition of the (squared) VIX into its three components,
$\mathrm{VIX^{2}}(M,\sigma^{2},\eta_{t},h_{t+1}^{\ast})=a_{1}(\eta_{t},M)+a_{2}(M)\sigma^{2}+a_{3}(\eta_{t},M)h_{t+1}^{*}$,
into their relative share of the squared VIX. Random variation in $\eta_{t}$
is responsible for the second term, $a_{2}(M)\sigma^{2}$, which is
monotonically increasing in maturity, $M$, while the third term, involving
$h_{t+1}^{*}$, is decreasing in $M$. The intuition behind this is that
$\eta_{t}$ is more persistent than volatility, such that the current level of
volatility, $h_{t+1}^{*}$, is less relevant at longer maturities than is the
stochastic innovation to $\eta_{t}$. The last row of Table 1 reports the
adjustment of $h_{t+1}^{*}$ in percent, as defined by the average value of
$\psi(\eta_{t},\sigma^{2},M)/h_{t+1}^{*}$. Note that this term is
approximately equal to the ratios between the values in the second and third
rows, because
$\psi(\eta_{t},\sigma^{2},M)=\frac{a_{2}(M)}{a_{3}(\eta_{t},M)}\sigma^{2}$.
The correction is a modest 1.5% with one month to maturity, and increases to a
whopping 49% at six month to maturity. This is again an implication of the
relative persistence of $\eta_{t}$ and $h_{t}^{\ast}$.

| Time to Maturity (Month)
---|---
|  $M=1$ |  $M=2$ |  $M=3$ |  $M=4$ |  $M=5$ |  $M=6$
$\mathrm{VIX}^{2}$ Components
$\mathbb{E}\left[\frac{a_{1}(\eta_{t},M)}{\mathrm{VIX^{2}}(M)}\right]$ |  42% |  60% |  69% |  74% |  78% |  80%
$\mathbb{E}\left[\frac{a_{2}(M)\sigma^{2}}{\mathrm{VIX^{2}}(M)}\right]$ |  1.3% |  3.1% |  4.3% |  5.2% |  5.7% |  6.1%
$\mathbb{E}\left[\frac{a_{3}(\eta_{t},M)h_{t+1}^{*}}{\mathrm{VIX^{2}}(M)}\right]$ |  57% |  37% |  27% |  21% |  17% |  14%
Volatility Adjustment
$\mathbb{E}\left[\frac{\psi(\eta_{t},\sigma^{2},M)}{h_{t+1}^{*}}\right]$ |  1.5% |  6.5% |  15% |  25% |  37% |  49%
Table 1: The average decomposition of $\mathrm{VIX}^{2}$ into the three
components in (8). The second term is the contribution ascribed to random
variation in $\eta_{t}$ and it ranges between 1.3% and 6.1%. The last row is
the average adjustment of $h^{\ast}$ that is required to compensate for the
random variation, when the true value is $\sigma=0.042$ (our point estimate
based on option prices). The compensation ranges from 1.5% at one month to
maturity to 49% at six months to maturity.

###  4 Observation-Driven Model for Variance Risk Aversion

In this section, we develop an observation-driven model for $\eta_{t}$, which
makes it possible to implement the results of the previous section. In
practice, we need to infer the current value of $\eta_{t}$ and estimate the
dynamic model for $\eta_{t}$. To this end, we adopt an intuitive observation-
driven model for $\eta_{t}$, inspired by the score-driven framework by Creal
et al., (2013). We use the first-order conditions for minimizing derivative
pricing errors to define innovations to $\eta_{t}$, such that $\eta$ is
adjusted to reduce pricing errors (on average).666Alternatively, one could
adopt state space approach. This is computationally more complicated without
necessarily providing any benefit. Even when the true model is a state space
model, observation-driven models are rarely found to be inferior, see Koopman
et al., (2016). The first-order condition of the log-likelihood function is
presented in the empirical section. At this point, it suffices to express the
log-likelihood function in its generic form,
$\sum_{t=1}^{T}\ell(R_{t},X_{t}|\mathcal{F}_{t-1})$, where $T$ is the sample
size, and observe that it relies on information from both $\mathbb{P}$ and
$\mathbb{Q}$. We factorize
$\ell(R_{t},X_{t}|\mathcal{F}_{t-1})=\ell(R_{t}|\mathcal{F}_{t-1})+\ell(X_{t}|R_{t},\mathcal{F}_{t-1})$,
where the latter can be expressed as $\ell_{t}(X_{t}|\mathcal{G}_{t})$. This
likelihood term measures how well the observed derivative prices are explained
by the statistical model for returns and the pricing kernel.

A score-driven model updates a parameter in the direction dictated by its
first-order conditions of the log-likelihood function, which is known as the
score.777The score-driven approach is locally optimal in the Kullback-Leibler
sense, see Blasques et al., (2015). The relevant score is
$\partial\ell(X_{t}|\mathcal{G}_{t})/\partial\eta_{t}$, because
$\ell(R_{t}|\mathcal{F}_{t-1})$ does not depend on $\eta$. The required
structure for $\varepsilon_{t}$, as specified in Assumption 3, motivates the
choice

$\varepsilon_{t}=\sigma s_{t-1}\qquad\text{where}\quad
s_{t}=\frac{\partial\ell(X_{t}|\mathcal{G}_{t})/\partial\eta_{t}}{\sqrt{\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\partial\ell(X_{t}|\mathcal{G}_{t})/\partial\eta_{t}\right)^{2}\right]}},$
(10)

where $s_{t}$ is $\mathcal{G}_{t+1}$-measurable. The score will satisfy
$\mathbb{E}^{\mathbb{P}}[\partial\ell(X_{t}|\mathcal{G}_{t})/\partial\eta_{t}|\mathcal{G}_{t}]=0$
when evaluated at the true parameter value. Importantly, we will show that
this moment condition also holds under $\mathbb{Q}$, so that
$\mathbb{E}^{\mathbb{Q}}[\varepsilon_{t}|\mathcal{G}_{t-1}]=0$ as required by
Assumption 3. Note that the normalization used in (10) is the Fisher
information under $\mathbb{Q}$, so that
$\mathrm{var}^{\mathbb{Q}}[s_{t}|\mathcal{G}_{t}]=1$. This ensures that
another condition in Assumption 3,
$\mathrm{var}^{\mathbb{Q}}[\varepsilon_{t}|\mathcal{G}_{t-1}]=\sigma^{2}$, is
satisfied. While $\varepsilon_{t}$ is homoskedastic under $\mathbb{Q}$
measure, this need not be the case under $\mathbb{P}$. From
$\eta_{t}=\theta\eta_{t-1}+(1-\theta)\zeta+\sigma s_{t-1}$ and the definition
of $s_{t}$ it follows that $\eta_{t}$ is $\mathcal{G}_{t}$-measurable.

This approach to modeling $\eta_{t}$ is analogous to the way the conditional
variance is modeled in GARCH models. For instance, the GARCH(1,1) model by
Bollerslev, (1986) has, $h_{t}=c+bh_{t-1}+a(r_{t-1}^{2}-h_{t-1})$, such that
changes in $h_{t}$ are driven by discrepancies between squared returns and the
conditional variance. The score model invokes a similar self-adjusting
property, where $\eta_{t}$ is updated whenever the derivative prices indicates
a that a different value of $\eta_{t}$ would result in smaller pricing errors.

This score-drive approach could be adopted with more general specifications
for $\eta_{t}$, than the AR(1) model we have used here. For instance, could
could include exogenous variables and lagged endogenous variables as
explanatory variables, or used a long-memory specification for $\eta_{t}$. A
main challenge to alternative specification is that the derivative prices
formulae would need to be established first. Such extensions might be possible
for the VIX, but would be challenging for option pricing in general. The
score-driven approach used here, can also be based on on first order
conditions for other criteria, than a log-likelihood function, such as a loss
function.

Table 2: Summary Statistics

A: S&P 500 returns and and CBOE VIX
---
|  |  Mean(%) |  Std(%) |  Skewness |  Kurtosis |  Obs.
Returns (annualized) |  5.26 |  17.51 |  -0.257 |  11.824 |  7,431
VIX |  |  19.22 |  7.76 |  2.11 |  10.84 |  7,431
|  |  |  |  |  |
B: Option Price Data
|  | Implied Volatility (%) | Average price ($) |  Obs.
All options |  | 19.25 | 81.61 |  30,633
|  |  |  |  |  |
Partitioned by Moneyness
Delta<0.3 |  | 13.92 | 8.02 |  5,235
0.3$\leq$Delta<0.4 |  | 16.20 | 21.19 |  1,691
0.4$\leq$Delta<0.5 |  | 17.53 | 34.12 |  2,335
0.5$\leq$Delta<0.6 |  | 19.05 | 47.17 |  4,644
0.6$\leq$Delta<0.7 |  | 19.97 | 65.12 |  4,019
0.7$\leq$Delta |  | 22.00 | 146.49 |  12,709
|  |  |  |  |  |
Partitioned By Maturity
DTM<30 |  | 19.09 | 58.54 |  4,597
30$\leq$DTM<60 |  | 19.36 | 58.45 |  6,526
60$\leq$DTM<90 |  | 19.20 | 76.22 |  6,844
90$\leq$DTM<120 |  | 19.72 | 96.27 |  4,731
120$\leq$DTM<150 |  | 19.20 | 106.41 |  3,659
150$\leq$DTM |  | 18.82 | 112.93 |  4,276
|  |  |  |  |  |
Partitioned By VIX Level
VIX<15 |  | 13.90 | 83.72 |  11,207
15$\leq$VIX<20 |  | 17.94 | 79.55 |  8,716
20$\leq$VIX<25 |  | 22.29 | 81.26 |  5,770
25$\leq$VIX<30 |  | 25.85 | 80.82 |  2,599
30$\leq$VIX<35 |  | 29.76 | 77.34 |  1,113
35$\leq$VIX |  | 39.44 | 84.09 |  1,228

Note: Summery statistics for close-to-close S&P 500 index returns, VIX and
option prices from January 1990 to June 2019. The reported statistics for S&P
500 and VIX index include the sample mean (Mean), standard deviation (Std),
skewness (Skew), kurtosis (Kurt), number of observations (Obs). Option prices
are based on closing prices of out-of-the-money call and put options. We
report the average Black-Scholes implied volatility (IV), average price, and
the number of option prices for different partitions of option prices.
“Moneyness” is defined by the Black-Scholes delta. DTM denotes the number of
calendar days to maturity. Data sources: S&P 500 returns from WRDS, VIX from
CBOE’s website. Option prices: Optsum data (1990–1995) and OptionMetrics
(1996–2019).

###  5 Empirical Analysis

####  5.1 Data

Our empirical analysis is based on daily returns for the S&P 500 index, the
CBOE VIX, and derivative prices on the exchange traded fund, SPX, that tracks
the S&P 500 index. Our sample period spans nearly 30 years (7,431 trading
days) from January 2nd, 1990 to June 28th, 2019. Returns are defined from cum-
dividend logarithmic transformed closing prices, downloaded from the Wharton
Research Data Services (WRDS). Daily VIX prices were downloaded from the
website of the Chicago Board Options Exchange (CBOE). The SPX option prices
were obtained from two sources. Prices for the first six year years
(1990–1995) are the so-called Optsum data, purchased from CBOEs website.
Option prices for the remaining 24 years are OptionMetrics data downloaded
from the WRDS database. Option prices were preprocessing following
Christoffersen et al., (2014) and Bakshi et al., (1997). Specifically, we
include out-of-the-money put and call options prices with positive trading
volume and with maturity between two weeks and six months, and put option
prices are converted to in-the-money call options using put-call parity. Much
of the existing literature use a single panel of observations per week,
typically Wednesday. The use of weekly data can make the data set more
manageable, and may be justified if there is little additional information in
daily data. Wednesdays are chosen because liquidity tends to be relatively
high on that day of the week. In our analysis, we use daily option prices,
because the score, $s_{t}$, is need in our equation that updates $\eta_{t}$.
We use six option prices per day, one for each of the maturities 1, 2, 3, 4, 5
and 6 months. In the event that there is not an option with the exact
maturity, then we use the option with the largest available maturity, which is
less than the desired maturity. For each of the six maturities, we retain the
most liquid option, as defined by traded volume.888The number of available
option prices has increased almost 50-fold over our sample period, initially
from about 38 daily option prices to well over 1,700. This results in a total
of 30,633 option prices in our data set. This number is not divisible by six,
because some days in the early part of the sample had fewer than six options
that meet the inclusion criteria.

Table 2 reports descriptive statistics for our data, S&P 500 returns and VIX
in Panel A and option prices in Panel B. The sample average of the VIX
(19.22%) is, as expected, larger than the standard deviation for annualized
returns (17.51%). Their difference reflect the (average) negative volatility
risk premium. The S&P 500 returns exhibit a small negative skewness and and a
high level of kurtosis, where as the skewness in positive for the VIX while
the kurtosis is also large. We also report summery statistics for option
prices, partitioned by moneyness, maturity, and the contemporaneous level of
VIX. Moneyness is defined by the Black-Scholes delta. Those with deltas less
than 0.5 are out-of-the-money call options, and those with deltas larger than
0.5 are based on the out-of-the-money put options. Deep out-of-the-money put
options (i.e. deltas larger than 0.7) are relatively expensive compared with
out-of-the-money calls, which illustrates the well-known stylized volatility
smirk. The implied volatility is relatively insensitive to the time to
maturity, resulting in a flat term structure. The bottom of Panel B has
options partitioned by the level of the contemporary VIX. The implied
volatility is, as expected, increasing in the level of the VIX.

####  5.2 The Log-Likelihood Function

Next we turn to the log-likelihood functions that is used to estimated model
parameters. Our observed data consist of returns, $R_{t}$, and a vector of
derivative prices, $X_{t}$, and we can, without loss of generality, factorize
the log-likelihood function as follows,

$\ell(R_{1},\ldots,R_{T},X_{1},\ldots,X_{T},\mathcal{F}_{0})=\sum_{t=1}^{T}\ell(R_{t},X_{t}|\mathcal{F}_{t-1})=\sum_{t=1}^{T}\ell(R_{t}|\mathcal{F}_{t-1})+\ell(X_{t}|R_{t},\mathcal{F}_{t-1}).$
(11)

This decomposition was also used in Christoffersen et al., (2013),.

The log-likelihood function for returns, $\ell(R_{t}|\mathcal{F}_{t-1})$, is
given from the HNG model and the assumption that return shocks are iid
Gaussian distributed, $z_{t}\sim iidN(0,1)$. Thus,

$\ell(R_{t}|\mathcal{F}_{t-1})=-\tfrac{1}{2}\left[\log(2\pi)+\log
h_{t}+(R_{t}-r-(\lambda-\tfrac{1}{2})h_{t})^{2}/h_{t}\right].$

To complete the model, we need to specify the log-likelihood function for the
vector of derivative prices, $X_{t}$,
$\ell(X_{t}|R_{t},\mathcal{F}_{t-1})=\ell(X_{t}|\mathcal{G}_{t})$. To this
end, we let $X_{t}^{m}\in\mathcal{G}_{t}$ denote the corresponding vector of
model-based prices. The derivative pricing errors are now given by

$e_{t}=X_{t}-X_{t}^{m},$

for which we make the following assumption.

######  Assumption 4.

The pricing errors, $e_{t}\sim iidN(0,\sigma_{e}^{2}\Omega_{N_{t}})$ where
$\Omega_{N_{t}}\equiv\left(\begin{array}[]{cccc}1&\rho&\cdots&\rho\\\
\rho&1&\rho&\vdots\\\ \vdots&\rho&\ddots&\rho\\\
\rho&\cdots&\rho&1\end{array}\right)$, where $N_{t}$ is the number of
derivative prices at times $t$.

Christoffersen et al., (2013) made the stronger assumption that $\rho=0$,
which is clearly violated in our data. The common correlation can be motivated
by a structure where the pricing error have a common term, $u_{t}$, and a
idiosyncratic term $v_{i,t}$, so that $e_{i,t}=u_{t}+v_{i,t}$. With all terms
being uncorrelated we have the structure in Assumption 4, where
$\rho=\sigma_{u}^{2}/(\sigma_{u}^{2}+\sigma_{v}^{2})$, and
$\sigma_{e}^{2}=\sigma_{u}^{2}+\sigma_{v}^{2}$.

It follows that the part of the log-likelihood function that is relate to
derivative prices, is given by

$\ell(X_{t}|\mathcal{G}_{t})=-\frac{1}{2}[N_{t}\log(2\pi\sigma_{e}^{2})+\log\det\Omega_{N_{t}}+\sigma_{e}^{-2}{\rm
tr}(\Omega_{N_{t}}^{-1}e_{t}e_{t}^{\prime})],$ (12)

which can be simplified further, see (C.1) in the Appendix. In our expression
of this part of the log-likelihood we utilized that $\eta_{t}$ is
$\mathcal{G}_{t}$-measurable, so it does not depend on the contemporaneous
derivative prices. A related virtue of the way we have introduced time-
variation in the pricing kernel, is that the model (for returns) under the
physical measure continues to be a time-homogeneous Heston-GARCH model, that
does not depend on derivative prices. So the parameters in the model could, if
needed, be estimated by two-stage estimation methods, where the parameters in
the GARCH model is estimated from returns in a first stage, while the
remaining parameters are estimated in a second stage using derivative
prices.999In our empirical analysis we used joint estimation because it was
straight forward and did not cause any computational issues. Two-stage
estimation can greatly reduce the computational burden, and is often used in
this context, see e.g. Broadie et al., (2007), Corsi et al., (2013),
Christoffersen et al., (2013) and Majewski et al., (2015), and two-stage
estimation could be implemented for our model, if computational considerations
required it.

In our empirical analysis, we use the Vega-weighted option price multiplied by
100, so that options are measures on the same scale as the VIX. This is
particularly important when both the VIX and option prices are used for
estimation to avoid a comparison of apples and oranges. Moreover, if the
elements of $X_{t}$ were measured in different units, the elements of $e_{t}$
would have different variances, which would violate Assumption 4.101010The
pricing errors for the Vega-weighted option prices is approximately the
difference in implied volatility, since the Vega is a measure of the
sensitivity of option price to changes in implied volatility.

####  5.3 Parameter estimation

We estimate HNG model with a time-varying variance risk aversion (SHNG) with
three different sets of derivatives, and the three corresponding
specifications of Christoffersen et al., (2013) (HNG) with constant $\eta$.
The estimated parameters for these six specifications are reported in Table 3.
Parameters are estimated by maximizing the joint log-likelihood function for
the full sample, which spans the period for the beginning of January 1990 to
the end of June 2019. The first row the table identifies the specification
being estimated in each column. HNG refer is the model proposed in
Christoffersen et al., (2013) and SHNG is the score-driven model proposed in
this paper. The entries in brackets, [VIX], [Opt], and [VIX+Opt], denotes the
derivatives used to estimate the models, the VIX, option prices, or both VIX
and option prices. Robust standard errors are provided in parentheses below
each parameter estimate.111111Following Christoffersen et al., (2013), we
impose $\omega$ = 0 in estimation because the non-negativity constraint, which
ensures positive variances, is binding. We also report the persistence of
volatility, $\pi^{\mathbb{P}}$ and $\pi^{\mathbb{Q}}$, under $\mathbb{P}$ and
$\mathbb{Q}$, respectively, where
$\pi^{\mathbb{Q}}=\mathbb{E}^{\mathbb{Q}}[\beta_{t}^{*}+\alpha_{t}^{*}\gamma_{t-1}^{*2}]$
for the SHNG models. We also report the value of the components of the log-
likelihood (that for returns, and that for derivative prices conditional of
contemporaneous returns). Some of these likelihood terms can be interpreted as
pseudo out-of-sample log-likelihood functions, and we have indicated these
with italic font. For instance the HNG[VIX] model is estimated using returns
and the VIX, but the estimated model also yields model-based option prices
that can be compared with actual option prices. The reported log-likelihood
for option prices is, in this case, computed with the resulting option pricing
errors, which involves estimated the sample variance for the HNG model, and
the sample covariance matrix, with constant diagonal elements, for the SHNG
model. We compute, similarly, the log-likelihood for VIX pricing errors for
the model estimated with options prices (and returns). The numbers in bold
indicate the largest log-likelihood within each row.

Table 3: Joint Estimation Results

Model |  HNG |  HNG |  HNG |  SHNG |  SHNG |  SHNG
---|---|---|---|---|---|---
|  [VIX] |  [Opt] |  [VIX+Opt] |  [VIX] |  [Opt] |  [VIX+Opt]
|  |  |  |  |  |
$\beta$ |  0.812 |  0.631 |  0.656 |  0.805 |  0.507 |  0.529
|  (0.009) |  (0.008) |  (0.006) |  (0.008) |  (0.015) |  (0.007)
|  |  |  |  |  |
$\alpha(\times 10^{-6})$ |  1.216 |  1.429 |  1.373 |  4.140 |  5.740 |  5.404
|  (0.088) |  (0.012) |  (0.051) |  (0.496) |  (0.256) |  (0.494)
|  |  |  |  |  |
$\gamma$ |  381.20 |  499.99 |  491.51 |  193.29 |  279.41 |  281.02
|  (13.95) |  (3.81) |  (9.09) |  (11.70) |  (10.28) |  (12.13)
|  |  |  |  |  |
$\lambda$ |  2.502 |  2.280 |  2.291 |  2.390 |  2.193 |  2.299
|  (1.225) |  (0.022) |  (0.039) |  (1.186) |  (0.644) |  (0.682)
|  |  |  |  |  |
$\eta,\zeta$ |  1.466 |  1.227 |  1.256 |  1.299 |  1.094 |  1.126
|  (0.041) |  (0.013) |  (0.022) |  (0.021) |  (0.013) |  (0.025)
|  |  |  |  |  |
$\theta$ |  |  |  |  0.983 |  0.996 |  0.994
|  |  |  |  (0.003) |  (0.001) |  (0.004)
|  |  |  |  |  |
$\sigma$ |  |  |  |  0.089 |  0.042 |  0.049
|  |  |  |  (0.006) |  (0.010) |  (0.002)
|  |  |  |  |  |
$\sigma_{e}$ |  3.715 |  4.084 |  4.056 |  1.007 |  1.892 |  1.860
|  (0.061) |  (0.024) |  (0.022) |  (0.032) |  (0.036) |  (0.031)
|  |  |  |  |  |
$\rho$ |  |  |  |  |  0.148 |  0.173
|  |  |  |  |  (0.010) |  (0.013)
|  |  |  |  |  |
$\pi^{\mathbb{P}}$ |  0.988 |  0.988 |  0.988 |  0.960 |  0.955 |  0.956
$\pi^{\mathbb{Q}}$ |  0.991 |  0.991 |  0.991 |  0.966 |  0.963 |  0.963
|  |  |  |  |  |
LogL |  |  |  |  |  |
$\ell(\rm{Ret})$ |  24426 |  24347 |  24363 |  24525 |  24214 |  24259
$\ell(\rm{VIX})$ |  13925 |  13456 |  13557 |  23630 |  18917 |  19839
$\ell(\rm{Opt})$ |  49749 |  54499 |  54442 |  69141 |  82245 |  82192
$\ell(\rm{VIX+Opt})$ |  63674 |  67955 |  67999 |  89097 |  102053 |  102630

Note: Estimation results for the full sample period (January 1990 to June
2019). The second row refers to the model specification: HNG is the model by
Christoffersen et al., (2013) and SHNG is the model introduced in this paper.
Contents in brackets denote the derivative prices used in the estimation,
which is either the VIX, option prices, or both VIX and option prices.
Parameter estimates are reported with robust standard errors (in parentheses),
$\pi^{\mathbb{P}}$ and $\pi^{\mathbb{Q}}$ refer to the volatility persistence
under $\mathbb{P}$ and $\mathbb{Q}$, respectively. The value each component of
the log-likelihood function is reported in the bottom of the table, where the
largest log-likelihood within each row is indicated with bold font.

There are a number of interesting observations to be made from the empirical
estimates in Table 3. First, the volatility process is found to be highly
persistent across all specifications, and the persistence is generally larger
under risk-neutral measure than under the physical measure. This is consistent
with the existing literature. Second, the equity premium parameter, $\lambda$,
is estimated to be positive an is significant in all specification. Third, the
variance risk ratio, $\eta_{t}$, is typically larger than one. It is constant
in the HNG model, where it is estimated to be significantly larger than one in
all specification. In the SHNG model the unconditional mean of $\eta_{t}$ is
given by, $\zeta,$ which is also found to be significantly larger than one in
all specification. The implication is that the risk neutral volatility,
$h^{\ast}$, is (on average) larger than the physical volatility, $h$. Fourth,
the parameter that define the leverage effects, $\gamma$, is also found to be
significant across all specification.

For the SHNG specification we note that the variance risk ratio, $\eta_{t}$,
is highly persistent, since $\theta$ is estimated to be close to one. In fact
is it estimated to be more persistent than $h_{t}$ and $h_{t}^{\ast}$, across
all SHNG specifications. Second, the coefficients, $\sigma$, that measures the
impact that $s_{t-1}$ has on $\eta_{t}$, is estimated to be positive and
significant. The estimated unconditional variance for $\eta_{t}$, implicitly
given by $\sigma^{2}/(1-\zeta^{2})$ , is similar across all SHNG
specifications, albeit it is slightly larger for the specification estimated
with the VIX, which is consistent with Figure 3 (below), where the estimated
path for $\eta_{t}$ tended to be slightly larger if estimated for the VIX than
if estimated from option prices. Third, all specifications yield similar
likelihood values for the returns, which is to be expected since all rely on
the Heston Nandi GARCH model for returns. The key difference between the HNG
and SHNG specifications is the pricing kernel, and we observed vast
improvements in the likelihood for derivative prices. This reflect a
substantial reduction in the pricing errors resulting from the pricing kernel
with time-varying variance risk aversion. Fourth, estimating the models with
VIX or option prices results in some differences across the estimated
parameters. This is to be expected, since the vector of option prices contains
more information (about future returns) than the VIX also. For instance, the
leverage parameter, $\gamma$, is estimates to be larger with option prices
than with the VIX. Fifth, the correlation between derivative errors is
estimated to be positive, at about 15%, and is very significant. Ignoring this
correlation (assuming it to be zero) would greatly underestimate the condition
variance of
$\partial\ell$($X_{t}$|$R_{t}$,$\mathcal{F}_{t-1}$)/$\partial\eta_{t}$, which
makes the variance of $s_{t}$ become larger than one.

Figure 3: This figure presents time series of the estimated variance risk
ratio, $\eta_{t}$. The blue line is based on the VIX and the red line on
option prices. The two time series are very similar, albeit the $\eta_{t}$
based on the VIX tend to be slightly larger than that based on option prices.

Figure 3 presents the estimated time series of the variance risk ratio,
$\eta_{t}=h_{t+1}^{\ast}/h_{t+1}$, based on the VIX (blue line) and option
prices (red line). The two have a high degree of comovement and, as we
discussed earlier, the unconditional variance of the two series is very
similar. Importantly, the variance risk ratio occasionally falls below one, as
is seen during the periods: 1993-1995, 2004-2007, and 2017. A value below one,
$\eta_{t}<1$, means that investors have an increased appetite for variance
risk, whereas a large value of $\eta$ implies that investors demand additional
compensation for variance risk. The latter was particularly the case during
the financial crises in the immediate aftermath of the Lehman collapse. The
conditional variance was unusually high during this period under physical
measure. The large value for $\eta_{t}$ during this period magnifies the
impact on $h_{t}^{\ast}$, so that it increased even more, and to a level that
was more than twice that of $h_{t}$. The level of $\eta_{t}$ based on VIX tend
to be slightly larger than that based on option prices. Since these VIX is
based on options with 30 days to maturity, whereas our option have maturity
with as much as 180 days to maturity. A possible explanation is that short-
term investors demanding demands larger compensation for variance risk than
long(er)-term investors. This would be consistent with the findings in
Eisenbach and Schmalz, (2016) and Andries et al., (2018).

Table 4: VIX and Option Pricing Performance (in sample)

Model |  HNG |  HNG |  HNG |  SHNG |  SHNG |  SHNG
---|---|---|---|---|---|---
|  [VIX] |  [Opt] |  [VIX+Opt] |  [VIX] |  [Opt] |  [VIX+Opt]
|  |  |  |  |  |
A: RMSE for VIX Pricing
Full Sample |  3.715 |  3.952 |  3.904 |  1.006 |  1.898 |  1.676
|  |  |  |  |  |
B: RMSE for Option Pricing
Full Sample |  4.128 |  3.793 |  3.793 |  2.855 |  1.877 |  1.888
|  |  |  |  |  |
Partitioned by moneyness
Delta<0.3 |  4.743 |  3.528 |  3.611 |  3.729 |  1.921 |  1.887
0.3$\leq$Delta<0.4 |  4.266 |  3.556 |  3.577 |  3.194 |  1.796 |  1.771
0.4$\leq$Delta<0.5 |  3.975 |  3.568 |  3.569 |  2.793 |  1.704 |  1.683
0.5$\leq$Delta<0.6 |  3.858 |  3.671 |  3.652 |  2.393 |  1.493 |  1.486
0.6$\leq$Delta<0.7 |  3.982 |  3.899 |  3.879 |  2.283 |  1.560 |  1.570
0.7$\leq$Delta |  4.003 |  3.978 |  3.960 |  2.712 |  2.110 |  2.153
|  |  |  |  |  |
Partitioned by maturity
DTM<30 |  4.222 |  4.078 |  4.082 |  2.611 |  2.600 |  2.444
30$\leq$DTM<60 |  4.114 |  3.972 |  3.964 |  2.485 |  1.976 |  1.892
60$\leq$DTM<90 |  4.109 |  3.833 |  3.835 |  2.695 |  1.546 |  1.616
90$\leq$DTM<120 |  3.971 |  3.545 |  3.543 |  2.947 |  1.603 |  1.695
120$\leq$DTM<150 |  4.157 |  3.588 |  3.589 |  3.169 |  1.704 |  1.824
150$\leq$DTM |  4.227 |  3.576 |  3.584 |  3.414 |  1.727 |  1.881
|  |  |  |  |  |
Partitioned by the level of VIX
VIX<15 |  4.029 |  2.773 |  2.875 |  2.841 |  1.611 |  1.628
15$\leq$VIX<20 |  3.503 |  3.215 |  3.218 |  2.635 |  1.731 |  1.736
20$\leq$VIX<25 |  3.999 |  4.451 |  4.365 |  2.689 |  2.034 |  2.036
25$\leq$VIX<30 |  4.360 |  4.831 |  4.785 |  3.043 |  2.110 |  2.136
30$\leq$VIX<35 |  4.589 |  4.865 |  4.783 |  3.076 |  2.394 |  2.382
35$\leq$VIX |  7.378 |  7.020 |  7.027 |  4.252 |  2.938 |  2.969

Note: This table reports the in-sample VIX and option pricing performance for
each model in Table 3. We evaluate the model’s option pricing ability through
the root of mean square errors of implied volatility (${\rm RMSE_{IV}}$). We
summarize the results by option moneyness, maturity and market VIX level.
Moneyness is measured by Delta computed from the Black-Scholes model. DTM
denotes the number of calendar days to maturity. The numbers in bold indicate
the minimum value within each row in each time period.

####  5.4 VIX and Option Pricing Performance (In Sample)

We evaluate and compare the performance of the models in terms of VIX and
option pricing, as measured by the RMSE. The VIX pricing error is defined as
the difference between the model-implied VIX and the observed VIX. We use the
same scale to to measure option pricing errors. This is achieved by converting
option prices to implied volatilities, and define these pricing errors as the
differences in between model-based implied volatilities for an option and the
implied volatility for the observed option price. The (in-sample) RMSEs for
these pricing errors are reported in Table 4 for each of the specifications
listed in Table 3.

We present the in-sample RMSEs for VIX pricing in Panel A. These are defined
by

${\rm
RMSE_{VIX}}=\sqrt{\frac{1}{T}\sum_{t=1}^{T}\text{$\left[\mathrm{VIX}_{t}^{m}-\mathrm{VIX}_{t}\right]^{2}$}}\times
100,$

where $\mathrm{VIX}_{t}^{m}$ is model-based and $\mathrm{VIX}_{t}$ is the
observed market VIX. Panel A shows that a time-varying variance risk ratio
result in substantially smaller average pricing errors, relative to the HNG
models which all assumes $\eta$ to be constant. The improvements in VIX prices
that are achieved by permitting time-variation in $\eta$ are substantial and
impressive. The specification that delivers the most accurate VIX pricing is
the SHNG model that is estimated with VIX data. Its RMSE is 73% smaller than
the best performing HNG specification. The other two SHNG specifications are
also substantially better at VIX pricing than any of the HNG specifications.
This is despite the fact that their parameters are estimated to also match
observed option pricing. The SHNG models estimated with option prices involves
an objective that must price options with maturities ranging from one to six
months. This entails a trade-off across maturities, whereas the VIX
specification specifically targets a one-month maturity.

In Panel B of Table 4 we report the in-sample performance in terms of option
pricing. We follow the literature and convert option prices to their
corresponding implied volatilities, as defined by the Black–Scholes formula.
We compare the model-based implied volatility, $\mathrm{IV}^{m}$, with the
market-based implied volatility, $\mathrm{IV}$, where the latter is defined by
the observed option price. The resulting RMSE,

${\rm
RMSE_{IV}}=\sqrt{\frac{1}{\sum_{t=1}^{T}N_{t}}\sum_{t=1}^{T}\text{$\sum_{i=1}^{N_{t}}\left[\mathrm{IV}_{t,i}^{m}-\mathrm{IV}_{t,i}\right]^{2}$}}\times
100,$

is reported the first row of Panel B in Table 4. We also compute the ${\rm
RMSE_{IV}}$ for option partitioned by the characteristics: Moneyness, Time to
maturity, and the contemporaneous volatility index level. The resulting RMSEs
can be used to identify shortcomings in a model, such as its inability to:
generate sufficient leverage effect; capture the dynamic properties; and
generate a proper variance risk premium. The smallest RMSE within each row is
highlighted with boldface.

The SHNG specifications clearly dominate the HNG specifications in terms of
option pricing. have an The average option pricing errors are substantially
better for than the HNG models., and the smallest RMSE is obtained by the
SHNG, which is estimated with option prices. Its RMSE is less than half that
of any of the HNG specifications. Impressively, the SHNG specifications
uniformly dominates all HNG specifications for both VIX and option pricing,
across all partitions, whenever the SHNG model is estimated with option
prices. Even the SHNG model that is estimated without option prices is, on
average, better at pricing option pricing performance than the all estimated
HNG models. This further supports the advantage of allowing for time variation
in the variance risk ratio. We note that the largest gains in pricing accuracy
are found during periods with high volatility and out-of-money put options,
which tends to be options with high variance risk premia.

While the SHNG estimated with option prices delivers the best in-sample option
pricing, the specification that use both VIX and option prices is a close
second. This is not unexpected since the parameter estimates for these two
specifications were very similar, see Table 3. We observe that the SHNG
specification that is estimated with both the VIX and option prices, yields
more accurate option pricing for low delta and at short maturities, and the
SHNG specification estimated with the VIX, delivers impressive option pricing
at for options with one month to maturity. The reason is that VIX is extracted
from option options with approximately 30 days to maturity, so when the VIX is
included in the estimation it tends to improve the pricing at shorter
maturities. This comes at the expense of slightly larger pricing errors for
options with longer maturities.

Interestingly, Christoffersen et al., (2013) also estimated an ad-hoc model
where the ratio of volatility under $\mathbb{P}$ and $\mathbb{Q}$ is not held
constant. This model was used an a benchmark for the HNG model they proposed
in their paper. The ad-hoc model consisted of two separate Heston-GARCH
models, one that was fitted to returns under $\mathbb{P}$,and one that was
estimated with option prices under $\mathbb{Q}$. This model is referred to as
ad-hoc, because it makes no attempt to uncover a pricing kernel that can
explain the differences between $\mathbb{P}$ and $\mathbb{Q}$. Moreover, their
ad-hoc model did not improve on option pricing relative to the HNG model with
a variance-dependent pricing kernel. Our Theorem 1 offers a theoretical
explanation for the poor performance of their ad-hoc models, because Theorem 1
shows a time-varying risk premium is incompatible with Heston-Nandi models
that have time-invariant parameters under both $\mathbb{P}$ and $\mathbb{Q}$.
121212Thus, the ad-hoc model is internally inconsistent, and interestingly it
only leads to minuscule improvements for option pricing, which follows form
the fact that their log-likelihood for option prices only improves by about
1/10.

In the next section we evaluate if the substantial in-sample improvements in
derivative pricing also holds out-of-sample.

Table 5: Out-of-sample VIX and Option Pricing (January 2008 - June 2019)

Model |  HNG |  HNG |  HNG |  SHNG |  SHNG |  SHNG
---|---|---|---|---|---|---
|  [VIX] |  [Opt] |  [VIX+Opt] |  [VIX] |  [Opt] |  [VIX+Opt]
|  |  |  |  |  |
A: RMSE for VIX Pricing
Out-of-sample |  4.355 |  4.235 |  4.194 |  1.286 |  1.873 |  1.692
|  |  |  |  |  |
B: RMSE for Option Pricing
Out-of-sample |  4.668 |  3.870 |  3.949 |  3.413 |  1.849 |  1.905
|  |  |  |  |  |
Partitioned by moneyness
Delta<0.3 |  6.094 |  3.980 |  4.295 |  4.789 |  1.818 |  1.863
0.3$\leq$Delta<0.4 |  5.088 |  3.845 |  3.974 |  4.093 |  1.791 |  1.832
0.4$\leq$Delta<0.5 |  4.497 |  3.458 |  3.568 |  3.480 |  1.496 |  1.551
0.5$\leq$Delta<0.6 |  4.337 |  3.856 |  3.848 |  2.920 |  1.371 |  1.409
0.6$\leq$Delta<0.7 |  4.545 |  4.254 |  4.239 |  2.709 |  1.574 |  1.612
0.7$\leq$Delta |  3.936 |  3.787 |  3.790 |  2.763 |  2.106 |  2.177
|  |  |  |  |  |
Partitioned by maturity
DTM<30 |  4.971 |  4.347 |  4.425 |  3.262 |  2.280 |  2.219
30$\leq$DTM<60 |  4.887 |  4.302 |  4.379 |  3.218 |  1.943 |  1.929
60$\leq$DTM<90 |  4.456 |  3.805 |  3.838 |  3.148 |  1.530 |  1.595
90$\leq$DTM<120 |  4.483 |  3.746 |  3.801 |  3.279 |  1.662 |  1.775
120$\leq$DTM<150 |  4.567 |  3.552 |  3.664 |  3.572 |  1.852 |  1.961
150$\leq$DTM |  4.754 |  3.501 |  3.644 |  3.974 |  1.920 |  2.038
|  |  |  |  |  |
Partitioned by the level of VIX
VIX<15 |  4.476 |  2.835 |  3.087 |  3.342 |  1.727 |  1.797
15$\leq$VIX<20 |  4.031 |  3.158 |  3.267 |  3.181 |  1.633 |  1.696
20$\leq$VIX<25 |  3.855 |  3.953 |  3.886 |  3.115 |  1.844 |  1.856
25$\leq$VIX<30 |  4.105 |  4.460 |  4.317 |  3.559 |  1.946 |  2.050
30$\leq$VIX<35 |  4.682 |  5.033 |  4.820 |  3.453 |  2.247 |  2.276
35$\leq$VIX |  9.141 |  8.472 |  8.536 |  5.084 |  2.947 |  2.984

Note: This table report the out-of-sample option pricing performance for each
model. We conduct our out-of-sample performance evaluation by splitting our
original dataset into two subsamples: the in-sample data consists of the years
before 2008 and the out-of-sample consists of the years 2008-2019, which spans
a 12-year period. The estimation for each model is done only once for the in-
sample data, and then price the out-of-sample data by the estimated
parameters. Therefore, the in-sample data is 1990-2007. We evaluate the
model’s option pricing ability through the root of mean square errors of
implied volatility (IVRMSE). We summarize the results by option moneyness,
maturity and market VIX level. Moneyness is measured by Delta computed from
the Black-Scholes model. DTM denotes the number of calendar days to maturity.
The numbers in shaded region indicate the minimum value within each row in
each time period..

####  5.5 Out-of-Sample Pricing Performance

In this section we turn to out-of-sample comparisons. The SHNG specification
with time-varying $\eta_{t}$ is obviously a more flexible model than the HNG
model, which as fewer parameters. Thus, the improved in-sample performance of
the SHNG model could potentially be driven by overfitting. It is therefore
interesting to investigate if the improved derivative pricing also holds in
out-of-sample comparisons. We therefore divide the full dataset into two
subsamples. Data from the years 1990-2007 (in-sample) are used to estimate
each of the six specifications, whereas data from the period 2008-2019 (out-
of-sample) are used to evaluate the estimated models. We proceed with the same
evaluation and compare models in terms of their average root mean square
pricing errors.

Figure 4: The autocorrelation function (ACF) for derivative pricing errors.
The ACFs for the VIX are in the upper panels and the ACFs for the average
option pricing error are in the lower panels. Left panels present results for
the HNG model and the right panels for the SHNG model.

The out-of-sample performance by the SHNG model is just as impressive as the
in-sample results. While the average RMSE tends to be larger for the HNG model
out-of-sample, the performance is similar for the SHNG model and is, in many
cases, better than the RMSE for the full sample. This is likely explained by
the time-varying volatility risk aversion being particularly beneficial during
this period that includes the Global Financial Crisis. The SHNG specifications
provide substantially more accurate derivative pricing than all the HNG
specifications. This holds for both the VIX and option prices, and it holds
uniformly across the three ways we partition option prices.

Analogous to the full sample results, the best VIX pricing is achieved with
SHNG model estimated with VIX only, and the best option pricing is achieved by
the SHNG model estimated with option prices only. While these two
specifications have similar point estimates and yield very similar path for
$\eta_{t}$, the differences are sufficiently large to influence derivative
pricing. Overall, the out-of-sample findings are consistent with the in-sample
results.

#####  5.5.1 Insight offered by Autocorrelations

Some key insight about the improvements offered by SHNG specifications can be
found in the autocorrelations of pricing errors. The substantial and
persistent time-variation in $\eta_{t}$ have implications for the HNG
specifications. By assuming $\eta$ to be constant, we should expect the HNG
will tend to make positive pricing errors when $\eta_{t}$ is small and and
negative pricing errors when $\eta_{t}$ is large. Moreover, since $\eta_{t}$
was found to be highly persistent this will induce autocorrelation in the
pricing errors of the HNG model. We present autocorrelation functions (ACF)
for pricing errors in Figure 4 that displays the autocorrelation functions
(ACF) for HNG pricing errors (left panels) and SHNG pricing errors (right
panels). The upper panels are the ACFs for the VIX pricing errors and the
lower panels are the ACFs for the average option pricing error. The model in
the the upper panels were estimated with the VIX and those in the lower panel
were estimated with option prices. Figure 4 demonstrates that a constant
$\eta$ does indeed induce a high degree of predictability in the pricing
errors. The ACF for the pricing errors in the SHNG model are substantially
closer to zero. However, for option prices, we do find the first
autocorrelation, at about 0.2, to be significantly different from zero. The
small autocorrelations observed for the SHNG specifications is testament to
the principle on which the SHNG is build. The constant adjustments of
$\eta_{t}$, as dictated by the local first order conditions, is the main
reason that the SHNG model improves derivative pricing. The score-driven
approach monitors the first order conditions and constantly makes adjustment
that are proportional to the signal, $s_{t}$. As a consequence, the first
order condition is closer to being satisfied throughout the sample, rather
than just on average over the full sample period.

Table 6: Economic Explanation for Time-Varying Variance Risk Aversion

Log Variance Risk Ratio (Monthly Average) |
---|---
|  |  |  |  |  |  |  |  |
Sentiment |  |  -0.118*** |  |  |  |  |  |  |  -0.099***
|  |  (0.018) |  |  |  |  |  |  |  (0.014)
Dispersion-UNEMP |  |  |  0.653*** |  |  |  |  |  |  0.254***
|  |  |  (0.115) |  |  |  |  |  |  (0.095)
Dispersion-GDP |  |  |  |  0.236*** |  |  |  |  |  0.079
|  |  |  |  (0.065) |  |  |  |  |  (0.053)
Economic Policy Uncertainty |  |  |  |  |  0.258*** |  |  |  |  0.203***
|  |  |  |  |  (0.044) |  |  |  |  (0.035)
Survey-based Uncertainty |  |  |  |  |  |  0.106*** |  |  |  0.065***
|  |  |  |  |  |  (0.025) |  |  |  (0.020)
Economic Uncertainty Index |  |  |  |  |  |  |  0.037*** |  |  -0.007
|  |  |  |  |  |  |  (0.008) |  |  (0.008)
Variance Risk Premium |  |  |  |  |  |  |  |  0.097*** |  0.074***
|  |  |  |  |  |  |  |  (0.008) |  (0.007)
|  |  |  |  |  |  |  |  |
Stock Market Volatility |  0.317*** |  0.318*** |  0.299*** |  0.300*** |  0.239*** |  0.297*** |  0.290*** |  0.202*** |  0.148***
|  (0.030) |  (0.028) |  (0.029) |  (0.030) |  (0.032) |  (0.030) |  (0.030) |  (0.026) |  (0.025)
Default Yield Spread |  0.162*** |  0.108** |  0.132*** |  0.143*** |  0.198*** |  0.093* |  0.080 |  0.160*** |  0.100**
|  (0.051) |  (0.048) |  (0.049) |  (0.050) |  (0.049) |  (0.052) |  (0.053) |  (0.041) |  (0.040)
Default Return Spread |  0.031*** |  0.031*** |  0.029*** |  0.028*** |  0.026*** |  0.029*** |  0.027*** |  0.024*** |  0.021***
|  (0.007) |  (0.006) |  (0.007) |  (0.007) |  (0.007) |  (0.007) |  (0.007) |  (0.006) |  (0.005)
Long Term Yield |  0.040*** |  0.046*** |  0.027*** |  0.027*** |  0.056*** |  0.021** |  0.035*** |  0.033*** |  0.033***
|  (0.009) |  (0.008) |  (0.008) |  (0.009) |  (0.009) |  (0.010) |  (0.008) |  (0.007) |  (0.008)
Dividend Price Ratio |  -0.158** |  -0.158** |  -0.118* |  -0.135* |  -0.233*** |  -0.116* |  -0.165** |  -0.174*** |  -0.179***
|  (0.070) |  (0.065) |  (0.067) |  (0.069) |  (0.068) |  (0.069) |  (0.068) |  (0.057) |  (0.052)
Term Spread |  -0.010 |  -0.021** |  -0.014 |  -0.020** |  -0.027*** |  -0.020** |  -0.025** |  0.004 |  -0.029***
|  (0.010) |  (0.010) |  (0.010) |  (0.010) |  (0.010) |  (0.010) |  (0.010) |  (0.008) |  (0.008)
|  |  |  |  |  |  |  |  |
Other Controls |  YES |  YES |  YES |  YES |  YES |  YES |  YES |  YES |  YES
|  |  |  |  |  |  |  |  |
$R^{2}$ |  0.649 |  0.693 |  0.681 |  0.663 |  0.683 |  0.667 |  0.670 |  0.770 |  0.829
$\rho(1)$ |  0.617 |  0.566 |  0.581 |  0.611 |  0.607 |  0.625 |  0.628 |  0.387 |  0.368

Note: The controlled economic variables are taken from Welch and Goyal, (2007)
and Baker and Wurgler, (2006). Unreported other controls including six macro
variables used in Baker and Wurgler, (2006) to control the macroeconomic
conditions and five variables from Welch and Goyal, (2007). The former
includes monthly growth rate of industrial production index, real durables
consumption, real nondurables consumption, real services consumption,
employment and NBER recession indicator. The latter includes earnings price
ratio, book to market ratio, net equity expansion, long term rate of return,
consumer price index. Most of the results for these unreported other controls
are insignificant.

###  6 Variance Risk Aversion and Economic Fundamentals

Time variation in the volatility risk aversion is the key characteristic of
the new pricing kernel. In this section, we seek economic explanation for this
time variation by relating $\eta_{t}$ to economic fundamentals and other
empirical quantities from the asset pricing literature.

We focus on seven core measures of uncertainty, disagreement, and sentiment.
The first set of variables are related to investor disagreements and we follow
Bollerslev et al., (2018) and use two types of proxies for disagreement. The
first concerns macroeconomy fundamentals, where we use the forecast dispersion
for both unemployment rate and GDP growth in the Survey of Professional
Forecasters.131313Bollerslev et al., (2018) only use the forecast dispersion
for one-quarter-ahead unemployment rate. The second is disagreement about
economic policy, where we adopt the economic policy uncertainty index by Baker
et al., (2016). We also include the _variance risk premium_ by Carr and Wu,
(2008), which Grith et al., (2017) argues is a proxy for market
uncertainty.141414Up-to-date VRP series was downloaded from Hao Zhou’s
website: https://sites.google.com/site/haozhouspersonalhomepage/ We also
include three uncertainty measures from existing literature: The _Economic
Uncertainty Index_ by Bali et al., (2014); the _Survey-based Uncertainty_
index by Ozturk and Sheng, (2018);151515This index is constructed from the
Consensus Forecasts publication by Consensus Economics Inc. Up-to-date data
were downloaded from the website:
https://www.american.edu/cas/faculty/sheng.cfm and the _Sentiment_ variable by
Baker and Wurgler, (2006). In addition to these seven core variables, we
include the economic variables previously used in Baker and Wurgler, (2006)
and Welch and Goyal, (2007), which we label as “control variables”. Most
variables are available at a monthly frequency. We will, for this reason, take
monthly average of $\eta_{t}$ and regress the logarithmically transformed
average on the variables listed above (and subsets of these variables).

The regression results are listed in Table 6, where the seven core variables
are listed in the top of the table, followed by the six control variables that
contributed the most to explaining variation in $\eta_{t}$. There are eleven
additional control variables, labelled “Other Controls” which were largely
insignificant, and to conserve space the results for these variables are
presented in the Web-Appendix.161616We follow Welch and Goyal, (2007) in the
definitions of these variables. Thus, the yield spread is the difference
between BAA- and AAA- rated corporate bond yields, default return spread is
the difference between the return on long-term corporate bonds and government
bonds, etc. see Welch and Goyal, (2007) for details. The Web-Appendix also
includes various diagnostics to rule out that our estimates are driven by
spurious regression results.171717Spurious regression results can arise in a
context, such as this one, with highly persistent variables. However, the unit
root hypothesis is rejected for all variable in the regression, and the
regression residual have little autocorrelation.

The first column in Table 6 is the specification that only includes the
control variables, and these can explain 64.9% of the variation in the monthly
average of $\eta$. The six most significant control variables are: stock
market volatility, the dividend price ratio, and four variables related to
interest rates and credit risk.

Columns two to eight present the results for the specification that include a
single core variable, in addition to the control variables. Every core
variable is significant in these regressions and the signs of their estimated
coefficients are consistent with expectation. The sentiment variables is
estimated to have a negative coefficient, such that low sentiment is
associated with high volatility risk aversion. The following five variables
are related to uncertainty, and they all have positive coefficients. So, high
uncertainty is associated with high volatility aversion. Finally, the variance
risk premium is an empirical measure of the difference between volatility
under $\mathbb{Q}$ and $\mathbb{P}$. It is therefore not surprising that this
variables has a positive coefficient and increases the $R^{2}$ by the largest
amount. The last column is a “kitchen sink regression” that includes all
explanatory variables. Most of the variables in Table 6 continue to be
significant in this regression and the $R^{2}$ increases to nearly 83%. Two of
core variables,Forecast Dispersion of GDP Growth and the Economic Uncertainty
Index, become insignificant in the kitchen sink regression. This could be due
to collinearity, because these two core variables have the largest
correlation, (65.7%), of all core variables.

###  7 Summery

We have introduced a novel approach to option pricing by introducing time-
varying risk aversion in the pricing kernel. An elegant and tractable
framework emerges when the pricing kernel is combined with the Heston-Nandi
GARCH model. In this framework, the variance risk ratio is found to be _the_
fundamental variable. This variable characterizes all time-variation in the
pricing kernel. It is functionally tied to the variance risk aversion
parameter that define the curvature of the pricing kernel, and we showed that
the framework can generate the shapes of the pricing kernel that are seen
empirically. This framework made it possible to derive closed-form expression
for the VIX and option prices. In the most general setting where the variance
risk aversion is driven by a stochastic process, we derive an option pricing
formula that is based on a novel approximation method. The approximation
method uses the insight provided by VIX formula, which is valid with
predetermined and random variance risk ratio, to determine an adjusted
volatility level that yields the correct VIX pricing, when the random
variation in the variance risk ratio is neglected. When applied to option
prices the resulting approximation works very well. Empirically, we have
documented that the proposed framework can greatly reduce derivative pricing
errors. The root mean square error of the derivative pricing errors is
typically reduced by more than 50% by allowing the variance risk ratio to be
time-varying. A reduction of this magnitude was found in-sample as well as
out-of-sample. We also found that that our key variable, the variance risk
ratio, is closely related to a range of well-known measures of sentiment,
disagreement, and uncertainty, as well as key economic variables.

The key idea in this paper is to model time variation in the pricing kernel
parameters, where the dynamic properties are driven by suitable first order
conditions. This idea can also be applied to other models, than the Heston-
Nandi GARCH model used in this paper. For instance, it could be used in
conjunction with other popular GARCH models, such as EGARCH, GJR-GARCH, and
NGARCH, see, Nelson, (1991), Glosten et al., (1993), and Engle and Ng, (1993).
The same idea can also be applied to modern volatility models that utilize
realized measures of volatility, such as Realized GARCH models by Hansen et
al., (2012) and Hansen and Huang, (2016), the GARV model by Christoffersen et
al., (2014) , and the LHARG, see Majewski et al., (2015). Closed-form option
pricing formula will likely be unattainable once time-varying variance risk
aversion is introduced to these model. However approximation methods, such as
those by Duan et al., (1999), may be used instead.

In this paper, we have introduced time-variation to the pricing kernel, while
preserving the model structure under the physical measure. An implication is
that the equity risk parameter, $\lambda$, which is a parameter in the Heston-
Nandi GARCH model, is constant. So the additional time-variation we have
introduced to the pricing kernel, get embedded in the compensation for
variance risk. However, the score-driven approach can also be applied to
$\lambda$, as well as other parameters in the physical part of the model. We
leave this for future work.

### References

  * Aït-Sahalia and Lo, (2000) Aït-Sahalia, Y. and Lo, A. W. (2000).  Nonparametric risk management and implied risk aversion.  Journal of Econometrics, 94:9–51.
  * Andries et al., (2018) Andries, M., Eisenbach, T. M., and Schmalz, M. C. (2018).  Horizon-dependent risk aversion and the timing and pricing of uncertainty.  FRB of New York Staff Report, (703).
  * Arrow and Debreu, (1954) Arrow, K. J. and Debreu, G. (1954).  Existence of an equilibrium for a competitive economy.  Econometrica, 22:265–290.
  * Baker and Wurgler, (2006) Baker, M. and Wurgler, J. (2006).  Investor sentiment and the cross-section of stock returns.  The journal of Finance, 61:1645–1680.
  * Baker et al., (2016) Baker, S. R., Bloom, N., and Davis, S. J. (2016).  Measuring economic policy uncertainty.  The quarterly journal of economics, 131:1593–1636.
  * Bakshi et al., (1997) Bakshi, G., Cao, C., and Chen, Z. (1997).  Empirical performance of alternative option pricing models.  The Journal of Finance, 52:2003–2049.
  * Bakshi and Madan, (2008) Bakshi, G. and Madan, D. (2008).  Investor heterogeneity and the non-monotonicity of the aggregate marginal rate of substitution in the market index.  Unpublished working paper, University of Maryland.
  * Bakshi et al., (2010) Bakshi, G., Madan, D., and Panayotov, G. (2010).  Returns of claims on the upside and the viability of u-shaped pricing kernels.  Journal of Financial Economics, 97:130–154.
  * Bali et al., (2014) Bali, T. G., Brown, S. J., and Caglayan, M. O. (2014).  Macroeconomic risk and hedge fund returns.  Journal of Financial Economics, 114:1–19.
  * Bali and Zhou, (2016) Bali, T. G. and Zhou, H. (2016).  Risk, uncertainty, and expected returns.  Journal of Financial and Quantitative Analysis, 51:707–735.
  * Barone-Adesi et al., (2008) Barone-Adesi, G., Engle, R. F., and Mancini, L. (2008).  A GARCH option pricing model with filtered historical simulation.  Review of Financial Studies, 21:1223–1258.
  * Beare and Schmidt, (2016) Beare, B. K. and Schmidt, L. D. (2016).  An empirical test of pricing kernel monotonicity.  Journal of Applied Econometrics, 31:338–356.
  * Bekaert et al., (2020) Bekaert, G., Engstrom, E. C., and Xu, N. R. (2020).  The time variation in risk appetite and uncertainty.  NBER.
  * Blasques et al., (2015) Blasques, F., Koopman, S. J., and Lucas, A. (2015).  Information-theoretic optimality of observation-driven time series models for continuous responses.  Biometrika, 102:325–343.
  * Bollerslev, (1986) Bollerslev, T. (1986).  Generalized autoregressive heteroskedasticity.  Journal of Econometrics, 31:307–327.
  * Bollerslev et al., (2011) Bollerslev, T., Gibson, M., and Zhou, H. (2011).  Dynamic estimation of volatility risk premia and investor risk aversion from option-implied and realized volatilities.  Journal of Econometrics, 160:235–245.
  * Bollerslev et al., (2018) Bollerslev, T., Li, J., and Xue, Y. (2018).  Volume, volatility, and public news announcements.  The Review of Economic Studies, 85:2005–2041.
  * Bollerslev et al., (2009) Bollerslev, T., Tauchen, G., and Zhou, H. (2009).  Expected Stock Returns and Variance Risk Premia.  The Review of Financial Studies, 22:4463–4492.
  * Broadie et al., (2007) Broadie, M., Chernov, M., and Johannes, M. (2007).  Model specification and risk premia: Evidence from futures options.  The Journal of Finance, 62:1453–1490.
  * Brown and Jackwerth, (2012) Brown, D. P. and Jackwerth, J. C. (2012).  The pricing kernel puzzle: Reconciling index option data and economic theory.  In Derivative Securities Pricing and Modelling, pages 155–183. Emerald Group Publishing Limited.
  * Campbell and Cochrane, (1999) Campbell, J. Y. and Cochrane, J. H. (1999).  By force of habit: A consumption-based explanation of aggregate stock market behavior.  Journal of political Economy, 107:205–251.
  * Carr and Wu, (2008) Carr, P. and Wu, L. (2008).  Variance Risk Premiums.  The Review of Financial Studies, 22:1311–1341.
  * Chabi-Yo, (2012) Chabi-Yo, F. (2012).  Pricing kernels with stochastic skewness and volatility risk.  Management Science, 58:624–640.
  * Chabi-Yo et al., (2008) Chabi-Yo, F., Garcia, R., and Renault, E. (2008).  State dependence can explain the risk aversion puzzle.  Review of Financial Studies, 21:973–1011.
  * Chabi-Yo and Song, (2013) Chabi-Yo, F. and Song, Z. (2013).  Recovering the probability weights of tail events with volatility risk from option prices.  Technical report, Working Paper, Ohio State University.
  * Chen and Ghysels, (2011) Chen, X. and Ghysels, E. (2011).  News – good or bad – and its impact on volatility predictions over multiple horizons.  Review of Financial Studies, 24:46–81.
  * Chernov and Ghysels, (2000) Chernov, M. and Ghysels, E. (2000).  A study towards a unified approach to the joint estimation of objective and risk neutral measures for the purpose of options valuation.  Journal of Financial Economics, 56:407 – 458.
  * Christoffersen et al., (2014) Christoffersen, P., Feunou, B., Jacobs, K., and Meddahi, N. (2014).  The economic value of realized volatility: Using high-frequency returns for option valuation.  Journal of Financial and Quantitative Analysis, 49:663–697.
  * Christoffersen et al., (2013) Christoffersen, P., Heston, S. L., and Jacobs, K. (2013).  Capturing option anomalies with a variance-dependent pricing kernel.  Review of Financial Studies, 26:1963–2006.
  * Corsi et al., (2013) Corsi, F., Fusari, N., and Vecchia, D. L. (2013).  Realizing smiles: Options pricing with realized volatility.  Journal of Financial Economics, 107:284–304.
  * Cox et al., (1985) Cox, J. C., Ingersoll, J. E., and Ross, S. A. (1985).  A theory of the term structure of interest rates.  Econometrica, 53:385–408.
  * Creal et al., (2013) Creal, D. D., Koopman, S. J., and Lucas, A. (2013).  Generalized autoregressive score models with applications.  Journal of Applied Econometrics, 28:777–795.
  * Cuesdeanu and Jackwerth, (2018) Cuesdeanu, H. and Jackwerth, J. C. (2018).  The pricing kernel puzzle in forward looking data.  Review of Derivatives Research, 21:253–276.
  * Duan, (1995) Duan, J.-C. (1995).  The GARCH option pricing model.  Mathematical finance, 5:13–32.
  * Duan et al., (1999) Duan, J.-C., Gauthier, G., and Simonato, J. (1999).  An analytical approximation for the GARCH option pricing model.  Journal of Computational Finance, 2:75–116.
  * Eisenbach and Schmalz, (2016) Eisenbach, T. M. and Schmalz, M. C. (2016).  Anxiety in the face of risk.  Journal of Financial Economics, 121:414–426.
  * Engle, (1982) Engle, R. F. (1982).  Autoregressive conditional heteroskedasticity with estimates of the variance of U.K. inflation.  Econometrica, 45:987–1007.
  * Engle and Ng, (1993) Engle, R. F. and Ng, V. K. (1993).  Measuring and Testing the Impact of News on Volatility.  Journal of Finance, 48:1749–78.
  * Feller, (1951) Feller, W. (1951).  Two singular diffusion problems.  Annals of mathematics, pages 173–182.
  * Glosten et al., (1993) Glosten, L. R., Jagannathan, R., and Runkle, D. E. (1993).  On the Relation between the Expected Value and the Volatility of the Nominal Excess Return on Stocks.  Journal of Finance, 48:1779–1801.
  * González et al., (2018) González, M., Nave, J., and Rubio, G. (2018).  Macroeconomic determinants of stock market betas.  Journal of Empirical Finance, 45:26–44.
  * Grith et al., (2017) Grith, M., Härdle, W. K., and Krätschmer, V. (2017).  Reference-dependent preferences and the empirical pricing kernel puzzle.  Review of Finance, 21:269–298.
  * Han, (2008) Han, B. (2008).  Investor sentiment and option prices.  The Review of Financial Studies, 21:387–414.
  * Hansen and Huang, (2016) Hansen, P. R. and Huang, Z. (2016).  Exponential GARCH modeling with realized measures of volatility.  Journal of Business & Economic Statistics, 34:269–287.
  * Hansen et al., (2012) Hansen, P. R., Huang, Z., and Shek, H. (2012).  Realized GARCH: A joint model of returns and realized measures of volatility.  Journal of Applied Econometrics, 27:877–906.
  * Hao and Zhang, (2013) Hao, J. and Zhang, J. E. (2013).  GARCH option pricing models, the CBOE VIX, and variance risk premium.  Journal of Financial Econometrics, 11:556–580.
  * Heston, (1993) Heston, S. L. (1993).  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  Review of Financial Studies, 6:327–343.
  * Heston and Nandi, (2000) Heston, S. L. and Nandi, S. (2000).  A closed-form GARCH option valuation model.  Review of Financial Studies, 13:585–625.
  * Huang et al., (2017) Huang, Z., Wang, T., and Hansen, P. R. (2017).  Option pricing with the realized GARCH model: An analytical approximation approach.  Journal of Futures Markets, 37:328–358.
  * Jackwerth, (2000) Jackwerth, J. C. (2000).  Recovering risk aversion from option prices and realized returns.  The Review of Financial Studies, 13:433–451.
  * Kiesel and Rahe, (2017) Kiesel, R. and Rahe, F. (2017).  Option pricing under time-varying risk-aversion with applications to risk forecasting.  Journal of Banking & Finance, 76:120–138.
  * Koopman et al., (2016) Koopman, S. J., Lucas, A., and Scharth, M. (2016).  Predicting time-varying parameters with parameter-driven and observation-driven models.  Review of Economics and Statistics, 98:97–110.
  * Li, (2007) Li, G. (2007).  Time-varying risk aversion and asset prices.  Journal of Banking & Finance, 31:243–257.
  * Lucas, (1978) Lucas, R. E. (1978).  Asset prices in an exchange economy.  Econometrica, 46:1429–1445.
  * Majewski et al., (2015) Majewski, A. A., Bormetti, G., and Corsi, F. (2015).  Smile from the past: A general option pricing framework with multiple volatility and leverage components.  Journal of Econometrics, 187:521–531.
  * Monahan, (2001) Monahan, J. F. (2001).  Numerical Methods of Statistics.  Cambridge University Press, Cambridge.
  * Nelson, (1991) Nelson, D. B. (1991).  Conditional heteroskedasticity in asset returns: A new approach.  Econometrica, 59:347–370.
  * Ozturk and Sheng, (2018) Ozturk, E. O. and Sheng, X. S. (2018).  Measuring global and country-specific uncertainty.  Journal of International Money and Finance, 88:276–295.
  * Polkovnichenko and Zhao, (2013) Polkovnichenko, V. and Zhao, F. (2013).  Probability weighting functions implied in options prices.  Journal of Financial Economics, 107:580–609.
  * Rosenberg and Engle, (2002) Rosenberg, J. V. and Engle, R. F. (2002).  Empirical pricing kernels.  Journal of Financial Economics, 64:341–372.
  * Rubinstein, (1976) Rubinstein, M. (1976).  The valuation of uncertain income streams and the pricing of options.  Bell Journal of Economics, 7:407–425.
  * Shefrin, (2001) Shefrin, H. (2001).  On kernels and sentiment.  Available at SSRN 288258.
  * Shefrin, (2008) Shefrin, H. (2008).  Risk and return in behavioral SDF-based asset pricing models.  Journal of Investment Management, 6:1–18.
  * Song and Xiu, (2016) Song, Z. and Xiu, D. (2016).  A tale of two option markets: Pricing kernels and volatility risk.  Journal of Econometrics, 190:176–196.
  * Welch and Goyal, (2007) Welch, I. and Goyal, A. (2007).  A comprehensive look at the empirical performance of equity premium prediction.  The Review of Financial Studies, 21:1455–1508.

###  Appendix A Appendix of Proofs

Proof of Theorem 1. The pricing kernel with time-varying risk parameters takes
the form

$M_{t+1,t}=\frac{\exp[\phi_{t}R_{t+1}+\xi_{t}h_{t+2}]}{\mathbb{E}_{t}^{\mathbb{P}}[\exp(\phi_{t}R_{t+1}+\xi_{t}h_{t+2})]},$

where the parameters $\phi_{t}$ and $\xi_{t}$ are
$\mathcal{G}_{t}$-measurable. The pricing kernel implies some dynamic
properties in risk neutral measure, which much satisfy the no-arbitrage
conditions:$\mathbb{E}_{t}^{\mathbb{P}}[M_{t+1,t}]=1$ and
$\mathbb{E}_{t}^{\mathbb{P}}[M_{t+1,t}\exp(R_{t+1})]=\exp(r)$. The first is
trivial, whereas the second condition reads

$\frac{\mathbb{E}_{t}^{\mathbb{P}}[\exp\\{(\phi_{t}+1)R_{t+1}+\xi_{t}h_{t+2}\\}]}{\mathbb{E}_{t}^{\mathbb{P}}[\exp(\phi_{t}R_{t+1}+\xi_{t}h_{t+2})]}=\exp(r).$
(A.1)

The moment generating function of $(R_{t+1},h_{t+2})$ under $\mathbb{P}$ is

$\displaystyle\text{$\mathbb{E}_{t}^{\mathbb{P}}$}\left[\exp\left(v_{1}R_{t+1}+v_{2}h_{t+2}\right)\right]$
$\displaystyle=$
$\displaystyle\exp\left[v_{1}r+v_{2}\omega+[v_{1}(\lambda-\tfrac{1}{2})+v_{2}(\beta+\alpha\gamma^{2})]h_{t+1}\right]$
$\displaystyle\times\text{$\mathbb{E}_{t}^{\mathbb{P}}$}\left[v_{2}\alpha
z_{t+1}^{2}+(v_{1}-2v_{2}\alpha\gamma)\sqrt{h_{t+1}}z_{t+1}\right]$
$\displaystyle=$
$\displaystyle\exp\left[v_{1}r+v_{2}\omega+[v_{1}(\lambda-\tfrac{1}{2})+v_{2}(\beta+\alpha\gamma^{2})]h_{t+1}\right.$
$\displaystyle\left.-\frac{1}{2}\log\left(1-2v_{2}\alpha\right)+\tfrac{(v_{1}-2v_{2}\alpha\gamma)^{2}}{2\left(1-2v_{2}\alpha\right)}h_{t+1}\right]$
$\displaystyle=$
$\displaystyle\exp\left[v_{1}r+v_{2}\omega-\tfrac{1}{2}\log\left(1-2v_{2}\alpha\right)\right.$
$\displaystyle\left.+\left(v_{1}(\lambda-\tfrac{1}{2})+v_{2}(\beta+\alpha\gamma^{2})+\tfrac{(v_{1}-2v_{2}\alpha\gamma)^{2}}{2\left(1-2v_{2}\alpha\right)}\right)h_{t+1}\right].$

Substituting $(\phi_{t+1}+1,\xi_{t})$ and $(\phi_{t+1},\xi_{t})$ for
$(v_{1},v_{2})$ and taking their ratio, reveals that (A.1) simplifies to
$[\lambda-\tfrac{1}{2}+\frac{1+2\phi_{t}-4\xi_{t}\alpha\gamma}{2(1-2\xi_{t}\alpha)}]h_{t+1}=0.$
Next, solving for $\lambda$ yields

$\lambda=\eta_{t}\left(\gamma-\phi_{t}-\tfrac{1}{2}\right)-\gamma+\tfrac{1}{2},\quad\text{where}\quad\eta_{t}\equiv\frac{1}{1-2\alpha\xi_{t}},$
(A.2)

and the results stated in part $(i)$ are established.

The risk-neutral dynamic is deduced from the moment generating function of
$z_{t+1}$ under the $\mathbb{\mathbb{Q}}$ measure is:

$\displaystyle\text{$\mathbb{E}_{t}^{\mathbb{Q}}$}\left[\exp\left(v_{1}z_{t+1}\right)\right]$
$\displaystyle=\exp\left[\left(\tfrac{\phi_{t}-2\xi_{t}\alpha\gamma}{1-2\xi_{t}\alpha}\sqrt{h_{t+1}}\right)v_{1}+\tfrac{1}{2\left(1-2\xi_{t}\alpha\right)}v_{1}^{2}\right]$
$\displaystyle=\exp\left[\left[-\left(\lambda+\tfrac{1}{2}\eta_{t}-\tfrac{1}{2}\right)\sqrt{h_{t+1}}\right]v_{1}+\tfrac{1}{2}\eta_{t}v_{1}^{2}\right]$
(A.3)

The last equality comes from (A.2). The expression (A.3) shows that $z_{t+1}$
is also normally distributed under $\mathbb{\mathbb{Q}}$, but with a different
mean and a different variance (than under $\mathbb{P}$). So we define

$z_{t+1}^{*}=\frac{z_{t+1}-\mathbb{E}_{t}^{\mathbb{Q}}(z_{t+1})}{\sqrt{\mathrm{var}^{\mathbb{Q}}(z_{t+1})}}=\frac{1}{\sqrt{\eta_{t}}}\left(z_{t+1}+\left(\lambda+\frac{1}{2}\eta_{t}-\tfrac{1}{2}\right)\sqrt{h_{t+1}}\right),\quad
t=1,2,\ldots,$ (A.4)

so that $z_{t}^{\ast}\sim iidN(0,1)$ under $\mathbb{\mathbb{Q}}$. Substituting
the risk-neutralized quantities, $z_{t+1}^{*}$ and
$h_{t+1}^{*}=\eta_{t}h_{t+1}$ into (1), we arrive at the expressions stated in
part $(ii)$. Finally, part $(iii)$ is an implication of the MGF for $z_{t+1}$
under $\mathbb{\mathbb{Q}}$, see (A.3). $\square$

Proof of Theorem 2. We rewrite the volatility dynamic under $\mathbb{Q}$
measure as

$h_{t+1}^{*}=(\omega_{t}^{*}+\alpha_{t}^{*})+(\beta_{t}^{*}+\alpha_{t}^{*}\gamma_{t-1}^{*2})h_{t}^{*}+\nu_{t},$

where
$\nu_{t}=\alpha_{t}^{*}((z_{t}^{*})^{2}-1-2\gamma_{t-1}^{*}z_{t}^{*}\sqrt{h_{t}^{*}})$
is a zero-mean variation condition on $\mathcal{F}_{t-1}$. We denote
$A_{t}\equiv\omega_{t}+\alpha_{t}$ and
$B_{t}\equiv\beta_{t}^{\ast}+\alpha_{t}\gamma_{t-1}^{\ast 2}$, then we have

$\displaystyle h_{t+k}^{*}$ $\displaystyle=$ $\displaystyle
A_{t+k-1}+B_{t+k-1}h_{t+k-1}^{*}+\nu_{t+k-1}$ $\displaystyle=$ $\displaystyle
A_{t+k-1}+B_{t+k-1}A_{t+k-2}+B_{t+k-1}B_{t+k-2}h_{t+k-2}^{*}+B_{t+k-2}\nu_{t+k-2}+\nu_{t+k-1}$
$\displaystyle=$ $\displaystyle
A_{t+k-1}+B_{t+k-1}A_{t+k-2}+\cdots+B_{t+k-1}\cdots
B_{t+2}A_{t+1}+B_{t+k-1}\cdots B_{t+1}h_{t+1}^{*}$
$\displaystyle+\sum_{i=1}^{t+k-1}B_{t+k-2}\cdots B_{i}\nu_{i}.$

In order to get its analytical formula, we replace $B_{t}$ with
$\tilde{B}_{t}\equiv\beta_{t}^{\ast}+\alpha_{t}\tilde{\gamma}_{t-1}^{\ast 2}$
where $\tilde{\gamma}_{t}^{\ast}=\tfrac{1}{\eta_{t}}(\gamma+\lambda)$.
Empirically $\tilde{\gamma}_{t}^{\ast}$ is indistinguishable from
$\gamma_{t}^{*}=\tfrac{1}{\eta_{t}}(\gamma+\lambda+\tfrac{1}{2}(\eta_{t}-1))$
because the magnitude of $\gamma+\lambda$ is much larger than
$\frac{1}{2}(\eta_{t}-1)$. Empirically we have $\gamma+\lambda\approx 280$
whereas $\frac{1}{2}(\eta_{t}-1)$ ranges between $-\tfrac{1}{4}$ and 1.5, with
an average value of about 0.1.

It follows that181818Recall that $\alpha_{t}^{\ast}=\alpha\eta_{t}\eta_{t-1}$
and $\beta_{t}^{*}=\beta\tfrac{\eta_{t}}{\eta_{t-1}}$.

$\tilde{B}_{t}=\beta_{t}^{\ast}+\alpha_{t}^{\ast}\gamma_{t-1}^{\ast
2}=\beta\frac{\eta_{t}}{\eta_{t-1}}+\alpha\eta_{t}\eta_{t-1}(\tfrac{\gamma+\lambda}{\eta_{t-1}})^{2}=\tilde{\beta}\frac{\eta_{t}}{\eta_{t-1}},$

where $\tilde{\beta}\equiv\beta+\alpha(\gamma+\lambda)^{2}$, so that
$\tilde{B}_{t+k}\tilde{B}_{t+k-1}\cdots\tilde{B}_{t+i}=\tilde{\beta}^{k-i+1}\frac{\eta_{t+k}}{\eta_{t+i-1}}$,
and

$\tilde{B}_{t+k}\cdots\tilde{B}_{t+i+1}A_{t+i}=\tilde{\beta}^{k-i}\frac{\eta_{t+k}}{\eta_{t+i}}(\omega\eta_{t+i}+\alpha\eta_{t+i}\eta_{t+i-1})=\tilde{\beta}^{k-i}(\omega\eta_{t+k}+\alpha\eta_{t+k}\eta_{t+i-1}).$

Using these results, we have

$\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})=\sum_{i=2}^{k}\tilde{\beta}^{k-i}\mathbb{E}_{t}^{\mathbb{Q}}(\omega\eta_{t+k-1}+\alpha\eta_{t+k-1}\eta_{t+i-2})+\tilde{\beta}^{k-1}\mathbb{E}_{t}^{\mathbb{Q}}(\frac{\eta_{t+k-1}}{\eta_{t}})h_{t+1}^{*}+\delta_{k},$

where $\delta_{k}$ is the term that results from substituting $\tilde{B}_{t}$
for $B_{t}$.

The $M$-days ahead VIX can be calculated as the annualized arithmetic average
of the expected daily variance over the following month under the risk-neutral
measure, which is

$\mathrm{VIX}_{t}=A\times\sqrt{\frac{1}{M}\sum_{k=1}^{M}\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})}.$

Where $A=100\sqrt{252}$ is the annualizing factor. The model implied VIX
formula is therefore given by,

$\sum_{k=1}^{M}\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})=\underset{(*)}{\underbrace{\sum_{k=2}^{M}\sum_{i=2}^{k}\tilde{\beta}^{k-i}\mathbb{E}_{t}^{\mathbb{Q}}(\omega\eta_{t+k-1}+\alpha\eta_{t+k-1}\eta_{t+i-2})}}+\underset{(**)}{\underbrace{\sum_{k=1}^{M}\tilde{\beta}^{k-1}\mathbb{E}_{t}^{\mathbb{Q}}(\frac{\eta_{t+k-1}}{\eta_{t}})}h_{t+1}^{*}}+\Delta,$
(A.5)

where $\Delta=\sum\delta_{k}$. We provide bounds for the approximation error
term $\Delta$ in Lemma A.1, and show that this term is negligible in practice.

From the AR(1) structure,
$\eta_{t}=(1-\theta)\zeta+\theta\eta_{t-1}+\epsilon_{t}$ with $|\theta|<1$, so
that $\mathbb{E}\eta_{t}=\zeta$, it follows that
$\mathbb{E}_{t}^{\mathbb{Q}}(\eta_{t+k}-\zeta)=\theta^{k}(\eta_{t}-\zeta)$,
and

$\mathbb{E}_{t}^{\mathbb{Q}}[(\eta_{t+k}-\zeta)^{2}]=(\theta^{2})^{k}(\eta_{t}-\zeta)^{2}+\tfrac{1-(\theta^{2})^{k}}{1-\theta^{2}}\sigma^{2},$

where $\sigma^{2}=\mathbb{E}_{t}^{\mathbb{Q}}(\varepsilon_{t+1}^{2})$.

Set $g(x)\equiv\sum_{i=1}^{M-1}x^{i}$, so that $g(x)=x(1-x^{M-1})/(1-x)$ if
$x\neq 1$. Then tedious algebraic calculus shows that the first term in (A.5),
takes the form

$(*)=\sum_{k=2}^{M}\sum_{i=2}^{k}\tilde{\beta}^{k-i}\mathbb{E}_{t}^{\mathbb{Q}}(\omega\eta_{t+k-1}+\alpha\eta_{t+k-1}\eta_{t+i-2})=a_{1}(\eta_{t},M)+a_{2}(M)\sigma^{2},$

where

$\displaystyle a_{1}(\eta_{t},M)$ $\displaystyle=$
$\displaystyle(\omega+\alpha\zeta)(\eta_{t}-\zeta)\tfrac{g(\theta)-g(\tilde{\beta}\theta)}{1-\tilde{\beta}}+\bar{\eta}(\omega+\alpha\zeta)\tfrac{g(1)-g(\tilde{\beta})}{1-\tilde{\beta}}$
$\displaystyle+\alpha\bar{\eta}(\eta_{t}-\zeta)\tfrac{g(\theta)-g(\tilde{\beta})}{\theta-\tilde{\beta}}+\alpha\left(\eta_{t}-\zeta\right)^{2}\tfrac{g(\theta^{2})-g(\tilde{\beta}\theta)}{\theta-\tilde{\beta}},$

and
$a_{2}(M)=\frac{\alpha}{1-\theta^{2}}\left(\frac{g(1)-g(\tilde{\beta}\theta)}{\theta^{-1}-\tilde{\beta}}-\frac{g(\theta^{2})-g(\tilde{\beta}\theta)}{\theta-\tilde{\beta}}\right)$.
Note that $a_{2}(M)$ does depend on the term structure for volatility risk
aversion.

For the second term in (A.5) we find,

$(**)=a_{3}(\eta_{t},M)=\sum_{k=1}^{M}\tilde{\beta}^{k-1}\mathbb{E}_{t}^{\mathbb{Q}}\left(\frac{\eta_{t+k-1}}{\eta_{t}}\right)=\tfrac{1-(\tilde{\beta}\theta)^{M}}{1-\tilde{\beta}\theta}+\left(\tfrac{1-\tilde{\beta}^{M}}{1-\tilde{\beta}}-\tfrac{1-(\tilde{\beta}\theta)^{M}}{1-\tilde{\beta}\theta}\right)\frac{\zeta}{\eta_{t}}.$

$\square$

######  Lemma A.1 (Bound for $\Delta$).

Suppose that $\eta_{t}\in[\eta_{L},\eta_{H}]$ is bounded where
$\eta_{L},\eta_{H}>0$. Then $\Delta\leq\delta(\tilde{\beta}_{{\rm
max}})-\delta(\tilde{\beta}_{{\rm min}})$, where

$\delta(x)=A\times\sqrt{\frac{1}{M}\left[\sum_{k=2}^{M}\sum_{i=2}^{k}x^{k-i}\mathbb{E}_{t}^{\mathbb{Q}}(\omega\eta_{t+k-1}+\alpha\eta_{t+k-1}\eta_{t+i-2})+\sum_{k=1}^{M}x^{k-1}\mathbb{E}_{t}^{\mathbb{Q}}(\frac{\eta_{t+k-1}}{\eta_{t}})h_{t+1}^{*}\right]},$

and $\tilde{\beta}_{{\rm
min}}=\beta+\alpha(\gamma+\lambda+\tfrac{1}{2}(\eta_{L}-1))^{2}$ and
$\tilde{\beta}_{{\rm
max}}=\beta+\alpha(\gamma+\lambda+\tfrac{1}{2}(\eta_{H}-1))^{2}$.

###### Proof.

We have $B_{t}=\beta_{t}^{\ast}+\alpha_{t}^{\ast}\gamma_{t-1}^{\ast
2}\in[\tilde{\beta}_{{\rm min}},\tilde{\beta}_{{\rm
max}}]\frac{\eta_{t}}{\eta_{t-1}}$. Since $\delta(x)$ is a monotone function
in $x$, there exists an $\bar{\beta}\in[\tilde{\beta}_{{\rm
min}},\tilde{\beta}_{{\rm max}}]$, such that
$\delta(\bar{\beta})=A\times\sqrt{\tfrac{1}{M}\sum_{k=1}^{M}\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})}$
(i.e. true value of VIX), and $\Delta^{{\rm
vix}}=|\delta(\bar{\beta})-\delta(\tilde{\beta})|$, then we have
$\Delta\leq\delta(\tilde{\beta}_{{\rm max}})-\delta(\tilde{\beta}_{{\rm
min}})$. Note that $\delta(x)$ can be expressed in terms of
$a_{1}(\eta_{t},M)$, $a_{2}(M)$, and $a_{3}(\eta_{t},M)$. ∎

######  Lemma A.2.

In the predetermined case,
$\eta_{\tau}=\mathbb{E}_{t}^{\mathbb{Q}}[\eta_{\tau}]$ for
$\tau=t+1,\ldots,T$. The moment generating function for cumulative returns has
the exponentially affine form:

$\displaystyle g_{t,T}^{*}(s)$
$\displaystyle=\mathbb{E}_{t}^{\mathbb{Q}}(\exp(s\sum_{\tau=t+1}^{T}R_{\tau}))=\exp\left(\text{$\mathcal{A}$}_{T}(s,M)+\mathcal{B}_{T}(s,M)h_{t+1}^{*}\right),$

where $\text{$\mathcal{A}$}_{T}(s,m)$ and $\mathcal{B}_{T}(s,m)$ are
recursively given by

$\displaystyle\text{$\mathcal{A}$}_{T}(s,m+1)$
$\displaystyle=\text{$\mathcal{A}$}_{T}(s,m)+sr+\mathcal{B}_{T}(m,s)\omega_{T-m}^{*}-\tfrac{1}{2}\log(1-2\alpha_{T-m}^{*}\mathcal{B}_{T}(m,s)),$
(A.6) $\displaystyle\mathcal{B}_{T}(s,m+1)$
$\displaystyle=s(\gamma_{T-m-1}^{*}-\tfrac{1}{2})-\tfrac{1}{2}\gamma_{T-m-1}^{*2}+\beta_{T-m}^{*}\mathcal{B}_{T}(s,m)+\tfrac{(s-\gamma_{T-m-1}^{*})^{2}}{2(1-2\alpha_{T-m}^{*}\mathcal{B}_{T}(s,m))},$
(A.7)

with initial values $\mathcal{A}_{T}(s,1)=sr$ and
$\mathcal{B}_{T}(s,1)=\tfrac{1}{2}\left(s^{2}-s\right)$, where all terms in
(A.6) and (A.7) are predetermined.

###### Proof.

For later use, we first note that the MGF for $(R_{t+1},h_{t+2}^{*})$ is given
by:

$\displaystyle\mathbb{E}_{t}^{\mathbb{Q}}(\exp(sR_{t+1}+uh_{t+2}^{*}))$
$\displaystyle=$
$\displaystyle\exp\Biggl{[}sr+u\omega_{t+1}^{*}-\tfrac{1}{2}\log(1-2u\alpha_{t+1}^{*})$
(A.8) $\displaystyle+$
$\displaystyle\left.\left(s(\gamma_{t}^{*}-\tfrac{1}{2})+u\beta_{t+1}^{*}-\tfrac{1}{2}\gamma_{t}^{*2}+\tfrac{(s-\gamma_{t}^{*})^{2}}{2(1-2u\alpha_{t+1}^{*})}\right)h_{t+1}^{*}\right],$

where we used Theorem 1. Next, we establish the exponentially affine form,
$\mathbb{E}_{t}^{\mathbb{Q}}(\exp(s\sum_{\tau=t+1}^{T}R_{\tau}))=\exp(\text{$\mathcal{A}$}_{T}(s,M)+\mathcal{B}_{T}(s,M)h_{t+1}^{*})$
where $M=T-t$. This is proven by backwards induction. First observe that for
$t=T-1$ we obviously have

$\mathbb{E}_{T-1}^{\mathbb{Q}}(\exp(sR_{T}))=\exp[sr+\tfrac{1}{2}(s^{2}-s)h_{T}^{\ast}].$

So $\mathcal{A}_{T}(s,1)=sr$ and $\mathcal{B}_{T}(s,1)=\tfrac{1}{2}(s^{2}-s)$
define the initial values for the recursive formulae we derive for
$\mathcal{A}_{T}$ and $\mathcal{B}_{T}$ below. Next, we establish the
induction, by showing that the exponentially affine form holds for $t=T-m$, if
it holds for $t=T-(m-1)$. Thus, consider

$\displaystyle\mathbb{E}_{T-m}^{\mathbb{Q}}(\exp(s\sum_{\tau=T-m+1}^{T}R_{\tau}))$
$\displaystyle=$
$\displaystyle\mathbb{E}_{T-m}^{\mathbb{Q}}(\exp(sR_{T-m+1}+s\sum_{\tau=T-m+2}^{T}R_{\tau}))$
$\displaystyle=$
$\displaystyle\mathbb{E}_{T-m}^{\mathbb{Q}}[\mathbb{E}_{t+1}^{\mathbb{Q}}\left\\{\exp(sR_{T-m+1}+s\sum_{\tau=T-m+2}^{T}R_{\tau})\right\\}]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{T-m}^{\mathbb{Q}}[\exp(sR_{T-m+1})\mathbb{E}_{T-m+1}^{\mathbb{Q}}\left\\{s\sum_{\tau=T-m+2}^{T}R_{\tau})\right\\}]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{T-m}^{\mathbb{Q}}(\exp(sR_{T-m+1}+\text{$\mathcal{A}$}_{T}(s,m-1)+\mathcal{B}_{T}(s,m-1)h_{T-m+2}^{*}).$

Now substitute $\mathcal{B}_{T}(s,m-1)$ for $u$ in (A.8) with $t=T-m$, then
the exponentially affine form holds if we set

$\displaystyle\text{$\mathcal{A}$}_{T}(s,m)$ $\displaystyle=$
$\displaystyle\text{$\mathcal{A}$}_{T}(s,m-1)+sr+\mathcal{B}_{T}(s,m-1)\omega_{T-m+1}^{*}-\tfrac{1}{2}\log(1-2\alpha_{T-m+1}^{*}\mathcal{B}_{T}(s,m-1))$
$\displaystyle\mathcal{B}_{T}(s,m)$ $\displaystyle=$ $\displaystyle
s(\gamma_{T-m}^{*}-\tfrac{1}{2})-\tfrac{1}{2}\gamma_{T-m}^{*2}+\beta_{T-m+1}^{*}\mathcal{B}_{T}(s,m-1)+\tfrac{(s-\gamma_{T-m}^{*})^{2}}{2(1-2\alpha_{T-m+1}^{*}\mathcal{B}_{T}(s,m-1))},$

and the result follows. ∎

Note that unlike the conventional formula for the moment generating function
in the Heston-Nandi GARCH model, the iterative terms
$\text{$\mathcal{A}$}_{T}(s,m)$ and $\mathcal{B}_{T}(m,s)$ are not only
function of days to maturity $m$, but also the maturity date $T$. This is
because key coefficients, $\omega^{\ast}$, $\alpha^{\ast}$, $\beta^{\ast}$,
and $\gamma^{\ast}$, are time-varying under $\mathbb{Q}$, as the depend on the
time-variation in $\eta_{t}$.

Proof of Theorem 3. From Christoffersen et al., (2013) we have the analytical
formula for a European call option price, and the only thing that remains is
to substitute the moment generating function for cumulative returns, that is
valid with time-varying, but predetermined variance risk ratio. The proper
expression is that for, $g_{t,T}^{*}(s)$, derived in Lemma A.2. $\square$

Proof of Proposition 1. The exponentially affine moment generating function
requires $\eta_{t}$ to be predetermined. Random variation in $\eta_{t}$ will
impact option prices, and our closed-form VIX formula, that holds with
predetermined $\eta_{t}$ and stochastic $\eta_{t}$, provides insight about the
extend to which derivative prices are influences by random variation in
$\eta_{t}$. From the VIX formula, we determine the hypothetical level of
volatility,
$\tilde{h}_{t+1}^{\ast}=h_{t+1}^{\ast}+\psi(\eta_{t},\sigma^{2},M)$, that
yields the correct VIX when the random variation is neglected, i.e.
$\sigma^{2}=0$ is assumed. Our approximate option pricing formula simply
amount to porting the adjusted volatility level, $\tilde{h}_{t+1}^{\ast}$,
into our expression for option pricing (that assumes $\sigma^{2}=0$). All
expressions are now evaluated using the predicted path for $\eta_{t}$ to
maturity, which influences many terms including the recursive formula for
$\mathcal{A}_{T}$ and $\mathcal{B}_{T}$. $\square$

###  Appendix B The Score for $\eta_{t}$

A key variable for the empirical implementation of the new model is the score,
$s_{t}$, as defined in (10). This expression requires the derivative to be
computed as well as the conditional variance of $s_{t}$, which we obtain by
numerical integration methods.

####  B.1 The Score for $\eta_{t}$

We observe $N_{t}$ derivative prices at time $t$, resulting in the vector of
pricing error, $e_{t}\in\mathbb{R}^{N_{t}}$. Thus the derivative for
$\eta_{t}$

$\displaystyle\nabla_{t}\equiv\tfrac{\partial\ell(X_{t}|\mathcal{G}_{t})}{\partial\eta_{t}}$
$\displaystyle=$ $\displaystyle-\tfrac{1}{2\sigma_{e}^{2}}\tfrac{\partial{\rm
tr}(\Omega_{N_{t}}^{-1}e_{t}e_{t}^{\prime})}{\partial\eta_{t}}$
$\displaystyle=$ $\displaystyle\tfrac{\partial{\rm
tr}(\Omega_{N_{t}}^{-1}e_{t}e_{t}^{\prime})}{\partial{\rm
vec}(e_{t}e_{t}^{\prime})^{\prime}}\tfrac{\partial{\rm
vec}(e_{t}e_{t}^{\prime})}{{\rm\partial
vech}(e_{t}e_{t}^{\prime})^{\prime}}\tfrac{{\rm\partial
vech}(e_{t}e_{t}^{\prime})}{\partial\eta_{t}}$ $\displaystyle=$
$\displaystyle{\rm
vec}(\Omega_{N_{t}}^{-1})^{\prime}D_{N_{t}}\tfrac{{\rm\partial
vech}(e_{t}e_{t}^{\prime})}{\partial\eta_{t}},$

where the expression for $D_{n}$ is given in appendix, and the elements of the
last term, are given by

$u_{ij}\equiv\tfrac{{\rm\partial
vech}(e_{i,t}e_{j,t})}{\partial\eta_{t}}=-\left[\left(X_{i,t}-X_{i,t}^{m}\right)\tfrac{\partial
X_{j,t}^{m}}{\partial\eta_{t}}+\left(X_{j,t}-X_{j,t}^{m}\right)\tfrac{\partial
X_{i,t}^{m}}{\partial\eta_{t}}\right],\quad i,j=1,\ldots,N_{t}.$

Next,

$\mathbb{E}_{t-1}^{\mathbb{Q}}[\nabla_{t}^{2}]=\frac{1}{4\sigma_{e}^{4}}{\rm
vec}(\Omega_{N_{t}}^{-1})^{\prime}D_{N_{t}}\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\tfrac{{\rm\partial
vech}(e_{i,t}e_{j,t})}{\partial\eta_{t}}\right)\left(\tfrac{{\rm\partial
vech}(e_{i,t}e_{j,t})}{\partial\eta_{t}}\right)^{\prime}\right]D_{N_{t}}^{\prime}{\rm
vec}(\Omega_{N_{t}}^{-1}).$

The terms in the conditional expectations are given by

$\mathbb{E}_{t-1}^{\mathbb{Q}}(u_{ij}u_{kl})=\sigma_{e}^{2}\left[\Theta(j,l)\mathbb{f}(i,k)+\Theta(j,k)\mathbb{f}(i,l)+\Theta(i,l)\mathbb{f}(j,k)+\Theta(i,k)\mathbb{f}(j,l)\right]$

where

$\displaystyle\Theta(i,j)$ $\displaystyle\equiv\begin{cases}1&i=j\\\
\rho&i\neq
j\end{cases},\quad\mathbb{f}(i,j)\equiv\mathbb{E}_{t-1}^{\mathbb{Q}}\left(\tfrac{\partial
X_{i,t}^{m}}{\partial\eta_{t}}\tfrac{\partial
X_{j,t}^{m}}{\partial\eta_{t}}\right).$

In the special case with a single derivative price, $N_{t}=1$, our expressions
simplify to

$\nabla_{t}=\frac{1}{\sigma_{e}^{2}}\left(X_{t}-X_{t}^{m}\right)\tfrac{\partial
X_{t}^{m}}{\partial\eta_{t}}\quad\text{and}\quad\mathbb{E}_{t-1}^{\mathbb{Q}}[\nabla_{t}^{2}]=\frac{1}{\sigma_{e}^{2}}\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\tfrac{\partial
X_{t}^{m}}{\partial\eta_{t}}\right)^{2}\right].$

To obtain the score, we need to derive $\tfrac{\partial
X_{t}^{m}}{\partial\eta_{t}}$, which must be done separately for the VIX and
option prices. This is done in Sections B.2 and B.3, respectively. For the
normalization of the score, we also need to determine $\mathbb{f}(i,j)$, which
is done by numerical integration in Section B.4.

####  B.2 Score for VIX

When $X_{i,t}$ is $VIX_{t}$, then there is only one derivative price (i.e.
$N_{t}=1$), it is the simplest case. As shown in the appendix, the analytical
formula for model implied $M$-days ahead VIX (denote as $VIX_{t}^{m}$) is:

$VIX_{t}^{m}=A\times\sqrt{\frac{1}{M}\sum_{k=1}^{M}\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})}$

$A=100\sqrt{252}$ is the annualizing factor. with

$\sum_{k=1}^{M}\mathbb{E}_{t}^{\mathbb{Q}}(h_{t+k}^{*})=\mathbb{f}_{1,t}(M)+\mathbb{f}_{2,t}(M)h_{t+1}^{*}$

where

$\mathbb{f}_{1,t}(M)\equiv\sum_{k=1}^{M}\tilde{\beta}^{k-1}\mathbb{E}_{t}^{\mathbb{Q}}(\frac{\eta_{t+k-1}}{\eta_{t}})=\tfrac{1-(\tilde{\beta}\theta)^{M}}{1-\tilde{\beta}\theta}+\left(\tfrac{1-\tilde{\beta}^{M}}{1-\tilde{\beta}}-\tfrac{1-(\tilde{\beta}\theta)^{M}}{1-\tilde{\beta}\theta}\right)\frac{\bar{\eta}}{\eta_{t}}.$
$\displaystyle\mathbb{f}_{2,t}(M)$
$\displaystyle\equiv\sum_{k=2}^{M}\sum_{i=2}^{k}\tilde{\beta}^{k-i}\mathbb{E}_{t}^{\mathbb{Q}}(\omega\eta_{t+k-1}+\alpha\eta_{t+k-1}\eta_{t+i-2})$
$\displaystyle=(\omega+\alpha\bar{\eta})(\eta_{t}-\bar{\eta})\frac{g(\theta)-g(\tilde{\beta}\theta)}{1-\tilde{\beta}}+\bar{\eta}(\omega+\alpha\bar{\eta})\frac{g(1)-g(\tilde{\beta})}{1-\tilde{\beta}}$
$\displaystyle\quad+\alpha\bar{\eta}(\eta_{t}-\bar{\eta})\frac{g(\theta)-f(\tilde{\beta})}{\theta-\tilde{\beta}}+\alpha\left(\eta_{t}-\bar{\eta}\right)^{2}\frac{g(\theta^{2})-g(\tilde{\beta}\theta)}{\theta-\tilde{\beta}}$
$\displaystyle\quad+\frac{\alpha\delta^{2}\sigma_{s}^{2}}{1-\theta^{2}}\left(\frac{g(1)-g(\tilde{\beta}\theta)}{\theta^{-1}-\tilde{\beta}}-\frac{g(\theta^{2})-g(\tilde{\beta}\theta)}{\theta-\tilde{\beta}}\right).$

with $g(x)\equiv\sum_{i=1}^{M-1}x^{i}$ and
$\tilde{\beta}\equiv\beta+\alpha(\gamma+\lambda)^{2}$. The $\sigma^{2}$ is the
variance of $s_{t}$, and equal to one by construction. Then we have

$\frac{\partial
VIX_{t}^{m}}{\partial\eta_{t}}=\frac{A}{2\sqrt{V_{t}(M)}}\frac{\partial
V_{t}(M)}{\partial\eta_{t}}=\frac{A}{2M\sqrt{V_{t}(M)}}\left[\frac{\partial\mathbb{f}_{1,t}(M)}{\partial\eta_{t}}+\frac{\partial\mathbb{f}_{2,t}(M)}{\partial\eta_{t}}h_{t+1}^{*}+\mathbb{f}_{2,t}(M)\frac{\partial
h_{t+1}^{*}}{\partial\eta_{t}}\right]$

where

$\displaystyle\frac{\partial\mathbb{f}_{1}(M,\eta_{t})}{\partial\eta_{t}}$
$\displaystyle=(\omega+\alpha\bar{\eta})\frac{g(\theta)-g(\tilde{\beta}\theta)}{1-\tilde{\beta}}+\alpha\bar{\eta}\frac{g(\theta)-g(\tilde{\beta})}{\theta-\tilde{\beta}}+2\alpha\left(\eta_{t}-\bar{\eta}\right)\frac{g(\theta^{2})-g(\tilde{\beta}\theta)}{\theta-\tilde{\beta}}$
$\displaystyle\frac{\partial\mathbb{f}_{2}(M,\eta_{t})}{\partial\eta_{t}}$
$\displaystyle=-\left(\frac{1-\tilde{\beta}^{M}}{1-\tilde{\beta}}-\frac{1-(\tilde{\beta}\theta)^{M}}{1-\tilde{\beta}\theta}\right)\frac{\bar{\eta}}{\eta_{t}^{2}}$
$\frac{\partial h_{t+1}^{*}}{\partial\eta_{t}}=h_{t+1}.$

Then

$\displaystyle\mathbb{E}_{t-1}^{\mathbb{Q}}[\nabla_{t}^{2}]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\tfrac{\mathrm{VIX}_{t}-\mathrm{VIX}_{t}^{m}}{\sigma_{e}^{2}}\right)^{2}\left(\tfrac{\partial\mathrm{VIX}_{t}^{m}}{\partial\eta_{t}}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\tfrac{\partial\mathrm{VIX}_{t}^{m}}{\partial\eta_{t}}\right)^{2}\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left.\left(\tfrac{\mathrm{VIX}_{t}-\mathrm{VIX}_{t}^{m}}{\sigma_{e}^{2}}\right)^{2}\right|\mathrm{VIX}_{t}^{m}\right]\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\tfrac{\partial\mathrm{VIX}_{t}^{m}}{\partial\eta_{t}}\right)^{2}\mathbb{E}_{t-1}^{\mathbb{P}}\left[\left.\left(\tfrac{\mathrm{VIX}_{t}-\mathrm{VIX}_{t}^{m}}{\sigma_{e}^{2}}\right)^{2}\right|\mathrm{VIX}_{t}^{m}\right]\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\sigma_{e}^{2}}\mathbb{E}_{t-1}^{\mathbb{Q}}\left[\left(\tfrac{\partial\mathrm{VIX}_{t}^{m}}{\partial\eta_{t}}\right)^{2}\right]$

####  B.3 Score for Options

So, the formula $\partial C^{m}/\partial\eta_{t}$ is

$\frac{\partial C^{m}(S_{t},K,r,T)}{\partial\eta_{t}}=S_{t}\frac{\partial
P_{1}(t)}{\partial\eta_{t}}-K\exp(-r(T-t))\frac{\partial
P_{2}(t)}{\partial\eta_{t}},$

with

$\displaystyle\frac{\partial P_{1}(t)}{\partial\eta_{t}}$
$\displaystyle=\text{$\frac{\exp(-r(T-t))}{\pi}\int_{0}^{\infty}{\rm
Re}\left[\frac{K^{-i\varphi}}{i\varphi S_{t}}\frac{\partial
g_{t,T}^{*}(i\varphi+1)}{\partial\eta_{t}}\right]d\varphi$,}$
$\displaystyle\frac{\partial P_{2}(t)}{\partial\eta_{t}}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}{\rm
Re}\left[\frac{K^{-i\varphi}}{i\varphi}\frac{\partial
g_{t,T}^{*}(i\varphi)}{\partial\eta_{t}}\right]d\varphi,$

where

$\frac{\partial
g_{t,T}^{*}(s)}{\partial\eta_{t}}=g_{t,T}^{*}(s)\left(\tfrac{\partial\text{$\mathcal{A}$}_{T}(s,M)}{\partial\eta_{t}}+\tfrac{\partial\mathcal{B}_{T}(s,M)}{\partial\eta_{t}}(h_{t+1}^{*}+\psi_{t}(M))+\mathcal{B}_{T}(s,M)(\tfrac{\partial
h_{t+1}^{*}}{\partial\eta_{t}}+\tfrac{\partial\psi_{t}(M)}{\partial\eta_{t}})\right),$

Here $\partial X/\partial\eta_{t}$ is denoted as $\dot{X}$, then
$\dot{\mathcal{A}}_{T}(s,m)$ and $\dot{\mathcal{\mathcal{B}}}_{T}(s,m)$ can be
calculated via iteration

$\displaystyle\dot{\mathcal{A}}_{T}(s,m+1)$
$\displaystyle=\dot{\mathcal{A}}_{T}(s,m)+\dot{\omega}_{T-m}^{*}\mathcal{B}_{T}(s,m)+\omega_{T-m}^{*}\dot{\mathcal{B}}_{T}(s,m)$
$\displaystyle\quad+\tfrac{\dot{\alpha}_{T-m}^{*}\mathcal{B}_{T}(s,m)+\alpha_{T-m}^{*}\dot{B}_{T}(s,m)}{1-2\alpha_{T-m}^{*}\mathcal{B}_{T}(s,m)},$
$\displaystyle\dot{\mathcal{B}}_{T}(s,m+1)$
$\displaystyle=(s-\gamma_{T-m-1}^{*})\dot{\gamma}_{T-m-1}^{*}+\beta_{T-m}^{*}\dot{\mathcal{B}}_{T}(s,m)+\dot{\beta}_{T-m}^{*}\mathcal{B}_{T}(s,m)$