1 Introduction

Two prominent approaches to model the term structure of interest rates are the classes of equilibrium and no-arbitrage models. Most equilibrium models concentrate on the dynamic of the short-rate—the instantaneous interest rate—and derive interest rates with longer maturities from it. Prominent candidates of this model class include the models of Cox et al. [3], Duffie and Kan [9] and Vasicek [18]. No-arbitrage models focus on exactly fitting the term structure at a specific point in time to prevent arbitrage possibilities. Representatives of this class are introduced by Heath et al. [11] and Hull and White [12].

Applications of these models often relate to pricing interest rate derivatives, which is the reason why they are directly defined under the risk neutral measure most of the time. A general form of a one-factor short-rate model under the risk neutral measure is, e.g. given by

$$\begin{aligned} dr(t) = \mu (t,r)dt + \sigma (t,r)dW(t), \end{aligned}$$

where \(\mu \) and \(\sigma \) are two functions, which can depend on time point t and the short-rate r, and W is a Brownian motion. A lot of advances in theoretic models and their estimation have been conducted in the last 30 years, but only in connection to pricing (see Diebold and Li [6]). Regarding these models little attention has been given to forecasting and risk management purposes (see Diebold and Li [6]). For these applications the corresponding model needs to be regarded under the real world measure. Under this measure the corresponding one factor short-rate model has the following dynamic

$$\begin{aligned} dr(t) = \bigg [\mu (t,r) + \lambda (t,r) \sigma (t,r) \bigg ]dt + \sigma (t,r)d\widetilde{W}(t), \end{aligned}$$

where \(\lambda \) is the market price of risk and can also depend on t and r. \(\widetilde{W}\) is a Brownian motion under the real world measure. The exact functional choice for \(\lambda \) completes the model specification under the real world measure. Dai and Singleton [5] as well as Jong [14] use a fixed multiple of the model’s variance for the market price of risk and investigate the in sample fit of specific short-rate models, but do not focus on forecasting. Duffee [8] concludes that the class of term structure models analysed in Dai and Singleton [5] fail in forecasting. He argues that a restriction for the market price of risk to be a fixed multiple of the variance reduces the flexibility of the model. Hull et al. [13] stress that the market price of risk for a model with few factors should be time-dependent. This results not from an economic interpretation but from a modelling issue because of an insufficient number of factors (see Hull et al. [13]). They estimated the market price of risk based on historical 3-month and 6-month interest rates and came to a similar result as Ahmad and Wilmott [1], Cox and Pedersen [4] and Stanton [17]. But they argue that this value is only valid in the short horizon. Keeping this market price of risk constant could lead to extreme risk premiums and interest rates in the long horizon.

In this paper we tackle exactly this problem for the Gauss2++ model. Instead of assuming a constant, we assume a time-varying function for the market price of risk. In contrast to Hull et al. [13], who estimate the market price of risk for each forecasting horizon individually, we propose two parametric functions. The step function is the easiest non-constant function, which allows to model a market price of risk valid in the short and one valid in the long horizon. The linear function assumes that the market price of risk in the short horizon converges linearly to a long-term level. With these simplified time-dependent functions it is possible to account for the problem mentioned by Hull et al. [13] and the functions can still be easily estimated by historical data or calibrated in a forward looking manner to interest rate forecasts.

In our backtest we use a very similar calibration approach as described in Korn and Wagner [15]. The framework illustrated in this monograph has been developed by the Fraunhofer ITWM on behalf of the Produktinformationsstelle Altersvorsorge GmbH (PIA) and is the industry standard to classify packaged retail and insurance based investment products (PRIIPs) into chance-risk classes. For the interest rate model they use a Gauss2++ model with a presumed constant market price of risk. Following their calibration procedure allows us to compare our results to real applications in the insurance industry.

The structure of the paper is as follows. In Sect. 2 we introduce the Gauss2++ model under the risk neutral and the real world measure in a very general framework. In Sect. 3 we propose the constant function for comparison reason as well as the step and the linear function to specify the change of measure and explain how they can be estimated. All three variants of the Gauss2++ model are applied to data and backtested for the last 3 years in Sect. 4. In the final section the results are summarized and concluded.

2 The Gauss2++ model in the risk neutral and the real world

Throughout this section a filtered probability space \((\Omega , \mathcal {F}, (\mathcal {F}_t)_{t \in [0,\mathcal {T}]}, \mathbb {M})\) is given, where \(\mathbb {M}\) is either the risk neutral measure \(\mathbb {Q}\) with respect to the bank account or the real world measure \(\mathbb {P}\). \(\mathcal {T}\) represents an appropriate modelling horizon. The bank account \((B(t))_{t \in [0,\mathcal {T}]}\) is given by

$$\begin{aligned} dB(t) = r(t)B(t)dt, \qquad B(0) = 1, \end{aligned}$$

where r(t) denotes the short-rate. We further adopted notations and descriptions of the Gauss2++ model from the relevant chapters in Brigo and Mercurio [2].

2.1 The Gauss2++ model under the risk neutral measure

Short-rate models differ in the underlying process for the short-rate. The Gauss2++ model assumes that the short-rate is given by a sum of two correlated normally distributed processes, \((x(t))_{t\in [0,\mathcal {T}]}\) and \((y(t))_{t\in [0,\mathcal {T}]}\), and a deterministic function \(\varphi \), which is well defined on the time interval \([0,\mathcal {T}]\):

$$\begin{aligned} r(t) = x(t) + y(t) + \varphi (t), \qquad r(0) = r_0, \end{aligned}$$

where \(r_0\) is the short-rate at time point 0. The processes \((x(t))_{t\in [0,\mathcal {T}]}\) and \((y(t))_{t\in [0,\mathcal {T}]}\) satisfy under the risk neutral measure \(\mathbb {Q}\) the following stochastic differential equations

