From the Black-Karasinski to the Verhulst model to accommodate the unconventional Fed's policy
aa r X i v : . [ q -f i n . C P ] J un From the Black-Karasinski to the Verhulst model toaccommodate the unconventional Fed’s policy.
Andrey Itkin Alexander Lipton and Dmitry Muravey Tandon School of Engineering, New York University, 1 Metro Tech Center, 10th floor, Brooklyn NY 11201, 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
June 23, 2020 I n this paper, we argue that some of the most popular short-term interest modelshave to be revisited and modified to reflect current market conditions better. Inparticular, we propose a modification of the popular Black-Karasinski model, whichis widely used by practitioners for modeling interest rates, credit, and commodities.Our adjustment gives rise to the stochastic Verhulst model, which is well-known in thepopulation dynamics and epidemiology as a logistic model. We demonstrate that theVerhulst model’s dynamics are well suited to the current economic environment and theFed’s actions. Besides, we derive new integral equations for the zero-coupon bond pricesfor both the BK and Verhulst models. For the BK model for small maturities up to2 years, we solve the corresponding integral equation by using the reduced differentialtransform method. For the Verhulst integral equation, under some mild assumptions,we find the closed-form solution. Numerical examples show that computationally ourapproach is significantly more efficient than the standard finite difference method. Introduction
The Global Economic Crisis (GEC) of 2008-2010 caused unprecedented changes in the way central banksin general, and the mighty Fed in particular, conduct their business. The Quantitative Easing (QE)resulted in central banks embracing the fractional reserve modus operandi . At the same time, commercialbanks switched to the narrow bank model, partly by choice and partly by necessity. The Federal Reservehas used short-term interest rates as the policy tool for achieving its macroeconomic goals. As a result,short rates were close to zero for much of the past decade, reflecting the effects of QE, low inflation causedby an aging population, and low productivity growth; see (Rudebusch, 2018). The current economicrecession due to the COVID-19 pandemic forces the Fed to push the short interest rates into extremelylow or outright negative territory. Given the unprecedented level of unemployment, the economic recession rom the Black-Karasinski to the Verhulst model ... is likely to pave the way for further use of the Fed’s unconventional monetary policy, resulting in meagershort rates.Typically, short-rate interest models use stochastic drivers governed by an Ornstein-Uhlenbeck (OU)process and transform these drivers into the actual rates via suitable mappings. For instance, the Vasicek-Hull-White model uses a linear mapping, while the Black-Karasinski (BK) uses an exponential mapping.The Cox-Ingersoll-Ross model is an exception, which uses a driver governed by a Feller process. Forshort-rate models driven by OU processes, the rate spends an approximately equal amount of time belowand above its equilibrium level. This assumption was valid for decades. However, as was mentionedearlier, it is no longer adequate due to the nontraditional interventions of central banks. Once the ratebecomes low, it tends to stay low for a very long time. Under these circumstances, we have to revisitshort-rate interest models and modify them to reflect prevailing current market conditions better. Withthis motivation in mind, we consider the popular BK model, which is widely used by practitioners formodeling interest rates, credit. A similar model, known as the Schwartz one-factor model, is often used tomodel commodities. The enduring popularity of the model is because, despite some lack of tractability,it is relatively simple and guarantees non-negativity of rates. Besides, one can calibrate it to a giventerm-structure of interest rates and prices or implied volatilities of caps, floors, or European swaptions,provided that the mean-reversion level and volatility are functions of time.For the reader’s convenience, we provide some stylized facts about the BK model in Appendix A.1.As can be seen from Eq. (A.1), the short-interest rate r t in this model is lognormal and positive. (Ifnecessary, it can be made negative by using a deterministic shift s ( t ).) Initially, this positivity was oneof the significant advantages of the BK model. However, in the current environment, this feature seemsto be less useful. Another problem with the model is in the lognormality of r t . Indeed, the lognormalitymeans the CDF of the distribution is right-skewed. Therefore, the time a typical path stays in the lowerrate area is short, because the short-rate quickly moves to the mean-reversion level. Accordingly, we needto choose a low mean-reversion speed to rectify this behavior, but the qualitative behavior of r t remainsthe same. In contrast, given the above discussion, we should design an interest rate model with fat tailsat the lower end.In addition to the structural drawbacks, the BK model is not sufficiently tractable, especially whenits coefficients are time-dependent. For instance, prices of zero-coupon bonds (ZCB) and highly liquidbarrier options are not known in the close form. We have to find these prices numerically by solvingthe corresponding partial differential equations (PDEs), see Eq. (B.1), either via finite differences orasymptotically. In this paper, we present an attractive alternative, by deriving a novel integral equationfor the ZCB price; see Appendix B. The corresponding integral equation can also be solved numerically.Moreover, for small maturities (up to 2 years or so), it can be solved by using the reduced differentialtransform method; see Appendix C. Numerical examples convincingly show that in this case, our methodis more efficient computationally than the standard approach of solving PDEs via finite differences.In this paper, we propose a modification of the BK model, which organically resolves the lower endfat tail issue, and improves the model tractability. We describe the model in the next Section. We alsopresent a new integral equation for the ZCB price for our model and provide its closed-form solution in aparticular case. We demonstrate that this solution significantly accelerates the computation of the ZCBprices and provides a basis for efficient calibration. Since the BK model doesn’t support fat tails at the lower end, besides not being analytically tractable,we introduce its modified version of the form dz t = k ( t )[¯ θ ( t ) − e z t ] dt + σ ( t ) dW t , (1) Page 2 of 20 rom the Black-Karasinski to the Verhulst model ... r t = s ( t ) + Re z t , R = r − s (0) . In other words, we modify the dynamics of the stochastic variable z t in Eq. (A.1) in the mean-reversionterm by replacing z t with e z t .In Eq. (1) r t is the short interest rate, t is the time, W t is the standard Brownian motion, κ ( t ) > θ ( t ) is the mean-reversion level, σ ( t ) is the volatility, R is someconstant with the same dimensionality as r t , eg., it can be R = r (0) − s (0). This model is similar to theHull-White model, but preserves positivity of r t by exponentiating the OU random variable z t . Becauseof that, usually practitioners add a deterministic function (shift) s ( t ) to the definition of r t to addresspossible negative rates and be more flexible when calibrating the term-structure of the interest rates.It can be seen, that at small t | z t | ≪
1, and so choosing ¯ θ ( t ) = 1 + θ ( t ) replicates the BK model inthe linear approximation on z t . Similarly, the choice ¯ θ ( t ) = e θ ( t ) replicates the BK model at z t close themean-reversion level θ ( t ). Thus, the modified BK model acquires the properties of the BK model whileis a bit more tractable as this will be seen below.By the Itô’s lemma and the Feynman–Kac formula any contingent claim written on the r t as theunderlying (for instance, the 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, (2)¯ r t = r t − s ( t ) r − s (0) = e z t , ˜ θ ( t ) = ¯ θ ( t ) + σ ( t )2 κ ( t ) . This equation should be solved subject to the same terminal and boundary conditions as in Eq. (B.2) F ( S, ¯ r ) = 1 , F ( t, ¯ r ) (cid:12)(cid:12)(cid:12) ¯ r →∞ = 0 . (3)It is worth noting that Eq. (2) is the stochastic Verhulst or stochastic logistic model, which are well-known in the population dynamics and epidemiology; see, eg., (Verhulst, 1838; Bacaer, 2011; Giet et al.,2015) and references therein. In the past, several authors attempted to use this model in finance; see, eg.,(Chen, 2010; Londono and Sandoval, 2015; Halperin and Feldshteyn, 2018). In our case, the stochasticVerhulst equation has the form d ¯ r t = k ( t )¯ r t [˜ θ ( t ) − ¯ r t ] dt + σ ( t )¯ r t dW t , (4)¯ r t = [ r t − s ( t )] /R, r ( t = 0) = r . Eq. (4) can be explicitly solved (for the time-homogeneous coefficients this is done, eg., in (Giet et al.,2015), Proposition 3.3). The following Proposition holds
Proposition 1.
The Eq. (4) admits a unique positive solution ¯ r t ¯ r t = 1 κ ( t ) S tS κ (0)¯ r + R t S q dq , t ≥ , (5) where S t solves the lognormal SDE dS t = µ ( t ) S t dt + σ ( t ) S t dW t , µ ( t ) = κ ′ ( t ) κ ( t ) + κ ( t )˜ θ ( t ) . (6) Also
Page 3 of 20 rom the Black-Karasinski to the Verhulst model ...1. The diffusion ¯ r t is recurrent if and only if q ≤ , where q = − κ ( t )˜ θ ( t ) σ ( t ) .2. If q < , the diffusion ¯ r t converges in law towards the unique stationary Gamma probability distri-bution γ (cid:0) − q, σ ( t ) / (2 κ ( t ) (cid:1) t →∞ .3. If q > , the diffusion goes a.s. to zero when time goes to infinity.Proof. The proof can be obtained by applying the Itô’s lemma to Eq. (5) and using Eq. (6). The secondpart follows from Proposition 3.3 in (Giet et al., 2015). It is interesting to note, that the condition q < / r t → ∞ .However, since, under the current market conditions, we are interested in modeling the lower end in thefirst place, the Verhulst model has a distinct advantage compared with the BK model. In other words,the probability of having lower rates for the Verhulst model is much higher than for the BK model.To illustrate this in a slightly different way, we produce a set of Monte-Carlo paths for both modelswhich have the same volatility and mean-reversion rate, while the mean -reversion level ¯ θ ( t ) in Eq. (1) ischosen as ¯ θ ( t ) = 1 + θ ( t ), so the dynamics Eq. (1) corresponds to the BK dynamics in Eq. (A.1) at small z t . The results obtained by using parameters given in Section 4 are presented in Fig. 1. Figure 1:
Typical paths of the difference r BK − r V h for the short-term BK and Verhulstinterest rates as a function of time. T -4 It can be seen that r BK is always higher than r V h , which confirms our theoretical observation inabove.Next, we aim to demonstrate that the Verhulst model is also more tractable than the BK model.
Page 4 of 20 rom the Black-Karasinski to the Verhulst model ...
In this section, we find the value of the ZCB by deriving and solving a Volterra integral equation of thesecond kind. We proceed with the elimination of the squared term in the drift in Eq. (2) via the followingchange of variables F ( t, ¯ r ) = e κ ( t ) σ t ) ¯ r + R tS s ( k ) dk W ( t, ¯ r ) . (7)This yields 0 = ∂W∂t + 12 σ ( t )¯ r ∂ W∂ ¯ r + κ ( t )˜ θ ( t )¯ r ∂W∂ ¯ r + " γ ( t )¯ r − κ ( t )2 σ ( t ) ¯ r W,γ ( t ) = − R + κ ( t )˜ θ ( t ) σ ( t ) + ddt (cid:18) κ ( t ) σ ( t ) (cid:19) ,W ( S, ¯ r ) = e − κ ( S ) σ S ) ¯ r , W ( S, ¯ r ) (cid:12)(cid:12)(cid:12) ¯ r →∞ = 0 . (8)Another change of variables W ( t, ¯ r ) = u ( τ, z ) , z = ln ¯ r + Z t κ ( k )˜ θ ( k ) dk, τ = 12 Z St σ ( κ ) dκ (9)transforms this PDE into the following one ∂u∂τ = ∂ u∂z + h α ( τ ) e z + β ( τ ) e z i u, (10) α ( τ ) = γ ( t ( τ )) e − R t ( τ )0 κ ( k )˜ θ ( k ) dk , β ( τ ) = − κ ( t ( τ ))2 σ ( t ( τ )) e − R t ( τ )0 κ ( k )˜ θ ( k ) dk . It is worth mentioning that, by the change of variables z → − z and τ → − i τ , this PDE transforms intothe time-dependent Schrödinger equation with the unsteady Morse potential.This PDE in Eq. (10) should be solved subject to the initial and boundary conditions (see discussionin (Itkin and Muravey, 2020)) u (0 , z ) = A ( S ) e − e z B ( S ) , u ( τ, z ) (cid:12)(cid:12)(cid:12) z →∞ = 0 , (11)where A ( S ) = e − R S κ ( k ) dk + κ ( S ) σ S ) R S κ ( k )˜ θ ( k ) dk , B ( S ) = κ ( S ) σ ( S ) . (12)Now applying the Fourier transform¯ u ( τ, ω ) = Z ∞−∞ u ( τ, z ) e − iωz dz, (13)to both parts of Eq. (10) yields the ordinary differential equation d ¯ udτ + ω ¯ u = g ( τ, ω ) , (14) g ( τ, ω ) = Z ∞−∞ u ( τ, z ) h α ( τ ) e z + β ( τ ) e z i e − iωz dz, ¯ u (0 , ω ) = Z ∞−∞ u (0 , z ) e − iωz dz. Page 5 of 20 rom the Black-Karasinski to the Verhulst model ...
The solution of this problem can be written as¯ u ( τ, ω ) = e − ω τ ¯ u (0 , ω ) + Z τ e − ω ( τ − s ) g ( s, ω ) ds (15)Now, applying the inversion formula u ( τ, z ) = 12 π Z ∞−∞ ¯ u ( τ, ω ) e iωz dω, we obtain the following representation for u ( τ, z ) u ( τ, z ) = 12 π Z ∞−∞ (cid:20) e − ω τ ¯ u (0 , ω ) + Z τ e − ω ( τ − s ) g ( s, ω ) ds (cid:21) e iωz dω (16)Substituting the explicit representations for ¯ u (0 , ω ) and g ( τ, ω ) into Eq. (16), and taking into accountthat the function e − ω τ is an even function, we obtain u ( τ, z ) = 1 π Z ∞ (cid:26) e − ω τ Z ∞−∞ u (0 , ξ ) cos[ ω ( z − ξ )] dξ (17)+ Z τ ds e − ω ( τ − s ) Z ∞−∞ u ( τ, ξ ) h α ( s ) e ξ + β ( s ) e ξ i cos[ ω ( z − ξ )] dξ (cid:27) dω. Applying the identity ((Gradshtein and Ryzhik, 2007)) Z ∞−∞ e − βx cos( bx ) dx = r πβ exp − b β ! , and changing the order of integration, we get the integral equation for u ( z, τ ) u ( τ, z ) = 12 √ π ( Z ∞−∞ u (0 , ξ ) √ τ e − ( z − ξ )24 τ dξ + Z ∞−∞ dξ Z τ u ( s, ξ ) √ τ − s e − ( z − ξ )24( τ − s ) h α ( s ) e ξ + β ( s ) e ξ i ds ) . (18)This is a two-dimensional Volterra equation of the second kind, see (Lipton, 2001; Lipton and Kaushan-sky, 2020; Itkin and Muravey, 2020; Carr et al., 2020) for the discussion. Here, we show that for some dependencies between the parameters of the model, the Cauchy problemEq. (8) can be solved explicitly in terms of the Gauss hypergeometric function, (Abramowitz and Stegun,1964).We start with the following change of variables x = φ ( t )¯ r, φ ( t ) = e R t [ C α σ ( k ) − κ ( k )˜ θ ( k ) ] dk . (19)This change of variables yields0 = ∂W∂t + 12 σ ( t ) x ∂ W∂x + C α σ ( t ) x ∂W∂x + " γ ( t ) φ ( t ) x − κ ( t )2 σ ( t ) φ ( t ) x W (20) W ( S, x ) = e − κ ( S ) σ S ) φ ( S ) x , W ( S, x ) (cid:12)(cid:12)(cid:12) x →∞ = 0 . Page 6 of 20 rom the Black-Karasinski to the Verhulst model ...
We assume κ ( t ), ˜ θ ( t ) and σ ( t ) satisfy the following conditions C γ = γ ( t ) κ ( t ) , C σ = σ ( t ) φ ( t )2 κ ( t ) , (21)where C γ , C σ > γ ( t ) , φ ( t ), and some algebra, we can showthat under these conditions, we have C γ = C α − Rκ ( t ) , C σ = σ (0)2 κ , ˜ θ ( t ) = C α σ κ + 2 κ σ ′ ( t ) σ ( t ) , φ ( t ) = σ (0) σ ( t ) . (22)Thus, in this case the mean reversion rate κ ( t ) = κ = const , and C γ depends on C α and κ , while C σ depends on κ and σ (0). It means that we can calibrate the MBK model as follows. First, we calibratethe volatility term structure to the market, together with the constant mean reversion rate of κ andthe constant C α . Second, we determine the time-dependent mean-reversion level ˜ θ ( t ) by using Eq. (22).Thus, in this version of the model, we have three calibration parameters: two of them - κ and C α areconstants, and the normal volatility σ ( t ) is time-dependent.Now, applying another change of variables to Eq. (20) z = xC σ , τ = 12 Z St σ ( m ) dm, W ( t, x ) = z − C α e C α (1 − C α ) τ u ( τ, z ) , (23)we obtain the PDE with the time-homogeneous coefficients0 = ∂ u∂z + (cid:20) −
14 + C γ z (cid:21) u − z ∂u∂τ . (24)This PDE should be solved subject to the initial and boundary conditions u (0 , z ) = z C α e − κ ( S ) σ S ) φ ( S ) C σ z = z C α e − z/ , u ( τ, z ) (cid:12)(cid:12)(cid:12) z →∞ = 0 . (25)Applying the Laplace transform ¯ u ( λ, z ) = Z ∞ e − λτ u ( z, τ ) dτ. (26)to Eq. (24) and introducing µ = p λ + 1 /
4, we obtain the following inhomogeneous ordinary differentialequation d ¯ udz + " −
14 + C γ z + 1 / − µ z ¯ u = − u (0 , z ) z , ¯ u ( τ, z ) (cid:12)(cid:12)(cid:12) z →∞ = 0 . (27)The corresponding homogeneous Eq. (27) is a Whittaker equation, which has two linearly independentsolutions (the Whittaker functions) M C γ ,µ ( z ) and W C γ ,µ ( z ), (Abramowitz and Stegun, 1964). A generalsolution of the problem Eq. (27) reads¯ u ( λ, z ) = C M C γ ,µ ( z ) + C W C γ ,µ ( z ) (28)+ Γ (1 / − C γ + µ )Γ (1 + 2 µ ) (cid:18) W C γ ,µ ( z ) Z z u (0 , ζ ) ζ M C γ ,µ ( ζ ) dζ − M C γ ,µ ( z ) Z z u (0 , ζ ) ζ W C γ ,µ ( ζ ) dζ (cid:19) = C M C γ ,µ ( z ) + C W C γ ,µ ( z )+ Γ (1 / − C γ + µ )Γ (1 + 2 µ ) (cid:18) W C γ ,µ ( z ) Z z u (0 , ζ ) ζ M C γ ,µ ( ζ ) dζ + M C γ ,µ ( z ) Z ∞ z u (0 , ζ ) ζ W C γ ,µ ( ζ ) dζ (cid:19) C = C + Z ∞ u (0 , ζ ) ζ W C γ ,µ ( ζ ) dζ. Page 7 of 20 rom the Black-Karasinski to the Verhulst model ...
Taking into account the asymptotic expressions for the Whittaker functions, (Abramowitz and Stegun,1964) z → M κ,µ ( z ) = z µ +1 / (1 + O ( z )) , W κ,µ ( z ) = Γ(2 µ )Γ(1 / µ − κ ) z / − µ + O (cid:16) z / −ℜ ( µ ) (cid:17) ,z → ∞ : M κ,µ ( z ) ∼ Γ(1 + 2 µ )Γ(1 / µ − κ ) e z/ z − κ , W κ,µ ( z ) ∼ e − z/ z κ , and the boundary condition in Eq. (27), we can set C = C = 0. Since the integrands in Eq. (28) havesingularities at the points z = 0 and z = ∞ we need to check that the function u ( τ, z ) is regular at thesepointslim z → u ( τ, z ) = lim z → (cid:20) µz µ − / Z z u (0 , ζ ) ζ − M C γ ,µ ( ζ ) dζ + Γ (1 / − C γ + µ )Γ (1 + 2 µ ) z − ( µ +1 / Z ∞ z u (0 , ζ ) ζ − W C γ ,µ ( ζ ) dζ (cid:21) = 12 µ lim z → (cid:20) u (0 , z ) µ − / u (0 , z ) µ + 1 / (cid:21) = 1 µ − / z → u (0 , z ) = ∞ , lim z →∞ u ( τ, z ) = lim z →∞ (cid:20) Γ (1 / − C γ + µ )Γ (1 + 2 µ ) e z/ z − C γ Z z u (0 , ζ ) ζ − M C γ ,µ ( ζ ) dζ + e z/ z − C γ Z ∞ z u (0 , ζ ) ζ − W C γ ,µ ( ζ ) dζ (cid:21) = lim z →∞ " e z/ z − C γ − u (0 , z ) e z/ z − C γ (1 / − C γ /z ) − e − z/ z C γ − u (0 , z ) e − z/ z C γ ( − / C γ /z ) = lim z →∞ u (0 , z ) z − C γ z = ∞ . Expression Eq. (28) can be further simplified by using the following formula, (Gradshtein and Ryzhik,2007) Z ∞ e − ( a + a ) t cosh( k ) coth ν (cid:18) k (cid:19) I µ ( t √ a a sinh( k )) dk = Γ (1 / µ − ν ) t √ a a Γ(1 + 2 µ ) W ν,µ ( a t ) M ν,µ ( a t ) , ℜ (1 / µ − ν ) > , ℜ ( µ ) > , a > a , (29)where I ν ( z ) is the modified Bessel function. Therefore, setting t = 1 , a = z, a = ξ , and perceiving thatEq. (29) is symmetric with respect to a and a (so that the integrals on ξ in Eq. (28) are complimentaryand sum up to a single integral from 0 to infinity), we obtain¯ u ( λ, z ) = Z ∞ Z ∞ u (0 , ζ ) r zζ e − ( z + ζ ) cosh( k ) coth C γ (cid:18) k (cid:19) I µ (cid:16)p zζ sinh( k ) (cid:17) dk dζ. (30)Using the inversion formula for the Laplace transform, we get u ( z, τ ) of the form u ( τ, z ) = 12 π i Z γ e λτ Z ∞ Z ∞ u (0 , ζ ) r zζ e − ( z + ζ ) cosh( k ) coth C γ (cid:18) k (cid:19) I µ (cid:16)p zζ sinh( k ) (cid:17) dk dζ dλ. (31)Here γ denotes any vertical line ℜ ( λ ) = γ in the complex plane such that all singularities of the integrandin Eq. (31) lie to the left of this line.Applying another identity, (Gradshtein and Ryzhik, 2007) Z ∞ x ν − e − αx I µ (cid:0) β √ x (cid:1) dx = Γ( ν + µ + 1 / µ + 1) β − e β α α − ν M − ν,µ β α ! , ℜ (cid:18) µ + ν + 12 (cid:19) > , to the internal integral of ζ , we obtain the following representation for u ( τ, z ) u ( τ, z ) = e − z/ π i Z γ e λτ Z ∞ k ) ( Υ( λ ) e − z (cosh( k ) − coth C γ (cid:18) k (cid:19) (cid:18) k )2 (cid:19) − C α (32) × M − C α , √ λ +1 / (cid:18) z k ) − (cid:19) ) dk dλ, Υ( λ ) = Γ (cid:16)p λ + 1 / C α − / (cid:17) Γ (cid:16) p λ + 1 / (cid:17) . Page 8 of 20 rom the Black-Karasinski to the Verhulst model ...
By changing the variable of integration k in Eq. (32) as2 log (cid:20) tanh (cid:18) k (cid:19)(cid:21) = − ν, k sinh( k ) = − dν, cosh( k ) −
12 = 1 e ν − , we get u ( τ, z ) = e − z/ π i Z γ e λτ Υ( λ ) Z ∞ e − z eν − + νC γ M − C α , √ λ +1 / (cid:18) ze ν − (cid:19) (cid:18) e ν e ν − (cid:19) − C α dνdλ. (33)This integral has poles at the points λ = λ k r λ k + 14 = µ k = 12 − C α − k, k = 0 , . . . , (cid:22) − C α (cid:23) , (34)since the Gamma function in the numerator of Υ( λ ) turns to complex infinity when its argument is anon-positive integer or zero. The number of poles depends on the value of C α : if C α ≥ /
2, the integrandfunction has no poles, if C α < / K = 1 + j − C α k , where ⌊ x ⌋ is the floorof x . The function p λ + 1 / λ , i.e., the point λ = − / Figure 2:
Contour of integration of Eq. (31) in the complex λ plane with poles at λ , λ , . . . , λ K . This picture corresponds to C α > , and C α < / . If C α < , all poles arepositive. Re λ Im λ − / γγ ε • • • • λ λ . . . λ K Γ Γ l : p λ + 1 / ωl : p λ + 1 / − i ωl l γ extending to the two big symmetric arcs Γ and Γ around the origin with the radius R ; connecting Page 9 of 20 rom the Black-Karasinski to the Verhulst model ... to two horizontal parallel lines segments l , l at Im( p λ + 1 /
4) = ± ω ; then extending to two verticalline segments l , l which end points are connected to the semi-arc γ ε with the radius ε around the pointRe( λ ) = − /
4. Using a standard technique, we take the limit ε → , R → ∞ , so in this limit theintegrals along the lines l and l are cancelled out. The integral along the contours Γ and Γ tendsto zero if R → ∞ due to Jordan’s lemma. Hence, according to the Cauchy residue theorem, the sumof integral along the vertical line γ and two integrals along the horizontal semi-infinite lines l and l isequal to the sum of residuals.Let us define the sum of residuals as R ( τ, z ) and the sum of integrals along the lines l and l with anegative sign as I ( τ, z ). We explicitly compute them in the next section. Using the well-known expressions for the poles of the Gamma function, and the connection formula forthe Whittaker functions M k,µ ( z ) and W k,µ ( z ), (Abramowitz and Stegun, 1964)Res x = − n Γ( x ) = ( − n n ! , W a + + n, a ( z ) = ( − n Γ( a + 1 + n )Γ( a + 1) M a + + n, a ( z ) , we obtain Res λ = λ k (cid:26) e λτ Υ( λ ) M − C α , √ λ +1 / (cid:18) ze ν − (cid:19)(cid:27) = 2 µ k k ! e ( µ k − ) τ Γ (cid:16) − C α + µ k (cid:17) W − C α ,µ k (cid:18) ze ν − (cid:19) . (35)Thus, the sum of the residuals after the substitutions 1 / ( e ν − ς, µ k = − k − C α reads R ( z, τ ) = e − z/ K X k =0 µ k k ! e ( µ k − ) τ Γ (cid:16) − C α + µ k (cid:17) Z ∞ W − C α ,µ k ( zς ) e − zς (1 + ς ) C γ − C α ς − − C γ dς. (36)The integrals of ς can be computed analytically Z ∞ W − C α ,µ ( zς ) e − zς (1 + ς ) C γ − C α ς − − C γ dς = z + µ Γ + C γ (i µ ) Γ − C γ (i µ )Γ( C α − C γ ) U (cid:18)
12 + µ − C γ , µ + 1 , z (cid:19) , Γ ± y ( ω ) = Γ (cid:18) − y ± i ω (cid:19) , ℜ ( C γ ) <
12 + ℜ ( µ ) , ℜ ( C γ ) < − ℜ ( µ ) , z > . (37)Here U ( a, b, z ) is a Kummer confluent hypergeometric function, (Abramowitz and Stegun, 1964).Using the relation between the Kummer and Whittaker functions W k,µ ( x ) = e − z z / µ U (cid:18)
12 + µ − k, µ, x (cid:19) , W k,µ ( x ) = W k, − µ ( x ) , we finally obtain from Eq. (36), Eq. (37) I ( z, τ ) = K X k =0 e ( µ k − ) τ k ! 2 µ k Γ (cid:16) − C α + µ k (cid:17) Γ + C γ (i µ k ) Γ − C γ (i µ k )Γ( C α − C γ ) W C γ ,µ k ( z ) . (38) Page 10 of 20 rom the Black-Karasinski to the Verhulst model ...
The integrals along the lines l and l read Z l e λτ Υ( λ ) M − C α , √ λ +1 / (cid:18) ze ν − (cid:19) dλ = 2 e − τ/ Z ∞ ωe − ω τ Γ (i ω + C α − / ω + 1) M − C α , i ω (cid:18) ze ν − (cid:19) dω, Z l e λτ Υ( λ ) M − C α , √ λ +1 / (cid:18) ze ν − (cid:19) dλ = − e − τ/ Z ∞ ωe − ω τ Γ ( − i ω + C α − / − ω + 1) M − C α , − i ω (cid:18) ze ν − (cid:19) dω. Applying the connection formula between the Whittaker functions M ν,µ ( z ) and W ν,µ ( z ), and theEuler’s reflection formulas, (Abramowitz and Stegun, 1964) W λ,µ ( z ) = Γ( − µ )Γ(1 / − µ − λ ) M λ,µ ( z ) + Γ(2 µ )Γ(1 / µ − λ ) M λ, − µ ( z )Γ( z )Γ(1 − z ) = π sin( πz ) , Γ(1 / z )Γ(1 / − z ) = π cos( πz )we get the following expression for the sum of the integrals along the lines l and l Z l + l e λτ Υ( λ ) M − C α , √ λ +1 / (cid:18) ze ν − (cid:19) dλ = 2 e − τ/ Z ∞ ωe − ω τ ω ( Γ (i ω + C α − / ω ) M − C α , i ω (cid:18) ze ν − (cid:19) + Γ ( − i ω + C α − / − ω ) M − C α , − i ω (cid:18) ze ν − (cid:19) ) dω = 2 e − τ/ Z ∞ ωe − ω τ ω Γ +1 − C α ( ω ) Γ − − C α ( ω )Γ (2i ω ) Γ ( − ω ) W − C α , i ω (cid:18) ze ν − (cid:19) dω = − π e − τ/ Z ∞ ωe − ω τ sinh (2 πω ) Γ +1 − C α ( ω ) Γ − − C α ( ω ) W − C α , i ω (cid:18) ze ν − (cid:19) dω. Substituting this expression into Eq. (33), we obtain I ( τ, z ) = e − z/ − τ/ π Z ∞ Z ∞ ωe − ω τ sinh (2 πω ) Γ +1 − C α ( ω ) Γ − − C α ( ω ) W − C α , i ω (cid:18) ze ν − (cid:19) e − z eν − (39) × e ν (1 − C α + C γ ) ( e ν − − C α dν dω = 1 π Γ( C α − C γ ) Z ∞ ωe − ( ω +1 / τ sinh (2 πω ) Γ − C γ ( ω )Γ + C γ ( ω ) × Γ − − C α ( ω )Γ +1 − C α ( ω ) W C γ , i ω ( z ) dω. Combining Eq. (39) and Eq. (38), we obtain u ( τ, z ) in closed-form u ( τ, z ) = K X k =0 e ( µ k − ) τ k ! 2 µ k Γ (cid:16) − C α + µ k (cid:17) Γ − C γ (i µ k ) Γ + C γ (i µ k )Γ( C α − C γ ) W C γ ,µ k ( z ) (40)+ 1 π Γ( C α − C γ ) Z ∞ ωe − ( ω +1 / τ sinh (2 πω ) Γ − C γ ( ω )Γ + C γ ( ω )Γ − − C α ( ω )Γ +1 − C α ( ω ) W C γ , i ω ( z ) dω. Since u ( z,
0) = z C α e − z/ , we obtain a previously unknown identity1 π Z ∞ ω sinh (2 πω ) Γ − C γ ( ω )Γ + C γ ( ω )Γ − − C α ( ω )Γ +1 − C α ( ω ) W C γ , i ω ( z ) dω (41)= z C α e − z/ Γ( C α − C γ ) − K X k =0 k ! 2 µ k Γ (cid:16) − C α + µ k (cid:17) Γ − C γ (i µ k ) Γ + C γ (i µ k ) W C γ ,µ k ( z ) Page 11 of 20 rom the Black-Karasinski to the Verhulst model ... which could be verified by numerical integration. Thus, Eq. (22), Eq. (40) can be represented in the form F ( τ,z ) e C a ( Ca − τ + R tS s ( k ) dk = 1 + e z/ z − C α K X k =0 µ k h e ( µ k − ) τ − i k !Γ (cid:16) − C α + µ k (cid:17) Γ − C γ (i µ k ) Γ + C γ (i µ k )Γ( R/κ ) W C γ ,µ k ( z ) (42)+ 1 π Γ( R/κ ) Z ∞ ω h e − ( ω +1 / τ − i sinh (2 πω ) Γ − C γ ( ω )Γ + C γ ( ω )Γ − − C α ( ω )Γ +1 − C α ( ω ) W C γ , i ω ( z ) dω ! . Also, it is known that Γ( a + b i)Γ( a − b i) ∈ R , (Cohen Jr., 1940). Since, (Temme, 1978) W k,µ (cid:18) z (cid:19) = 2 − k √ z ∞ X n =0 n ! Γ( + 2 µ + n )Γ( + 2 µ − n ) D k − n − / ( z )(2 z ) n , D ν ( x ) = 2 ν/ U (cid:18) − ν, , x (cid:19) , we have Γ (cid:16) + 2 µ + n (cid:17) Γ (cid:16) + 2 µ − n (cid:17) = Γ (cid:18)
12 + 2 µ + n (cid:19) cos( π (2 µ − n )) π Γ (cid:18) − µ + n (cid:19) ∈ R , µ = i ω, and so W C γ , i ω ( z ) ∈ R . Thus, in Eq. (40) Γ − C γ ( ω )Γ + C γ ( ω )Γ − − C α ( ω )Γ +1 − C α ( ω ) W C γ , i ω ( z ) ∈ R . To validate our analytical solution, we compute the ZCB prices by using numerical integration of Eq. (42),where the explicit form of the parameter σ ( t ) is σ ( t ) = σ a + σ b t + σ c , (43)where σ a , σ b , σ c are constants. We also assume that s ( t ) = 0, and so R = r . With these assumptions,we have˜ θ = C α ( σ a ( σ c + t ) + σ b ) κ ( σ c + t ) − σ b k ( σ c + t )[ σ a ( σ c + t ) + σ b ] , τ = 12 (cid:20) σ a t + σ b log (cid:18) t + σ c σ c (cid:19)(cid:21) (44)Since we are interested in the positive values of ˜ θ ( t ), it implies σ b < T ∈ [1 / , . , . , , , , , , ,
50] years. We show the ZCB prices computed in our numerical experimentin Fig 3. As a benchmark, we use the numerical solution of Eq. (2) obtained by the FD method, Itkin(2017), and the solution of Eq. (1) obtained by Monte Carlo. To accelerate the FD solution, instead ofEq. (2), we solve the corresponding forward equation for the density and then find the prices by integratingthe payoff with thus found density. The FD solver runs on a non-uniform grid with 100 nodes in space z and 200 steps in time t . However, for long maturities, more nodes in space might be necessary. TheMonte Carlo method uses 500,000 paths and 500 steps in time. We perform all calculations in Matlab. Table 1:
Parameters of the test. r κ σ a σ b σ c C α Page 12 of 20 rom the Black-Karasinski to the Verhulst model ...
Figure 3:
ZCB prices computed in the test by using Analytic, FD and Monte Carlo solu-tions. T F AnalyticMBKFD-BKFD-MBKMC-MBK
The results obtained by all methods coincide with high accuracy. We compare the correspondingrelative errors of those results in Table 2.Also, we compute the ZCB prices in the BK model at small T , so that e z t ≈ z t . We choose¯ θ ( t ) BK = ¯ θ ( t ) MBK −
1. We show the result in Fig 3 as well. It can be seen that the BK ZCB prices agreewith the corresponding MBK ZCB prices for
T <
3, which is due to the fact that z is small.As far as the performance of the methods is concerned, the elapsed time for computing all 10 ZCBbond prices by using the FD methods is 130 msec. For the analytics, since the only term under the integralin Eq. (42), which depends on τ is e − ( ω +1 / τ , the other terms, including complex-valued Gamma andWhittaker functions, can be computed just once and then re-used. We calculate the integral by usingthe Simpson rule with 75 nodes; the elapsed time for getting all 10 ZCB prices is 55 msec. We note thatMatlab uses Simulink to compute the Whittaker functions, so it is a bit slow. However, the performanceof our method is on par with that of the forward FD solver, while the accuracy is higher. Moreover, wecan further improve the accuracy of the integration by using higher-order quadratures while keeping theelapsed time similar. At the same time, for the FD method, this is problematic (but perhaps can be doneby using Radial Basis Functions methods). Table 2:
Relative error in bps between Analytic, FD and MC solutions in the test.
T, yearsRel error 0.0833 0.3 0.5 1 2 5 10 20 30 50Anal - FD -0.0385 -0.0625 -0.0713 -0.0929 -0.2377 -0.2818 -0.5767 -1.4387 -2.8153 -8.7802Anal - MC -0.0844 -0.2447 -0.2803 -0.6167 -1.1802 -2.5238 -5.8173 -11.7015 -18.0272 -36.9372
Page 13 of 20 rom the Black-Karasinski to the Verhulst model ...
In this paper, we have shown that in the current market environment, it is necessary to update theclassical short-rate models. We introduced a useful extension of the popular BK model, which naturallyproduces prolonged periods of low rates. We have derived several complementary numerical and analyticalmethods for its efficient solution.
References
R. Abazari and A. Kilicman. Numerical Study of Two-Dimensional Volterra Integral Equations by RDTMand Comparison with DTM.
Abstract and Applied Analysis , (929478), 2013.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.A. Antonov and M. Spector. General short-rate analytics.
Risk , pages 66–71, 2011.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.P. Carr and A. Itkin. Semi-closed form solutions for barrier and American options written on a time-dependent Ornstein Uhlenbeck process, March 2020. Arxiv:2003.08853.P. Carr, A. Itkin, and D. Muravey. Semi-closed form prices of barrier options in the time-dependent cevand cir models, May 2020. URL https://arxiv.org/abs/2005.05459 .S. Chen. decisions Modeling the dynamics of commodity prices for investement decisions under uncer-tainty . PhD thesis, University of Waterloo, Ontario, Canada, 2010.A.C. Cohen Jr. The numerical computation of the product of conjugate imaginary gamma functions.
Ann. Math. Statist. , 11(2):213–218, 1940.J.S. Giet, P. Vallois, and S. Wantz-Mezieres. The logistic sde.
Theory of Stochastic Processes , 20(36):28–62, 2015.I.S. Gradshtein and I.M. Ryzhik.
Table of Integrals, Series, and Products . Elsevier, 2007.I. Halperin and I. Feldshteyn. Market self-learning of signals, impact and optimal trading:invisible handinference with free energy, May 2018. SSRN:3174498.
Page 14 of 20 rom the Black-Karasinski to the Verhulst model ...
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.A. Itkin.
Pricing Derivatives Under Lévy Models. Modern Finite-Difference and Pseudo-DifferentialOperators Approach. , volume 12 of
Pseudo-Differential Operators . Birkhauser, 2017.A. Itkin and D. Muravey. Semi-closed form prices of barrier options in the Hull-White model, April 2020.Arxiv:2004.09591.A. Lipton.
Mathematical Methods For Foreign Exchange: A Financial Engineer’s Approach . WorldScientific, 2001.A. Lipton and V. Kaushansky. On three important problems in mathematical finance.
The Journal ofDerivatives. Special Issue , 2020.J.A. Londono and J. Sandoval. A new logistic-type model for pricing european options.
SpringerPlus ,(4):762, December 2015.G.D. Rudebusch. A review of the fed’s unconventional monetary policy.
FRBSF Economic Letter , (27),December 2018.N.M. Temme. Uniform asymptotic expansions of confluenthypergeometric functions.
J. Inst. MathsApplies , 22:215–223, 1978.S.M. Torabi, A. Tari, and S. Shahmorad. Two-step collocation methods fortwo-dimensional volterraintegral equationsof the second kind.
Journal of Applied Analysis , 25(1):1–11, 2019.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 Stylized facts about the BK model
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-rate model of the form dz t = k ( t )[ θ ( t ) − z t ] dt + σ ( t ) dW t , r ∈ R , t ≥ , (A.1) r t = s ( t ) + Re z t , r ( t = 0) = r . Here t is the time, r t is the short interest rate, κ ( t ) > θ ( t ) isthe 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), W t is the standard Brownian motion. This model is similar to the Hull-White modelbut preserves the positivity of r t by exponentiating the OU random variable z t . Frequently, practitionersadd a deterministic function (shift) s ( t ) to the definition of r t to address possible negative rates and bemore flexible when calibrating the term-structure of the interest rates. Page 15 of 20 rom the Black-Karasinski to the Verhulst model ...
By the 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 . (A.2)This SDE can be explicitly integrated. Let 0 ≤ s ≤ t ≤ S , with S being the maturity of a ZCB. Then r t can be represented as, (Brigo and Mercurio, 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) , (A.3)and thus, conditionally on filtration F s is lognormally distributed and always stays positive. The expec-tation E Q [ r t |F s ] and variance V[ r t |F s ] can be found analytically, (Brigo and Mercurio, 2006) E Q [¯ r t |F s ] = exp (cid:20) e − k ( t − s ) log ¯ r s + k Z ts e − k ( t − u ) θ ( u ) du + 12 Z ts σ ( u ) e − k ( t − u ) du (cid:21) , (A.4)V[¯ r t |F s ] = exp (cid:20) e − k ( t − s ) log ¯ r s + 2 k Z ts e − k ( t − u ) θ ( u ) du (cid:21) h e I ( s,t ) − e I ( s,t ) i . However, in the BK model, the price F ( t, ¯ r ) of a ZCB is not known in the closed form, since this modelis not affine. Multiple good approximations have been developed in the literature using asymptoticexpansions of various flavors; see, e.g., (Antonov and Spector, 2011; Capriotti and Stehlikova, 2014;Horvath et al., 2017), and also survey in (Turfus, 2020). B An integral equation for the ZCB price in the BK model
It is known that, written in terms of z , the corresponding PDE for the ZCB price F ( t, z ) reads, (Andersenand Piterbarg, 2010) 0 = ∂F∂t + 12 σ ( t ) ∂ F∂z + κ ( t )[ θ ( t ) − z ] ∂F∂z − [ s ( t ) + Re z ] V, (B.1)This equation should be solved subject to the terminal and boundary conditions, (Andersen and Piterbarg,2010) (see also discussion in (Carr and Itkin, 2020)) F ( S, z ) = 1 , F ( t, z ) (cid:12)(cid:12)(cid:12) z →∞ = 0 . (B.2)Let us make the change of variables x = a ( t ) z + b ( t ) , τ = 12 Z tν σ ( m ) e R mν κ ( m ) dm dm, v ( τ, x ) = e − R tν s ( t ) dt F ( t, z ) , (B.3) a ( t ) = e R tν κ ( m ) dm , b ( t ) = − Z tν θ ( m ) κ ( m ) e R mν κ ( m ) dm dm, . where ν = const . As is discussed below, this constant can be chosen to simplify the final expressions.With this change Eq. (B.1) can be transformed to0 = ∂v∂τ + ∂ v∂x − β ( t ( τ )) e ¯ a ( t ( τ )) x v, (B.4) β ( t ) = 2 Rσ ( t ) e − R tν κ ( m ) dm − b ( t ) a ( t ) , ¯ a ( t ) = 1 a ( t ) v (0 , x ) = 1 , τ = 12 Z Sν σ ( m ) e R mν κ ( m ) dm dm. Page 16 of 20 rom the Black-Karasinski to the Verhulst model ...
Next, we apply the Fourier transform w ( ω, τ ) = Z ∞−∞ v ( τ, x ) e i ωx dx, to Eq. (B.4) to get the following problem w τ − ω w = β ( τ ) Z ∞−∞ v ( τ, x ) e ¯ a ( t ) x e i ωx dx, w ( τ , ω ) = 2 πδ ( ω ) , (B.5)where δ ( ω ) is the Dirac delta function. The solution of Eq. (B.5) reads w = e ω ( τ − τ ) w ( τ , ω ) + Z ττ e ω ( τ − k ) β ( k ) Z ∞−∞ F ( k, ζ ) e ¯ a ( k ) ζ − R t ( τ ) ν ( τ ) s ( t ) dt e i ωζ dζ dk. (B.6)Applying the inverse transform v ( τ, x ) = 12 π Z ∞−∞ w ( τ, ω ) e − i ωx dω yields v ( τ, x ) = 1 + 12 π Z ∞−∞ Z ττ e ω ( τ − k ) β ( k ) Z ∞−∞ v ( k, ζ ) e ¯ a ( k ) ζ e i ω ( ζ − x ) dζ dk dω = 1 − √ π Z ∞−∞ Z ττ β ( k ) e ζ ¯ a ( k ) √ τ − k e − ( x − ζ )24( τ − k ) v ( k, ζ ) dζdk, (B.7)where in the last line the change of variables k → − k, τ → − τ was made.Thus, the ZCB price solves the following two-dimensional Volterra integral equation of the secondkind v ( τ, x ) = 1 − √ π Z ∞−∞ Z ττ β ( t ( k )) e ¯ a ( t ( k )) ζ √ τ − k e − ( x − ζ )24( τ − k ) v ( k, ζ ) dζdk. (B.8) C Methods for solving Eq. (B.8)
Eq. (B.8) is a two dimensional integral Volterra equation of the second kind. Various authors haveproposed efficient numerical methods for solving this type of equations. These methods include the block-by-block method, collocation and iterated collocation methods, the differential transform method (DTM),Galerkin and spectral Galerkin methods, multi-step collocation methods, and several other, see (Torabiet al., 2019) and references therein. However, the complexity of the numerical methods (excluding theDTM) is at least O ( N ), where N is the number of computational nodes. On the other hand, N could betaken relatively small compared, e.g., with the corresponding finite-difference (FD) method, if the highorder quadrature rules are used when approximating the integrals.Also, when applying all the methods mentioned above, the infinite interval should be replaced witha finite one. Another change of variables can do this, e.g., ζ tanh( ζ ). Then another class of methodscan be used where the unknown function v ( τ, x ) is expanded into series on some basis. This basis couldbe a set of orthogonal functions, or Taylor series, etc.However, a quick estimation of the solution of Eq. (B.8) can be obtained along the lines of the reduceddifferential transform method (RDTM), (Abazari and Kilicman, 2013). The RDTM can be consideredas an asymptotic solution of Eq. (B.8) around some time t = t . It can be constructed with arbitraryprecision. It is worth mentioning that the RDTM can not be directly applied to Eq. (B.8) as the kernelin Eq. (B.8) depends on τ itself. Therefore, we propose a modification of the RDTM suitable to handlethis situation as well. Page 17 of 20 rom the Black-Karasinski to the Verhulst model ...
Next, we briefly present basic definitions of the RDTM and some theorems from (Abazari and Kilic-man, 2013) necessary to use this method for solving Eq. (B.8).Consider a function of two variables w ( t, x ), and suppose that it can be represented as a product oftwo single-variable functions w ( x, t ) = f ( x ) g ( t ). The function w ( t, x ) can be represented as w ( x, t ) = ∞ X i =0 F ( i ) x i ∞ X j =0 G ( j ) t j = ∞ X i =0 ∞ X j =0 W j ( i )( i, j ) x i t j , (C.1)where W ( i, j ) = F ( i ) G ( j ) is called the spectrum of w ( x, t ).To start with, we briefly present basic definitions of the RDTM and some theorems from (Abazariand Kilicman, 2013) necessary to use this method for solving Eq. (B.8).Consider a function of two variables w ( t, x ), and suppose that it can be represented as a product oftwo single-variable functions w ( x, t ) = f ( x ) g ( t ). The function w ( t, x ) can be represented as.If the double sum in Eq. (C.1) is truncated to the N terms in each variable, this expressions is thePoisson series of the input expression w ( x, t ) with respect to the variables ( x, t ) to order N using thevariable weights W ( i, j ).If w ( x, t ) is an analytic function in the domain of interest, then the spectrum function W k ( x ) = 1 k ! " ∂ k ∂t k w ( x, t ) t = t (C.2)is called the reduced transformed function of w ( x, t ). We use the notation where the lowercase w ( x, t )denotes the original function while the uppercase W k ( x ) stands for the reduced transformed function.The differential inverse transform of W k ( x ) is defined as w ( x, t ) = ∞ X k =0 W k ( x )( t − t ) k . (C.3)Combining Eq. (C.2) and Eq. (C.3) one can get w ( x, t ) = ∞ X k =0 k ! " ∂ k ∂t k w ( x, t ) t = t ( t − t ) k . (C.4)To proceed, we need the following fragment of Theorem 7 in (Abazari and Kilicman, 2013) Theorem 1.
Assume that U k ( x ) , H k ( x ) and W k ( x ) are the reduced differential transforms of the functions u ( x, t ) , h ( x, t ) and w ( x, t ) , respectively. If w ( x, t ) = Z tt Z xx h ( y, z ) u ( y, z ) dydz, (C.5) then W k ( x ) = 1 k Z xx k − X ν =0 H ν ( y ) U k − ν − ( y ) ! dy, k = 1 , , . . . . (C.6) Proof.
See (Abazari and Kilicman, 2013).Now one can observe that Eq. (B.8) actually has the form of Eq. (C.5) with x = ∞ , x = −∞ , t = 0,and thus u ( k, ζ ) = − β ( t ( k )) e ζ ¯ a ( t ( k )) p π ( τ − k ) e − ( x − ζ )24( τ − k ) . (C.7) Page 18 of 20 rom the Black-Karasinski to the Verhulst model ...
Our modification of the RDTM consists in eliminating the definition in Eq. (C.2) for the function u ( k, ζ ). Then, based on the definition of the reduced differential transform in Eq. (C.2), we get U ( s, ζ ) = − β ( t ( s )) e ζ ¯ a ( t ( s )) p π ( τ − s ) e − ( x − ζ )24( τ − s ) , (C.8) U ( s, ζ ) = − Z s β ( t ( k )) e ζ ¯ a ( t ( k )) p π ( s − k ) e − ( x − ζ )24( s − k ) dk, and so on.Let us denote the reduced differential transform of v ( τ, x ) as W k ( x ). From Eq. (B.8) it follows that W ( x ) = 1, as the double integral vanishes at τ = 0, and the following properties of the RDT hold w k ( x ) = u k ( x ) ± v k ( x ) W k ( x ) = U k ( x ) ± V k ( x ) ,w ( x, t ) = A = const, W ( x ) = 1 , W k ( x ) = 0 , k > , Then from Eq. (C.6), and Eq. (C.8) we have W ( x ) = Z τ Z ∞−∞ W ( s, ζ ) U ( s, ζ ) dζds = Z τ e ¯ a ( s )(¯ a ( s )( τ − s )+ x ) ds. (C.9)The next iteration reads W ( x ) = Z τ Z ∞−∞ [ W ( ζ ) U ( ζ ) + W ( ζ ) U ( ζ )] dζ ds = I + I , (C.10) I = Z τ Z s β ( k ) β ( s ) e a ( k )(2 a ( s )( τ − s )+ x )+ a ( k ) ( τ − k )+ a ( s )( a ( s )( τ − s )+ x ) dk ds,I = − Z τ Z s β ( t ( k )) e ζ ¯ a ( t ( k )) p π ( s − k ) e − ( x − ζ )24( s − k ) dk ds, etc.Once all the terms W k ( x ) , k = 1 , . . . , N are found, the final representation of the solution followsfrom the inverse formula Eq. (C.3) changed according to our modification of the RDTM v ( τ, x ) = N X k =0 W k ( x ) . (C.11)The time-integrals in Eq. (C.9), Eq. (C.10) can be computed either numerically, or analytically if functions¯ a ( τ ) , β ( τ ) could be expanded into series on τ around some τ . In the latter case the method becomesalmost identical to the original RDTM. When doing so, one has to remember that derivatives of a ( τ ) , β ( τ )are the derivatives on τ while the definitions of these functions in Eq. (B.3) are given in terms of t = t ( τ ).The latter map is also given in Eq. (B.3). C.1 Numerical example
To test the RDTM as applied to our problem, we solve Eq. (B.8) by using the modified RDTM describedin Section C. Here we use the following explicit form for κ ( t ) , ˜ θ ( t ) , σ ( t ) κ ( t ) = κ , ˜ θ ( t ) = ˜ θ e θ t , σ ( t ) = σ e − σ t , (C.12) This expression now contains an integral in time since we eliminated the Taylor series expansion in Eq. (C.2).Page 19 of 20 rom the Black-Karasinski to the Verhulst model ... where κ , θ , σ , θ , σ are constants. We also assume s ( t ) = 0, and R = 1. With these definitions, wehave a ( t ) = e − k ( t − ν ) , b ( t ) = (cid:16) e θ ν − e t ( κ + θ ) − κ ν (cid:17) κ θ κ + θ , (C.13) β ( t ) = 2 σ exp (cid:20) θ κ + θ (cid:16) e θ t − e ν ( θ + κ ) − κ t (cid:17) + 2 t ( σ − κ ) + 2 κ ν (cid:21) ,t ( τ ) = log (cid:16) e − σ ν − τ ( κ − σ ) σ (cid:17) + κ νκ − σ . Table 3:
Parameters of the test. r κ θ σ θ σ Table 4:
Prices of ZCB bonds with different maturities computed by using the modifiedRDTM with 1 and 2 terms, and the FD difference method.
ZCB price, $T 0.0833 0.3 0.5 1.0 2.0 5.0FD
RDTM-1
RDTM-2
Difference with the FD method, %T 0.0833 0.3 0.5 1.0 2.0 5.0FD
RDTM-1
RDTM-2
We present the model parameters for this test in Table 3. We run the test for a set of maturities T ∈ [1 / , . , . , , , T >
2, one or two terms inthe expansion are insufficient to get the correct price. Therefore, for larger T , the Eq. (B.8) has to eitherbe solved numerically or more terms should be taken in the RDTM.Note, that there are at least two choices for ν in Eq. (B.3): ν = 0 and ν = T . We found that forsmall T the choice ν = 0 provides slightly better results, while for T > ν = T .Also, note that this method, in some sense, is similar to that in (Capriotti and Stehlikova, 2014).However, we solve an integral equation instead of a PDE in (Capriotti and Stehlikova, 2014). Besides,there is a difference in parametrization, since we assume that all the parameters are time-dependent. Incontrast, in (Capriotti and Stehlikova, 2014), all model parameters are constant.As far as the performance of the RDTM is concerned, we compared it with the performance of the FDmethod applied to the forward equation (the forward analog of Eq. (B.1). The elapsed time of gettingthe ZCB prices by solving such the equation is 40 msec while using the RDTM even with the numericalcomputation of all integrals in Eq. (C.10) takes 13 msec. Therefore, this method allows fast calculationof ZCB prices in the time-dependent BK model for T <2.