American Option Pricing Using Particle Filtering Under Stochastic Volatility Correlated Jump Model

2.1 Optimal stopping problem under partial information

Assume that we have a hidden Markov model {(L_n, S_n), n = 0, 1, ⋯, T} as described in (1)∼(3) under a finite dimensional time space 𝓣 : eqq{0, 1, ⋯, T}. The filtration generated by {(L_n, S_n)} is 𝓕_n = σ {(L_i, S_i);i = 0, 1, ⋯, n}, and the observable filtration generated by observable data is FnS = σ (S₀, S₁, ⋯, S_n). A random variable τ : Ω → 𝓣 is a FnS -stopping time if {τ ≤ n} ∈ FnS for every n ∈ 𝓣. Also define 𝓣^S as the set of FnS -stopping times τ ∈ 𝓣. At time 0, the initial price S₀ is known and the latent factor L₀ follows a known distribution π₀ = p(l₀), which is derived from the historical data up to S₀.

Let u₀ = u₀ (S₀, π₀) denotes the price of an American option on S₀. Under certain regular conditions, the fundamental theorem of asset pricing shows that the arbitrage-free price of the option with maturity T is a finite-horizon partially observable optimal stopping problem:

where g : 𝓢 × 𝓛 × 𝓣 → ℝ is the payoff function at time τ. For a strike price K, the payoff of a call option is g = g(S_τ) = (S_τ − K)⁺; and for a put option, g(S_τ) = (K − S_τ)⁺. In addition, 𝔼^ℚ [⋅|⋅] stands for the conditional expectation under the chosen equivalent martingale measure ℚ with respect to ℙ. In this setting, the decision maker has only access to S_n at time n, so that the decision is made only relying on FnS.

According to [12], the above partially observable problem can be transformed to an equivalent fully observable form by introducing a new state “filtering distribution”, denoted as Π_n = p(L_n|s_0:n), which can be estimated by sequential Monte Carlo techniques. Given (S_n, Π_n), the optimal pricing problem is equivalent to

where

Theoretically, we can solve it following the dynamic programming recursion:

where C_n(S_n, Π_n, n) is the continuation value at time t defined as

Then the optimal stopping time is

The above recursion show that (S_n, Π_n) are sufficient statistics that determine the optimal stopping time. However, it is often impossible to solve the problem exactly following (7). The difficulty inside comes first from the fact that the filtering distribution Π_n is infinite dimensional, and that a lack of accurate estimation of continuation value C_n(S_n, Π_n, n). To overcome these two problems, Monte Carlo technique can be used to provide approximation. This is essentially done by firstly using particle filter to approximate the filtering distribution Π_n, and summarize it within a finite-dimensional vector that describes the statistical property of it. Then a least-square Monte Carlo can be applied to estimate C_n(S_n, Π_n, n). The detail algorithm is introduced in the next section.

2.2 Generic particle filter based American options pricer

We combine the dynamic programming algorithm with sequential Monte-Carlo techniques to pricing American options by particle filtering. The key to solve (7) is to provide a finite dimensional summary vector πnm of from m filtered particles at time n. The filtering distribution vary from time to time. For stochastic volatility model, the posterior distribution is usually near to the previous one. In such a case, measures of location and scale would be sufficient. The posterior distribution of volatility of the jump models can be highly skewed and may contains multiple peaks, especially when a jump in volatility and price occurs. In such a case, one may think of taking two modes of the filtered particles.

Assumes that under our algorithm, we have M scenarios indexed by i, with N time steps indexed by j, and m particles in each step indexed by k. For sake of simplicity, we present the SIR version here. Since the filtering step can be separated from the backward decision step, it is straightforward to implement the corresponding APF version.

Assume that under the physical measure ℙ, the asset price S_t is driven by two latent factors, namely, the latent variance and the latent jump process. Specifically, we assume that they are the solutions of the following stochastic differential equations:

