Approximate Factor Models for Functional Time Series

Approximate Factor Models for Functional Time Series

\doublespacing

{\it Keywords:} Curve data, functional data analysis, identification, information criterion, \\ forecasting, yield curve modeling

{\it JEL Classification:} C32, C38, C55, E43

Introduction

Over the past decades, approximate and dynamic factor models have emerged as powerful tools for analyzing high-dimensional datasets, offering a parsimonious framework that captures the essential underlying structure while filtering out irrelevant parts. These models, first conceptualized by chamberlain1983 and expanded upon by forni2000, stock2002, stock2002b, and bai2003, have become fundamental tools in areas such as economic forecasting, monetary policy analysis, psychology, environmental and social sciences. For comprehensive reviews, see bai2008, breitung2013, and stock2016.

Building on the success of these models in handling high-dimensional data, the generalization of factor modeling techniques to functional or infinite-dimensional data structures has gained significant attention. This transition to functional data analysis (FDA) is a natural extension, as this field focuses on datasets consisting of continuous curves or functions, as discussed extensively in ramsay2005, horvath2012, hsing2015, and kokoszka2017. FDA is particularly valuable in economic applications, where data often take functional forms, such as term structures of bond yields, credit default swaps, income profiles, and inflation expectations, for which conventional multivariate methods may prove overly restrictive or inadequate. In this paper, we contribute to the literature on factor models for functional data by proposing a novel approximate factor model that explicitly models the dynamics of time-dependent functional data (or functional time series).

The study of factor models in FDA has seen significant progress in recent years. Early contributions focused on generalizing classical factor models to the functional setting. For instance, hays2012 and liebl2013 introduced functional factor models with a discrete error component, while hormann2022, hormann2023 and ofner2024 demonstrated how discretely observed functional data follow an approximate factor model structure. These works primarily utilized the functional counterpart of principal component analysis (PCA), building on directions of maximal variability for factor identification. The same PCA-based approach has been extended to functional time series. For example, hyndman2009 and aue2015 applied PCA techniques to model and forecast functional time series, while further extensions have addressed challenges related to nonstationarity and long-range dependence (see, e.g., Chang2016, Li2020, and salish2019). More recent contributions, such as hallin2023 and tavakoli2023, have refined these methods, identifying factors asymptotically through panel structures of large functional time series datasets. In addition, Bayesian methods have also offered a rich set of tools for modeling factors in functional datasets both cross-sectional and time series (e.g., Montanga2012, kowal2017, and Kowal2023).

The identification and estimation of factors in functional data analysis typically relies on selecting directions of highest variability, typically obtained through the functional adaptation of PCA. While this approach effectively reduces dimensionality, the resulting factors do not necessarily capture all temporal dependencies or reflect predictability, both of which are integral to functional time series analysis. As argued by forni2000, panaretos2013 and hormann2015, focusing solely on variability can lead to an incomplete representation of the data's dynamic structure, potentially overlooking important serial dependencies. Our paper addresses this limitation by complementing and extending existing methodologies to develop a factor model tailored specifically for functional time series. In this framework, factors are identified based on their contribution to the dynamic and predictive structure of the data, rather than directions of maximal variability.

Similar to existing factor models, our approach decomposes a functional time series into two distinct components: (i) a low-dimensional factor component and (ii) a remainder term, referred to as an error component. The key distinction lies in how the factor component is derived in our framework. Rather than relying on directions of maximal variability identified through PCA, our method extracts the factor component from the autocovariance structure of the functional time series to reflect serial dependencies and predictability. This departure from conventional factor analysis enables a better characterization of the temporal dynamics inherent in functional time series. Our methodology is related to those proposed by pena1987, bathia2010, lam2012, and diks2023. Consequently, the factors and error components in our model can also be interpreted as predictable and unpredictable, respectively, reflecting their roles in characterizing the process dynamics. The term “approximate” in the model’s name stems from the allowance for correlations between values at different domain points of the error term, consistent with the terminology introduced by chamberlain1983.

Our construction also changes the role and structure of the error component. The primary requirement to the error component is to not carry any relevant information about dynamics of the original functional time series. This, in turn, allows us to relax the typical restrictions imposed in conventional factor models on the error term, allowing for more flexibility. First, the eigenvalues of both the factor and error components can be of comparable magnitude. Second, as an approximate factor model, it permits non-negligible correlations across different domain points of the error function, even asymptotically. This stands in contrast to conventional models, which assume lower magnitude errors and only weak correlations, crucial for separating the factor component from the error one. Finally, in our framework, the factor and error components may exhibit weak cross-correlations, accommodating various forms of nonstationarity and heteroskedasticity.

Following the development of this new factor modeling framework for functional time series, the second main contribution of this paper is to ensure both theoretical rigor and practical applicability. We develop a consistent estimation procedure for the model's primitives, such as the autocovariance structure, and address the challenge of estimating the factor component’s dimension, which is typically unknown and must be inferred from the data. To tackle this challenge, we propose an information criterion, drawing from traditional solution in factor modeling literature (e.g., bai2002; hallin2007). Our criterion is based on the prediction error curve and assumes that the factor component follows a vector autoregressive (VAR) process with an unknown number of lags. We incorporate a carefully designed penalty term in the criterion to mitigate over-selection, ensuring consistent estimates for both the number of factors and the lag order under mild theoretical assumptions. All proposed estimation procedures are available in our R-package\footnote{https://github.com/ottosven/dffm}.

Finally, to illustrate the practical utility of our model, we apply it to mortality and yield curve modeling and forecasting. In both applications, we compare our approach with functional PCA, which has arguably become the main applied workhorse model in functional time series analysis since the seminal work of aue2015. Additionally, for yield curves, we benchmark our model against the widely used dynamic Nelson-Siegel model, as reviewed in diebold2013. Our main findings highlight the importance of a careful structural representation of process dynamics -- that is, reducing the dimensionality of the original process in a meaningful way by extracting factors that drive the dynamics rather than merely capturing variability. For instance, in the case of mortality rates, while the first three factors in our framework closely resemble those identified by PCA, the remaining factors differ. This distinction proves important, as selecting factors based on their role in driving dynamics rather than explaining variability leads to improved forecasting performance. In the case of yield-curve modeling, our results indicate that the number of factors required to capture yield curve dynamics varies over time, with more factors needed during periods of economic uncertainty than in stable periods. These findings suggest that yield curves have a richer and more adaptive dynamic structure than typically considered in the literature.

