Calibrating Local Volatility Models with Stochastic Drift and Diffusion
CCalibrating Local Volatility Models with Stochastic Drift andDiffusion
Orcan ¨Ogetbil ∗ , Narayan Ganesan † , and Bernhard Hientzsch ‡ Corporate Model Risk, Wells Fargo Bank
Abstract
We propose Monte Carlo calibration algorithms for three models: local volatilitywith stochastic interest rates, stochastic local volatility with deterministic interest rates,and finally stochastic local volatility with stochastic interest rates. For each model, weinclude detailed derivations of the corresponding SDE systems, and list the requiredinput data and steps for calibration. We give conditions under which a local volatilitycan exist given European option prices, stochastic interest rate model parameters, andcorrelations. The models are posed in a foreign exchange setting. The drift term for theexchange rate is given as a difference of two stochastic short rates, domestic and foreign,each modeled by a G1++ process. For stochastic volatility, we model the variance forthe exchange rate by a CIR process. We include tests to show the convergence and theaccuracy of the proposed algorithms.
In quantitative finance, considerable amount of research focuses on modeling the market-observed smile of the implied volatility surface. There are competing approaches to tacklethis problem. Most notably, stochastic volatility and local volatility models are often ap-plied in practice to imitate the market-observed smile. Stochastic volatility models aim tocapture the volatility dynamics observed in the market. While such models often capturethe implied volatilities in certain tenor and strike ranges well, there are often ranges thatare not repriced well. Moreover, the parametric structure of these models makes proper cal-ibration computationally challenging. The local volatility models are simpler due to theirnon-parametric nature; and they are constructed to fit completely arbitrage-free impliedvolatility surfaces. However, they fail to capture the proper dynamics of implied volatilities,as they construct a static local volatility surface, which inherently assumes deterministicspot volatility for the underlying process. These shortcomings lead practitioners to combinethem in unified frameworks, called stochastic local volatility models, with the intention ofcombining the advantages of both approaches. The combined model is of special interest ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ q -f i n . M F ] S e p o the financial industry, as it would provide a methodology to work with a more completeset of risk factors associated with exotic instruments, such as options on exchange rates,foreign stocks, quantoes and baskets.Dupire’s local volatility model [1, 2] constructs a unique diffusion process that is con-sistent with the European vanilla option market quotes. Being calibrated to vanilla optionquotes, the local volatility surface has become a standard tool to capture the risk-neutralmarginal distributions of the underliers implied by market European option price quotes,and is being utilized by practitioners for pricing and risk-assessment of more exotic instru-ments. In its original formulation, the local volatility model assumes deterministic drift anddiffusion terms. In our work, we relax both of these constraints in turn, and study the gen-eralizations of the local volatility model first under a drift that is the difference between twostochastic short rates, then under stochastic diffusion, and eventually under both stochasticdrift and stochastic diffusion. While we base our set up in foreign exchange context withthe foreign exchange rate, the domestic and base short rates as the risk drivers, our resultscan be applied in equity or commodity contexts. The theoretical setup and algorithmspresented in the paper can be simplified with minimal effort to the cases where the drift ismodeled more simplistically with a single short rate driver.The theoretical framework for extending the local volatility model with a single stochas-tic rate was presented in [3]. Calibration of this model using finite difference and MonteCarlo (MC) methods was discussed in [4, 5]. As a further extension, [6] studies the modelwith two stochastic interest rates.Stochastic volatility models embody a spot volatility or variance process that is cor-related to the underlying asset. Popular choices are the Heston model [7] with variancefollowing a Cox-Ingersoll-Ross (CIR) process [8], and the Stein-Stein model [9, 10] withvolatility following an arithmetic Ornstein-Uhlenbeck process [11]. Both models admitclosed form solutions, at least in their constant parameter formulations, for pricing Euro-pean vanilla options. Atlan’s work proposes an extension to stochastic volatility modelswith a stochastic interest rate [3]. A more common extension is to join the local volatilitywith the stochastic volatility or variance in the diffusion term of the stochastic differentialequation (SDE) to form a stochastic local volatility model. In the literature, the theoreti-cal framework for stochastic local volatility has been studied in several forms with varyingdegrees of model complexity. Dupire’s “unified theory of volatility” imposes dynamics onlocal variances [12], whereas Alexander-Nogueira’s work focuses on modeling the stochasticevolution of the local volatility surface rather than of the spot volatility [13]. Lipton pro-poses a further extension to incorporate jumps in the model [14]. Calibration algorithmsfor several forms of stochastic local volatility models have been introduced in [15, 16, 17].We use a number of acronyms for the models under consideration for easy reference.The standard Dupire local volatility model where the drift is the difference between twodeterministic rates is denoted by LV2DR. The stochastic rates extension of this model isLV2SR, whereas the stochastic local volatility extension is SLV2DR. The full model withstochastic rates and stochastic local volatility is referred to as SLV2SR.While this paper reviews and makes use of existing findings and algorithms in part,to our knowledge, the following are our novel contributions: SLV2SR model; completederivations of the SDE systems that are used during Monte Carlo calibration for all threemodels; presentation of Monte Carlo calibration algorithms for the SLV2DR and SLV2SR2odels; and convergence and accuracy studies of all calibration algorithms presented.The structure of this paper is as follows. In Section 2, we give a detailed descriptionof the algorithm to calibrate the local volatility model subject to stochastic interest rates(LV2SR), including derivations of the relevant formulas and SDEs in forward measure. Sec-tion 3 describes the two-fold calibration algorithm of the stochastic local volatility model(SLV2DR): First the underlying Heston model is calibrated to match near at-the-money for-ward (ATMF) market quotes; then the leverage function is calibrated to the entire impliedvolatility surface. Section 4 blends the algorithms from the previous sections to developa calibration algorithm for our ultimate model: stochastic local volatility with stochasticinterest rates (SLV2SR). As the models being discussed are theoretically built on top of eachother, the calibration algorithm for the SLV2SR model references and uses the results of thealgorithms for calibrating the simpler models in earlier sections. In this paper stochasticinterest rates are modeled by single factor G1++ short rate processes. Other short ratemodels can be adapted to our framework as well, as our methodology is not tightly coupledwith the choice for the type of the short rate process. Each of these sections are accompa-nied by tests that measure the convergence and the accuracy of the calibrated models. InSection 5 we summarize our findings and exhibit a study to compare the calibrated models.Appendix A proves some identities used in the measure transformations prescribed in thispaper. Let S t be the exchange rate, that is the amount of domestic currency needed to buy oneunit of foreign currency. In the domestic risk neutral measure Q DRN the exchange rate isassumed to follow the LV2SR local volatility model, dS t = (cid:104) r dt − r ft (cid:105) S t dt + σ ( S t , t ) S t dW S (DRN) t , (2.1)where σ ( S t , t ) is the state dependent diffusion coefficient that is commonly referred to as local volatility , and S t is denominated in domestic currency. The domestic short rate r dt and the foreign short rate r ft follow G1++ processes. In particular, the domestic short rateevolves in domestic risk neutral measure as r dt = x dt + φ dt ,dx dt = − a dt x dt dt + σ dt dW d (DRN) t , (2.2)whereas the foreign short rate evolves in foreign risk neutral measure Q FRN as r ft = x ft + φ ft ,dx ft = − a ft x ft dt + σ ft dW f (FRN) t . (2.3)Here φ it are the shift functions that are calibrated to market yield curves; a it are the meanreversion coefficients, and σ it are the volatility coefficients, with i = d, f .3ur derivations and computations assume constant coefficients of correlation betweenthe returns of the underlying assets; however we note that it is straightforward to generalizeour findings to more advanced models with time-dependent or even stochastic coefficientsof correlation. The first step is to derive the evolution of the foreign short rate in the domestic risk neutralmeasure. For now assume that the evolution of the foreign short rate in domestic riskneutral measure has the form dx ft = g ( · , t ) dt + σ ft dW f (DRN) t (2.4)for some drift function g ( · , t ) of the underlying assets of the SDE system and time that weare going to determine.For any asset V t denominated in domestic currency, the discounted asset price is amartingale under the domestic risk neutral measure. Defining the domestic money marketaccount as B dt = exp[ (cid:82) t r du du ], we have V B d = V = E Q DRN (cid:20) V t B dt (cid:21) . Likewise, the discounted value of V t S t , that is the price of the asset V t denominated in theforeign currency, is a martingale under the foreign risk neutral measure. Defining the foreignmoney market account as B ft = exp[ (cid:82) t r fu du ], we have V S B f = V S = E Q FRN (cid:34) V t S t B ft (cid:35) . Therefore the Radon-Nikodym derivative [18] writes d Q FRN d Q DRN = V t B dt V S V V t S t B ft = S t B ft S B dt . (2.5)The exchange rate process (2.1) is a geometric Brownian motion SDE which has the solution S t = S exp (cid:20)(cid:90) t (cid:18) r du − r fu − σ ( S u , u )2 (cid:19) du + (cid:90) t σ ( S u , u ) dW S (DRN) u (cid:21) = S B dt B ft exp (cid:20) − (cid:90) t σ ( S u , u ) du + (cid:90) t σ ( S u , u ) dW S (DRN) u (cid:21) . Thus, the Radon-Nikodym derivative (2.5) becomes d Q FRN d Q DRN = exp (cid:20) − (cid:90) t σ ( S u , u ) du + (cid:90) t σ ( S u , u ) dW S (DRN) u (cid:21) . (2.6)4hen, according to the Girsanov theorem dW S (FRN) t = dW S (DRN) t − σ ( S t , t ) dt (2.7)is a Brownian motion under the foreign risk neutral measure Q FRN .Let ρ Sf be the coefficient of correlation between the Brownian motions W S (FRN) t and W f (FRN) t , that is d (cid:10) W S (FRN) , W f (FRN) (cid:11) t = ρ Sf dt. Following Lemma A.1, the foreign shortrate process (2.4) evolves in foreign risk neutral measure Q FRN as dx ft = (cid:104) g ( · , t ) + ρ Sf σ ft σ ( S t , t ) (cid:105) dt + σ ft dW f (FRN) t . (2.8)Comparing (2.3) and (2.8), we find that g ( · , t ) = − a ft x ft − ρ Sf σ ft σ ( S t , t ) . Collectively, we can write the three processes in the domestic risk neutral measure as dS t = (cid:104) r dt − r ft (cid:105) S t dt + σ ( S t , t ) S t dW S (DRN) t ,dx dt = − a dt x dt dt + σ dt dW d (DRN) t , r dt = x dt + φ dt ,dx ft = (cid:104) − a ft x ft − ρ Sf σ ft σ ( S t , t ) (cid:105) dt + σ ft dW f (DRN) t , r ft = x ft + φ ft . (2.9) T -forward measure For computational ease, that is to decouple the discounting terms from expectations, wewould like to transform (2.9) to the domestic T -forward measure. We take the zero couponbond P d ( t, T ) maturing at time T as the num´eraire. We have P d ( T, T ) = 1. Under themeasure Q T defined by this num´eraire, the discounted price of an asset is a martingale, V P d (0 , T ) = E Q T (cid:20) V T P d ( T, T ) (cid:21) = E Q T [ V T ] . We arrive at the following Radon-Nikodyn derivative, d Q T d Q DRN = V T B dT V P d (0 ,T ) V V T = 1 B dT P d (0 , T ) = exp (cid:104) − (cid:82) T r du du (cid:105) P d (0 , T ) . (2.10)Since the domestic short rate follows a G1++ process, this expression can be written as(see Lemma A.2 for proof) d Q T d Q DRN = exp (cid:20) − (cid:90) T σ du b d ( u, T ) dW d (DRN) u − (cid:90) T (cid:16) σ du b d ( u, T ) (cid:17) du (cid:21) , (2.11)with b d ( t, T ) ≡ (cid:90) Tt e − (cid:82) vt a dz dz dv. dW d (T) t = dW d (DRN) t + b d ( t, T ) σ dt dt (2.12)is a Brownian motion under the domestic T -forward measure Q T . This allows us to writedown the domestic short rate process (2.2) as dx dt = (cid:104) − a dt x dt − b d ( t, T )( σ dt ) (cid:105) dt + σ dt dW d (T) t . (2.13)Let ρ Sd be the coefficient of correlation between the Brownian motions W S (DRN) t and W d (DRN) t , that is d (cid:10) W S (DRN) , W d (DRN) (cid:11) t = ρ Sd dt. Following Lemma A.1, the exchange rateprocess (2.1) evolves in domestic T -forward measure Q T as dS t = (cid:104) r dt − r ft − ρ Sd b d ( t, T ) σ dt σ ( S t , t ) (cid:105) S t dt + σ ( S t , t ) S t dW S (T) t . (2.14)Finally, let ρ df be the coefficient of correlation between the Brownian motions W d (DRN) t and W f (DRN) t , that is d (cid:10) W d (DRN) , W f (DRN) (cid:11) t = ρ df dt. Following Lemma A.1, the foreignshort rate process from (2.9) evolves in domestic T -forward measure Q T as dx ft = (cid:104) − a ft x ft − ρ Sf σ ft σ ( S t , t ) − ρ df b d ( t, T ) σ dt σ ft (cid:105) dt + σ ft dW f (T) t . (2.15)Collecting everything, dS t = (cid:104) r dt − r ft − ρ Sd b d ( t, T ) σ dt σ ( S t , t ) (cid:105) S t dt + σ ( S t , t ) S t dW S (T) t ,dx dt = (cid:104) − a dt x dt − b d ( t, T )( σ dt ) (cid:105) dt + σ dt dW d (T) t , r dt = x dt + φ dt ,dx ft = (cid:104) − a ft x ft − ρ Sf σ ft σ ( S t , t ) − ρ df b d ( t, T ) σ dt σ ft (cid:105) dt + σ ft dW f (T) t , r ft = x ft + φ ft (2.16)describe the evolutions of the exchange rate, domestic short rate, and foreign short rateprocesses under the domestic T -forward measure Q T . Note that the above SDE system isdifferent than what is given in [6]. The standard formulation of the local volatility model [1, 2] with deterministic interest rateshas been studied extensively in the literature. The local volatility surface can be computedfrom a call option price surface that can be constructed by market quotes of call optionprices C = C ( K, T ) as [19] σ = ∂C∂T + ( r dT − r fT ) K ∂C∂K + r fT C K ∂ C∂K . (2.17)When the interest rates have stochastic dynamics, the above equation generalizes to [20] σ = ∂C∂T − P d (0 , T ) E Q T [( r dT K − r fT S T ) { KK (cid:105) ∂C BS ∂w (cid:20) − yw ∂w∂y + ∂ w∂y + (cid:16) ∂w∂y (cid:17) (cid:16) − − w + y w (cid:17)(cid:21) , (2.20)where the Black-Scholes model price C BS = C BS ( T, y, w ) and its derivatives are given by C BS ( T, y, w ) = P d (0 , T ) F T [ N ( d ) − e y N ( d )] ,∂C BS ∂w = 12 P d (0 , T ) F T e y N (cid:48) ( d ) w − ,∂C BS ∂y = − P d (0 , T ) F T e y N ( d ) ,∂C BS ∂T = − f f (0 , T ) C BS + ∂C BS ∂w ∂w∂T + (cid:18) ∂C BS ∂y + ∂C BS ∂w ∂w∂y (cid:19) ( f f (0 , T ) − f d (0 , T )) . Here the instantaneous forward rate is defined as f i (0 , T ) ≡ − ∂ log P i (0 ,T ) ∂T = − P i (0 ,T ) ∂P i (0 ,T ) ∂T ,with i = d, f . Moreover N ( · ) is the cumulative Gaussian probability distribution function;7 = − yw − + w , and d = d − w . r dT and r fT are the time T values of the domestic andforeign short rates. Similar to the call price surface formulation, equation (2.20) reduces to(2.19) in the deterministic interest rates limit. Analogous to the call price formulation case,the evaluation of the local volatility requires the construction of the total implied variancesurface .Dupire’s equations (2.18) and (2.20) use the local volatility on both sides since thecomputation of the expectation on the right hand side is under the dynamics that involvesthe local volatility function. This can be used to define iterative approaches. We found thatthe bootstrapping approach presented below typically yields satisfactory results without anyiterative refinement. Inputs for calibration
Our calibration routine expects the following quantities as inputfor local volatility calibration: • Spot FX rate S • Market implied volatility Σ(
K, t ) for FX rate • Market yield curves P d (0 , t ) and P f (0 , t ) • For both domestic and foreign rates, G1++ model parameters mean reversion, volatil-ity and shift function calibrated to market data. See [21] for example calibrationmethods for both constant and time dependent cases. • Coefficients of correlation between the underlying assets: the FX rate, the domesticand foreign short rates
Steps for calibration
In our framework, we calibrate the local volatility surface timeslice by time slice, in a bootstrapping fashion. Let t i ; i = 1 , . . . , n be the increasing sequenceof (positive) times where we will perform the calibration.1. Using the market implied volatility Σ( K, t ), generate a vanilla call option price surface C ( K, t ) interpolator or a total implied variance surface w ( y, t ) interpolator. Theinterpolator must be able to compute the partial derivatives appearing in the localvolatility expressions.2. For the first time slice t , evaluate the deterministic equation (2.17) or (2.19) tocompute the FX local volatilities for a predetermined range of strikes. This steprequires no Monte Carlo simulation. As a result, obtain local volatility values to beused in the simulation until time t in the subsequent calibration steps.3. For each of the subsequent time slices t j , j >
1, simulate the SDE system (2.16) up totime t j . Compute the Monte Carlo estimate for the expectation appearing in (2.18)or (2.20) for a predetermined range of strikes. Use these equations to obtain thelocal volatility values. These local volatility values will be used during subsequentsimulation steps from time t j to time t j +1 . This step is first performed with j = 2and is then repeated for the remaining time slices.8he strike grid can be chosen to be uniform across all calibration time slices. In thisapproach, however, the strike grid needs to be sufficiently large to cover attainable values ofthe FX rate at long expiries. This in turn would result in unreliable local volatility valuesat short expiries and strikes in the far wings. To overcome this problem, we suggest using amore adequate strike grid at each calibration time slice, e.g. one that spans a predeterminednumber of standard deviations away from the ATMF strike value. In order to test the validity of the calibrated local volatility surface, one needs to use it forpricing. The prices generated by Monte Carlo simulation by this method, however, havetwo sources of Monte Carlo errors. First, the estimation of the expectations appearing in(2.18) or (2.20) is subject to Monte Carlo error. Second, the evaluation of the Monte Carloaverage of the payoffs computed during the pricing introduces an additional source of error.Keeping the number of paths high in one of these two steps will allow us to study theconvergence of the other.For both calibration and pricing, we also generate the antithetic conjugate paths toreduce the variance. Therefore when we talk about a simulation with N as the number ofpaths, the actual total number of paths simulated is 2 N .We use EURUSD market data as of 2020-04-30 in our tests. The domestic T -forwardmeasure SDE system (2.16) is simulated via forward Euler discretization in both calibrationand pricing steps. We vary the number of paths used during calibration while we fix thenumber of paths for the pricing simulation at 100000.In this test, we use the total implied variance formulation (2.20). The Monte Carloestimate for the local volatility is given byˆ σ = ∂C BS ∂T − P d (0 , T ) ˆ E ∂C BS ∂w (cid:20) − yw ∂w∂y + ∂ w∂y + (cid:16) ∂w∂y (cid:17) (cid:16) − − w + y w (cid:17)(cid:21) , where ˆ E is the Monte Carlo estimate for the expectation appearing in (2.20). In thisformulation, the Monte Carlo error δ ˆ E of this estimate translates to the error in localvolatility as δ ˆ σ LV = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P d (0 , T ) δ ˆ E σ LV ∂C BS ∂w (cid:20) − yw ∂w∂y + ∂ w∂y + (cid:16) ∂w∂y (cid:17) (cid:16) − − w + y w (cid:17)(cid:21) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . For every point of the calibrated surface, we estimate the Monte Carlo error in the localvolatility values. In the end we price with three volatility surfaces: the original calibratedsurface, the original calibrated surface bumped down by 2 Monte Carlo errors, and theoriginal calibrated surface bumped up by 2 Monte Carlo errors.The total implied variance interpolator has a time slice every 0.05 years, and 100 log-moneynesses per slice. The log-moneyness points are spanned uniformly over 3.5 standarddeviations from the ATMF strike value for each slice. Similarly, the local volatility surface is9alibrated to have a time slice every 0.05 years. The strike grid of the local volatility surfacehas 200 points, spanned uniformly over 3 standard deviations from the ATMF strike value.The maximum simulation time step is set at 0.01 years for both calibration and pricing.We price a set of vanilla call options expiring at 10 years. Since we have the impliedvolatility surface data, we can compare the Monte Carlo prices to analytical Black-Scholesvanilla call option prices.Figure 2.1: LV2SR: Repricing of vanilla call options with local volatility surfaces calibratedwith varying number of paths. The pricing simulation is done with 100000 paths to keepthe simulation Monte Carlo error small. The pricing is done with the original calibratedsurface as well as bumped down and bumped up surfaces according to Monte Carlo errorintroduced during the calibration. Comparison to Black-Scholes prices is made by observingthe differences between the Monte Carlo and the analytical prices. One can observe that thepricing differences become comparable in magnitude to the simulation Monte Carlo errorsat around 1000 calibration paths.In Figure 2.1 we see that the convergence is achieved quickly with a relatively low numberof paths at 1000, where the pricing differences begin to be comparable in magnitude to thesimulation Monte Carlo errors. The relatively low number of paths allow fast calibration.For the convergence of the calibrated local volatility surface, two surfaces calibrated with100,000 and 50,000 calibration runs respectively, were studied by looking at the difference10etween the surfaces in Figure 2.2. For the surfaces the x-axis represents the spot valuesof the underlier, the y-axis the time and z-axis the corresponding local volatility values.It can be seen that difference between the surfaces is roughly 1% of the magnitude of thelocal volatilities and the convergence of the surface is achieved for relatively small numberof runs.Figure 2.2: LV2SR: Calibrated local surfaces with 100000 calibration runs (left), 50000calibration runs (middle) and difference between the surfaces (right). The difference betweenthe surfaces is roughly 1% of the magnitude of the local volatilities.In order to study the effects of repricing, a local volatility surface calibrated with 100000paths was used for all subsequent tests presented in this section. The effect of number ofrepricing simulation runs for the given local volatility surface is shown in Figure 2.3 againstthe call option price and the corresponding MC error. As the MC error decreases with thenumber of simulation paths, so does the difference between the analytical and Monte-Carlopriced call option values. 11igure 2.3: LV2SR: Repricing call options at 100 uniformly spaced strikes at maturityT=9.95 years, each with 100000 simulation Monte-Carlo runs (lower-right), 50000 MC sim-ulation runs (lower-left), 10000 MC simulation runs (upper-right) and 1000 MC simulationruns (upper-left), using the local volatility surface calibrated with 100,000 MC simulationruns. The MC errors and the difference between the analytical and MC computed pricesdecreases with the number of simulation runs.Furthermore, the calibrated local volatility surface was used to reprice the call optionsat multiple maturities and various strikes in the strike-grid to generate the so-called callprice surface. The 100 maturies that are uniformly spaced between T=0 and 9.95 years,and 100 strikes per maturity were used to generate the call price surface shown in Figure2.4a. The difference between the Monte-Carlo repriced call option values and analyticalBlack-Scholes call option prices assuming constant interest rates and volatilty is shown inFigure 2.4b, which is found to be less than 0.1% of the call option price.12 a) Repriced call option surface at all strikes andmaturities in the grid (b) Difference in analytical call surface andrepriced call surface
Figure 2.4: LV2SR: (a) Call options repriced at 100 uniformly spaced maturities betweenT=0 to 9.95 years and 100 strikes per maturity and (b) the difference between Monte-Carloand Black-Scholes analytical price.Next, the market implied volatility and the implied volatility recovered from the Monte-Carlo repriced options, by inverting the option price to evaluate Black-Scholes impliedvolatility, are compared in Figure 2.5 at maturity T=9.95. Additionally the implied volatil-ity recovered from MC prices with ± ± ± ± The standard local volatility model with deterministic interest rates (LV2DR) dS t = (cid:104) r dt − r ft (cid:105) S t dt + σ LV ( S t , t ) S t dW S (DRN) t (3.1)can be extended to incorporate a stochastic nature in the diffusion term by replacing σ LV ( S t , t ) with L ( S t , t ) √ U t where L ( S t , t ) is the leverage function and U t is the varianceprocess. A common choice for U t is the Cox-Ingersoll-Ross (CIR) process [8]. With thischoice the SLV2DR SDE system is dS t = (cid:104) r dt − r ft (cid:105) S t dt + L ( S t , t ) (cid:112) U t S t dW S (DRN) t ,dU t = κ t ( θ t − U t ) dt + ξ t (cid:112) U t dW U (DRN) t , (3.2)with coefficient of correlation ρ SU between the two Brownian drivers. This model can beseen as an augmentation of the Heston model [7], with κ t , θ t , and ξ t representing the time-dependent mean reversion, long term variance, and vol-of-vol parameters. Together withthe initial variance V , they form the set of Heston model parameters that will be calibrated14o market data as we will describe below. The leverage function L ( S t , t ) is to be calibratedto recover market option prices.The standard local volatility σ LV ( S t , t ) is related to the leverage function L ( S t , t ) as[20] σ LV ( x, t ) = L ( x, t ) E Q (DRN) [ U t | S t = x ] . (3.3)The main idea is to have the Heston parameters recover market vanilla option pricesnear ATMF strikes. The leverage function will then serve as a correction factor at the wingsof the volatility surface. In the limit where the calibration function is set to L ( S t , t ) = 1, the stochastic local volatilitymodel (3.2) reduces to the Heston model with time dependent coefficients, dS t = (cid:104) r dt − r ft (cid:105) S t dt + (cid:112) U t S t dW S (DRN) t ,dU t = κ t ( θ t − U t ) dt + ξ t (cid:112) U t dW U (DRN) t . (3.4)Here we assume constant correlation between the two Brownian motions, (cid:68) dW S (DRN) , dW U (DRN) (cid:69) t = ρdt. To improve calibration accuracy, one can trivially extend the model to admit time dependentcorrelation. However, we found that our simpler setup is sufficient for our purposes. The setof parameters we need to calibrate are the mean reversion κ t , the long term variance θ t , thevol-of-vol ξ t , the coefficient of correlation ρ , and the initial variance U . Several methodshave been studied in literature to calibrate these parameters, including an asymptotic ap-proximation [22], or a semi-analytical approach computing the characteristic function andusing control variates to regularize the numerical integration [23]. While these methods canbe directly applied to our setup, we choose a simpler approach of calibrating the parametersusing a PDE solver. Inputs for calibration
Our calibration routine expects the following quantities as inputfor Heston model with piecewise constant coefficients calibration: • Spot FX rate S • Market implied volatility Σ(
K, t ) for FX rate • Market yield curves P d (0 , t ) and P f (0 , t ). Since the rates are deterministic, we have r it = f i (0 , t ) , i = d, f where the instantaneous forward rate can be computed from themarket yield curves, f i (0 , t ) ≡ − ∂ log P i (0 ,t ) ∂t . We note that the expectation in [20] is given in domestic T -forward measure Q T , which is equal todomestic risk neutral measure Q (DRN) under deterministic rates. teps for calibration The calibration is done in a bootstrapping fashion. Let t i ; i =1 , . . . , n be the increasing sequence of (positive) times where we will perform the calibration.At each calibration step, we make sure the Feller condition 2 κ t i θ t i > ξ t i is satisfied [8, 24]in order to keep the variance U t strictly positive, by adding a suitable penalty function tothe objective function. We choose the simplex algorithm [25] to do the optimization. Theparameters κ t , θ t , and ξ t are assumed to be piecewise constant. • For all calibration time slices t i and predetermined ranges of strikes K j , compute agrid of market vanilla call or put option prices using the input market implied volatilityΣ( K j , t i ). • Using the PDE pricer, solve for all five parameters κ t , θ t , ξ t , ρ , and U to matchthe market vanilla option prices expiring at t , generated in the first step. • For the subsequent slices t i ; i >
1, using the results of the previous slices in the PDEpricer, solve for the three parameters κ t i , θ t i , and ξ t i to match the market vanillaoption prices expiring at t i .Since our primary goal of calibrating the Heston model parameters is to recover marketquotes around ATMF strikes, we choose the strike grid to cover a narrower range than inthe subsequent calibration routine of the leverage function. The strike grid is chosen uniformly over 1 standard deviation away in both directions fromthe ATMF strike value for each slice. The calibrated Heston parameters obtained by theprocedure outlined above are visualized in Figure 3.1.Figure 3.1: Heston model: Time-Series of calibrated model parametersThe figure shows that the calibrated parameters are mostly within their expected upperand lower boundaries. These parameters are used to calibrate the leverage function for thestochastic local volatility model as outlined in the subsequent section. We test the validityof the calibrated Heston parameters, which corresponds to the SLV2DR model (3.2) withthe leverage function L ( S t , t )) set to 1.0, by pricing vanilla options by simulation. Thuswe can examine these call option prices close to ATMF. The repriced call option values atmaturity and the corresponding recovered implied volatilities are shown in Figure 3.2. Itcan be seen that the recovered call option price as well as the implied volatility are within2MC error for strikes close to ATMF, and near market quoted prices and implied volatilitiesfor strikes within the calibration region. 16igure 3.2: Heston model: Difference in repriced vs analytical call options (left), impliedvolatility from repriced call options vs market implied volatility at maturity(right) (T=9.95).The dashed lines indicate the calibration range.We perform an implied volatility recovery test at various maturities by inverting thesimulated Heston model prices to compute the Black-Scholes implied volatility. We comparethis with the market implied volatility in Figure 3.3. The recovered implied volatility isfound to be in good agreement with the market implied volatility for strikes within thecalibration range, as indicated in the plot. There are various methods proposed in the literature to estimate the conditional expectationin (3.3). The standard approach involves solving a forward Kolmogorov PDE that describesthe forward evolution of the probability density function of the underlier; e.g. the FX rate.Here we introduce two methods that can be implemented by Monte Carlo simulation.
Suppose we simulate N Monte Carlo paths ( S i , U i ) , i = 0 , ..., N − t . At time t , we sort these paths by S i . Let us denote by ( ˆ S i , ˆ U i ) the sorted pairs. We divide thesepairs into M bins, each bin containing N/M pairs. We compute the bin averages as˜ S j = MN NM − (cid:88) k =0 ˆ S NM j + k , ˜ U j = MN NM − (cid:88) k =0 ˆ U NM j + k . with j = 0 , ..., M −
1. By computing the interpolation function for ˜ S j against ˜ U j , we canestimate the expectation in (3.3) for a given S t .17igure 3.3: Heston model: Implied volatility computed from repriced call options vs themarket implied volatility at various maturities. The dashed lines indicate the calibrationrange at each maturity. The idea is to linearly regress the variance values U t against basis functions f n ( · ) of theunderlying spot rate process S t . After simulating N Monte Carlo paths ( S i , U i ) , i =0 , ..., N − t , we compute the regression coefficients a n by solving the leastsquares problem ˆ U t = (cid:88) n a n f n ( S t ) . (3.5)Standard monomials or orthogonal polynomials with appropriate limits can be used as basisfunctions. For example, if we use a constant term and the first two orders of monomials,we need to solve ˆ U t = a + a S t + a S t . (3.6)After computing the regression coefficients a n , we can use this regression equation to eval-uate the expected value of U t for a given S t , which gives us the conditional expectation in(3.3). 18 nputs for calibration Our calibration routine expects the following quantities as inputfor stochastic local volatility calibration: • Spot FX rate S • Market implied volatility Σ(
K, t ) for FX rate • Market yield curves P d (0 , t ) and P f (0 , t ). Since the rates are deterministic, we have r it = f i (0 , t ) , i = d, f where the instantaneous forward rate can be computed from themarket yield curves, f i (0 , t ) ≡ − ∂ log P i (0 ,t ) ∂t . • Heston model parameters κ t , θ t , ξ t , ρ , and U calibrated to market vanilla option pricesas in Section 3.2. The strike grid is chosen to be near ATMF strikes. In particular ithas to be smaller than (e.g. one third of) the local volatility strike range. Steps for calibration
The calibration is done in a bootstrapping fashion. Let t i ; i =1 , . . . , n be the increasing sequence of (positive) times where we will perform the calibration.After computing the leverage function values for a time slice t j , the values are used duringthe subsequent simulations to estimate the leverage function values at the next time slice.1. Using the market implied volatility Σ( K, t ), generate a vanilla call option price surface C ( K, t ) interpolator or a total implied variance surface w ( y, t ) interpolator. Theinterpolator must be able to evaluate the partial derivatives appearing in the localvolatility expressions.2. Compute the deterministic local volatility σ LV ( K, t ) by (2.17) or (2.19) for all calibra-tion time slices t i and predetermined ranges of strikes.3. Simulate the SDE system (3.2) up to time t j . Compute the Monte Carlo estimate forthe conditional expectation appearing in (3.3) for the same range of strikes from theprevious step by using one of the approaches described in the previous sections. Usethis equation to obtain the leverage function values L ( K, t ). These leverage functionvalues will be used during subsequent simulation steps from time t j to time t j +1 .The strike grid can be chosen to be uniform across all calibration time slices. In thisapproach, however, the strike grid needs to be sufficiently large to cover attainable valuesof the FX rate at long expiries. This in turn would result in unreliable local volatilityvalues at short expiries and strikes in the far wings. To overcome this problem, we suggestusing a more adequate strike grid at each calibration time slice, e.g. one that spans apredetermined number of standard deviations away from the ATMF strike value. For theHeston parameters calibration, we suggest using a similar grid with a smaller range, e.g.one third of the local volatility grid range. In order to test the calibration convergence, the differences between the calibrated leveragefunctions with 100,000 and 50,000 MC paths are visualized in Figure 3.4 across the strikegrid at various time slices. As the figure shows, the differences are minor, indicating theachievement of convergence. 19igure 3.4: SLV2DR: Difference in leverage function with respect to the calibration Monte-Carlo runsThe leverage function calibrated with 100,000 MC paths was used for all subsequenttests including repricing and implied volatility recovery. The call option was repriced forstrikes at maturity with the above leverage function, with various number of simulation MCpaths as shown in Figure 3.5.Figure 3.5: SLV2DR: Repricing the call options with different simulation Monte-Carlo runsusing leverage function calibrated with 100,000 Monte-Carlo runs20e perform a simulation to reprice vanilla call options across the strike grid at variousmaturities to obtain the so called call surface. The difference between the analytical callprice surface and repriced surface using 100,000 MC paths at each strike and maturity isvisualized in Figure 3.6. The differences appear to be small compared to option prices. (a) SLV2DR: Repriced call surface at all strikesand maturities in the grid (b) SLV2DR:Difference between analytical andMC repriced call surfaces
Figure 3.6: SLV2DR: Call options repriced at all strikes and maturities and the differencefrom BS analytical priceThe implied volatility for the repriced call options at maturity was obtained by invertingoption values with Black-Scholes formula. The recovered implied volatility is comparedto market implied volatility as shown in Figure 3.7. This is also performed for optionprices obtained by bumping them by ± This model can be seen an extension to both local volatility model with stochastic inter-est rates (LV2SR) and stochastic local volatility model with deterministic interest rates(SLV2DR). Considering an FX rate process with stochastic local volatility and stochasticinterest rates, the SLV2SR SDE system can be written in the domestic risk neutral measure22
DRN as dS t = (cid:104) r dt − r ft (cid:105) S t dt + L ( S t , t ) (cid:112) U t S t dW S (DRN) t ,dU t = κ t ( θ t − U t ) dt + ξ t (cid:112) U t dW U (DRN) t ,r dt = x dt + φ dt ,dx dt = − a dt x dt dt + σ dt dW d (DRN) t ,r ft = x ft + φ ft ,dx ft = (cid:104) − a ft x ft − ρ Sf σ ft L ( S t , t ) (cid:112) U t (cid:105) dt + σ ft dW f (DRN) t , (4.1)with coefficients of correlation ρ SU (FX rate/FX variance), ρ Sd (FX rate/domestic interestrate), ρ Sf (FX rate/foreign interest rate), ρ Ud (FX variance/domestic interest rate), ρ Uf (FX variance/foreign interest rate), and ρ df (domestic interest rate/foreign interest rate)between the respective Brownian motions.One can relate the standard local volatility σ LV ( S t , t ) to the leverage function L ( S t , t )by [20] σ LV ( x, t ) = L ( x, t ) E Q T [ U t | S t = x ] . (4.2)Since our calibration algorithm presented in the next section demands Monte Carlo estima-tion of the above expectation in the domestic T -forward measure Q T , we need to formulatethe SDE system in this measure. The T -forward measure SDEs for the exchange rate,domestic short rate and domestic short rate are derived the same way as in Section 2.1.2.The remaining computation is for the stochastic variance SDE, which we present here.Let ρ Ud is coefficient of correlation between the Brownian motions W U (DRN) t and W d (DRN) t ,that is d (cid:10) W U (DRN) , W d (DRN) (cid:11) t = ρ Ud dt. The domestic risk neutral measure Q DRN to do-mestic T -forward measure Q T transformation is characterized by (2.12), Following LemmaA.1, the variance process (4.1) evolves in domestic T -forward measure Q T as dU t = (cid:104) κ t ( θ t − U t ) − ρ Ud b d ( t, T ) σ dt ξ t (cid:112) U t (cid:105) dt + ξ t (cid:112) U t dW U (T) t . (4.3)Collecting everything, dS t = (cid:104) r dt − r ft − ρ Sd b d ( t, T ) σ dt L ( S t , t ) (cid:112) U t (cid:105) S t dt + L ( S t , t ) (cid:112) U t S t dW S (T) t ,dU t = (cid:104) κ t ( θ t − U t ) − ρ Ud b d ( t, T ) σ dt ξ t (cid:112) U t (cid:105) dt + ξ t (cid:112) U t dW U (T) t ,r dt = x dt + φ dt ,dx dt = (cid:104) − a dt x dt − b d ( t, T )( σ dt ) (cid:105) dt + σ dt dW d (T) t ,r ft = x ft + φ ft ,dx ft = (cid:104) − a ft x ft − ρ Sf σ ft L ( S t , t ) (cid:112) U t − ρ df b d ( t, T ) σ dt σ ft (cid:105) dt + σ ft dW f (T) t (4.4)describe the evolutions of the exchange rate, exchange rate variance, domestic short rate,and foreign short rate processes under the domestic T -forward measure Q T .23 .2 Calibration of the Leverage Function The standard forward Kolmogorov PDE approach to solve the conditional expectation in(4.2) suffers from the curse of dimensionality, as this is now a 4D problem. Similarly, thebinning approach utilized in Section 3.4.1 for the model with deterministic interest ratescan not be applied directly, at least without any simplifying assumptions, as the sorting ofthe underliers becomes nontrivial.As in the case of stochastic local volatility model with deterministic interest rates, thecalibration is done in a bootstrapping fashion, after computing the leverage function valuesfor a time slice t , the values are used during the subsequent simulation to estimate theleverage function values at the next time slice. The Heston model parameters are assumedto be calibrated to match an appropriate subset of market data.The idea is to linearly regress the variance values U t against basis functions f n ( · ) ofthe underlying spot rate process S t , and the two interest rate processes r dt and r ft . Aftersimulating N Monte Carlo paths ( S i , r d,i , r f,i , U i ) , i = 0 , ..., N − t , we computethe regression coefficients a n by solving the least squares problemˆ U t = (cid:88) n a n f n ( S t , r dt , r ft ) . (4.5)Standard monomials or orthogonal polynomials can be used as basis functions. For example,if we use a constant term and the first two orders of monomials for all underliers, we needto solve ˆ U t = a + a S t + a S t + a x dt + a x dt + a x ft + a x ft . (4.6)Note that in this example we used the x dt and x ft as the basis functions instead of theshort rates r dt and r ft ; the deterministic parts φ dt and φ ft of the latter can be absorbed intoother coefficients. After computing the regression coefficients a n , we can use this regressionequation to evaluate the expected value of U t for given S t , r dt , and r ft , which gives us theconditional expectation in (4.2). Inputs for calibration
Our calibration routine expects the following quantities as inputfor local volatility calibration: • Spot FX rate S • For both domestic and foreign rates, G1++ model parameters mean reversion, volatil-ity and shift function calibrated to market data. See [21] for some calibration methodsfor both constant and time dependent cases. • Coefficients of correlation between all underlying assets: the FX rate, its variance, thedomestic and foreign short rates • Local volatility (with stochastic rates) surface data, as calibrated in Section 2.2 • Heston model (with deterministic rates) parameters, as calibrated in Section 3.224 teps for calibration
In our framework, we calibrate the leverage function surface timeslice by time slice, in a bootstrapping fashion. Let t i ; i = 0 , . . . , n be the increasing sequenceof times where we will perform the calibration. The first time slice is t = 0.1. For t , (4.2) simplifies to σ LV ( K, = L ( K, U , where σ LV ( K, t ) is the localvolatility with stochastic rates. Use this equation to evaluate the leverage functionfor the t slice for a predetermined range of strikes.2. For each of the subsequent positive time slices t j , j ≥
1, simulate the SDE system(4.4) up to time t j . Compute the Monte Carlo estimate for the expectation appearingin (4.2) for a predetermined range of strikes. Use this equation to obtain the lever-age function values. These leverage function values will be used during subsequentsimulation steps from time t j to time t j +1 .The strike grid for the leverage function L ( K, t ) is chosen to be the same as the strike gridfor the local volatility σ LV ( K, t ) to avoid any inaccuracy introduced by interpolation.
In order to calibrate the SLV2SR model, the Heston model parameters computed during theSLV2SR model calibration and the local volatility function calibrated for the LV2SR modelwith 100,000 MC paths were used. The procedure outlined above was followed in orderto iteratively compute the leverage function for each time slice. In order to evaluate theexpectation appearing in (4.2), 100,000 MC paths followed by regression with monomials oforder 2 for each of the underlying regressors ( S t , x dt , x ft ) along with the constant coefficient a as in (4.6) were used. The results for the calibration and repricing tests are presentedbelow.The effect of number of repricing simulation runs for the given local volatility surface isshown in Figure 4.1 against the call option price and the corresponding MC error. As theMC error decreases with the number of simulation paths, so does the difference betweenthe analytical and Monte-Carlo priced call option values.25igure 4.1: SLV2SR: Repricing call options with various simulation Monte-Carlo runs us-ing leverage function calibrated for SLV2DR with 100,000 MC paths and multi-regressionapproach with 100,000 MC paths.Furthermore, the local volatility surface was used to reprice call options at multiplematurities and various strikes in the strike-grid to generate the so-called call price surface.100 maturies uniformly spaced between T=0 and 9.95 years, and 100 strikes per maturitywere used to generate the call price surface shown in Figure 4.2a. The difference betweenthe Monte-Carlo repriced call option values and analytical Black-Scholes call option pricesassuming constant interest rate and volatilty is shown in Figure 4.2b, which is found to beless than 1% of the call option price. 26 a) SLV2SR: Repriced call surface at all strikesand maturities in the grid (b) SLV2SR: Difference between analytical andMC repriced call surfaces Figure 4.2: SLV2SR: Call options repriced at all strikes and maturities and the differencefrom BS analytical priceNext, the market implied volatility and the implied volatility recovered from the MonteCarlo repriced options, by inverting the Black-Scholes formula, are compared in Figure 4.3at maturity T=9.95. In addition the implied volatility recovered from MC prices ± ± ± We studied the convergence and the vanilla option repricing accuracy of the LV2DR,SLV2DR and SLV2SR models calibrated with the proposed algorithms. While all threemodels perform decently with recovering market implied volatilities, we found that as themodels get more complex, e.g. when they have higher number of parameters, one needs toincrease the number of simulation paths to maintain the accuracy, which in turn results inincreased calibration time.We are now in a position to compare the pricing inaccuracies for the three main modelswe considered in this paper. With fixed number of calibration and simulation paths, weobserve that the SLV2SR model (see Figures 4.1, 4.2, 4.3) reprices market quotes slightlymore accurately than the SLV2DR model (see Figures 3.5, 3.6, 3.7). In the meantime, wesee that our simplest model, LV2SR (see Figures 2.3, 2.4, 2.5) gives the smallest repricingerrors among the three models.Recovering market quotes is clearly a property desired for any pricing model. In ourcase, the market quotes are vanilla option prices or an implied volatility surface. Yet the28isk neutral price of a vanilla option depends on the terminal distribution of the underlier,for which there is a closed form solution. Therefore the three models we are considering arenot strictly mandatory for pricing vanilla options. What about the instruments for whichthe payoff depends on the intermediate values of the underliers?To study the pricing of path dependent options, we consider an up-and-out barriercall option with 5-year maturity struck at ATMF strike. Without the barrier feature, theinstrument becomes a plain European vanilla call option, which has an analytical solutionunder the Black-Scholes model. The models under consideration price the option withvaluation date 2020-04-30 and 100,000 simulation paths as
Model Price Error
Analytical (BS) 0.084379LV2DR 0.084553 0.000293LV2SR 0.084633 0.000286SLV2DR 0.084336 0.000293SLV2SR 0.084651 0.000289We note that the analytical price is within the Monte Carlo error for all the four models,in consistency with the findings of the previous sections.Now we turn the up-and-out barrier on and set the barrier position at 1.25 times theATMF strike. The barrier is active throughout the lifetime of the option. The standardBlack-Scholes model admits an analytical solution for such barriers when the interest ratesare constant [26]. The models under consideration price the option with the same valuationdate and number of simulation paths as
Model Price Error
Analytical (BS) 0.028555LV2DR 0.028705 0.000108LV2SR 0.028487 0.000107SLV2DR 0.029990 0.000114SLV2SR 0.030084 0.000112We find that the LV2DR and LV2SR model prices are within 2 Monte Carlo errors ofthe analytical price. However the SLV2DR and SLV2SR model prices are observed not toconverge to the analytical BS price. This shows us the impact of stochastic volatility onthe barrier option price.While the stochasticity of the local volatility has a clear impact on the price of pathdependent instruments, the stochasticity of interest rates have little effect under standardmarket conditions, e.g. when the interest rate volatilities are low. We expect this effect tobe more prominent in stressed environments with higher interest rate volatilities, which iswhat we study next.Consider the stochastic rates extension of the Black Scholes model (BS2SR) dS t = (cid:104) r dt − r ft (cid:105) S t dt + σ S S t dW S (DRN) t ,dx dt = − a dt x dt dt + σ dt dW d (DRN) t , r dt = x dt + φ dt ,dx ft = (cid:104) − a ft x ft − ρ Sf σ ft σ S (cid:105) dt + σ ft dW f (DRN) t , r ft = x ft + φ ft . (5.1)29his model can be seen as a special case of LV2SR with flat local volatility, which allows usto incorporate results from Section 2. Using (A.9) and Itˆo’s lemma one can write the SDEsfor the domestic and foreign zero coupon bonds in domestic risk neutral measure as dP d ( t, T ) P d ( t, T ) = r dt dt − σ dt b d ( t, T ) dW d (DRN) t ,dP f ( t, T ) P f ( t, T ) = (cid:104) r ft + ρ Sf σ S σ ft b f ( t, T ) (cid:105) dt − σ ft b f ( t, T ) dW f (DRN) t . (5.2)Using (2.11) and Lemma A.1, these can be written in the domestic T -forward measure as dP d ( t, T ) P d ( t, T ) = (cid:20) r dt + (cid:16) σ dt b d ( t, T ) (cid:17) (cid:21) dt − σ dt b d ( t, T ) dW d (T) t ,dP f ( t, T ) P f ( t, T ) = (cid:104) r ft + σ ft b f ( t, T ) (cid:16) ρ Sf σ S + ρ df σ dt b d ( t, T ) (cid:17)(cid:105) dt − σ ft b f ( t, T ) dW f (T) t . (5.3)Similarly, the exchange rate process written in the domestic T -forward measure is given by(c.f. (2.14)) dS t = (cid:104) r dt − r ft − ρ Sd b d ( t, T ) σ dt σ S (cid:105) S t dt + σ S S t dW S (T) t . (5.4)The forward value of the exchange rate is F ( t, T ) ≡ E Q (T) [ S T | F t ] = S t P f ( t, T ) P d ( t, T ) , (5.5)which is a martingale under the T -forward measure Q (T) . Its SDE can be computed from(5.3), (5.4), and application of Itˆo’s lemma, dF ( t, T ) F ( t, T ) = σ S dW S (T) t + σ dt b d ( t, T ) dW d (T) t − σ ft b f ( t, T ) dW f (T) t . (5.6)Since the diffusion process above is a linear combination of correlated Brownian motions,we can extract the total implied variance easily,Σ T =( σ S ) T + 2 σ S (cid:90) T (cid:104) ρ Sd σ dt b d ( t, T ) − ρ Sf σ ft b f ( t, T ) (cid:105) dt + (cid:90) T (cid:20)(cid:16) σ dt b d ( t, T ) (cid:17) − ρ df σ dt b d ( t, T ) σ ft b f ( t, T ) + (cid:16) σ ft b f ( t, T ) (cid:17) (cid:21) dt. (5.7)With this quantity one can write down the Black formula for the time zero value of a vanillacall option as C ( K, T ) = P d (0 , T ) P f (0 , T ) (cid:104) F (0 , T ) N ( ˜ d ) − KN ( ˜ d ) (cid:105) , (5.8)with ˜ d = log F (0 ,T ) K + Σ T Σ √ T and ˜ d = ˜ d − Σ √ T .The integrals in (5.7) can be evaluated numerically. Therefore, for a given marketimplied volatility Σ at a given maturity T and strike K , one can solve this quadratic equation30o find the BS2SR volatility σ S that will reproduce the market quotes. Conversely, the righthand side of (5.7) dictates the lower bound of total implied variance for which there is asolution for the BS2SR model. The lower bound can be evaluated, by taking the derivativewith respect to σ S and setting it to zero,Σ T = (cid:90) T (cid:20)(cid:16) σ dt b d ( t, T ) (cid:17) − ρ df σ dt b d ( t, T ) σ ft b f ( t, T ) + (cid:16) σ ft b f ( t, T ) (cid:17) (cid:21) dt − T (cid:18)(cid:90) T (cid:104) ρ Sd σ dt b d ( t, T ) − ρ Sf σ ft b f ( t, T ) (cid:105) dt (cid:19) . (5.9)If the market total implied variance at a given maturity T and strike K is lower than this,the BS2SR model will not have a solution. As a consequence, the local volatility extensionof this model (LV2SR) will not be calibratable, that is the evaluation of (2.18) during theapplication of the calibration algorithm will lead to imaginary values for local volatilitygiven the already fixed parameters of the interest rate models and the correlations. Thissignifies that no real and positive local volatility exists that would reproduce the marketquotes for vanilla options.As a study, we take the market data as of 2020-04-30, but vary the interest rate modelparameters and corrrelations, and compare the minimum total implied variance admittedby the BS2SR model for a number of values of interest rate volatilities and correlations. Inparticular for a mean reversion a dt = a ft = 0 .
01, and volatilities σ dt = σ ft = 0 . , .
05 we lookat the minimum total implied variances admitted by the BS2SR model for various values of ρ df . Figure 5.1 shows that for interest rate volatilities set to 0.02, the market total impliedFigure 5.1: Minimum total implied variances (TIV) allowed by the BS2SR model for variousvalues of interest rate volatilities and correlations compared to market TIV. ρ Sd = 0 . ρ Sf = 0 . ρ df is below -0.4. In case the interest rate volatilities31re set to 0.05, the market total implied variance is attainable only when the correlation ρ df is quite high, around 0.8. A Supplementary computations
The following finding is utilized in measure transformations throughout the paper.
Lemma A.1.
Let dX t = a ( X t , t ) dt + b ( X t , t ) dW X (A) t , (A.1) with P (cid:18)(cid:90) t (cid:0) | a ( X s , s ) | + (cid:12)(cid:12) b ( X s , s ) (cid:12)(cid:12)(cid:1) ds < ∞ (cid:19) = 1 , ∀ t ≥ , be an SDE describing the evolution of asset X t under measure Q A whose Brownian motion W X (A) t is correlated to another Brownian motion W Y (A) t with a coefficient of correlation ρ XY , that is d (cid:10) W X (A) , W Y (A) (cid:11) t = ρ XY dt . Let W Y (B) t be a Brownian motion under anequivalent measure Q B characterized by the transformation d Q B d Q A = exp (cid:20) − (cid:90) t c ( · , s ) ds − (cid:90) t c ( · , s ) dW Y (A) s (cid:21) , (A.2) such that dW Y (B) t = dW Y (A) t + c ( · , t ) dt, (A.3) with P (cid:16)(cid:82) t | c ( · , s ) | ds < ∞ (cid:17) = 1 , ∀ t ≥ , and c ( · , t ) is a function of underlying assets of theSDE system, and time. Then the evolution of X t under measure Q B is described by dX t = [ a ( X t , t ) − ρ XY b ( X t , t ) c ( · , t )] dt + b ( X t , t ) dW X (B) t , (A.4) where W X (B) t is a Brownian motion under Q B .Proof. We can decompose the Brownian motion W X (A) t into W Y (A) t and an independentBrownian motion Z t , that is d (cid:10) W Y (A) , Z (cid:11) t = 0 , as dW X (A) t = ρ XY dW Y (A) t + (cid:113) − ρ XY dZ t . (A.5)We note that, as a result of the multi-dimensional Girsanov theorem, Z t is a Brownianmotion under Q B .Now we can use (A.3) and (A.5) to write the process (A.1) as dX t = a ( X t , t ) dt + b ( X t , t ) (cid:20) ρ XY dW Y (A) t + (cid:113) − ρ XY dZ t (cid:21) = a ( X t , t ) dt + b ( X t , t ) (cid:20) ρ XY (cid:16) dW Y (B) t − c ( · , t ) dt (cid:17) + (cid:113) − ρ XY dZ t (cid:21) = [ a ( X t , t ) − ρ XY b ( X t , t ) c ( · , t )] dt + b ( X t , t ) dW X (B) t , (A.6)32ith dW X (B) t = ρ XY dW Y (B) t + (cid:113) − ρ XY dZ t = dW X (A) t + ρ XY c ( · , t ) dt. (A.7)We use the following result during the change to T -forward measure in Section 2.1.2.The limiting case with constant coefficients was investigated in [27]. Here we study thegeneral case with time dependent coefficients. Lemma A.2.
For the G1++ model (2.2) describing the evolution of the short rate r t , r t = x t + φ t ,dx t = − a t x t dt + σ t dW t , (A.8) where φ t is the deterministic shift function that is calibrated to market yield curve; a t is themean reversion coefficient, and σ t is the volatility coefficient, the following identity holds, exp (cid:20) − (cid:90) Tt r s ds (cid:21) = P ( t, T ) exp (cid:20) − (cid:90) Tt σ v b ( v, T ) dW v − (cid:90) Tt σ v b ( v, T ) dv (cid:21) . (A.9) Here, P ( t, T ) ≡ E (cid:104) e − (cid:82) Tt r s ds (cid:12)(cid:12) F t (cid:105) is the time t value of a zero coupon bond maturing at time T , and b ( t, T ) ≡ (cid:82) Tt e − (cid:82) vt a z dz dv .Proof. The integral on the left hand side of (A.9) can be split into − (cid:82) Tt r s ds = − (cid:82) Tt x s ds − (cid:82) Tt φ s ds . We start with integrating the first of these integrals by parts, (cid:90) Tt x s ds = sx s (cid:12)(cid:12) Tt − (cid:90) Tt sdx s =( T − t ) x t + (cid:90) Tt ( T − v )( − a v x v dv + σ v dW v ) . (A.10)We compute x s by evaluating the following integral (cid:90) st d u (cid:16) x u e (cid:82) ut a z dz (cid:17) = (cid:90) st a u e (cid:82) ut a z dz x u du + (cid:90) st e (cid:82) ut a z dz ( − a u x u du + σ u dW u ) , which leads to x s = x t e − (cid:82) st a z dz + (cid:90) st e − (cid:82) su a z dz σ u dW u . (A.11)We plug this into (A.10) to get (cid:90) Tt x s ds = ( T − t ) x t + (cid:90) Tt ( T − v ) (cid:20) − a v (cid:18) x t e − (cid:82) vt a z dz + (cid:90) vt e − (cid:82) vu a z dz σ u dW u (cid:19) dv + σ v dW v (cid:21) . − (cid:90) Tt ( T − v ) a v e − (cid:82) vt a z dz dv =( T − v ) e − (cid:82) vt a z dz (cid:12)(cid:12) Tv = t + (cid:90) Tt e − (cid:82) vt a z dz dv = − ( T − t ) + b ( t, T ) . The second one evaluates − (cid:90) Tt ( T − v ) a v (cid:90) vt e − (cid:82) vu a z dz σ u dW u dv = − (cid:90) Tt (cid:18)(cid:90) vt e (cid:82) ut a z dz σ u dW u (cid:19) d v (cid:18)(cid:90) vt a y ( T − y ) e − (cid:82) yt a z dz dy (cid:19) = − (cid:20)(cid:90) vt e (cid:82) ut a z dz σ u dW u (cid:90) vt a y ( T − y ) e − (cid:82) yt a z dz dy (cid:21) Tv = t + (cid:90) Tt e (cid:82) vt a z dz σ v dW v (cid:90) vt a y ( T − y ) e − (cid:82) yt a z dz dy = − (cid:90) Tt e (cid:82) vt a z dz σ v dW v (cid:90) Tv a y ( T − y ) e − (cid:82) yt a z dz dy = (cid:90) Tt [ − ( T − v ) + b ( v, T )] σ v dW v . This leads to (cid:90) Tt x s ds = x t b ( t, T ) + (cid:90) Tt σ v b ( v, T ) dW v . (A.12)Hence E (cid:20)(cid:90) Tt x u du (cid:21) = x t b ( t, T ) Var (cid:20)(cid:90) Tt x u du (cid:21) = (cid:90) Tt σ v b ( t, T ) dv. Since E (cid:2) e − X (cid:3) = e − µ + σ where X ∼ N ( µ, σ ), we have P ( t, T ) = E (cid:104) e − (cid:82) Tt ( φ s + x s ) ds (cid:12)(cid:12) F t (cid:105) = e − (cid:82) Tt φ s ds − x t b ( t,T )+ (cid:82) Tt σ v b ( v,T ) dv = e (cid:82) Tt x s ds e − (cid:82) Tt r s ds e − x t b ( t,T )+ (cid:82) Tt σ v b ( v,T ) dv . Plugging in (A.12) into this expression gives the desired result.
Acknowledgments
The authors are grateful to Dooheon Lee and Kisun Yoon for numer-ous enlightening discussions and guidance about the theoretical foundations of this paper.The authors would also like to thank Agus Sudjianto for supporting this research, and Vi-jayan Nair for suggestions, feedback, and discussion regarding this work. Any opinions,findings and conclusions or recommendations expressed in this material are those of theauthors and do not necessarily reflect the views of Wells Fargo Bank, N.A., its parentcompany, affiliates and subsidiaries. 34 eferences [1] Bruno Dupire. Pricing with a Smile.
Risk Magazine , pages 18–20, January 1994.[2] Emanuel Derman and Iraj Kani. Riding on a Smile.
Risk Magazine , pages 32–39,February 1994.[3] Marc Atlan. Localizing Volatilities, 2006.[4] Bing Hu. Local Volatility Model with Stochastic Interest Rate. Master’s thesis, YorkUniversity, August 2015.[5] Julien Hok and Shih-Hau Tan. Calibration of Local Volatility Model with StochasticInterest Rates by Efficient Numerical PDE Method, 2018.[6] Griselda Deelstra and Gr´egory Ray´ee. Local Volatility Pricing Models for Long-datedFX Derivatives. Papers 1204.0633, arXiv.org, April 2012.[7] Steven L. Heston. A Closed-Form Solution for Options with Stochastic Volatilitywith Applications to Bond and Currency Options.
The Review of Financial Studies ,6(2):327–343, 04 2015.[8] John C. Cox, Jonathan E. Ingersoll, and Stephen A. Ross. A Theory of the TermStructure of Interest Rates.
Econometrica , 53(2):385–407, 1985.[9] Elias M. Stein and Jeremy C. Stein. Stock Price Distributions with Stochastic Volatil-ity: An Analytic Approach.
Review of Financial Studies , 4:727–752, 1991.[10] Rainer Sch¨obel and Jianwei Zhu. Stochastic Volatility With an Ornstein-UhlenbeckProcess: An Extension.
Review of Finance , 3(1):23–46, 04 1999.[11] G. E. Uhlenbeck and L. S. Ornstein. On the Theory of the Brownian Motion.
Phys.Rev. , 36:823–841, Sep 1930.[12] Bruno Dupire. A unified theory of volatility. In
Discussion paper Paribas CapitalMarkets, reprinted in ”Derivatives Pricing: The Classic Collection”, edited by PeterCarr, 2004 , 1996.[13] Carol Alexander and Leonardo M. Nogueira. Stochastic local volatility. In
Proceedingsof the Second IASTED Int. Conf. on Financial Engineering and Applications , 2004.[14] Alexander Lipton. The vol smile problem.
Risk Magazine , 15(2):61–66, February 2002.[15] Michal Jex, R. C. W. Henderson, and Desheng Wang. Pricing Exotics Under the Smile.In
Pricing Exotics Under the Smile , 1999.[16] Y. Ren, D. Madan, and Qian M. Q. Calibrating and Pricing with Embedded LocalVolatility Models.
Risk , pages 138–143, September 2007.[17] Yuri F. Saporito, Xu Yang, and Jorge P. Zubelli. The Calibration of Stochastic-LocalVolatility Models - An Inverse Problem Perspective, 2017.3518] S.E. Shreve.
Stochastic Calculus for Finance II: Continuous-Time Models . Number v.11 in Springer Finance Textbooks. Springer, 2004.[19] Jim Gatheral.
The Volatility Surface , chapter 1, pages 1–14. John Wiley & Sons, Ltd,2012.[20] Orcan Ogetbil. Extensions of Dupire Formula: Stochastic Interest Rates and StochasticLocal Volatility, 2020.[21] S´ebastien Gurrieri, Masaki Nakabayashi, and Tony Wong. Calibration methods ofHull-White model.
Available at SSRN 1514192 , 2009.[22] Sergei Mikhailov and Ulrich N¨ogel. Heston’s Stochastic Volatility Model Implementa-tion, Calibration and Some Extensions.
Willmott Magazine , 1:74–79, 2005.[23] Daniel Guterding and Wolfram Boenkost. The Heston stochastic volatility model withpiecewise constant parameters - efficient calibration and pricing of window barrieroptions.
Journal of Computational and Applied Mathematics , 343:353 – 362, 2018.[24] William Feller. Two Singular Diffusion Problems.
Annals of Mathematics , 54(1):173–182, 1951.[25] J. A. Nelder and R. Mead. A Simplex Method for Function Minimization.
The Com-puter Journal , 7(4):308–313, 01 1965.[26] Paul Wilmott.
Paul Wilmott on Quantitative Finance . John Wiley & Sons, 2013.[27] D. Brigo and F. Mercurio.