Sequential Monte Carlo pricing of American-style options under stochastic volatility models\thanksrefTITL1

# Sequential Monte Carlo pricing of American-style options under stochastic volatility models\thanksrefTitl1

[ [    [ [ U.S. Department of the Treasury and Carnegie Mellon University Office of the Comptroller of the Currency
U.S. Department of the Treasury
250 E Street SW
Washington, DC 20219
USA
Department of Statistics
Carnegie Mellon University
232 Baker Hall
Pittsburgh, Pennsylvania 15213
USA
\smonth1 \syear2008\smonth6 \syear2009
\smonth1 \syear2008\smonth6 \syear2009
\smonth1 \syear2008\smonth6 \syear2009
###### Abstract

We introduce a new method to price American-style options on underlying investments governed by stochastic volatility (SV) models. The method does not require the volatility process to be observed. Instead, it exploits the fact that the optimal decision functions in the corresponding dynamic programming problem can be expressed as functions of conditional distributions of volatility, given observed data. By constructing statistics summarizing information about these conditional distributions, one can obtain high quality approximate solutions. Although the required conditional distributions are in general intractable, they can be arbitrarily precisely approximated using sequential Monte Carlo schemes. The drawback, as with many Monte Carlo schemes, is potentially heavy computational demand. We present two variants of the algorithm, one closely related to the well-known least-squares Monte Carlo algorithm of Longstaff and Schwartz [The Review of Financial Studies 14 (2001) 113–147], and the other solving the same problem using a “brute force” gridding approach. We estimate an illustrative SV model using Markov chain Monte Carlo (MCMC) methods for three equities. We also demonstrate the use of our algorithm by estimating the posterior distribution of the market price of volatility risk for each of the three equities.

\kwd
\doi

10.1214/09-AOAS286 \volume4 \issue1 2010 \firstpage222 \lastpage265 \setattributecopyrightownerIn the Public Domain \newproclaimassumption[theorem]Assumption \newproclaimremarkRemark

\runtitle

Sequential Monte Carlo pricing of American options \thankstextTITL1Supported in part by equipment purchased with NSF SCREMS Grants DMS-05-32083 and DMS-04-22400.

{aug}

A]\fnmsBhojnarine R. \snmRambharat\thanksreft1label=e1]ricky.rambharat@occ.treas.gov\corref and B]\fnmsAnthony E. \snmBrockwelllabel=e2]abrock@stat.cmu.edu \thankstextt1The views expressed in this paper are solely those of the authors and do not represent the opinions of the U.S. Department of the Treasury or the Office of the Comptroller of the Currency. All remaining errors are the authors’ responsibility.

Optimal stopping \kwddynamic programming \kwdarbitrage \kwdrisk-neutral \kwddecision \kwdlatent volatility \kwdvolatility risk premium \kwdgrid \kwdsequential \kwdMonte Carlo \kwdMarkov chain Monte Carlo.

## 1 Introduction

American-style option contracts are traded extensively over several exchanges. These early-exercise financial derivatives are typically written on equity stocks, foreign currency and some indices, and include, among other examples, options on individual equities traded on The American Stock Exchange (AMEX), options on currency traded on the Philadelphia Stock Exchange (PHLX) and the OEX index options on the S&P 100 Index traded on the Chicago Board Options Exchange (CBOE). As with any other kind of option, methods for pricing are based on assumptions about the probabilistic model governing the evolution of the underlying asset. Arguably, stochastic volatility models are the most realistic models to date for underlying equities, but existing methods for pricing American-style options have mostly been developed using less realistic models, or else assuming that volatility is observable. In this paper we develop a new method for pricing American-style options when the underlying process is governed by a stochastic volatility model, and the volatility is not directly observable. The method yields near-optimal solutions under the model assumptions, and can also formally take into account the market price of volatility risk (or volatility risk premium).

It follows from the fundamental theorem of arbitrage that an option price can be determined by computing the discounted expectation of the payoff of the option under a risk-neutral measure, assuming that the exercise decision is made so as to maximize the payoff. While this is simple to compute for European-style options, as illustrated in the celebrated papers of Black and Scholes (1973) and Merton (1973), the pricing problem is enormously more difficult for American-style options, due to the possibility of early exercise. For American-style options, the price is in fact the supremum over a large range of possible stopping times of the discounted expected payoff under a risk-neutral measure. A range of methods has been developed to find this price, or equivalently, to solve a corresponding stochastic dynamic programming problem. Glasserman (2004) provides a thorough review of American-style option pricing with a strong emphasis on Monte Carlo simulation-based procedures. Due to the difficulty of the problem, certain assumptions are usually made. For instance, a number of effective algorithms [including those developed in Brennan and Schwartz (1977), Broadie and Glasserman (1997), Carr, Jarrow, and Myneni (1992), Geske and Johnson (1984), Longstaff and Schwartz (2001), Rogers (2002), and Sullivan (2000)] are based on the assumption that the underlying asset price is governed by a univariate diffusion process with a constant and/or directly observable volatility process.

As recognized by Black and Scholes (1973) and others, the assumption of constant volatility is typically inconsistent with observed data. The volatility “smile” (or “smirk”) is one example where empirical data show evidence against constant volatility models. The smile (smirk) effect arises when option contracts with different strike prices, all other contract features being equivalent, result in different implied volatilities (i.e., the volatility required to calibrate to market observed option prices). A variety of more realistic models has subsequently been developed for asset prices, with stochastic volatility models arguably representing the best models to date. This has led researchers to develop pricing methods for European-style options when the underlying asset price is governed by stochastic volatility models [e.g., Fouque, Papanicolaou, and Sircar (2000), Heston (1993), Hull and White (1987) and Stein and Stein (1991)]. However, work on pricing of American-style options under stochastic volatility models is far less developed. A number of authors [including Clarke and Parrott (1999), Finucane and Tomas (1997), Fouque, Papanicolaou, and Sircar (2000), Guan and Guo (2000), Tzavalis and Wang (2003) and Zhang and Lim (2006)] have made valuable inroads in addressing this problem, but most assume that volatility is observable. Fouque, Papanicolaou, and Sircar (2000) provide an approximation scheme based on the assumption of fast mean-reversion in the volatility process, and they use a clever asymptotic expansion method to correct the constant volatility option price to account for stochastic volatility. The correction involves parameters estimated from the implied volatility surface and they derive a pricing equation that does not depend directly on the volatility process. Tzavalis and Wang (2003) use an analytic-based approach whereby they compute the optimal stopping boundary using Chebyshev polynomials to value American-style options in a stochastic volatility framework. They derive an integral representation of the option price that depends on both the share price and level of volatility.

