Calibration of the Ueno’s Shadow Rate Model of Interest Rates

Abstract

Shadow rate models of interest rates are based on the assumption that the interest rates are determined by an unobservable shadow rate. This idea dates back to Fischer Black, who understood the interest rate as an option that cannot become negative. Its possible zero values are consequences of negative values of the shadow rate. In recent years, however, the negative interest rates have become a reality. To capture this behavior, shadow rate models need to be adjusted. In this paper, we study Ueno’s model, which uses the Vasicek process for the shadow rate and adjusts its negative values when constructing the short rate. We derive the probability properties of the short rate in this model and apply the maximum likelihood estimation method to obtain the parameters from the real data. The other interest rates are—after a specification of the market price of risk—solutions to a parabolic partial differential equation. We solve the equation numerically and use the long-term rates to fit the market price of risk.

1. Introduction

Interest rates are fundamental to the economy. Understanding their dynamics is essential in situations where either discounting of cash flows is performed or the interest rate influences the variables of interest (investing in bonds, interest rate derivatives, etc.). Therefore, their mathematical modeling is necessary as a part of the modeling of the original problems. Also, since interest rates are an important indicator of the state of the economy, information about their behavior and tendencies is of interest by itself. We refer the reader to books [1,2] for a complex coverage of interest rate modeling.

Short rate models of interest rates are a class of models in which the instantaneous interest rate (called short rate) is modeled by a stochastic differential equation. A concept of shadow rate dates back to Black [3], who argued that the nominal interest rate cannot become negative, but its zero values can be observed for a certain time interval. This can be explained by the existence of a so-called shadow rate. In the case it is positive, the short rate coincides with it. However, if the shadow rate becomes negative, the short rate will be equal to zero.

At the time when negative interest rates were considered to be unrealistic, it was this zero lower bound that characterized shadow rate models. This property can be achieved in different ways. Usually, a particular choice of the process in Black’s framework is given. Often, it is an affine combination of Gaussian factors; see [4,5,6,7]. There are also other possible approaches: the sum of squares of the factors and the sum of exponentials were considered in [8]. Paper [9] studies a shadow rate version of the model from [10], which is based on the Nelson–Siegel model [11]. A nonzero lower bound can be found in some of the papers, either as a note about a possible generalization or as a particular question of interest; see [8,9]. When interest rates in some countries fell to negative values, a modification of Black’s model became necessary. This motivated different ways of incorporating this feature into modeling. Dividing the time period into subperiods with different lower bounds was considered in [12], and several ways of estimating the lower bound from the data were implemented in [13]. The stochastic character of the shadow rate was proposed in [14,15], while Ref. [15] also assumed dependence of the shadow rate on an exogenous variable.

Since the factors are not observable, a typical approach to statistical analysis of interest rate time series is the Kalman filter and the extended Kalman filter. Pricing bonds, which are used to evaluate interest rates with different maturities, are complicated, and several different methods can be found in the literature. They include a numerical solution of the partial differential equation in [8], an approximation based on options in [6,9], the computation of forward rates in [15], and eigenfunction expansion in [16].

Because of the complicated calculations in shadow rate models, the simple Vasicek model [17] with a Gaussian short-term rate and an explicit formula for the bond prices is often employed to allow negative interest rates when the interest rate model is not the main question but one of the building blocks of an analysis. Such complex problems include a pension planning analysis [18,19], a portfolio optimization [20], and a stochastic optimization arising from a reinsurance–investment problem [21]. Therefore, there is undoubtedly a demand for uncomplicated interest rate models that still capture important features of the observed interest rates.

Ueno in [14] considered a rather complicated model and subsequently its simplified versions. One of the models for the short rate $r_t$ and the one that we study in this paper was given by the Vasicek-type (cf. [17] for the original short rate model) shadow rate $s_t$ and the transformation rule $r_t=max(s_t,k{s}_t)$ for $k\in (0,1]$. We show that this particular model has the potential to serve as a computationally tractable model that, at the same time, allows different dynamics for positive and negative interest rates.

The paper is organized as follows. In Section 2, we describe the model under consideration. Subsequently, in Section 3, we derive the probability density of the short rate in this model, and in Section 4, we derive the conditional expected value of the short rate and note some of its properties. Before the calibration of the model, we introduce the data in Section 5. Firstly, we use the density to construct the likelihood function in Section 6, and we obtain the maximum likelihood estimates for the market data of short-term interest rates. Afterward, in Section 7, we use the method of lines to numerically compute the solution of the partial differential equation for the bond prices. Using this computation, we estimate the market price of risk of the model by fitting the long-term interest rates. A discussion of the results and further considerations are given in Section 8.

