Asymptotic behavior of bond yields and volatilities for the extended 3/2 model under the real-world measure

Abstract

The extended $3/2$ short rate model is a mean reverting model of the short rate which, for suitably chosen parameters, permits a sensible term structure of bond yields and closed-form valuation formulae of zero-coupon bonds and options. This article supplies proofs of the formulae for the expected present values of future cash flows under the real-world probability measure, known as actuarial valuation. Finally, we give formulae for asymptotic levels of bond yields and formulae for bond option prices for the extended $3/2$ model, under particular conditions on its parameters.

1. Introduction

The choice of discount rate that is employed in actuarial valuations of pension funds and life insurance portfolios depends upon the duration of cashflows, rates assumed in previous actuarial reviews and, of course, actuarial judgement. Nowadays, actuaries use stochastic asset models, a key component of which is a stochastic interest rate term structure model. Much of the research on stochastic interest rate models has focussed on interest rate derivatives pricing, varying in their complexity from single-factor short rate models, as described in Rebonato (1998) and Brigo & Mercurio (2006), to multi-factor forward rate models that cater for non-linear behaviors, as pioneered by Heath et al. (1992) and Brace et al. (1997). However, less research has been devoted to stochastic models in actuarial valuations, with some work on long-term yields being done by Dufresne (1990), Carriere (1999) and Thomson & Gott (2009).

Single factor term structure models have been comprehensively studied in the late 1990’s and early 2000’s, particularly by Yacine Aït-Sahalia who analyzed many relevant properties for a wide variety of single factor models, including the 3/2-model. Nonlinear diffusions of the short rate are examined in Aït-Sahalia (1999), whereby the transition density of the nonlinear SDE is approximated by a closed-form expression, thus allowing maximum likelihood estimation of the parameters. Such an approach is generalised to discretely sampled processes in Aït-Sahalia (2002). A nonparametric approach is described in Aït-Sahalia (1996a), whereby the forms of the drift and diffusion terms are inferred from the data, and the associated PDE for the price of the derivative, e.g., bond option or zero-coupon bond (ZCB), is solved.

One advantage of single-factor short rate models for modeling interest rates is their parsimony, yet this is also a disadvantage when pricing derivatives as they lack the flexibility in being able to capture nonparallel shifts in yield curves. While the short-term interest rate is important in valuations of future cashflows pertaining to short-tail insurance policies, its effect on long-term rates is not clear because the long-term rates that are implied by a mathematical model of the instantaneous, continuously compounded deposit rate are dependent upon the parameters of the short rate model. Importantly for actuaries, the long-term bond yield implied by the model should be a finite constant and our aim in this paper is to examine conditions on parameters of the well known $3/2$ short rate model, derived in Platen (1999) and studied by Ahn & Gao (1999), which guarantee this. Previously, considerations of this kind for Vasicek and CIR models have been undertaken in Fergusson (2017) and (2019), respectively.

In Fig. 1, the asymmetry of the distribution of annual changes in the US short rate is shown, which is symptomatic of the leverage effect in bond markets. While this effect is captured by the $3/2$ model, it is a short-term effect which is less pronounced when analyzing the asymptotic behavior of bond yields.

Figure 1 Comparison of empirical probability density function of annual change in short rate with that of the fitted normal distribution (US one-year cash rates 1871–2023).

Unlike the Vasicek model of Vasicek (1977) and the CIR model of Cox et al. (1985), the $3/2$ model is a nonlinearly mean reverting stochastic model. The mean-reverting aspect ensures that interest rates adhere to a long-run reference level. Working on a filtered probability space $(\Omega, {\underline{\mathcal{A}}}, ({\mathcal{A}}_t)_{t\ge 0}, P)$ , the SDE for the $3/2$ model is

(1.1) \begin{align} dr_t = (\bar{p} r_t + q r_t^2 ) dt + \sigma r_t^{3/2} dZ_t, \end{align}

where

(1.2) \begin{align} \bar{p}&\gt 0 \\ q &\lt \frac{1}{2}\sigma ^2 \notag \\ \sigma &\gt 0,\notag \end{align}

