Data driven modeling of multiple interest rates with generalized Vasicek-type models

Data driven modeling of multiple interest rates with generalized Vasicek-type models

\address{Aalto University School of Science, Department of Mathematics and Systems Analysis, PO Box 11100, 00076 Aalto, Finland} \email{[email removed]}

\address{Aalto University School of Science, Department of Mathematics and Systems Analysis, PO Box 11100, 00076 Aalto, Finland} \email{[email removed]}

\address{University of Vaasa, School of Technology and Innovations, P.O. Box 700, Vaasa, FIN-65101, Finland} \address{Taras Shevchenko National University of Kyiv, Department of Probability, Statistics and Actuarial Mathematics, Volodymyrska St., Kyiv, 01601, Ukraine}

\email{[email removed]}

\address{University of Vaasa, School of Technology and Innovations, P.O. Box 700, Vaasa, FIN-65101, Finland}

\email{[email removed]}

\address{Aalto University School of Business, Department of Information and Service Management, PO Box 11110, 00076 Aalto, Finland} \email{[email removed]}

\subjclass[2020]{91G30, 62M09, 60G10}

\allowdisplaybreaks

Introduction

Interest rates are important for economics as they affect the consumption of households and businesses. Moreover, they can be used as economic indicators since their values often change in accordance with economic shifts. However, the changes in interest rates are influenced by multiple different factors, and thus interest rates themselves are considered to be stochastic and can be modeled as random processes. The first model utilizing this idea was the one-factor short-rate interest model introduced by both Merton merton1973theory and Black and Scholes black1973pricing. In this model, the interest rate is described using a stochastic differential equation featuring a Brownian motion with drift. Building on this model, Vasicek vasicek1977equilibrium modeled interest rates as an Ornstein-Uhlenbeck process and included mean reversion in the model.

Inspired by black1973pricing, merton1973theory and vasicek1977equilibrium, many extensions and different versions of interest rate models have been proposed in the literature. Langetieg langetieg1980multivariate proposed a multivariate version of the Vasicek model with a vector of independent Brownian motions as the noise process. Cox, Ingersoll and Ross (CIR) cox1985theory used a similar form to Vasicek but defined the interest rate as a square-root process, disallowing negative values. The CKLS model chan1992empirical combines all these models into one, in which the parameters can be chosen to suit the data at hand. Hull and White hull1990pricing extended the Vasicek and CIR models to allow time-dependent parameters. This model is sometimes referred to as an extended Vasicek model.

Most generalizations of the Vasicek model are univariate. For example, the univariate fractional Vasicek model with univariate fractional Brownian motion as the noise process was introduced already in 1998 by Comte comte1998long and its statistical properties were studied in xiao2019frac1 and xiao2019frac2. Nourdin and Diu Tran nourdin2019statistical generalized the univariate fractional Vasicek model further by using Hermite processes as noise.

In this article, we study the generalized multivariate Vasicek model

where $X$ is a random process with stationary increments and $r$ is a vector of interest rates. This model is inspired by the vector-valued generalized Ornstein-Uhlenbeck process studied in Marko1, which can be obtained from (ref) by setting $b=0$ and $r_0 = \int_{-\infty}^0 e^{\Theta s}dX_s$. The generalized multivariate Vasicek model provides a framework to simultaneously model multiple interest rates. The proposed model is very flexible. As shown in Marko1, every stationary process with continuous sample paths can be viewed as a vector-valued generalized Ornstein-Uhlenbeck process. As such, the generalized Vasicek model (ref) covers all interest rates that are (almost) stationary around a mean level $b$.

The rest of the article is organized as follows. In Section (ref), we introduce the model and present our main results. In particular, we derive estimators for parameters and show that one obtains consistency and asymptotic distributions under rather minimal assumptions. We illustrate the practical applicability of our approach by presenting examples with simulated data in Section (ref) and a real data example in Section (ref). Technical proofs are presented in Section (ref).

Main results

