2. Actuarial valuation

We denote by $B_t$ the dollar amount in a savings account at time $t$ , compounding interest at an instantaneous rate $r_t$ and having an initial amount of one dollar deposited at time $0$ . We have, for $t\ge 0$ ,

(2.1) \begin{align} B_t = B_0 \mathrm{exp} \bigg \{ \int _0^t r_s\, ds \bigg \} . \end{align}

In actuarial valuations of liabilities consisting of cash flows $CF_i$ at times $t_i$ , $i\in \{1,2,\ldots, n\}$ , the present $V_t$ value at time $t$ , where $t\lt t_1\lt \ldots \lt t_n$ , is computed as

(2.2) \begin{align} V_t = \sum _{i=1}^n E\bigg ( \frac{B_t}{B_{t_i}} CF_i \bigg | \mathcal{A}_t\bigg ), \end{align}

where $E({\cdot} |\mathcal{A}_t)$ denotes the conditional expectation at time $t$ with respect to the real-world probability measure $P$ . As per Delbaen & Schachermayer (Reference Delbaen and Schachermayer1994), classical risk-neutral assumptions are that in a market which prohibits arbitrage, in the sense of no free lunches with vanishing risk (NFLVR), discounted price processes are martingales with respect to an equivalent risk-neutral probability measure $Q$ , and we have the risk neutral valuation

(2.3) \begin{align} V_t^* = \sum _{i=1}^n E^Q\bigg ( \frac{B_t}{B_{t_i}} CF_i \bigg | \mathcal{A}_t\bigg ). \end{align}

When the random variables $\frac{B_t}{B_{t_i}} CF_i$ are independent of the corresponding Radon–Nikodym derivatives $\frac{dQ}{dP}|_{\mathcal{A}_{t_i}}$ for $i\in \{1,2,\ldots, n\}$ and these form a martingale, (2.2) arises from (2.3) and vice versa. The existence of an equivalent risk-neutral probability measure $Q$ is often assumed for most models in the literature, but it may not exist for more general models. In such situations, we can perform valuations under the real-world probability measure using Platen’s benchmark approach, described in Platen (Reference Platen2002) and Platen & Heath (Reference Platen and Heath2006), where the actuarial pricing formula and the risk-neutral pricing formula emerge as particular cases.

When the cash flows are independent of the short rate, (2.2) can be written as

(2.4) \begin{align} V_t = \sum _{i=1}^n E(CF_i|\mathcal{A}_t) \times G_{t_i}(t), \end{align}

where $G_{t_i}(t)$ denotes the actuarial value of a $t_i$ -maturity ZCB at time $t$ , computed as

(2.5) \begin{align} G_{t_i}(t) = E\bigg (\mathrm{exp} \bigg \{ - \int _t^{t_i} r_s\, ds \bigg \} \bigg | \mathcal{A}_t\bigg ) . \end{align}

For certain cash flows, the actuarial valuation formula (2.4) is often otherwise written as

(2.6) \begin{align} V_t = \sum _{i=1}^n CF_i \times (1+ R)^{-(t_i-t)}, \end{align}

where $R$ is the internal rate of return, which can be inferred from the valuations of $G_{t_i}(t)$ supplied via the $3/2$ short rate model. In practice, a valuation actuary commences with an assumed value for $R$ , called the discount rate, and then computes the present value of a series of certain cashflows using (2.6). The assumed value of $R$ is based typically on the government bond yield having commensurate term to maturity with those of the cashflows, as well as having the same currency of denomination.

Our aim in this paper is to determine the level of the long-term discount rate $R$ when using the $3/2$ model of the short rate. In forthcoming sections, we determine formulae for $G_{t_i}(t)$ , and, for appropriately chosen parameters, determine the level of the asymptotic bond yield under the $3/2$ model. In our penultimate section, we determine formulae for bond option prices and provide estimates of the asymptotic implied volatility.

3. Explicit formula

Setting $R_t = 1/r_t$ in (1.3) we obtain the SDE

(3.1) \begin{align} dR_t = (\sigma ^2 - q - p_t\, R_t ) dt - \sigma \sqrt{R_t} dZ_t, \end{align}

which shows that $R_t = 1/r_t$ follows a square root process. This fact makes the derivation of the explicit solution to the SDE of the $3/2$ model in the following lemma, and its transition density function straightforward. The conditions in (1.2) ensure that the coefficients in (3.1) satisfy the Feller condition, which ensures that $R_t$ is strictly positive and nonexplosive. Under such conditions, the explicit solution permits exact simulations of the short rate to be performed for the extended $3/2$ short rate model.

Lemma 1. A solution to the SDE (1.3), for $p_t\gt 0$ , and $q$ and $\sigma$ satisfying (1.2), is

(3.2) \begin{align} r_t = \mathrm{exp} \bigg(\int _0^t p_u\, du\bigg ) / \sum _{i=1}^\nu (\lambda ^{(i)} + Z_{\varphi _t}^{(i)} )^2, \end{align}

where $\nu$ is an integer such that

(3.3) \begin{align} \nu = \frac{4 (\sigma ^2 - q)}{ \sigma ^2 }, \end{align}

$\lambda ^{(1)},\ldots, \lambda ^{(\nu )}\in{\mathrm{I\!R}}$ are chosen such that $r_0 = 1/\sum _{i=1}^\nu (\lambda ^{(i)})^2$ , where

(3.4) \begin{align} \varphi _t = \varphi _0 + \frac{1}{4} \sigma ^2 \int _0^t \mathrm{exp}\,\! \bigg (\int _0^s p_u\, du\bigg)\, ds \end{align}

and $Z^{(1)},\ldots, Z^{(\nu )}$ are independent Brownian motions.

Proof. From (3.1) $R_t = 1/r_t$ follows a square root process and, therefore, using the explicit solution for the CIR model given in Section 8.7 of Platen & Heath (Reference Platen and Heath2006), the explicit formula for $R_t$ is

(3.5) \begin{align} R_t = \mathrm{exp} \big (-\int _0^t p_u\, du\big ) \sum _{i=1}^\nu (\lambda ^{(i)} + Z_{\varphi _t}^{(i)} )^2, \end{align}

where $\lambda ^{(1)},\ldots, \lambda ^{(\nu )}$ are chosen such that $R_0 = \sum _{i=1}^\nu (\lambda ^{(i)})^2$ and where $\varphi _t$ is given in (3.4) and $\nu$ is given in (3.3). The formula for $r_t$ is simply the reciprocal of (3.5).

Remark 1. The solution to the SDE (1.3) can be extended to the case where $\nu$ is a positive noninteger, as indicated in Chapter 3 of Baldeaux & Platen (Reference Baldeaux and Platen2013). Let $\nu$ be a nonnegative real number satisfying (3.3), and let $X$ be a strong solution to the SDE $dX_t = \nu \,dt + 2\sqrt{X_t}\, dZ_t$ , with initial condition $X_0 = 1/r_0$ . Then 3.2 emerges as a particular case of the general solution

(3.6) \begin{align} r_t = \mathrm{exp} \big (\int _0^t p_u\, du\big ) / X_{\varphi _t} \end{align}

of (1.3) for nonnegative $\nu$ .

4. Transition density

The transition density function of the extended $3/2$ short rate model is easily deduced from the transition density of the CIR short rate model as specified in Section 4.4 of Platen & Heath (Reference Platen and Heath2006) because the SDE of the reciprocal of a $3/2$ -process short rate obeys the CIR SDE (4.3), as demonstrated in (3.1).

Lemma 2. The transition density function of the short rate process in (1.3) is

(4.1) \begin{align} p_r(s,r_s, t, r_t ) &= \frac{r_t^{-2}}{2(\varphi _t -\varphi _s ) \mathrm{exp} \big (-\int _0^t p_u\, du\big )} \left ( \frac{r_t^{-1} \mathrm{exp} \big (\int _0^t p_u\, du\big )}{r_s^{-1} \mathrm{exp} \big ( \int _0^s p_u\, du\big )} \right )^{\frac{1}{2} (\frac{\nu }{2} - 1)}\end{align}

\begin{align} &\times \mathrm{exp} \left ( -\frac{1}{2} \frac{r_s^{-1} \mathrm{exp} \big (\int _0^s p_u\, du\big ) + r_t^{-1} \mathrm{exp} \big (\int _0^t p_u\, du\big )}{ ( \varphi _t-\varphi _s)} \right ) \nonumber\\ \notag &\times I_{\frac{\nu }{2}-1} \left ( \frac{\sqrt{r_s^{-1} r_t^{-1}\mathrm{exp} \big (\int _0^s p_u\, du\big )\mathrm{exp} \big (\int _0^t p_u\, du\big )}}{ (\varphi _t-\varphi _s)} \right ) \nonumber \end{align}

where $\varphi _t$ is given in (3.4), $\nu$ is given in (3.3) and

(4.2) \begin{align} I_\nu (x) = \sum _{i=0}^\infty \frac{1}{i! \Gamma (i+\nu + 1) } \left ( \frac{x}{2} \right )^{2i+\nu } \end{align}

