A New Nonparametric Estimate of the Risk-Neutral Density with Applications to Variance Swaps

We develop a new nonparametric approach for estimating the risk-neutral density of asset prices and reformulate its estimation into a double-constrained optimization problem. We evaluate our approach using the S\&P 500 market option prices from 1996 to 2015. A comprehensive cross-validation study shows that our approach outperforms the existing nonparametric quartic B-spline and cubic spline methods, as well as the parametric method based on the Normal Inverse Gaussian distribution. As an application, we use the proposed density estimator to price long-term variance swaps, and the model-implied prices match reasonably well with those of the variance future downloaded from the CBOE website.


Introduction
A financial derivative, such as option, swap, future, or a forward contract, is an asset that is contingent on an underlying asset. Its fair price can be obtained by calculating the expected future payoff under a risk-neutral probability distribution of the future payoff of the underlying asset. Therefore, the problem of pricing a derivative can be addressed via estimating the risk-neutral density of the future payoff of the underlying asset. On the other hand, the market prices of the derivatives traded in a financial market reveal information about the risk-neutral density. Breeden and Litzenberger (1978) was among the first to use option prices to estimate the risk-neural probability distribution of the future payoff of the underlying asset. Among the financial products that can be used to recover the risk-neutral density, European options are the most common ones, which give the investors rights to trade assets at a pre-agreed strike price at the maturity date. Among all the underlying assets that options are written on, Standard & Poor's 500 Index (S&P 500) is a popular one, which aggregates the values of stocks of 500 large companies traded on American stock exchanges and provides a credible view of American stock market for investors.
There are numerous approaches towards recovering risk-neutral density functions in the literature (see, for example, Bliss and Panigirtzoglou (2002) for an extensive review), including parametric and nonparametric ones. Parametric approaches typically specify a parametric statistical model for the risk-neutral density and then recover the parameters by solving an optimization problem. For examples, Jarrow and Rudd (1982) used a lognormal distribution; Melick and Thomas (1997) considered a mixture of lognormal distributions proposed by Ritchey (1990); Sherrick et al. (1992) employed a three-parameter Burr distribution, called the Burr family, which covers a broad range of shapes, including distributions similar to gamma, lognormal, and J-shaped beta. Another commonly used probability distribution in the literature of derivative pricing is the Generalized Hyperbolic distribution that contains Variance Gamma, Normal Inverse Gaussian, and t distributions as special cases (see, for instance, Ghysels and Wang (2014)).
A nonparametric procedure, by contrast, is free from distribution assumptions on the underlying asset and thus achieves more flexibility than parametric methods. For example, Monteiro et al. (2008) used cubic spline functions to model the unknown riskneutral density. An estimated density is numerically obtained by solving a quadratic programming problem with a convex objective function and non-negativity constraints. They deliberately chose more knots than option strikes for higher flexibility. Lee (2014) approximated the risk-neutral cumulative distribution function using Quartic B-splines with power tails and the minimum number of knots that meet zero bid-ask spread. Their estimation was based on out-of-the-money option prices.
In this paper, we propose a simpler but more powerful nonparametric solution using piecewise constant functions to estimate the risk-neutral density. It is easy to implement since the estimating problem is formulated as a weighted least squared procedure. It is more powerful since our method can recover the risk-neutral density more effectively with all available option market prices without screening. Furthermore, our solution provides a practical way to explore profitable investment opportunities in financial markets by comparing the estimated price and the corresponding market price.
The rest of this paper is structured as follows. Section 2 introduces the proposed nonparametric approach after reviewing cubic splines, Quartic B-splines, and the Normal Inverse Gaussian (NIG) parametric approaches in the literature. In Section 3, we run comprehensive cross-validation studies using real data to compare different methods. In Section 4, we apply the proposed nonparametric approach to price variance swap, which is challenging in practice. We conclude in Section 5. The proofs and more tables are collected in the Appendix.

Methodologies
In this section, we first briefly review the cubic splines, quartic B-splines, and NIG approaches in the literature for recovering the risk-neutral density (RND). Then we propose our piecewise constant (PC) nonparametric approach with least square and weighted least square procedures.

