Dynamically Consistent Analysis of Realized Covariations in Term Structure Models

Dynamically Consistent Analysis of Realized Covariations in Term Structure Models

Introduction

We present a nonparametric method to measure covariations in a general arbitrage-free term structure setting in the spirit of bjork1999. A motivation is to determine the number of statistically relevant random drivers needed to describe bond market dynamics under minimal assumptions and in a dynamically consistent manner. That is, by working coherently in the abstract setting of bjork1999, we circumenvent the well-known consistency problems of finding arbitrage-free finite-dimensional dynamics that reflect the empirical observations.

In the bond market the notion of a zero coupon bond is fundamental. A zero coupon bond guarantees its holder at time $t$ a fixed amount of money at some time $t+x$ in the future. The cost $P_t(x)$ of entering this contract at time $t$ is the price of this bond, which depends on the time to maturity $x$, implying a price curve $x\mapsto P_t(x)$ at each time $t\geq 0$ called the discount curve. We assume to observe bond price or yield curve data, potentially derived by smoothing as in FPY2022 or liu2021, such that we can recover log bond prices $$P_{i,j}^n:=\log P_{i\Delta_n}(j\Delta_n)\quad j=0,1,...,\lfloor M/\Delta_n\rfloor, \,\, i=0,1,...,\lfloor T\Delta_n\rfloor$$ for a resolution $\Delta_n=1/n$ and with different maturities $j\Delta_n $ and on different time points $i\Delta_n $. Here $M>0$ is some maximal time to maturity (e.g. $M=10$ or $M=30$ when time is measured in years) and $T$ is the time until which the data are observed or considered.

Often, risk factor analyses are conducted on the basis of transformations of the discount curve, such as yield differences or excess returns, which typically suggests that three factors explain a large amount of variation in bond market dynamics, c.f. Litterman1991. Recently Crump2022, raised the concern that these low-dimensional factor structures are obtained irrespectively of the data generating process due to the high correlation of bond prices with close maturities. To remedy this effect, dimension reduction could be based on difference returns, which are the returns of the trading strategy of buying an $x+\Delta$-maturity bond and shorting an $x$-maturity bond. Precisely, difference returns are defined for $i=1,...,\ensuremath{\lfloor T/\Delta_n\rfloor}-1$ and $j=1,...,\lfloor M/\Delta_n\rfloor-1$ by

In this article we develop an asymptotic econometric theory for the realized covariations of these difference returns, that is, for $j_1,j_2=1,...,\lfloor M/\Delta_n\rfloor$ we analyse the covariations $$\hat q_T^n(j_1,j_2):= \sum_{i=1}^{\ensuremath{\lfloor T/\Delta_n\rfloor}} d_i(j_1)d_i(j_2).$$

Importantly, while $\ensuremath{\lfloor T/\Delta_n\rfloor}^{-1} q_T^n$ is the empirical covariance of the data $d_1^n,...d_{\ensuremath{\lfloor T/\Delta_n\rfloor}-1}^n$, assuming w.l.o.g. $\mathbb E[d_i^n]=0$, we do not consider it as an estimator of the population covariance of difference returns. Such an interpretation is not invariant with respect to the resolution $\Delta_n$ and requires the restrictive assumption that difference returns are i.i.d. or at least covariance stationary and ergodic. More importantly, it is not clear how dimension reduction can be conducted without entailing arbitrage opportunities. General arbitrage-free term structure models in the sense of bjork1997 and FilipovicTappeTeichmann2010b require that forward rates $ f_t(x)=-\partial_x\log(P_t(x))$ for $ x,t\geq 0,$ satisfy dynamics of the form

where the equation holds in an appropriate function space. The latent process $X$ is a possibly infinite-dimensional It{\^o} semimartingale $$X_t:= \int_0^t \alpha_s ds+\int_0^t\sigma_s dW_s+J_t. $$ where $\alpha$ is a curve-valued drift process, $\sigma$ is the (in general operator-valued) volatility, $W$ is an (in principle infinite-dimensional) Wiener process and $J$ is a jump process that we assume to model rare extreme events (c.f. Section (ref) for the details). Under a risk-neutral measure, the drift $\alpha$ can in addition be characterized as a deterministic function of the volatility $\sigma$ and characteristics of the jump process $J$ (c.f. bjork1997, FilipovicTappeTeichmann2010b). Parametrizations of forward cuves are then required to be viable in the dynamic setting (ref) to avoid the introduction of arbitrage opportunities to the model. This is known to be an intricate problem in term structure modelling (see e.g. bjork1999, bjork2001, filipovic2003, Filipovic2000, filipovic2000b ) and some frequently employed parametrizations of forward curves are incompatible with arbitrage-free dynamics or induce restrictive additional conditions (see e.g. filipovic1999).

The concerns on realized covariations of difference returns can be resolved when we consider infill asymptotics ($n\to \infty$). Precisely, we show without imposing further assumptions and interpreting $\hat q_T^n$ as a piecewise constant kernel that

where the limits hold in $L^2([0,M]^2)$. The right hand side describes the quadratic covariation of the latent driver $X$, which always exists (see e.g. Schroers2024) and can equivalently be described as the limit of covariance operators in $L^2(\mathbb R_+)$ by

where $\langle \cdot, \cdot\rangle$ is the inner product of $L^2(\mathbb R_+)$ and $\Delta_i^n X:= X_{i\Delta_n}-X_{(i-1)\Delta_n}$ denotes the $i$'th increment of $X$.