Additionally, the approach in Clarke and Parrott (1999) uses a multi-grid technique where both the asset price and volatility are state variables in a two-dimensional parabolic partial differential equation (PDE). Pricing options in a stochastic volatility framework using PDE methods are feasible once we assume that volatility itself is an observed state variable. Typically, there is a grid in one dimension for the share price and another dimension for the volatility. The final option price is a function of both the share price and volatility. Guan and Guo (2000) derive a lattice-based solution to pricing American options with stochastic volatility. Following the work in Finucane and Tomas (1997), they construct a lattice that depends directly on the asset price and volatility and they illustrate an empirical study where they back out the parameters from the stochastic volatility model using data on American-style S&P 500 futures options. The valuation algorithm that they develop, however, involves an explicit dependence on both the share price and volatility state variables. Additionally, Zhang and Lim (2006) propose a valuation approach that is based on a decomposition of American option prices, however, volatility is a variable in the pricing result.

The approach we develop, in contrast with the aforementioned methods, combines the optimal decision-making problem with the volatility estimation problem. We assume that the asset price follows a stochastic volatility model, that observations are made at discrete points in time and that exercise decisions are made immediately after each observation. It could be argued that volatility should be considered observable since one could simply compute “Black and Scholes type” implied volatilities from observed option prices. However, implied volatilities are based on the assumption of a simple geometric Brownian motion model (or some other simplified diffusion process) and, thus, their use would defeat the purpose of developing pricing methods with more realistic models. The implied volatility calculation is not as straightforward in a stochastic volatility (multivariate) modeling framework. Renault and Touzi (1996) approximate a Black and Scholes type implied volatility quantity in a Hull and White (1987) setting. The analysis in Renault and Touzi (1996) computes filtered volatilities using an iterative approach and illustrates applications to hedging problems. On the other hand, we aim to compute the posterior distribution of volatility conditional on observed data. Our pricing scheme is based on two key observations. First, we use a sequential Monte Carlo (also referred to as “particle filtering”) scheme to perform inference on the unobserved volatility process at any given point in time. Second, conditional distributions of the unobserved volatility at a given point in time, given current and past observations of the price process, which are necessary for finding an exact solution to the dynamic programming problem, can be well approximated by a summary vector or by a low-dimensional parametric family of distributions.

Inference on the latent volatility process for the purpose of option pricing is an application of the more general methodology that addresses partially observed time series models in a dynamic programming/optimal control setting [see, e.g., Bertsekas (2005, 2007) and the references therein]. Among earlier work, Sorenson and Stubberud (1968) provide a method based on Edgeworth expansions to estimate the posterior density of the latent process in a nonlinear, non-Gaussian state-space modeling framework. Our objective in this paper is to illustrate an algorithm that allows an agent (a holder of an American option) to optimally decide the exercise time assuming that both the share price and volatility are stochastic variables. In our partially observed setting, only the share price is observable; volatility is a latent process. The main challenge for the case of American-style options is determining the continuation value so that an optimal exercise/hold decision can be made at each time point.

Several researchers have applied Monte Carlo methods to solve the American option pricing problem. Most notable among these include Longstaff and Schwartz (2001), Tsitsiklis and Van Roy (2001), Broadie and Glasserman (1997) and Carrière (1996). An excellent summary of the work done on Monte Carlo methods and American option pricing is presented in Chapter 8 of Glasserman (2004). The least-squares Monte Carlo (LSM) algorithm of Longstaff and Schwartz (2001) has achieved much popularity because of its intuitive regression-based approach to American option pricing. The LSM algorithm is very efficient to price American-style options since as long as one could simulate observations from the pricing model, then a regression-based procedure could be employed along with the dynamic programming algorithm to price the options. Pricing American-style options in a stochastic volatility framework is straightforward using the methodology of Longstaff and Schwartz (2001) as long as draws from the share price and volatility processes can be obtained. That is, at each time point, , the decision of whether or not to exercise an American option will be a function of the time share price, (and possibly some part of its recent history , ), and the time volatility, . Our pricing framework, however, needs to accommodate a latent volatility process.

We propose to combine the Longstaff and Schwartz (2001) idea with a sequential Monte Carlo step whereby at each time point , we estimate the conditional (posterior) distribution of the latent volatility given the observed share price data up to that time. We propose a Monte Carlo based approach that uses a summary vector to capture the key features of this conditional distribution. As an alternative, we also explore a grid-based approach, studied in Rambharat (2005), whereby we propose a parametric approximation to the conditional distribution that is characterized by the summary vector. Therefore, our exercise decision at a given time point is based on the observed share price and the computed summary vector components at that time. Although our pricing approach is computationally intensive, as it combines nonlinear filtering with the early-exercise feature of American-style option valuation, it provides a way to solve the associated optimal stopping problem in the presence of a latent stochastic process. We compare our approach to a basic LSM method whereby at a given time point, we use a few past observations to make the exercise decision in order to price American-style options in a stochastic volatility framework. We also compare our method to the LSM approach using the current share price and an estimate of the current realized volatility as a proxy for the true volatility. In order to assess precision, we compare LSM with past share prices, LSM with realized volatility, and our proposed valuation technique to the LSM-based American option price assuming that share price and volatility can be observed (the full information state). The method closest to the full information state benchmark would be deemed the most accurate.

