Estimating Option Pricing Models Using a Characteristic Function-Based Linear State Space Representation

-1.0cm Estimating Option Pricing Models Using a Characteristic Function-Based Linear State Space Representation

{ Keywords: Options; Characteristic Function; Affine Jump-Diffusion; State Space Representation.}\\ { JEL Classification: Primary: C13; C58; G13; Secondary: C32; G01.}

Introduction

Over the past decades, explosive growth in the trading of option contracts has attracted the attention of academics and practitioners to the development and estimation of increasingly sophisticated option pricing models. The building blocks of many continuous-time option pricing models are semimartingale stochastic processes that govern the dynamics of the underlying asset. These processes are often latent with stochastic diffusive volatility as the prototypical example, as in the classical \citeA{heston1993} model. The literature also suggests the need to allow for a discontinuous jump component, both in the asset price dynamics and in its volatility process, potentially with a time-varying stochastic jump intensity.

An important econometric challenge lies in estimating the parameters of these continuous-time models and in filtering their unobserved and time-varying components, since option prices are highly nonlinear functions of the state vector. This stands in contrast to, for instance, term structure models, where bond yields can be represented as a linear function of the states, at least within the affine framework (see, e.g., \citeNP{piazzesi2010affine}, for a review of the affine term structure literature). To evaluate option prices as a function of the state vector, one typically needs to apply either Fourier-based methods or simulation-based approaches, in both cases at a substantial computational cost. This is one of the reasons why in much of the empirical research on option pricing, only a subset of the available option price data is used, such as at-the-money contracts or weekly (typically Wednesday) options data.

In this paper, we develop a new latent state filtering and parameter estimation procedure for option pricing models governed by general affine jump-diffusion processes. Our procedure leverages the linear relationship between the logarithm of an option-implied, model-free spanning formula for the conditional characteristic function of the underlying asset return on the one hand, and the state vector induced by parametric model specification on the other hand. From this relationship, we formally derive a linear state space representation, and establish the asymptotic properties of the corresponding measurement errors. Linearity of the measurement and state updating equations that make up the state space representation, with coefficient and variance matrices that are (semi-)closed-form functions of the parameters, allows us to exploit Kalman filtering techniques. The proposed estimation procedure is fast and easy to implement, circumventing the typical computational burden in conducting inference on option pricing models.

Exploiting the option-spanning formula of \citeA{carr2001optimal} for European-style payoff functions, we replicate the risk-neutral conditional characteristic function (CCF) of the underlying log-asset price at the expiration date in a completely model-independent way. In other words, we imply information about the CCF from the option prices without imposing any parametric assumptions on the underlying asset price dynamics. A similar option-spanning approach for the CCF is used by \citeA{todorov2019nonparametric} to develop an option-based nonparametric spot volatility estimator. On the other hand, a large stream of literature is devoted to parametric option pricing models belonging to the general affine jump-diffusion (AJD) family; canonical examples are \citeA{heston1993}, \citeA{DPS2000}, \citeA{pan2002}, and \citeA{bates2006maximum}.\footnote{See also, e.g., \citeA{broadie2007model}, \citeA{ait2015contagion}, \citeA{AFT2017} and the references therein.} The defining property of the AJD class is the exponential-affine joint CCF, which is available in semi-closed form. By comparing the two option pricing representations---model-free and model-implied---we can obtain a linear relation between the logarithm of the option-implied CCF and the model-dependent CCF within the affine framework.

The state vector in AJD option pricing models typically contains both observable processes and latent factors. We address the filtering of such latent factors by developing a linear state space representation for this model class. The development includes an asymptotic analysis of the measurement error components, consisting of observation, truncation and discretization errors, under a double asymptotic scheme in the moneyness dimension. The state space representation allows us to employ suitably modified Kalman filtering techniques to learn about the unobserved intrinsic components of the model and estimate the model parameters using quasi-maximum likelihood (QML). QML approaches based on Kalman filtering are often used in the affine term structure literature, where the yields themselves are linear functions of the state vector (see, e.g., \citeNP{duffee1999estimating}, \citeNP{jong2000time}, \citeNP{driessen2005default}). Besides the possibility to exploit Kalman filtering and QML estimation techniques, another advantage of our approach is that, once the model-free CCF has been obtained from the data, no further numerical option pricing methods, such as the FFT approach of \citeA{CM1999} or simulation-based methods, are needed. Therefore, our method reduces computational costs considerably relative to many existing approaches in the option pricing literature. We note that, whereas the parametric CCF is used to price options in Fourier-based methods, here we use the CCF to directly learn about the latent factors and model parameters.

We analyze the developed estimation procedure in Monte Carlo simulations based on several AJD specifications. We consider a one-factor AJD option pricing model, with the stochastic volatility and jump intensity both being affine functions of a single latent process, and a two-factor AJD model specification with an observable exogenous factor. We find good finite-sample performance in both cases, notwithstanding the challenging nature of the econometric problem.

Finally, we illustrate our new filtering and estimation approach in an empirical application to S&P 500 index options. In particular, we filter and estimate the latent volatility and jump intensity from a stochastic volatility model with co-jumps in returns and volatility. We also investigate the impact of the Covid-19 propagation rate on the stock market within this model, by embedding the associated reproduction number as an exogenous factor into the volatility and jump intensity dynamics. Our results show that while the reproduction number has only a mild effect on total diffusive volatility, it contributes substantially to the likelihood of jumps. By contrast, when we consider an Economic Policy Uncertainty index as exogenous factor, the jump intensity process is not affected, but the exogenous factor contributes significantly to diffusive volatility.

Various estimation and filtering strategies for option pricing models have been developed in the literature. These include the (penalized) nonlinear least squares methods in, for instance, \citeA{bakshi1997empirical}, \citeA{broadie2007model}, \citeA{andersen2015parametric}; the efficient method of moments of \citeA{gallant1996moments} as applied in \citeA{chernov2000study} and \citeA{andersen2002empirical}; the implied-state methods initiated by \citeA{pan2002} and further analyzed by \citeA{santa2010crashes}; the Markov Chain Monte Carlo method in \citeA{eraker2004stock} and \citeA{eraker2003impact}; and the particle filtering method, see \citeA{johannes2009optimal} and \citeA{bardgett2019inferring}. Most of these estimation methods use as inputs option prices or a monotonic transformation thereof, such as implied volatilities. By contrast, we propose an estimation procedure based on the prices of spanning option portfolios that by the bijection between CCFs and conditional distributions, in principle, contain all probabilistic information about the stochastic process governing the dynamics of the underlying asset.

