A moment matching method for option pricing under stochastic interest rates
AA moment matching method for option pricing understochastic interest rates
F. Antonelli ∗ , A. Ramponi † , S. Scarlatti ‡ May 29, 2020
Abstract
In this paper we present a simple, but new, approximation methodology forpricing a call option in a Black & Scholes market characterized by stochastic in-terest rates. The method, based on a straightforward Gaussian moment matchingtechnique applied to a conditional Black & Scholes formula, is quite general and itapplies to various models, whether affine or not. To check its accuracy and com-putational time, we implement it for the CIR interest rate model correlated withthe underlying, using the Monte Carlo simulations as a benchmark. The method’sperformance turns out to be quite remarkable, even when compared with analogousresults obtained by the affine approximation technique presented in [9] and by theexpansion formula introduced in [11], as we show in the last section.
Keywords : Option pricing, Stochastic interest rates, Moment matching, Non-affine models, Cox-Ingersoll-Ross model.
Since the appearance of the seminal Black & Scholes/Merton option pricing funda-mental formula, there has been an intensive effort to incorporate in the market modeladditional stochastic factors, such as the volatility and/or the interest rates, the latteralready discussed by Merton himself in [15]. Along the years, a huge field of researchdeveloped, leading to a very rich literature on stochastic volatility models, while fewerpapers aimed at the inclusion of a dynamic term structure into the valuation of deriva-tives, e.g. [17], [1], [20], [19], [11], [7], [18].Nowadays, the improvement in the performances of option pricing formulas obtainedby adding these risk factors is widely recognized in the empirical literature (see e.g [2],[3]), indeed in ([12]) the author remarked that even including solely stochastic interestrates in the model does affect the pricing formula, especially for longer-dated options,in a noticeable manner. ∗ DISIM, University of L’Aquila, [email protected] † Dept. Economics & Finance, University of Roma - Tor Vergata, [email protected] ‡ Dept. Enterprise Engineering, University of Roma - Tor Vergata, [email protected] a r X i v : . [ q -f i n . C P ] M a y f course this generalization implies a higher degree of mathematical complexityand the search for efficient pricing techniques, able to provide accurate answers in ashort computational time (as opposed to Monte Carlo methods) has been relentless,even more so in modern quantitative finance where a huge amount of data allows toconsider strategies that call for on real-time model calibration. Hence, computationalefficiency has become one of the primary concerns of risk managers and this requirementessentially restricted the choice of models to the affine class (see [7]).Indeed, when the interest rates are modeled in a Gaussian processes framework, as inthe very popular Hull-White / Vasicek models, even analytical prices can be obtained.These models are appropriate for modeling periods that admit positive probability ofnegative rates, such as the current one, but this feature becomes a drawback in usualperiods of positive rates. The most popular model used to avoid this drawback isthe Cox-Ingersoll-Ross (CIR) one, which guarantees the rate’s strict positivity underFeller’s condition. Its popularity comes from the fact that falls into the so called affinemodels, that can exploit a very efficient and fast Fourier transform technique to pricethe bonds.Unfortunately, the affinity of the model is lost when the interest rate is coupled,with correlation, with a risky asset’s dynamics, making the search for efficient approx-imations of risk-neutral pricing formulas very challenging.Here we present a simple, but new, approximation methodology for pricing a Eu-ropean call option in a market model given by a linear diffusion dynamics for theunderlying (a Black & Scholes (BS) framework) coupled with a stochastic short termrisk-free rate. The problem is a classical one and the novelty lies on the fact that wepropose a quite straightforward moment matching (MM) technique, easy to implementand leading to very efficient approximations.In building our procedure a few issues have to be addressed and we first provide,by appropriate conditioning, a representation formula for the claim’s price in terms ofthe BS formula, then we exploit a Gaussian approximation by properly matching thefirst two moments of the involved random variables, that allows to use the propertiesof the Normal cumulative distribution function (c.d.f.) (see Lemma (1)). When apply-ing the method to the affine models, we also employ a change-of-numeraire technique(introducing the T -forward measure as in [5]) to partially disentangle the contributionsdue to the underlying and to the interest rate, exploiting the explicit expressions of thebond’s price in an affine framework. To keep computations as simple as possible, anytime a quantity is computable, it is stored and treated as a constant in the sequel. Thisleads to an efficient mixed use of the risk free probability and the T -forward measureto evaluate the separate quantities.The paper is organized as follows. In Section 2 we derive a representation formulafor the call option’s price in Black & Scholes market with stochastic risk-free short rates,while in Section 3 the Moment Matching method is fully described. Finally, in Section 4we restrict to the affine models and we apply our technique with a CIR interest rate. Inthe same Section, we briefly introduce other two techniques, the affine approximation,inspired by Grzelak and Oosterlee [9] and the expansion method proposed by Kim andKunimoto [11] alternative to prices obtained by Monte Carlo simulations. Hence we run2 numerical study comparing those methods with ours, using Monte Carlo evaluationas a benchmark. The underlying problem we are concerned with is the pricing of a European call option,whose payoff is given by the function f ( x ) = (e x − e κ ) + for some κ ∈ R , when stochasticinterest rates come into play.Thus, given a finite time interval [0 , T ] and a complete probability space (Ω , F , Q ),endowed of a filtration {F t } { t ∈ [0 ,T ] } satisfying the “usual hypotheses” (see [16]), themarket model is defined by the log-price of a risky asset and a risk-free interest ( X t , r t ),whose joint dynamic for any initial condition ( t, x, r ) ∈ [0 , T ] × R × R and ∀ s ∈ [ t, T ] isgiven by (cid:40) X s = X t + (cid:82) st ( r v − σ ) dv + σ (cid:104) ρ ( B s − B t ) + (cid:112) − ρ ( B s − B t ) (cid:105) , X t = xr s = r t + (cid:82) st µ ( v, r v ) dv + (cid:82) st η ( v, r v ) dB v , r t = r, (1)where ( B , B ) is a two dimensional standard Brownian motion and ρ ∈ ( − , µ ( · , · ) and η ( · , · ) are in a class thatensures the existence and uniqueness of a strong solution of (1) (see e.g. [10]) and that Q is some risk neutral probability selected by the market.Under these assumptions, the pair ( X t , r t ) is Markovian, whence the arbitrage-freeoption’s price is a deterministic function of the state variables, given by u ( t, x, r, T ; ρ ) = E (e − (cid:82) Tt r s (e X T ( ρ ) − e κ ) + ds | X t = x, r t = r ) , (2)provided that the coefficients µ and η are chosen to guarantee the exponential integra-bility of X T and (cid:82) T | r s | ds . Here we wrote X T ( ρ ), to stress the prices’ dependence onthe correlation parameter.If u ( t, x, r, T ; ρ ) is regular enough in t, x, r , Feymann-Kac’s formula implies that itis a classical solution of the following two-dimensional parabolic problem (cid:40) ∂u∂t + L ρ u = 0 u ( T, x, r, T ; ρ ) = (e x − e κ ) + , (3)where L ρ = L + A , with L := (cid:18) σ ∂ ∂x + ( r − σ ∂∂x − r (cid:19) + (cid:18) η ( t, r )2 ∂ ∂r + µ ( t, r ) ∂∂r (cid:19) (4) A := ρση ( t, r ) ∂ ∂x∂r . (5)In what follows to keep notation easy, we take t = 0 and we omit the dependence on t in the pricing function. The general case may be readily obtained substituting T inthe final formulas with the time to maturity T − t .3y conditioning internally with respect to F T = σ ( { B s : 0 ≤ s ≤ T } ), we have u ( x, r, T ; ρ ) = E (cid:16) e − (cid:82) T r s ds (e X T ( ρ ) − e κ ) + (cid:17) = E (cid:16) e − (cid:82) T r s ds E (cid:0) (e X T ( ρ ) − e κ ) + |F T (cid:1)(cid:17) . (6)But X T ( ρ ) (cid:12)(cid:12) F T ∼ N ( M T , Σ T ), where M T = x + (cid:90) T ( r s − σ ds + σρB T , and Σ T = σ (1 − ρ ) T. so we obtain E (cid:16) (e X T ( ρ ) − e κ ) + |F T (cid:17) =e M T + Σ T N (cid:18) M T − κ + Σ T Σ T (cid:19) − e κ N (cid:18) M T − κ Σ (cid:19) =e x + (cid:82) T ( r s − σ ) ds + σρB T + σ (1 − ρ ) T N ( d ( ρ )) − e κ N ( d ( ρ )) , where we define d ( ρ ) = x − κ + (cid:82) T r s ds + σρB T + σ T − σ ρ Tσ (cid:112) − ρ √ T (7) d ( ρ ) = x − κ + (cid:82) T r s ds + σρB T − σ Tσ (cid:112) − ρ √ T (8)and N denotes the cumulative distribution function of the standard Gaussian. It isconvenient to introduce the following notationsΛ T = (cid:90) T r s ds, β ( T, ρ ) = ρ (1 − ρ ) / √ T , γ ( T, ρ ) = 1 σ (1 − ρ ) / √ Tα ( x, T, ρ ) = x − κ + σ T − σ ρ Tσ (1 − ρ ) / √ T , α ( x, T, ρ ) = x − κ − σ Tσ (1 − ρ ) / √ T , so that d i ( x, T, ρ ) = α i ( x, T, ρ ) + β ( T, ρ ) B T + γ ( T, ρ )Λ T , i = 1 , . (9)Setting S T = e − Λ T , we can finally write u ( x, r, T ; ρ ) = e x e − σ ρ T E (cid:16) e σρB T N (cid:0) d ( ρ ) (cid:1)(cid:17) − e κ E (cid:0) S T N (cid:0) d ( ρ ) (cid:1)(cid:1) . (10)In the forthcoming section we shall introduce the moment matching approximationprocedure. The main idea of this section is to replace the r.v.’s d i ( ρ ), i = 1 ,
2, defined by (9), withGaussian r.v.’s D i ( ρ ) matching the first and second moments of d i ( ρ ).4e define D i ( ρ ) := α i ( x, T, ρ ) + ˆ β ( T, ρ ) B T + γ ( T, ρ ) E (Λ T ) , i = 1 , , consequently E (cid:0) D i ( ρ ) (cid:1) = α i ( x, T, ρ ) + γ ( T, ρ ) E (Λ T ) = E (cid:0) d i ( ρ ) (cid:1) (11)and the new coefficient ˆ β > (cid:0) D i ( ρ ) (cid:1) = T ˆ β ( T, ρ ) = var( d ( ρ )) = var( d ( ρ ))= β ( T, ρ ) T + γ ( T, ρ )var(Λ T ) + 2 β ( T, ρ ) γ ( T, ρ ) E ( B T Λ T ) (12)with var(Λ T ) = E (cid:18)(cid:16) (cid:90) T r s ds (cid:17) (cid:19) − (cid:104) E (cid:18)(cid:90) T r s ds (cid:19) (cid:105) , (13) E ( B T Λ T ) = E (cid:18) B T (cid:90) T r s ds (cid:19) . (14)The moment matching method with Gaussian r.v’s may be motivated by looking atthe empirical distributional properties of the random variables d i in some well-knownrate models: see as examples Figs (1), (2) and (3).We finally introduce a call price approximation u appr ( x, r, T ; ρ ) :=e x e − σ ρ T E (cid:16) e σρB T N (cid:0) D ( ρ ) (cid:1)(cid:17) − e κ E (cid:0) S T N (cid:0) D ( ρ ) (cid:1)(cid:1) =:e x e − σ ρ T F ( ρ ) − e κ G ( ρ ) . (15)The function F can be evaluated in closed form by means of the following Lemma 1
Let p ∈ R and X ∼ N ( µ, ν ) , ( µ, ν ) ∈ R × R + , then E (e pX N ( X )) = e pµ + ( pν )22 N (cid:18) µ + pν √ ν (cid:19) . Proof:
See [21] for p = 0, the general case follows by a “completing the squares”argument. (cid:3) Since B T = (cid:2) D ( ρ ) − α ( x, T, ρ ) − γ ( T, ρ ) E (Λ T ) (cid:3) ˆ β ( T, ρ ) − , we may rewrite F as F ( ρ ) = E (cid:16) e σρ ( D ( ρ ) − α ( x,T,ρ ) − γ ( T,ρ ) E (Λ T )) ˆ β ( T,ρ ) − N ( D ( ρ )) (cid:17) =e − σρ [ α ( x,T,ρ )+ γ ( T,ρ ) E (Λ T )] ˆ β ( T,ρ ) − E (cid:16) e σρD ( ρ ) ˆ β ( T,ρ ) − N ( D ( ρ )) (cid:17) =e − σρ [ α ( x,T,ρ )+ γ ( T,ρ ) E (Λ T )] ˆ β ( T,ρ ) − e σρ E ( D ( ρ )) ˆ β ( T,ρ ) − + ˆ β ( T,ρ ) − σ ρ D ρ ))2 × N (cid:16) E ( D ( ρ )) + σρ var( D ( ρ )) ˆ β ( T, ρ ) − (cid:112) D ( ρ )) (cid:17) . F ( ρ ) = e σ ρ T N α ( x, T, ρ ) + σρ ˆ β ( T, ρ ) T + γ ( T, ρ ) E (Λ T ) (cid:113) β ( T, ρ ) T . (16)If E (Λ T ) and (13), (14) can be computed, then F is totally explicit. From now on, wedenote λ ( T ) := E (Λ T ) to point out this is a known constant.On the contrary, the function G cannot be evaluated in such a straightforwardmanner, as it involves a detailed knowledge of the joint distribution of Λ T and B T andnot only of their moments and covariance. So, to represent G , we suggest applying a achange-of-numeraire technique that allows us to exploit the bond pricing theory.Let us define P ( s, T ) := E (cid:16) (e − (cid:82) Ts r v dv |F s (cid:17) , (17)the Zero Coupon Bond price. Again, since r. is a Markov process, P ( s, T ) is a determin-istic function of the state variable, say g ( s, r s ), which we assume to be C , ([0 , T ] × R + ).For 0 ≤ s ≤ T , we define the F s -martingale (we remark that is a true martingale thanksto the exponential integrability of Λ T ) L s = E (e − (cid:82) T r v dv |F s ) P (0 , T ) = S s P ( s, T ) P (0 , T ) = S s g ( s, r s ) g (0 , r ) , r = r. (18)By applying Itˆo’s formula, we have the dynamic of LdL s = S s g (0 , r ) (cid:104) ∂g∂t ( s, r s ) + 12 η ( s, r s ) ∂ g∂r ( s, r s ) + µ ( s, r s ) ∂g∂r ( s, r s ) − r s g ( s, r s ) (cid:105) ds + S s g (0 , r ) η ( s, r s ) ∂g∂r ( s, r s ) dB s = S s g (0 , r ) η ( s, r s ) ∂g∂r ( s, r s ) dB s , L = 1and we may define the T -forward measure on every A ∈ F by Q T ( A ) := E ( L T A ) (see[4] for the method and [6] for a similar application). Under Q T , we get G ( ρ ) = E ( S T N ( D ( ρ ))) = P (0 , T ) E Q T ( N ( D ( ρ ))) , (19)and by Girsanov theorem, by setting ξ s := (cid:90) s η ( v, r v ) g ( v, r v ) ∂g∂r ( v, r v ) dv, we have that the process ˜ B s := B s − ξ s is a Q T − Brownian motion. When choosing an in-terest rate model that allows an explicit expression of the bond’s price, E Q T ( N ( D ( ρ )))will be the last quantity to compute. Under Q T , D ( ρ ) has the expression D ( ρ ) = α ( x, T, ρ ) + ξ T ˆ β ( T, ρ ) + ˆ β ( T, ρ ) ˜ B T + γ ( T, ρ ) λ ( T ) , whence its distribution is no longer known.6o compute the final expectation, we replace D ( ρ ) by the r.v.¯ D ( ρ ) := α ( x, T, ρ ) + E ( ξ T ) ˆ β ( T, ρ ) + ˆ β ( T, ρ ) ˜ B T + γ ( T, ρ ) λ ( T ) , where we are taking (cid:15) ( T ) := E ( ξ T ) under the original probability Q , so that ¯ D ( ρ ) is aGaussian r.v. and we may apply Lemma 1 once again to obtain E Q T ( N ( ¯ D ( ρ ))) = N E Q T ( ¯ D ( ρ )) (cid:113) Q T ( ¯ D ( ρ )) , with E Q T ( ¯ D ( ρ )) = α ( x, T, ρ ) + (cid:15) ( T ) ˆ β ( T, ρ ) + γ ( T, ρ ) λ ( T ) , var Q T ( ¯ D ( ρ )) = ˆ β ( T, ρ ) T. Hence we shall denote by ¯ G ( ρ ) := P (0 , T ) E Q T ( N ( ¯ D ( ρ )))the approximation of G ( ρ ) and we may define the final approximation of the call optionprice u ( x, r, T ; ρ ) as¯ u ( x, r, T ; ρ ) := e x − σ ρ T F ( ρ ) − e κ ¯ G ( ρ )= e x N α ( x, T, ρ ) + σρ ˆ β ( T, ρ ) T + γ ( T, ρ ) λ ( T ) (cid:113) β ( T, ρ ) T (20) − e κ P (0 , T ) N α ( x, T, ρ ) + (cid:15) ( T ) ˆ β ( T, ρ ) + γ ( T, ρ ) λ ( T ) (cid:113) β ( T, ρ ) T . As a conclusion, we summarize the key requirements to make the approximation (20)explicitly computable and hopefully efficient1. the distributions of d i ( ρ ) , i = 1 , P ( t, T ) should be theoretically computable. Moreover one canexploit the observed (today) bond price for P (0 , T ) in (18) and for calibrationpurposes;3. the quantities E (Λ T ), var(Λ T ) and E (Λ T B T ) and/or their approximations, shouldbe easily computable;4. the change of numeraire technique (Girsanov’s theorem) should be applicable.The performance of this approximation needs to be compared with Monte-Carlosimulated prices and then with other methods present in the literature. This will bedone in the next section. 7 Numerics and comparison with other methodologies
In this section we employ an affine model for the interest rate. This choice provides anexplicit expression for the the ZCB’s price (17). So our market model is given by X s = X t + (cid:90) st ( r v − σ dv + σ (cid:104) ρ ( B s − B t ) + (cid:112) − ρ ( B s − B t ) (cid:105) , X t = xr s = r t + (cid:90) st [ a ( v ) r v + b ( v )] dv + (cid:90) st [ c ( v ) r v + d ( v )] / dB v , r t = r, (21)where a, b, c, d : [0 , T ] −→ R are bounded functions. In this framework, we have for r t = r P ( t, T ) = g ( t, r ) = A ( t, T )e − rB ( t,T ) , for suitable deterministic functions A ( · , T ) and B ( · , T ). Two very classical models fallinto this setting(Vasicek) a ( v ) = − γ, b ( v ) = γθ, c ( v ) = 0 , d ( v ) = η (CIR) a ( v ) = − γ, b ( v ) = γθ, c ( v ) = η , d ( v ) = 0 γ, θ > , for which A ( t, T ) and B ( t, T ) are explicitly known ([5]), the same being true also forthe Hull-White / Vasicek and Hull-White / CIR models, considering time dependentcoefficients.The functions A and B are usually characterized by the solution of a Riccati systemof ODE’s. Unfortunately, when in presence of correlation, the same procedure cannotbe applied to the pair ( X., r. ), since its diffusion matrix σ ( v, x, r ) σ ( v, x, r ) T = (cid:18) σ ρσ [ c ( v ) r + d ( v )] / ρσ [ c ( v ) r + d ( v )] / c ( v ) r + d ( v ) (cid:19) (22)may haves entries which are non-linear in the state variables, so that the joint diffusionis no longer affine, as it happens for the CIR model. Hence, in this context it makessense to apply the approximation procedure presented in the previous section. Asbefore we consider t = 0.In this case (see e.g. [5]), setting δ = (cid:112) γ + 2 η , we have A (0 , T ) = e γθη δ e γ + δT δ − γ + ( δ + γ )e δT , B (0 , T ) = 2(e δT − δ − γ + ( δ + γ )e δT , and let us proceed to the computation of E (Λ T ), var(Λ T ) and E (Λ T B T ).1. Computation of E (Λ T ). It is straightforward to see E (Λ T ) = (cid:90) T E ( r s ) ds = (cid:90) T (cid:2) ( r − θ )e − γs + θ (cid:3) ds = θT + ( r − θ ) 1 − e − γT γ . Computation of var(Λ T ) . Taking into account the first point, we only have tocompute the second moment E (cid:18)(cid:16)(cid:90) T r s ds (cid:17) (cid:19) = E (cid:18)(cid:90) T (cid:90) T r s r v dsdv (cid:19) = (cid:90) T (cid:90) T E ( r s r v ) dsdv = (cid:90) T (cid:90) t E ( r s r v ) dsdv + (cid:90) T (cid:90) Tt E ( r s r v ) dsdv = (cid:90) T (cid:90) t E ( r s r v ) dsdv + (cid:90) T (cid:90) s E ( r s r v ) dvds = 2 (cid:90) T (cid:90) s E ( r s r v ) dvds. By the independence of the increments of the process r , for v < s we have E ( r s r v ) = E (cid:0) ( r s − r v ) r v + r v (cid:1) = E ( r s − r v ) E ( r v ) + E ( r v )= E (cid:16) θ ( s − v ) + ( r v − θ ) 1 − e − γ ( s − v ) γ (cid:17) E ( r v ) + E ( r v )= θ (cid:104) ( s − v ) − − e − γ ( s − v ) γ (cid:105) E ( r v ) + 1 − e − γ ( s − v ) γ [ E ( r v )] + var( r v ) + [ E ( r v )] Since var( r s ) = r η γ (e − γs − e − γs ) + θη γ (1 − e − γs ) , all the integrals appearing inthe second moment can be calculated analytically.3. Computation of E ( B T Λ T ). By Itˆo’s integration-by-parts formula, we get E ( B T Λ T ) = (cid:90) T E ( B s r s ) ds and again by integration by parts we have B s r s = (cid:90) s B v dr v + (cid:90) s r v dB v + (cid:104) B , r (cid:105) s = (cid:90) s B v γ ( θ − r v ) dt + η (cid:90) s B t √ r v dB v + (cid:90) s r v dB v + η (cid:90) s √ r v dv, so that E ( B s r s ) = − γ (cid:90) s E ( B v r v ) dv + η (cid:90) s E ( √ r v ) dv. Solving this linear ODE for h ( s ) := E ( B s r s ), since h (0) = 0, we obtain E ( B s r s ) = η (cid:90) s e − γ ( s − v ) E ( √ r v ) dv, E ( B T Λ T ) = η (cid:90) T (cid:90) s e − γ ( s − v ) E ( √ r v ) dvds. Thus the final crucial point is computing E ( √ r v ), which is rather delicate (see[8]). No explicit expression can be provided and we employ the approximationproposed in [9], that we are going to present in the next subsection E ( √ r v ) ≈ a + b e − cv (23)9here the parameters a, b and c are obtained by an ad hoc matching procedure,which proved to be numerically very efficient.Finally, given the above three points, ˆ β ( T, ρ ) is easily computed from (12).As a last step, we have to approximate E ( ξ T ), which is readily done, given the lastremark, since E ( ξ T ) = E (cid:16) − η (cid:90) T B ( s, T ) √ r s ds (cid:17) = E (cid:16) − η (cid:90) T B ( s, T ) E (cid:0) √ r s (cid:1) ds (cid:17) ≈ − η (cid:90) T B ( s, T )( a + b e − cs ) ds. Here and in the next subsection, for completeness, we briefly describe the two approx-imation techniques, we are going to compare with.The GO approximation consists simply in modifying the A operator given in (5) byreplacing the state variable in the coefficient with a constant, namely we define A GO u ( s, x, r ) := ρσ E ( η ( s, r t s ) ∂ u∂x∂r . In the case of the CIR model, η is time-homogeneous and this operator becomes A GO u ( s, x, r ) := ρση E ( √ r s )) ∂ u∂x∂r ≈ ρση ( a + b e − cs ) ∂ u∂x∂r Once this replacement has been made then the Fourier transform methods apply, henceit is possible to compute approximated prices of the call option. We shall denote thisapproximation by u GO ( t, x, r, T ; ρ ). To evaluate the accuracy of this approximation acomparison with the prices of the (non-affine) true model, obtained by MC simulations,must be performed.Once again we specialize the formulas for t = 0 for a direct comparison with ourresults, so X = x and r = r .The discounted transform, for ζ ∈ C , (see [7]) for the affine approximation is φ ( ζ, x, r, T ) := E (cid:16) e − (cid:82) T r s ds e ζX T (cid:17) = e A ( ζ,T )+ B ( ζ,T ) x + C ( ζ,T ) r , where the functions A, B, C satisfy a system of solvable ODE’s, that give B ( ζ, T ) = ζ,C ( ζ, T ) = 1 − e − dT η (1 − g e − dT ) , d = (cid:112) γ + 2 η (1 − ζ ) , g = γ − dγ + d ,A ( ζ, T ) = − σ T ζ (1 + ζ ) + γ − dη (cid:90) T (cid:104) γθη + ρσζr sqs (cid:105) − e − ds (1 − g e − ds ) ds, where r sqs = E ( √ r s ) is approximated as in (23).10inally, by L´evy inversion formula as in [7] or Fourier inversion as in [13], one getsan integral representation for the price function: in our implementation we use theFourier inversion u GO ( x, r, T ; ρ ) = e νγ π (cid:90) + ∞ R (cid:18) e − i ζγ ν − ν − ζ + i ζ (1 − ν ) φ ( ζ, x, r, T ) (cid:19) dζ, (24)where ν < R ( z ) is the real part for z ∈ C . Kim and Kunimoto, in [11], consider a Taylor expansion of the process r s in powers of η around η = 0. Considering the first order polynomial and setting ϕ ( s ) = r exp( − γs ) + θ (1 − exp( − γs )), they obtain r s = ϕ ( s ) + η (cid:90) s e − γ ( s − v ) (cid:112) ϕ ( v )( ρdB v + (cid:112) − ρ dB v ) + o ( η ) . (25)Inserting the approximation (25) in the evaluation formula for the call option, aftersome manipulations one can approximate the option’s price as u KK ( x, r, T ; ρ ) =e x N ( d ) − e κ − (cid:82) T ϕ ( s ) ds N ( d )+ ηC (cid:104) d e x N (cid:48) ( d ) − d e κ − (cid:82) T ϕ ( s ) ds N (cid:48) ( d ) (cid:105) (26)where C = − ρσT √ θ (cid:2) (1 + 2e γT ) √ r − γ K (cid:3) + (cid:2) r − θ (1 + 2 e γT ) (cid:3) λ K γT γ √ θ ,d = x − κ + θT + ( r − θ )(1 − e − γT ) /γ + σ T / √ σ T , d = d − σ √ T , being γ K = e γT/ (cid:112) r − θ (1 − e γT ) and λ K = log (cid:16) ( √ r + √ θ ) r − θ (1 − γT )+2 γ K √ θ ) (cid:17) . We compare the results of the different approximations with the benchmark MonteCarlo method, applied to the price (10). In particular this means that we only haveto simulate the rate process to get samples from d ( ρ ) and d ( ρ ). The simulationwas implemented by means of the Euler discretization with full truncation algorithm(see [14]). In our numerical experiments we generated M = 10 sample paths with atime step discretization equal to 10 − for all the maturities. All the algorithms wereimplemented in MatLab (R2019b) and ran on an Intel Core i7 2.40GHZ with 8GBRAM, by using the available building-in functions, in particular for the computationof all the integrals involved. The average time to compute one price was (in secs) 32 . .
055 (GO), 0 .
005 (KK) and 0 .
009 (MM).11e chose different set of parameters ( κ, θ, η ) and volatility scenarios: a low volatility σ L = 0 . σ L = 0 .
4; hence we varied the correlation ρ , the ratevolatility η and the maturity of the contract T . The initial price of the underlyingwas set to 100 as well as the strike price K . Numerical results are summarized inTables (1) - (8). At least in the CIR model, the numerical results show that the MMmethod produces the best approximations with respect to the benchmark Monte Carloevaluation in most scenarios. References [1] Amin, K. and R. Jarrow,
Pricing options on risky assets in a stochastic interestrate economy , Mathematical Finance 2, (1992) 217 – 237.[2] Bakshi, G., Cao, C., Chen, Z.,
Empirical performance of alternative option pricingmodels , The Journal of finance, 52(5), (1997) 2003–2049.[3] Bakshi, G., Cao, C., Chen, Z.
Pricing and hedging long-term options , Journal ofeconometrics, 94(1-2), (2000) 277–318.[4] T. Bjork,
Arbitrage Theory in Continuous Time , Oxford University Press 2009.[5] D. Brigo, F. Mercurio,
Interest Rates Models , Springer 2013.[6] D. Brigo, F. Vrins,
Disentangling wrong-way risk: pricing credit valuation ad-justment via change of measures , European Journal of Operational Research 269(2018) 1154–1164.[7] D. Duffie, J. Pan, K. J. Singleton,
Tranform analysis and asset pricing for affinejum-diffusions , Econometrica, Vol. 68, No. 6, November, (2000), 1343–1376.[8] D. Dufresne,
The integrated square root process , Research Paper, University ofMelbourne, (2001).[9] L.A. Grzelak, C.W. Oosterlee,
On the Heston model with stochastic interest rates ,SIAM Journal on Financial Mathematics 2 (1) (2011) 255-286.[10] I. Karatzas, S: Shreve
Brownian Motion and Stochastic Calculus , Graduate Textsin Mathematics 113, Springer-Verlag New Yprk (1998)..[11] Y. J. Kim, N. Kunimoto,
Pricing Options under Stochastic Interest Rates: A NewApproach.
Asia-Pacific Financial Markets, 6, 49–70 (1999).[12] Y.J. Kim,
Option Pricing under Stochastic Interest Rates: An Empirical Investi-gation , Asia-Pacific Financial Markets 9, (2002) 23 – 44.[13] R. W. Lee,
Option Pricing by Transform Methods: Extensions, Unification, andError Control , Journal of Computational Finance, 7, 51–86, (2004).1214] R. Lord, R. Koekkoek, D. Van DijK,
A comparison of biased simulation schemesfor the stochastic volatility models , Quantitative Finance, 10 (2) (2010) 177–194.[15] Merton, R.,
The theory of rational option pricing , Bell J. Econom. Managt Sci. 4,(1973) 14 – 183.[16] P. Protter,
Stochastic Integration and Differential Equations , Stochastic Modellingand Applied Probability, 21, Springer-Verlag Berlin Heidelberg (2005).[17] Rabinovitch, R.
Pricing stock and bond options when the default-free rate isstochatic , J. Financ. Quantitat. Anal. 24, (1989) 447 – 457.[18] A. Ramponi, S. Scarlatti, S.
Option pricing in a hidden Markov model of theshort rate with application to risky debt evaluation , Int. J. Risk Assessment andManagement, Vol. 11, (1/2), (2009) 88–103.[19] Rindell, K.
Pricing of index options when interest rates are stochastic: an empiricaltest , J. Banking Finance 19, (1995) 785 – 802.[20] Shimko D.C., Tejima N., Van Deventer D.
The pricing of risky debts when interestrates are stochastic , Journal of Fixed Income Vol. 3, No.2, (1993), pp. 58–65.[21] S. Zacks,
Parametric Statistical Inference , Pergamon Press 1981.13 -5 0 501000200030004000 d -5 0 501000200030004000-5 0 5 Standard Normal Quantiles -10-50510 Q uan t il e s o f I npu t S a m p l e qq-plot d -5 0 5 Standard Normal Quantiles -10-50510 Q uan t il e s o f I npu t S a m p l e qq-plot d d -5 0 501000200030004000 d -5 0 501000200030004000-5 0 5 Standard Normal Quantiles -10-50510 Q uan t il e s o f I npu t S a m p l e qq-plot d -5 0 5 Standard Normal Quantiles -10-50510 Q uan t il e s o f I npu t S a m p l e qq-plot d Figure 1:
The histograms of d and d for ρ = 0 . , T = 1 (left) and T = 5 (right), incomparison with the standard normal law (in red) and related qq-plot, CIR dynamic: dr t = κ ( θ − r t ) dt + η √ r t dB t . d -1 0 101000200030004000 d -1 -0.5 0 0.5 101000200030004000-5 0 5 Standard Normal Quantiles -1012 Q uan t il e s o f I npu t S a m p l e qq-plot d -5 0 5 Standard Normal Quantiles -2-1012 Q uan t il e s o f I npu t S a m p l e qq-plot d d d Standard Normal Quantiles Q uan t il e s o f I npu t S a m p l e qq-plot d -5 0 5 Standard Normal Quantiles Q uan t il e s o f I npu t S a m p l e qq-plot d Figure 2:
The histograms of d and d for ρ = 0 . , T = 1 (left) and T = 5 (right),in comparison with the standard normal law (in red) and related qq-plot, ExponentialVasicek dynamic: dr t = r t ( θ − a ln( r t )) dt + ηr t dB t . -1 -0.5 0 0.5 101000200030004000 d -1 -0.5 0 0.5 101000200030004000-5 0 5 Standard Normal Quantiles -1012 Q uan t il e s o f I npu t S a m p l e qq-plot d -5 0 5 Standard Normal Quantiles -2-101 Q uan t il e s o f I npu t S a m p l e qq-plot d d d Standard Normal Quantiles -10123 Q uan t il e s o f I npu t S a m p l e qq-plot d -5 0 5 Standard Normal Quantiles -1012 Q uan t il e s o f I npu t S a m p l e qq-plot d Figure 3:
The histograms of d and d for ρ = 0 . , T = 1 (left) and T = 5 (right), incomparison with the standard normal law (in red) and related qq-plot, Dothan dynamic: dr t = ar t dt + ηr t dB t . ρ -0.9 -0.6 -0.3 0.0 0.3 0.6 0.9PricesMC 8.1543 8.1799 8.2055 8.2314 8.2574 8.2832 8.3085(0.0225) (0.0137) (0.0064) (0.0003) (0.0069) (0.0142) (0.0230)GO 8.1192 8.1568 8.1943 8.2315 8.2686 8.3055 8.3423KK 8.1361 8.1677 8.1993 8.2309 8.2625 8.2941 8.3258MM 8.146 8.1745 8.2029 8.2313 8.2595 8.2877 8.3157ErrorsGO 0.0351 0.0231 0.0113 -0.0001 -0.0113 -0.0223 -0.0338KK 0.0182 0.0121 0.0062 0.0005 -0.0052 -0.0109 -0.0172MM 0.0083 0.0053 0.0026 0.0001 -0.0022 -0.0045 -0.0072Rel. Err.GO 0.0043 0.0028 0.0013 1.4e-05 0.0014 0.0027 0.0041KK 0.0022 0.0015 0.0008 5.9e-05 0.0006 0.0013 0.0021MM 0.0010 0.00065 0.0003 1.7e-05 0.0003 0.0005 0.0009Table 1: Results of the approximations for the parameters κ = 0 . , θ = 0 . , η = 0 . , r = 0 . and σ L . The time to maturity is T = 1 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. -0.9 -0.6 -0.3 0.0 0.3 0.6 0.9PricesMC 19.8443 20.1287 20.4125 20.6936 20.9705 21.2425 21.5086(0.0595) (0.0351) (0.0153) (0.0026) (0.0202) (0.0404) (0.0649)GO 19.6375 19.9974 20.3492 20.6936 21.0308 21.3614 21.6856KK 19.7487 20.0582 20.3678 20.6773 20.9869 21.2964 21.606MM 19.7747 20.085 20.3892 20.6875 20.981 21.269 21.5522ErrorsGO 0.2067 0.1313 0.0632 5.2e-07 -0.0603 -0.1189 -0.1769KK 0.0956 0.0705 0.0447 0.0162 -0.0164 -0.0540 -0.0973MM 0.0695 0.0436 0.0232 0.0061 -0.0104 -0.0265 -0.0435Rel. Err.GO 0.0104 0.0065 0.0031 2.5e-08 0.0029 0.0056 0.0082KK 0.0048 0.0035 0.0022 0.0008 0.0008 0.0025 0.0045MM 0.0035 0.0022 0.0011 0.0003 0.0005 0.0012 0.0020Table 2: Results of the approximations for the parameters κ = 0 . , θ = 0 . , η = 0 . , r = 0 . and σ L . The time to maturity is T = 5 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. ρ -0.9 -0.6 -0.3 0.0 0.3 0.6 0.9PricesMC 16.0337 16.0533 16.073 16.0933 16.1141 16.1351 16.156(0.0504) (0.0301) (0.0139) (0.0002) (0.0144) (0.0306) (0.0509)GO 15.9831 16.0199 16.0567 16.0934 16.1300 16.1665 16.2030KK 15.9997 16.0309 16.062 16.0932 16.1243 16.1555 16.1866MM 16.0094 16.0374 16.0654 16.0933 16.1211 16.1489 16.1767ErrorsGO 0.0506 0.0333 0.0162 -0.0001 -0.0159 -0.0314 -0.0469KK 0.0339 0.0224 0.0109 0.0001 -0.0102 -0.0203 -0.0306MM 0.0242 0.0159 0.0076 1.8e-05 -0.0071 -0.0138 -0.0207Rel. Err.GO 0.0032 0.0021 0.0010 7.0e-06 0.0010 0.0019 0.0029KK 0.0021 0.0013 0.0007 6.9e-06 0.0006 0.0012 0.0019MM 0.0015 0.0009 0.0005 1.1e-06 0.0004 0.0009 0.0013Table 3: Results of the approximations for the parameters κ = 0 . , θ = 0 . , η = 0 . , r = 0 . and σ H . The time to maturity is T = 1 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. -0.9 -0.6 -0.3 0.0 0.3 0.6 0.9PricesMC 36.1379 36.3566 36.5912 36.8358 37.0875 37.3439 37.6015(0.0405) (0.0102) (0.0141) (0.0008) (0.0163) (0.0324) (0.0571)GO 35.8574 36.1877 36.5138 36.8358 37.154 37.4683 37.7789KK 35.9641 36.2539 36.5437 36.8335 37.1233 37.4132 37.703MM 35.9876 36.2725 36.5543 36.8329 37.1089 37.3819 37.6520ErrorsGO 0.2805 0.1690 0.0774 2.2e-06 -0.0664 -0.1244 -0.1773KK 0.1739 0.1028 0.04748 0.0023 -0.0358 -0.0693 -0.1014MM 0.1504 0.0842 0.0368 0.0029 -0.0214 -0.0379 -0.0504Rel. Err.GO 0.0078 0.0045 0.0021 6.1e-08 0.0018 0.0033 0.0047KK 0.0048 0.0028 0.0013 6.2e-05 0.0010 0.0018 0.0027MM 0.0041 0.0023 0.0010 7.8e-05 0.0006 0.0010 0.0013Table 4: Results of the approximations for the parameters κ = 0 . , θ = 0 . , η = 0 . , r = 0 . and σ H . The time to maturity is T = 5 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. η Results of the approximations for the parameters κ = 0 . , θ = 0 . , ρ = 0 . , r = 0 . and σ L . The time to maturity is T = 1 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. Results of the approximations for the parameters κ = 0 . , θ = 0 . , ρ = 0 . , r = 0 . and σ H . The time to maturity is T = 1 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. η Results of the approximations for the parameters κ = 0 . , θ = 0 . , ρ = 0 . , r = 0 . and σ L . The time to maturity is T = 5 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation. Results of the approximations for the parameters κ = 0 . , θ = 0 . , ρ = 0 . , r = 0 . and σ H . The time to maturity is T = 5 and K = 100 . In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation.. In parenthesis theconfidence interval of the Monte Carlo (MC) estimates. The error is defined as thedifference between the MC price and the related approximation.