$$\begin{aligned} dx(t)&= -ax(t)dt + \sigma dW^1(t), \qquad x(0) = 0,\\ dy(t)&= -by(t)dt + \eta dW^2(t), \qquad y(0) = 0,\\ \rho dt&= dW^1(t)dW^2(t), \end{aligned}$$

where a, b, \(\sigma \), \(\eta \) are non-negative constants and \(-1 \le \rho \le 1\) is the instantaneous correlation between the two Brownian motions \(W^1\) and \(W^2\).

Short-rate models derive spot rates via prices of zero-coupon bonds. As the short-rate in the Gauss2++ model is normally distributed, there exists an analytic solution for a zero-coupon bond price, P(t, T), at time point t and maturity T:

$$\begin{aligned} P(t,T) = e^{-\int _t^T \varphi (s) ds - B(a,t,T) x(t) - B(b,t,T) y(t) + \frac{1}{2} V(t,T)}, \end{aligned}$$
(1)

where

$$\begin{aligned} B(z,t,T) = \frac{1-e^{-z(T-t)}}{z} \end{aligned}$$

and

$$\begin{aligned} V(t,T)&= \frac{\sigma ^2}{a^2} \left[ (T-t) + \frac{2}{a}e^{-a(T-t)} - \frac{1}{2a}e^{-2a(T-t)} - \frac{3}{2a} \right] \nonumber \\&\quad + \frac{\eta ^2}{b^2} \left[ (T-t) + \frac{2}{b}e^{-b(T-t)} - \frac{1}{2b}e^{-2b(T-t)} - \frac{3}{2b} \right] \nonumber \\&\quad + 2\rho \frac{\sigma \eta }{ab} \left[ (T-t) + \frac{e^{-a(T-t)}-1}{a} + \frac{e^{-b(T-t)}-1}{b} - \frac{e^{-(a+b)(T-t)}-1}{a+b} \right] . \end{aligned}$$

A derivation can be found in Brigo and Mercurio [2]. With formula (1) for the zero-coupon bond price under the risk neutral measure spot rates can be directly derived via

$$\begin{aligned} r(t,T) = \frac{-ln(P(t,T))}{T-t}, \end{aligned}$$
(2)

where r(t, T) represents the spot rate at time point t and a maturity of T.

The financial market we actually model consists of a bank account and a set of zero-coupon bonds, P(t, T), which differ in the maturity T. The dynamic of a zero-coupon bond price can be derived from the bond price formula in (1) by applying Ito’s formula and is given by

$$\begin{aligned} dP(t,T) = P(t,T) \bigg [r(t)dt - \sigma B(a,t,T) dW^1(t) - \eta B(b,t,T)dW^2(t)\bigg ]. \end{aligned}$$

A detailed derivation can be found in “Appendix 1”. Note that all assets have the same drift as it is the case in the risk neutral world.

2.2 The Gauss2++ model under the real world measure

To calculate performance scenarios and risk indicators the Gauss2++ model must be regarded under the real world measure \(\mathbb {P}\).

2.2.1 The change of measure

By specifying the Gauss2++ model under the risk neutral measure, we implicitly assume an arbitrage free market. Therefore, we can make the transition to a real world measure \(\mathbb {P}\) by defining the change of measure according to Girsanov, who states that a progressive and square-integrable process \(\varvec{\Phi }=\left( \Phi ^1(t), \Phi ^2(t),\ldots ,\Phi ^d(t)\right) _{t\in [0,\mathcal {T}]}\) determines a new probability measure \(\mathbb {P}\) such that if \((\varvec{\widehat{W}}(t))_{t\in [0,\mathcal {T}]}\) is a standard d-dimensional \((\mathcal {F}_t)_{t \in [0,\mathcal {T}]}\)-Brownian motion under \(\mathbb {Q}\), then

$$\begin{aligned} \varvec{\breve{W}}(t) := \varvec{\widehat{W}}(t) + \int _0^t \varvec{\Phi }(s) ds \end{aligned}$$

defines a standard d-dimensional \((\mathcal {F}_t)_{t \in [0,\mathcal {T}]}\)-Brownian motion under \(\mathbb {P}\) (see Girsanov [10]).

The Gauss2++ model is a two-factor model and \(\varvec{\Phi }\) is therefore two-dimensional. Its components can be interpreted as the market price of risk for each factor in the model. We will represent \(\varvec{\Phi }\) such that the resulting processes \((x(t))_{t\in [0,\mathcal {T}]}\) and \((y(t))_{t\in [0,\mathcal {T}]}\) still belong to the class of Ornstein–Uhlenbeck processes under \(\mathbb {P}\)

$$\begin{aligned} \varvec{\Phi }(t) = \left( \begin{array}{c} \Phi ^1(t) \\ \Phi ^2(t) \end{array} \right) = \left( \begin{array}{c} -\frac{ad_x(t)}{\sigma } \\ -\frac{bd_y(t)}{\eta \sqrt{1-\rho ^2}} + \frac{\rho ad_x(t)}{\sigma \sqrt{1-\rho ^2}} \end{array} \right) . \end{aligned}$$
(3)

Note that we restrict \(\varvec{\Phi }\) to be a function of time. By this the change of measure only changes the mean reversion level. More general measure change specification can be applied. For example, Diez and Korn [7] introduce a measure change for the 1-Factor Vasicek model, which influences the mean reversion level as well as the mean reversion speed.

The conditions for the Girsanov theorem translate directly to the functions \(d_x(t)\) and \(d_y(t)\). In the following we will specify the change of measure via \(d_x(t)\) and \(d_y(t)\). An appropriate interpretation of these functions will be given in Sect. 2.2.2.

2.2.2 The dynamics under the real world measure \(\varvec{\mathbb {P}}\)

With the representation of \(\varvec{\Phi }\) as in (3) the dynamics of the processes x and y in the Gauss2++ model change according to Girsanov to

