Semi-closed form prices of barrier options in the time-dependent CEV and CIR models
aa r X i v : . [ q -f i n . C P ] M a y Semi-closed form prices of barrier options in thetime-dependent CEV and CIR models
Peter Carr, 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
May 13, 2020 W e continue a series of papers where prices of the barrier options written on theunderlying, which dynamics follows some one factor stochastic model with time-dependent coefficients and the barrier, are obtained in semi-closed form, see(Carr and Itkin, 2020, Itkin and Muravey, 2020). This paper extends this methodology tothe CIR model for zero-coupon bonds, and to the CEV model for stocks which are usedas the corresponding underlying for the barrier options. We describe two approaches.One is generalization of the method of heat potentials for the heat equation to the Besselprocess, so we call it the method of Bessel potentials. We also propose a general schemehow to construct the potential method for any linear differential operator with time-independent coefficients. The second one is the method of generalized integral transform,which is also extended to the Bessel process. In all cases, a semi-closed solution meansthat first, we need to solve numerically a linear Volterra equation of the second kind,and then the option price is represented as a one-dimensional integral. We demonstratethat computationally our method is more efficient than both the backward and forwardfinite difference methods while providing better accuracy and stability. Also, it is shownthat both method don’t duplicate but rather compliment each other, as one providesvery accurate results at small maturities, and the other one - at high maturities. Introduction
This paper continues a series of papers where prices of the barrier options written on the underlying, whichdynamics follows some one factor stochastic model with time-dependent coefficients and the barrier, areconstructed in semi closed form, see (Carr and Itkin, 2020; Itkin and Muravey, 2020). Here we extendour approach for two additional models: the Cox–Ingersoll–Ross (CIR) model, (Cox et al., 1985), and thetime-dependent constant elasticity of variance (CEV) model, (Cox, 1975). Both models are very popularamong practitioners and used for pricing various derivatives such asset classes as Equities, Fixed Income,Commodities, FX, etc. 1 emi-closed form solutions for barrier options ...
For pricing the time-dependent barrier options we develop in parallel two analytic methods, bothbased on the notion of generalized integral transform, (Carr and Itkin, 2020; Itkin and Muravey, 2020).The first method is the method of heat potentials which, as applied to the models considered in thispaper, is discussed in detail in Section 3. We extend this method, since the partial differential equation(PDE) we need to solve cannot be transformed to the heat equation, but rather to the Bessel PDE. Thisapproach is new and has not been developed yet in the literature. However, once this is done, the samemethod could be used for solving some other problems implicitly related to pricing of barrier options, e.g.,pricing American options, (Carr and Itkin, 2020), analyzing the stability of a single bank and a group ofbanks in the structural default framework, (Kaushansky et al., 2018), calculating the hitting time density,(Alil et al., 2005; Lipton and Kaushansky, 2020a), finding an optimal strategy for pairs trading, (Liptonand de Prado, 2020), etc. Also, the method could be used for solving various problems in physics whereit was originally developed for the heat equation, (Kartashov, 2001; Friedman, 1964.) and referencestherein.The other method is the method of generalized integral transform was actively developed by theRussian mathematical school to solve parabolic equations at the domain with moving boundaries, seee.g., (Kartashov, 1999) and references therein. However, again for the Bessel PDE this approach hasnot been developed yet in the literature in full (i.e., up to the final formula), despite some comments onpossible ways of achieving this could be found in (Kartashov, 1999). It is also worth mentioning that so farthe only known problem solved by using this method is the heat equation with the boundary y ( t ) movingin time t , so the solution is defined at the domain [0 , y ( t )]. In (Itkin and Muravey, 2020) this approachwas extended to the domain [ y ( t ) , ∞ ). Also, the problem at the domain [ y ( t ) , z ( t )], which emerges, eg,for the time-dependent double barrier options where the underlying follows a time-dependent Hull-Whitemodel, was also constructed in (Itkin and Muravey, 2020).Going back to the CIR and CEV models with time-dependent coefficients, the prices of the barrieroptions in these models are not known in closed form yet. Instead, various numerical methods are usedto compute them. This obviously produces a computations burden which could be excessive when thesenumerical procedures are used as a part of calibration process. Therefore, our closed form solutions couldbe of importance for practitioners. By a semi-closed solution we mean that first, one needs to solvenumerically a linear Volterra equation of the second kind, and then the option price is represented asa one-dimensional integral of thus found solution. We demonstrate that computationally our methodis more efficient than both the backward and forward finite difference methods while providing betteraccuracy and stability.Overall, our contribution to the existing literature is twofold. First, we solve the problem of pricingbarrier options in the CIR and CEV models in semi-closed form, and provide the resulting expressionsnot known yet in the literature. Second, we solve these problems by two methods which the extension ofthe existing methods and are developed by the authors.The rest of the paper is organized as follows. Section 1 shortly describes the CEV model and showshow to transform the pricing equation to the Bessel PDE. In Section 2 we do same for the CIR model.Note, that we use the CEV model to price barrier options written on Equities, while the CIR modelis used to price barrier options written on Zero-coupon bonds. Nevertheless, the transformed PDE issame for both models, while solution could be defined at same or different domains. In Section 3 wedevelop a method of Bessel potentials which is an extension of the method of heat potentials, and obtainthe solution of our problems using this approach. In Section 4 same program is fulfilled for the methodof generalized integral transform. In Section 5 the results of some numerical experiments are presentedwhich compare the performance and accuracy of our method with the finite-difference method used tosolve the forward Kolmogorov equation (this is currently the standard way to price barrier options in thetime-dependent CIR and CEV models). The last Section concludes.Also, we discovered that the potential method could be constructed for any PDE where the space oper- Page 2 of 32 emi-closed form solutions for barrier options ... ator is a linear differential operator with constant coefficients. We propose and discuss this generalizationof the heat potential method in Appendix.
The time-dependent constant elasticity of variance (CEV) model is a one-dimensional diffusion processthat solves a stochastic differential equation dS t = µ ( t ) S t dt + σ ( t ) S β +1 t dW t , S t =0 = S . (1)Here t ≥ S t is the stochastic stock price, µ ( t ) is the drift, σ ( t ) is the volatility and β is theelasticity parameter such that β < , β = { , − } , W t is the standard Brownian motion. It is known,that under risk-neutral measure µ ( t ) = r ( t ) − q ( t ) where r ( t ) is the deterministic short interest rate,and q ( t ) is the continuous dividend. We assume that all parameters of the model are known either as acontinuous functions of time t ∈ [0 , ∞ ), or as a discrete set of N values for some moments t i , i = 1 , . . . , N .The CEV model with constant coefficients has been introduced in (Cox, 1975) as an alternative tothe geometric Brownian motion for modeling asset prices. Despite some sophistication as compared withthe Black-Scholes model, the model is still analytically tractable, and prices of the European optionscan be obtained in closed-form. That is because the CEV process with constant coefficients is related tothe Bessel process, (Linetsky and Mendoza, 2010). Also, as mentioned in that reference, the elasticityparameter β controls the steepness of the skew (the larger the | β | - the steeper the skew), while thevolatility (scale) parameter σ fixes the at-the-money volatility level. This ability to capture the skew hasmade the CEV model popular in equity options markets.For the standard CEV process with constant parameters it is known that change of variable z t =1 / ( σ | β | ) S − βt reduces the CEV process without drift ( µ = 0) to the standard Bessel process of order 1 / β (see (Revuz and Yor, 1999; Davydov and Linetsky, 2001)). Then the continuous part of the risk-neutraldensity of S t , conditional on S = S , is obtained from the well-known expression for transition density ofthe Bessel process. If µ = 0, using the result of (Goldenberg, 1991) this CEV process could be obtainedfrom the process without drift via a scale and time change S µt = e µt S τ ( t ) , τ ( t ) = 12 µβ (cid:16) e µβt − (cid:17) . Also, as shown in (Davydov and Linetsky, 2001), at β > S t = 0 is a natural boundary, and infinity is an entrance boundary.Below we show that a similar connection with the Bessel process can be established for the time-dependent version of the model in Eq. (1). Let us consider an Upper-and-Out barrier Call option C ( t, S )written on the underlying process S t . By the Feynman–Kac formula this price solves the following partialdifferential equation (PDE) ∂C∂t + 12 σ ( t ) S β +2 ∂ C∂S + [ r ( t ) − q ( t )] S ∂C∂S = r ( t ) C. (2)This equation should be solved subject to the terminal condition at the option maturity t = TC ( T, S ) = ( S − K ) + , (3)where K is the option strike, and the boundary conditions C ( t,
0) = 0 , C ( t, H ( t )) = 0 , (4)where H ( t ) is the upper barrier, perhaps time-dependent. Then the following Proposition holds. In case β = 0 this model is the Black-Scholes model, while for β = − emi-closed form solutions for barrier options ... Proposition 1.1.
The PDE in Eq. (2) can be transformed to ∂u∂τ = 12 ∂ u∂z + bz ∂u∂z , (5) where b is some constant, u = u ( τ, z ) is the new dependent variable, and ( τ, z ) are the new independentvariables. The Eq. (5) is the PDE associated with the one-dimensional Bessel process, (Revuz and Yor,1999) dX t = dW t + bX t dt. (6) Proof.
This transformation can be done in two steps. First, we make a change of variables S = ( − xβ ) − /β , C ( t, S ) → u ( t, x ) e R t r ( k ) dk , φ = Z Tt σ ( k ) dk. (7)This reduces the PDE in Eq. (2) to the form ∂u∂φ = 12 ∂ u∂x + (cid:18) xf ( t ) + bx (cid:19) ∂u∂x , (8) f ( t ) = β r ( t ) − q ( t ) σ ( t ) , b = β + 12 β , t = t ( φ ) . The function t ( φ ) is the inverse map of the last term in Eq. (7). It can be computed for any t ∈ [0 , T ] bysubstituting it into the definition of φ , then finding the corresponding value of φ ( t ), and finally inverting.The Eq. (8) is a known type of PDE. Therefore, following (Polyanin, 2002), at the second step wemake a new change of variables z = xF ( φ ) , τ = Z φ F ( k ) dk, F ( φ ) = e R φ f ( k ) dk . (9)After this change the final PDE takes the form of Eq. (5) which finalizes the proof.Note, that Carr and Linetsky in (Carr and Linetsky, 2006) extend the time-dependent CEV modelconsidered in this paper by allowing a jump to zero. They also reduce their stock price process to a timehomogeneous Bessel process. We will discuss this extension as applied to the barrier options further inthis paper.As far as the terminal (now the initial ) condition in Eq. (3) and the boundary conditions in Eq. (4)in the new variables is concerned, we must distinguish two cases, which are determined by the sign of β .If − < β <
0, the domain of definition for z is z ∈ [0 , y ( τ )], where y ( τ ) = − β H − β ( t ( τ )) F ( φ ( τ )) > . (10)Accordingly, the initial condition now reads u (0 , z ) = e − R T r ( k ) dk "(cid:18) − βzF ( φ (0)) (cid:19) − /β − K + , (11)and the boundary conditions are u ( τ,
0) = u ( τ, y ( τ )) = 0 . (12) Page 4 of 32 emi-closed form solutions for barrier options ...
However, if 0 < β <
1, the left boundary goes to −∞ . Therefore, in this case it is convenient toredefine x → ¯ x = − x . This also redefines z → ¯ z = − z . Then the domain of definition for ¯ z becomes¯ z ∈ [ y ( τ ) , ∞ ) where y ( τ ) = 1 β H − β ( t ( τ )) F ( φ ( τ )) > . (13)The initial condition now reads u (0 , ¯ z ) = e − R T r ( k ) dk "(cid:18) β ¯ zF ( φ (0)) (cid:19) − /β − K + , (14)and the boundary conditions are u ( τ, ¯ z ) (cid:12)(cid:12)(cid:12) ¯ z →∞ = u ( τ, y ( τ )) = 0 . (15)Also, in the case 0 < β < z variables. It can be seen,that in this case the Up-and-out barrier options transforms to the Down-and-out counterpart.We interrupt here dealing with the time-dependent CEV model and postpone construction of thesolution of problems in Eq. (5), Eq. (11), Eq. (12) (or Eq. (5), Eq. (14), Eq. (15)) till Section 3. Instead, inthe next Section we consider the time-dependent CIR model, and barrier options written on a Zero-couponbond (ZCB) which follows the CIR model. This is done, because, as we show below, the correspondingPDE could also be transformed to Eq. (5). Thus, our proposed method could be applied uniformly toboth models. The Cox–Ingersoll–Ross (CIR) model has been invented in (Cox et al., 1985) for modeling interest rates. Inthis model the instantaneous interest rate r t is stochastic variable which follows the stochastic differentialequation (SDE), also named the CIR process. For the time-dependent version of the model this SDEreads dr t = κ ( t )[ θ ( t ) − r t ] dt + σ ( t ) √ r t dW t , r t =0 = r. (16)Here κ ( t ) > θ ( t ) is the mean-reversion level. The CIR modelis an extension of the Hull-White model that we analyzed in (Itkin and Muravey, 2020) by making thevolatility proportional to √ r t . This, on the one hand, allows avoiding the possibility of negative interestrates when the Feller condition 2 κ ( t ) θ ( t ) /σ ( t ) > F ( r, t, S )for this model is known in closed form. Here S is the bond expiration time. It is known, that F ( r, t, S )under a risk-neutral measure solves a linear parabolic partial differential equation (PDE), (Privault, 2012) ∂F∂t + 12 σ ( t ) r ∂ F∂r + κ ( t )[ θ ( t ) − r ] ∂F∂r = rF. (17)It should be solved subject to the terminal condition F ( r, S, S ) = 1 , (18)and the boundary condition F ( r, t, S ) (cid:12)(cid:12)(cid:12) r →∞ = 0 . (19)The second boundary condition is necessary in case the Feller condition is violated, so the interest rate r t can hit zero. Otherwise, the PDE in Eq. (17) itself at r = 0 serves as the second boundary condition. Page 5 of 32 emi-closed form solutions for barrier options ...
The ZCB price can be obtained from Eq. (17) assuming that the solution of can be represented in theform F ( r, t, S ) = A ( t, S ) e B ( t,S ) r . (20)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 = 1 + κ ( t ) B ( t, S ) − σ ( t ) B ( t, S ) , (21) ∂A ( t, S ) ∂t = − A ( t, S ) B ( t, S ) θ ( t ) κ ( t ) . To obey the terminal condition Eq. (18), the first PDE in Eq. (21) should be solved subject to the terminalcondition B ( S, S ) = 0, and the second one - to A ( S, S ) = 1.The first equation in Eq. (21) is the Riccati equation. It this general form it cannot be solvedanalytically for arbitrary functions κ ( t ) , σ ( t ), but can be efficiently solved numerically. Also, in somecases it can be solved approximately (asymptotically), see e.g., an example in (Carr and Itkin, 2020).Once the solution is obtained, the second equation in Eq. (21) can be solve analytically to yield A ( t, S ) = e − R tS B ( m ) θ ( m ) κ ( m ) dm . (22)When coefficients κ ( t ) , θ ( t ) , σ ( t ) are constants, it is known that the solution B ( t, S ) can be obtained inclosed form and reads, (Andersen and Piterbarg, 2010) B ( t, S ) = − S − t ) h ) − h + ( θ + h )[exp(( S − t ) h ) − , h = p θ + 2 σ . (23)Thus, B ( t, S ) < t < S . Therefore, F ( r, t, S ) → r → ∞ . In other words, the solution inEq. (20) satisfies the boundary condition at r → ∞ . In case when all the parameters of the model aredeterministic functions of time, and B ( t, S ) solves the first equation in Eq. (21), this also remains to betrue. This can be checked as follows. Since κ ( t ) > , σ ( t ) > B ( t, S ) = κ ( t ) σ ( t ) − σ ( t ) q κ ( t ) + 2[1 − B ′ ( t, S )] . Therefore, B ( t, S ) < B ′ ( t, S ) <
1. On the other hand, if B ( t, S ) < B ′ ( t, S ) < Let us consider a Down-and-Out barrier Call option written on a ZCB. Under a risk-neutral measure theoption price C ( t, r ) solves the same PDE as in Eq. (17), (Andersen and Piterbarg, 2010). ∂C∂t + 12 σ ( t ) r ∂ C∂r + κ ( t )[ θ ( t ) − r ] ∂C∂r = rC. (24)The terminal condition at the option maturity T ≤ S for this PDE reads C ( T, r ) = ( F ( r, T, S ) − K ) + , (25)where K is the option strike. Page 6 of 32 emi-closed form solutions for barrier options ...
By a standard contract, the lower barrier L F ( t ) (which we assume to be time-dependent as well) isset on the ZCB price, and not on the underlying interest rate r . This means that it can be written inthe form C ( t, r ) = 0 if F ( r, t, S ) = L F ( t ) . (26)However, since the ZCB price F ( r, t, S ) can be expressed in closed form in Eq. (20), this condition canbe translated into the r domain by solving the equation F ( r, t, S ) = A ( t, S ) e B ( t,S ) r = L F ( t ) , with respect to r . Denoting the solution of this equation as L ( t ) we find L ( t ) = 1 B ( t, S ) log (cid:18) L F ( t ) A ( t, S ) (cid:19) > , (27)where it is assumed that L F > A ( t, S ). Accordingly, in the r domain the boundary condition to Eq. (24)reads C ( t, L ( t )) = 0 . (28)The second boundary can be naturally set at r → ∞ . As at r → ∞ the ZCB price tends to zero, seeEq. (20), the Call option price also vanishes in this limit. This yields C ( t, r ) (cid:12)(cid:12)(cid:12) r →∞ = 0 . (29)The PDE in Eq. (24) can also be transformed to that for the Bessel process in Eq. (6). Proposition 2.1.
The Eq. (24) can be transformed to ∂u∂τ = 12 ∂ u∂z + bz ∂u∂z , (30) where b is some constant, u = u ( τ, z ) is the new dependent variable, and ( τ, z ) are the new independentvariables, if κ ( t ) θ ( t ) σ ( t ) = m , (31) where m ∈ [0 , ∞ ) is some constant. The Eq. (30) is the PDE associated with the one-dimensional Besselprocess in Eq. (6). Proof.
First make a change of variables C ( t, r ) = u ( t, z ) e a ( t ) r + R t a ( s ) κ ( s ) θ ( s ) ds , z = g ( t ) √ r, (32) g ( t ) = exp (cid:20) Z t (cid:16) κ ( s ) − a ( s ) σ ( s ) (cid:17) ds (cid:21) (33)where a ( t ) solves the Riccati equation da ( t ) dt = − σ ( t ) a ( t )2 + κ ( t ) a ( t ) + 1 . (34)This reduces the PDE in Eq. (24) to the form4 k ( t ) θ ( t ) − σ ( t )2 z ∂u∂z + 12 σ ( t ) ∂ u∂z + 4 g ( t ) ∂u∂t = 0 . (35) Page 7 of 32 emi-closed form solutions for barrier options ...
Now make a change of time τ ( t ) = 14 Z Tt g ( s ) σ ( s ) ds, (36)which transforms Eq. (35) to (cid:18) k ( t ) θ ( t ) σ ( t ) − (cid:19) z ∂u∂z + 12 ∂ u∂z = ∂u∂τ , t = t ( τ ) . (37)The function t ( τ ) is the inverse map of Eq. (36). It can be computed for any t ∈ [0 , T ] by substituting itinto the definition of τ , then finding the corresponding value of τ ( t ), and inverting.Finally, as by assumption k ( t ) θ ( t ) /σ ( t ) = m − const , we set b = m − /
2. Thus, the final PDE takesthe form of Eq. (5) which finalizes the proof.As follows from Proposition 2.1, for the time-dependent CIR model the transformation from Eq. (24)ro Eq. (30) cannot be done unconditionally. However, from practitioners’ points of view the conditionEq. (31) seems not to be too restrictive, Indeed, the model parameters already contain the independentmean-reversion rate κ ( t ) and volatility σ ( t ). Since m is an arbitrary constant, it could be calibrated tothe market data together with κ ( t ) and σ ( t ). Therefore, in this form the model should be capable forcalibration to the term-structure of interest rates.Also, according to the change of variables made in Proposition 2.1, the terminal condition Eq. (25)in new variables reads u (0 , z ) = e − a ( t (0)) z /g ( t (0)) − R T a ( s ) κ ( s ) θ ( s ) ds (cid:16) A ( T, S ) e B ( T,S ) z /g ( T ) − K (cid:17) + . (38)And the boundary conditions in Eq. (28) and Eq. (29) transform to C ( τ, y ( τ )) = 0 , C ( τ, z ) (cid:12)(cid:12)(cid:12) z →∞ = 0 , (39) y ( τ ) = 1 B ( t ( τ ) , S ( τ )) log (cid:18) L F ( t ( τ )) A ( t ( τ ) , S ( τ )) (cid:19) . For convenience of notation, further let us call as the CEV problem the PDE in Eq. (5) that has to besolved subject to the initial condition in Eq. (14) and the boundary conditions in Eq. (15). Also, we call the CIR problem the PDE in Eq. (30) solved subject to the initial condition Eq. (38) and the boundaryconditions in Eq. (39). Both problems can be considered simultaneously, as the PDE in Eq. (5) coincideswith that in Eq. (30), and the boundary conditions in Eq. (15) coincide with that in Eq. (39). Accordingly,both solutions are defined at the domain z ∈ [ y ( τ ) , ∞ ], thought the definitions of y ( τ ) in Eq. (13) andEq. (39) differ. We will describe our method for this domain in Section 3.1. Another type of the CEVproblem where ¯ z ∈ [0 , y − ( τ )] will be considered separately in Section 3.2.In (Carr and Itkin, 2020; Itkin and Muravey, 2020) a similar problem but for the heat equationwas solved by using two approaches. The first one is a method of generalized integral transform, activelyelaborated on by the Russian mathematical school to solve parabolic equations at the domain with movingboundaries, see (Kartashov, 1999, 2001) and references therein. These kind of problems are known inphysics for a long time and arise in the field of nuclear power engineering and safety of nuclear reactors;in studying combustion in solid-propellant rocket engines; in the theory of phase transitions (the Stefanproblem and the Verigin problem (in hydromechanics)); in the processes of sublimation in freezing andmelting; in the kinetic theory of crystal growth; etc. In (Carr and Itkin, 2020; Kartashov, 1999) thismethod was applied to the domain z ∈ [0 , y ( τ )], and in (Itkin and Muravey, 2020) for the first time the Page 8 of 32 emi-closed form solutions for barrier options ... solution was obtained for the semi-infinite domain z ∈ [ y ( τ ) , ∞ ). We will further develop this method tosolve the CEV and CIR problems in Section 4.The second method used to solve the same problems in (Carr and Itkin, 2020; Itkin and Muravey, 2020)is the method of heat potentials, see, e.g., (Tikhonov and Samarskii, 1963; Friedman, 1964.; Kartashov,2001) and references therein. 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 dePrado, 2020) (also see references therein). However, the CEV and CIR problems which we deal with inthis paper, cannot be reduced to the heat equation, but rather to the Bessel PDE. Therefore, in thisSection we propose generalization of the method for this type of equations. Accordingly, we call thisgeneralization as the method of Bessel potentials.Note, that the potential method could be constructed for any PDE where the space operator is a lineardifferential operator with time-independent coefficients. We propose and discuss this generalization ofthe heat potential method in Appendix. Thus, the heat and Bessel potentials are just two particularcases of this general scheme. y ( τ ) ≤ z < ∞ . Since both the CEV and CIR problems have the inhomogeneous initial condition, our first step is toreduce them to the alternative problems with a homogeneous initial condition. Since the Green functionof the Bessel equation at the infinite domain is known in closed form, (Polyanin, 2002), this can beachieved by representing u ( τ, z ) in the form u ( τ, z ) = q ( τ, z ) + Z ∞ y (0) u (0 , ξ ) q τ ( z, ζ, b ) dζ. (40)Here q τ ( z, ζ, b ) is the fundamental solution (or the transition density, or the Green function) of Eq. (30)at the domain z ∈ [0 , ∞ ). This density can be obtained assumed that the Bessel process stops when itreaches the origin. But since the domain of definition of z is z ∈ [ y ( τ ) , ∞ ), we moved the left boundaryfrom 0 to y (0).By the definition of b in Eq. (8), for the CEV model b = (1 + β ) / (2 β ). Since in this case 0 < β < b >
1. It is known, (Lawler, 2018; Linetsky and Mendoza, 2010), that in case b ≥ / q τ ( z, ζ, b ) is a good density with no defect of mass, i.e., it integrates into 1. The explicit representationreads, (Cox, 1975; Emanuel and Macbeth, 1982) q τ ( z, ζ, b ) = √ zζτ (cid:18) ζz (cid:19) b e − z ζ τ I b − / (cid:18) zζτ (cid:19) . (41)Here I ν ( x ) is the modified Bessel function of the first kind, (Abramowitz and Stegun, 1964).For the CIR model b = m − / m ∈ [0 , ∞ ) can be found by calibration. Therefore, if m > q τ ( z, ζ, b )is given by Eq. (41). Otherwise, 0 < m < − / < b < /
2. Then by another change of variables,(Polyanin, 2002) w ( τ, z ) = z − m ) u ( τ, z ) , the Eq. (30) transforms to the same equation with respect to w ( τ, z ) but now with b = (3 − m ) / < m < b > /
2, Therefore, again the Green function is represented byEq. (41).
Page 9 of 32 emi-closed form solutions for barrier options ...
The function q ( x, τ ) solves the problem ∂q ( τ, z ) ∂τ = 12 ∂ q ( τ, z ) ∂z + bz ∂q ( τ, z ) ∂z , (42) q (0 , z ) = 0 , y (0) < z < ∞ ,q ( τ, z ) (cid:12)(cid:12)(cid:12) z →∞ = 0 , q ( τ, y ( τ )) = ς ( τ ) ,ς ( τ ) = − Z ∞ y (0) u (0 , ζ ) q τ ( y ( τ ) , ζ, b ) dζ. This problem is like that in Eq. (30), Eq. (14), Eq. (15) Eq. (38), Eq. (39)), but now with a homogeneousinitial condition. Therefore, following the general idea of the method of heat potentials, we represent thesolution in the form of a generalized potential for the Bessel PDE q ( τ, z ) = Z τ Ψ( k ) ∂∂ξ " √ zξτ − k (cid:18) ξz (cid:19) b e − z ξ τ − k ) I b − / (cid:18) zξτ − k (cid:19) ξ → y ( k ) dk, (43)where Ψ( k ) is the potential density. It can be seen that q ( τ, z ) solves Eq. (30) as the derivative of theintegral on the upper limit is proportional to the Delta function which vanishes due to z = y ( τ ). Thesolution in Eq. (43) also satisfies the initial condition at τ = 0, and the vanishing condition at z → ∞ .For the large values of argument zξ/ ( τ − k ) we propose to use the following approximation : q ( τ, z ) ≈ √ π Z τ Ψ( k ) ∂∂ξ " √ τ − k (cid:18) ξz (cid:19) b e − ( z − ξ )22( τ − k ) ξ → y ( k ) dk (44)At the barrier z = y ( τ ) function q ( τ, z ) is discontinuous. Following a similar approach for the heatpotentials method, (Tikhonov and Samarskii, 1963)), it can be shown that the limiting value of q ( τ, z ) at z = y ( τ ) + 0 is equal to ς ( τ ): ς ( τ ) = Ψ( τ ) + Z τ Ψ( k ) ∂∂y ( k ) " p y ( τ ) y ( k ) τ − k (cid:18) y ( k ) y ( τ ) (cid:19) b e − y τ )+ y k )2( τ − k ) I b − / (cid:18) y ( τ ) y ( k ) τ − k (cid:19) dk. (45)The Eq. (45) is a linear Volterra equations of the second kind, (Polyanin and Manzhirov, 2008). Since ς ( τ )is a continuously differentiable function, Eq. (45) has a unique continuous solution for Ψ( τ ). The Volterraequation can be efficiently solved numerically, see (Itkin and Muravey, 2020) for the discussion on variousapproached to the numerical solution of this type of equations. In brief, for instance, the integral in theRHS is approximated using some quadrature rule with N nodes in k space, and the solution is obtainedat N nodes in the τ space. Thus, obtained matrix equation can be solved with the complexity O ( N )since the matrix is lower triangular. Since N could be small ( N ≈ − τ ) is found, the final solution reads u ( τ, z ) = Z τ Ψ( k ) ∂∂y ( k ) " p zy ( k ) τ − k (cid:18) y ( k ) z (cid:19) b e − z y k )2( τ − k ) I b − / (cid:18) zy ( k ) τ − k (cid:19) dk + Z ∞ y (0) u (0 , ξ ) q τ ( z, ζ, b ) dζ. (46) Page 10 of 32 emi-closed form solutions for barrier options ... < z < y ( τ ) . The construction of the solution in this case is similar to that described in the previous Section. Again,to obtain a PDE with a homogeneous initial condition, we represent the solution in the form u ( τ, z ) = q ( τ, z ) − ς ( τ ) + Z y (0)0 u (0 , ξ ) q τ ( z, ζ, b ) dζ, (47) ς ( τ ) = − Z y (0)0 u (0 , ζ ) q τ (0 , ζ, b ) dζ, where, (Abramowitz and Stegun, 1964) q τ (0 , ζ, b ) = 2 / − b ζ b τ b +1 / Γ (cid:16) b + (cid:17) e − ζ τ , and Γ( x ) is the Euler Gamma function.Then the function q ( x, τ ) solves the problem ∂q ( τ, z ) ∂τ = 12 ∂ q ( τ, z ) ∂z + bz ∂q ( τ, z ) ∂z , (48) q (0 , z ) = 0 , < z < y (0) ,q ( τ,
0) = 0 , q ( τ, y ( τ )) = ς ( τ ) + ς ( τ ) . We search for the solution in the form of the Bessel potential in Eq. (43). The potential density Ψ( τ )solves the following Volterra equation of the second kind ς ( τ ) + ς ( τ ) = Ψ( τ ) + Z τ Ψ( k ) ∂∂y ( k ) " p y ( τ ) y ( k ) τ − k (cid:18) y ( k ) y ( τ ) (cid:19) b e − y τ )+ y k )2( τ − k ) I b − / (cid:18) y ( τ ) y ( k ) τ − k (cid:19) dk. (49)Once this function is found, the final solution reads u ( τ, z ) = Z τ Ψ( k ) ∂∂y ( k ) " p zy ( k ) τ − k (cid:18) y ( k ) z (cid:19) b e − z y k )2( τ − k ) I b − / (cid:18) zy ( k ) τ − k (cid:19) dk + Z y (0)0 u (0 , ξ ) q τ ( z, ζ, b ) dζ. (50) Double barrier options for time-dependent models can be also priced by using the method of potentials. Inparticular, in (Itkin and Muravey, 2020) this is demonstrated for the time-dependent Hull-White model.Here we use a similar approach and apply the idea proposed in (Itkin and Muravey, 2020) to constructionof the semi-closed form solutions for double barrier options for the CIR and CEV models.Let us provide the explicit formulae just for the CEV model as for the CIR model this can be done inthe exact same way. Suppose we need the price of a double barrier Call option with the lower barrier L ( t )and the upper barrier H ( t ) > L ( t ). After doing transformation to the Bessel PDE as this is described inSection 1, this implies solving the following problem ∂u∂τ = 12 ∂ u∂z + bz ∂u∂z , (51) u ( τ = 0 , z ) = u (0 , z ) , y (0) < x < h (0) ,u ( y ( τ ) , τ ) = u ( h ( τ ) , τ ) = 0 , Page 11 of 32 emi-closed form solutions for barrier options ... where for − < β < y ( τ ) = − β L ( t ( τ )) − β F ( φ ( τ )) , h ( τ ) = − β H ( t ( τ )) − β F ( φ ( τ )) , (52)and for 0 < β < y ( τ ) = 1 β L ( t ( τ )) − β F ( φ ( τ )) , h ( τ ) = 1 β H ( t ( τ )) − β F ( φ ( τ )) . (53)Thus, in this case the solution is defined at the z -domain with two moving (time-dependent) boundaries.Since this problem has an inhomogeneous initial condition, the method of potentials cannot be directlyapplied. Therefore, similar to Eq. (40) we represent the solution in the form u ( τ, z ) = q ( τ, z ) + Z h (0) y (0) u (0 , ξ ) q τ ( z, ζ, b ) dζ. (54)Now the function q ( x, τ ) solves a similar problem but with the homogeneous initial condition ∂q∂τ = 12 ∂ q∂z + bz ∂q∂z , (55) q (0 , z ) = 0 , y (0) < x < h (0) ,q ( τ, y ( τ ) = − φ ( τ ) , q ( τ, h ( τ )) = − ψ ( τ ) ,φ ( τ ) = − Z h (0) y (0) u (0 , ξ ) q τ ( y ( τ ) , ζ, b ) dζ, ψ ( τ ) = − Z h (0) y (0) u (0 , ξ ) q τ ( h ( τ ) , ζ, b ) dζ, . Based on the method of (Itkin and Muravey, 2020), we construct the solution of Eq. (55) in the formof a generalized Bessel potential q ( τ, z ) = Z τ ( Ψ( k ) ∂∂ξ " √ zξτ − k (cid:18) ξz (cid:19) b e − z ξ τ − k ) I b − / (cid:18) zξτ − k (cid:19) ξ → y ( k ) (56)+Φ( k ) ∂∂ξ " √ zξτ − k (cid:18) ξz (cid:19) b e − z ξ τ − k ) I b − / (cid:18) zξτ − k (cid:19) ξ → h ( k ) ) dk. Here Ψ( k ) , Φ( k ) are the Bessel potential densities to be determined. Using the boundary conditions inEq. (55), and the fact that the expression in square brackets at τ = k is the Dirac delta function, onecan find that they solve a system of two Volterra equations of the second kind φ ( τ ) = Ψ( τ ) + Z τ ( Ψ( k ) ∂∂ξ " p y ( τ ) ξτ − k (cid:18) ξy ( τ ) (cid:19) b e − y τ )+ ξ τ − k ) I b − / (cid:18) y ( τ ) ξτ − k (cid:19) ξ → y ( k ) (57)+Φ( k ) ∂∂ξ " p y ( τ ) ξτ − k (cid:18) ξy ( τ ) (cid:19) b e − y τ )+ ξ τ − k ) I b − / (cid:18) y ( τ ) ξτ − k (cid:19) ξ → h ( k ) ) dk.ψ ( τ ) = Φ( τ ) + Z τ ( Ψ( k ) ∂∂ξ " p h ( τ ) ξτ − k (cid:18) ξh ( τ ) (cid:19) b e − h τ )+ ξ τ − k ) I b − / (cid:18) h ( τ ) ξτ − k (cid:19) ξ → y ( k ) (58)+Φ( k ) ∂∂ξ " p h ( τ ) ξτ − k (cid:18) ξh ( τ ) (cid:19) b e − h τ )+ ξ τ − k ) I b − / (cid:18) h ( τ ) ξτ − k (cid:19) ξ → h ( k ) ) dk. This system can be solved by various numerical methods with complexity O ( N ) (see the discussion afterEq. (45)). Once this is done, the solution of the double barrier problem is found. Page 12 of 32 emi-closed form solutions for barrier options ...
In this Section we solve the same problem but using the method of generalized integral transform. Asapplied to finance this method was successfully used in (Carr and Itkin, 2020; Itkin and Muravey, 2020)to price barrier options in the time-dependent OU model and American option for equites, and the Hull-White models for interest rates. The method was borrowed from physics, where it was used to solvethe Stefan problem and other heat and mass transfer problems with a moving boundary (or movinginterphase boundary), see (Kartashov, 1999, 2001) and references therein. In particular, in (Itkin andMuravey, 2020) the authors extended this approach to an infinite domain where the solution was notknown yet. Below we extend this approach and apply it to the CEV and CIR problems.Note, that so far, the method was elaborated on and used just for getting a semi closed form solutionof the heat equation. But in this paper we deal with the Bessel PDE. It (Kartashov, 1999) it is proposedto construct the direct integral transform for this equation by using Bessel functions. However, except thisrecommendation the explicit solution has not been presented. Moreover, our analysis showed that usingthe form of the transform proposed in (Kartashov, 1999)doesn’t give rise to the solution, as constructionof the inverse transform faces various technical problems.Therefore, our method presented in this Section is i) completely original and ii) gives rise to theclosed-form solution of the problem. In other words, we solve the CIR and CEV problems by using themethod of generalized integral transform to the very end, and to the best of the authors knowledge thisis done for the first time in the literature. As such, this approach could be also very useful in physicsfor solving various problems. As mentioned in (Kartashov, 1999) those problems appear (but not limitedto) in the field of nuclear power engineering and safety of nuclear reactors; in studying combustion insolid-propellant rocket engines; in laser action on solids; in the theory of phase transitions (the Stefanproblem and the Verigin problem (in hydromechanics)); in the processes of sublimation in freezing andmelting; in the kinetic theory of crystal growth; etc., see (Kartashov, 1999) and references therein. < z < y ( τ ) Recall, that based on the description in Section 1, this problem emerges when we consider the CEV prob-lem with β <
0. Since the Laplace transform of Eq. (30) gives rise to the Bessel ODE, (Abramowitz andStegun, 1964), it would be natural seeking for the general integral transform in the class of Bessel func-tions. Therefore, by analogy with (Kartashov, 2001; Carr and Itkin, 2020) we introduce the generalizedintegral transform of the form ¯ u ( τ, p ) = Z y ( τ )0 z ν +1 u ( τ, z ) J | ν | ( zp ) dz, (59)where p = a + i ω is a complex number, J ν ( x ) is the Bessel function of the first kind, and ν = 1 / (2 β ) < β <
0. Next, let us multiply both parts of Eq. (30) by z ν +1 J | ν | ( zp ) and integrate on z from 0 to y ( τ ). For the LHS this yields Z y ( τ )0 z ν +1 ∂u∂τ J | ν | ( zp ) dz = ∂ ¯ u∂τ − y ′ ( τ )[ y ( τ )] ν +1 u ( τ, y ( τ )) J | ν | ( y ( τ ) p ) . (60)The last term in the RHS of Eq. (60) vanishes due to the boundary condition in Eq. (12). Page 13 of 32 emi-closed form solutions for barrier options ...
Accordingly, for the RHS of Eq. (30) we have b = ν + 1 /
2, and J = 12 Z y ( τ )0 z ν +1 ∂ u∂z J | ν | ( zp ) dz = 12 z ν +1 ∂u∂z J | ν | ( zp ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y ( τ )0 − u ( τ, z ) ∂∂z (cid:16) z ν +1 J | ν | ( zp ) (cid:17) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y ( τ )0 (61)+ 12 Z y ( τ )0 u ( τ, z ) ∂ ∂z (cid:16) z ν +1 J | ν | ( zp ) (cid:17) dz,J = Z y ( τ )0 z ν +1 ν + 1 / z ∂u∂z J | ν | ( zp ) dz = ( ν + 1 / z ν J | ν | ( zp ) u ( τ, z ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y ( τ )0 − ( ν + 1 / Z y ( τ )0 u ( τ, z ) ∂∂z (cid:16) z ν J | ν | ( zp ) (cid:17) dz. Due to the boundary conditions in Eq. (12), the sum J + J can be represented as J + J = y ν +1 ( τ ) J | ν | ( y ( τ ) p )Ψ( τ ) + 12 Z y ( τ )0 u ( τ, z ) z " − νz ∂∂z (cid:16) z ν J | ν | ( zp ) (cid:17) + ∂ ∂z (cid:16) z ν J | ν | ( zp ) (cid:17) dz, Ψ( τ ) = ∂u∂z (cid:12)(cid:12)(cid:12) z = y ( τ ) . (62)It can be checked by using the theory of cylinder functions, (Bateman and Erdélyi, 1953) that theBessel function J | ν | ( zp ) also solves the following ordinary differential equation (ODE) d Zdz + 1 − νz dZdz + p Z = 0 , Z ( z ) = C z ν J | ν | ( zp ) + C z ν Y | ν | ( pz ) . (63)Here Y | ν | ( z ) denotes the Bessel function of the second kind (also known as the Neumann or Weberfunction) which is linearly independent of J | ν | ( z ). Assuming C = 1 , C = 0, from Eq. (60), Eq. (61) weobtain the following Cauchy problem for ¯ ud ¯ u ( τ, p ) dτ = 12 h − p ¯ u ( τ, p ) + y ν +1 ( τ ) J | ν | ( y ( τ ) p )Ψ( τ ) i , (64)¯ u ( p,
0) = Z y (0)0 z ν +1 J | ν | ( zp ) u (0 , z ) dz. The solution of this problem reads¯ u = e − p τ/ (cid:20) ¯ u (0 , p ) + 12 Z τ e p s/ Ψ( s ) y ν +1 ( s ) J | ν | ( y ( s ) p ) ds (cid:21) . (65)By analogy with (Carr and Itkin, 2020), we can obtain the Fredholm equation of the first type forthe function Ψ( τ ). For doing so, let us set p = iλ, λ ∈ R , tend τ to infinity and apply the formula J | ν | (i x ) = e i νπ/ I ν ( x ) which connects the Bessel function J ν ( x ) with the modified Bessel function I ν ( x ).This yields Z ∞ e − λ s/ Ψ( s ) y ν +1 ( s ) I ν ( y ( s ) λ ) ds = − Z y (0)0 q ν +1 I ν ( qλ ) u (0 , q ) dq. (66)The solution of this integral equation Ψ( τ ) can be found numerically on a grid by solving a systemof linear equations, see e.g., (Carr and Itkin, 2020) for the discussion on this subject and a numericalexample. Once the function Ψ( τ ) is found, it has to be substituted into Eq. (64) to obtain the generalizedtransform of u ( τ, z ) in the explicit form. Then, if this transform can be inverted back, we solved theproblem of pricing Up-and-Out barrier Call options for the CEV model with β < Page 14 of 32 emi-closed form solutions for barrier options ...
As it was already mentioned, it is reasonable to seek the solution of the CEV problem in the class of theBessel functions. Therefore, we represent the solution in the form u ( τ, z ) = z − ν ∞ X n =1 α n ( τ ) J | ν | (cid:18) µ n zy ( τ ) (cid:19) (67)Here µ n is an ordered sequence of the positive zeros of the Bessel function J | ν | ( µ ): J | ν | ( µ n ) = J | ν | ( µ m ) = 0 , µ n > µ m > , n > m. Note, that the definition in Eq. (67) automatically respects the vanishing boundary conditions for u ( τ, z ).We assume that this series converges absolutely and uniformly ∀ z ∈ [0 , y ( τ )] for any τ > z → ˆ z = zy ( τ ) yields¯ u ( τ, p ) y ( τ ) = ∞ X n =1 α n ( τ ) Z ˆ zJ | ν | ( µ n ˆ z ) J | ν | ( py ( τ )ˆ z ) d ˆ z. (68)The set of functions J | ν | ( α ˆ z ) with α ∈ µ n , n = 1 , . . . , forms an orthogonal basis in the space C [0 ,
1] withthe scalar product h J | ν | ( αz ) , J | ν | ( βz ) i = 2 Z zJ | ν | ( αz ) J | ν | ( βz ) dzJ | ν | +1 ( α ) J | ν | +1 ( β ) = ( , α = β, , α = β (69)Therefore, the explicit formula for each coefficient α n ( τ ) is straightforward α n ( τ ) = 2 ¯ u ( τ, µ n /y ( τ )) y ( τ ) J | ν | +1 ( µ n ) . (70)Thus, the final solution for u ( τ, z ) reads u ( τ, z ) = 2 z − ν y ( τ ) ∞ X n =1 (cid:20) Z y (0)0 u (0 , s ) s ν +1 e − µ n y τ ) τ J | ν | ( µ n s/y ( τ )) J | ν | ( µ n z/y ( τ )) J | ν | +1 ( µ n ) ds (71)+ 12 Z τ y ν +1 ( s )Ψ( s ) e − µ n y τ ) ( τ − s ) J | ν | ( µ n y ( s ) /y ( τ )) J | ν | ( µ n z/y ( τ )) J | ν | +1 ( µ n ) ds (cid:21) . This expression can be also re-written in the form u ( τ, z ) = 2 z − ν y ( τ ) h Z y (0)0 s u (0 , s )Θ | ν | √ τy ( τ ) , sy ( τ ) , zy ( τ ) ! ds (72)+ 12 Z τ y ( τ )Ψ( s )Θ | ν | √ τ − sy ( τ ) , y ( s ) y ( τ ) , zy ( τ ) ! ds i , where we introduced a new functionΘ | ν | ( θ, x , x ) = ∞ X n =1 e − µ nθ ( x x ) ν J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n ) . (73)The function Θ | ν | ( θ, x , x ) is an analog of the Jacobi theta function which is a periodic solution of theheat equation. Indeed, in (Carr and Itkin, 2020) the solution of a similar problem for the time-dependent Page 15 of 32 emi-closed form solutions for barrier options ...
OU model (so β = − , ν = − / | ν | = 1 /
2) with moving boundaries but for the heat equation hasbeen obtained in terms of the theta functions. It can be checked that, if ν = 1 /
2, we haveΘ / √ τy ( τ ) , sy ( τ ) , zy ( τ ) ! = y ( τ )2 (cid:20) θ (cid:18) e − π θ , nπ ( s − z )4 y ( τ ) (cid:19) − θ (cid:18) e − π θ , nπ ( s + z )4 y ( τ ) (cid:19)(cid:21) , (74)where θ ( ω, x ) is the Jacobi theta function, (Mumford et al., 1983). Accordingly, function Θ | ν | ( θ, x , x )is a periodic solution of the Bessel equation.Alternatively to the Fredholm equation of the first kind in Eq. (66) which is ill-posed and requiresspecial methods to solve it, see (Carr and Itkin, 2020), we can use a trick proposed in (Itkin and Muravey,2020) and instead derive the Volterra equation of the second kind for the function Ψ( τ ). For doing that,one needs to differentiate Eq. (71) on z , and then let z = y ( τ ). This yieldsΨ( τ ) = − y ν ( τ ) ∞ X n =1 ( µ n + ν ) (cid:20) Z y (0)0 u (0 , s ) s ν +1 e − µ nτ y τ ) J | ν | ( µ n s/y ( τ )) J | ν | +1 ( µ n ) ds (75)+ 12 Z τ y ν +1 ( s )Ψ( s ) e − µ n ( τ − s )2 y τ ) J | ν | ( µ n y ( s ) /y ( τ )) J | ν | +1 ( µ n ) ds (cid:21) . This equation has to be solved numerically, again see (Itkin and Muravey, 2020) for the discussion andexamples.
In some cases, Eq. (75) can be solved asymptotically. For instance, one can apply the following approxi-mations J | ν | ( z ) + J | ν | +1 ( z ) ≈ πz , x ≫ ν, J | ν | +1 ( µ n ) ≈ πµ n , µ n ≫ ν (76) J | ν | ( z ) ∼ r πz cos (cid:18) z − | ν | + 14 π (cid:19) , z → ∞ µ n ≈ π (cid:18) n + 2 | ν | + 14 (cid:19) , n → ∞ . Then the infinite sum in Eq. (71) can be truncated up to keep first N terms. The reminder (the error ofthis method) reads R ( N, τ, z ) = 1 z ν +1 / y ( τ ) ∞ X n = N +1 (cid:26) (77) Z y (0)0 s ν +1 / u (0 , s ) e − π λ n y τ ) τ (cid:20) cos (cid:18) πλ n ( s + z ) y ( τ ) − πλ n (cid:19) + cos (cid:18) πλ n ( s − z ) y ( τ ) (cid:19)(cid:21) ds + 12 Z τ y ν +1 / ( s )Ψ( s ) e − π λ n y τ ) ( τ − s ) (cid:20) cos (cid:18) πλ n ( y ( s ) + z ) y ( τ ) − πλ n (cid:19) + cos (cid:18) πλ n ( y ( s ) − z ) y ( τ ) (cid:19)(cid:21) ds (cid:27) . Here λ n = n + (2 ν − /
4A simple assessment shows that, since − < β <
0, from Eq. (8) at typical values of the modelparameters we have f ( t ) ≈ f ≪
1. Assume that H = const . Hence, in Eq. (9) φ = log(2 βf τ + 1)2 βf , y ( τ ) = − H − β √ βf τ + 1 β . Then it can be checked that the expression − µ n τ y ( τ ) Page 16 of 32 emi-closed form solutions for barrier options ... rapidly drops down for − < β < , τ >
0, unless τ ≪ H ≫
1. Therefore, in this case just firstfew terms in the sum in Eq. (75) would be a good approximation.Another approximation can be proposed to compute function Θ | ν | ( θ, x , x ) defined in Eq. (73). Theidea is that, as mentioned in above if τ ≪ H = O (1), the first argument θ of this function is small, θ ≪
1. Then the function c = e − µ nθ is small at large n , and µ n θ is small at small n . Therefore, forsmall n we represent c by using the Padé approximation ( k, n we replace the smallvalues of c with the same Padé approximation. To estimate how accurate is this trick, let us pick, forinstance, k = 2, so c ≈
11 + x (cid:18) − x (cid:19) + O ( x ) , x = µ n θ . (78)Thus, expanding the parenthesis into two terms, the function ( x x ) − ν Θ | ν | ( θ, x , x ) can be representedas ( x x ) − ν Θ | ν | ( θ, x , x ) = A ( θ, x , x ) + A ( θ, x , x ) . (79)Now observe that, (Bateman and Erdélyi, 1953) A ( θ, x , x ) = ∞ X n =1 J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n ) 11 + µ n θ = − θ ∞ X n =1 J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n )( γ − µ n ) (80)= πJ | ν | ( x γ ) θ J | ν | ( γ ) h J | ν | ( x γ ) Y | ν | ( γ ) − J | ν | ( γ ) Y | ν | ( x γ ) i = 2 I | ν | ( x ¯ γ ) θ I | ν | (¯ γ ) h K | ν | ( x ¯ γ ) I | ν | (¯ γ ) − K | ν | (¯ γ ) I | ν | ( x ¯ γ ) i , γ = 2i /θ, ¯ γ = 2 /θ. Here we again used the formulae J ν ( ix ) = e i νπ/ I ν ( x ) , Y ν ( ix ) = e ( ν +1) π i / I ν ( x ) − /πe − νπ i / K ν ( x ) , which connect the Bessel functions J ν ( x ) and Y ν ( x ) with the modified Bessel functions I ν ( x ) and K ν ( x ).To compute the next term in Eq. (78) A ( θ, x , x ) = − θ ∞ X n =1 J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n ) µ n µ n θ (81)let us differentiate the LHS of Eq. (80) by θ , so ∂∂θ ∞ X n =1 J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n ) 11 + µ n θ = − θ ∞ X n =1 J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n ) µ n (cid:16) µ n θ (cid:17) (82) ≈ − θ ∞ X n =1 J | ν | ( µ n x ) J | ν | ( µ n x ) J | ν | +1 ( µ n ) µ n µ n θ . Therefore, the second term in Eq. (78) takes the form A ( θ, x , x ) = − ¯ γ ∂∂ ¯ γ ( ¯ γ I | ν | ( x ¯ γ ) I | ν | (¯ γ ) h K | ν | ( x ¯ γ ) I | ν | (¯ γ ) − K | ν | (¯ γ ) I | ν | ( x ¯ γ ) i) , ¯ γ = 2 √ /θ, (83) ∂I | ν | ( γ ) ∂γ = 12 (cid:16) I | ν |− ( γ ) + I | ν | +1 ( γ ) (cid:17) , ∂K | ν | ( γ ) ∂γ = 12 (cid:16) K | ν |− ( γ ) + K | ν | +1 ( γ ) (cid:17) . Page 17 of 32 emi-closed form solutions for barrier options ...
A similar approximation can be developed for Eq. (75) by using the identity, (Bateman and Erdélyi,1953) ∞ X n =1 µ n J | ν | ( µ n x )( µ n − k ) J | ν | +1 ( µ n ) = J | ν | ( kx )2 J | ν | ( k ) . (84)As an example, let consider the CEV model with constant parameters given in Table 1. In Fig. 1 we Table 1:
Parameters of the test. r q σ H N x /y ( τ ) x /y ( τ )0.02 0.01 0.5 0.2 100 0.5 0.5present the difference of two values of the function G ν ( x , x ) = ( x x ) − ν Θ | ν | ( θ, x , x ): one computedby the definition in Eq. (73) using the first N terms in the sum; and the other computed by using theapproximation in Eq. (79). It can be seen that the latter approximation provides an accuracy about 25%except the area where both τ and β are simultaneously kind of large, and hence, the value of G ν ( x , x )is very small. Figure 1:
Comparison of the exact and approximated values of the function G ν ( x , x ) . This approach can be further developed by increasing k in the Padé approximations ( k, θ , and then express this derivative via the RHS of Eq. (80), etc. Let us consider the following stopping moment T yx = inf { τ ≥ , X τ ≥ y ( τ ) } , Page 18 of 32 emi-closed form solutions for barrier options ... where X τ is the Bessel process defined by Eq. (6) and originated from X = x . From standard results inprobability theory the p.d.f. ρ yx ( τ ) of the moment T yx can be found via the Fokker–Planck–Kolmogorovequation associated with the process X τ :12 ∂ F x ( τ, x ) ∂x − bx ∂F x ( τ, x ) ∂x + bx F x ( τ, x ) = ∂F x ( τ, x ) ∂τ (85) F x ( τ, y ( τ )) = 0 , F x ( τ,
0) = δ ( x − x )here δ ( x − x ) is the Dirac measure at the point x . The density ρ yx ( τ ) has the following representation: ρ yx ( τ ) = 12 ∂F x ( τ, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = y ( τ ) . (86)The solution of the problem Eq. (30) and the function Ψ( τ ) can be represented in terms of the function F : u ( z, τ ) = z − b Z y (0)0 s b u (0 , s ) F s ( τ, z ) ds, Ψ( τ ) = 2 y ( τ ) − b Z y (0)0 s b u (0 , s ) ρ ys ( τ ) ds (87)For the boundary moving linearly in time y ( τ ) ≈ α + βτ it is possible to propose the following approxi-mation for the function Ψ( τ )Ψ( τ ) ≈ y ( τ ) − (2 ν +1) Z y (0)0 s ν +1 u ( s, e β α ( α − s ) + b τ ( α + βτ ) ν +2 ∞ X n =1 s − ν µ n J | ν | ( µ n z/α ) α − ν J | ν | +1 ( µ n e − µ n τ α ( α + βτ ) ds. (88)This is due to the fact that the p.d.f ρ yx ( τ ) of the first hitting time of the line aτ + b for the Bessel processis known explicitly, (Alili and Patie, 2010).Note, that a similar problem for the OU process with both constant and time-dependent coefficientshas been studied in (Lipton and Kaushansky, 2020a). The authors considered the first hitting timedensity to a moving boundary for a diffusion process, which satisfies the Cherkasov condition, and hence,can be reduced to a standard Wiener process. They give two complementary (forward and backward)formulations of this problem and provide semi-analytical solutions for both by using the method of heatpotentials. z > y ( τ ) . To recall, this problem occurs in both the CEV model with 0 < β < u ( τ, p ) = Z ∞ y ( τ ) z ν +1 W ( τ, p, z ) u ( τ, z ) dz (89) u ( τ, z ) = z − ν Z ∞ pW ( τ, p, z ) V ( τ, p ) ¯ u ( τ, p ) dp. The kernel W ( a, b ) and the function V ( p ) are defined as follows, (Bateman and Erdélyi, 1953) W ( τ, a, b ) = J | ν | ( ab ) Y | ν | ( ay ( τ )) − Y | ν | ( ab ) J | ν | ( ay ( τ )) , (90) V ( τ, p ) = J | ν | ( py ( τ )) + Y | ν | ( py ( τ )) . The definitions in Eq. (90) are generalizations of the Pythagorean and Angle sum identities for trigono-metric functions to the case of cylinder functions J | ν | and Y | ν | . The functions W ( τ, a, b ) as the functionsof the second argument a also form an orthogonal basis in the space C [ y ( τ ) , ∞ ) for all τ > Page 19 of 32 emi-closed form solutions for barrier options ...
However, we cannot apply this transform directly to the Bessel equation due to fact that the kernelis time-dependent. Therefore, we propose to represent ¯ u as the weighted sum of the following transforms¯ u J ( p, τ ) = Z ∞ y ( τ ) z ν +1 J | ν | ( zp ) u ( τ, z ) dz, ¯ u Y ( p, τ ) = Z ∞ y ( τ ) z ν +1 Y | ν | ( zp ) u ( τ, z ) dz, (91)so ¯ uτ, p ) = ¯ u J ( τ, p ) Y | ν | ( y ( τ ) p ) − ¯ u Y ( τ, p ) J | ν | ( y ( τ ) p ) (92)The explicit formulae for ¯ u J and ¯ u Y read¯ u J ( τ, p ) = e − p τ/ (cid:20) ¯ u J (0 , p ) + 12 Z τ e p s/ Ψ( s ) y ν +1 ( s ) J | ν | ( y ( s ) p ) ds (cid:21) , (93)¯ u Y ( τ, p ) = e − p τ/ (cid:20) ¯ u Y (0 , p ) + 12 Z τ e p s/ Ψ( s ) y ν +1 ( s ) Y | ν | ( y ( s ) p ) ds (cid:21) , where Ψ( τ ) is defined in Eq. (62). Using the inversion formula from Eq. (89) we immediately get theexplicit formula for u ( τ, z ) : u ( z, τ ) = z − ν Z ∞ Z ∞ y (0) s ν +1 u (0 , s ) e − p τ W ( τ, p, z ) W ( τ, p, s ) V ( τ, p ) p dp ds (94)+ z − ν Z ∞ Z τ y ν +1 ( s )Ψ( s ) e − p ( τ − s ) W ( τ, p, z ) W ( τ, p, y ( s )) V ( τ, p ) p dp ds. By analogy with the previous section, the function Ψ( τ ) solves the Fredholm equation of the first kind Z ∞ e − λ s/ Ψ( s ) y ν +1 ( s ) K ν ( y ( s ) λ ) ds = − Z ∞ y (0) q ν +1 K ν ( qλ ) u (0 , q ) dq, (95)or the Volterra equation of the second kindΨ( τ ) = y − ν ( τ ) (cid:20) Z ∞ Z ∞ y (0) q ν +1 u (0 , q ) e − p τ Q ( τ, p ) W ( τ, p, q ) V ( τ, p ) dp dq (cid:21) (96)+ 12 Z ∞ Z τ y ν +1 ( s )Ψ( s ) e − p ( τ − s ) Q ( τ, p ) W ( τ, p, y ( s )) V ( τ, p ) dp ds. Here Q ( τ, p ) = J ν +1 ( py ( τ )) Y ν ( py ( τ )) − Y ν +1 ( py ( τ )) J ν ( py ( τ )) . (97)In some cases, Eq. (94) can be solved asymptotically. For instance, one can apply the followingapproximations J | ν | ( z ) + Y | ν | ( z ) ≈ πz ∞ X k =0 (2 k − k z k Γ( | ν | + k + 1 / k !Γ( | ν | − k + 1 / , J | ν | ( z ) + Y | ν | ( z ) ∼ πz , z → ∞ (98) Y | ν | ( z ) ∼ r πz sin (cid:18) z − | ν | + 14 (cid:19) , z → ∞ , W ( τ, a, b ) ∼ a [ b − y ( τ )]) πa p by ( τ ) , a → ∞ . Then, the outer integral in Eq. (94) can be truncated from above and approximated by the integral overthe domain [0 , P ] , < P < ∞ . The reminder (the error of this method) reads R ( P, τ, z ) = z − ν − / Z ∞ P (cid:26) Z ∞ y (0) s ν +1 / u (0 , s ) e − p τ sin ( p [ z − y ( τ )]) sin ( p [ s − y ( τ )]) ds (99)+ 12 Z τ y ν +1 / ( s )Ψ( s ) e − p ( τ − s ) sin ( p [ z − y ( τ )]) sin ( p [ y ( s ) − y ( τ )]) ds (cid:27) dp. Page 20 of 32 emi-closed form solutions for barrier options ...
Introducing the new function Υ(
P, t, η )Υ(
P, t, η ) = e − η t (cid:20) erfc (cid:18) P t + i η √ t (cid:19) + erfc (cid:18) P t − i η √ t (cid:19)(cid:21) , and taking into account the identity Z ∞ P e − p t cos( ηp ) dp = 12 r π t Υ( P, t, η ) , we obtain the following explicit representation for R ( P, τ, z ) R ( P, τ, z ) = z − ν − / √ π √ (cid:26) Z ∞ y (0) s ν +1 / u (0 , s ) √ τ [Υ( P, τ, z − s ) − Υ( P, τ, z + s − y ( τ ))] ds (100)+ 12 Z τ y ν +1 / ( s )Ψ( s ) √ τ − s [Υ( P, τ − s, z − y ( s )) − Υ( P, τ − s, z + y ( s ) − y ( τ ))] ds (cid:27) . Now we show that under the assumptions Z ∞ y (0) q ν +1 / u (0 , q ) dq ≤ M , Z ∞ y (0) q ν +1 u (0 , q ) dq ≤ M , M , M − const, (101)we can set the upper limit of integration P such that | R ( P, τ, z ) | < ǫ for any ǫ >
0. Indeed, using thefollowing inequalities (see Eq. (87)):1 √ t Υ( P, t, η ) ≤ erfc( P ) , Ψ( τ ) ≤ y ( τ ) − (1+2 ν ) Z ∞ y (0) q ν +1 u (0 , q ) dq, (102)yields | R ( P, τ, z ) | ≤ z − ν − / √ π erfc( P )2 (cid:26) Z ∞ y (0) q ν +1 / u (0 , q ) dq + Z τ y − ν − / ( s ) ds Z ∞ y (0) q ν +1 u (0 , q ) dq (cid:27) . (103)Since the second integral is bounded for any τ , we obtain the following inequality: | R ( P, τ, z ) | ≤ z − ν − / √ π erfc( P )2 ( M + M M ( τ )) , M ( τ ) = Z τ y − ν − / ( s ) ds. (104) Similar to (Itkin and Muravey, 2020), to check performance and accuracy of the proposed methods weconstruct the following test. We consider Up-and-Out Barrier Call option written on the underlyingstock which follows the CEV process with 0 < β <
1. This case is described at the end of Section 1. Torecall, after the change of variables proposed in that Section is done, the problems is transformed to thesolution of the Bessel PDE at the domain ¯ z ∈ [ y ( τ ) , ∞ ) with the boundary and initial conditions givenin Eq. (15), Eq. (14) In this test we use the explicit form of parameters r ( t ) , q ( t ) , σ ( t ) r ( t ) = r − r k ( a + t ) , q ( t ) = q − q k ( a + t ) , σ ( t ) = σ √ a + t, (105) Hence, in new variables the Up-and-Out option transforms to the Down-and-Out option. Page 21 of 32 emi-closed form solutions for barrier options ... where r , q , σ , r k , q k , σ k are constants. We also assume r = q , and H − const . With these definitionsone can find φ ( t ) = − σ ( t − T )(2 a + t + T ) , φ ( τ ) = σ β ( q k − r k ) log σ + 2 β ( q k − r k ) τσ ! , (106) F ( φ ) = q βτ ( q k − r k ) + σ σ , y ( τ ) = F ( φ ) H − β /β. We approach pricing the Up-and-Out barrier Call option in the CEV model twofold. First, as abenchmark we solve the PDE in Eq. (2) by using a finite-difference (FD) scheme of the second order inspace and time. We use the Crank-Nicolson scheme with few first Rannacher steps on a non-uniform gridcompressed close to the barrier level, see (Itkin, 2017). Accordingly, our domain in S space is S ∈ [0 , H ] . Alternatively, we apply the method of Bessel potentials (BP) developed in Section 3.1 to solve theBessel PDE Eq. (5). For doing so, first we solve the Volterra equation in Eq. (45) where the kernel isapproximated on a rectangular grid M × M , and the integral is computed using the trapezoidal rule.This implies solving the following system of linear equations k ς k = ( I + P ) k Ψ k . (107)Here k Ψ k is the vector of discrete values of Ψ( τ ) , τ ∈ (cid:20) , τ ( t ) (cid:12)(cid:12)(cid:12) t =0 (cid:21) on a grid with M nodes, k ς k is asimilar vector of ς ( τ ), I is the unit M × M matrix, and P is the M × M matrix of the kernel values onthe same grid. Note, that the matrix P is lower triangular. Therefore, solution of Eq. (107) can be donewith complexity O ( M ).As the kernel (and so the matrix P ) doesn’t depend of strikes K , but only the function ς ( τ ), Eq. (107)can be solved simultaneously for all strikes by inverting the matrix I + P with the complexity O ( M ),and then multiplying it by vectors k ς k k , k = 1 , . . . , ¯ k , ¯ k is the total number of strikes. Therefore, thetotal complexity of this step remains O (¯ kM ), but this operation, however, can be vectorized in k . TheVolterra equation could also been solved by iterative methods, but with almost the same complexity, seediscussion in (Itkin and Muravey, 2020). Table 2:
Parameters of the test. r q σ r k q k σ k a H S T ∈ [1 / , . , . ,
1] and strikes K ∈ [59 , , , , , The same results computed by using the BP method are displayed in Fig. 2 for the option prices. Also,in Fig. 3 the percentage difference between the prices obtained by using the BP and FD methods is A similar approach can be developed for the Down-and-Out options with S ∈ [ L, ∞ ). Then instead of truncating theinfinite semi-interval, one can transform it to the fixed interval [ − ,
1) or to [0 ,
1) and solve a modified PDE on the newinterval. The boundary behavior of the solution can be obtained using Fichera theory and/or Greeen’s integral formula, see(Wilmott et al., 2014) where this was done in many cases and proved that this approach works well. Page 22 of 32 emi-closed form solutions for barrier options ... presented as a function of the option strike K and maturity T . Here to provide a comparable accuracywe run the FD solver with 101 nodes in space S and 100 steps in time t . Otherwise the quality of theFD solution is not sufficient. Table 3:
Up-and-Out barrier Call option prices computed by using the BP and FD methods. BP FD Difference %K\T 0.0833 0.3 0.5 1.0 0.0833 0.3 0.5 1.0 0.0833 0.3 0.5 1.059 Figure 2:
Up-and-Out barrier Call option price computed by using the BP method.
It can be seen that the agreement of both methods is good (less than 1%) if the option price is not toosmall which happens when the strike K is close to the barrier or at high maturities. In this case, as thisis seen from Table 3, the relative difference becomes large, but the absolute difference of two methods isabout one cent, which is almost insignificant. Obviously, such cases are a challenge for any FD method,as at t = T there is a jump in the initial condition at the boundary, and the first derivative of the solutiondoesn’t exists in this point.As far as performance of both methods is concerned, to decrease the elapsed time for the FD methodinstead of Eq. (2) we solve the corresponding forward PDE. Therefore, prices of all options for a givenset of strikes and maturities could be obtained within one sweep. This also requires ¯ m × ¯ k integrations Page 23 of 32 emi-closed form solutions for barrier options ...
Figure 3:
Percentage difference of Up-and-Out barrier Call option prices computed by usingthe BP and FD methods. of the product of thus found density function with the payoff function, where ¯ m is the total number ofmaturities, and ¯ k is the total number of strikes. In this test the elapsed time for the FD method is, onaverage, 140 msec.For the BP method, since the expression for ς ( τ ) in Eq. (42) is not known in closed form, we computethis integral numerically by using the Simpson quadratures. Nevertheless, to make the results accurate,we need to increase the number of the grid nodes M . As compared with (Itkin and Muravey, 2020), wherea similar expression for the Hull-White model could be computed in closed form, and M = 20 provideda sufficient accuracy, here we need to take M = 100. Nevertheless, the elapsed time, on average, is 70msec, i.e. twice faster than the FD method for the forward equation.Decreasing M almost doesn’t impact the accuracy of the method at large values of the options prices,while slightly deteriorates the quality of the feed at large maturities and large strikes. Changing M from100 to 70 drops down the elapsed time to 50 msec. However, for the FD scheme decreasing the gridto 50 ×
50 drops down the elapsed time to 30 msec while almost twice increasing the error for smallmaturities. Overall, we can conclude that the method of BP demonstrates, at least, same performanceas the forward FD solver.
Here we solve the same problem by using the GIT method developed in Section 3.1. Our numericalscheme is similar to that for the BP method: first we solve the Volterra equation Eq. (96) and thencompute the value of the integrals in Eq. (94) using a trapezoidal rule. The inner integrals in the firstsummand in Eq. (96) and Eq. (94) can be computed explicitly for the payoff Eq. (11) via the followingformulas (Gradshtein and Ryzhik, 2007) Z x ν +1 J ν ( ax ) dx = a − J ν +1 ( a ) , ℜ ( ν ) > − , Page 24 of 32 emi-closed form solutions for barrier options ... Z x − ν J ν ( ax ) dx = a ν − ν − Γ( ν ) − a − J ν − ( a ) , ℜ ( ν ) < , Z x ν +1 Y ν ( ax ) dx = a − Y ν +1 ( a ) + 2 ν +1 a − ν − Γ( ν + 1) , ℜ ( ν ) > − , Z x − ν Y ν ( ax ) dx = a ν − cot( νπ )2 ν − Γ( ν ) − a − Y ν − ( a ) , ℜ ( ν ) < . However, for the second summands we have to numerically compute two-dimensional integrals containinga lot of special functions.We run the same test described in above, and the results of this numerical experiment are presentedin Table 4 and also in Fig. 3 which depicts the percentage difference between the prices obtained by usingthe GIT and FD methods. For this test we use M = 10 steps in time. This algorithm was implementedin python. Table 4:
Up-and-Out barrier Call option prices computed by using the GIT and FD methods.
GIT
FD Difference %K\T 0.0833 0.3 0.5 1.0 0.0833 0.3 0.5 1.0 0.0833 0.3 0.5 1.059 Figure 4:
Percentage difference of Up-and-Out barrier Call option prices computed by usingthe FD and GIT methods. -4 -2
75 0.8
Percentage difference of the GIT and FD solutions K T Page 25 of 32 emi-closed form solutions for barrier options ...
It can be seen that this method produces very accurate results at high strikes and maturities (i.e.where the option price is relatively small) in contrast to the BP method. This can be verified by lookingat the exponents in Eq. (94) which are proportional to the time τ . Contrary, when the price is higher(short maturities, low strikes) the GIT method is slightly less accurate than the BP method, as in Eq. (45)the exponent is inversely proportional to τ . Obviously, the accuracy of the GIT method increases when M increases.This situation is well investigated for the heat equation with constant coefficients. As applied topricing double barrier options, it is described in (Lipton, 2002). There exist two representation of thesolution: one - obtained by using the method of images, and the other one - by the Fourier series. Despiteboth solutions are equal as the infinite series, their convergence properties are different. In particular,the Fourier method is superior when the difference between the upper H and lower L barriers is small.and the time is relatively large. And the image expansion should be used otherwise.In this paper we come to a similar principle for the time-dependent problems, and not just for theheat equation but also for the Bessel one. Thus, it is important that both the BP and GIT methodsdon’t duplicate but rather compliment each other.The speed of our python implementation is a bit slower than that for BP (approx. 0.158 sec). However,the latter method was implemented in Matlab. It is known that linear algebra in python (numpy) is almost3 times slower than that in Matlab. Therefore, performance of both the BP and GIT methods is roughlysame.In addition, one can find that the main computational time is spent by a lot of calls to the routinecomputing the values of the integrands. The integrand in Eq. (94) is a four-dimensional function of z , p , y ( τ ) and y ( s ). Let us denote this functions as Z ( z, p, v, w ), where v = y ( τ ) and w = y ( s ). We can gainthe speed by the following trick: first pre-compute the values of Z on a regular four-dimensional grid,and then design this computational routine as interpolation.It is also worth mentioning that in many situations the parameters of the model are such that theboundary y ( τ ) changes slowly with time, i.e. y ( τ ) is almost const. Then the first integral in Eq. (94) isa good approximation of the price. Accordingly, one don’t need to solve the Volterra equation Eq. (96)that makes the algorithm about 2.5 faster. We have already mentioned that both the method of heat potentials (which in our case is extended tothe method of Bessel potentials) and the method of generalized integral transforms were first developedin physics, see (Kartashov, 1999, 2001; Tikhonov and Samarskii, 1963; Friedman, 1964.) and referencestherein. This was done to solve various problems of heat and mass transfer which are widely present inphysics, chemistry, energetic, nuclear engineering, geology and many other areas of science and engineer-ing. It turns out that many of those problems mathematically can be formulated in terms of stationaryand non-stationary heat transfer. This includes such problems as diffusion, sedimentation, viscosity flowsaccompanied by various kinetic processes, astronomy, atomic physics, absorbtion, combustion, phasetransitions and many others.As an example in this Section we consider the Stefan problem which is a particular kind of theboundary value problem for the heat equation adapted to the case in which a phase boundary can movewith time. It was introduced by I. Stefan in 1889 (see a detailed review in (Lyubov, 1978)). The classicalStefan problem describes the temperature distribution in a homogeneous medium undergoing a phasechange, for example ice passing to water. This is accomplished by solving the heat equation imposing theinitial temperature distribution on the whole medium, and a particular boundary condition, the Stefancondition, on the evolving boundary between its two phases. Note that this evolving boundary is anunknown hyper-surface; hence, Stefan problems are examples of free boundary problems. However, a
Page 26 of 32 emi-closed form solutions for barrier options ... temperature gradient at this boundary is supposed to be known.As such, treating it in financial terms, one can immediately recognize this as the pricing problem forthe American option where the exercise boundary is also a free boundary, i.e. is not known. However, theoption Delta ∂u∂z at the boundary z = y ( τ ) is known (it follows from the conditions ∂C A ∂S | S = S B ( t ) = 1 and ∂P A ∂S | S = S B ( t ) = −
1. Also the boundary condition for the American Call and Put at the exercise boundary(the moving boundary) is set as 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 Put. This problem was solved in (Carr and Itkin, 2020) by using the method of generalizedintegral transform for the time-dependent OU model, and in (Lipton and Kaushansky, 2020b) for theBlack-Scholes model with constant coefficients by using the method of heat potentials.On contrary, we want to emphasize that many of the problems considered in this paper have not beensolved yet in physics, and are mentioned in (Kartashov, 2001) as yet unsolved problems. Therefore, theresults obtained in this paper also make a contribution to physics as they can be easily re-formulated interms of the above mentioned physics problems.Another connection of our results with physics is about the first passage time (FPT) problem consid-ered in Section 4.1.3. As mentioned in (Ding and Rangarajan, 2004) (see also references therein), the FPTproblem finds applications in many areas of science and engineering . A sampling of these applicationsincludes (but not limited to) • statistical physics (study of anomalous diffusion) • neuroscience (analysis of neuron firing models) • civil and mechanical engineering (analysis of structural failure) • chemical physics (study of noise assisted potential barrier crossings) • hydrology (optimal design of dams) • imaging (study of image blurring due to hand jitter)In particular, in (Redner, 2001) the author analyses the fundamental connection between the first-passage properties of diffusion and electrostatics. Basic questions of first passage include where a diffusingparticle is absorbed on a boundary and when does this absorption event occur. These are time-integratedattributes, obtained by integration of a time-dependent observable over all time. For example, to deter-mine when a particle is absorbed, we should compute the first-passage probability to the boundary andthen integrate over all time to obtain the eventual hitting probability. However, it is more elegant toreverse the order of calculation and first integrate the equation of motion over time and then compute theoutgoing flux at the boundary. This first step transforms the diffusion equation to the simpler Laplaceequation. Then, in computing the flux, the exit probability is just the electric field at the boundary point.Thus, there is a complete correspondence between a first-passage problem and an electrostatic problemin the same geometry. This mapping is simple yet powerful, and can be adapted to compute relatedtime-integrated properties, such as the splitting probabilities and the moments of the exit time.Other connections to physics problems such as kinetics of spin systems, first passage in composite andfluctuating systems, hydrodynamic transport. reaction-rate theory, etc. could also be found in (Redner,2001).With the hope that we managed to convince the reader about a strong connection between the subjectof this paper and physics, we stop here this excursus into a wonderful world of physics leaving curiousreaders to extend it themselves. In Sections 1, 2 we constructed semi-closed form solutions for the prices of Up-and-Out barrier Call option C uao . Despite same could be done for the Down-and-Out options, alternatively we can use the parityfor barrier options, (Hull, 1997). Then the price of the Down-and-Out barrier Call option C dao can be Page 27 of 32 emi-closed form solutions for barrier options ... found as C dao = C v − C uao , where C v is the price of the European vanilla Call option. For the modelsconsidered in this paper the latter is known in closed form, (Andersen and Piterbarg, 2010). The doublebarrier case was also considered in Section 3.3.Despite in our test we assumed the barrier H to be constant in time, the whole framework is developedfor the genera case where the barrier is some arbitrary function of time.From the computational point of view the proposed solution is very efficient as this is shown inSection 5. 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 z space can be easily increasedby using high order quadratures. On the other hand doing same for the FD method is not easy (i.e., itsignificantly increases the complexity of the method, e.g., see (Itkin, 2017)).Another advantage of the approach advocated in this paper is computation of option Greeks. Sincethe option prices in both the BP and GIT methods are represented in closed form via integrals, theexplicit dependence of prices on the model parameters is available and transparent. Therefore, explicitrepresentations of the option Greeks can be obtained by a simple differentiation under the integrals.This means that the values of Greeks can be calculated simultaneously with the prices almost with noincrease in time. This is because differentiation under the integrals slightly changes the integrands, andthese changes could be represented as changes in weights of the quadrature scheme used to numericallycompute the integrals. Since the major computational time has to be spent for computation of densitieswhich contain special functions, they can be saved during the calculation of the prices, and then reusedfor computation of Greeks.Note, that the FD method also provides the values of Delta, Gamma and Theta on the FD grid, while,for instance, for Vega one need to bump the model volatility and rerun the whole scheme. But for the BPand GIT methods computation of Delta or Vega is done uniformly. Also, the ability of fast computationof Greeks is important for model calibration. Therefore, one can efficiently calibrate the CIR and CEVmodels to the market data by using the BP and GIT methods, since the semi-explicit nature of the finalexpressions allows quasi-analytical formulae for the gradient of the loss function. Acknowledgments
We are grateful to Daniel Duffy and Alex Lipton for some helpful comments. Dmitry Muravey acknowl-edges support by the Russian Science Foundation under the Grant number 20-68-47030.
References
M. Abramowitz and I. Stegun.
Handbook of Mathematical Functions . Dover Publications, Inc., 1964.L. Alil, P. Patie, and J. L. Pedersen. Representations of the first hitting time density of an ornstein-uhlenbeck process.
Stochastic Models , 21(4):967–980, 2005.L. Alili and P. Patie. Boundary-crossing identities for diffusions having the time-inversion property.
Journal of Theoretical Probability , 23(1):65–85, 2010.L.B.G. Andersen and V.V. Piterbarg.
Interest Rate Modeling . Number v. 2 in Interest Rate Modeling.Atlantic Financial Press, 2010. ISBN 9780984422111.
Page 28 of 32 emi-closed form solutions for barrier options ...
H. Bateman and A. Erdélyi.
Higher Transcendental Functions , volume 1 of
Bateman Manuscript ProjectCalifornia Institute of Technology . McGraw-Hill, 1953.P. Carr and A. Itkin. Semi-closed form solutions for barrier and American options written on a time-dependent Ornstein Uhlenbeck process, March 2020. Arxiv:2003.08853.P. Carr and V. Linetsky. A jump to default extended CEV model: an application of Bessel processes.
Finance and Stochastics , 10:303–330, 2006.J. Cox. Notes on option pricing i. constant elasticity of variance diffusions. Technical report, StanfordUniversity working paper, 1975.J.C. Cox, J.E. Ingersoll, and S.R. Ross. A theory of the term structure of interest rates.
Econometrica ,53(2):385–408, 1985.D. Davydov and V. Linetsky. Pricing and hedging path-dependent options under the CEV process.
Management Science , 47(7):949–965, 2001.M. Ding and G. Rangarajan. First Passage Time Problem: A Fokker-Planck Approach. In Luc T. Wille,editor,
New Directions in Statistical Physics: Econophysics, Bioinformatics, and Pattern Recognition ,chapter 3, pages 31–46. Springer-Verlag, Berlin Heidelberg, 2004.D. Emanuel and J. Macbeth. Further results on the constant elasticity of variance call option pricingmodel.
Journal of Financial and Quantitative Analysis , 17:533–554, 1982.A. Friedman.
Partial Differential Equations of Parabolic Type . Prentice-Hall, New Jersey„ 1964.D. Goldenberg. A unified method for pricing options on diffusion-processes.
Journal of Financial Eco-nomics , 29:3–34, 1991.I.S. Gradshtein and I.M. Ryzhik.
Table of Integrals, Series, and Products . Elsevier, 2007.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.A. Itkin and D. Muravey. Semi-closed form prices of barrier options in the Hull-White model, April 2020.Arxiv:2004.09591.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.V. Kaushansky, A. Lipton, and C. Reisinger. Numerical analysis of an extended structural default modelwith mutual liabilities and jump risk.
Journal of Computational Science , 24(218–231), 2018.Gregory F. Lawler. Notes on the bessel process. 2018. URL . Corpus ID: 52200396.V. Linetsky and R. Mendoza. Encyclopedia of quantitative finance. In
Constant Elasticity of Variance(CEV) Diffusion Model . John Wiley & Sons, 2010. ISBN 9780470061602.
Page 29 of 32 emi-closed form solutions for barrier options ...
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 for a reducible diffusion process.
Quan-titative Finance, , 2020a. published online.A. Lipton and V. Kaushansky. On three important problems in mathematical finance.
The Journal ofDerivatives. Special Issue , 2020b.B. Ya. Lyubov.
Kinetic Theory of Phase Transformations . Amerind Publishing, 1978. ASIN:B0000EGG4O.D. Mumford, C. Musiliand M. Nori, E. Previato, and M. Stillman.
Tata Lectures on Theta . Progress inMathematics. Birkhäuser Boston, 1983. ISBN 9780817631093.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.S. Redner.
A guide to first-passage processes . Cambridge University Press, Cambridge, 2001.D. Revuz and M. Yor.
Continuous Martingales and Brownian Motion . Springer, Berlin, Germany, 3rdedition, 1999.A.N. Tikhonov and A.A. Samarskii.
Equations of mathematical physics . Pergamon Press, Oxford, 1963.P. Wilmott, A.L. Lewis, and D.J. Duffy. Modeling volatility and valuing derivatives under anchoring.
Wilmott , (73):48–57, September 2014.
Appendices
A General construction of the potential method.
In this Section we generalize the construction of the potential method originally proposed for the heatequation. For convenience, we follow the notation of (Tikhonov and Samarskii, 1963).Consider a PDE ∂V ( t, x ) ∂t = L ( V ( t, x )) (A.1)where the operator L is a linear differential operator with time-independent coefficients. An example ofsuch the equation is the heat equation and the Bessel equation in Eq. (5). Suppose that the fundamentalsolution (or the Green function) of Eq. (A.1) G ( x, t | ξ, τ ) is known in closed form. Page 30 of 32 emi-closed form solutions for barrier options ...
Suppose we need to solve Eq. (A.1) subject to the homogeneous initial condition V (0 , x ) = 0 . (A.2)Note that this doesn’t restrict our consideration, as by the change of variables in Eq. (54) the problemwith inhomogeneous initial condition can be transformed to the problem with V (0 , x ) = 0.We are interesting in solving the first boundary-value problem assuming that the boundaries are somefunctions of time. Consider, for simplicity just a semi-bounded region with x ≥ y ( t ) , x < ∞ , t >
0, withthe following boundary conditions V ( t, x ) (cid:12)(cid:12)(cid:12) x →∞ = 0 , V ( t, y ( t )) = φ ( t ) . (A.3)The assumption that the problem has just one time-dependent boundary can be easily relaxed as this isdemonstrated in Section 3.3, and is not a restriction of the method.Let us introduce the single layer potential, (Tikhonov and Samarskii, 1963)Π( x, t ) = Z t Ψ( τ ) ∂G∂ξ ( x, t | ξ, τ ) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ, (A.4)where Ψ( t ) the potential density. The single layer potential is continuous and twice differentiable functionin x , and continuous and differentiable in t . Then the below proposition follows Proposition A.1.
The potential function in Eq. (A.4) solves Eq. (A.1). The potential density is deter-mined by the boundary condition in Eq. (A.3) and solves the Volterra equation of second kind φ ( t ) = b Ψ( t ) + Z t Ψ( τ ) ∂G∂ξ ( y ( t ) , t | y ( τ ) , τ ) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ, (A.5) where b is a constant. Proof.
First, it can be checked that substituting the definition in Eq. (A.4) into Eq. (A.1) we get ∂ Π( x, t ) ∂t = Z t Ψ( τ ) ∂∂ξ ∂G∂t (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ + Ψ( t ) ∂G∂y ( τ ) ( x, y ( t ) | t, t ) dτ = Z t Ψ( τ ) ∂∂ξ ∂G∂t (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ, (A.6) L (Π( t, x )) = Z t Ψ( τ ) ∂∂ξ L ( G ) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ. Combining these two expressions yields Z t Ψ( τ ) ∂∂ξ (cid:20) ∂G∂t − L ( G ) (cid:21) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ = 0 . (A.7)The first line follows from the fact that the Green function solves Eq. (A.1), and G ( x, y ( τ ) | t, τ ) (cid:12)(cid:12)(cid:12) τ = t = G ( x, y ( τ ) | t − τ ) (cid:12)(cid:12)(cid:12) τ = t = δ ( x − y ( t ) = 0 as x = y ( t ). The second line is a consequence of time-independenceof the operator coefficients.To summarize what we got: the potential function satisfies the PDE for x ≥ y ( t ), is bounded atinfinity, and has a zero initial value for any choice of Ψ( t ). Thus, it is the solution of Eq. (A.1), i.e. V ( t, x ) = Z t Ψ( τ ) ∂G∂ξ ( x, t | ξ, τ ) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ. (A.8) Page 31 of 32 emi-closed form solutions for barrier options ...
Now using the boundary condition at x = y ( t ) we obtain from Eq. (A.8) φ ( t ) = Z t Ψ( τ ) ∂G∂ξ ( y ( t ) , t | ξ, τ ) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ. (A.9)However, as shown in (Tikhonov and Samarskii, 1963), for x = y ( t ) the RHS is discontinuous, but withthe finite limiting value at x = y ( t ) + 0. The limiting value could be represented as b Ψ( t ) + Z t Ψ( τ ) ∂G∂ξ ( y ( t ) , t | y ( τ ) , τ ) (cid:12)(cid:12)(cid:12) ξ = y ( τ ) dτ. The constant b depends on the particular form of the operator L . In particular, if L is a second orderparabolic operator with the diffusion coefficient a , then b = 1 / (2 a ), (Tikhonov and Samarskii, 1963).Thus, for instance, for a = 1 / b = 1.= 1.