The paper is structured as follows: Section (ref) develops the approximate functional factor model, providing detailed discussions on all necessary assumptions and the identification of its components. In Section (ref), we introduce the estimator for the functional component, discuss the information criterion for jointly estimating the number of factors and their dynamics, and show their consistency. We also provide guidance on practical implementation. Section 4 presents a Monte Carlo simulation to assess the model's performance and the proposed estimation methods in finite samples. In Section 5, we apply the method to mortality rate and yield curves. Finally, Section 6 concludes the paper.

To ease the readability of the paper, we collect the key notations below. Consider $H=L^2([a,b])$, the space of square-integrable functions $f:[a,b] \rightarrow \mathbb{R}$ satisfying $\int_a^b f(r)^2 dr < \infty$ for $a < b$. $H$ forms a Hilbert space equipped with the inner product $\langle f, g \rangle = \int_a^b f(r)g(r) dr$ and the norm $\|f\| = \sqrt{\langle f, f \rangle}$, for $f,g \in H$. Any square-integrable kernel function $\rho: [a,b] \times [a,b] \rightarrow \mathbb{R}$ defines an integral operator $\mathcal{R}: H \rightarrow H$, $f(\cdot) \mapsto \int_a^b \rho(\cdot,s)f(s) ds$ with squared Hilbert-Schmidt norm $\Vert\mathcal R\Vert_{\mathcal{S}}^2 = \int_a^b \int_a^b \rho(r,s)^2 \,\mathrm{d} s \,\mathrm{d} r < \infty$. The image space of $\mathcal R$ is $Im(\mathcal R) = \{ g \in H: g(r) = \int_a^b \rho(r,s) f(s) \,\mathrm{d} s \ \text{for some} \ f \in H\}$, and the rank of $\mathcal R$ is the dimension of its image space. An eigenvalue-eigenfunction pair $(\xi, v)$ of $\mathcal R$ satisfies $\int_a^b \rho(r,s)v(s) ds = \xi v(r)$, for all $r \in [a,b]$. If $\rho(r,s)$ is symmetric and positive semi-definite, all eigenvalues are real, and eigenfunctions associated with distinct eigenvalues are orthogonal. The adjoint operator $\mathcal R^*$ of $\mathcal R$ is the integral operator with kernel function $\rho^*(r,s) = \rho(s,r)$. See hsing2015, Sections 3 and 4, for a detailed exposition of the relevant operator theory. To clarify vector and matrix norms, we use $\|\cdot\|_2$ for the Euclidean vector norm and $\|\cdot\|_M$ for the compatible Frobenius matrix norm.

The Approximate Functional Factor Model

We consider a time series of curves $Y_1(r), \ldots, Y_T(r)$ defined on the domain $r \in [a,b]$. Our main goal is to develop a framework that is capable of capturing the low-dimensional dynamic behavior of the given functional time series. To achieve this, we employ a general factor model framework:

Here, $F_t = (F_{1,t}, \ldots, F_{K,t})'$ represents the $K\times 1$ vector of factors driving the dynamic part of $Y_t(r)$, while the error term $\epsilon_t(r)$ carries no relevant serial dependence signal, and $\mu(r)$ is the standard intercept function. The vector $\Psi(r) = (\psi_1(r), \ldots, \psi_K(r))'$ collects $K$ loading functions determining how each factor contributes to the curve series. All components of this model, including the number of factors, $K$, the vector of loading functions, $\Psi(r)$, and the factors, $F_t$, are unobserved. To properly identify their role and estimate them, an additional set of conditions is necessary, which we discuss in detail below.

We begin by separating the factor component from the error term, which is a fundamental step in our modeling approach. Central to this step is the concept of serial dependence, which allows us to isolate the low-dimensional dynamic factor component from the noise. This concept is formalized through a set of restrictions outlined in the following assumption:

By Assumption (ref)(a), the error term $\epsilon_t$ has no correlation with lagged $Y_t$, implying that it cannot be predicted from the past observations of $Y_t$. In contrast, the factors globally correlate with at least one lag of the original process, $Y_t$, as specified in Assumption (ref)(b). To be more specific, the positive definiteness of the matrix $M$, implies that $$ \lim_{T\to \infty} \sum_{\tau=1}^{q_0} \int_a^b \Big( \frac{1}{T} \sum_{t=1}^T E[F_{l,t} Y_{t-\tau}(r)] \Big)^2 \,\mathrm{d} r > 0, $$ ensuring that the global covariance between $F_{l,t}$ and $Y_{t-\tau}$ is non-zero for any $l=1, \ldots, K$ for at least one lag $\tau=1,\ldots,q_0$, and there is no cancellation of this covariance by the integration. To ensure broad applicability of our framework, we adopt the concept of global covariance, as defined by white2001, (i.e., $\lim_{T \to \infty} T^{-1} \sum_{t=\tau+1}^TE[F_{t} Y_{t-\tau}(r)]$ is used instead of $E[F_{t} Y_{t-\tau}(r)]$). This concept accommodates potential heteroskedasticity and local nonstationarities within the original process $Y_t$, which are common in economic applications. The role of $q_0$ is similar to that of portmanteau tests for autocorrelation, ensuring each factor correlates with at least one of the first $q_0$ lags of $Y_t$. lam2012 and zhang2019 argue that a small $q_0$ is sufficient in practice, and in non-seasonal setups, $q_0 = 1$ can be chosen. It is important to note that $q_0$ serves only to identify the factor model parameters and does not specify the type or features of temporal dependencies imposed on the process $Y_t$. Finally, from Assumption (ref)(a) and (b), the intercept $\mu(r)$ is identified as the mean function $E[Y_t(r)]$. This renders the first two terms in (ref), expressed as

predictable from the original process $Y_t(r)$, justifying the terminology introduced in the introduction section for the corresponding parts as predictable and unpredictable.