Nonnegative cubic spline estimate for RND
Given the current trading date t and the expiration date T of European options, let [K 1 , K q ] be the range of strike prices of all available options traded in the market. Monteiro et al. (2008) considered n + 1 equally spaced knots for a cubic spline with K 1 = x 1 < x 2 < x 3 < · · · < x n < x n+1 = K q . These knots are not necessarily a subset of the available strikes. Nevertheless, the closer these knots are to the strikes, the better their solution is. Monteiro et al. (2008) also claimed that the number of knots should not be too much larger than the number of distinct strikes.
In order to guarantee the estimated RND nonnegative everywhere, Monteiro et al. (2008)'s solution is much more complicated and computationally expensive than the usual cubic spline estimates. For comparison purpose, we keep only the constraints that ensure the nonnegativity of the density function on knots in their optimization procedure. In terms of consistency between the estimated fair prices and the market prices of options, if our approach performs better than the cubic spline estimate with less constraints, then our approach is better than the nonnegative RND estimate of Monteiro et al. (2008).
In order to apply their RND estimate for practice, Monteiro et al. (2008) suggested eliminating option prices that led to potential arbitrage opportunities according to bidask interval, put-call parity, monotonicity, and strict convexity. They also generated "fake" call option prices using put-call parity to eliminate "artificial" arbitrage opportunities. Our comprehensive studies in Section 3 show that their screening and cleaning procedure may result in substantial information loss.

Quartic B-spline estimate
Lee (2014) adopted a uniform quartic B-spline to estimate the risk-neutral cumulative distribution function (CDF). They used power tails to extrapolate the area outside the strike price range. Lee (2014) suggested using only the out-of-the-money (hereafter OTM) options to estimate the CDF, including OTM call option whose strike is higher than the underlying asset price, and OTM put option whose strike is lower than the underlying asset price. The opposite concept is in-the-money (ITM) options. OTM options are typically cheaper than ITM options and are considered to be more liquid as well. Nevertheless, our case studies in Section 3 show that ITM options may help recover the risk-neutral distribution as well.
Due to fewer parameters, the quartic B-spline estimate is computationally more efficient than the nonnegative cubic spline approach. Lee (2014) chose the number of knots needed as the minimum number that satisfies zero bid-ask pricing spread. They also suggested an elimination of options subject to monotonicity and strict convexity constraints.

NIG parametric approach
For comparison purpose, we choose one parametric approach for approximating the riskneutral density, as suggested by Eriksson et al. (2004). It is based on the Normal Inverse Gaussian (NIG) distribution, which belongs to the Generalized Hyperbolic class and can be characterized by its first four moments, i.e., mean, variance, skewness, and kurtosis. According to Bakshi et al. (2003a), these four moments can be estimated by the OTM European call and put options. One major issue with NIG density estimate is that, as shown in Ghysels and Wang (2014), the feasibility of NIG approach drops down as the time to maturity increases, since more estimated skewness and kurtosis pairs fall out of the feasible domain of NIG distribution.

The proposed piecewise constant nonparametric approach
The piecewise constant (PC) approach that we propose in this paper is a simpler but more efficient nonparametric approach. Let S t and S T stand for the current price on day t and the future price on day T of an equity. In order to estimate the risk-neutral density function f Q of log(S T ) conditional on the information up to day t, we propose to use a piecewise constant function, or a step function, to approximate f Q , using all distinct strike prices as knots. The constants in the step function can be estimated by solving an optimization problem with certain constraints. By forcing the constants nonnegative, the nonnegativity of the estimated risk-neutral density can be guaranteed easily.
To be precise, suppose that we have a collection of market prices of European put and call options that are traded on date t and expire on date T . Let {K 1 , K 2 , . . . , K q } stand for the distinct strikes in ascending order, and C be the collection of indices for call options and P for put options. Then C ∪ P = {1, 2, . . . , q}. Let m = |C| and n = |P| be the numbers of call options and put options, respectively. Then m + n ≥ q .
Given the risk-neutral density f Q , the fair prices of put option and call option with strike K i are respectively, where R tT stands for the cumulative risk-free interest rate from t to T , that is, $1 on day t ends for sure with e R tT dollars on day T . We denote by r t the risk-free interest rate for period [t, t + 1], which is obtained from risk-free zero-coupon bond, and clearly R tT = T −1 j=t r j . To account for the risk-neutral density outside the range [K 1 , K q ], we add K 0 = K 1 /c K and K q+1 = c K K q , where c K > 1 is a predetermined constant that can be chosen by means of cross-validation or prior knowledge. See details in Section 3. We then use a piecewise constant function f ∆ to approximate f Q , that is, and zero elsewhere. Here ∆ = {log K 1 , . . . , log K q } stands for the collections of distinct strikes in log scale, and {a l , l = 1, . . . , q + 1} are nonnegative constants satisfying q+1 l=1 a l log K l K l−1 = 1 (2) due to the requirement´+ ∞ −∞ f ∆ (y)dy = 1. Given the approximate risk-neutral density f ∆ , the estimated put and call prices with respectively, which are essentially linear functions of a 1 , . . . , a q .
Proposition 2.1. Given a l ≥ 0, l = 1, . . . , q + 1 satisfying (2), the estimated prices for put and call options with strike K i satisfy where X The proof of Proposition 2.1 is relegated into A.
In order to estimate the unknown parameters a 1 , . . . , a q+1 , a natural solution is to minimize the least square (LS) objective function subject to a l ≥ 0, l = 1, 2, . . . , q + 1, and Equation (2), whereC i andP i are market prices of call option and put option, respectively, with strike K i . When there exists a risk-neutral density f Q , we have C i =C i , i ∈ C and P i =P i , i ∈ P. That is, the market prices are fair when there is no arbitrage in the financial market.
From an investment point of view, the more expensive the options are, the less liquid they are. An alternative approach to determining a 1 , . . . , a q+1 is to minimize a weighted least square (WLS) objective function Since typically OTM options are less expensive and more liquid, the WLS estimate is in favor of OTM options over ITM options.