The limit (ref) implies that it is possible to infer on the number of random drivers in the bond market on the basis of discrete bond price data without further assumptions on moments, stationarity and ergodicity. In fact, if the quadratic variation of $X$ is $d$-dimensional for $d\in \mathbb N$, then $X$ is $d$-dimensional and its state space (until time $T$) is spanned by the eigenvectors corresponding to the nonzero eigenvalues of $[X,X]_T$. The eigenvalues of $[X,X]_T$ indicate the amount of variation that the corresponding random factor explained of $X$ up to time $T$. This is in contrast to the explained variation of factors derived from a covariance of $f$ or its increments, which is a priori not informative on the number of random drivers. In fact, it is possible that $f$ is infinite-dimensional, while $X$ is a one-dimensional process (c.f. Example 4.1 in Schroers2024). The most striking advantage of the interpretation of the realized covariation of difference returns in (ref) is, however, that exchanging $X$ by an arbitrary finite-dimensional semimartingale (with the correct form of the drift) in the formulation of the dynamics in (ref) does not affect the capability of the model to be free of arbitrage, such that an investigation of the number of relevant factors can be conducted independently of further consistency conditions. This underlines that our abstract infinite-dimensional setting relaxes the analysis of the term structure when compared to models in which state spaces are assumed to be finite-dimensional a priori. An additional advantage of quadratic variations is that they are naturally interpreted as time-varying objects enabling their temporal analysis.

In practice, the distortion of the measurements due to ouliers can bias the analysis of the relevant factors. For instance, in the context of sudden interest rate movements during an economic crisis it is possible that a single outlying difference return impacts the measurement of covariations and, thus, the measured dimensionality of the driver $X$ substantially. For this reason, besides $[X,X]$, the continuous part $[X^C,X^C]$ of the quadratic variation where $X^C_t= \int_0^t \alpha_s ds+\int_0^t \sigma_s dW_s$ is central for the task of identifying the statistically relevant number of random drivers, considering the jump part to model outlying events. We describe how estimation of the continuous quadratic covariation is possible by a truncation technique, which sets outlying difference returns in the realized covariation on the left of (ref) to $0$. We derive rates of convergence and a central limit theorem for these estimators and also show how the long-time limit of $[X^C,X^C]_T/T$ as $T\to \infty$ can be estimated, if it exists. Our limit theory holds under weak assumptions, which mainly reflect those for finite-dimensional semimartingales, although, $f$ does not need to be a semimartingale.

We conduct an emprical study on the relevant drivers in the market via this covariation estimates based on real bond market data. The procedure is numerically equivalent to a principal component analysis based on the (truncated) empirical covariance of difference returns with a daily resolution and mean zero. However, the classification of jumps takes into account the abstract setting (ref). By investigation of a truncated version of the realized covariations $q_j^*(x,y):=\Delta_n^{-2} \sum_{i=\lfloor ( j-1)/\Delta_n \rfloor}^{\lfloor j/\Delta_n \rfloor} d_i(\lfloor x/\Delta_n\rfloor)d_i(\lfloor y/\Delta_n\rfloor)$ for any year $j$ from 1990 to 2022, we find evidence for the dimension of the driver to vary from year to year but also to be consistently high (in each year more than $8$ drivers are needed to explain at least $99\% $ of the variation). We further observe that quadratic variations vary in shape over time and not just their level. We provide Monte-Carlo evidence for the validity of the limit theory in the context of sparse and noisy data.

Formal validity of our method is guaranteed by relating difference returns via cross-sectional and temproal discretization to the abstract setting (ref) and then apply the results from the article Schroers2024. The truncation procedure is inspired from the truncated realized variation estimators of Mancini2001, Mancini2004, Mancini2009 and Jacod2008 for finite-dimensional semimartingales. Here, we also provide a data-driven variant of the truncation rule, to account for functional outliers in a similar way as the trimmed least squares method in ren2017. Naturally, our asymptotic theory also allows for nonparametric estimation of characteristics of infinite-dimensional volatility models in continuous time employed for term structure modeling (c.f. BenthRudigerSuss2018, BenthSimonsen2018, BenthSgarra2021, Cox2020, Cox2021, Petersson2022 and cox2023).

The article is structured as follows: Section (ref) describes the general bond market setting that we consider for this article. Section (ref) presents the estimation theory for the central application of term structure models. Identification for the quadratic variations of $X$, $X^C$ and $J$ on the basis of difference return variations can be found in Section (ref), while rates of convergence for estimating $[X^C,X^C]$ and a central limit theorem can be found in Section (ref). Section (ref) discusses long-time asynmptotics for estimation of a stationary mean of $[X^C,X^C]_T/T$. Section (ref) contains practical considerations on smoothing of discrete bond price data and presents a data-driven truncation rule for robust estimation. Section (ref) provides a simulation study. Finally, we apply our theory to bond market data in Section (ref). Technical proofs of our results along with further remarks on the simulation scheme and additional empirical results can be found in the appendix.

Technical preliminaries and notation

Let $I$ be an interval in $\mathbb R$. We write $\langle h,g\rangle= \int_I h(x)g(x) dx$ for the $L^2$-scalar product of two elements $h,g\in L^2(I)$ as well as $\|h\|=\sqrt{\langle h,h\rangle}$ for the norm. We write $L_{\text{HS}}(L^2(I))$ for the Hilbert space of Hilbert-Schmidt operators from $L^2(I)$ into itself and $\|\mathcal T_k\|_{\text{HS}}$ for the Hilbert-Schmidt norm of an $\mathcal T_k\in L_{\text{HS}}(L^2(I))$. Recall, that a Hilbert-Schmidt operator $\mathcal T_k:L^2(I)\to L^2(I)$ can be uniquely associated to a kernel $k\in L^2(I^2)$ such that $\|\mathcal T_k\|_{\text{HS}}=\|k\|_{L^2(I^2)}$ and

Importantly, the operator $h\otimes g:= \langle h,\cdot\rangle g$ is Hilbert-Schmidt for two elements $h,g \in L^2(I)$. We shortly write $h^{\otimes 2}=h\otimes h$. Finally, for $L^2(I)$-valued processes $X^n, n\in \mathbb N$, $X$, we write $X^n\stackrel{u.c.p.}{\longrightarrow}{X}\quad \text{ as }n\to\infty$ for the convergence uniformly on compacts in probability, i.e. it is $\mathbb P[\sup_{t\in [0,T]} \|X^n(t)-X(t)\|>\epsilon]\to 0$ for all $\epsilon, T>0$.

General Arbitrgae-Free Bond Market-Dynamics

