Semi-closed form prices of barrier options in the Hull-White model
aa r X i v : . [ q -f i n . C P ] A p r Semi-closed form prices of barrier options in theHull-White model
Andrey Itkin and Dmitry Muravey Tandon School of Engineering, New York University, 1 Metro Tech Center, 10th floor, Brooklyn NY 11201, USA Moscow State University, Moscow, Russia
April 22, 2020 I n this paper we derive semi-closed form prices of barrier (perhaps, time-dependent)options for the Hull-White model, ie., where the underlying follows a time-dependentOU process with a mean-reverting drift. Our approach is similar to that in (Carrand Itkin, 2020) where the method of generalized integral transform is applied to pricingbarrier options in the time-dependent OU model, but extends it to an infinite domain(which is an unsolved problem yet). Alternatively, we use the method of heat potentialsfor solving the same problems. By semi-closed solution we mean that first, we need tosolve numerically a linear Volterra equation of the first kind, and then the option price isrepresented as a one-dimensional integral. Our analysis shows that computationally ourmethod is more efficient than the backward and even forward finite difference methods(if one uses them to solve those problems), while providing better accuracy and stability. Introduction
The Hull-White model since it was invented in (Hull and White, 1990a) became to be very popular amongpractitioners for modeling interest rates and credit. That is because it is relatively simple and allows fornegativity prices (while for a long time this behaviour was treated as a deficiency of the model, nowadaysthis became its advantage). The model could be calibrated to the given term-structure of interest ratesand to the prices or implied volatilities of caps, floors or European swaptions since the mean-reversionlevel and volatility are functions of time. Under the Hull-White model the prices of Zero-coupon bondsand European Vanilla options are known in closed form, (Andersen and Piterbarg, 2010). However, forexotic options, e.g., highly liquid barrier options, these prices are not known yet in closed form. Therefore,various numerical methods are used to obtain such prices.In this paper we present an analytical solution of this problem, and demonstrate that it significantlyaccelerates computation of the prices and, accordingly, calibration of the model. Our contribution istwofold. First, we solve the problem of pricing Barrier options in semi-closed form and provide theresulting expressions not known yet in the literature. Second, we solve this problem by two methods.1 emi-closed form solutions for barrier options ...
The first one is the method of heat potentials, which came to mathematical finance due to A. Lipton whoborrowed it from mathematical physics (see references in the next section). The other method is a methodof generalized integral transform, also known in physics, and introduced to mathematical finance in (Carrand Itkin, 2020). However, this method solves the problem where the underlying is defined at the domain S ∈ [0 , y ( t )] with S being the stock price, and y ( t ) being the time-dependent barrier. Contrary, in thispaper we consider a complimentary domain r ∈ [ y ( t ) , ∞ ) where the solution is not known yet. Therefore,our paper fills this gap, and the constructed solution can be applied to a wide class of problems, both inphysics and finance. In this section we consider a one-factor short interest rate model, first introduced in (Hull and White,1990a), and named after the authors as the Hull-White model. The model assumes dynamics of the shortinterest rate r t to follow the Ornstein-Uhlenbeck (OU) process with time-dependent coefficients dr t = κ ( t )[ θ ( t ) − r t ] dt + σ ( t ) dW t , r t =0 = r. (1)Here t ≥ κ ( t ) > θ ( t ) is the mean-reversion level, σ ( t ) is the volatility of the process, W t is the standard Brownian motion under the risk-neutral measure.This model is also popular for modeling prices of the plain vanilla and exotic options. In particular, inthis paper we consider a Down-and-Out barrier option with the time-dependent lower barrier L ( t ) wherethe underlying is a zero-coupon bond with maturity S and price F ( r, t, S ). We assume that the short interest rate evolves in time as in Eq. (1). It is known that F ( r, t, S ) under arisk-neutral measure solves a linear parabolic partial differential equation (PDE), (Privault, 2012) ∂F∂t + 12 σ ( t ) ∂ F∂r + κ ( t )[ θ ( t ) − r ] ∂F∂r = rF. (2)This equation should be solved subject to the terminal condition F ( r, S, S ) = 1 , (3)and the boundary conditions F (0 , t, S ) = g ( t ) , F ( r, t, S ) (cid:12)(cid:12)(cid:12) r →∞ = 0 , (4)where g ( t ) is some function of the time t . See, for instance, (Zhang and Yang, 2017) and referencestherein.To find g ( t ), recall that according to (Ekström and Tysk, 2011) for single-factor models that predictnonnegative short rates, no second boundary condition is required for Eq. (2) if r t > , ∀ t ≥ r t = 0 is not attainable. Otherwise, the PDE itself at this point serves as the boundarycondition. In particular, as applied to our OU process, it reads ∂F∂t + 12 σ ( t ) ∂ F∂r + κ ( t ) θ ( t ) ∂F∂r ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r =0 = 0 . (5)However, since Eq. (1) is the OU process, it allows zero or even negative interest rates, i.e. in this model r ∈ R . Therefore, the PDE Eq. (2) must be solved subject to the second boundary condition either at Page 2 of 15 emi-closed form solutions for barrier options ... r → −∞ , or at some artificial left boundary r = r min <
0. It is not obvious, however, how to setup thisboundary condition at r min .On the other hand, since the Hull-White model belongs to the class of affine models, (Andersen andPiterbarg, 2010), the solution of Eq. (2) can be represented in the form F ( r, t, S ) = A ( t, S ) e B ( t,S ) r . (6)Substituting this expression into Eq. (2) and separating the terms proportional to r , we obtain twoequations to determine A ( t, S ) , B ( t, S ) ∂B ( t, S ) ∂t = B ( t, S ) κ ( t ) + 1 , (7)2 ∂A ( t, S ) ∂t = − A ( t, S ) B ( t, S ) h θ ( t ) κ ( t ) + B ( t, S ) σ ( t ) i . To obey the terminal condition Eq. (3), the first PDE in Eq. (7) should be solved subject to the terminalcondition B ( S, S ) = 0, and the second one - to A ( S, S ) = 1. The solution reads B ( t, S ) = e R t κ ( x ) dx Z tS e − R x κ ( q ) dq dx, (8) A ( t, S ) = exp (cid:20) − Z tS B ( x, S ) (cid:16) θ ( x ) κ ( x ) + B ( x, S ) σ ( x ) (cid:17) dx (cid:21) . It can be seen that B ( t, S ) < t < S . Therefore, F ( r, t, S ) → r → ∞ . Let us consider a Down-and-Out barrier Call option. It is known that the price of the option C ( r, t )written on this bond under a risk-neutral measure solves the same PDE as in Eq. (2) ∂C∂t + 12 σ ( t ) ∂ C∂r + κ ( t )[ θ ( t ) − r ] ∂C∂r = rC. (9)This PDE should be solved subject to the terminal condition at the option maturity T ≤ S , and someboundary conditions provided to guarantee a uniqueness of the solution. The terminal condition reads C ( r, T ) = ( F ( r, T, S − K )) + , (10)where K is the option strike.By the contract definition, the lower barrier L F is set on the Zero coupon bond (ZCB) price. In otherwords, at the barrier we have the following condition C ( r, t ) = 0 if F ( r, t, S ) = L F . (11)Since the ZCB price F ( r, t, S ) is known in closed form in Eq. (6), the above condition can be reformulatedin the r domain by solving the equation F ( r, t, S ) = A ( t, S ) e B ( t,S ) r = L F , with respect to r . This yields the following equivalent barrier L ( t ) in the r domain L ( t ) = 1 B ( t, S ) log (cid:18) L F A ( t, S ) (cid:19) > , (12) Page 3 of 15 emi-closed form solutions for barrier options ... where it is assumed that L F > A ( t, S ). Thus, the boundary condition to Eq. (9) now reads C ( L ( t ) , t ) = 0 . (13)Hence, in this case the barrier becomes time-dependent. Accordingly, this allows a natural generalizationof Eq. (11) to make the barrier L F being also dependent of time, i.e. L F ( t ).At the second boundary as r → ∞ the ZCB price tends to zero, see Eq. (6), and, therefore, the Calloption price tends to zero. Thus, we set C ( r, t ) (cid:12)(cid:12)(cid:12) r →∞ = 0 . (14)Our goal now is to build a series of transformations to transform Eq. (9) to the heat equation. To transform the PDE Eq. (2) to the heat equation we first rewrite Eq. (9) in the form ∂C∂t = − σ ( t ) ∂ C∂r − [ − κ ( t ) r + κ ( t ) θ ( t )] ∂C∂r + rC. (15)This equation belongs to the type of equations considered in (Polyanin, 2002), Section 3.8.7. It is shownthere that by transformation C ( r, t ) = exp[ α ( t ) r + β ( t )] u ( x, τ ) , τ = φ ( t ) , x = rψ ( t ) + ξ ( t ) , (16)Eq. (15) can be reduced to the heat equation ∂u∂τ = ∂ u∂x . (17)Here, ψ ( t ) = C exp (cid:18)Z t κ ( q ) dq (cid:19) , (18) φ ( t ) = 12 Z Tt σ ( q ) ψ ( q ) dq + C ,α ( t ) = ψ ( t ) Z t ψ ( q ) dq + C ψ ( t ) ,β ( t ) = − Z t α ( q ) h κ ( q ) θ ( q ) + σ ( q ) α ( q ) i dq + C ,ξ ( t ) = − Z t h κ ( q ) θ ( q ) + σ ( q ) α ( q ) i ψ ( q ) dq + C , where C , . . . , C are some constants. In our case we can choose C = 1, C = C = C = C = 0.The Eq. (17) should be solved subject to the initial condition u ( x,
0) = exp (cid:20) − α ( T ) ψ ( T ) ( x − ξ ( T )) − β ( T ) (cid:21) (cid:16) ¯ F ( x, T, S ) − K (cid:17) + , (19)¯ F ( x, T, S ) = A ( T, S ) exp (cid:20) B ( T, S ) ψ ( T ) ( x − ξ ( T )) (cid:21) , as it follows from Eq. (10), Eq. (6) and Eq. (16). Page 4 of 15 emi-closed form solutions for barrier options ...
The boundary conditions should be set at the new domain x ∈ Ω : [ y ( τ ) , ∞ ), where y ( τ ) = L ( t ( τ )) ψ ( t ( τ )) + ξ ( t ( τ )), and t ( τ ) is the inverse map t → τ . The latter can be found explicitly bysolving the second equation in Eq. (18) τ = 12 Z Tt σ ( q ) ψ ( q ) dq, (20)with respect to t.Accordingly, the conditions at the boundaries of Ω can be obtained from Eq. (14), Eq. (13) and read u ( x, τ ) (cid:12)(cid:12)(cid:12) x →∞ = 0 , u ( y ( τ ) , τ ) = 0 . (21) As applied to equities, the problem of pricing barrier options, where the underlying follows a time-dependent OU process, was considered in (Carr and Itkin, 2020). There the authors utilized and extendeda method of generalized integral transform actively elaborated on by the Russian mathematical school tosolve parabolic equations at the domain with moving boundaries, see eg., (Kartashov, 1999) and referencestherein. However, in (Carr and Itkin, 2020; Kartashov, 1999) those problems were solved at the domain x ∈ [0 , y ( τ )] while in this paper we have to deal with the infinite domain Ω. For that kind of domains theabove method is not fully elaborated yet since there is a problem with constructing the inverse transform.Some examples where the inverse transform is not necessary can be found in (Kartashov, 1999, 2001).Therefore, in this Section we will use the method of heat potentials (while the method of generalizedintegral transform is considered in Section 2.3). This method is known in the theory of heat equation fora long time, see, eg., (Tikhonov and Samarskii, 1963; Friedman, 1964.; Kartashov, 2001) and referencestherein. The first use of this method in mathematical finance is due to (Lipton, 2002) for pricing path-dependent options with curvilinear barriers, and more recently in (Lipton and Kaushansky, 2018; Liptonand de Prado, 2020) (also see references therein).Recall, that we have to solve the heat equation in Eq. (17) with the initial condition in Eq. (19) andthe boundary conditions in Eq. (21). Since this problem has an inhomogeneous initial condition, themethod of heat potentials cannot be directly applied. Therefore, let us represent u ( x, τ ) in the form u ( x, τ ) = q ( x, τ ) + 12 √ πτ Z ∞ y (0) u ( x ′ , e − ( x − x ′ )24 τ dx ′ . (22)Here the second term in the RHS is the solution of the heat equation in Eq. (17) at the infinite domainwith the initial condition in Eq. (10). But due to domain of definition x ∈ Ω, we moved the left boundaryfrom −∞ to y (0).The function q ( x, τ ) solves the problem ∂q ( x, τ ) ∂τ = ∂ q ( x, τ ) ∂x , (23) q ( x,
0) = 0 , y (0) < x < ∞ ,q ( x, τ ) (cid:12)(cid:12)(cid:12) x →∞ = 0 , q ( y ( τ ) , τ ) = φ ( τ ) ,φ ( τ ) = − √ πτ Z ∞ y (0) u ( x ′ , e − ( y ( τ ) − x ′ )24 τ dx ′ . This problem is similar to that in Eq. (17), Eq. (10), Eq. (21), but now with a homogeneous initialcondition. Therefore, the solution can be found by the method of heat potentials.
Page 5 of 15 emi-closed form solutions for barrier options ...
Following the general idea of this method, we are looking for the solution of Eq. (23) in the form ofa generalized heat potential q ( x, τ ) = 14 √ π Z τ Ψ( k ) x − y ( k ) p ( τ − k ) e − ( x − y ( k ))24( τ − k ) dk, (24)where Ψ( k ) is the heat potential density. It is easy to check that thus defined function q ( x, τ ) solves thefirst line of Eq. (23), and satisfies the initial condition and the vanishing condition at x → ∞ . Also, fromEq. (24) at the barrier x = y ( τ ) we must have, (Tikhonov and Samarskii, 1963)2 φ ( τ ) = Ψ( τ ) + 12 √ π Z τ Ψ( k ) y ( τ ) − y ( k ) p ( τ − k ) e − ( y ( τ ) − y ( k ))24( τ − k ) dk, (25)since for x = y ( τ ) function q ( x, τ ) is discontinuous, and its limiting value at x = y ( τ ) + 0 is equal to φ ( τ ).The Eq. (25) is a linear Volterra equations of the second kind, (Polyanin and Manzhirov, 2008). Since φ ( τ ) is a continuously-differentiable function, Eq. (25) has a unique continuous solution for Ψ( τ ).The Eq. (25) can be attacked twofold. First, it can be efficiently solved numerically. By a standardapproach, the integral in the RHS is approximated using some quadrature rule with N nodes in k , andthe solution is obtained at M nodes in τ . Since the kernel is proportional to Gaussian, the discrete sumapproximating the integral can be computed with linear complexity O ( N + M ) using the Fast GaussTransform, see eg., (Spivak et al., 2010). The final solution can be obtained by using Pickard’s iterations.Another approach is discussed in Section 4.Second, if in Eq. (20) τ (0) is small, we can approximate a curvilinear boundary y ( τ ) by a linearfunction y ( τ ) = a + bτ, a = y (0) , b = y ( τ (0)) − y (0) y ( τ (0)) . (26)Then the integral kernel becomes a function of the difference τ − k , and so the integrand is a convolutionalfunction. Thus, Eq. (25) can be solved by using the Laplace or G- transforms. For instance, with allowancefor Eq. (26) let us re-write Eq. (25) in the form, (Kartashov, 2001) φ ( τ ) = Ψ ( τ ) + b √ π Z τ Ψ ( k ) √ τ − k dk, (27) φ ( τ ) = 2 φ ( τ ) e b τ/ , Ψ ( k ) = Ψ( k ) e b k/ . This is the Abel integral equation of the second kind, (Polyanin and Manzhirov, 2008), which can besolved in closed form by using the Laplace transform. Since φ (0) = 0, the solution readsΨ ( τ ) = F ( τ ) + b Z τ e τ − k F ( k ) dk, F = φ ( τ ) − b √ π Z τ φ ( k ) √ τ − k dk. In case where the linear approximation is too crude, this expression can be used as a smart initial guessfor the function Ψ( τ ) which is needed by the iterative numerical method described in above.Once Eq. (25) is solved and the function Ψ( τ ) is found, the final solution reads u ( x, τ ) = 14 √ π Z τ Ψ( k ) x − y ( k ) p ( τ − k ) e − ( x − y ( k ))24( τ − k ) dk + 12 √ πτ Z ∞ y (0) u ( x ′ , e − ( x − x ′ )24 τ dx ′ . (28) Page 6 of 15 emi-closed form solutions for barrier options ...
The second integral can be further simplified with allowance for the definition of u ( x,
0) in Eq. (10). Thisyields 12 √ πτ Z ∞ y (0) u ( x ′ , e − ( x − x ′ )24 τ dx ′ (29)= 12 e − β ( T ) " e A x + A A ( T, S ) (cid:18) Erf (cid:18) x − y (0) + 2 τ A √ τ (cid:19) − Erf (cid:18) x − K + 2 τ A √ τ (cid:19)(cid:19) − Ke B x + B (cid:18) Erf (cid:18) x − y (0) + 2 τ B √ τ (cid:19) − Erf (cid:18) x − K + 2 τ B √ τ (cid:19)(cid:19) ,A = A ψ ( T ) [ τ ( B ( T, S ) − α ( T )) − ξ ( T ) ψ ( T )] , A = B ( T, S ) − α ( T ) ψ ( T ) ,B = B ψ ( T ) [ − τ α ( T ) − ξ ( T ) ψ ( T )] , B = − α ( T ) ψ ( T ) ,K = max n ξ ( T ) + ψ ( T ) B ( T, S ) log (cid:18) KA ( T, S ) (cid:19) , y (0) o . Also, by definition in Eq. (23), the function φ ( τ ) can be represented in closed form just by substituting x = y ( τ ) into Eq. (29) and multiplying by -1. In this Section we solve the same problem but using the method of generalized integral transform. Thismethod was invented by the Russian mathematical school in the 20th century starting from A.V. Luikov,and then by B.Ya. Lyubov, E.M. Kartashov, and some others, see a detailed survey in (Kartashov, 1999).However, as mentioned in (Kartashov, 2001), the solution for a semi-infinite domain is not known yet,while some recommendations were given on how one can try to proceed. Therefore, to the best of authors’knowledge, the solution presented in this Section is new and compliments the method of heat potentialspresented in Section 2.2.We attack this problem by introducing the following integral transform¯ u ( p, τ ) = Z ∞ y ( τ ) u ( x, τ ) e −√ px dx, (30)where p = a + iω is a complex number with Re( p ) = β >
0, and − π/ < arg( √ p ) < π/
4. Let us multiplyboth parts of Eq. (17) by e − x √ p and then integrate on x from y ( τ ) to ∞ : ∂∂τ Z ∞ y ( τ ) u ( x, τ ) e −√ px dx + u ( y ( τ ) , τ ) e −√ py ( τ ) y ′ ( τ ) = (31)= ∂u∂x e −√ px (cid:12)(cid:12)(cid:12) x → + ∞ x = y ( τ ) + √ pu ( x, τ ) e −√ px (cid:12)(cid:12)(cid:12) x → + ∞ x = y ( τ ) + p Z ∞ y ( τ ) u ( x, τ ) e −√ px dx. By taking into account the boundary conditions in Eq. (21), Eq. (31) can be represented as thefollowing Cauchy problem d ¯ udτ − p ¯ u = Ψ( τ ) e −√ py ( τ ) , (32)¯ u ( p,
0) = Z ∞ y (0) u ( x, e −√ px dx, Ψ( τ ) = − ∂u ( x, τ ) ∂x (cid:12)(cid:12) x = y ( τ ) , Page 7 of 15 emi-closed form solutions for barrier options ... where u ( x,
0) us given in Eq. (19). The solution of this problem reads¯ u ( p, τ ) = e pτ Z τ Ψ( s ) e − ps e −√ py ( s ) ds + e pτ Z ∞ y (0) u ( z, e −√ pz dz. (33)By analogy with (Carr and Itkin, 2020), the function Ψ( τ ) solves the Fredholm equation of the first type Z ∞ Ψ( τ ) e − pτ −√ pτ dτ = F ( p ) , (34)with F ( p ) = − e α ( T ) ξ ( T ) ψ ( T ) − β ( T ) " A ( T, S ) e B ( T,S ) ξ ( T ) ψ ( T ) e − f ( p,T ) y − e − f ( p,T ) K f ( p, T ) − K e − f ( p,T ) y − e − f ( p,T ) K f ( p, T ) (35) f ( p, T ) = √ p + α ( T ) ψ ( T ) − β ( T ) ψ ( T ) , f ( p, T ) = √ p + α ( T ) ψ ( T ) , where K is defined in Eq. (29).To construct the inverse transform, recall that the solution of the heat equation L u ( x, τ ) = 0, L = ∂/∂τ − ∂ /∂x in the half-plane domain x ∈ (0 , ∞ ) can be expressed via the Fourier sine integral, (Cannonand Browder, 1984) Z ∞ α ( ξ ) e − ξ τ sin( ξx ) dξ. Therefore, by analogy we look for the inverse transform of ¯ u (or for the solution of Eq. (32) in terms of u ) to be an oscillatory integral of the form u ( x, τ ) = Z ∞ α ( ξ, τ ) sin[ ξ ( x − y ( τ ))] dξ, (36)where α ( ξ, τ ) is some function to be determined. Note, that this definition automatically respects the van-ishing boundary conditions for u ( x, τ ). We assume that this integral converges absolutely and uniformly ∀ x ∈ [ y ( τ ) , ∞ ) for any τ > u ( p, τ ) = Z ∞ y ( τ ) e −√ px Z ∞ α ( ξ, τ ) sin( ξ ( x − y ( τ ))) dξdx = e −√ py ( τ ) Z ∞ α ( ξ, τ ) ξdξξ + p . (37)Replacing ¯ u ( p, τ ) with the solution found in Eq. (33), we obtain Z ∞ α ( ξ, τ ) ξdξξ + p = Z τ Ψ( s ) e p ( τ − s ) e √ p ( y ( τ ) − y ( s )) ds + e pτ Z ∞ y (0) u ( z, e √ p ( y ( τ ) − z ) dz. (38)Now, similar to inverse operator methods, like the inverse Laplace transform, we need an analytic contin-uation of the transform parameter p into the complex plane. Let us integrate both sides of Eq. (38) on p along the so-called keyhole contour presented in Fig. 1. In more detail, this contour can be describedas follows. It starts with a big symmetric arc Γ around the origin with the radius R ; extending to twohorizontal line segments l , l (a cut around the line Im p = 0 , Re p > γ ε around the origin with the radius ε ≪
1; then extending to two vertical line segments up topoints Im p = ± ξ ; then again two horizontal parallel line segments l , l at Im p = ± ξ , which end pointsare connected to the arc Γ with a cut at Im p = − ξ (it consists of two vertical line segments and twosemi-circles γ r with the radius ε ), such that the whole contour is continuous.Using a standard technique, we take a limit ε → , R → ∞ , so in this limit the contour takes theform as depicted in Fig, 2. It has a horizontal cut along the positive real line with point p = 0 excluded Page 8 of 15 emi-closed form solutions for barrier options ... Re p Im p γ ε γ r ••− ξ Γ l : √ p = i ξl : √ p = − i ξ l l Figure 1:
Contour of integration of Eq. (38) in a complex plane of p . from the area inside the contour; another vertical cut at Re p = − ξ with the point p = − ξ lying insidethe contour; and a branch cut l , l of the multivalued function √ p at p = − ξ . Also, in this limit l → , l →
0, but in Fig. 2 we left them as it is for a better readability.Now we are ready to compute the integrals in Eq. (38). That one in the LHS is regular everywhereinside this contour except a single pole p = − ξ . By the residue theorem, we obtain I γ (cid:18)Z ∞ α ( ξ, τ ) ξdξξ + p (cid:19) dp = − πi Z ∞ ξα ( ξ, τ ) dξ. (39)The integral in the RHS doesn’t have any singularity inside the contour γ , however, it has severalcuts. As can be easily checked, the integrals along the segments l and l cancel out, as well as thosealong l and l , and those along l and l . The integral along the contour Γ tends to zero if R → ∞ dueto Jordan’s lemma. Hence, the only remaining integrals are those along the horizontal semi-infinite lines l and l . They read Z l Z τ Ψ( s ) e p ( τ − s ) e √ p ( y ( τ ) − y ( s )) ds + e pτ Z ∞ y (0) u ( z, e √ p ( y ( τ ) − z ) dz ! dp (40)= − Z ∞ ξ Z τ Ψ( s ) e − ξ ( τ − s ) e iξ ( y ( τ ) − y ( s )) ds + e − ξ τ Z ∞ y (0) u ( z, e iξ ( y ( τ ) − z ) dz ! dξ, Z l Z τ Ψ( s ) e p ( τ − s ) e √ p ( y ( τ ) − y ( s )) ds + e pτ Z ∞ y (0) u ( z, e √ p ( y ( τ ) − z ) dz ! dp = 2 Z ∞ ξ Z τ Ψ( s ) e − ξ ( τ − s ) e − iξ ( y ( τ ) − y ( s )) ds + e − ξ τ Z ∞ y (0) u ( z, e − iξ ( y ( τ ) − z ) dz ! dξ. (41) Page 9 of 15 emi-closed form solutions for barrier options ... Re p Im p •− ξ l l l l l l l l γ Figure 2:
Contour of integration γ of Eq. (38) in a complex plane of p at ε → , R → ∞ . The RHS of Eq. (38) is equal to a sum of these integrals4 i Z ∞ ξ "Z τ Ψ( s ) e − ξ ( τ − s ) sinh ( iξ ( y ( τ ) − y ( s ))) ds + e − ξ τ Z ∞ y (0) u ( z,
0) sin ( ξ ( y ( τ ) − z ) dz dξ. Equating the LHS and RHS provides an explicit representation of α ( ξ, τ ) α ( ξ, τ ) = − π "Z τ Ψ( s ) e − ξ ( τ − s ) sin ( ξ ( y ( τ ) − y ( s ))) ds + e − ξ τ Z ∞ y (0) u ( z,
0) sin ( ξ ( y ( τ ) − z ) dz (42)Finally, we substitute Eq. (42) into Eq. (38) and take into account the identity, (Gradshtein and Ryzhik,2007) Z ∞ e − βx sin ( ax ) sin ( bx ) dx = 14 r πβ (cid:18) e − ( a − b )24 β − e − ( a + b )24 β (cid:19) , β > , which yields u ( x, τ ) = 12 √ π Z τ Ψ( s ) √ τ − s e − ( x − y ( s ))24( τ − s ) − e − ( x − y ( τ )+ y ( s ))24( τ − s ) ! ds (43)+ 12 √ πτ Z ∞ y (0) u ( z, (cid:18) e − ( x − z )24 τ − e − ( x − y ( τ )+ z )24 τ (cid:19) dz. Thus, we obtained another representation of the solution which reads a bit different from that in Eq. (28),despite the general ansatz looks similar (the solution is a sum of two integrals, one in time τ and theother one in space x of the initial condition with a Gaussian weight). The difference can be attributed to Page 10 of 15 emi-closed form solutions for barrier options ... the different definitions of function Ψ( k ), as in Eq. (24) it is the heat potential density, while in Eq. (32)this is the gradient of the solution at x = y ( τ ). The first one is determined by the solution of the Volterraequation of the second kind Eq. (25), and the second one - by the solution of the Fredholm equation ofthe first kind in Eq. (34). However, the latter can either be transformed to the Volterra equation of thesecond kind. For doing that, one needs to differentiate Eq. (43) on x , and then let x = y ( τ ). This yieldsΨ( τ ) = 12 √ π Z τ Ψ( s ) y ( τ ) − y ( s )( τ − s ) / e − ( y ( τ ) − y ( s ))24( τ − s ) ds + 12 √ πτ Z ∞ y (0) u ( z, y ( τ ) − z ) e − ( z − y ( τ ))24 τ dz. (44) Similar to Section 2.2, the method of heat potentials can be used to price double barrier options with thelower barrier L F and the upper barrier H F . In this case, we have the following problem to solve ∂u ( x, τ ) ∂τ = ∂ u ( x, τ ) ∂x , (45) u ( x, τ = 0) = u ( x, , y (0) < x < z (0) ,u ( y ( τ ) , τ ) = u ( z ( τ ) , τ ) = 0 , where z ( τ ) = H ( t ( τ )) ψ ( t ) + ξ ( t ) is the moving upper boundary, H ( t ( τ ) is defined in Eq. (12) by replacing L F with H F , and u ( x,
0) is defined in Eq. (19).Similar to Eq. (22) we represent the solution in the form u ( x, τ ) = q ( x, τ ) + 12 √ πτ Z z (0) y (0) u ( x ′ , e − ( x − x ′ )24 τ dx ′ , (46)so function q ( x, τ ) solves a problem with homogeneous initial condition ∂q ( x, τ ) ∂τ = ∂ q ( x, τ ) ∂x , (47) q ( x,
0) = 0 , y (0) < x < z (0) ,q ( y ( τ ) , τ ) = − φ ( τ ) , q ( z ( τ ) , τ ) = − ψ ( τ ) ,φ ( τ ) = − √ πτ Z z (0) y (0) u ( x ′ , e − ( y ( τ ) − x ′ )24 τ dx ′ , ψ ( τ ) = − √ πτ Z z (0) y (0) u ( x ′ , e − ( z ( τ ) − x ′ )24 τ dx ′ . Again, we are looking for the solution of Eq. (47) in the form of a generalized heat potential q ( x, τ ) = 14 √ π Z τ p ( τ − k ) ( x − y ( k ))Ψ( k ) e − ( x − y ( k ))24( τ − k ) + ( x − z ( k ))Φ( k ) e − ( x − z ( k ))24( τ − k ) ! dk, (48)where Ψ( k ) , Φ( k ) are the heat potential densities. They solve a system of two Volterra equations of thesecond kind2 φ ( τ ) = Ψ( τ ) + 12 √ π Z τ Ψ( k ) y ( τ ) − y ( k )( τ − k ) / e − ( y ( τ ) − y ( k ))24( τ − k ) + Φ( k ) y ( τ ) − z ( k )( τ − k ) / e − ( y ( τ ) − z ( k ))24( τ − k ) ! dk, (49)2 ψ ( τ ) = Φ( τ ) + 12 √ π Z τ Ψ( k ) z ( τ ) − y ( k )( τ − k ) / e − ( z ( τ ) − y ( k ))24( τ − k ) + Φ( k ) z ( τ ) − z ( k )( τ − k ) / e − ( z ( τ ) − z ( k ))24( τ − k ) ! dk, where functions ψ ( τ ) , φ ( τ ) can be expressed in closed form, similar to Eq. (29). This system, again canbe solved by the Variational Iteration Method (VIM), see (Wazwaz, 2011) with a linear complexity byusing the FGT. Once this is done, the solution of our double barrier problem is found. Page 11 of 15 emi-closed form solutions for barrier options ...
To test performance and accuracy of our method we run a test where the explicit form of parameters κ ( t ) , θ ( t ) , σ ( t ) is chosen as κ ( t ) = κ , θ ( t ) = θ e − θ k t , σ ( t ) = σ e − σ k t , (50)with κ , θ , σ , θ k , σ k being constants. With these definitions all functions in Eq. (18) can be found inclosed form.We approach pricing of Down-and-Out barrier Call option written on a ZCB in two ways. First, wesolve the PDE in Eq. (9) by using a finite-difference scheme of the second order in space and time. Weuse the Crank-Nicolson scheme with few first Rannacher steps on a non-uniform grid compressed close tothe barrier level at t = 0, in more detail, see (Itkin, 2017). Second, we use the method of heat potentials(HP) and solve the same problem as this is described in Section 2.2. To solve the Volterra equation inEq. (25) we approximate the kernel on a rectangular grid M × M , and the integral using the trapezoidalrule, which gives rise to the following system of linear equations k φ k = ( I + P ) k Ψ k . (51)Here k Ψ k is the vector of discrete values of Ψ( τ ) , τ ∈ [0 , τ (0)] on a grid with M nodes, k φ k is a similarvector of φ ( τ ), I is the unit M × M matrix, and P is the M × M matrix of the kernel values on thesame grid. Due to the specific structure of the Gaussian kernel, matrix P is lower triangular. Therefore,solution of Eq. (51) can be done with complexity O ( M ).An important point here is that the kernel (and so the matrix P ) doesn’t depend of strikes K , butonly the function φ ( τ ). Therefore, Eq. (51) can be solved simultaneously for all strikes by inverting thematrix I + P with the complexity O ( M ), and then multiplying it by vectors k φ k k , k = 1 , . . . , ¯ k , ¯ k is thetotal number of strikes. Hence, the total complexity is O (¯ kM ).We emphasize, that this algorithm of solving the Volterra equation Eq. (25) doesn’t require iterations,as that which makes use of the FGT. To compare both algorithms, note that complexity of the matrixalgorithm is O ( M ) versus O (2 nM ) - the complexity of the iterative algorithm which requires n iterationsto converge. Therefore, if M is relatively small, eg., M = 20, and the number of iterations is near 10,both algorithms are of the same complexity.For this test parameters of the model are presented in Table 1. r κ θ σ θ k σ k H F S Table 1:
Parameters of the test.
We run the test for a set of maturities T ∈ [1 / , . , . ,
1] and strikes K ∈ [0 . , . , . , . , . , . K and maturity T . Here to provide a comparableaccuracy we run the FD solver with 200 nodes in space r and 201 steps in time t . Otherwise the qualityof the FD solution is poor.A relatively large error about 8% at simultaneously high maturities and strikes is due to the verylow computed price of the option, which is 1.99 cents for the FD method, and 2.17 cents for the HPmethod, respectively. Otherwise, the error is about 1-3%. Also, since in this test we have cosen a fixedbarrier in the ZCB price space, the corresponding barrier in the r space is moving down. Therefore, the Page 12 of 15 emi-closed form solutions for barrier options ...
Figure 3:
Down-and-Out barrier Call option price computed by using the HP method.
Figure 4:
The relative difference in % of the Down-and-Out barrier Call option pricescomputed by using the HP method and the FD solver with 201 nodes in space r and time t . computed option price could vanish for strikes close to the barrier at short maturities, but contrary havesome positive value for the same strikes at longer maturities.Obviously, since the FD solver needs a high number of nodes in both space and time to achievea reasonable accuracy, the cost for this is speed. Suppose that this solver uses N r nodes in space r , Page 13 of 15 emi-closed form solutions for barrier options ... and M t nodes in time t . Then the total complexity of solving the forward PDE to simulatneously getoption prices for all given strikes K i , i = 1 , . . . , ¯ k and maturities T j , j = 1 , . . . , ¯ m is O ( N r × M t ).This should be compared with the complexity of the HP method which is O (¯ k ¯ mM ). Since, as we saw, N r = M t = 200 , ¯ k = 6 , ¯ m = 4 , M = 20, the HP method should be four times faster than the FDmethod. In reality, our test with 24 points in the K × T space shows that the elapsed time of the HPmethod is 20 mls, while for the FD method this is 146 mls which is 7 times slower. Decreasing the sizeof the FD grid to 100 ×
100 nodes also decreases the elapsed time to 21 mls, but by the cost of increasinga relative error up to ±
15% for a wide range of maturities and strikes. Thus, overall, the method of HPdemonstrates, at least same performance as the forward FD solver.
In Sections 2.2,2.3 we constructed semi-closed form solutions for the prices of Down-and-Out barrier Calloption C dao where the underlying is a Zero-Coupon Bond, and the interest rate dynamics follows the Hull-White model. Obviously, using the parity for barrier options, (Hull, 1997), the price of the Up-and-Outbarrier Call option C uao can be found as C dao = C van − C uao , where C van is the price of the Europeanvanilla Call option in the Hull-White model. Since this model allows closed-form solutions for Europeanoptions on Zero-coupon bonds, (Andersen and Piterbarg, 2010), our solution also provides a closed formsolution for C uao . The double barrier case was also considered in Section 3.As shown in Section 2.2, this solution also covers the case when the barrier is some arbitrary functionof time, as this just slightly modifies the definition of function y ( τ ).From the computational point of view the proposed solution is very efficient as this is shown inSection 4. Using theoretical analysis justified by a test example we conclude that our method is, at least,of the same complexity, or even faster than the forward FD method. On the other hand, our approachprovides high accuracy in computing the options prices, as this is regulated by the order of a quadraturerule used to discretize the kernel. Therefore, the accuracy of the method in x space can be easily increasedby using high-order quadratures. For instance, using the Simpson instead of the trapezoid rule doesn’taffect the complexity of our method, while increasing the accuracy for the FD method is not easy (i.e.,it significantly increases the complexity of the method, e.g., see (Itkin, 2017)). References
L.B.G. Andersen and V.V. Piterbarg.
Interest Rate Modeling . Number v. 2 in Interest Rate Modeling.Atlantic Financial Press, 2010. ISBN 9780984422111.J.R. Cannon and F.E. Browder.
The One-Dimensional Heat Equation . Encyclopedia of Mathematicsand its Applications. Cambridge University Press, 1984. ISBN 9780521302432.P. Carr and A. Itkin. Semi-closed form solutions for barrier and American options written on a time-dependent Ornstein Uhlenbeck process, March 2020, https://arxiv.org/abs/2003.08853 .E. Ekström and J. Tysk. oundary conditions for the single-factor term structure equation.
The Annalsof Applied Probability , 21(1):332–350, 2011.A. Friedman.
Partial Differential Equations of Parabolic Type . Prentice-Hall, New Jersey„ 1964.I.S. Gradshtein and I.M. Ryzhik.
Table of Integrals, Series, and Products . Elsevier, 2007.
Page 14 of 15 emi-closed form solutions for barrier options ...
J. Hull and A. White. Pricing interest-rate-derivative securities.
Review Of Financial Studies , 3:573–592,1990a.J.C. Hull.
Options, Futures, and Other Derivatives . Prentice Hall, 3rd edition, 1997.A. Itkin.
Pricing Derivatives Under Lévy Models. Modern Finite-Difference and Pseudo-DifferentialOperators Approach. , volume 12 of
Pseudo-Differential Operators . Birkhauser, 2017.E. M. Kartashov. Analytical methods for solution of non-stationary heat conductance boundary problemsin domains with moving boundaries.
Izvestiya RAS, Energetika , (5):133–185, 1999.E.M. Kartashov.
Analytical Methods in the Theory of Heat Conduction in Solids . Vysshaya Shkola,Moscow, 2001.A. Lipton. The vol smile problem.
Risk , pages 61–65, February 2002.A. Lipton and M.L. de Prado. A closed-form solution for optimal mean-reverting trading strategies, 2020.available at https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3534445 .A. Lipton and V. Kaushansky. On the first hitting time density of an ornstein-uhlenbeck process, October2018. URL https://arxiv.org/pdf/1810.02390.pdf .A.D. Polyanin.
Handbook of linear partial differential equations for engineers and scientists . Chapman& Hall/CRC, 2002.P. Polyanin and A.V. Manzhirov.
Handbook of Integral Equations: Second Edition . Handbooks of math-ematical equations. CRC Press, 2008. ISBN 9780203881057.N. Privault.
An Elementary Introduction to Stochastic Interest Rate Modeling . Advanced series onstatistical science & applied probability. World Scientific, 2012. ISBN 9789814390859.M. Spivak, S.K. Veerapaneni, and L. Greengard. The fast generalized gauss transform.
SIAM Journalon Scientific Computing , 32(5):3092–3107, 2010.A.N. Tikhonov and A.A. Samarskii.
Equations of mathematical physics . Pergamon Press, Oxford, 1963.A. M. Wazwaz.
Linear and Nonlinear Integral Equations . Higher Education Press, Beijing and Springer-Verlag GmbH Berlin Heidelberg, 2011.K. Zhang and X.Q. Yang. Pricing European options on zero-coupon bonds with a fitted finite volumemethod.
International Journal of Numerical Analysis and Modeling , 14(3):405–418, 2017., 14(3):405–418, 2017.