Pricing European Options
In this section, we use S&P 500 European option market prices to evaluate the performances of different RND estimators by regenerating the fair prices of options.

S&P 500 European option data
The option market price data considered here are European calls and puts written on the S&P 500 indices from January 2, 1996 to August 31, 2015 in the US. The expiration dates of the options are the third Saturday of the delivery month. Following Carr et al. (2012), we keep only the options with positive bid prices, positive volumes, and with expiration more than seven days in our analysis. Similar to Ghysels and Wang (2014)

Comprehensive comparisons with existing methods
In this section, we use the S&P 500 European options to evaluate the performances of the parametric NIG estimate, quartic B-spline (Bspline) estimate, the nonnegative cubic spline estimates with either least square criterion (Cubic + LS) or weighted least square criterion (Cubic + WLS), as well as the proposed piecewise constant estimate with either least square or weighted least square objective function using OTM options only (PC + LS + OTM or PC + WLS + OTM) or using all available options (PC + LS + ALL or PC + WLS + ALL). All the comparisons are based on the performance of recovering option market prices. For each of the seven categories listed in Table 1, we randomly select 200 pairs of (t, T ). For each pair, the market prices of calls and puts are collected. The aforementioned approaches are applied to estimate the RND of the underlying asset at time T . We then use the estimated conditional density to price the OTM and ITM options. Deviations between the market prices and the estimated prices are assessed by calculating the absolute error L a and relative error L r which are defined below In Table 2 and Table 3, we report L a and L r that are averaged over 200 randomly selected pairs of (t, T ) from each time-to-maturity category for each estimation approach. The columns labeled "200" show the actual number of pairs that yield a valid RND estimate. The higher the count is, the more effective the method is. As explained in Section 2.3, the NIG approach is very picky in selecting calls and puts. For B-spline and Cubic methods, following the same filtering procedures as in Monteiro et al. (2008) and Lee (2014), we observe that fewer options become available as the time-to-maturity increases, which results in substantial information loss. On the contrary, our PC methods with LS or WLS are feasible for almost all cases, especially when using both ITM and OTM options.
In terms of the averaged absolute error L a and relative error L r for different combinations of time-to-maturities and RND estimates, our PC estimates are more stable and accurate than the other three existing approaches. As illustrated in Tables 2 and  3, the proposed PC methods always yield the lowest L a or L r , regardless of the type of options used. Nevertheless, it should be noted that we do not include the scenario when solely using ITM options. This is because its performance is not as good as the one with all options included. Thus, when it comes to practical implementation, we would recommend using all available option prices, including both ITM and OTM options.
For practical uses, if the goal is to obtain the most precise price, we recommend "PC+LS+ALL", which controls absolute error L a the best; if the goal is for best investment, we recommend "PC+WLS+ALL", which controls relative error L r the best.