The present analysis addresses the problem of pricing an American-style option, once a good model has already been found. It is worth noting, however, that since neither the sequential Monte Carlo scheme nor the gridding approach used in our pricing technique are tied to a particular model, the method is generalizable in a straightforward manner to handle a fairly wide range of stochastic volatility models. Thus, it could be used to perform option pricing under a range of variants of stochastic volatility models, such as those discussed in Chernov et al. (2003). Although the focus of our paper is not model estimation/selection methodology, we do implement a thorough statistical exercise using share price history to estimate model parameters from a stochastic volatility model. [See Chernov and Ghysels (2000), Eraker (2004), Gallant, Hsieh, and Tauchen (1997), Ghysels, Harvey, and Renault (1996), Jacquier, Polson, and Rossi (1994), Kim, Shephard, and Chib (1998) and Pan (2002) for examples of work that mainly address model estimation/selection issues]. Our approach will be to employ a sequential Monte Carlo based procedure to estimate the log-likelihood of our model [see, e.g., Kitagawa and Sato (2001) and the references therein] and then use a Markov chain Monte Carlo (MCMC) sampler to obtain posterior distributions of all model parameters. Conditional on a posterior summary measure of the model parameters (such as the mean or median), we then estimate the (approximate) posterior distribution of the market price of volatility risk. Our analysis is illustrated in the context of three equities (Dell Inc., The Walt Disney Company and Xerox Corporation).

The paper is organized as follows. In Section 2 we formally state the class of stochastic volatility models with which we work. In Section 3 we review the dynamic programming approach for pricing American-style options and demonstrate how it can be transformed into an equivalent form, and introduce (i) a sequential Monte Carlo scheme that yields certain conditional distributions, and (ii) a gridding algorithm that makes use of the sequential Monte Carlo scheme to compute option prices. Section 4 describes the pricing algorithms and presents some illustrative numerical experiments. In Section 5 we describe (i) the MCMC estimation procedure for the stochastic volatility model parameters, and (ii) the inferential analysis of the market price of volatility risk. Section 6 contains posterior results of our empirical analysis with observed market data. Section 7 provides concluding remarks. Finally, we present additional technical details in the Appendix. We also state references to our computing code and data sets in supplementary.

## 2 The stochastic volatility model

Let be a probability space, and let be a stochastic process defined on , describing the evolution of our asset price over time. Time will be referred to as the “current” time, and we will be interested in an option with expiry time , with being some positive integer, and some positive real-valued constant. Assume that we observe the process only at the discrete time points , and that exercise decisions are made immediately after each observation. (We would typically take the time unit here to be one year, and to be , representing one trading day. However, both and the time units can be chosen arbitrarily subject to the constraints mentioned above.)

Assume that, under a risk-neutral measure, the asset price evolves according to the Itô stochastic differential equations (SDEs)333Under the statistical (or real-world) measure, the asset price evolves on another probability space. Under the real-world measure, the drift term in equation (1) is replaced by the physical drift and the term does not appear in the drift of equation (3). The change of measure between real-world and risk-neutral is formalized through a Radon–Nikodym derivative.

 dS(t) = rS(t)dt+σ(Y(t))S(t)[√1−ρ2dW1(t)+ρdW2(t)], (1) σ(Y(t)) = exp(Y(t)), (2) dY(t) = [α(β−λγα−Y(t))]dt+γdW2(t), (3)

where represents the risk-free interest rate (measured in appropriate time units), is referred to as the “volatility,” measures the co-dependence between the share price and volatility processes, (volatility mean reversion rate), (volatility mean reversion level), and (volatility of volatility) are constants with , , and are assumed to be two independent standard Brownian motions, and is a constant referred to as the “market price of volatility risk” or “volatility risk premium” [Bakshi and Kapadia (2003a), Melenberg and Werker (2001) and Musiela and Rutkowski (1998)]. If we set

 dW∗1(t)=[√1−ρ2dW1(t)+ρdW2(t)],

we can more clearly see that is the correlation between the Brownian motions and . The parameter quantifies the so-called “leverage effect” between share prices and their volatility.

Observe that is not uniquely determined in the above system of SDEs. Since we are working in a stochastic volatility modeling framework, markets are said to be “incomplete” because volatility is not a traded asset and cannot be perfectly hedged. This is to be contrasted with the constant volatility Black and Scholes (1973) framework where a unique pricing measure exists and all risks can be perfectly hedged away. There is a range of possible risk-neutral measures when pricing under stochastic volatility models, each one having a different value of . In fact, there are also risk-neutral measures under which varies over time. However, for the sake of simplicity, we will assume that is a constant. Later in this paper, we illustrate how to estimate .

Equations (1)–(3) represent a stochastic volatility model that accommodates mean-reversion in volatility and we will use it to illustrate our American option valuation methodology. It is an example of a nonlinear, non-Gaussian state-space model. Scott (1987) represents one of the earlier analyses of this stochastic volatility model. This same type of model has also been studied in a Bayesian context by Jacquier, Polson, and Rossi (1994) and, more recently, in Jacquier, Polson, and Rossi (2004) and Yu (2005). It should be noted that our methodology is not constrained to a specific stochastic volatility model. The core elements of our approach would apply over a broad spectrum of stochastic volatility models such as, for instance, the Hull and White (1987) and Heston (1993) stochastic volatility models.

Since our observations occur at discrete time-points we will make extensive use of the discrete-time approximation to the solution of the risk-neutral stochastic differential equations (1)–(3) given by

 St+1 = St⋅exp{(r−σ2t+12)Δ+σt+1√Δ[√1−ρ2Z1,t+1+ρZ2,t+1]}, (4) σt+1 = exp(Yt+1), (5) Yt+1 = (6)

