The Global Economic Crisis (GEC) of 2008-2010 caused unprecedented changes in the way central banks in 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, commercial banks switched to the narrow bank model, partly by choice and partly by necessity. The Federal Reserve has 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 caused by an aging population, and low productivity growth; see (Rudebusch, 2018) . The current economic recession due to the COVID-19 pandemic forces the Fed to push the short interest rates into extremely low or outright negative territory. Given the unprecedented level of unemployment, the economic recession
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
r t = s(t) + Re zt , R = r 0 − s(0).
In other words, we modify the dynamics of the stochastic variable z t in Eq. (A.1) in the mean-reversion term by replacing z t with e zt . In Eq. (1) r t is the short interest rate, t is the time, W t is the standard Brownian motion, κ(t) > 0 is the constant speed of mean-reversion,θ(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 R = r(0) − s(0). This model is similar to the Hull-White model, but preserves positivity of r t by exponentiating the OU random variable z t . Because of that, usually practitioners add a deterministic function (shift) s(t) to the definition of r t to address possible 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 in the linear approximation on z t . Similarly, the choiceθ(t) = e θ(t) replicates the BK model at z t close the mean-reversion level θ(t). Thus, the modified BK model acquires the properties of the BK model while is 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 the underlying (for instance, the price F (r, t, S) of a Zero-coupon bond (ZCB) with maturity S) solves the following partial differential equation
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)
It is worth noting that Eq.
(2) is the stochastic Verhulst or stochastic logistic model, which are wellknown 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 stochastic Verhulst equation has the form
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 solutionr t
where S t solves the lognormal SDE
Also 1. The diffusionr t is recurrent if and only if q ≤ 0, where q = 1 2 − κ(t)θ(t) σ 2 (t) . 2. If q < 0, the diffusionr t converges in law towards the unique stationary Gamma probability distribution γ −2q, σ 2 (t)/(2κ(t) t→∞ . 3. If q > 0, 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 second part follows from Proposition 3.3 in (Giet et al., 2015) . It is interesting to note, that the condition q < 1/2 is precisely the Feller condition for the famous CIR model, (Andersen and Piterbarg, 2010) .
Thus, the stationary distribution for the Verhulst model is the Gamma distribution. It is easy to check that as compared with the mean-reversion lognormal model (the BK model) with the same parameters, the former has much fatter tails at the lower end, while the latter has the fatter tails whenr t → ∞. However, since, under the current market conditions, we are interested in modeling the lower end in the first 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 models which have the same volatility and mean-reversion rate, while the mean -reversion levelθ(t) in Eq. (1) is chosen 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 . It can be seen that r BK is always higher than r V h , which confirms our theoretical observation in above.
Next, we aim to demonstrate that the Verhulst model is also more tractable than the BK model.
In this section, we find the value of the ZCB by deriving and solving a Volterra integral equation of the second kind. We proceed with the elimination of the squared term in the drift in Eq.
(2) via the following change of variables
This yields
Another change of variables
transforms this PDE into the following one
It is worth mentioning that, by the change of variables z → −z and τ → −iτ , this PDE transforms into the 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 discussion in )
where
Now applying the Fourier transform
to both parts of Eq. (10) yields the ordinary differential equation
The solution of this problem can be written as
Now, applying the inversion formula
we obtain the following representation for u(τ, z)
Substituting the explicit representations forū(0, ω) and g(τ, ω) into Eq. (16), and taking into account that the function e −ω 2 τ is an even function, we obtain
Applying the identity ( (Gradshtein and Ryzhik, 2007) )
and changing the order of integration, we get the integral equation for u(z, τ )
This is a two-dimensional Volterra equation of the second kind, see (Lipton, 2001; Lipton and Kaushansky, 2020; for the discussion.
Here, we show that for some dependencies between the parameters of the model, the Cauchy problem Eq. (8) can be solved explicitly in terms of the Gauss hypergeometric function, (Abramowitz and Stegun, 1964) .
We start with the following change of variables
This change of variables yields
From the Black-Karasinski to the Verhulst model ...
We assume κ(t),θ(t) and σ(t) satisfy the following conditions
where C γ , C σ > 0 are some constants. Using the definitions of γ(t), φ(t), and some algebra, we can show that under these conditions, we have
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 calibrate the volatility term structure to the market, together with the constant mean reversion rate of κ and the 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 α are constants, and the normal volatility σ(t) is time-dependent. Now, applying another change of variables to Eq. (20)
we obtain the PDE with the time-homogeneous coefficients
This PDE should be solved subject to the initial and boundary conditions
Applying the Laplace transformū
to Eq. (24) and introducing µ = λ + 1/4, we obtain the following inhomogeneous ordinary differential equation
The corresponding homogeneous Eq. (27) is a Whittaker equation, which has two linearly independent solutions (the Whittaker functions) M Cγ ,µ (z) and W Cγ ,µ (z), (Abramowitz and Stegun, 1964) . A general solution of the problem Eq. (27) reads
Taking into account the asymptotic expressions for the Whittaker functions, (Abramowitz and Stegun, 1964) 
and the boundary condition in Eq. (27), we can set C 3 = C 2 = 0. Since the integrands in Eq. (28) have singularities at the points z = 0 and z = ∞ we need to check that the function u(τ, z) is regular at these points
Expression Eq. (28) can be further simplified by using the following formula, (Gradshtein and Ryzhik, 2007 
where I ν (z) is the modified Bessel function. Therefore, setting t = 1, a 1 = z, a 2 = ξ, and perceiving that Eq. (29) is symmetric with respect to a 1 and a 2 (so that the integrals on ξ in Eq. (28) are complimentary and sum up to a single integral from 0 to infinity), we obtain
Using the inversion formula for the Laplace transform, we get u(z, τ ) of the form
Here γ denotes any vertical line ℜ(λ) = γ in the complex plane such that all singularities of the integrand in Eq. (31) lie to the left of this line. Applying another identity, (Gradshtein and Ryzhik, 2007 
to the internal integral of ζ, we obtain the following representation for u(τ, z)
From the Black-Karasinski to the Verhulst model ...
By changing the variable of integration k in Eq. (32) as
This integral has poles at the points λ = λ k
since the Gamma function in the numerator of Υ(λ) turns to complex infinity when its argument is a non-positive integer or zero. The number of poles depends on the value of C α : if C α ≥ 1/2, the integrand function has no poles, if C α < 1/2 the number of the poles is K = 1 + 1 2 − C α , where ⌊x⌋ is the floor of x. The function λ + 1/4 is a multivalued function of λ, i.e., the point λ = −1/4 is a branching point. Therefore, let us construct the contour of integration in Eq. (31) as the so-called keyhole contour presented in Fig. 2 . In more detail, this contour can be described as follows. It starts with a vertical line 
Re λ
γ extending to the two big symmetric arcs Γ 1 and Γ 2 around the origin with the radius R 1 ; connecting to two horizontal parallel lines segments l 1 , l 2 at Im( λ + 1/4) = ±ω; then extending to two vertical line segments l 3 , l 4 which end points are connected to the semi-arc γ ε with the radius ε around the point Re(λ) = −1/4. Using a standard technique, we take the limit ε → 0, R 1 → ∞, so in this limit the integrals along the lines l 3 and l 4 are cancelled out. The integral along the contours Γ 1 and Γ 2 tends to zero if R 1 → ∞ due to Jordan's lemma. Hence, according to the Cauchy residue theorem, the sum of integral along the vertical line γ and two integrals along the horizontal semi-infinite lines l 1 and l 2 is equal to the sum of residuals.
Let us define the sum of residuals as R(τ, z) and the sum of integrals along the lines l 1 and l 2 with a negative 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 for the Whittaker functions M k,µ (z) and W k,µ (z), (Abramowitz and Stegun, 1964) Res
Thus, the sum of the residuals after the substitutions 1/
The integrals of ς can be computed analytically
Here U (a, b, z) is a Kummer confluent hypergeometric function, (Abramowitz and Stegun, 1964) . Using the relation between the Kummer and Whittaker functions
we finally obtain from Eq. (36), Eq. (37)
The integrals along the lines l 1 and l 2 read
Applying the connection formula between the Whittaker functions M ν,µ (z) and W ν,µ (z), and the Euler's reflection formulas, (Abramowitz and Stegun, 1964 )
we get the following expression for the sum of the integrals along the lines l 1 and l 2
Substituting this expression into Eq. (33), we obtain
Combining Eq. (39) and Eq. (38), we obtain u(τ, z) in closed-form
Since u(z, 0) = z Cα e −z/2 , we obtain a previously unknown identity
which could be verified by numerical integration. Thus, Eq. (22), Eq. (40) can be represented in the form
Also, it is known that Γ(a + bi)Γ(a − bi) ∈ R, (Cohen Jr., 1940) . Since, (Temme, 1978 )
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
where σ a , σ b , σ c are constants. We also assume that s(t) = 0, and so R = r 0 . With these assumptions, we havẽ
Since we are interested in the positive values ofθ(t), it implies σ b < 0. We present the model parameters for this test in Table 3 . We run the test for a set of maturities T ∈ [1/12, 0.3, 0.5, 1, 2, 5, 10, 20, 30, 50] years. We show the ZCB prices computed in our numerical experiment in 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 of Eq.
(2), we solve the corresponding forward equation for the density and then find the prices by integrating the 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. The Monte Carlo method uses 500,000 paths and 500 steps in time. We perform all calculations in Matlab. The results obtained by all methods coincide with high accuracy. We compare the corresponding relative errors of those results in Table 2 .
Also, we compute the ZCB prices in the BK model at small T , so that e zt ≈ 1 + z t . We choosē θ(t) BK =θ(t) M BK − 1. We show the result in Fig 3 as well. It can be seen that the BK ZCB prices agree with 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 ZCB bond prices by using the FD methods is 130 msec. For the analytics, since the only term under the integral in Eq. (42), which depends on τ is e −(ω 2 +1/4)τ , the other terms, including complex-valued Gamma and Whittaker functions, can be computed just once and then re-used. We calculate the integral by using the Simpson rule with 75 nodes; the elapsed time for getting all 10 ZCB prices is 55 msec. We note that Matlab uses Simulink to compute the Whittaker functions, so it is a bit slow. However, the performance of our method is on par with that of the forward FD solver, while the accuracy is higher. Moreover, we can further improve the accuracy of the integration by using higher-order quadratures while keeping the elapsed time similar. At the same time, for the FD method, this is problematic (but perhaps can be done by using Radial Basis Functions methods). 5238 -5.8173 -11.7015 -18.0272 -36.9372 
In this paper, we have shown that in the current market environment, it is necessary to update the classical short-rate models. We introduced a useful extension of the popular BK model, which naturally produces prolonged periods of low rates. We have derived several complementary numerical and analytical methods for its efficient solution.
The Black-Karasinski (BK) model was introduced in (Black and Karasinski, 1991) , see also (Brigo and Mercurio, 2006) for a more detailed discussion. The BK is a one-factor short-rate model of the form
Here t is the time, r t is the short interest rate, κ(t) > 0 is the constant speed of mean-reversion, θ(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), W t is the standard Brownian motion. This model is similar to the Hull-White model but preserves the positivity of r t by exponentiating the OU random variable z t . Frequently, practitioners add a deterministic function (shift) s(t) to the definition of r t to address possible negative rates and be more flexible when calibrating the term-structure of the interest rates.
By the Itô's lemma the short rater t = (r t − s(t))/R in the BK model solves the following stochastic differential equation (SDE)
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) 
and thus, conditionally on filtration F s is lognormally distributed and always stays positive. The expectation E Q [r t |F s ] and variance V[r t |F s ] can be found analytically, (Brigo and Mercurio, 2006 )
However, in the BK model, the price F (t,r) of a ZCB is not known in the closed form, since this model is not affine. Multiple good approximations have been developed in the literature using asymptotic expansions of various flavors; see, e.g., (Antonov and Spector, 2011; Capriotti and Stehlikova, 2014; Horvath et al., 2017) , and also survey in (Turfus, 2020) .
It is known that, written in terms of z, the corresponding PDE for the ZCB price F (t, z) reads, (Andersen and Piterbarg, 2010) 
This equation should be solved subject to the terminal and boundary conditions, (Andersen and Piterbarg, 2010 ) (see also discussion in )
Let us make the change of variables
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 to
Next, we apply the Fourier transform
to Eq. (B.4) to get the following problem
where δ(ω) is the Dirac delta function. The solution of Eq. (B.5) reads
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 second kind
Eq. (B.8) is a two dimensional integral Volterra equation of the second kind. Various authors have proposed efficient numerical methods for solving this type of equations. These methods include the blockby-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 (Torabi et al., 2019) and references therein. However, the complexity of the numerical methods (excluding the DTM) is at least O(N 3 ), where N is the number of computational nodes. On the other hand, N could be taken relatively small compared, e.g., with the corresponding finite-difference (FD) method, if the high order quadrature rules are used when approximating the integrals. Also, when applying all the methods mentioned above, the infinite interval should be replaced with a finite one. Another change of variables can do this, e.g., ζ → tanh(ζ). Then another class of methods can be used where the unknown function v(τ, x) is expanded into series on some basis. This basis could be 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 reduced differential transform method (RDTM), (Abazari and Kilicman, 2013) . The RDTM can be considered as an asymptotic solution of Eq. (B.8) around some time t = t 0 . It can be constructed with arbitrary precision. It is worth mentioning that the RDTM can not be directly applied to Eq. (B.8) as the kernel in Eq. (B.8) depends on τ itself. Therefore, we propose a modification of the RDTM suitable to handle this situation as well.
Next, we briefly present basic definitions of the RDTM and some theorems from (Abazari and 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 of two single-variable functions w(x, t) = f (x)g(t). The function w(t, x) can be represented as
To start with, we briefly present basic definitions of the RDTM and some theorems from (Abazari and 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 of two 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 the Poisson series of the input expression w(x, t) with respect to the variables (x, t) to order N using the variable weights W (i, j).
If w(x, t) is an analytic function in the domain of interest, then the spectrum function
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
Combining Eq. (C.2) and Eq. (C.3) one can get
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
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 0 = −∞, t 0 = 0, and thus u(k, ζ) = −β(t(k))e ζā(t(k))
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
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 0 (x) = 1, as the double integral vanishes at τ = 0, and the following properties of the RDT hold
Then from Eq. (C.6), and Eq. (C.8) we have 1
The next iteration reads
Once all the terms W k (x), k = 1, . . . , N are found, the final representation of the solution follows from the inverse formula Eq. (C.3) changed according to our modification of the RDTM v(τ, x) = N 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 τ 0 . In the latter case the method becomes almost 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).
To test the RDTM as applied to our problem, we solve Eq. (B.8) by using the modified RDTM described in Section C. Here we use the following explicit form for κ(t),θ(t), σ(t) κ(t) = κ 0 ,θ(t) =θ 0 e θ 1 t , σ(t) = σ 0 e −σ 1 t , (C.12)
where κ 0 , θ 0 , σ 0 , θ 1 , σ 1 are constants. We also assume s(t) = 0, and R = 1. With these definitions, we have a(t) = e −k(t−ν) , b(t) = e θ 1 ν − e t(κ 0 +θ 1 )−κ 0 ν κ 0 θ 0 κ 0 + θ 0 , (C.13) β(t) = 2 σ 2 0 exp 2θ 0 κ 0 + θ 1 e θ 1 t − e ν(θ 1 +κ 0 )−κ 0 t + 2t(σ 1 − κ 0 ) + 2κ 0 ν , t(τ ) = log e −2σ 1 ν − 2τ (κ 0 −2σ 1 ) We present the model parameters for this test in Table 3 . We run the test for a set of maturities T ∈ [1/12, 0.3, 0.5, 1, 2, 5]. We show the corresponding ZCB prices in Table 4 . As a benchmark, we use the solution of Eq. (B.1) obtained by using the FD method described above. The modified RDTM provides reasonable accuracy for small maturities (up to 2 years), while for T > 2, one or two terms in the expansion are insufficient to get the correct price. Therefore, for larger T , the Eq. (B.8) has to either be 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 for small T the choice ν = 0 provides slightly better results, while for T > 1 it is better to use ν = 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. In contrast, 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 FD method applied to the forward equation (the forward analog of Eq. (B.1). The elapsed time of getting the ZCB prices by solving such the equation is 40 msec while using the RDTM even with the numerical computation of all integrals in Eq. (C.10) takes 13 msec. Therefore, this method allows fast calculation of ZCB prices in the time-dependent BK model for T < 2.