is the power series expansion of the modified Bessel function of the first kind.

Proof. It is known, for example from Corollary 5.2 of Fergusson (Reference Fergusson2019), that the transition density function of the short rate process

(4.3) \begin{align} dR_t = \kappa _t (\bar{r}_t - R_t)\, dt + \sigma _t \sqrt{R_t}\,{d}Z_t, \end{align}

where $\bar{r}_t$ , $\kappa _t$ and $\sigma _t$ are positively valued functions satisfying

(4.4) \begin{align} \bar{r}_t\kappa _t \ge \frac{1}{2}\sigma _t^2, \end{align}

is

(4.5) \begin{align} p_R(s,R_s, t, R_t ) =\,& \frac{1}{2(\varphi _t -\varphi _s ) \mathrm{exp} \big (-\int _0^t \kappa _u \, du \big )} \left ( \frac{R_t \mathrm{exp} \big (\int _0^t \kappa _u \, du \big )}{R_s \mathrm{exp} \big (\int _0^s \kappa _u \, du \big )} \right )^{\frac{1}{2} (\frac{\nu }{2} - 1)} \\ \notag &\times \mathrm{exp} \left ( -\frac{1}{2} \frac{R_s \mathrm{exp} \big (\int _0^s \kappa _u \, du \big ) + R_t \mathrm{exp} \big (\int _0^t \kappa _u \, du \big )}{ ( \varphi _t-\varphi _s)} \right ) \\ \notag &\times I_{\frac{\nu }{2}-1} \left ( \frac{\sqrt{R_s R_t\mathrm{exp} \big (\int _0^s \kappa _u \, du \big )\mathrm{exp} \big (\int _0^t \kappa _u \, du \big )}}{ (\varphi _t-\varphi _s)} \right ), \end{align}

where, in this case,

(4.6) \begin{align} \varphi _t = \varphi _0 + \frac{1}{4} \int _0^t \sigma _s^2 \mathrm{exp} \big ( \int _0^s \kappa _u\, du\big ) \, ds, \end{align}

$\nu = 4\kappa _t\bar{r}_t/\sigma _t^2$ and $I_\nu$ is given in (4.2). Note that, because the formula for $I_{\nu }$ given in (4.2) can be analytically continued to values of $\nu \in \mathbb{R}-\{ \ldots, -2,-1 \}$ , the transition density formula (4.5) is valid for noninteger values of $\nu \gt 2$ .

Comparing coefficients of (3.1) with those in (4.3) gives the correspondences

(4.7) \begin{align}{\bar{r}_t} & \leftrightarrow \frac{\sigma ^2 - q}{p_t} \\ \notag \kappa _t & \leftrightarrow p_t \\ \notag \sigma _t & \leftrightarrow \sigma \end{align}

and making these substitutions in (4.5) gives

(4.8) \begin{align} p_R(s,R_s, t, R_t ) =\,& \frac{1}{2(\varphi _t -\varphi _s ) \mathrm{exp} \big (-\int _0^t p_u \, du \big )} \left ( \frac{R_t \mathrm{exp} \big (\int _0^t p_u \, du \big )}{R_s \mathrm{exp} \big (\int _0^s p_u \, du \big )} \right )^{\frac{1}{2} (\frac{\nu }{2} - 1)} \\ \notag &\times \mathrm{exp} \left ( -\frac{1}{2} \frac{R_s \mathrm{exp} \big (\int _0^s p_u \, du \big ) + R_t \mathrm{exp} \big (\int _0^t p_u \, du \big )}{ ( \varphi _t-\varphi _s)} \right ) \\ \notag &\times I_{\frac{\nu }{2}-1} \left ( \frac{\sqrt{R_s R_t\mathrm{exp} \big (\int _0^s p_u \, du \big )\mathrm{exp} \big (\int _0^t p_u \, du \big )}}{ (\varphi _t-\varphi _s)} \right ) . \end{align}

Finally, the transition density function of $r_t = 1/R_t$ is

(4.9) \begin{align} p_r(s,r_s, t, r_t ) &= \frac{d}{d r_t} \int _0^{r_t} p_r(s,r_s, t, w ) dw \\ \notag &= \frac{d}{d r_t} \int _{R_t}^\infty p_R(s,R_s, t, x ) dx \\ \notag &= -\frac{d R_t}{d r_t} \times p_R(s, R_s, t, R_t ) \\ \notag &= r_t^{-2} p_R(s, r_s^{-1}, t, r_t^{-1} ) \end{align}

and using (4.8) gives the result.

The following corollary explicitly gives the distribution of the reciprocal of the extended $3/2$ short rate.

Corollary 1. For $t\gt s$ and for the short rate process $r$ defined in (1.3), the conditional random variable

(4.10) \begin{align} \frac{\mathrm{exp} \big (\int _0^t p_u\, du\big ) }{ r_t(\varphi _t-\varphi _s)} \end{align}

given $r_s$ has a non-central chi-squared distribution with $\nu$ degrees of freedom and non-centrality parameter $\lambda = \mathrm{exp} \big (\int _0^s p_u\, du\big )/(r_s(\varphi _t - \varphi _s ))$ , namely

(4.11) \begin{align} \frac{\mathrm{exp} \big (\int _0^t p_u\, du\big ) }{ r_t(\varphi _t-\varphi _s) } \sim \chi ^2_{\nu, \mathrm{exp} \big (\int _0^s p_s\, du\big )/( r_s(\varphi _t-\varphi _s))}, \end{align}

given $r_s$ , where $\varphi _t$ is given in (3.4) and $\nu$ is given in (3.3).

Proof. We know from Corollary 5.4 of Fergusson (Reference Fergusson2019) that for $t\gt s$ and for the CIR short rate process $r$ in (4.3), the conditional random variable

(4.12) \begin{align} \frac{\mathrm{exp} (\int _0^t \kappa _u \, du ) }{ \varphi _t-\varphi _s} r_t \end{align}

given $r_s$ has a non-central chi-squared distribution with $\nu = 4\kappa _t{\bar{r}_t}/\sigma _t^2$ degrees of freedom and non-centrality parameter $\lambda = r_s \mathrm{exp} ( \int _0^s \kappa _u \, du)/(\varphi _t - \varphi _s )$ , namely

(4.13) \begin{align} \frac{\mathrm{exp} (\int _0^t \kappa _u \, du) }{ \varphi _t-\varphi _s } r_t \sim \chi ^2_{\nu, r_s \mathrm{exp} (\int _0^s \kappa _u \, du )/( \varphi _t-\varphi _s)}, \end{align}

where $\varphi _t$ is given in (4.6). Adapting this to the extended $3/2$ model gives the result.

We augment Robert Shiller’s annual series of one-year deposit rates from 1871 to 2011, given in Shiller (Reference Shiller2022), with the series IR3TCD01USM156N of certificates of deposit rates 2011–2023, given in Federal Reserve Bank of St. Louis (2023). Using maximum likelihood estimation, as described in Fergusson (Reference Fergusson2021), the fitted parameters of the $3/2$ model are given by

(4.14) \begin{align} \bar{p} &= 0.296974\,(0.074103 )\\ \notag q &= 6.188698\,(2.450063 ) \\ \notag \sigma &= 4.930868\,(0.383273 ), \end{align}

where the standard errors are shown in brackets.

In Fig. 2 the transition density of the $3/2$ model is plotted, illustrating the positivity of $r_t$ .

Figure 2 $3/2$ transition density of US cash rates.

The mean reverting level of $R_t = 1/r_t$ in (3.1) is given by

(4.15) \begin{align} \frac{\sigma ^2 - q}{\bar{p}} . \end{align}

However, the inverse of this level is not the mean reverting level of $r_t$ , as the following lemma shows.

Lemma 3. Let $r_t$ obey (1.1). Then the limiting distribution of $\frac{1}{4}\sigma ^2 r_t/\bar{p}$ as $t\to \infty$ is an inverse chi-squared distribution with $\nu$ degrees of freedom and the mean of $r_t$ as $t\to \infty$ is

(4.16) \begin{align} \lim _{t\to \infty } E(r_t|\mathcal{A}_0) = \frac{2\bar{p}}{\sigma ^2 - 2q} . \end{align}

Proof. The Kolmogorov forward equation for the transition density function $p_r(s,x,t,y)$ , given in (4.1), is

(4.17) \begin{align} \frac{\partial }{\partial t}p_r(s,x,t,y) +\frac{\partial }{\partial y}\big \{a(t,y)\,p_r(s,x,t,y)\big \} -\frac{1}{2} \frac{\partial ^2}{\partial y^2}\big \{b^2(t,y)\,p_r(s,x,t,y)\big \} = 0, \end{align}

where $a(t,y)=a(y) = \bar{p} y+qy^2$ is the drift term and $b(t,y)=b(y)=\sigma y^{3/2}$ is the diffusion term in (1.1). Since we are seeking the stationary density $\bar{p}_r$ , which is independent of the starting and ending times and starting value, the stationary, or time-independent, forward equation for $\bar{p}_r$ is