The set of loading functions forms a basis for the factor space $H_F := \text{span}(\psi_1, \ldots, \psi_K)$, where $Y_t$ exhibits autocovariances. The linear independence of the loading functions, as dictated by Assumption (ref)(c), implies that $H_F$ is $K$-dimensional. Essentially, this indicates the presence of $K$ distinct directions along which $Y_t$ displays temporal dependence. Consequently, identifying $H_F$ plays a pivotal role in understanding the autocovariance structure of $Y_t$, leading to the subsequent step in our model: identifying the loading functions and determining their number.

The structure of the global autocovariance of $Y_t$ depends solely on the structure of the autocovariance of the factors and the directions represented by the loadings. According to our model (ref), the global $\tau$-th order autocovariance function is expressed as:

Here, $c_\tau(r,s)$ is the kernel function of the global integral autocovariance operator $C_\tau$. Heuristically, each autocovariance $C_\tau$ for $\tau=1,...,q_0$ captures some, but not necessarily all, of the directions of interest across which the original series exhibits serial dependence. This fact follows from (ref), as these directions are linear combinations of $\psi_1, \ldots, \psi_K$. Specifically, the directions captured by $C_\tau$ are represented by its image, $Im(C_\tau)$, the space spanned by the right-singular functions of $C_\tau$. Following standard results on singular value decomposition, these functions are the orthonormal eigenfunctions of the positive semi-definite operator $C_\tau C_\tau^*$, where $C_\tau^*$ denotes the adjoint operator of $C_\tau$. Since $Im(C_\tau)$ may only capture a subset of the relevant directions (i.e., $Im(C_\tau) \subset H_F$), a combined analysis of such operators for $\tau=1,...,q_0$ becomes necessary.

We adopt the approach used in bathia2010, introducing a cumulative operator $D = \sum_{\tau=1}^{q_0} C_\tau C_\tau^*$, with its kernel function given by:

The cumulative autocovariance operator $D$ is symmetric with $Im(D) = Im(D^*) = H_F$ and $rank(D) = K$, as implied by Assumption (ref)(b) and (c). By construction, its collection of ordered orthogonal eigenfunctions, denoted here as $d_1,...,d_K$, span the entire space $H_F$, thereby capturing all directions in which $Y_t$ exhibits serial dependence. However, the loading functions $\psi_1, \ldots, \psi_K$ also span $H_F$, indicating that they are linear combinations of $d_1,...,d_K$. Hence, while $H_F$ can be properly identified through $d_1,...,d_K$, the loadings are identified up to some rotation of these eigenfunctions.

To achieve the exact identification, we resort here to the solution routinely used in conventional factor analysis, where loadings are assumed to be orthonormal and factors have a diagonal covariance matrix (see, e.g., stock2002, and bai2013). Unlike classical factor analysis, we impose a diagonal structure on the positive definite matrix $M$, as our model aims to identify factors based on the cumulative autocovariance rather than the highest variability. This leads us to the next set of restrictions outlined in the following assumption:

These conditions fix the rotation of the loading functions in the factor space $H_F$, ensuring their exact identification. To see this, by Assumption (ref)(b), the kernel function of the operator $D$ satisfies

implying that $\lambda_1, \ldots, \lambda_K$ are the descendingly ordered nonzero eigenvalues of $D$. Then the loadings $\psi_1, \ldots, \psi_K$ are identified as the eigenfunctions $d_1,...,d_K$ of operator $D$ up to a sign change and their number as the rank of this operator. It's noteworthy that $\psi_1, \ldots, \psi_K$ as eigenfunctions satisfy certain optimality properties given in the proposition below:

According to Proposition (ref), the projection of $Y_t - \mu$ onto $H_F$, as defined by

captures all components of the functional time series that correlate with its past $q_0$ lags. The projection coefficient $\langle Y_t - \mu, \psi_1\rangle$ is optimal in the sense that there are no other projection coefficients that have a higher dependency with these lags. The second projection coefficient $\langle Y_t - \mu, \psi_2 \rangle$ is optimal among all projections orthogonal to $\psi_1$, and this pattern continues with subsequent coefficients.

As the loadings $\psi_1,\ldots,\psi_K$ are deterministic, the dynamic nature of $Y_t$ can be equivalently represented by the $K\times 1$ vector $F_t^* := (F_{1,t}^*, \ldots, F_{K,t}^*)'$, where $F_{l,t}^* := \langle Y_t - \mu, \psi_l \rangle$ are the scores of projection (ref) for $l=1,\ldots,K$. Further from model (ref), we deduce:

This implies that the factors $F_{l,t}$ are partially identified as $F_{l,t}^*$, up to an unpredictable noise error term $\langle \epsilon_t, \psi_l \rangle$. In other words, under Assumptions (ref) and (ref), the factor model (ref) cannot be distinguished from its orthogonalized form: $$ Y_t(r) = \mu(r) + (\Psi(r))'F_t^* + \epsilon_t^*(r), $$ where $ \epsilon_t^*(r) := \epsilon_t(r) - \sum_{l=1}^K \langle \epsilon_t, \psi_l \rangle \psi_l(r)$.

Achieving complete factor identification may require imposing further restrictions on the model. For example, one could assume that $\epsilon_t$ solely takes values in $H_F^\perp$ (the orthogonal complement of $H_F$), implying $F_t^* = F_t$ for all $t$. Alternatively, a weaker version of this restriction, where the variance of $\epsilon_t$ in $H_F$ is asymptotically negligible, could be considered, as often done in conventional factor literature. Additionally, adopting smoothness versus roughness type of restrictions could help to separate factors from the part of the error term in $H_F$, as proposed in descary2019. We refrain from imposing additional restrictions in this study to maintain the generality of our framework and proceed with partially identified factors $F_{l,t}^*$ in the subsequent analysis. Furthermore, all results presented in the subsequent sections regarding the estimation of the model's primitives, the subspace $H_F$ and the suggested information criterion for the number of factors, $K$, do not require complete factor identification.

