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.
68,432 characters · 15 sections · 59 citation commands
Macroeconomic Forecasting with Fractional Factor Models
\spacingset{1.45}
\thispagestyle{empty} \setcounter{page}{0} \paragraph{\bf Abstract.} We combine high-dimensional factor models with fractional integration methods and derive models where nonstationary, potentially cointegrated data of different persistence is modelled as a function of common fractionally integrated factors. A two-stage estimator, that combines principal components and the Kalman filter, is proposed. The forecast performance is studied for a high-dimensional US macroeconomic data set, where we find that benefits from the fractional factor models can be substantial, as they outperform univariate autoregressions, principal components, and the factor-augmented error-correction model.
\paragraph{\bf Keywords.} Fractional integration, state space model, principal components, long memory, Kalman filter
\paragraph{\bf JEL-Classification.}
C32, C38, C51, C53
At least since the seminal work of ForHalLi2000 and StoWat2002, factor models have become a popular tool for forecasting macroeconomic dynamics, as they handle covariation in the cross-section efficiently by condensing it to a typically small number of common latent factors. Regardless their applicability to large data sets, the major drawback of standard factor models is an inefficient use of longitudinal information: in contrast to e.g.\ VARMA models, the vast majority of factor models requires stationarity. Consequently, key features of macroeconomic data, such as nonstationary trends and cointegration, are not captured adequately by standard factor models but rather differenced away. Over-differencing of latent processes poses an additional risk, since model selection criteria and model specification tests for the number of factors are likely to miss these components, as eigenvalues corresponding to over-differenced series converge to zero.
A more flexible setup is suggested by a young strand of the factor model literature that adds unit roots to the model PenPon2006, Eic2009, ChaMilPa2009, BanMarMa2014, BanMarMa2016, BarLipLu2016. But these models come with the drawback of requiring a priori assumptions about the degree of persistence, and typically the series under study are assumed to be $I(1)$. This makes an endogenous treatment of the (unknown) long-run dynamic characteristics of observable time series impossible. Statistical inference about the degree of persistence of an observable variable is then limited to prior unit root testing, ignoring the non-standard behavior of many economic series that are fractionally integrated. Misspecifying the integration orders of the observable variables may bias the factor estimates, can yield wrong inference about the number of common factors, and is likely to deteriorate the forecast performance.
To address these problems, semiparametric methods that are robust to fractional integration have been proposed by LucVer2015 for a single fractionally integrated factor and by Erg2017 for pervasive fractionally integrated nuisance. Allowing for a wide range of persistence and an endogenous treatment of integration orders, HarWei2018 derive a parametric fractionally integrated factor model and apply it to realized covariance matrices.
In macroeconomics, fractionally integrated factor models have not played a role so far, although there is comprehensive evidence for long memory and fractional cointegration in the data HasWol1995, Bai1996, GilRob1997, TscWebWe2013.
Tackling this issue, this paper aims to provide insights on whether fractional integration techniques have merit at least for a relevant fraction of the numerous and heterogeneous macroeconomic variables typically under study. By elaborating fractionally integrated factor models, we construct setups where cross-sectional covariation in the data in levels is driven by fractionally integrated latent factors that may impose cointegration relations. In detail, we propose three different factor models that generalize the aforementioned factor models to fractionally integrated processes. The first model introduces ARFIMA processes in the nonstationary factor model setup of BarLipLu2016, while the second model distinguishes between purely fractionally integrated factors that impose cointegration relations and $I(0)$ factors that model common short-run behavior of the data. Finally, our third model generalizes the pre-differencing of the data for standard $I(0)$ factor models by taking fractional differences.
As standard factor models they are applicable to high-dimensional data, but bear several advantages: the fractional factor models allow for a joint modelling of data of different persistence, do not require prior assumptions about the degree of persistence of the data but treat the integration orders endogenously, they capture cointegration via the common fractionally integrated factors and are more robust to over-differencing.
For the estimation of the latent factors we introduce a two-stage estimator, where initial factor estimates are obtained via principal components, until the model is cast in state space form such that the Kalman filter and smoother is applicable. For the latter to be computationally feasible, we use ARMA approximations for fractionally integrated processes as suggested in HarWei2018a. Estimation of the unknown model parameters and the latent factors is then carried out jointly via an expectation maximization algorithm.
In a pseudo out-of-sample forecast experiment for a high-dimensional US macroeconomic data set of McCNg2015, we study the forecast accuracy of the fractional factor models. We provide a guided choice among the different models by considering the forecast performance for $112$ macroeconomic variables. Finally, we find comprehensive evidence that adequately combining fractional integration techniques and factor models can improve forecasts substantially compared to standard factor models and other benchmarks.
The remaining paper is organized as follows. Section (ref) details the construction of fractional factor models. The two-stage estimator for the factors and model parameters is discussed in section (ref). Section (ref) compares the forecast performance of the fractional factor models to different benchmarks in a pseudo out-of-sample forecast experiment, until section (ref) concludes.
To begin with, consider the general form of a high-dimensional factor model for possibly fractionally integrated data
where $\bm y_t=\left(y_{1,t},...,y_{N,t}\right)' $ is an $N$-dimensional observable time series with entries $y_{i,t} \sim I(d^*_i)$ that are integrated of order $d_i^*$, $d^*_i \in \mathbb{R}_{\geq 0}$. An integration order $d_i^*$ implies that the fractional difference of a series is $I(0)$, i.e.\ $\Delta^{d_i^*}y_{i,t}\sim I(0)$, $i=1,...,N$. The vector $\bm \chi_t$ is $r$-dimensional and accounts for common short- and long-run dynamics among the $\bm y_t$, and $\bm u_t= (u_{1, t},...,u_{N, t})'$ holds the $N$ idiosyncratic errors and has a diagonal variance.
The fractional difference operator $\Delta^d$ is defined as
and a $+$-subscript amounts to a truncation of an operator at $t \leq 0$, e.g.\ for an arbitrary stochastic process $z_t$, $\Delta^d_+ z_t = \sum_{j=0}^{t-1}\pi_j(d) L^j z_t$ Joh2008. For $d \in \mathbb{N}_0$ fractionally integrated processes nest the standard integer integrated specifications (e.g.\ $I(0)$, $I(1)$, and $I(2)$ processes), whereas $d \in \mathbb{R}_{\geq 0}$ adds flexibility to the weighting of past shocks. Throughout the paper, we adopt the type II definition of fractional integration MarRob1999 that assumes zero starting values for all fractional processes, and, as a consequence, allows for a seamless treatment of the asymptotically stationary ($d < 1/2$) and the nonstationary ($d \geq 1/2$) case. Due to the type II definition the inverse fractional difference $\Delta_+^{-d}z_t$ exists.
Standard factor models as those considered in ForHalLi2000, BaiNg2002, and StoWat2002, are special cases of (ref). They extract $r$ common factors of a data set in first and second differences, implying that the series in $\bm y_t$ are $I(1)$ and $I(2)$. The common factors in $f(\bm \chi_t)$ then correspond to the common trends in the Granger representation theorem for cointegrated data BarLipLu2016.
To give an intuition on how fractional integration affects the long-run properties of a time series we note the following. For positive $d$ the autocovariance function of an $I(d)$ process decays at a hyperbolic rate, implying that a shock has a persistent impact on the $I(d)$ process, and a greater $d$ implies a more persistent impact of a shock. While an $I(d=1)$ process is an unweighted sum of past shocks, an $I(d)$ process in general can be interpreted as a weighted sum of past shocks, where the weights depend on $d$ via (ref). Furthermore, if a linear combination of a vector $I(d)$ process exists that is integrated of order $b < d$, then the series are cointegrated. Cointegration implies common (fractionally) integrated trends, which our models capture via $f(\bm \chi_t)$. For a discussion of cointegration relations in a fractionally integrated factor model setup we refer to HarWei2018.
We introduce three different fractionally integrated factor models in the next sections that are nested in (ref) and differ in the functional relation between $\bm \chi_t$ and $\bm y_t$. Section (ref) generalizes nonstationary factor models BarLipLu2016 to fractionally integrated processes. In section (ref) we distinguish between fractionally integrated factors, that account for long-run co-movements in $\bm y_t$, and $I(0)$ factors, that allow common short-run dynamics. Finally, section (ref) generalizes the pre-differencing of standard factor models to fractional differencing.
Consider a simple multivariate unobserved components model
where $f(\bm \chi_t)=\bm \varLambda \bm f_t$ in (ref), $\bm f_t = (f_{1,t}, ..., f_{r,t})'$ holds the $r$ common factors, $\bm \varLambda$ is a $N \times r$ matrix of factor loadings that is assumed to have full column rank, and the errors $\bm u_t$ account for idiosyncratic dynamics. The latent factors are assumed to follow $r$ fractionally integrated autoregressive processes
where ${B}_j(L)= 1-\sum_{k =1}^{p}{B}_{j,k}L^{k}$ is a stable lag polynomial. For the pervasive shocks that drive $\bm f_t$ we assume $(\zeta_{1,t},...,\zeta_{r,t})'= \bm{\zeta}_{t}\sim \mathrm{NID}(\bm{0}, \bm Q)$, where $\bm Q$ is diagonal. A matrix formulation of (ref) follows directly by defining $\bm d = (d_1,...,d_r)'$, the lag polynomials $\bm D(\bm d) = \mathrm{diag}(\Delta^{d_1}_+,...,\Delta^{d_r}_+)$ and $\bm B(L) = \mathrm{diag}(B_1(L),...,B_r(L))$, such that $\bm B(L)\bm D(\bm d)\bm f_t = \bm \zeta_t$.
The errors ${u}_{i,t}$ are assumed to be mutually independent and are allowed to be autocorrelated
where $\rho_{i}(L)=\sum_{k=0}^{p_i}\rho_{i,k}L^k$ is a stationary autoregressive lag polynomial.
As a consequence, the model may explain various degrees of common persistence that characterize the data by common components with long memory. For $d_1 =...=d_r =0$, the model nests the approximate dynamic factor model of StoWat2002, while $d_j \in \{0, 1\}$, $j=1,...,r$, yields a nonstationary dynamic factor model with $I(1)$ factors as considered in e.g.\ BarLipLu2016. Therefore, the model can be interpreted as a fractional generalization that neither requires prior differencing of the data, nor a prior assumptions about the integration orders.
A more parsimonious factor model is proposed by HarWei2018. Their model distinguishes between $r_1$ purely fractional factors $\bm f_t^{(1)}=(f_{1, t}^{(1)},...,f_{r_1, t}^{(1)})'$, that establish cointegration relations among the $\bm y_t$, and $r_2$ stationary autoregressive components $\bm f_t^{(2)}=(f_{1, t}^{(2)},...,f_{r_2, t}^{(2)})'$, that account for common short-run behavior. We consider a slightly more general modification that allows for autocorrelated idiosyncratic errors. The general framework for the dynamic orthogonal fractional components model is then given by
for all $t = 1,...,T$ and $r=r_1 + r_2 \leq N$. The $N$ idiosyncratic shocks $\bm{\xi}_t = (\xi_{1,t},...,\xi_{N,t})'$ are assumed to follow independent Gaussian white noise processes $\bm{\xi}_t \sim \mathrm{NID}(\bm{0}, \bm{H})$. For the pervasive shocks $\bm \zeta_t^{(1)} = (\zeta_{1, t}^{(1)}, ..., \zeta_{r_1, t}^{(1)})'$, $\bm \zeta_t^{(2)} = (\zeta_{1, t}^{(2)}, ..., \zeta_{r_2, t}^{(2)})'$ we assume $\mathrm{vec}(\bm{\zeta}^{(1)}_{t}, \bm{\zeta}^{(2)}_t) \sim \mathrm{NID}(0, \bm Q)$ where $\bm Q$ is diagonal. In addition, we assume that the errors $\bm u_t$ are independent of the common components $\bm f_t$.
Define the polynomials $\bm B^{(2)}(L) = \operatorname{diag} (B_1^{(2)}(L),...,B_{r_2}^{(2)}(L))$, $\bm D^{(1)}(\bm d) = \operatorname{diag} (\Delta_+^{d_1},...,\Delta_+^{d_{r_1}})$. Then, the model can be shown to be nested in the setup of section (ref) for $\bm f_t=\mathrm{vec}(\bm f_t^{(1)}, \bm f_t^{(2)})$, $\bm B(L) = \mathrm{diag}(\bm{I}, \bm B^{(2)}(L))$, and $\bm D(\bm d)=\mathrm{diag}(\bm D^{(1)}(\bm d), \bm{I})$. In terms of (ref) the model specifies $f(\bm \chi_t)=\bm \varLambda^{(1)} \bm f_t^{(1)} + \bm \varLambda^{(2)} \bm f_t^{(2)}$.
Note that the NID assumption on $\bm \zeta_t$ together with $\bm Q$ diagonal yields $r$ orthogonal factors $\bm f_t$. This common feature of many unobserved components models, which also applies to the models in sections (ref) and (ref), reduces estimation uncertainty of the loadings and makes the framework very attractive for forecasting. Since $\bm u_t$, $\bm \zeta_t$ are assumed to be independent, any correlation among the variables in $\bm y_{t}$ stems from the common long- and short-run components $\bm f^{(1)}_{t}$ and $\bm f^{(2)}_{t}$.
A third model that completes our toolbox of fractionally integrated factor models takes fractional differences of the observable variables to arrive at a short memory model, where all components are at most $I(0)$. Hence, we contrast our two models from sections (ref) and (ref) with an additional approach that excludes fractional integration from the factors. For this purpose we define
As before, $\bm y_t=(y_{1,t},...,y_{N,t})'$ are the observable variables, $\bm{\Lambda} = [\bm{\Lambda}_{1}', ..., \bm{\Lambda}_{N}']'$ holds the factor loadings, and $\bm f_{t}=(f_{1,t},...,f_{r, t})'$ contains the $r$ latent factors. In the notation of (ref) this implies $f(\bm \chi_t) = \bm D(-\bm d^*)\bm \varLambda \bm f_t$ and $\bm u_t = \bm D(-\bm d^*)\bm \xi_t$ with $\bm \xi_t = (\xi_{1, t},..., \xi_{N, t})'$, and $\bm d^* = (d_1^*,...,d_N^*)'$.
By defining $\bm B(L) = \mathrm{diag}({B}_{1}(L), ..., {B}_{r}(L))$ as in sections (ref) and (ref) the factors $\bm f_t$ can be written as a diagonal VAR process $\bm B(L)\bm f_t = \bm \zeta_t$, where $\bm \zeta_t = (\zeta_{1, t}, ..., \zeta_{r, t})'$. The idiosyncratic and pervasive shocks are assumed to be orthogonal and to follow independent Gaussian white noise processes $\bm{\xi}_t \sim \mathrm{NID}(\bm{0}, \bm H)$ and $\bm{\zeta}_t \sim \mathrm{NID}(\bm{0}, \bm Q)$.
By taking fractional differences prior to estimating a factor model, our approach generalizes the pre-differencing of standard factor models to the fractional domain. In fractional differences, our model is an approximate dynamic factor model and, therefore, it nests the model of StoWat2002 for $d_1^*,...,d_N^* \in \mathbb{N}_0$.
Taking fractional differences of order $d_i^{*}$ ensures for each $\Delta_+^{d_i^*}y_{i,t}$ that the common and idiosyncratic components are at most $I(0)$. Note that fractional differences are less sensitive to over-differencing compared to integer differences, since the method ensures that the fractional difference of the most persistent factor that loads on $y_{i,t}$ is $I(0)$ $\forall i =1,...,N$.
In this section we discuss both, the estimation of the latent factors $f(\bm \chi_t)$ in (ref) for the three different factor models proposed in sections (ref) to (ref), and the estimation of the unknown model parameters. The expectation-maximization (EM) algorithm is a natural choice for the estimation of parametric factor models JunKoo2015 and has been derived for fractionally integrated factor models in HarWei2018a. In the E-step, the latent factors are estimated given a set of parameters via the Kalman filter. The M-step then updates the parameter vector by maximizing the likelihood function given the factor estimates from the E-step. Therefore, the EM-algorithm allows for a joint estimation of factors and model parameters.
Since the EM-algorithm is a parametric estimator, it requires starting values for the unknown model parameters in sections (ref) to (ref). We tackle this problem by proposing a two-stage estimator. The first stage is described in section (ref). We estimate the latent factors via the nonparametric method of principal components (PC) and propose estimators for the unknown model parameters. We include a consistency proof for the PC estimator for fractionally integrated factors with integration orders in $\mathbb{R}_{\geq 0}$, since consistency of the PC estimator has so far only been shown in more restrictive settings.
The second stage is considered in section (ref). We derive an approximate state space formulation for each of the factor models in sections (ref) to (ref), so that the Kalman filter can be applied to estimate the latent factors. Finally, we discuss the joint estimation of the model parameters and the latent factors via the EM algorithm.
Sufficient conditions for a consistent estimation of $f(\bm \chi_t)$ in (ref) via PC were derived in BaiNg2002 for stationary processes, in Bai2004 for $I(1)$ common components and in BaiNg2004 for $\bm y_t \sim I(1)$, where nonstationarity may also stem from the idiosyncratic components. For a single fractionally integrated factor and fractional integration orders in $[0, 1]$ LucVer2015 have shown that the methods of BaiNg2004 are also applicable. We generalize their results to non-negative integration orders and multiple fractionally integrated factors by showing consistency of the PC estimator for $f(\bm \chi_t)$ in (ref).
Since PC are estimated via an eigendecomposition of $\mathrm{Var}(\bm y_t)$, the applicability of the PC estimator depends crucially on the stability of the variance. For $\mathrm{max}(d_i^*)<0.5$ all $\bm y_t$ are asymptotically stationary, and consequently the variance of $\bm y_t$ converges as $t \to \infty$. Therefore, the PC estimator satisfies the assumptions of BaiNg2002, where assumption A postulates boundedness of $\mathrm{plim}_{T \to \infty} T^{-1} \sum_{t=1}^T \bm f_t \bm f_t' = \bm \varSigma_f < \infty$.
For $\mathrm{max}(d_i^*)\geq 0.5$ assumption A of BaiNg2002 is violated. Nonetheless, under a suitable scaling the PC estimator is still consistent. We report an updated set of assumptions for consistency of the PC estimator for nonstationary data in appendix (ref). Following BaiNg2002 and Bai2004, for $d_1 = ... = d_r$ we show that there exists a matrix $\bm H$ such that the factors $ \bm f_t$ are estimated consistently up to a rotation by PC
Expressions for $\hat{\bm f}_t$ and $\bm H$, together with a detailed proof, are given in appendix (ref).
Whenever there is at least one $d_j \neq d_k$, $j, k = 1,...,r$, direct estimation of all common fractionally integrated factors is not feasible, since, depending on the scaling of the PC, either the contribution of the least persistent factors to the covariance of $\bm y_t$ converges to zero, or the contribution of the most persistent factors diverges. In this case, one needs to separate $\bm y_t$ into blocks of equal persistence. Starting with the most persistent block, latent factors are estimated via PC and projected out. The adjusted variables are then added to the next block of $\bm y_t$, and the procedure repeats, until a stationary set of variables is obtained.
Having established consistency of PC for the estimation of the latent fractionally integrated factors, we turn to the estimation of the dynamic parameters for the common factors. Since the dynamic properties differ among the three frameworks discussed in section (ref), we consider them separately in the following.
\paragraph{ARFI factors} The common components of the model in section (ref) are assumed to follow $r$ independent autoregressive fractionally integrated processes. Therefore, we rotate the PC estimates via the method of MatTsa2011 to obtain dynamic orthogonal components. The parameters in (ref) are estimated by maximizing the likelihood function for a multivariate fractionally integrated process Nie2004 that is given by
with $\bm B = \mathrm{vec}(\bm B_{1}, ..., \bm B_{p})= \mathrm{vec}(\bm B(L))$, and $\bm B(L)$, $\bm D(\bm d)$ as defined in section (ref). Plugging in the first-order condition $\bm{\bm Q} = \bm{\bm Q}(\bm d, \bm B) = T^{-1}\sum_{t=1}^{T}\left(\bm B(L)\bm D(\bm d)\bm f_{t}\right)\left(\bm B(L)\bm D(\bm d)\bm f_{t}\right)'$, and dropping the constant terms gives $ l^{*}(\bm d, \bm B)=-\frac{T}{2}\mathrm{log}\vert\bm Q(\bm d, \bm B)\vert, $ which we maximize to get estimates for the unknown parameters $d_{1}, ..., d_{r}$, $\bm B_{1}, ..., \bm B_{p}$. For some data sets the assumption of orthogonal factors may be violated. Then, the diagonal assumption on $\bm B(L)$ can be dropped, which does not affect the identification of the fractional factor VAR but increases the number of unknown parameters in (ref). Factor loadings $\bm \varLambda$ in (ref) are estimated via ordinary least squares (OLS).
\paragraph{FI and AR factors} To derive an estimator for the dynamic parameters of the model in section (ref), we first need to distinguish between the space spanned by the purely fractional factors and the stationary autoregressive components. We identify the two factor subspaces of ${\bm f}_t^{(1)}$ and ${\bm f}_t^{(2)}$ up to a rotation by estimating the fractional cointegration subspace and its orthogonal complement via the semiparametric method of CheHur2006, who use eigenvectors of an averaged periodogram matrix of the first $m$ Fourier frequencies to estimate the fractional cointegration subspace. Finally, orthogonal series within the fractional and non-fractional factors are obtained by applying the decorrelation method of MatTsa2011. The resulting fractional and non-fractional factor estimates are denoted as $\hat{\bm f}_{t}^{(1)}$ and $\hat{\bm f}_{t}^{(2)}$ respectively.
Given the factor estimates $\hat{\bm f}_t^{(1)}$ and $\hat{\bm f}_t^{(2)}$ together with the observable variables $\bm y_t$, we estimate the factor loadings $\bm{\Lambda}$ in (ref) and the AR coefficients of (ref) via OLS. Estimates for the fractional integration orders of the common components in (ref) are obtained by maximizing the likelihood of the $r_1$ ARFIMA(0, $d_j$, 0) processes, $j=1,...,r_1$. \paragraph{AR factors} Due to the stationary representation of the model in section (ref) the PC estimator of BaiNg2002 is directly applicable. The factors are again decorrelated by means of dynamic orthogonal components of MatTsa2011. For a discussion of the consequences when a diagonal representation of the common factors is not feasible we refer to the ARFI case. The dynamic coefficients for the $r$ common factors in (ref) together with their factor loadings in (ref) are estimated via OLS. \paragraph{AR errors} An estimate for the idiosyncratic errors is obtained via $\hat{\bm u}_t = \bm y_t - \hat{\bm \varLambda}\hat{\bm f}_t$. Since the errors are assumed to follow $N$ independent autoregressive processes, the AR parameters are estimated via OLS.
The second stage of our estimator combines factor estimation for a given set of parameters via the Kalman filter and smoother together with parameter optimization via maximum likelihood (ML) in an EM algorithm. For the Kalman filter to be applicable, the different components of our fractional factor models are cast in state space form. Note that for a given sample size $T$ a finite state space representation of a type II fractionally integrated process exists but requires a state vector of dimension $T-1$, as (ref) shows. Since the Kalman filter sequentially inverts the $(T-1) \times (T-1)$ autocovariance matrix for each factor, a full representation of a fractionally integrated process can be very costly from a computational perspective, in particular for long time series. Therefore, section (ref) discusses finite approximations that resemble the dynamic properties of fractionally integrated processes well and are computationally feasible. Section (ref) derives the state space representation and section (ref) considers parameter estimation.
The literature has considered a variety of approximations for long memory processes: Pal2007 suggests truncated AR approximations, whereas ChaPal1998 study truncated MA approximations. In a simulation study, HarWei2018a find that small ARMA($v, w$) models with $v, w \in \{3, 4\}$ outperform pure AR and MA approximations even if a high number of lags enter the latter models. In addition, the ML estimator for the integration order is found to be more precise when an ARMA approximation is used. As their simulation studies show, the ML estimates for an approximate representation of a fractionally integrated process converge to the ML estimates of the exact state space representation as $T \to \infty$. For the latter, consistency is proven in HarTscWeb2019.
Following the suggestions of HarWei2018a, an ARMA($4, 4$) process is used to approximate the purely fractional factors of section (ref). For ARFIMA processes, whose dynamic properties stem not only from the fractional differencing operator, the approximation quality of ARMA processes is not clear. Therefore, we use pure AR($5$) processes to resemble the properties of the fractional differencing operator in the ARFI-case of section (ref). For an arbitrary integration paramter $b$ the approximations are given by
where $m_{k}(b)$ are the MA coefficients, $k = 1, ..., w$, and $a_{l}(b)$ are the AR parameters, $l = 1,...,v$, $(v=4, w=4)$ for purely fractional factors as in (ref), and $(v=5, w=0)$ for ARFI-factors as in (ref).
The ARMA parameters are chosen beforehand for a given sample size $T$ and fractional integration order $b$ by minimizing the distance between the generic process $x_t = \Delta_+^{-b}z_t=\sum_{j=0}^{t-1}\pi_j(-b)z_{t-j}$ and its approximation $\tilde{x}_t = [m(L, b)a(L, b)^{-1}]_+z_t=\sum_{j=0}^{t-1}\tilde{\psi}_j(-b) z_{t-j}$, $z_t \sim NID(0, 1)$, over all $t=1,...,T$, where $\tilde{\psi}_j(-b)$ is the $j$-th coefficient of the ARMA Wold representation, and $\pi_j(-b)$ is its counterpart from (ref). We use the mean squared error over $t=1,...,T$ as the distance measure $ MSE_{T}^{b} = \frac{1}{T} \sum_{t=1}^{T}\sum_{j=0}^{t-1}\left(\tilde{{\psi}}_j(-b) - {\pi}_j(-b)\right)^2. $
For a given sample size $T$ and integration order $b$, we collect the ARMA coefficients in a $(v+w)$-vector $\bm \varphi_T(b)=(a_1(b),...,a_v(b), m_1(b),...,m_w(b))'$. The ARMA coefficient estimates are then defined via $\hat{\bm \varphi}_T(b) = \arg \min \bm \varphi MSE_T^b$. Following HarWei2018a, for a given $T$ optimization is carried out for each value on a grid for $b$. The ARMA coefficients are smoothed using cubic regression splines, such that a continuous, differentiable function $\bm \varphi_T(b)$ in $b$ is obtained. The technical details and several simulation studies are contained in HarWei2018a.
With a smooth function $\bm \varphi_T(b)$ in $b$ at hand, parameter optimization for the models in sections (ref) and (ref) can be conducted over the low-dimensional vector of fractional integration orders $\bm d$, which keeps the dimension of the parameter vector within the optimization procedure manageable and independent of the length of the ARMA approximations.
With these approximations at hand, we can turn to the state space representation of our fractional factor models. A general representation of a state space model is given by
where $\bm a_{t|T}=\mathrm{E}(\bm{\alpha}_{t}|\tilde{\bm y}_1,...,\tilde{\bm y}_{T}, \bm{\theta})$, $\bm P_{t|T}=\mathrm{Var}(\bm{\alpha}_t|\tilde{\bm y}_1,...,\tilde{\bm y}_T, \bm{\theta})$. The covariance matrices of the disturbances $\bm Q = \mathrm{Var}(\bm{\zeta}_{t})$, $\bm H = \mathrm{Var}(\bm{\xi}_{t})$ are diagonal $\forall t = 1,..., T$. Without loss of generality, we set $\bm Q = \bm I$ for all fractional factor models in state space form to distinguish between the factor loadings $\bm \varLambda$ and the variance of the factor innovations $\bm Q$. The matrices $\bm Z$, $\bm T$, $\bm R$, and the states $\boldsymbol{\alpha}_t$ differ for the three fractional factor models and are derived separately in the following. \paragraph{ARFI factors} By approximating the fractional difference operator of our first model that is given in section (ref), equation (ref) becomes
Let $\bm A(L, \bm d) = \bm I - \sum_{j=1}^v \bm A_j(\bm d)L^j$, $\bm A_j(\bm d) = \operatorname{diag}(a_{j}(d_1),...,a_{j}(d_r))$, $j=1,...,r$. Then $\bm B(L)\bm A(L, \bm d)=\sum_{k=0}^{p+v}\sum_{l=0}^{k}\bm B_{l}\bm A_{k-l}(\bm d)L^k$ where $\bm A_{0}(\bm d)=\bm B_{0}=-\bm I$, $\bm A_l(\bm d)=\boldsymbol{0} \ \forall l > v$, and $\bm B_l = \boldsymbol{0} \ \forall l > p$.
JunKoo2015 suggest to eliminate autocorrelation in the idiosyncratic errors $\bm u_t$ via the observations equation instead of accounting for them via the state equation. We follow their suggestion and manipulate the observable variables
For the state space representation we collect the adjusted observable variables in $\tilde{\bm y_t}=(\tilde{y}_{1,t},..., \tilde{y}_{N,t})'$ and define $\bm{\Psi}_j = \mathrm{diag}(\rho_{1,j}, ..., \rho_{N, j})$ such that $\tilde{\bm y}_t=\bm y_t-\sum_{j=1}^{\mathrm{max}(p_i)}\bm{\Psi}_j\bm y_{t-j}.$
A state space representation of (ref), (ref), and (ref) follows directly by defining the system matrices $\bm T$, $\bm Z$, $\bm R$, together with the state vector $\bm \alpha_t$ in (ref) as follows. $\bm T$ depends on $\bm d$ and $\bm B(L)$, whereas $\bm Z$ depends on $\bm \varLambda$ and ${\rho_i}(L)$, $i=1,...,N$,
$\boldsymbol{\alpha}_t=(\bm f_{t}',...,\bm f_{t-u+1}')'$ holds the states, $\bm R=[\bm I, \bm{0}]'$ is a selection matrix and $u$ is defined as $\mathrm{max} (p+v, \mathrm{max} (p_i)+1)$. To distinguish between the $r$ factors, we restrict the first $r$ rows of $\bm{\Lambda}$ to form a lower triangular matrix. \paragraph{FI and AR factors} Starting with the ARMA approximations of the purely fractional factors in (ref), an approximate representation of the latent fractionally integated factors is given by $\bm f_t^{(1)}\overset{a}{=}\bm M(L, \bm d)\bm A(L, \bm d)^{-1}_+\boldsymbol{\zeta}_t^{(1)}$, where the matrix AR and MA polynomials are $\bm M(L, \bm d)=\bm I + \bm M_1(\bm d)L + ... + \bm M_w(\bm d)L^w$, $\bm M_j(\bm d)=\operatorname{diag} (m_j(d_1),...,m_j(d_{r_1}))$, $\bm A(L, \bm d)=\bm I - \bm A_1(\bm d)L - ... - \bm A_v(\bm d)L^v$, $\bm A_j(\bm d)=\operatorname{diag} (a_j(d_1),...,a_j(d_{r_1}))$, and $\bm M_j(\bm d) = \bm{0}$ $\forall j > w$, $\bm A_j(\bm d)=\bm{0}$ $\forall j > v$.
Regarding $\bm u_t$, we again eliminate autocorrelation from the idiosyncratic errors by manipulating $\bm y_t$ as in (ref), i.e.\ $\tilde{\bm y}_t = \bm y_t - \sum_{j=1}^{\mathrm{max}(p_i)} \bm \varPsi_j \bm y_{t-j}= \bm \varPsi(L)\bm y_t$. For the latent fractionally integrated factors, this implies $\bm \varPsi(L)\bm \varLambda^{(1)}\bm f_t^{(1)}\overset{a}{=}\bm \varPsi(L)\bm \varLambda^{(1)}\bm M(L, \bm d)\bm A(L, \bm d)^{-1}_+\boldsymbol{\zeta}_t^{(1)}$.
The state space form (ref) of the model is then obtained by imposing a block diagonal structure on $\bm T=\mathrm{diag}(\bm T^{(1)}, \bm T^{(2)})$, where the first block $\bm T^{(1)}$ solely depends on $\bm d$, whereas the second block $\bm T^{(2)}$ depends on $\bm B(L)$
where $u_1=\mathrm{max} (v, w+\mathrm{max} (p_i)+1)$, $u_2=\mathrm{max}(p, \mathrm{max} (p_i)+1)$, and $\mathrm{max}(p_i)$ is the maximum lag order of the idiosyncratic errors $\bm u_t$ in (ref). $\bm T^{(1)}$ accounts for the dynamic properties of the fractionally integrated factors, whereas $\bm T^{(2)}$ models the stationary variation of the $\bm f_t^{(2)}$.
The two blocks for $\bm Z=
$ depend on $\bm{\Lambda}^{(1)}$, $\bm{\Lambda}^{(2)}$, $\bm d$, and $\boldsymbol{\rho}(L)$, and are given by
whereas the state vector is given by
Note that $\bm{\Psi}_j = \bm{0} \ \forall j > \mathrm{max}(p_i)$ and $\bm B_j^{(2)}=\bm{0}\ \forall j>p$. Finally, the selection matrices are given by $\bm R = \mathrm{diag}(\bm R^{(1)}, \bm R^{(2)})$, $\bm R^{(1)} = [\bm I, \bm{0}]'$, $\bm R^{(2)} = [\bm I, \bm{0}]'$, whereas the disturbances in the state equation are $\bm{\zeta}_{t+1} = (\bm{\zeta}_{t+1}^{(1)'}, \bm{\zeta}_{t+1}^{(2)'})'$, $\bm{\zeta}_{t} \sim \mathrm{NID}(\bm{0}, \bm Q)$ for all $t = 1,...,T$. Note that for the fractionally integrated factors the observations equation yields $\bm Z^{(1)}\bm \alpha_t^{(1)} = \bm \varPsi(L)\bm \varLambda^{(1)} \bm M(L, \bm d) \bm \alpha_t^{(1)} =\bm \varPsi(L)\bm \varLambda^{(1)} \bm M(L, \bm d)\bm A(L, \bm d)_+ \bm \zeta_{t}^{(1)} \overset{a}{=} \bm \varPsi(L)\bm \varLambda^{(1)}\bm f_t^{(1)}$, whereas for the stationary AR factors it gives $\bm Z^{(2)}\bm \alpha_t^{(2)} = \bm \varPsi(L)\bm \varLambda^{(2)}\bm f_t^{(2)}$. \\ The $r_1$ independent fractional factors are identified by imposing a block triangular structure on $\bm{\Lambda}^{(1)}$ while sorting the observations $\bm y_t$ with respect to their order of fractional integration in ascending order. As a consequence, the first block of variables in $\bm y_t$ is driven by the least persistent factor $f_{1t}$, the second block of variables depends on $f_{1t}$ and $f_{2t}$ whereas the $r_1$-th block with the highest order of fractional integration is allowed to be influenced by all fractional factors. In addition, the first $r_2$ rows of $\bm{\Lambda}^{(2)}$ form a lower triangular matrix to identify the $I(0)$ factors $\bm f_t^{(2)}$.
\paragraph{AR factors} Since the factors of our third model (ref) are stationary autoregressive processes, a state space representation as in (ref) follows immediately by defining $\tilde{\bm y}_t = (\Delta_+^{d_1^*}y_{1, t},...,\Delta_+^{d_N^*}y_{N, t})'$. The factors enter the state vector directly, whereas their dynamic coefficients in (ref) are contained in $\bm T$. Furthermore, the factor loadings are modelled via $\bm Z$, and $\bm R$ is again a selection matrix
For identification of the factors, we restrict the first $r$ rows of $\bm{\Lambda}$ to be lower triangular.
We collect the unknown parameters in $\bm d$, ${\bm \varLambda}$, $\bm B_1, ..., \bm B_p$, $\rho_{1,1}, ..., \rho_{N,p_N}$, and $\bm H$, that enter the system matrices of the state space model $\bm T$, $\bm Z$, and $\bm H$, in a parameter vector $\boldsymbol{\theta}$. To estimate ${\bm \theta}$ we adopt the approach of HarWei2018a, who derive an analytical solution to the optimization problem of the expected complete Gaussian likelihood function of the state space model, together with a computationally fast combination of the EM algorithm and gradient-based optimization.
In the expectation step of the EM algorithm, we estimate the smoothed states and disturbances, together with the corresponding covariance matrices for a given set of parameters $\hat{\bm \theta}_j$ via the Kalman filter and smoother. The M-step then maximizes the likelihood given the Kalman filter and smoother estimates to obtain $\hat{\boldsymbol{\theta}}_{j+1}$.
After either a convergence criterion is satisfied, or a predefined number of iterations $m$ is reached, the resulting parameter estimates from the EM algorithm $\hat{\bm{\theta}}$ are used as starting values for the maximum likelihood estimation via the BFGS algorithm, which uses the analytical solution for the score vector of HarWei2018a, since the EM algorithm was found to be slow around the optimum.
In case of the stationary factor model in fractional differences, the matrices $\bm T$ and $\bm Z$ are functions of two disjoint parameter spaces and, therefore, a simplification of the EM algorithm is obtained directly by solving the score vector for $\mathrm{vec}(\bm T)$ and $\mathrm{vec}(\bm Z)$ JunKoo2015.
Forecasts are obtained by shifting the system one period ahead and plugging in the smoothed factor estimates from the Kalman filter.
Having discussed the estimation of $f(\bm \chi_t)$ together with the unknown parameters for the three fractionally integrated factor models in sections (ref)--(ref), we investigate their predictive accuracy when neither the DGP, nor the starting values, nor the number of factors, are known to the researcher. For this purpose we study the forecast performance of our three models in a pseudo out-of-sample forecast experiment with an underlying data set for the United States of America that consists of 112 macroeconomic variables and spans from January 1960 to December 2016 McCNg2015.
To compare the forecast performance of the different factor models, we report the resulting mean squared prediction errors (MSPE) for a selected subset of economic variables that represent different segments of the economy.
All forecast models allow for seven common factors in the data, which is suggested by the $PC_{(p3)}$ criterion of BaiNg2002 after deterministic terms have been eliminated from the fractionally differenced data set. Lag lengths of the different AR polynomials for the common factors and idiosyncratic components are chosen via the Bayesian Information Criterion (BIC). For the first forecasting period we obtain starting values for the Kalman filter from the principal components estimator, as described in section (ref). In all subsequent periods the optimized parameters from the preceding step are used as starting values. Finally, the number of iterations of the EM algorithm is set to ten. To distinguish between the first stage and the second stage estimator, we denote the principal components forecasts as PC and the Kalman filter forecasts as KF. Abbreviations for the three fractional factor models are: Dynamic fractional factor model (DFFM) in section (ref), dynamic orthogonal fractional components (DOFC) in section (ref), and dynamic factor model in fractional differences (DFFD) in section (ref).
Forecasts are conducted for horizons $h=1,...,12$ in a recursive window forecast experiment, where the first forecast period is January 2000, whereas the last is December 2016, leading to 204 forecasts for 112 variables and 12 horizons.
The DOFC model introduced in section (ref) includes $r_1=3$ fractionally integrated factors, since a higher number was not found to increase the forecast precision substantially. As a consequence, the number of remaining $I(0)$ factors is set to $r_2=4$. The latter restriction is confirmed by the $PC_{(p3)}$ criterion of BaiNg2002, which suggests four factors after the fractionally integrated factors have been projected out.
A stationary data set for the DFFD model in section (ref) is obtained by estimating the integration order of each $y_{i,t}$ via the exact local Whittle estimator of ShiPhi2005 with a tuning parameter of $0.5$ and taking fractional differences.
In addition, we include four benchmark models to evaluate the forecast performance of the fractional factor models relative to widely used alternatives. The first benchmark is an autoregressive model (AR) where the AR lag order is chosen via the Akaike Information Criterion for each $y_{i,t}$. The second benchmark is a standard approximate dynamic factor model StoWat2002 that is estimated via principal components (PC) based on a pre-differenced data set, i.e. $\Delta^{k_i}y_{i,t+h}=\bm{\Lambda}_i\bm f_{t+h}+\xi_{i,t+h}$, $\bm{\phi}(L)\bm{f}_{t+h}= \bm{\zeta}_{t+h}$, where $\xi_{i,t}$, $\zeta_{j,t}$ are mutually independent and white noise $\forall t=1,...,T$ and $k_i$ is an integer that is taken from McCNg2015. Our third model adds lagged dependent variables to the approximate dynamic factor model. It is given by $\bm{\phi}(L)\bm{f}_{t+h}= \bm{\zeta}_{t+h}$, $c_i(L)\Delta^{k_i}y_{i,t+h}=\bm{\Lambda}_i\bm f_{t+h}+\xi_{i,t+h}$ where $\xi_{i,t}$, $\zeta_{j,t}$ are again mutually independent and white noise $\forall t=1,...,T$. We denote it as PCAR. Finally, the last benchmark is the so-called factor-augmented error-correction model (FECM), which separates the observable variables into two disjoint samples $\bm y = (\bm y^{(1)'}, \bm y^{(2)'})'$ and shrinks the latter sample via principal components to $\hat{\bm f}$. A vector error-correction model is then estimated for $(\bm y^{(1)'}, \hat{\bm f}')'$. Details on the forecast properties are found in BanMarMa2014. Since we only obtain predictions for $\bm y^{(1)}$, the $\textbf{FECM}$ results are only reported in tables (ref) and (ref).
Table (ref) shows for a given forecast horizon $h$ how often each specification leads to the smallest MSPE for all $112$ variables. Hence, it illustrates how frequently fractional factor models are able to outperform widely used forecast methods like autoregressive models and principal components of integer differences. To draw inference on the extent of forecast improvement from the fractional factor models, the tables (ref) and (ref) report the relative MSPE for the twelve depicted variables and for $h = 1, 2, 3, 6, 9, \text{ and }12$. Consequently, they also show how large the forecast accuracy fluctuates for each specification and highlight the robustness of the forecast results when a model is not chosen to be the best one.
\linespread{0.5}{
} We find that fractional factor models tend to outperform classical autoregressive models, pre-differenced principal components models and mixtures of these two model classes. Over all 1344 conducted forecasts, the benchmarks only exhibit a smaller MSPE than the fractional factor models in 357 cases (26.6%), as table (ref) shows. Hence, for the remaining 987 forecasts (73.4 %) the smallest MSPE is achieved by one of the six fractional factor models. Within the benchmarks, one often finds that principal components come with the smallest MSPE. Nonetheless, they are often beaten by one of the fractional factor models. Among those, the dynamic orthogonal fractional components model in state space form produces the best predictions for forecast horizons up to 9 months most frequently.
In addition to the good performance of the DOFC-KF specification, the DFFD models complement the predictive power of fractional factor models. Whenever the DOFC-KF model does not provide the best forecasts, the fractionally differenced models are likely to exhibit the smallest MSPE. Furthermore, principal components are found to perform relatively well at least for smaller forecast horizons when the data is in fractional differences, whereas they are typically beaten by the state space formulation in the DOFC framework. This might be a result of the additional structure that is imposed on the DOFC-KF model via the block-triangular identification of the fractional factors relative to the DOFC-PC case, whereas only little additional structure is imposed on the DFFD-KF specification relative to principal components. For larger forecast horizons, the forecast performance of the DFFD-KF model improves, leading to the highest amount of best predictions for $h = 10, 11, 12$.
\linespread{0.5}{
} We are able to uncover more details about the forecast performance of fractional factor models by having a closer look at the tables (ref) and (ref) that visualize the relative MSPEs for selected variables and forecast horizons $h$. We find that gains from the fractional factor models can be large, relative to the four benchmarks. In many cases, fractional factor models can reduce the MSPE relative to the AR benchmark by more than 25%. For some target variables, the MSPE is cut by half when fractional factor models are used, and reductions of more than 80% are possible.
Within the class of fractional factor models, we find the DOFC-KF specification to perform best. For $h=1,2,3$, the most accurate predictions for the consumer price index, personal consumption index and average hourly earnings are obtained from the DOFC-KF specification, which reduces the MSPE relative to the AR benchmark by more than 50%. In addition, the DOFC-KF specification exhibits the smallest MSPE for the St. Louis adjusted monetary base, total reserves of depository institutions and the S&P500 frequently. The stable and reliable performance of the DOFC-KF forecasts is illustrated by the fact that their largest relative MSPE is 1.29, whereas the smallest relative MSPE is 0.17.
The good performance of the DOFC-KF specification is complemented by the DFFD model. For the industrial production index, the DFFD-PC specification exhibits the smallest MSPE for any forecast horizon. In addition, the DFFD-KF specification produces accurate predictions for the S&P500, average hourly earnings and the US / UK foreign exchange rate. Furthermore, its forecast performance is almost as stable as the DOFC-KF prediction quality.
Finally, the DFFM model, which serves as the most general framework as it nests the two remaining fractional factor model formulations, cannot compete with the other fractional factor models, as its predictive power fluctuates largely. Nonetheless, for larger forecast horizons, the DFFM-KF formulation produces accurate forecasts for the consumer price and personal consumption index.
Note that the only difference between the benchmark PC model and the DFFD-PC specification is the pre-differencing. As one can see, the two models coincide regarding their relative performance to the AR benchmark. The advantages over the AR model are therefore likely to result from cross-sectional dependencies that are detected by the common factors. In addition, the better performance of the DFFD-PC model can be explained by the sensitivity of standard PC methods to spurious coefficients, as Fra2017 argue.
Turning to the DOFC-KF specification, which explicitly models fractional cointegration relations instead of eliminating them as in the DFFD model, we note that the forecast quality of the two models is similar for many predictions. Nonetheless, gains from the DOFC-KF specification relative to the DFFD model can be large, especially in situations where the latter produces a relative MSPE $> 1$. Consider e.g.\ the forecasts for the adjusted monetary base (AMBSL) and the total reserves of depository institutions (TOTRESNS) in tables (ref) and (ref), where the DOFC-KF and the FECM model perform well, wheres the DFFD-KF model yields large MSPEs. While the former two models take cointegration into account, the DFFD-KF model eliminates long-run components by prior differencing and is likely to produce over-differenced short-run components. Hence, the better performance of the DOFC-KF model over the DFFD-KF model is likely to result from cointegration relations and over-differencing of additive short-run factors.
Finally, we want to draw inference on the performance of the fractional factor models during the world financial crisis. By studying the predictive power of the fractional factor models during this period, we shed light on the behavior of this model class when the economy is hit by a large shock and pushed out of its equilibrium growth path. For this purpose, figure (ref) sketches the three step ahead predictions for the twelve selected target variables and the two best performing fractional factor models together with the AR benchmark and the realization of the target variable from January 2007 to December 2011. As the graphs show, the forecast performance of the fractional factor models is not systematically affected by the global financial crisis relative to the AR benchmark. Instead, the forecasts converge towards the realizations of the observable variables rapidly after the crisis. The DOFC-KF forecasts seem to be the least affected by the large shock, as they converge faster towards the observable variables. Furthermore, the AR and DFFD-KF predictions for the adjusted monetary base and total reserves of depository institutions seem to be polluted by the crisis until the end of 2009, which substantiates the relative robustness of the DOFC-KF specification.
We have derived three different fractional factor models that allow the joint modelling of data of different persistence. A two-stage estimator for the fractional factors and model parameters was derived. In a macroeconomic forecast experiment, it was shown that incorporating fractional integration into the class of factor models improves forecast performance substantially. \\ Future research could examine whether a combination of the DOFC model in state space form and a factor model in fractional differences can improve the predictive power of fractional factor models. Furthermore, one could combine principal components and the Kalman filter analogous to BraeKoo2014 by reducing the dimension of a subset of observable variables via principal components in order to speed up the estimation of the parameters. Additionally, fractional factor models could be used to explore common trends and cycles in macroeconomic variables and to identify cointegrated blocks. Finally, future research could address the predictive power of fractional factor models for other data sets and economies. If gains are of similar size as for the US, we are confident that fractional factor models have the potential to become a widely used tool for predicting macroeconomic dynamics.