(4.18) \begin{align} \frac{d}{d y}\big \{a(y)\,\bar{p}_r(y)\big \} -\frac{1}{2} \frac{d^2}{d y^2}\big \{b^2(y)\,\bar{p}_r(y)\big \} = 0. \end{align}

Letting

(4.19) \begin{align} H(y)=\big \{a(y)\,\bar{p}_r(y)\big \} -\frac{1}{2} \frac{d}{d y}\big \{b^2(y)\,\bar{p}_r(y)\big \}, \end{align}

we see from (4.18) that $H(y)$ must be a constant, and in the light of $\bar{p}_r(y)\to 0$ and $\frac{d}{d y}\bar{p}_r(y)\to 0$ as $y\to \infty$ , this constant must be zero. The equation $H(y)=0$ boils down to

(4.20) \begin{align} \big \{ a(y) - \frac{1}{2}\frac{d}{d y}b^2(y)\big \} \, \bar{p}_r(y) -\frac{1}{2} b^2(y) \frac{d}{d y}\bar{p}_r(y) =0, \end{align}

and solving this linear ordinary differential equation gives

(4.21) \begin{align} \bar{p}_r(y) = \frac{C}{b^2(y) } \mathrm{exp} \bigg \{ 2\int \frac{a(y)}{b^2(y)}\, dy \bigg \} = \frac{C}{\sigma ^2 } y ^{2q/\sigma ^2 - 3} \mathrm{exp} \bigg \{ - \frac{2\bar{p}}{\sigma ^2 y} \bigg \} . \end{align}

upon making the substitutions $a(y) = \bar{p}y+qy^2$ and $b^2(y)=\sigma ^2 y^3$ . We recognize this as the probability density of the reciprocal gamma distribution, having parameters $\alpha =2-2q/\sigma ^2$ and $\lambda = 2\bar{p}/\sigma ^2$ . We know that the mean of the reciprocal gamma distribution is $\lambda / (\alpha - 1)$ , and the result of the lemma is deduced to be

(4.22) \begin{align} \frac{4\bar{p}/\sigma ^2}{\nu -2} = \frac{4\bar{p}}{2\sigma ^2 - 4q} . \end{align}

Figure 3 The fitted reverting levels under the $3/2$ model and the extended $3/2$ models ( $N=1,2,3$ ).

In Fig. 3, the graph of the short rate is shown along with the implied reverting level, the mean, of $0.0498$ of the $3/2$ model.

As an illustration of the extended $3/2$ model, we fit three versions of (1.3) to our augmented data set, with

\begin{align*}p_t = \bar {p} + \sum _{j=1}^N \left ( a_j\cos 2\pi t/\lambda _j + b_j\sin 2\pi t/\lambda _j\right ), \end{align*}

for each $N\in \{1,2,3\}$ . The maximum likelihood parameter estimates are shown in Table 1, where the standard errors are shown in brackets. The values of the log-likelihood, AIC and BIC for each of the models indicate the improvement in fit of the model as $N$ is increased. The likelihood ratio test (LRT) is employed to compare the fit of each model relative to the extended $3/2$ model with $N=3$ . The LRT statistic of a model is computed as $\lambda _{LR} = 2(\ell ^* - \ell )$ , where $\ell$ is the model’s log-likelihood and $\ell ^*$ is the log-likelihood of the $N=3$ model. It is approximately distributed as chi-squared with the degrees of freedom equal to the different in the number of estimated parameters. Despite the models with $N\lt 3$ being rejected at the 99% level of significance, there is potential for overfitting as $N$ is increased.

Table 1. Parameter estimates and model diagnostics in respect of the extended $3/2$ models

Shown in Fig. 4 are the probability–probability plots of the four models corresponding to $N=0,1,2,3$ . The extent to which the P-P plot adheres to the line $y=x$ demonstrates the model’s ability to fit the distribution of the one-year-ahead short rate. It is evident in Fig. 4 that the extended $3/2$ models, i.e., models with $N\in \{1,2,3\}$ , provide a better fit to the data.

For comparison with the $3/2$ model, Fig. 3 also depicts the implied reverting levels of the extended $3/2$ models with $N=1,2,3$ , calculated as $2\bar{p}/(\sigma ^2-2q)=0.0404$ , $0.0327$ , and $0.0327$ , respectively.

5. ZCB price

We calculate the $\bar{T}$ -maturity ZCB price $G_{\bar{T}} (t)$ . The following lemma gives a formula for the Laplace transform of the random variable $\int _t^T r_s ds$ when $r_t$ obeys (1.3). A variation of this result involving the joint Fourier-Laplace transform is proven in Theorem 3 in Carr & Sun (Reference Carr and Sun2007) and proofs involving Lie symmetries of partial differential equations are provided in Craddock & Lennox (Reference Craddock and Lennox2009) and Section 5.5 of Baldeaux & Platen (Reference Baldeaux and Platen2013). A proof of the following lemma which uses the Lie algebra of differential operators is given in Lo (Reference Lo2013) and requires some astute observations on eigenfunctions of differential operators and a clever modification of the form of the payoff function. In contrast, the original proof of the ZCB price formula which we supply here, makes use of a transformation of the two variables $x$ and $t$ in the PDE that reduces it to an ordinary differential equation (ODE), which then is readily identified as the confluent hypergeometric differential equation. This is similar to the method of proof of Theorem 3 in Carr & Sun (Reference Carr and Sun2007), where this transformation or form of the solution is assumed and subsequently verified.

Lemma 4. Let $r_t$ obey (1.3). Then the conditional Laplace transform of $\int _t^T r_s ds$ satisfies

(5.1) \begin{align}{\mathrm{E}} \bigg ( \mathrm{exp} (-u \int _t^T r_s ds) \bigg |{\mathcal{A}}_t \bigg ) = \frac{\Gamma (\gamma _u - \alpha _u )}{\Gamma (\gamma _u )} \bigg ( \frac{2}{\sigma ^2 y(t,r_t)} \bigg )^{\alpha _u} M(\alpha _u, \gamma _u, \frac{-2}{\sigma ^2 y(t,r_t)} ) \end{align}

where

(5.2) \begin{align} y(t,r_t) &= r_t \int _t^T \mathrm{exp} \big (\int _t^\xi p_u\, du \big ) \, d\xi \\ \notag \alpha _u &= -\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2} \bigg ) + \sqrt{\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2 } \bigg )^2 + \frac{2u}{\sigma ^2}} \\ \notag \gamma _u &= 2\bigg ( \alpha _u +1 - \frac{q}{\sigma ^2} \bigg ) \end{align}

and $M(\alpha, \gamma, z)$ is the confluent hypergeometric function, given by

(5.3) \begin{align} M(\alpha, \gamma, z) = \frac{\Gamma (\gamma )}{\Gamma (\alpha ) \Gamma (\gamma - \alpha )} \int _0^1 \mathrm{exp} (zu) u^{\alpha - 1} (1-u)^{\gamma - \alpha - 1} du \end{align}

or the alternative expression

(5.4) \begin{align} M(\alpha, \gamma, z) = \sum _{n=0}^\infty \frac{(\alpha )_n}{(\gamma )_n} \frac{z^n}{n!}, \end{align}

where $(\alpha )_n := \alpha (\alpha +1)\ldots (\alpha + n-1)$ is the Pochhammer notation.

Furthermore, we have

(5.5) \begin{align} E\bigg (\mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} \bigg | r_{\bar{T}}, \,{\mathcal{A}}_t\bigg ) = \frac{I_{\sqrt{8/\sigma ^2+(\nu /2-1)^2}} (\sqrt{xy}/\tau )}{I_{\nu /2-1} (\sqrt{xy}/\tau )}, \end{align}

where $ y = \mathrm{exp} (\int _0^{\bar{T}} p_u\,du) /r_{\bar{T}}$ and $ x = \mathrm{exp} (\int _0^t p_u\,du) /r_{t}$ .

Figure 4 Graphs of quantiles of the CDFs of one-year-ahead short rates under the $3/2$ model ( $N=0$ ) and the extended $3/2$ models ( $N=1,2,3$ ).

Proof. From Section 9.7 of Platen & Heath (Reference Platen and Heath2006), the Feynman-Kac Theorem states that if

(5.6) \begin{align} u(x,t) &={\mathrm{E}} \bigg ( \int _t^T \mathrm{exp} \bigg \{-\int _t^s V(X_\tau, \tau )\,d\tau \bigg \} f(X_s, s)\,ds \\ \notag &+ \mathrm{exp} \bigg \{-\int _t^T V(X_\tau, \tau )\,d\tau \bigg \}\psi (X_T) \bigg | X_t = x\bigg ), \end{align}

where $X=\{X_t:t\ge 0\}$ is a stochastic process satisfying

(5.7) \begin{align} dX_t = \mu (X_t,t)dt+\sigma (X_t,t)dW_t \end{align}

then $u(x,t)$ satisfies the partial differential equation (PDE)

(5.8) \begin{align} \frac{\partial u}{\partial t} + \mu (x,t) \frac{\partial u}{\partial x} + \frac{1}{2}\sigma (x,t)^2\frac{\partial ^2 u}{\partial x^2} - V(x,t)\,u(x,t) - f(x,t)=0, \end{align}