Consistency of PC estimates for fair prices
Given the distinct strike prices K 1 < K 2 < · · · < K q with the corresponding market prices {C i , i ∈ C} of call options and {P i , i ∈ P} of put options traded on date t with expiration date T , if a risk-neutral density f Q of log S T exists, theñ That is, the market prices (C i ,P i ) agree with the fair prices (C i , P i ) .
Recall that the proposed PC approach provides an approximation   for the RND f Q with (a 1 , . . . , a q+1 ) minimizing absolute error L(a 1 , . . . , a q+1 ) or relative error W (a 1 , . . . , a q+1 ). The estimated fair prices based on f ∆ are (P i ,Ĉ i ) determined by equations (3) and (4), respectively. A natural question is that how close f ∆ is to f Q . Unfortunately, the question is not valid since f Q is typically not unique in practice. In other words, the market is not complete enough to determine a unique f Q .
An alternative question is whether the fair prices based on f ∆ could recover the market prices well. From the comprehensive studies reported in Tables 2 and 3, the answer seems to be yes. The following theorem justifies the answer with a large number of options and dense strike prices.
Theorem 3.1. Suppose there exits a continuous risk-neutral density f Q of log S T satisfy-ing´∞ 0 e x f Q (x)dx < ∞. Let ∆ = {log K 1 , . . . , log K q } be the collection of distinct strike prices in log scale with both call and put option market prices available. Then as K 1 → 0, K q → ∞, q → ∞, and |∆| := max 1≤i<q log(

is a necessary and sufficient condition forC i < ∞.
Remark 3.2. The proof for Theorem 3.1 is relegated to B. It shows the existence of (a 1 , . . . , a q+1 ) such that max 1≤i≤q |Ĉ i −C i | < , max 1≤i≤q |P i −P i | < for any given > 0 when K 1 , |∆| are sufficiently small and K q , q are sufficiently large. In other words, |Ĉ i −C i |, |P i −P i |, i = 1, . . . , q, can be uniformly small.

Detecting profitable opportunities
Theorem 3.1 guarantees the consistency of the proposed PC method given the assumption of the existence of a continuous risk-neutral density. Nevertheless, the PC method can still be applied even when there is any arbitrage opportunity in the market. In this case, a significant difference between some market price and its estimated fair price is expected.
With a given set of market option prices, our nonparametric method can recover a fair option price with any strike price. From an investment point of view, we are able to detect options on the markets that are under or over priced. It may not be adequate to claim arbitrage opportunities due to the lack of guarantee to earn and there is a mature market system designed to catch such kind of difference among the option prices. Nevertheless, we can still report profitable investment opportunities for investors.
In Figure 1, we provide an illustrative example using m + n = 95 available market prices of options traded on 04/14/2014 with expiration 5/9/2014. For each of the 95 options, we obtain its fair price based our PC+LS method using the market prices of the rest m + n − 1 = 94 options. Then we compare the market price and the leave-one-out fair price, known as leave-one-out cross validation. Figure 1 market prices in dots and leave-one-out fair prices in solid line against the corresponding strike prices. It seems that they match each other very well.
In order to check the difference between market price and fair price, we plot the relative difference, that is, (market price -fair price)/fair price, against strike price in Figure 1(b). Then the sign of the difference between the market price and the fair price indicates whether the option is under or over priced. More than that, in order to check if the difference between a market price and its fair price is statistically significant, we bootstrap the rest market prices 50 times to obtain a 95% confidence interval of the fair price. The dash lines in Figure 1(b) show the upper and lower ends of the bootstrap confidence intervals. When a market price falls outside its bootstrap confidence interval, we may report to investors that the corresponding option is significantly under/over priced compared with the market prices of the other options.

Pricing Variance Swap
With an estimated risk-neutral density, one can calculate the fair price of any derivative whose payoff is a function of S T . In this section, we apply the proposed method to price variance swaps. Our study shows that our fair prices match the market prices of long-term variance swaps reasonably well.
A variance swap is a financial product that allows investors to trade realized variance against current implied variance of log returns. More specifically, let S t stand for the closing price of the underlying asset on day t, t = 0, 1, . . . , T , and let R t = log(S t /S t−1 ) represent the tth daily log return. The annualized realized variance over T trading days is defined as where A is the number of trading days per year, which on average is 252. The payoff of a variance swap is defined as where the variance notional N var and variance strike σ 2 strike are specified before the sale of a variance swap contract.
Variance swap provides investors with pure exposure to the variance of the underlying asset without directional risk. It is notably liquid across major equities, indices, and stock markets, and is growing across other markets. Historical evidence indicates selling variance systematically is profitable.
There are numerous methods in the literature of pricing variance swap, including both analytical and numerical methods (see Zhu and Lian (2010) for an extensive review). Nevertheless, a pricing formula or procedure based on a certain stochastic process assumption, such as Lévy process, may not fit the real data well due to the validity of the model assumption (Carr et al., 2012)).
In this section, we propose a moment-based method based on our PC risk-neutral density estimate to price variance swaps, which is free of model assumption.
Assuming the existence of a risk-neutral measure Q, the fair price V S t,T of a variance swap on day t is the discounted expected payoff where R tT is the cumulative risk-free interest rate from t to T defined in Section 2.4. To proceed, we further assume that Assumption 4.1. The increments of the process log S t are independent, that is, log(S t+1 /S t ) is independent of S 0 , . . . , S t , t = 0, . . . , T − 1.
Consequently, the fair price of a variance swap can be represented by a sequence of the risk neutral moments of the underlying asset.
Proposition 4.1. Assuming the existence of a risk-neutral measure Q and Assumption 4.1 is fulfilled, the fair price of variance swap is The proof of Proposition 4.1 is relegated to C.