where , , is an independent and identically distributed (i.i.d.) sequence of random variables with standard normal [] distributions,

 β∗=β−λγα,

and all other parameters are as previously defined. Thus, and represent approximations, respectively, to and . [The expression for is obtained directly from the exact solution to (3), while the expression for is the solution to (1) that one would obtain by regarding to be constant on successive intervals of length . The approximation for becomes more accurate as becomes smaller.]

It will sometimes be convenient to express (4) in terms of the log-returns , as

 Rt+1=(r−σ2t+12)Δ+σt+1√Δ(√1−ρ2Z1,t+1+ρZ2,t+1). (7)

To complete the specification of the model, we can assign a normal distribution,

 Y0∼N(β∗,γ22α). (8)

This is simply the stationary (limiting) distribution of the first-order autoregressive process . However, for practical purposes, it will usually be preferable to replace this distribution by the conditional distribution of , given some historical observed price data Another reasonable starting value for is a historical volatility based measure (i.e., the log of the standard deviation of a few past observations). Additionally, there are examples of stochastic volatility models where exact simulation is not feasible. In such cases, we must resort to numerical approximation schemes such as the Euler–Maruyama method (or any other related approach). Kloeden and Platen (2000) illustrate several numerical approximation schemes that could be applied to simulate from a stochastic volatility model where no exact simulation methodology exists.

## 3 Dynamic programming and option pricing

The arbitrage-free price of an American-style option is

 supτ∈TERN[exp(−rτ)g(Sτ)], (9)

where is a random stopping time at which an exercise decision is made, is the set of all possible stopping times with respect to the filtration defined by

 Ft=σ(S0,…,St),

represents the expectation, under a risk-neutral probability measure, of its argument, and denotes the payoff from exercise of the option when the underlying asset price is equal to . For example, a call option with strike price has payoff function , and a put option with strike price has payoff function . [The analysis in this paper is in the context of American put options since these options typically serve as canonical examples of early-exercise derivatives; see Karatzas and Shreve (1991, 1998) and Myneni (1992) for key mathematical results concerning American put options.] Since is a stopping time, the event must be -measurable, or equivalently, the decision to exercise or hold at a given time must be made only based on observations of the previous and current values of the underlying price process. To allow for the possibility that the option is never exercised, we adopt the convention that if the option is not exercised at or before expiry, along with the convention that In order to price the American option, we need to find the stopping time at which the supremum in (9) is achieved. While it is not immediately obvious how one might search through the space of all possible stopping times, this problem is equivalent to a stochastic control problem, which can be solved (in theory) using the dynamic programming algorithm, which was originally developed by Bellman (1953). Ross (1983) also describes some key principles of stochastic dynamic programming and a thorough treatment is presented, in the context of financial analysis, by Glasserman (2004).

Our objective is to find the optimal stopping rule (equivalently, the optimal exercise time) while taking into account the latent stochastic volatility process. The key difficulty arises because the price itself is not Markovian. We would like to get around this by using the fact that the bivariate process, composed of both price and volatility, is Markovian, but unfortunately we can only observe one component of that process. It is known that one can still find the optimal stopping rule if one can determine (exactly) the conditional distribution of the unobservable component of the Markov process, given current and past observable components [see, e.g., DeGroot (1970) and Ross (1983)]. In such a case, we can use algorithms like the ones described in Brockwell and Kadane (2003) to find the optimal decision rules. These algorithms effectively integrate the utility function at each point in time over the distribution of unknown quantities given observed quantities. Unfortunately, in the context of this paper, we cannot obtain the required conditional distributions exactly, but we can find close approximations to them. We will therefore approach our pricing problem by using numerical algorithms in the style of Brockwell and Kadane (2003), in conjunction with close approximations to the required distributions. In doing so, we make the assumption (stated later in this paper) that our distributional approximations are close to the required conditional distributions.

### 3.1 General method

The dynamic programming algorithm constructs the exact optimal decision functions recursively, working its way from the terminal decision point (at time ) back to the first possible decision point (at time ). In addition, the procedure yields the expectation in the expression (9), which is our desired option price. The algorithm works as follows.

Let denote the decision made immediately after observation of , either to exercise () or hold () the option. While either decision could be made, only one is optimal, given available information up to time . (In the event that both are optimal, we will assume that the agent will exercise the option.) We denote the optimal decision, as a function of the available observations, by

 d∗t(s0,…,st)∈{E,H}.

Here and in the remainder of the paper, we adopt the usual convention of using (upper case) to denote the random variable representing the equity price at time , and (lower case) to denote a particular possible realization of the random variable. Next, let

and for , let denote the discounted expected payoff of the option at time , assuming that decision is made, and also assuming that optimal decisions are made at times It is obvious that, at the expiration time,

 d∗T(s0,…,sT)=argmaxdT∈{E,H}uT(s0,…,sT,dT).

The optimal decision functions can then be obtained by defining

 u∗t(s0,…,st)=ut(s0,…,st,d∗t(s0,…,st)),t=0,…,T, (11)

and using the recursions

 ut(s0,…,st,dt) (12)
 d∗t(s0,…,st)=argmaxdt∈{E,H}ut(s0,…,st,dt). (13)

These recursions are used sequentially, for and yield the (exact) optimal decision functions (Each is optimal in the space of all possible functions of historical data ) The corresponding stopping time is simply

 τ=min({t∈{0,…,T}|dt=E}∪{∞}).

Furthermore, the procedure also gives the risk-neutral option price, since

 u∗0(s0)=supτ∈TERN[exp(−rτ)g(Sτ)].

In practice, it is generally not possible to compute the optimal decision functions, since each needs to be computed and stored for all (infinitely many) possible combinations of values of its arguments However, in what follows we will develop an approach which gives high-quality approximations to the exact solution. The approach relies on exploiting some key features of the American-style option pricing problem.