with the boundary condition

(5.9) \begin{align} u(x,T)=\psi (x) . \end{align}

Applying this theorem to the problem of computing the Laplace transform

(5.10) \begin{align} u(r_t,t) ={\mathrm{E}} \bigg ( \mathrm{exp} \bigg \{-s\int _t^T r_\tau \,d\tau \bigg \} \bigg |{\mathcal{A}}_t \bigg ), \end{align}

where $r=\{r_t\,:\,t\ge 0\}$ is the stochastic process given in (1.1), we must solve the PDE (5.8). We have

(5.11) \begin{align} \mu (x,t) &= p_t x+qx^2\\ \notag \sigma (x,t) &= \sigma x^{3/2} \\ \notag \psi (x) &= 1 \\ \notag V(x,t) &= s\, x\\ \notag f(x,t) &= 0 \end{align}

and (5.8) becomes

(5.12) \begin{align} \frac{\partial u}{\partial t} + (p_t x+qx^2) \frac{\partial u}{\partial x} + \frac{1}{2}\sigma ^2x^3\frac{\partial ^2 u}{\partial x^2} - s\,x\,u(x,t) =0. \end{align}

The boundary conditions for the solution to (5.12) are

(5.13) \begin{align} u(x,T) &= 1 \\ \notag u(0,t) &= 1 \\ \notag \lim _{x\to \infty }u(x,t) &= 0. \end{align}

Time-Dependent Coefficient of $r_t$ in the Drift Term

We now consider the general case $p_t$ being time-dependent, where $p_t\gt 0$ for all $t\gt 0$ . We combine $x$ and $t$ into the single variable

(5.14) \begin{align} y = F(t ) x, \end{align}

where $F(t)$ is some function to be determined, and we set

(5.15) \begin{align} w(y) := u(x,t ) . \end{align}

A consequence of this transformation is that $w(0) = u(0,t )$ and $w(F(T)x) = u(x,T)$ . So that the boundary conditions in (5.13) are satisfied we insist that $F(T)=0$ and that the boundary conditions for $w(y)$ be

(5.16) \begin{align} w(0) &= 1 \\ \notag \lim _{y\to \infty }w(y) &= 0. \end{align}

We have

(5.17) \begin{align} \frac{\partial u}{\partial t} &= w^\prime (y)\,F^\prime (t)\, x\\ \notag \frac{\partial u}{\partial x} &= w^\prime (y)\,F(t)\\ \notag \frac{\partial ^2 u}{\partial x^2} &= w^{\prime \prime } (y)\,F(t)^2 = w^{\prime \prime } (y) \frac{y^2}{x^2} \end{align}

and making these substitutions in (5.12), we obtain

(5.18) \begin{align} w^\prime (y)\, (F^\prime (t)\, x+p_t\, F(t) \, x+q\,x\,y) + \frac{1}{2}\sigma ^2 x y^2 w^{\prime \prime }(y) - s\,x \, w(y) =0. \end{align}

Now if we can choose $F(t)$ such that

(5.19) \begin{align} F^\prime (t)+p_t\, F(t) = -1, \end{align}

then (5.18) simplifies to

(5.20) \begin{align} w^\prime (y)\, (-1+q\,y) + \frac{1}{2}\sigma ^2 y^2 w^{\prime \prime }(y) - s\, w(y) =0. \end{align}

We see immediately that

(5.21) \begin{align} F(t) = \int _t^T \mathrm{exp} \bigg (\int _t^\xi p_\zeta \, d\zeta \bigg )\, d\xi \end{align}

solves (5.19), with the boundary condition $F(T)=0$ .

Analytic Solution for $w(y)$ at Point at Infinity

Instead of seeking an analytic solution for $w$ at $y=0$ and then trying to satisfy the boundary condition at infinity, we seek an analytic solution at the point at infinity and then attempt fulfilling the boundary condition at $y=0$ . Hence, let us write $w(y)=z^\alpha h(z)$ where $z= 1/y$ , $h$ is an analytic function of $z$ around $z=0$ , and $\alpha$ is some positive real number to be determined. Positivity of $\alpha$ ensures that at $z=0$ , and hence as $y\to \infty$ , $w(y) = 0$ , which is the second boundary condition in (5.16). The first and second derivatives of $w$ are

(5.22) \begin{align} w^\prime (y) &= - \alpha \, z^{\alpha +1}h(z)- z^{\alpha +2}h^\prime (z) \\ \notag w^{\prime \prime } (y) &= \alpha (\alpha +1) z^{\alpha +2}h(z)+(2\alpha +2) z^{\alpha +3}h^\prime (z) + z^{\alpha +4}h^{\prime \prime } (z) \end{align}

and substituting into (5.20) gives

(5.23) \begin{align} &\frac{1}{2}\sigma ^2 z^{\alpha +2} h^{\prime \prime }(z) + \bigg (\frac{1}{2}(2\alpha +2)\sigma ^2z^{\alpha +1}-z^{\alpha +1}(-z+q )\bigg ) h^\prime (z) \\ \notag &+\bigg ( - s\, z^\alpha - \alpha \, z^\alpha (-z+q )+\frac{1}{2}\alpha (\alpha +1)\sigma ^2z^\alpha \bigg ) h(z) =0. \end{align}

Dividing both sides by $z^{\alpha }$ gives

(5.24) \begin{align} &\frac{1}{2}\sigma ^2 z^{2} h^{\prime \prime }(z) + \bigg (\frac{1}{2}(2\alpha +2)\sigma ^2z- z(-z+q)\bigg ) h^\prime (z) \\ \notag &+\bigg ( - s + \alpha z-\alpha q+\frac{1}{2}\alpha (\alpha +1)\sigma ^2 \bigg ) h(z) =0. \end{align}

Choosing $\alpha$ such that

(5.25) \begin{align} - s -\alpha q+\frac{1}{2}\alpha (\alpha +1)\sigma ^2 =0 \end{align}

and dividing both sides of (5.24) by $\frac{1}{2}\sigma ^2 z$ gives

(5.26) \begin{align} z h^{\prime \prime }(z) + \bigg ((2\alpha +2)+\frac{2}{\sigma ^2 }z-\frac{2}{\sigma ^2}q\bigg ) h^\prime (z) +\frac{2\alpha }{\sigma ^2 } h(z) =0. \end{align}

Making the change of variable $\zeta = -2z/\sigma ^2$ , letting $m(\zeta ) = h(z)$ and letting

(5.27) \begin{align} \gamma = 2\alpha + 2 - \frac{2q}{\sigma ^2}, \end{align}

we obtain

\begin{align*} \zeta \, m^{\prime \prime }(\zeta ) + (\gamma - z) m^\prime (\zeta ) - \alpha \, m(\zeta ) = 0, \end{align*}

which we recognize as the confluent hypergeometric differential equation, whose two independent solutions are given by $M(\alpha, \gamma, \zeta )$ for $\alpha$ satisfying the quadratic equation (5.25); see Chapter 4 of Watson (Reference Watson1966) for details. Hence, $w(y)$ is some linear combination of

\begin{align*} \left ( 1/ y \right )^{\alpha _+} M(\alpha _+,\gamma _+,-2/(\sigma ^2 y))\quad \mathrm{and}\quad \left ( 1/ y \right )^{\alpha _-} M(\alpha _-,\gamma _-,-2/(\sigma ^2 y)), \end{align*}

where $\alpha _+$ and $\alpha _-$ are the upper and lower roots of (5.25), and $\gamma _+$ and $\gamma _-$ are the corresponding values of $\gamma$ given by (5.27).

Satisfying the Boundary Conditions in (5.16)

We seek a choice of $w(y)$ such that the first boundary condition in (5.16) is satisfied. First, choosing

\begin{align*} w(y) = h_0\left ( 1/ y \right )^\alpha M(\alpha, \gamma, -2/(\sigma ^2 y)), \end{align*}

where $h_0$ is some constant and $\alpha$ equals the upper root of (5.25), that is,

(5.28) \begin{align} \alpha = \bigg(\frac{q}{\sigma ^2}-\frac{1}{2}\bigg)+\sqrt{\bigg (\frac{q}{\sigma ^2}-\frac{1}{2}\bigg )^2+\frac{2s}{\sigma ^2}}, \end{align}

ensures that $\alpha \gt 0$ , and, therefore, as $y\to \infty$ we have $w(y) \to 0$ . Second, using the integral form of the confluent hypergeometric function

(5.29) \begin{align} M(\alpha, \gamma, \zeta ) = \frac{\Gamma (\gamma )}{\Gamma (\alpha ) \Gamma (\gamma - \alpha )} \int _0^1 \mathrm{exp} (\zeta u) u^{\alpha - 1} (1-u)^{\gamma - \alpha - 1} du \end{align}

we have, using the transformation in the integral $t=-\zeta u$ , that as $\zeta \to -\infty$ ,