$$\begin{aligned} dx(t)= &\, {} a(d_x(t)-x(t))dt + \sigma d\widetilde{W}^1(t), \qquad x(0) = 0, \end{aligned}$$
(4)
$$\begin{aligned} dy(t)= &\, {} b(d_y(t)-y(t))dt + \eta d\widetilde{W}^2(t), \qquad y(0) = 0, \end{aligned}$$
(5)

where \(\widetilde{W}^1\) and \(\widetilde{W}^2\) are two correlated Brownian motions under \(\mathbb {P}\). The derivation can be found in “Appendix 2”. We observe that x and y are still Ornstein–Uhlenbeck processes with the solutions

$$\begin{aligned} x(t)&= \int _0^t e^{-a(t-u)}ad_x(u)du + \sigma \int _0^t e^{-a(t-u)} d\widetilde{W}(u), \end{aligned}$$
(6)
$$\begin{aligned} y(t)&= \int _0^t e^{-b(t-u)}bd_y(u)du + \eta \int _0^t e^{-b(t-u)} d\widetilde{W}(u). \end{aligned}$$
(7)

The mean reversion level of each process at time point t amounts to \(d_x(t)\) and \(d_y(t)\), respectively. Recall that the sum of x(t) and y(t) and a deterministic function \(\varphi (t)\) under the risk neutral measure adds up to the instantaneous return rate r(t) of a risk free investment. Changing the measure changes the mean reversion level at time point t from 0 to \(d_x(t)\) for the process x and to \(d_y(t)\) for the process y. Therefore, \(d_x(t)+d_y(t)\) can be interpreted as the local long run risk premium of the short-rate—the amount, which is added in the real world to the risk neutral short-rate in the long run, if \(d_x(t)+d_y(t)\) would stay constant over time. If this amount is negative, future bond prices increase in expectation compared to the risk neutral world and a risk averse investor, therefore, gets compensated for the risk of investing in a risky bond. This means in contrast to equity prices, in a market where investors are risk averse, future interest rates tend to be lower in the real world than in the risk neutral world (see, e.g. Hull et al. [13]). Therefore, \(d_x(t)\) and \(d_y(t)\) can be interpreted as the local long run risk premium the corresponding risk factor is mean reverting to at time point t.

In the following we will specify the change of measure by these two functions instead of the market prices of risk. The market price of risk of each risk factor is then directly defined by these two functions.

$$\begin{aligned}&\text {Market price of risk of risk factor 1:} -\frac{ad_x(t)}{\sigma } \\&\text {Market price of risk of risk factor 2:} -\frac{bd_y(t)}{\eta \sqrt{1-\rho ^2}} + \frac{\rho ad_x(t)}{\sigma \sqrt{1-\rho ^2}} . \end{aligned}$$

If we assume a step or a piecewise linear function for \(d_x(t)\) and \(d_y(t)\) the functional form of the individual market prices of risk are the same.

The dynamics of a zero-coupon bond with maturity T under \(\mathbb {P}\) has the following form

$$\begin{aligned} dP(t,T) =&P(t,T)\left[ r(t) - B(a,t,T)ad_x(t) - B(b,t,T)bd_y(t)\right] dt \nonumber \\&- P(t,T)B(a,t,T)\sigma d\widetilde{W}^1(t) - P(t,T)B(b,t,T)\eta d\widetilde{W}^2(t) \end{aligned}$$
(8)

The derivation can be found in “Appendix 3”.

2.2.3 The bond price formula under the real world measure

The price of a zero-coupon bond under \(\mathbb {P}\) is obtained by the same analytic formula as in (1). The only difference is that the x- and the y-process are now regarded under the real world measure (see Diez and Korn [7]). In the following we will shortly explain why the formula does not change under this new measure.

To calculate the price of a zero-coupon bond under the real world measure we use the following conditional expectation

$$\begin{aligned} \frac{P(t,T)}{X_{P(t,T)}(t)} = E^{\mathbb {P}}\left[ \frac{P(T,T)}{X_{P(t,T)}(T)} \bigg \vert \mathcal {F}_t\right] , \end{aligned}$$

where \(X_{P(t,T)}\) represents the cash flow, with which we have to discount the zero-coupon bond such that the discounted price process is a martingale under \(\mathbb {P}\). The dynamic of \(X_{P(t,T)}\) coincides with the deterministic part of the zero-coupon bond price dynamic in (8) and is therefore specified by the change of measure:

$$\begin{aligned} dX_{P(t,T)}(t) = X_{P(t,T)}(t)\left[ r(t) - B(a,t,T)ad_x(t) - B(b,t,T)bd_y(t)\right] dt, \qquad X_{P(t,T)}(0) = 1. \end{aligned}$$

A short proof can be found in “Appendix 4”. The solution of this dynamic is given by

$$\begin{aligned} X_{P(t,T)}(t) = e^{\int _0^t \left( r(u) - B(a,u,T) a d_x(u) - B(b,u,T) bd_y(u) \right) du}. \end{aligned}$$

The price of a zero-coupon bond at time point t is therefore given by

$$\begin{aligned} P(t,T) = E^{\mathbb {P}} \left[ \frac{X_{P(t,T)}(t)}{X_{P(t,T)}(T)} \bigg \vert \mathcal {F}_t \right] . \end{aligned}$$

The ratio in the expectation amounts to

$$\begin{aligned} \frac{X_{P(t,T)}(t)}{X_{P(t,T)}(T)} = e^{-\int _t^T \left( r(u) - B(a,u,T) a d_x(u) - B(b,u,T) bd_y(u) \right) du}. \end{aligned}$$

To determine the distribution of this ratio, we first derive the distribution of the integral in the exponent, i.e.

$$\begin{aligned} I(t,T) {:}{=}\int _t^T \left( r(u) - B(a,u,T) a d_x(u) - B(b,u,T) bd_y(u) \right) du. \end{aligned}$$

