Calibration of Local Volatility Model with Stochastic Interest Rates by Efficient Numerical PDE Method
RRevision HistoryRevision Date Author(s) Description a r X i v : . [ q -f i n . M F ] M a r alibration of Local Volatility Model with Stochastic Interest Rates byEfficient Numerical PDE Method Julien Hok Credit Agricole CIB, Broadwalk House, 5 Appold StLondon, EC2A 2DA, United [email protected]
Shih-Hau Tan Cuemacro, London, United Kingdom
Abstract
Long maturity options or a wide class of hybrid products are evaluated using a local volatility typemodelling for the asset price S ( t ) with a stochastic interest rate r ( t ). The calibration of the localvolatility function is usually time-consuming because of the multi-dimensional nature of the problem. Inthis paper, we develop a calibration technique based on a partial differential equation (PDE) approachwhich allows an efficient implementation. The essential idea is based on solving the derived forwardequation satisfied by P ( t, S, r ) Z ( t, S, r ), where P ( t, S, r ) represents the risk neutral probability densityof ( S ( t ) , r ( t )) and Z ( t, S, r ) the projection of the stochastic discounting factor in the state variables( S ( t ) , r ( t )). The solution provides effective and sufficient information for the calibration and pricing. ThePDE solver is constructed by using ADI (Alternative Direction Implicit) method based on an extension ofthe Peaceman-Rachford scheme. Furthermore, an efficient algorithm to compute all the corrective termsin the local volatility function due to the stochastic interest rates is proposed by using the PDE solutionsand grid points. Different numerical experiments are examined and compared to demonstrate the resultsof our theoretical analysis. Keywords: local volatility model; stochastic interest rates; hybrid, calibration; forward Fokker-Plancktype equation; alternating direction implicit (ADI) method
1. Introduction
In quantitative finance, the local volatility type model as introduced in [1, 2, 3] is widely used tomodel the price of underlying in order to capture the market volatility skew or smile in equity or foreignexchange market. It is known that with deterministic rates, the local volatility function can be obtainedwith the Dupire formula (see equation (13)) by using the European call and put option prices. Corresponding author Quantitative analyst consultant
Preprint submitted to Elsevier March 13, 2018 or long maturity options (e.g pure equity autocall derivative) or some hybrid products with a payoffinvolving interest rate and underlying asset like a best-of interest rate-equity which pays coupons of theform max (cid:20)
LIBOR, a (cid:18) S t S − (cid:19)(cid:21) , (1)where the interest rates is potentially needed to be modelled as stochastic. It is then natural to extendthe local volatility model to incorporate stochastic interest rates. This modeling framework is widelyused in the financial industry (see e.g [4, 5, 6, 7]). To perfectly match the market implied volatility, thelocal volatility type formula can be derived and is given by equation (12). We observe on top of theDupire local volatility function, there is an additional corrective term taking into account the covariancebetween the equity price and short rate. Unfortunately, this extension of the Dupire formula is not easilyapplicable for calibration over the market since there seems no immediate way to link the expectationterm with the European option prices or other liquid products.The main challenge for the implementation of the model consists in the calibration (see also discus-sion in [8]). First and foremost, it corresponds to a two factor models. Also the formula (12) requiresa correction term on top of
Dupire expression for each strike point K per maturity in the definitionof the local volatility function. Accuracy, robustness in the calibration and pricing associated with anefficient implementation are required for execution in real time. Indeed, the model will be used not onlyfor the pricing but also to compute all sensitivity factors for hedging purpose (e.g delta, gamma and vega).Many research results appeared in the last decade in different areas of quantitative finance on thelocal volatility model with stochastic rates. A brief introduction is given as following: • Theoretical results about the local volatility function and its calibration were exposed and discussedin [4, 9, 10]. • For the option pricing, in [6, 11], the authors developed expansion formulas by applying the per-turbation method using a proxy introduced by [12]. A pricing framework via Partial DifferentialEquation (PDE) approach was studied in [13] and the Crank-Nicolson scheme also the AlternatingDirection Implicit (ADI) method were applied to build the PDE solver. • In terms of model calibration, for pricing the Power Reverse Dual Currency (PRDC) derivatives,Piterbarg in [8] modeled the local volatility function for the forward foreign exchange rate usingthe constant elasticity variance (CEV) dynamic as a parametric form. A fast calibration proce-dure based on the so-called Markovian projection method was developed and the skew averagingtechnique were discussed in [14, 15]. The calibration essentially captures the slope of the impliedvolatility surface but does not exactly fit its convexity. In [7], the authors proposed a PDE calibra-tion method which bootstraps the local volatility function for each expiry T i by solving a forward3DE for the joint distribution of ( S, r ) under the forward measure Q T i . Depending on the numberof expiries in the local volatility calibration, the algorithm is potentially very intensive in compu-tations. Calibration by Monte Carlo approach using McKeans particle method was studied in [16].The authors in [17] used Malliavin calculus to derive an equation of the local volatility function,which is then solved by using fixed-point method. However the calibration quality deterioratessignificantly for far out of the money or long maturities as discussed in [16]. • For models calibration in quantitative finance, different numerical solvers were constructed to solvesimilar types of high dimensional PDEs. For example, an idea of using operator splitting methodwas provided by Ren [18] for the calibration of local volatility with stochastic volatility. A Heston-like term-structure model in FX market was examined by Tian [19] with a modified Douglas schemeand Wyns [20] with a modified Craig-Sneyd scheme. These ADI-type schemes were discussed indetailed for dealing with convection-diffusion equations in terms of the stability and second orderconvergence in [21].In this paper, we assume a Markovian setting with a stochastic differential equation (SDE) given by(2), i.e., the underlying S follows a local volatility type diffusion and the short rate r ( t ) is represented bya general form dynamic. This model covers a particular case, a local volatility diffusion for S associatedwith a Gaussian Hull-White dynamic for r (see [22]), which is widely used for pricing hybrid products(see e.g [6, 11, 7]). We focus on the model calibration and propose a methodology for an efficient im-plementation using a PDE approach. The idea starts with introducing the projection of the stochasticdiscounting factor in the state variable ( S ( t ) , r ( t )), quantity defined by the functional Z ( t, S ( t ) , r ( t )) inequation (18). Using a martingale argument, we derive the forward partial differential equation satisfiedby the product P ( t, S, r ) Z ( t, S, r ), where P ( t, S, r ) is the probability density of ( S ( t ) , r ( t )) (see our mainresult in proposition 3.3). The solution of the PDE can substantially help us to obtain the correctiveterms and perform the pricing. For getting the numerical PDE solution efficiently, a simple and intuitiveADI scheme based on an extension of the Peaceman-Rachford scheme in [23] is proposed which leadsto solve the discretized linear system with a reasonable matrix size. Moreover, we suggest an efficientalgorithm to compute all the corrective terms sequentially in strike per maturity using the solutions andpoints in the PDE grid (see section 4.2). Generally in the two dimensional case, we believe our method-ology is able to make PDE approach very efficient w.r.t Monte Carlo method (see discussion in [16]).The paper is organized as follows. Section 2 describes the hybrid equity-interest rates modelling andthe derivation of the local volatility function is explained. The calibration framework is described insection 3 to obtain the forward equation satisfied by P ( t, S, r ) Z ( t, S, r ). For its resolution, an ADI-typemethod extending Peaceman-Rachford scheme is constructed in section 4. Section 5 is devoted to thenumerical tests. Finally the conclusions and discussions are given in section 6.4 . Hybrid equity interest rates model Let’s consider the 2-d stochastic differential equations, describing the spot price S ( t ) and short rates r ( t ) under the risk neutral probability Q , defined by dS ( t ) S ( t ) = r ( t ) dt + σ ( t, S ( t )) dW ( t ) , S (0) = S ,dr ( t ) = µ ( t, r ( t )) dt + α ( t, r ( t ))( ρdW ( t ) + (cid:112) − ρ dW ( t )) , r (0) = r , (2)where ( W ( t )) t ≥ is a standard Brownian motion in R on a filtered probability space (Ω , F , ( F t ) t ≥ , Q )with the usual assumptions on the filtration ( F ( t )) t ≥ . We assume the existence and uniqueness of thesolution for (2) (see e.g theorem 3.5.5 in [24]).This model corresponds to an extension of the Dupire local volatility model (e.g [1, 2, 3]) which allowsstochastic interest rates. The local volatility function σ ( t, S ) allows the model to calibrate to the surfaceof European call price C ( T, K ) where K represents the strike and T is the maturity. For the sake ofcompleteness, we provide its derivation by following the methodology in [25]. Let’s note by Z ( t ) := e − (cid:82) t r ( u ) du , (3)and applying Tanaka’s formula to the convex but non-differentiable function Z ( t )( S ( t ) − K ) + leads to Z ( t )( S ( t ) − K ) + = ( S (0) − K ) + − (cid:90) T r ( u ) Z ( u )( S ( u ) − K ) + du + (cid:90) T Z ( u )1 S ( u ) ≥ K dS ( u ) + 12 (cid:90) T Z ( u ) dL Ku ( S ) , (4) where L Ku ( S ) is the local time of S , and ( S (0) − K ) + = max( S (0) − K, S is a continuoussemimartingale, then almost-surely (see e.g [26]) L ( t ) K ( S ) = lim (cid:15) (cid:38) (cid:15) (cid:90) t [ K,K + (cid:15) ]( S ( u )) d u . (5)Using (2), it comes Z ( t )( S ( t ) − K ) + = ( S (0) − K ) + − (cid:90) T r ( u ) Z ( u )( S ( u ) − K ) + du + (cid:90) T Z ( u )1 S ( u ) ≥ K r ( u ) S ( u ) du (6)+ (cid:90) T Z ( u )1 S ( u ) ≥ K S ( u ) σ ( u, S ( u )) dW u + 12 (cid:90) T Z ( u ) dL Ku ( S )= ( S (0) − K ) + + (cid:90) T KZ ( u )1 S ( u ) ≥ K r ( u ) du + (cid:90) T Z ( u )1 S ( u ) ≥ K S ( u ) σ ( u, S ( u )) dW u + 12 (cid:90) T Z ( u ) dL Ku ( S )= ( S (0) − K ) + + (cid:90) T KZ ( u )1 S ( u ) ≥ K ( r ( u ) − f (0 , u )) du + (cid:90) T KZ ( u )1 S ( u ) ≥ K f (0 , u ) du + (cid:90) T Z ( u )1 S ( u ) ≥ K S ( u ) σ ( u, S ( u )) dW u + 12 (cid:90) T Z ( u ) dL Ku ( S ) , (7)5here f (0 , u ) = − ∂∂u log( ZC (0 , u )) means the forward rate at time 0 for investing at time u and ZC ( t, T )is the zero coupon price at t for maturity T . We introduce the forward rate in the last equality in orderto compare below the expression of the local volatility when interest rate is stochastic, σ ( T, K ), andwhen interest rate is deterministic σ Dup ( T, K ).Assuming that the function Z ( u )1 S ( u ) ≥ K S ( u ) σ ( u, S ( u )) is a member of the class V (see def 3.1.4in [27]), namely the measurable and adapted functions f s.t E [ (cid:82) t f ( s ) ds ] < ∞ , and by taking theexpectation C ( T, K ) = C (0 , K ) + K (cid:90) T E [ Z ( u )1 S ( u ) ≥ K ( r ( u ) − f (0 , u ))] du + K (cid:90) T E [ Z ( u )1 S ( u ) ≥ K f (0 , u )] du (8)+ 12 K (cid:90) T E [ Z ( u ) δ ( S ( u ) − K ) σ ( u, K )] , with δ ( . ) the delta function at 0. Differentiating w.r.t T leads to C T ( T, K ) = KE [ Z ( T )1 S ( T ) ≥ K ( r ( T ) − f (0 , T ))] + KE [ Z ( T )1 S ( T ) ≥ K f (0 , T )] (9)+ 12 K E [ Z ( T ) δ ( S ( T ) − K ) σ ( T, K )] . Standard computations give E [ Z ( T )1 S ( T ) ≥ K ] = − ∂C ( T, K ) ∂K , (10) E [ Z ( T ) δ ( S ( T ) − K )] = ∂ C ( T, K ) ∂K . (11)Using the last two equations in (9), we obtain the expression of the local volatility σ ( T, K ) in termsof call prices C ( T, K ) σ ( T, K ) = σ Dup ( T, K ) − E [ Z ( T )( r ( T ) − f (0 , T ))1 S ( T ) >K ] K ∂ C ( T,K ) ∂K , (12)with σ Dup ( T, K ) = ∂C ( T,K ) ∂T + Kf (0 , T ) ∂C ( T,K ) ∂K K ∂ C ( T,K ) ∂K . (13) σ Dup ( T, K ) represents the
Dupire local volatility function when interest rates are deterministic.Equation (12) shows the corrections to employ on the tractable
Dupire local volatility surface in or-der to obtain the local volatility surface which takes into account the effect of stochastic interest rates. E [ Z ( T )( r ( T ) − f (0 , T ))1 S ( T ) >K ] represents the numerator of the extra term in the local volatility ex-pression. No closed form solution exists and it is not directly related to European call prices or other6iquid products. Its calculations need to be estimated. For sanity check, when interest rate becomes de-terministic, we have r ( T ) = f (0 , T ) (see equation (16)) and σ ( T, K ) reduces to σ Dup ( T, K ) in equation(12).
Remark 2.1.
Under T − forward measure Q T , the extra term can be written as E [ Z ( T )( r ( T ) − f (0 , T ))1 S ( T ) >K ] = ZC (0 , T ) E T [( r ( T ) − f (0 , T ))1 S ( T ) >K ] . (14)In the HeathJarrowMorton (HJM) framework under Q T (see e.g [28]), the forward rate f ( t, T ) is amartingale, with a vector of volatility σ ( s, T ) which can be written as f ( t, T ) = f (0 , T ) + (cid:90) t σ ( s, T ) .dW Ts . (15)Using the fact that r ( T ) = f ( T, T ), it becomes r ( T ) = f (0 , T ) + (cid:90) T σ ( s, T ) .dW Ts , (16)and with (16), the expression of the extra term in (14) shows clearly the impact of stochastic rates.By assuming (cid:82) T σ ( s, T ) ds < + ∞ , (cid:82) T σ ( s, T ) .dW Ts is a Q T martingale and E T [( r ( T ) − f (0 , T ))1 S ( T ) >K ] = Cov T [( r ( T ) − f (0 , T )) , S ( T ) >K ] , (17)which provides an interpretation to the corrective term as a covariance, under Q T , between r ( T ) − f (0 , T )and 1 S ( T ) >K .
3. Calibration
Before using a model to price any derivatives, we usually calibrate it on the vanilla market whichmeans that it is able to price vanilla options with the concerned model and the resulting implied volatili-ties match the market-quoted ones. More precisely, it is necessary to determine all parameters presentingin the different stochastic processes which define the model. In such a way, all the European option pricesderived in the model are as consistent as possible with the corresponding market ones.The calibration procedure for the two-factor model with local volatility can be decomposed into threesteps: • Parameters present in the one-factor dynamic for the interest rates, µ ( t, r ( t )) and α ( t, r ( t )), arechosen to match European swaption / cap-foors values. Methods for doing so are well developed inthe literature (see e.g [28]). 7 The correlation parameter ρ is typically chosen either by historical estimation or from occasionallyobserved prices of hybrid product involving interest rate and the underlying spot (see discussion in[11]). • After these two steps, the calibration problem consists in finding the local volatility function σ ( t, S ( t )) which is consistent with its associated implied volatility surface.Here we focus on the third step of the calibration and propose to use some martingales properties ofthe model to allow an efficient implementation.Let’s introduce the projection of the discount factor on the state variables ( S ( t ) , r ( t )), defined as Z ( t ) := E [ Z ( t ) | S ( t ) , r ( t )] = Z ( t, S ( t ) , r ( t )) , (18)with Z ( t ) given in (3) and assume that Z ∈ C , on [0 , T ] × R .Our objective is to determine the product functional P ( t, S, r ) Z ( t, S, r ) where P ( t, S, r ) representsthe joint distribution of ( S ( t ) , r ( t )). Indeed, we write E [ e − (cid:82) T r ( s ) ds ( r ( T ) − f (0 , T ))1 S ( T ) >K ] = E [ E [ e − (cid:82) T r ( s ) ds ( r ( T ) − f (0 , T ))1 S ( T ) >K /S ( t ) , r ( t )]] (19)= E [ Z ( T, S ( t ) , r ( t ))( r ( T ) − f (0 , T ))1 S ( T ) >K ] (20)= (cid:90) ( r − f (0 , T )) 1 S ≥ K ( P Z )( T, S, r ) dSdr. (21)We can then compute the corrective term (12) at least numerically. As ( S ( t ) , r ( t )) are the model statevariables, P ( t, S, r ) Z ( t, S, r ) can also be used for option pricing.For any fix T > h ( S, r ) a Borel-measurable function, let’s define the function f ( t, S, r ) = E t,S,r [ e − (cid:82) Tt r ( s ) ds h ( S ( T ) , r ( T ))] , (22)where we assume E t,S,r | h ( S ( T ) , r ( T )) | < + ∞ for all t, S, r .Using martingale argument as for the discounted Feynman-Kac theorem (see theorem 6.4.3 in [29]),we derive the partial differential equation satisfied by f ( t, S, r ): f t + rSf s + µf r + 12 S σ f ss + 12 α f rr + ρσαSf sr − rf = 0 , (23)with the terminal condition f ( T, S, r ) = h ( S, r ) ∀ ( S, r ) . (24)8n the risk neutral pricing framework, we have the following result Proposition 3.1.
For Z ( t ) and f defined, respectively, as in (18) and in (22), Z ( t ) f ( t, S ( t ) , r ( t )) is amartingale.Proof. Using the martingality of Z ( t ) f ( t, S ( t ) , r ( t )) and the Markov property of solutions in (2), for t ≥ s ,we write E [ Z ( t ) f ( t, S ( t ) , r ( t )) / F s ] = Z s f ( s, S s , r s ) , (25) E [ E [ Z ( t ) f ( t, S ( t ) , r ( t )) / F t ] / F s ] = E [ Z s f ( s, S s , r s ) / F s ] , (26) E [ E [ Z ( t ) f ( t, S ( t ) , r ( t )) /S ( t ) , r ( t )] / F s ] = E [ Z s f ( s, S s , r s ) /S s , r s ] , (27) E [ Z ( t ) f ( t, S ( t ) , r ( t )) / F s ] = Z s f ( s, S s , r s ) . (28)Since the martingale Z ( t, S ( t ) , r ( t )) f ( t, S ( t ) , r ( t )) is an Itˆo process, it must have zero drift. Calculatingthe drift term using Itˆo formula and setting it to zero give the following corollary: Corollary 3.2. Z ( t, S, r ) f ( t, S, r ) is a solution of the following partial differential equation ( Z f ) t + rS ( Z f ) s + µ ( Z f ) r + 12 S σ ( Z f ) ss + 12 α ( Z f ) rr + ρσαS ( Z f ) sr = 0 . (29)Next proposition provides the forward equation satisfied by P ( t, S, r ) Z ( t, S, r ). Proposition 3.3.
Assume the probability density P ( t, S, r ) and its derivatives decay fast enough to forlarge | ( S, r ) | to preclude boundary terms, i.e., for | ( S, r ) | → ∞ µP Z = σ s P Z = α P Z = sσαP Z = 0 , ∂ ( σ s P Z ) ∂S = ∂ ( α P Z ) ∂r = ∂ ( sσαP Z ) ∂r = ∂ ( sσαP Z ) ∂s = 0 . (30) Then P Z satisfies the following forward equation with the Dirac delta function as the initial condition ( P Z ) t + ( rSP Z ) s + ( µP Z ) r − ( s σ P Z ) ss − ( α P Z ) rr − ρ ( ασSP Z ) sr + r ( P Z ) = 0 , ( P Z )(0 , S, r ) = δ ( S − S , r − r ) . (31) or ( P Z ) t + (cid:2) S ( r − σ − ρσα r ) − σσ S S (cid:3) ( P Z ) s + [ µ − αα r − αρ ( σ + Sσ s )] ( P Z ) r − S σ ( P Z ) ss − α ( P Z ) rr − ραSσ ( P Z ) sr + (cid:2) r + µ r − ( σ + 4 Sσσ s + S ( σ ss σ + σ s )) − ( αα rr + α r ) − ρα r ( σ + Sσ s ) (cid:3) ( P Z ) = 0 , ( P Z )(0 , S, r ) = δ ( S − S , r − r ) . (32)9 roof. Let’s note by C ( R ) = { h ∈ C ( R ) with compact support } (33)and consider f ( t, S, r ) defined as in (22) for h ∈ C ( R ).Using the martingale property of Z ( t ) f ( t, S ( t ) , r ( t )) , we write for any t ≥ f (0 , S (0) , r (0)) = E [ Z ( t ) f ( t, S ( t ) , r ( t ))] (34)= (cid:90) [ P ( t, S, r ) Z ( t, S, r )] f ( t, S, r ) dSdr. (35)By taking derivative w.r.t t in (35), using (23), performing integration by parts and using zero boundaryconditions (30), we get (cid:90) f [( P Z ) t + ( rSP Z ) s + ( µP Z ) r −
12 ( s σ P Z ) ss −
12 ( α P Z ) rr − ρ ( ασsP Z ) sr + r ( P Z )] dsdr = 0 . (36)Note that when t approaches T , the function f ( t, S, r ) approaches h ( S, r ) and the last equation becomes (cid:90) h [( P Z ) t + ( rSP Z ) s + ( µP Z ) r −
12 ( s σ P Z ) ss −
12 ( α P Z ) rr − ρ ( ασsP Z ) sr + r ( P Z )] dsdr = 0 . (37)Since h ∈ C ( R ) is arbitrary, we obtain the forward equation (31).The Dirac delta initial condition means that at time t = 0 we are sure that the spot S (0) equals S ,the interest rate r (0) equals r and Z (0) = 1. Remark 3.4.
We know (see e.g proposition 11.5 in [30]) that P ( t, S, r ) satisfies the forward Fokker-Planckequation given by ( P ) t + ( rSP ) s + ( µP ) r − ( s σ P ) ss − ( α P ) rr − ρ ( ασSP ) sr = 0 , ( P Z )(0 , S, r ) = δ ( S − S , r − r ) . (38)In comparison, the equation for P Z in (31) has an additional reaction term. Also by definition, wehave for all t ≥ (cid:90) P ( t, S, r ) dSdr = 1 , (39) (cid:90) ( P Z )( t, S, r ) dSdr = E [ Z ( t, S ( t ) , r ( t ))] = ZC (0 , t ) . (40)
4. Numerical methods
The targeted equation to solve is (32) with a Dirac delta function as initial condition. For the calibra-tion, the space variables (
S, r ) can potentially be defined in an unbounded domain typically S ≥ r ∈ R (see section 5). This unbounded domain is actually truncated by ( S, r ) ∈ [ Smin, Smax ] × [ rmin, rmax ]10n the finite difference spatial discretization where Smin is close to 0, | Smax | , | rmin | , | rmax | are suffi-ciently large in the numerical experiments (see [31]).For the boundary values, we opt for Dirichlet type conditions. The values depend on the modelconsidered. In our numerical experiments at section 5, r t has Gaussian distribution for all t and Q ( S t > , ∀ t ≥
0) = 1. The values for Z ( t, S, r ) are expected to be around 1 (see Figures 13 and 14). Then,reasonably, the boundary values for P Z ( t, S, r ) are set to be 0. For the initial condition, we propose toapproximate the Dirac delta function by the Gaussian kernels i.e a family of Gaussian functions γ N ( S, r ) = 12 π √ det Σ exp − (cid:32) S − S r − r (cid:33) (cid:48) Σ − (cid:32) S − S r − r (cid:33) (41)parametrized by the parameter N > N N (see e.g [32] and [33]).Constructing fully implicit scheme for solving equation (31) may result some difficulties as the de-scritized linear system becomes banded structure and requires efficient linear solvers. Besides if theemployed grid numbers are large, the initialization of program can cause problem as the size of matrix istoo large. In the following, the details of doing time discretization of solving equation (31) is explainedby using Peaceman-Rachford scheme which can resolve these issues. Let N S , N r be the number of grid points for doing spatial discretizations on the direction S and r ,respectively, and N t be the number of grid points employed on temporal discretization. The essentialidea of Peaceman-Rachford scheme is to split the calculation in the time-marching scheme into severalsteps with respect to different spatial variables. More precisely, the scheme used for evaluating ( P Z ) n +1 from the known value ( P Z ) n can be specified as the following two steps. ADI Step 1
In the first step of ADI, the finite difference spatial discretization on S is treated implicitly and the restterms are explicitly for time step from n to n + , namely( P Z ) n + ij − ( P Z ) nij ∆ t/ C ( P Z ) n + i +1 ,j − ( P Z ) n + i − ,j S + C ( P Z ) ni,j +1 − ( P Z ) ni,j − r + C ( P Z ) n + i +1 ,j − P Z ) n + ij + ( P Z ) n + i − ,j (∆ S ) + C ( P Z ) ni,j +1 − P Z ) nij + ( P Z ) ni,j − (∆ r ) + C ( P Z ) ni +1 ,j +1 + ( P Z ) ni − ,j − − ( P Z ) ni − ,j +1 − ( P Z ) ni +1 ,j − S ∆ r + C ( P Z ) n + ij = 0 , (42)11here C = ( r j S i − S i σ i − S i ) σ i ( σ S ) i − ρσ i S i ( α r ) j ) , C = ( µ j − ρσ i α j − ρ ( σ S ) i S i α j − α j ( α r ) j ) ,C = (cid:18) − S i σ i (cid:19) , C = (cid:32) − α j (cid:33) , C = ( − ρσ i S i α j ) ,C = (2 r j +( µ r ) j − σ i − S i σ i ( σ S ) i − (( σ S ) i ) S i − σ i ( σ SS ) i S i − (( α r ) j ) − α j ( α rr ) j − ρ ( σ S ) i S i ( α r ) j − ρσ i ( α r ) j ) , here ( P Z ) nij is the discretized notation of ( P Z )( n ∆ t, i ∆ S, j ∆ r ) for i = 1 , ..., N S , j = 1 , ..., N r and n + is a dummy time step. Note the notations ( σ S ) i , ( σ SS ) i present that the evaluation of the sensitivitiesis at the discretized point S i . Similarly the notations ( α r ) j , ( α rr ) j show the values evaluated with thediscretized point r j . Then it leads to solve H ( P Z ) n + = f (( P Z ) n ) . (43)for each j , where H is a tri-diagonal matrix with size ( N S × N S ) and entries equal to a i = − C S + C (∆ S ) , b i = 2∆ t − C (∆ S ) + C , c i = C S + C (∆ S ) , where a i , b i , c i are the entries for lower, main and upper diagonal for i = 1 , ..., N S . The right hand sidevalue f (( P Z ) n ) can be evaluated with the known value ( P Z ) n at current time step and is given by f (( P Z ) n ) = ( P Z ) nij (cid:18) t + 2 C (∆ r ) (cid:19) − ( P Z ) ni,j +1 (cid:18) C r + C (∆ r ) (cid:19) + ( P Z ) ni,j − (cid:18) C r − C (∆ r ) (cid:19) − C (cid:0) ( P Z ) ni +1 ,j +1 + ( P Z ) ni − ,j − − ( P Z ) ni − ,j +1 − ( P Z ) ni +1 ,j − (cid:1) S ∆ r . (44) ADI Step 2
Once the solution at the dummy time step ( P Z ) n + is obtained, the second step of ADI is to considerthe finite difference spatial discretization on r which is treated implicitly and the rest terms are explicitlyfor time step from n + to n + 1 as follows( P Z ) n +1 ij − ( P Z ) n + ij ∆ t/ C ( P Z ) n + i +1 ,j − ( P Z ) n + i − ,j S + C ( P Z ) n +1 i,j +1 − ( P Z ) n +1 i,j − r + C ( P Z ) n + i +1 ,j − P Z ) n + ij + ( P Z ) n + i − ,j (∆ S ) + C ( P Z ) n +1 i,j +1 − P Z ) n +1 ij + ( P Z ) n +1 i,j − (∆ r ) + C ( P Z ) n + i +1 ,j +1 + ( P Z ) n + i − ,j − − ( P Z ) n + i − ,j +1 − ( P Z ) n + i +1 ,j − S ∆ r + C ( P Z ) n +1 ij = 0 , (45)for i = 1 , ..., N S , j = 1 , ..., N r and n + 1 is the time step of the pursued solution. Then it leads to solve H ( P Z ) n +1 = f (( P Z ) n + ) , (46)12or each i , where H is a tri-diagonal matrix with size ( N r × N r ) and entries equal to d j = − C r + C (∆ r ) , e j = 2∆ t − C (∆ r ) + C , f j = C r + C (∆ r ) . where d j , e j , f j are the entries for lower, main and upper diagonal for j = 1 , ..., N r . The right hand sidevalue f (( P Z ) n + ) again can be evaluated with the calculated solution ( P Z ) n + at the dummy timestep and is given by f (( P Z ) n + ) = ( P Z ) n + ij (cid:18) t + 2 C (∆ S ) (cid:19) − ( P Z ) n + i +1 ,j (cid:18) C S + C (∆ S ) (cid:19) + ( P Z ) n + i − ,j (cid:18) C S − C (∆ S ) (cid:19) − C (cid:16) ( P Z ) n + i +1 ,j +1 + ( P Z ) n + i − ,j − − ( P Z ) n + i − ,j +1 − ( P Z ) n + i +1 ,j − (cid:17) S ∆ r . (47) Remark 4.1. • In the ADI steps 1 and 2, the coefficients C i , i = 1 , .., t n . Itcorresponds to some approximations and allows simplifications. • Certain PDE schemes for solving the forward Fokker-Planck equation for P ( t, S, r ) lead to somenegative probabilities and are source of instability (see e.g [34, 35]). In our numerical examples insection 5, we have not seen these issues when solving the equation (32) for P Z . • Using the concept of generator and M -matrix in [34], Itkin proposed PDE spliting scheme s.t theintegral of the discrete solution P ( t, S, r ) is equal to 1. Here, our solution for ( P Z )( t, S, r ) shouldtheoretically satisfy the relation (40). In our experiments in section 5, integrating the numericalsolution can give slightly different results. In that case, we normalize the numerical solution for( P Z )( t, S, r ) s.t relation (40) is satisfied. Algorithm 1:
Alternative Direction Implicit Scheme
Input: initial condition ( P Z ) = δ ( S − S , r − r ), parameters Output: ( P Z ) N t for n = 0 : N t − do
1. ADI step 1: for j = 1 : N r − do Solve H ( P Z ) n + = f (( P Z ) n ) to get ( P Z ) n + ; end
2. ADI step 2: for i = 1 : N S − do Solve H ( P Z ) n +1 = f (( P Z ) n + ) to get ( P Z ) n +1 ; endend .2. Efficient implementation of the corrective terms Let’s note by K < K < ...., < K N the strike grid where we want to compute the corrective terms E [ Z ( T )( r ( T ) − f (0 , T ))1 S ( T ) >K ] at each maturity date T . We observe the corrective terms for twoconsecutive strikes are highly correlated. Indeed for K i < K i +1 , we have Adj ( K i ) = Adj ( K i +1 ) + E [ Z ( T )( r ( T ) − f (0 , T ))1 K i +1 ≥ S ( T ) >K i ] , (48)with Adj ( K ) = E [ Z ( T )( r ( T ) − f (0 , T ))1 S ( T ) >K ].So an efficient algorithm consists to calculate Adj ( K N ) first and then to use the relation (48) tocompute sequentially the other terms. With this method, computing all corrective terms for a givenmaturity consists essentially to perform one numerical integration in the whole domain of discretizationfor S and r , which allows to speed up significantly the computation time. In our numerical experiments,the integrals for calibration and pricing are estimated numerically using the PDE solutions and gridpoints. Algorithm 2:
Calibration Algorithm
Input:
Necessary parameters for using the ADI method.
Output: σ ( T i , K j ), i = 1 , ..., N T , j = 1 , ..., N K . for i = 1 : N T do Solve equation (31) to get ( P Z ) ; for j = 1 : N K do
1. Do numerical integration using equation (48) to obtain the extra term;2. Calculate the sensitivities to get σ Dup ;3. Evaluate σ ( T i , K j ) ; endend5. Numerical experiments Let’s consider the Black-Scholes economy with Hull-White stochastic interest rates model: dS ( t ) S t = r ( t ) dt + σ dW ( t ) , S (0) = S ,dr ( t ) = a ( θ ( t ) − r ( t )) dt + σ ( ρdW ( t ) + (cid:112) − ρ dW ( t )) , r (0) = r . (49) a represents the speed of reversion, θ the long term mean level. From [28], the zero coupon ZC (0 , T )price and the forward rates are given respectively by14 C (0 , T ) = A (0 , T ) e − B (0 ,T ) r (0) , (50) f (0 , T ) = − σ a + θ − ( θ − σ a − r e − aT − σ a e − aT , (51)with B (0 , T ) = a (1 − e − aT ) and A (0 , T ) = e (cid:20) ( θ − σ a )( B (0 ,T ) − T ) − σ a B (0 ,T ) (cid:21) . Remark 5.1.
We can exactly fit the term structure of interest rates being observed in the market byconsidering θ = θ ( t ) . From [28], the formula is given by θ ( t ) = 1 a ∂f (0 , t ) ∂t + f (0 , t ) + 12 (cid:16) σ a (cid:17) (1 − e − at ) . (52)Using the martingale method, the European call option price C ( T, K ) with maturity T and strike K is derived in [36] which we summarize in the next proposition. Proposition 5.2. C ( T, K ) = S N ( d ) − KZC (0 , T ) N ( d ) , (53) where n ( t ) = √ π e − t , N ( x ) = (cid:82) x −∞ n ( t ) dt , k = ln ( K ) , d = log ( S K ) − log( ZC (0 ,T ))+ (cid:82) T ˆ σ ( t ) dt √ (cid:82) T ˆ σ ( t ) dt , d = d − (cid:113)(cid:82) T ˆ σ ( t ) dt , ˆ σ ( t ) = (cid:112) σ + 2 ρσ σ X ( t ) + σ X ( t ) and X ( t ) = − ZC (0 ,T ) ∂ZC (0 ,T ) ∂r . From (50), we have X ( t ) = B (0 , t ) and (cid:90) T ˆ σ ( t ) dt = σ T + 2 ρσ σ a (cid:20) T + e − aT − a (cid:21) + σ a (cid:20) T − a (3 − e − aT + e − aT ) (cid:21) . (54)The following corollary provides analytical formulas for C T , C K and C KK used in the local volatilitycalibration expression (12). Corollary 5.3.
Using the definitions in proposition (5.2) and g ( T ) = (cid:82) T ˆ σ ( t ) dt , we have ∂C∂T ( T, K ) := C T ( T, K ) = S n ( d )2 ˆ σ ( T ) (cid:112) g ( T ) + KZC (0 , T ) f (0 , T ) N ( d ) , (55) ∂C∂K ( T, K ) := C K ( T, K ) = − ZC (0 , T ) f (0 , T ) N ( d ) , (56) ∂ C∂K ( T, K ) := C KK ( T, K ) = ZC (0 , T ) n ( d ) K (cid:112) g ( T ) . (57)The proof is given in appendix. 15his model offers tractability and we can compute the function Z ( t, y, r ) analytically whose expressionis given in the next proposition. Proposition 5.4.
Let’s define Y ( T ) := log ( S ( T )) , R ( T ) := (cid:82) T r ( s ) ds and assume the matrix Σ yr ,defined in (65), invertible. Then we have E [ e − R ( T ) | Y ( T ) , r ( T )] = exp (cid:40) − µ R − Σ tyr,R Σ − yr . (cid:32) Y ( T ) − µ y r ( T ) − µ r (cid:33) + Σ R − Σ tyr,R Σ − yr Σ yr,R (cid:41) , (58) with µ y = log ( S (0)) + µ R − σ T, (59) µ r = r (0) e − aT + θ (1 − e − aT ) , (60) µ R = r (0) (1 − e − aT ) a + θT − θa (1 − e − aT ) , (61)Σ y = Σ R + σ T + 2 σ σ ρa (cid:20) T − a (1 − e − aT ) (cid:21) , (62)Σ r = σ a (1 − e − aT ) , (63)Σ R = (cid:16) σ a (cid:17) (cid:20) T + 12 a (1 − e − aT ) − a (1 − e − aT ) (cid:21) , (64)Σ yr = Σ y ρσ σ a (1 − e − aT ) ρσ σ a (1 − e − aT ) Σ r , (65)Σ yr,R = ρσ σ a [ T − a (1 − e − aT )] (cid:0) σ a (cid:1) ( − e − aT + e − aT ) . (66) Proof.
Standard computations give Y ( T ) = log( S (0)) + R ( T ) − σ T + σ W ( T ) , (67) r ( T ) = r (0) e − aT + θ (1 − e − aT ) + σ (cid:90) T e − a ( T − t ) (cid:16) ρdW ( t ) + (cid:112) − ρ dW ( t ) (cid:17) , (68) R ( T ) = 1 a (1 − e − aT )( r (0) − θ ) + θT + σ a (cid:90) T (1 − e − a ( T − t ) ) (cid:16) ρdW ( t ) + (cid:112) − ρ dW ( t ) (cid:17) . (69)Then Y ( T ) r ( T ) R ( T ) is a Gaussian vector with mean µ = µ y µ r µ R and covariance matrix Σ yr Σ yr,R Σ tyr,R Σ R .It is well-known (see e.g chap 2.3 in [37]) the conditional distribution of R ( T ) given by Y ( T ) r ( T ) is anormal random variable with mean µ R | y,r and variance σ R | y,r given respectively by16 R | y,r = µ R + Σ tyr,R Σ − yr . Y ( T ) − µ y r ( T ) − µ r , (70) σ R | y,r = Σ R − Σ tyr,R Σ − yr Σ yr,R . (71)Using the moment generating function for a standard normal random variable N , E [ e uN ] = e u , ∀ u ∈ R , (72)then we obtain the expression in equation (58). For the numerical tests, we consider two sets of parameters:parameters set 1 set 2 S r
2% 2% σ
20% 20% σ
4% 4% ρ − a . . θ
2% 2% T σ for the dynamic of interest rates andcorrelations values ρ to better measure the impact of stochastic rates process.For each set of model parameters, we solve the PDE (32) in an uniform grid using the ADI method ex-plained in the previous section. For set 1 and set 2 respectively, the sizes of doing temporal and spatial dis-cretizations are given by ds = 0 . , dr = 0 . , dt = 0 . ds = 0 . , dr = 0 . , dt = 0 . P Z are shown in Figure 1, numerical solutions in Figure 2 andtheir discrepancy in Figure 3. For set 2, similar results are illustrated in Figure 7, Figure 8 and Figure9. For both sets, we observe the discrepancies between analytic formula and numerical solutions on thePDE grid are reasonable small . For a more accurate assessment, we perform option pricing with thenumerical solution P Z at maturity T for each set. We evaluate numerically the European call option17rices for various strikes. Numerical and analytical prices solutions are illustrated respectively in Figure4 and Figure 10 for set 1 and 2. The discrepancies are shown respectively in Figure 5 and Figure 11.For both cases, we observe very good pricing accuracy as the differences are within couple of basis points(10 − ).Also we quantify the impact of the corrective term in expression (12) for stochastic interest rates. Toavoid any scaling, we show the analytic formula in the numerator i.e E [ Z ( T )( r ( T ) − f (0 , T ))1 S ( T ) >K ] onthe grid ( T, K ). The results are illustrated in Figure 6 and Figure 12 respectively for set 1 and 2. In thefirst case with positive correlation parameter ρ = 0 .
4, we observe positive corrective terms where highervalues are concentrated around the moneyness 1. It is expected as discussed in the remark of section 2where the corrective term measures the covariance between interest rates and equity spot. Similarly, forset 2 where the correlation parameter is negative ρ = − .
4, we obtain negative values with higher levelsconcentrated around the moneyness 1.Finally, we draw the projected discounted factor Z ( T, S, r ) in Figure 13 and Figure 14 for set 1 and2 respectively. In both cases, the forms are similar and their levels vary around 1.
Figure 1: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . T = 1 . igure 2: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . T = 1 . S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . T = 1 . .5 1 1.500.20.40.60.8 K C a ll p r i c e Analytic FormulaPDE
Figure 4: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . T = 1 . −5 K D i sc r epan cy Figure 5: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . T = 1 . igure 6: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . T = 2 . igure 8: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . T = 2 . S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . T = 2 . .5 1 1.500.20.40.60.8 K C a ll p r i c e Analytic FormulaPDE
Figure 10: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . T = 2 . −4 K D i sc r epan cy Figure 11: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . T = 2 . igure 12: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = 0 . a = 0 . θ = 0 . T = 1 . igure 14: Parameters S = 1 . r = 0 . σ = 0 . σ = 0 . ρ = − . a = 0 . θ = 0 . T = 2 . In our second example, we consider a skew model given by dS ( t ) S t = r ( t ) dt + σ H ( S t ) dW ( t ) , S (0) = s ,dr ( t ) = a ( θ ( t ) − r ( t )) dt + σ ( ρdW ( t ) + (cid:112) − ρ dW ( t )) , r (0) = r , (73)where σ H ( S t ) = ν (cid:110) (1 − β + β ) β + ( β − βS t (cid:0)(cid:113) S t + β (1 − S t ) − β (cid:1)(cid:111) , (74)with ν > β ∈ [0 ,
1] shows the skew parameter.This model introduced in [39] behaves closely to the CEV model and has been used for numericalexperiments as in [40, 41]. It presents the advantage to avoid zero to be an attainable boundary andthen allows to avoid some numerical instabilities as seen in the CEV model when the underlying S isclose to 0 (see e.g [42]). It corresponds to the Black-Scholes model for β = 1 and exhibits a skew for thevolatility surface when β (cid:54) = 1. Figure 15 illustrates the impact of the parameter β on the skew of thevolatility surface. We observe that the skew increases significantly with decreasing value of β . For ex-ample with ν = 0 . , β = 0 .
2, the difference in volatility between strikes at 50% and at 100% is about 20%.25e run two family of tests by considering respectively negative correlation ρ = −
30% and positivecorrelation ρ = 30%. Other model parameters were chosen in Table 2. For both tests, we solve the PDEin equation (31) for P Z up to maturity T in an uniform grid with ds = 0 . , dr = 0 . , dt = 0 . T for various strike by numerical integrationusing PDE grid points and solution. We also run a Monte Carlo simulation for doing pricing by usingEuler discretisation for the SDE (73) with dt = and one million of paths (see e.g [37] or [43]).For set 1 (respectively for set 2) the pricing results are shown in Figure 16 (respectively in Figure 18)and their discrepancies in Figure 17 (respectively in Figure 19). We obtain very accurate results as thedifferences of prices given by using PDE and Monte Carlo methods are all within a couple of basis pointsfor all strikes in the range [0 , S r . ν β . σ a . T Table 2: Model parameters for the hyperbolic local volatility Hull-White model.
6. Conclusions and discussions
In this paper, we have proposed a new PDE-based technique for calibration on the local volatilitymodel with stochastic interest rate. The main results are the derivation of the forward equation satis-fied by P ( t, S, r ) Z ( t, S, r ) and the constructed PDE solver based on ADI scheme which leads to a moreefficient calibration method. Besides, some techniques of accelerating the calibration algorithm are alsointroduced which make the model practical for real time execution. The numerical experiments com-plement our theoretical analysis and show consistent tests results. Furthermore, the discussed model isactually a general case which can cover most of the well-known extension of local volatility models usedthese days. Therefore this calibration framework is useful for various problems in different asset classmarkets like equity, exchange or inflation.We suggest a couple of interesting avenues of research: • Here, we have focused our numerical experiments on the resolution of the forward equation (31)with two widely used hybrid models in quantitative finance: The Black-Scholes Hull-White and the26 H y pe r bo li c v o l a t ili t y β = 0.2 β = 0.4 β = 0.5 β = 0.6 β = 0.8 Figure 15: Impact of the value β on the hyperbolic local volatility σ H for fixed volatility level ν = 0 . C a ll p r i c e Monte CarloPDE
Figure 16: Parameters S = 1 . r = 3 . ν = 20%, β = 0 . σ = 4%, ρ = − . a = 0 . T = 1 . −4 K D i sc r epan cy Figure 17: Parameters S = 1 . r = 3 . ν = 20%, β = 0 . σ = 4%, ρ = − . a = 0 . T = 1 . C a ll p r i c e Monte CarloPDE
Figure 18: Parameters S = 1 . r = 3 . ν = 20%, β = 0 . σ = 4%, ρ = 0 . a = 0 . T = 1 . −4 K D i sc r epan cy Figure 19: Parameters S = 1 . r = 3 . ν = 20%, β = 0 . σ = 4%, ρ = 0 . a = 0 . T = 1 . Hyperbolic Local Volatility Hull-White models. It would be interesting to complement the testingof the PDE calibration procedure with live market data. • Through the numerical tests, our simple ADI scheme gives good convergence results and showsrobustness w.r.t strong skew and high correlations parameters. Performing numerical analysis andproviding comparisons w.r.t modern ADI schemes are parts of our future research.
7. Appendix
For the proof of corollary (5.3), we provide the following useful lemma
Lemma 7.1.
Using the definitions in proposition (5.2) and corollary (5.3), we have S n ( d ) = KZC (0 , T ) n ( d ) (75)29 roof. d − d = ( d − d )( d + d ) (76)= − (cid:112) g ( T )( d + d ) (77)= − (cid:112) g ( T ) (cid:16) d − (cid:112) g ( T ) (cid:17) (78)= − (cid:18) log (cid:18) S K (cid:19) − log ZC (0 , T ) (cid:19) (79)log (cid:18) n ( d ) n ( d ) (cid:19) = log (cid:18) K log ZC (0 , T ) S (cid:19) (80)From the last expression, we deduce directly (75).For expression (55), we write C T ( T, K ) = S n ( d ) d ,T − K [ ZC T (0 , T ) N ( d ) + ZC (0 , T ) n ( d ) d ,T ] (81)= S n ( d )( d ,T − d ,T ) − KZC T (0 , T ) N ( d ) (82)= S n ( d )2 ˆ σ ( T ) (cid:112) g ( T ) + KZC (0 , T ) f (0 , T ) N ( d ) (83)where we have used (75) in the second equality, d ,T − d ,T =
12 ˆ σ ( T ) √ g ( T ) and ZC T (0 , T ) = − ZC (0 , T ) f (0 , T )to obtain the third expression. C K ( T, K ) = S n ( d ) d ,K − ZC (0 , T ) [ N ( d ) + Kn ( d ) d ,K ] (84)Using d ,K = d ,K and result of lemma (75), we get expression (56). Finally, we obtain formula (57)by deriving (56) w.r.t K and using d ,K = − K √ g ( T ) .30 eferenceseferences