(5.30) \begin{align} M(\alpha, \gamma, \zeta ) &= \frac{\Gamma (\gamma )}{\Gamma (\alpha ) \Gamma (\gamma - \alpha )} \int _0^{-\zeta } \mathrm{exp} (-t ) (-t/\zeta )^{\alpha - 1} (1+t/\zeta )^{\gamma - \alpha - 1} \frac{1}{-\zeta }dt \end{align}

(5.31) \begin{align} &= \frac{\Gamma (\gamma )}{\Gamma (\alpha ) \Gamma (\gamma - \alpha )} \bigg ( -\frac{1}{\zeta }\bigg )^\alpha \int _0^{-\zeta } \mathrm{exp} (-t ) t^{\alpha - 1} dt \bigg ( 1 + O(1/\zeta )\bigg )\\[-10pt]\nonumber \end{align}

(5.32) \begin{align} &\sim \frac{\Gamma (\gamma )}{\Gamma (\alpha ) \Gamma (\gamma - \alpha )} \bigg ( -\frac{1}{\zeta }\bigg )^\alpha \Gamma (\alpha ) \\[10pt]\nonumber . \end{align}

Thus, as $z\to \infty$ , we have that $y\to 0^+$ and

(5.33) \begin{align} M(\alpha, \gamma, -2z/\sigma ^2) \sim \frac{\Gamma (\gamma )}{\Gamma (\gamma -\alpha )} (2z/\sigma ^2)^{-\alpha }. \end{align}

Hence, as $z\to \infty$ ,

(5.34) \begin{align} w(y) &=z^\alpha \, h_0\, M(\alpha, \gamma, -\frac{2z}{\sigma ^2}) \\ \notag &\sim z^\alpha \, h_0\,\frac{\Gamma (\gamma )}{\Gamma (\gamma -\alpha )} \bigg (\frac{2z}{\sigma ^2}\bigg )^{-\alpha } \\ \notag &= h_0\,\frac{\Gamma (\gamma )}{\Gamma (\gamma -\alpha )} \bigg (\frac{2}{\sigma ^2}\bigg )^{-\alpha } \end{align}

and setting

(5.35) \begin{align} h_0 = \frac{\Gamma (\gamma -\alpha )}{\Gamma (\gamma )} \bigg (\frac{2}{\sigma ^2}\bigg )^{\alpha } \end{align}

makes $w(y) = h_0\, z^\alpha M(\alpha, \gamma, -2z/\sigma ^2 )\to 1$ as $z\to \infty$ .

Solution for $u(x,t)$

Thus, from (5.15) we have

(5.36) \begin{align} u(x,t) = \frac{\Gamma (\gamma -\alpha )}{\Gamma (\gamma )} \bigg ( \frac{2}{\sigma ^2 y}\bigg )^\alpha M\big (\alpha, \gamma, -\frac{2}{\sigma ^2 y}\big ), \end{align}

where $y$ is given by (5.14), $\alpha$ is given by (5.28), and $\gamma$ is given by (5.27).

Finally, we now evaluate

\begin{align*} E\bigg (\mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} \bigg | r_{\bar{T}}, \,{\mathcal{A}}_t\bigg ) . \end{align*}

Proposition 6.2.1 on p. 164 of Baldeaux & Platen (Reference Baldeaux and Platen2013) gives

(5.37) \begin{align} E\bigg (\mathrm{exp} \bigg \{-\frac{b^2}{2}\int _0^\tau \frac{ds}{\tilde{X}_s} \bigg \} \bigg | \tilde{X}_\tau = y,\, \tilde{X}_0=x \bigg ) = \frac{I_{\sqrt{b^2+(\nu /2-1)^2}} (\sqrt{xy}/\tau )}{I_{\nu /2-1} (\sqrt{xy}/\tau )}, \end{align}

where $\tilde{X} = \{\tilde{X}_t |t\ge 0\}$ is given by the SDE

\begin{align*} d\tilde{X}_t = \nu \, dt + 2\sqrt{\tilde{X}_t}\, dW_t, \end{align*}

with $\tilde{X}_0=x\gt 0$ and $\nu \ge 2$ . Observe that making the change of variable $\xi = \varphi ^{-1}(s+ \varphi _t )$ in the integral in (5.37), noting that ${\bar{T}}=\varphi ^{-1}(\tau + \varphi _t)$ , and employing (3.6), with $\tilde{X}_s = X_{s+\varphi _t}$ , we have

\begin{align*} \int _0^\tau \frac{ds}{\tilde{X}_s} = \int _{t}^{\bar{T}} \frac{\frac{1}{4}\sigma ^2 \mathrm{exp} \left (\int _0^\xi p_u\, du\right )}{\tilde{X}_{\varphi _\xi - \varphi _t} } \, d\xi = \int _{t}^{\bar{T}} \frac{\frac{1}{4}\sigma ^2 \mathrm{exp} \left (\int _0^\xi p_u\, du\right )}{X_{\varphi _\xi } } \, d\xi = \frac{1}{4}\sigma ^2\int _{t}^{\bar{T}} r_\xi \, d\xi . \end{align*}

Thus,

\begin{align*} &E\bigg (\mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} \bigg | r_{\bar{T}}, \,{\mathcal{A}}_t\bigg ) \\ &= E\bigg (\mathrm{exp} \bigg \{-\frac{4}{\sigma ^2}\int _0^\tau \frac{ds}{\tilde{X}_s} \bigg \} \bigg | \tilde{X}_\tau = y,\, \tilde{X}_0=x \bigg ) \\ &= \frac{I_{\sqrt{8/\sigma ^2+(\nu /2-1)^2}} (\sqrt{xy}/\tau )}{I_{\nu /2-1} (\sqrt{xy}/\tau )}, \end{align*}

where $ y = \mathrm{exp} (\int _0^{\bar{T}} p_u\,du) /r_{\bar{T}}$ and $ x = \mathrm{exp} (\int _0^t p_u\,du) /r_{t}$ .

Equivalently, we have that the moment generating function of the random variable $L=\int _t^T r_s ds$ is

(5.38) \begin{align} MGF_L(-u) = \frac{\Gamma (\gamma - \alpha )}{\Gamma (\gamma )} \bigg ( \frac{2}{\sigma ^2 y(t,r_t)} \bigg )^\alpha M(\alpha, \gamma, \frac{-2}{\sigma ^2 y(t,r_t)} ), \end{align}

where the variables are defined as in the above lemma.

The following corollary follows straightforwardly from the above lemma.

Corollary 2. Let $r_t$ obey the SDE (1.3). For time $t \in [ 0, \bar{T} ]$ , the $\bar{T}$ -maturity ZCB price is

(5.39) \begin{align} G_{\bar{T}} (t) = \frac{\Gamma (\gamma _1 - \alpha _1)}{\Gamma (\gamma _1 )} \bigg ( \frac{2}{\sigma ^2 y(t,r_t)} \bigg )^{\alpha _1} M(\alpha _1, \gamma _1, \frac{-2}{\sigma ^2 y(t,r_t)} ), \end{align}

where

(5.40) \begin{align} y(t,r_t) &= r_t\, \int _t^T \mathrm{exp} \big (\int _t^\xi p_u\, du \big ) \, d\xi \\ \notag \alpha _1 &= -\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2} \bigg ) + \sqrt{\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2 } \bigg )^2 + \frac{2}{\sigma ^2}} \\ \notag \gamma _1 &= 2\bigg ( \alpha _1 +1 - \frac{q}{\sigma ^2} \bigg ) . \end{align}

6. Bond yields and forward rates

We employ the notation $h_{\bar{T}}(t)$ for the $\bar{T}$ -maturity ZCB yield, as given by

(6.1) \begin{align} h_{\bar{T}}(t) = -\frac{1}{\bar{T}-t} \log G_{\bar{T}}(t). \end{align}

The asymptotic value of the ZCB yield is given in the following corollary.

Theorem 1. For the short rate $r_t$ obeying SDE (1.3), with $p_t$ satisfying Condition 1, the long ZCB yield is

(6.2) \begin{align} h_\infty (t) = \alpha _1 \bar{p}, \end{align}

where

(6.3) \begin{align} \alpha _1 = -\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2} \bigg ) + \sqrt{\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2 } \bigg )^2 + \frac{2}{\sigma ^2}} \end{align}

and $\bar{p}$ is given in Condition 1.

Proof. We write the logarithm of (5.39) as

(6.4) \begin{align} \log G_{\bar{T}} (t) = \log \Gamma (\gamma _1 - \alpha _1) - \log \Gamma (\gamma _1 ) + \alpha _1 \log z + \log M(\alpha _1, \gamma _1, -z ), \end{align}

where $z = \frac{2}{\sigma ^2 y(t,r_t)}$ and $y(t,r_t) = r_t \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi p_u\, du \}\, d\xi$ . Let $\epsilon \gt 0$ be arbitrarily small and, using Condition1, let $\bar{T}_0$ be such that for all $\bar{T}\gt \bar{T}_0$ ,

(6.5) \begin{align} \sup _{\xi \in [t,\bar{T}]} \bigg | \int _t^\xi (p_u -\bar{p})\, du \bigg | \lt \epsilon \, (\bar{T} - t) . \end{align}

Writing