Let $(\Omega,\mathcal F,(\mathcal F_t)_{t\geq 0},\mathbb P)$ be a filtered probability space with right-continuous filtration. From here on, we assume that the forward rate process $(f_t)_{t\geq 0}$ is an $L^2(\mathbb R_+)$-valued stochastic process that is the mild solution to the stochastic partial differential equation (ref), defined on $(\Omega,\mathcal F,(\mathcal F_t)_{t\geq 0},\mathbb P)$. That is,

where $ \mathcal S(t)f(x)=f(x+t)$ for $t\geq 0 $ and $f\in L^2(\mathbb R_+)$ defines the left-shift operator semigroup. We relegate all further technical discussions on $X$, $\alpha$, $\sigma$, $W$ and $J$ and various related technical assumptions that we need in for the validity of our limit theory to Sections (ref) in the appendix. We remark, that under arbitrage-free dynamics, that is, under an equivalent local martingale measure the drift is necessarily a deterministic function of $\sigma$ and $\gamma$ (c.f. bjork1997), which was in the continuous case the original inside leading to the popular Heath-Jarrow-Morton framework of HJMoriginal for pricing bonds and interest rate sensitive contingent claims. However, this will not be of particular importance for our purposes, as the drift later vanishes asymptotically in our limit theory. We will discuss however some important examples subsequently. Before that, we make a remark on the choice $L^2(\mathbb R_+)$ as the state space of forward rate curves.

Let us now discuss some important simple examples.

The term structure setting (ref) contains the vast majority of existing arbitrage-free term structure models considered in the literature. Among them is the class of affine term structure models, which are widely appreciated for their parsimony.

Many simple ways to model term structures lead to nonaffine dynamics, such as

Example (ref) (and also Example 3.16 in BSV2022) shows that the process $(f_t)_{t\geq 0}$ is in general not an $L^2(\mathbb R_+)$-valued semimartingale. However, all implied bond prices $(P_t(T-t))_{0\leq t\leq T}$ are semimartingales for all $T>0$, which is necessary to guarantee the absence of arbitrage in the bond market. Moreover, while the quadratic covariations of $f$ must not be convergent, we show in the next section, in which we present our main results, that the realized covariation of difference returns measures quadratic covariations of the latent driver asymptotically and without further conditions.

Estimation of Quadratic Covariations

In this section we present our asymptotic theory for estimation of quadratic variations. We start with the identifiability of $[X,X]$, $[X^C,X^C]$ and $[J,J]$.

Identification of the quadratic covariation of the latent semimartingale

We rely on infill asymptotics $\Delta_n\to 0$ as $n\to \infty$ and recall the definition of the realized covariation $(\hat q_t^n)_{t\geq 0}$ as a piecewise constant kernel. That is, for $x\in [(j_1-1)\Delta_n,j_1\Delta_n] $ and $y\in [(j_2-1)\Delta_n,j_2\Delta_n]$, $j_1,j_2\in\mathbb N$ and $t\geq 0$ we define

For each $n\in \mathbb N$ and $t\geq 0$, we have that $\hat q_t^n\in L^2(\mathbb R_+^2)$, which follows from the Assumption that forward curves are elements in $L^2(\mathbb R_+)$ (c.f. Remark (ref) below). We now state the general identifiability result for the quadratic covariation of $X$.

As in the case of finite-dimensional semimartingales, we do not have to impose any further conditions to identify the quadratic covariations of the driving semimartingale $X$ although the observable process $(f_t)_{t\geq 0}$ is not necessarily an $L^2(\mathbb R_+)$-valued semimartingale. This is due to the relation of difference returns to semigroup-adjusted forward rate returns, which where shown in Schroers2024, Benth2022 and BSV2022 to be well-suited for volatility estimation for processes of the form (ref). This relationship is made clear in the next remark.

Next, we examine how to identify the continuous part of the quadratic covariation.

Identification of $[X^C,X^C]$ and $[J,J]$ via truncated covariation estimators

We now turn to the estimation of the continuous part of the quadratic covariation. We will derive jump robust estimators by a truncated form of $\hat q_t^{n}$ defined by the piecewise constant kernel $\hat q^{n,-}_t$ given for $x\in [(j_1-1)\Delta_n,j_1\Delta_n]$, $y\in [(j_2-1)\Delta_n, j_2\Delta_n]$, $j_1,j_2\in \mathbb N$ and $t\geq 0$ by

for $u_n=\alpha \Delta_n^{w}$, with $w\in (0,1/2)$ and $\alpha>0$ and a particular sequence of truncation functions $g_{n}$ that takes into account only the discrete data $d^i_n(j)$ for $j\in \mathbb N$. Precisely, the corresponding sequence of truncation functions $g_n: l^2\to \mathbb R_+$ must satisfy for constants $c,C>0$ and for all $f,h\in l^2$ and all $n\in \mathbb N$

While the particular choice of the functions $g_n$ will not play a role for the asymptotic behavior of $\hat q^{n,-}_t$, it is important to modify it in practice. For the moment, one can take in mind the legitimate choice $g_n=\|\cdot\|_{l^2}$ for all $n\in \mathbb N$ for which we have with $\Pi_{n,\infty}$ defined as in (ref) for $M=\infty$ that $\|d_i^n/\Delta_n\|_{l^2}= \| \Pi_{n,\infty} \tilde \Delta_i^n f\|_{L^2(\mathbb R_+)}.$ We will discuss a data-driven specification of $g_n$ and the truncation level in Section (ref).

The next result states that $\hat q^{n,-}$ consistently estimates the quadratic covariation of $X^C$.

Let us make a remark on the feasibility of the estimator.

Observe that we can also define an upward truncated estimator $ \Delta_n^{-2}\hat q^{n,+}$ given for $t\geq 0$, $x\in [(j_1-1)\Delta_n,j_1\Delta_n]$, $y\in [(j_2-1)\Delta_n, j_2\Delta_n]$ and $j_1,j_2\in \mathbb N$ by