We consider the generalized multivariate Vasicek model for interest rates $r^{\Theta, b} = \left(r^{\Theta, b}_1,\ldots,r^{\Theta, b}_d\right)$ given by

The initial value $r_0\in \mathbb{R}^d$ is assumed to be a square-integrable random vector. The vector $b\in \mathbb{R}^d$ represents the long-term means of the interest rates and $\Theta \in \mathbb{R}^{d\times d}$ is a positive definite matrix modeling the mean-reversion effect and interactions between the interest rates. The vector $X =(X_t)_{t\geq 0}\in \mathbb{R}^d$ is a square-integrable centered $d$-dimensional noise process with stationary increments with the vector $X_0=0$. The matrix $\sigma \in \mathbb{R}^{d\times d}$ is a positive definite diagonal matrix describing the size of the noise.

In the sequel, we assume that the vector $X$ has uncorrelated components and finite fourth moments. Moreover, we use the notation $\Vert \cdot \Vert$ for the $L^2$-norm. In addition, we drop the superscript from $r^{\Theta, b}$ and simply write $r$ whenever confusion cannot arise, and we denote the stationary solution corresponding to $b=0$ and $r_0 = \int_{-\infty}^0 e^{\Theta s}dX_s$ by $U=(U_t)_{t\geq 0}$, given by

where the integral can be defined via integration-by-parts, see Marko1. The covariance matrix of $U_t$ with time-lag $t$ is denoted by $\gamma(t) = \mathbb{E}[U_tU_0^\top] \in \mathbb{R}^{d\times d}$. Note that $\gamma(-t) = \gamma(t)^\top$. For the cross-covariance matrix $\operatorname{Cov}(X_t)$, we use the short notation $V(t)$.

As in Marko1, we use the notations

Note that we can express $\operatorname{Cov}(U_t-U_0)$ with $\gamma$ as $\operatorname{Cov}(U_t-U_0) = 2\gamma(0)-\gamma(t)-\gamma(-t)$. It was shown in Marko1 that the matrix $\Theta$ satisfies, for any $t\geq 0$, the continuous-time algebraic Riccati equation (CARE) $$ B_t^\top \Theta + \Theta B_t - \Theta C_t \Theta + D_t = 0. $$

In the sequel, we assume that a sample path of the process $r$ is observed continuously on the interval $[0,T]$. Our objective is to estimate the unknown parameters $\Theta$, $b$, and $\sigma$. We assume that the covariance matrix $V(t)$ of the process $X$ is known on $[0,t]$ for some fixed $t$. Note that this is not restrictive, as the uncertainty of the noise is encoded in the unknown $\sigma$.

Estimation of $\sigma$ requires densely spaced discrete observations over a fixed interval. From a theoretical perspective, under a continuous-time observation scheme, one effectively has infinitely many such observations. However, in practice, data are always available only in discrete form. For this reason, we prefer to retain the explicit dependence on the mesh size in our estimators; see Remark 2.6 for further details.

The estimators for $b$ and $\gamma$ depend on the length of the observed time series $T$ (or $T+s$), and the estimates get more accurate if a longer series is observed. Note that the estimators $\hat{B}_t,\hat{C}_t$, and $\hat{D}_t$ depend on a chosen time instance $t$. On a theoretical level, the choice of $t$ does not play a role as long as we know the value $V(t)$. However, in practice, the choice of $t$ affects the performance of the estimators. For small $t$, the matrices $B_t$ and $C_t$ are effectively zero matrices, and this may cause inaccuracies. For large $t$, the estimation of the matrices becomes computationally more expensive and, at the same time, more inaccurate since then one needs to estimate $\gamma(s)$ for larger lags $s$. The estimator of $\sigma \sigma^\top$ depends on the number of observations on a fixed time interval, and for simplicity, we use the $[0, 1]$ interval here. To get a better estimate of $\sigma \sigma^\top$, a larger number of observations on this interval is required so that the contribution of the drift term in (ref) vanishes. Therefore, a good estimate for $\Theta$ requires a long time series and densely observed data.

For obtaining asymptotical results, we make the following assumption.

