Semi-closed form solutions for barrier and American options written on a time-dependent Ornstein Uhlenbeck process
SSemi-closed form solutions for barrier and Ameri-can options written on a time-dependent OrnsteinUhlenbeck process
Peter Carr and Andrey Itkin
Tandon School of Engineering, New York University, 1 Metro Tech Center, 10th floor, Brooklyn NY 11201, USA
March 31, 2020 I n this paper we develop a semi-closed form solutions for the barrier (perhaps, time-dependent) and American options written on the underlying stock which follows atime-dependent OU process with a lognormal drift. This model is equivalent to thefamiliar Hull-White model in FI, or a time dependent OU model in FX. Semi-closed formmeans that given the time-dependent interest rate, continuous dividend and volatilityfunctions, one need to solve numerically a linear (for the barrier option) or nonlinear(for the American option) Fredholm equation of the first kind. After that the optionprices in all cases are presented as one-dimensional integrals of combination of the abovesolutions and Jacobi theta functions. We also demonstrate that computationally ourmethod is more efficient than the backward finite difference method used for solvingthese problems, and can also be as efficient as the forward finite difference solver whileproviding better accuracy and stability. Introduction
The Orstein-Uhlenbeck (OU) process with time-dependent coefficients is very popular among practitionersfor modeling interest rates and credit. That is because it is relatively simple, allows negative interestrates (which recently has become a hot feature) , and could be calibrated to the given term-structureof interest rates and to the prices or implied volatilities of caps, floors or European swaptions since themean-reversion level and volatility are functions of time. The most known among this class are Hull-Whiteand Vasicek models, see (Brigo and Mercurio, 2006) and references therein.The Hull-White model is a one-factor model for the stochastic short interest rate r t of the form dr t = k [ θ ( t ) − r t ] dt + σ ( t ) dW t , (1)where t is the time, k > θ ( t ) is the mean-reversion level, σ ( t ) isthe volatility of the process, W t is the standard Brownian motion under the risk-neutral measure. This1 a r X i v : . [ q -f i n . P R ] M a r emi-closed form solutions for barrier and American options ... model can also be used for pricing Equity or FX derivatives if one assumes that the mean-reversion levelvanishes, while the mean-reversion rate is replaced either by q ( t ) − r ( t ) for Equities, or by r f ( t ) − r d ( t ) forFX, where r ( t ) , q ( t ) are the deterministic interest rate and continuous dividends, and r d ( t ) , r f ( t ) are thedeterministic domestic and foreign interest rates.Without loss of generality in this paper we are concentrated on the Equity world. Since the processin Eq. (1) is Gaussian the model is tractable for pricing European plain vanilla options. However, forexotic options, e.g., highly liquid barrier options, or for American options these prices are not known yetin closed form. Therefore, various numerical methods are used to obtain them, that sometimes couldbe computationally expensive. In this paper we attack this problem by constructing a semi-closed formsolutions for the prices of barrier and American options written on the process Eq. (1). The resultsobtained in the paper are new. Our approach to a certain degree is similar to that in (Mijatovic, 2010), who,however, used a different underlying process (the lognormal model with local spot-dependent volatility, andconstant interest rates and dividends, but time-dependent barriers). Therefore, our model is more generalin a sense that all parameters of the model are time-dependent, while adding time-dependent barrierscan be naturally done within our approach. Also, as compared with (Mijatovic, 2010), we don’t use aprobabilistic argument, rather a theory of partial differential equations (PDE). At the end we demonstratethat computationally our method is more efficient than the backward finite difference method used tosolve these problems, and can also be as efficient as the forward finite difference solver while providingbetter accuracy and stability. We start by specifying the dynamics of the underlying stock price S t to be dS t = [ r ( t ) − q ( t )] S t dt + σ ( t ) dW t , S t =0 = S , (2)where now r ( t ) is the deterministic short interest rate, and S t is the stock price . Here we don’t specifythe explicit form of r ( t, q ( t ) , σ ( t ) but assume that they are known either as a continuous functions of time t ∈ [0 , ∞ ), or as a discrete set of N values for some moments t i , i = 1 , . . . , N .Further in this section we consider a contingent claim written on the underlying process S t in Eq. (2)which is the Up-and-Out barrier Call option. It is known that by the Feynman-Kac formula, (Klebaner,2005) one can obtain a parabolic ( linear ) PDE which solution gives the Up-and-Out barrier Call optionprice C ( S, t ) conditional on S t =0 = S , which reads ∂C∂t + 12 σ ( t ) ∂ C∂S + [ r ( t ) − q ( t )] S ∂C∂S = r ( t ) C. (3)This equation should be solved subject to the terminal condition at the option maturity t = TC ( S, T ) = ( S − K ) + , (4)and the boundary conditions C (0 , t ) = 0 , C ( H, t ) = 0 , (5)where H is the upper barrier.Our goal now is to build a series of transformations to transform Eq. (3) to the heat equation. It is easy to show that this model is equivalent to the Hull-White model. Page 2 of 16 emi-closed form solutions for barrier and American options ...
To transform the PDE Eq. (3) to the heat equation we first make a change of the dependent andindependent variables as follows: S → x/g ( t ) , C ( S, t ) → e f ( x,t ) u ( x, t ) , (6)where new functions f ( x, t ) , g ( t ) has to be determined in such a way, that the equation for u is the heatequation. This can be done by substituting Eq. (6) into Eq. (3) and providing some tedious algebra. Theresult reads f ( x, t ) = k ( t ) − g ( t ) + g ( t )( r ( t ) − q ( t ))2 g ( t ) σ ( t ) x , (7) k ( t ) = 12 log (cid:18) g ( t ) g (0) (cid:19) + 12 Z t [3 r ( s ) − q ( s )] ds. The function g ( t ) solves the following ordinary differential equation (ODE)0 = b ( t ) g ( t ) − g ( t ) + 2 g ( t ) σ ( t ) σ ( t ) + 2 g ( t ) g ( t ) , (8) b ( t ) = 2( r ( t ) − q ( t )) σ ( t ) σ ( t ) − h ( r ( t ) − q ( t )) + r ( t ) − q ( t ) i . The Eq. (8) by substitution g ( t ) → e R t w ( s ) ds (9)can be further transformed to the Riccati equation w ( t ) = b ( t ) + w ( t ) + 2 w ( t ) σ ( t ) σ ( t ) . (10)This equation cannot be solved analytically for arbitrary functions r ( t ) , q ( t ) , σ ( t ), but can be efficientlysolved numerically. Also, in some cases it can be solved in closed form. For instance, if | r ( t ) − q ( t ) | t = ε (cid:28) b ( t ) can be reduced to b ( t ) = 2( r ( t ) − q ( t )) σ ( t ) /σ ( t ).Then assuming in the first approximation on ε | w ( t ) | (cid:29) | r ( t ) − q ( t ) | , (11)we obtain the solution w ( t ) = σ ( t ) D − R t σ ( s ) ds , (12)where D is an integration constant. Thus, Eq. (11) can be re-written asV( t ) (cid:20) V ¯ V ε (cid:21) (cid:29) D ( r ( t ) − q ( t )) , where V( t ) = σ ( t ) is the normal variance, and ¯ V ( t ) = t R t σ ( s ) ds is the average normal variance. Thus,our solution in Eq. (12) is correct if ε (cid:28) V ε/ ¯ V (cid:28)
1, because then D can always be chosen toobey the inequality V( t ) (cid:29) D ( r ( t ) − q ( t )) , ∀ t ∈ [0 , T ].With these explicit definitions Eq. (3) transforms to the form12 σ ( t ) e R t w ( s ) ds ∂ u∂x + ∂u∂t = 0 . (13) Page 3 of 16 emi-closed form solutions for barrier and American options ...
The next step is to make a change of the time variable τ ( t ) → Z Tt σ ( s ) e R s w ( m ) dm ds, (14)so Eq. (13) finally takes the form of the heat equation ∂u∂τ = ∂ u∂x . (15)The Eq. (15) should be solved subject to the terminal condition u ( x,
0) = ( xe − R T w ( s ) ds − K ) + e − f ( x,T ) , (16)and the boundary conditions u (0 , τ ) = 0 , u ( y ( τ ) , τ ) = 0 , y ( τ ) = Hg ( t ( τ )) . (17)These conditions directly follow from Eq. (4), Eq. (5), while y ( τ ) is now a time-dependent upper barrier .The function t ( τ ) is the inverse map of Eq. (14). It can be computed for any t ∈ [0 , T ] by substituting itinto Eq. (14), then finding the corresponding value of τ ( t ), and finally inverting. The PDE in Eq. (15), Eq. (16), Eq. (17) is a parabolic equation which solution should be found atthe domain with moving boundaries. These kind of problems are known in physics for a long time.Similar problems arise in the field of nuclear power engineering and safety of nuclear reactors; in studyingcombustion in solid-propellant rocket engines; in laser action on solids; in the theory of phase transitions(the Stefan problem and the Verigin problem (in hydromechanics)); in the processes of sublimation infreezing and melting; in the kinetic theory of crystal growth; etc., see (Kartashov, 1999) and referencestherein. Analytical solutions of these problems require non-traditional, and sometimes sophisticatedmethods. Those methods were actively elaborated on by the Russian mathematical school in the 20thcentury starting from A.V. Luikov, and then by B.Ya. Lyubov, E.M. Kartashov, and many others.As applied to mathematical finance, one of these methods - the method of heat potentials - wasactively utilized by A. Lipton and his co-authors who solved various problems of mathematical financeby using this approach, see (Lipton, 2002; Lipton and de Prado, 2020) and references therein. Anothermethod that we use in this paper is the method of a generalized integral transform. Below we closelyfollow (Kartashov, 2001) when give an exposition of the method.We start by introducing an integral transform of the form¯ u ( p, τ ) = Z y ( τ )0 u ( x, τ ) sinh( x √ p ) dx, (18)where p = a + i ω is a complex number with R ( p ) ≥ β >
0, and − π < arg (cid:0) √ p (cid:1) < π . Let us multiply bothparts of Eq. (15) by sinh( x √ p ) and then integrate on x from zero to y ( τ ): Z y ( τ )0 ∂u∂τ sinh( x √ p ) dx = Z y ( τ )0 ∂ u∂x sinh( x √ p ) dx. Therefore, we can also naturally solve the same problem with the time-dependent upper barrier H = H ( t ) as this justchanges the definition of y ( τ ). Page 4 of 16 emi-closed form solutions for barrier and American options ... Integrating by parts, we obtain ∂∂τ Z y ( τ )0 u ( x, τ ) sinh( x √ p ) dx − u ( y ( τ ) , τ ) sinh( y ( τ ) √ p ) y ( τ )= ∂u ( x, τ ) ∂x sinh( x √ p ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y ( τ )0 + √ pu ( x, τ ) cosh( x √ p ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y ( τ )0 + p Z y ( τ )0 u ( x, τ ) sinh( x √ p ) dx. With allowance for the boundary conditions in Eq. (17) and the definition in Eq. (18) we obtain thefollowing Cauchy problem ∂ ¯ u∂τ − p ¯ u = Ψ( τ ) sinh( y ( τ ) √ p ) , (19)¯ u ( p,
0) = Z y (0)0 u ( x,
0) sinh( x √ p ) dx, Ψ( τ ) = ∂u ( x, τ ) ∂x (cid:12)(cid:12)(cid:12) x = y ( τ ) . The Eq. (19) can be solved explicitly, assuming that Ψ( τ ) is known. The solution reads¯ ue − pτ = Z τ Ψ( k ) e − pk sinh( y ( k ) √ p ) dk + Z y (0)0 u ( x,
0) sinh( x √ p ) dx. (20)As R ( p ) ≥ β >
0, and ¯ u ( x, τ ) < ∞ , the function ¯ ue − pτ → τ → ∞ . Therefore, letting τ tend to ∞ ,we obtain an equation which makes connection between the moving boundary y ( τ ) and Ψ( τ ): Z ∞ Ψ( τ ) e − pτ sinh( y ( τ ) √ p ) dτ = − Z y (0)0 u ( x,
0) sinh( x √ p ) dx. (21)Using the definitions in Eq. (16) and Eq. (7), the integral in the RHS of Eq. (21) can be represented as F ( p ) ≡ − e − R T w ( s ) ds Z y (0) K ( x − K ) sinh( x √ p ) e k ( T ) − a ( T ) x dx (22)= e − R T w ( s ) ds e k ( T ) − x [ a ( T ) x + √ p ]8 a / ( T ) ( q a ( T ) (cid:16) e √ px − (cid:17) − √ πe [ a ( T ) x + √ p ] a ( T ) · " [ √ p − a ( T ) K ] erf a ( T ) x − √ p p a ( T ) ! + [ √ p + 2 a ( T ) K ] erf a ( T ) x + √ p p a ( T ) ! y (0) K ,K = Ke R T w ( s ) ds , a ( t ) = g ( t ) + g ( t )( r ( t ) − q ( t ))2 g ( t ) σ ( t ) . Thus, Eq. (21) takes the form Z ∞ Ψ( τ ) e − pτ sinh( y ( τ ) √ p ) dτ = F ( p ) , (23)where F(p) is known from Eq. (22).The Eq. (23) is a linear Fredholm integral equation of the first kind, (Polyanin and Manzhirov, 2008).The solution Ψ( τ ) can be found numerically on a grid by solving a system of linear equations. In otherwords, given functions r ( t ) , q ( t ) , σ ( t ) we can compute first w ( t ), then g ( t ), and finally τ ( t ) (or t ( τ )), thusdetermining the moving boundary y ( τ ). Next we can solve Eq. (23) for Ψ( τ ) and substitute it into Eq. (20)to obtain the generalized transform of u ( x, τ ) in the explicit form. Therefore, if this transform can beinverted back, we solved the problem of pricing Up-and-Out barrier Call options. Page 5 of 16 emi-closed form solutions for barrier and American options ...
In this Section the description of inversion is borrowed from (Kartashov, 2001). Since that book hasnever been translated into English, we provide a wider exposition of the method. Also, the book containsvarious typos that are fixed here.As known from a general theory of the heat equation, the solution of the heat equation L u ( x, τ ) =0 , L = ∂/∂τ − ν∂ /∂x at the space domain 0 < x < l = const can be expressed via Fourier series of theform, (Polyanin, 2002) u ( x, τ ) = ∞ X n =1 α n e − νγ n t sin (cid:18) nπxl (cid:19) , where ψ ( x ) = sin( nπx/l ) are the eigenfunctions of the heat operator L , and γ n = nπ/l are its eigenvalues.Therefore, by analogy we look for the inverse transform of ¯ u , or for the solution of Eq. (20) in terms of u ( x, τ ), to be a generalized Fourier series of the form, (Kartashov, 2001) u ( x, τ ) = ∞ X n =1 α n ( τ ) e − (cid:0) nπy ( τ ) (cid:1) τ sin (cid:18) nπxy ( τ ) (cid:19) . (24)where α ( τ ) are some functions to be determined. Note, that this definition automatically respects thevanishing boundary conditions for u ( x, τ ). We assume that this series converges absolutely and uniformly ∀ x ∈ [0 , y ( τ )] for any τ > ∞ X n =1 ( − n +1 nα n ( τ ) e − ( nπ/y ( τ )) τ p + ( nπ/y ( τ )) = y ( τ ) π sinh ( √ py ( τ )) ¯ u ( p, τ ) . (25)The LHS of this equation is regular everywhere except simple poles on the negative semi-axis, see Fig, 1 p n = − (cid:18) nπy ( τ ) (cid:19) , n = 1 , , ... Let us sequentially integrate both sides of Eq. (25) on p along contours γ , γ , . . . . The contour γ n consistsof the vertical line γ >
0, the half-round of radius R n = [ π / (2 y ( τ )](2 n + 2 n + 1) (the contour γ n crossesthe Re ( p ) axis in the middle point between p n and p n +1 with the center in the origin), and two horizontallines Y = ± [ π / (2 y ( τ ))](2 n + 2 n + 1). It means that the circle R n doesn’t hit any pole of the LHS ofEq. (25). Then by the Cauchy’s residual theorem, (Mitrinovic and Keckic, 1984) the integral taken alongthe contour γ n is equal to 2 π i times the sum of residuals of the LHS of Eq. (25) that lie inside γ n As poles are simple, and the function under the integral in the LHS of Eq. (25) has the form F ( p ) /F ( p ),the residual of such a function is, (Mitrinovic and Keckic, 1984)Res[ F ( p ) /F ( p ); p k ] = F ( p ) /F ( p ) (cid:12)(cid:12)(cid:12) p = p k The above analysis is the basis for running a residual machinery to calculate all the coefficients α n ( τ ). Let us denote via I k the following contour integral I k = 12 π i I γ k ¯ u ( p, τ )sinh (cid:0) √ py ( τ ) (cid:1) dp. Below we show that all coefficients α n , n = 1 , . . . , ∞ can be expressed via these integrals. Page 6 of 16 emi-closed form solutions for barrier and American options ... • p γ • p γ • p γ • p γ Re p Im p γ . . . Figure 1:
Contours of integration in a complex plane: p n , n = 1 , ... are simple poles of the LHS of theEq. (25), γ n - the integration contours. Figure 2:
Typical exercise boundaryfor the American Call option.
1. Coefficient α ( τ ) . Integrating Eq. (25) along the contour γ gives α ( τ ) e − ( π/y ( τ )) τ I γ p + ( π/y ( τ )) dp + ∞ X n =2 ( − n +1 nα n ( τ ) e − ( nπ/y ( τ )) τ I γ p + ( nπ/y ( τ )) dp = y ( τ ) π I γ ¯ u ( p, τ ) dp sinh ( √ py ( τ )) . Observe, that I γ p + ( π/y ( τ )) dp = 2 π i , I γ p + ( nπ/y ( τ )) dp = 0 , n ≥ , where the second result is due to the Cauchy integral theorem, (Mitrinovic and Keckic, 1984). Then α ( τ, y ) = y ( τ ) π e ( π/y ( τ )) τ I . (26)
2. Coefficient α ( τ ) . By analogy, integrating the second equation of Eq. (25) along the contour γ weobtain α ( τ ) e − ( π/y ( τ )) τ I γ dpp + ( π/y ( τ )) − α ( τ ) e − (2 π/y ( τ )) τ I γ dpp + (2 π/y ( τ )) + ∞ X n =3 ( − n +1 α n ( τ ) e − ( nπ/y ( τ )) τ I γ dpp + ( nπ/y ( τ )) = y ( τ ) π I γ ¯ u ( p, τ )sinh ( √ py ( τ )) dp, whence using again the residual theorem and Eq. (26) we find α ( τ, y ) = − y ( τ )2 π e (2 π/y ( τ )) τ [ I − I ] . Page 7 of 16 emi-closed form solutions for barrier and American options ...
3. Coefficient α n ( τ ) . Proceeding in a similar manner, we obtain a general formula for the coefficients α n , n ≥ α n ( τ ) = ( − n +1 y ( τ ) nπ e ( nπ/y ( τ )) τ [ I n − (1 − δ n, ) I n − ] , (27)where δ n, is the Kronecker symbol. To calculate the integrals in the RHS of Eq. (27), we rewrite them in the explicit form by using thesolution for ¯ u ( p, τ ) previously found in Eq. (20) I j = 12 π i I γ k e pτ sinh( y ( τ ) √ p ) "Z τ Ψ( s ) e − ps sinh( y ( s ) √ p ) ds + Z y (0)0 u ( x,
0) sinh( x √ p ) dx dp. As sinh( x ) is a periodic complex function with the period πk/i , the RHS of this equation is regulareverywhere except simple poles , where sinh (cid:0) √ py ( τ ) (cid:1) vanishes. It is easy to checks that these poles areexactly p i , i = 1 , . . . , k . Therefore, we again can directly apply the Cauchy residual theorem. Computingresiduals, after some algebra we obtain α n ( τ ) = 2 y ( τ ) " Z y (0)0 u ( x,
0) sin (cid:18) nπxy ( τ ) (cid:19) dx + Z τ e ( nπ/y ( τ )) s Ψ( s ) sin (cid:18) nπy ( s ) y ( τ ) (cid:19) ds. (28)Thus, from Eq. (24) and Eq. (28) we find the final solution u ( x, τ ) = 2 y ( τ ) " ∞ X n =1 e − (cid:0) nπy ( τ ) (cid:1) τ sin (cid:18) nπxy ( τ ) (cid:19) Z y (0)0 u ( x,
0) sin (cid:18) nπxy ( τ ) (cid:19) dx + ∞ X n =1 sin (cid:18) nπxy ( τ ) (cid:19) Z τ e ( nπ/y ( τ )) ( s − τ ) Ψ( s ) sin (cid:18) nπy ( s ) y ( τ ) (cid:19) ds . This can also be re-written as u ( x, τ ) = 2 y ( τ ) " Z y (0)0 dz u ( z, ∞ X n =1 e − (cid:0) nπy ( τ ) (cid:1) τ sin (cid:18) nπxy ( τ ) (cid:19) sin (cid:18) nπzy ( τ ) (cid:19) (29)+ Z τ ds Ψ( s ) ∞ X n =1 e ( nπ/y ( τ )) ( s − τ ) sin (cid:18) nπy ( s ) y ( τ ) (cid:19) sin (cid:18) nπxy ( τ ) (cid:19) . We proceed with the observation that the sums in Eq. (29) could be expressed via Jacobi theta functionsof the third kind, (Mumford et al., 1983). By definition θ ( z, ω ) = 1 + 2 ∞ X n =1 ω n cos(2 nz ) . (30)Therefore, ∞ X n =1 e − (cid:0) nπy ( τ ) (cid:1) τ sin (cid:18) nπxy ( τ ) (cid:19) sin (cid:18) nπzy ( τ ) (cid:19) = 14 [ θ ( φ − ( z ) , ω ) − θ ( φ + ( z ) , ω )] , (31) ∞ X n =1 e ( nπ/y ( τ )) ( s − τ ) sin (cid:18) nπy ( s ) y ( τ ) (cid:19) sin (cid:18) nπxy ( τ ) (cid:19) = 14 [ θ ( φ − ( y ( s )) , ω ) − θ ( φ + ( y ( s )) , ω )] ,ω = e − (cid:16) π √ τy ( τ ) (cid:17) , ω = e (cid:16) π √ s − τy ( τ ) (cid:17) , φ − ( z ) = π ( x − z )2 y ( τ ) , φ + ( z ) = π ( x + z )2 y ( τ ) . Page 8 of 16 emi-closed form solutions for barrier and American options ...
A well-behaved theta function must have parameter | ω | <
1, (Mumford et al., 1983). This conditionholds at any τ > u ( x, τ ) = 12 y ( τ ) " Z y (0)0 dz u ( z, θ ( φ − ( z ) , ω ) − θ ( φ + ( z ) , ω )] (32)+ Z τ ds Ψ( s )[ θ ( φ − ( y ( s )) , ω ) − θ ( φ + ( y ( s )) , ω )] . The RHS of Eq. (29) depends on x via functions φ − , φ + . Since the theta function θ ( z, ω ) is even in z ,the boundary condition at x = 0 is satisfied. At x = y ( τ ) it is also satisfied as follows from Eq. (31) if onereads it from right to left.The result in Eq. (32) to some extent is not a surprise, as it is known that the Jacobi theta function isthe fundamental solution of the one-dimensional heat equation with spatially periodic boundary conditions,(Ohyama, 1995). We recall that an American option is an option that can be exercised at anytime during its life. Americanoptions allow option holders to exercise the option at any time prior to and including its maturity date,thus increasing the value of the option to the holder relative to European options, which can only beexercised at maturity. The majority of exchange-traded options are American. For a more detailedintroduction, see (Detemple, 2006; Hull, 1997).It is known that pricing American (or Bermudan) options requires solution of a linear complimentaryproblem. Various efficient numerical methods have been proposed for doing that. For instance, when theunderlying stock price S t follows the time-dependent Black-Scholes model these (finite-difference) methodsare discussed in (Itkin, 2017) (see also references therein).Another approach, elaborated e.g., in (Andersen et al., 2016) for the Black-Scholes model with constantcoefficients, uses a notion of the exercise boundary S B ( t ). The boundary is defined in such a way, that, e.g.,for the American Put option P A ( S, t ) at S ≤ S B ( t ) it is always optimal to exercise the option, therefore P A ( S, t ) = K − S . For the complimentary domain S > S B ( t ) the earlier exercise is not optimal, and in thisdomain P A ( S, t ) obeys the Black-Scholes equation. This domain is called the continuation (holding) region.The problem of pricing American options lies in the fact that S B ( t ) is not known in advance. Instead, weonly know the price of the American option at the boundary. For instance, for the American Put we have P A ( S B ( t ) , t ) = K − S B ( t ), and for the American Call - C A ( S B ( t ) , t ) = S B ( t ) − K . A typical shape of theexercise boundary for the Call option obtained with the parameters K = 100 , r = 0 . , q = 0 . , σ = 0 . S B ( t ) by numerically solvingan integral (Volterra) equation for S B ( t ). The resulting scheme is straightforward to implement andconverges at a speed several orders of magnitude faster than existing approaches.In terms of this paper, the continuation region is a domain with the moving boundary where theoption price solves the corresponding PDE. In case of our model in Eq. (2) this is the PDE in Eq. (3).Therefore, this problem is, by nature, similar to that for the barrier options considered in Section 1, butthe difference is as follows. • For the barrier option pricing problem the moving boundary (the time-dependent barrier) is known,as this is stated in Eq. (17). But the Option Delta ∂u∂x at the boundary x = y ( τ ) is not, and shouldbe found by solving the linear Fredholm equation Eq. (23). Also, the problem is solved subject tothe vanishing condition at the barrier (the moving boundary) for the option value. Page 9 of 16 emi-closed form solutions for barrier and American options ... • For the American option pricing problem the moving boundary is not known. However, the optionDelta ∂u∂x at the boundary x = y ( τ ) is known (it follows from the conditions ∂C A ∂S | S = S B ( t ) = 1 and ∂P A ∂S | S = S B ( t ) = − x and y ( τ ) according to their definitions in Section 1).Also the boundary condition for the American Call and Put at the exercise boundary (the movingboundary) differs from that for the Up-and-Out barrier option, namely: it is C A ( S B ( t ) , t )) = S B ( t ) − K for the Call, and P A ( S B ( t ) , t )) = K − S B ( t ) for the PutBecause of the similarity of these two problems, it turns out that the American option problem can besolved for the continuation region together with the simultaneous finding of the exercise boundary, byusing the same approach that we proposed for solving the barrier option pricing problem. However, dueto the highlighted differences, some equations slightly change. Since the PDE we need to solve is the same as in Eq. (3), we do same transformations as in Section 1, andcome up to the same heat equation as in Eq. (15). It should be solved subject to the terminal condition u ( x,
0) = (cid:18) xe − R T w ( s ) ds − K (cid:19) + e − f ( x,T ) , and the boundary conditions u (0 , τ ) = 0 , u ( y ( τ ) , τ ) ≡ ψ ( τ ) = y ( τ ) − K, Ψ( τ ) ≡ ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y ( τ ) = e − f ( y ( τ ) ,t ) g ( t ) [1 + a ( t ) y ( τ )( y ( τ ) − K )] , a ( t ) = g ( t ) + g ( t )( r ( t ) − q ( t )) g ( t ) σ ( t ) , where t = t ( τ ). We underline once again that the function y ( τ ) here is not known yet, while Ψ( τ ) isknown. These problems with the free boundaries are also well known in physics.We proceed by using the same transformation in Eq. (18), and by analogy with Eq. (19) obtain thefollowing Cauchy problem ∂ ¯ u∂τ − p ¯ u = Ψ( τ ) sinh( x √ p ) + ψ ( τ ) √ p, ¯ u ( p,
0) = Z y (0)0 u ( x,
0) sinh( x √ p ) dx, This problem can be solved explicitly to yield, (Kartashov, 2001)¯ ue − pτ = Z τ Ψ( τ ) e − pτ sinh( y ( τ ) √ p ) dτ + Z y (0)0 u ( x,
0) sinh( x √ p ) dx + √ p Z τ e − pτ ψ ( τ ) dτ. Accordingly, instead of Eq. (21) we obtain Z ∞ e − pτ " Ψ( τ ) sinh( y ( τ ) √ p ) √ p + y ( τ ) dτ = Kp − √ p Z y (0)0 u ( x,
0) sinh( x √ p ) dx = Kp + F ( p ) √ p . This is a nonlinear Fredholm equation of the first kind, but now with respect to the function y ( τ ). It canalso be solved numerically (iteratively).The next step is to reduce our problem to that with homogeneous boundary conditions. This can bedone by change of the dependent variable u ( x, τ ) = W ( x, τ ) + Θ( x, τ ) , Θ( x, τ ) = (1 − x/y ( τ )] ψ ( τ ) . Page 10 of 16 emi-closed form solutions for barrier and American options ...
The function W ( x, τ ) solves the same heat equation with the same terminal condition, and with thehomogeneous boundary conditions. Therefore, it can be solved by using the method of generalized integraltransform described in Section 1.2. The solution reads u ( x, τ ) = Θ( x, τ ) + ∞ X n =1 α n ( τ ) e − (cid:0) nπy ( τ ) (cid:1) τ sin (cid:18) nπxy ( τ ) (cid:19) ,α n ( τ ) = 2 y ( τ ) " Z y (0)0 [ u ( z, − Θ( z, (cid:18) nπzy ( τ ) (cid:19) dz + Z τ e ( nπ/y ( τ )) s (cid:20) Ψ( s ) + ψ ( s ) y ( s ) (cid:21) sin (cid:18) nπy ( s ) y ( τ ) (cid:19) ds . Again, using the definition of the Jacobi theta function in Eq. (30), this can be finally re-written as u ( x, τ ) = Θ( x, τ ) + 12 y ( τ ) " Z y (0)0 dz [ u ( z, − Θ( z, θ ( φ − ( z ) , ω ) − θ ( φ + ( z ) , ω )]+ Z τ ds (cid:20) Ψ( s ) + ψ ( s ) y ( s ) (cid:21) [ θ ( φ − ( y ( s )) , ω ) − θ ( φ + ( y ( s )) , ω )] . To test performance and accuracy of our method in this Section we provide a numerical example where aparticular time dependence of r ( t ) , q ( t ) , σ ( t ) is chosen as r ( t ) = r e − r k t , q ( t ) = q , σ ( t ) = σ e − σ k t . (33)Here r , q , σ , r k , σ k are constants. With this model Eq. (10) can be solved analytically to yield w ( t ) = q − r e − r k t . (34)Accordingly, from Eq. (9) we find g ( t ) = exp (cid:20) q t + r r k (cid:16) e − r k t − (cid:17)(cid:21) , (35)and from Eq. (7) f ( x, t ) = k ( t ) = 12 (cid:20) log (cid:18) g ( t ) g (0) (cid:19) − q t + 3 r r k (cid:16) − e − r k t (cid:17)(cid:21) . (36)The algorithm described in Section 1 was implemented in python. We did it for two reasons. First,we found neither any standard implementation of the Jacobi theta functions in Matlab, nor any customgood one. Surprisingly, this is also not a part of numpy or scipy packages in python. However, theyare available as a part of the python package mpmath which is a free (BSD licensed) Python library forreal and complex floating-point arithmetic with arbitrary precision, see (Johansson, 2007). It has beendeveloped by Fredrik Johansson since 2007, with help from many contributors.Also, we didn’t find any standard implementation of solver for the Fredholm integral equation of thefirst kind in both python and Matlab. Therefore, we implemented a Tikhonov regularization method asthis is described in (Fuhry, 2001). In particular, with the model used in this Section, the function F ( p )reads F ( p ) = e − k ( T ) p [ −√ p ( K − y ) cosh ( √ py ) + sinh ( K √ p ) − sinh ( √ py )] . (37)Finally to validate the results provided by our method, we implemented a FD solver for pricingUp-and-Out barrier options. This solver is based on the Crank–Nicolson scheme with a few Rannacher Page 11 of 16 emi-closed form solutions for barrier and American options ... first steps, and uses a non-uniform grid, in more detail see, eg., (Itkin, 2017). We implemented two solvers:one for the backward PDE, and the other one - for the forward PDE. But logically, since in this paper wesolved the backward PDE, it does make sense to compare our method with the backward solver. Thisimplementation has been done in Matlab.In our particular test we choose parameters of the model as they are presented in Table 1. r q σ r k σ k H S · H Table 1:
Parameters of the test.
We recall, that here σ ( t ) is the normal volatility. Therefore, we choose its typical value by multiplyingthe log-normal volatility by the barrier level.We run the test for a set of maturities T ∈ [1 / , . , . ,
1] and strikes K ∈ [50 , , , , , , T K
50 55 60 65 70 75 80 C ( K , T ) Figure 3:
Up-and-Out barrier Call option price computed by our method in the test.
In Fig. 4 the relative errors between the Up-and-Out barrier Call option prices obtained by using ourmethod and the FD solver are presented as a function of the option strike K and maturity T . Here toprovide a comparable accuracy we run the FD solver with 101 nodes in space S and the time step ∆ t =0.01.It can be seen that the quality of the FD solution is not sufficient. Therefore, we reran it by using 201nodes in space S and the time step ∆ t = 0.001. The relative errors between our semi-analytic and the FDsolutions in this case are presented in Fig. 5.It can be seen that the agreement becomes better, so the relative error decreases. However, the cost forthis improvement of the FD method is speed. In Table 2 we compare the elapsed time of both methods.The column "no Ψ" has the following meaning. Since the volatility and the interest rate change with timerelatively slow, contribution of the second integral in Eq. (32) to the option price is negligible. Therefore, Page 12 of 16 emi-closed form solutions for barrier and American options ... T K
50 55 60 65 70 75 80 E rr o r ( K , T ) , % -0.82-0.47-0.120.220.570.911.261.601.952.29 −0.50.00.51.01.52.0 Figure 4:
The relative error in the Up-and-Out barrier Call option prices obtained by usingour method and the FD solver with 101 nodes in space S and time step ∆ t = 0.01. T K E rr o r ( K , T ) , % -0.84-0.56-0.29-0.020.250.520.791.061.341.61 −0.50.00.51.01.5 Figure 5:
The relative error in the Up-and-Out barrier Call option prices obtained by usingour method and the FD solver with 201 nodes in space S and time step ∆ t = 0.001. in this particular case we can find the option price by computing only the first integral in Eq. (32).Accordingly, we don’t need to solve the Fredholm equation Eq. (23) that almost halves the elapsed time.Finally, it is known that linear algebra in python (numpy) is almost 3 times slower than that in Matlab. Page 13 of 16 emi-closed form solutions for barrier and American options ...
Test semi-analytic semi-analytic, no Ψ FD-101,∆ t = 0 .
01 FD-201,∆ t = 0 .
01 FD-201,∆ t = 0 . Table 2:
Elapsed time in secs of various tests.
Therefore, given the same accuracy, our method is about 30-40 times faster than the backward FD solver.Of course, the forward FD solver is by an order of magnitude faster than the backward one if weneed to simultaneously price multiple options of various strikes and maturities, but written on the sameunderlying. However, for barrier options this approach requires a very careful implementation, whichoften is not universal and with a lot of tricks involved.
Our attention in Section 1 was drawn to the Up-and-Out barrier Call option C uao . Obviously, using thebarriers parity, (Hull, 1997), the price of the Down-and-Out barrier Call option C dao can be found as C dao = C van − C uao , where C van is the price of the European vanilla Call option. It is known that thelatter is given by the corresponding formula for the process with constant coefficients, where those efficientconstant coefficients ¯ r, ¯ q, ¯ σ are defined as T ¯ r = Z T r ( s ) ds, T ¯ q = Z T q ( s ) ds, T ¯ σ = Z T σ ( s ) ds. Second, as shown in Section 1.2, the barriers also could be some arbitrary functions of time, as thischanges just the definition of function y ( τ ). And our method provides the full coverage of this case withno changes.Third, and perhaps the most important point is about computational efficiency of our method.In addition to what was presented in Section 3, let’s look at this problem from a theoretical pint ofview. Suppose the barrier pricing problem is attacked by solving the forward PDE for a set of strikes K i , i = 1 , . . . , ¯ k and a set of maturities T j , j = 1 , . . . , ¯ m numerically by some FD method on a gridwith N nodes in the space domain S ∈ [0 , H ], and M nodes in the time domain t ∈ [0 , T ¯ m ]. Then thecomplexity of this method is known to be O ( M N + 4 N ). This should be compared with the complexityof our approach.Let’s assume that the Riccati equation in Eq. (10) can be solved either analytically, or, at least,approximately, as this is discussed in Section 1.1. Then the first computational step consists in solvingthe linear Fredholm equation in Eq. (23). This can be done on a rarefied grid with M < M nodes andcomplexity O ( M ). The intermediate values in t can be found (if necessary) by interpolation with thecomplexity O ( M ). As the integral kernel doesn’t depend on strikes K i , this calculation can be donesimultaneously for all strikes still preserving the complexity O ( M ).The final solution of the pricing problem is provided in the form of two integrals in Eq. (32). Therefore,if we need the option price at a single value of S (same as when solving the forward PDE), but forall strikes and maturities, the complexity is O (2¯ kL ( M + N ), where N is the number of points inthe x space, and O ( L ) is the complexity of computing the Jacobi theta function θ ( z, ω ). Normally, M ≤ N, L (cid:28)
N, N (cid:28) N for the typical values of N in the FD method (about 50–100 or even more).Thus, the total complexity of our method is fully determined by the solution of the Fredholm equation.Therefore, our method is slower than the corresponding FD method if M > ( M N ) / . For the Americanoption this situation is worse since instead of solving a linear Fredholm equation we need to solve anonlinear equation. This can be done iteratively, e.g., using k iterations until the method converges to the Page 14 of 16 emi-closed form solutions for barrier and American options ... given tolerance. Then the total complexity becomes O ( kM ). However, our experiments show that usingjust M = 10 points in p space could be sufficient, while further increase of M doesn’t change the results.Also, the accuracy of the method in x can be increased if one uses high-order quadratures for computingthe final integrals. For instance, one can use the Simpson instead of the trapezoid rule that doesn’t affectthe complexity of our method. While increasing the accuracy for the FD method is not easy (i.e., itsignificantly increases the complexity of the method, e.g., see (Itkin, 2017)). References
L. Andersen, M. Lake, and D. Offengenden. High-performance american option pricing.
Journal ofComputational Finance , 20(1):39–87, 2016.D. Brigo and F. Mercurio.
Interest Rate Models – Theory and Practice with Smile, Inflation and Credit .Springer Verlag, 2nd edition, 2006.J. Detemple.
American-Style Derivatives: Valuation and Computation . Financial Mathematics Series.Chapman & Hall/CRC, Boca Raton, London, New York, 2006.M. Fuhry. A new tikhonov regularization mathod. mathesis, Kent State University Honors College, May2001.John 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.F. Johansson. Mpmath library. URL http://mpmath.org/ , 2007.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.F. Klebaner.
Introduction to stochastic calculus with applications . Imperial College Press, London, UK.,2005.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. Mijatovic. Local time and the pricing of time-dependent barrier options.
Finance and Stochastics , 14(1):13–48, 2010.D.S. Mitrinovic and J.D. Keckic.
The Cauchy method of residues: Theory and Applications . Mathematicsand its Applications. Springer, Netherlands, 1984. ISBN 978-90-277-1623-1.D. Mumford, C. Musiliand M. Nori, E. Previato, and M. Stillman.
Tata Lectures on Theta . Progress inMathematics. Birkhäuser Boston, 1983. ISBN 9780817631093.
Page 15 of 16 emi-closed form solutions for barrier and American options ...
Y. Ohyama. Differential relations of theta functions.
Osaka Journal of Mathematics , 32(2):431–450, 1995.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 ofmathematical equations. CRC Press, 2008. ISBN 9780203881057.. Handbooks ofmathematical equations. CRC Press, 2008. ISBN 9780203881057.