Obviously, $\hat q^n=\hat q^{n,-}+\hat q^{n,+} $ and $\mathcal T_{\hat q^n}=\mathcal T_{\hat q^{n,-}}+\mathcal T_{\hat q^{n,+} }$. Then, combining Theorem (ref) and Theorem (ref), we also obtain

This result shows that the quadratic covariation corresponding to the jump part is identifiable in the context of general bond market models. However, the finer analysis of jumps is not part of this paper and relegated to future work. Instead, we derive convergence rates for the estimation of the continuous part of the quadratic covariation in the next section.

Convergence rates and central limit theorem for estimation for $\hat q^{n,-}$

In order to derive rates of convergence and a central limit theorem for estimating the continuous part of the quadratic variation, we need to impose further regularity Assumptions, which depend on the smoothness of the kernel corresponding to the operators $[X^C,X^C]_t$. For the error bounds, this is Assumption (ref), which is discussed in detail in Section (ref). We discuss these Assumptions in the context of Example (ref) right below the subsequent abstract result.

The rate implied by Theorem (ref) is at most $\mathcal O_p(\Delta_n^{ 1/2})$, which is achieved if Assumptions (ref)(r) and Assumption (ref)($\gamma$) hold for $\gamma = 1/2$ and $r<1$. Let us now discuss Theorem (ref) and Assumptions (ref)($\gamma$) and (ref)($r$) in the context of Example (ref).

If we assume a slightly stronger Assumption than (ref)($1/2$), which can also be found in Section (ref), we can even obtain a stable CLT in the next result, where stable convergence in law is denoted by $\overset{st.}{\longrightarrow}$ \footnote{ Recall that a sequence of random variables $(X_n)_{n\in\mathbb N}$ defined on a probability space $(\Omega, \mathcal F,\mathbb P)$ and with values in a Hilbert space $H$ converges stably in law to a random variable $X$ defined on an extension $(\tilde{\Omega}, \tilde{\mathcal F},\tilde{\mathbb P})$ of $(\Omega, \mathcal F,\mathbb P)$ with values in $H$, if for all bounded continuous $f:H\to \mathbb R$ and all bounded random variables $Y$ on $(\Omega,\mathcal F)$ we have $\mathbb E[Y f(X_n)]\to\tilde {\mathbb E}[Y f(X)]$ as $n\to\infty$, where $\tilde {\mathbb E}$ denotes the expectation w.r.t. $\tilde{\mathbb P}$.}.

The partial derivative $\Sigma_t:=\partial_t [X,X]_t$ has to be interpreted as a Frechet-derivative and does always exists, due to the Assumption on $X$ being an It{\^o} semimartingale (c.f. Section (ref)). Let us derive the form of $\mathfrak Q_t$ in the context of Example (ref):

So far we have discussed limit theorems for infill asymptotics leaving $T$ fixed. In the next section, we outline how it is possible under additional assumptions and letting $T\to\infty$ to make use of all available data to estimate the stationary instantaneous covariance for difference returns.

Long-time volatility estimation

The truncated estimation procedure described previously enables estimations of a time series of the integrated volatilities $\int_{i}^{i+1}\Sigma_s ds=[X^C,X^C]_i-[X^C,X^C]_{i-1} $ for $i\in \mathbb N$. If the aim is to derive a time-invariant mean for the volatility, we have to impose further conditions, which are described in detail in Section (ref). These Assumptions are much stricter than the ones we considered in the previous section for the infill asymptotics on finite intervals and in particular imply that the mean $$\mathcal C:=\frac 1 T\mathbb E[[X,X]_T]$$ is independent of $T$. However, they allow us to derive a stationary mean of $\Sigma_t$ via large $T$ asymptotics.

Assumption (ref) does not impose very strong Assumptions on the dynamics of the volatility and is satisfied by most stochastic volatility models. To verify this condition for particular models for the infinite-dimensional volatility process one might investigate the vast literature for ergodic properties of Hilbert space-valued processes and, in particular, SPDEs (c.f. DPZ2014 or PZ2007). For the existence of invariant measures for term structure models, we further mention vargiolu1999, Tehranchi2005, marinelli2010, rusinek2010, and FFRS2020. Recently, friesen2022 examined the long-time behavior of infinite-dimensional affine volatility processes. Here, we only review the validity of the Assumptions employed in Theorem (ref) in the context of our running Example (ref):

Practical considerations

In this section, we discuss some practical complications on the implementations of the estimator. Namely, we present a data-driven truncation rule and comment on the use of nonparametrically smoothed yield or bond price curve data.

We start with a data driven choice of the truncation function $g_n$ and the tuning parameters $\alpha$ and $w$.

Truncation in practice

While the asymptotic theory of Section (ref) justifies the use of truncated estimators, the choice of the truncation level and functions remains a practical issue. Even in finite dimensions, this can be challenging and we refrain from finding optimal choices. However, we outline how the truncation rule can be reasonably implemented.

Truncation rules in finite dimensions often necessitate preliminary estimators for the average realized variance in the corresponding interval of interest (c.f. Mancini2012, page 418, for an overview of some truncation rules). One sorts out a large amount of data first, to obtain a preliminary estimator of $[X^C,X^C]_T$. As this can be interpreted as the average covariation of the increments in the interval $[0,T]$, one then chooses truncation levels in terms of multiples of standard deviations as measured by the preliminary estimate. In our infinite-dimensional setting, we mimic this procedure, but it is harder to distinguish typical increments and outliers as we cannot argue componentwise. While the choice $g_n\equiv \|\cdot\|$ leads to consistent estimators in terms of the limit theory developed in Section (ref), it is not necessarily a good choice in the context of finite data since the continuous martingale might vary considerably more in one direction than another .