(6.6) \begin{align} \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi p_u\, du \}\, d\xi = \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi \bar{p}\, du \}\mathrm{exp} \{\int _t^\xi (p_u-\bar{p})\, du \}\, d\xi \end{align}

and making use of (6.5) gives

(6.7) \begin{align} &\frac{1}{\bar{p}} \big (\mathrm{exp} \{ \bar{p} (\bar{T} - 1)\} - 1\big )\, \mathrm{exp} \{ -\epsilon (\bar{T} - t) \} \\ \notag &= \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi \bar{p}\, du \}\, d\xi \, \mathrm{exp} \{ -\epsilon (\bar{T} - t) \} \\ \notag &\lt \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi p_u\, du \}\, d\xi \\ \notag &\lt \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi \bar{p}\, du \}\, d\xi \, \mathrm{exp} \{ \epsilon (\bar{T} - t) \} \\ \notag &= \frac{1}{\bar{p}} \big (\mathrm{exp} \{ \bar{p} (\bar{T} - 1)\} - 1\big )\, \mathrm{exp} \{ \epsilon (\bar{T} - t) \} . \end{align}

Taking logarithms and dividing by $\bar{T} - t$ gives

(6.8) \begin{align} &\frac{-\log \bar{p} + \log \big (\mathrm{exp} \{ \bar{p} (\bar{T} - t)\} - 1\big ) -\epsilon (\bar{T} - t) }{\bar{T}-t} \\ \notag &\lt \frac{\log \int _t^{\bar{T}} \mathrm{exp} \{\int _t^\xi p_u\, du \}\, d\xi }{\bar{T} - t}\\ \notag & \lt \frac{-\log \bar{p} + \log \big (\mathrm{exp} \{ \bar{p} (\bar{T} - t)\} - 1\big ) +\epsilon (\bar{T} - t) }{\bar{T}-t} . \end{align}

Thus,

(6.9) \begin{align} \bar{p}-\epsilon \le \lim _{{\bar{T}}\to \infty } \frac{\log y(t,r_t)}{{\bar{T}}-t} \le \bar{p}+\epsilon . \end{align}

Since $\epsilon$ is arbitrary, we have

\begin{align*} \lim _{{\bar {T}}\to \infty } \frac {\log y(t,r_t)}{{\bar {T}}-t} = \bar {p} .\end{align*}

The long ZCB yield is given by the formula

(6.10) \begin{align} h_\infty (t) &= -\lim _{{\bar{T}}\to \infty } \frac{1}{{\bar{T}}-t} \log G_{\bar{T}} (t) \\ \notag &= \lim _{{\bar{T}}\to \infty } \frac{-1}{{\bar{T}}-t} ( \log \Gamma (\gamma _1 - \alpha _1) - \log \Gamma (\gamma _1 ) + \alpha _1 \log z + \log M(\alpha _1, \gamma _1, -z ) ) \\ \notag &= \lim _{{\bar{T}}\to \infty } \frac{-\alpha _1 \log z}{{\bar{T}}-t} = \lim _{{\bar{T}}\to \infty } \frac{\alpha _1 \log y(t,r_t)}{{\bar{T}}-t} \\ \notag &= \alpha _1 \bar{p} \end{align}

as required.

In Fig. 5, the continuously compounded yield curve as of January 1871 is plotted corresponding to the start time in Fig. 3. For the $3/2$ model with constant $p_t$ , we have an inverted yield curve, which portends an economic recession because decreasing forward rates indicate expectations of low inflation and low economic growth; see for example Harvey (Reference Harvey1991). In Fig. 6, the corresponding curves for the extended $3/2$ model with $N=2$ are plotted, where we also observe a flat yield curve for about 40 years, beyond which the yield curve inverts for 35 years, after which the yields increase for maturities up to 110 years.

Figure 5 The zero-coupon yield curve and instantaneous forward rate curve under the $3/2$ model as at 1871.

Figure 6 The zero-coupon yield curve and instantaneous forward rate curve under the extended $3/2$ model, with $N=2$ , as at 1871.

Employing the notation $g_{\bar{T}} (t)$ for the instantaneous $\bar{T}$ -forward rate, based at time $t$ and given by

(6.11) \begin{align} g_{\bar{T}} (t) = -\frac{\partial }{\partial \bar{T}} \log G_{\bar{T}}(t), \end{align}

we calculate $g_{\bar{T}} (t)$ in the following lemma.

Lemma 5. For time $t \in [ 0, \bar{T} ]$ the $\bar{T}$ -forward rate at time $t$ is computed to be

(6.12) \begin{align} g_{\bar{T}} (t) = \alpha _1 \frac{1}{\int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi } \bigg ( 1 - \frac{z }{\gamma _1} \frac{M(\alpha _1 +1, \gamma _1 +1, -z)}{M(\alpha _1, \gamma _1, -z) } \bigg ) \end{align}

where

(6.13) \begin{align} z = \frac{2}{\sigma ^2 y(t,r_t)} = \frac{2}{\sigma ^2 r_t} \frac{1}{\int _t^{\bar{T}} \mathrm{exp} \{ \int _t^\xi p_u\, du \} \, d\xi }. \end{align}

Proof. We write the logarithm of (5.39) as

(6.14) \begin{align} \log G_{\bar{T}} (t) = \log \Gamma (\gamma _1 - \alpha _1) - \log \Gamma (\gamma _1 ) + \alpha _1 \log z + \log M(\alpha _1, \gamma _1, -z ) \end{align}

where $z = \frac{2}{\sigma ^2 y(t,r_t)}$ and $y(t,r_t) ={r_t}\int _t^{\bar{T}} \mathrm{exp} \{ \int _t^\xi p_u\, du \} \, d\xi$ . We note that

(6.15) \begin{align} \frac{\partial z}{\partial{\bar{T}}} = \frac{-z}{\int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi } . \end{align}

Also, we have straightforwardly

(6.16) \begin{align} \frac{\partial M(\alpha _1, \gamma _1, z)}{\partial z} = \frac{\alpha _1}{\gamma _1} M(\alpha _1 +1, \gamma _1 +1, z) . \end{align}

From (6.11) and the above relations, we have

(6.17) \begin{align} g_{\bar{T}} (t) &= - \frac{\partial }{\partial{\bar{T}} } \log G_{\bar{T}} (t) \\ \notag &= z \frac{1}{\int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi } \bigg ( \frac{\alpha _1}{z} - \frac{\alpha _1}{\gamma _1} \frac{M(\alpha _1 +1, \gamma _1 +1, -z)}{M(\alpha _1, \gamma _1, -z) } \bigg ) \end{align}

and simplifying gives the result.

The asymptotic value of the forward rate is given by the following lemma.

Theorem 2. For the extended $3/2$ short rate model, with $p_t$ fulfilling Condition 2, the asymptotic instantaneous forward rate is

(6.18) \begin{align} g_\infty (t) = \alpha _1 \bar{p} \end{align}

where

(6.19) \begin{align} \alpha _1 = -\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2} \bigg ) + \sqrt{\bigg ( \frac{1}{2} - \frac{q}{\sigma ^2 } \bigg )^2 + \frac{2}{\sigma ^2}} . \end{align}

Proof. Let $\epsilon \gt 0$ be an arbitrarily small positive number and let $T= (1/\epsilon )^{1/\delta }\gt 0$ such that $|p_u-\bar{p}|\lt \epsilon$ for all $u\ge T$ . Also, let $\bar{T} = 2T$ . Then

(6.20) \begin{align} &\int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi \\ \notag &= \int _t^{T} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi + \int _T^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi \\ \notag &= \int _t^{T} \mathrm{exp} \{ - \int _\xi ^T p_u\, du \}\mathrm{exp} \{ - \int _T^{\bar{T}} p_u\, du \} \, d\xi + \int _T^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi \end{align}

Condition2 gives

(6.21) \begin{align} \mathrm{exp} \{ -K_2 (T-\xi )\} \le \mathrm{exp} \{ - \int _\xi ^T p_u\, du \} \le \mathrm{exp} \{ -K_1 (T-\xi )\} \end{align}

and, integrating with respect to $\xi$ over the interval $[t,T]$ , we have

(6.22) \begin{align} \frac{1}{K_2}\big ( 1-\mathrm{exp} \{ -K_2 (T-t )\}\big ) \le \int _t^{T}\mathrm{exp} \{ - \int _\xi ^T p_u\, du \}\, d\xi \le \frac{1}{K_1}\big ( 1-\mathrm{exp} \{ -K_1 (T-t )\}\big ) . \end{align}

Also

(6.23) \begin{align} \mathrm{exp} \{ -(\bar{T} - T) (\bar{p} +\epsilon )\} \le \mathrm{exp} \{ - \int _T^{\bar{T}} p_u\, du \} \le \mathrm{exp} \{ -(\bar{T} - T) (\bar{p} -\epsilon )\} . \end{align}

Finally, we have