2. Shadow Rate Model of the Short Rate

The first step in the analysis of the model is its calibration using the maximum likelihood method; therefore, we define it the under the physical (instead of risk-neutral) probability measure. Since it is a model considered in [14] as a special case of the presented general framework, we refer to it as Ueno’s model.

Definition 1.

Let $(\Omega ,\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t\in [0,T]},\mathbb{P})$ be a complete probability space and $w_t$ be a Wiener process defined on this probability space. Here, ${\left\{{\mathcal{F}}_t\right\}}_{t\in [0,T]}$ is the filtration generated by the Wiener process $w_t$, and it represents the information up to time t. The short rate $r_t$ in the Ueno’s model is defined as

$$r_t=max(s_t,k{s}_t),d{s}_t=\kappa (\theta -s_t)dt+\sigma d{w}_t,$$

(1)

where $\kappa ,\sigma >0$, $\theta \in \mathbb{R}$, $k\in (0,1]$ are constant parameters.

If $k=1$, the model becomes the original classical Vasicek model [17]. However, if $k\in (0,1)$, the short rate equals the shadow rate governed by the Vasicek model only for its positive values. When the shadow rate becomes negative, the short rate is equal to only a k-multiple of the shadow rate. Since $k\in (0,1)$, it is “less negative” than the shadow rate. This behavior is illustrated in Figure 1 for two different values of the parameter k.

3. Conditional Probability Distribution of the Short Rate

The derivation of the probability distribution of the short rate is based on the following two observations: Firstly, for the given value of k and the value of the short rate $r_t$, we can reconstruct the value of the shadow rate $s_t$. Secondly, when evaluating the probability of the inequality that determines the value of the cumulative distribution function of the short rate in the argument x, it can be transformed into an inequality for the shadow rate s.

Theorem 1.

Let $r_t$ be a short rate process given by Definition 1. Then, the conditional density of $r_{t+\Delta t}$, given $r_t$, is