We Therefore present a method that is based on a measure of functional outlyingness in the spirit of ren2017: Assuming that $\Sigma$ is independent of the driving Wiener process and does not vary too wildly on the interval $[0,T]$, and that no jumps exist, we have approximately that $\Sigma_t \approx \frac 1T\int_0^T \Sigma_s ds$ for $t\in [0,T]$ and $\tilde \Delta_i^n f/\sqrt{\Delta_n}|\Sigma\sim N(0, \frac 1T\int_0^T \Sigma_s ds)$. If the largest $d$ eigenvalues $e_1,...,e_d$ of $\int_0^T \Sigma_s ds/T$ account for a large amount of the variation as measured by the summed eigenvalues of $\int_0^T \Sigma_s ds/T$ (e.g. 90 percent), we know that $P_d \tilde \Delta_t^n f$ with $P_d=\sum_{i=1}^d e_i^{\otimes 2}$ is a linearly optimal approximation of $\tilde \Delta_t^n f$ in the $L^2(\mathbb R_+)$-norm. We can also define $P_d(\int_0^T \Sigma_s ds/T)^{-1}P_d =\sum_{i=1}^d(1/\lambda_i)e_i^{\otimes 2}$ and define $g^d(h):=\| P_d(\int_0^T \Sigma_s ds/T)^{-1}P_d h\|^2$. This distance resembles the measure proposed in ren2017, however, it is not a valid truncation function, since (ref) cannot hold. Further, if a truncation at level $d$ is made, outliers impacting the higher eigenfactors might be overlooked. We Therefore propose an adjusted method defining $g(x)=g^d(x)+ \frac{\|(I-P_d)x\|^2}{\sum_{i=d+1}^{\infty} \lambda_i} $. Then, for $p_n:l^2\to L^2(\mathbb R_+)$ given by $p_n x:= \sum_{j=1}^{\infty} x_j \ensuremath{\mathbb{I}}_{[(j-1)\Delta_n,j\Delta_n]}$ we define the sequence of truncation functions $(g_n)_{n\in \mathbb N}$ via

where the index $d$ can be chosen freely as long as $\lambda_{d+1}>0$. E.g. we can choose $d$ such that the first $d$ eigenfactors for $\mathcal C$ explain $90$ % of the variation. It is then with $\lambda_i^d:=\lambda_i$ for $i\leq d$ and $\lambda_i^d:=\sum_{j=d+1}^{\infty} \lambda_j$ for $i\geq d$ and with $\Delta_n$ small

In practice, we do not know the eigenvalues and eigenfunctions of $\int_0^T \Sigma_s ds/T$ and derive them from a preliminary estimate. We suggest a simple truncation procedure in two steps:

Arguably, there can be many other methods for deriving truncation rules, which however have to take into account the infinite dimensionality of the data and deal with the subtlety of functional outliers. The simulation study in Section (ref) shows the good performance of our method.

Presmoothing bond market data

It is rarely the case that term structure data are observed in the same resolution in time as in the maturity dimension. For bond market data, points on the discount curve $x\mapsto P_{i\Delta_n}(x)$ are observed irregularly with a lower resolution than daily along the maturity dimension. Additionally, information on the discount curve is sometimes latent as bonds are often coupon-bearing and assumed to be corrupted by market microstructure noise. To account for these difficulties and in accordance with the classical “smoothing first, then estimation" procedure for functional data analysis advocated in Ramsay2005 we pursue the simple yet effective approach of presmoothing the data. We derive smoothed yield or discount curves, as described, for instance, in FPY2022, liu2021 or Linton2001, which allows us to derive approximate zero coupon bond prices for any desired maturity and for which the impact of market microstructure noise is mitigated. While some theoretical guarantees in terms of asymptotic equivalence of discrete and noisy to perfect curve observation schemes could be derived for certain smoothing techniques and the task of estimating means and covariances of i.i.d. functional data (c.f. Zhang2007), in our case they would depend on the respective smoothing technique, the volatility and the semigroup as well as the magnitude of distorting market microstructure noise. A detailed theoretical analysis in that regard is beyond the scope of this article and instead, we showcase the robustness of our approach in the context of sparse, irregular and noisy bond price data within a simulation study in Section (ref).

Simulation study

In our simulation study we examine the performance of the truncated estimator $\mathcal T_{\hat q^{n,-}_{\cdot}}$ defined in (ref) as a measure of the continuous part $[X^C,X^C]$ of the quadratic covariation of the latent driver. As an important application of our theory is the identification of the number of statistically relevant drivers, we also examine how reliable the estimator can be used to determine the effective dimensions of $X^C$. In this context, we also want to assess the robustness of our estimator concerning three important aspects: First, we need to confirm the robustness of the truncated estimator to jumps. Second, we study the effect of the common practice of presmoothing sparse, noisy, and irregular bond price data on the estimator's performance. Moreover, we examine how the routine of projecting these data onto a small finite set of linear factors (c.f. for instance Litterman1991 or the survey piazzesi2010), influences conclusions on the quadratic covariation.

For that, we simulate log bond prices for some sampling size $m\leq 1000$, and $n=100$ time points, that is, $$ P_{i,l}:=\log P_{i\Delta_n}(j_{i,l}\Delta_n)+\epsilon_{i,l}, \qquad i=1,...,100,\quad l=1,...,m,$$ where $(\epsilon_{i,j})_{i,j=1,...,100}$ are i.i.d centered Gaussian errors with variance $\sigma_{\epsilon }>0$ and the $j_{i,r}$ are drawn randomly from $\{0,...,1000\}$ without replacement for $r=1,...,m$. We distinguish two cases: First, as a benchmark, we observe the data densely and without noise such that $\sigma_{\epsilon}=0$ and $m=1000$ and, second, we observe the data with noise $\sigma_{\epsilon}=0.01$ and sparsely with $m=100$. In this case, the prices for all maturities $j\Delta_n, j=1,...,1000$ are recovered by quintic spline smoothing. The roughness penalty for the smoothing splines is for each date chosen by a Bayesian information criterion and implemented via the ss-function from the npreg package in R. Using quintic splines and a Bayesian information criterion induces smooth implied forward curves.