(6.24) \begin{align} &\frac{1}{\bar{p} + \epsilon } \big ( 1-\mathrm{exp} \{ -(\bar{p}+\epsilon ) (\bar{T} - T )\}\big ) = \int _T^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} (\bar{p} + \epsilon ) \, du \} \, d\xi \\ \notag &\le \int _T^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi \le \int _T^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} (\bar{p} - \epsilon ) \, du \} \, d\xi \\ \notag &= \frac{1}{\bar{p} - \epsilon } \big ( 1-\mathrm{exp} \{ -(\bar{p}-\epsilon ) (\bar{T} - T )\}\big ) . \end{align}

Applying the inequalities (6.22), (6.23) and (6.24) to (6.20) gives

(6.25) \begin{align} &\frac{1}{K_2}\big ( 1-\mathrm{exp} \{ -K_2 (T-t )\}\big ) \mathrm{exp} \{ -(\bar{T} - T) (\bar{p} +\epsilon )\} \\ \notag &+ \frac{1}{\bar{p} + \epsilon } \big ( 1-\mathrm{exp} \{ -(\bar{p}+\epsilon ) (\bar{T} - T )\}\big ) \\ \notag &\lt \int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi \\ \notag &\lt \frac{1}{K_1}\big ( 1-\mathrm{exp} \{ -K_1 (T-t )\}\big ) \mathrm{exp} \{ -(\bar{T} - T) (\bar{p} -\epsilon )\} \\ \notag &+\frac{1}{\bar{p} - \epsilon } \big ( 1-\mathrm{exp} \{ -(\bar{p}-\epsilon ) (\bar{T} - T )\}\big ) . \end{align}

Hence,

(6.26) \begin{align} \frac{1}{\bar{p} + \epsilon } \lt \lim _{\bar{T}\to \infty } \int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi \lt \frac{1}{\bar{p} - \epsilon } . \end{align}

Since $\epsilon$ is arbitrary, we have

(6.27) \begin{align} \lim _{\bar{T}\to \infty } \int _t^{\bar{T}} \mathrm{exp} \{ - \int _\xi ^{\bar{T}} p_u\, du \} \, d\xi ={\frac{1}{\bar{p}}} \end{align}

and, inserting this into (6.12), we see that as ${\bar{T}}\to \infty$

(6.28) \begin{align} g_{\bar{T}}(t) = \alpha _1 \bar{p} \bigg ( 1 + o(1) \bigg ) \bigg ( 1 + O(z ) \bigg ) . \end{align}

Because $z\to 0$ as ${\bar{T}}\to \infty$ , we have the result.

Remark 3. Note that we may deduce Theorem 2 alternatively as follows. Applying L’Hôpital’s rule to (6.10) we obtain

(6.29) \begin{align} \alpha _1 \bar{p} = h_\infty (t) &= \lim _{{\bar{T}}\to \infty } \frac{\alpha _1 \log y(t,r_t)}{{\bar{T}}-t} \\ \notag &= \lim _{{\bar{T}}\to \infty } \frac{\frac{\partial }{\partial \bar{T}}(\alpha _1 \log y(t,r_t))}{\frac{\partial }{\partial \bar{T}}({\bar{T}}-t)} \\ \notag &= \alpha _1 \lim _{{\bar{T}}\to \infty } \frac{\partial }{\partial \bar{T}} \log y(t,r_t) . \end{align}

From the definition of $g_{\bar{T}}(t)$ in (6.11) and taking the partial derivative of both sides of (6.14) with respect to $\bar{T}$ , we have

\begin{align*} g_{\bar {T}}(t)= -\alpha _1 \frac {\partial }{\partial \bar {T}}\log z - \frac {\partial z}{\partial \bar {T}}\, \frac {\partial z}{\partial z}\log M(\alpha _1, \gamma _1, -z ) .\end{align*}

From (6.13) and (6.15) we have

\begin{align*} g_{\bar {T}}(t) =\alpha _1 \frac {\partial }{\partial \bar {T}}\log y(t,r_t) +z \frac {1}{\int _t^{\bar {T}} \mathrm{exp} \{ - \int _\xi ^{\bar {T}} p_u\, du \} \, d\xi }\, \frac {\partial }{\partial z}\log M(\alpha _1, \gamma _1, -z ) .\end{align*}

Taking limits of both sides as $\bar{T}\to \infty$ and, using $\lim _{\bar{T}\to \infty } z = 0$ and (6.27), gives

\begin{align*}g_{\infty }(t) = \alpha _1 \lim _{{\bar {T}}\to \infty } \frac {\partial }{\partial \bar {T}}\log y(t,r_t),\end{align*}

which, in tandem with (6.29), gives the result.

In Fig. 5 the instantaneous forward rate $g_{\bar{T}}$ is plotted and can be seen to be asymptotic to $g_\infty (t) = 0.0392$ based upon the parameters in (4.14).

7. Historical simulation and distribution of bond yields

The bond pricing formula in respect of a ten-year semi-annual coupon bond, with unit face value, coupon rate $g$ and quoted yield $y$ , is

(7.1) \begin{align} P(y) = \frac{1}{(1+y/2)^{20}} + \sum _{i=1}^{20} \frac{g/2}{(1+y/2)^i} . \end{align}

On the other hand, the same ten-year bond is valued at time $t$ under the $3/2$ short rate model as

(7.2) \begin{align} V_t \equiv V_t(r_t) = G_{t+10}(t) + \sum _{i=1}^{20} \frac{g}{2} G_{t+i/2}(t), \end{align}

and the implied ten-year bond yield under the $3/2$ model is found by solving $P(y) = V_t$ for $y$ . This allows us to compare the actual ten-year bond yield with that which is implied by the $3/2$ model.

Figure 7 Comparison of the distribution of annual changes in 10-year bond yields with those implied under the $3/2$ model using 10-year bond yield data 1871–2023.

In Fig. 7, the distribution of annual changes in 10-year bond yields is compared with those implied under the $3/2$ model and extended $3/2$ model using a 10Y bond yield data set for period 1871–2023. Our 10Y bond yield data set is obtained by augmenting Shiller’s long-term bond data 1871–2011, given in Shiller (Reference Shiller2022), with 10Y bond yield data for the period 2011–2023, given in YahooFinance (2023). It is evident that the $3/2$ short rate model is able to capture the peaked shape of the empirical density near zero, which is not the case for short rate models that induce normally distributed changes in ten-year bond yields.

In Fig. 8, the time series of actual ten-year bond yields is plotted alongside those for the $3/2$ model and extended $3/2$ models for $N=1,2,3$ . Note that the time series for $N=2$ and $N=3$ are visually indistinguishable. The summary statistics of the annual changes in ten-year bond yields are shown in Table 2, where the correlation of the annual changes in short rate with those of ten-year bond yields is closely matched by those corresponding to the extended $3/2$ models, i.e., for $N\ge 1$ . The ten-year bond yield for the $3/2$ model is highly correlated with the short rate, which is a consequence of the fixed parameters. The time-dependence of the parameters in the extended $3/2$ models permit this correlation with the short rate to be consistent with the actual correlation. However, the standard deviation of annual changes in ten-year bond yields is best matched by the $3/2$ model, with those for the extended $3/2$ models being dampened by the deterministic formulae for $p_t$ . The serial correlations of the $3/2$ model more closely match reality than those of the extended $3/2$ model. These stylized features of the extended $3/2$ models need to be borne in mind when used for modeling the yield curve.

8. Bond options and implied volatilities

In this section, we provide formulae for prices of bond options and provide estimates of the implied volatilities under each of the $3/2$ models fitted in Section 4.

Table 2. Summary statistics of annual changes in ten-year bond yields

Figure 8 Time series of actual 10-year bond yields and those implied by the $3/2$ model and extended $3/2$ models for $N=1,2,3$ over the period 1871–2023. Note that the time series for $N=2$ and $N=3$ are visually indistinguishable.

For a $\bar{T}$ -expiry call option, $0 \lt{\bar{T}} \le T \lt \infty$ , on a $T$ -maturity ZCB with strike price $K$ , we denote by $c_{{\bar{T}},T,K}(t)$ its actuarial price at time $t\in [0,T]$ . We have

(8.1) \begin{align} c_{{\bar{T}},T,K}(t) = E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} \big ( G_T({\bar{T}}) - K \big )^{+} \bigg |{\mathcal{A}}_t\bigg ) . \end{align}

Also, formulae for prices at time $t$ of a put option, asset-or-nothing call option, asset-or-nothing put option, cash-or-nothing call option, and cash-or-nothing put option, each on a $T$ -maturity ZCB, are given by

(8.2) \begin{align} p_{{\bar{T}},T,K}(t) &= E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} \big ( K-G_T({\bar{T}}) \big )^{+} \bigg |{\mathcal{A}}_t\bigg ) \\ \notag A^{+}_{{\bar{T}},T,K}(t) &= E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} G_T({\bar{T}}){\bf{1}}_{G_T({\bar{T}}) \gt K } \bigg |{\mathcal{A}}_t\bigg ) \end{align}

