Stochastic Volatility Surface Estimation - Stanford Mathematics
Stochastic Volatility Surface Estimation Suhas Nayak and George Papanicolaou May 9, 2006
Abstract We propose a method for calibrating a volatility surface that matches option prices using an entropy-inspired framework. Starting with a stochastic volatility model for asset prices, we cast the estimation problem as a variational one and we derive a Hamilton-Jacobi-Bellman (HJB) equation for the volatility surface. We study the asymptotics of the HJB equation assuming that the stochastic volatility model exhibits fast mean-reversion. From the asymptotic solution of the HJB equation we get an estimate of the stochastic volatility surface. We also incorporate uncertainty in quoted derivative prices through a penalty term, i.e. by softening the constraints in the HJB equation. We present numerical solutions of our estimation scheme. We find that, depending on the softness of the constraints, certain parameters of the volatility surface related to the implied volatility smile can be calibrated so that they are stable over time. These parameters are essentially the ones found in previous fast mean-reversion asymptotics papers by Fouque, Papanicolaou and Sircar. We find that our procedure provides a natural way of interpolating between the prior parameters and the parameters of Fouque, Papanicolaou and Sircar.
1
Introduction and Review
Parameter identification for systems governed by partial differential equations is a well-studied inverse problem. In mathematical finance, finding the volatility of a risky asset from options with multiple strikes and maturities is an example. The information contained in the market is not enough to identify a pricing model and so many sets of model parameters and many types of models could potentially be compatible with the observed option prices. We study the problem of estimating a volatility surface from option prices in an incomplete market setting.
1
1.1
Complete Markets
In the complete markets case, there have been several approaches to estimating volatilities from observed option prices. One may try to use the Black-Scholes partial differential equation 1 Ct + σ(S, t)2 CSS + rCS − rC = 0, t < T 2 C(S, T ) = h(S) directly to estimate σ(S, t). This was done in Andersen and Brotherton-Ratcliffe ([2], [3]). Here, h(S) is the payoff of the option. Penalization criteria like smoothness norms were used by Lagnado and Osher [18] and Jackson et. al. [16] to regularize the volatility extraction procedure. Another possible approach in the complete market setting is to use a dual of the option-pricing PDE. Dupire’s equation (see [10], [11]) determines the dependence of option prices on strikes and expiration dates. It has the form 1 ∂2C ∂C ∂T C − σ 2 (T, K)K 2 + rK = 0, T > t 2 ∂K 2 ∂K C(t, K) = h∗ (K) Here, h∗ (K) has the form as h(S). Since this equation has derivatives in K and not in S, it is easier to handle because derivative prices are quoted for certain expiration times and strikes. Achdou and Pironneau [1] used this equation in combination with least-squares and a regularization (penalty) term to solve the inverse problem of estimating σ(T, K). Their objective was to minimize X J(σ) = |C(ti , Ki ) − ci |2 + Jr (σ) i
over σ subject to C solving Dupire’s equation. Here the {ci } are observed derivative prices and Jr (σ) is an appropriate Tychonoff regularization functional that involves the L 2 -norm of σ and the L2 -norm of its derivatives with respect to K and T .
1.2
Entropy-based methods
Regularization may also be achieved through the use of entropy. Entropy minimization for calibrating one-period asset pricing models was used by Buchen and Kelly [7], by Gulko [15], by Jackwerth and Rubinstein [17] and by Platen and Rebolledo [19]. Relative entropy-based methods can be motivated as follows. There is a range of strategies available to hedge an option. If the hedger is fully confident that the option will follow dynamics given by a known model, then it would be worthwhile to ignore observed market volatilities and hedge using the model. If, on the other hand, the hedger believes that current prices of the option fully determine the future evolution of the stock, then it would be a good idea to hedge according to 2
these current option prices. Relative entropy-based methods help bridge the gap between these two possible strategies. They provide estimates for model parameters that are close to prior information while still matching current option prices. Within the class of relative entropy methods there are two general approaches. In Avellaneda [4], a probability law is found for the risky asset that satisfies certain moment constraints, namely that it matches observed market prices of options, and is close to a a prior probability law, which can arise from historical or other econometric information. The closeness to the prior is measured by the relative entropy, which is given by dP P H(P |P0 ) = E ln dP0 dP is the Radon-Nikodym derivative of P with respect to where P0 is the prior distribution, and dP 0 P0 . For H(P |P0 ) to be finite, the measure P has to be absolutely continuous with respect to P 0 . So, in particular, if under P0 , we assume our asset prices follow
dXi (t) =
ν X
(0)
(0)
σij dBj (t) + µi dt
i=1
where {Bj } are standard Brownian motions, then, by Girsanov’s theorem, we must have that under P , the asset follows the price process dXi (t) =
ν X
(0)
σij dBj (t) + µi dt.
i=1
P (0) (0) Here µi = µi + j σij mj and mj is a market price of risk. In other words, once the volatility of the asset under the prior probability law is specified, the volatility of the asset under the new measure must be the same. Only the market price of risk can be adjusted to give a consistent pricing measure. In this framework we assume that many different, perhaps correlated, shocks drive the evolution of each asset. Moreover, we assume that there is a good prior that describes the effect of fluctuations (0) through the volatilities σij . It is unlikely that we will have an accurate assessment of these volatilities. Yet once they are fixed in the model we are unable to deviate from them in the new probability law because of the absolute continuity that is assumed. In order for the market price of risk to correct for possible misspecification of volatilities, wild swings in its value may be necessary. The work of Carmona and Xu [8], who introduce a stochastic volatility model in an entropy framework, suffers from the same problems. There, the market price of risk associated with the stochastic volatility process is the only degree of flexibility and all other parameters are fixed prior to the estimation procedure. It is desirable, therefore, to use a different relative entropy approach, one that was introduced in [5]. They considered processes for equity of the form dSt = σ0,t dBtP0 + µPt 0 dt, under P0 St 3
and
dSt = σt dBtP + µPt dt, under P. St
Unless σ0 = σ under P , the relative entropy of the two measures is infinite, since P and P 0 are then mutually singular. Avellaneda et. al. in [5] extended the concept of relative entropy to probability measures under which the processes do not necessarily have the same volatility. They considered the most singular part of the relative entropy using a time discretization that is based on trinomial trees. With this discretization and a small σ − σ0 expansion, they found that to highest order the relative entropy looks like (σ 2 − σ02 )2 over each time step. Both approaches to relative entropy calibration have the same general objective. Once the form of the relative entropy is determined the objective is to minimize it over all possible P subject to the constraint that the prices of options under P match observed prices.
1.3
Stochastic volatility surface estimation
Finding volatilities across strikes and expiration dates for incomplete markets is a very difficult task. We focus our attention on stochastic volatility models. These models have a large number of parameters that need to be known for pricing purposes and options can be quite sensitive to them. There are, however, several papers that deal with this issue. Broadie et al. [6], for example, fit the parameters of various stochastic volatility models. They first obtain parameters using a long-run time series of underlying asset returns. They then use the information found in option prices to estimate volatility and risk premia. This second stage involves the minimization of an objective function that is just the sum of the squares of the difference in the model-derived Black-Scholes implied volatilities and the implied volatilities that correspond to the data. Their method forces consistency in parameters between the data for the underlying asset and the data for the option prices. The problem of option pricing in a stochastic volatility setting was tackled in a series of papers by Fouque, Papanicolaou and Sircar (see [12] and the references therein). They developed a method that reduced the number of parameters that were needed for pricing and hedging to just three. Two were derived directly from the smile associated with options prices, and the other was some estimate of the underlying’s volatility. They introduced a type of model where the stochastic volatility factor followed a fast mean-reverting process. The data set of options they then considered was chosen so that the asymptotics would be valid. For example, the strikes of the options were near at-the-money and the options had expiration dates that were neither close nor far away. With this data set they found that their parameters were stable over time. Importantly, the option pricing equations they obtained within the fast mean-reverting regime were independent of the stochastic volatility driving factor. This was not the case in the entropy-based estimation procedure of Carmona and Xu [8]. A subsequent paper of Fouque et al. [14] extended the domain of applicability by introducing structure to the maturity cylces of options. We will not pursue that method here. In this paper, we introduce a method for calibrating volatility surfaces in a stochastic volatility 4
environment. We use an entropy-inspired framework that allows flexibility in volatilities when passing from the prior probability law to the pricing one. Because we are unable to observe the underlying stochastic volatility process, we adopt the framework of Fouque, Papanicolaou and Sircar and consider a regime where fast mean-reversion is apparent. We also relax the pricing constraints so that our volatility surface is less sensitive to unreliable out-of-the-money option prices. We formulate and solve a Hamilton-Jacobi-Bellman (HJB) equation associated with the variational problem. An HJB equation in a fast mean-rversion setting was studied by Sircar and Zariphopoulou [20] in the context of portfolio optimization. Our problem is different and results in a different HJB equation. From this equation, we obtain a volatility surface that consistently incorporates the Fouque-Papanicolaou-Sircar stochastic volatility parameters. These parameters are found to be stable over time. Our method, in fact, permits a natural interpolation between the model parameters of Fouque, Papanicolaou and Sircar and our own prior. It is an interpolation that depends on the level of pricing uncertainty. The rest of this paper is organized as follows. In Section 2 we formulate the problem and introduce a method for obtaining the volatility surface. We review the fast mean-reversion regime and simplify the corresponding HJB equation in Section 3. We reduce the complexity of the problem by studying fast mean-reversion asymptotics in Section 4. In Section 5 we describe the iterative procedure that was used to estimate the surface, while in Section 6 we detail the numerical methods we used to solve the resulting partial differential equations. In Section 7 we present the results of our method. We discuss the sensitivity of the method to the various parameters in Section 8. We end with a brief summary and conclusion.
2 2.1
Stochastic volatility in the relative entropy framework A simplified form for the relative entropy
In the relative entropy framework that we follow, we allow the functional form of the volatility to change from P0 , our prior probability law, to P , the probability law we wish to determine, because otherwise we would only have freedom to pick a risk premium. We therefore adapt the method of Avellaneda et. al. to the case of stochastic volatility. In order to do this, we discretize in time the stock and stochastic volatility processes so that we may determine the most singular part of the relative entropy. We assume the following dynamics under P0 : dSt = rdt + σ0 (St , Yt )dWt St dYt = γ0,t dt + κdZˆt where, as before, we assume the correlation between the two Brownian motions is ρ. In other words, dZˆt dWt = ρdt. 5
Here, r is the riskless interest rate process and σ0 is our prior estimate of the volatility surface. Yt is the stochastic volatility driving factor, which has a drift of γ0 and has its own volatility κ. Under P , we assume that the dynamics are given by: dSt = rdt + σ(St , Yt )dWt∗ St dYt = γt dt + κdZˆt∗
(1)
Here, γ and γ0 are functions that may depend on both t and Yt , and again we assume the same level of correlation between the Brownian motions. Details of the procedure used to get the relative entropy given the dynamics above are provided in Appendix A. The result is that the relative entropy generated per unit of time is just η(σ) = (σ 2 − σ02 )2
when σ 2 is not very far from σ02 . This particular form is an unsurprising consequence of approximating a distance function.
2.2
Formulation of the variational problem
Since we wish to minimize the relative entropy between our prior and a model that fits the data, we consider the optimization problem Z T P η(σs )ds sup −E σ
0
subject to the constraints:
P −rT E e i hi (ST ) − Ci ≤ βi , i
for each i and where βi ≥ 0. Here, hi (STi ) is the payoff of option i, which has expiration date Ti . Ci is the price of this option that is quoted in the market. The constants βi may be considered to be a measure of our confidence in the market data. We note that in previous works involving relative entropy minimization, these constants were taken to be zero (hard price constraints). By allowing deviation from market prices under our new measure P , we have softened the pricing constraints, since it is not practical for us to consider the market data as being exact. Some data points may be the result of only one trade and hence would not be reliable. Imposing hard constraints also leads to kinks in the smile near the data points. We simplify the objective above by forming a single Lagrangian. The following analysis is similar to the one in [9] for a different problem. Our objective may be rewritten as Z T inf E P e−rs η(σs )ds subject to σ 0 (λ+ ∀i : E P e−rTi h(STi ) − Ci ≤ βi i ) (λ− ∀i : Ci − E P e−rTi h(STi ) ≤ βi i ). 6
Here, Lagrange multipliers are written next to the constraints. From this, we obtain the augmented Lagrangian which we try to minimize over both sets of λ: X Z T X −rT − + P −rs P i (2) )(E e h(S ) − C ) + βi (λ+ − λ e η(σs )ds + (λ− sup −E i T i i + λi ). i i σ
0
i
i
− Without loss of generality, we may assume that, for each i, at most one of λ + i and λi is nonzero. This is because we are minimizing over λi , and so we could otherwise decrease both λi ’s by an equal but positive amount. We would then obtain a new expression where the third term in Equation (2) is decreased but the other two are unchanged.
Now suppose that λ+ i is nonzero (at the minimum). Then our solution is the same as the one for an optimization with the first constraint replaced by equality for that value of i. In this case, the second constraint cannot simultaneously be an equality and the solution must therefore be an interior one. Similarly, if the λ− i is nonzero, then the second constraint is an equality. These are the − only two possible cases. If we now write for each i, λi = λ+ i − λi , we may consider the objective: sup −E σ
P
Z
T
e
−rs
η(σs )ds +
0
X i
X λi (E P e−rTi h(STi ) − Ci ) + βi |λi |.
(3)
i
If we have to minimize this over λi (for each i), then there are two cases. First, the λi , for some i, may be positive, in which case we have solved the original optimization problem with the first inequality replaced by equality. Or, the λi , for that same i, is negative, in which case the original problem with the second constraint replaced by equality has been solved. Hence, the two problems given in Equations (2) and (3) are equivalent, as they both reduce to the same cases once minimization over λi has been performed. We may write an indirect value function Z T X i h λ P −r(s−t) V (t, s, y) = sup −Et e η(σs )ds + λi EtP e−r(Ti −t) h(STi ) . σ
t
i
We note that our value function discounts the entropy to make the subsequent expressions simpler. Our objective in Equation (3) amounts to minimizing the expression X X V λ (0, S, y) − λi C i + βi |λi | (4) i
i
over λ. We know that this V λ solves the HJB equation that is derived in Appendix B and has the form X 1 1 λi δ(t − Ti )hi (S). Vt − rV + rSVS + γVy + κ2 Vyy + Φ( S 2 VSS , ρκSVSy ) = − 2 2 i
V (T, S, y) = 0
7
Here, Φ(X1 , X2 ) = sup X1 σ 2 + X2 σ − η(σ). σ
Before describing how to obtain a volatility surface from these equations, it is worthwhile noting some features of the HJB equation. With ρ = 0, the HJB equation becomes much simpler because Φ(X, Y ) may be directly evaluated. The equation decouples and the value function would solve the same PDE that was found in Avellaneda et al. [5]. The volatility surface may be derived from the value function as follows. Once we solve for the λ = λ∗ that minimizes Equation (4), we may write the desired volatility surface as 1 λ∗ 2 λ∗ σ ∗ (S, y, t) = arg sup S 2 VSS σ + ρκSVSy . σ 2
(5)
We have thus described a way to incorporate stochastic volatility in a relative entropy framework. We did this by first determining an appropriate form for the relative entropy. Once we did this, however, encoding our objective into a value function is standard, and this value function must solve an HJB equation. Apart from the form of the relative entropy, our other main contribution so far is the incorporation of pricing uncertainty through the parameters β i . Although this methodology seems reasonable, there are some issues that need to be addressed. The HJB equation as written is highly dependent on many prior parameters. We need to know γ, κ and ρ accurately in order to proceed. Since these are all parameters associated with the unobservable stochastic volatility driving factor, determining them is a difficult task. The volatility surface, moreover, is dependent on the actual value of the stochastic volatility driving factor. The next section develops a way to get around these difficulties by considering a regime in which the stochastic volatility driving factor is fast mean-reverting. This represents the point of departure from the standard HJB equation methodology.
3 3.1
The fast mean-reversion setting A special form for the dynamics of the stochastic volatility driving process
We have so far dealt with the presence of stochastic volatility in a fairly general way. We now consider some special models for the dynamics of the volatility driving factor, Y . One feature that most models of stochastic volatility incorporate is mean reversion. This refers to the tendency of the process to go back to its invariant or long-run distribution. A particularly tractable model that exhibits mean reversion is the Ornstein-Uhlenbeck process dYt = α(m − Yt )dt + κdZˆt . 8
Here, α is the rate of mean reversion, while m is the long-run mean of the process. We suppose that α, κ and m are constants. Moreover, Zˆt is a Brownian motion that is correlated to the Brownian motion that drives the stock price process, Bt , with correlation ρ. The invariant distribution for this process is a Y0 that satisfies E[Lg(Y0 )] = 0 for any smooth and bound g, where L = α(m − y)
∂ 1 ∂2 + κ2 2 . ∂y 2 ∂y
If we let Ψ(y) be the density function of the invariant distribution, it is easily seen that (y − m)2 1 exp − Ψ(y) = √ , 2ν 2 2πν 2 where ν2 =
(6)
κ2 . 2α
This is exactly the normal density with mean m and variance ν 2 , which tells us that the parameter ν controls the size of the equilibrium fluctuations from the mean, m. The fast mean reversion limit was considered in Fouque, Papanicolaou and Sircar [12]. They found empirical evidence for fast mean reversion in options data. Analytically, fast mean reversion corresponds to the limit α → ∞. Care, however, must be taken to ensure that the invariant distribution remains the same as we take this limit. In other words, we would like our model to have fluctuations of the same size regardless of how quickly the volatility reverts to its mean. We therefore take √ κ = ν 2α for some constant ν. If this is the specification of Y under the physical measure, then to price derivatives we need to consider the process under the risk-neutral law. Suppose that under some physical probability law the stock and volatility driving process follow dSt = µSt dt + σ0 (St , Yt )dB0,t √ dYt = α(m − Yt )dt + ν 2αdZˆ0,t . Then, under a risk-neutral law, the processes must follow, by Girsanov’s theorem, ∗ dSt = rSt dt + σ0 (Yt )dB0,t p √ √ √ µ−r ∗ 2 − ν 2αγ0,t 1 − ρ dt + ν 2αdZˆ0,t dYt = α(m − Yt ) − ν 2αρ . σ0
9
(7)
Here, γt is a risk-premium factor or the market price of volatility risk that parametrizes the space of equivalent martingale measure. By taking, γ0,t = γ0 (t, St , Yt ), we specifically restrict ourselves to a Markovian setting. Equation (7) is the specification of the dynamics of the processes S and Y under our prior risk-neutral probability law. We now just need to find the relative entropy between another risk-neutral probability law and our prior. We assume that under P the price of the stock evolves as dSt = rSt dt + σ(St , Yt )dBt∗ p √ √ √ µ−r 2 dYt = α(m − Yt ) − ν 2αρ − ν 2αγt 1 − ρ dt + ν 2αdZˆt∗ . σ
(8)
With these particular assumptions on the dynamics of Yt , the form of η(σ) to first order in ρ is unchanged from the analysis in Appendix A. It is worth noting one important aspect of the dynamics under the prior probability law. We have assumed that σ0 (St , Yt ) = σ0 (Yt ). In other words, our prior estimate of the volatility is dependent only on the stochastic volatility driving factor and not on the stock price. With this choice of the functional form σ 0 now corresponds to function f in the work of Fouque, Papanicolaou and Sircar [12]. Although this particular form seems like a heavy restriction, it is not. The reason is that the fast mean-reversion setting does not need an explicit functional form for the prior volatility estimate. Instead, this setting replaces explicit functional forms with observable quantities and averages. This will be made clear later.
3.2
Simplifying the HJB equation
The HJB equation we need to solve is now a little different to the one we derived earlier. Because of the appearance of the σ in the denominator of one of the diffusive terms for Y t in Equation (8), Φ must now be a function of three variables. Specifically, we define Φ as: 1 2 Φ(X1 , X2 , X3 ) = sup σ X1 + ρX2 σ + ρX3 − η(σ) . (9) σ σ The HJB equation thus becomes: p √ Vt − rV + rSVS + (α(m − y) − ν 2αγ 1 − ρ2 )Vy + ν 2 αVyy X √ √ 1 2 +Φ λi δ(t − Ti )hi (S) S VSS , ν 2αSVSy , −ν 2α(µ − r)Vy = − 2
(10)
i
Before we carry out the asymptotic analysis we simplify Φ for σ close to σ 0 and for ρ small. We therefore expand σ in ρ (a small ρ expansion). Let: σ=σ ˜ + Aρ 10
where A is a coefficient to be determined. Here, σ ˜ is taken to be the maximizer when ρ = 0, which is just given by σ ˜ 2 = X1 /2 + σ02 . The first-order condition for a maximum gives: 2σX1 + ρX2 − ρ
X3 − 4σ(σ 2 − σ02 ) = 0. σ2
Upon substituting for σ, we find that the coefficient for ρ is 2AX1 + X2 −
X3 X3 − 4A(˜ σ 2 − σ02 ) − 8A˜ σ 2 = X2 − 2 − 8A˜ σ2. σ ˜2 σ ˜
Equating the coefficient of ρ to zero yields: X2 −
X3 − 8A˜ σ 2 = 0, σ ˜2
which implies that 3 X2 − X σ ˜2 . A= 4X1 + 8σ02
So, to first order in ρ: Φ(X1 , X2 , X3 ) =
X12 X3 + σ02 X1 + ρX2 σ ˜+ρ 4 σ ˜
which, upon expanding σ ˜ , leads to: X2 Φ(X1 , X2 , X3 ) = 1 + σ02 X1 + ρX2 4
X1 + σ02 2
1
2
+ ρX3
X1 + σ02 2
− 1 2
.
Given this form of Φ, we note that the maximizing value of σ is clearly approximated by σ ˜. We remark also that the small ρ expansion may be justified a posteriori. We will find that the y-dependence in our value function only occurs at the α1 scale. This implies that, since we find √ X2 ∝ αVxy , X2 scales like √1α , which is small for large α. A similar argument holds for X3 . Because ρ always premultiplies X2 or X3 , the small ρ expansion of Φ is reasonable. We now have a suitable form of the HJB equation, which uses the small ρ approximation for Φ. 2 p √ 1 1 1 2 2 2 Vt − rV + rSVS + (α(m − y) − ν 2αγ 1 − ρ )Vy + ν αVyy + S VSS + σ02 S 2 VSS 4 2 2 1 − 1 √ √ 2 2 1 2 1 2 + ρν 2αSVSy S VSS + σ02 − ρν 2α(µ − r)Vy S VSS + σ02 4 4 X =− λi δ(t − Ti )hi (S) i
V (0, S, y) = 0
(11)
11
4
Fast mean-reversion asymptotics applied to the HJB equation
In order to apply fast mean-reversion asymptotics, we suppose that α is large. For book-keeping purposes, we replace α by αε . We consider the value function corresponding to this situation in the limit as ε → 0. We first let x = log S (which means that 21 S 2 VSS = 12 (Vxx − Vx )). Substituting this into Equation (11) yields p √ 1 1 Vt − rV + (r − σ02 )Vx + σ02 Vxx + (α(m − y) − ν 2αγ 1 − ρ2 )Vy + ν 2 αVyy 2 2 1 − 1 √ √ 2 2 1 Vxx − Vx Vxx − Vx 2 2 2 + σ0 + σ0 + (Vxx − Vx ) + ρν 2αVxy − ρν 2α(µ − r)Vy 16 4 4 X =− λi hi (x)δ(t − Ti ). (12) i
Here, σ0 = σ0 (Yt ). We now replace α by αε so that we may apply some asymptotic techniques. We therefore write our equation as √ 1 1 (Vxx − Vx )2 1 + √ ρν 2αVxy (g(Vxx , Vx , y) − 1) L2 + √ L1 + L0 V + ε 16 ε ε X √ 1 1 − √ ρν 2αVy −1 − λi hi (x)δ(t − Ti ) (13) g(Vxx , Vx , y) ε i
where ∂2 ∂ 1 ∂ 1 − r + (r − σ0 (y)2 ) + σ0 (y)2 2 ∂t 2 ∂x 2 ∂x p √ √ √ ∂2 ∂ ∂ − ν 2αγ 1 − ρ2 − ρν 2α(µ − r) L1 = ρν 2α ∂x∂y ∂y ∂y 2 ∂ ∂ L0 = α(m − y) + ν 2α 2 ∂y ∂y 1 2 Vxx − Vx . g(Vxx , Vx , y) = + σ02 4 L2 =
V
(14)
−V
If we suppose that 0,xx4 0,x is small in relation to σ0 , which is reasonable as long as our prior is somewhat close to the actual σ, then we may approximate g as just g(y) = σ0 (y).
(15)
This approximation is consistent with the other approximations we have made in deriving a form of the relative entropy. We note that this approximation for g is also what we would have if we had happened to pick a prior that matched the prices within the tolerance desired. 12
We also expand V in powers of ε V = V0 +
4.1
√ √ εV1 + εV2 + ε εV3 + ....
The leading order term
In order to solve for the V0 term we equate the expansions of the HJB equation at the first few scales. At the 1ε scale L0 V0 = 0, which implies that V0 = V0 (t, x) (it is independent of y). At the
√1 ε
scale
√ √ L0 V1 + L1 V0 + ρν 2αV0,xy (g(V0,xx , V0,x , y) − 1) − ρν 2α(µ − r)V0,y (g(V0,xx , V0,x , y)) = 0 but since V0 is independent of y, the equation becomes L 0 V1 = 0 and this implies, as before, V1 = V1 (t, x). At the order 1 scale L 2 V0 + L 1 V1 + L 0 V2 +
X 1 (V0,xx − V0,x )2 = − λi hi (x)δ(t − Ti ), 16 i
where, again, we have used that V1 and V0 is independent of y to get rid of the Vxy part of the nonlinearity. Similarly, we may get rid of the L1 V1 = 0. Solvability for V2 requires that the other terms have average 0. But since V0 is independent of y, this means that the equation that determines solvability yields the following equation for V0 X 1 LBS (¯ σ0 )V0 + (V0,xx − V0,x )2 = − λi hi (x)δ(t − Ti ). (16) 16 i
Here, LBS is the Black-Scholes operator after the log transformation has been performed.The terminal condition is V0 (T, x) = 0. This is exactly the equation that is solved in Avellaneda et. al. [5]. It appears here as the zero order approximation in our method.
4.2
The correction term
Once we have solved the equation for V0 , we may remove the average (which we have set to zero) to get an equation for V2 2 1 ∂ V0 ∂V0 2 2 L0 V2 + (σ0 (y) − σ = 0, ¯ ) − 2 ∂x2 ∂x 13
which implies that 1 V2 (t, x, y) = − L−1 (σ0 (y)2 − σ ¯2) 2 0 We may therefore write V2 as 1 V2 (t, x, y) = − (φ(y) + c(t, x)) 2
∂ 2 V0 ∂V0 − ∂x2 ∂x
∂ 2 V0 ∂V0 − ∂x2 ∂x
.
,
(17)
where φ(y) is a solution of the Poisson equation L0 φ = σ0 (y)2 − < σ02 > . √ Moving on to the next scale ( ε), the relevant equation is √ 1 L2 V1 + L1 V2 + L0 V3 + (V0,xx − V0,x )(V1,xx − V1,x ) + ρν 2αV2,xy (g(V0,xx , V0,x , y) − 1) 8 √ − ρν 2αV2,y (1/g(V0,xx , V0,x , y) − 1) = 0. Solvability for V3 implies √ 1 ∂ 2 V2 LBS (¯ σ0 )V1 + (V0,xx − V0,x )(V1,xx − V1,x ) + ρν 2α g(V0,xx , V0,x , y) 8 ∂x∂y p √ √ ∂V2 1 ∂V2 − 2α 1 − ρ2 γ = 0. − ρν 2α(µ − r) g(V0,xx , V0,x , x) ∂y ∂y
(18)
The averages in the last three terms of Equation (18) may be further expanded into terms with derivatives in V0 , using Equation (17). Specifically, we obtain the following relations 3 D E 2 ∂ 2 V2 g(V0,xx , V0,x , y) ∂x∂y = 12 hgφ0 i ∂∂xV30 − ∂∂xV20 E D E 2 D ∂V0 ∂V2 ∂ V0 1 1 0 1 = (19) φ − 2 2 g ∂x g(V0,xx ,V0,xD ,x) ∂y E 2∂x 2 0 γ ∂V = 21 hγφ0 i ∂∂xV20 − ∂V ∂y ∂x So our equation for the correction V1 becomes
2 1 ∂ 3 V0 ˜2 ∂ V0 + A˜1 ∂V0 = 0, LBS (¯ σ0 )V1 + (V0,xx − V0,x )(V1,xx − V1,x ) + A˜3 + A (20) 8 ∂x3 ∂x2 ∂x which is a linear equation with a source term. In particular, the coefficients A˜1 , A˜2 and A˜3 are given by
√ 1 p √ 1 1 0 ˜ φ + ν α√ 1 − ρ2 γφ0 A1 = ρν α √ (µ − r) g 2 2
√ 1 0 √ 1 √ 1 p 1 0 ˜ A2 = −ρν α √ gφ − ρν α √ (µ − r) 1 − ρ2 γφ0 φ − ν α√ g 2 2 2
√ 1 A˜3 = ρν α √ gφ0 2
14
We note that, with this method, and as it stands now, we need the correlation ρ and the various parameters of the stationary distribution of the stochastic volatility driving factor.
4.3
Determination of the group parameters
With the approximation given in Equation (15), our parameters may be written as
√ √ 1 p √ 1 1 0 A1 = ε ρν α √ (µ − r) φ + ν α√ 1 − ρ2 γφ0 σ0 2 2
√ √ 1
√ 1 √ 1 p 1 0 0 0 A2 = − ε ρν α √ σ0 φ + ρν α √ (µ − r) φ + ν α√ 1 − ρ 2 σ0 φ σ0 2 2 2 A3 =
√ 1
√ ερν α √ σ0 φ0 . 2
(21)
These parameters, which we call our group parameters, are the same as the A˜i but they are now √ scaled by ε. We may determine these parameters as combinations of parameters that are observable from the smile. To do this, we outline the results in Fouque, Papanicolaou and Sircar [12]. In their work, they found that the price of an option, C may be expanded as ∂CBS ∂ 2 CBS ∂ 3 CBS C(t, x; T, K) = CBS (t, x; T, K) − (T − t) (2V3 − V2 ) , + (V2 − 3V3 ) + V3 ∂x ∂x2 ∂x3 (22) where CBS is the Black-Scholes option price and x = log(S). Here,
√ √ 1 √ 1 p √ 1 1 0 0 0 2 √ √ √ V2 = ε ν α 1 − ρ σ0 φ φ −ν α 2ρ σ0 φ − ρν α (µ − r) σ0 2 2 2 √ 1
√ (23) V3 = ερν α √ σ0 φ0 2 We deduce, then, that our approximate Ai ’s given by Equation (21 )are exactly A1 = 2V3 − V2
A2 = V2 − 3V3
A3 = V 3 .
The parameters V2 and V3 may be determined from the smile. By expanding the price of the option in powers of ε, we may write C(t, x; T, K; I) = CBS (t, x; T, K; σ0 ) + 15
√ ∂CBS εI1 , σ
(24)
where I1 is the correction to the implied volatility. By equating Equations (22) and (24), we find that the implied volatility I is just: I = σ0 +
3 V3 V2 V3 log(K/S) (r + σ02 ) − − 3 . 2 2 σ0 σ0 T − t σ0
The implied volatility is therefore an affine function of the log-moneyness-to-maturity-ratio LMMR =
log(K/S) , T −t
which means that we may fit the parameters aF P S and bF P S below to the smile curve determined by options prices I = aF P S LMMR + bF P S , where aF P S = − bF P S
V3 σ03
V3 = σ0 + 3 σ0
3 r + σ02 2
−
V2 . σ0
This means that the scaled parameters Ai which are required to determine the correction term √ εV1 are explicitly expressed as combinations of the parameters aF P S and bF P S . Moreover, given any such parameters a and b determined from a smile curve, the parameters A i are given by 3 A1 = −2aσ03 − σ0 ((σ0 − b) − a(r + σ02 )) 2 3 2 A2 = σ0 ((σ0 − b) − a(r + σ0 )) + 3aσ03 2 A3 = −aσ03
(25)
We have thus expressed the parameters needed for the correction term as simple expressions involving observable quantities. This is the power of the fast mean-reversion setting. Individual estimates of the parameters of the stochastic volatility models are not needed. Rather, group parameters like the Ai are all that are needed to not only price options as in Fouque, Papanicolaou and Sircar [12] but also to determine minimal entropy volatility surfaces as in this paper.
4.4
Summary of our analysis
What we have achieved now is a set of equations that incorporates the objectives we wanted 1. Uncertainty over the prior in the form of stochastic volatility that is independent of the specific form of the volatility since we only really need the group parameters, A 1 , A2 and A3 . 16
2. A value function that does not depend on the stochastic volatility driving factor because of the use of fast mean-reversion. 3. Uncertainty in the observed market prices The asymptotics helped us isolate the leading two terms of the value function that solved the HJB equation. These two terms were, indeed, independent of the stochastic volatility driving factor. The first term, moreover, solved the same equation that was solved in Avellaneda et al. [5].
5
Description of the estimation procedure
Our analysis permits a simplified approach to the estimation of the stochastic volatility surface. There are two possible approaches from here. We may estimate the parameters A i once or we may try to self-consistently determine the same parameters in an iteration scheme. In the first approach, we may determine A1 , A2 and A3 by following the smile-fitting procedure given in Section 4. We may then solve for V0 using Equation (16). With this V0 and the coefficients √ {Ai }, we may solve for εV1 for each value of λ using search P(20). Using some gradient P Equation √ algorithm, we may then minimize (over λ) V (0, x0 )− i λi Ci + i |λi |βi , where V = V0 + εV1 and x0 is the log of the current stock price. Finally, the appropriate derivative of V gives a volatility surface as described by Equation (5). We choose to iterate this scheme. We may obtain A1 , A2 and A3 from the prices as before, solve for V0 and V1 and obtain a volatility surface. But instead of stopping there, we price the same set of options in our data set again using the newly derived volatility surface and use these prices to get adjusted estimates for A1 , A2 and A3 . Repeating this procedure a few times then leads to a converged set of parameters. The fact that our procedure leads to a converged set of parameters tells us the group parameters are intrinsic to the volatility surface as a whole and that they may be self-consistently determined. They thus have the same role in volatility surface estimation as they do in option pricing. We call this method the volatility surface iteration procedure. It is used in our numerical calculations and is summarized succinctly as 1. Fit the implied volatility smile to an affine function of the log-moneyness-to-maturity ratio. This affine fit gives us the parameters aF P S and bF P S from which we obtain A1 , A2 and A3 . √ 2. Solve for V0 and εV1 . 3. Obtain a volatility surface. 4. Price options again using the volatility surface in Step 3. We do this using the Black-Scholes PDE with volatility σ. 17
5. Derive the implied volatility smile from the prices obtained in Step 4. Fit this smile to the log-moneyness-to-maturity ratio to obtain parameters a and b. Obtain new estimates of the parameters A1 , A2 and A3 from the fitted a and b. . 6. Stop if adjusted parameters differ from the previous iteration by less than δ, a small parameter (which we take to be of order 10−5 in the l2 norm). Otherwise, repeat the procedure from Step 2.
6 6.1
Numerical Methods The highest-order equation
Supposing that the stock pays dividend of d (constant), then the function V0 (x, t) satisfies the equation: 1 2 2 2 2 sup V0,t − rV0 + (r − d)V0,x + σ (V0,xx − V0,x ) − (σ − σ ¯0 ) 2 σ X =− λi hi (ex )δ(t − Ti ) i
Here, we choose σ ∈ [σmin , σmax ] in order to preserve the stability of our numerical scheme. Now, using the substitution V˜0 = e−rt V0 , yields: 1 ¯02 )2 V˜0,t + (r − d)V˜0,x + σ 2 (V˜0,xx − V˜0,x ) − e−rt (σ 2 − σ sup 2 σ∈[σmin ,σmax ] X =− λi hi (ex )e−rt δ(t − Ti ) i
We discretize the interior equation (spatial h, temporal ∆t), and supposing that we have the optimal σ at time t, we have the equation (after dropping the tildes): V0 (x, t) − V0 (x, t − ∆t) 1 2 V0 (x + h, t) + V0 (x − h, t) − 2V0 (x, t) + σ ∆t 2 h2 1 V0 (x + h, t) − V0 (x − h, t) + (r − d − σ 2 ) − e−rt (σ 2 − σ ¯02 )2 2 2h X =− e−rt λi hi (ex )δ(t − Ti ) i
In order to regularize the problem, we replace each option payoff e −rt hi (ex )δ(t − Ti ) by Pi = C(ex , ∆t, σ)1t≤Ti
Abstract We propose a method for calibrating a volatility surface that matches option prices using an entropy-inspired framework. Starting with a stochastic volatility model for asset prices, we cast the estimation problem as a variational one and we derive a Hamilton-Jacobi-Bellman (HJB) equation for the volatility surface. We study the asymptotics of the HJB equation assuming that the stochastic volatility model exhibits fast mean-reversion. From the asymptotic solution of the HJB equation we get an estimate of the stochastic volatility surface. We also incorporate uncertainty in quoted derivative prices through a penalty term, i.e. by softening the constraints in the HJB equation. We present numerical solutions of our estimation scheme. We find that, depending on the softness of the constraints, certain parameters of the volatility surface related to the implied volatility smile can be calibrated so that they are stable over time. These parameters are essentially the ones found in previous fast mean-reversion asymptotics papers by Fouque, Papanicolaou and Sircar. We find that our procedure provides a natural way of interpolating between the prior parameters and the parameters of Fouque, Papanicolaou and Sircar.
1
Introduction and Review
Parameter identification for systems governed by partial differential equations is a well-studied inverse problem. In mathematical finance, finding the volatility of a risky asset from options with multiple strikes and maturities is an example. The information contained in the market is not enough to identify a pricing model and so many sets of model parameters and many types of models could potentially be compatible with the observed option prices. We study the problem of estimating a volatility surface from option prices in an incomplete market setting.
1
1.1
Complete Markets
In the complete markets case, there have been several approaches to estimating volatilities from observed option prices. One may try to use the Black-Scholes partial differential equation 1 Ct + σ(S, t)2 CSS + rCS − rC = 0, t < T 2 C(S, T ) = h(S) directly to estimate σ(S, t). This was done in Andersen and Brotherton-Ratcliffe ([2], [3]). Here, h(S) is the payoff of the option. Penalization criteria like smoothness norms were used by Lagnado and Osher [18] and Jackson et. al. [16] to regularize the volatility extraction procedure. Another possible approach in the complete market setting is to use a dual of the option-pricing PDE. Dupire’s equation (see [10], [11]) determines the dependence of option prices on strikes and expiration dates. It has the form 1 ∂2C ∂C ∂T C − σ 2 (T, K)K 2 + rK = 0, T > t 2 ∂K 2 ∂K C(t, K) = h∗ (K) Here, h∗ (K) has the form as h(S). Since this equation has derivatives in K and not in S, it is easier to handle because derivative prices are quoted for certain expiration times and strikes. Achdou and Pironneau [1] used this equation in combination with least-squares and a regularization (penalty) term to solve the inverse problem of estimating σ(T, K). Their objective was to minimize X J(σ) = |C(ti , Ki ) − ci |2 + Jr (σ) i
over σ subject to C solving Dupire’s equation. Here the {ci } are observed derivative prices and Jr (σ) is an appropriate Tychonoff regularization functional that involves the L 2 -norm of σ and the L2 -norm of its derivatives with respect to K and T .
1.2
Entropy-based methods
Regularization may also be achieved through the use of entropy. Entropy minimization for calibrating one-period asset pricing models was used by Buchen and Kelly [7], by Gulko [15], by Jackwerth and Rubinstein [17] and by Platen and Rebolledo [19]. Relative entropy-based methods can be motivated as follows. There is a range of strategies available to hedge an option. If the hedger is fully confident that the option will follow dynamics given by a known model, then it would be worthwhile to ignore observed market volatilities and hedge using the model. If, on the other hand, the hedger believes that current prices of the option fully determine the future evolution of the stock, then it would be a good idea to hedge according to 2
these current option prices. Relative entropy-based methods help bridge the gap between these two possible strategies. They provide estimates for model parameters that are close to prior information while still matching current option prices. Within the class of relative entropy methods there are two general approaches. In Avellaneda [4], a probability law is found for the risky asset that satisfies certain moment constraints, namely that it matches observed market prices of options, and is close to a a prior probability law, which can arise from historical or other econometric information. The closeness to the prior is measured by the relative entropy, which is given by dP P H(P |P0 ) = E ln dP0 dP is the Radon-Nikodym derivative of P with respect to where P0 is the prior distribution, and dP 0 P0 . For H(P |P0 ) to be finite, the measure P has to be absolutely continuous with respect to P 0 . So, in particular, if under P0 , we assume our asset prices follow
dXi (t) =
ν X
(0)
(0)
σij dBj (t) + µi dt
i=1
where {Bj } are standard Brownian motions, then, by Girsanov’s theorem, we must have that under P , the asset follows the price process dXi (t) =
ν X
(0)
σij dBj (t) + µi dt.
i=1
P (0) (0) Here µi = µi + j σij mj and mj is a market price of risk. In other words, once the volatility of the asset under the prior probability law is specified, the volatility of the asset under the new measure must be the same. Only the market price of risk can be adjusted to give a consistent pricing measure. In this framework we assume that many different, perhaps correlated, shocks drive the evolution of each asset. Moreover, we assume that there is a good prior that describes the effect of fluctuations (0) through the volatilities σij . It is unlikely that we will have an accurate assessment of these volatilities. Yet once they are fixed in the model we are unable to deviate from them in the new probability law because of the absolute continuity that is assumed. In order for the market price of risk to correct for possible misspecification of volatilities, wild swings in its value may be necessary. The work of Carmona and Xu [8], who introduce a stochastic volatility model in an entropy framework, suffers from the same problems. There, the market price of risk associated with the stochastic volatility process is the only degree of flexibility and all other parameters are fixed prior to the estimation procedure. It is desirable, therefore, to use a different relative entropy approach, one that was introduced in [5]. They considered processes for equity of the form dSt = σ0,t dBtP0 + µPt 0 dt, under P0 St 3
and
dSt = σt dBtP + µPt dt, under P. St
Unless σ0 = σ under P , the relative entropy of the two measures is infinite, since P and P 0 are then mutually singular. Avellaneda et. al. in [5] extended the concept of relative entropy to probability measures under which the processes do not necessarily have the same volatility. They considered the most singular part of the relative entropy using a time discretization that is based on trinomial trees. With this discretization and a small σ − σ0 expansion, they found that to highest order the relative entropy looks like (σ 2 − σ02 )2 over each time step. Both approaches to relative entropy calibration have the same general objective. Once the form of the relative entropy is determined the objective is to minimize it over all possible P subject to the constraint that the prices of options under P match observed prices.
1.3
Stochastic volatility surface estimation
Finding volatilities across strikes and expiration dates for incomplete markets is a very difficult task. We focus our attention on stochastic volatility models. These models have a large number of parameters that need to be known for pricing purposes and options can be quite sensitive to them. There are, however, several papers that deal with this issue. Broadie et al. [6], for example, fit the parameters of various stochastic volatility models. They first obtain parameters using a long-run time series of underlying asset returns. They then use the information found in option prices to estimate volatility and risk premia. This second stage involves the minimization of an objective function that is just the sum of the squares of the difference in the model-derived Black-Scholes implied volatilities and the implied volatilities that correspond to the data. Their method forces consistency in parameters between the data for the underlying asset and the data for the option prices. The problem of option pricing in a stochastic volatility setting was tackled in a series of papers by Fouque, Papanicolaou and Sircar (see [12] and the references therein). They developed a method that reduced the number of parameters that were needed for pricing and hedging to just three. Two were derived directly from the smile associated with options prices, and the other was some estimate of the underlying’s volatility. They introduced a type of model where the stochastic volatility factor followed a fast mean-reverting process. The data set of options they then considered was chosen so that the asymptotics would be valid. For example, the strikes of the options were near at-the-money and the options had expiration dates that were neither close nor far away. With this data set they found that their parameters were stable over time. Importantly, the option pricing equations they obtained within the fast mean-reverting regime were independent of the stochastic volatility driving factor. This was not the case in the entropy-based estimation procedure of Carmona and Xu [8]. A subsequent paper of Fouque et al. [14] extended the domain of applicability by introducing structure to the maturity cylces of options. We will not pursue that method here. In this paper, we introduce a method for calibrating volatility surfaces in a stochastic volatility 4
environment. We use an entropy-inspired framework that allows flexibility in volatilities when passing from the prior probability law to the pricing one. Because we are unable to observe the underlying stochastic volatility process, we adopt the framework of Fouque, Papanicolaou and Sircar and consider a regime where fast mean-reversion is apparent. We also relax the pricing constraints so that our volatility surface is less sensitive to unreliable out-of-the-money option prices. We formulate and solve a Hamilton-Jacobi-Bellman (HJB) equation associated with the variational problem. An HJB equation in a fast mean-rversion setting was studied by Sircar and Zariphopoulou [20] in the context of portfolio optimization. Our problem is different and results in a different HJB equation. From this equation, we obtain a volatility surface that consistently incorporates the Fouque-Papanicolaou-Sircar stochastic volatility parameters. These parameters are found to be stable over time. Our method, in fact, permits a natural interpolation between the model parameters of Fouque, Papanicolaou and Sircar and our own prior. It is an interpolation that depends on the level of pricing uncertainty. The rest of this paper is organized as follows. In Section 2 we formulate the problem and introduce a method for obtaining the volatility surface. We review the fast mean-reversion regime and simplify the corresponding HJB equation in Section 3. We reduce the complexity of the problem by studying fast mean-reversion asymptotics in Section 4. In Section 5 we describe the iterative procedure that was used to estimate the surface, while in Section 6 we detail the numerical methods we used to solve the resulting partial differential equations. In Section 7 we present the results of our method. We discuss the sensitivity of the method to the various parameters in Section 8. We end with a brief summary and conclusion.
2 2.1
Stochastic volatility in the relative entropy framework A simplified form for the relative entropy
In the relative entropy framework that we follow, we allow the functional form of the volatility to change from P0 , our prior probability law, to P , the probability law we wish to determine, because otherwise we would only have freedom to pick a risk premium. We therefore adapt the method of Avellaneda et. al. to the case of stochastic volatility. In order to do this, we discretize in time the stock and stochastic volatility processes so that we may determine the most singular part of the relative entropy. We assume the following dynamics under P0 : dSt = rdt + σ0 (St , Yt )dWt St dYt = γ0,t dt + κdZˆt where, as before, we assume the correlation between the two Brownian motions is ρ. In other words, dZˆt dWt = ρdt. 5
Here, r is the riskless interest rate process and σ0 is our prior estimate of the volatility surface. Yt is the stochastic volatility driving factor, which has a drift of γ0 and has its own volatility κ. Under P , we assume that the dynamics are given by: dSt = rdt + σ(St , Yt )dWt∗ St dYt = γt dt + κdZˆt∗
(1)
Here, γ and γ0 are functions that may depend on both t and Yt , and again we assume the same level of correlation between the Brownian motions. Details of the procedure used to get the relative entropy given the dynamics above are provided in Appendix A. The result is that the relative entropy generated per unit of time is just η(σ) = (σ 2 − σ02 )2
when σ 2 is not very far from σ02 . This particular form is an unsurprising consequence of approximating a distance function.
2.2
Formulation of the variational problem
Since we wish to minimize the relative entropy between our prior and a model that fits the data, we consider the optimization problem Z T P η(σs )ds sup −E σ
0
subject to the constraints:
P −rT E e i hi (ST ) − Ci ≤ βi , i
for each i and where βi ≥ 0. Here, hi (STi ) is the payoff of option i, which has expiration date Ti . Ci is the price of this option that is quoted in the market. The constants βi may be considered to be a measure of our confidence in the market data. We note that in previous works involving relative entropy minimization, these constants were taken to be zero (hard price constraints). By allowing deviation from market prices under our new measure P , we have softened the pricing constraints, since it is not practical for us to consider the market data as being exact. Some data points may be the result of only one trade and hence would not be reliable. Imposing hard constraints also leads to kinks in the smile near the data points. We simplify the objective above by forming a single Lagrangian. The following analysis is similar to the one in [9] for a different problem. Our objective may be rewritten as Z T inf E P e−rs η(σs )ds subject to σ 0 (λ+ ∀i : E P e−rTi h(STi ) − Ci ≤ βi i ) (λ− ∀i : Ci − E P e−rTi h(STi ) ≤ βi i ). 6
Here, Lagrange multipliers are written next to the constraints. From this, we obtain the augmented Lagrangian which we try to minimize over both sets of λ: X Z T X −rT − + P −rs P i (2) )(E e h(S ) − C ) + βi (λ+ − λ e η(σs )ds + (λ− sup −E i T i i + λi ). i i σ
0
i
i
− Without loss of generality, we may assume that, for each i, at most one of λ + i and λi is nonzero. This is because we are minimizing over λi , and so we could otherwise decrease both λi ’s by an equal but positive amount. We would then obtain a new expression where the third term in Equation (2) is decreased but the other two are unchanged.
Now suppose that λ+ i is nonzero (at the minimum). Then our solution is the same as the one for an optimization with the first constraint replaced by equality for that value of i. In this case, the second constraint cannot simultaneously be an equality and the solution must therefore be an interior one. Similarly, if the λ− i is nonzero, then the second constraint is an equality. These are the − only two possible cases. If we now write for each i, λi = λ+ i − λi , we may consider the objective: sup −E σ
P
Z
T
e
−rs
η(σs )ds +
0
X i
X λi (E P e−rTi h(STi ) − Ci ) + βi |λi |.
(3)
i
If we have to minimize this over λi (for each i), then there are two cases. First, the λi , for some i, may be positive, in which case we have solved the original optimization problem with the first inequality replaced by equality. Or, the λi , for that same i, is negative, in which case the original problem with the second constraint replaced by equality has been solved. Hence, the two problems given in Equations (2) and (3) are equivalent, as they both reduce to the same cases once minimization over λi has been performed. We may write an indirect value function Z T X i h λ P −r(s−t) V (t, s, y) = sup −Et e η(σs )ds + λi EtP e−r(Ti −t) h(STi ) . σ
t
i
We note that our value function discounts the entropy to make the subsequent expressions simpler. Our objective in Equation (3) amounts to minimizing the expression X X V λ (0, S, y) − λi C i + βi |λi | (4) i
i
over λ. We know that this V λ solves the HJB equation that is derived in Appendix B and has the form X 1 1 λi δ(t − Ti )hi (S). Vt − rV + rSVS + γVy + κ2 Vyy + Φ( S 2 VSS , ρκSVSy ) = − 2 2 i
V (T, S, y) = 0
7
Here, Φ(X1 , X2 ) = sup X1 σ 2 + X2 σ − η(σ). σ
Before describing how to obtain a volatility surface from these equations, it is worthwhile noting some features of the HJB equation. With ρ = 0, the HJB equation becomes much simpler because Φ(X, Y ) may be directly evaluated. The equation decouples and the value function would solve the same PDE that was found in Avellaneda et al. [5]. The volatility surface may be derived from the value function as follows. Once we solve for the λ = λ∗ that minimizes Equation (4), we may write the desired volatility surface as 1 λ∗ 2 λ∗ σ ∗ (S, y, t) = arg sup S 2 VSS σ + ρκSVSy . σ 2
(5)
We have thus described a way to incorporate stochastic volatility in a relative entropy framework. We did this by first determining an appropriate form for the relative entropy. Once we did this, however, encoding our objective into a value function is standard, and this value function must solve an HJB equation. Apart from the form of the relative entropy, our other main contribution so far is the incorporation of pricing uncertainty through the parameters β i . Although this methodology seems reasonable, there are some issues that need to be addressed. The HJB equation as written is highly dependent on many prior parameters. We need to know γ, κ and ρ accurately in order to proceed. Since these are all parameters associated with the unobservable stochastic volatility driving factor, determining them is a difficult task. The volatility surface, moreover, is dependent on the actual value of the stochastic volatility driving factor. The next section develops a way to get around these difficulties by considering a regime in which the stochastic volatility driving factor is fast mean-reverting. This represents the point of departure from the standard HJB equation methodology.
3 3.1
The fast mean-reversion setting A special form for the dynamics of the stochastic volatility driving process
We have so far dealt with the presence of stochastic volatility in a fairly general way. We now consider some special models for the dynamics of the volatility driving factor, Y . One feature that most models of stochastic volatility incorporate is mean reversion. This refers to the tendency of the process to go back to its invariant or long-run distribution. A particularly tractable model that exhibits mean reversion is the Ornstein-Uhlenbeck process dYt = α(m − Yt )dt + κdZˆt . 8
Here, α is the rate of mean reversion, while m is the long-run mean of the process. We suppose that α, κ and m are constants. Moreover, Zˆt is a Brownian motion that is correlated to the Brownian motion that drives the stock price process, Bt , with correlation ρ. The invariant distribution for this process is a Y0 that satisfies E[Lg(Y0 )] = 0 for any smooth and bound g, where L = α(m − y)
∂ 1 ∂2 + κ2 2 . ∂y 2 ∂y
If we let Ψ(y) be the density function of the invariant distribution, it is easily seen that (y − m)2 1 exp − Ψ(y) = √ , 2ν 2 2πν 2 where ν2 =
(6)
κ2 . 2α
This is exactly the normal density with mean m and variance ν 2 , which tells us that the parameter ν controls the size of the equilibrium fluctuations from the mean, m. The fast mean reversion limit was considered in Fouque, Papanicolaou and Sircar [12]. They found empirical evidence for fast mean reversion in options data. Analytically, fast mean reversion corresponds to the limit α → ∞. Care, however, must be taken to ensure that the invariant distribution remains the same as we take this limit. In other words, we would like our model to have fluctuations of the same size regardless of how quickly the volatility reverts to its mean. We therefore take √ κ = ν 2α for some constant ν. If this is the specification of Y under the physical measure, then to price derivatives we need to consider the process under the risk-neutral law. Suppose that under some physical probability law the stock and volatility driving process follow dSt = µSt dt + σ0 (St , Yt )dB0,t √ dYt = α(m − Yt )dt + ν 2αdZˆ0,t . Then, under a risk-neutral law, the processes must follow, by Girsanov’s theorem, ∗ dSt = rSt dt + σ0 (Yt )dB0,t p √ √ √ µ−r ∗ 2 − ν 2αγ0,t 1 − ρ dt + ν 2αdZˆ0,t dYt = α(m − Yt ) − ν 2αρ . σ0
9
(7)
Here, γt is a risk-premium factor or the market price of volatility risk that parametrizes the space of equivalent martingale measure. By taking, γ0,t = γ0 (t, St , Yt ), we specifically restrict ourselves to a Markovian setting. Equation (7) is the specification of the dynamics of the processes S and Y under our prior risk-neutral probability law. We now just need to find the relative entropy between another risk-neutral probability law and our prior. We assume that under P the price of the stock evolves as dSt = rSt dt + σ(St , Yt )dBt∗ p √ √ √ µ−r 2 dYt = α(m − Yt ) − ν 2αρ − ν 2αγt 1 − ρ dt + ν 2αdZˆt∗ . σ
(8)
With these particular assumptions on the dynamics of Yt , the form of η(σ) to first order in ρ is unchanged from the analysis in Appendix A. It is worth noting one important aspect of the dynamics under the prior probability law. We have assumed that σ0 (St , Yt ) = σ0 (Yt ). In other words, our prior estimate of the volatility is dependent only on the stochastic volatility driving factor and not on the stock price. With this choice of the functional form σ 0 now corresponds to function f in the work of Fouque, Papanicolaou and Sircar [12]. Although this particular form seems like a heavy restriction, it is not. The reason is that the fast mean-reversion setting does not need an explicit functional form for the prior volatility estimate. Instead, this setting replaces explicit functional forms with observable quantities and averages. This will be made clear later.
3.2
Simplifying the HJB equation
The HJB equation we need to solve is now a little different to the one we derived earlier. Because of the appearance of the σ in the denominator of one of the diffusive terms for Y t in Equation (8), Φ must now be a function of three variables. Specifically, we define Φ as: 1 2 Φ(X1 , X2 , X3 ) = sup σ X1 + ρX2 σ + ρX3 − η(σ) . (9) σ σ The HJB equation thus becomes: p √ Vt − rV + rSVS + (α(m − y) − ν 2αγ 1 − ρ2 )Vy + ν 2 αVyy X √ √ 1 2 +Φ λi δ(t − Ti )hi (S) S VSS , ν 2αSVSy , −ν 2α(µ − r)Vy = − 2
(10)
i
Before we carry out the asymptotic analysis we simplify Φ for σ close to σ 0 and for ρ small. We therefore expand σ in ρ (a small ρ expansion). Let: σ=σ ˜ + Aρ 10
where A is a coefficient to be determined. Here, σ ˜ is taken to be the maximizer when ρ = 0, which is just given by σ ˜ 2 = X1 /2 + σ02 . The first-order condition for a maximum gives: 2σX1 + ρX2 − ρ
X3 − 4σ(σ 2 − σ02 ) = 0. σ2
Upon substituting for σ, we find that the coefficient for ρ is 2AX1 + X2 −
X3 X3 − 4A(˜ σ 2 − σ02 ) − 8A˜ σ 2 = X2 − 2 − 8A˜ σ2. σ ˜2 σ ˜
Equating the coefficient of ρ to zero yields: X2 −
X3 − 8A˜ σ 2 = 0, σ ˜2
which implies that 3 X2 − X σ ˜2 . A= 4X1 + 8σ02
So, to first order in ρ: Φ(X1 , X2 , X3 ) =
X12 X3 + σ02 X1 + ρX2 σ ˜+ρ 4 σ ˜
which, upon expanding σ ˜ , leads to: X2 Φ(X1 , X2 , X3 ) = 1 + σ02 X1 + ρX2 4
X1 + σ02 2
1
2
+ ρX3
X1 + σ02 2
− 1 2
.
Given this form of Φ, we note that the maximizing value of σ is clearly approximated by σ ˜. We remark also that the small ρ expansion may be justified a posteriori. We will find that the y-dependence in our value function only occurs at the α1 scale. This implies that, since we find √ X2 ∝ αVxy , X2 scales like √1α , which is small for large α. A similar argument holds for X3 . Because ρ always premultiplies X2 or X3 , the small ρ expansion of Φ is reasonable. We now have a suitable form of the HJB equation, which uses the small ρ approximation for Φ. 2 p √ 1 1 1 2 2 2 Vt − rV + rSVS + (α(m − y) − ν 2αγ 1 − ρ )Vy + ν αVyy + S VSS + σ02 S 2 VSS 4 2 2 1 − 1 √ √ 2 2 1 2 1 2 + ρν 2αSVSy S VSS + σ02 − ρν 2α(µ − r)Vy S VSS + σ02 4 4 X =− λi δ(t − Ti )hi (S) i
V (0, S, y) = 0
(11)
11
4
Fast mean-reversion asymptotics applied to the HJB equation
In order to apply fast mean-reversion asymptotics, we suppose that α is large. For book-keeping purposes, we replace α by αε . We consider the value function corresponding to this situation in the limit as ε → 0. We first let x = log S (which means that 21 S 2 VSS = 12 (Vxx − Vx )). Substituting this into Equation (11) yields p √ 1 1 Vt − rV + (r − σ02 )Vx + σ02 Vxx + (α(m − y) − ν 2αγ 1 − ρ2 )Vy + ν 2 αVyy 2 2 1 − 1 √ √ 2 2 1 Vxx − Vx Vxx − Vx 2 2 2 + σ0 + σ0 + (Vxx − Vx ) + ρν 2αVxy − ρν 2α(µ − r)Vy 16 4 4 X =− λi hi (x)δ(t − Ti ). (12) i
Here, σ0 = σ0 (Yt ). We now replace α by αε so that we may apply some asymptotic techniques. We therefore write our equation as √ 1 1 (Vxx − Vx )2 1 + √ ρν 2αVxy (g(Vxx , Vx , y) − 1) L2 + √ L1 + L0 V + ε 16 ε ε X √ 1 1 − √ ρν 2αVy −1 − λi hi (x)δ(t − Ti ) (13) g(Vxx , Vx , y) ε i
where ∂2 ∂ 1 ∂ 1 − r + (r − σ0 (y)2 ) + σ0 (y)2 2 ∂t 2 ∂x 2 ∂x p √ √ √ ∂2 ∂ ∂ − ν 2αγ 1 − ρ2 − ρν 2α(µ − r) L1 = ρν 2α ∂x∂y ∂y ∂y 2 ∂ ∂ L0 = α(m − y) + ν 2α 2 ∂y ∂y 1 2 Vxx − Vx . g(Vxx , Vx , y) = + σ02 4 L2 =
V
(14)
−V
If we suppose that 0,xx4 0,x is small in relation to σ0 , which is reasonable as long as our prior is somewhat close to the actual σ, then we may approximate g as just g(y) = σ0 (y).
(15)
This approximation is consistent with the other approximations we have made in deriving a form of the relative entropy. We note that this approximation for g is also what we would have if we had happened to pick a prior that matched the prices within the tolerance desired. 12
We also expand V in powers of ε V = V0 +
4.1
√ √ εV1 + εV2 + ε εV3 + ....
The leading order term
In order to solve for the V0 term we equate the expansions of the HJB equation at the first few scales. At the 1ε scale L0 V0 = 0, which implies that V0 = V0 (t, x) (it is independent of y). At the
√1 ε
scale
√ √ L0 V1 + L1 V0 + ρν 2αV0,xy (g(V0,xx , V0,x , y) − 1) − ρν 2α(µ − r)V0,y (g(V0,xx , V0,x , y)) = 0 but since V0 is independent of y, the equation becomes L 0 V1 = 0 and this implies, as before, V1 = V1 (t, x). At the order 1 scale L 2 V0 + L 1 V1 + L 0 V2 +
X 1 (V0,xx − V0,x )2 = − λi hi (x)δ(t − Ti ), 16 i
where, again, we have used that V1 and V0 is independent of y to get rid of the Vxy part of the nonlinearity. Similarly, we may get rid of the L1 V1 = 0. Solvability for V2 requires that the other terms have average 0. But since V0 is independent of y, this means that the equation that determines solvability yields the following equation for V0 X 1 LBS (¯ σ0 )V0 + (V0,xx − V0,x )2 = − λi hi (x)δ(t − Ti ). (16) 16 i
Here, LBS is the Black-Scholes operator after the log transformation has been performed.The terminal condition is V0 (T, x) = 0. This is exactly the equation that is solved in Avellaneda et. al. [5]. It appears here as the zero order approximation in our method.
4.2
The correction term
Once we have solved the equation for V0 , we may remove the average (which we have set to zero) to get an equation for V2 2 1 ∂ V0 ∂V0 2 2 L0 V2 + (σ0 (y) − σ = 0, ¯ ) − 2 ∂x2 ∂x 13
which implies that 1 V2 (t, x, y) = − L−1 (σ0 (y)2 − σ ¯2) 2 0 We may therefore write V2 as 1 V2 (t, x, y) = − (φ(y) + c(t, x)) 2
∂ 2 V0 ∂V0 − ∂x2 ∂x
∂ 2 V0 ∂V0 − ∂x2 ∂x
.
,
(17)
where φ(y) is a solution of the Poisson equation L0 φ = σ0 (y)2 − < σ02 > . √ Moving on to the next scale ( ε), the relevant equation is √ 1 L2 V1 + L1 V2 + L0 V3 + (V0,xx − V0,x )(V1,xx − V1,x ) + ρν 2αV2,xy (g(V0,xx , V0,x , y) − 1) 8 √ − ρν 2αV2,y (1/g(V0,xx , V0,x , y) − 1) = 0. Solvability for V3 implies √ 1 ∂ 2 V2 LBS (¯ σ0 )V1 + (V0,xx − V0,x )(V1,xx − V1,x ) + ρν 2α g(V0,xx , V0,x , y) 8 ∂x∂y p √ √ ∂V2 1 ∂V2 − 2α 1 − ρ2 γ = 0. − ρν 2α(µ − r) g(V0,xx , V0,x , x) ∂y ∂y
(18)
The averages in the last three terms of Equation (18) may be further expanded into terms with derivatives in V0 , using Equation (17). Specifically, we obtain the following relations 3 D E 2 ∂ 2 V2 g(V0,xx , V0,x , y) ∂x∂y = 12 hgφ0 i ∂∂xV30 − ∂∂xV20 E D E 2 D ∂V0 ∂V2 ∂ V0 1 1 0 1 = (19) φ − 2 2 g ∂x g(V0,xx ,V0,xD ,x) ∂y E 2∂x 2 0 γ ∂V = 21 hγφ0 i ∂∂xV20 − ∂V ∂y ∂x So our equation for the correction V1 becomes
2 1 ∂ 3 V0 ˜2 ∂ V0 + A˜1 ∂V0 = 0, LBS (¯ σ0 )V1 + (V0,xx − V0,x )(V1,xx − V1,x ) + A˜3 + A (20) 8 ∂x3 ∂x2 ∂x which is a linear equation with a source term. In particular, the coefficients A˜1 , A˜2 and A˜3 are given by
√ 1 p √ 1 1 0 ˜ φ + ν α√ 1 − ρ2 γφ0 A1 = ρν α √ (µ − r) g 2 2
√ 1 0 √ 1 √ 1 p 1 0 ˜ A2 = −ρν α √ gφ − ρν α √ (µ − r) 1 − ρ2 γφ0 φ − ν α√ g 2 2 2
√ 1 A˜3 = ρν α √ gφ0 2
14
We note that, with this method, and as it stands now, we need the correlation ρ and the various parameters of the stationary distribution of the stochastic volatility driving factor.
4.3
Determination of the group parameters
With the approximation given in Equation (15), our parameters may be written as
√ √ 1 p √ 1 1 0 A1 = ε ρν α √ (µ − r) φ + ν α√ 1 − ρ2 γφ0 σ0 2 2
√ √ 1
√ 1 √ 1 p 1 0 0 0 A2 = − ε ρν α √ σ0 φ + ρν α √ (µ − r) φ + ν α√ 1 − ρ 2 σ0 φ σ0 2 2 2 A3 =
√ 1
√ ερν α √ σ0 φ0 . 2
(21)
These parameters, which we call our group parameters, are the same as the A˜i but they are now √ scaled by ε. We may determine these parameters as combinations of parameters that are observable from the smile. To do this, we outline the results in Fouque, Papanicolaou and Sircar [12]. In their work, they found that the price of an option, C may be expanded as ∂CBS ∂ 2 CBS ∂ 3 CBS C(t, x; T, K) = CBS (t, x; T, K) − (T − t) (2V3 − V2 ) , + (V2 − 3V3 ) + V3 ∂x ∂x2 ∂x3 (22) where CBS is the Black-Scholes option price and x = log(S). Here,
√ √ 1 √ 1 p √ 1 1 0 0 0 2 √ √ √ V2 = ε ν α 1 − ρ σ0 φ φ −ν α 2ρ σ0 φ − ρν α (µ − r) σ0 2 2 2 √ 1
√ (23) V3 = ερν α √ σ0 φ0 2 We deduce, then, that our approximate Ai ’s given by Equation (21 )are exactly A1 = 2V3 − V2
A2 = V2 − 3V3
A3 = V 3 .
The parameters V2 and V3 may be determined from the smile. By expanding the price of the option in powers of ε, we may write C(t, x; T, K; I) = CBS (t, x; T, K; σ0 ) + 15
√ ∂CBS εI1 , σ
(24)
where I1 is the correction to the implied volatility. By equating Equations (22) and (24), we find that the implied volatility I is just: I = σ0 +
3 V3 V2 V3 log(K/S) (r + σ02 ) − − 3 . 2 2 σ0 σ0 T − t σ0
The implied volatility is therefore an affine function of the log-moneyness-to-maturity-ratio LMMR =
log(K/S) , T −t
which means that we may fit the parameters aF P S and bF P S below to the smile curve determined by options prices I = aF P S LMMR + bF P S , where aF P S = − bF P S
V3 σ03
V3 = σ0 + 3 σ0
3 r + σ02 2
−
V2 . σ0
This means that the scaled parameters Ai which are required to determine the correction term √ εV1 are explicitly expressed as combinations of the parameters aF P S and bF P S . Moreover, given any such parameters a and b determined from a smile curve, the parameters A i are given by 3 A1 = −2aσ03 − σ0 ((σ0 − b) − a(r + σ02 )) 2 3 2 A2 = σ0 ((σ0 − b) − a(r + σ0 )) + 3aσ03 2 A3 = −aσ03
(25)
We have thus expressed the parameters needed for the correction term as simple expressions involving observable quantities. This is the power of the fast mean-reversion setting. Individual estimates of the parameters of the stochastic volatility models are not needed. Rather, group parameters like the Ai are all that are needed to not only price options as in Fouque, Papanicolaou and Sircar [12] but also to determine minimal entropy volatility surfaces as in this paper.
4.4
Summary of our analysis
What we have achieved now is a set of equations that incorporates the objectives we wanted 1. Uncertainty over the prior in the form of stochastic volatility that is independent of the specific form of the volatility since we only really need the group parameters, A 1 , A2 and A3 . 16
2. A value function that does not depend on the stochastic volatility driving factor because of the use of fast mean-reversion. 3. Uncertainty in the observed market prices The asymptotics helped us isolate the leading two terms of the value function that solved the HJB equation. These two terms were, indeed, independent of the stochastic volatility driving factor. The first term, moreover, solved the same equation that was solved in Avellaneda et al. [5].
5
Description of the estimation procedure
Our analysis permits a simplified approach to the estimation of the stochastic volatility surface. There are two possible approaches from here. We may estimate the parameters A i once or we may try to self-consistently determine the same parameters in an iteration scheme. In the first approach, we may determine A1 , A2 and A3 by following the smile-fitting procedure given in Section 4. We may then solve for V0 using Equation (16). With this V0 and the coefficients √ {Ai }, we may solve for εV1 for each value of λ using search P(20). Using some gradient P Equation √ algorithm, we may then minimize (over λ) V (0, x0 )− i λi Ci + i |λi |βi , where V = V0 + εV1 and x0 is the log of the current stock price. Finally, the appropriate derivative of V gives a volatility surface as described by Equation (5). We choose to iterate this scheme. We may obtain A1 , A2 and A3 from the prices as before, solve for V0 and V1 and obtain a volatility surface. But instead of stopping there, we price the same set of options in our data set again using the newly derived volatility surface and use these prices to get adjusted estimates for A1 , A2 and A3 . Repeating this procedure a few times then leads to a converged set of parameters. The fact that our procedure leads to a converged set of parameters tells us the group parameters are intrinsic to the volatility surface as a whole and that they may be self-consistently determined. They thus have the same role in volatility surface estimation as they do in option pricing. We call this method the volatility surface iteration procedure. It is used in our numerical calculations and is summarized succinctly as 1. Fit the implied volatility smile to an affine function of the log-moneyness-to-maturity ratio. This affine fit gives us the parameters aF P S and bF P S from which we obtain A1 , A2 and A3 . √ 2. Solve for V0 and εV1 . 3. Obtain a volatility surface. 4. Price options again using the volatility surface in Step 3. We do this using the Black-Scholes PDE with volatility σ. 17
5. Derive the implied volatility smile from the prices obtained in Step 4. Fit this smile to the log-moneyness-to-maturity ratio to obtain parameters a and b. Obtain new estimates of the parameters A1 , A2 and A3 from the fitted a and b. . 6. Stop if adjusted parameters differ from the previous iteration by less than δ, a small parameter (which we take to be of order 10−5 in the l2 norm). Otherwise, repeat the procedure from Step 2.
6 6.1
Numerical Methods The highest-order equation
Supposing that the stock pays dividend of d (constant), then the function V0 (x, t) satisfies the equation: 1 2 2 2 2 sup V0,t − rV0 + (r − d)V0,x + σ (V0,xx − V0,x ) − (σ − σ ¯0 ) 2 σ X =− λi hi (ex )δ(t − Ti ) i
Here, we choose σ ∈ [σmin , σmax ] in order to preserve the stability of our numerical scheme. Now, using the substitution V˜0 = e−rt V0 , yields: 1 ¯02 )2 V˜0,t + (r − d)V˜0,x + σ 2 (V˜0,xx − V˜0,x ) − e−rt (σ 2 − σ sup 2 σ∈[σmin ,σmax ] X =− λi hi (ex )e−rt δ(t − Ti ) i
We discretize the interior equation (spatial h, temporal ∆t), and supposing that we have the optimal σ at time t, we have the equation (after dropping the tildes): V0 (x, t) − V0 (x, t − ∆t) 1 2 V0 (x + h, t) + V0 (x − h, t) − 2V0 (x, t) + σ ∆t 2 h2 1 V0 (x + h, t) − V0 (x − h, t) + (r − d − σ 2 ) − e−rt (σ 2 − σ ¯02 )2 2 2h X =− e−rt λi hi (ex )δ(t − Ti ) i
In order to regularize the problem, we replace each option payoff e −rt hi (ex )δ(t − Ti ) by Pi = C(ex , ∆t, σ)1t≤Ti
