Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
34,442 characters · 11 sections · 38 citation commands
Inference for Model Misspecification in Interest Rate Term Structure using Functional Principal Component Analysis
\doublespacing
The dynamic evolution of the interest rate term structure is indispensable for many tasks, including pricing interest rate derivatives, discounting future payoff for asset prices, explaining return predictability, conducting firm evaluations, managing financial risk, making monetary and fiscal policy elenev2021can,elenev2022can,krishnamurthy2011effects,woodford2012methods, and most importantly, forecasting rare events like financial crises and predicting future growth hu2013noise,koijen2017cross.
There seems to be a linear structure in the interest rate term structure most of the times. Eigendecomposition of the covariance matrix of yield changes attributes over 97% of the variance to three principal components, known as the level, the slope, and the curvature, seemingly universal across different units of maturity chapman2001recent,litterman1991common,piazzesi2010affine. Using the daily zero-coupon yield curve of bonds maturing in 1 year to 30 years gurkaynak2007us from November 25, 1985\footnote{The data starts in 1985 after dropping the missing and incomplete data, so that all cross-sectional data points have observations from hypothetical securities with 1-year to 30-year maturity. The missing data before was mainly due to lack of actual securities for model calibration.}, to December 9, 2022, we obtain the first three principal components of the entire dataset in Figure (ref).
On the other hand, a dispersion in actual yields around the smooth term structure has been used to measure liquidity, whose friction spikes in crises while remaining small otherwise hu2013noise. Following fama1973risk regression and construction of a factor mimicking portfolio ang2006cross, a significant liquidity premium is revealed so that the noise measure passes the cross-sectional test. However, the decent returns gained from trading illiquidity do not necessarily suggest that the interest-rate model is well-specified, nor does it justify using a smooth three-factor term structure model as a reference model (as in hu2013noise) from which dispersion is measured.
To put it in another way, even if the first three principal components can explain a high variance in the yield curves (which also means that reconstruction error is low), as shown in Figure (ref), and even if the full-sample stylistic features (in Figure (ref)) of first three principal components have an immaculate interpretation of their shapes, the three principal components do not always correspond to the level, slope, and curvature factors from a more granular perspective - in some years they may not have rich economic interpretations at all lord2007level. Indeed, from Figure (ref), the level-slope-curvature seems to be occasionally violated, as sometimes the first principal component (which “ought to be" the level) is not flat (e.g., in 2007 and 2017), and the second principal component (which “ought to be" the slope) is not monotonic (e.g., in 1990, 2008, and 2021).
With such economic interpretations being invalidated occasionally, some standard models of interest rate term structure based on these interpretations, such as nelson1987parsimonious, might problematically impose a specific smooth structure of yield curves, and the dispersion from a smooth yield curve that hu2013noise benchmarks against could be a biased benchmark.
As pointed out in previous literature, model misspecification (also known as model uncertainty or ambiguity) should be acknowledged in macroeconomic theory hansen2001acknowledging,hansen2001robust; market participants probably have a preference towards robustness anderson2000robustness,hansen1999robust,maenhout2004robust; the monetary policies could have differed a lot when the actual model describing the environment appears ambiguous to the central bank giannoni2002does,hansen2010wanting,kasa2002model,onatski2002robust,woodford2010robustly. It is also likely that interest rate models also suffer from misspecification brenner1996another, mainly when market crashes or jumps are observed in continuous time johannes2004statistical. Despite the overall satisfactory exponentially affine term structure models, the moments of latent variables present in the affine structure could also be ambiguous liu2018return,liu2022volatility,zhang2022optimal, leading to a misspecified affine model.
In light of hu2013noise, this paper examines the model misspecification of interest rate term structures both on an omniscient scope and a granular perspective. Figure (ref) presents the functional data structure we have, i.e., the time series of term structure as a function, which motivates our methodology different from the inadequate conventional one. We utilize recent advances in functional data analysis, particularly a nonparametric hypothesis test of functional principal component analysis (FPCA), and compare its results implied by classical principal component analysis (PCA), which highly relies on the homogeneity of term structure throughout time. Instead of directly constructing a stochastic partial differential equation to model the interest rate dynamics as a stochastic process in the space of functions cont2005modeling, we employ the nonparametric and distribution-free method in this research to lean on minimal assumptions.
Our codes are released on \url{https://github.com/Kaiwen-Hou-KHou/termStructureFPCA}.
In this section, we introduce the essential statistical and econometric tools to deal with functional data such as the time series of yield curves, including a generalization of classical principal component analysis and eigendecomposition into infinite-dimension space $L^2([0,1],d\xi)$ (the space of square-integrable functions on $[0,1]$, where the inner product weighted by the constant function $1$ is defined as $\langle f,g \rangle := \int\limits_0^1 fg d\xi$), the assumption of smoothness and smoothing techniques, followed by some asymptotics that helps examine FPCA method through statistical inference.
A generalized PCA method, FPCA is developed to investigate the dominant modes of functional data variation, such as modeling sparse longitudinal data or conducting functional classification and regression hall2006properties,li2010uniform,yao2005functional,yao2005functional2.
Let $y_i$ be i.i.d. samples of a random function $y: (\Omega, {\mathcal{F}}, P)\to (C^0[0,1], {\mathcal{B}}_{C[0,1]})$ with $\mathbb{E}||y||^2<\infty$. Denote by $$ m=\mathbb{E} y\in L^2([0,1],d\xi) $$ the mean function, and by $$ G(x,x')=Cov(y(x),y(x'))\in L^2([0,1]^2,d\xi d\xi') $$ the covariance function\footnote{Since covariance is symmetric (and continuous here), then $G$ is a positive semi-definite kernel, so Mercer's theorem applies.}. By Mercer's theorem, there exists a sequence of eigenvalues $\lambda_k \downarrow 0$ such that $\sum\limits_{k=1}^\infty \lambda_k <\infty$ and $$ G(x,x') = \sum\limits_{k=1}^{\infty}\lambda_k\psi_k(x)\psi_k(x'), $$ where the orthonormal basis functions $(\psi_k)_{k=1}^\infty$ are the eigenfunctions satisfying $$ \int\limits_0^1 G(x,x')\psi_k(x')dx' = \lambda_k \psi_k(x). $$ By the Kosambi–Karhunen–Lo\`eve theorem kosambi2016statistics, one can write $$ \xi_i = y_i - m = \sum\limits_{k=1}^\infty \zeta_{ik}\psi_k, $$ where the principal component\footnote{$\zeta_{ik}$ is also known as the rescaled functional principal component scores.} associated with the $k$-th eigenfunction $\psi_k$ is $\zeta_{ik} = \langle \xi_i, \psi_k \rangle$ with the following moment properties
for all $1\leq i\leq n,$ and $ k\neq k'\in{\mathbb{N}}^+$.
It is also justified in cont2005modeling and kusuoka2000term\footnote{This is a work suggested by Professor Darrell Duffie.} that the term structure could be modeled to have a two-dimensional dynamics (time and maturity), even if the dynamics on maturity is slightly anti-intuitive. Therefore, the interest rate term structure fits in this stochastic framework.
Note that a crucial assumption of FPCA is that data are modeled as sample paths of smooth stochastic processes, although they are often observed discretely and noisily. If sample paths are not smooth, then $G\notin L^2([0,1]^2,d\xi d\xi')$. Fortunately, this assumption can be relaxed by restricting $G$ to a subset to model some diffusion processes. For example, standard Brownian motion has non-differentiable covariance $G(x,x')=\min(x,x')$ on the diagonal, but it becomes $G(x,x')=x$, which is infinitely differentiable when restricted on the upper-triangular block
With similar restriction, the covariance functions of common diffution processes such as Brownian bridge, geometric Brownian motion, and Ornstein-Uhlenbeck process jouzdani2021functional would also lie in $L^2([0,1]^2,d\xi d\xi')$.
For discrete observations, methods like local linear or spline smoothing can be used to obtain the functional form of each $y_i$ yao2005functional. In particular, if observations are evenly-spaced at $\left(\frac{j}{N}\right)_{j=1}^N$, one can adopt the B-spline estimator by projecting the discrete-valued function $y_i$ to the space ${\mathcal{H}}^{(p-2)}[0,1]$ spanned by B-splines of order $p$:
Correspondingly, the smoothed and centered process is estimated by $$ \hat{\xi}_i = \hat{y}_i - \hat{m} = \hat{y}_i - \frac{1}{n}\sum\limits_{j=1}^n \hat{y_j}. $$ Getting rid of the diagonal components in $G$ by restrictions in Eq. (ref) also ignores the measurement errors after smoothing staniswalis1998nonparametric.
In dealing with interest rate term structures, their smoothness has been a reasonable assumption as widely used in interpolating yield curves, which rely on piecewise smooth functions that are also smoothly joined at selected knots christensen2011affine,diebold2006forecasting,diebold2008global,fisher1995fitting,mcculloch1975tax,nelson1987parsimonious,svensson1994estimating.
Given an orthonormal set $(\psi_{0k})_{k=1}^{\infty}$, song2022hypotheses propose a statistical test for the hypothesis $H_0$ that the eigenfunctions $(\psi_k)_{k=1}^\infty$ form the same set up to permutations and change of signs. Fixing a hyper-parameter $\kappa_n$, the test statistic is
where $$ \hat{Z}_{kk'} = \left\langle \hat{G}(x,x'), \psi_{0k}(x)\psi_{0k'}(x') \right\rangle_{L^2([0,1]^2,d\xi d\xi')} := \int\limits_0^1\int\limits_0^1 \hat{G}(x,x')\psi_{0k}(x)\psi_{0k'}(x') dxdx', $$ with the two-step estimate of the covariance
Note that $\hat{G}$ is essentially a low-rank approximation of the bi-infinite covariance matrix, as $y_i$ has been projected onto ${\mathcal{H}}^{(p-2)}[0,1]$. This na\"ive sub-block estimator is genuinely a choice out of simplicity; one could also adopt stochastic estimation such as the Russian roulette estimator to attain unbiasedness hou2022spectral.
It has also been shown in song2022hypotheses that the asymptotic distribution of the test statistic under $H_0$ is an infinite Gaussian quadratic form, and a finite sample estimator of its quantile $\hat{Q}_{1-\alpha}$, which is the $100(1-\alpha)$-th percentile of
is consistent. Here $$ \hat{\lambda}_k = \frac{1}{n}\sum\limits_{i=1}^n \hat{\zeta}_{ik}^2, $$ where $$ \hat{\zeta}_{ik} = \left\langle \hat{\xi}_i,\psi_{0k} \right\rangle. $$ The distribution of the Gaussian quadratic form $\bar{S}_n$ in Eq. (ref) is essentially a weighted sum of $\chi^2(1)$ random functions, where the weights are the eigenvalue pairs of $G$. In the special case where these eigenvalues are 1, the distribution becomes $\chi^2\left( \frac{\kappa_n(\kappa_n-1)}{2}\right)$, but this does not hold in general, as the scale parameter of a gamma distribution changes when multiplied by a constant and $\hat{\lambda}_k\hat{\lambda}_{k'}\chi_{kk'}^2(1)$ might be no longer of Chi-squared distribution. Unfortunately, $\bar{S}_n$ does not have a neat closed-form representation of its density or distribution function, so its quantiles should probably be computed numerically bodenham2016comparison,davies1980algorithm.
The Nelson-Siegel model nelson1987parsimonious fits the yield curve with a simple functional form
where the parameters are $\Theta=(\beta_0, \beta_1, \beta_2, \theta)$ and $\tau$ is the maturity in months\footnote{The model is invariant with respect to the unit of maturity, as $\theta$ absorbs the scale difference between units. For example, transforming $\tilde{\theta}=12\theta$ and $\tilde{\tau}=\frac{\tau}{12}$ in years, Eq. (ref) is of the same form.}.
It is acknowledged (e.g., by diebold2006forecasting,diebold2008global) that $\beta_0,\beta_1$, and $\beta_2$ represent loadings on the level factor, the slope factor, and the curvature factor, respectively, so Eq. (ref) can be generalized to the dynamic Nelson-Siegel Model in the following Eq. (ref)
through adding the auto-regressive state variables $L_t$ for level, $S_t$ for slope, and $C_t$ for curvature.
We now proceed to test whether the orthonormal basis beneath the three factors $1, \frac{1-e^{-\theta\tau}}{\theta\tau}$ and $\left( \frac{1-e^{-\theta\tau}}{\theta\tau} - e^{-\theta\tau} \right)$ is equivalent to the actual principal component basis to see if the level, slope, and curvature factors always explain the major variation in the interest rate term structures, and offer insights on the dispersion between the two bases.
Typically the tradable bond maturity in months is bounded both from above and from below, i.e., $0<\underaccent{\bar}{\tau}\leq\tau\leq\bar{\tau}$. We can reproduce Eq. (ref) in the following theorem so that $\tau$ is rescaled in $[0,1]$.
For the closed-form expressions of the orthonormal basis $(\psi_{0k})_{k=1}^3$, refer to Proof of Theorem (ref).
Letting $\theta$ vary from a reasonable range, we observe from Figure (ref) that $\psi_{01}$ is always the level, $\psi_{02}$ the slope, and $\psi_{03}$ the curvature. Furthermore, with an increase in positive $\theta$ increases, the slope and the curvature change more drastically at short-maturity bonds; a decrease in negative $\theta$ pushes the changes of the slope and curvature towards long-term yields. Hence, just like the coefficients $(\beta_k)_{k=0}^2$ in Eq. (ref), the coordinates in the new basis $(\psi_{0k})_{k=1}^3$ can still be interpreted as the loadings on the level/slope/curvature factors, if they truly are these factors.
If we view $(\psi_{0k})_{k=1}^3$ as the factors with the most significant loadings among an extended orthonormal set $(\psi_{0k})_{k=1}^{\infty}$, we now proceed to test if $(\psi_{0k})_{k=1}^{\infty}$ is equivalent to the actual eigenfunctions $(\psi_k)_{k=1}^\infty$ up to permutations and changes of sign, i.e., whether the orthonormal basis underlying the Nelson-Siegel model is well-specified to match the actual principal component basis of the covariance kernel $G$.
We discretize $x$ and $x'$ on an $N\times N$ grid and estimate a rank-$N$ approximation of the bi-infinite covariance matrix in Eq. (ref). The entries of the low-rank sample covariance matrix are shown in Figure (ref) for the entire data. The positive entries demonstrate the strong cross-sectional co-movements of yield curves over time. Such co-movements are much stronger on short-term bond yields as concentrated at the upper-left block of Figure (ref), which apparently have higher liquidity and respond to the market more effectively than long-term bond yields.
However, such a statement is only true when there are no significant frictions in the market. Notice the estimates $\hat{G}$ are heterogeneous in different years shown in Figure (ref). Unlike hu2013noise's noise measure of liquidity and arbitrage capital that equally weighed the pricing errors of bond yields with different maturities, we argue that the cross-sectional difference in covariance is also indicative for different sources of liquidity. In 1986, for example, the short-term covariation and long-term covariation were both high, while the covariance between short-term yields and long-term ones stayed low, which indicates that the frictions that drove the illiquid yields might have some components uncorrelated with everyday frictions that mainly impacted the short-term yields, and further gave insights on the 1987 market crash. For another instance, a different but unique covariance structure was observed in both 1992 and 2008, representing the strong and correlated co-movements on both short-term and long-term yields, which could be impacted by factors from the external financial market which tilted the yield curves at their most extreme values, such as the exchange rate risk exposure in the 1992 UK currency crisis and the overall exposure to the global financial crisis in 2008.
To show the interpretations above are robust with respect to the choice of B-spline subspace and the number of uniform knots $N$, so that the projection in Eq. (ref) is not overfitting the yield curves, we demonstrate highly similar patterns revealed by Figure (ref) using linear B-spline basis and by Figure (ref) where $N=20$. This also justifies the robustness of the estimation of covariance function $G$.
Based on the covariance estimate $\hat{G}$, we first straightforwardly implement FPCA, the results given in Figure (ref). The shapes of the first three functional principal components and their explained variance are almost identical to the PCA results in Figure (ref).
To better illustrate the effects of the obtained functional principal components, we add and subtract the components to the mean function $\hat{m}$. We can then observe that the first component (FPC1) shows the parallel shift in the mean $\hat{m}$, the second component (FPC2) changes the slope of yield curves by tilting around one knot, the third (FPC3) embeds rich curvature by tilting the yield curve around two fixed knots, FPC4 tilts around 3 points, and so on so forth. Thus, similar interpretations on the level, slope, and curvature to PCA could be given to the functional principal components. Figure (ref) offers an intuition of how the flexibility of yield curves can be traced back by their functional principal components. Formally, such intuition is a natural result of the following theorem.
Eq. (ref) is also the svensson1994estimating specification of yield curves on all maturities, also used in our data after 1980\footnote{Prior to 1980, the Nelson-Siegel model with fewer parameters was used to fit the yield curve, as there were not enough Treasury securities to fit the Svensson model, which has six parameters.} gurkaynak2007us and widely accepted as smoothed yields on hypothetical Treasury securities\footnote{Refer to the nominal yield curve (\url{https://www.federalreserve.gov/data/nominal-yield-curve.htm}).}. Regarding the parameters, $\theta_i=\frac{1}{\tau_i}$ for $i=1,2$, and $(\beta_k)_{k=0}^3$ are the same as the Svensson model.
Theorem (ref) also implies that if $|Ker(\psi_k)|$ deviates a lot from $k-1$, either the regularity conditions imposed on $\epsilon$ fail to be satisfied, or the sample size is not large enough, or the data-generating process is misspecified. We argue that the former two scenarios are explicitly testable, both involving tests (whether asymptotic or not) on the distribution of $\epsilon$. So it would be possible that the model is misspecified after eliminating the former two scenarios. For instance, it might not be appropriate to name the corresponding functional principal components in Figure (ref) according to level, slope, and curvature anymore, and the fact that $|Ker(\psi_k)|\neq k-1$ for $k=2,3$ is evidence in favor of model misspecification.
We now formally compute the test statistic $\hat{S}_n$ and simulate $\bar{S}_n$, the estimated asymptotic distribution of the test statistic under $H_0$, in the full sample. The hyper-parameter $\kappa_n=3$, since the eigenvalues $\lambda_k \ (k>3)$ are negligible under $H_0$, i.e., $(\psi_{0k})_{k=1}^3$ given in Theorem (ref) already explains most variation in the yield curves.
Comparing Figure (ref) and Figure (ref), we see the score $\zeta_{i1}$ on the first proposed eigenfunction $\psi_{01}$ captures the macroscopic variation in yield curves of all maturities, as a natural result of the definition of functional principal component analysis.
Directly applying the test in song2022hypotheses for a set of possible parameters $\theta\in\{\pm 0.001, \pm 0.01, \pm 0.1, 0.2, 0.5, 1\}$ for the Nelson-Siegel orthonormal basis, $H_0$ is rejected both on the entire sample and on almost each yearly data with $p << 0.05$ shown in Figure (ref). This means that the Nelson-Siegel basis in general cannot be interpreted as the principal components, and there is strong evidence in the model being misspecified.
To further gain some insights on how confident we are when regarding the Nelson-Siegel orthonormal basis as the principal components, we loosen the above test statistic in Theorem (ref).
Now we proceed to conduct the hypothesis test in Theorem (ref) on the full sample. The $p$-values reported in Table (ref) suggest that the Nelson-Siegel basis $(\psi_{0k})_{k=1}^3$ is not highly sensitive enough to the parameter $\theta$ to influence the test conclusion within a reasonable range (e.g., $|\theta|\in[0.01, 0.2]$). Hence there is probably rich flexibility embedded by $\theta$ in the Nelson-Siegel model.
However, the heterogeneity of the test results each year in demonstrated by Figure (ref). We notice that when the highest $p$-values attained are low (so that there is evidence in favor of $H_1$, and the level-slope-curvature argument is suspicious), there was a recession or some influential adjustments on the interest rate policies. In 2014, for instance, the Federal Reserve ceased its expansion in holding longer-term securities through open market purchases, originally aiming at putting downward pressure on longer-term interest rates and supporting economic activities through making more accommodating financial conditions. 2007 marked that the Federal Reserve cut the target on short-term interest rate by half of a percentage point. The recession and RTC in 1991, and the UK currency crisis in 1992 causes explainable deviations to the reference model as well. These are all possible shocks that lead to interest rate model deviating from the level-slope-curvature component structure.
This paper investigates the model misspecification of interest rate term structure using functional principal component analysis. We first show by orthonormalization in the function space that the Nelson-Siegel model, as a standard reference model, can be interpreted as a three-factor interest rate term structure involving the level, the slope, and the curvature. This justifies the presence of the Nelson-Siegel orthonormal basis in the null hypothesis to test whether the actual data-driven principal components coincide with the level/slope/curvature factors, i.e., whether (and when) the three-factor model is misspecified.
Two notions have been used to determine whether two bases are the same. The stricter one in Section (ref) offers a non-rigorous way to test whether the functional principal components match the level, the slope, and the curvature in the exact order. We conclude from the entire sample that the components appear in the desired order, which is not the case for specific individual years. For instance, in 2008, a curvature factor potentially explained more variance than a slope factor. The same conclusion is drawn based on the other less strict notion in Section (ref), where permutations and sign changes in the basis are accepted. To sum up, this paper characterizes the heterogeneity of misspecification in the interest rate term structure when rare events or relevant shocks in policies occur, even if the model is accepted overall.
Future research could generalize the proposed methodology to test the reference models other than the Nelson-Siegel model, which lacks financial intuition compared to other affine models such as vasicek1977equilibrium, and to include calibration from the actual market data in bond yields rather than merely hypothetical securities.
\onehalfspacing \setlength\bibsep{0pt}
\doublespacing