To analyze the impact of the customary procedure of projecting the bond price data onto a low-dimensional linear subspace, we conduct our experiments in two scenarios. Scenario 1 in which we do not project the log bond prices and Scenario 2 in which we project the log bond prices onto the first three eigenvalues of their covariance $c_{logbond}\equiv\frac 1{100}\sum_{i=2}^{100} (P_{i,\cdot}-P_{i-1,\cdot})(P_{i,\cdot}-P_{i-1,\cdot})'$ before we calculate $\hat q^{n,-}_t$. Indeed, as usual for bond market data, the first three eigenvectors of $c_{logbond}$ explain over 99% the variation in the log bond prices.

The log bond prices are derived from simulated instantaneous forward rates $F_{i,j}:=\langle \ensuremath{\mathbb{I}}_{[(j-1)\Delta_n, j\Delta_n]}, f_{i\Delta_n}\rangle_{L^2(\mathbb R_+)}$ for $n=100$, $i=1,...,100$ and $j=1,...,1000$ from a forward rate process driven by a semimartingale $X$. Precisely, we define $X_t= \int_0^t \sqrt{\Sigma_s}dW_s+ J_t$ where

Here $x$ is a univariate mean-reverting square root process $$dx(t)=1.5\left(0.058- x(t)\right)dt+ 0.05 \sqrt{x(t)}d\beta(t), \quad t\geq 0,\,\,x(0)=0.058$$ and $Q_{a}$ is a covariance operator on $L^2(\mathbb R_+)$ such that the corresponding covariance kernel $q_a$ (s.t. $Q_a=\mathcal T_{q_a}$) restricted to $[0,M]$ is a Gaussian covariance kernel $q_{a}(x,y) \propto \exp(-a(x-y)^2)$ for some $a>0$ and $\|q_{a}\big|_{[0,M]^2}\|_{L^2([0,M]^2)}=1$. The jumps are specified by two Poisson processes $N^1,N^2$ with intensities $\lambda_1, \lambda_2>0$ and jump distributions $\chi_i^2 \sim N(0,\rho_2 Q_{0.01})$ and $\chi_i^1 \sim N(0,\rho_1K)$ for $\rho_1,\rho_2 \geq 0$ and where $K$ is another covariance operator with kernel $k(x,y)\propto e^{-(x+y)}$ and $\|k\big|_{[0,M]^2}\|=1$.

We specify the corresponding parameters of this infinite-dimensional model as follows: We choose $a=50$ reflecting a high dimensional setting since the decay rate of the eigenvalues of $Q_a$ is slow (10 eigencomponents are needed to explain 99% of the variation of $X^C$). The mean reversion level $0.058$ of the square root process corresponds to the Hilbert-Schmidt norm of the long-time estimator of volatility derived from bond market data as discussed in the next section. Jumps corresponding to the first component are considered large and rare outliers reflected by a high $\rho_1=0.0116$ and low $\lambda_1=1$. The second component describes outliers which are more frequent and smaller in norm reflected by a lower $\rho_2=0.0029$ and higher $\lambda_2=4$ but correspond to changes of the shape of the forward curves. Both jump processes are chosen such that their Hilbert-Schmidt norm accounts for approximately 10 % of the quadratic variation. We also consider cases in which no jumps are present (corresponding to the parameter choices $\lambda_1=\lambda_2=0$).

In each considered scenario, we compute the estimator $\Delta_n^{-2}\hat q_1^{n,10,-}$ (as defined in Remark (ref)) for $q_a\cdot \int_0^1 x(s)ds$. In the cases in which jumps are present, we consider the truncated estimator via the data-driven truncation rule of Section (ref) with different values of $l=3,4,5$ for the truncation level $u_n=l\sqrt{d+1}\Delta_n^{0.49}$. Here $d$ is chosen as the smallest value such that the first $d$ eigencomponents account for 90% of the variation as measured by the preliminary covariance estimator for which we truncate at the 0.75-quantile of the sequence of difference return curves as measured in their $l^2$ norm. Models M1 and M2 for which no jumps are present serve as benchmarks for the truncated estimators and no truncation is conducted ($l=\infty$).

We assess the performance of the respective estimator in the context of two criteria of which each reflects an important application of our estimator. First, we measure the relative approximation error $rE(\Delta_n^{-2} \hat q_1^{n,10,-}, IV)$ where for $x\in [(j_1-1)\Delta_n,j_1\Delta_n] $ and $y\in [(j_2-1)\Delta_n,j_2\Delta_n]$ the $\Delta_n$-resolution of the integrated volatility is $IV(x,y):=\int_0^1x(s)ds \cdot n^2\int_{(j_1-1)\Delta_n}^{j_1\Delta_n}\int_{(j_2-1)\Delta_n}^{j_2\Delta_n}q_{a}(z_1,z_2) dz_1,dz_2$ and $$rE(k_1,k_2) : \frac{\| k_1- k_2 \|_{L^2([0,10]^2)}}{\| k_2 \|_{L^2([0,10]^2)}}\qquad k_1,k_2\in L^2([0,M]^2).$$ Second, we will investigate how reliably the estimator can be used to determine the number of factors needed to explain certain amounts of variation of the latent driving semimartingale. For that, we define

for a symmetric positive nuclear operator $ C$ and an orthonormal basis $e=(e_i)_{i\geq 0}\subset L^2([0,M])$. Let $\hat e=(\hat e_i)_{i\in \mathbb N}$ denote the eigenfunctions of $\mathcal T_{\hat q_1^{n,10,-}}$ ordered by the magnitude of the respective eigenvalues. We report the numbers $D_C^e$ for $C=\mathcal T_{\hat q_1^{n,10,-}}$, $e=\hat e$ and $p=0.85$, $0.90$, $0.95$ and $0.99$, which are the numbers of factors needed to explain respectively $85\%$, $90\%$, $95\%$ and $99\%$ of the variation.