Where W=(Wts)t≥0 and W=(Wtv)t≥0 are two Brownian motions under probability measure ℙ, with parameters κ, θ, η and ρ are the same as Heston model. In addition, we assume that the jump sizes are Zts∼N(μs,σs2),Ztv∼exp⁡(μv) with correlation ρ_z. The jump times N_t ∼ Poi(λ_t) is a Poisson process with intensity λ_t for both spot and variance. This implicates the spot and variance would jump simultaneously, which is coherent with listed option data. For sake of simplicity, we assume the independence of jump and spot volatility, i.e. the correlation ρ_z = 0, because this parameter is difficult to estimate, as argued by [3] and [2].

2.3 Particle filtering in stochastic volatility jump model

The only difference between SVCJ and SV model is to simultaneously sample two latent factors and calculate the importance density in a more careful way. In SIR algorithm, the derivation is more straghtforward. The APF algorithm is more challenging because by first resampling, posterior distributions become some normal mixture.

2.3.3 Modification of APF in SVCJ model

APF algorithm can overcome sampling impoverishment problem by first resampling. However, in the resampling step of spot variance in APF, the conditional expectation as the predictive sample, which may again suffers from such a problem. The problem arises in SVCJ model when the jumps in volatility occur. When there is a jump in volatility, the distribution of post-jump volatility changes dramatically comparing to the pre-jump volatility, and there is a probability that a finite-dimensional sample space does not contain the post-jump particles. We consider the following example.

Consider the case that at time n, we have a spot volatility V_n = 4, and a jump occurs with Zn+1v = 10, resulting in a post-jump volatility to be V_n+1 = 14 and observed return y_n+1 = −15. A typical view of such a case would looks like the first panel in Figure 1.

Figure 1

Illustration of predictive particles calculated via conditional expectation in APF

In Figure 1, we have some initial particles {Vn(k)} as histogramed in red. The corresponding probability distribution of these particles is non-central χ²-distribution, which looks similar with normal distribution, as plotted in red dash-dot line. Now a jump occured at a size of Zn+1v = 10, leads to the updated true volatility value to be V_n+1 = 14, plotted in a blue line. The corresponding probability density function conditioned on the jump occured is plotted in blue dash-dot line. By Bayesian formula, we can update the initial particles incorporating the new information y_n+1 = −15, using first the conditional expectation (19) to “esimate” the volatilities, then calculate the importance weights using (20), and resampling according to the weights. The posterior resampled particles {Vn+1i(k)}k=1m becomes then the reweighted version of initial particles.

However, none of the resampled particles {Vn+1i(k)}k=1m contains, or even be near to, the true post-jump value of volatility. Although we give a large weight for the maximum value of initial particle, and obtain some particles with a value around 8.25, we cannot capture any post-jump particle. The interpretation is quite straghtforward: the prediction step using conditional expectation (19) just shifts forward the particle a bit, and the resampling step is just to reweight the probability of the existed particles and thus does not change the sample state space of the particles. As a result, if a jump does occur, neither do the two steps would create reasonable particles for the post-jump volatility.

To overcome this problem, we can think of a scheme to diversify the particles by adding some particles conditioning on the case when we have jumps. Assume that we have a probability of jump λ, we can assign some of the particles to be having some jumps. Precisely, we can assign a number k_λ with the form

kλ=max1,[λ⋅m](24)

where [x] stands for taking the floor of x to make it an integer, and compute the predictive particles by

V^n+1(k)=Vn(k)+κ(θ−Vn(k))+Zn+1v,(k)(25)

for k = 1, ⋯, k_λ randomly chosen, and the jumps Zn+1v,(k) ∼ exp(μ_v). The intuition behind (24) is quite straightforward: we add at least one particle with jump, and add more if we have enough total particles. With this, we can diversify the particles by manually assigning some particles with jumps. In our example here, the result can be shown in panel 2 of figure 1.

In the panel 2 of figure 1, we first sample V^n+1(k) from (25), and compute the importance weights as in (20) and resample. Not surprisingly, we add some diversification on particles, in the sense that they contains the region of “true post-jump volatility”. Notice that the resampled particles do not typically concentrate on some single point, even at the point where the true value stays. This is a well-known effect shown by many Bayesian statistical problems, that the historical data would tend to reject the posterior distribution and make a compromise in between. In the example here, since the probability of jumps is quite low, the model made a compromise in such a way that it increases almost all the numbers of with-jump particles (comparing to one particle in before), but with a highly non-normal posterior density.