In general, estimation strategies based on the transform space of conditional characteristic functions are, of course, not new to the literature. For instance, \citeA{carrasco2000generalization} develop a generalized method of moments (GMM) estimator with a continuum of moment conditions based on the CCF; see also \citeA{singleton2001estimation}, \citeA{carrasco2007}. In applications to option prices, \citeA{BLL2015} and \citeA{boswijk2021jump} propose to imply the latent state vector from a panel of options and then estimate the model via GMM with a continuum of moments. \citeA{bates2006maximum} develops maximum likelihood estimation and filtering using CCFs. In particular, he proposes a recursive likelihood evaluation by updating the CCF of a latent variable conditional upon observed data. However, unlike our approach, these methods require numerical integration over the dimension of the state vector, thus suffering from a `curse of dimensionality'.

Our work is also related to \citeA{feunou2018risk}, who exploit the linear relation between the first four risk-neutral cumulants of the log-asset price and latent factors. They obtain these cumulants via a portfolio of options and employ the Kalman filter to estimate the latent factors. The main difference with our approach is that we exploit the CCF, and the corresponding state space representation we develop, instead of the first four moments. The CCF contains much richer information, leading to more efficient inference. Another difference is in dimension reduction: \citeA{feunou2018risk} use a two-step principal components analysis (PCA) to reduce the dimension of the risk-neutral cumulants observed at different maturities. Instead, we use a modified version of a so-called collapsed Kalman filtering approach, originally developed by \citeA{jungbacker2015likelihood}, which does not suffer from information losses relative to the full-dimensional setting.

The paper is organized as follows. Section (ref) provides the theoretical framework for aligning the option-implied and model-implied CCFs. In Section (ref), we develop the state space representation, and establish the main result about the orders of measurement errors, under a double asymptotic scheme. This allows us to next develop the filtering approach and corresponding estimation procedure. Section (ref) presents the Monte Carlo simulation results. We describe the data in Section (ref) and the empirical applications in Section (ref). Conclusions are in Section (ref). In supplementary material, four appendices provide details on ($i$) the proof of Proposition (ref), ($ii$) the computation of conditional moments, ($iii$) the inter- and extrapolation scheme for option prices and the measurement errors in the CCF replication, and ($iv$) additional simulation and empirical results.

Theoretical Framework

In this section, we provide the theoretical framework for our approach. We start with extracting information about the CCF from option prices allowing for general underlying dynamics. Next, we consider the CCF within the AJD class, which is exponentially affine in the model's state variables. Finally, we discuss how to align the two CCFs---option-implied model-free and AJD model-implied---in order to conduct inference about the model parameters and the latent state variables.

Option-implied CCF

Throughout the paper, we fix a filtered probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq 0}, \mathbb{P})$. On this probability space, we consider the dynamics of an arbitrage-free financial market. The no-arbitrage assumption guarantees the existence of a risk-neutral probability measure $\mathbb{Q}$, locally equivalent to $\mathbb{P}$. Since we are interested in exploiting information from options, we formulate the model dynamics under $\mathbb{Q}$.

Let us denote by $F_t$ the futures price at time $t$ for a stock or an index futures contract with some fixed maturity. The absence of arbitrage implies that the futures price process is a semimartingale. In this subsection, we assume the following general dynamics for $F_t$ under $\mathbb{Q}$:

where $v_t$ is an adapted, locally bounded, but otherwise unspecified stochastic volatility process; $W_t$ is a standard Brownian motion; $\mu$ is a counting random measure with compensator $\nu_t(\mathrm{d} x)\mathrm{d}t$ such that $\tilde{\mu}(\mathrm{d}t, \mathrm{d} x) := \mu(\mathrm{d}t, \mathrm{d} x) - \nu_t(\mathrm{d} x)\mathrm{d}t $ is the associated martingale measure and $\int (x^2 \wedge 1)\nu_t(\mathrm{d} x) < \infty$.

We further denote out-of-the-money (OTM) European-style option prices at time $t$ with time-to-maturity $\tau > 0$ and strike price $K>0$ by $O_t(\tau, K)$. Under the no-arbitrage assumption, the option prices equal the risk-neutral conditional expectations of the corresponding discounted payoff functions:

The OTM price $O_t(\tau, K)$ is a call option price if $K > F_{t}$ and a put option price if $K\leq F_{t}$. For simplicity, we assume a constant interest rate $r$.

Following \citeA{carr2001optimal}, any twice continuously differentiable European-style payoff function $g(F_{t+\tau})$, with first and second derivatives $g_{F}$ and $g_{FF}$, can be spanned via a position in risk-free bonds, futures (or stocks) and options with a continuum of strikes, as follows:

Here, $x\in\mathbb{R}^{+}$, the first and second terms on the right-hand side correspond to risk-free bonds and futures positions, and the third and fourth terms correspond to OTM options. Taking conditional expectations under the risk-neutral measure for $x=F_t$, we find that the price at time $t$ of a contingent claim with payoff function $g(F_{t+\tau})$ can be expressed as a weighted portfolio of a risk-free bond and OTM options:

This general spanning result lies behind the construction of one of most popular `fear' indices---the VIX index, when $g(F_{t+\tau}) = \log(F_{t+\tau}/F_t)$. Some other applications of the spanning formula (ref) include the calculation of the option-implied skewness and kurtosis bakshi2003stock and of the corridor implied volatility andersen2007construction.

Applying this result to the complex-valued payoff function $g(x) = e^{\mathrm{i }u \log(x/F_t)}$ yields that the discounted CCF of log returns can be spanned as

where $m = \log (K/F_t)$ is the log-moneyness of an option with strike price $K$.\footnote{With slight abuse of notation, we use the same symbol $O_t$ for the option value as a function of $(\tau,K)$ and as a function of $(\tau,m)$.}

It is important to emphasize that the spanning of the CCF in equation (ref) is exact and is furthermore completely model independent akin to the VIX construction. Therefore, the CCF of log returns over a particular horizon $\tau$ can be replicated in a model-free way given a single cross-sectional slice of liquid option prices with all strikes (and the same maturity $\tau$). A similar approach of CCF spanning is taken by \citeA{todorov2019nonparametric} to nonparametrically estimate spot volatilities from option prices (considering the limit as $\tau \downarrow 0$).

The expression for $\phi_t(u,\tau)$ in (ref) cannot be computed in reality as we do not observe option prices for a continuum of strikes. Nevertheless, as we detail in Section (ref), the expression in (ref) is easy to approximate using a limited number of observable option prices. When developing our estimation procedure, we take both the resulting approximation errors as well as the observation errors in option prices, and hence in the CCF approximation, into account. Henceforth, we denote by $\widehat{\phi}_t(u,\tau)$ the computationally feasible counterpart of the option-implied CCF; it is explicitly defined in (ref) below. In our simulation experiments and empirical applications, we further employ an interpolation-extrapolation scheme to improve the reliability of the approximation.

Affine jump-diffusion CCF

Whereas the CCF in (ref) is model independent, the CCF of log returns of the underlying asset is often considered under some parametric assumptions on the return dynamics. A model-implied CCF depends on the model parameters, which we generally do not know, and potentially on the dynamics of other latent processes, which affect the distribution of returns. Therefore, by suitably aligning the model-free and parameter-dependent CCFs, we may learn about the model parameters and the unobservable state dynamics.

We restrict our attention to the broad class of AJD models defined in \citeA{DPS2000}. The main attraction of the AJD class is that the Laplace transform has a semi-closed-form expression and is of the exponential-affine form. Suppose that $X_t$ is a Markov process representing an $d_X$-dimensional state vector in $D \subset \mathbb{R}^{d_X}$ with the first component being the log price of an asset. We assume that under the physical and risk-neutral probability measures, the state vector $X_t$ solves the following stochastic differential equation:

where $W_t$ is a standard Brownian motion in $\mathbb{R}^{d_W}$; $\mu{:}\ D \to \mathbb{R}^{d_X}$ and $\sigma{:}\ D \to \mathbb{R}^{d_X \times d_W}$ are the drift and diffusion functions; $N_{i,t}$ is a pure jump process with intensity $\{\lambda^{i}(X_t;\theta){:}\ t\geq 0\}$, $\lambda^{i}{:}\ D \to \mathbb{R}^{+}$; $\{J_{i,t}\}_{t\geq 0}$ constitutes a sequence of jump sizes with generic conditional distribution $\nu^{i}$ on $\mathbb{R}^{d_X}$ for $i=1,\dots,d_J$; and $\theta$ is a vector of unknown parameters that governs the model for $X_t$. We note that we allow for multiple jump types each arriving with their own intensity process as in the generalized AJD class in Appendix B of \citeA{DPS2000}. The specification (ref) can be extended further, e.g., to include a time-dependent structure and infinite activity jumps; see \citeA{DPS2000}, in particular Appendix B, and \citeA{duffie2003affine} for more details on the AJD class formulation.

Following \citeA{DPS2000}, the drift $\mu(x)$, diffusive variance $\sigma(x)\sigma(x)'$ and jump intensities $\lambda^{i}(x)$ are assumed to be affine on $D$:

where $x_{j}$ is the $j$-th element of a vector $x$ and $H_1^{(j)}$ for $j{=}1,\dots,d_X$ form a $d_X {\times} d_X {\times} d_X$ tensor $H_1$ by stacking matrices along a new dimension. The joint regularity conditions on $(D, \mu, \sigma, \lambda, \nu)$ that guarantee a unique solution to the SDE (ref) are discussed in \citeA{duffie1996yield} and \citeA{dai2000specification}. These joint conditions put constraints on the parameter vector $\theta$. Therefore, we consider a model from the AJD class indexed by $\theta$ in a parameter space $\Theta$ containing such admissible parameter values, on which there is a unique solution to (ref) that remains in $D$. For instance, in the case of the stochastic volatility component, the admissible parameter values in $\Theta$ ensure that the volatility process remains nonnegative, by satisfying Feller's condition; see also the discussion of the admissibility problem in \citeA[Chapter 5]{singleton2009empirical}.

\citeA{DPS2000} show that the affine dependence of the functions $\mu(x)$, $\sigma(x)\sigma(x)'$ and $\lambda(x)$ implies an exponential-affine form of the CCF of the state vector $X_t$. Specifically, the discounted joint CCF of $X_{t+\tau}$ conditional on $\mathcal{F}_t$ with $\tau > 0$ is given by

where $\mathbf{u} \in \mathbb{R}^{d_X}$ is an argument vector and $\alpha(\mathbf{u}, \tau; \theta)$ and $\beta(\mathbf{u}, \tau; \theta)$ are solutions to the following complex-valued system of ordinary differential equations (ODEs) in time:

with initial conditions $\beta(\mathbf{u}, 0) = \mathrm{i} \mathbf{u}$ and $\alpha(\mathbf{u}, 0) = 0$. Here, $\chi^{i}(c) = \int_{\mathbb{R}^n}\exp(c \cdot z) \mathrm{d} \nu^{i}(z),\ c\in \mathbb{C}^{d_X}$, are jump transforms, which determine the conditional jump-size distributions. The ODE for $\beta$ is known as a generalized Riccati equation, whereas the solution for the second ODE can be obtained by simply integrating the right-hand side expression over time.