Table (ref) shows the results of the simulation study based on $500$ Monte-Carlo iterations. For Scenario S1, reflecting our proposed fully infinite-dimensional estimation procedure, the log bond prices are not projected onto a finite-dimensional subspace a priori. In this scenario, at least if jumps are truncated at a low level ($l=3$), the medians of relative errors are of a comparable magnitude when using either nonparametrically smoothed data (M2,M4) or perfect observations (M1,M3). While some jumps are overlooked by the truncation rule, the medians of the relative errors in the cases with jumps (M3,M4) just moderately increased compared to the respective cases in which no jumps appeared (M1,M2), at least for a low truncation level. The reported dimensions needed to explain the various levels of explained variation are estimated quite reliably (observe that the true thresholds are respectively $5,6,7$ and $10$). In the noisy and irregular settings (M2, M4) the estimators tend to add a dimension compared to the perfectly observed settings in the median but have low interquartile ranges, which contain the correct dimension. We conclude that the measurement of dimensions can be conducted accurately under realistic conditions.

Comparing scenarios S1 and S2 it becomes evident that the customary finite-dimensional projections of log bond prices (S2) affected the estimator's performance significantly. All of the medians of relative errors are significantly higher compared to the case in which no projection was conducted, while for the practically important case in which data were smoothed from irregular sparse and noisy observations ($M4$) and jumps were truncated at level $l=3$, the error more than doubled. For all considered thresholds of explained variation ($85\%$, $90\%$, $95\%$ and $99\%$) the reported dimension is constantly $3$, where we just reported the results for the $99\%$ threshold in the table. This is not surprising, since we started from a three-factor model for the log bond prices, but it demonstrates, that the common practice of projecting price or yield curves onto a few linear factors can disguise statistically important information, despite their high explanatory power for the variation of log bond prices, which is in line with Crump2022.

Empirical analysis of bond market data

In this section, we apply our theory to bond market data. In particular, we investigate the influence of jumps on the estimators and the dimensions of the integrated volatility, that is, the continuous part of the quadratic covariation, to determine how many random drivers are statistically relevant.

We consider nonparametrically smoothed yield curve data from FPY2022. For constructing smooth curves on each day, the authors of FPY2022 use a kernel ridge-regression approach based on the theory of reproducing kernel Hilbert spaces. We measure time in years and the data are available for approximately $250$ trading days in each year, yielding $n\approx 250$, using a day count convention in trading days, with a daily resolution in the maturity direction, where we consider a maximal time of $M=10$ years to maturity. The data are given as yields, which we first transform to zero coupon bond prices and then derive the difference returns $d_i^n(j)$ for $i=1,...,\ensuremath{\lfloor T/\Delta_n\rfloor}$ and $j=1,...,\lfloor M/\Delta_n\rfloor$ by formula (ref). We consider data from the first trading day of the year $1990$ ($i=1$) to the last trading day of the year $2022$ ($i=33$). We then derive the estimators $\hat q_{i}^n\big|_{[0,10]},$ and $\hat q_{i}^{n,10,-}$ (as defined in Remark (ref)) for $i=1,...,33$ and derive the yearwise covariation kernels $$\hat q_{i}^{*,-}:=\Delta_n^{-2}(\hat q_{i}^{n,10,-}-\hat q_{i-1}^{n,10,-})\quad \text{ and }\quad \hat q_{i}^*:=\Delta_n^{-2}(\hat q_{i}^n\big|_{[0,10]}-\hat q_{i-1}^n\big|_{[0,10]})$$ for $i=1,2,...,33$ and $q_{0}^n=0$. The data-driven truncation rule described in Section (ref) is applied for a preliminary truncation at the $0.75$-quantile of the sequence of difference return curves as measured in their $l^2$ norm and $d$ is chosen as the smallest value such that $d$ eigencomponents explain 90% of the variation of the preliminary estimator. Importantly, the truncation rule is conducted for different $l=3,4,5$ and for each year separately and only takes into account data within the respcetive year. We also consider an estimator for a potential long-term volatility given by $$\hat q_{long}^*:= \frac 1{33}\sum_{i=1}^{33} q_i^{*,-}.$$ Under the Assumption of Section (ref), this is an estimator for a stationary volatility kernel. The results suggest that quadratic covariations in each year are rather complex in the sense that they exhibit a slow relative eigenvalue decay, unveil a varying shape and magnitude over time, and often differ quite substantially from measured quadratic covariations due to the existence of jumps. Subsequently, we provide a thorough discussion of these observations. A table containing all results of the analysis is contained in the appendix

Impact of jumps

On one hand, jumps that have a moderate impact on the magnitude of the overall quadratic covariation can visually distort the shape of the volatility. Figure (ref) depicts plots of the graphs of the estimated truncated kernels $q_i^{*,-}$ (with the truncation level $l=3$) and nontruncated kernels $q^*_i$ for the years 2005, 2006 and 2007. In 2006, which is also a year in which the yield curve inverted before the financial crisis in 2007, two jumps had a visible impact on the shape of the measured quadratic covariation kernels although they together accounted for less than 3 % of the magnitude of the quadratic covariation. This is due to a higher emphasis on the variation in difference returns with short maturities where one should note the different scalings in the plots. Removing these two jumps leads to a more time-homogeneous shape of the integrated volatility surfaces in the sense of the relation of the variation in the shorter maturities to the variation in higher maturities. On the other hand, jumps influence the magnitude of the quadratic variation. For instance, nine increments in the year 2020 (Covid-19 outbreak) sorted out by the truncation rule for $l=3$ accounted for more than 50% of the overall quadratic variation in the data as measured by its norm. The statistics for jumps in all years can be found in the supplment to this article. Interestingly, our measurements suggest that jumps tend to cluster.

Dimensionality of the continuous part of the quadratic variations

We now examine the statistically relevant number of random processes that are driving the continuous part of the forward curve dynamics by investigating the dimensionality of the continuous quadratic covariation of the latent driving semimartingale via the estimator $\mathcal T_{\hat q_i^{*,-}}$.