so as to avoid explosive values of $r_t$ . The $3/2$ power law for the diffusion coefficient in a short rate model was shown in Chan et al. (1992) to be the best fitting power law. Also the nonlinear drift term of this model could not be rejected in Aït-Sahalia (1996b). As for the CIR model, but unlike the Vasicek model, the $3/2$ model avoids the possibility of negative interest rates.

The extended $3/2$ model has the SDE

(1.3) \begin{align} dr_t = (p_t r_t + q r_t^2 ) dt + \sigma r_t^{3/2} dZ_t, \end{align}

where $p_t$ is a positively valued function, and $q$ and $\sigma$ satisfy (1.2). The risk-neutral ZCB pricing formula for this model is given in Lo (2013), and later in this article we supply an original proof of the pricing formula for real-world short rate dynamics. Moreover, we state Conditions1 and 2 on $p_t$ which ensure that the long-term bond yield is a finite constant, with Condition2 also ensuring that the long-term instantaneous forward rate is a finite constant.

Condition 1. There exist a constant $\bar{p}\gt 0$ and a positive function $K(x)$ , satisfying $\lim _{x\to \infty }{} K(x) = 0$ , such that for all $s$ and $t$ , with $0\lt s\lt t$ ,

\begin{align*} \sup _{\xi \in [s,t]}\bigg | \int _s^\xi (p_u - \bar {p})\, du \bigg | \lt (t-s)\, K(t-s) .\end{align*}

Condition 2. There exist constants $K_1\gt 0$ , $K_2\gt 0$ and $\bar{p}\gt 0$ such that $\lim _{t\to \infty } p_t = \bar{p}$ and for all $s$ and $t$ , with $0\lt s\lt t$ ,

\begin{align*} K_1 (t-s)\le \int _s^t p_u \, du \le K_2 (t-s).\end{align*}

Note that if $p_t$ satisfies Condition2, then it automatically satisfies Condition1; see the appendix for a proof. Furthermore, the constant $\bar{p}$ appearing in each condition is necessarily unique. Clearly, the $3/2$ model, with $p_t = \bar{p}\gt 0$ , fulfils Conditions1 and 2. However, the extended $3/2$ model, with $p_t = \bar{p} \big ( 1+\frac{1}{2} \sin \frac{2\pi t}{14}\big )$ accommodating a 14-year interest rate cycle, fulfils only Condition1.

In Section 2, we give some background on actuarial valuations. In Section 3, we provide an explicit solution to the SDE (1.1). In Section 4, we calculate the transition density function. In Section 5, we provide an original proof of the closed-form formula for the ZCB price. In Section 6, we determine the long-term bond yield and instantaneous forward rate. Then in Section 7, we compare the distribution of annual changes in the 10-year bond yield with those implied by the $3/2$ model. In Section 8, we provide formulae for prices of bond options and compute the implied volatilities under each of the $3/2$ models. Finally, in Section 9, we conclude. Because our valuations of ZCBs are under the real-world measure, we supply all steps of the derivations, despite the fact that some of these are well known under risk-neutral assumptions.

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 (1994), 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 (2002) and Platen & Heath (2006), 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 (2006), 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 (2013). 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 (2006) 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 (2019), 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 (2019) 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 (2022), 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 (2021), 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.

Remark 2. The dimension of the square root process $1/r_t$ is here estimated as $\nu = 4(\sigma ^2-q)/\sigma ^2 \approx 2.9818$ , which is reasonably close to three.

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 (2007) and proofs involving Lie symmetries of partial differential equations are provided in Craddock & Lennox (2009) and Section 5.5 of Baldeaux & Platen (2013). A proof of the following lemma which uses the Lie algebra of differential operators is given in Lo (2013) 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 (2007), 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 (2006), 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 (1966) 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 (2013) 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}

Proof. Substituting $u=1$ in Lemma4 gives the result.

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 (1991). 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 (2022), 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.

9. Conclusions

We have demonstrated how ZCBs and bonds, and options on ZCBs, can be valued in an actuarial and real-world probabilistic framework for the extended $3/2$ short rate model and have determined formulae for the asymptotic levels of bond yields. The examples illustrate that sensible long-term bond yields result when conditions on $p_t$ , $q$ , and $\sigma$ are imposed and when implied from market data. In the light of the model’s suitability for short rate modeling, as described in Chan et al. (1992) and Ahn & Gao (1999), we have highlighted its potential for modeling long-term rates in actuarial valuation models.