While Assumptions (ref) and (ref) establish the necessary restrictions to identify the role of each unobserved component in our model, additional constraints are required to specify time dependencies allowed for the process $Y_t$ as well as factors. These restrictions are required to establish the asymptotic properties of the estimators proposed in the subsequent sections.

Vector autoregressions, as described in Assumption (ref)(a), offer a general and convenient framework to model the time-dependencies of factors and allows us to estimate the loadings and the number of factors with parametric rates. Under Assumption (ref), the functional time series $Y_t$ is not required to be strictly or covariance stationary. However, in combination with Assumption (ref)(a), $Y_t$ is globally covariance stationary. Alternative weak dependence assumptions may be based on strong mixing conditions or $L^4$-$m$-approximability (see hoermann2010). Weak dependencies between factors and lagged errors, similar to those in bai2003, are allowed by Assumption (ref)(b) as long as they become sufficiently small in the limit.

Estimation

The previous section's identification results demonstrate that all unobserved components in the model can be expressed using the global autocovariances of the functional time series $Y_t$. Most of these components can be estimated using method of moments estimators by substituting population moments with sample equivalents. Section (ref) elaborates on the consistency of these estimators. However, estimating the number of factors presents a more challenging task, which is addressed in Section (ref), where we introduce an information criterion for consistently determining the number of factors. Finally, Section (ref) provides practical implementation guidelines, including an estimation and prediction algorithm.

Estimation of parameter functions