\begin{align*} A^{-}_{{\bar{T}},T,K}(t) &= E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} G_T({\bar{T}}){\bf{1}}_{G_T({\bar{T}}) \le K } \bigg |{\mathcal{A}}_t\bigg )\\ B^{+}_{{\bar{T}},T,K}(t) &= E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \}{\bf{1}}_{G_T({\bar{T}}) \gt K } \bigg |{\mathcal{A}}_t\bigg ) \\ \notag B^{-}_{{\bar{T}},T,K}(t) &= E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \}{\bf{1}}_{G_T({\bar{T}}) \le K } \bigg |{\mathcal{A}}_t\bigg ), \end{align*}

respectively. The indicator function is denoted above by ${\bf{1}}_{X\gt K}$ , which is equal to one if the random variable $X$ exceeds the value $K$ , and zero otherwise.

By virtue of the relations,

(8.3) \begin{align} c_{{\bar{T}},T,K}(t) &= A^{+}_{{\bar{T}},T,K}(t) - K\, B^{+}_{{\bar{T}},T,K}(t) \\ \notag p_{{\bar{T}},T,K}(t) &= K\, B^{-}_{{\bar{T}},T,K}(t) - A^{-}_{{\bar{T}},T,K}(t) \\ \notag A^{-}_{{\bar{T}},T,K}(t) &= G_T(t) - A^{+}_{{\bar{T}},T,K}(t) \\ \notag B^{-}_{{\bar{T}},T,K}(t) &= G_{\bar{T}}(t) - B^{+}_{{\bar{T}},T,K}(t), \end{align}

it suffices to determine formulae for $A^{+}_{{\bar{T}},T,K}(t)$ and $B^{+}_{{\bar{T}},T,K}(t)$ .

We have the following theorem, which gives an original formula for bond options under the extended $3/2$ short rate model.

Theorem 3. For the short rate $r_t$ obeying SDE (1.3), the formulae for prices of the asset-or-nothing call option and the cash-or-nothing call option are given by

(8.4) \begin{align} A^{+}_{{\bar{T}},T,K}(t) =& \int _{y^*}^\infty \frac{1}{2\tau } \left ( \frac{y}{x}\right )^{(\nu /2-1)/2} \mathrm{exp} \{ -\frac{1}{2}(x+y)/\tau \} \\ \notag &\times I_{\sqrt{(\nu /2-1)^2 + 8/\sigma ^2}} (\sqrt{xy}/\tau ) \, G(y)\, dy, \\[6pt] \notag B^{+}_{{\bar{T}},T,K}(t) =& \int _{y^*}^\infty \frac{1}{2\tau } \left ( \frac{y}{x}\right )^{(\nu /2-1)/2} \mathrm{exp} \{ -\frac{1}{2}(x+y)/\tau \} \\ \notag &\times I_{\sqrt{(\nu /2-1)^2 + 8/\sigma ^2}} (\sqrt{xy}/\tau ) \, dy \end{align}

where

\begin{align*} \tau &= \varphi _{\bar{T}} - \varphi _t = \frac{1}{4}\sigma ^2 \int _t^{\bar{T}} \mathrm{exp} \left (\int _0^s p_u\,du\right )\, ds, \\ G(y) &= \frac{\Gamma (\gamma _1 - \alpha _1)}{\Gamma (\gamma _1 )} z_y^{\alpha _1} M(\alpha _1, \gamma _1, -z_y ), \\ z_y &= \frac{2y}{\sigma ^2 \int _{\bar{T}}^T \mathrm{exp} (\int _0^s p_u\,du)\, ds }, \\ y^* &= \frac{1 }{2}z^* \sigma ^2 \int _{\bar{T}}^T \mathrm{exp} \left (\int _0^s p_u\,du\right )\, ds \end{align*}

and $z^*$ solves the equation

\begin{align*}\frac {\Gamma (\gamma _1 - \alpha _1)}{\Gamma (\gamma _1 )} {z^*}^{\alpha _1} M(\alpha _1, \gamma _1, -z^* ) = K.\end{align*}

Proof. We have from (8.2),

(8.5) \begin{align} A^{+}_{{\bar{T}},T,K}(t) &= E\bigg ( \mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} G_T({\bar{T}}){\bf{1}}_{G_T({\bar{T}}) \gt K } \bigg |{\mathcal{A}}_t\bigg ) \\ \notag &= E\bigg ( G_T({\bar{T}}){\bf{1}}_{G_T({\bar{T}}) \gt K } E\bigg (\mathrm{exp} \bigg \{-\int _t^{\bar{T}} r_s\, ds \bigg \} \bigg | r_{\bar{T}},\,{{\mathcal{A}}_t}\bigg ) \bigg |{\mathcal{A}}_t\bigg ) . \end{align}

Lemma4 gives

\begin{align*} E\bigg (\mathrm{exp} \bigg \{-\int _t^{\bar {T}} r_s\, ds \bigg \} \bigg | r_{\bar {T}}, \, {\mathcal {A}}_t\bigg ) = \frac {I_{\sqrt {8/\sigma ^2+(\nu /2-1)^2}} (\sqrt {xy}/\tau )}{I_{\nu /2-1} (\sqrt {xy}/\tau )},\end{align*}

where $ y = \mathrm{exp} (\int _0^{\bar{T}} p_u\,du) /r_{\bar{T}}$ and $ x = \mathrm{exp} (\int _0^t p_u\,du) /r_{t}$ .

Hence,

\begin{align*} &A^{+}_{{\bar{T}},T,K}(t) \\ &=E\bigg ( G_T({\bar{T}}){\bf{1}}_{G_T({\bar{T}}) \gt K } \frac{I_{\sqrt{8/\sigma ^2+(\nu /2-1)^2}} (\sqrt{xy}/\tau )}{I_{\nu /2-1} (\sqrt{xy}/\tau )} \bigg |{\mathcal{A}}_t\bigg ) \\ &= \int _0^\infty p_r(t,r_t,\bar{T},r_{\bar{T}})\, G_T({\bar{T}}){\bf{1}}_{G_T({\bar{T}}) \gt K } \frac{I_{\sqrt{8/\sigma ^2+(\nu /2-1)^2}} (\sqrt{xy}/\tau )}{I_{\nu /2-1} (\sqrt{xy}/\tau )} \, dr_{\bar{T}} . \end{align*}

Making the change of variable $ y = \mathrm{exp} (\int _0^{\bar{T}} p_u\,du) /r_{\bar{T}}$ in the integral above gives the result. The proof for $B^{+}_{{\bar{T}},T,K}(t)$ is similar.

To illustrate how derivative prices vary across the extended $3/2$ models, we compute the implied volatilities of at-the-money-forward 30-year ZCB call options as at 1871, for option expiries beyond this value date. These are shown in Fig. 9.

Figure 9 Plot of implied volatilities in respect of at-the-money-forward call options on thirty-year ZCBs for the $3/2$ model ( $N=0$ ) and extended $3/2$ models for $N=1,2$ , having value date 1871 and having expiry times in the period 1872–1978.

Figure 10 Plot of implied volatilities in respect of at-the-money-forward call options on thirty-year ZCBs, having ten years to expiry, for the $3/2$ model ( $N=0$ ) and extended $3/2$ models for $N=1,2$ , as at value dates in the period 1871–1979.

The calculation of the implied volatility $\sigma ^{BS}$ involves firstly computing the price of an at-the-money-forward 30-year ZCB call option, expiring at time $\bar{T}$ , using the formula

\begin{align*} c_{\bar{T},\bar{T}+30,K} (t) =& \int _{y^*}^\infty \frac{1}{2\tau } \left ( \frac{y}{x}\right )^{(\nu /2-1)/2} \mathrm{exp} \{ -\frac{1}{2}(x+y)/\tau \} \\ &\times I_{\sqrt{(\nu /2-1)^2 + 8/\sigma ^2}} (\sqrt{xy}/\tau ) \, \left ( G(y)- K\right )\, dy, \end{align*}

where the variables are given in Theorem3, with $K = G_{T}(t)/G_{\bar{T}}(t)$ . Then, the implied volatility is computed by solving

\begin{align*} c_{\bar {T},\bar {T}+30,K} (t) = G_T(t) \, \Phi (d_1) - K\, G_{\bar {T}}(t)\, \Phi (d_2) \end{align*}

for $\sigma ^{BS}$ . Here, $\Phi (x)$ is cumulative distribution function of the standard normal distribution, $T = \bar{T} + 30$ ,

\begin{align*}d_1 = \{ \log (G_T(t)/(G_{\bar {T}}(t)\, K)) + \frac {1}{2} (\sigma ^{BS})^2 (\bar {T}-t) \} / (\sigma ^{BS} \sqrt {\bar {T}-t})\end{align*}

and $d_2 = d_1 - \sigma ^{BS} \sqrt{\bar{T}-t}$ .

In Fig. 10 are shown implied volatilities of call options on forward thirty-year ZCBs, with ten years to expiry, for the various $3/2$ models. We note that the implied volatilities for the ZCB call option under the $3/2$ model are approximately 0.09, which for a thirty-year ZCB, corresponds to a standard deviation of annual changes in the thirty-year ZCB yield of 0.09/30 = 0.003. This is around half of the standard deviation of annual changes in actual ten-year bond yields, which is shown in Fig. 7 as 0.007.