The affine dependence of the characteristic exponent $\alpha(\mathbf{u},\tau ; \theta) + \beta(\mathbf{u}, \tau; \theta){\cdot} X_t$ on the current state $X_t$ is even the defining property of the AJD class under some regularity conditions (see \citeNP{duffie2003affine}). In other words, the AJD class can be defined as a class in which characteristic exponents of $X_{t+\tau}$ given $X_{t}$ are affine functions of $X_t$. In fact, this is a key property in our estimation procedure. While it is also possible to obtain the CCF for some non-affine models, the exponential-affine form allows us to use linear Kalman filtering techniques in the estimation procedure. This is the main motivation why we restrict our attention to the parametric models of the AJD class.\footnote{The considered AJD class could, in principle, be broadened further to the linear-quadratic jump-diffusion class by augmenting the state vector (see \citeNP{cheng2007linear}, for more details).}

Unlike the option-implied CCF (ref), the CCF in (ref) is fully parametric, that is, it requires parametric AJD model dynamics of the state vector $X_t$. Although the AJD class is more restrictive than the general dynamics of $F_t$ in (ref), it includes a myriad of popular option pricing models such as those in \citeA{heston1993}, \citeA{DPS2000}, \citeA{pan2002}, \citeA{bates2006maximum}, \citeA{broadie2007model}, \citeA{BLL2015}, and \citeA{AFT2017} among many others.

The state process $X_t$ often includes both observed and unobserved state variables that affect the dynamics of the log futures price $\log F_t$. In our empirical application, we consider the presence of both. Therefore, it is convenient to partition the state vector as $X_t' = (w_t', x_t')$, where $w_t$ represents the observable component and $x_t$ includes $d < d_X$ latent state variables. Then, the dynamics of $X_t$ given by equation (ref), can be rewritten as

where $\mu^{w}{:}\ D \to \mathbb{R}^{d_X-d},\ \mu^{x}{:}\ D \to \mathbb{R}^{d},\ \sigma^w{:}\ D \to \mathbb{R}^{(d_X-d) \times d_J},\ \sigma^x{:}\ D \to \mathbb{R}^{d \times d_J}$ and $J_{i,t}^w$ and $J_{i,t}^x$ are marginal jump sizes of $J_{i,t}$ associated with $w_t$ and $x_t$, respectively. In the simplest case, the observable component includes only the log futures prices, that is, $w_t = \log F_t$. In more general settings, the stochastic volatility is often a main latent driver of the log returns dynamics, as e.g., in \citeA{heston1993}.

Marrying the two CCFs

Given the two CCFs (ref) and (ref), we can now align them to conduct inference about the model parameters and the unobservable state variables. For that purpose, first note that the CCF in (ref) is joint for the state vector $X_t$. We assume, without loss of generality, that the first component of the state vector $X_t$ is the log futures price. Therefore, we can easily obtain its marginal CCF by plugging in an argument vector of the form $\mathbf{u_1} := (u, 0,\dots, 0)' \in \mathbb{R}^{d_X}$ with $u\in \mathbb{R}$. To obtain the marginal CCF of log returns, we further subtract the term $\mathrm{i}u \log F_t$ in the exponent. That is, the marginal CCF of log returns under the AJD specification is aligned to that in (ref) as follows:

where $\tilde\beta(\mathbf{u_1},\tau; \theta) := \beta(\mathbf{u_1},\tau; \theta) - \mathrm{i}\mathbf{u_1}$; i.e., the first component of $\tilde\beta(\mathbf{u_1},\tau; \theta)$ differs from that of $\beta(\mathbf{u_1},\tau; \theta)$, since we are interested in the CCF of log returns rather than that of log prices.

