Multilayer heat equations: application to finance
aa r X i v : . [ q -f i n . C P ] F e b Multilayer heat equations: application to finance
Andrey Itkin Alexander Lipton and Dmitry Muravey Tandon School of Engineering, New York University, New York, USA The Jerusalem School of Business Administration, The Hebrew University of Jerusalem, Jerusalem, Israel;Connection Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA, USA Moscow State University, Moscow, Russia
February 17, 2021 I n this paper, we develop a Multilayer (ML) method for solving one-factor parabolicequations. Our approach provides a powerful alternative to the well-known finitedifference and Monte Carlo methods. We discuss various advantages of this approach,which judiciously combines semi-analytical and numerical techniques and provides a fastand accurate way of finding solutions to the corresponding equations. To introducethe core of the method, we consider multilayer heat equations, known in physics for arelatively long time but never used when solving financial problems. Thus, we expandthe analytic machinery of quantitative finance by augmenting it with the ML method.We demonstrate how one can solve various problems of mathematical finance by usingour approach. Specifically, we develop efficient algorithms for pricing barrier optionsfor time-dependent one-factor short-rate models, such as Black-Karasinski and Verhulst.Besides, we show how to solve the well-known Dupire equation quickly and accurately.Numerical examples confirm that our approach is considerably more efficient for solvingthe corresponding partial differential equations than the conventional finite differencemethod by being much faster and more accurate than the known alternatives. Introduction
The problem of solving partial differential equations (PDEs) with moving boundaries appears naturallyin various areas of science and technology. As mentioned in (Kartashov, 2001), such problems have beenknown in physics for a long time. They arise in several fields, such as (a) nuclear power engineering andsafety of nuclear reactors; (b) combustion in solid-propellant rocket engines; (c) laser action on solids;(d) the theory of phase transitions (the Stefan problem and the Verigin problem); (e) the processes ofsublimation in freezing and melting; (f) in the kinetic theory of crystal growth; etc., see (Kartashov,1999) and references therein. Analytical solutions to these problems often require rather sophisticatedmethods., which were actively developed by the Russian mathematical school in the 20th century startingfrom A.V. Luikov, and then by B.Ya. Lyubov, E.M. Kartashov, and many others.1 ultilayer heat equations: application to finance
As applied to mathematical finance, one of these methods - the method of heat potentials (HP) -was actively utilized by A. Lipton and his co-authors to solve various mathematical finance problems,see (Lipton, 2001; Lipton and de Prado, 2020) and references therein. A complementary method of ageneralized integral transform (GIT) is developed in (Carr and Itkin, 2021; Itkin and Muravey, 2020a;Carr et al., 2020) to price barrier and American options in the semi-closed form. These authors studiedthe time-dependent Ornstein-Uhlenbeck (OU), Hull-White, CIR, and CEV models. An extension of themethod of heat potentials for the Bessel process called the method of Bessel potentials is developed by(Carr et al., 2020), who also describe a general scheme of how to construct the potential method for anylinear differential operator with time-independent coefficients. Finally, they also extended the method ofgeneralized integral transform to the Bessel process. In all cases, a semi-analytical (or semi-closed form)solution means that first, one needs to solve a linear Volterra equation of the second kind. Then theoption price is represented as a one-dimensional integral.(Carr and Itkin, 2021; Itkin and Muravey, 2020a; Carr et al., 2020) show that the new method iscomputationally more efficient than the existing ones, such as the backward and forward finite differencemethods while providing better accuracy and stability. Also, the heat potential and GIT methods donot duplicate but rather complement each other. The former provides very accurate results for shortmaturities, and the latter for long maturities.Even though many new problems have been solved in the above-cited papers, some of the financialmodels are hard to solve by using these methods directly. For instance, this is the case for the Black-Karasinski model, popular among practitioners. Another problem is the calibration of the local (orimplied) volatility surface in various one-factor models. Almost all popular analytic and semi-analyticalmethods approach the solution of this problem by doing it term-by-term, which, doubtless, producescomputational errors. For more details, see (Itkin, 2020) and references therein.In this paper, we attack this class of problems (some of them unsolved in the semi-analytical form)by using another method, which we call the method of multilayer (ML) heat equation. An alternativeapproach is given in (Dias), where an innovative technique of recursive images is presented to obtainsolutions to the transient diffusion equation in a N -layered material based on the superposition of Greenfunctions for a semi-infinite material. The solution is initially built for a single layer over a substrateby constructing a sequential sum of reflected image functions. These functions are chosen to satisfy insequence the boundary conditions, first at the front interface, then at the back interface, then again atthe front interface, and so on until the added functions’ magnitude becomes negligible.Based on this so-called "1-layer" algorithm, the author also constructs a "2-layer" algorithm by sequen-tial application of the "1-layer" algorithm first to layer 1, then to layer 2, then again to layer 1, and so on.The sequential application of the N − N -layer algorithm. This schemeworks for the first and second kind boundary conditions but does not apply to the case where there is acontact resistance between layers or the convective heat transfer at the end interfaces.Note that this algorithm as applied to the local volatility calibration problem is similar to the approachused in (Lipton and Sepp, 2011; Itkin and Lipton, 2018; Carr and Itkin, 2020, 2019).Since the ML method splits the whole (possibly infinite) domain in the space variable into a sequenceof sub-domains, one could extend it naturally to solving parabolic equations with coefficients beingfunctions of time t and location x . At every sub-interval, the corresponding parabolic operator could beeither approximated by the operator with the space-homogeneous coefficients or, possibly, reduced to theheat equation by a series of transformations. After either approximation or reduction, the ML methodcan be applied.Moreover, the method could be extended further to deal with non-linear volatility, drift, and killingterm. Again, piecewise approximations of these terms lead to the parabolic equations at every sub-interval that could be transformed to the heat equation. Then, the application of the ML method solvesthe problem. Page 2 of 36 ultilayer heat equations: application to finance
The main idea of this paper is to combine the ML method with the method of heat potentials andthe GIT method. Since both provide a semi-analytical solution for sub-interval problems, a combinationof these solutions within the ML heat equation method results in the problem’s full solution, expressedexplicitly via one-dimensional integrals. At each layer, these integrals depend either on the yet unknownpotential density (in the HP method) or on the solution gradient at the layer’s boundaries (the GITmethod). These unknown functions solve the interconnected systems of the integral Volterra equationsof the second kind derived in the paper. Once this solution is found (either numerically or, sometimes,analytically), the whole problem is solved. Note that one can transform the system of integral equationsto linear equations on a time-space grid, which is lower banded (in our case, block lower triangular).Therefore, the corresponding system can be solved with complexity O ( M N ) where N is the number oflayers, and M is the number of time steps, see (Itkin and Muravey, 2020a) in more detail.We also propose a particular construction of the layers’ internal boundaries, which allows the repre-sentation of every integral in the Volterra equation as convolution. Applying the Laplace transform, weobtain a system of linear equations with a block-tridiagonal matrix (it contains four blocks). This systemcan be efficiently solved numerically (with complexity O ( N )). In some cases, it can be solved analytically.After this system’s solution is found, we use the Gaver-Stehfest method to compute the inverse Laplacetransform, also with linear complexity in the number of layers N . This algorithm solves the system ofthe Volterra equations and thus the whole problem.We illustrate these novel ideas by representing several significant financial problems in the formsuitable for solving them by the ML method. These problems include pricing barrier options in the time-dependent Black-Karasinski and modified Black-Karasinski (Verhulst) models, see (Itkin et al., 2020), aswell as the solution of the Dupire equation. We also provide several numerical examples to demonstrateour method’s high speed and accuracy compared with standard finite-difference (FD) methods.To the best of our knowledge, all the paper results are new and contribute to the existing financialand physics literature. It is interesting to note that our method is capable of solving similar problemsthat appear in medicine and biology in addition to finance. For instance, our technique is well-suitedfor studying (a) the growth of diffusive brain tumors, which considers the brain tissue’s heterogeneity,(Asvestas et al., 2014); (b) the transdermal drug release from an iontophoretic system, (Pontrelli et al.,2016); (c) and many other similar problems. It is imperative to emphasize that our method allows solvingthe ML problems with time-dependent boundaries and time- and space-dependent diffusion coefficients.In contrast, the method of (Carr and March, 2018) and all other known approaches operate only withconstant boundaries (possibly with time-dependent boundary conditions) and spatially piecewise constantdiffusion coefficients. Moreover, their setting corresponds to one of our numerical examples in Section 5.Since in (Carr and March, 2018) the solution is obtained by using spectral (eigenvector) series, while weapply the Laplace transform method, our approach is about 1000 times faster.The rest of the paper is organized as follows. In Section 1 we construct the solution of the ML heatequation by using the method of heat potentials. In Section 2, we solve this equation by using the GITmethod. In Section 3.1.1 we describe the pricing of barrier options in the time-dependent Black-Karasinski(BK) model and also in our modification of this model, which was introduced in (Itkin et al., 2020) and iscalled the Verhulst model. In particular, we demonstrate how to reduce the pricing PDEs for both modelsto the ML heat equation. Also, in Section 3.1.1 we provide a generalization of this approach for someother models. In Section 3.2 we apply the results of Section 2.3 to investigate the case of space-dependentvolatility σ ( x ) in conjunction with solving the Dupire equation. Section 4 is dedicated to the solutionof the Volterra equations. In particular, we describe a construction of the internal boundaries, whichallows a transformation of the Volterra equations of the second kind to Abel equations. We solve thelatter equations via the Laplace transform. Section 5 describes some numerical experiments with the MLmethod. The final section concludes. More general potential methods, e.g., the method of Bessel potentials, can also be used in such a scheme. Page 3 of 36 ultilayer heat equations: application to finance
Let us consider the following initial-boundary problem L u ( τ, x ) = 0 , ( x, τ ) ∈ Ω : h y − ( τ ) , y + ( τ ) i × R + , (1) u (0 , x ) = f ( x ) , y − (0) < x < y + (0) ,u ( τ, y − ( τ )) = χ − ( τ ) , u ( τ, y + ( τ )) = χ + ( τ ) . Here the operator L is a partial differential operator of the parabolic type L = − ∂∂τ + ∂∂x (cid:18) σ ( τ, x ) ∂∂x (cid:19) + µ ( τ, x ) ∂∂x + ν ( τ, x ) , (2) σ ( τ, x ) , µ ( τ, x ) , ν ( τ, x ) are some known functions, Ω is the spatial-temporal domain with curvilinear tem-poral boundaries, and χ − ( τ ) , χ + ( τ ) are known functions of time (the boundary conditions).Similar to (Itkin and Muravey, 2020a), we represent the solution in the form u ( x, τ ) = q ( x, τ ) + Z y + (0) y − (0) f ( ξ ) G ( x, ξ, τ ) dξ, (3)where G ( x, ξ, τ ) is Green’s function of the problem. Then the function q ( x, τ ) solves a problem similar toEq. (1) but with the homogeneous initial condition L q ( τ, x ) = 0 , ( x, τ ) ∈ Ω : h y − ( τ ) , y + ( τ ) i × R + , (4) q (0 , x ) = 0 , y − (0) < x < y + (0) ,q ( τ, y − ( τ )) = χ − ( τ ) − Z y + (0) y − (0) f ( ξ ) G ( y − ( τ ) , ξ, τ ) dξ = φ − ( τ ) ,q ( τ, y + ( τ )) = χ + ( τ ) − Z y + (0) y − (0) f ( ξ ) G ( y + ( τ ) , ξ, τ ) dξ = φ + ( τ ) . If the Green function G ( x, ξ, τ ) is known, the problem in Eq. (4) can be solved via the HP method,(Itkin and Muravey, 2020b). Otherwise, one can apply the ML method as this is described below.To use the ML method, suppose the domain Ω could be split into N layers: Ω = S Ni =1 Ω i , where eachlayer is a curvilinear stripΩ i = [ y i ( τ ) , y i +1 ( τ )] × R + , y i ( τ ) < y i +1 ( τ ) , ∀ τ > , ∀ i = 1 , . . . , N, (5) y ( τ ) = y − ( τ ) , y N +1 ( τ ) = y + ( τ ) . Let us seek for the solution of the problem Eq. (4) in the form u ( τ, x ) = N X i =1 u i ( τ, x ) x − y i ( τ ) y i +1 ( τ ) − x , x = ( , x ≥ , x < , (6)and request that both u ( τ, x ) and its flux are continuous functions of x . Using these conditions atevery boundary y i ( τ ) , i = 2 , . . . , N together with the boundary conditions yields the following system ofequations u i ( τ, y i +1 ( τ )) = u i +1 ( τ, y i +1 ( τ )) , (7) These conditions are natural in physics if by u ( t, x ) we assume, e.g., the temperature and interpret σ ∂ x u ( t, x ) as theheat flux. Therefore, it is standard to require continuity of the heat flux rather than the first derivative ∂ x u ( t, x ), (LienhardIV and Lienhard V, 2019) Page 4 of 36 ultilayer heat equations: application to finance σ i ( τ, y i +1 ( τ )) ∂u i ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i +1 ( τ ) = σ i +1 ( τ, y i +1 ( τ )) ∂u i +1 ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i +1 ( τ ) , i = 1 , . . . , N − ,u ( τ, y − ( τ )) = χ − ( τ ) , u N +1 ( τ, y + ( τ )) = χ + ( τ ) . The first condition means a continuity of the function u at every boundary y i , i = 1 , . . . , N −
1. Thesecond condition is a continuity of the heat flux at the same boundary. The last line follows from theboundary conditions in Eq. (4).Also, let us define the operator L for the whole domain Ω as follows L = N X i =1 L i x − y i ( τ ) y i +1 ( τ ) − x , (8)where L i = − ∂∂τ + ∂∂x (cid:18) σ i ( τ, x ) ∂∂x (cid:19) + µ i ( τ, x ) ∂∂x + ν i ( τ, x ) , (9)and L i u i = 0.The idea of the ML method is to assume that Green’s function G i ( x, τ | ξ, s ) associated with theoperator L i can be obtained in the closed form. For an arbitrary dependencies σ i ( τ, x ) , µ i ( τ, x ) , ν i ( τ, x )this is not the case, but for various specific forms of these functions this can be done. For instance,when µ ( τ, x ) = ν ( τ, x ) = 0 and σ ( τ, x ) = σ ( τ ) or σ ( τ, x ) = σ ( x ), etc., (Polyanin, 2002). Otherwise,the functions σ i ( τ, x ) , µ i ( τ, x ) , ν i ( τ, x ) can be approximated at every layer, e.g., by piecewise constant orlinear function in x and an arbitrary function of τ , or by piecewise constant functions in τ and piecewiselinear functions in x , etc. These approximations make the ML method somewhat similar to the FDmethod, however, with some critical distinctions, see Section 6.It is important to mention, that the operator L i in Eq. (9), while natural for physics where a divergentform of the parabolic equation (e.g., the heat equation) is commonly accepted, is just rarely used inmathematical finance. Instead, in finance it is natural to consider a non-divergent (non-conservative)form, which for the heat equation reads L i = − ∂∂τ + σ i ( τ, x ) ∂ ∂x . (10)Obviously, when σ i = σ i ( τ ) , ∀ i , i.e. σ i ( τ, x ) is a straight line at given τ , both operators in Eq. (9) andEq. (10) coincide. However, if one solves Eq. (10) by the ML method, it can be unclear what continuitycondition should be used. As shown in (Lejay, 2006), for the divergent heat/diffusion equation withdrift this condition remains the same, i.e. this is a continuity of flux over the boundary. Obviously, anon-divergent heat equation can be represented in this form, i.e. the divergent diffusion part plus drift.Therefore, the continuity condition is still represented by the equality of fluxes over the boundary, butthe equation now includes an extra drift term.As applied to the ML method, this can be seen as follows. Suppose we apply the ML method to someparabolic equation, and approximate all coefficients in the drift and killing terms by piecewise constantfunction at every interval. Then, at the i -th interval this equation reads ∂u i ( τ, x ) ∂τ = σ ( τ, x ) ∂ u i ( τ, x ) ∂x + α i ∂u i ( τ, x ) ∂x + β i u i ( τ, x ) , (11)where α i = const, β i = const, i = 1 , N . By transforming it to a divergent form we obtain ∂u i ( τ, x ) ∂t = ∂∂x (cid:18) σ ( τ, x ) ∂u i ( τ, x ) ∂x (cid:19) + [ α i − σ ( τ, x ) σ x ( τ, x )] ∂u i ( τ, x ) ∂x + β i u i ( τ, x ) . (12) Page 5 of 36 ultilayer heat equations: application to finance
Hence, again, the continuity condition for this equation is given in Eq. (7). Further, Eq. (12) by aseries of transformations can be reduced to a non-divergent heat equation in Eq. (10). Accordingly, thesetransformations should be applied to Eq. (7) as well to obtain the correct continuity conditions.When the external boundaries are constant, i.e. y − ( t ) = χ − ( t ) = const, y + ( t ) = χ + ( t ) = const one may use an alternative where a non-divergent heat equation can be reduced to a divergent one by achange of variables x y = g ( x ), where g ( x ) is some function which depends on σ ( x ). In more detailthis is shown in Appendix A. Accordingly, the operator Eq. (10) transforms to L i = − ∂∂τ + ∂∂y (cid:18) Ξ i ( τ, y ) ∂∂y (cid:19) , (13)where Ξ ( τ, y ) is a new diffusion coefficient, which can be expressed via σ ( τ, x ), again see Appendix A.In what follows, we provide our analysis for Eq. (9); Eq. (10) can be analyzed similarly, as explainedabove. For simplicity and without loss of generality, we give an exposition of the HP method assuming µ i ( τ, x ) = ν i ( τ, x ) = 0, and either σ i ( τ, x ) = σ i ( τ ), or σ i ( τ, x ) is a piecewise constant function of x forevery layer. In this case each equation L i u i = 0 by some change of variables τ ¯ τ , x ¯ x can betransformed to the heat equation Eq. (63) with σ i (¯ τ , ¯ x ) = 1 , (Polyanin, 2002), and the correspondingGreen function G (¯ x, ξ, ¯ τ ) reads G (¯ x, ξ, ¯ τ ) = 12 √ π ¯ τ e − (¯ x − ξ )24¯ τ . (14)Also, these transformations modify the boundary y ( τ ) ¯ y (¯ τ ). Some examples of such transforma-tions are presented in Section 3. In Appendix B we also provide some recipes on how to proceed if oneneeds to generalize this approach by considering a general case σ = σ ( τ, x ).To the end of this section, for easiness of reading let us drop the bar over new variables. Now, followingthe general idea of the method of heat potentials for pricing double barrier options, (Itkin and Muravey,2020a; Carr et al., 2020), we represent each function q i ( τ, x ) as q i ( τ i , x ) = Z τ i ( Ψ i ( k ) ∂G ( x, ξ, τ i − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i ( k ) + Φ i ( k ) ∂G ( x, ξ, τ i − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i +1 ( k ) ) dk. (15)In Eq. (15) the second integral is a sum of two single layer potentials with the potential densities Ψ i ( τ ) andΦ i ( τ ). By writing Eq. (15) we take into account that according to Eq. (45) and Eq. (72), the transformedtime τ might differ for each interval, therefore, the notation τ i is used. However, e.g., for the problemdescribed in Section 3.1.1 all new times are equal, i.e. τ i = τ, i = 1 , . . . , N + 1.Since the domain Ω consists of N layers, there are 2 N unknown density functions Ψ i ( τ i ) , Φ i ( τ i ) , i =1 , . . . , N . To determine them one need to plug the representation of q i ( τ i , x ) in Eq. (15) into Eq. (7), andthen solve thus obtained system of the integral Volterra equations of the second kind.However, it is known, (Tikhonov and Samarskii, 1963) that the integral in Eq. (15) for x = y i ( τ i ) and x = y i +1 ( τ i ) is discontinuous, but with the finite value at x = y i ( τ i ) + ε, ∀ i = 1 , . . . , N + 1 when ε → q i ( τ, y i ( τ )) = 12 σ i ( y i ( τ )) Ψ i ( τ ) (16)+ Z τ ( Ψ i ( k ) ∂G ( y i ( τ ) , ξ, τ − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i ( k ) + Φ i ( k ) ∂G ( y i ( τ ) , ξ, τ − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i +1 ( k ) ) dk, τ = τ i , Of course, there exist other possible representations of σ i which give rise to the heat equation, or e.g., to the Besselequation, Carr et al. (2020). Page 6 of 36 ultilayer heat equations: application to finance q i ( τ, y i +1 ( τ )) = − σ i ( y i +1 ( τ )) Φ i ( τ )+ Z τ ( Ψ i ( k ) ∂G ( y i +1 ( τ ) , ξ, τ − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i ( k ) + Φ i ( k ) ∂G ( y i +1 ( τ ) , ξ, τ − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i +1 ( k ) ) dk. τ = τ i +1 . The gradients of q i ( τ, x ) for the heat equation with σ i = σ = const have been derived first in (Liptonand Kaushansky, 2020a,b), and later in (Itkin and Muravey, 2020b) by using a different method . Theresult read ∂q i ( τ, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i ( τ ) = − Ψ i ( τ )2 σ (cid:18) √ πτ + y ′ i ( τ ) σ (cid:19) + Z τ Ψ i ( k ) e − ( yi ( τ ) − yi ( k ))24 σ τ − k ) − Ψ i ( τ )4 σ p π ( τ − k ) dk (17) − Z τ Ψ i ( k ) ( y i ( τ ) − y i ( k )) e − ( yi ( τ ) − yi ( k ))24 σ τ − k ) σ p π ( τ − k ) dk − Z τ Φ i ( k ) ∂ G ( x, ξ, σ ( τ − k )) ∂ξ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i +1 ( k ) x = y i ( τ ) dk, τ = τ i ,∂q i ( τ, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i +1 ( τ ) = − Z τ Ψ i ( k ) ∂ G ( x, ξ, σ ( τ − k )) ∂ξ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = y i ( k ) x = y i +1 ( τ ) dk − Φ i ( τ )2 σ √ πτ − y ′ i +1 ( τ ) σ ! + Z τ Φ i ( k ) e − ( yi +1( τ ) − yi +1( k ))24 σ τ − k ) − Φ i ( τ )4 σ p π ( τ − k ) dk − Z τ Φ i ( k ) ( y i +1 ( τ ) − y i +1 ( k )) e − ( yi +1( τ ) − yi +1( k ))24 σ τ − k ) σ p π ( τ − k ) dk, τ = τ i +1 , where G ( x, ξ, τ ) is given in Eq. (14). An alternative method to construct the ML problem solution is generalized integral transform (GIT). TheGIT method is used in physics, (Kartashov, 1999, 2001), but was unknown in finance until its first use in(Carr and Itkin, 2021). The previously known solution to the heat equation, using the GIT method, wasobtained only for the domain S ∈ [0 , y ( t )]. For other domains, the solution was unknown even in physics.(Itkin and Muravey, 2020a) were the first to construct the GIT solution for the domain S ∈ [ y ( t ) , ∞ ).The latter technique was extended further for the CIR and CEV models, (Carr et al., 2020), the Black-Karasinski model, (Itkin et al., 2020), and finally for double barrier options in (Itkin and Muravey, 2020b).The latter problem deals with the spatial domain determined by two moving in time boundaries, andboundary conditions, which are arbitrary functions of time.The GIT and HP methods are similar but have an essential difference. In the HP method, the solutionis represented in the form of heat potential with the unknown potential density function Ψ( τ ) which solvesthe corresponding Volterra equation of the second kind, see Section 1. In the GIT method, similar to,e.g., the Fourier method, we start with applying some integral transform to the PDE under consideration.The transform has to be such that the transformed equation with x p is solvable analytically in time.The second step is to construct an inverse transform, which could be computed analytically using thecomplex analysis. If this is possible, then the solution can be represented as an explicit integral of somekernel multiplied by the unknown function Ω( τ ). Hence, this looks pretty similar to the HP method.However, the function Ω( τ ) has now a transparent meaning - this is the gradient of the solution at themoving boundary. It turns out that this gradient also solves a Volterra equation of the second kind. As These results can be naturally generalized for the case σ = σ ( x ). Page 7 of 36 ultilayer heat equations: application to finance mentioned, explicit construction of such forward and inverse transforms is performed in (Carr and Itkin,2021; Itkin and Muravey, 2020a; Carr et al., 2020; Itkin et al., 2020; Itkin and Muravey, 2020b) for variousmodels and spatial domains. Also, the authors show that the performance of both methods is the same.Both HP and GIT methods are faster than the finite-difference approach and provide higher accuracy.As mentioned in (Itkin and Muravey, 2020b), it is not unreasonable to ask why we need two methods -the HP and GIT, which are used to solve the same problem and demonstrate the same performance. Theanswer is interesting. As shown in (Carr et al., 2020), the GIT method produces very accurate results athigh strikes and maturities (i.e., when the option price is relatively small), in contrast to the HP method,which struggles under these circumstances. One can verify this fact by looking at the exponents underthe GIT solution integral proportional to the time τ . Contrary, when the price is higher (short maturities,low strikes), the GIT method is slightly less accurate than the HP method since the exponents in the HPsolution integral are inversely proportional to τ .Thus, the GIT and HP methods complement each other for the CIR and CEV models. For othermodels reducible to the heat equation, the same conclusion holds; see (Itkin and Muravey, 2020a). Thisstatement is true because GIT integrals contain the difference of two exponents, which becomes small atlarge τ . On the contrary, the HP exponent tends to one at large τ . Therefore, the convergence propertiesof the two methods are different at large τ , so they complement each other.This situation is well known for the heat equation with constant coefficients, (Lipton, 2001). Thereexist two representations of the solution: one - obtained by using the method of images, and the otherone - by the Fourier series. Although both solutions are equal when considered as infinite series, theirconvergence properties are different. To apply the GIT method to the solution of Eq. (1), we can use the results obtained in (Itkin and Muravey,2020a). There it is assumed that L i is the heat equation operator L i = − ∂∂τ + ∂ ∂x , (18)and then the solution of Eq. (1) can be represented in the form u i ( τ, x ) = ∞ X n = −∞ ( Z y i +1 (0) y i (0) u i (0 , ξ )Υ n,i ( x, τ | ξ, dξ + Z τ h Ω i ( s ) + χ + i ( s ) y ′ i +1 ( s ) i Υ n,i ( x, τ | y i +1 ( s ) , s ) ds, + Z τ h Θ i ( s ) − χ − i ( s ) y ′ i ( s ) i Υ n,i ( x, τ | y i ( s ) , s ) ds (19)+ Z τ χ − i ( s )Λ n,i ( x, τ | y i ( s ) , s ) − χ + i ( s )Λ n,i ( x, τ | y i +1 ( s ) , s ) ds ) , Υ n,i ( x, τ | ξ, s ) = 12 p π ( τ − s ) " e − (2 nli ( τ )+ x − ξ )24( τ − s ) − e − (2 nli ( τ )+ x + ξ − yi ( τ ))24( τ − s ) , Λ n,i ( x, τ | ξ, s ) = x − ξ + 2 nl i ( τ )4 p π ( τ − s ) e − (2 nli ( τ )+ x − ξ )24( τ − s ) + x + ξ − y i ( τ ) + 2 nl i ( τ )4 p π ( τ − s ) e − (2 nli ( τ )+ x + ξ − yi ( τ ))24( τ − s ) . Here χ − i ( τ ) , χ + i ( τ ) are the boundary conditions at the left and right boundaries of the i -th interval,and l i ( τ ) = y i +1 ( τ ) − y i ( τ ) , τ = τ i , (20) Page 8 of 36 ultilayer heat equations: application to finance Ω i ( τ ) = − ∂u i ( τ, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i ( τ ) Θ i ( τ ) = ∂u i ( τ, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i +1 ( τ ) . The functions Ω( τ ) , Θ( τ ) for the heat equation in Eq. (18) can be found by solving a system of theVolterra equations of the second kind. As applied to our problem for the i -th interval with i ∈ [1 , N ], itreads, (Itkin and Muravey, 2020b) − Ω i ( τ ) = Z y i +1 (0) y i (0) u (0 , ξ ) υ − i ( τ | ξ, dξ (21) − χ − i ( τ ) √ πτ + Z τ χ − i ( s ) − χ − i ( τ )2 p π ( τ − s ) ds + Z τ h χ − i ( s ) d (cid:16) η − i ( τ | y i ( s ) , s ) (cid:17) − χ + i ( s ) d (cid:16) η − i ( τ | y i +1 ( s ) , s ) (cid:17)i − Z τ Ω i ( s ) y i ( τ ) − y i ( s )2 p π ( τ − s ) e − ( yi ( τ ) − yi ( s ))24( τ − s ) ds + Z τ h Θ i ( s ) υ − i ( τ | y i +1 ( s ) , s ) + Ω i ( s ) υ − ,i ( τ | s ) i ds Θ i ( τ ) = Z y i +1 (0) y i (0) u (0 , ξ ) υ + i ( τ | ξ, dξ + χ + i ( τ ) √ πτ − Z τ χ + i ( s ) − χ + i ( τ )2 p π ( τ − s ) ds + Z τ h χ − i ( s ) d (cid:16) η + i ( τ | y i ( s ) , s ) (cid:17) − χ + i ( s ) d (cid:16) η + i ( τ | y i +1 ( s ) , s ) (cid:17)i − Z τ Θ i ( s ) y i +1 ( τ ) − y i +1 ( s )2 p π ( τ − s ) e − ( yi +1( τ ) − yi +1( s ))24( τ − s ) ds + Z τ h Θ i ( s ) υ +0 ,i ( τ | s ) + Ω i ( s ) υ + i ( τ | y i ( s ) , s ) i ds. Here the following notation is used η − i ( τ | ξ, s ) = − δ ξ,y i ( s ) p π ( τ − s ) + 1 p π ( τ − s ) ∞ X n = −∞ e − ( yi ( τ ) − ξ +2 nli ( τ ))24( τ − s ) , (22) η + i ( τ | ξ, s ) = − δ ξ,y i +1 ( s ) p π ( τ − s ) + 1 p π ( τ − s ) ∞ X n = −∞ e − ( yi ( τ ) − ξ +(2 n +1) li ( τ ))24( τ − s ) ,υ − i ( τ | ξ, s ) = − ∞ X n = −∞ y i ( τ ) − ξ + 2 nl i ( τ )2 p π ( τ − s ) e − ( yi ( τ ) − ξ +2 nli ( τ ))24( τ − s ) ,υ + i ( τ | ξ, s ) = − ∞ X n = −∞ y i ( τ ) − ξ + (2 n + 1) l i ( τ )2 p π ( τ − s ) e − ( yi ( τ ) − ξ +(2 n +1) li ( τ ))24( τ − s ) ,υ − ,i ( τ | s ) = υ − i ( τ | y i ( s ) , s ) + y i ( τ ) − y i ( s )2 p π ( τ − s ) e − ( yi ( τ ) − yi ( s ))24( τ − s ) ,υ +0 ,i ( τ | s ) = υ + i ( τ | y i +1 ( s ) , s ) + y i +1 ( τ ) − y i +1 ( s )2 p π ( τ − s ) e − ( yi +1( τ ) − yi +1( s ))24( τ − s ) , where δ ξ,x is the Kronecker symbol. It is worth emphasizing that all summands in Eq. (21) are regular.The integrands in the first two lines have weak (integrable) singularities, while other summands areregular.At the boundaries of the domain where the solution of our problem is defined, we have χ − = χ − , χ − N = χ + , (23)where N is the number of intervals.In (Itkin and Muravey, 2020b), an alternative system of the Volterra equations of the second kind isalso obtained, which has the form of Eq. (21), but with a different definition of the coefficients. We have η − i ( τ | ξ, s ) = − δ ξ,y i ( s ) p π ( τ − s ) + 1 l i ( τ ) θ ( φ i ( ξ, y i ( τ )) , ω i ) , (24) Page 9 of 36 ultilayer heat equations: application to finance η + i ( τ | ξ, s ) = − δ ξ,y i +1 ( s ) p π ( τ − s ) + 1 l ( τ ) θ ( φ i ( ξ + l i ( τ ) , y i ( τ )) , ω i ) ,υ − i ( τ | ξ, s ) = − π l i ( τ ) θ ′ ( φ i ( ξ, y i ( τ )) , ω i ) ,υ + i ( τ | ξ, s ) = − π l i ( τ ) θ ′ ( φ i ( ξ + l i ( τ ) , y i ( τ )) , ω i ) . Here θ ( z, ω ) is the Jacobi theta function of the third kind, (Mumford et al., 1983), which is defined asfollows: θ ( z, ω ) = 1 + 2 ∞ X n =1 ω n cos (2 nz ) , (25)Also, in Eq. (24) the following notation is used ω i = e − π τ − s ) l i ( τ ) , φ i ( x, ξ ) = π ( x − ξ )2 l i ( τ ) , (26) ∂θ ( z, ω ) ∂z = θ ′ ( z, ω ) = − ∞ X n =1 nω n sin (2 nz ) . Formulas Eq. (24) and Eq. (22) are complementary. Since the exponents in the definition of the thetafunctions in Eq. (24) are proportional to the difference τ − s , the Fourier series Eq. (24) converge fast if τ − s is large. Contrary, the exponents in Eq. (22) are inversely proportional to τ − s . Therefore, theseries Eq. (22) converge fast if τ − s is small. σ is piecewise constant Here we assume that σ i ( x ) = σ i , i = 1 , . . . , N , i.e. the volatility is a piecewise constant function of x .For instance, this is true for the problem described in Section 3.2. As shown there, the pricing PDE canbe transformed to Eq. (9) instead of Eq. (71). According to the transformation in Eq. (72) the clock willrun differently at each ML interval, which is inconvenient. Therefore, instead of a change of temporalvariable, below we use a transformation of the spatial variable x . This transformation allows using thesame time at each ML interval. To achieve this, we change the definition of τ ( t ) in Eq. (72) to τ = 12 Z T e − R s [ r ( s ) − q ( s )] dk ds, (27)where r ( t ) , q ( t ) for the problem considered in Section 3.2 are the deterministic interest rate and continuousdividends. This converts the problem in Section 3.2 and the PDE Eq. (70) to ∂U ( τ, x ) ∂τ = σ ( x ) ∂ U ( τ, x ) ∂x , (28) U (0 , x ) = U ( x ) = ( x − S ) + ,U ( τ,
0) = 0 , U ( τ, x ) (cid:12)(cid:12)(cid:12) x →∞ = x − e − R T ( r ( s ) − q ( s )) ds S. And, according to Eq. (69) σ ( T, K ) = v i ( T ) , K ∈ [ K i , K i +1 ] . (29)The boundary and initial conditions in Eq. (28) are the direct translation of those conditions in Eq. (66)and Eq. (67). Page 10 of 36 ultilayer heat equations: application to finance
Again, as shown in Appendix A, the problem in Eq. (28) can be transformed to ∂U ( τ, ˆ x ) ∂τ = ∂∂ ˆ x (cid:18) Ξ (ˆ x ) ∂U ( τ, ˆ x ) ∂ ˆ x (cid:19) , ˆ x = ˆ x ( x ) , (30) U (0 , ˆ x ) = U (ˆ x ) = ( x (ˆ x ) − S ) + ,U ( τ, ˆ x (0)) = 0 , U ( τ, x (ˆ x )) (cid:12)(cid:12)(cid:12) x (ˆ x ) →∞ = x (ˆ x ) − e − R T ( r ( s ) − q ( s )) ds S. Also, since based on Eq. (72) U ( τ, x (ˆ x )) = P ( T, K ) e − R T q ( s ) ds , the continuity of the solution and its flux at all internal boundaries can be expressed as χ + i ( τ ) = χ − i +1 ( τ ) . (31)Ξ i +1 Ω i +1 ( τ ) = − Ξ i Θ i ( τ ) . To use the results of the previous Section, we proceed by applying the following transformation toEq. (19) ¯ x = Ξ i ˆ x, ¯ y ( τ ) = Ξ i y ( τ ) , ¯ ξ = Ξ i ξ, ¯ l ( τ ) ≡ Ξ i l ( τ ) . (32)Note, that the last equality in Eq. (32) is actually the definition of ¯ l ( τ ). Another complication whichcomes due to this transformation is that in new variables ¯ x the layers stop to be continuous. In otherwords, the upper boundary of the i -th layer ¯ y + i ( τ ) and the lower boundary of the ( i + 1)-th layer ¯ y − i +1 ( τ )are now not equal. Therefore, in what follows to avoid any confusion we will explicitly use this notation,i.e. the left and right boundaries of the i -th layer are denoted as y − i ( τ ) and y − i +1 ( τ ).Also, per these transformations, we have u i (0 , ¯ ξ ) = u i (0 , Ξ i , ξ ) , u i (0 , Ξ i , ξ ) dξ = 1Ξ i u i (0 , ¯ ξ ) d ¯ ξ, (33) Z y i +1 (0) y i (0) u (0 , ξ ) υ − i ( τ | ξ, dξ = 1Ξ i Z ¯ y i +1 (0)¯ y i (0) u (0 , ¯ ξ )¯ υ − i ( τ | ¯ ξ, d ¯ ξ, ¯ υ − i ( τ | ¯ ξ,
0) = − ∞ X n = −∞ ¯ y − i ( τ ) − ¯ ξ + 2 n ¯ l i ( τ )2Ξ i p π ( τ − s ) e − (¯ y − i ( τ ) − ¯ ξ +2 n ¯ li ( τ ))24Ξ2 i ( τ − s ) . It is easy to check that this transformation leaves φ i (ˆ x, ξ )) and ∂ ˆ x ( φ i (ˆ x, ξ )) invariant, but with thenew definition ¯ ω i = e − π i ( τ − s )¯ l i ( τ ) . Finally, let us redefine the partial derivatives¯Ω i = − ∂u ( τ, ¯ x ) ∂ ¯ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ¯ x =¯ y − i ( τ ) , ¯Θ i = ∂u ( τ, ¯ x ) ∂ ¯ x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ¯ x =¯ y + i ( τ ) (34)instead of their definitions in Eq. (20), i.e. Ω i = Ξ i ¯Ω i , Θ i = Ξ i ¯Θ i . Then the continuity conditions inEq. (83) change to ¯ χ + i ( τ ) = ¯ χ − i +1 ( τ ) , (35)Ξ i +1 ¯Ω i +1 ( τ ) = − Ξ i ¯Θ i ( τ ) . Page 11 of 36 ultilayer heat equations: application to finance
To simplify notation, we omit bars from all new variables assuming this doesn’t bring any confusion.Then Eq. (19) transforms to u i ( τ, x ) = ∞ X n = −∞ ( i Z y i +1 (0) y i (0) u i (0 , ξ )Υ n,i ( x, τ | ξ, dξ (36)+ Z τ (cid:20) i Ω i ( s ) + 1Ξ i +1 χ + i ( s ) y + ′ i ( s ) (cid:21) Υ n,i ( x, τ | y + i ( s ) , s ) ds, + 1Ξ i Z τ h Θ i ( s ) − χ − i ( s ) y − ′ i ( s ) i Υ n,i ( x, τ | y − i ( s ) , s ) ds + Z τ χ − i ( s )Λ n,i ( x, τ | y − i ( s ) , s ) − χ + i ( s )Λ n,i ( x, τ | y + i ( s ) , s ) ds ) , Υ n,i ( x, τ | ξ, s ) = 12 p π ( τ − s ) e − (2 nli ( τ )+ x − ξ )24Ξ2 i ( τ − s ) − e − (2 nli ( τ )+ x + ξ − y − i ( τ ))24Ξ2 i ( τ − s ) , Λ n,i ( x, τ | ξ, s ) = x − ξ + 2 nl i ( τ )4Ξ i p π ( τ − s ) e − (2 nli ( τ )+ x − ξ )24Ξ2 i ( τ − s ) + x + ξ − y − i ( τ ) + 2 nl i ( τ )4Ξ i p π ( τ − s ) e − (2 nli ( τ )+ x + ξ − y − i ( τ ))24Ξ2 i ( τ − s ) . By analogy, the modified Volterra equations can be obtained from Eq. (21).In Eq. (21) the unknown variables are [ χ − ( τ ) , χ +1 ( τ ) , Ω ( τ ) , Θ ( τ ) , . . . , χ − N ( τ ) , χ + N ( τ ) , Ω N ( τ ) , Θ N ( τ )],so that there are 4 N unknowns in total. The boundary conditions in Eq. (23) and the continuity conditionsin Eq. (83) reduce the number of unknown variables to 2 N −
2, because χ + i ( τ ) , Θ i ( τ ) can be expressed via χ − i ( τ ) , Ω i ( τ ) and substituted into Eq. (21). Thus, the GIT method provides a significant simplificationof the system of Volterra equations as compared with the HP method. In this section, we consider several models that are frequently used in mathematical finance. We providea short description of each model and demonstrate how to reduce the corresponding pricing problem tothe form suitable for solving it by the ML method.
As the first example, we consider one-factor short interest rate (IR) models. Although these modelswere developed a long time ago, they are still essential and widely used by practitioners. While onecan price zero-coupon bonds (ZCB) and European options on the ZCB and swaptions for many of themanalytically, this is not true for exotic options. For instance, pricing of barrier options when the barriersare time-dependent and could pay time-dependent rebates has to be done numerically. The same is truefor American options.However, as mentioned in the Introduction, one can find the solution to these problems semi-analyticallyusing the HP and GIT methods for some one-factor models, including the time-dependent OU (Vasicek)model in (Carr and Itkin, 2021; Lipton and Kaushansky, 2020a,b), for the Hull-White model in (Itkinand Muravey, 2020a), for the CEV and CIR models in (Carr et al., 2020), and then in a general formfor any model that can be reduced to the heat equation - in (Itkin and Muravey, 2020b). In other words,solving these problems doesn’t require the ML method. Therefore, below we consider some other modelsfor which the barrier pricing problems cannot be directly solved by the HP or GIT methods but can besolved by using the ML method.
Page 12 of 36 ultilayer heat equations: application to finance
The Black-Karasinski (BK) model was introduced in (Black and Karasinski, 1991), see also (Brigo andMercurio, 2006) for a more detailed discussion. The BK is a one-factor short interest rate model of theform dz t = k ( t )[ θ ( t ) − z t ] dt + σ ( t ) dW t , r ∈ R , t ≥ , (37) r t = s ( t ) + Re z t , r ( t = 0) = r . Here κ ( t ) > θ ( t ) is the mean-reversion level, σ ( t ) is the volatility, R is some constant with the same dimensionality as r t , eg., it can be 1/(1 year). This model is similar tothe Hull-White model but preserves the positivity of r t by exponentiating the Ornstein-Uhlenbeck (OU)random variable z t . Because of that, usually, practitioners add a deterministic function (shift) s ( t ) to thedefinition of r t to address possible negative rates and be more flexible when calibrating the term-structureof the interest rates.By Itô’s lemma the short rate ¯ r t = ( r t − s ( t )) /R in the BK model solves the following stochasticdifferential equation (SDE) d ¯ r t = [ kθ ( t ) + 12 σ ( t ) − k log ¯ r t ]¯ r t dt + σ ( t )¯ r t dW t . (38)This SDE can be explicitly integrated. Let 0 ≤ s ≤ t ≤ T , Then r t can be represented as, (Brigo andMercurio, 2006)¯ r t = exp (cid:20) e − k ( t − s ) log ¯ r s + k Z ts e − k ( t − u ) θ ( u ) du + Z ts σ ( u ) e − k ( t − u ) dW ( u ) (cid:21) , (39)and thus, conditionally on filtration F s is lognormally distributed and always positive.However, in the BK model, the price P ( t, T ) of a (ZCB) with the maturity T is not known in closedform since this model is not affine. Multiple good approximations have been developed in the literatureusing asymptotic expansions of various flavors, see, e.g., (Antonov and Spector, 2011; Capriotti andStehlikova, 2014; Horvath et al., 2017), and also survey in (Turfus, 2020).Despite this lack of tractability, the BK model is widely used by practitioners for modeling interestrates and credit and is also known in commodities as the Schwartz one-factor model. The BK model isattractive because it is relatively simple, guarantees non-negativity of the prices (which could be a badfeature in the current environment). It could also be calibrated to the given term-structure of interestrates and the prices or implied volatilities of caps, floors, or European swaptions since the mean-reversionlevel and volatility are functions of time. However, for exotic options, e.g., highly liquid barrier options,these prices are not known yet in closed form. Therefore, various numerical methods are used to obtainthem.Here we describe how one can reduce the pricing problem for the ZCB to the ML heat equation. Sincethis problem is defined at a semi-infinite domain, the corresponding ML heat equation is also definedat a semi-infinite interval. Thus, the number of layers could be infinite. Therefore, truncation of thesemi-infinite interval to a finite is needed. Of course, the impact of the remainder should be assessedappropriately.Along the BK model lines, consider a model where the dynamics of the underlying stochastic variable z t is the OU process defined in Eq. (37). We assume that the interest rate r t is some deterministicfunction of z t r t = s ( t ) + f ( t, z t ) , z = 0 . (40)In particular, according to Eq. (37) for the BK model we have f ( t, z t ) = Re z t , and so R = r − s . Page 13 of 36 ultilayer heat equations: application to finance
In terms of z , the corresponding PDE for the ZCB price F ( t, r ) in Eq. (49) and for the option price C ( t, r ) in Eq. (52) reads0 = ∂V∂t + 12 σ ( t ) ∂ V∂z + κ ( t )[ θ ( t ) − z ] ∂V∂z − [ s ( t ) + f ( t, z )] V, (41)where V = V ( t, z ) is either F ( t, z ) or C ( t, z ). This equation should be solved subject to the terminaland boundary conditions. For the ZCB price they are given in Eq. (50) and Eq. (51), and for theDown-and-out barrier Call option price - in Eq. (53) and in Eq. (55), Eq. (56). Note, that solvingEq. (41) for F ( t, z ) assumes that z ∈ (cid:0) f − ( t, −∞ ) , f − ( t, ∞ ) (cid:1) , while for C ( t, z ) the domain of definitionis z ∈ [ L z ( t ) , f − ( t, ∞ ))., where L z ( t ) = f − ( t, L ( t )), and f − ( t, r ) is the inverse function.To apply the ML method to Eq. (41), for instance, when solving the barrier option pricing problem,we truncate the interval [ L ( t ) , ∞ ) from above at z = z max to make it [ L z ( t ) , z max ]. The reason thisis possible lies in the fact that when z increases, the ZCB price tends to zero based on the boundarycondition. Therefore, the Call option price in this limit vanishes as well. Thus, the contribution of theregion [ z max , f − ( t, ∞ )) to the Call option price becomes negligible .Now we split the interval [ L z ( t ) , z max ] into N > z i , z i +1 ] , i =1 , . . . , N assume that f ( t, z ) = a i ( t ) + b i ( t ) z ). Accordingly, at every interval i, i = 1 , . . . , N the PDEEq. (41) takes the form0 = ∂V∂t + 12 σ ( t ) ∂ V∂z + κ ( t )[ θ ( t ) − z ] ∂V∂z − [ s ( t ) + a i ( t ) + b i ( t ) z ] V. (42)This PDE can be transformed to the heat equation ∂U∂τ = ∂ U∂x , (43)by the change of variables, (Polyanin, 2002; Itkin and Muravey, 2020a; Lipton and Kaushansky, 2020b) V ( t, z ) = exp[ α i ( t ) z + β i ( t )] U ( τ, x ) , τ = φ ( t ) , x = zψ ( t ) + ̺ ( t ) , (44)where ψ ( t ) = C exp (cid:18)Z tS κ ( q ) dq (cid:19) , φ ( t ) = 12 Z St σ ( q ) ψ ( q ) dq + C , (45) α i ( t ) = ψ ( t ) Z tS b i ( q ) ψ ( q ) dq + C ψ ( t ) , ̺ i ( t ) = − Z tS h κ ( q ) θ ( q ) + σ ( q ) α i ( q ) i ψ ( q ) dq + C ,β i ( t ) = − Z tS α i ( q ) h κ ( q ) θ ( q ) + σ ( q ) α i ( q ) i dq + Z tS [ s ( q ) + a i ( q )] dq + C , where C , . . . , C are some constants. In our case we can choose C = 1, C = C = C = C = 0.One of the advantages of such an approach is that the new time τ doesn’t depend on the specificinterval i , i.e. the time τ runs in sync for all intervals [ z i , z i +1 ] , i = 1 , . . . , N .Thus, the main trick here is in using the approximation f ( t, z ) = a i ( t ) + b i ( t ) z ). This approximationprovides the second order of accuracy in the length of the interval (similar to the finite-difference methodof the second order), and allows reduction of the PDE at each interval to the heat equation (while theoriginal PDE doesn’t hold this property).At the end of this Section, note that Eq. (42), if used for pricing ZCB, doesn’t need the boundarycondition at the left boundary z → −∞ , as this is discussed in (Itkin and Muravey, 2020a) with a For some choices of the functions f ( t, z ) the value f − ( t, ∞ ) could be finite which eliminates the need for truncation.Page 14 of 36 ultilayer heat equations: application to finance reference to Fichera theory, (Oleinik and Radkevich, 1973). However, the price of the ZCB at some fixedleft boundary z min , i.e. V ( t, z min ) can be found having in mind that the transformed PDE in Eq. (42) isaffine, which yields V ( z, t, S ) = A ( t, S ) e B ( t,S ) Re z . (46)With allowance for the terminal condition in Eq. (50), the solution reads, (Itkin and Muravey, 2020a) B ( t, S ) = e R t κ ( m ) dm Z tS b i ( m ) e − R m κ ( q ) dq dm, (47) A ( t, S ) = exp (cid:20)Z tS (cid:18) a i ( m ) + s ( m ) − B ( m, S ) (cid:16) θ ( m ) κ ( m ) + B ( m, S ) σ ( m ) (cid:17)(cid:19) dm (cid:21) . It can be seen that B ( t, S ) < t < S . Therefore, F ( r, t, S ) → z → ∞ . Since the BK model is not fully tractable, in (Itkin et al., 2020) we introduced a slightly modified versionof the model as follows dz t = k ( t )[¯ θ ( t ) − e z t ] dt + σ ( t ) dW t , (48) r t = s ( t ) + Re z t , z = 0 , R = r − s (0) . It can be seen, that at small t | z t | ≪
1, and so choosing ¯ θ ( t ) = 1 + θ ( t ) replicates the BK model in thelinear approximation on z t . Similarly, the choice ¯ θ ( t ) = e θ ( t ) replicates the BK model at z t close themean-reversion level θ ( t ). Thus, this model acquires the properties of the BK model while is a bit moretractable as this will be seen below.It is worth noting that if by using Itô’s lemma we re-write Eq. (48) for the stochastic variable r t , theresulting dynamics can be recognized as the stochastic Verhulst or stochastic logistic model, which arewell-known in the population dynamics and epidemiology; see, eg., (Verhulst, 1838; Bacaer, 2011; Gietet al., 2015) and references therein. For more information, see (Itkin et al., 2020).By the Itô’s lemma and the Feynman–Kac formula any contingent claim written on the r t as theunderlying (for instance, price F (¯ r, t, S ) of a Zero-coupon bond (ZCB) with maturity S ) solves thefollowing partial differential equation0 = ∂F∂t + 12 σ ( t )¯ r ∂ F∂ ¯ r + κ ( t )¯ r [Θ( t ) − ¯ r ] ∂F∂ ¯ r − ( s ( t ) + R ¯ r ) F, (49)¯ r t = r t − s ( t ) r − s (0) = e z t , Θ( t ) = ¯ θ ( t ) + 12 σ ( t ) . This equation should be solved subject to the terminal condition F (¯ r, S, S ) = 1 , (50)and the boundary condition F (¯ r, t, S ) (cid:12)(cid:12)(cid:12) ¯ r →∞ = 0 , (51)see, eg., (Andersen and Piterbarg, 2010).In the sequel we will also consider a Down-and-Out barrier Call option written on the ZCB. It isknown, (Andersen and Piterbarg, 2010), that under a risk-neutral measure the option price C ( t, ¯ r ) solvesthe same PDE as in Eq. (49),0 = ∂C∂t + 12 σ ( t )¯ r ∂ C∂ ¯ r + κ ( t )¯ r [Θ( t ) − ¯ r ] ∂C∂ ¯ r − ( s ( t ) + R ¯ r ) C. (52) Page 15 of 36 ultilayer heat equations: application to finance
The terminal condition at the option maturity T ≤ S for this PDE reads C ( T, ¯ r ) = ( F (¯ r, T, S ) − K ) + , (53)where K is the option strike.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 ) . (54)This condition can be translated into the ¯ r domain by solving the equation F (¯ r, t, S ) = L F ( t ) , with respect to ¯ r . Denoting the solution of this equation as L ( t ) we find that Eq. (54) in the ¯ r domainreads C ( t, L ( t )) = 0 . (55)The second boundary can be naturally set at ¯ r → ∞ . As at ¯ r → ∞ the ZCB price tends to zero, the Calloption price also vanishes in this limit. This yields C ( t, ¯ r ) (cid:12)(cid:12)(cid:12) ¯ r →∞ = 0 . (56)Accordingly, Eq. (52) has to be solve at ¯ r ∈ [ L ( t ) , ∞ ). As we have already mentioned, in contrast to other similar one-factor models like the time-dependentOrnstein-Uhlenbeck, Hull-White, CIR and CEV models which have been considered in (Carr and Itkin,2021; Itkin and Muravey, 2020a; Carr et al., 2020), the solution of the pricing problem for the BK modelis not known in closed form . Therefore, we propose an approximation that gives rise to a semi-analyticalsolution for the barrier Call option price. This approximation is inspired by the ML heat equations whichare discussed in Section .Since our problem in Eq. (52) is defined at the semi-infinite domain ¯ r ∈ [ L ( t ) , ∞ ), using the MLapproximation is time-consuming, as we need to split this semi-infinite interval into a fixed number ofsub-intervals. Therefore, it is feasible first to make a change of variables C ( t, ¯ r ) = V ( t, x ) e R t s ( k ) dk , x = a ( t )¯ r , a ( t ) = e R t ( κ ( m )Θ( m ) − σ ( m ) ) dm , (57)so the problem to solve in new variables reads0 = ∂V∂t + 12 σ ( t ) x ∂ V∂x + a ( t ) κ ( t ) ∂V∂x − Ra ( t ) Vx . (58)Thus, now our problem is defined at a fixed domain x ∈ [0 , /L ( t )], where the upper boundary is timedependent. Accordingly, in the new variables the Down-and-Out barrier Call option becomes the Up-and-Out barrier Call.Doing the second change of the dependent variable V ( t, x ) = u ( t, x ) e d ( t ) /x , d ( t ) = Re − R t σ ( m ) dm Z t e R y κ ( m )Θ( m ) dm dy (59) Some approaches for doing that are discussed in (Itkin et al., 2020). Page 16 of 36 ultilayer heat equations: application to finance yields the equation 0 = ∂u∂t + 12 σ ( t ) x ∂ u∂x + g ( t ) ∂u∂x − f ( t ) ux , (60) f ( t ) = 12 d ( t ) (cid:16) a ( t ) κ ( t ) − d ( t ) σ ( t ) (cid:17) , g ( t ) = a ( t ) κ ( t ) − d ( t ) σ ( t ) . Accordingly, in the new variables the initial and boundary conditions read u ( T, x ) = exp − d ( T ) x − Z T s ( k ) dk ! ( F ( x, T, S ) − K ) + , (61) u ( t,
0) = 0 , u ( t, y ( t )) = 0 , y ( t ) = a ( t ) /L ( t ) . The problem in Eq. (60), Eq. (61) cannot be solved in closed form. Therefore, we proceed by borrowingthe idea from the ML approach in physics which is described in Sections 1,2. This approach impliesthat the interval x ∈ [0 , y ( t )] we approximate the function ζ ( x ) = x by using a piecewise constantapproximation. In more detail, we split the interval [0 , y ( t )] into N > x i , x i +1 ] , i = 1 , . . . , N assume that x ≈ ν i ( t ) . For instance, one can choose the middle valueof the function ζ ( x ) at each sub- , so ν i ( t ) = y ( t ) ( i + 1 / N . With allowance for this approximation, at every i -th interval Eq. (60) takes the form0 = ∂u∂t + 12 σ ( t ) ν i ( t ) ∂ u∂x + g ( t ) ∂u∂x − f ( t ) ν i ( t ) u. (62)The Eq. (62) can be transformed to the heat equation ∂U∂τ = ∂ U∂ς , (63)using the transformation, (Polyanin, 2002) u ( t, x ) = U ( τ, ς ) exp (cid:18) − Z t f ( t ) ν i ( t ) (cid:19) , l = x − Z t g ( k ) dk, τ = 12 Z Tt σ ( k ) ν i ( k ) dk. (64)Note, that under this approximation the new time τ also becomes a function of the interval i . Calibration of the local volatility model (constructed by using a one-factor Geometric Brownian motionprocess) to a given set of option prices is a classical problem of mathematical finance. It was considered inmultiple papers, and various solutions were proposed; see, e.g., a survey in (Itkin, 2020; Itkin and Lipton,2018) and references therein. In particular, in (Itkin and Lipton, 2018) an analytical approach to solvingthe calibration problem is developed. This approach extends the method in (Lipton and Sepp, 2011)by replacing a piecewise constant local variance construction with a piecewise linear one and allowingnon-zero interest rates and dividend yields. This approach remains analytically tractable as it combinesthe Laplace transform in time with an analytical solution of the resulting spatial equations in terms ofKummer’s degenerate hypergeometric functions. Since the upper boundary of the whole interval y ( t ) is the function of time, we need to put ν i = ν i ( t ) Page 17 of 36 ultilayer heat equations: application to finance A similar problem could be formulated not just for the Black-Scholes model but also for other models.For instance, in (Carr and Itkin, 2020, 2019) two extensions of the Local Variance Gamma model proposedinitially in (Carr and Nadtochiy, 2017) were developed. The first new model (ELVG) considers a Gammatime-changed arithmetic
Brownian motion with drift and the local variance to be a piecewise linearfunction of the strike. The second model (GLVG) is a geometric version of the ELVG with drift. Italso treats various cases by introducing three piecewise linear models: the local variance as a function ofstrike, the local variance as a function of log-strike, and the local volatility as a function of strike (so, thelocal variance is a piecewise quadratic function of the strike). For all these extensions, the authors derivean ordinary differential equation for the option price, which plays the role of Dupire’s equation for thestandard local volatility model. Moreover, it can be solved in closed form.In (Itkin and Lipton, 2018; Carr and Itkin, 2020, 2019) all models were calibrated to the marketquotes term-by-term. Therefore, various types of no-arbitrage interpolation were proposed to guaranteeno-arbitrage while keeping the model analytically tractable on the other hand; further details are givenin (Itkin, 2020).Two advantages of the semi-analytical approach, which are essential for calibration the model, shouldbe emphasized. First, the option prices can be found analytically in a semi closed form. Here "semi"means that the analytic solution requires an additional inverse transform to be applied to get the finalprices; see (Itkin and Lipton, 2018). However, in (Carr and Itkin, 2020, 2019) since the Dupire-likeequation is an ODE and not a PDE, this step is eliminated. Nevertheless, all these models are calibratedterm-by-term.This idea can be extended by constructing a semi-analytical solution of the ML heat equation, whichis analytic in time. Thus, the term-by-term calibration could be eliminated, and the quotes for all strikesand maturities can be used simultaneously. Therefore, this approach allows a further acceleration of thecalibration process.For brevity, let us consider European options, for instance, a Put option on a stock. It is well-knownthat the price P ( T, K ) of the option written on the underlying stock price S t as a function of the optionmaturity T and strike K solves Dupire’s equation, (Dupire, 1994) ∂P∂T = 12 σ ( T, K ) ∂ P∂K − ( r − q ) K ∂P∂K − q ( T ) P. (65)with σ = σ ( t, S ) being the local volatility, and q ( t ) is the dividend yield. This PDE should be solvedsubject to the initial condition at T = 0 P (0 , K ) = ( K − S ) + , (66)and natural boundary conditions for the put option price that read, (Hull, 2011) P ( T, K ) = 0 , K → ,P ( T, K ) = D ( t, T )( K − S ) ≈ D ( t, T ) K, K → ∞ , (67)where the discount factor D ( t, T ) is defined as D ( t, T ) = exp − Z Tt r ( k ) dk ! . (68)To proceed further, we use the idea in (Lipton and Sepp, 2011; Itkin and Lipton, 2018) and approx-imate the local variance using some piecewise approximation in the strike space. However, in contrastto (Lipton and Sepp, 2011; Itkin and Lipton, 2018) we make this approximation a function of time. Fur-ther, suppose that for each trading maturity T j , j ∈ [1 , M ] the market quotes are provided at a set Page 18 of 36 ultilayer heat equations: application to finance of strikes K i,j , i = 1 , . . . , n j where the strikes are assumed to be sorted in the increasing order. Let usconstruct a finite grid in the strike space G ( K ) : K ∈ [min( K i,j ) , max( K i,j )] , j = 1 , . . . , M, i = 1 , . . . , n j ,by splitting the whole interval K ∈ [min( K i,j ) , max( K i,j )] into n sub-intervals. At every interval[ K i , K i +1 ] , i = 1 , . . . , N , we approximate the local variance function σ ( T, K ) by a piecewise constantfunction in K as follows: σ ( T, K ) = v i ( T ) , K ∈ [ K i , K i +1 ] . (69)This approximation is not continuous, so the local variance σ ( T, K ) experience a finite jump at everypoint K i . However, it is continuous in the maturity T .Accordingly, at every interval [ K i , K i +1 ] Eq. (65), V ( T, K ) takes the form ∂P i ∂T = 12 v i ( T ) ∂ P i ∂K − [ r ( T ) − q ( T )] ∂P i ∂K − q ( T ) P i . (70)This equation can be transformed to the heat equation ∂U i ∂τ = ∂ U i ∂x , (71)by a change of variables P i ( T, K ) = U i ( τ, x ) e − R T ( q ( s ) ds , x = e − R T ( r ( s ) − q ( s ) ds K, τ = 12 Z T v i ( s ) e − R s ( r ( s ) − q ( s )) dk ds. (72)It is important to note that the new time τ runs differently at every interval K ∈ [ K i , K i +1 ] as it dependson the local variance value v i ( s ) at this interval. An efficient solution of the derived systems of Volterra equations is a problem that requires some attentionand extended description. Therefore, it will be published elsewhere. Instead, here we show that aparticular choice of the internal boundaries can reduce this problem to a linear system with a tridiagonalmatrix allowing the inverse Laplace transform. In this section, we explain this approach in detail.We start with Eq. (21). Using the definitions of η ± , υ ± , υ ± in Eq. (22), we observe that under allintegrals in Eq. (21) the functions f i ( s ) , Θ i ( s ) , Ω i ( s ) are functions of s , while functions η ± , υ ± , υ ± arefunctions of t − s and y i ( t ) − y i ( s ) and z i ( t ) − z i ( s ). Recall that functions y ( t ) = χ − , y N ( t ) = χ + definethe external boundaries of the computational domain, while functions y i ( t ) , i = 2 , . . . , N − y i ( s ) = a i s + b i s + c i , where a i , b i , c i are some constants. Then y i ( t ) − y i ( s )can also be represented as a certain function g ( t − s ). Indeed y i ( t ) − y i ( s ) = − a i ( t − s ) + ( b i + 2 a i t )( t − s ) . (73)A similar representation can be obtained for y i +1 ( s ) − y i ( s ) y i +1 ( s ) − y i ( s ) = ( a i +1 − a i ) s + ( b i +1 − b i ) s + ( c i +1 − c i ) = A ( t − s ) + B ( t − s ) + C, (74) A = a i +1 − a i , B = 2 At + b i +1 − b i , C = c i − c i +1 + t [ b i − b i +1 + ( a i − a i +1 ) t ] . Same can be done for a polynomial of any degree.All coefficients a i , b i , c i can be precomputed given the external boundaries. An example of this con-struction is given in Fig. 1. Page 19 of 36 ultilayer heat equations: application to finance • y N ( t ) y ( t ) y i ( t ) y i +1 ( t ) y i − ( t ) ...... ...... Figure 1:
Internal layers constructed for the given external boundaries y ( t ) and y N ( t ) , andthe number of layers N , by using 3 points for each boundary y i ( t ) and polynomial curves. In more detail, suppose that for the given functions y ( t ) = χ − ( t ) , y N +1 = χ + ( t ), we want to have N layers. We use N uniform nodes to split the interval [ χ − (0) , χ + (0)] into N subintervals. We do thesame for the interval [ χ − ( t ) , χ + ( t )]. If the boundaries χ − ( t ) , χ + ( t ) are smooth enough, we can connectpoints y i (0) , y i ( t ) , i = 2 , . . . , N by straight lines in such a way that all boundaries don’t cross each other.Suppose this is not possible because the external boundaries are too convex or concave. In that case,we can find some s = τ where the distance between the external boundaries is minimal and put N nodes there. Then, we can connect all points y i (0) , y i ( τ ) , y i ( t ) by parabolas, and again check that all theboundaries don’t cross each other. We can continue this process by using polynomials of a higher degreeto provide the final representation of the boundaries.Thus, we can find a polynomial of the necessary degree to guarantee that all boundaries don’t crosseach other unless the external curves have a very peculiar shape, which does not happen in the context offinancial applications. Otherwise, we need to use a general approach to the computation of the integralsin Eq. (21), which will be published elsewhere.Provided that all the boundaries are constructed in such a way, one can observe that all integrals inEq. (21) are convolutions. Therefore, we can apply the Laplace transform L ( f | s, λ ) L ( f | s, λ ) = Z ∞ e − λs f ( s ) ds (75)to both parts of each equation in Eq. (21). Taking into account that L ( f ∗ g ) = L ( f ) L ( g ) , L ( f ′ ) = λ L ( f ) − f (0) , (76)and obtain from Eq. (21) a linear system of equations for L ( χ + i ) , L ( χ − i ) , L (Ω i ) , L (Θ i ). Using the conditionsat the internal boundaries (such as in Eq. (35)), this system can be reduced to a linear system for only Page 20 of 36 ultilayer heat equations: application to finance L ( χ − i ) , L (Ω i ). One can check that the resulting system is block-diagonal, with all blocks being tridiagonalmatrices. Once this system is solved, the functions χ − i ( t ) , Ω i ( t ) is found by using the inverse Laplacetransform.Consider the heat equation in a strip with a piecewise constant thermal conductivity coefficient toillustrate our approach. Consider the following problem for the heat equation with a piecewise constant thermal conductivitycoefficient (the ML problem): ∂∂x (cid:18) σ ( x ) ∂U∂x (cid:19) = ∂U∂t , ( x, t ) ∈ [ y , y N ] × R + (77) U ( t, y ) = 0 , U ( t, y N ) = 0 .U (0 , x ) = δ ( x − x ) . Here the thermal conductivity coefficient σ ( x ) is a piecewise constant function of x , which changesfrom layer to layer σ ( x ) = N X i =1 ( y i − < x ≤ y i ) σ i , (78) y < y < y < ... < y i < ... < y N . and δ ( x ) denotes the Dirac delta function. As before y i are the boundaries of the layers in x space.Without loss of generality, we assume that x ∈ [ y j − , y j ) , < j < N .Due to the initial condition in Eq. (77) the solution of this problem is Green’s function for Eq. (77).We represent the solution U ( t, x ) in the form U ( t, x ) = N X i =1 ( y i − < x ≤ y i )[ U i ( t, x ) + H i ( t, x )] , (79)where the functions U i ( t, x ) and H i ( t, x ) solve the following problems ∂∂x (cid:18) σ i ∂U i ∂x (cid:19) = ∂U i ∂t , ( x, t ) ∈ ( y i − , y i ] × R + , (80)lim x → y i − U i ( t, x ) = χ − i ( t ) , U i ( t, y i ) = χ + i ( t ) ,U (0 , x ) = 0 , and ∂∂x (cid:18) σ i ∂H i ∂x (cid:19) = ∂H i ∂t , ( x, t ) ∈ ( y i − , y i ] × R + , (81)lim x → y i − H i ( t, x ) = 0 , H i ( t, y i ) = 0 ,H (0 , x ) = δ ( x − x ) . A well-known physical argument shows that the solution and its flux must be continuous at the layers’boundaries. The first condition yields U ( t, y ) = 0 , (82) Page 21 of 36 ultilayer heat equations: application to finance lim x → y i U i ( t, x ) = U i +1 ( t, y i ) , i = 1 ...N − , lim x → y N U N ( t, x ) = 0 . According to (Itkin and Muravey, 2020b), the function H ( t, x ) = 0 only at that interval which containsthe point x , i.e. [ y j − , y j ). Therefore, the flux continuity conditions could be written aslim x → y i σ i ∂U i ∂x ( t, y i ) = σ i +1 ∂U i +1 ∂x ( t, y i ) , i = j − , j, (83)lim x → y j − σ j − ∂U j − ∂x ( t, y j − ) = σ j ∂U j ∂x ( t, y j − ) + σ j ∂H j ∂x ( t, y j − ) , lim x → y j (cid:20) σ j ∂U j ∂x ( t, y j ) + σ j ∂H j ∂x ( t, y j ) (cid:21) = σ j +1 ∂U j +1 ∂x ( t, y j ) . It follows from Eq. (82) that χ − ( t ) = χ + N ( t ) = 0 , χ + i ( t ) = χ − i +1 ( t ) , i = 1 , . . . , N − . (84)To simplify the notation let us introduce new functions f i ( t ), such as χ + i ( t ) = χ − i +1 ( t ) = f i ( t ) , i = 0 , . . . , N, (85)so obviously f ( τ ) = f N +1 ( τ ) = 0. Using Eq. (21) (or Eq. 3.33 in (Itkin and Muravey, 2020b)), one canget an explicit representation for the derivatives of U ( t, x ) at each interval ∂U i ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i = f i ( t ) σ i √ πt − Z t f i ( s ) − f i ( t )2 σ i q π ( t − s ) ds + Z t h f i − ( s ) λ + i ( t | y i − , s ) − f i ( s ) λ +0 ,i ( t | y i , s ) i ds, (86) ∂U i +1 ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i = − f i ( t ) σ i +1 √ πt + Z t f i ( s ) − f i ( t )2 σ i +1 p π ( t − s ) ds + Z t h f i ( s ) λ − ,i +1 ( t | y i , s ) − f i +1 ( s ) λ − i +1 ( t | y i +1 , s ) i ds. Here λ − i ( t | ξ, s ) = ∞ X n = −∞ e − ( yi − ξ +2 nli )24 σ i ( t − s ) σ i p π ( t − s ) " − ( y i − ξ + 2 nl i ) σ i ( t − s ) = − π l i θ ′′ (cid:20) π ( y i − ξ )2 l i , q i ( s ) (cid:21) , (87) λ + i ( τ | ξ, s ) = ∞ X n = −∞ e − ( yi − ξ +(2 n +1) li )24 σ i ( t − s ) σ i p π ( t − s ) " − ( y i − ξ + (2 n + 1) l i ) σ i ( t − s ) = − π l i θ ′′ (cid:20) π ( y i + l i − ξ )2 l i , q i ( s ) (cid:21) ,λ − ,i ( t | ξ, s ) = λ − i ( t | ξ, s ) − e − ( yi − − ξ )24 σ i ( t − s ) σ i p π ( t − s ) " − y i − − ξ σ i ( t − s ) ,λ +0 ,i ( t | ξ, s ) = λ + i ( t | ξ, s ) − e − ( yi − ξ + li )24 σ i ( t − s ) σ i p π ( t − s ) " − ( y i − ξ + l i ) σ i ( t − s ) .q i ( s ) = e − π σ i ( t − s ) /l i . Here θ i ( z, p ) , i = 2 , θ ′′ ( z, p ) is the second derivative of θ ( z, p ) on the first argument. Page 22 of 36 ultilayer heat equations: application to finance
Substituting ξ = y i − , y i , y i +1 into Eq. (87) we obtain λ − i +1 ( t | y i +1 , s ) = − π l i +1 θ ′′ (0 , q i +1 ( s )) , (88) λ + i ( τ | y i − , s ) = − π l i θ ′′ ( π, q i ( s )) ,λ − ,i +1 ( t | y i , s ) = ∞ X n = −∞ n =0 e − (2 nli +1)24 σ i +1( t − s ) σ i +1 p π ( t − s ) " − (2 nl i +1 ) σ i +1 ( t − s ) ,λ +0 ,i ( t | y i , s ) = − π l i θ ′′ ( π, q i ( s )) − e − l i σ i ( t − s ) σ i p π ( t − s ) " − l i σ i ( t − s ) . Since at s → t we have q ( s ) →
1, all RHSs in Eq. (88) are regular in this limit and vanish. Thelatter is due to the fact that lim q → θ ′′ (0 , q ) = θ ′′ ( π, q ) = 0. Therefore, all integral kernels in Eq. (86)are regular. We also assume that functions f i ( s ) are smooth enough, so that in the limit all the integralsvanish.Applying integration by parts to Eq. (86), we get the following simplified system ∂U i ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i − = − Z t h η eveni ( t, s ) d ( f i − ( s )) − η oddi ( t, s ) d ( f i ( s )) i (89) − f i − (0) σ i √ πt − f i − (0) η eveni ( t,
0) + f i (0) η oddi ( t, ,∂U i ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i = − Z t h η oddi ( t, s ) d ( f i − ( s )) − η eveni ( t, s ) d ( f i ( s )) i + f i (0) σ i √ πt + f i (0) η eveni ( t, − f i − (0) η oddi ( t, , where η eveni ( t − s ) = 1 σ i p π ( t − s ) ∞ X n = −∞ e − (2 nli )24 σ i ( t − s ) , η oddi ( t − s ) = 1 σ i p π ( t − s ) ∞ X n = −∞ e − ((2 n +1) li )24 σ i ( t − s ) . (90)Indeed, since η − i ( t | y i − , t ) = 0 , η − i ( t | y i , t ) = 0 , f i − (0) = 0 , f i (0) = 0 we have − f i − ( t ) σ i √ πt + Z t f i − ( s ) − f i − ( t )2 σ i p π ( t − s ) + Z t h f i − ( s ) d (cid:16) η − i ( t | y i − , s ) (cid:17) − f i ( s ) d (cid:16) η − i ( t | y i , s ) (cid:17)i (91)= − f i − ( t ) σ i √ πt − f i − ( t ) σ i p π ( t − s ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ts =0 + f i − ( s ) σ i p π ( t − s ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ts =0 − Z t σ i p π ( t − s ) d ( f i − ( s )) ++ f i − ( s ) η − i ( t | y i − , s ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ts =0 − f i ( s ) η − i ( t | y i , s ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ts =0 − Z t h η − i ( t | y i − , s ) d ( f i − ( s )) − η − i ( t | y i , s ) d ( f i ( s )) i = − Z t " η − i ( t | y i − , s ) + 1 σ i p π ( t − s ) ! d ( f i − ( s )) − η − i ( t | y i , s ) d ( f i ( s )) = − Z t h η eveni ( t, s ) d ( f i − ( s )) − η oddi ( t, s ) d ( f i ( s )) i . Page 23 of 36 ultilayer heat equations: application to finance
In turn, as shown in (Itkin and Muravey, 2020b), Eq. 3.31, the function H j ( t, x ) can be representedas follows H j ( t, x ) = 12 σ j √ πt ∞ X n = −∞ e − (2 nlj + x − x σ j t − e − (2 nlj + x + x − yj − σ j t . (92)Hence, the gradients at the boundaries are ∂H j ( t, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y j − ≡ υ − ( t | x ,
0) (93)= ∂∂x ( l j θ " π ( x − x )2 l j , q j (0) − l j θ " π ( x + x − y j − )2 l j , q j (0) x = y j − = 14 l j ( θ ′ " π ( y j − − x )2 l j , q j (0) − θ ′ " π ( x − y j − )2 l j , q j (0) ,∂H j ( t, x ) ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y j = 14 l j ( θ ′ " π ( y j − x )2 l j , q j (0) − θ ′ " π ( x − y j )2 l j , q j (0) . Thus, we obtain the following system of Volterra equations Z t " − σ i η oddi ( t − s ) d ( f i − ( s )) + (cid:16) σ i η eveni ( t − s ) + σ i +1 η eveni +1 ( t − s ) (cid:17) d ( f i ( s )) (94) − σ i +1 η oddi +1 ( t − s ) d ( f i +1 ( s )) = h i ( t ) ,h i ( t ) = 0 , i = j − , j, h j − ( t ) = σ j υ − ( t | x , , h j ( t ) = − σ j υ + ( t | x , . Since the kernels depend only on t − s one can rewrite the above equations as a convolution (cid:16) − σ i η oddi ( · ) ∗ f ′ i − ( · ) + h σ i η eveni ( · ) + σ i +1 η eveni +1 ( · ) i ∗ f ′ i ( · ) − σ i +1 η oddi +1 ( · ) ∗ f ′ i +1 ( · ) (cid:17) ( t ) = h i ( t ) . (95) Applying the Laplace transform to Eq. (95), we get √ λ (cid:16) − σ i L ( η oddi ) L ( f i − ) + (cid:16) σ i L ( η eveni ) + σ i +1 L ( η eveni +1 ) (cid:17) L ( f i ) − σ i +1 L ( η oddi +1 ) L ( f i +1 ) (cid:17) = L ( h i ) √ λ , (96)or, in the matrix form Mg = σ j √ λ h L ( υ − ( t | x , j − − L ( υ + ( t | x , j i . (97)Here j denotes the indicator vector, i.e. j = , , .. | {z } j − , , , ... ⊤ , (98)the vector g is the column vector g = ( L f , . . . , L f N − ) ⊤ , (99) Page 24 of 36 ultilayer heat equations: application to finance and the matrix M is a symmetric tridiagonal matrix M = D − β − β D − β − β . . . . . .. . . . . . − β N − − β N − D N − , (100)Coefficients of the matrix M have the form β i = √ λσ i +1 L ( η oddi +1 ) , D i = √ λ h σ i L ( η eveni ) + σ i +1 L ( η eveni +1 ) i , (101)and can be found explicitly, see Appendix C L ( η eveni ) = 1 σ i √ λ coth √ λl i σ i ! , L ( η oddi ) = 1 σ i √ λ (cid:16) √ λl i σ i (cid:17) , (102) L ( υ − ( t | x , σ j sinh (cid:18) ( y j +1 − x ) √ λσ j (cid:19) sinh (cid:18) l j √ λσ j (cid:19) , L ( υ + ( t | x , − σ j sinh (cid:18) ( x − y j ) √ λσ j (cid:19) sinh (cid:18) l j √ λσ j (cid:19) . Finally, introducing the notation ω i = l i σ i , γ = y j +1 − x l j , γ = x − y j l j , (103)the system Eq. (97) can be represented as Mg = 1 √ λ sinh (cid:16) γ ω j √ λ (cid:17) sinh (cid:16) ω j √ λ (cid:17) j − + sinh (cid:16) γ ω j √ λ (cid:17) sinh (cid:16) ω j √ λ (cid:17) j , (104) D i = σ i coth (cid:16) ω i √ λ (cid:17) + σ i +1 coth (cid:16) ω i +1 √ λ (cid:17) , β i = σ i +1 sinh (cid:16) ω i +1 √ λ (cid:17) . In this section, we solve the problem in Eq. (77) by using the ML method. Note that such problemsappear both in physics and in finance. A simple example of a financial problem is finding Green’s functionfor pricing double barrier options written on the underlying S t with local volatility dS t = σ ( S t ) dW t , (105)where σ ( S ) is a piecewise constant function. Below we describe two numerical experiments. σ i To start with, we assume that σ i = const, i = 1 , . . . , N , and hence, σ i in Eq. (77) can be pulled out of thederivative in x . This problem has an analytic solution, see (Lipton, 2001) and references therein, which Page 25 of 36 ultilayer heat equations: application to finance can be represented as the Fourier series. Re-writing it by using the definition of Jacobi theta functionsyields U ( T, y ) = 12 l (cid:20) θ (cid:18) π ( y − x )2 l , q (cid:19) − θ (cid:18) π ( y + x − y )2 l , q (cid:19)(cid:21) , (106) l = y N − y , q = e − π σ l T . To solve this problem, we need first to solve the linear system in Eq. (104) numerically, and then usethe corresponding Laplace images to find the function f i = f ( y i ) , i = 1 , . . . , N by applying an inverseLaplace transform. For the latter step, we use the Gaver-Stehfest method f ( T, y ) = Λ [ m ] X s =1 St ( m ) s L ( f i (Λ)) , Λ = log 2
T . (107)This algorithm was widely studied (see, e.g., (Kuznetsov, 2013) and references therein), and, providedthat the resulting function is non-oscillatory, converges very quickly. For instance, choosing m = 12 termsin the series representing the solution is usually sufficient. The coefficients St s can be found explicitly inadvance. Table 1:
Parameters of the test. y y N σ T N m -1.0 1.0 0.5 1.0 20 16The model parameters for this test are given in Table 1, and the results are depicted in Fig. 2a. Here,the left vertical axis shows the values of y i ( T ), and the right vertical axis shows the relative error (inpercent) of the solution compared with the analytic one. -1 -0.5 0 0.5 1 y V a li e o f t he G r een F un c t i on -0.015-0.01-0.00500.0050.010.015 R e l d i ff e r en c e , % AnalyticILTDiff (a) -1 -0.5 0 0.5 1 y V a l ue o f t he G r een F un c t i on -0.08-0.06-0.04-0.0200.020.04 R e l d i ff e r en c e , % The ML vs FD (N=41, M=40) methods
AnalyticFDILTDifILTDifFD (b)
Figure 2:
Comparison of the Analytic and ML solutions (a), and Analytic, ML and FDsolutions (grid with × nodes) (b) for σ i = 0 . , T = 1 . Here Analytic denotes the analyticsolution of the problem, ILT - the ML solution, FD - the FD solution, DiffILT - the relativeerror of the ML solution with respect to the analytic one, DifFD - same for the FD method. Page 26 of 36 ultilayer heat equations: application to finance
The ML solution coincides with the analytic one with high accuracy. The elapsed time of the experi-ment is 8.8 ms (we run our code, written in Matlab, on a PC with two Quadcore CPU Intel i7-4790 3.60GHz). The elapsed time doesn’t depend on the option maturity T , so the calculation is fast even for longmaturities. Note that the ML solution’s computation takes only 2.3 ms, while the remaining time is usedfor computing the Gaver-Stehfest coefficients (but those can be precomputed if so desired). Since thematrix M is tridiagonal, the ML method’s complexity is O ( mN ). By comparison, the complexity of thefinite difference (FD) method is O ( M N ), M is the number of time steps. Obviously, for long maturities M ≫ m , so the FD method is slower.To validate this, we also implemented an FD method to solve the same problem. The FD solver runson a uniform grid and is a Crank-Nicolson scheme after four steps, while for the first four steps it uses animplicit Euler scheme. In other words, we start with four Ranacher steps, see (Itkin, 2017) and referencestherein.To have the same spatial approximation in y , we need to run both the ML and FD methods withthe same number of nodes. While the ML method provides an accurate result even at N = 10, the FDmethod fails and needs at least 40 nodes to converge to the solution. Therefore, we choose N = 41. Thesame is true for the time step, so the minimal number of FD time steps in our experiment is M = 40,while the ML method provides the solution at t = T just at once. The results of the comparison of bothmethods with the analytic solution are given in Fig. 2b. Again, the left vertical axis shows the valuesof y i ( T ), and the right vertical axis is the relative error (in percent) of the solution compared with theanalytic one.It is clear that the accuracy of the FD method is worse than that of the ML method. By increasing N and M , one can improve the FD method’s accuracy, but it takes time. The elapsed time for the FDmethod with N = 41 , M = 40 is 41 ms, so it is about 18 times slower than the ML method. Obviously,for longer maturities, more time-steps are necessary, so the FD method becomes even slower.It is also known that for small T and volatilities, the FD method’s error increases. To illustrate thisfact, we run the same experiment, but now with T = 0 . , σ = 0 .
3. The results are presented in Fig. 3. -1 -0.5 0 0.5 1 y V a l ue o f t he G r een F un c t i on -1-0.8-0.6-0.4-0.200.20.40.60.8 R e l d i ff e r en c e , % The ML vs FD (N=41, M=40) methods
AnalyticFDILTDifILTDifFD (a) -1 -0.5 0 0.5 1 y V a l ue o f t he G r een F un c t i on -0.35-0.3-0.25-0.2-0.15-0.1-0.0500.05 R e l d i ff e r en c e , % The ML vs FD (N=151, M=150) methods
AnalyticFDILTDifILTDifFD (b)
Figure 3:
Comparison of the Analytic, ML and FD solutions for σ i = 0 . , T = 0 . . HereAnalytic denotes the analytic solution of the problem, ILT - the ML solution, FD - the FDsolution, DiffILT - the relative error of the ML solution with respect to the analytic one,DifFD - same for the FD method. Page 27 of 36 ultilayer heat equations: application to finance
Fig. 3a shows that the accuracy of the ML method is still good with the same number of internallayers N = 40, while the error of the FD method with N = 41 , M = 40 is quite significant. The error canbe reduced by running the FD method with N = 151 , M = 150 (see Fig. 3b); however, the correspondingelapsed time is 86 ms. Thus, for the same accuracy, the ML method is 37 times faster. σ i Here σ ( x ) is a piecewise constant function defined in Eq. (78). In this case, there is no analytic solution ofthe problem , hence as a benchmark, we use an FD method, namely the same FD solver as in the previousexperiment. However, since we solve the problem in Eq. (77), an FD scheme has to be implemented forthe conservative heat equation. While such an implementation is possible, we prefer to rewrite Eq. (80)and Eq. (81) in a non-divergent form. For instance, ∂U i ∂t = Ξ i ( x ) ∂ U i ∂x + ∂ Ξ ( x ) ∂x ∂U i ∂x , (108)Ξ i ( x ) = σ i [Θ H ( x i +1 − x ) − Θ H ( x i − x )] , where Θ H ( x ) is the Heaviside theta function, (Abramowitz and Stegun, 1964) with Θ H (0) = 1. Accord-ingly, on the interval x ∈ ( y i , y i +1 ] we have ∂ Ξ ( x ) ∂x = σ i +1 [ δ ( x i − x ) − δ ( x i +1 − x )] . (109)At the point x = x i +1 this gives ∂ Ξ ( x ) /∂x = − σ i +1 δ (0). In turn, δ (0) can be numerically approximatedas δ (0) = 2 y N − y , (110)which provides the correct normalization of the Dirac delta function. Indeed, the integral over the interval[ y , y n ] of the test function equal to 1 at x = y i +1 and 0 otherwise computed by using a trapezoidal rule isequal to ( y N − y ) /
2. Therefore, we need to use Eq. (110) to provide the correct numerical normalization.
Table 2:
Parameters of the second experiment. y y N T N m
M-1.0 4.0 2.0 50 16 100In this experiment, we use parameters of the model given in Table 2, and the piecewise constantvolatility σ i , which is defined as follows σ i ( s ) = e − i/N , s ∈ ( y i , y i +1 ] , i = 1 , . . . , N. (111) Based on Eq. (104) it is possible to derive an explicit series representation of the solution. It will be published elsewhere.Page 28 of 36 ultilayer heat equations: application to finance -1 0 1 2 3 4 y V a l ue o f t he G r een F un c t i on -16-14-12-10-8-6-4-2024 R e l d i ff e r en c e , % The ML vs FD (N=51, M=100) methods
FDILTDif (a) -1 0 1 2 3 4 y V a l ue o f t he G r een F un c t i on -0.5-0.4-0.3-0.2-0.100.10.2 R e l d i ff e r en c e , % The ML vs FD (N=201, M=100) methods
FDILTDif (b)
Figure 4:
Comparison of the ML and FD solutions for a piecewise constant σ ( x ) . Here ILTdenotes the ML solution, FD - the FD solution, Dif - the relative error of the FD solutionwith respect to the ML one. The results of the test are presented in Fig. 4a. Again, the number of nodes for the ML and FDmethods is the same. The difference between the two solutions reaches 16% at the right external boundary,4% at the left external boundary, and changes in this range in between. The elapsed time is 2.2 ms forthe ML and 50 ms for the FD methods.To check the convergence of the solution, we rerun the calculation with N = 200. The results arepresented in Fig. 4b. The relative error ε drops down to be in ε ∈ [ − . , .
2] percent. The elapsed timeis 5 ms for the ML and 88 ms for the FD methods. Note that in this case, to reduce the error, we alsoneed to increase the number of layers for the ML method. This effect is explained in Section 6; it is dueto the non-smoothness of σ in this experiment.The physics meaning of the obtained results is as follows. Suppose we consider diffusion rather thanheat conduction. According to Eq. (111) the diffusion coefficient σ ( x ) is a decreasing function of x whenmoving from y to y N . Since we request continuity of the flux at the internal boundaries, the gradientof the solution increases with x when moving from left to right. The maximum of the solution, which islocated at x = x when t = 0, travels to the right when t increases. Recall that the solution is the Greenfunction of our problem. This behavior was also observed in (Lançon et al., 2001), where the authorsstudied particles trapped between two nearly parallel walls making their confinement position dependent.They not only measured a diffusion coefficient which depended on the particles’ position but also reportedand explained a new effect: a drift of the particles individual positions (so change in concentration) inthe direction of the diffusion coefficient gradient, in the absence of any external force or concentrationgradient. In the previous section, we have demonstrated that the ML method’s complexity is linear in N . The sameis true in the general case because, as mentioned at the end of Section 4, from Eq. (21) we obtain a linearsystem of equations for L ( χ − i ) , L (Ω i ) which has a block diagonal matrix with all blocks being tridiagonalmatrices. Therefore, the ML method’s complexity remains linear and approximately is O (4 mN ), so it Page 29 of 36 ultilayer heat equations: application to finance doesn’t depend on T . Hence, if 4 m is of the order of M , the ML and the FD methods have the samecomplexity. For typical values m = 12 we have 4 m = 48. Therefore, for short maturities T < x (once the values at the layers’ boundaries are found). Second, the Greeks, i.e.,derivatives of the solution, can be expressed semi-analytically by differentiating the solution with respectto x or some parameter of the model and performing numerical integration, provided that the valuesat the internal boundaries are found. For the FD method, the Greeks can be found only numerically.Moreover, to compute the Vega, a new run of the FD method is required, while for the ML method, allGreeks can be calculated in one go, as described.As far as an approximation with respect to x is concerned, the following observation holds. Using theML method, we obtain an analytical solution at every interval i, i = 1 , . . . , N . However, to do this, weneed to approximate the corresponding coefficient, e.g., σ ( x ) over layers by piecewise constant or linearfunctions. For the linear approximation, the solution’s accuracy is O ((∆ x ) ), i.e., same as for the FDmethod of the second order. Therefore, it seems that the spatial accuracies of both ML and FD methodsare the same.On the other hand, the error of both methods is also proportional to the second derivative. For theFD method, this is the second derivative of the solution; for the ML method - the second derivative ofthe coefficient, e.g., σ x,x ( x ). If the latter is smaller than the second derivative of the solution (say, theoption Gamma), then the number of layers N can be decreased while providing the same accuracy. Thisreduction provides an additional speedup of our method as compared with the FD method. This factis illustrated by our first experiment where function σ ( x ) is smooth, so even a small number of layersis sufficient to obtain a very accurate solution. In the second experiment, σ ( x ) jumps at the layer’sboundaries, and, therefore, one needs to increase the number of layers to provide the same accuracy.Note that for the FD method, the difficulties caused by sharp gradients can be alleviated by usingnonuniform grids where the nodes are condensed in the area where gradients are high. The same approachcould be applied to the construction of internal layers in the ML method.Overall, we can conclude that the new ML method proposed in this paper is significantly faster thanthe FD method, provides better accuracy, and represents the solution in a semi-analytical form. Themethod’s speed is close to that for the Radial Basis Functions (RBF) approach, (Hon and Mao, 1999;Fasshauer et al., 2004; Pettersson et al., 2008), while other properties listed above are superior to theRBF. Acknowledgments
We are grateful to Peter Carr for some fruitful discussions. Dmitry Muravey acknowledges support bythe 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.B.G. Andersen and V.V. Piterbarg.
Interest Rate Modeling . Number v. 2 in Interest Rate Modeling.Atlantic Financial Press, 2010. ISBN 9780984422111.
Page 30 of 36 ultilayer heat equations: application to finance
A. Antonov and M. Spector. General short-rate analytics.
Risk , pages 66–71, 2011.M Asvestas, A.G Sifalakis, E.P Papadopoulou, and Y.G Saridakis. Fokas method for a multi-domainlinear reaction-diffusion equation with discontinuous diffusivity.
Journal of Physics: Conference Series ,490(012143), 2014.N. Bacaer.
A short history of mathematical population dynamics , chapter 6, pages 35–39. Springer-Verlag,London, 2011. ISBN 978-0-85729-114-1.F. Black and P. Karasinski. Bond and option pricing when short rates are lognormal.
Financial AnalystsJournal , pages 52–59, 1991.D. Brigo and F. Mercurio.
Interest Rate Models – Theory and Practice with Smile, Inflation and Credit .Springer Verlag, 2nd edition, 2006.L. Capriotti and B. Stehlikova. An Effective Approximation for Zero-Coupon Bonds and Arrow-DebreuPrices in the Black-Karasinski Model.
International Journal of Theoretical and Applied Finance , 17(6):1650017, 2014.E.J. Carr and N.G. March. Semi-analytical solution of multilayer diffusion problems with time-varyingboundary conditions and general interface conditions.
Applied Mathematics and Computation , 333(15):286–303, 2018.P. Carr and A. Itkin. Geometric local variance gamma model. 27(2):7–30, 2019.P. Carr and A. Itkin. An expanded local variance gamma model. 2 2020. doi: 10.1007/s10614-020-10000-w.P. Carr and A. Itkin. Semi-closed form solutions for barrier and American options written on a time-dependent Ornstein Uhlenbeck process.
Journal of Derivatives , Fall, 2021.P. Carr and S. Nadtochiy. Local Variance Gamma and explicit calibration to option prices.
MathematicalFinance , 27(1):151–193, 2017.P. Carr, A. Itkin, and D. Muravey. Semi-closed form prices of barrier options in the time-dependent cevand cir models.
Journal of Derivatives , 28(1):26–50, 2020.M. Craddock. Fundamental solutions, transition densities and the integration of Lie symmetries.
Journalof Differential Equations , 246:2538–2560, 2009.C.J. Dias. A method of recursive images to solve transient heat diffusionin multilayer materials. 85:1075–1083.Bruno Dupire. Pricing with a smile.
Risk , 7:18–20, 1994.G. E. Fasshauer, A. Q. M. Khaliq, and D. A. Voss. Using meshfree approximation for multi-asset Americanoption problems.
J. Chinese Inst. Engrs. , 27(4):563–571, 2004.J.S. Giet, P. Vallois, and S. Wantz-Mezieres. The logistic sde.
Theory of Stochastic Processes , 20(36):28–62, 2015.Y. C. Hon and X. Z. Mao. A radial basis function method for solving options pricing model.
FinancialEngineering , 8(1):31–49, 1999.
Page 31 of 36 ultilayer heat equations: application to finance
B Horvath, A. Jacquier, and C. Turfus. Analytic option prices for the black-karasinski short rate model,2017. URL https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3253833 . SSRN: 3253833.J.C. Hull.
Options, Futures, and Other Derivatives . Prentice Hall, 8rd edition, 2011.A. Itkin.
Pricing derivatives under Lévy models . Number 12 in Pseudo-Differential Operators. Birkhauser,Basel, 1 edition, 2017.A. Itkin.
Fitting Local Volatility: Analytic and Numerical Approaches in Black-Scholes and Local VarianceGamma Models . Number 11623. World Scientific Publishing Co. Pte. Ltd., 2020.A. Itkin and A. Lipton. Filling the gaps smoothly.
Journal of Computational Sciences , 24:195–208, 2018.A. Itkin and D. Muravey. Semi-closed form prices of barrier options in the Hull-White model.
Risk ,December 2020a.A. Itkin and D. Muravey. Semi-analytic pricing of double barrier options with time-dependent barriersand rebates at hit, September 2020b. URL https://arxiv.org/abs/2009.09342 .A. Itkin, A. Lipton, and D. Muravey. From the black-karasinski to the verhulst model to accommodatethe unconventional fed’s policy, June 2020. URL https://arxiv.org/abs/2006.11976 .E. M. Kartashov. Analytical methods for solution of non-stationary heat conductance boundary problemsin domains with moving boundaries.
Izvestiya RAS, Energetika , (5):133–185, 1999.E.M. Kartashov.
Analytical Methods in the Theory of Heat Conduction in Solids . Vysshaya Shkola,Moscow, 2001.A. Kuznetsov. On the convergence of the Gaver-Stehfest algorithm.
SIAM J. Numerical Analysis , 51(6):2984–2998, 2013.P Lançon, G Batrouni, L Lobry, and N Ostrowsky. Drift without flux: Brownian walker with a space-dependent diffusion coefficient.
Europhysics Letters (EPL) , 54(1):28–34, 2001.A. Lejay. On the constructions of the skew brownian motion.
Probability Surveys , 3:413–466, 2006.J.H. Lienhard IV and J.H. Lienhard V.
A Heat Transfer Textbook . Phlogiston Press, Cambridge, MA,5th edition, 8 2019.A. Lipton.
Mathematical Methods For Foreign Exchange: A Financial Engineer’s Approach . WorldScientific, 2001.A. Lipton and M.L. de Prado. A closed-form solution for optimal mean-reverting trading strategies.
Risk ,June 2020.A. Lipton and V. Kaushansky. On the first hitting time density for a reducible diffusion process.
Quan-titative Finance, , 5, 2020a. published online.A. Lipton and V. Kaushansky. On three important problems in mathematical finance.
The Journal ofDerivatives. Special Issue , 28(2), 2020b.A. Lipton and A. Sepp. Filling the gaps.
Risk Magazine , pages 86–91, 10 2011.D. Mumford, C. Musiliand M. Nori, E. Previato, and M. Stillman.
Tata Lectures on Theta . Progress inMathematics. Birkhäuser Boston, 1983. ISBN 9780817631093.
Page 32 of 36 ultilayer heat equations: application to finance
O. A. Oleinik and E. V. Radkevich.
Second order equations with non-negative characteristic form . KluwerAcademic Publishers, 1973.Ulrika Pettersson, Elisabeth Larsson, Gunnar Marcusson, and Jonas Persson. Improved radial basisfunction methods for multi-dimensional option pricing.
J. Comput. Appl. Math. , 222(1):82–93, 2008.doi: 10.1016/j.cam.2007.10.038.A.D. Polyanin.
Handbook of linear partial differential equations for engineers and scientists . Chapman& Hall/CRC, 2002.G. Pontrelli, M. Lauricella, J.A. Ferreira, and G. Pena. Iontophoretic transdermal drug delivery: Amulti-layered approach.
Mathematical Medicine and Biology , 00:1–18, 2016.A.N. Tikhonov and A.A. Samarskii.
Equations of mathematical physics . Pergamon Press, Oxford, 1963.C. Turfus. Analytic swaption pricing in the black-karasinski model, February 2020. URL https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3253866 . SSRN: 3253866.P.F. Verhulst. Notice sur la loi que la population suit dans son accroisseement.
Correspondance mathe-matique et physique , 10:113–121, 1838.
Appendices
A Transformation of a non-divergent heat equation to a divergent form
Consider the PDE in Eq. (30) which is a divergent form of the heat equation ∂U ( t, x ) ∂t = ∂∂x (cid:18) Ξ ( x ) ∂U ( t, x ) ∂x (cid:19) . (A.1)In this Section we show how to transform it to a non-divergent form as in Eq. (28) when the externalboundaries are constant, i.e. y ( t ) = χ − ( t ) = const, y N ( t ) = χ + ( t ) = const . We start with making achange of variables x z = f ( x ) with f ( x ) = c + c Z x ( k ) dk, (A.2)where c , c are some constants. This transformation reduces Eq. (30) to ∂U ( t, z ) ∂t = σ ( z ) ∂ U ( t, z ) ∂z , (A.3) σ ( z ) = Ξ( x ( z )) x ′ ( z ) = c Ξ( x ( z )) . The Eq. (A.3) is a non-divergent form of the heat equation. The only thing which remains to be doneis finding the dependence x ( z ). Obviously, it solves the equation z = f ( x ) = c + c Z x ( k ) dk. (A.4) Page 33 of 36 ultilayer heat equations: application to finance
Given Ξ( x ), it can be solved either numerically (so this dependence can be precomputed), or in somecases analytically. As an example, assume that Ξ( x ) = e − ax , a = const = 0, and also let c = 0. Then, x = 1 a log (cid:18) ac z (cid:19) , (A.5) σ ( z ) = c ( c + az ) . Reverting these steps, we obtain the inverse transformation from a non-divergent heat equation to adivergent one.Also, the second continuity condition for Eq. (A.1) (an equality of fluxes over the boundary) is givenby Eq. (7) which, by using our notation in this Section, can be re-written asΞ i ( y i +1 ) ∂U i ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i +1 = Ξ i +1 ( y i +1 ) ∂U i +1 ∂x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i +1 , i = 1 , . . . , N − . (A.6)Using Eq. (A.3) this can be transformed to ∂U i ∂z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z = z ( y i +1 ) = ∂U i +1 ∂z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z = z ( y i +1 ) , i = 1 , . . . , N − , (A.7) z ( y i +1 ) = c + c Z y i +1 ( k ) dk = c + 1 c Z y i +1 σ ( k ) dk. This is the continuity condition for Eq. (A.3).
B Multilayer method for time-inhomogeneous coefficients and the do-main
In this section we generalize the ML method to the case σ = σ ( τ, x ). We again consider the initial-boundary problem Eq. (4) for the differential operator L i of the form Eq. (8) where now each operator L i reads L i = − ∂∂τ + σ i ( τ, x ) ∂ ∂x . (B.1)As before, we look for the solution of the problem Eq. (4) in the form Eq. (6) such that the conditionsEq. (7) still hold. Our goal is to show that under certain assumptions the problem Eq. (4) with time andspace dependent volatility σ ( τ, x ) can be reduced to the corresponding time-homogeneous problem.Suppose that for each sub-domain Ω i we can construct a map M i transforming Eq. (B.1) into PDEof the form Eq. (9). Then the solution of the transformed PDE can be represented in the form of theheat potential Eq. (15), and then transformed back by inverting the map M i . More precisely, consider acollection of maps {M i } Ni =1 acting on triplets ( τ, x, u i ( τ, x ))( τ, x, u i ( τ, x )) M i ( T i ( τ ) , X i ( τ, x ) , U i ( T i , X i )) , T i (0) = 0 , (B.2)such that the function U i ( T i , X i ) solves the following PDE with time-independent coefficients − ∂ U i ∂ T i + A ( X i ) ∂ U i ∂ X i = 0 . (B.3)Also, let us denote the inverse map as Υ i ( τ, x ), such that the following representation holds u i ( τ, x ) = Υ i ( τ, x ) U i ( X i ( τ, x ) , T i ( τ )) . (B.4) Page 34 of 36 ultilayer heat equations: application to finance
The map M i transforms the sub-domain Ω i to the sub-domain Ξ i Ξ i : h Y − i ( T i ) , Y + i ( T i ) i × R + bounded by the curves Y − i ( T i ) and Y + i ( T i ) which are defined as Y − i ( T i ) = X i ( λ i ( T i ) , y i ( λ i ( T i ))) , Y + i ( T i ) = X i ( λ i ( T i ) , y i +1 ( λ i ( T i ))) . Here λ i is the inverse map T − i , i.e. λ i ( T i ) = τ ( T i ) = T − i . Since the new time variables T i are differentfor each layer Ξ i the transformed boundaries are different as well, i.e., Y + i ( T i ) = Y − i +1 ( T i +1 ). Also, theinitial value function f ( x ) is transformed to the function F i f ( x ) M i −−→ F i ( X i ) , F i ( X i ) = f ( η i ( X i )) / Υ i ( η i ( X i )) , where η i ( x ) solves the equation η i ( x ) : X i (0 , η i ( x )) = x. Since the equations Eq. (B.3) are time-homogeneous, we can represent their solutions in the form ofEq. (15) U i ( T i , X i ) = Z T i ( Φ i ( k ) ∂G ( X i , ξ, T i − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = Y − i ( k ) + Ψ i ( k ) ∂G ( X i , ξ, T i − k ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = Y + i ( k ) ) dk. (B.5)Then making the inversion in Eq. (B.4), applying the chain rule ∂u i ∂x = U i ( T i ( τ ) , X i ( τ, x )) ∂ Υ( τ, x ) ∂x + Υ( τ, x ) ∂ X i ( τ, x ) ∂x ∂ U i ( T i ( τ ) , X ) ∂ X (cid:12)(cid:12)(cid:12)(cid:12) X = X i ( τ,x ) , and taking into account the discontinuity of the layer potentials on the boundaries, we arrive at thesystem of Volterra equations in Eq. (7).The map Eq. (B.2) can be explicitly found via two different approaches. The first is by applicationof Lie symmetry analysis. It is well known, that if Eq. (B.3) has six or four independent groups ofsymmetries, it can be reduced to the heat or Bessel PDE, see (Craddock, 2009).Another method is based on the theory of diffusion processes. Since any PDE of the form Eq. (B.1)and Eq. (B.3) can be associated with some diffusion process, say X = { X t , t ≥ } for Eq. (B.1) and Y = { Y t , t ≥ } for Eq. (B.3), the map in Eq. (B.2) can be found via reduction methods, see (Lipton andKaushansky, 2020a) and references therein. The terms T ( τ ) , X i ( τ, x ) and Υ i ( τ, x ) are interpreted as ascale, time and measure changes. C Coefficients of Eq. (97)
By using the definitions of coefficients of Eq. (97) given in Eq. (90) and Eq. (24) and tables of Laplacetransforms we find L ( η eveni ) = L ∞ X n = −∞ e − (2 nli )24 σ i t σ i √ πt = 1 σ i √ λ ∞ X n = −∞ e − √ λ | nli | σi = 1 σ i √ λ ∞ X n =1 e − √ λnliσi ! (C.1) Page 35 of 36 ultilayer heat equations: application to finance = 1 σ i √ λ e − √ λliσi − e − √ λliσi = 1 σ i √ λ coth √ λl i σ i ! , L ( η oddi ) = L ∞ X n = −∞ e − ((2 n +1) li )24 σ i t σ i √ πt = 1 σ i √ λ ∞ X n = −∞ e − √ λ | (2 n +1) li | σi = 2 e − √ λnliσi σ i √ λ ∞ X n =0 e − √ λnliσi = 1 σ i √ λ e − √ λliσi − e − √ λliσi = 1 σ i √ λ (cid:16) √ λl i σ i (cid:17) , L ( υ − ( t | x , − σ j ∞ X n = −∞ ( y j − x + 2 nl i ) e − √ λσj | y j − x +2 nl j | = − σ j ∞ X n =1 e − √ λσj ( y j − x ) − √ λσj nl j + 1 σ j ∞ X n =0 e √ λσj ( y j − x ) − √ λσj nl j = 1 σ j " e √ λσj ( y j − x ) − e − √ λσj ( y j − x ) − √ λσj l j ∞ X n =0 e − √ λσj nl j = 1 σ j e √ λσj ( y j − x ) − e − √ λσj ( y j − x ) − √ λσj l j − e − √ λσj l j = 1 σ j sinh (cid:18) ( y j − x + l ) √ λσ j (cid:19) sinh (cid:18) l j √ λσ j (cid:19) = 1 σ j sinh (cid:18) ( y j +1 − x ) √ λσ j (cid:19) sinh (cid:18) l j √ λσ j (cid:19) , L ( υ + ( t | x , − σ j ∞ X n = −∞ ( y j +1 − x + 2 nl i ) e − √ λσj | y j +1 − x +2 nl j | = − σ j ∞ X n =0 e − √ λσj ( y j +1 − x ) − √ λσj nl j + 1 σ j ∞ X n =1 e √ λσj ( y j +1 − x ) − √ λσj nl j = 1 σ j " e √ λσj ( y j +1 − x ) − √ λσj l j − e − √ λσj ( y j +1 − x ) ∞ X n =0 e − √ λσj nl j = 1 σ j e √ λσj ( y j +1 − x ) − e − √ λσj ( y j +1 − x ) − √ λσj l j − e − √ λσj l j = 1 σ j sinh (cid:18) ( y j +1 − x − l ) √ λσ j (cid:19) sinh (cid:18) l j √ λσ j (cid:19) = − σ j sinh (cid:18) x − y j ) √ λσ j (cid:19) sinh (cid:18) l j √ λσ j (cid:19) ..