Closed-form approximations with respect to the mixing solution for option pricing under stochastic volatility
aa r X i v : . [ q -f i n . M F ] D ec CLOSED-FORM APPROXIMATIONSWITH RESPECT TO THE MIXING SOLUTIONFOR OPTION PRICING UNDER STOCHASTIC VOLATILITY
KAUSTAV DAS † AND NICOLAS LANGREN´E ‡ Abstract.
We consider closed-form approximations for European put option prices withinthe Heston and GARCH diffusion stochastic volatility models with time-dependent parameters.Our methodology involves writing the put option price as an expectation of a Black-Scholesformula and performing a second-order Taylor expansion around the mean of its argument. Thedifficulties then faced are simplifying a number of expectations induced by the Taylor expansion.Under the assumption of piecewise-constant parameters, we derive closed-form pricing formulasand devise a fast calibration scheme. Furthermore, we perform a numerical error and sensitivityanalysis to investigate the quality of our approximation and show that the errors are well withinthe acceptable range for application purposes. Lastly, we derive bounds on the remainder termgenerated by the Taylor expansion.Keywords: Stochastic volatility, Closed-form expansion, Closed-form approximation, Heston,GARCH Introduction
We consider the European put option pricing problem within the Heston and GARCH dif-fusion stochastic volatility models [19, 30, 36]. Our goal is to study how European put optionprices can be approximated in each of these frameworks via expansion of the so-called mixingsolution, which will be detailed later in this article. We find that through a second-order Taylorexpansion of the mixing solution, we can derive accurate approximations to European put op-tion prices in the aforementioned models. In addition, our methodology works naturally withtime-dependent parameters. This is seen as a major positive compared to other methodologies,for example transform methods, which cannot handle time-dependent parameters well. Fur-thermore, the approximation formulas are written solely in terms of iterated integrals, whichwe show obey convenient recursive properties when parameters are assumed to be piecewise-constant. This results in the approximation being essentially closed-form, and moreover leadsto a fast calibration scheme. Our method builds upon Drimus [14], in which the Heston modelis considered with constant parameters. A sensitivity analysis is performed in order to assessour approximation numerically. Furthermore, we present mathematical bounds on the errorterm in terms of higher order moments of the underlying variance process.It has been well established that volatility is highly dependent on the strike and maturity ofEuropean option contracts. This phenomenon is called the volatility smile or skew, an attributethe well known Black-Scholes model fails to take into account [7]. In response, there have beena number of frameworks proposed to explain the volatility smiles and skews observed in themarket. In particular, stochastic volatility models have been introduced, where the volatilityitself is a stochastic process possibly correlated with the spot. However with this added com-plexity, often option prices cannot be computed in a closed-form fashion. This is detrimental, † School of Mathematics, Monash University, Victoria, 3800 Australia. ‡ CSIRO Data61, Victoria, 3008 Australia.
E-mail addresses : [email protected], [email protected] .Research supported by an Australian Government Research Training Program (RTP) Scholarship. Generalised AutoRegressive Conditional Heteroskedasticity. For example, local volatility models, stochastic volatility models, and local-stochastic volatility models. as closed-form solutions lead to rapid option pricing, an attribute necessary for fast calibrationof financial models. Without a closed-form solution, option pricing must be done numericallyvia Monte Carlo or PDE methods, both methods being computationally costly.If we assume that the characteristic function of the log-spot is known explicitly, then the optionprice can be computed quasi-explicitly , albeit under the restrictive assumption of constant orpiecewise-constant parameters [9, 19, 29]. One class of models where this occurs is the class ofaffine models such as the Heston and Sch¨obel-Zhu model. However, for non-affine models, thecharacteristic function of the log-spot is rarely known explicitly, and such a procedure will notbe effective. Non-affine models, although usually intractable compared to their affine counter-parts, are often far more realistic. This has been shown in a number of studies, see for exampleChristoffersen et al. [11], Gander and Stephens [16], Kaeck and Alexander [22]. For these rea-sons, numerical procedures such as PDE and Monte Carlo methods have been substantiallydeveloped in the literature [3, 34].Closed-form approximations are an alternative methodology for option pricing, where the optionprice is approximated by a closed-form expression. The main advantages are that the optionprice can be computed rapidly and since transform methods are not used, time-dependent pa-rameters can usually be handled well. One motivation for quick option pricing formulas iscalibration, where the option price must be computed several times within an optimisation pro-cedure.There have been a plethora of results on closed-form expansions in the literature. For example,Lorig et al. [28] derive a general closed-form expression for the price of an option via a PDEapproach, as well as its corresponding implied volatility. Hagan et al. [18] use singular pertur-bation techniques to obtain an explicit approximation for the option price and implied volatilityin their SABR model. Al`os [2] show that from the mixing solution, one can approximate theput option price by decomposing it into a sum of two terms, one being completely correlationindependent and the other dependent on correlation. However, neither terms are explicit. Fur-thermore, similar to our work, Antonelli et al. [4], Antonelli and Scarlatti [5] show that underthe assumption of small correlation, an expansion can be performed with respect to the mixingsolution, where the resulting expectations can be computed using Malliavin calculus techniques.Similarly, in the case of the time-dependent Heston model, Benhamou et al. [6] consider themixing solution and expand around vol-of-vol, performing a combination of Taylor expansionsand computing the resulting terms via Malliavin calculus techniques.Stochastic volatility models usually either model the volatility directly, or indirectly via thevariance process. A critical assumption is that volatility or variance has some sort of meanreversion behaviour, and this is supported by empirical evidence, see for example Gatheral [17].Specifically, for modelling the variance, a large class of stochastic volatility models is given by d S t = S t (( r dt − r ft )d t + p V t d W t ) , S , d V t = κ t ( θ t V ˆ µt − V ˜ µt )d t + λ t V µt d B t , V = v , d h W, B i t = ρ t d t, Quasi-explicit meaning in terms of at most one-dimensional complex integrals, where the integrands are explicitfunctions. Affine stochastic volatility models are such that ln E (cid:16) e iu ln( S t ) (cid:17) is affine in ln S , i.e., the log of the characteristicfunction of the log-spot is an affine function in ln S [1, 15]. Our model formulation here is for FX market purposes, but can be easily adjusted for equity and fixed incomemarkets. whereas for modelling the volatility, this class is of the formd S t = S t (( r dt − r ft )d t + V t d W t ) , S , d V t = κ t ( θ t V ˆ µt − V ˜ µt )d t + λ t V µt d B t , V = v , d h W, B i t = ρ t d t, for some ˜ µ, ˆ µ and µ ∈ R . In this article, we will focus on this class of stochastic volatilitymodels. Some popular models in the literature include:Model Variance/Volatility Dynamics of V ˆ µ ˜ µ µ Heston [19] Variance d V t = κ t ( θ t − V t )d t + λ t √ V t d B t V t = κ t ( θ t − V t )d t + λ t d B t V t = κ t ( θ t − V t )d t + λ t V t d B t V t = κ t ( θ t − V t )d t + λ t V t d B t V t = κ t ( θ t V t − V t )d t + λ t V / t d B t [10, 26] Volatility d V t = κ t ( θ t V t − V t )d t + λ t V t d B t • Section 2 details preliminary calculations, where we express the put option price as the mixingsolution. Once done, a second-order Taylor expansion is performed, expressing the approxi-mation formula in terms of a number of expectations of functionals of the underlying varianceprocess. • Section 3 details how to derive more convenient expressions for these expectations obtainedin Section 2 through change of measure techniques. Specifically, we rewrite the spot S T as aconvenient expression so as to construct a term which is a Dol´eans-Dade exponential, therebydefining a Radon-Nikodym derivative. This term allows us to change measure, allowing formore convenient calculations. • Section 4 introduces specific models. As precise dynamics are assumed, the objective is toderive expressions in terms of specific iterated integrals for the expectations from Section 3.In particular, we consider the Heston and GARCH diffusion models. • In Section 5, under the assumption of piecewise-constant parameters, we obtain closed-formexpressions for our pricing formulas. Moreover, we design a fast calibration scheme. Specif-ically, we show that the iterated integrals in which the pricing formulas in Section 4 are interms of satisfy some convenient recursive properties. There exist other classes of stochastic volatility models, for example the exponential Ornstein-Uhlenbeck modelWiggins [35] is not a part of either of these classes. Also known as the Logistic or XGBM model. • Section 6 is dedicated to a numerical error and sensitivity analysis for the Heston and GARCHdiffusion models. • In Section 7 we perform an error analysis, bounding the error in the expansion in terms ofhigher order moments of the underlying variance process.2.
Preliminaries
Suppose the spot S with variance σ follows the dynamicsd S t = S t (( r dt − r ft )d t + √ σ t d W t ) , S , d σ t = α ( t, σ t )d t + β ( t, σ t )d B t , σ , d h W, B i t = ρ t d t, where W and B are Brownian motions with deterministic, time-dependent instantaneous corre-lation ( ρ t ) ≤ t ≤ T , defined on the filtered probability space (Ω , F , ( F t ) ≤ t ≤ T , Q ). Here T is a finitetime horizon, where ( r dt ) ≤ t ≤ T and ( r ft ) ≤ t ≤ T are the deterministic, time-dependent domesticand foreign interest rates respectively. In addition, ρ t ∈ [ − ,
1] for any t ∈ [0 , T ]. Furthermore,( F t ) ≤ t ≤ T is the filtration generated by ( W, B ) which satisfies the usual assumptions. In thefollowing, E ( · ) denotes the expectation under Q , where Q is a domestic risk-neutral measurewhich we assume to be chosen. We assume that the drift and and diffusion coefficients of σ aresuch that σ has a pathwise unique strong solution. Remark 2.1 (Geometric Brownian motion process) . Let ˜ B be an arbitrary Brownian motionprocess. We call Y a GBM( y ; µ t , ν t ) if it solves the SDEd Y t = µ t Y t d t + ν t Y t d ˜ B t , Y = y , where µ and ν are adapted to the Brownian filtration and satisfy some regularity conditions(for example, µ and ν bounded on [0 , T ] is sufficient).We decompose the Brownian motion W as W t = R t ρ u d B u + R t p − ρ u d Z u , where Z isa Brownian motion under Q such that B and Z are independent. Then, noticing S is aGBM( S ; r dt − r ft , √ σ t ), we obtain the expression S T = S ξ T exp (cid:26)Z T ( r dt − r ft )d t − Z T σ t (1 − ρ t )d t + Z T q σ t (1 − ρ t )d Z t (cid:27) , where ξ t := exp (cid:26)Z t ρ u √ σ u d B u − Z t ρ u σ u d u (cid:27) . Assumption A. (i) sup t ∈ [0 ,T ] lim sup x →∞ ρ t β ( t,x ) √ x + α ( t,x ) x < ∞ . (ii) β ( t, x ) ≤ K (1 + x ) . We will enforce Assumption A for the rest of this article.
Lemma 2.1. ξ is a martingale. Proof.
The proof is inspired by arguments made in Theorem 2.4 (i) of Lions and Musiela [27].First notice that d ξ t = ρ t √ σ t ξ t d B t , ξ = 1 . (2.1) Meaning that ( F t ) ≤ t ≤ T is right continuous and augmented by Q null-sets. Let τ n := inf { t ≥ σ t > n } be the first time σ crosses the level n . Then τ n ↑ ∞ Q a.s.. Define ξ nt := ξ t ∧ τ n . Then by solving eq. (2.1) up to τ n , this yields the pathwise unique strong solution ξ nt = exp (cid:18)Z t ρ u √ σ u { u ≤ τ n } d B u − Z t ρ u σ u { u ≤ τ n } d u (cid:19) . Since the integrand of the above Itˆo integral is bounded, ( ξ nt ) t is a martingale for each n ∈ N .Furthermore, ξ is a non-negative local-martingale. Hence by Fatou’s lemma, we have that ξ isa non-negative supermartingale. Our goal now is to show thatsup n E ( ξ nT log( ξ nT )) < ∞ since then, by the Vall´ee Poussin theorem, ( ξ nT ) n is uniformly integrable and thus E ( ξ T ) = 1.Combined with the fact that ξ is a non-negative supermartingale, this will then imply that ξ isa martingale. Now, ξ nT log( ξ nT ) = ξ nT (cid:18)Z T ρ t √ σ t { t ≤ τ n } d B t − Z T ρ t σ t { t ≤ τ n } d t (cid:19) . Clearly ξ nT is a Radon-Nikodym derivative which defines a measure ˆ Q . Hence by Girsanov’stheorem, ˆ B t = B t − R t ρ u √ σ u { u ≤ τ n } d u is a Brownian motion under ˆ Q . Denote the expectationunder ˆ Q by ˆ E . Then E ( ξ nT log( ξ nT )) = E (cid:18) ξ nT (cid:20)Z T ρ t √ σ t { t ≤ τ n } d B t − Z T ρ t σ t { t ≤ τ n } d t (cid:21)(cid:19) = ˆ E (cid:18)Z T ρ t √ σ t { t ≤ τ n } d ˆ B t + 12 Z T ρ t σ t { t ≤ τ n } d t (cid:19) = 12 Z T ˆ E ( ρ t σ t { t ≤ τ n } )d t ≤ Z T ˆ E ( σ t )d t. So it suffices to determine a bound on ˆ E ( σ t ) for t ≤ T . Under ˆ Q , we haved σ t = h ρ t β ( t, σ t ) √ σ t { t ≤ τ n } + α ( t, σ t ) i d t + β ( t, σ t )d ˆ B t . Now as a consequence of item (i) in Assumption A, this implies that there exists some constant
M > ρ t β ( t, x ) √ x + α ( t, x ) ≤ M (1 + x ) , (2.2)for all x ≥
0, uniformly in t ≤ T . Utilising eq. (2.2) yieldsd σ t ≤ M (1 + σ t )d t + β ( t, σ t )d ˆ B t . Hence ˆ E ( σ t ) ≤ ˆ E (cid:18) σ + M Z t (1 + σ u )d u + Z t β ( u, σ u )d ˆ B u (cid:19) ≤ σ + 4 M ˆ E (cid:18)Z t (1 + σ u )d u (cid:19) + 2ˆ E (cid:18)Z t β ( u, σ u )d ˆ B u (cid:19) ≤ σ + 4 M Z t ˆ E (1 + σ u ) d u + 2 Z t ˆ E (cid:0) β ( u, σ u ) (cid:1) d u ≤ σ + 4 M Z t ˆ E h (2 + 2 σ u ) + K M (1 + σ u ) i d u, where we have used that ( a + b ) ≤ a + 2 b , (cid:16)R t f ( u )d u (cid:17) ≤ R t f ( u )d u (Jensen’s inequality)and item (ii) in Assumption A. Define m t := ˆ E ( σ t ). After redefining constants, we obtain thefollowing integral inequality, m t ≤ c (1 + t ) + c Z t m u d u. Then by Gronwall’s inequality, we obtain m t ≤ c (1 + t ) e ct . (cid:3) Let Put := e − R T r dt d t E ( K − S T ) + be the price of a put option on S . Denote by ( F Bt ) ≤ t ≤ T the filtration generated by the Brownianmotion B , as well as N ( · ) and φ ( · ) the standard Normal distribution and density functionsrespectively. Proposition 2.1.
Put can be expressed asPut = E n e − R T r dt d t E (cid:2) ( K − S T ) + |F BT (cid:3)o = E (cid:18) Put BS (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19)(cid:19) , (2.3)where Put BS ( x, y ) := Ke − R T r dt d t N ( − d − ) − xe − R T r ft d t N ( − d + ) ,d ± ( x, y ) := d ± := ln( x/K ) + R T ( r dt − r ft )d t √ y ± √ y. (2.4) Proof.
This is a consequence of the mixing solution methodology, which is detailed in Appen-dix A. (cid:3)
Proposition 2.2 (Second-order put option price approximation) . The second-order put optionprice approximation, denoted by Put (2) , is given byPut (2) = Put BS (ˆ x, ˆ y )+ 12 ∂ xx Put BS (ˆ x, ˆ y ) S E ( ξ T − + 12 ∂ yy Put BS (ˆ x, ˆ y ) E (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + ∂ xy Put BS (ˆ x, ˆ y ) S E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19)(cid:27) , (2.5)where (ˆ x, ˆ y ) := ( S , R T (1 − ρ t ) E ( σ t )d t ). Proof.
Notice that Put BS is a composition of smooth functions. Hence, Put BS is smooth on( R ; R + ). We Taylor expand Put BS around the mean of (cid:16) S ξ T , R T σ t (1 − ρ t )d t (cid:17) up to second-order. By Lemma 2.1, the expansion point is (ˆ x, ˆ y ) = ( S , R T (1 − ρ t ) E ( σ t )d t ). ThusPut BS (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19) ≈ Put BS (ˆ x, ˆ y )+ ∂ x Put BS (ˆ x, ˆ y ) S ( ξ T −
1) + ∂ y Put BS (ˆ x, ˆ y ) (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + 12 ∂ xx Put BS (ˆ x, ˆ y ) S ( ξ T − + 12 ∂ yy Put BS (ˆ x, ˆ y ) (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + ∂ xy Put BS (ˆ x, ˆ y ) S ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) . Taking expectation yields a second-order approximation to the put option price, that is, Put (2) .Notice that Put BS (ˆ x, ˆ y ) is a deterministic quantity, thus the first-order terms will vanish. (cid:3) The partial derivatives of Put BS are analogous to the second, third and fourth order Black-Scholes Greeks, which are explicit and provided in Appendix B. In order to simplify Put (2) ,what remains to be done is the calculation of each of the expectations from Proposition 2.2,which are E ( ξ T − , (2.6) E (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) , (2.7) E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19)(cid:27) . (2.8) Remark 2.2 (Call options) . When considering the pricing of call options, the approximationmethodology is essentially the same as that of put options. It is clear that the mixing solutionmethodology from Appendix A can be easily adapted to the case of call options, in which casewe obtain Call := e − R T r dt d t E ( S T − K ) + = E n e − R T r dt d t E (cid:2) ( S T − K ) + |F BT (cid:3)o = E (cid:18) Call BS (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19)(cid:19) , where Call BS ( x, y ) := xe − R T r ft d t N ( d + ) − Ke − R T r dt d t N ( d − ) . Moreover, Put-Call parity states thatCall BS ( x, y ) − Put BS ( x, y ) = xe − R T r ft d t − Ke − R T r dt d t . Thus the second-order Call partial derivatives will be the same as the second-order Put partialderivatives. Comparing the mixing solution expressions for Call and Put, one can see that theonly difference between the call and put option second-order approximation is the zero-th orderterm. For this reason, we need only consider the pricing of put options.
Remark 2.3 (Greeks) . One can obtain second-order approximation of put option Greeks bysimply differentiating the expression for Put (2) . These expressions are provided in Appendix C. Calculation of expectations
Lemma 3.1.
Let ( Q n ) n ≥ be a sequence of probability measures equivalent to Q , defined bythe Radon-Nikodym derivativesd Q n +1 d Q n := ξ ( n ) T := exp (cid:26)Z T ρ u √ σ u d B nu − Z T ρ u σ u d u (cid:27) , ξ (0) T := ξ T , n ≥ , where Q := Q and B := B . Under Q n , B nt := B n − t − R t ρ u √ σ u d u is a Brownian motion. Proof.
This is a clear consequence of Girsanov’s theorem. (cid:3)
Remark 3.1.
Through Lemma 3.1, we have the following relationships: ξ ( n ) T = ξ ( n − T e − R T ρ u σ u d u , n ≥ , (3.1)and E Q n ( X ) = E Q n − ( Xξ ( n − T ) , E Q n − ( X ) = E Q n X ξ ( n − T ! . (3.2)These relationships allow for alternative and often more convenient calculations of expectationsunder Q .Using Remark 3.1, we can now give alternative expressions for the expectations seen in eqs. (2.6)to (2.8).3.1. E ( ξ T − . First, expanding eq. (2.6) gives E ( ξ T − = E ( ξ T ) − . This second moment can be dealt with a number of changes of measures. E ( ξ T ) = E Q ( ξ T ) = E Q ( ξ (1) T e R T ρ t σ t d t )= E Q ( e R T ρ t σ t d t ) . (3.3)Under the assumption of constant parameters we may calculate eq. (3.3) explicitly via theLaplace transform for certain processes σ . However to our knowledge, there exists no explicitsolution when parameters are time-dependent, see Hurd and Kuznetsov [21]. Instead, we ap-proximate eq. (3.3) by expanding the exponential around the mean of R T ρ t σ t d t to second-order. E Q ( e R T ρ t σ t d t ) ≈ E Q ( e R T ρ t E Q ( σ t )d t " Z T ρ t ( σ t − E Q ( σ t )) d t + 12 (cid:18)Z T ρ t ( σ t − E Q ( σ t )d t ) (cid:19) = e R T ρ t E Q ( σ t )d t ( E Q "(cid:18)Z T ρ t ( σ t − E Q ( σ t )d t ) (cid:19) = e R T ρ t E Q ( σ t )d t (cid:26) Z T ρ t Z t ρ s Cov Q ( σ s , σ t )d s d t (cid:27) , where we have used the fact that (cid:16)R T f ( t )d t (cid:17) = 2 R T f ( t ) (cid:16)R t f ( s )d s (cid:17) d t .3.2. E (cid:16)R T (1 − ρ t )( σ t − E ( σ t ))d t (cid:17) . To calculate eq. (2.7), we use the same approach fromSection 3.1. E (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) = 2 Z T (1 − ρ t ) (cid:18)Z t (1 − ρ s )Cov( σ s , σ t )d s (cid:19) d t. That is, α ( t, σ t ) = α ( σ t ), β ( t, σ t ) = β ( σ t ) and ρ t = ρ . E n ( ξ T − (cid:16)R T (1 − ρ t )( σ t − E ( σ t ))d t (cid:17)o . Calculation of the mixed expectation eq. (2.8)gives E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19)(cid:27) = Z T (1 − ρ t ) (cid:0) E ( ξ T σ t ) − E ( σ t ) (cid:1) d t = Z T (1 − ρ t ) (cid:0) E Q ( σ t ) − E ( σ t ) (cid:1) d t. Pricing for specific models
We now introduce specific dynamics for both the spot and its underlying variance process. FromSection 3, it is apparent that simplifying the expression for Put (2) will largely depend on thetractability of the variance process σ under the original measure Q , as well as the artificialmeasures Q and Q . Remark 4.1 (Notation for key processes) . Let (˜ κ t ) ≤ t ≤ T , (˜ θ t ) ≤ t ≤ T and (˜ λ t ) ≤ t ≤ T be time-dependent, deterministic, strictly positive and bounded on [0 , T ]. Let ˜ B be an arbitrary Brow-nian motion. We will utilise the following terminology:(1) If ˜ V solves d ˜ V t = ˜ κ t (˜ θ t − ˜ V t )d t + ˜ λ t q ˜ V t d ˜ B t , ˜ V = ˜ v , then we call ˜ V a CIR(˜ v ; ˜ κ t , ˜ θ t , ˜ λ t ).(2) If ˜ V solves d ˜ V t = ˜ κ t (˜ θ t − ˜ V t )d t + ˜ λ t ˜ V t d ˜ B t , ˜ V = ˜ v , then we call ˜ V an IGa(˜ v ; ˜ κ t , ˜ θ t , ˜ λ t ).We point the reader towards Appendix D for further information on these processes.4.1. Heston model.
Suppose the spot S with variance V follows the Heston dynamicsd S t = S t (( r dt − r ft )d t + p V t d W t ) , S , d V t = κ t ( θ t − V t )d t + λ t p V t d B t , V = v , d h W, B i t = ρ t d t, (4.1)where ( κ t ) ≤ t ≤ T , ( θ t ) ≤ t ≤ T and ( λ t ) ≤ t ≤ T are time-dependent, deterministic, strictly positiveand bounded on [0 , T ]. It is clear that the variance process V in eq. (4.1) is a CIR( v ; κ t , θ t , λ t ). Remark 4.2.
Note that the Heston model satisfies Assumption A, sincesup t ∈ [0 ,T ] lim sup x →∞ ρ t λ t √ x √ x + κ t ( θ t − x ) x = sup t ∈ [0 ,T ] lim sup x →∞ (cid:18) ρ t λ t − κ t + κ t θ t x (cid:19) < ∞ . Hence by Lemma 2.1, ξ is a martingale. Lemma 4.1.
Let ( Q n ) n ≥ be a sequence of probability measures equivalent to Q , defined bythe Radon-Nikodym derivativesd Q n +1 d Q n := ξ ( n ) T := exp (cid:26)Z T ρ u p V u d B nu − Z T ρ u V u d u (cid:27) , ξ (0) T := ξ T , n ≥ , where Q := Q and B := B . Under Q n , B nt := B n − t − R t ρ u √ V u d u is a Brownian motion. For n ≥
0, the dynamics of V under the measure Q n ared V t = ( κ t − nλ t ρ t ) (cid:18) θ t κ t κ t − nλ t ρ t − V t (cid:19) d t + λ t p V t d B nt , which is a CIR( v ; κ t − nλ t ρ t , θ t κ t κ t − nλ t ρ t , λ t ). Proof.
This lemma is simply obtained through Lemma 3.1, then expressing V under the newmeasures Q n . (cid:3) Thus the variance process V is a CIR process under all measures considered.4.1.1. Pricing under the Heston framework.
The second-order approximation of the put optionprice in the Heston framework is given by the following theorem.
Theorem 4.1 (Second-order Heston put option price) . Let ˆ x = S and ˆ y = R T (1 − ρ t ) E ( V t )d t = R T (1 − ρ t ) n v e − R t κ z d z + R t e − R tu κ z d z κ u θ u d u o d t . The second-order approximation to the putoption price in the Heston model, denoted by Put (2)H , isPut (2)H = Put BS (ˆ x, ˆ y )+ 12 ∂ xx Put BS (ˆ x, ˆ y ) S E ( ξ T − + 12 ∂ yy Put BS (ˆ x, ˆ y ) E (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19) + ∂ xy Put BS (ˆ x, ˆ y ) S E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19)(cid:27) , (4.2)where E ( ξ T − ≈ e R T ρ t E Q ( V t )d t (cid:26) Z T ρ t Z t ρ s Cov Q ( V s , V t )d s d t (cid:27) − , E Q ( V t ) = v e − R t κ z − λ z ρ z d z + Z t e − R tu κ z − λ z ρ z d z κ u θ u d u, Cov Q ( V s , V t ) = e − R ts κ z − λ z ρ z d z · Z s λ u e − R su κ z − λ z ρ z d z (cid:20) v e − R u κ z − λ z ρ z d z + Z u e − R up κ z − λ z ρ z d z κ p θ p d p (cid:21) d u, E (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19) = 2 Z T (1 − ρ t ) Z t (1 − ρ s ) · (cid:20) e − R ts κ z d z Z s λ u e − R su κ z d z (cid:26) v e − R u κ z d z + Z u e − R up κ z d z κ p θ p d p (cid:27) d u (cid:21) d s ! d t, and E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19)(cid:27) = Z T (1 − ρ t ) ( v (cid:16) e − R t κ z − λ z ρ z d z − e − R t κ z d z (cid:17) + Z t (cid:16) e − R tu κ z − λ z ρ z d z − e − R tu κ z d z (cid:17) κ u θ u d u ) d t. Proof.
Use Proposition 2.2 and adapt Section 3 to the Heston framework. The moments of V correspond to the moments of a CIR process, which are obtained from Appendix D.1. (cid:3) GARCH diffusion model.
Suppose the spot S with variance V follows the GARCHdiffusion dynamics d S t = S t (( r dt − r ft )d t + p V t d W t ) , S , d V t = κ t ( θ t − V t )d t + λ t V t d B t , V = v , d h W, B i t = ρ t d t, (4.3)where ( κ t ) ≤ t ≤ T , ( θ t ) ≤ t ≤ T and ( λ t ) ≤ t ≤ T are time-dependent, deterministic, strictly positiveand bounded on [0 , T ]. It is evident that the variance process V in eq. (4.3) is an IGa( v ; κ t , θ t , λ t ). Remark 4.3.
Unlike the Heston model, the GARCH diffusion model does not satisfy Assump-tion A, sincesup t ∈ [0 ,T ] lim sup x →∞ ρ t λ t x √ x + κ t ( θ t − x ) x = sup t ∈ [0 ,T ] lim sup x →∞ (cid:18) ρ t λ t √ x − κ t + κ t θ t x (cid:19) = ∞ . Thus, there is no guarantee that ξ is indeed a martingale. Instead, we will now assume that ξ is a martingale and show that, even if it is indeed a martingale, it does not help us. Lemma 4.2.
Let ( Q n ) n ≥ be a sequence of probability measures equivalent to Q , defined bythe Radon-Nikodym derivativesd Q n +1 d Q n := ξ ( n ) T := exp (cid:26)Z T ρ u p V u d B nu − Z T ρ u V u d u (cid:27) , ξ (0) T := ξ T , n ≥ , where Q := Q and B := B . Under Q n , B nt := B n − t − R t ρ u √ V u d u is a Brownian motion. For n ≥
0, the dynamics of V under the measure Q n ared V t = κ t (cid:18) θ t − V t + nλ t ρ t κ t V / t (cid:19) d t + λ t V t d B nt . Proof.
This lemma is simply obtained through Lemma 3.1, then expressing V under the newmeasures Q n . (cid:3) Remark 4.4.
Let Q n be defined as in Lemma 4.2. Under the measures Q n , n ≥ V has noknown expression for its solution, nor known expression for its moments. Validity of Remark 4.4.
The SDE in Lemma 4.2 is a linear diffusion type SDE. From Appen-dix E, it is known that if an explicit solution exists, it is given by V t = Y t /F t , where F is a GBM(1; λ t , − λ t ) and Y is the solution to the integral equation (written in differ-ential form) d Y t = (cid:18) κ t θ t F t − κ t Y t + nλ t ρ t κ t Y / t F − / t (cid:19) d t. Define A t := κ t θ t F t and C t := nλ t ρ t κ t F − / t . Note that A t and C t are both non-differentiable in t . Thus d Y t = (cid:16) A t − κ t Y t + C t Y / t (cid:17) d t. As far as we know, there is no explicit solution to these types of integral equations in theliterature, even when A and C are constants. As for explicit moments, it is unclear how toapproach this problem. There seems to be no approach to this problem in the literature,especially in the case of time-dependent parameters, see for example Kloeden and Platen [23]Chapter 4.4 for a comprehensive list of explicitly solvable SDEs. Furthermore, as an explicitsolution does not exist, we cannot use the method of approximating moments via the SDE’ssolution. (cid:3) Pricing under the GARCH diffusion framework: ρ = 0 . In the GARCH diffusion model,Assumption A is violated, and thus there is no guarantee ξ is a martingale. Additionally,assuming ξ is a martingale is not helpful, as the change of measure technique gives an intractabledynamic for V ; we cannot appeal to it for calculating expectations. However, in the case of ρ = 0 a.e., this implies ξ T = 1 Q a.s., and one will notice that the terms in the expansionrequiring a change of measure will disappear. Of course, the cost is the additional restrictiveassumption that spot and volatility movements are uncorrelated. We hope to mitigate thisissue in future work by combining this approach with small correlation expansion methods, seeAntonelli et al. [4], Antonelli and Scarlatti [5]. Theorem 4.2 (Second-order GARCH put option price) . Assume ρ = 0 a.e.. Let ˆ x = S and ˆ y = R T E ( V t )d t . Then the second-order put option price in the GARCH diffusion model,denoted by Put (2)GARCH , isPut (2)GARCH = Put BS (ˆ x, ˆ y ) + ∂ yy Put BS (ˆ x, ˆ y ) Z T (cid:18)Z t Cov( V s , V t )d s (cid:19) d t. (4.4)Both Cov( V s , V t ) and E ( V t ) are given in Appendix D.2. Proof.
Use Proposition 2.2 under the assumption of ρ = 0 a.e.. Then notice that E (cid:18)Z T ( V t − E ( V t ))d t (cid:19) = 2 Z T (cid:18)Z t Cov( V s , V t )d s (cid:19) d t. (cid:3) Explicit formulas and fast calibration
In this section, under the assumption of piecewise-constant parameters, we provide closed-formformulas for the price of a put option as well as a fast calibration scheme. To do this, weappeal to the results from Section 4 on the approximation of prices of put options in the variousmodels. We recognise that these prices are in terms of specific iterated integrals. Our goal isto show that these iterated integrals obey a convenient recursive property when parameters arepiecewise-constant. To this end, define the integral operator ω ( k,l ) T := Z T l u e R u k z d z d u. (5.1)In addition, we define the n -fold integral operator using the following recurrence: ω ( k ( n ) ,l ( n ) ) , ( k ( n − ,l ( n − ) ,..., ( k (1) ,l (1) ) T := ω (cid:0) k ( n ) ,l ( n ) w ( k ( n − ,l ( n − , ··· , ( k (1) ,l (1)) · (cid:1) T , n ∈ N . (5.2) The rest of the section is dedicated to making these integral operators explicit when parametersare piecewise-constant.Let T = { T , T , . . . , T N − , T N = T } , where T i < T i +1 be a collection of maturity dates on[0 , T ], with ∆ T i := T i +1 − T i and ∆ T ≡
1. When the dummy functions are piecewise-constant,that is, l ( n ) t = l ( n ) i on t ∈ [ T i , T i +1 ) and similarly for k ( n ) , then we can recursively calculate theintegral operators eq. (5.1) and eq. (5.2). Define e ( k ( n ) ,...,k (1) ) t := e R t P nj =1 k ( j ) z d z ,ϕ ( k,p ) T i ,t := Z tT i γ pi ( u ) e R uTi k z d z d u, where γ i ( u ) := ( u − T i ) / ∆ T i and p ∈ N ∪ { } . In addition, define n -fold extension of ϕ ( · , · ) T i , · recursively by ϕ ( k ( n ) ,p n ) ,..., ( k (2) ,p ) , ( k (1) ,p ) T i ,t := Z tT i γ p n i ( u ) e R uTi k ( n ) z d z ϕ ( k ( n − ,p n − ) ,..., ( k (2) ,p ) , ( k (1) ,p ) T i ,u d u, For example ω ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T = Z T l (3) u e R u k (3) z d z (cid:18)Z u l (2) u e R u k (2) z d z (cid:18)Z u l (1) u e R u k (1) z d z d u (cid:19) d u (cid:19) d u . where p n ∈ N ∪ { } . With the assumption that the dummy functions are piecewise-constant,we can obtain the integral operator at time T i +1 expressed by terms at time T i . ω ( k (1) ,l (1) ) T i +1 = ω ( k (1) ,l (1) ) T i + l (1) i e ( k (1) ) T i ϕ ( k (1) , T i ,T i +1 ,ω ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i +1 = ω ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (2) i e ( k (2) ) T i ϕ ( k (2) , T i ,T i +1 ω ( k (1) ,l (1) ) T i + l (2) i l (1) i e ( k (2) ,k (1) ) T i ϕ ( k (2) , , ( k (1) , T i ,T i +1 ,ω ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i +1 = ω ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (3) i e ( k (3) ) T i ϕ ( k (3) , T i ,T i +1 ω ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (3) i l (2) i e ( k (3) ,k (2) ) T i ϕ ( k (3) , , ( k (2) , T i ,T i +1 ω ( k (1) ,l (1) ) T i + l (3) i l (2) i l (1) i e ( k (3) ,k (2) ,k (1) ) T i ϕ ( k (3) , , ( k (2) , , ( k (1) , T i ,T i +1 ,ω ( k (4) ,l (4) ) ,..., ( k (1) ,l (1) ) T i +1 = ω ( k (4) ,l (4) ) , ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (4) i e ( k (4) ) T i ϕ ( k (4) , T i ,T i +1 ω ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (4) i l (3) i e ( k (4) ,k (3) ) T i ϕ ( k (4) , , ( k (3) , T i ,T i +1 ω ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (4) i l (3) i l (2) i e ( k (4) ,k (3) ,k (2) ) T i ϕ ( k (4) , , ( k (3) , , ( k (2) , T i ,T i +1 ω ( k (1) ,l (1) ) T i + l (4) i l (3) i l (2) i l (1) i e ( k (4) ,k (3) ,k (2) ,k (1) ) T i ϕ ( k (4) , , ( k (3) , , ( k (2) , , ( k (1) , T i ,T i +1 ,ω ( k (5) ,l (5) ) ,..., ( k (1) ,l (1) ) T i +1 = ω ( k (5) ,l (5) ) , ( k (4) ,l (4) ) , ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (5) i e ( k (5) ) T i ϕ ( k (5) , T i ,T i +1 ω ( k (4) ,l (4) ) , ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (5) i l (4) i e ( k (5) ,k (4) ) T i ϕ ( k (5) , , ( k (4) , T i ,T i +1 ω ( k (3) ,l (3) ) , ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (5) i l (4) i l (3) i e ( k (5) ,k (4) ,k (3) ) T i ϕ ( k (5) , , ( k (4) , , ( k (3) , T i ,T i +1 ω ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i + l (5) i l (4) i l (3) i l (2) i e ( k (5) ,k (4) ,k (3) ,k (2) ) T i ϕ ( k (5) , , ( k (4) , , ( k (3) , , ( k (2) , T i ,T i +1 ω ( k (1) ,l (1) ) T i + l (5) i l (4) i l (3) i l (2) i l (1) i e ( k (5) ,k (4) ,k (3) ,k (2) ,k (1) ) T i ϕ ( k (5) , , ( k (4) , , ( k (3) , , ( k (2) , , ( k (1) , T i ,T i +1 . The only terms here that are not explicit are the functions e ( · ,..., · ) · and ϕ ( · , · ) ,..., ( · , · ) T i , · . For t ∈ ( T i , T i +1 ], we can derive the following: e ( k ( n ) ,...,k (1) ) t = e ( k ( n ) ,...,k (1) ) T i e ∆ T i γ i ( t ) P nj =1 k ( j ) i = e P i − m =0 ∆ T m P nj =1 k ( j ) m e ∆ T i γ i ( t ) P nj =1 k ( j ) i , For example ϕ ( k (3) ,p ) , ( k (2) ,p ) , ( k (1) ,p ) T i ,t = Z tT i γ p i ( u ) e R u Ti k (3) z d z (cid:18)Z u T i γ p i ( u ) e R u Ti k (2) z d z (cid:18)Z u T i γ p i ( u ) e R u Ti k (1) z d z d u (cid:19) d u (cid:19) d u . In general ω ( k ( n ) ,l ( n ) ) ,..., ( k (2) ,l (2) ) , ( k (1) ,l (1) ) T i +1 = n +1 X m =1 ω ( k ( n − m +1) ,l ( n − m +1) ) ,..., ( k (1) ,l (1) ) T i m − Y j =0 l ( n − j ) i ! e ( k ( n − m +2) ,...,k (1) ) T i · ϕ ( k ( n − m +2) , ,..., ( k (1) , T i ,T i +1 , where whenever the index goes outside of { , . . . , n } , then that term is equal to 1. where e ( k ( n ) ,...,k (1) )0 = 1. Using integration by parts and basic integration properties, we obtain ϕ ( k,p ) T i ,t = k i (cid:16) γ pi ( t ) e k i ∆ T i γ i ( t ) − p ∆ T i ϕ ( k,p − T i ,t (cid:17) , k i = 0 , p ≥ , k i (cid:0) e k i ∆ T i γ i ( t ) − (cid:1) , k i = 0 , p = 0 , p +1 ∆ T i γ p +1 i ( t ) , k i = 0 , p ≥ . In addition, for n ≥ ϕ ( k ( n ) ,p n ) ,..., ( k (1) ,p ) T i ,t = k ( n ) i (cid:16) γ p n i ( t ) e k ( n ) i ∆ T i γ i ( t ) ϕ ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t − p n ∆ T i ϕ ( k ( n ) ,p n − , ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t − ϕ ( k ( n ) + k ( n − ,p n + p n − ) , ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t (cid:17) , k ( n ) i = 0 , p n ≥ , k ( n ) i (cid:16) e k ( n ) i ∆ T i γ i ( t ) ϕ ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t − ϕ ( k ( n ) + k ( n − ,p n − ) , ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t (cid:17) , k ( n ) i = 0 , p n = 0 , ∆ T i p n +1 (cid:16) γ p n +1 i ( t ) ϕ ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t − ϕ ( k ( n − ,p n + p n − +1) , ( k ( n − ,p n − ) ,..., ( k (1) ,p ) T i ,t (cid:17) , k ( n ) i = 0 , p n ≥ . Explicit prices.
We now express the second-order put option prices under the Heston andGARCH diffusion models in terms of the integral operators eq. (5.1) and eq. (5.2). These aregiven by the following propositions. By doing so, we essentially obtain a closed-form expressionfor the put option prices when parameters are piecewise-constant. This is a consequence ofthe relationships derived in this section. Since we are essentially expressing the results fromSection 4 with new notation, we omit the proofs for these propositions.
Proposition 5.1 (Second-order Heston explicit put option price) . The second-order put optionprice in the Heston model can be written asPut (2)H = Put BS (ˆ x, ˆ y )+ 12 ∂ xx Put BS (ˆ x, ˆ y ) S E ( ξ T − + 12 ∂ yy Put BS (ˆ x, ˆ y ) E (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19) + ∂ xy Put BS (ˆ x, ˆ y ) S E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19)(cid:27) , where E ( ξ T − ≈ exp ( v ω ( − ( κ − λρ ) ,ρ ) T + ω ( − ( κ − λρ ) ,ρ ) , ( κ − λρ,κθ ) T )! · ( v ω ( − ( κ − λρ ) ,ρ ) , ( − ( κ − λρ ) ,ρ ) , ( κ − λρ,λ ) T + ω ( − ( κ − λρ ) ,ρ ) , ( − ( κ − λρ ) ,ρ ) , ( κ − λρ,λ ) , ( κ − λρ,κθ ) T ) − . E (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19) = 2 v ω ( − κ, − ρ ) , ( − κ, − ρ ) , ( κ,λ ) T + 2 ω ( − κ, − ρ ) , ( − κ, − ρ ) , ( κ,λ ) , ( κ,κθ ) T , E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( V t − E ( V t ))d t (cid:19)(cid:27) = v (cid:16) ω ( − ( κ − λρ ) , − ρ ) T − ω ( − κ, − ρ ) T (cid:17) + ω ( − ( κ − λρ ) , − ρ ) , ( κ − λρ,κθ ) T − ω ( − κ, − ρ ) , ( κ,κθ ) T . Furthermore ˆ x = S and ˆ y = v ω ( − κ, − ρ ) T + ω ( − κ, − ρ ) , ( κ,κθ ) T . Proposition 5.2 (Second-order GARCH explicit put option price) . The second-order put op-tion price in the GARCH diffusion model can be written asPut (2)GARCH = Put BS (ˆ x, ˆ y ) + ∂ yy Put BS (ˆ x, ˆ y ) Z T (cid:18)Z t Cov( V s , V t )d s (cid:19) d t, where Z T (cid:18)Z t Cov( V s , V t )d s (cid:19) d t = (cid:16) v ω ( − κ, , ( − κ, , ( λ ,λ ) T + 2 v ω ( − κ, , ( − κ, , ( λ ,λ ) , ( − ( λ − κ ) ,κθ ) T + 2 ω ( − κ, , ( − κ, , ( λ ,λ ) , ( − ( λ − κ ) ,κθ ) , ( κ,κθ ) T (cid:17) . Furthermore ˆ x = S and ˆ y = Z T E ( V t )d t = v ω ( − κ, T + ω ( − κ, , ( κ,κθ ) T . Remark 5.1 (Fast calibration) . Since the integral operators obey a recursion, we can exploitdynamic programming to greatly speed up the calibration procedure. To implement our fastcalibration scheme, one executes the following algorithm. Let µ t ≡ µ = ( µ (1) , µ (2) , . . . , µ ( n ) ) bean arbitrary set of parameters and denote by ω t an arbitrary integral operator. • Calibrate µ over [0 , T ) to obtain µ . This involves computing ω T . • Calibrate µ over [ T , T ) to obtain µ . This involves computing ω T which is in terms of ω T ,the latter already being computed in the previous step. • Repeat until time T N .6. Numerical tests and sensitivity analysis
We test our approximation method by considering the sensitivity of our approximation withrespect to one parameter at a time. Specifically, for an arbitrary parameter set ( µ , µ , . . . , µ n ),we vary only one of the µ i at a time and keep the rest fixed. Then, we compute impliedvolatilities via our approximation method as well as the Monte Carlo for strikes correspondingto Put 10, 25 and ATM deltas. Specifically,Error( µ ) = σ IM − Approx ( µ, K ) − σ IM − Monte ( µ, K )for K corresponding to Put 10, Put 25 and ATM strikes.For all our simulations, we use 2,000,000 Monte Carlo paths, and 24 time steps per day. Thisis to reduce the Monte Carlo and discretisation errors sufficiently well. Heston sensitivity analysis.
We consider maturity times T ∈ { / , / , / , } . Westart from a ‘safe’ parameter set, which are parameters calibrated by Bloomberg USD/JPY FXoption price data, Bloomberg [8]. The safe parameter set is( S , v , r d , r f ) = (100 . , . , . , κ, θ, λ, ρ ) = (5 . , . , . , − . , T = 1 / , (5 . , . , . , − . , T = 3 / , (5 . , . , . , − . , T = 6 / , (5 . , . , . , − . , T = 1 . In our analysis, we vary one of the ( κ, θ, λ, ρ ) and keep the rest fixed. Varying κ . We vary κ over the values { , , , , , , , } . Table 6.1. κ : Heston error for ATM implied volatilities in basis points κ Table 6.2. κ : Heston error for Put 25 implied volatilities in basis points κ Table 6.3. κ : Heston error for Put 10 implied volatilities in basis points κ The Feller condition is 2 κθ > λ . We will use red text to indicate when the Feller condition is not satisfied. Note that in application, this conditionis almost always violated. That is, parameters calibrated from market data almost always violate the Fellercondition, see for example Clark [12], Da Fonseca and Grasselli [13], Ribeiro and Poulsen [31]. In fact, this is areason why practitioners now favour non-affine models. Varying θ . We vary θ over the values { , , , , , , , } . Table 6.4. θ : Heston error for ATM implied volatilities in basis points θ Table 6.5. θ : Heston error for Put 25 implied volatilities in basis points θ Table 6.6. θ : Heston error for Put 10 implied volatilities in basis points θ Varying λ . We vary λ over the values { . , . , . , . , . , . , . , . } . Table 6.7. λ : Heston error for ATM implied volatilities in basis points λ Table 6.8. λ : Heston error for Put 25 implied volatilities in basis points λ Table 6.9. λ : Heston error for Put 10 implied volatilities in basis points λ Varying ρ . We vary ρ over the values {− . , − . , − . , − . , − . , − . , − . , } . Table 6.10. ρ : Heston error for ATM implied volatilities in basis points ρ -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 01M 52.52 24.09 9.93 0 2.59 0 -1.25 -3.20 -4.08 -4.283M 53.60 21.85 4.84 -4.89 -10.63 -13.97 -15.71 -16.176M 52.02 20.13 2.84 -7.24 -13.33 -17.01 -18.99 -19.551Y 48.20 19.17 4.09 -4.33 -9.24 -12.11 -13.60 -13.91 Table 6.11. ρ : Heston error for Put 25 implied volatilities in basis points ρ -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 01M 55.43 26.08 11.36 3.660 -0.42 -2.54 -3.55 -3.853M 57.58 24.55 6.73 -3.55 -9.66 -13.29 -15.24 -15.896M 54.73 21.62 3.61 -6.90 -13.27 -17.11 -19.19 -19.831Y 49.97 19.92 4.28 -4.47 -9.55 -12.50 -14.02 -14.34 Table 6.12. ρ : Heston error for Put 10 implied volatilities in basis points ρ -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 01M 56.28 26.32 11.26 3.370 -0.79 -2.92 -3.91 -4.173M 60.03 26.03 7.61 -3.05 -9.42 -13.20 -15.28 -16.016M 58.74 24.43 5.63 -5.45 -12.23 -16.40 -18.75 -19.611Y 54.35 23.14 6.70 -2.62 -8.14 -11.45 -13.27 -13.87The above sensitivity analysis is consistent with what we expect. For example, for large maturity T , large vol-of-vol λ or large correlation | ρ | , the component-wise variance of the difference in theexpansion and evaluation point increases. Thus, when these parameters are large, we expect alower approximation accuracy. As we can see, this indeed occurs. For realistic parameter valueswe see that the magnitude of error is around 10-50bps, which can be deemed reasonable forapplication purposes.6.2. GARCH sensitivity analysis.
We start from the same ‘safe’ parameter set from Sec-tion 6.1, albeit with ρ = 0. Varying κ . We vary κ over the values { , , , , , , , } . Table 6.13. κ : GARCH error for ATM implied volatilities in basis points κ Table 6.14. κ : GARCH error for Put 25 implied volatilities in basis points κ Table 6.15. κ : GARCH error for Put 10 implied volatilities in basis points κ Varying θ . We vary θ over the values { , , , , , , , } . Table 6.16. θ : GARCH error for ATM implied volatilities in basis points θ Table 6.17. θ : GARCH error for Put 25 implied volatilities in basis points θ Table 6.18. θ : GARCH error for Put 10 implied volatilities in basis points θ Varying λ . We vary λ over the values { . , . , . , . , . , . , . , . } . Table 6.19. λ : GARCH error for ATM implied volatilities in basis points λ Table 6.20. λ : GARCH error for Put 25 implied volatilities in basis points λ Table 6.21. λ : GARCH error for Put 10 implied volatilities in basis points λ ρ is assumed to be0 always. Otherwise, the approximation behaves as we expect, with errors being larger for largematurity T and vol-of-vol λ , as the variance of the difference in the expansion and evaluationpoint will grow with these parameters.7. Error analysis
We present an explicit bound on the error term in our expansion in terms of higher ordermoments of the corresponding variance process. Specifically, this means bounding the remainderterm in the second-order expansion of the function Put BS , and for the case when ρ = 0, theerror term associated with the expansion of e R T ρ u σ u d u . We will need explicit expressions for the error terms. These are given by Taylor’s theorem,which we will present here to fix notation. We only consider the results up to second-order.
Theorem 7.1 (Taylor’s theorem for f : R → R ) . Let A ⊆ R , B ⊆ R and f : A → B be a C function in a closed interval about the point a ∈ A . Then the Taylor series of f around thepoint a is given by f ( x ) = f ( a ) + f ′ ( a )( x − a ) + 12 f ′′ ( a )( x − a ) + R ( x ) , where R ( x ) = 12 Z xa ( x − u ) f ′′′ ( u )d u = 12 ( x − a ) Z (1 − u ) f ′′′ ( a + u ( x − a ))d u. We will prefer the integration bounds to be from 0 to 1 rather than a to x , as a and x willcorrespond to random variables. Theorem 7.2 (Taylor’s theorem for g : R → R ) . Let A ⊆ R , B ⊆ R and g : A → B be a C function in a closed ball about the point ( a, b ) ∈ A . Then the Taylor series of g around thepoint ( a, b ) is given by g ( x, y ) = g ( a, b ) + g x ( a, b )( x − a ) + g y ( a, b )( y − b )+ 12 g xx ( a, b )( x − a ) + 12 g yy ( a, b )( y − b ) + g xy ( a, b )( x − a )( y − b ) + R ( x, y ) , where R ( x, y ) = X | α | =3 | α | α ! α ! E α ( x, y )( x − a ) α ( y − b ) α ,E α ( x, y ) = Z (1 − u ) ∂ ∂x α ∂y α g ( a + u ( x − a ) , b + u ( y − b ))d u, with α := ( α , α ) and | α | := α + α .7.1. Expression for error term.
The representation for the total error generated by theexpansion can be summarised by the following theorem.
Theorem 7.3 (Total expansion error) . As a functional of the underlying variance process σ ,let E BS ( σ ) and ˜ E ( σ ) correspond to the error induced by the expansion of Put BS and e R T ρ u σ u d u respectively. The error generated by Taylor expansions for a general variance process σ is givenby E ( σ ) = E BS ( σ ) + ˜ E ( σ ) , where E BS ( σ ) = X | α | =3 | α | α ! α ! E α (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19) S α ( ξ T − α (cid:18)Z T (1 − ρ u )( σ u − E ( σ u ))d u (cid:19) α ,E α (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19) = Z (1 − u ) ∂ ∂x α ∂y α Put BS ( F ( u ) , G ( u )) d u,F ( u ) := S + uS ( ξ T − ,G ( u ) := Z T (1 − ρ t ) E ( σ t )d t + u (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) , and ˜ E ( σ ) = 14 ∂ xx Put BS (ˆ x, ˆ y ) S ξ T e − R T ρ u σ u d u (cid:18)Z T ρ u ( σ u − E Q ( σ u ))d u (cid:19) · Z (1 − u ) e R T ρ m E Q ( σ m )d m e u R T ρ m ( σ m − E Q ( σ m ))d m d u. Proof.
First, we deal with the error term associated with the function Put BS , that is, E BS ( σ ).Recall the expansion of Put BS around the point (ˆ x, ˆ y ) := ( S , R T E ( σ t )(1 − ρ t )d t ) evaluated at( S ξ T , R T σ t (1 − ρ t )d t ) for a general variance process σ :Put BS (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19) = Put BS (ˆ x, ˆ y )+ ∂ x Put BS (ˆ x, ˆ y ) S ( ξ T −
1) + ∂ y Put BS (ˆ x, ˆ y ) (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + 12 ∂ xx Put BS (ˆ x, ˆ y ) S ( ξ T − + 12 ∂ yy Put BS (ˆ x, ˆ y ) (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + ∂ xy Put BS (ˆ x, ˆ y ) S ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + E BS ( σ ) . Using Theorem 7.2 for the function Put BS , this gives the error term as E BS ( σ ) = X | α | =3 | α | α ! α ! E α (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19) S α ( ξ T − α (cid:18)Z T (1 − ρ u )( σ u − E ( σ u ))d u (cid:19) α ,E α (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19) = Z (1 − u ) ∂ ∂x α ∂y α Put BS ( F ( u ) , G ( u )) d u,F ( u ) := S + uS ( ξ T − ,G ( u ) := Z T (1 − ρ t ) E ( σ t )d t + u (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) . We now investigate the error term associated with the calculation of E ξ T , that is, ˜ E ( σ ). Let uslook at this term without the expectation. ξ T = (cid:16) ξ T e − R T ρ u σ u d u (cid:17) e R T ρ u σ u d u . We expand e R T ρ u σ u d u around the expectation of the exponential’s argument under Q . Notethat ξ T e − R T ρ u σ u d u is the Radon-Nikodym derivative which changes measure from Q to Q .Expanding to second-order gives e R T ρ u σ u d u = e R T ρ u E Q ( σ u )d u Z T ρ u ( σ u − E Q ( σ u ))d u + 12 (cid:18)Z T ρ u ( σ u − E Q ( σ u ))d u (cid:19) ! + 12 (cid:18)Z T ρ u ( σ u − E Q ( σ u )d u (cid:19) Z (1 − u ) e R T ρ m E Q ( σ m )d m e u R T ρ m ( σ m − E Q ( σ m ))d m d u. Finally, the coefficient in front of ξ T in the pricing formula is ∂ xx Put BS (ˆ x, ˆ y ) S . Thus, theerror term ˜ E ( σ ) can be written as˜ E ( σ ) = 14 ∂ xx Put BS (ˆ x, ˆ y ) S ξ T e − R T ρ u σ u d u (cid:18)Z T ρ u ( σ u − E Q ( σ u ))d u (cid:19) · Z (1 − u ) e R T ρ m E Q ( σ m )d m e u R T ρ m ( σ m − E Q ( σ m ))d m d u. (cid:3) Corollary 7.1 (Total expansion error: ρ = 0) . The error due to Taylor expansions for a generalvariance process σ with ρ = 0 a.e., denoted by E ( σ ), is given by E ( σ ) = 12 (cid:18)Z T ( σ t − E ( σ t ))d t (cid:19) Z (1 − u ) ∂ yyy Put BS (cid:16) S , ˜ G ( u ) (cid:17) d u, ˜ G ( u ) := Z T E ( σ t )d t + u (cid:18)Z T ( σ t − E ( σ t )d t (cid:19) , which is just E BS ( σ ) when ρ = 0 a.e..7.2. Bounding error term.
The hope now is to be able to bound E ( E ( σ )) in terms of thehigher order moments of the variance process σ . To do this, we will need to show that thepartial derivatives ∂ Put BS ∂x α ∂y α ( F ( u ) , G ( u )) for u ∈ (0 ,
1) where α + α = 3, are functions of T and K which are bounded. Lemma 7.1.
Consider the third-order partial derivatives of Put BS , ∂ Put BS ∂x α ∂y α , where α + α = 3.Let ˆ f ( x ) := ln( x/K ) + R T (cid:16) r dt − r ft (cid:17) d t . Thenlim y ↓ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Put BS ∂x α ∂y α (cid:12)(cid:12)(cid:12)(cid:12) ˆ f ( x )=0 = ∞ . Furthermore, this is the only case where the partial derivatives explode.
Proof.
In the following, we will repeatedly let F be an arbitrary polynomial of some degree,as well as A to be an arbitrary constant. That is, they may be different on each use. FromAppendix B, it can be seen that as a function of x and y , the third-order partial derivatives areof the form A φ ( d + ) x n y m/ F ( d + , d − , √ y ) , n ∈ Z , m ∈ N . (7.1)Recall d ± = d ± ( x, y ) = ln( x/K ) + R T (cid:16) r dt − r ft (cid:17) d t √ y ± √ y,φ ( x ) = 1 √ π e − x / . Written in this form eq. (7.1), it is evident that the partial derivatives could only blow up ifeither x or y tend to 0 or infinity. We need only look at these limits independently of the othervariable. In the following, we shall say f = o ( g ) if and only if lim x →∞ f ( x ) g ( x ) = 0 and f = o ( g ) ifand only if lim x ↓ f ( x ) g ( x ) = 0.(1) For fixed x : From eq. (7.1) the partial derivatives are of the form A φ ( d + ) F ( d + ,d − , √ y ) y m/ . It canbe shown that φ ( d + ) = Ae − D y − D y , where D = ( ˆ f ( x )) and D = 1 /
8. Hence both D and D are non-negative. However,there will be two cases to consider, when D > D = 0.(a) Suppose D >
0, then ˆ f ( x ) = 0. As F is a polynomial in d + , d − , and √ y , we can saythat F ( d + , d − , √ y ) = o (1 /y M / ) and F ( d + , d − , √ y ) = o ( y M/ ) for some M, M ∈ N .Thus (cid:12)(cid:12)(cid:12)(cid:12) φ ( d + ) F ( d + , d − , √ y ) y m/ (cid:12)(cid:12)(cid:12)(cid:12) = | A | e − D y − D y o (1 /y M / ) y m/ and also (cid:12)(cid:12)(cid:12)(cid:12) φ ( d + ) F ( d + , d − , √ y ) y m/ (cid:12)(cid:12)(cid:12)(cid:12) = | A | e − D y − D y o ( y M/ ) y m/ . Then as y ↓ y → ∞ , the partial derivatives tend to 0. (b) Suppose D = 0, then ˆ f ( x ) = 0. Thus d + = √ y and φ ( d + ) = Ae − D y . Evidently, F ( d + , d − , √ y ) = P Ni =0 C i y i/ for some N ∈ N and constants C , . . . , C N . Thus (cid:12)(cid:12)(cid:12)(cid:12) φ ( d + ) F ( d + , d − , √ y ) y m/ (cid:12)(cid:12)(cid:12)(cid:12) = | A | e − D y | P Ni =0 C i y i/ | y m/ . This quantity tends to 0 as y → ∞ , as the exponential decay makes the polynomialgrowth/decay irrelevant. However, when y ↓
0, then this limit depends on the poly-nomial F . If N > m and if one of the C , C , . . . , C m are non-zero then this quantitytends to ∞ as y ↓
0. If
N < m , then this quantity tends to ∞ as y ↓
0. For each ofthe partial derivatives, it can be shown that either of these cases are satisfied. Thusthe partial derivatives tend to ∞ when y ↓ x , if ˆ f ( x ) = 0, then the partial derivatives tend to 0 if y → ∞ , andto ∞ if y ↓
0. When ˆ f ( x ) = 0, the partial derivatives tend to 0 if y ↓ y → ∞ .(2) For fixed y : From eq. (7.1), the partial derivatives are of the form φ ( d + ) F ( d + ,d − ) x n . It can beshown that φ ( d + ) = Ax − E ln( x ) − E , where E > E ∈ R . Furthermore, as F is a polynomial in d + and d − , then | F ( d + , d − ) | ≤ N X i =0 | C i | (cid:12)(cid:12) ln i ( x ) (cid:12)(cid:12) , for some N ∈ N and constants C , . . . C N . It is clear F = o ( x ) and F = o (ln N +1 ( x )). Thus (cid:12)(cid:12)(cid:12)(cid:12) φ ( d + ) F ( d + , d − ) x n (cid:12)(cid:12)(cid:12)(cid:12) = | A | x − E ln( x ) − E − n o ( x ) . Then as x → ∞ , this quantity tends to 0. In addition (cid:12)(cid:12)(cid:12)(cid:12) φ ( d + ) F ( d + , d − ) x n (cid:12)(cid:12)(cid:12)(cid:12) = | A | x − E ln( x ) − E − n o (ln N +1 ( x )) . Then as x ↓ y , the partial derivatives tend to 0 as x ↓ x → ∞ . (cid:3) We will however, be concerned with the behaviour of u ∂ ∂x α ∂y α Put BS ( F ( u ) , G ( u )), meaningwe will have to consider both arguments simultaneously, as they are both linear functions of u . Lemma 7.2.
Consider the third-order partial derivatives of Put BS , ∂ Put BS ∂x α ∂y α , where α + α = 3as well as the linear functions h , h : [0 , → R + such that h ( u ) = u ( d − c ) + c and h ( u ) = u ( d − c ) + c . Assume there exists no point a ∈ (0 ,
1) such thatlim u → a ln( h ( u ) /K ) + R T ( r dt − r ft )d t p h ( u ) = 0 and lim u → a h ( u ) = 0 . Then there exists functions M α bounded on R such thatsup u ∈ (0 , (cid:12)(cid:12)(cid:12)(cid:12) ∂ Put BS ∂x α ∂y α ( h ( u ) , h ( u )) (cid:12)(cid:12)(cid:12)(cid:12) = M α ( T, K ) . Furthermore, the behaviour of M α for fixed K and T is characterised by the functions ζ and η respectively, where ζ ( T ) = ˆ Ae − R T r ft d t e − E ˜ r ( T ) e − E ˜ r ( T ) n X i =0 c i ˜ r i ( T ) , with ˜ r ( T ) := R T ( r dt − r ft )d t and E > E ∈ R , ˆ A ∈ R , n ∈ N and c , . . . , c n are constants, and η ( K ) = ˜ AK − D ln( K )+ D N X i =0 C i ( − i ln i ( K ) , with D > , D ∈ R , ˜ A ∈ R , N ∈ N and C , . . . , C N are constants. Proof.
Under the assumptions of h and h , this supremum will be bounded by a direct appli-cation of Lemma 7.1. Next, we need to show that M α are bounded on R and behave as ζ and η for fixed K and T respectively.In the following, A denotes an arbitrary constant and F an arbitrary polynomial of some degree.They may be different on each use.(1) Behaviour in T : Fix all variables as constant except T . Then we can write the partialderivatives as Ae − R T r ft d t φ ( d + ) F ( d + , d − ) . Expanding and collecting terms with T , we can write the partial derivatives in the form Ae − R T r ft d t e − E ˜ r ( T ) e − E ˜ r ( T ) n X i =0 c i ˜ r i ( T ) = ζ ( T ) , where E > E ∈ R , n ∈ N and c , . . . , c n are constants. As ζ is a composition of poly-nomials and exponentials of ˜ r ( T ), then it is bounded for any closed interval not containing0. Now since sup t ∈ [0 ,T ] ( | r dt − r ft | ) =: R < | ˜ r ( T ) | < RT . Thus ˜ r i ( T ) = o ( T i ) and˜ r i ( T ) = o ( T i ). Hence ζ tends to 0 as T ↓ T → ∞ . Thus ζ is bounded on R + .(2) Behaviour in K : Now fix all variables as constant except K . Then the partial derivativescan be written as Aφ ( d + ) F ( d + , d − ) . Expanding and collecting terms with K , the partial derivatives can be written in the form AK − D ln( K )+ D F (ln(1 /K )) = η ( K ) , where D > D ∈ R . Then writing out the polynomial explicitly η ( K ) = AK − D ln( K )+ D N X i =0 C i ( − i ln i ( K ) , where N ∈ N and C , . . . , C N are constants. Thus | η ( K ) | ≤ | A | K − D ln( K )+ D N X i =0 | ln i ( K ) | .η is bounded for any closed interval not containing 0 since it is a composition of exponentialsand logarithms. Then as ln i ( K ) = o ( K ) and ln i ( K ) = o (ln N +1 ( K )), η tends to 0 as K ↓ K → ∞ . Thus η is bounded on R + . (cid:3) Proposition 7.1.
There exists functions M α as in Lemma 7.2 such thatsup u ∈ (0 , (cid:12)(cid:12)(cid:12)(cid:12) ∂ ∂x α ∂y α Put BS ( F ( u ) , G ( u )) (cid:12)(cid:12)(cid:12)(cid:12) ≤ M α ( T, K ) Q a.s. . Proof.
Since F and G are linear functions, then from Lemma 7.2, this claim is immediately trueif we can show that G is strictly greater than 0 for any u ∈ (0 , Q a.s.. Recall G ( u ) = (1 − u ) (cid:18)Z T (1 − ρ t ) E ( σ t )d t (cid:19) + u Z T (1 − ρ t ) σ t d t. G corresponds to the linear interpolation of R T (1 − ρ t ) E ( σ t )d t and R T (1 − ρ t ) σ t d t . It is clearsup t ∈ [0 ,T ] (1 − ρ t ) >
0. As σ corresponds to the variance process, in application this is alwayschosen to be a non-negative process such that the set { t ∈ [0 , T ] : σ t > } has non-zero Lebesguemeasure. Thus, these integrals are strictly positive and hence G is strictly larger than 0 for any u ∈ (0 , Q a.s.. (cid:3) We obtain the following corollary.
Corollary 7.2.
There exists a function M as in Lemma 7.2 such thatsup u ∈ (0 , (cid:12)(cid:12)(cid:12) ∂ yyy Put BS (cid:16) S , ˜ G ( u ) (cid:17)(cid:12)(cid:12)(cid:12) ≤ M ( T, K ) Q a.s. . Proof.
Recall ˜ G ( u ) = (1 − u ) (cid:18)Z T E ( σ t )d t (cid:19) + u Z T σ t d t. Then by the same argument in the proof of Proposition 7.1, ˜ G is strictly greater than 0 Q a.s..Hence by Lemma 7.2, the claim is true. (cid:3) Theorem 7.4 (Error bounds for general σ ) . The error term in the pricing formula is boundedas | E ( E ( σ )) | ≤ X | α | =3 C α M α ( T, K ) T α − S α (cid:0) E ( ξ T − α (cid:1) / (cid:26)Z T (1 − ρ u ) α E | σ u − E ( σ u ) | α d u (cid:27) / + C ˜ M ( T, K ) S T / e R T ρ m E Q ( σ m )d m (cid:16) E Q e R T ρ m | σ m − E Q ( σ m ) | d m (cid:17) / · (cid:18)Z T ρ u E Q | σ u − E Q ( σ u ) | d u (cid:19) / , where ˜ M ( T, K ) = ∂ xx Put BS (ˆ x, ˆ y ) is bounded on R and C = 1 / , C α = | α | α ! α ! are constants,the latter depending on α . Proof.
First, by Proposition 7.1, we have that E α ( S ξ T , R T σ t (1 − ρ t )) ≤ M α ( T, K ). ByTheorem 7.3, the error is decomposed as E ( σ ) = E BS ( σ ) + ˜ E ( σ ). We will make use of the wellknown integral inequality (cid:18)Z T | f ( u ) | d u (cid:19) p ≤ T p − Z T | f ( u ) | p d u, p ≥ . (7.2)For the term E BS ( σ ), we have |E BS ( σ ) | ≤ X | α | =3 C α M α ( T, K ) S α | ξ T − | α T α − (cid:18)Z T (1 − ρ u ) α | σ u − E ( σ u ) | α (cid:19) , where we have used the integral inequality eq. (7.2). Applying the Cauchy-Schwarz inequality,we obtain E |E BS ( σ ) | ≤ X | α | =3 C α M α ( T, K ) S α T α − (cid:0) E | ξ T − | α (cid:1) / ( E (cid:18)Z T (1 − ρ u ) α | σ u − E ( σ u ) | α (cid:19) ) / ≤ X | α | =3 C α M α ( T, K ) S α T α − (cid:0) E | ξ T − | α (cid:1) / T / · (cid:26) E (cid:18)Z T (1 − ρ u ) α | σ u − E ( σ u ) | α (cid:19)(cid:27) / , where we have used the integral inequality eq. (7.2) for the second inequality. For the term ˜ E ( σ ), notice that E ( ˜ E ( σ )) = E Q (cid:16) ˜ E ( σ ) ξ − T e R T ρ u σ u d u (cid:17) = 14 ∂ xx Put BS (ˆ x, ˆ y ) S E Q ( (cid:18)Z T ρ u ( σ u − E Q ( σ u ))d u (cid:19) · Z (1 − u ) e R T ρ m E Q ( σ m )d m e u R T ρ m ( σ m − E Q ( σ m ))d m d u ) . Now for u ∈ (0 , e u R T ρ m ( σ m − E Q ( σ m ))d m ≤ e u R T ρ m | σ m − E Q ( σ m ) | d m . Thussup u ∈ (0 , e u R T ρ m ( σ m − E Q ( σ m ))d m ≤ e R T ρ m | σ m − E Q ( σ m ) | d m . Hence Z (1 − u ) e R T ρ m E Q ( σ m )d m e u R T ρ m ( σ m − E Q ( σ m ))d m d u ≤ e R T ρ m ( E Q ( σ m )+ | σ m − E Q ( σ m ) | ) d m . Thus E | ˜ E ( σ ) | ≤ C∂ xx Put BS (ˆ x, ˆ y ) S e R T ρ m E Q ( σ m )d m · E Q ( (cid:18)Z T ρ u ( σ u − E Q ( σ u ))d u (cid:19) e R T ρ m | σ m − E Q ( σ m ) | d m ) . Finally, using the Cauchy-Schwarz inequality and the integral inequality eq. (7.2), we obtain E | ˜ E ( σ ) | ≤ C∂ xx Put BS (ˆ x, ˆ y ) S e R T ρ m E Q ( σ m )d m T / Z T ρ u E Q | σ u − E Q ( σ u ) | d u ! / · (cid:16) E Q (cid:16) e R T ρ m | σ m − E Q ( σ m ) | d m (cid:17)(cid:17) / . Furthermore, notice that ∂ xx Put BS (ˆ x, ˆ y ) = ˜ M ( T, K ), where ˜ M ( T, K ) is a function which be-haves like M ( T, K ). (cid:3) Corollary 7.3 (Error bounds for general σ : ρ = 0) . For ρ = 0 a.e., the error term in the pricingformula is bounded as | E ( E ( σ )) | ≤ CM ( T, K ) T Z T E | σ t − E ( σ t ) | d t, where C = 1 / Proof.
From Corollary 7.1, we have E ( σ ) = 12 (cid:18)Z T ( σ t − E ( σ t ))d t (cid:19) Z (1 − u ) ∂ yyy Put BS (cid:16) S , ˜ G ( u ) (cid:17) d u. Notice that ( R T | f ( t ) | d t ) ≤ T R T | f ( t ) | d t and use Corollary 7.2. Then the result is immediate. (cid:3) Conclusion
We have derived closed-form approximation formulas for European put option prices in thecontext of stochastic volatility models with time-dependent parameters. Our method involvesa second-order Taylor expansion of the mixing solution, followed by simplification of a numberof expectations via use of change of measure techniques. Such a method has been consideredby Drimus with respect to the Heston model with constant parameters [14]. We extend thismethod as we consider time-dependent parameters, which requires an additional approximationof one of the expectations. Furthermore, we obtain the expression for the error term inducedby the expansion. We give general bounds on the induced error term in terms of higher ordermoments of the underlying variance process. Under the Heston framework, we find that allterms induced by the expansion are able to be expressed as iterated integrals. Furthermore, we attempt to generalise this method to the GARCH diffusion model. We show that obtainingan expression for the solution or moments of the resulting variance process after the change ofmeasure is non-trivial to achieve. By assuming ρ = 0 a.e. in the GARCH diffusion model, weare able to work around this problem, albeit with the added assumption of uncorrelated spotand volatility movements. Additionally, we show that the iterated integrals obey a convenientrecursive property when parameters are piecewise-constant. By doing so, we make the approxi-mation formulas closed-form. In addition, we devise a fast calibration scheme which exploits therecursive property of our integral operators. Lastly, we perform a numerical error and sensitiv-ity analysis to investigate the quality of our approximation in the Heston and GARCH diffusionmodels. We find that the error is well within the range appropriate for application purposesand behaves as we expect for certain parameter values, such as long maturity, large vol-of-voland large correlation. The purely probabilistic mixing solution approach, which is the backboneof our expansion method, is very appealing due to its generality and ability to handle time-dependent parameters. Further research would be needed to combine it, with no correlationrestriction, with the type of non-affine stochastic volatility models favoured by practitioners. References [1] Alfonsi, A. et al. (2015).
Affine diffusions and related processes: simulation, theory andapplications , volume 6. Springer.[2] Al`os, E. (2006). A generalization of the Hull and White formula with applications to optionpricing approximation.
Finance and Stochastics , 10(3):353–365.[3] Andersen, L. B. (2007). Efficient simulation of the Heston stochastic volatility model.
Avail-able at SSRN 946405 .[4] Antonelli, F., Ramponi, A., and Scarlatti, S. (2010). Exchange option pricing under sto-chastic volatility: a correlation expansion.
Review of Derivatives Research , 13(1):45–73.[5] Antonelli, F. and Scarlatti, S. (2009). Pricing options under stochastic volatility: a powerseries approach.
Finance and Stochastics , 13(2):269–303.[6] Benhamou, E., Gobet, E., and Miri, M. (2010). Time-dependent Heston model.
SIAMJournal on Financial Mathematics , 1(1):289–325.[7] Black, F. and Scholes, M. (1973). The pricing of options and corporate liabilities.
TheJournal of Political Economy , 81(3):637–654.[8] Bloomberg, L. P. (2018). USD/JPY FX option price data: expiries 1,3,6,12 months, notional1000000. Retrieved July 09, 2018 from Bloomberg Terminal.[9] Carr, P. and Madan, D. (1999). Option valuation using the fast Fourier transform.
Journalof Computational Finance , 2(4):61–73.[10] Carr, P. and Willems, S. (2019). A lognormal type stochastic volatility model with quadraticdrift. arXiv preprint arXiv:1908.07417 .[11] Christoffersen, P., Jacobs, K., and Mimouni, K. (2010). Volatility dynamics for the S&P500:evidence from realized volatility, daily returns, and option prices.
The Review of FinancialStudies , 23(8):3141–3189.[12] Clark, I. J. (2011).
Foreign exchange option pricing: a practitioner’s guide . John Wiley &Sons.[13] Da Fonseca, J. and Grasselli, M. (2011). Riding on the smiles.
Quantitative Finance ,11(11):1609–1632.[14] Drimus, G. G. (2011). Closed-form convexity and cross-convexity adjustments for Hestonprices.
Quantitative Finance , 11(8):1137–1149.[15] Duffie, D., Filipovi´c, D., Schachermayer, W., et al. (2003). Affine processes and applicationsin finance.
The Annals of Applied Probability , 13(3):984–1053.[16] Gander, M. P. and Stephens, D. A. (2007). Stochastic volatility modelling in continuoustime with general marginal distributions: Inference, prediction and model selection.
Journalof Statistical Planning and Inference , 137(10):3068–3081.[17] Gatheral, J. (2011).
The volatility surface: a practitioner’s guide , volume 357. John Wiley& Sons. [18] Hagan, P. S., Kumar, D., Lesniewski, A. S., and Woodward, D. E. (2002). Managing smilerisk. The Best of Wilmott , 1:249–296.[19] Heston, S. L. (1993). A closed-form solution for options with stochastic volatility withapplications to bond and currency options.
Review of Financial Studies , 6(2):327–343.[20] Hull, J. and White, A. (1987). The pricing of options on assets with stochastic volatilities.
The Journal of Finance , 42(2):281–300.[21] Hurd, T. R. and Kuznetsov, A. (2008). Explicit formulas for Laplace transforms of sto-chastic integrals.
Markov Processes and Related Fields , 14(2):277–290.[22] Kaeck, A. and Alexander, C. (2012). Volatility dynamics for the S&P 500: further evidencefrom non-affine, multi-factor jump diffusions.
Journal of Banking & Finance , 36(11):3110–3121.[23] Kloeden, P. E. and Platen, E. (2013).
Numerical solution of stochastic differential equations ,volume 23. Springer Science & Business Media.[24] Langren´e, N., Lee, G., and Zhu, Z. (2016). Switching to non-affine stochastic volatility: Aclosed-form expansion for the Inverse-Gamma model.
International Journal of Theoretical andApplied Finance , 19(5):1650031.[25] Lewis, A. (2002). Option valuation under stochastic volatility with Mathematica code.
International Review of Economics & Finance , 11(3):331–333.[26] Lewis, A. L. (2019). Exact solutions for a GBM-type stochastic volatility model having astationary distribution.
Wilmott magazine , 2019(101):20–41.[27] Lions, P.-L. and Musiela, M. (2007). Correlations and bounds for stochastic volatilitymodels.
Annales de l’Institut Henri Poincar´e C, Analyse Non Lin´eaire , 24(1):1–16.[28] Lorig, M., Pagliarani, S., and Pascucci, A. (2017). Explicit implied volatilities for multi-factor local-stochastic volatility models.
Mathematical Finance , 27(3):926–960.[29] Mikhailov, Sergei and N¨ogel, Ulrich (2004).
Heston’s stochastic volatility model: Imple-mentation, calibration and some extensions . John Wiley and Sons.[30] Nelson, D. B. (1990). ARCH models as diffusion approximations.
Journal of Econometrics ,45(1-2):7–38.[31] Ribeiro, A. and Poulsen, R. (2013). Approximation behoves calibration.
QuantitativeFinance Letters , 1(1):36–40.[32] Romano, M. and Touzi, N. (1997). Contingent claims and market completeness in a sto-chastic volatility model.
Mathematical Finance , 7(4):399–412.[33] Sch¨obel, R. and Zhu, J. (1999). Stochastic volatility with an Ornstein–Uhlenbeck process:an extension.
Review of Finance , 3(1):23–46.[34] Van der Stoep, A. W., Grzelak, L. A., and Oosterlee, C. W. (2014). The Heston stochastic-local volatility model: Efficient Monte Carlo simulation.
International Journal of Theoreticaland Applied Finance , 17(07):1450045.[35] Wiggins, J. B. (1987). Option values under stochastic volatility: Theory and empiricalestimates.
Journal of financial economics , 19(2):351–372.[36] Willard, G. (1997). Calculating prices and sensitivities for path-independent derivativesecurities in multifactor models.
Journal of Derivatives , 5(1):45–61.
Appendix A. Mixing solution
In this appendix, we give a derivation of the result referred to as the mixing solution by Hulland White [20]. This result is crucial for the expansion methodology implemented in Section 2.Hull and White first established the expression for the case of independent Brownian motions W and B . Later on, this was extended for the correlated Brownian motions case, see [32, 36].Under a domestic risk-neutral measure Q , suppose that the spot S with variance σ follows thedynamics d S t = S t (( r dt − r ft )d t + √ σ t d W t ) , S , d σ t = α ( t, σ t )d t + β ( t, σ t )d B t , σ , d h W, B i t = ρ t d t. Theorem A.1 (Mixing solution) . Let Put = e − R T r dt d t E ( K − S T ) + . ThenPut = E n e − R T r dt d t E (cid:2) ( K − S T ) + |F BT (cid:3)o = E (cid:18) Put BS (cid:18) S ξ T , Z T σ t (1 − ρ t )d t (cid:19)(cid:19) , where Put BS is given in eq. (2.4). Proof.
By writing the driving Brownian motion of the spot as W t = R t ρ u d B u + R t p − ρ u d Z u ,where Z is a Brownian motion under Q which is independent of B , this gives the explicit pathwiseunique strong solution of S as S T = S ξ T exp (cid:26)Z T ( r dt − r ft )d t − Z T σ t (1 − ρ t )d t + Z T q σ t (1 − ρ t )d Z t (cid:27) ,ξ t := exp (cid:26)Z t ρ u √ σ u d B u − Z t ρ u σ u d u (cid:27) . First, notice that both σ and ξ are adapted to the filtration ( F Bt ) ≤ t ≤ T . Thus, it is evident that S T |F BT will have a log-normal distribution, namely S T |F BT ∼ LN (cid:0) ˜ µ ( T ) , ˜ σ ( T ) (cid:1) , ˜ µ ( T ) := ln( S ξ T ) + Z T ( r dt − r ft )d t − Z T σ t (1 − ρ t )d t, ˜ σ ( T ) := Z T σ t (1 − ρ t )d t. Hence, the calculation of e − R T r dt d t E (( K − S T ) + |F BT ) will result in a Black-Scholes like formula. e − R T r dt d t E (( K − S T ) + |F BT )= Ke − R T r dt d t N (cid:18) ln( K ) − ˜ µ ( T )˜ σ ( T ) (cid:19) − e − R T r dt d t e ˜ µ ( T )+ ˜ σ ( T ) N (cid:18) ln( K ) − ˜ µ ( T ) − ˜ σ ( T )˜ σ ( T ) (cid:19) = Ke − R T r dt d t N ln( K ) − ˜ µ ( T ) − ˜ σ ( T )˜ σ ( T ) + 12 ˜ σ ( T ) ! − S ξ T e − R T r ft d t N ln( K ) − ˜ µ ( T ) − ˜ σ ( T )˜ σ ( T ) −
12 ˜ σ ( T ) ! = Ke − R T r dt d t N ln( K/S ξ T ) − R T ( r dt − r ft )d t ˜ σ ( T ) + 12 ˜ σ ( T ) ! − S ξ T e − R T r ft d t N ln( K/S ξ T ) − R T ( r dt − r ft )d t ˜ σ ( T ) −
12 ˜ σ ( T ) ! . It is immediate that e − R T r dt d t E (( K − S T ) + |F BT ) = Put BS (cid:0) S ξ T , ˜ σ ( T ) (cid:1) . (cid:3) Appendix B. Put BS partial derivatives This appendix contains some partial derivatives for the Black-Scholes put option formula Put BS ,given in eq. (2.4). One can think of these partial derivatives as being analogous to Black-ScholesGreeks. However, these are slightly different as our Black-Scholes formulas are parametrisedwith respect to integrated variance rather than volatility. B.1.
First-order
Put BS . ∂ x Put BS = e − R T r fu d u ( N ( d + ) − ,∂ y Put BS = xe − R T r fu d u φ ( d + )2 √ y . B.2.
Second-order
Put BS . ∂ xx Put BS = e − R T r fu d u φ ( d + ) x √ y ,∂ yy Put BS = xe − R T r fu d u φ ( d + )4 y / ( d − d + − ,∂ xy Put BS = ( − e − R T r fu d u φ ( d + ) d − y . B.3.
Third-order
Put BS . ∂ xxx Put BS = ( − e − R T r fu d u φ ( d + ) x y ( d + + √ y ) ,∂ xxy Put BS = e − R T r fu d u φ ( d + )2 y ( d − d + − ,∂ xyy Put BS = ( − e − R T r fu d u φ ( d + )2 y (cid:18) d − d + − d + − d − (cid:19) ,∂ yyy Put BS = xe − R T r fu d u φ ( d + )8 y / (cid:0) ( d − d + − − ( d − + d + ) + 2 (cid:1) . B.4.
Fourth-order
Put BS . ∂ xxxx Put BS = e − R T r fu d u φ ( d + ) x y / ( d + 3 d + √ y + 2 y + 1) ,∂ xxxy Put BS = e − R T r fu d u φ ( d + )2 xy / ( d − (1 − d )) ,∂ xxyy Put BS = ( − e − R T r fu d u φ ( d + )2 xy / (cid:18)
12 ( d − + d + ) + d − d + (cid:18) − d − d + (cid:19) − (cid:19) ,∂ xyyy Put BS = e − R T r fu d u φ ( d + )8 y (cid:16) ( √ y − d + ) (cid:2) ( d − d + − − ( d − + d + ) + 2 (cid:3) + 4 [ d + ( d − − d + ) − d − − d + ] (cid:17) ,∂ yyyy Put BS = xe − R T r fu d u φ ( d + )8 y /
12 ( d − d + − ( d − d + − − ( d − d + − d − + d + ) −
12 ( d − + d + ) ( d − d + −
7) + ( d − d + − ! . Appendix C. Greeks
Expressions for second-order approximations of put option Greeks can be obtained via simplepartial differentiation of Put (2) . The Put Delta approximation is obtained via partial differen-tiation of Put (2) with respect to the underlying S . ∂ S Put (2) = ∂ x Put BS (ˆ x, ˆ y )+ 12 (cid:2) S ∂ xx Put BS (ˆ x, ˆ y ) + S ∂ xxx Put BS (ˆ x, ˆ y ) (cid:3) E ( ξ T − + 12 ∂ xyy Put BS (ˆ x, ˆ y ) E (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + [ ∂ xy Put BS (ˆ x, ˆ y ) + S ∂ xxy Put BS (ˆ x, ˆ y )] E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19)(cid:27) . The Put Gamma approximation is obtained via partial differentiation of the Put Delta approx-imation with respect to the underlying S . ∂ S S Put (2) = ∂ xx Put BS (ˆ x, ˆ y )+ 12 (cid:2) ∂ xx Put BS (ˆ x, ˆ y ) + 2 S ∂ xxx Put BS (ˆ x, ˆ y ) + S ∂ xxxx Put BS (ˆ x, ˆ y ) (cid:3) E ( ξ T − + 12 ∂ xxyy Put BS (ˆ x, ˆ y ) E (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19) + [2 ∂ xxy Put BS (ˆ x, ˆ y ) + S ∂ xxxy Put BS (ˆ x, ˆ y )] · E (cid:26) ( ξ T − (cid:18)Z T (1 − ρ t )( σ t − E ( σ t ))d t (cid:19)(cid:27) . Notice that the above expectations are the same as those from eqs. (2.6) to (2.8). Theseexpectations are made explicit in Section 3.
Appendix D. Calculation of moments
In this appendix, we derive expressions for some of the moments, mixed moments and covari-ances of the CIR and IGa processes utilised in this article. The results for the CIR process arewell known, however, we have not found a source for the IGa moments with time-dependentparameters in the literature.D.1.
Moments of the CIR process.
Let V be a CIR( v ; κ t , θ t , λ t ). It satisfies the SDEd V t = κ t ( θ t − V t )d t + λ t p V t d B t , V = v , where we assume ( κ t ) ≤ t ≤ T , ( θ t ) ≤ t ≤ T and ( λ t ) ≤ t ≤ T are time-dependent, deterministic, strictlypositive and bounded on [0 , T ]. For s < t , it can be integrated to obtain V t = V s e − R ts κ z d z + Z ts e − R tu κ z d z κ u θ u d u + Z ts e − R tu κ z d z λ u p V u d B u . (D.1)In particular, for s = 0, V t = v e − R t κ z d z + Z t e − R tu κ z d z κ u θ u d u + Z t e − R tu κ z d z λ u p V u d B u . (D.2) Proposition D.1. V has the following moments: E ( V nt ) = e − R t nκ z d z (cid:18) v n + Z t e R u nκ z d z (cid:18) nκ u θ u + 12 n ( n − λ u (cid:19) E ( V n − u )d u (cid:19) , Var( V t ) = Z t λ u e − R tu κ z d z (cid:26) v e − R u κ z d z + Z u e − R up κ z d z κ p θ p d p (cid:27) d u, Cov( V s , V t ) = e − R ts κ z d z Z s λ u e − R su κ z d z (cid:26) v e − R u κ z d z + Z u e − R up κ z d z κ p θ p d p (cid:27) d u, E ( V ms V nt ) = e − R t nκ z d z (cid:18) E ( V m + ns ) + Z ts e R u nκ z d z (cid:18) nκ u θ u + 12 n ( n − λ u (cid:19) E ( V ms V n − u )d u (cid:19) , Cov( V ms , V nt ) = E ( V ms V nt ) − E ( V ms ) E ( V nt ) , all for m, n ≥ s < t . Proof.
We give an outline for obtaining Var( V t ) and Cov( V s , V t ). The other terms follow asimilar methodology. Notice that Var( V t ) = E ( V t − E ( V t )) . Then using eq. (D.2) and E ( V t ),Var( V t ) = E (cid:18)Z t e − R tu κ z d z λ u p V u d B u (cid:19) = Z t e − R tu κ z d z λ u E ( V u )d u. Assume s < t . Using the representation of V t in terms of V s eq. (D.1), we haveCov( V s , V t ) = Cov (cid:18) V s , V s e − R ts κ z d z + Z ts e − R tu κ z d z κ u θ u d u + Z ts e − R tu κ z d z λ u p V u d B u (cid:19) = e − R ts κ u d u Var( V s ) , where we have used that V s is independent of the Itˆo integral R ts e − R tu κ z d z λ u √ V u d B u . (cid:3) D.2.
Moments of the IGa process.
Let V be an IGa( v ; κ t , θ t , λ t ). It satisfies the SDEd V t = κ t ( θ t − V t )d t + λ t V t d B t , V = v , where we assume ( κ t ) ≤ t ≤ T , ( θ t ) ≤ t ≤ T and ( λ t ) ≤ t ≤ T are time-dependent, deterministic, strictlypositive and bounded on [0 , T ]. Then for s < t , V has the explicit strong solution V t = V s Y t Y s v + R t κ u θ u /Y u d uv + R s κ u θ u /Y u d u ! . In particular, for s = 0, V t = Y t (cid:18) v + Z t κ u θ u Y u d u (cid:19) . Proposition D.2. V has the following moments: E ( V nt ) = e R t n ( n − λ z − nκ z d z (cid:18) v n + n Z t κ u θ u e − R u n ( n − λ z − nκ z d z E ( V n − u )d u (cid:19) , Var( V t ) = e − R t κ z d z Z t λ u E ( V u ) e R u κ z d z d u, Cov( V s , V t ) = Var( V s ) e − R ts κ z d z , E ( V ms V nt ) = e R t n ( n − λ z − nκ z d z (cid:16) E ( V m + ns ) e − R s n ( n − λ z − nκ z d z + n Z ts κ u θ u e − R u n ( n − λ z − nκ z d z E ( V ms V n − u )d u (cid:17) , Cov( V ms , V nt ) = E ( V ms V nt ) − E ( V ms ) E ( V nt ) , all for m, n ≥ s < t . Proof.
We show how to obtain E ( V ns V mt ). The other terms follow a similar methodology. Weconsider the differential of V n .d( V nt ) = (cid:18) nκ t θ t V n − t + (cid:18) n ( n − λ t − nκ t (cid:19) V nt (cid:19) d t + nλ t V nt d B t = ⇒ V nt = V ns + Z ts nκ u θ u V n − u + (cid:18) n ( n − λ u − nκ u (cid:19) V nu d u + Z ts nλ u V nu d B u . Multiplying both sides by V ms and taking expectation yields E ( V ms V nt ) = E ( V n + ms ) + Z ts nκ u θ u E ( V ms V n − u ) + (cid:18) n ( n − λ u − nκ u (cid:19) E ( V ms V nu )d u. Differentiating both sides in t and letting M m,ns ( t ) := E ( V ms V nt ), thendd t M m,ns ( t ) = nκ t θ t M m,n − s ( t ) + (cid:18) n ( n − λ t − nκ t (cid:19) M m,ns ( t ) . This is a first order ODE, which can be solved with the integrating factor method by integratingfrom s to t . (cid:3) Appendix E. Solutions to SDEs with linear diffusion
Suppose the diffusion U solves the SDEd U t = f ( t, U t )d t + ν t U t d ˜ B t , U = u , (E.1)where ( ν t ) ≤ t ≤ T is adapted to the Brownian filtration and f and ν satisfy some regularityconditions so that a pathwise unique strong solution for U exists. For example, f Lipschitz in x uniformly in t , and ν bounded on [0 , T ] is sufficient. Proposition E.1.
The solution to eq. (E.1) can be expressed as U t = Y t /F t , where F is a GBM(1; ν t , − ν t ), That is,d F t = ν t F t d t − ν t F t d B t , F = 1 , = ⇒ F t = exp (cid:26)Z t ν u d u − Z t ν u d ˜ B u (cid:27) , and Y solves the integral equation (written in differential form)d Y t = F t f (cid:18) t, Y t F t (cid:19) d t, Y = u . (E.2) Proof.
We essentially verify that this form of U satisfies the SDE eq. (E.1).d (cid:18) Y t F t (cid:19) = d (1 /F t ) Y t + 1 F t d Y t + d (1 /F t ) d Y t = (cid:18)(cid:26) ν t F t d B t − ν t F t d t (cid:27) + ν t F t d t (cid:19) Y t + f ( t, Y t /F t ) d t + 0= Y t F t ν t d B t + f ( t, Y t /F t ) d t = ν t U t d B t + f ( t, U t )d t.t.