Consider the sample mean function $$ \widehat \mu(r) = \frac{1}{T} \sum_{t=1}^T Y_t(r), $$ the $\tau$-th order sample autocovariance function $$ \widehat c_\tau(r,s) := \frac{1}{T} \sum_{t=\tau+1}^T (Y_t(r) - \widehat \mu(r) (Y_{t-\tau}(s) - \widehat \mu(s)), $$ and the sample counterpart of the cumulative autocovariance function $d(r,s)$ given as

The integral operators with kernel functions $\widehat c_\tau(r,s)$ and $\widehat d(r,s)$ are denoted as $\widehat C_\tau$ and $\widehat D$, respectively. Let $\widehat \lambda_1 \geq \ldots \geq \widehat \lambda_T \geq 0$ be the eigenvalues of $\widehat D$, and let $\widehat \psi_1, \ldots, \widehat \psi_T$ be corresponding orthonormal eigenfunctions. In practice, the eigenequation and the integral in (ref) are computed by numerical integration. A ready-to-use implementation is provided in our accompanying R package.

Theorem (ref) implies that the parameter functions in model (ref) are consistently estimated with parametric rates of convergence. To estimate the factors themselves, given the estimated intercept and loading functions, we use the sample equivalents of the factors in (ref), defined as $$\widehat F_{l,t} := \langle Y_t - \widehat \mu, \widehat \psi_l \rangle, \quad l=1, \ldots, K,\ t= 1, \ldots, T.$$ Using this projection coefficient as our factor estimator is further justified by the least squares principle, as it optimizes the model fit by minimizing $$ \Big\| Y_t - \widehat \mu - \sum_{k=1}^K f_{k,t} \widehat \psi_k \Big\|^2 = \sum_{k=1}^K \big( f_{k,t}^2 - 2 f_{k,t} \langle Y_t - \widehat \mu, \widehat \psi_k \rangle \big) + \big\| Y_t - \widehat \mu \|^2, $$ where the minimum is attained when $f_{l,t} = \widehat F_{l,t}$.

Theorem (ref)(e) highlights the relevance of the sign $s_l = \text{sign}(\langle\widehat \psi_l, \psi_l \rangle)$ for our theoretical analysis, since the signs of the eigenfunctions are unidentified. Conditional on the chosen signs for the eigenfunctions of $\widehat D$, the sign-adjusted quantities $s_l \psi_l$ and $s_l F_{l,t}^* = \langle Y_t - \mu, s_l \psi_l \rangle$ serve as the population equivalents of $\widehat \psi_l$ and $\widehat F_{l,t}$. However, flipping the signs of the loadings also flips the signs of the factors, so the products $\widehat F_{l,t} \widehat \psi_l$ and $F_{l,t}^* \psi_l$ remain invariant to sign changes. Together with Theorem (ref), a direct consequence is that $$ \frac{1}{T} \sum_{t=1}^T \bigg\| \sum_{l=1}^K \widehat F_{l,t} \widehat \psi_l - F_{l,t}^* \psi_l \bigg\| = O_P(T^{-1/2}), $$ since, for any given $t=1,...,T$, the estimation error $|\widehat F_{l,t} - s_l F_{l,t}^*|$ is bounded from above by $\|\mu - \widehat \mu \| + \|Y_t - \mu \| \|\widehat \psi_l - s_l \psi_l\|$.

Estimation of the number of factors

To determine the number of factors $K$, we exploit their dynamic VAR structure, defined in Assumption (ref), and construct an information criterion for its selection. The advantage of this approach is that it not only aids in the estimation of K but also provides the opportunity to estimate the number of lags of the VAR structure, which is an essential step for empirical applications.

The temporal dynamics of the curve process $Y_t$ and its latent factors are characterized by the $K\times pK$ matrix of autoregressive coefficients denoted by $\bm{A} = [A_1, \ldots, A_p]$, as defined in Assumption (ref)(a). We employ the standard conditional least squares (LS) estimator to estimate $\bm A$. For a selected number of factors $J$ and lags $m$, the unknown $K\times 1$ vectors $F_t^*$ are replaced with the $J \times 1$ vectors of sample scores $\widehat F_t^{(J)} = (\widehat F_{1,t}, \ldots, \widehat F_{J,t})'$. The LS estimator is given by

with $\widehat\Gamma_{(J,m)}= T^{-1}\sum_{t=m+1}^T \widehat{F}_t^{(J)}(\bm{\widehat x}_{t-1}^{(J,m)})'$ and $\widehat\Sigma_{(J,m)}= T^{-1}\sum_{t=m+1}^T \bm{\widehat x}_{t-1}^{(J,m)}(\bm{\widehat x}_{t-1}^{(J,m)})'$, where the stacked vector of lagged sample scores is $\bm{\widehat x}_{t-1}^{(J,m)}= ((\widehat F_{t-1}^{(J)})', \ldots, (\widehat F_{t-m}^{(J)})')'$. Conditional on the selected number of factors and lags, the one-step ahead curve predictor is expressed as

and the corresponding mean squared error (MSE) is given by

Given that $MSE_T(J,m)$ depends on both the number of selected factors and the number of lags, it can be used to construct a consistent information criterion. To achieve this, we first must examine how $MSE_T$ behaves with respect to $J$ and $m$, considering two sources of uncertainty: one from estimating the model parameter functions and the other from estimating the vector autoregression itself. The model parameter functions that enter $MSE_T$ (i.e., $\mu$ and $\psi_l$) are invariant with respect to $J$ and $m$ and, according to Theorem (ref), are estimated with parametric rates of convergence. Hence, it remains crucial to understand how the asymptotic properties of the conditional LS estimator and therefore $MSE_T$ are impacted by a misspecified number of factors and lags.

The population coefficient matrix, $\bm{A}$, and the conditional LS estimator matrix, $\bm{\widehat A}_{(J,m)}$, are of different dimensions. To align the $K \times Kp$ matrix $\bm{A}$ with the $J \times Jm$ matrix $\bm{\widehat A}_{(J,m)} = [\widehat A_1^{(J)}, \ldots, \widehat A_m^{(J)}]$, we transform them into matrices of order $J^* \times J^*m^*$ with $J^*=\max\{J,K\}$ and $m^* = \max\{m,p\}$ by inserting zeros where their dimensions do not match. Using the completion matrix $$ \bm R_{J,K} =

$$ where $\bm 0_{J,K}$ is the $J \times K$ matrix of zeros and $\bm I_J$ is the identity matrix, we define the aligned LS estimator

Given that the loadings are only identified and correctly estimated up to a sign change, we condition our notation on the selected signs and consider the sign-adjusted matrices $\widetilde A_i = \bm S A_i \bm S$ for $i=1, \ldots, K$, with the sign transformation matrix $\bm S = \text{diag}(s_1, \ldots, s_K)$, where $s_l = \text{sign}(\langle\widehat \psi_l, \psi_l \rangle)$. The VAR process can be written as $\bm S F_t^* = \sum_{i=1}^p \widetilde A_i \bm S F_{t-i}^* + \bm S \eta_t$ since $\bm S \bm S = \bm I_K$, and the aligned sign-adjusted stacked population coefficient matrix is

Theorem (ref) shows that the consistency of the LS estimator hinges on the condition that $J\geq K$ and $m\geq p$. If either $J$ or $m$ is smaller than the actual values, the vector autoregression cannot be consistently estimated using the conditional LS estimator. This underscores the importance of the simultaneous selection of $K$ and $p$ when employing the LS estimator. For example, if the selected number of factors exceeds $K$ and the chosen lags are $m<p$, the LS estimator is biased. However, it is consistent when $m\geq p$.

The central insight from Theorem (ref) is that the MSE is asymptotically minimized when $J \geq K$ and $m \geq p$. Specifically, a model estimated with $K+j$ factors and $p+i$ lags for $i,j>0$ cannot asymptotically fit worse than a model with $K$ factors and $p$ lags. Once the threshold with the true $K$ and $p$ is met, an increase in the number of selected factors and lags does not impact the asymptotic MSE, but it may lead to parameter proliferation and a loss of efficiency. Consequently, we propose an MSE-based information criterion for estimating $K$ and $p$ of the form

where $g_T(J,m)$ serves as a penalty term for overfitting the model, and $f(\cdot)$ is a strictly increasing function. Then, the numbers of factors and lags are estimated as $$ (\widehat K, \, \widehat p ) = \operatorname*{argmin}_{\substack{J = 1, \ldots, K_{max} \\ m = 1, \ldots, p_{max}}} \mathrm{CR}_T(J,m). $$

The results of Theorem (ref) indicate that penalized MSE-based information criteria select both the correct number of factors and lags with probability 1. The crucial element for the consistent estimation of $K$ and $p$ is a penalty term that vanishes at an appropriate rate to ensure that an overparameterized model is not chosen. Commonly employed penalty terms from established information criteria in multivariate time series analysis, such as the Bayesian Information Criterion (BIC) and the Hannan-Quinn Criterion (HQC), meet the conditions outlined in Theorem (ref). Furthermore, it is also standard practice to use a logarithmic transformation to put all terms of the generic information criterion (ref) onto the same scale. These practical considerations naturally lead us to propose two types of information criteria for estimating $K$ and $p$. First, a BIC-type estimator is formulated as

where $J m$ is the number of estimated parameters in the model, and $T^{-1} \log(T)$ is the penalization rate. Second, the HQC-type estimator employs a lower penalization rate and is formulated as

Both (ref) and (ref) satisfy the conditions from Theorem (ref) and are therefore provide consistent estimators for $K$ and $p$. In practice, the minimization problem can be solved by grid search, where $K_{max}$ and $p_{max}$ must be selected.

Practical guidance

This section details the practical implementation of our estimation method and information criterion. Our primary objective is to present a procedure that can be easily executed using existing software. Building upon the theoretical foundations established in Section 3, we offer two approaches for implementing the information criterion: one based on the analytical representation of the expression in (ref), and the other based on a graphical representation. Both methods require numerical integration for computing empirical eigenfunctions and eigenvalues. Additionally, we provide an estimation and prediction algorithm that outlines how one can execute an empirical analysis of a factor models of functional time series. Our accompanying R package facilitates the execution of all proposed steps.

\paragraph{Analytical representation.} Given the selected number of factors $J$ and lags $m$, the fitted factor and error components are $$ \widehat \chi_t^{(J)}(r) = \widehat \mu(r) + \sum_{l=1}^J \widehat F_{l,t} \widehat \psi_l(r), \quad \hat \epsilon_t^{(J)}(r) = Y_t(r) - \widehat \chi_t^{(J)}(r), $$ where $Y_t(r) = \widehat \chi_t^{(J)}(r) + \hat \epsilon_t^{(J)}(r)$. The one-step ahead curve predictor can be written as $\widehat{Y}_{t|t-1}^{(J,m)}(r) = \widehat{\mu}(r) + \sum_{l=1}^{J} \widehat{F}_{l,t|t-1}\widehat{\psi}_l(r)$ with $\widehat{F}_{t|t-1}^{(J)} = (\widehat{F}_{1,t|t-1}, \ldots, \widehat{F}_{J,t|t-1})' = \bm{\widehat A}_{(J,m)} \bm{\widehat x}_{t-1}^{(J,m)}$, and the functional forecast error is

where $\hat{\eta}_{l,t}=\widehat{F}_{l,t}-\widehat{F}_{l,t|t-1}$ are the VAR residuals. From the orthonormality of the estimated loading functions, the MSE given in (ref) simplifies to

The advantage of the expression (ref) over the MSE in (ref) is that all components can be easily computed. In particular, $\widehat{\Sigma}_{\eta}^{(J,m)}$ is the least squares estimator of ${\Sigma}_{\eta}$ obtained by fitting a VAR($m$) model based on the $J$-variate time series $\widehat{F}_t^{(J)}$. The integral residual sample variance $\int_a^b \frac{1}{T} \sum_{t=1}^T (\hat \epsilon_t^{(J)}(r))^2 \,\mathrm{d} r$ equals the sum of all eigenvalues of the sample covariance operator of the residual curves $\hat \epsilon_t^{(J)}$.

\paragraph{Graphical representation.}

A careful inspection of the proof of Theorem (ref) shows that the MSE reaches its asymptotic minimum when $J \geq K$ and $m \geq p$. This result can be used to select $(K,p)$ graphically, similar to the concept of the scree plot. More precisely, one can plot $MSE_T(J,m)$ for various combinations of $J$ and $m$ and choose the minimum vertex of a rectangular surface with respect to $J$ and $m$ for which the MSE remains “flat”. For this purpose, expression (ref) can be used. Figure (ref) shows an example illustrating an MSE surface. This figure suggests that $\widehat{K}=4$ and $\widehat{p}=4$ should be selected.

The graphical approach has an advantage over the analytical expressions presented in (ref) and (ref), since it does not require the specification of the penalty term. However, it cannot be automated when it comes to multiple model selection (for instance, in Monte Carlo simulations). Furthermore, it often comes to a subjective decision of a researcher where the smallest point of the MSE rectangular “flat” area is since the estimated MSE will also fluctuate in this area in finite samples.

\paragraph{Estimation and prediction algorithm.} \

\noindentStep 1: Estimation of parameter functions. Compute the sample mean function $\widehat \mu(r)$ and the cumulative autocovariance function $\widehat d(r,s)$ from the observed curves $Y_1(r), \ldots, Y_T(r)$. Fix $K_{max}$ large enough and compute eigencomponents $\{(\widehat \lambda_l, \widehat \psi_l)\}_{l=1}^{K_{max}}$ and sample scores $\widehat{F}_{l,t} = \langle Y_t - \widehat{\mu}, \widehat{\psi}_l \rangle$, $l=1,\ldots,K_{max}$ as estimates for the factors.

\noindentStep 2: Estimation of $K$, $p$, and the factor dynamics. Fix $p_{max}$ large enough, compute $MSE_T(J,m)$ from (ref) for $J=0, \ldots, K_{max}$ and $m=0, \ldots, p_{max}$, and select $K$ and $p$ according to the BIC or HQC criterion in (ref)--(ref). Estimate a VAR($\widehat p$) model by the LS estimator given in (ref), yielding $[\widehat A_1^{(\widehat K)}, \ldots, \widehat A_{\widehat p}^{(\widehat K)}] = \widehat{\bm{A}}_{(\widehat K, \widehat p)}$.

\noindentStep 3: Fitted curves and forecasting. Compute the fitted curves for $t=1, \ldots, T$ as follows: $\widehat Y_t(r) = \widehat \mu(r) + \sum_{l=1}^{\widehat K} \widehat F_{l,t} \widehat \psi_l(r)$. Compute the VAR($\widehat p$) factor forecasts as $\widehat F_{T+1|T}^{(\widehat{K})} = \sum_{i=1}^{\widehat{p}} \widehat A_i^{(\widehat{K})} \widehat{F}_{T+1-i|T}^{(\widehat{K})}$ with $\widehat F_{T+j|T}^{(\widehat K)} = \widehat F_{T+j}^{(\widehat K)}$ for $j \leq 0$. The $h$-step curve forecast is given by $\widehat Y_{T+h|T}^{(\widehat{K},\widehat{p})}(r) = \widehat \mu(r) + \big(\widehat \Psi^{(\widehat K)}(r)\big)' \widehat F_{T+h|T}^{(\widehat{K})}$.

Simulations

We evaluate the finite sample properties of our estimators through Monte Carlo simulations. Functional time series are generated in the space spanned by the first 20 Fourier basis functions, given as $v_1(r) = 1$, $v_{2j}(r) = \sqrt 2 \sin(2 j \pi r)$, and $v_{2j+1}(r) = \sqrt 2 \cos(2 j \pi r)$ for $j=1, \ldots, 20$. We use a subset $\mathcal I$ of these Fourier functions to specify the loading functions as $(\psi_1, \ldots, \psi_K) = (v_j, \ j \in \mathcal I)$, while the remaining functions are used to form the error component. This leads to the following data generating procedure:

where $\mathcal I^c = \{1, \ldots, 20\} \setminus \mathcal I$. Here, the first term represents the factor component (ref), while the second term represents the error component. Their stochastic nature is generated through a random score vector $e_t = (e_{1,t}, \ldots, e_{20,t})' \sim \mathcal N(0,\text{diag}(1, 2^{-1}, \ldots, 20^{-1}))$, which is independent across $t=1, \ldots, T$. The elements with indices $\mathcal I^c$ of this vector are used to form the error term, as given in (ref), and the factors are defined as follows: $$ F_t = (F_{1,t}, \ldots, F_{K,t})' = A(L)\eta_t, \quad \eta_t' = (\eta_{1,t}, \ldots, \eta_{K,t}) = (e_{j,t}, \ j \in \mathcal I). $$

We consider five scenarios, M1--M5, for the dynamics of the factors, represented by $A(L)$, and subsets $\mathcal I$ of Fourier basis functions that form the factor space and indicate the number of the factors, detailed in Table (ref). It is noteworthy that the scenarios we investigate encompass empirically relevant and challenging cases where the error component can have a higher variance compared to the factors. The single-factor models M1--M3 have different loading functions and factor variances, model M4 uses a two-factor setup, and M5 is similar to the three-factor setup used in aue2015. Figure (ref) and Table (ref) confirm our theoretical findings and illustrate the consistency of $\widehat \psi_l$ as well as the BIC-type and HQ-type information criteria from equations (ref) and (ref). The estimators $\widehat K_{bic}$, $\widehat p_{bic}$, $\widehat K_{hqc}$, and $\widehat p_{hqc}$ provide a good approximation of the true parameters for reasonable sample sizes.

Empirical applications

We apply our methodology to two datasets that naturally align with a functional time series perspective: yearly mortality rate curves and monthly bond yield curves (see Figure (ref)). The first dataset consists U.S.\ log mortality rate curves of the male population, defined as $Y_t(r) = \log(D_t(r)/P_t(r))$, where $D_t(r)$ represents the number of deaths in calendar year $t$ for individuals aged $r$, and $P_t(r)$ denotes the corresponding population of age $r$. This dataset is available from the Human Mortality Database\footnote{Accessible via the R package demography, though website registration at \url{https://www.mortality.org} is required.}, comprising $T=90$ yearly curves, each observed at 111 equidistant points spanning ages 0 to 110.

The second dataset, from liu2021\footnote{Data source: \url{https://sites.google.com/view/jingcynthiawu/yield-data}.}, consists of reconstructed annualized continuously-compounded zero-coupon U.S.\ Treasury yield curves. Each curve is given by $Y_t(r)$, where $t$ denotes the calendar month and $r$ the time to maturity in months. The dataset spans November 1985 to December 2023 ($T=458$), with each curve containing 360 equidistant observations corresponding to maturities from 1 to 360 months.

Mortality curve forecasting

Mortality rates, when viewed across age groups over time, naturally take the form of functional curves, making them well-suited for functional time series analysis. U.S. mortality rate curves have been widely used in the literature, including in kokoszka2017, to illustrate functional PCA-based forecasting methods developed in hyndman2007 and aue2015. Their approach relies on a truncated Karhunen-Loève expansion: $$ Y_t(r) \approx \widehat \mu(r) + \sum_{l=1}^J g_{lt} \widehat \phi_l(r), \quad g_{lt} = \langle Y_t - \widehat \mu, \widehat \phi_l \rangle, $$ where $\widehat{\phi}_l$ are the sample functional principal components of $Y_1, \ldots, Y_T$. The forecasted curve is $\widehat{\mu}(r) + \sum_{l=1}^{J} \widehat{g}_{l,T+h|T} \widehat{\phi}_l(r)$, where $\widehat{g}_{l,T+h|T}$ are predicted scores obtained from a multivariate time series model for the vectors of sample functional principal component scores $(\widehat g_{1t}, \ldots, \widehat g_{Jt})'$. To determine the optimal truncation parameter $J$, aue2015 propose a functional final prediction error (fFPE) information criterion. While this approach is widely used and serves as a proven forecasting method, its key limitation is that the first $J$ functional principal components do not necessarily align with the most predictable components of the functional time series.

Using the full sample, we obtain $\widehat{K}_{\text{bic}} = 6$, $\widehat{K}_{\text{hqc}} = 7$, and $\widehat{p}_{\text{bic}} = \widehat{p}_{\text{hqc}} = 1$, while the fFPE criterion suggests using $J = 8$ functional principal component scores with one lag. The in-sample mean squared forecast errors (MSFE) for $\widehat{K}_{\text{hqc}} = 7$ factors is lower than for $J=8$ principal components (as suggested by the fFPE criterion), which indicates that our approximate functional factor model provides a more parsimonious dynamic representation than the PCA-based approach. Furthermore, the left panel of Figure (ref) confirms that factor-based predictions yield uniformly lower in-sample MSFE than PCA-based predictions for the same number of included components ranging from 1 to 8. The difference becomes noticeable when higher-order components are included. This is because the first three factor loadings closely resemble the first three functional principal components, whereas the fourth and fifth factor loadings deviate significantly from the PCA counterparts (see Figure (ref)).

This difference becomes more prominent when we evaluate out-of-sample forecasts using a rolling window approach with a window size of 50 observations. Model parameters are estimated within the rolling training sample, and $h$-step out-of-sample MSFEs are computed. The second plot of Figure (ref) shows that factor-based predictions more clearly outperform PCA-based predictions for 1-step forecasts in the VAR(1) model when ranging the number of included factors from 1 to 8. Furthermore, when turning to $h$-step ahead forecasts where $h=1,...,10$, we find another argument why structural representation of the dynamics is important. The right panel of Figure (ref) demonstrates that factor-based predictions using a VAR($p$) model with $K$ factors, where $K$ and $p$ are selected via the BIC or HQC criterion, increasingly outperform forecasts from a VAR($p$) model with $J$ functional principal component scores, where $J$ and $p$ are chosen using the fFPE criterion of aue2015 with increasing $h$. In conclusion, by focusing on factors that capture the underlying dynamics of the process, we obtain better forecasting performance compared to the traditional PCA-based approach.

Yield curve modeling

The dynamic Nelson-Siegel (DNS) framework introduced first by nelson1987 and further developed by diebold2006 has emerged as a workhorse model in the financial econometrics literature and has been the basis for many modifications and extensions (see svensson1995, christensen2009, lengwiler2010, and diebold2013). Functional factor models also find their role in the analysis of the term structure of bond yields. For instance, functional data models for yield curves have been explored (see hays2012, bardsley2017, sen2019, and horvath2022).

Central to the models of the Nelson-Siegel class is the assumption that the bond yield $Y_t(r)$ with time to maturity $r \in [a, b]$ at point in time $t$ follows a strict factor model framework that incorporates an additive discrete white noise component. The curve process is represented as $Y_t(r_i) = \widetilde{\chi}_t(r_i) + \widetilde{\epsilon}_{i,t}$, where $\widetilde{\chi}_t$ denotes a finite-dimensional factor component, $\widetilde{\epsilon}_{i,t}$ is white noise across $i$ and $t$, and $r_1, \ldots, r_N$ forms the discrete grid of available maturities. The shape of the loading functions and the number of factors are treated as pre-specified parameters. Specifically, the DNS model assumes the three-factor structure

with diebold2006 suggesting the decay parameter value $\xi = 0.0609$. The factors, $F_{1,t}, F_{2,t}, F_{3,t}$, are estimated through ordinary least squares at each time $t$ using the available maturities.

lengwiler2010 and nielsen2024 have already highlighted that the three-factor DNS model fails to fully capture the dynamics of yield curve data, suggesting the need for additional or more complex loading functions. Our information criteria support these findings, estimating the number of factors as $\widehat{K}_{\text{bic}} = 5$ and $\widehat{K}_{\text{hqc}} = 10$, with the number of lags estimated as $\widehat{p}_{\text{bic}} = \widehat{p}_{\text{hqc}} = 1$. In Figure (ref), we compare the first six estimated loading functions with the pre-specified loadings from the DNS model. The first three estimated loadings show similarities to the DNS loadings in terms of magnitude and curvature and have similar economic interpretations: the first factor represents long-term effects, the second short-term, and the third medium-term. The fourth, fifth and sixth factors mediate between short, long, and medium-term effects.

However, a key finding of our analysis is that the estimated number of factors and lags varies significantly depending on the time period. A rolling window analysis with a 120-month (10-year) window reveals that during periods of economic stability, 3-4 factors suffice, whereas during economic crises, substantially more factors are required. Figure (ref) shows that the BIC criterion suggests around four factors until the mid-1990s, a period that includes the 1990–1991 recession and the 1994 bond market crisis. During the relatively stable mid-2000s, the estimated number of factors drops to three. Following the 2007 housing bubble collapse, it increases to seven before fluctuating between three and five in the late 2010s and early 2020s. These findings highlight the sensitivity of factor specifications to prevailing economic conditions. Rather than fixing the number of factors a priori, we strongly recommend that practitioners adopt a data-driven approach to determine the appropriate specification.

To assess the impact of additional factors on forecasting accuracy, we employ a rolling window approach similar to diebold2006, performing sequential monthly out-of-sample yield curve forecast comparisons. At each step, we use a fixed rolling window of $w=120$ and $w=240$ for the training period. Specifically, the $h$-step-ahead forecast for time $t$ is based on data from periods $t-h-w$ to $t-h$.

In addition to comparing our forecasts with those from the dynamic DNS model, we include as a benchmark the naive random walk forecast, which is simply the observed curve of $h$ periods before. diebold2006 and caldeira2025 previously noted that DNS yield curve forecasts offer little to no improvement over the naive benchmark. Their results suggest that any potential improvements from DNS over the random walk are limited to short- and medium-term interest rates, with no significant gains observed for long-term rates. Therefore, when evaluating the MSFEs, we consider short-term interest rates (up to 12 months to maturity), medium-term interest rates (13 months to 24 months to maturity), and long-term interest rates (25 months to 360 months to maturity) separately.

Table (ref) confirms that the DNS model rarely outperforms the random walk forecast. In contrast, functional factor-based forecasts, where $K$ and $p$ are selected using the BIC or HQC criterion, consistently outperform the random walk for short-term interest rate prediction. However, as illustrated in Figure (ref), the selected values of $K$ and $p$ can become excessively large in certain periods, leading to overparameterization. Consequently, this can result in suboptimal forecasts, particularly for PCA-based predictions using the fFPE criterion, which often tends to select a large number of components and lags.

To address this issue, we also include forecasts from an L1 shrinkage-estimated VAR model, where $K$ and $p$ are selected based on the HQC criterion. The shrinkage parameter is tuned within each rolling training sample using 10-fold cross-validation, as implemented by default in the glmnet R package. Our results indicate that lasso estimators for the factor-based VAR model outperforms all other models, including the random walk, particularly for longer forecasting horizons.

Conclusion

This paper provides an in-depth study of the factor model for functional time series, including its identification, estimation, and prediction. From a practical point of view, the approximate functional factor model is an attractive modeling framework for infinitely-dimensional temporal data, as it allows analyses and predictions via a low-dimensional factor component of the data. Our results are useful for a broad range of applications in which the number of factors is unknown, and the error component potentially has strong cross-correlation and is weakly correlated with the common component. We have developed a simple-to-use novel method, yielding consistent estimates of the number of factors and their dynamics. A Monte Carlo study and an empirical illustration of yield curves show that our method provides an attractive modeling and predictive framework.

Several methodological problems await further analysis. The first is to develop the distributional and inferential theory for the estimators beyond the consistency results obtained in this paper. For instance, in the empirical illustration of yield curves, it might be interesting to provide confidence bands or test some restrictions on the loading functions. The second is to go beyond the weakly stationary assumption on the factors, for instance, by allowing some factors to have short memory while others are permitted to have long memory (persistence). Finally, the third is to develop a predictive methodology for the factors using semiparametric or nonparametric models.

Acknowledgments

We thank Jörg Breitung, Juan Carlos Escanciano, Joachim Freyberger, Tobias Hartl, Justus Henseler, Alois Kneip, Malte Knüppel, Dominik Liebl, Alexander Mayer, Daan Opschoor, and Luis Winter for their valuable comments and suggestions, and Justin Franken for his assistance with software implementations. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under project number 511905296. Moreover, the first author received financial support from the University of Bonn's Argelander Grants, while the second author was funded by Juan de la Cierva Incorporación, grant number IJC2019-041742-I. Additionally, we acknowledge the use of the CHEOPS HPC cluster for parallel computing.

Supporting Information

An accompanying R package is available at \url{https://github.com/ottosven/dffm}.

\setcounter{equation}{0} \setcounter{lemma}{0}