$$f_{r_{t+\Delta t}}\left(x\right|r_t)=\left\{\begin{array}{cc}f\left(x;\mu =r_t{e}^{-\kappa \Delta t}+\theta (1-e^{-\kappa \Delta t}),{\Sigma }^2\right)\hfill & \mathrm{if}{r}_t\ge 0,x>0,\hfill \\ f\left(x;\mu ={\displaystyle \frac{1}{k}}{r}_t{e}^{-\kappa \Delta t}+\theta (1-e^{-\kappa \Delta t}),{\Sigma }^2\right)\hfill & \mathrm{if}{r}_t<0,x>0,\hfill \\ {\displaystyle \frac{1}{k}}f\left({\displaystyle \frac{1}{k}}x;\mu =r_t{e}^{-\kappa \Delta t}+\theta (1-e^{-\kappa \Delta t}),{\Sigma }^2\right)\hfill & \mathrm{if}{r}_t\ge 0,x<0,\hfill \\ {\displaystyle \frac{1}{k}}f\left({\displaystyle \frac{1}{k}}x;\mu ={\displaystyle \frac{1}{k}}{r}_t{e}^{-\kappa \Delta t}+\theta (1-e^{-\kappa \Delta t}),{\Sigma }^2\right)\hfill & \mathrm{if}{r}_t<0,x<0,\hfill \end{array}\right.$$

(2)

where f is the probability density of normal distribution with specified expected value μ and variance ${\Sigma }^2={\displaystyle \frac{{\sigma }^2}{2\kappa }}\left(1-e^{-2\kappa \Delta t}\right)$.

Proof. 

Clearly,

$$F_{r_{t+\Delta t}}(x|r_t)=\mathbb{P}(r_{t+\Delta t}<x|r_t)=\left\{\begin{array}{cc}\mathbb{P}(s_{t+\Delta t}<x|s_t=r_t)\hfill & \mathrm{if}{r}_t\ge 0,x\ge 0,\hfill \\ \mathbb{P}(s_{t+\Delta t}<x|s_t={\displaystyle \frac{1}{k}}{r}_t)\hfill & \mathrm{if}{r}_t<0,x\ge 0,\hfill \\ \mathbb{P}(s_{t+\Delta t}<{\displaystyle \frac{1}{k}}x|s_t=r_t)\hfill & \mathrm{if}{r}_t\ge 0,x<0,\hfill \\ \mathbb{P}(s_{t+\Delta t}<{\displaystyle \frac{1}{k}}x|s_t={\displaystyle \frac{1}{k}}{r}_t)\hfill & \mathrm{if}{r}_t<0,x<0.\hfill \end{array}\right.$$

(3)

We recall that the distribution of the shadow rate $s_t$, governed by the Vasicek process, is given by [17]:

$$s_{t+\Delta t}|s_t\sim \mathcal{N}\left(s_t{e}^{-\kappa \Delta t}+\theta \left(1-e^{-\kappa \Delta t}\right),{\displaystyle \frac{{\sigma }^2}{2\kappa }}\left(1-e^{-2\kappa \Delta t}\right)\right).$$

(4)

By substituting probabilities given by the distribution (4) into (3) and differentiating the resulting cumulative distribution function, we obtain the conditional distribution of the short rate, as expressed in (2). □

We note that the cumulative distribution function is not differentiable for $x=0$ and the density is not continuous at this point. Examples of the densities are shown in Figure 2.

4. Conditional Expected Value of the Short Rate

It is possible to derive the expected value of the short rate, given its current value, by integrating the density obtained in the previous section or by integrating the normal density. Alternatively, its derivation can be simplified by using known probability distributions, as we do in the proof of the following theorem.

Theorem 2.

Let $r_t$ be a short rate process given by Definition 1. Then,

1. 

The conditional expected value $r_t$, given $r_0$, is

$$\mathbb{E}\left(r_t\right|r_0)={\displaystyle \frac{1}{2}}\left((1-k)\left[\sqrt{{\displaystyle \frac{2}{\pi }}}\Sigma e^{-{\displaystyle \frac{{\mu }^2}{2{\Sigma }^2}}}+\mu \left(1-2\Phi \left(-{\displaystyle \frac{\mu }{\Sigma }}\right)\right)\right]+(1+k)\mu \right),$$

(5)

where Φ is the cumulative distribution function of the normalized normal distribution and

$$\mu =s_0{e}^{-\kappa t}+\theta \left(1-e^{-\kappa t}\right),{\Sigma }^2={\displaystyle \frac{{\sigma }^2}{2\kappa }}\left(1-e^{-2\kappa t}\right),s_0=\left\{\begin{array}{cc}{r}_0\hfill & for{r}_0\ge 0,\hfill \\ {\displaystyle \frac{1}{k}}{r}_0\hfill & for{r}_0<0.\hfill \end{array}\right.$$

2. 

As time t approaches infinity, the conditional expected value (5) converges to a value independent of $r_0$, which equals

$$\underset{t\to \infty }{lim}\mathbb{E}\left(r_t\right|r_0)={\displaystyle \frac{1}{2}}\left((1-k)\left[{\displaystyle \frac{\sigma }{\sqrt{\pi \kappa }}}{e}^{-{\displaystyle \frac{{\theta }^2\kappa }{{\sigma }^2}}}+\theta \left(1-2\Phi \left(-{\displaystyle \frac{\sqrt{2\kappa }\theta }{\sigma }}\right)\right)\right]+(1+k)\theta \right).$$

Proof. 

Using a basic property $max(a,b)={\displaystyle \frac{1}{2}}\left(|a-b|+a+b\right)$ and the definition of the short rate $r_t$ we can write

$$r_t=max(s_t,k{s}_t)={\displaystyle \frac{1}{2}}\left(|s_t-k{s}_t|+s_t+k{s}_t\right)={\displaystyle \frac{1}{2}}\left((1-k)|s_t|+(1+k)s_t\right).$$

For the given value of the short rate $r_0$, the corresponding value of the shadow rate $s_0$ can be reconstructed, and we can use (4) to obtain the conditional expected value of $s_t$. Furthermore, if a random variable Y follows a normal $\mathcal{N}(\mu ,{\Sigma }^2)$ distribution, then its absolute value has the so-called folded normal distribution [22], and its expected value is known to be

$$\mathbb{E}\left(\right|Y\left|\right)=\sqrt{{\displaystyle \frac{2}{\pi }}}\Sigma e^{-{\displaystyle \frac{{\mu }^2}{2{\Sigma }^2}}}+\mu \left(1-2\Phi \left(-{\displaystyle \frac{\mu }{\Sigma }}\right)\right).$$

(6)

Substituting the distribution of the shadow rate (4) into (6) gives the conditional expected value of $|s_t|$. These expected values imply the Formula (5) for the conditional expected value of the short rate. Its limit for time approaching infinity is then a direct consequence. □

Examples of the expected values for selected initial values $r_0$ and parameters are depicted in Figure 3. It is important to note that, unlike in the case of the Vasicek model, the expected value can be non-monotone. Furthermore, the convergence of the conditional expected values for time approaching infinity is illustrated in Figure 4.

5. Data Used in Calibration

We use the OECD monthly data for short-term [23] and long-term [24] interest rates in Denmark (DNK), Euro Area (EA), and Sweden (SWE). The time period was chosen to be 2011–2020, which spans a period of 10 years; these datasets include a sufficient number of both positive and negative values, which is necessary for estimating the model under consideration. By focusing on this time period, we intended to investigate the alterations in the dynamics of interest rates that occurred after their transition into negative values.

According to OECD, the short-term interest rates are the rates at which short-term borrowings are effected between financial institutions or the rate at which short-term government paper is issued or traded in the market [23]. The long-term interest rate refers to government bonds maturing in ten years [24]. We consider the short-term rates as a proxy to the short rate (therefore, we refer to them as short rates), and the particular maturity of the long-term rates (long rates hereafter) is important when we fit them using the computed bond prices. Finally, we convert the interest rates from the original percentage points to decimal numbers since this representation is assumed in the derivation of the partial differential equation for the bond prices (see [25] for the derivation).

The evolution of the data over time is given in Figure 5.

6. Maximum Likelihood Estimation

The likelihood function for the given time series of the short rate is given as a product of the densities (2). Its maximization with respect to $\kappa$, $\theta$, $\sigma$, and k gives the maximum likelihood estimates (MLEs) of the parameters. It is important to note that when maximizing the likelihood function for a fixed value of k, the optimal values of the remaining parameters can be found in a closed form using the formulae from the Vasicek model estimation (see, for example, [1] for the expressions) since the evolution of the shadow rate can be reconstructed. Therefore, the optimization, which is used to compute the MLE of the shadow rate model, is essentially one-dimensional. For any given value of k, the optimal value of the likelihood is computed, and the one-dimensional optimization is performed with respect to k.

We summarize the results in

Table 1. Results of maximum likelihood estimation.

	DNK	EA	SWE
Estimate of k	0.67051	0.38895	0.68963
Estimate of $\kappa$	0.32729	0.06004	0.17251
Estimate of $\theta$	−0.00438	−0.04235	−0.01023
Estimate of $\sigma$	0.00320	0.00250	0.00312
p-Value for $H_0:k=1$	0.00540	0.00000	0.00427
Interval estimate of k	$(0.50086,0.88908)$	$(0.30041,0.50042)$	$(0.52968,0.89018)$

. For every dataset, we provide the estimates of the parameters $\kappa$, $\theta$, $\sigma$, and k. We test the null hypothesis $k=1$ (which corresponds to the Vasicek model) using the likelihood ratio test; it is rejected for all cases. Furthermore, since the parameter k determines the shadow rate, which is our point of interest, we provide an interval estimate of this parameter. In particular, we compute the interval of the values of k, which are not rejected by the likelihood ratio test on the usual 5 percent significance level. Figure 6 shows the estimated shadow rate and the shadow rates corresponding to boundary values of the interval estimates of k.

7. Pricing Bonds and Estimating the Market Price of Risk

A discount bond is a security that pays its holder a unit amount of money at a specified time T, called the maturity of the bond. If the short rate $r=r\left(x\right)$ is defined via a stochastic process x governed by a stochastic differential equation

$$dx=\mu (x,t)dt+\sigma (x,t)dw,$$

(7)

then, after the specification of the so-called market price of risk $\lambda =\lambda (x,t)$, the price of a discount bond $P=P(x,t)$ is a solution of a parabolic partial differential equation

$${\displaystyle \frac{\partial P}{\partial t}}+(\mu (x,t)-\lambda (x,t)\sigma (x,t)){\displaystyle \frac{\partial P}{\partial x}}+{\displaystyle \frac{{\sigma }^2(x,t)}{2}}{\displaystyle \frac{{\partial }^{2}P}{\partial x^2}}-r\left(x\right)P=0,$$

where $t\in (0,T)$ and the range of x is given by the range of values attainable by the stochastic process (7). The terminal condition is given by $P(x,T)=1$ for all x. Let us note that this approach to pricing bonds (see [25], [Chapter 7.2.1]) is based on the construction of a portfolio of bonds, for which its stochastic part vanishes. As a consequence, the partial differential Equation (7) emerges as a condition satisfied by the bond price. Another approach leading to Equation (7) is based on the so-called risk-neutral probability measure $\mathbb{Q}$, which is used to express the bond price as the expected value ${\mathbb{E}}_t^{\mathbb{Q}}\left[e^{-{\int }_t^{T}r\left(x_s\right)ds}\right]$. The Feynman–Kac formula then enables us to obtain the partial differential equation satisfied by the expected value and, therefore, the bond price. The change in the measure and the specification of market prices of risk are closely related; we refer the reader to [25] for more details regarding pricing bonds in interest rate models.

Let T be a fixed maturity. The interest rate $R(x,t)$ at time $t\in [0,T]$ when the stochastic factor equals x is connected with the price of a discount bond through the relation $P(x,t)=e^{-R(x,t)(T-t)}$. This allows us to compute the bond prices by solving the partial differential equation, transform them into interest rates, and compare them with market data.

We consider Ueno’s shadow rate model for two choices of the market price of risk:

• Constant market price of risk $\lambda$;
• Continuous piecewise linear market price of risk, taking ${\lambda }_1$ for a negative shadow rate and ${\lambda }_2$ for a shadow rate greater than 0.01 (i.e., one percent).

We note that these boundary values, when the market price of risk changes from its constant value, correspond to the same values of the observed short rate since they are both nonnegative. However, the values of zero and one percent were chosen only for illustration purposes and can also be taken as parameters that are being estimated.

We use the methodology from [8]: we consider the boundary conditions given by vanishing second derivative with respect to x at $\pm \infty$ and the method of lines for solving the partial differential equation for bond prices. In order to make the paper self-contained, we review the method of lines according to [8] while performing necessary modifications in order to apply it to our model. The discretization is achieved in space variable x, and the solution is computed on grid points $x_1,\cdots ,x_J$, uniformly distributed over the interval $[x_0,x_J]$ with uniform step $\delta x$. The time t is transformed into the time remaining to maturity $\tau =T-t$, where T is the maturity of the bond. Then, after denoting the vector of values at the grid points by $u\left(\tau \right)={(u_1\left(\tau \right),\cdots ,u_J\left(\tau \right))}^{\prime }$, the discretization that we use takes the form of the system of the ordinary differential equations

$${\displaystyle \frac{\partial u_j}{\partial \tau }}={\displaystyle \frac{{\sigma }^2}{2}}{\displaystyle \frac{u_{j+1}-2u_j+u_{j-1}}{{\left(\delta x\right)}^2}}+(\kappa (\theta -x_j)-\lambda \left(x_j\right)\sigma ){\displaystyle \frac{u_{j+1}-u_{j-1}}{2\delta x}}-max(x_j,k{x}_j)u_j=0$$

for the interior points. Taking the boundary conditions into account, the system can be written in the matrix form as

$${\displaystyle \frac{\partial u}{\partial \tau }}=\mathbb{M}u,u\left(0\right)=1_J,$$

where $1_J$ is a J-dimensional vector with all elements equal to 1, ${\displaystyle \frac{\partial u}{\partial \tau }}$ is computed by differentiating the components of the vector u, and

$$\mathbb{M}=\left(\begin{array}{ccccccc}2A_1+B_1& C_1-A_1& 0& ..& 0& 0& 0\\ A_2& B_2& C_2& ..& 0& 0& 0\\ ..& ..& ..& ..& 0& 0& 0\\ 0& 0& 0& ..& A_{J-1}& B_{J-1}& C_{J-1}\\ 0& 0& 0& ..& 0& A_J-C_J& B_J+2C_J\end{array}\right)$$

with

$$A_j={\displaystyle \frac{1}{2}}\left({\left({\displaystyle \frac{\sigma }{\delta x}}\right)}^2-{\displaystyle \frac{\kappa (\theta -x_j)-\lambda \left(x_j\right)\sigma }{\delta x}}\right),C_j={\displaystyle \frac{1}{2}}\left({\left({\displaystyle \frac{\sigma }{\delta x}}\right)}^2+{\displaystyle \frac{\kappa (\theta -x_j)-\lambda \left(x_j\right)\sigma }{\delta x}}\right),$$

$$B_j=-max(x_j,k{x}_j)-{\left({\displaystyle \frac{\sigma }{\delta x}}\right)}^2.$$

The solution can be written in terms of a matrix exponential as

$$u\left(\tau \right)=e^{\mathbb{M}\tau }.$$

To evaluate the matrix exponential, we use Padé’s computation, as implemented in the expm package [26] for R. We use the same number of discretization points (i.e., 200) as in [8], but we use a smaller interval for x, namely, $(-0.2,0.2)$ instead of $(-1,1)$. We test the accuracy on the Vasicek model, which is a special case for $k=1$ and which has an explicit solution for the bond prices; see [17]. We consider the parameter values estimated for the Vasicek model by the maximum likelihood method for EA data, which are equal to $\kappa =0.14271$, $\theta =-0.01033$, and $\sigma =0.00181$. The market price of risk is taken to be zero. Since we are going to fit the 10-year interest rates, we consider the absolute error of the 10-year interest rates when using the numerical method for bond prices, as described above, and computing the corresponding interest rates.

Table 2. Characteristics of absolute error of the interest rates.

	Test 1	Test 2	Test 3
Minimum	2.839 × 10−9	2.839 × 10−9	1.906 × 10−8
First quartile	1.798 × 10−7	1.656 × 10−8	3.967 × 10−8
Median	3.495 × 10−7	4.371 × 10−8	6.027 × 10−8
Mean	4.920 × 10−7	4.550 × 10−8	6.027 × 10−8
Third quartile	5.230 × 10−7	7.086 × 10−8	8.088 × 10−8
Maximum	3.666 × 10−6	9.801 × 10−8	1.015 × 10−7

provides the statistics about the error in the following cases:

• Test 1. Interest rates computed for all the discretization points.
• Test 2. Interest rates computed for the discretization points withing the range of EA short rate.
• Test 3. Using linear interpolation of the computed interest rates for the vector of x values with the step 10⁻⁶ (the precision used to quote the EA data), with minimum and maximum values equal to those of the EA short rate.

As we can see in Table 2, the numerical results are already more precise than the precision of the data. Furthermore, in the paper [8], the author refers to a good precision of this method when applied to Black’s model with the short rate given by $r_t=max(s_t,0)$, where $s_t$ follows a process of the Vasicek type.

We compare the calibration results for the proposed shadow rate model and the classical Vasicek model (i.e., using $k=1$ to obtain maximum likelihood estimates of $\kappa$, $\theta$, and $\sigma$) with its usual choice of a constant market price of risk.

For every dataset, we define the objective function equal to the sum of squares of differences between observed long rates and the long rates implied by the model for a particular choice of market price of risk. The latter are computed by linear interpolation of the interest rates obtained by the numerical method described above, while the remaining parameters of the model (i.e., $\kappa$, $\theta$, $\sigma$, and k) are taken from the maximum likelihood estimation from the previous section. In the case of a constant market price of risk, the resulting optimization problem is one-dimensional. If the market price of risk is switching its value in the way specified above, the optimization is two-dimensional, and as an initial approximation, we use the point corresponding to the optimal constant market price of risk. We use the L-BFGS-B method to solve the resulting optimization problems. For comparison, we compute the optimal constant market price of risk in the Vasicek model, in which case the objective function is quadratic.

When interpreting the quality of fit, we need to keep in mind the difference between moving from the Vasicek model to a shadow rate model with a constant market price of risk and moving from a constant market price of risk to switching in the shadow rate model. While the shadow rate model is a generalization of the Vasicek model when both have constant market prices of risk (the Vasicek model corresponds to $k=1$), it is no longer the case at the moment of fitting the long-term rates in our approach. At this moment, the parameters from the maximum likelihood step are fixed, and therefore, it is not necessarily true that the shadow rate model cannot perform worse than the Vasicek model. On the other hand, when specifying the market price of risk in the shadow rate model, the fixed parameters are the same. Therefore, since considering the same market prices of risk in the switching case reduces the model to a constant market price of risk, the switching model cannot lead to a worse fit.

The optimal values of the objective function are presented in

Table 3. Optimal values of the objective function when fitting the long-term interest rates with different specifications of the market price of risk (MPR).

	DNK	EA	SWE
Vasicek model (constant MPR)	0.00755	0.01235	0.00438
Shadow rate model with constant MPR	0.00742	0.00730	0.00417
Shadow rate model with switching MPR	0.00222	0.00254	0.00299

and Figure 7; the estimated market values of risk are given in

Table 4. Estimated market prices of risk (MPR).

	DNK	EA	SWE
Constant MRP in the Vasicek model	−1.86103	−3.81415	−1.84068
Constant MPR in the shadow rate model	−1.68033	−2.90624	−1.67634
${\lambda }_1$ in the shadow rate model with switching MPR	−0.21014	−1.73524	−1.02982
${\lambda }_2$ in the shadow rate model with switching MPR	−3.38341	−3.80948	−2.12516

. The quality of fit in each case is shown in Figure 8, where the observed data of the long rates are compared with the fitted values. We can observe different behaviors: In the case of Denmark and Sweden, the Vasicek model and the shadow rate model with a constant market price of risk provide a very similar fit. However, allowing the market price of risk to switch its value lowers the objective function (and, therefore, makes the fit better) more pronouncedly in the case of Denmark. In the case of the Euro Area, the improvement is nontrivial both when changing from the Vasicek model to the shadow rate model with a constant market price of risk and when moving from a constant to a switching market price of risk in the shadow rate model.

8. Conclusions

In this paper, we studied one of the shadow rate models considered in the literature and provided a methodology for estimating its parameters. Firstly, the parameters from the time evolution of the short rate were estimated by means of the maximum likelihood method. The form of the model made it possible to reduce the maximization of the likelihood function to a one-dimensional optimization problem. In the second step, the market price of risk was estimated. We considered a constant market price or risk and a special form of switching market price of risk; for comparison, we also assumed the classical Vasicek model.

The first major implication of our study is a clear rejection of the Vasicek model in favor of the proposed shadow rate model based on the short rate evolution in all of the datasets. This approach uses only one new parameter, and the shadow rate can be reconstructed after its specification, which simplifies the use of the model. Therefore, it can be seen as a suitable alternative when modeling interest rates is only a part of a more complicated setting. In such a case, the interest rate model has to be kept sufficiently simple so that the whole analysis is trackable.

The second important result is the importance of a suitable specification of the market price of risk. In two out of three cases, the shadow rate model with a constant market price of risk was able to produce only practically the same fit of the long-term rates as the Vasicek model. Also, in two out of three datasets, a nonconstant market price of risk in the shadow rate model decreased the objective function by factors of 2.9 and 3.3, respectively, although the assumed form of the market price of risk was rather restricted by fixing some of its parameters and estimating only two constants. We consider this to be a very promising result, and we see a potential to improve the fit of the model by a proper choice of the market price of risk. This can be accomplished by introducing new parameters (in a generalization of our example, we may also be estimating the values for which the market price of risk switches its values) or by considering different classes of functions for its specification. Depending on the number of parameters and the behavior of the objective function, more consideration may be necessary regarding the choice of the algorithm used for its optimization. Furthermore, the optimization of the likelihood of the short rates and fitting the long rates do not need to be performed sequentially. While it simplifies the optimization problems, an equally interesting path of possible future research lies in a bicriteria optimization, i.e., looking for Pareto optimal parameters with respect to these two criteria. Furthermore, it is possible to look for a suitable market price of risk using the techniques for solving inverse problems instead of using a parametric form for its specification.

We estimated the model using data from 2011 to 2020, a period marked by negative short-term rates in its later years. Our results, therefore, offer insights into how interest rate behavior changes in such an environment. It can be, for example, noticed that the estimated equilibrium level of the shadow rate is negative across all cases. Obviously, this prevents the model from predicting subsequent growth of the interest rates. However, with sufficient amounts of “post-negative period” data, the model can be reestimated. In principle, a model with a positive equilibrium short rate can accommodate temporary negative rates before their reversion to the equilibrium and return to the positive region. During this period, the model would constrain the extent of negative rates. A practical evaluation of this concept can provide valuable new insights.

Finally, we note that the idea of a modified behavior when the short-term rate is negative can also be applied to any other model with negative values, which would take the role of the Vasicek process. It may be suitably modified also for the models with exogenous variables and model their different effects at the times with positive and negative short-term rates.