It can be shown that I(t, T) is normally distributed with mean

$$\begin{aligned} M(t,T) = \int _t^T \varphi (u)du + \frac{1-e^{-a(T-t)}}{a} x(t) + \frac{1-e^{-b(T-t)}}{b} y(t) \end{aligned}$$
(9)

and variance

$$\begin{aligned} V(t,T)&= \frac{\sigma ^2}{a^2} \left[ (T-t) + \frac{2}{a}e^{-a(T-t)} - \frac{1}{2a}e^{-2a(T-t)} - \frac{3}{2a} \right] \nonumber \\&\quad + \frac{\eta ^2}{b^2} \left[ (T-t) + \frac{2}{b}e^{-b(T-t)} - \frac{1}{2b}e^{-2b(T-t)} - \frac{3}{2b} \right] \nonumber \\&\quad + 2\rho \frac{\sigma \eta }{ab} \left[ (T-t) + \frac{e^{-a(T-t)}-1}{a} + \frac{e^{-b(T-t)}-1}{b} - \frac{e^{-(a+b)(T-t)}-1}{a+b} \right] . \end{aligned}$$
(10)

The variance is the same as in the risk neutral world as the change of measure does not influence the variance of the processes. Note that also the mean has the same form as in the risk neutral case as the terms \(B(a,u,T)ad_x(u)\) and \(B(b,u,T)bd_y(u)\) in I(t, T) cancel out in the calculations. The derivations can be found in “Appendix 5”.

The expression \(e^{-I(t,T)}\) is therefore log-normally distributed and the zero-coupon bond price under \(\mathbb {P}\) is given by the same analytic formula as under \(\mathbb {Q}\):

$$\begin{aligned} P(t,T)&= E^{\mathbb {P}}\left[ e^{-\int _t^T r(u) - B(a,u,T) a d_x(u) - B(b,u,T) bd_y(u) du} \mid \mathcal {F}_t\right] \\&= e^{-M(t,T) + \frac{1}{2} V(t,T)} \\&= e^{-\int _t^T \varphi (u)du - \frac{1-e^{-a(T-t)}}{a} x(t) - \frac{1-e^{-b(T-t)}}{b} y(t) + \frac{1}{2}V(t,T)}. \end{aligned}$$

3 Local long run risk premium functions—specification and calibration

In the following three different types of functions for \(d_x(t)\) and \(d_y(t)\) are introduced: the constant, the step and the linear function. Following the interpretation in Sect. 2.2.2 these functions represent the long run risk premium for each risk factor at a specific time point t in the Gauss2++ model. The functional equations of the three types are

$$\begin {aligned} {\text{Constant:}}\quad d_x(t) &= d_x \\ d_y(t) &= d_y \\ {\text{Step:}}\quad d_x(t) &= \mathbb{1}_{t\leq\tau}d_x + \mathbb{1}_{t>\tau}l_x \\ d_y(t) &= \mathbb {1}_{t\le \tau }d_y + \mathbb {1}_{t>\tau }l_y \\ {\text{Linear:}}\quad d_x(t) &= \mathbb {1}_{t\le \tau }(1-m_xt)d_x + \mathbb {1}_{t>\tau }l_x \\ d_y(t) &= \mathbb {1}_{t\le \tau }(1-m_yt)d_y + \mathbb {1}_{t>\tau }l_y\end {aligned}$$

where \(d_x\), \(l_x\), \(m_x\) and \(d_y\), \(l_y\), \(m_y\) are real valued constants and \(\mathbb {1}_{A}\) represents the indicator function of a subset A.

The constant function assumes that the local long run risk premium is constant for the whole modelling horizon. The latter two functions distinguish between a local long run risk premium valid in the short and in the long horizon, seperated at time point \(\tau \). As mentioned in Sect. 2.2.2 the same holds for the market price of risk, respectively. Hull et al. [13] argue that a time-varying market price of risk is necessary to account for unobserved risk factors and to prevent unrealistic interest rate forecasts in the long horizon. They therefore estimate an individual market price of risk for each forecasting horizon. We use a more parsimonious function with regard to the number of parameters. The step function we propose is the simplest time-varying function that expects that the local long run risk premium differs in the short and the long horizon but is still constant in each period. The linear function implements the property that the local long run risk premium in the short horizon approaches the long term level linearly. The simplicity of these functions allows a straight forward calibration to interest rate forecasts.

Because of the distributional properties of the Gauss2++ model the expected values for interest rates under the real world measure \(\mathbb {P}\) for any future time point can be calculated:

$$\begin{aligned} E^{\mathbb {P}}[r(t,T)]&= E^{\mathbb {Q}}[r(t, T)] + \frac{B(a,t,T)}{T-t} RP_x(t) + \frac{B(b,t,T)}{T-t} RP_y(t), \end{aligned}$$
(11)

where \(RP_x(t)\) and \(RP_y(t)\) represent the actual risk premium of the short-rate at time point t for each risk factor and are given by the first integral in (6) and (7)

$$\begin{aligned} RP_x(t)&{:}{=}\int _0^t e^{-a(t-u)}ad_x(u)du, \\ RP_y(t)&{:}{=}\int _0^t e^{-b(t-u)}bd_y(u)du. \end{aligned}$$

For the constant, the step and the linear function these integrals can be easily calculated. To get the risk premium for longer maturities the functions \(RP_x(t)\) and \(RP_y(t)\) are weighted by a loading function, which accounts for the different riskiness of the corresponding zero-coupon bonds

$$\begin{aligned} \frac{B(a,t,T)}{T-t} \qquad \text {and} \qquad \frac{B(b,t,T)}{T-t}. \end{aligned}$$