This would introduce new challenges in the American option pricing problem, in which we have to summarize particles in each step.

2.3.4 Numerical experiment

We now show some numerical results on the stochastic volatility model with jumps. The parameters are as follows:

The parameters in Table 1 are based on [3], and a larger observation period T = 2000 can capture more jumps in the simulation. We show an comparison of performances for all models by RMSE and MAE. For each model, we implemented the SIR algorithm, the original APF algorithm (APFo) as proposed by [13] and [5], and the modification of APF algorithm (APFm). The main results are shown in Table 2.

Table 1

Parameters of SVCJ model in the numerical example

κ	θ	η	ρ	λ	μS	σS	μV	T
0.02	0.9	0.15	0	0.006	−2.5	4	2	2000

Table 2

Simulation results for SVCJ model

m	SIR		APFo		APFm
	R[1]	M	R	M	R	M
Variances Vn
100	12.45	8.15	13.01	8.39	11.52	7.41
1000	13.15	6.93	10.90	6.85	9.77	6.26
10000	11.12	6.70	10.96	6.90	9.58	6.29
Jumps in price Jn⋅Zns
100	3.54	0.27	2.34	0.26	2.74	0.29
1000	3.40	0.28	2.16	0.25	2.12	0.23
10000	3.19	0.28	2.05	0.22	2.06	0.22
Jumps in volatility Jn⋅Znv
100	2.44	0.23	2.29	0.25	2.00	0.21
1000	3.11	0.28	1.94	0.22	1.94	0.22
10000	1.82	0.20	1.91	0.21	1.88	0.21

In Table 2, in general, APFm outperforms the other two. Comparing the RMSE of filter variance, for example, we can observe that when a particle number of 10000 is used, which we think can eliminate the Monte Carlo approximation error, the RMSE is near to 1. However, we can find that the APF algorithm is not always performing better than SIR, especially when the jump intensity is large, say 0.2. This is because we are using an approximation of Bernoulli distribution rather than Binomial to compute the importance weights of jumps, and thus would fail when the jump occurs twice. But for daily return, the event is so rare that it is almost negligible.

In Figure 2, panel 1∼4 show the daily return, spot volatility, jumps in price and jumps in volatility respectively. Both SIR and the modified APF algorithm are implemented. The price process is simulated from SVCJ model, which makes it more difficult to “disentangle” the jumps from the volatility comparing to SVJ model.

Figure 2

A representative sample path of daily returns, volatility, jumps in prices and jumps in volatility with corresponding filtered estimators. Parameters are taken from Table 2 with m = 10000, n_h = 100 and T = 5000.

Generally, both SIR and the APFm algorithm follow the true volatility trend. The APFm outperforms SIR in general. It would provide some possibly unreasonable jump particles and thus predict small jumps in volatility. Figure 3 shows the evolution of particles over time. The algorithm is run with APFm with m = 10000 particles. At time n = 325, a jump in both price and volatility occurs, and the posterior distribution shifts dramatically. The histogram provides the posterior density at this time. Initially, the particles are centered around 2.8, and the distribution looks like that from stochastic volatility model. When a jump occurred and is detected, the importance weights favors towards the new post-jump value, with a resistency from history. In such a scenario, none of any traditional statistical measure provides an accurate estimation.

Figure 3

Illustrative filtering distributions over time

3.1 Applications with S&P 500 returns

we utilize pure particle filter over SV, SVJ and SVCJ models under historical measure ℙ. The data set we consider are from 01/02/1986 to 01/02/2013. We use data from 01/02/1986 to 12/31/2004 of S&P 500 returns to calibrate the parameter, and implement particle filter with SV, SVJ and SVCJ models for the period during 01/02/1990 to 12/31/2012 since VIX index are available only from 01/02/1990. We calibrate the model using SA algorithm for SV, SVJ and SVCJ model, the estimated parameters are given as follows:

Then we utilize particle filtering to the subsequent observed return, and the result is shown as in Figure 4.

Figure 4

Estimated latent factors for SV, SVJ and SVCJ model for S&P 500 return, a modified APF algorithm with m = 10000 particles is used.