### 3.2 Equivalent formulation of the dynamic programming problem

First, it follows from the Markov property of the bivariate process that we can transform the arguments of the decision functions so that they do not increase in number as increases. Let us define

 πt(yt)dt=P(Yt∈dyt|S0=s0,…,St=st),t=0,…,T, (14)

where denotes the observed value of , so that denotes the conditional density of the distribution of (with respect to the Lebesgue measure), given historical information . Then we have the following result.

###### Lemma 3.1

For each , can be expressed as a functional of only , and , that is,

 ut(s0,…,st,dt)=~ut(st,πt,dt).

Consequently, for each , and can also be expressed, respectively, as functionals

 d∗t(s0,…,st) = ~d∗t(st,πt), u∗t(s0,…,st) = ~u∗t(st,πt).

This is a special case of the well-known result [see, e.g., Bertsekas (2005, 2007)] on the sufficiency of filtering distributions in optimal control problems. A proof is given in the Appendix. It is important to note here that the argument to the functions , and is a function itself.

Lemma 3.1 states that each optimal decision function can be expressed as a functional depending only on the current price and the conditional distribution of (the latent process that drives volatility) given observations of prices In other words, we can write the exact equivalent form of the dynamic programming recursions,

where and

### 3.3 Summary vectors and sequential Monte Carlo

In order to implement the ideal dynamic programming algorithm as laid out in Section 3.2, we would need (among other things) to be able to determine the filtering distributions Unfortunately, since these distributions are themselves infinite-dimensional quantities in our modeling framework, we cannot work directly with them. We can, however, recast the dynamic programming problem in terms of -dimensional summary vectors

 Qt=⎡⎢ ⎢⎣f1(πt)⋮fl(πt)⎤⎥ ⎥⎦, (17)

where are some functionals. The algorithms we introduce can be used with any choice of these functionals, but it is important that they capture key features of the distribution. As typical choices, one can use the moments of the distribution. In the examples in this paper, we take , and Adding components to this vector typically provides more comprehensive summaries of the required distributions, but incurs a computational cost in the algorithms.

Ultimately, we need to be able to compute and thus . One way to do this is with the use of sequential Monte Carlo simulation [also known as “particle filtering,” see, e.g., Doucet, de Freitas, and Gordon (2001), for discussion, analysis and examples of these algorithms]. The approach yields samples from (arbitrarily close approximations to) the distributions , and thus allows us to evaluate components of . A full treatment of sequential Monte Carlo methods is beyond the scope of this paper, but the most basic form of the method, in this context, appears in Algorithm 1.

This algorithm yields collections of particles, with the property that for each , can be regarded as an approximate sample of size from the distribution . The algorithm has the convenient property that as increases, the empirical distributions of the particle collections converge to the desired distributions . Typically is chosen to be as large as possible subject to computational constraints; for the algorithms in this paper, we choose to be around . There are several ways to improve the sampling efficiency of Algorithm 1, namely, ensuring that the weights in equation (18) do not lead to degeneration of the particles. [For more details and various refinements of this algorithm, refer to, among others, Kitagawa and Sato (2001), Liu and West (2001) and Pitt and Shephard (1999).]

### 3.4 Approximate dynamic programming solution

In order to motivate our proposed American option pricing methodology, we state an assumption that describes how relates to past price history captured in .

{assumption}

The summary vector is “close to sufficient,” that is, it captures enough information from the past share price history, so that is close to .

(If the summary vector was sufficient, then the dynamic programming algorithm would yield exact optimal decision rules. Of course, even in this ideal case, the numerical implementation of a dynamic programming algorithm introduces some small errors.)

It is beyond the scope of this paper to quantify the “closeness” between and . One could, however, theoretically do so using standard distribution distance measures and perform an analysis of the propagation of the error through the algorithms discussed in this paper.

Combining the equivalent form of the dynamic programming recursions (15), (16) along with Assumption 3.4, we can approximate the dynamic programming recursions by

with

 ^d∗t(st,Qt)=argmaxdt∈{E,H}^ut(st,Qt,dt) (21)

and

 ^u∗t(st,Qt)=^ut(st,Qt,^d∗t(st,Qt)), (22)

where is the vector summarizing [cf. equation (17)]. The recursion (3.4) is convenient because the required expectation can be approximated using the core of the sequential Monte Carlo update algorithm.

Assumption 3.4 motivates a practical, approximate solution to the ideal formulation of the dynamic programming problem described in Section 3.2. In the special case of linear Gaussian state-space models, the vector would form a sufficient statistic since would be Gaussian, and our approach would reduce to the standard Kalman filtering procedure [see Chapter 12 of Brockwell and Davis (1991) for details]. For the case of nonlinear, non-Gaussian state-space models, such as the illustrative one in this paper, is not summarized by a finite-dimensional sufficient statistic. Assumption 3.4 permits an approximate solution to our American option pricing problem and, in particular, motivates the conditional expectation in equation (3.4) by using to summarize information about the latent volatility using the past price history. We next introduce two algorithms for pricing American-style options, making use of the summary vectors described in Section 3.3, and we illustrate how our algorithms perform through a series of numerical experiments.

## 4 Pricing algorithms

### 4.1 A least-squares Monte Carlo based approach

The popular least-squares Monte Carlo (LSM) algorithm of Longstaff and Schwartz (2001) relies essentially on approximating the conditional expectation in (3.1) by a regression function, in which one or several recent values are used as explanatory variables. Since the conditional expectation is in fact (exactly) a functional of the filtering distribution , we might expect to obtain some improvement by performing the regression using summary features of as covariates instead. This variant of the LSM algorithm, describing the simulation component of our pricing methodology, is stated in Algorithm 2.

{remark*}