To calibrate the local long run risk premium functions, \(d_x(t)\) and \(d_y(t)\), the parameters of the functions are chosen in such a way that the model meets specific interest rate forecasts in expectation. For the constant type two interest rate forecasts are needed. For the other two types four interest rate forecasts are necessary—two short term and two long term forecasts. The time parameter \(\tau \), which determines the separation between the short and the long term local long run risk premium must lie between the forecasting horizons of the two short and the two long term forecasts.

In Fig. 1 the three types of local long run risk premium functions have been exemplary calibrated. \(\tau \) has been set to 24 months, which is the forecasting horizon of the short term interest rate forecasts.

Fig. 1
figure 1

Local long run risk premium functions

Full size image

In the following subsections the calibration procedures for all three types of local long run risk premium functions, which are applied in this paper, are described.

3.1 The constant function

The constant functions represented in Fig. 1a implement a constant local long run risk premium for the whole modelling horizon, which can amount to 40 years or more for actual applications in the insurance industry, e.g. to classify certified pension contracts into risk classes. The absolute risk premiums, \(RP_x(t)\) and \(RP_y(t)\), are given by:

$$\begin{aligned} RP_x(t)&= (1-e^{-at})d_x, \\ RP_y(t)&= (1-e^{-bt})d_y. \end{aligned}$$

Note that if \(t \rightarrow \infty \), \(RP_x(t)\) and \(RP_y(t)\) indeed converge to \(d_x\) and \(d_y\), the long run risk premiums, respectively. To calibrate the parameters of the constant functions two interest rate forecasts, \( \hat{r}(t_1, T_1) \) and \(\hat{r}(t_2, T_2)\), are used. Plugging the absolute risk premium functions, \(RP_x(t)\) and \(RP_y(t)\), into (11) and setting the expectations equal to the interest rate forecasts results in the following two equations

$$\begin{aligned} \begin{array}{c l} \hbox {(I)}&{} \hat{r}(t_1, T_1) \overset{!}{=} E^Q[r(t_1, T_1)] + \frac{B(a,t_1,T_1)}{(T_1-t_1)} (1-e^{-at_1})d_x + \frac{B(b,t_1,T_1)}{(T_1-t_1)} (1-e^{-bt_1})d_y ,\\ \hbox {(II)}&{} \hat{r}(t_2, T_2) \overset{!}{=} E^Q[r(t_2, T_2)] + \frac{B(a,t_2,T_2)}{(T_2-t_2)} (1-e^{-at_2})d_x + \frac{B(b,t_2,T_2)}{(T_2-t_2)} (1-e^{-bt_2})d_y . \end{array} \end{aligned}$$

As the expectations are linear functions in \(d_x\) and \(d_y\), the two parameters can be easily determined.

The constant function for the local long run risk premium in the Gauss2++ model and this calibration procedure is a standard approach in the insurance industry. As the values for \(d_x\) and \(d_y\) determine the risk premium for the whole modelling horizon, their calibration is crucial for the model’s interest rate distribution. Especially if the interest rate forecasts used for the calibration have a short forecasting horizon, the resulting distribution in the long horizon is very sensitive to these forecasts. For example if the interest rate forecasts and the forward rates—calculated from the current yield curve—are very different, to reach the forecasts a huge risk premium is necessary, which might be valid in the short horizon, but produces extreme interest rates in the long horizon. The next two functions account for this problem by representing a time-varying local long run risk premium.

3.2 The step function

The step functions represented in Fig. 1b take the same value as the corresponding constant function up to time \(\tau \) as the same interest rate forecasts have been used for the short horizon, but then they jump to a different level to account for the risk premium in the long horizon. Similar to the constant function the absolute risk premium functions can easily be calculated and amount to

$$\begin{aligned} RP_x(t)&= \left( e^{-a(t-\min (t,\tau ))}-e^{-at}\right) d_x + \left( 1-e^{-a (t - \min (t,\tau ))}\right) l_x, \\ RP_y(t)&= \left( e^{-b(t-\min (t,\tau ))}-e^{-bt}\right) d_y + \left( 1-e^{-b (t- \min (t,\tau ))}\right) l_y. \end{aligned}$$

Note that if \(t \rightarrow \infty \), \(RP_x(t)\) and \(RP_y(t)\) now converge to \(l_x\) and \(l_y\), respectively. To calibrate the four parameters of the step function two short term and two long term interest rate forecasts are used resulting in the following equations:

$$\begin{aligned} \begin{array}{c l} \hbox {(I)}&{} \hat{r}(t_1, T_1) \overset{!}{=} E^Q[r(t_1, T_1)] + \frac{B(a,t_1,T_1)}{(T_1-t_1)} RP_x(t_1) + \frac{B(b,t_1,T_1)}{(T_1-t_1)} RP_y(t_1) ,\\ \hbox {(II)}&{} \hat{r}(t_2, T_2) \overset{!}{=} E^Q[r(t_2, T_2)] + \frac{B(a,t_2,T_2)}{(T_2-t_2)} RP_x(t_2) +\frac{B(b,t_2,T_2)}{(T_2-t_2)} RP_y(t_2) , \\ \hbox {(III)}&{} \hat{r}(t_3, T_3) \overset{!}{=} E^Q[r(t_3, T_3)] + \frac{B(a,t_3,T_3)}{(T_3-t_3)} RP_x(t_3) + \frac{B(b,t_3,T_3)}{(T_3-t_3)} RP_y(t_3) ,\\ \hbox {(IV)}&{} \hat{r}(t_4, T_4) \overset{!}{=} E^Q[r(t_4, T_4)] + \frac{B(a,t_4,T_4)}{(T_4-t_4)} RP_x(t_4) + \frac{B(b,t_4,T_4)}{(T_4-t_4)} RP_y(t_4) , \end{array} \end{aligned}$$

where \(t_1 \le t_2 < t_3 \le t_4 \). \(\tau \) must lie between \(t_2\) and \(t_3\), i.e. \(t_2 \le \tau < t_3\).