The first part of Assumption (ref) (stating that the true covariance of the process $U$ decays to zero), together with Lemma (ref), allows us to consistently estimate $b$ from observations of $r$. This is also required for covariance estimation. The second part of Assumption (ref) (uniform consistency of the covariance estimator), together with Lemma (ref), is needed for the estimators of the coefficient matrices $B, C$ and $D$ to be consistent.

The third and the fourth part of Assumption (ref) are related to the estimation of $\sigma \sigma^\top$. The third part ensures that we can estimate the quadratic variation of the random process $X$ with the quadratic variation of $r$. Note that the third part implies that $\mathbb{E}[X_{1/N,i}^2]$ decays, for all $i=1,\ldots, d$, to zero slower than $N^{-2}$, meaning that none of the components $X_i$ are smooth (in the mean-square sense). This implies that the drift term in (ref) is asymptotically negligible in the computation of the quadratic variation. Finally, the fourth part ensures that the quadratic variation of $X$ is asymptotically deterministic and behaves as its expectation.

We now present our main results. Theorem (ref) is related to consistency of the estimators and Theorem (ref) is related to limiting distributions.

Examples and discussion

Suppose first that $X$ is a Gaussian process. Then $U$ is Gaussian as well. Consider first the case $\Vert\gamma(t)\Vert \sim t^{2H-2}$ for $0<H<1$ as $t\to \infty$. Then Assumption (ref) parts (i)-(ii) are valid, see Marko1. For the rate functions $l_1(T)$ and $l_2(T)$ in Theorem (ref), note that

and hence the rate function $l_1(T)$ for $\hat b_{U,T}$ is

Similarly,

and thus the rate function $l_2(T)$ for $\hat\gamma_{U,T}(t)$ is

Therefore, $l_1^2 > l_2$ for $0<H\leq\frac{3}{4}$ and $l_1^2 = l_2$ for $\frac{3}{4}<H<1$. The case $l_1^2(T)<l_2(T)$ is not possible. If the covariance function decays faster than $t^{2H-2}$, e.g., exponentially fast, we obtain the standard rate $l_1(T) = l_2(T) = \sqrt{T}$ for both $\hat{b}_{U,T}$ and $\hat\gamma_{U,T}$.

For Assumption (ref) parts (iii)-(iv), note that typically (iii) is automatically satisfied whenever $X$ is not smooth. In particular, (iii) is automatically satisfied whenever the components of $X$ are merely H\"older continuous for some $H\in(0,1)$, see azmoodeh2014necessary, and then $\Vert V(t)\Vert \sim t^{2H}$ for some $H\in(0,1)$. For (iv), note that, in the Gaussian case, the matrix $f(j,k,N)$ is a diagonal matrix with diagonal entries of the form $2\eta^2_i+\sum_{v=1,v\neq i}^d \eta_i\eta_v$, where $\eta_v = \mathrm{Cov}(\Delta X^{(v)}_{t_k},\Delta X^{(v)}_{t_j})$. If $X$ is the standard $d$-dimensional fractional Brownian motion (fBm) with Hurst index $H\in(0,1)$, this behaves as $N^{-4H}|j-k|^{4H-4}$. This yields (iv) directly. More generally, in the Gaussian case, the behavior of $\mathrm{Cov}(\Delta X^{(v)}_{t_k},\Delta X^{(v)}_{t_j})$ is closely related to the decay of $\gamma(t)$ determined by the behavior of function $V(t)$ for $t\to \infty$ while $V(1/N)$ is related to the function $V(t)$ for $t\to 0$.