Note that the log of the (joint) CCF (also known as cumulant generating function) is linear in the state vector $X_t$. Therefore, under a correctly specified AJD model we obtain a simple linear relation between the log of the option-spanned CCF\footnote{Although the logarithm of a complex number is a multivalued function, here, the ambiguity is resolved given the fact that $\phi(0) = 1$ and the CCF is a continuous function. In fact, in practice we ensure that the logarithm of the CCF does not have `jumps' by taking the logarithm sequentially with respect to $u$, starting from the origin.} of log returns and the model's state vector:

Replacing the cumulant generating function on the left-hand side with its computationally feasible counterpart $\widehat \phi_t(u, \tau)$, which we will explicitly define in Section (ref), we obtain the following equation, which will play a central role in our estimation procedure:

Here, $\xi_t(u, \tau)$ is the measurement error, which is related to the observation, truncation and discretization errors in the CCF-spanning option portfolios. We elaborate in detail on the relation between the computable counterpart of the CCF and the source of the measurement errors in the next section.

Equation (ref) is the key relation in our analysis and a few remarks shall be made here regarding it. First, (ref) is essentially a functional linear model since this equation holds for any argument variable of the CCF, $u \in \mathbb{R}$. Furthermore, the functions $\alpha(\mathbf{u_1}, \tau; \theta)$ and $\tilde\beta(\mathbf{u_1}, \tau; \theta)$ are parameter-dependent and solutions to the system of Riccati ODEs (ref). Therefore, if the state vector $X_t$ is observable, then the model parameters can be estimated by solving a continuum version of a non-linear least-squares problem.

Second, in the case in which the state vector is (partially) unobservable, (ref) represents a linear latent factor model with a continuum of linear relations. The factors are given by the state components of the AJD model. Therefore, one could apply, e.g., a (functional) principal component analysis to learn about the unobserved factors. In this paper, we utilize a (suitably modified) Kalman filtering technique to conduct inference about the model parameters and the latent factors.

In other words, (ref) reveals that, using the present approach, AJD models become amenable to filtering and estimation using approaches from the rich literature on linear factor and state space models. This is reminiscent of the term structure literature, where in affine term structure models (see \citeNP{piazzesi2010affine}, for a review of this class of models) bond yields themselves are assumed to be linear functions of the state vector. For instance, \citeA{duffee1999estimating}, \citeA{jong2000time}, \citeA{driessen2005default} use the Kalman filter in their estimation of affine term structure models.

Furthermore, another advantage of this approach is that it does not require evaluating option prices given a certain parametric model. Therefore, our estimation procedure is computationally more appealing than many alternative approaches, which often involve the Carr-Madan FFT pricer (\citeNP{CM1999}) or the COS method (\citeNP{fang2008novel}) to price options. This also implies that the usage of the characteristic function is different: with the FFT or COS methods one needs a model-dependent CCF only to evaluate option prices, while here we use the CCF to directly learn about the latent factors and the model parameters.

Finally, given the partition of the state vector into observable and unobservable components, the linear relation between the option-implied and model-implied CCFs in (ref) can be rewritten as

where $\beta^{w}(\mathbf{u_1}, \tau; \theta) \in \mathbb{C}^{d_X-d}$ and $\beta^{x}(\mathbf{u_1}, \tau; \theta) \in \mathbb{C}^{d}$ are such that $\tilde\beta' = (\beta^{w\prime}, \beta^{x\prime})$ is the solution to the ODE system (ref). Representation (ref) serves as the basis for an observation (or measurement) equation in our estimation procedure.

Estimation Procedure

In this section, we develop our filtering approach and corresponding estimation procedure for the general class of AJD models under consideration. First, we provide the formal state space representation for the defined class of models. Then, we describe our estimation strategy, which uses the collapsed Kalman filter.

State space representation

As discussed in the previous section, we restrict our attention to the parametric models of the AJD class due to their exponential-affine form of the characteristic function. This form will allow us to exploit a linear Kalman filter in the estimation procedure. In the following, we summarize the assumptions we impose on the parametric model:

Assumption (ref) guarantees the existence of a unique solution to the SDE (ref) within the AJD class. As discussed in Section (ref), admissible values $\theta \in \Theta$ reflect the regularity conditions imposed on the model such that there is a unique solution to (ref), with, e.g., non-negative volatilities and jump intensities. Such admissibility conditions will need to be checked in a case-by-case model analysis. Assumption (ref)(i) also presumes the technical conditions required to represent the AJD process, defined via the affine dependence of its drift, diffusive variance and jump intensities on the state vector, through the exponential-affine characteristic function. For a detailed analysis of the AJD theory, we refer to \citeA{DPS2000} and \citeA{duffie2003affine}. Note that Assumption (ref) does not require the state process to be stationary. Stationarity of the latent state variables $x_t$ is reasonable but not essential for the results to follow; the observed state variables $w_t$ (often including the log-forward price) are typically non-stationary.

In our estimation procedure, we discretize the continuous-time model along two dimensions: with respect to time and with respect to the argument of the CCF. The former naturally follows from the discrete sampling times of financial data, which we denote by the integer indices $t=1,\dots,T$. The latter allows us to rely on the existing literature about filtering techniques. For that, let us denote the collection of discretely sampled arguments by a set $\mathcal{U} \subseteq \mathbb{R}$ with cardinality $q \in \mathbb{N}$. We further consider options with $k \in \mathbb{N}$ different maturities $\tau$ and $n \in \mathbb{N}$ different log-moneyness values $m$ on each day.

Since the input of our estimation procedure is a portfolio of option prices, we need to take into account the measurement errors in these option portfolios. For that purpose, we assume an observation error scheme on the option prices that constitute the portfolios. The measurement errors will be defined on the common probability space $(\Omega, \mathcal{F}, \mathbb{P})$, but in what follows, the filtration $\{\mathcal{F}_t\}_{t\ge 0}$ is generated by the state process $\{X_t\}_{t\ge 0}$ only. Note that the theoretical option prices $O_t(\tau, m)$ are $\mathcal{F}_t$-measurable, and hence the same applies to functionals of the option prices such as the (theoretical) Black-Scholes implied volatility (BSIV) and vega.

The additive error assumption is commonly imposed in the option pricing literature. For instance, \citeA{andersen2015parametric} and \citeA{todorov2019nonparametric} use additive error assumptions for option prices quoted in terms of BSIV and dollar amount, respectively. Additive observation errors are also often implicitly assumed when calibrating an option pricing model to market-observed prices, since the calibration is often performed using non-linear least squares as in, e.g., \citeA{broadie2007model}.

Assumption (ref)(i) excludes in particular dependence of the observation errors across strikes and is also often imposed in the literature (see, for instance, \citeNP{christoffersen2010volatility}, \citeNP{andersen2015parametric} and \citeNP{todorov2019nonparametric}). This assumption can be relaxed by introducing a spatial dependence as in \citeA{andersen2021spatial}. This would, however, result in more complex expressions for the covariance terms in the measurement errors that we derive below. Furthermore, \citeA{andersen2021spatial} find evidence of limited dependence in the observation errors for S&P 500 index options. They also show that this dependence declined sharply for short-dated options in recent years, due to improved liquidity. Since in our empirical application we consider S&P 500 index options with short tenors focusing on the past three years, the independence assumption will play a secondary role for the estimation procedure.

The conditional mean zero Assumption (ref)(ii) is crucial for our main result. Assumption (ref)(iii) asserts the standard deviation of the observation errors to be proportional to the product of the option's BSIV and vega. The motivation for this structure is as follows. Let $\widehat\kappa(m_j)$ and $\kappa(m_j)$ denote the error-distorted and true BSIV of an option, and assume that the relative volatility errors $\varkappa_j = (\widehat\kappa(m_j) - \kappa(m_j))/\kappa(m_j)$ are homoskedastic across the strikes, such that $\mathbb{E}[\varkappa_j^2|\mathcal{F}_t] = \sigma_\varkappa^2$. A Taylor-series expansion of the Black-Scholes pricing function $O^{BS}(\widehat\kappa(m_j), m_j)$ around $\kappa(m_j)$ then gives $\widehat{O}(m_j) = O^{BS}(\widehat\kappa(m_j), m_j) \approx O(m_j) + \nu(m_j) \kappa(m_j) \varkappa_j $, with $\nu(m_j) = \partial O^{BS}(\kappa(m_j), m_j) / \partial \kappa(m_j)$ the theoretical Black-Scholes vega. Homoskedastic errors in relative implied volatilities are also assumed by \citeA{christoffersen2012dynamic} and \citeA{du2019pricing} in their MLE based on the particle filter and the unscented Kalman filter, respectively.

Finally, to assess the error sizes of the CCF approximation specified below, we impose the following assumption on the existence of moments for the underlying asset and on the log-moneyness grid that allows nonequidistant sampling in the moneyness dimension:

Using $n > 1$ observable option prices with time-to-maturity $\tau >0$ and log-moneyness values $\{m_j\}_{j=1}^n$, we may approximate the CCF $\phi_t(u,\tau)$ given in (ref) by replacing the theoretical option prices by their observed counterparts, and the integral by a Riemann sum:

where we use the notation $u_t := (u^2 + \mathrm{i} u)/F_t$, and where $\widehat{O}_t(\tau, m_j)$ satisfies Assumption (ref).

The deviation of the option-spanned CCF from its theoretical counterpart, $\zeta_t^{\phi}(u,\tau) := \widehat{\phi}_t(u, \tau) - \phi_t(u, \tau)$, stems from observation, truncation and discretization errors, where truncation refers to the fact that the integration interval $[\underline{m},\overline{m}]$ does not cover the entire real line. The truncation and discretization errors also arise in VIX calculations and depend on the availability of option prices. They will be shown to be of smaller order than the observation errors, and can further be efficiently reduced by using an interpolation-extrapolation scheme (see, e.g., \citeNP{jiang2005model,jiang2007extracting}, \citeNP{chang2012option}, and Appendix (ref)). Appendix (ref) illustrates the impact of the three different types of measurement errors on the CCF approximation, and the effectiveness of the interpolation-extrapolation scheme.\footnote{The interpolation-extrapolation scheme may induce some cross-sectional dependence in the observation errors $\zeta_t(\tau,m)$. This is in deviation from Assumption (ref), which is only realistic when referring to the errors before application of the interpolation-extrapolation scheme. We will not consider this effect explicitly in Proposition (ref) that follows; it would lead to a more complicated expression for the covariance matrix of the measurement errors, but, importantly, would not affect the main result otherwise.}

From the preceding analysis, the functional measurement equation (ref) is then obtained using the following log-linearization:

where the log-linearized observation errors $\xi_t^{(1)}(u,\tau)$ are defined by $\zeta_t^{(1)}(u,\tau)/\phi_t(u, \tau)$, with

and where $r_t(u,\tau)$ is a remainder term that collects the log-linearized truncation and discretization errors as well as the higher-order terms from the required Taylor-series expansion. (The superscript $^{(1)}$ refers to the first, and prime, source of the measurement errors, the observation errors; see also the detailed decomposition in equation (ref).)

To formulate the main result, we turn the complex-valued functional measurement equation (ref) into a real vector measurement equation, as usual in state space model formulations. First, we stack the log CCF and the corresponding measurement errors along $q$ values $u_1,\ldots, u_q$ for the CCF argument $u \in \mathcal{U}$, for a fixed expiration period $\tau_i$:

In a similar way, we denote by $a_{t,i}$, $b^w_{t,i}$ and $b^x_{t,i}$ the stacked outputs\footnote{Here, we attribute these elements (and system matrices $\tilde{d}_t, W_t$ and $Z_t$ in equation (ref)) with an additional time index although the coefficient functions are assumed to be time-invariant in the exposition. This is because in practice we can have different expiration periods for different days.} of the functions $\alpha(u,\tau_i)$, $\beta^w(u,\tau_i)$ and $\beta^x(u,\tau_i)$, respectively. Next, to tackle the complex-valued measurement equation (ref), we stack the real and imaginary parts, as well as $k$ maturities:

where $p=2qk$. Stacking the real and imaginary parts of the measurements is a natural approach when the state vector is real-valued;\footnote{See \citeA{singleton2001estimation} and \citeA{chacko2003spectral}, who use this approach in a GMM estimation setting based on the empirical characteristic function.} a complex-valued state vector would have required a complex Kalman filter based on the so-called widely linear complex estimator, as in \citeA{dini2012class}. The stacked observation equation (ref) links all available information from option prices with several tenors at time $t$ to the state vectors $w_t$ and $x_t$ in a linear way.

To complete the state space model, we need to augment the measurement equation (ref) by a transition equation for the unobservable state vector $x_t$. This is a linear, discrete-time dynamic system, to be derived from the continuous-time stochastic differential equation. An Euler discretization of the state process (ref) would converge to the true transition dynamics as the discretization step $\Delta t \to 0$. However, the maximum likelihood (ML) estimator based on the Euler discretization is, in general, inconsistent for fixed non-zero $\Delta t$ lo1988maximum, because the discretization has conditional moments different from those of the true process piazzesi2010affine. Fortunately, the AJD assumption under $\mathbb{P}$ implies that the first and second conditional moments of $x_{t+1}$ given $\mathcal{F}_t$ are linear and available in semi-closed form (possibly requiring the solution of a system of ODEs):

where $Q_t{:}\ \mathbb{R}^d \to \mathbb{R}^{d\times d}$ is an affine function in $x_t$. The finiteness of the conditional moments is ensured by Assumption (ref)(i). Both conditional moments will in general be linear in both the observed state $w_t$ and the latent state $x_t$; but because the former does not need filtering, we absorb its effect in the time-varying intercept $c_t$, and similarly in the intercept of the affine function $Q_t$.\footnote{The transition matrix $T_t$ will not be time-varying in stationary AJD processes with equidistant observations, but we do not impose this time-constancy in the notation, also to avoid confusion with the sample size $T$.}

In Appendix (ref), we show how these transition coefficients can be computed for the AJD model. Using this approach, which will in principle be model-dependent and hence has to be applied case by case, we obtain a discrete-time transition equation with the same conditional mean and variance as the true continuous-time process (but possibly different higher-order moments). Quasi-maximum likelihood (QML) estimation based on conditionally normally distributed measurement and transition errors in the state space representation yields consistent estimation results fisher1996estimating. A similar approach has been adopted in the term structure literature (see, e.g., \citeNP{jong2000time}, \citeNP{duffee2002term}).

We summarize the development of the state space representation, and analyze properties of the errors, in the following proposition. The main result contains a remainder term in the measurement equation that collects the truncation and discretization errors in the construction of $\log \widehat{\phi} (u,\tau)$ and higher-order terms in the log-linearization. This term vanishes under an asymptotic scheme, where $\overline{m} = \max_{1\le j \le n} m_j \to \infty$, $\underline{m} = \min_{1\le j \le n} m_j \to -\infty$ and $\Delta m \to 0$. We also denote the corresponding smallest and largest strike prices by $\underline{K}$ and $\overline{K}$, and express the asymptotic orders with respect to the number of option prices $n$ with fixed maturity.

The proof is given in Appendix (ref). The orders indicate that the remainder term goes to zero faster than the observation term given some minimum non-zero requirements for $\underline{\alpha}$ and $\overline{\alpha}$. In the sequel, we assume that $(\underline{\alpha} \wedge \overline{\alpha}) > \frac{1}{4}$ and neglect the remainder term in the estimation and filtering procedures. The system matrices $Z_t, T_t, Q_t(x_t)$ and system vectors $d_t$ and $c_t$ are known up to a parameter vector $\theta$, assumed to lie in the interior of a compact parameter space $\Theta$ by Assumption (ref)($ii$). Similarly, the system matrix $H_t$ depends on the data and $\theta$ (via $u_t, \phi_t, \kappa_t$ and $\nu_t$), and an additional unknown parameter $\sigma_{\varkappa}^2$. Note that $(d_t,Z_t,H_t)$ are derived from the $\mathbb{Q}$-dynamics of (ref), whereas $(c_t, T_t, Q_t(\cdot))$ correspond to the $\mathbb{P}$-dynamics. Therefore, possible deviations between $\mathbb{P}$ and $\mathbb{Q}$, reflecting the presence of factor risk premia, will require an extension of the parameter vector; we discuss this possibility further in Section (ref) and Appendix (ref). Estimation of $\theta$ and filtering of the latent state vector via (versions of) the Kalman filter is considered in the next sub-section.

Modified and collapsed Kalman filter

Consider the state space representation (ref)--(ref), where from now on we will ignore the remainder term $r_{t,n}$, and hence assume that the set of strike prices $\{m_j\}_{j=1}^n$ on each day is rich enough to make this term negligible. Define the dataset $Y_t = \{y_1, \ldots, y_t\}$, and linear projections (denoted by $\widehat \mathbb{E}$) of the latent state vector conditional on the data: $\widehat x_{t|t} = \widehat \mathbb{E}[x_t|Y_t]$ and $\widehat{x}_{t|t-1} = \widehat \mathbb{E}[x_{t}| Y_{t-1}]$, with corresponding mean square error matrices $P_{t|t} = \mathbb{E} [(x_t-\widehat{x}_{t|t})(x_t-\widehat{x}_{t|t})']$ and $P_{t|t-1} = \mathbb{E} [(x_t-\widehat{x}_{t|t-1})(x_t-\widehat{x}_{t|t-1})']$. Then a modified version of the Kalman filter reads as follows:

for $t=1,\dots,T$. If the latent state process $x_t$ is stationary, the initial conditions $\widehat{x}_{1|0}$ and $P_{1|0}$ for the filter can be set to the unconditional mean and variance, respectively.

In traditional homoskedastic Gaussian state space models, where the distribution of the vector $(\varepsilon_t^{\prime}, \eta_{t+1}^{\prime})^{\prime}$, conditional on $x_t$, is Gaussian with a constant variance matrix, the filtered state $\widehat{x}_{t|t}$ is the conditional expectation of the true process $x_t$ given the observations up to time $t$. When the errors are non-Gaussian homoskedastic, the filtered state represents the linear projection (or minimum mean square error linear predictor) instead of the conditional expectation. This property can be used to prove that quasi-maximum likelihood (QML) estimation based on the Gaussian likelihood still yields consistent and asymptotically normal parameter estimates hamilton1994time. In general AJD models, on the other hand, the distribution of the errors will be non-Gaussian with a conditional variance $Q_t(x_t)$ that is an affine function of the true latent state vector $x_t$. Therefore, the Kalman filter recursions have been modified by using $Q_t(\widehat{x}_{t|t})$ instead of the unobserved $Q_t(x_{t})$. A similar modification is used in, e.g., \citeA{jong2000time}, \citeA{monfort2017staying} and \citeA{feunou2018risk}. Although consistency of QML based on this modification has not been proved, Monte Carlo simulation results in these articles suggest that the method works well in practice.

Given the large dimension of the observation vector $p=2qk$, the Kalman filter and its QML estimation will be computationally challenging if not infeasible. In fact, an important caveat with this approach is that one needs a non-singular innovation variance matrix $G_t$. Since our CCF approximation is based on common option price data for $q$ different arguments $u$ and fixed time-to-maturity $\tau$, this matrix is likely to be (near-)singular for large $q$. Furthermore, with large cross-sectional dimension, the computation of the inverse matrix for each time $t$ adds a significant computational burden to the estimation procedure. To overcome these issues, we consider the collapsed Kalman filter, originally developed by \citeA{jungbacker2015likelihood}, which we describe below. We modify their method to allow for a (near-)singular variance matrix $H_{t}$, using generalized inverses.

The idea of the collapsed Kalman filter is to transform the observation vector $y_t$ into an uncorrelated pair of vectors $y_t^*$ and $y_t^+$ such that $y_t^*$ depends on the state vector $x_t$ and has dimension $d \times 1$, whereas $y_t^+$ does not depend on $x_t$ and has dimension $(p-d) \times 1$. Such a transformation can be done using, for instance, the projection matrices $A_t^* = (Z_t'H_t^- Z_t)^{-1} Z_t'H_t^-$ and $A_t^+ = L_t H_t^-(I_p - Z_t A_t^*)$, where $L_t$ is chosen such that $A_t^+$ has full row rank and where $H^-$ is the generalized inverse of $H$ and $I_p$ is the identity matrix of size $p$. Since $A_t^* Z_t = I_p$ and $A_t^+ Z_t = 0$, the observation equation is then transformed into

with $d_t^* = A_t^* d_t$, $d_t^+ = A_t^+ d_t$, $\varepsilon_t^* = A_t^* \varepsilon_t$ and $\varepsilon_t^+ = A_t^+ \varepsilon_t$. Using $H^-HH^-=H^-$, we have

In the preceding display, it has been assumed that $\mathrm{rank}(Z_t'H_t^{-}Z_t) =d$; this is not very restrictive, given that the dimension $d$ of the state vector will typically be much smaller than the dimension $p$ of the observation vector. We also require that the matrix $A_t = [A_{t}^{*\prime},A_{t}^{+\prime} ]'$ is non-singular, such that the transformation $A_t y_t$ does not lead to a loss of information.

The representation (ref) shows that information about the state vector $x_t$ is contained in the observation equation for $y_t^*$; thus we may ignore the second equation with $y_t^+$ and focus on the collapsed state space model:

Let us emphasize that the collapsing transformation into a lower-dimensional state space form is also valid for the Moore-Penrose inverse covariance matrix $H_t^{-}$. Therefore, we can collapse a high-dimensional data vector into a lower-dimensional vector even when the covariance system matrix of disturbances is (near-)singular.

The logarithm of the Gaussian likelihood function of the data vector $Y_T=(y_1', \dots, y_T')'$ is given by

where $p_{\theta}(y_t|Y_{t-1})$ is the (misspecified) Gaussian distribution of $y_t$ conditional on $Y_{t-1}$ (and $w_1,\ldots,w_{t-1}$), which can be evaluated via the prediction error decomposition based on the original state space representation (ref)--(ref). Given the assumption of a full rank transformation matrix $|A_t|$, the collapsed transformation allows to decompose the log-likelihood function $l(Y_T; \theta)$ into three parts to ease computation:

where $Y_T^*$ and $Y_T^+$ are stacked vectors of $y_t^*$ and $y_t^+$ over $t=1,\dots,T$, respectively.

The first term in (ref) is the quasi-loglikelihood evaluated by the Kalman filter applied to the collapsed state space system (ref)--(ref):

where $\omega_t^{*}$ are the prediction errors and $G_t^*$ are their mean square error matrices from the Kalman filter.

Since $y_t^+$ does not depend on the state vector $\alpha_t$ and $|H_t^+|=1$ may be imposed without loss of generality, the second term in (ref) is given by

Fortunately, the last term in the expression above can be calculated without construction of the matrix $A_t^+$:

where $M_Z = I - Z_t(Z_t'H_t^-Z_t)^{-1}Z_t'H_t^- = I - Z_t A_t^*$,\ $J_t^+ = A_t^{+\prime}(A_t^{+}H_t A_t^{+})^{-1} A_t^{+} H_t$ and $e_t = M_Z (y_t - d_t)$, that is, these are the generalized least squares (GLS) residuals from the observation vector $y_t$ with the covariate matrix $Z_t$ and variance matrix $H_t$. For derivation details,\footnote{ The derivation in \citeA{jungbacker2015likelihood} is based on the invertible covariance matrix $H_t$, but the same result and the same derivation are valid when using the pseudo-inverse matrix $H_t^-$. } see \citeA{jungbacker2015likelihood}.

Finally, the third term in (ref), $|A_t|$, can be found from the relation

which follows from the fact that the covariance matrix $A_t H_t A_t'$ is block diagonal given the uncorrelated error terms $\varepsilon_t^*$ and $\varepsilon_t^+ $ and using again $|H_t^+|=1$.

Given the measurement error structure as implied by Proposition (ref), the single scale parameter $\sigma_\varkappa^2$ of the covariance matrix can be factored out as $H_t = \sigma_\varkappa^2 \cdot \widetilde{H}_t$. The matrix $\widetilde{H}_t$ has a block-diagonal structure; although its blocks depend on the state vector and parameters via the theoretical BSIV $\kappa_t(\tau, m)$ and vega $\nu_t (\tau,m)$, we estimate these quantities directly from the data, hence they are not updated as we optimize over $\theta$. Thus, we have from (ref) that

Therefore, the log-likelihood (ref) is proportional to

Note that the inversions and determinants of the matrices $G_t^*$ and $H_t^*$ can be computed efficiently since they have small dimensions $d\times d$. This eases maximization of the log-likelihood function (ref) substantially.

The quasi maximum-likelihood parameter estimates $\widehat{\theta}$ are obtained by maximizing (ref) over the model parameter space $\Theta$, where we implicitly assume that the parameter vector $\theta$ has been extended to include the additional parameter $\sigma_{\varkappa}^2$. Its asymptotic properties are analogous to QML estimation based on the (modified) Kalman filter, as discussed at the beginning of this sub-section. In cases in which the conditional covariance matrix $Q_t$ does not depend on the latent state vector $x_t$, and the latent state process $x_t$ is stationary, QML based on the Kalman filter will yield consistent and asymptotically normal estimators. When $Q_t$ is affine in $x_t$, then QML based on the modified Kalman filter appears to have comparable properties in Monte Carlo simulations, but no formal consistency proof is available.

Monte Carlo Study

In this section, we study the finite-sample performance of our estimation procedure. In particular, we consider two AJD specifications: a one-factor model and two versions of an option pricing model with two factors.

SVCDEJ

As a starting point, we illustrate the developed estimation procedure based on a modification of the widely used `double-jump' stochastic volatility model of \citeA{DPS2000}. The modification is due to using double-exponential (rather than Gaussian) jump sizes in returns as in \citeA{kou2002jump} and \citeA{andersen2015parametric}, and a stochastic (rather than constant) jump intensity that is a multiple of the stochastic variance as in \citeA{pan2002}. We label this specification as `SVCDEJ' for stochastic volatility model with co-jumps in volatility and double-exponential jumps in returns.

In particular, we assume the following data-generating process for the log forward price under both the $\mathbb{P}$ and $\mathbb{Q}$ probability measures:

where the two standard Brownian motions $W_1$ and $W_2$ are assumed to be correlated with coefficient $\rho\in[-1,1]$, and $N_t$ is a Poisson jump process with jump intensity proportional to the stochastic variance, $\lambda_t = \delta v_t$, $\delta>0$. We further assume that $J_t$ is a double-exponentially distributed jump size with generic probability density function

where $p^+$ and $p^-$ are probabilities of positive and negative jumps, respectively, and $\eta^+$ and $\eta^-$ are the corresponding conditional means of the jump sizes. We assume that all of these parameters are positive, $p^+ + p^- =1$ and $\eta^+<1$. Given the jump size distribution, the expected relative jump size in returns is

We allow the volatility to co-jump only with negative jumps in returns, with exponentially distributed jump sizes $J_t^v$ with mean $\mu_v>0$. Finally, we assume $\kappa$, $\bar{v}$ and $\sigma$ to be positive and impose Feller's condition $2 \kappa \bar{v} > \sigma^2$ and the covariance stationarity condition $\kappa > p^- \delta \mu_v$.

The model in (ref)--(ref) belongs to the AJD class and exhibits all important ingredients of option pricing models: stochastic volatility, jump components in returns and volatility, time-varying jump intensity and a self-excitation feature (because a negative jump in returns is associated with a positive jump in volatility, which increases the volatility and hence the jump intensity). Furthermore, this specification assumes a double-exponential jump size distribution in returns, which has recently been advocated in the literature (see, e.g., \citeNP{kou2002jump}, \citeNP{ait2015contagion}, \citeNP{andersen2015parametric} and \citeNP{bardgett2019inferring}).

Our developed estimation and filtering approach uses information from option prices, and is agnostic about equity risk premia. Indeed, the measurements are constructed as portfolios of options rather than the underlying asset. On the other hand, since the transition equation in the state space representation reflects the dynamics of the latent components (under $\mathbb{P}$), it is, in principle, possible to learn about the risk premia associated with the latent processes (for instance, the variance risk premium). However, additional simulation results, reported in Appendix (ref), suggest that the $\mathbb{Q}$-information in the option prices largely dominates the $\mathbb{P}$-information, making the identification of risk premium parameters weak. A similar difficulty of identifying the physical dynamics arises in the term structure literature (see, e.g., \citeNP{kim2012term}). Therefore, we assume no variance (or state related) risk premia, that is, the latent components have the same dynamics under both probability measures. Importantly, the results in Appendix (ref) suggest that estimation of the $\mathbb{Q}$-parameters is hardly affected by imposing this (possibly invalid) restriction.

The discounted marginal CCF of the log forward prices in the SVCDEJ model can be derived using the results in \citeA{DPS2000} and is given by

where $\alpha(u, \tau)$ and $\beta(u, \tau)$ are solutions to the complex-valued ODE system in time:

with initial conditions $\beta_1(u,0) = \mathrm{i} u,\ \beta_2(u,0) = 0$ and $\alpha(u,0) = 0$. Here the `jump transform' takes the form

The CCF of the log price in (ref) is used to price options. For the state space representation, we turn it into the CCF of log returns as described in Section (ref). Using the fact that the solution to the ODE system satisfies $\beta_1(u,\tau)=\mathrm{i}u$, the linear relation between the log of the option-implied CCF and the state vector is given by

where $\widehat \phi_t(u, \tau)$ is the option-implied CCF, $\tau > 0$ is the time-to-maturity of available options and $\xi_t(u, \tau)$ is the measurement error term due to observation and approximation errors in the option-implied CCF. We use this linear relation to construct the measurement equation as discussed in Section (ref).

Following Appendix (ref), the conditional mean and variance of the latent stochastic volatility process are given by

with $g_0 = \kappa \bar{v}$ and $g_1 = -\kappa + p^-\delta \mu_v$. Equations (ref)--(ref) are then used to define the state updating equation:

where $c_t =\frac{g_0}{g_1}\left(e^{g_1\Delta t} - 1 \right),\ T_t = e^{g_1\Delta t}$ and $\mbox{Var}(\eta_{t+1}| \mathcal{F}_t) =\mbox{Var}(v_{t+1}|\mathcal{F}_t) =: Q_t(v_t)$.

The model specification has nine parameters of interest and one additional parameter that characterizes the observation errors. We note that the parameter $\delta$ often enters as a multiple of $p^-$, which can possibly cause identification issues in the estimation procedure. Therefore, to avoid these identification issues, we fix the probability of negative jumps to be $p^- = 0.7$. This is consistent with findings in \citeA{ait2015contagion} and our empirical results for the unrestricted model provided in Appendix (ref), where we also assess the robustness of our empirical results to fixing $p^- = 0.7$.

In the simulation study, we use $T=500$ time points with $\Delta t =1/250$. The time-series of the log prices and true spot volatilities are simulated using an Euler scheme applied to the specification (ref)--(ref). The initial values are set to $F_0 = 100$ and $v_0=0.015$. The options data are generated using the COS method of \citeA{fang2008novel} based on the CCF, specified in (ref). The true model parameters are displayed in Table (ref).

In the simulations, we consider three tenors for options equal to 10, 30 and 60 days. For each tenor, we simulate a finite number of options with log-moneyness between $\underline{m} = -10 \cdot \sigma_{ATM, \tau} \sqrt{\tau}$ and $\overline{m} = 4 \cdot \sigma_{ATM, \tau} \sqrt{\tau}$, where $\sigma_{ATM,\tau}$ is the BSIV of the ATM option with time-to-maturity $\tau$. Furthermore, the strikes are generated equidistantly with $\Delta K = 0.01\cdot F_t$. Finally, we distort the options data by adding the observation errors to the option prices for each tenor $\tau$ and each log-moneyness level $m$ as specified in Assumption (ref), i.e.,

where $\epsilon$ is an i.i.d.\ standard normal random variable and $\sigma_{\varkappa}=0.02$. The distorted option prices, in terms of total implied variance, are then interpolated using a cubic spline and extrapolated linearly in log-moneyness, as described in Appendix (ref).

The covariance matrix of the errors in the measurement equation is calculated according to equation (ref). To calculate the pseudo-inverse of the $2q \times 2q$ covariance matrix $\widetilde{H}_{t,i}$ for each of the maturities $i=1,\dots,k$, we set the following level of the threshold for the singular values:

with $\bar{s} = 10^{-7}$ and where the maximum is taken over all singular values $s_j$ of $\widetilde{H}_{t,i}$. We also analyze the robustness of our results to the choice of $\bar{s}$.

With these specifications, we take the number of replications $N$ to be $N=300$, thus running the estimation procedure of Section (ref) 300 times. Table (ref) provides the Monte Carlo results for the SVCDEJ model, for six different ranges of the argument set $\mathcal{U}$. The results in general show a good finite-sample performance. We notice that for smaller ranges of the CCF argument, the estimates exhibit biases for some model parameters. This is expected since the smaller ranges provide coarser information on which we build the filtering and parameter estimation procedures. On the other hand, we also notice that the variance of some parameter estimates increases when using a very large range of arguments (in particular, $u=1,\dots,30$). This is likely due to an increased variance in the CCF approximation for large arguments $u$.

To explore the robustness to the choice of the truncation level in the pseudo-inversion of the covariance matrix, we also consider other values of $\bar{s}$. In particular, we run $N=300$ simulations for each level of $\bar{s}$ using the same parameter values as in Table (ref), and construct the root mean square percentage error (RMSPE) metrics, defined as the square root of $N^{-1} \sum_{i=1}^N \sum_{j=1}^{d_{\theta}} \left ( (\widehat{\theta}_{i,j} - \theta_{0,j})/\theta_{0,j} \right ) ^2 $, with $d_{\theta}$ the dimension of $\theta$. Figure (ref) plots the resulting RMSPEs for different levels of $\bar{s}$ and three different ranges of the argument $u$. As we can see, the levels $\bar{s}$ in between $10^{-7}$ and $10^{-6}$ yield the smallest RMSPE. In the following simulations and empirical applications, we therefore set $\bar{s} = 10^{-7}$.

We end this subsection by noting that we have also conducted simulation studies for some related alternative one-factor specifications. In particular, in Appendix (ref), we provide additional simulation results for the `SVCJ' model with a Gaussian jump size distribution, and the `SVCEJ' model with two separate counting processes for positive and negative jumps. The former shows a very good finite-sample performance, while the latter, a richer model specification, shows reasonable results, gradually reaching the limits of what can be identified using the present input data and design.

SVCDEJ with external factors

Now we extend the one-factor specification by adding an external factor. This modification can be seen as a two-factor specification, but we will assume that the second factor is observable. The motivation comes from the fact that in some situations we might have an understanding of possible drivers of the risks in the market. Therefore, we would like to embed exogenous variables into the model's risk factors and quantify their impact.

In particular, next to the stochastic volatility component we introduce the exogenous factor $h_t$, which affects the intensity of jumps and the diffusive component. The model reads as follows:

where $V_t = v_t + q^2 h_t$ is the total diffusive variance of the process and the jump intensity process $\lambda_t$ is also affected by $h_t$ with $\lambda_t = \delta v_t + \gamma h_t$, $q,\delta,\gamma>0$. We assume that $W_{3,t}$ and $W_{4,t}$ are independent standard Brownian motions, jointly independent of $(W_{1,t},W_{2,t})$. The process $h_t$ is exogenous to the SVCDEJ dynamics, meaning that the dynamics of $\log F_t$ and $v_t$ do not affect the dynamics of $h_t$. In turn, the exogenous factor $h_t$ affects the intensity of jumps and the diffusive component of the log return dynamics. This specification is similar to the two-factor model in \citeA{andersen2015parametric}, which includes short- and long-term stochastic volatility components. The difference is that here the exogenous process $h_t$ is observable, although its parameters are unknown.

In the Monte Carlo simulations, we consider two possible estimation approaches. In the first approach, we assume a correct specification of the dynamics of $h_t$ with known true parameters $\kappa_h,\ \bar{h}$ and $\sigma_h$. In practice, these parameters can be pre-estimated given the observed path of the exogenous process. In the second approach, we estimate the misspecified model in which the contribution of $h_t$ is constant throughout the maturity of an option. In other words, under this approach we ignore the dynamics of $h_t$ when pricing options, but let $h_t$ still affect the level of the jump intensity and of the total variance. The motivation is that when the exogenous process is persistent and smooth relative to $v_t$, its dynamics can be neglected when pricing options with short expiration periods. In a similar way, interest rates are often assumed to enter option prices in a deterministic way. Moreover, the true parametric specification for an exogenous variable is likely unknown in practice, but if its dynamics are persistent and smooth, we can find its effect on option prices via this approach. Therefore, in the Monte Carlo experiment, we simulate $h_t$ with a mean-reversion rate that is smaller than that in $v_t$, mimicking the specification we will explore in the empirical application.

For the rest, the Monte Carlo setting for the SVCDEJ model with an external factor is the same as for the SVCDEJ specification in the previous subsection. The parameters of the external factor are set to $\kappa_h = 1$, $\bar{h}=1$ and $\sigma_h=0.1$. The simulation results are provided in Table (ref). The parameters of the SVCDEJ model exhibit similar good performance under both estimation approaches. Importantly, the parameters related to the external factors, $\gamma$ and $q$, also show similar good performance in the correctly specified model as in the misspecified setting. We emphasize that this is achieved due to simulating a relatively smooth and persistent exogenous process $h_t$ and using short-dated options in the estimation procedure.

Data

This section describes the data and the data selection process, which we use in our empirical application. Since our estimation procedure utilizes option-implied CCFs, we also pay attention to the construction of these objects in this section. Further details are in Appendix (ref).

Data description

In this paper, we use options data on the S&P 500 stock market index obtained from the Chicago Board Options Exchange (CBOE). We focus on the period from May 1, 2017, to April 1, 2021, covering in particular the turbulent period in the stock market due to the outbreak of the Covid-19 pandemic. The CBOE provides end-of-day option quotes and a snapshot at 3:45 pm ET, 15 minutes prior to the market closure. We use the latter to calculate mid-quotes since it is considered to be a more accurate representation than the former in view of market liquidity. The data contain both the `standard' AM-settled SPX options and Weeklys and End-of-Months PM-settled SPXW products. The settlement value for the SPX options is based on the opening level of the S&P 500 index on the settlement day, whereas for the SPXW options it is based on the closing prices of the index.

Given that we need a reliable and wide coverage of option prices for each tenor, we use a fairly generous set of filters. In particular, we retain option observations that satisfy the following criteria: ($i$) bid price is strictly positive and ask-to-bid ratio is less than a factor 10; ($ii$) the maturity is larger than or equal to 2 calendar days, but less than or equal to 365 calendar days; ($iii$) it is not an early-closure day. The first criterion filters out illiquid observations and the second one limits our consideration in terms of options' maturity. The third criterion rules out shortened trading sessions, which in total constitute 10 days in our sample.

For each tenor, we determine the moneyness based on the forward index level, $F_t(\tau)$. For that, we use the put-call parity to calculate the forward price for close to at-the-money (ATM) options. Specifically, we use up to 5 option pairs with the smallest absolute difference between the call and put prices. The median of their forward-implied prices is taken as the forward index level for the corresponding tenor. The risk-free rates are obtained by interpolating the LIBOR rates to any particular tenor. Finally, given the calculated forward prices and moneyness levels, we retain only out-of-the-money (OTM) options for further exploitation. Descriptive statistics of the S&P 500 index options data sample are provided in Table (ref). We observe that the largest portion of the trading volume is due to trades of OTM contracts and options with time-to-maturity less than 60 calendar days. Figure (ref) plots the frequency of tenors up to 70 calendar days.

CCF-spanning option portfolios

The construction of the CCF-spanning option portfolios requires reliable option slices with wide coverage of strikes. Given that most of the trading volume is concentrated in option contracts with time-to-maturity of less than 60 days, our empirical application relies on the use of short-dated option slices with expiration period of no more than 2 months. In particular, on each trading day, we keep the six tenors closest to 8, 15, 22, 29, 36 and 61 days\footnote{The first five of these tenors are the most representative in the sample, see Figure (ref).} from below with the largest trading volume and number of quoted OTM option contracts. Specifically, starting with the option slices closest to the indicated tenors, we compare them with every next shorter maturity option slice, and prefer the next one if it has a larger trading volume and a larger number of quoted contracts for OTM options. Table (ref) provides the descriptive statistics for each of the six selected tenors over the considered time span. We notice the wide coverage of strikes, since the average minimum put and call prices are close to the tick size of \$0.05, especially for very short-dated options. We also mention that in the selected option sample, each option slice at each trading day contains at least 55 different quoted contracts. Therefore, no additional filters on the minimum number of contracts are imposed. In total, we have 978 trading days, with six different tenors at each one of them, resulting in a total number of 1,158,059 contracts in the sample.

The inputs of our estimation procedure are option portfolios representing CCFs rather than BSIVs that are commonly used in the literature. Therefore, we pay careful attention to the construction of the option-implied CCF. As discussed in Section (ref), we use a Riemann sum approximation to obtain a computationally feasible counterpart of the CCF spanning (ref). However, in order to reduce the truncation and discretization errors, we further employ an interpolation-extrapolation technique. In particular, we interpolate option prices using cubic splines with carefully selected knot sequences and extrapolate beyond the observable range of strike prices using a parametrization that satisfies the asymptotic results of \citeA{lee2004moment}. The details of the interpolation-extrapolation scheme are provided in Appendix (ref).

The calculation of the option-implied CCF then uses the Riemann sum approximation (ref) applied to the result of the interpolation-extrapolation scheme. The construction is conducted for each day and for each maturity separately. In particular, for equation (ref), we set $\Delta m = 0.0001$ with a sufficiently wide range of log-moneyness between $\underline{m} = -6$ and $\overline{m} = 2$.

To conclude this section we emphasize again that, contrary to what is common in many existing approaches, the option prices---or a monotonic transformation thereof---are not used as inputs in our developed estimation procedure. Instead, we use the option portfolios that replicate the CCF of log returns. Furthermore, unlike in many other papers, our option dataset is daily and utilizes the information from short-dated options with maturities between two days and two months.

Empirical Applications

Having thus constructed the dataset of option-implied CCFs for S&P 500 index options, we now illustrate our estimation procedure in two empirical applications, without and with an external factor.

SVCDEJ

We start with estimating the SVCDEJ model specified in Section (ref), (ref)--(ref), using the CCF-spanning option portfolios with six short-term tenors described in Section (ref). Table (ref) provides the parameter estimates. Informed by the Monte Carlo results, the estimates are based on the range of CCF arguments $u=1,\dots,20$, a singular value threshold $\bar{s}=10^{-7}$, and a fixed parameter $p^- = -0.7$. Standard errors are calculated using the familiar sandwich form covariance matrix.

The parameter estimates in Table (ref) are meaningful, intuitive and broadly consistent with the literature. For instance, \citeA{andersen2015parametric} find the mean jump sizes to be 1.71% and 5.33% for positive and negative jumps in their three-factor model specification. (They use, however, only the Wednesday options with a different sample period, from 1996 to 2010.)

We note that the leverage parameter $\rho$ is estimated close to its boundary value of $-1$, implying almost perfectly correlated diffusive components in returns and volatility. The empirical literature suggests that $\rho$ is negative and large in absolute value. The estimate of $\rho$ being nearly equal to its boundary value might be due to the use of short-dated options that typically exhibit steep implied volatility slopes. Indeed, \citeA{AFT2017} also find this correlation to be close to $-1$ in their dataset dominated by option contracts with maturities of less than 2 months.

We also note that the estimated measurement standard error $\sigma_{\varkappa}$ corresponds to a standard deviation of about 25% of the implied volatility. This is somewhat larger than what one might expect of measurement errors in option prices only, and might be interpreted to indicate e.g., missing state variables. In agreement with this, some of the extensions of the SVCDEJ model considered below and in Appendix (ref) show a slightly lower estimate of $\sigma_{\varkappa}$.

Figure (ref) plots the filtered volatility (i.e., the square root of the filtered state $\widehat{x}_{t+1|t}$) given the parameter estimates of the SVCDEJ model. As is clearly visible, the filtered volatility exhibits a relatively stable volatility regime prior to 2020 and jumps up in March 2020 at the outbreak of the Covid-19 pandemic.

SVCDEJ with external factors

Now we turn to model specifications with embedded external factors. In some situations, we might have specific information on possible drivers of the risks in the market, and would like to quantify their effect.

An example is the recent Covid-19 crisis. The Covid-19 pandemic has dramatically affected our lives. It has also had a tremendous impact on the world's economy and financial markets. The beginning of the pandemic, in particular, was associated with a spike in uncertainty. This uncertainty surrounded many aspects: the contagiousness and lethality of the virus, the time required to develop vaccines, the effectiveness of measures, the work-from-home policies, travel bans, and so on. In this application, we explore the impact of the Covid-19 pandemic on the stock market through the lens of option prices. In particular, we consider how the spread of the virus affected the likelihood of jump events and the volatility in the U.S. stock market.

Figure (ref), Panel (a), plots the daily cases of Covid-19 infections around the world obtained from the World Health Organization (WHO). The reported number of daily cases, however, does not represent well the contagiousness of the virus. Therefore, Panel (b) of Figure (ref) displays the so-called reproduction number $R_t$, according to two measures: the first one is taken from the website `Our World in Data' and the second one is calculated as the ratio $R_t = I_t/I_{t-7}$, where $I_t$ is the number of infected people in day $t$ and $7$ is the reported serial interval for Covid-19. The former is based on the parametric methodology of \citeA{arroyo2021tracking} and is smoothed over time.\footnote{In fact, \citeA{arroyo2021tracking} use a Kalman smoother.} The latter is non-smoothed and based on the assumption that the serial interval is $7$ days, which is consistent with the recent epidemiology literature (see, e.g., \citeNP{maier2020effective}, \citeNP{prem2020effect}, \citeNP{flaxman2020estimating}, \citeNP{arroyo2021tracking}). We will use the latter non-parametric and non-smoothed measure as the reproduction number in our application.

To quantify the effect of Covid-19 propagation on the financial market, we embed the reproduction number as an external factor into the (time-varying) levels of the stochastic volatility and jump intensity processes, as described in Section (ref), equations (ref)--(ref), with $h_t$ replaced by $R_t$. Given that the reproduction number constitutes a relatively persistent process, we will treat it as a deterministic process when pricing options; in other words, we follow the second estimation approach described in Section (ref). In a similar way, the risk-free rate and dividend yields are often assumed to be deterministic in the option pricing literature. This allows us to be agnostic about the parametric dynamics of the reproduction number. Furthermore, given the short-dated options under consideration, the errors due to this deterministic treatment are likely to be negligible.\footnote{Similarly, \citeA{AFT2017} and \citeA{boswijk2021jump} consider an approximation of the return process with `freezed' spot volatility when estimating their option pricing models with short-dated options.}

Table (ref) provides the parameter estimates for the SVCDEJ model with Covid-19 reproduction numbers as an exogenous factor. With $q$ estimated at $0.0003$, the results indicate that the reproduction number dynamics have no substantial effect on the total diffusive volatility. A one unit increase in reproduction numbers, however, leads to an increase in the intensity of jumps by $\gamma$ which is estimated at $2.64$. In other words, the reproduction number contributes substantially to the likelihood of jumps. Figure (ref) illustrates the dynamics of the jump intensity without and with the added effect of reproduction numbers.

It is also possible to investigate the contribution of other external factors to the diffusive volatility and jump intensity. As an example, we provide in Table (ref) estimation results for the SVCDEJ model with the Economic Policy Uncertainty (EPU) index embedded as an external factor. The EPU index, developed by \citeA{baker2016measuring}, reflects policy-related economic uncertainty based on newspaper coverage frequency. The estimation results indicate that, unlike the reproduction number, the EPU index has no effect on the jump intensity process, but contributes significantly to the total diffusive volatility of the model, with $q$ estimated at $0.0369$; see also Figure (ref). Thus, whereas Covid-19 reproduction numbers contribute substantially to the jump intensity dynamics, the policy uncertainty index EPU contributes significantly to the total diffusive volatility.

Conclusion

We have proposed a novel state filtering and parameter estimation procedure for option pricing models that belong to the affine jump-diffusion class. Our procedure utilizes the log of the option-implied and model-free conditional characteristic function and the model-implied conditional log-characteristic function, which is functionally affine in the model's state vector. We have developed a linear state space representation for the considered class of option pricing models, which allows us to exploit suitably modified collapsed Kalman filtering techniques. Our estimation procedure is fast and easy to implement, circumventing the typical computational burden when working with option pricing models. We have demonstrated the applicability of our procedure in two empirical illustrations that analyze S&P 500 index options and the impact of exogenous variables capturing Covid-19 reproduction and economic policy uncertainty data.

Although we have focused on Gaussian QML estimation based on Kalman filtering techniques, which delivers good results in our Monte Carlo simulations, the same state space formulation can also be analyzed by more refined methods such as those based on particle filters; see, e.g., \citeA{johannes2009optimal}, \citeA{christoffersen2014particle} and \citeA{bardgett2019inferring}. Such methods could exploit the non-Gaussianity and heteroskedasticity in the data to obtain more efficient estimates, at the cost of some increased computational complexity. We note that such extensions would still not require option price evaluation by the FFT or COS methods, and thus retain an important advantage of our approach.

Our proposed estimation procedure in principle allows for identification and estimation of factor risk premium parameters, by combining the risk-neutral parameters entering the measurement equation with the objective parameters entering the transition equation. Monte Carlo simulation results suggest, however, that option price data are not very informative about such risk premia, which is why we have concentrated on the case where the objective and risk-neutral measures coincide. Fortunately, the simulation results also suggest that inference on the risk-neutral parameters is quite robust with respect to deviations from this assumption. For more focused inference on (volatility) risk premium parameters, it may be possible to combine the information in daily option prices as considered in this paper with realized measures based on high-frequency returns on the underlying. We intend to explore this in future research.