For the stochastic volatility model used in this analysis, measures of center and spread will suffice to capture the key features of the distribution. Therefore, our summary vector in Algorithm 2 describes the mean and standard deviation of the filtering distribution.

{remark*}

The summary vector in Algorithm 2 can include as many key measures of the filtering distribution as needed to accurately describe it. Other types of stochastic volatility models may require additional measures that capture potential skewness, kurtosis or modalities in the filtering distribution. One can learn of the need for such measures by doing an empirical analysis on historical data using Algorithm 1 to gain insights into the behavior of the filtering distribution .

Next, we illustrate in Algorithm 3 the implementation of the LSM regression step using our summary vector of .

{remark*}

The regression in Step 3 of Algorithm 3 uses basis functions of the share price and the summary vector to form the explanatory variables. We choose Laguerre functions as basis functions. Our summary vector consists of the mean and standard deviation of the filtering distribution . The design matrix used in the regression at time point consists of the first two Laguerre functions in , and , and a few cross-terms of these covariates. Specifically, our covariates used in the regression at time (in addition to the intercept term) are

 L0(Sk),L1(Sk),L0(μk),L1(μk),L0(ζk),L1(ζk), L0(Sk)∗L0(μk),L0(Sk)∗L0(ζk),L1(Sk)∗L1(μk), L1(Sk)∗L1(ζk),L0(μk)∗L0(ζk),L1(μk)∗L1(ζk),

where and and, in general, . Other choices of basis functions, such as Hermite polynomials or Chebyshev polynomials, are also reasonable alternatives that could be used to implement the least-squares Monte Carlo algorithm of Longstaff and Schwartz (2001).

{remark*}

Longstaff and Schwartz (2001) actually adjust Step 4 of Algorithm 3 as follows:

in order to avoid computing American option prices with a slight upward bias due to Jensen’s inequality; we follow their recommendation and use this adjustment. They also suggest using only paths where (i.e., the “in-the-money” paths) as a numerical improvement. However, we could use all paths since the convergence of the algorithm also holds in this case [see Clement, Lamberton, and Protter (2002)]. Therefore, we use all paths in our implementation of Algorithm 3, as this produces similar results.

### 4.2 A grid-based approach

We now present a grid-based algorithm for determining approximate solutions to the dynamic programming problem. This algorithm is based on a portion of the research work in Rambharat (2005). The approach uses the vectors summarizing the filtering distributions as arguments to the decision and value functions. In contrast to the Monte Carlo based approach of Section 4.1 where we directly use a summary vector in the LSM method, the grid-based technique requires a distribution to approximate . This distribution will typically be parameterized by the summary vector when used to execute the grid-based algorithm. In our illustrative pricing model, we use a Gaussian distribution to approximate the filtering distribution, , at each time point . Kotecha and Djuric (2003) provide some motivation for using Gaussian particle filters, namely, Gaussian approximations to the filtering distributions in nonlinear, non-Gaussian state-space models. However, any distribution that approximates reasonably well could be used for our purposes. Such a choice will depend on the model and an empirical assessment of .

Since the steps of the sequential Monte Carlo algorithm are designed to make the transition from a specified to the corresponding distribution , we can combine a standard Monte Carlo simulation approach with the use of Steps 1, 2 and 3 of Algorithm 1 to compute the expectation. To be more specific, the next algorithm computes the expectation on the right-hand side of equation (3.4).

Algorithm 4 works by drawing pairs from the conditional distribution of , given , and using these to compute a Monte Carlo estimator of the required conditional expectation. Hence, we are able to evaluate at various points, given knowledge of , and will thus form a key component of the backward induction step. Note that this algorithm also relies on Assumption 3.4 in its use of the summary vector .

Since we will be interested in storing the functions and , we next introduce some additional notation. Let

 G={gi∈Rdq+1,i=1,2,…,G} (24)

denote a collection of grid points in , where denotes the dimensionality of These are points at which we will evaluate and store the functions and . We will typically take

 G=G1×G2×⋯×Gdq+1, (25)

where is a grid of possible values for the share price, and , , is a grid of possible values for the st component of .

We state our grid-based pricing routine in Algorithm 5. [This is a standard gridding approach to solving the dynamic programming equations, as described, e.g., in Brockwell and Kadane (2003).]

{remark*}

As is the case for the Monte Carlo based approach that we describe in this paper, the grid-based scheme stated in Algorithm 5 also gives an option price which assumes that no information about volatility is available at time . In the absence of such information, we just assume that the initial log-volatility can be modeled as coming from the limiting distribution of the autoregressive process and take as the summary measure of this distribution (e.g., its mean and variance if the limiting distribution is Gaussian). However, in most cases, it is possible to estimate log-volatility at time using previous observations of the price process . In such cases, would represent the mean and conditional variance of , given “previous” observations for some . These could be obtained in a straightforward manner by making use of the sequential Monte Carlo estimation procedure described in Algorithm 1 (appropriately modified so that time becomes time ). We could also base on a historical volatility measure (i.e., the standard deviation of a few past observations).

{remark*}

As with any quadrature-type approach, grid ranges must be chosen with some care. In order to preserve quality of approximations to the required expectations in Algorithm 5, it is necessary for the ranges of the marginal grids to contain observed values of the respective quantities with probability close to one. Once the stochastic volatility model has been fit, it is typically relatively easy to determine such ranges.

{remark*}

The evaluation of in the previous algorithm is performed making use of the Monte Carlo approximation given by (23). Obviously the expression relies on knowledge of , but since we have only evaluated at grid points , it is necessary to interpolate in some manner. Strictly speaking, one could simply choose the nearest grid point, and rely on sufficient grid density to control error. However, inspection of the surface suggests that local linear approximations are more appropriate. Therefore, in our implementations, we use linear interpolation between grid points.

### 4.3 Numerical experiments