Table (ref) reports the number of eigenfunctions of $\mathcal T_{\hat q_i^{*,-}}$ that are needed to explain resp. $85\%$, $90\%$, $95\%$, and $99\%$ of the continuous covariation in the years 2005, 2006 and 2007 showing that to explain 99% of the variation at least 10 factors are needed in each year. The situation looks similar for all other years from 1990-2022, while the detailed results were relegated to the appendix. We find that the complexity of the covariation seems to have decreased over the years, indicating a time-varying pattern of the volatility term structure that goes beyond its overall level. It is noteworthy that in almost every year (30 out of 33), the number of linear factors needed to explain at least 99% of the truncated variation of the data is at least $10$. These dimensions even increase if we employ the static factors $(e_j^{long})_{j\in \mathbb N}$ of eigenvectors of $\mathcal T_{\hat q_{long}^*}$ and do not update them in each year. In that case, in 27 out of 33 years at least 12 factors are needed to explain at least $99\%$ of the variation in each year.

Importance of higher-order factors for short term trading strategies

A natural question is if the higher-order factors indicated by the analysis of real bond market data in Section (ref) are of economic significance beyond capturing variation in difference returns. Therefore, we investigate whether the high dimensionality of the continuous quadratic varitations indicated by the estimators $\hat q_i^*$ and $q^*_{long}$ are important for other short term trading strategies than difference returns. Precisely, Define the daily return $(d_L)_{i}^n(j)$ of the trading strategy of buying an $(j+L)\Delta_n$ bond and shorting an $j\Delta_n$ bond $$(d_L)_{i}^n(j):= \sum_{l=0}^{L}d_i^n(j+l)=\tilde\Delta^n_{i\Delta_n} \log P((j+L)\Delta_n)-\tilde\Delta_{i\Delta_n}^n \log P(j\Delta_n)$$ where $\tilde\Delta_t^n \log P(x)=\log P_{t+\Delta_n}(x)-\log P_{t+\Delta_n}(x-\Delta_n)$. Evidently, we can derive them as linear functionals of either log bond price returns or difference returns.

We want to determine the adequacy of approximation of these higher-order difference returns when they are derived either from approximated log price curves, which are projected onto its leading principal components or when they are derived from difference returns, which are projected onto the leading eigencomponents of the long-term volatility estimator. For that, we calculate the relative mean absolute error ($RMAE$) for a set $\mathcal V=\{i_1 \Delta_n^{val},...,i_{825}^{val}\Delta_n\}$ of dates (where in each year from 1990 to 2022 we randomly draw 25 dates making a total of 825 validation dates). That is, defining the piecewise constant kernels $(\tilde d_L)_{i}^n=\sum_{j=1}^{\lfloor M/\Delta_n\rfloor}(d_L)_i(j)\ensuremath{\mathbb{I}}_{[j\Delta_n,j\Delta_n)}$ we calculate $$RMAE_L(f_1,...,f_d):=\frac 1{825}\sum_{l=1}^{825}\frac{\left\|(\tilde d_L)_{i_l}^n-\mathcal P_{f_1,...,f_d}(\tilde d_L)_{i_l}^n \right\|_{L^2(0,10)}}{\left\|(\tilde d_L)_{i_l}^n\right\|_{L^2(0,10)}}. $$ where $\mathcal P_{f_1,...,f_d}:= \sum_{i=1}^d f_i^{\otimes 2}$. The factors are derived in two different ways. In the first scenario (S1), the factors $f_1,...,f_d$ correspond to principal components of the empirical covariance of log-price differences $\tilde \Delta \log P_{i\Delta_n}$ for $i\notin \mathcal V$ and in the second scenario (S2) the factors $f_1,...,f_d$ correspond to the leading eigenfunctions of the estimated stationary volatility kernel $\hat q_{long}^*$ where as before the truncation of jumps is conductcted yearwise with truncation level $l=3$ according to the truncation procedure described in Section (ref).

We compare the results for lags of 7, 30, 90 and 180 days, since they approximately correspond to the returns of buying a bond and shorting another bond with time to maturity that is resp. a week, a month, a quarter or half a year higher. The $RMAE$s can be found in Table (ref).

It can be observed that a high number of factors is needed to approximate the lagged difference returns precisely and that approximations based on low factor structures as indicated by the covariance of difference returns imply high approximation errors. While it is not surprising that the approximation gets better if we use more factors, the high discrepancy of the approximation errors is noteworthy. The errors for a typically chosen three factor model based on log price differences (the factors correspond to level, slope and curvature), which explain more than 99,7% of the variation in log-price returns is for all lags higher than $0.4$, whereas for the approximation error for 14 factors, which we would need to explain 99% of the variation in difference returns as measured by $\hat q_{long}^*$ is never higher than $0.11$. Interestingly, for all lags, choosing the factors equal to the leading eigenfunctions of the long term volatility $\hat q_{long}^*$ instead of the ones indicated by log-price differences can reduce the error for the higher-order approximations quite significantly and for $d=16$ and for lags not higher than 90 days by almost 50 %. Higher-order factors of volatility can, thus, not easily be ignored and might carry important economic information.

Concluding remarks on the empirical study

We conclude that the reported dimensions are overall quite high compared to the few factors needed to explain a large amount of the variation in yield and discount curves. This suggests that low-dimensional factor models are not able to capture all statistically relevant codependencies of bond prices. Still, exact magnitudes of explained variations of the higher order components have to be interpreted cautiously and conditional on the smoothing technique that was employed to derive yield or discount curves. However, higher-order factors seem to be economically relevant for capturing variations in short term trading strategies as indicated by the out-of-sample study of Section (ref).

Underestimation of the number of statistically relevant random drivers can have undesirable effects. For instance, Crump2022 showcase the potential economic impact on mean-variance optimal portfolio choices and hedging errors. At the same time, not every model that is parsimonious in its parameters needs to entail a low-dimensional factor structure such as the simple volatility model of Section (ref). It seems desirable to derive parsimonious models that match the empirical observation of high or infinite-dimensional covariations and reflect the characteristics of their dynamic evolution.

Acknowledgements

I would like to thank Dominik Liebl, Fred Espen Benth, Alois Kneip and Andreas Petersson for helpful comments and discussions. Funding by the Argelander program of the University of Bonn is gratefully acknowledged.