Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
110,084 characters · 20 sections · 0 citation commands
iCOS: Option-Implied COS Method
{ Keywords: Derivatives; Option pricing; Risk-neutral distribution; Option Greeks; COS method.}\\ { JEL Classification: C13; C58; G13.}
Option prices play a crucial role in financial markets, providing investors and researchers with valuable insights into market expectations, risks, and dynamics associated with underlying instruments. This information is essential for effective risk management, portfolio optimization, and academic research aimed at understanding market dynamics. However, accurate extraction of this valuable information is a challenging task due to the complex nature of option prices and the presence of various sources of uncertainty.
Traditional option pricing models, commonly used to extract embedded information from options, rely on parametric assumptions about market dynamics and the distribution of asset returns. While these models have been widely used and have provided valuable insights, they often fail to capture the complexity and nuances of real-world market behavior and are often subject to model misspecification. To address these limitations, non-parametric methods have gained increasing attention in option pricing research. These methods aim to extract information about underlying dynamics directly from market data without making explicit assumptions about the underlying asset price dynamics.
In this paper, we propose a unified non-parametric estimation procedure for the risk-neutral density (RND), option prices, and option sensitivities, which are the key objects of interest in the option pricing literature. Our approach leverages the Fourier-based cosine technique, the COS method, proposed by \citeA{fang2008novel}, in a model-free way by implying information from observed option contracts. Therefore, we refer to this method as the (option-)implied COS method, or, in short, iCOS. The proposed estimation method is fully non-parametric and does not require any optimization routines, offering a flexible and computationally appealing alternative to traditional techniques.
Fourier-based methods are widely used numerical techniques for option evaluation that exploit the relation between the probability density function (PDF) and the characteristic function (CF). Examples of these methods include those proposed by \citeA{CM1999}, \citeA{lewis2001simple}, and \citeA{fang2008novel}. The popularity of these methods is due to the CFs -- the Fourier transforms of the PDFs -- that can be obtained in a (semi-)closed form for a large class of parametric models, such as the affine jump-diffusion class as defined in \citeA{DPS2000}. As such, the main requirement for these methods is a fully parametric specification of the CF, often induced from a parametrized asset price dynamics.
In contrast, the method proposed in this paper only requires the availability of observed option prices, which is a natural setup when working with option data. This allows us to bypass the unnecessarily restrictive requirement of a parametric CF and then utilize the flexibility of the COS machinery. In particular, we exploit the option spanning result of \citeA{carr2001optimal} to estimate the cosine series coefficients -- the projection coefficients of the orthonormal Fourier-cosine basis. These implied cosine coefficients are the building blocks for constructing non-parametric estimators of the RND, option prices with strikes that are not observed in the market, and option sensitivities.
In particular, by utilizing the spanning result of \citeA{carr2001optimal} and the COS machinery, we provide portfolio representations for the RND, option prices and option deltas. These portfolios consist of option contracts with a continuum of strike prices and offer model-free representations of the objects of interest. Since, in practice, we observe only a finite number of option contracts subject to observation errors, we develop non-parametric estimators based on these spanning results. We derive the asymptotic properties of these non-parametric estimators in an asymptotic setting in which the mesh of the strike grid shrinks to zero, while maturity and largest and smallest strike prices remain fixed. The resulting limiting distributions allow constructing confidence intervals for our estimates.
The existing semi- and non-parametric methods for the estimation or interpolation of option prices include parametric curve fitting (\citeNP{shimko1993bounds}, \citeNP{gatheral2014arbitrage}), local polynomial estimators (\citeNP{ait2003nonparametric}), (penalized) cubic splines (\citeNP{bliss2002testing}, \citeNP{malz2014simple}), and kernel-based methods (\citeNP{ait1998nonparametric}, \citeNP{grith2012nonparametric}) among many others. These methods are also often used to obtain the RND via the famous \citeA{breeden1978prices} formula. In comparison, our iCOS approach does not rely on (local) parametric representations and can be considered as a `global' non-parametric method. In other words, the developed estimators utilize all available option prices via the portfolio spanning results, while kernel-based approaches or local polynomial regression use only limited information. Furthermore, the proposed procedure provides a unified framework for estimating the RND, option prices and option sensitivities, without the need to take derivatives of estimated option functions.
There are also various methods that estimate the RND directly without invoking the \citeA{breeden1978prices} formula. Examples of such methods include the mixture of distributions (\citeNP{melick1997recovering}, \citeNP{gemmill2000useful}), the Gram–Charlier approximation (\citeNP{jarrow1982approximate}, \citeNP{rompolis2007retrieving}), and the Hermite expansion (\citeNP{xiu2014hermite}, \citeNP{lu2021sieve}), among many others. For a comprehensive review of various methods, see \citeA{figlewski2018risk}. In this paper, we utilize the Fourier-cosine expansion of the RND and imply the cosine expansion coefficients directly from the observed option prices. This offers a flexible, computationally appealing, and model-free alternative to traditional techniques. A similar direction is taken in \citeA{cui2021model}, \citeA{cui2022new}, and \citeA{bossu2022static}, who also imply the expansion coefficients from the observed option prices. In fact, \citeA{cui2021model} use the Fourier cosine method in a model-free way to extract the RND. However, our paper goes beyond these studies by proposing a unified framework to estimate the RND together with the option prices and option sensitivities. Furthermore, we explicitly control the truncation of the RND to a finite interval, allow for observation errors in options, and derive asymptotic results for the proposed estimators.
Our paper also contributes to the literature on model-free Greeks estimation (\citeNP{bates2005hedging}, \citeNP{alexander2007model}). Similar to these studies, our approach enables non-parametric estimation of option sensitivities, such as option deltas. However, unlike the methods proposed by \citeA{bates2005hedging} and \citeA{alexander2007model}, our approach does not require knowledge of the derivative of option prices with respect to strike and, hence, does not involve fitting the implied volatility curve. Instead, we estimate the option-implied deltas using an approach similar to the estimation of the RND and option prices, utilizing the portfolio spanning result based on the Fourier expansion. This eliminates the need for optimization and calibration to the market, thereby reducing the impact of model misspecification and calibration errors.
Our paper is also closely related to several studies that propose non-parametric approaches to estimate various risk measures from short-dated options, such as the Levy density in \citeA{qin2019nonparametric}, spot volatility in \citeA{todorov2019nonparametric}, and jump variation in \citeA{todorov2022nonparametric}. Like these papers, our approach is fully non-parametric and accounts for errors in the observed option prices. Our analysis differs in terms of the specific information extracted from the option contracts, and our methodology is not restricted to short-dated options, making it more broadly applicable.
We conduct extensive simulation experiments to assess the finite-sample properties of the proposed non-parametric estimators. We consider the classical and well-understood \citeA{black1973pricing} model and the more realistic `double-jump' stochastic volatility model of \citeA{DPS2000} as data generating processes. We find good finite-sample performance of all three non-parametric estimators for different maturities and show superiority of our procedure to kernel-based smoothing methods.
Finally, in our empirical application, we demonstrate the effectiveness of the iCOS method across various market settings. For that, we analyze the performance of the method in the highly liquid and well-studied market of S&P 500 index (SPX) options. We also consider Amazon (AMZN) options during an Earnings Announcement Day (EAD) in a unique, high-volatility conditions with short maturities and a distinct bimodal density pattern due to the EAD effect. Additionally, we apply our methodology to dissect and quantify errors in the VIX index, one of most popular measure of market volatility. We find that observation errors in the VIX, that arise from the imperfect observation of option prices, are centered around zero and have a small magnitude, reaching up to 0.04 percentage points of the index value. In contrast, discretization errors lead to positive biases, reaching up to 0.7 percentage points during periods of high volatility.
The rest of the paper is organized as follows. In Section (ref), we describe the option-implied COS method. The non-parametric estimators for the RND, option prices, and option deltas along with their asymptotic properties are discussed in Section (ref). Section (ref) provides the Monte Carlo simulation results. The empirical applications are in Section (ref). Section (ref) concludes the paper. The proofs of the propositions are collected in Appendix (ref) and some additional results are in Appendix (ref).
In this section, we start with the discussion of the COS method, proposed by \citeA{fang2008novel}. Then, we show how the option-implied information can be incorporated into this machinery. The option-implied COS allows us to get the expansion expressions for the RND, option prices and option sensitivities that serve as the building blocks for the non-parametric estimators introduced in the next section.
The COS method introduced by \citeA{fang2008novel} is based on the idea that the conditional density function $f(y)$ on an interval $[a,b] \subset \mathbb{R}$ can be represented via its Fourier cosine series expansion:
where $\sideset{}{'}\sum$ indicates the sum with the first term weighted by one-half, $u_m := \frac{m \pi}{b-a}$, and the cosine coefficients
\citeA{fang2008novel} showed that the cosine coefficients (ref) can be calculated via the (`truncated') characteristic function (CF). In fact, let us denote the CF of the density function restricted to the interval $[a,b]$ by
Then, premultiplying $\phi^{[a,b]}(u_m) $ by $e^{-\mathrm{i} u_m a}$ and taking the real part we find that
Working with the regular CF of the density $\phi(u_m)$ turns out to be more convenient than with its `truncated' version $\phi^{[a,b]}(u_m)$ since the CF is often available in a semi-closed form for many parametric option pricing models. Fortunately, the truncated version is well approximated by its infinite counterpart $\phi(u_m)$ for sufficiently wide interval $[a,b]$, and so are the cosine coefficients:
Therefore, the COS method allows for efficient pricing of options under any parametric model when a closed-form CF is available. In particular, consider a European-style option with a general payoff $v(y,T)$ as a function of the state variable $y$ at a maturity time $T$. Define the cosine series coefficients of the payoff function $v(y,T)$ as
Then, the COS formula to evaluate the price of this contract at time $t=0$ is derived by plugging the Fourier cosine series expansion of the conditional density into the risk-neutral valuation. That is, under the risk-neutral measure $\mathbb{Q}$ with the deterministic interest rate $r$, we have:
where by $\stackrel{\scriptscriptstyle\mathrm{(i)}}{\approx}$ we denote the subsequent numerical approximation. Note that the product of the two functions $v(y,T)$ and $f(y)$ is represented by the product of their Fourier-cosine series coefficients $A_m$ and $H_m$. The coefficients $H_m$ can be calculated analytically for many types of options. In Appendix (ref), we provide the analytic formulas for call options.
The COS method allows for fast option evaluation using the Fourier cosine expansion. There are three numerical approximations involved: (1) truncation of the integration range in the risk-neutral expectation, (2) usage of the CF $\phi(u_m)$ (and hence, $\widetilde{A}_m$) instead of the truncated counterpart $\phi^{[a,b]}(u_m)$, and (3) the cosine series truncation. Moreover, what is more important, this method requires a parametric model assumption, which is unknown a priori.
In contrast to the traditional COS method, the approach proposed in this paper does not rely on a parametric specification of the CF. Instead, we use a finite number of plain vanilla option prices observed in the market. Given these observable prices, we can extract the risk-neutral density, price European options with strike prices that are not observed in the market (i.e.\ perform interpolation), and compute the option sensitivities. Furthermore, as we show in the next subsection, our approach does not require a proper choice of the interval $[a,b]$, i.e.,\ it entirely eliminates the first two numerical approximation errors.
Let us denote by $S_t$ the underlying price at time $t$ for a stock or an index under consideration and by $F_t$ the futures price at time $t$ for this underlying asset with some fixed maturity. Let us further denote by $C_0(K)$ and $P_0(K)$ the call and put option prices at time $t=0$ maturing at date $T>0$ with a strike price $K$. Assuming the existence of an arbitrage-free financial market and denoting with $\mathbb{Q}$ the risk-neutral measure, the prices of out-of-the-money (OTM) options at time $t=0$ are given as the discounted risk-neutral conditional expectations of the corresponding payoff functions:
where $r$ is a deterministic interest rate. Throughout the paper we consider option contracts with a fixed maturity, thus we eliminate the dependence on $T$ in our notations.
The payoff spanning result of \citeA{carr2001optimal} allows expressing the value of any European-style contingent contract with a general payoff function $v(S_T)$ as a weighted portfolio of a risk-free bond and plain vanilla OTM options:
where $f_S(\cdot)$ is the risk-neutral density of the future underlying price, and $F:=F_0$ is the futures price at time $t=0$. This result has been extensively used in the literature to, e.g.\ construct the VIX index (\citeNP{CBOE}), extract the risk-neutral expectations (\citeNP{bakshi2003stock}), and imply the characteristic function in a model-free way (\citeNP{todorov2019nonparametric}, \citeNP{blv}).
For our purposes, it turns out to be convenient to modify this spanning formula to find values of contingent claims restricted to a finite interval $[\alpha, \beta] \subset \mathbb{R}$. For that, let us denote the risk-neutral valuation on a finite interval $[\alpha, \beta]$ as \[ v_0^{[\alpha, \beta]} = e^{-rT}\mathbb{E}^{\mathbb{Q}}[v(S_T) \mathbf{1}_{\{\alpha \leq S_T \leq \beta\}} ] = e^{-rT} \int_\alpha^\beta v(S_T) f_S(S_T) \mathrm{d} S_T. \] Then we have the following spanning result for $v_0^{[\alpha, \beta]}$.
The proof follows from the application of the general payoff spanning formula (ref) to a function $v(S_T) \mathbf{1}_{\{\alpha \leq S_T \leq \beta\}}$ and using properties of the Dirac Delta function. Like the spanning result of \citeA{carr2001optimal}, Proposition (ref) allows to replicate $v_0^{[\alpha, \beta]}$ by constructing a portfolio of risk-free asset and plain vanilla OTM options with strike prices from the interval $ [\alpha, \beta]$. Due to the restriction to the interval, the weights at the boundary option contracts, $C_0(\beta)$ and $P_0(\alpha)$, and risk-free asset position are adjusted. Note that with $\alpha \to 0$ and $\beta \to \infty$, the value $v_0^{[\alpha, \beta]}$ converges to the unrestricted contract value $ v_0$.
A natural choice for the truncated interval is the range of observable strike prices, i.e., we can set $\alpha$ to be the smallest observable strike price $\underline{K}$ and $\beta$ to be the largest observable strike $\overline{K}$. Alternatively, we can restrict the estimation to the interval with the most liquid options, i.e., $(\alpha, \beta) \subset (\underline{K}, \overline{K})$, which can be practically more appealing. This choice prevents the truncation errors, which are inevitable in the standard Carr-Madan spanning result (ref).
Now we can find the replicating portfolio for the (discounted) cosine coefficients $A_m$. For that, we consider the transformed\footnote{The motivation for this transformation comes from the COS method, where $x$ is the option's strike price. } variable $y= \log\frac{S_T}{x}$ (and, thus, $a= \log\frac{\alpha}{x}$ and $b= \log\frac{\beta}{x}$) with some $x>0$ and notice that
Applying Proposition (ref) to the function $v(S_T) = \cos\left(u_m \log \frac{S_T}{\alpha} \right)$ and denoting the second-order derivative of this function with respect to $S_T$ as
we get the (discounted) option-implied cosine coefficients as a portfolio of options:
Here, the term $b_m$ adjusts $D_m$ to account for the restriction to a finite interval, and with $\alpha \to 0$ and $\beta \to \infty$, the adjustment $b_m \to 0$. Therefore, we will also refer to $D_m$ as the cosine coefficient. It is important to emphasize that the cosine coefficient $A_m$ in equation (ref) is completely model-free. This is in contrast to the standard COS method, where one needs to specify a parametric assumption on the dynamics of the underlying asset to get the cosine expansion coefficients via parametric CF. Furthermore, the option-implied coefficients $A_m$ are exact, thus, approximation (2) in equation (ref) is avoided.
After having implied the cosine coefficients, we can obtain the risk-neutral density (RND) using equation (ref). For instance, setting $x=1$ gives us the RND of the log future price $\log S_T$:
where $ \nu_f := \tfrac{2e^{rT}}{\log\left(\beta/\alpha \right) }$. Equation (ref) is an important representation of a portfolio of option prices with strike prices within the interval $[\alpha, \beta]$. Like the cosine coefficients $A_m$, the spanning result for the RND in (ref) is exact and model-free. This serves as a basis for our non-parametric estimator of the RND, which we discuss in Section (ref).
It is worth noting that equation (ref) provides the values of the RND for any $y \in [a,b]$, although the density itself may have support on $\mathbb{R}$. This restriction to the finite interval is coherent with the availability of option data: if there are no options traded with strike prices $K < \alpha$, then it is difficult to infer information about the density accurately for $y < \log \alpha$ without making further (often parametric) assumptions.
After having extracted the option-implied cosine coefficients $A_m$, the same COS machinery can be used to price options but in a model-free way. While there is generally no need to price options that are already observed in the market, the developed iCOS method can be used to further evaluate option prices with the strikes that are not listed in the market. In other words, this approach allows interpolating plain vanilla options within the interval $[\alpha, \beta]$ in a completely model-free way. Additionally, as we discuss in Section (ref), option evaluation can be helpful in estimating unobserved quantities and constructing a feasible limiting distribution of the estimated RND and option sensitivities.
For the accurate pricing/interpolation, it is important to take into account the truncation levels, i.e., we shall separate the information available in the interval $[\alpha, \beta]$ from the information outside of this range. For that, we can decompose the risk-neutral valuation of a contract with the general payoff function $v(S_T)$ as follows:
That is, we can express the price of the contract as a sum of values over the three non-overlapping intervals. The motivation for this is to approximate the infinite counterpart by a value on $[\alpha, \beta]$, while possibly taking into account the values outside this interval. In fact, the COS method of \citeA{fang2008novel} assumes that the value of a contract on $[\alpha, \beta]$, $v_0^{[\alpha, \beta]}$, represents the contract value $v_0$ well, i.e., it assumes that the values $v_0^{(0,\alpha)}$ and $v_0^{(\beta,\infty)}$ are negligible.
It turns out that for the plain vanilla options, the values outside this finite interval $[\alpha, \beta]$ can be well controlled, completely eliminating the integration range truncation errors, represented by approximation (1) in equation (ref). In particular, for a call option with a strike price $x \in [\alpha, \beta]$, the value on the interval $(\beta,\infty)$ is given by
where $C_0(\beta)$ is the call price with the strike $\beta$ and $C_K'(\beta)$ is its derivative with respect to the strike price evaluated at $\beta$. $C_0^{(\beta,\infty)} (x)$ is the price of a so-called gap call option with a strike price $x$ and a trigger price $\beta$. A similar relation can be found with the gap put option for the value of a put option evaluated on the interval $(0,\alpha)$. See Appendix (ref) for the details about the put options.
Therefore, we have the following relations for the call and put option prices with strike prices $x$ such that $\alpha \leq x \leq \beta$:
These decompositions allow us to take into account the truncation of the integral in the risk-neutral valuation. In particular, by setting $\alpha$ and $\beta$ to the smallest and largest observable strike prices $\underline{K}$ and $\overline{K}$ respectively, we can account for the price of the gap call option using information from the observable range of strike prices. This means that our method is not affected by the choice of $[\alpha, \beta]$, and, hence, $[a,b]$, unlike the COS method.
Therefore, to interpolate call options, we compute $C_0^{[\alpha, \beta]} (x)$ using the COS machinery with the option-implied information from the corresponding interval and add the price of the gap call option. The value of the call contract truncated to the interval $[\alpha, \beta]$ is obtained using the COS formula as
where $H_m(x)$ are the cosine series coefficients specific to the call payoff function with the strike price $x$. The closed-form expression for $H_m(x)$ is provided in Appendix (ref).
Therefore, the price of a call option with strike price $x \in [\alpha, \beta]$ can be represented as
where we additionally denote $\theta_c := C_K'(\beta) $ and $ \theta_p := P_K'(\alpha) $. Equation (ref) represents the price of a call option as a portfolio of a continuum of OTM contracts with strike prices from the interval $[\alpha, \beta]$, with additional hedging terms due to the truncation on the finite interval. It is important to note that the decomposition (ref) is exact, i.e., it does not involve any numerical approximations and integration range truncation errors thanks to the gap options. In practice, however, we only observe a finite number of OTM options and truncate the cosine series expansion with a finite number of terms $N$. We address these issues in the next section.
Although this decomposition is circular (to find the price of a single option we need to know the prices of a continuum of options), it is essential in practice, where we observe only a finite number of option prices but might be interested in pricing options with strikes that are not observed in the market, i.e., we use (ref) to perform the interpolation.
Finally, the call and put price first-order derivatives with respect to the strike price, $\theta_c$ and $\theta_p$, are not directly observable in the market. However, we can approximate them using, e.g., the finite-difference approach. Alternatively, as we show in the next section, we can estimate them from the observed cross-section of option contracts as they are linearly loaded on the call prices. After having estimated these derivatives, we can use them to non-parametrically estimate the RND and option sensitivities.
In a similar model-free way, we can replicate some option sensitivities such as the option delta, $\delta$, the first derivative of the option price with respect to the underlying asset price $S_0$. For that, we first note that
where we additionally assume that $\frac{\partial S_T}{\partial S_0} = \frac{S_T}{S_0}$, i.e.\ the solution to the stochastic process $S_T$ is homogenous of degree one as a function of initial stock price $S_0$.
Let us further denote $B_m := e^{-rT} \mathbb{E}^{\mathbb{Q}}\left[ \sin\left(u_m \log \frac{S_T}{\alpha} \right) \mathbf{1}_{\{\alpha \leq S_T \leq \beta\}} \right]$. Then, $B_m$ can be replicated similar to the cosine coefficients terms $A_m$ using Proposition (ref):
with
Hence, the delta of the call option restricted to the interval $[\alpha, \beta]$ is given by
Furthermore, the delta of the gap call option for $x < \beta$ can be expressed as
which does not depend on $x$. Therefore, combining two deltas, the implied delta for a European call is given as following:
As the RND expansion (ref) and option evaluation (ref), the replication result for the delta (ref) is exact and model-free for all strike values $x \in [\alpha, \beta]$. A similar spanning result can be derived for option gamma, the second order sensitivity. However, one can also easily obtain gamma by properly scaling the option-implied RND (ref) (see, e.g., \citeNP{bates2005hedging} and \citeNP{alexander2007model}).
It is important to emphasize that all three spanning results (given by equations (ref), (ref), and (ref)) are exact, and the restriction to the finite interval does not lead to the associated truncation errors due to the usage of the gap call option. Furthermore, the choice of the interval is data-driven. This is in contrast to the traditional COS method, where one has to select wide interval to minimize the integration range truncation error.
All three option-implied quantities depend on a continuum number of option contracts. In practice, however, we observe only a finite number of option contracts. Nevertheless, all three expressions are easy to approximate using a limited number of observable option prices. Additionally taking into account the observation errors, in the next section, we develop feasible non-parametric estimators for the RND, option prices, and deltas.
In this section, we introduce the computationally feasible estimators based on the option-implied COS valuation method and derive their asymptotic properties.
Unlike the COS method, our approach does not require a parametric specification of model dynamics under the risk-neutral measure. Instead, we use a cross-section of option prices observed in the market. However, these options are observed across a finite set of strike prices and prone to observation errors due to, e.g., bid-ask spread, tick sizes of quotes, and liquidity issues. Therefore, we first describe the observation option scheme.
Our data consists of $n$ OTM option prices observed at time $t=0$ and expiring at a fixed time $T>0$ with a deterministic sequence of strike prices: \[ 0< \underline{K}:= K_1 < K_2 <\dots < K_n =: \overline{K} <\infty. \] For the asymptotic analysis developed below, we assume that the smallest and the largest strike prices, $\underline{K}$ and $\overline{K}$, are fixed, and the number of options $n$ with strike prices between them goes to infinity by shrinking the strike mesh. We could also consider the joint asymptotic scheme as in \citeA{todorov2019nonparametric} and \citeA{blv}, where $\underline{K} \to 0$ and $ \overline{K} \to \infty$ at certain rates as $n\to \infty$. However, we fix $\alpha=\underline{K}$ to be the smallest strike price and $\beta = \overline{K}$ to be the largest strike price since the choice of the interval $[\alpha, \beta]$ in the spanning results does not introduce an integration range truncation error, as discussed in the previous section.
Assumption (ref) can be relaxed to allow for a non-equidistant grid. However, using an equidistant grid simplifies our analysis and enables us to use various numerical approximation methods under one umbrella. Table (ref) summarizes several popular numerical integration methods for an equidistant grid. Later, we comment on our results with a non-equidistant grid of strike prices.
We emphasize again that we do not require an increasing range of strike prices since we account for the restriction to a finite interval. This is different from other methods that estimate RNDs using expansion series, such as in \citeA{lu2021sieve} and \citeA{cui2021model}. Fixing the interval $[\alpha, \beta]$ comes as a great advantage in practice since it allows us to choose different intervals of strike prices without introducing additional approximation errors. For instance, we can consider an interval with only actively traded options.
As common in the option pricing literature, Assumption (ref) imposes an additive error structure form with independent but possibly heteroskedastic error terms (see, e.g.,\ \citeNP{andersen2015parametric}, \citeNP{todorov2019nonparametric}, \citeNP{blv}). The independence assumption can be further relaxed by considering a spatial dependence as in \citeA{andersen2021spatial} at the cost of more complex expressions for the limiting distributions. This would, however, play a secondary role in the developed estimation procedure.
Note that we drop the null index to denote the observed OTM prices. Furthermore, due to the put-call parity, the same observation errors translate into the counterpart in-the-money contracts, i.e., both $ C(K_i) = C_0(K_i) + \varepsilon_i$ and $P(K_i) = P_0(K_i) + \varepsilon_i$ for the call and put contracts with the same strike price $K_i$ and error term $\varepsilon_i$.
Using $n$ observable option prices, we can estimate (part of the) cosine coefficients $D_m$ defined in equation (ref), by using a numerical approximation of the integral:
where $w_i$ are the coefficients of a chosen numerical integration method, as listed in Table (ref).
The deviation of the estimated cosine expansion coefficient $\widehat{D}_m$ from its true value $D_m$ stems from the observation and discretization errors. These errors also arise in the VIX calculation (see, e.g., \citeNP{jiang2005model} and \citeNP{jiang2007extracting}). \citeA{todorov2019nonparametric} and \citeA{blv} also analyze these errors in their estimation procedures along with the truncation errors that arise due to integration over a finite interval. However, in our setting, there are no truncation errors for the cosine coefficients $D_m$ since we take this truncation further into account. See the discussion in Section (ref).
Given the fixed smallest and largest strike prices $\alpha$ and $\beta$, we have the following asymptotic result for the cosine coefficients.
The proof of Proposition (ref) is provided in Appendix (ref). The proposition states that although the estimator $\widehat{D}_m$ with fixed $m$ based on a finite number of option prices is biased, it is asymptotically unbiased as the number of option prices $n$ increases. This serves as a building block for the non-parametric estimators introduced below.
Next, we introduce the computationally feasible option-implied call price estimator $\widehat{\overline{C}}(x)$ of the error-free counterpart $\overline{C}_0(x)$, defined in equation (ref), with a strike $x$ and the payoff restricted to the interval $[\alpha, \beta]$. It can be expressed as a linear combination of asymptotically unbiased estimators $\widehat{D}_m$ with $m= 1, \dots, N-1$ as follows:
where $N$ is the number of expansion terms in the Fourier-cosine expansion and $\widehat{D}_0 = e^{-rT}$.
Unlike its error-free counterpart $\overline{C}_0(x)$, the option-implied call price estimator $\widehat{\overline{C}}(x)$ is based on a finite number of noisy option prices and is prone to three types of errors. These errors can be decomposed as follows:
where $\xi(x),\ \zeta(x),$ and $\overline{\eta}(x)$ are observation, discretization and series truncation errors, respectively, all formally defined in Appendix (ref). The series truncation error $\overline{\eta}(x)$ refers to the truncation of the cosine expansion to the finite number of terms $N$. To get an order of this truncation error, we additionally impose the assumption on the smoothness of the RND.
The proof of Proposition (ref) is provided in Appendix (ref).
The evaluation of options introduces two types of biases: the discretization error bias $\zeta(x)$ and the series truncation bias $\overline{\eta}(x)$, which arises from the truncation of the cosine expansion as in the original COS method. The former vanishes with an increase in the number of option prices $n$, while the latter decreases with an increase in the number of expansion terms $N$ as in the COS method. The joint asymptotic for $n$ and $N$, with $N$ increasing slower than $\sqrt{n}$, guarantees that the option-implied call price estimator $\widehat{\overline{C}}(x)$ is asymptotically unbiased. However, in a finite setting with noisy option prices, an increase in expansion terms can lead to a higher variance of the estimators. We address this bias-variance trade-off in the next subsection by choosing an optimal number of expansion terms $N^*$.
We also note that, like in the COS method, for a sufficiently large interval $[\alpha, \beta]$, the estimator $\widehat{\overline{C}}(x)$ gives a good approximation for the call price with a payoff unrestricted to this interval, $C_0(x)$. However, to accurately evaluate call options, we also need the first-order derivatives of call and put options $\theta_c$ and $\theta_p$ evaluated at the boundaries of this interval (see equation (ref)). Although we could use finite differences to estimate the first-order derivatives non-parametrically, here we use a simple linear relation of observed option prices on these derivatives instead. In particular, for the observed call price with the strike price $K_i$, we can get the following decomposition:
where the second equality follows from equation (ref) and the third one from the decomposition (ref). Therefore, given $n$ observed option prices, the first derivatives of call and put options $\theta_c$ and $\theta_p$ can be estimated using a simple linear regression of $ C(K_i) - \widehat{\overline{C}}(K_i) - C(\beta)$ on $Z_{c}^N(K_i)$ and $Z_{p}^N(K_i)$ with an intercept $\bar{\theta}$, where $Z_{c}^N(K_i)$ and $Z_{p}^N(K_i)$ are the partial sum counterparts of $Z_{c}(K_i)$ and $Z_{p}(K_i)$, respectively, defined in (ref). Adding an intercept into this regression reduces finite-sample biases due to the discretization and truncation errors.
Finally, the option-implied call price estimator for any strike price $x \in [\alpha, \beta]$ is given by
where $\widehat{\bar{\theta}},\ \widehat\theta_c$ and $\widehat\theta_p$ are the OLS estimates of the aforementioned regression. We emphasize here that while the OLS estimates are obtained using a finite number of observed option prices, the estimator (ref) can be obtained for any strike price $x \in [\alpha, \beta]$. Therefore, the call price estimator $\widehat{C}(x)$ can be seen as an interpolation-approximation method, for which we have the following asymptotic result.
The proof of Proposition (ref) can be found in Appendix (ref). Like the estimator $\widehat{\overline{C}}(x)$, the call price estimator $\widehat{C}(x)$ is asymptotically unbiased when the number of expansion terms grows slower than $\sqrt{n}$.
The developed call price estimator is related to non-parametric kernel smoothing methods that are widely used in the literature (see, e.g., \citeNP{ait1998nonparametric}, \citeNP{grith2012nonparametric}, \citeNP{dalderop2020nonparametric}), but it can be considered as a `global' smoother. While kernel methods are typically local smoothers (bandwidth parameters control the locality of these estimators), our call price estimator uses all available option prices via the portfolio spanning result (ref) discussed in Section (ref). This difference allows our method to provide a more flexible approximation of option prices. In the simulation section, we compare these two approaches and demonstrate the superiority of our method.
To demonstrate the `global' nature of our approach, in Figure (ref) we display the weights of option portfolios for the at-the-money call option for different number of expansion terms $N$. The illustration is based on the Black-Scholes model and the simulation set-up is outlined in Section (ref). As shown in the figure, our call price estimator for the strike price $K=F_0$ utilizes all available option contracts instead of restricting information to a local neighborhood. When the number of expansion terms increases, the weights concentrate more around the target strike price while still incorporating information from all available contracts.
After having estimated the cosine coefficients $\widehat{D}_m$ and the first order derivatives $\widehat{\theta}_c$ and $\widehat{\theta}_p$, we can get the non-parametric estimator for the RND of the log price:
where $\nu_f = \frac{2e^{rT}}{\log\left(\beta/\alpha \right) }$. The RND of the future price $S_T$ is obtained by the appropriate transformation of the log price density (ref). A similar asymptotic result carries over to the non-parametric RND estimator.
The proof of Proposition (ref) can be found in Appendix (ref) as well. Unlike the option-implied call price estimator, here we require the number of expansion terms $N$ to grow slower than $n^{-1/6}$ due to the different cosine function. Nevertheless, when this condition is met, the option-implied RND remains asymptotically unbiased.
We emphasize again that this estimator is not for the truncated density, but for the full RND evaluated at any $y$ within the interval $[\log\alpha, \log\beta]$. If one wishes to estimate the RND outside of this interval, additional, often parametric assumptions have to be made about the behavior of the density of options in areas where no option prices are observed. For instance, one possible approach to estimating the RND outside of this interval is to extrapolate option prices beyond the observable range of strike prices using a parametric form based on no-arbitrage conditions. This extrapolated data can then be used to estimate the RND based on the same estimator (ref). In practice, for sufficiently liquid options, the observed range of strike prices covers almost an entire distribution.
The non-parametric estimator for the option delta can also be derived in a similar way using the spanning result (ref) given a finite number of option prices:
where
An analogous asymptotic result holds for the delta estimator (ref), but with an additional assumption.
Assumption (ref) is required to derive the expansion for the delta given in equation (ref). The same assumption is imposed for the other non-parametric delta estimators (see, e.g., \citeA{bates2005hedging} and \citeA{alexander2007model}).
It is worth noting, that in the current formulations, all limiting results are self-scaling. This implies that equidistant strike price Assumption (ref) can be easily relaxed without affecting the limiting distributions.
As discussed earlier, the number of expansion terms $N$ controls the bias-variance tradeoff in the developed non-parametric estimators. An increase in the number of terms reduces the bias resulting from the series truncation error, but increases the variance of the estimators. To find the optimal number of expansion terms $N$ for the Fourier-cosine expansion in the iCOS method, we consider the expansion for the RND\footnote{Depending on the purposes, one could also determine the optimal $N$ that minimizes the difference between the observed and estimated option prices. However, since option prices are observed with noise, this approach can potentially lead to severe arbitrage violations.}.
To assess the impact of truncation on the Fourier-cosine expansion, it is convenient to consider the fit of the density based on the Mean Integrated Squared Error (MISE), defined as
where $\widehat{f}(y)$ is the density estimate based on $N$ expansion terms. Following \citeA{leitao2018data} and \citeA{kronmal1968estimation}, the MISE can be decomposed as follows:
where $A_m$ is the $m$-th Fourier-cosine coefficient, and $\widehat{A}_m$ is its estimate. Since the discretization errors contribute to the asymptotically vanishing bias, we consider the Asymptotic MISE (AMISE), where the second moment equals the variance of $A_m$. Hence, the optimal number of expansion terms $N$ trades off the bias, given by the first part, and the variance of the estimator.
We can derive a recursive relationship in $N$ as follows:
from which we can see that if $A_N^2 - \mbox{Var}\left( \widehat{A}_N \right) >0$, then $\mbox{AMISE}_{N+1} < \mbox{AMISE}_N$. We can use this inequality as a rule to determine the optimal number of expansion terms.
The variance of the cosine coefficient estimators depends on the estimates $\widehat{D}_N$ and $\widehat{\theta}$, and is provided in a closed-form in Appendix (ref). A true value of $A_N$ is, however, unknown a priori. Furthermore, due to the presence of discretization error bias in finite settings, this inequality may not accurately represent the MISE. Nevertheless, we can operationalize this inequality to obtain a rule-of-thumb for the optimal number of expansion terms $N$ given the feasible estimates $\widehat{A}_N$. We provide such a rule-of-thumb algorithm in Appendix (ref).
Finally, we note that the MISE given by equation (ref) is for the Fourier-cosine expansion. Hence, the delta estimator, which essentially utilizes the Fourier-sine expansion, may require a different\footnote{The sine expansion is known to have a slower rate of convergence than the cosine series. In fact, this is the main reason for popularity of the Fourier cosine expansions rather than the Fourier or sine series.} optimal choice of expansion terms $\widetilde{N}$. In this case, a similar rule can be applied but with the sine coefficients $\widehat{B}_N$ instead.
In this section, we investigate the finite-sample performance of the developed non-parametric estimators. We consider two models to generate data: the \citeA{black1973pricing} model and the `double-jump' stochastic volatility model of \citeA{DPS2000}. The former has closed-form solutions for the true quantities of interests, while the latter offers a more realistic depiction of options data that features two stylized facts -- stochastic volatility and jump components in returns and volatility.
For each model, we set the initial spot price $S_0=4000$, the interest rate $r=0$, and the strike prices between 85% and 110% of the spot price with equidistant increments of 5, similar to the available S&P 500 index option data. This results in $n=201$ option contracts for each maturity.
We distort the true option prices of each model by adding homoskedastic observation errors, i.e.,
where $\epsilon$ is an i.i.d.\ standard normal random variable. This error structure roughly matches the dispersion of errors in empirical applications, where option errors are typically within the tick size of \$0.05. Note that the smallest and the largest strike prices are fixed and the corresponding OTM options are strictly positive.
First, we consider the Black-Scholes model with short (30 days) and long (1 year) maturities, and set the volatility parameter $\sigma=0.3$. The true option prices are generated via the Black-Scholes formula and then distorted with the additive error terms as described above. Table (ref) provides the simulation results of the estimated option call prices for a selection of strike prices, along with the estimates of $\boldsymbol{\widehat{\theta}}$. The latter includes the intercept $\bar{\theta}$ and the first-order derivatives $\theta_c$ and $\theta_p$, which have closed-form solutions in the Black-Scholes model. The number of expansion terms is set to $N=14$ for short maturity options and to $N=7$ for long maturity options. This choice is motivated by the rule-of-thumb discussed in Section (ref) and Appendix (ref). The numerical integration scheme is set to Simpson's 1/3 rule throughout the simulations.
Although we estimate option prices from the set of OTM prices themselves, the simulation results indicate the convergence of the considered approach. Furthermore, the Monte Carlo results show no significant bias for all levels of strikes and maturities and a reduction in the variance of option prices. In fact, the Monte Carlo standard deviations of the estimated prices are smaller than the standard deviations used to simulate option errors, indicating the smoothing effect of the estimation procedure. The asymptotic standard deviations, defined as the square root of the average estimated asymptotic variance, roughly correspond to the Monte Carlo standard errors, which indicates the validity of the constructed standard errors.
The estimated parameters $\boldsymbol{\widehat{\theta}}$ also exhibit good finite-sample performance. The estimated intercept $\bar{\theta}$, which collects the average of finite-sample biases, indicates that these errors are of rather small order. The good finite-sample performance of the first-order derivatives $\theta_c$ and $\theta_p$ is crucial for the RND and delta estimators considered below.
The Monte Carlo simulation results for the RND estimates are reported in Table (ref). Note that the RND estimator is for the log price $\log S_T$ and evaluated at a few strike levels after the corresponding log transformation. The true RND for the Black-Scholes model is the normal distribution with mean $\log S_0 - \tfrac{1}{2}\sigma^2 T$ and variance $\sigma^2 T$. Similar to the option price estimates, the estimated RNDs exhibit good finite-sample performance. The asymptotic standard errors roughly match the Monte Carlo standard errors, and both tend to increase towards the bounds of the considered interval.
Finally, the simulation results for the delta estimates are reported in Table (ref). The true values are obtained in the closed-form under the Black-Scholes assumptions. The number of expansion terms for the delta estimator is set to $\widetilde{N}=25$ since the sine series exhibits a slower convergence rate (see discussion in Section (ref)). Unlike the RND and call price estimates, the estimation of the delta exhibits small bias terms, which are economically likely to be negligible.
The `double-jump' stochastic volatility model of \citeA{DPS2000}, labeled as SVCJ, allows for stochastic volatility and jumps in returns and volatility, and under the risk-neutral measure $\mathbb{Q}$ is given by the following system of stochastic differential equations:
where two Brownian motions $W_1$ and $W_2$ are correlated with coefficient $\rho$, and $N_t$ is a Poisson jump process with intensity $\lambda$. The jump sizes in returns are Gaussian, $J \sim \mathcal{N}(\mu_j, \sigma_j^2)$ with the expected relative jump size in returns $\mu = \exp(\mu_j + \frac{1}{2}\sigma_j^2) - 1$, while the co-jump sizes in volatility are exponentially distributed, $J^v \sim \exp(1/\mu_v)$, and independent of jump sizes in returns. We choose the following parameter values for the simulation:
Since there is no closed-form solution for option prices under the SVCJ model, we simulate them using the COS method. We use the analytic solution for the CCF and a large number of expansion terms $N=1024$ with $[a,b] = [-4\sqrt{T}, 4\sqrt{T}]$. We then add observation errors to these option prices as described previously.
Figure (ref) displays the simulated option prices from the SVCJ model on BSIV and log OTM spaces for $T=30$ days. The chosen model parameters generate the so-called implied volatility `smile', which is commonly observed in the market, particularly for short-dated options. Capturing such pronounced `smiles' can be challenging for many parametric and non-parametric methods since they require the methods to be rather flexible. As a consequence, these methods often fail to accurately capture option prices, RND, and deltas.
For the SVCJ model, we compare the simulation results of the developed approach with the closest non-parametric and widely-used alternative, the kernel smoother. In fact, kernel smoothing methods are also model-free and do not require any optimization routines. In particular, we consider the Nadaraya–Watson kernel estimator with the Gaussian kernel applied to the BSIV space, as in, e.g., \citeA{ait1998nonparametric} and \citeA{grith2012nonparametric}. After smoothing BSIV observations, we convert them into price dimension to obtain call price estimates. We then calculate the second-order derivatives to obtain the RND estimates due to \citeA{breeden1978prices}. Fitting option prices on implied volatility space is commonly used in practice (see, e.g., \citeA{ait1998nonparametric}, \citeA{andersen2015parametric} among many others).
Finding the bandwidth parameter $h$ is crucial for the kernel smoothing methods as it controls the bias-variance tradeoff. Since we are interested in estimating both option prices and the RND, we consider the kernel bandwidths, as in \citeA{ait1998nonparametric}, of the following form:
with $p=2$ and some constant $c>0$. In practice, one can use the cross-validation to find the optimal bandwidth, but in simulations we vary the constant $c$ to illustrate the bias-variance tradeoff.
Table (ref) provides simulation results for the call price estimates under the simulated SVCJ model. We compare the results obtained using our proposed iCOS method with the widely-used kernel smoothing approach with different bandwidth parameters. First, we observe that the simulation results for the call price estimates under the SVCJ model using our approach are similar to the results based on the Black-Scholes model discussed in the previous subsection. Second, as expected, the biases for the kernel smoother decline as the bandwidth parameter decreases, but this comes at the cost of increased variance. Notably, the biases for the kernel smoother are especially pronounced at the moneyness level of 1.05, which roughly corresponds to the `turning' point of the smile depicted in Figure (ref). Finally, when comparing two methods, we notice that only the results with the parameter $c=0.03$ for the kernel smoother are comparable to the iCOS approach in terms of biases. However, such a small bandwidth value results in a non-smooth RND as we discuss below.
Table (ref) provides the Monte Carlo simulation results for the RND of log future prices. We again compare our approach with the kernel smoothing method using the same set of bandwidth parameters. The developed iCOS approach demonstrates a good finite-sample performance, similar to the Black-Scholes case. In contrast, the results for the kernel smoothing approach are often biased or exhibit large variances, which is an indication of non-smooth and, hence, arbitrage-violating RND estimates.
Analogously, Table (ref) provides the Monte Carlo results for the call delta estimates under the SVCJ model. The non-parametric iCOS method yields insignificant biases but larger variances than the kernel smoothing method. The latter, however, again exhibits biases at the moneyness level of 1.05 except for the parameter $c=0.03$, which corresponds to a non-smooth RND.
Overall, when comparing two approaches, we observe that the kernel smoothing method fails to fully capture the shape of the observed option data, which results in the biased estimates of the call prices, RND, and deltas. Decreasing the bandwidth parameter reduces the biases but at the cost of a non-smooth RND with large variance. In contrast, the iCOS method is able to simultaneously capture the shape of the observed option data, RND and deltas with insignificant biases and small variances.
In this section, we first demonstrate the proposed estimation procedures in two empirical illustrations. Then, we apply the developed method to analyse errors in the VIX index.
In the first empirical illustration, we consider options on the S&P 500 stock market index (SPX) obtained from the Chicago Board Options Exchange (CBOE), which are commonly used in the literature. We consider a snapshot of options at 3:45 pm ET with a time-to-maturity $T = 29$ days traded on April 1, 2021. The forward price implied from the put-call parity is $F = \$ 4008.5$. For the developed estimation procedure, we utilize the mid-quote prices of OTM contracts. Additionally, we compare the pricing accuracy of our approach with the bid-ask spread of the corresponding contracts, which allows us to assess the performance of our method in the real-world options market context.
We focus our analysis on the interval $[\alpha, \beta] = [2950, 4400]$, i.e., we use options with strikes from this interval. Although there are a few option contracts with strike prices beyond this range, they tend to be less liquid and spaced further apart\footnote{ The distance between strike prices outside the considered interval is \$50 or \$100, while the distance for close to ATM options is only \$5.}. We do not filter out any options, except for those with zero bid prices. This results in a total number of 239 OTM options with non-equidistant strike prices.
To approximate the integrals in the option-implied cosine coefficients (ref), we use Simpson's 1/3 rule. Figure (ref) displays the estimated coefficients $\widehat{D}_m$ and $\widehat{A}_m$ plotted against the term number $m$. The coefficients $\widehat{A}_m$ are displayed after applying a logarithmic transformation, along with their corresponding standard deviations $\sigma_A$. As expected, these coefficients converge towards zero up to a certain point but then exhibit a divergent pattern due to the discretization errors and increased variance. The rule-of-thumb algorithm, detailed in Appendix (ref), selects $N^{*}=23$ as the optimal number of terms. This number roughly corresponds to the point where the cosine coefficients stabilize and their standard deviations surpass the magnitude of the coefficients themselves. We use this number of terms for the subsequent analysis.
Figure (ref) depicts the option-implied call price estimates $\widehat{C}(x)$ given by equation (ref) for the considered data. The figure displays the market prices in terms of the BSIV and the logarithm of OTM prices, although the estimation is performed in terms of dollar-amount OTM prices. As observed, the estimates closely capture the shape of the implied volatility smile and the log OTM prices. Note that the solid line does not pass exactly through all market prices, but provides an approximation of them.
To investigate the accuracy of our estimation procedure, we plot in Figure (ref)(a) the pricing errors as the difference between the call option price estimates $\widehat{C}(x)$ and the market observed call prices $C(x)$ for $x\in\{K_1, \dots, K_n\}$. Most of these differences range between \$-0.05 and \$0.05, corresponding to two ticks of size \$0.05 in this dataset. Notably, these differences exhibit a nearly homoskedastic pattern with respect to strike prices. In other words, the error terms do not vary much with the price level of the options and strike prices, contrary to what is often assumed in the literature. It suggests that the primary source of observation errors for these highly liquid SPX options can be attributed to the minimum tick size.
Panel (b) in Figure (ref) displays the same pricing differences divided by the half-spread, which is defined for each option contract with strike price $x$ as \[ HS(x) = \frac{AskO(x) - Bid O(x)}{2}, \] where $AskO(x)$ and $Bid O(x)$ are the ask and bid prices for the OTM option with strike price $x$, respectively. The half-spread is calculated based on put option quotes if $x < F$ and on call options otherwise. Since we use mid-quote prices $O(x) = \frac{AskO(x) + BidO(x)}{2}$ for our estimation, these pricing errors indicate the percentage distance between the mid-quote price and the bid (if negative) or ask (if positive) prices. As shown in Figure (ref)(b), most pricing errors fall within half of the bid-ask spread, indicating a very good pricing performance of the method.
Finally, Figure (ref) plots the non-parametric estimates for the RND of future asset price and call deltas for the SPX options data. The RND estimate of a future price is obtained from the RND estimate of a log future price given by equation (ref). Given the derived asymptotic theory (and the appropriate delta-rule), the confidence interval for the RND is displayed in a gray area. It is very narrow in the main part of the distribution and slightly widens towards the ends of the considered interval. We also note that we do not impose any arbitrage-free conditions on our RND and option price estimators. Thus, the RND estimates have negative values for some values of strikes. However, such minor arbitrage violations are unlikely to have any practical implications since all corresponding call price estimates fall within the minimum tick size and bid-ask spread.
The estimated call deltas $\widehat{\delta}$ are displayed alongside deltas based on the Black-Scholes model, $\delta_{BS}$. As observed, the Black-Scholes deltas can substantially underestimate the in-the-money call deltas. This might potentially result in hedging errors as discussed in \citeA{bates2005hedging} and \citeA{alexander2007model}.
In our second application, we examine equity options on Amazon with a very short time-to-maturity of $T = 1$ day. These options are traded on the Earning Announcement Day (EAD) of April 26, 2018, prior to the announcement itself. Compared to the SPX options, Amazon options are less liquid and are prone to larger observation errors due to their very short maturity. Moreover, the EAD introduces extra uncertainty about the stock price at the expiry.
Similar to the SPX options, we use mid-quote prices of OTM contracts and filter out only zero-bid contracts. We concentrate our analysis on the interval $[\alpha, \beta] = [1250, 1760]$, which corresponds to approximately 18% below and 16% above the underlying spot price of \$1518.96 on this EAD. The availability of such a wide interval for short-dated options is attributed to the information uncertainty surrounding the EAD. Based on the rule-of-thumb for the optimal number of expansion terms, we set $N = 13$.
Figure (ref) presents the estimation result for the option-implied call prices displayed on BSIV and log OTM spaces. Notably, the BSIVs are exceptionally high, reaching approximately 130% for ATM options expiring in just one day. Furthermore, these BSIVs exhibit a distinctive W-shaped pattern, which is atypical for implied volatility curves\footnote{Note that most parametric curves amd models commonly used in the literature would fail to capture this pattern, leading to large estimated errors.}. The W-shape arises from the anticipation of a significant stock price jump following the earnings announcement release. \citeA{alexiou2021pricing} document frequent concave patterns in implied volatilities prior to the EAD for equity options.
Additionally, we observe a large dispersion of option prices. However, our estimation procedure effectively smoothes out the noisy data, resulting in accurate price estimates. For this example, in Figure (ref) we also display the 95% confidence interval around the estimated call prices, obtained by applying the appropriate delta rules for the derived asymptotic results. We emphasize that this confidence interval reflects the uncertainty around the estimates and not the observation errors in option prices. Therefore, it does not and need not cover the observed prices.
Figure (ref) displays the pricing errors for call price estimates. As with the SPX options, we plot the pricing errors and the errors relative to the half-spread. Consistent with Figure (ref), the pricing errors are larger than those for the SPX options but are still centered around zero.
Finally, Figure (ref) displays the estimated RND and deltas for Amazon options. The consequence of the W-shaped implied volatility curve is a bimodal RND, reflecting the market's anticipation of two possible outcomes. The two modes of the estimated RND are at \$1442 and \$1590, which corresponds to around 5% down and 4.7% up from the spot underlying stock price, respectively. After the announcement, the next day's opening price for Amazon was \$1634 and it closed at \$1574 (3.62% up from the spot price). The bimodality of the RND is also reflected in the estimated deltas. As shown in Figure (ref)(b), the Black-Scholes deltas overestimate the deltas for the strikes around the first mode and underestimate them for strike prices close to the second mode.
The developed methodology allows us to analyse errors embedded in the VIX index. As noted by \citeA{jiang2007extracting}, the construction of the VIX is prone to several types of approximation errors, including truncation and discretization errors. The former arises from truncating the real line to the range of observed strike prices, and the latter is due to the discreteness of strike prices. On top of that, option prices used in the VIX construction are subject to observation errors since the true prices are not observed perfectly, as argued in Section (ref). This results in observation errors in the VIX index. Our methodology enables us to estimate and disentangle observation and discretization errors in the VIX.
In particular, the CBOE calculates the VIX index as\footnote{For simplicity of notation, the exposition is based on a single maturity of 30 days. The CBOE averages (in total variances) the two VIX measures constructed using the near-term and the next-term options. In our empirical application, we follow the same procedure.}
where $K_0$ is the largest strike price below the forward level $F$, $\Delta K_i = \frac{1}{2}(K_{i+1} - K_{i-1})$ for $i=2,\dots,n{-}1$, and $T=30$ days. For more details, see the \citeA{CBOE} white paper. Since the OTM option prices are observed with noise, the VIX itself contains measurement error. Given the consistent estimator of option prices $\widehat{O}(K_i)$, we can (re-)construct the VIX using these estimates and obtain estimates of the observation errors in the VIX. That is, we calculate
and define $\widehat{\xi}_{\mbox{vix}} := \mbox{VIX} - \widehat{\mbox{VIX}}$ as the estimator of the observation error in the VIX index.
On the other side, the VIX is developed to approximate the model-free implied volatility. Since our method allows us to further consistently interpolate between observed strike prices, we can construct the measure of the model-free corridor implied volatility (CIV) as
The VIX can be seen as a measure of CIV with barriers fixed at the lowest and highest strike prices that the CBOE uses for calculating the index (\citeNP{andersen2015exploring}). Therefore, we define $\widehat{\zeta}_{\mbox{vix}} := \widehat{\mbox{VIX}} - \widehat{\mbox{CIV}}$ as the estimator of the discretization error in the VIX index.
To estimate the observation and discretization errors in the VIX, we consider the SPX options obtained from the CBOE from January 3, 2017 until April 1, 2021. We follow the exact same procedure for the construction of the index as outlined in their white paper (\citeNP{CBOE}). To reduce the finite-sample bias in the iCOS procedure due to the discreteness of the observed option strikes, for each tenor, we interpolate option prices using cubic splines applied to implied volatilities. This can be seen as a bias-reduction procedure as motivated in \citeA{blv} in the context of option-implied CCFs.
Figure (ref) displays the time series plots of the estimated observation and discretization errors in the VIX over the period of more than four years. Figure (ref) complements it with histograms of the percentage errors over the same time period. We notice that the observation errors are centered around zero, while the discretization errors are mainly positive. This is expected since observation errors in option prices do not introduce biases in the VIX, while the discreteness of the strikes leads to a finite-sample bias in the constructed index. In fact, the sample average of the percentage observation errors is nearly zero, and the average of the percentage discretization errors is estimated at around 0.135%. Furthermore, the magnitude of observation errors is rather low, reaching in absolute terms up to 0.04 percentage points. The discretization errors, on the other hand, can result in a substantial overestimation of the index, with the differences up to 0.7 percentage points during high volatility periods.
In this paper, we proposed a non-parametric estimation procedure for option prices, RND, and option sensitivities. This method is based on the combination of Fourier-based cosine technique and the option spanning result of \citeA{carr2001optimal}. This combination allows for a flexible and accurate estimation of the density, option prices and option sensitivities without imposing parametric assumptions on the dynamics of underlying asset and on the shape of implied volatility surface. We have also established the asymptotic properties of the proposed estimators and demonstrated the finite sample properties through the Monte Carlo simulations.
The usage of the proposed method is illustrated in empirical applications using options data on the S&P 500 stock market index and Amazon equity options on the Earning Announcement Day. The empirical analysis demonstrates the effectiveness of the iCOS method in accurately estimating option prices and capturing important market features in different market conditions. Additionally, we demonstrated the usefulness of our methodology to dissect and quantify errors in the VIX index, one of most popular measure of market volatility. We found that observation errors in the VIX are centered around zero and have a small magnitude, while discretization errors can lead to positive and substantial biases in the VIX index.