The pricing algorithms outlined in Sections 4.1 and 4.2 demonstrate how to price American-style options in a latent stochastic volatility framework. These pricing algorithms are computationally intensive, however, their value will depend on how accurately they price American options in this partial observation setting. We next illustrate the applicability of our pricing algorithms through a series of numerical experiments. We assess the accuracy of our valuation procedure by pricing American-style put options using the following methods. (All methods use the current share price in the LSM regression, however, the difference in each method is how volatility is measured.)

• Method A (basic LSM). This method simulates asset prices according to the model (4)–(6), however, the regression step in the LSM algorithm uses a few past observations () as a measure of volatility in lieu of the summary vector . This procedure is most similar to the traditional LSM approach of Longstaff and Schwartz (2001).

• Method B (realized volatility). This procedure simulates asset prices according to the model (4)–(6), however, for each time point () and for each path (), we compute a measure of realized volatility

 RVk,l=1kk∑j=1R2j,l, (27)

where is the return at time along path . The LSM regression step then proceeds to use the realized volatility measure at each time point as a measure of volatility.

• Method C (MC/Grid). This method uses the algorithms described in Sections 4.1 and 4.2 to price American-style options in a latent stochastic volatility framework either using (i) a pure simulation-based Monte Carlo (MC) approach or (ii) a Grid-based approach. Along with the simulated share prices , this approach makes extensive use of the summary vectors that capture key features of the volatility filtering distribution .

• Method D (observable volatility). This approach simulates the asset prices and the volatility , however, it assumes that both asset price and volatility are observable. This is the full observation case that we will use as a benchmark. Whichever of methods A, B or C is closest to method D will be deemed the most accurate.

Figure 1 presents an illustrative example of the difference in American put option prices between method A and method D for several types of option contracts. One could think of this illustration as reporting the difference in pricing results for two extremes: the minimum observation case (method A or basic LSM) and the full observation case (method D or observable volatility). This figure uses parameter settings where stochastic volatility is prevalent. The differences in option prices indeed show that stochastic volatility matters when computing American option prices, especially when volatility of volatility is high [i.e., when in equation (3) is large].

In order to demonstrate the value of our proposed approach, method C (MC/Grid), we must illustrate that it produces more accurate American option pricing results than either of the simpler (and faster) methods (A and B). We experimented with various model and option parameter settings and found situations where simpler methods work just as well as our approach. However, we also found model/option parameter settings where our approach outperforms the other pricing methods. Table 1 describes the settings of our numerical experiments. This table outlines selected values of the model parameters and the American option pricing inputs for use in the pricing algorithms ( is the strike price, is the expiration in days, is the interest rate, is the initial share price, and is a fixed initial volatility).

Table 2 reports the American option pricing results (and standard errors) for methods A (basic LSM) and B (realized volatility) and Table 3 reports the pricing results (and standard errors) for methods C (MC/Grid) and D (observable volatility). The computation times for all methods are also reported. For methods A, B and D, we used LSM paths in the pricing exercise. When running method C, we observed negligible differences between the MC-based and grid-based approaches. Hence, we report the pricing results (and standard errors) for method C using Algorithms 2 and 3 with LSM paths. (The MC-based approach also permits comparisons to the other approaches in terms of standard errors.)

There are some notable observations to be made from the results of the numerical experiments in Tables 2 and 3. First, when the effect of volatility is weak/moderate and the effect of mean reversion is moderate/strong (experiments 1 and 2), all methods result in similar American option pricing results. There are cases in the literature that discuss fast mean reversion [see, e.g., Fouque, Papanicolaou, and Sircar (2000)], and in these cases it would certainly make sense to use a faster pricing method. However, in the situations where we experiment with dominant volatility effects (experiments 3 through 9), although not often encountered in the empirical stochastic volatility literature, but pertinent to volatile market scenarios, method C comes closest to method D. One should note that in cases of dominant volatility method B does better than method A, as it uses a more accurate measure of volatility. Method C, however, which actually makes use of the filtering distributions , comes within standard error of method D (observable volatility case). Additionally, we observed when stochastic volatility is dominant and the American option has a long maturity and is at/in-the-money, the difference in pricing results is more pronounced.

Method C (MC/Grid) is a computationally intensive approach, although this feature is shared by many Monte Carlo and grid-based techniques. Clearly, the proposed pricing algorithm using method C is not competitive in terms of computation time. The simpler methods (A and B) are much faster and also accurate under strong mean reversion and weak stochastic volatility. It should be observed, however, that the pricing accuracy is much higher using method C, as it is always within standard error of method D (observable volatility). The accuracy of method C holds in all model/option parameter settings. One does not require special cases (high mean reversion, low volatility) or special option contract features (long/short maturity, in/at/out-of-the money) for our approach to be competitive in terms of accuracy. We also ascertain the robustness of our approach by computing (and storing) the exercise rule from each of the four methods (A through D). We use the rule to revalue the American put options on an independent, common set of paths. The results are described in Tables 4 and 5; note that the respective pricing results are similar to those stated in Tables 2 and 3, hence leading to similar conclusions.

The differences in pricing results, as noted in Tables 2 and 3 (or Tables 4 and 5), are relevant when thinking about how some option transactions take place in practice. For example, on the American Stock Exchange (AMEX), American-style options have a minimum trade size of one contract, with each contract representing 100 shares of an underlying equity.444Source: http://www.amex.com. Hence, the discrepancies between methods A and B and our proposed method C could be magnified under such trading scenarios. Our proposed approach would be especially useful when constructing risk management strategies during volatile periods in the market.

We also repeated the numerical experiments in Table 1 using a first-order Euler discretization of the model as opposed to exact simulation. The Euler-based results are reported in Tables 12 and 13 of Section .4. Upon inspection, one observes that the corresponding results for each of methods A through D are virtually identical. Thus, if a model does not permit exact simulation, a numerical procedure such as the first-order Euler scheme or any other scheme [see, e.g., Kloeden and Platen (2000)] should suffice for the purposes of executing our pricing algorithm.