Instead of interest rate forecasts direct forecasts of the absolute risk premium of the short-rate can be used. This approach is applied by Hull et al. [13], who estimate risk premiums for each forecasting horizon from historical data, but they also scale their result to a long term short-rate forecast. Another possible approach is to take the ultimate forward rate (UFR) from Solvency II as a long term target, which is reached at a future time point with a specific percentage (e.g. \(95\%\) of the UFR in 40 years) and to \(100\%\) in the limit, i.e. \(t \rightarrow \infty \).

3.3 The linear function

The linear functions represented in Fig. 1c avoid the sudden jump as it is the case in the step functions and converge in the short term linearly to a long term level. The absolute risk premiums at time point t can be calculated as before and amount to

$$\begin{aligned} RP_x(t) &= \left ( \left(e^{-a(t-\min(t,\tau))}-e^{-at}\right)\left(1+\frac{m_x}{a}\right) - e^{-a (t-\min(t,\tau))}m_x\min(t,\tau)\right) d_x \\ &+ \left(1-e^{-a (t-\min(t,\tau))}\right)l_x, \\ RP_y(t) &= \left(\left(e^{-b(t-\min(t,\tau))}-e^{-bt}\right)\left(1+\frac{m_y}{b}\right) - e^{-b (t-\min(t,\tau))}m_y\min(t,\tau)\right)d_y \\ &+ \left(1-e^{-b (t-\min(t,\tau))}\right)l_y. \end{aligned}$$

Note again that if \(t \rightarrow \infty \), \(RP_x(t)\) and \(RP_y(t)\) converge to \(l_x\) and \(l_y\), the long term risk premiums, respectively. To calibrate \(d_x\), \(l_x\), \(d_y\) and \(l_y\) four interest rate forecasts as for the step function are used. By imposing that the absolute risk premium functions, \(RP_x(t)\) and \(RP_y(t)\), are differentiable at the forecasting horizon \(\tau \) to prevent a kink in the absolute risk premium function, two further conditions are incorporated to specify \(m_x\) and \(m_y\):

$$\begin{aligned} \begin{array}{c l} \hbox {(V)}&{} (RP_x)'_-(\tau ) = (RP_x)'_+(\tau ) ,\\ \hbox {(VI)}&{} (RP_{y})'_-(\tau ) = (RP_{y})'_+(\tau ), \end{array} \end{aligned}$$

where \((RP_{\cdot })'_-(\tau )\) and \((RP_{\cdot })'_+(\tau )\) denote the derivative from the left and from the right, respectively. Solving the equations for \(m_x\) and \(m_y\) leads to the following closed form solutions reducing the number of free parameters to four:

$$\begin{aligned} m_x&= \frac{d_x - l_x}{d_x \tau }, \\ m_y&= \frac{d_y - l_y}{d_y \tau }. \end{aligned}$$

Note that with this condition the same number of interest rate forecasts as for the step function are needed to calibrate \(d_x(t)\) and \(d_y(t)\).

4 Results

In this section the calibration results of three variants of our framework for the Gauss2++ model are presented. The variants differ in the assumption about the local long run risk premium functions, which determine the change from the risk neutral to the real world measure. Variant 1 assumes a constant, variant 2 a step and variant 3 a linear local long run risk premium function for the risk factors. In the first subsection the three variants of the Gauss2++ model are compared if calibrated at the same valuation date. In Sect. 4.2 we show with a backtest over the last three years that variant 2 and 3 produce much more stable interest rate scenarios for the long forecasting horizon over this time period. This stability would transfer to performance scenarios and risk measures of, e.g. an interest rate sensitive fonds.

4.1 Calibration at one valuation date

The calibration process of the Gauss2++ model can be split into two steps. In the first step the model is calibrated under the risk neutral measure. This step does not depend on the choice of the local long run risk premium function and is therefore the same for all modelling cases. In the second step the change of measure is calibrated. The choice of the local long run risk premium function plays an important role and leads to different interest rate scenarios, performance measures and risk indicators.

To calibrate the model at a specific valuation date under the risk neutral measure the term structure of interest rate swaps and swaption volatilities at this date are used. In the Gauss2++ model the deterministic and time-dependent function \(\varphi \) ensures the market consistency regarding the current term structure by being defined as follows:

$$\begin{aligned} \varphi (t) = f^M(0,t) + \frac{\sigma ^2}{2a}(1-e^{-at})^2 + \frac{\eta ^2}{2b}(1-e^{-bt})^2 + \rho \frac{\sigma \eta }{ab}(1-e^{-at})(1-e^{bt}). \end{aligned}$$

\(f^M(0,t)\) represents the instantaneous forward rate at time point 0 for a maturity t, i.e. \(f^M(0,t) = \frac{\partial P^M(0,t)}{\partial T}\), where \(\frac{\partial P^M}{\partial T}\) denotes the partial derivative with respect to the second argument and \(P^M(0,t)\) is the market zero-coupon bond price. For the derivation and further information the reader is referred to Brigo and Mercurio [2]. The parameters \(a,b,\sigma ,\eta \) and \(\rho \) of the model are chosen in such a way that the model prices of the considered swaptions coincide with the market prices. For this the downhill simplex algorithmFootnote 1 is used to minimize the root mean squared error (RMSE):

$$\begin{aligned} RMSE = \sqrt{\sum _{i=1}^N \left( C_{Model,i}(a,b,\sigma , \eta , \rho ) - C_{Market,i}\right) ^2}, \end{aligned}$$