In the non-Gaussian case, Assumption (ref) (i) and (ii) are also typically valid provided that the memory vanishes as the lags grow without a limit. One can again compute the rates $l_1(T)$ and $l_2(T)$ explicitly from $\left\Vert\mathrm{Cov}(\hat b_{U, T})\right\Vert$ and $ \left\Vert\mathrm{Cov}(\operatorname{vec}(\hat \gamma_{U,T}(t))) \right\Vert$. However, unlike in the Gaussian case, also $l_1^2(T)<l_2(T)$ is possible, see lietzen2021modeling for examples. For Assumption (ref) parts (iii)-(iv), part (iii) is again related to the smoothness of $X$. For example, under hypercontractivity for the moments, one has $\Vert V(t)\Vert \sim t^{2H}$, exactly as in the Gaussian case, whenever the components are H\"older continuous, see nummi2024necessary. Finally, for part (iv), as long as the components of $X$ are independent, the matrix $f(j,k,N)$ is diagonal with $i$:th entry given by $$ \mathrm{Cov}\left((\Delta X^{(v)}_{t_k})^2,(\Delta X^{(v)}_{t_j})^2\right) + \sum_{v=1,v\neq i}^d \eta_i\eta_v. $$ Here the asymptotic behavior of $\mathrm{Cov}\left((\Delta X^{(v)}_{t_k})^2,(\Delta X^{(v)}_{t_j})^2\right)$ is linked to $ \left\Vert\mathrm{Cov}(\operatorname{vec}(\hat \gamma_{U,T}(t))) \right\Vert$ while the asymptotic behavior of the terms $\eta_v$ is linked to the asymptotic $\left\Vert\mathrm{Cov}(\hat b_{U, T})\right\Vert$, and hence any of the terms can be the dominating one.

As a concrete non-Gaussian example, we mention the $d$-dimensional Hermite process, a non-Gaussian generalization of the fBm. In this case, the behavior of all of the terms is similar to the case of the fBm, and hence all the assumptions are valid. Finally, we mention suitable Levy processes (such as Poisson processes or $\alpha$-stable processes) as an example of non-continuous drivers $X$. In this case, due to the independence of the increments, the covariances decay rapidly, and hence the rates are the standard $\sqrt{T}$. In the presence of jumps, one does not have smooth paths. This provides (iii). Finally, the independence of the increments gives (iv) under suitable assumptions on the Levy measure. The broad range of processes that fulfill the required assumptions allows constructing, e.g., interest rate models that incorporate linear combinations of these processes that enables simultaneous modeling of jumps and various memory effects.

Simulated examples

In order to illustrate finite sample behavior and the applicability of the provided modeling approach, we present two different simulated examples. We provide one example where $\Theta$ is a diagonal matrix and one with a non-diagonal $\Theta$. In both examples, two-dimensional interest rates with fBm as the random process $X$ were simulated using the Euler-Maryama method. The model parameters and the Hurst indices $H$ for the simulations are given in Table (ref). The processes were simulated with 1000 different realizations of fBm, resulting in 1000 sample paths for each set of parameters in both examples. The length of the time series in every sample was 4000, consisting of 10000 observations.

After simulating, the model parameters $b, \Theta$ and $\sigma$ were estimated from each sample. The upper bound used in the integrals for estimating the matrices $B, C$ and $D$ was 5. Figures (ref) and (ref) display, in the case of diagonal $\Theta$, the componentwise differences between the estimates $\hat \Theta$ and the true value $\Theta$ as histograms for $H=0.35$ and $H=0.8$, respectively. Based on the histograms, it seems that the parameter $\Theta$ can be estimated rather well from the simulated data and, as expected, the differences of each component seem normally distributed. However, as $H$ increases, the estimation error increases as the estimates for $\sigma$ become more inaccurate, see Figure (ref). The histograms in Figure (ref) also hint to that direction, as the centers of the distributions corresponding to the diagonal elements of $\Theta$ are slightly shifted from zero. This phenomenon is expected as well, since when $H$ increases, one has longer memory inducing slower rate of convergence.

In the case of non-diagonal $\Theta$ and Brownian motion as the noise, Figure (ref) displays the componentwise values of the differences $\hat \Theta-\Theta$ and Figure (ref) displays the Frobenius norms of the differences. The histograms in Figure (ref) seem normal. In this case, due to the independence of the increments, one has fast convergence towards normal distribution with mean zero. It follows that then the Frobenius norm is essentially a weighted sum of $\chi^2$-distributed random variables. This is visible in Figure (ref).