We now illustrate a statistical application of our proposed pricing methodology. We first demonstrate how to estimate model parameters. Next, we make inference, using observed American put option prices, on the market price of volatility risk. Since method C outperforms either A or B in all model/option settings, we will use it as our primary tool for statistical analysis.

## 5 Statistical estimation methodology

### 5.1 Model parameter estimation

The stochastic volatility model, outlined in equations (1), (2) and (3), has been analyzed extensively in previous work, namely, Jacquier, Polson, and Rossi (1994, 2004) and Yu (2005). These authors analyze the model under the statistical (or “real-world”) measure as described in a Section 2 footnote. Model estimation of share prices under the real-world measure would only require data on price history. Estimation of model parameters in a risk neutral setting, however, is a bit more involved as we need both share price data on the equity as well as option data. Pan (2002) illustrates how to use both share and option price data to jointly estimate parameters under the real-world measure and risk-neutral measure. Eraker (2004) also implements joint estimation methodology for share and option price data. However, both Pan (2002) and Eraker (2004) only deal with the case of European-style options.

We propose a two-step procedure to estimate our illustrative stochastic volatility model in an American-style derivative pricing framework. In the first step, we use share price data to estimate the statistical model parameters. Define the parameter vector555Strictly speaking, estimation using only share prices (i.e., under the physical measure) involves the physical drift rate in the parameter vector. However, since we do not use the physical drift rate in risk-neutral pricing, we do not present its estimation results.

 θ=(ρ,α,β,γ).

Conditional on estimated model parameter values, we then estimate the market price of volatility risk parameter . Estimation of requires data on both share and American option prices. Although it would be comprehensive to do a joint statistical analysis of both share and option prices, this problem is quite formidable in the American option valuation setting. (The full joint estimation problem is left for future analysis.) We adopt a Bayesian approach to parameter estimation. The first objective is to estimate the posterior distribution of ,

 p(θ|r1,…,rn)=p(r1,…,rn|θ)⋅p(θ)∫p(r1,…,rn|θ)⋅p(θ)dθ. (28)

Since , , , and , we reparameterize a sub-component of the vector so that each component will have its domain in . This reparameterization facilitates exploration of the parameter space. Let us define

 ~ρ=tan(ρπ2),~α=log(α),and~γ=log(γ).

We assign independent standard normal priors to the components of the (reparameterized) vector

 ~θ=[~ρ,~α,β,~γ].

The next step is the evaluation of the likelihood (or log-likelihood) for , however, an analytical expression for the likelihood is not available in closed-form. We employ Kitagawa’s algorithm [see, e.g., Kitagawa (1987, 1996) and Kitagawa and Sato (2001) and the references therein] to estimate the log-likelihood of our model for share prices under the real-world measure. Kitagawa’s algorithm, used to compute the log-likelihood for nonlinear, non-Gaussian state-space models, employs the fundamental principles of particle-filtering. The essence of Kitagawa’s approach uses the weights described in equation (18) of Algorithm 1 to provide a Monte Carlo based approximation to the log-likelihood. The details of this log-likelihood approximation are available in Kitagawa and Sato (2001). We provide Kitagawa’s log-likelihood estimation approach for nonlinear, non-Gaussian state-space models in Algorithm 6.

{remark*}

The approximation to the log-likelihood value associated with a particular parameter value in equation (29) of Algorithm 6 becomes more accurate as the number of particles gets large. (We used in our estimation exercise.) When resampling using the weights in equation (18), we sometimes work with the shifted log-weights as this leads to improved sampling efficiency.

Once we obtain an estimate of the log-likelihood for the model parameters, we then combine it with the (log) priors and use a random-walk Metropolis algorithm to estimate the (log) posterior distribution of (or, equivalently, ). That is, we estimate the posterior distribution for and then transform back to the original scale to calculate the posterior distribution of as stated in equation (28). Our Markov chain Monte Carlo (MCMC) algorithm utilizes a Gaussian proposal density in order to facilitate the estimation of the posterior distribution in the parameter space. The details of our MCMC sampler are described in Algorithm 7, which is found in Section .2.

Given the estimates of the model parameters under the real-world measure using the share prices, we now describe how to use both share and option price data to estimate the market price of volatility risk. We work with a summary measure of , namely, the posterior mean, although one could use another measure such as the posterior median. The estimation of the market price of volatility risk will be computed conditional on this posterior summary measure. We implement this estimation exercise using the algorithms outlined in Section 4.

### 5.2 Volatility risk estimation

The estimation of the market price of volatility risk, , in equation (3) presents a computational challenge, particularly in the American option valuation framework. We aim to use the algorithms described in this paper to propose methodology that will facilitate statistical inference of in the American option pricing setting. To the best of our knowledge, this is the first analysis to compute and make posterior inference on the volatility risk premium for American-style (early-exercise) options. In many applications of option-pricing in a stochastic volatility framework, it is often the case that is set to a prespecified value [see, e.g., Heston (1993), Hull and White (1987) and Pastorello, Renault, and Touzi (2000)]. One convenient approach is to set . This is known as the “minimal martingale measure” in some strands of the option-pricing literature [Musiela and Rutkowski (1998)].

Both Pan (2002) and Eraker (2004) estimate stochastic volatility model parameters, including the market price of volatility risk, for the case of European options. However, the early-exercise feature of American-style options adds further complexities to the estimation problem. One method to estimate the market price of volatility risk for American-style options would be to set up the following nonlinear regression model [similar in spirit to what Eraker (2004) does in his analysis of European options]. Let

 Ui=Piθ∗(λ)+εi, (30)

where is the observed American option price, is the model predicted American option price conditional on the mean, , of the posterior distribution in equation (28), and is the market price of volatility risk.666Although it is not explicitly stated in equation (30), both