where \(C_{Model,i}\) represents the model price of swaption i of the Gauss2++ model and \(C_{Market,i}\) is the market price of that swaption. The swaptions considered in the calibration process differ with respect to their tenor and maturity combination, which is denoted by the subscript i. N represents the number of considered swaptions. The analytic formula for the price of a swaption in the Gauss2++ model can be found in Brigo and Mercurio [2]. Table 1 shows the result of a calibration at the 31.12.2019. We used at-the-money receiver swaptions with a maturity and tenor combination of \(\{5,7,10,12,15,20\} \times \{5,7,10,12,15,20\}\), i.e. in total \(N=36\) swaption prices. The RMSE amounts to 0.0619. In the optimization we further restricted \(\rho \) to lie between \(-\,1\) and 1 as well as all other parameters to be \(> 0\).

Table 1 Parameters of the Gauss2++ model calibrated at 31.12.2019
Full size table

These parameters together with the current interest rate curve determine the dynamics of the Gauss2++ model under the risk neutral measure.

In the second step the local long run risk premium functions, which determine the change of measure, are calibrated to interest rate forecasts as described in Sects. 3.1–3.3. For the short term interest rate forecasts we use forecasts published by the OECD for a 3-month and a 10-year interest rate. To take the OECD forecasts we have been inspired by the framework developed by the Frauenhofer ITWM on behalf of PIA to classify PRIIPs into chance-risk classes (see Korn and Wagner [15]). Their model represents the industry standard for PRIIP calculations. The latest OECD forecasts regarding the 31.12.2019 for the longest horizon, which is the fourth quarter of 2021, amount to \(-\, 0.4\%\) and \(0.4\%\), respectively.Footnote 2 For the long term interest rate forecasts, which are needed to calibrate the step and the linear function, we take the average of monthly 3-month and 10-year interest rates over the last 15 years also published by the OECD. This is a valid approach if interest rates follow a stationary process, because in this case historical data can be considered as a random sample from the corresponding interest rate distribution. Hull et al. [13] point out that this approach is questionable if monetary and fiscal policies are expected to be materially different from those in the past. Nevertheless any other model based on historical data would be questionable and the user of the model can alternatively provide personal estimates or an expert judgment. The historical average amounts to \(1.08\%\) for the 3-month and \(1.84\%\) for the 10-year interest rate and as we assume these forecasts to be a long run average we set the forecasting horizon to 40 years—the modelling horizon. We further set \(\tau \) to 24 months, which is the forecasting horizon of the short term OECD forecasts.

Table 2 shows the calibration results for the three local long run risk premium function types.

Table 2 Parameters of the local long run risk premium functions
Full size table

The values of \(d_x\) and \(d_y\) coincide for the constant and the step function as the same interest rate forecasts have been used in the calibration process. But in contrast to the step function, which takes the values of \(l_x\) and \(l_y\) after 24 months, the constant function stays constant for the whole modelling horizon. It also appears that the step and the linear function take the same values for \(l_x\) and \(l_y\). But there is a slight difference as their functional forms differ in the first two years, which influences the absolute risk premium in future time points. This influence decreases in time, such that the difference is negligible as we calibrated \(l_x\) and \(l_y\) to forecasts with an forecasting horizon of 40 years.

Figures 2, 3 and 4 visualize for the three calibrated variants of the Gauss2++ model the development of the expectation of the short-rate, the 10-year and the 20-year interest rate for forecasting horizons of up to 40 years. The solid line represents the expectation under the risk neutral measure, the dashed line shows the expected values under the real world measure.

Fig. 2
figure 2

Constant function: expected values of the short-rate, the 10-year and the 20-year interest rate under the risk neutral and the real world measure

Full size image
Fig. 3
figure 3

Step function: expected values of the short-rate, the 10-year and the 20-year interest rate under the risk neutral and the real world measure

Full size image
Fig. 4
figure 4

Linear function: expected values of the short-rate, the 10-year and the 20-year interest rate under the risk neutral and the real world measure

Full size image

For the variant of the Gauss2++ model, which uses the constant function as the local long run risk premium function, the expected real world interest rates lie above the risk neutral expectation. This means, that a risk seeking behaviour of the investors is assumed for the whole modelling period, because an investor accepts a lower expected return for a corresponding bond if the interest rates are expected to be higher in the real world compared to the risk neutral world. Ahmad and Wilmott [1] show that there have been time periods where investors seem to have historically behaved in this way. But in general investors are assumed to be risk averse and therefore interest rates should be lower in the real world than in the risk neutral world, which is an opposite behaviour to equity prices (see, e.g. Hull et al. [13]). For the other two variants of the Gauss2++ model the expected real world interest rates lie also above the risk neutral interest rates in the short horizon but below in the long horizon. This assumption of risk seeking behaviour in the short horizon stems from the quite high forecasts of the OECD for the short horizon, but it might be valid in the current market situation. In contrast to the constant case, which keeps this risk seeking behaviour assumption for the whole modelling horizon, in the long run the other two variants of the Gauss2++ model assume in this calibration a risk averse behaviour. The difference between the step and the linear function is only visible in the short horizon. While the step function has a kink in the expectation after \(\tau \) years, the linear function is smoother due to its condition that the derivative of the absolute risk premium function exists at this time point.

Furthermore, the absolute difference in the risk neutral and real world expectations decreases for interest rates with longer maturities. This results from the less variation of interest rates with longer maturities, which is an implicit model characteristic of the Gauss2++ model and is supported by historical data as well. A risk premium is therefore higher (less negative) for a risk averse and lower (less positive) for a risk seeking investor in an arbitrage free market.

Figure 5 shows the absolute risk premium functions of the short-rate for all three modelling types.

Fig. 5
figure 5

Absolute risk premium function for the variants of the Gauss2++ model

Full size image

It can be observed that for the constant and the step function the absolute risk premium is the same up to year 2. After that year the Gauss2++ variant with the step function has a kink in the absolute risk premium as the local long run risk premium changes to a different level, while the modelling case with the constant function continuous to approach the long term risk premium determined by the short term interest rate forecasts. The modelling case with the linear function results in a different risk premium for the first 2 years, but approaches—without a kink—the same long term risk premium as the step function.