Comparing the spot volatility, SV model provides a biggest filtered value, because the other two will “smooth” it by introducing jumps. For SVJ and SVCJ model, we observe that the volatility in SVCJ model would generally be larger than the SVJ model. This can be interpreted by the fact that the jumps in volatility would generally make the contribution of volatility to the price change larger. Numerically, since the jump of volatility can only be positive, the volatility in SVCJ should also be larger.

For the option pricing problem, we consider the stochastic volatility model for S&P 500 index future returns and index future options, because the algorithm is quite slow and is impractical to estimate any parameters from a large data set. The S&P 500 index future options are American style.

we take the most active S&P 500 future contract defined by Bloomberg with ticker being “SP1 Index”, the time period is from the first trading day of 2005 to the first trading day of 2013 in total 2,000 observations. The data are obtained from Bloomberg®. For the options data, we utilize a scheme to fetch the data. Precisely, we search for those options that are nearest to be at-the-money, and with the nearest time-to-maturity being 30 days. These options are typically most actively traded. For the rates, we use one-month USD Libor rate as the risk-free rate to match the maturities, and Bloomberg best expectation of dividends to approximate the dividend yields.

We calibrate the constant volatility risk premium introduced by option markets. we fix all the parameters for the stochastic volatility model, and calibrate the volatility risk-premium through a constant λ^ℚ, by minimizing the normalized Mean Squared Error (MSE) of the model price and market observed price:

This method takes a long time to be done, as we have 2000 observations, and to price a single American option using limited samples and particles, e.g. M = 100, N = 100, m = 100, takes already 10 seconds in Matlab. After calibration, we obtain a volatility risk premium being λ = −4.2496%. The result is shown in Figure 5.

Figure 5

Calibrated risk-neutral filtered volatility v.s. VIX index

In panel 1 of Figure 5, we plotted: 1. the calibrated filtered volatility with volatility premium of λ = −4.2496%, denoted by SVQ; 2. the pure particle filtered volatility, denoted by SV; 3. the VIX index. After calibration, SVQ is generally closer to VIX index, especially during the financial crisis, though bias still exists especially because we have only one parameter changed. The RMSE of SVQ is 20.7853 compared to 21.6994 of SV.

In panel 2 of Figure 5, we plotted: 1. the Black-Scholes implied volatility of predicted price, computed via calibrated volatility using our algorithm; 2. the Black-Scholes implied volatility of the market prices, which is different from deriving directly from an European option; 3. the VIX index. We also plot the VIX index. To the S&P 500 index future options, it seems that the volatility risk premium has changed over period, because we observe an inconsistency through the time, especially before and after the global financial crisis at 2008, of which the sign being Lehman’s bankruptcy on September 15, 2008. In order to study further the volatility risk premium, we adopt another set of most liquidly traded American options written on Apple Inc.

3.2 Application to American options written on equities

In this section, we apply the particle filter to some other more liquidly traded options. We use daily closing price of Apple Inc. that we obtain also from Bloomberg®. Apple Inc. is one of the largest companies in the world and its options are among in the most liquid. The data we used for estimating the stochastic volatility under historical measure is from 01/03/2012 to 10/30/2012. We choose this period because it goes through a relatively complete cycle.

We also obtain the American put option prices from Bloomberg®. During that period, the most actively traded options are 𝓚 = {550, 555, 560, 565, 570} with maturities 𝓣 = {12/22/2012, 01/19/2013, 02/16/2013}, so we have 15 options at each day. The filtered volatility of Apple Inc. is given in Figure 6.

Figure 6

Apple daily closing prices and filtered volatility under stochastic volatility model

Once the volatility is estimated from the underlying prices, we may start to calibrate the volatility premium from market option prices. We also implemented a gradient-based approach to minimize the normalized Mean Squared Error (MSE) of the model price and market observed price. Since the pricing algorithm is quite time-consuming. In our empirical study, we fit only one day option prices, namely 15 options from 5 different strikes and 3 time to maturities.

Figure 7 shows the calibrated result. We observe that the volatility premium λ^ℚ is −5.882%, which we think is reasonable because the more negative λ^ℚ is, the bigger the price would be since Vega is usually positive, which incorporate our assumptions that traders requires additional premium from the uncertainty from volatility.

Figure 7

Model calibrated option prices and observed prices