Moments calculation
In view of Equation (12), we need to estimate the first and second moments of log S i under the risk-neutral measure. One option is a moment-based method described by Bakshi et al. (2003b). In this section, we employ an alternative way of calculating the moments by the proposed nonparametric approach, which provides similar estimates for our purpose.
Recall that the step function f ∆ defined in (1) or (10) provides an approximation to the risk-neutral density f Q of log S T . We use all the available market prices of options to estimate f ∆ , then the moments calculated from f ∆ serve as the estimates of required moments. Since f ∆ is piecewise constant, it can be verified that the first and second moments of log(S T ) are Note that there is no market price available for options that expire on a day that is other than the third Saturday of the delivery month. Therefore, we would have to interpolate the mean and standard deviation of log(S i ) for t < i < T , and this is done via linear interpolation in this paper. Detailed procedure is described in D.

Replicating by variance futures
In order to evaluate the fair price of variance swap that we obtain in Section 4.1, we replicate variance swap using available market prices of variance futures. Variance future is a financial contract that is traded over the counter. As stated by Biscamp and Weithers (2007), variance swap and variance future are essentially the same since they both trade the difference of variance and one can replicate a variance swap by the corresponding variance future. As a matter of fact, if variance future and variance swap share the same expiration date, then at the start point of the observation period, there is no difference between trading a variance future and trading a variance swap with $50 variance notional. The formula for a fair price of variance swap induced from variance future is where M is the number of observed days to date, N e is the expected number of trading days in the observation period, IU G is the square of market implied volatility given by

Results
In order to assess the accuracy of our estimated fair price of variance swap, we compare three relevant quantities: 1. OP: Fair price of variance swap based on our moment-based nonparametric approach, using option market prices till day t; 2. VF: Induced market price of variance swap from CBOE traded variance future till day t; 3. True: Realized price of variance swap at expiration day T with known S 0 , S 1 , . . . , S T .
We present three ratios, OP/True, VF/True, and OP/VF, in Figure 2, against the remaining calendar days of variance swaps. Compared with "True" price based on realized underlying asset prices, Figure 2(a) and Figure 2 (b) show that OP and VF have a similar increasing pattern along with the remaining calendar days. It might be due to the uncertainty in the estimate of the variance, which increases along with the number of days to expiration. On the other hand, Figure 2(c) and Figure 2(d) show that the fair price OP based on our proposed method matches the market price VF pretty well on variance swaps with expiration between 365 days and 800 days. For variance swaps expiring in less than 365 days (not shown here), OP and VF do not match well. A reasonable explanation is that long-term option data are more reasonable and stable, which are less likely to be affected by external factors or noises. For variance swaps longer than 800 days, the relatively low VF might indicate underpriced variance futures.

Discussion
The proposed nonparametric approach is data-driven. It is not built on any model assumption about the underlying data generating process of asset prices. It only assumes the existence of a risk-neutral density and independence of increments of log return for pricing variance swap. That is why it can capture the market price very well.
Compared with other nonparametric methods, such as cubic spline and B-spline, our method is much simpler but fit the real data better. We choose only distinct strikes as knots and assume constant values between knots to avoid the overfitting issue. By sacrificing the continuity of estimated risk-neutral density, the non-negativity requirement for a density function is easy to meet.
On the other hand, the proposed approach utilizes market prices of all options, not just OTM option prices. In our opinion, ITM options, although not so liquid as OTM options, still contain market information and should be incorporated when estimating

VF/True (less than 5) vs Remaining Cdays
Remaining calender days of contract  Figure 2: Comparison OP, VF and True variance swap prices a risk-neutral density. Our comprehensive analysis shows that it recovers OTM option prices better by incorporating ITM option prices. Pricing variance swap is a difficult job when involving real data. We display in Figure 2 only the cases where the ratio OP/True is less than 5. There are cases where OP and VF disagree significantly. Overall, our OP prices work better for variance futures that expire in the last four months of 2015, which are also the last four months available in our data.

C Proof of Proposition 4.1
Since E Q