Other parameter combinations yield similar results. Histograms under all our parameter combinations are provided in an electronic supplementary material of the article, see repo. Moreover, the repository provides a Python implementation for the proposed approach and the simulation results as numerical tables.

Real data

To further elaborate the practical usability of our proposed approach, we used it to model bivariate real data. The data consist of 1-month Euribor interest rate euribor and the Federal Funds Effective Rate (DFF) dff with 6765 daily observations from the time period 4.1.1999 - 6.5.2025.

Based on visualizing the data, the paths of the interest rates are somewhat smooth, indicating that their Hurst indices should be larger than $\frac12$ (see Figure (ref)). This guided us to use $V(t) \sim t^{2H} I$ for $H\in\left(\frac12,1\right)$. We tried different values within the $\left(\frac12,1\right)$ range, and they all produced similar estimates. Therefore, estimation of the Hurst index does not seem to play a large role here. In the sequel, we report the results for the case $H=0.7$.

The estimated model parameters were

where the first dimension corresponds to the Euribor rate and the second dimension corresponds to the DFF rate. Figure (ref) displays the paths of both interest rates and their respective estimated mean levels $\hat b$.

The components of $\hat \sigma$ can be interpreted as the estimated volatility of the components of $r$. The value corresponding to DFF is a bit larger than the value corresponding to Euribor. This can be seen when looking at the first half of the time series in Figure (ref): the magnitude of fluctuations of DFF is slightly larger than that of the Euribor rate.

The diagonal components of $\hat \Theta$ can be interpreted as the estimated mean reversion rates of the components of $r$. A larger value pushes the interest rate more strongly towards its mean, whereas a smaller value allows the rate to go and stay further away from the mean. Therefore, a larger value corresponds to a smaller variance of the interest rate.

The off-diagonal component of $\hat \Theta$ describes the interaction between the components of $r$. A negative value indicates a positive correlation, meaning that the rates are often on the same side of their means. Each rate tries to push the other one to the same direction as it itself is: down if the rate itself is below its mean and up if the rate is above its mean.

In the model with Euribor and DFF, the off-diagonal value is negative, and the interest rates seem visibly positively correlated as they are often moving in the same direction and on the same side of their means, see Figure (ref). Between 2016 and 2020, DFF increases while Euribor stays low. However, even though DFF is increasing during this period, it still stays below the mean level. Therefore, it is pushing Euribor downwards compensating the mean reversion effect.

A positive off-diagonal value would indicate a negative correlation with similar but opposite interactions between interest rates. A value close to zero compared to the diagonal values would indicate no correlation between interest rates.

The model assumes stationary increments of the random process $X$. The estimated increments of $X$ for the Euribor-DFF data are visualized in Figure (ref). Based on the figure, the stationarity assumption is not violated too badly: both increment series seem fat tailed, but they do not show trends or clear heteroskedasticity.

Finally, to assess the model performance, we provide one-step predictions of the interest rates for the last fifth of the time series (1353 last observations). The predictions are visualized in Figure (ref). The predictions seem to be very good.

Proofs

All the technical proofs are presented in this section. The section is divided into three subsections. Section (ref) contains the statements and proofs of lemmas that are used in proving the main results. Section (ref) provides the proof of Theorem (ref) and Section (ref) contains the proof of Theorem (ref).

Throughout this section, we use the $O_P$ notation for stochastic boundedness: a random object $Z_n$ satisfies $Z_n= O_P(a_n)$ as $n \to \infty$, if for any $\epsilon >0$ there exist finite $M > 0$ and $N>0$ such that

for all $n>N$.

Auxiliary lemmas

The first two lemmas in this section, Lemma (ref) and (ref), are fundamental results used throughout the section. Lemmas (ref), (ref) and (ref) are used in the proof of Theorem (ref). The remaining two lemmas, Lemma (ref) and (ref), are used in the proof of Theorem (ref).

Proof of Theorem (ref)

Proof of Theorem (ref)