All three functions intersect after 2 years as this is the forecasting horizon of the short term interest rate forecasts, which were used for the calibration. The absolute risk premium at this time point must be the same for all modelling cases such that the expected interested rates of the model coincide with the forecasts.

We further investigated the resulting yield curve shapes of the three variants of the Gauss2++ model. The variant, which uses a constant function, represents the industry standard regarding PRIIP calculations (see Korn and Wagner [15]). An unpleasent feature of this model is the unrealistic high frequency of inverse yield curves with growing time (see Diez and Korn [7]). In their paper they show that for the 2-Factor Vasicek model the number of inverse yield curves can be reduced by assuming a negative risk premium. The share of inverse yield curves in our calibration of the three variants were investigated in a simulation study. We simulated 10,000 yield curve paths with each calibrated model and counted the number of yield curves, which have a higher 1-year interest rate than a 30-year interest rate. The result is visualized in Fig. 6. We can see a similar behaviour as described in the paper of Diez and Korn [7]. The variant with the constant function, which has a positive risk premium over the modelling horizon, shows an unrealistic high share of inverse yield curves. The other two variants have a negative risk premium and decrease the number of inverse yield curves in the long run compared to the risk neutral case. Using the step or the linear function for the risk premium results therefore not only in more realistic interest rates but also in more realistic yield curve shapes in the long horizon.

Fig. 6
figure 6

The share of inverse yield curves for the Gauss2++ model under the risk neutral measure and under the real world measure using a constant, a step and a linear function for the market price of risk

Full size image

4.2 Backtest

In this subsection the different variants of the Gauss2++ model calibrated on a quarterly basis over the last 3 years are compared.

As in Sect. 4.1 interest rate swaps and swaption volatilities have been used for the risk neutral calibration of the Gauss2++ model. To calibrate the parameters of the local long run risk premium functions in the second calibration step short term interest rate forecasts published by the OECD and a long term average have been used. The forecasts are shown in Table 3. The calibration results of the parameters of the Gauss2++ model under the risk neutral measure and of the local long run risk premium function for each variant of the Gauss2++ model can be found in Tables 4, 5, 6 and 7 in “Appendix 6”.

Table 3 Interest rate forecasts of the OECD and historical average of the 3-month and the 10-year interest rate
Full size table

For each calibration the absolut risk premium function of the short-rate and the development of the expected 10-year interest rate have been calculated and visualised in Figs. 7 and 8.

Fig. 7
figure 7

Absolute risk premium functions

Full size image
Fig. 8
figure 8

Development of the expectation of the 10-year interest rate over the modelling horizon for all three variants of the Gauss2++ model

Full size image

The absolute risk premium function of the short-rate for the Gauss2++ model, which uses the constant function for the local long run risk premium, depends highly on the risk neutral calibration results and the forecasts of the OECD. An unfavorable combination of market data and interest rate forecasts can lead to a high value for the local long run risk premium. This value might be reasonable to meet the short term forecasts used for the calibration, but as it stays constant over time it is the value the absolute risk premium is converging to. Therefore, this problem can strike through if the modelling horizon is much longer than the forecasting horizon of the interest rates used for the calibration. In this case a time-varying local long run risk premium function, which can be calibrated to a short and a long term forecast, is more convenient to regularize the risk premium. As it can be seen in Fig. 7 the variants of the Gauss2++ model, which use the step or the linear function for the local long run risk premium, produce more stable risk premiums in the long horizon. In each calibration the absolute risk premium is positive in the first years, which presumes a risk seeking behaviour of the investors, but in the long horizon the absolute risk premium lies between \(-\, 0.5\) and \(-2.5\%\) representing a risk averse market. Also the interest rate distribution in the long horizon is more stable. Figure 8b, c show that the expectation of the 10-year interest rate in the long horizon change only little in each calibration according to the historical average, which was used for the long term interest rate forecast.

5 Conclusion

As the Gauss2++ model is often used for pricing purposes, the focus in the literature lies on the evolution of interest rates under the risk neutral measure \(\mathbb {Q}\). But regarding risk management and forecasting applications the model under the real world measure is needed. In this paper we introduced a framework to apply the model under both measures in a consistent manner. This framework first conducts a calibration under the risk neutral measure and then determines the change of measure such that it is possible to switch between the risk neutral and the real world. We showed that according to Girsanov this change of measure can be specified by any time-dependent function without loosing the analytic tractability of, e.g. zero-coupon bond prices. Hull et al. [13] argue that because of unobserved risk factors, which are not included in the model, a time-varying function should be used, because otherwise unrealistic interest rates in the long forecasting horizon could be reached. We therefore compared the industry standard, which uses a constant function to model the change of measure, with two variants, which use either a step or a linear function. These functions are the simplest extensions of the constant function to a time-varying function without increasing the computational effort much. By accounting for different risk premiums in the short and in the long horizon the time-varying functions result in much more stable interest rate forecasts in the long run if calibrated at different valuation dates. From a macroeconomical point of view it makes sense that current market fluctuations should not influence interest rate forecasts in the long horizon, e.g. in 40 years, much. This would also imply that risk measures calculated with the Gauss2++ model, which uses one of the time-varying functions for the change of measure, would be more consistent if estimated at different valuation time points.

We further investigated the yield curve shapes by conducting a simulation study. The result is in line with the findings of Diez and Korn [7] for the 2-Factor Vasicek model. Assuming a positive risk premium—as it was the case in our calibration for the constant function—the number of inverse yield curves increases compared to the risk neutral case. This also replicates the problem of too many inverse yield curves in the insurance industry for PRIIP calculations (see Diez and Korn [7]). The other two variants represented in this paper, which apply a time-varying function for the market price of risk, assume a negative risk premium in the long run and have a much lower amount of inverse yield curves. Using a step or a linear function for the market price of risk, therefore, not only leads to more realistic interest rates in the long run, but also creates more realistic yield curve shapes.