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.
137,582 characters · 13 sections · 32 citation commands
The Forecasting Performance of Factor Models with Martingale Difference errors
\affil[1]{ University of Rome “Tor Vergata", Italy} \affil[2]{ University of L'Aquila, Italy}
\noindentKeywords: Multivariate Time Series, Factor models, Nonlinear dependence, Martingale difference hypothesis, Factor selection.
This paper investigates the forecasting performance of factor models with martingale difference error (FMMDE, \citeA{lee2018martingale}), a recent advancement in factor modeling literature. The main features of FMMDE are outlined below. Under the assumption that the idiosyncratic component is a martingale difference sequence - a more stringent condition than just assuming that it follows a white noise process - it is possible to show that there exists a contemporaneous linear transformation of the original series such that the resulting series are segmented into two groups, one of being conditionally mean independent upon past information. This property of the model has important implications in terms of optimal prediction, as it reduces the modeling effort to the subset of factors that exhibit some form of serial dependence in the conditional mean. Moreover, employing FMMDE is particularly advantageous when analyzing time series that exhibit nonlinear dependencies, as the model relies upon measures of statistical dependence that are suitable for detecting nonlinear relationships in the data. In this regard, we provide simulated examples that highlight the good performance of FMMDE in terms of forecast accuracy within nonlinear contexts.
We contribute to the existing literature in two respects. First, we introduce a novel technique for determining the dimension of the latent factor process. Our method, in particular, rests upon a sequential testing procedure to identify those transformed series for which the martingale difference hypothesis can be rejected, i.e., those exhibiting some form of conditional mean dependence upon the past. Through an extensive Monte Carlo simulation exercise, we show an improved accuracy in the estimation of the true number of factors as we compare our method to the eigenvalue ratio procedure outlined in Lee and Shao's original paper. The evidence is confirmed in high-dimensional settings when the number of variables greatly exceeds the number of observations.
Second, we assess the forecasting performance of FMMDE to establish whether it can represent a valid alternative when predicting macroeconomic aggregates. For this purpose, we compare its empirical properties to those of the factor model defined in Stock and Watson ( SW hereafter), see \citeA{stock2002forecasting}, \citeA{stock2002macroeconomic}. Furthermore, we provide a comparison with the model proposed by lam2011estimation, hereafter referred to as LYB, which represents an analog of FMMDE in the case where the idiosyncratic term is white noise by construction, and the factors are assumed to exhibit serial linear dependence. The two models, FMMDE and LYB, share a significant number of features, as they both rely on the information obtained from specific autocovariance matrices at non-zero lags. For this class of models, the factors capture all the non-trivial dynamics of the data, while cross-sectional dependence may be explained by both the factors and the idiosyncratic components. This approach differs from the one followed by more traditional techniques, which rely upon the information contained in the variance-covariance matrix. In this case, in fact, the factors capture most of the cross-sectional dependence, while the idiosyncratic terms, not subjected to serial independence constraints, may contain some non-trivial dynamics.
The estimation of FMMDE rests upon the definition of martingale difference divergence matrix, which can be regarded as a generalization of the standard autocovariance matrix based upon measures of nonlinear serial dependence. In particular, the martingale difference divergence matrix encodes the information about the number and form of linear combinations of the original vector of time series that are conditional mean independent of past information. We can retrieve this set of transformations by performing principal component analysis on the sample martingale difference divergence matrix. The estimation involves selecting the value of a parameter to determine the lag at which we truncate the serial information used by the autocovariances. Since the original paper provides no optimal way to determine the parameter, in our empirical analysis, we resort to a cross-validation procedure to determine the most appropriate value to be used at each time interval, following a squared forecast error minimization criterion.
Using a comprehensive US macroeconomic dataset, the FRED-MD \footnote{The dataset, introduced by mccracken2016fred, consists of 123 monthly macroeconomic time series observed from January 1959 to December 2019. We remind to the enclosed Appendix for more information about the nature and composition of the series in the dataset.}, we assess the forecasting performance of the aforementioned techniques by conducting a pseudo-real-time exercise with forecast horizon at 1, 6, and 12 months ahead. In addition to the usual target variables such as Industrial Production, Unemployment Rate, All Employees, and Consumer Price Index, the prediction exercise is also conducted on the remaining variables composing the dataset. In this way, it is possible to get a more general picture of the different sectors beyond the forecasting performance on the individual macroeconomic aggregates.
Comparing results obtained for FMMDE using the sequential testing and eigenvalue ratio procedures for factor selection, our empirical evidence suggests that the former performs better for the Output and Labour Market sectors. On the other hand, the eigenvalue ratio procedure leads to better results in predicting the Prices sector, especially for the medium to long horizons. We attribute this difference in the results to the generally more considerable number of factors retrieved by the sequential testing method, which gives an advantage in predicting real economic activity measures. The previous results are confirmed when looking at the forecast performance of key macroeconomic indicators such as Industrial Production, Unemployment Rate, and Consumer Price Index. Comparing the performance of FMMDE with the factor models introduced by SW, it emerges that the two generally exhibit similar performance for short-term forecasting horizons. However, for the medium to long horizons, using FMMDE with the sequential testing procedure offers significant advantages.
The paper is structured as follows. We introduce the main theoretical properties of the FMMDE and the nonlinear dependence measures used in its definition in Section (ref). In particular, we discuss the martingale difference divergence matrix, clarifying its role as a generalized covariance measure. Section (ref) details the sequential testing methodology employed to estimate the number of factors in the model. We compare the finite sample performance of our method to that of the original eigenvalue ratio methodology through simulation examples. Section (ref) briefly discusses the SW model and its main properties. Additionally, it provides an introduction to the LYB model, detailing its close relationship with FMMDE. We study the models' predictive capabilities in both linear and nonlinear settings. The composition of the dataset, estimation methodology, and forecasting procedure are treated in Section (ref). In the same section, we define the cross-validation strategy used for selecting the value of the lag parameter used by the martingale difference divergence matrix in the estimation of FMMDE. Finally, the forecasting performance of the competing models is discussed in Section (ref), where we report results obtained for both individual series and the entire dataset as an aggregate.
Let $\mathbf{x}_{t}=(x_{1t},x_{2t},...,x_{nt})^{\prime}\in\mathbb{R}^{n}$ be a strictly stationary time series, and $ \mathcal{F} _{t-1}=\sigma\left( \mathbf{x}_{t-1},...,\mathbf{x}_{1},...\right) $ be the past information set. It is well known that the best predictor of $\mathbf{x}_t$, in the mean squared error sense, coincides with the conditional mean $E\left(\mathbf{x}_t| \mathcal{F}_{t-1} \right) $. In the original paper, lee2018martingale demonstrate how it is possible to identify a contemporaneous linear transformation, denoted as $\mathbf{M} $ $\in$ $R^{n\times n}$, such that $\mathbf{Mx}_{t}$ $=$ $\left[ \mathbf{F}_{a,t}^{\prime}\mathbf{,F}_{b,t}^{\prime}\right] ^{\prime}$, where $E\left( \mathbf{F}_{a,t}| \mathcal{F} _{t-1}\right) \neq E\left( \mathbf{F}_{a,t}\right) $ and $E\left( \mathbf{F}_{b,t}| \mathcal{F} _{t-1}\right) $ $=E\left( \mathbf{F}_{b,t}\right) $ with $\mathbf{F}_{a,t} $ $\in$ $\mathbb{R}^{r}$ and $\mathbf{F}_{b,t}\in\mathbb{R}^{n-r}$. In other words, the resulting series can be separated into two parts, with one part being conditionally mean dependent upon $ \mathcal{F} _{t-1}$, and the other being conditionally mean independent upon the past information set. The modeling task for the entire series $\mathbf{x}_{t}$ will be therefore reduced to that for the lower dimensional series $\mathbf{F}_{a,t}$ without any loss in terms of predictive accuracy.
Such a problem can be formulated equivalently in a factor model framework: under the assumption that the variables have a latent factor structure; the approximate factor model is defined as:
where $\mathbf{F}_{t}=(F_{1t},...,F_{rt})'$ are $r<n$ unobserved common factors, $\mathbf{\Lambda}\in\mathbb{R}^{n \times r}$ is the matrix of factor loadings and $\boldsymbol{\zeta}_{t}=(\zeta_{1t},...,\zeta _{nt})^{\prime}\in\mathbb{R}^{n}$ represents the idiosyncratic disturbance not explained by the factors. The factor model with martingale difference error can be expressed as
with $E(\mathbf{x}_{t}| \mathcal{F} _{t-1})=\mathbf{\Lambda}\mathbf{F}_{t}$, where, by construction, the idiosyncratic component $\boldsymbol{\zeta }_{t}$ is a vector martingale difference sequence. Defining $\mathbf{M}$ as the partitioned matrix $\mathbf{M}=\left[ \mathbf{\Lambda },\mathcal{B}\right] ^{\prime }$, for $\mathcal{B} \in \mathbb{R}^{n\times \left( n-r\right) }$ being a matrix for which the relationship $\mathbf{\Lambda }^{\prime }\mathcal{B=}$ $0$ holds true, we can retrieve an estimate of the factor process as a result of obtaining an estimate of the matrix $\widehat{\mathbf{M}}$ as the transformation $\mathbf{\widehat{\mathbf{F}}}_{t}= \widehat{\mathbf{\Lambda}}^{\prime}\mathbf{x}_{t}.$ The identification of $\mathbf{M}$ is made feasible by the introduction of martingale difference divergence matrix, which allows us to express $\mathbf{x}_{t}$ as a nonsingular linear transformation of a $r$-dimensional dynamically dependent common factor process, and a martingale difference error component of dimension $(n-r)$. Given matrix can be regarded as a multivariate extension of the notion of martingale difference divergence, introduced in shao2014martingale as a measure of conditional mean independence between random vectors and defined as (the non-negative square root of)
where $G_{h}\left( \mathbf{s}\right) =$ cov$\left( x_{h,t},\exp i\left\langle \mathbf{s},\mathbf{x}_{t-j} \right\rangle \right) $ and $c_{n}=\pi ^{\left( 1+n\right) /2}/\Gamma \left( \left( 1+n\right) /2\right) $.\footnote{Here we denote $i=\sqrt{-1}$. For a complex-valued function $f(.)$ the complex conjugate of $f$ is denoted by $f^{\ast}$ and $|f|^{2}=ff^{\ast}.$} \footnote{Martingale difference divergence can be regarded as an extension of the concept distance covariance (szekely2007measuring) using a similar weighting function under the integral.} In our case, martingale difference divergence can be used to measure the departure from the conditional mean independence relationship $E\left( x_{h,t}|\mathbf{x}_{t-j}\right) =E\left( x_{h,t}\right)$ almost surely, for a fixed $j\in\mathbb{N}^{+}$, and for each $h$-th element in the vector $\mathbf{x}_{t}$. In particular, martingale difference divergence is defined so that conditional mean independence is verified if and only if $ MDD\left( {x}_{h,t}|\mathbf{x}_{t-j}\right)^{2} =0$. The function $G_{h}\left( \cdot \right) $ can be interpreted as a generalized covariance function between $h$-th element in $\mathbf{x}_{t}$ and the complex exponential of its lag: this allows to express dependence conditions in terms of characteristic functions of the variable in question. As discussed in hong1999hypothesis, characteristic functions can be used to check for a variety of dependence conditions different from simple correlation, allowing, in particular, to assess the existence of nonlinear dependence relations. It is well known that, in the case of non-Gaussian and nonlinear time series, the autocorrelation function cannot fully characterize dependence structures. In this sense, MDD allows us to assess the existence of conditional mean dependence in the absence of linear relationships when only nonlinear dependence is present.
In the case of multivariate time series objects, the martingale difference divergence matrix is defined as the $n\times n$ matrix
where $\mathbf{G}\left( \mathbf{s}\right) =$ cov$\left( \mathbf{x}_{t},\exp i\left\langle \mathbf{s},\mathbf{x}_{t-j}\right\rangle \right) =\left( G_{1}\left( \mathbf{s}\right) ,...,G_{n}\left( \mathbf{s}\right) \right) ^{\prime}$ for $\mathbf{s}\in \mathbb{R}^{n}$. For $h=1,...,n$, elementwise, we have\ $G_{h}\left( \mathbf{s}\right) =$ cov$\left( x_{ht},\exp i\left\langle \mathbf{s},\mathbf{x}_{t-j} \right\rangle \right) $. Accordingly, its $(i,i)$th entry equals to $MDD(x_{it}|\mathbf{x}_{t-j})^{2}$.
Assuming that $E\left( |\mathbf{x}_{t}|_{n} ^{2}+|\mathbf{x}_{t-j}|_{n}^{2}\right) $ $<\infty$, the matrix is real, symmetric, and positive semi-definite, so that, for any real matrix $\mathbf{\Lambda}\in \mathbb{R}^{n\times r}$ we observe the relationship $ \mathbf{MDDM}\left( \mathbf{\Lambda}^{\prime}\mathbf{x}_{t}|\mathbf{x}_{t-j}\right) =\mathbf{\Lambda}^{\prime}\mathbf{MDDM}(\mathbf{x}_{t}|\mathbf{x}_{t-j})\mathbf{\Lambda }$. There exist, subsequently, $\left(n-r\right)$ linearly independent combinations that are conditionally mean independent of $\mathbf{x}_{t-j}$, i.e., $MDD\left( \mathcal{B}_{i}^{\prime}\mathbf{x}_{t}|\mathbf{x}_{t-j}\right) =0$, $i=1,...,n-r$, if and only if Rank$\left( \mathbf{MDDM}\left( \mathbf{x} _{t}|\mathbf{x}_{t-j}\right) \right) =r$. In this sense, the martingale difference divergence matrix can be employed to determine the set of linear combinations that are conditional mean independent of past information. Let $\left\{ \lambda_{i},\gamma_{i}\right\} _{i=1}^{n}$ \ be the $n$ pairs of eigenvalues and eigenvectors of $\mathbf{\Gamma}_{k_{0}}$, the matrix $\mathbf{M}=\left[ \mathbf{\Lambda},\mathcal{B}\right] ^{\prime}$ can be obtained by the spectral analysis of $\mathbf{\Gamma}_{k_{0}}$ as the matrix of its eigenvectors, so that $\mathbf{M=}\left[ \gamma _{1},...,\gamma _{n}\right] $. In particular, $\mathbf{\Lambda}$ will be the matrix collecting eigenvectors corresponding to the $r$ factors that are not martingale difference sequences; conversely, $\mathcal{B}$ will correspond to the remaining $n-r$ eigenvectors for which $MDD\left( \mathcal{B}_{i}^{\prime }\mathbf{x}_{t}|\mathbf{x}_{t-j}\right) =0$, for every $i=r+1,...,n,$ and for any $j \in \mathbb{N}^{+}$.
For a given sequence of observations $\{\mathbf{x}_t \}^{T}_{t=1}$, defining $\mathbf{z}_{t}=\mathbf{x}_{t-j}$, for $j=1,...,T-1$, sample martingale difference divergence matrix can be expressed as
From a practical standpoint, as we have only a limited number of observations available, we need to approximate the past information set,$\ \mathcal{F}_{t-1} $, with its finite sample analogue $\mathcal{F}_{t-1,t-k_{0}}=\sigma\left( \mathbf{x}_{t-1} ,...\mathbf{x}_{t-k_{0}}\right)$, with $k_{0}$ being a pre-specified integer. For this purpose, the original paper introduced the concept of cumulative martingale difference divergence matrix to quantify the conditional mean independence of $\mathbf{x}_{t}$ on its recent past $\mathcal{F}_{t-1,t-k_{0}}$, defined as
This reduces the analysis to that of the matrix $\widehat{\mathbf{M}}$ $=\left[ \widehat{\mathbf{\Lambda}},\mathcal{\widehat {\mathcal{B}}}\right] ^{\prime}$ of the eigenvectors of $\mathbf{\widehat{\Gamma}}_{k_{0}}$. From now on, we mainly refer to the vector of transformed variables $\widehat{\mathbf{P}}_t = \widehat{\mathbf{M}} \mathbf{x}_t$. By construction, such linear transformation gives rise to a separation between the estimated factors and a number of MDS components, so that $\widehat{\mathbf{P}}_t = \left[ \widehat{\mathbf{F}}_t^{\prime},\widehat {\mathbf{E}}_t^{\prime} \right] ^{\prime}\in$ $\mathbb{R}^{n}$, with $\widehat{\mathbf{F}}_{t} $ $\in$ $\mathbb{R}^{r}$ and $\widehat{\mathbf{E}}_{t}\in\mathbb{R}^{n-r}$.
The methodology for factor selection adopted in the original paper is based upon the eigenvalue ratio estimator defined in lam2011estimation, and lam2012factor. Let $\left\{ \widehat{\lambda} _{i},\widehat{\gamma}_{i}\right\} _{i=1}^{n}$ be the $n$ pairs of eigenvalues and eigenvectors obtained by the spectral analysis of $\mathbf{\widehat{\Gamma}}_{k_{0}}$ rearranged in order of decreasing eigenvalues. It is possible estimate $r$, i.e., the rank of $\mathbf{\Gamma}_{k_{0}}$, using the ratio based estimator
In the original paper, no formal rule is given for selecting the value of $R$. In general, such value can be defined as a fraction of the total number of variables $n$.
We propose a novel approach to determining the number of factors in the model. In particular, we exploit the nature of FMMDE, which allows obtaining separation of the unobservable series into different groups according to whether they are conditionally mean dependent upon past information or not. Specifically, we test the dynamic properties of the individual transformed series $\widehat{\mathbf{P}}_{t}=\left(\widehat{P}_{1,t},...,\widehat{P}_{n,t}\right)$ in a sequential scheme, aiming to identify the subset of those that are not a martingale difference sequence. In performing this operation, we start from the first few elements in $\widehat{\mathbf{P}}_t$, i.e., those associated with the largest estimated eigenvalues of $\mathbf{\widehat{\Gamma}}_{k_{0}}$. The testing sequence stops as soon as we are able to identify the first transformed series, $\widehat{P}_{i,t}$, for which the null hypothesis of conditional mean independence upon the past cannot be rejected. In this case we have $\widehat{r}=i-1$, with $\widehat{P}_{i-1,t}$ being the last element in $\widehat{\mathbf{F}}_t$ and $\widehat{P}_{i,t}$ being the first element in $\widehat{\mathbf{E}}_t$. If all transformed series are a martingale difference sequence, the testing process stops at the first iteration, $\widehat{r} = 0$, and $\widehat{\mathbf{F}}_t=\mathbf{0_v}$.
A similar approach to determining the number of factors is discussed in gao2021modeling, which considers a generalization of the model of lam2011estimation. In particular, estimation of the loading matrix associated with the common factors is based upon the eigenanalysis of certain autocovariance matrices, which, similarly to what happens with MDDM, allow for obtaining a separation between a set of serially dependent factors (defined in this case by a non-null degree of autocorrelation) and a white noise idiosyncratic component. The sequential testing procedure described by gao2021modeling differs from ours for various reasons. Firstly, they investigate an alternative null hypothesis, namely the absence of serial correlation in the idiosyncratic term. In this regard, the consideration of the martingale difference hypothesis can be viewed as more comprehensive, since it allows for the possibility of observing white noise series that are not strictly martingale difference sequences. Consequently, even in the absence of linear dependence, such a series may exhibit some form of nonlinear dependence in the conditional mean. Furthermore, unlike the test statistics they utilize, the chosen test statistic in our study does not have a known limiting null distribution. As a result, we need to employ a bootstrap procedure to estimate its distribution.
We make use of the test statistic for the martingale difference hypothesis introduced by wang2022testing, defined as
where $\omega_{j}=\left( N-j+1\right) /\left( Nj^{2}\right) $ $\ $and $\left\vert \left\vert \cdot\right\vert \right\vert _{F}$ indicates the Frobenius norm.\footnote{Among all the different test specifications proposed in the original paper, we chose $\widehat{T}_{wn}^{F}$ as our statistic of reference as it represents the best compromise in terms of size and power properties.} For the rest of this section, to avoid using too heavy notation, we prefer, from now on, not to index the test $\widehat{T}_{wn}^{F}$ at $i$, as it is the case in Equation (ref), even though we are testing the $i$-th element in a sequential scheme. As the asymptotic null distribution of $\widehat{T}_{wn}^{F}$ is non-pivotal, the critical values of the test statistics need to be approximated by an appropriate bootstrap methodology. Following the original paper, we adopt a fixed-design wild bootstrap procedure, see wu1986jackknife and liu1988bootstrap, where, given a sequence $\left\{ w_{t}^{\ast }\right\} _{t=1}^{T}$ of iid auxiliary random variables, and a random sample $\left\{z_t\right\}^{T}_{t=1}$, the bootstrap sample $\left\{z^{*}_t\right\}^{T}_{t=1}$ is generated as $z^{*}_t=z_tw^{*}_t$ for each $t=1,...,T$. In our specific case, we assume $\left\{ w_{t}^{\ast }\right\} _{t=1}^{T}$ to be a sequence of iid Bernoulli variates with distribution $ P\left[ w_{t}^{\ast }={\scriptstyle \left(1-\sqrt{5}\right)/2}\right] ={\scriptstyle\left(\sqrt{5}+1\right)/2 \sqrt{5},}$ and $ P\left[w_{t}^{\ast }={\scriptstyle\left(1+\sqrt{5}\right)/2}\right] ={\scriptstyle\left(\sqrt{5}-1\right)/2 \sqrt{5}.} $ Analogously to the methodology implemented in escanciano2006goodness and wang2022testing for testing the presence of nonlinear dependence in the residuals of a time series model, our methodology relies upon re-estimating the factor model for each of the bootstrap replicates using the bootstrap sample. The resulting vector, $\widehat{\mathbf{P}}^{*}_{t}$, will be then used to compute the bootstrap test statistic $\widehat{T}^{*}_{wn,b}$, so that, for a total of $B$ bootstrap replicates, we obtain $\left\{\widehat{T}^{*}_{wn,b}\right\}^{B}_{b=1}$, which allows us to approximate the test distribution.
To describe the step-by-step functioning of our testing methodology, we start by illustrating the initial operation of testing the transformed series $\widehat{P}_{1,t}$, to then consider the generic $i$-th series $\widehat{P}_{i,t}$. We denote as $\mathbf{\widehat{P}}_{i:n,t}$ the vector collecting the elements of $\mathbf{\widehat{P}}_t$ from the $i$-th to the $n$-th element. The sequential testing procedure can be summarized as follows
In this section, we report the results of some simulation experiments to assess the performance of the proposed factor selection methodology. Results for the testing procedure are compared to those obtained for the original eigenvalue ratio methodology, defined in Equation (ref). The baseline model for all simulations will be
Following gao2021modeling, we assume that the factor process is characterized by an autoregressive structure
where $\mathbf{\Pi}$ is a $r\times r$ diagonal matrix with its diagonal elements being drawn independently from $U\left( 0.5,0.9\right) ,\mathbf{\eta} _{t}\sim N\left( 0,I_{r}\right) $. The entries of the factor loading matrix $\mathbf{\Lambda}$ are all drawn independently from $N\left( 0,1\right)$. Throughout the analysis we consider three data-generating processes (DGPs) for the idiosyncratic component similar to those in bai2002determining and alessi2010improved.
Tables from (ref) to (ref) report values of the empirical probabilities computed, for each methodology, in the process of determining the correct number of factors. These probabilities are obtained considering 1000 independently simulated samples, counting the number of successful factor identification attempts. Results are displayed for different values of $\theta$, $DGPs$, and for a true number of factors ranging from 2 to 5. Regarding our sequential testing procedure, we assume the fixed truncation lag $M$ - the argument used by the statistic $\widehat{T}_{wn}^{F}\left( M\right)$ - to be equal to 6. Empirical $p$-values are computed over 499 bootstrap replications for a significance level $\alpha=5\%$. In the case of the eigenvalue ratio methodology, the value of the parameter $R$ is set to $n/3$ in all simulations. Moreover, we set $k_{0}=1$ for simplicity: our evidence suggests that results are not affected by selecting alternative values for $k_{0}$ when considering this class of models. All results are obtained after standardizing the series. The benefits of utilizing the sequential testing methodology are particularly evident when analyzing $DGP_3$ with cross-sectional correlations among idiosyncratic components and when considering higher values of theta. The most significant results, relative to the eigenvalue ratio methodology, are in general attained for smaller values of $T$ and $n$, as it is apparent from Table (ref) where we display results obtained for $n=50$, $T=50$ (upper part of the table), and $T=100$ (lower part of the table). It is worth noticing that better results are also obtained for $n=100$ and $T=100$ (lower part of Table (ref)): these values of $T$ and $n$ lie approximately in the same range as those that we select for our empirical study.
We extend our simulation study by considering the model reported in Example 6 of lee2018martingale's paper. The three-dimensional factor space, with $\mathbf{F}_{t}=\left( F_{1t},F_{2t},F_{3t}\right) ^{\prime}$, is generated from the nonlinear process $\omega_{t}$, defined by
with
Concerning the factor loading matrix $\mathbf{\Lambda}$, the first $n/2$ elements of each column are iid $U\left( -2,2\right) $ and are kept fixed once generated and the other elements are set to be zero. The parameter $\theta$ is now assumed to be always equal to one, while the error term $\boldsymbol{\zeta}\mathbf{_{t}}$ is a random sample of $N\left( 0,0.25\mathbf{\Sigma}\right) $ independent of $\mathbf{F}_{t}$. Finally, the covariance matrix $\mathbf{\Sigma=}\left( \sigma_{i,j}\right) _{i,j=1}^{n}$ is generated as
for Hurst exponent $H=0.9$. As mentioned in the original paper, the above assumptions about the dependence structure of error term make it possible to address the large $n$ case, where $n$ can exceed the number of available observations $T$. Table (ref) displays results for values $k_0 \in \left\{1,10,25\right\}$. In particular, the parameter has a non-negligible influence on the factor selection process: increasing $k_0$ improves the performance of the sequential testing method while worsening the performance of the eigenvalue ratio at the same time. Still, at the two extremes, the methods provide comparable results, with the eigenvalue ratio performing best for $k_0=1$, and the sequential testing method for $k_0=25$.
In this section, we investigate the forecasting performance of the FMMDE model in a simulated environment. The exercise is based on a systematic comparison with the factor model defined in Stock and Watson, see stock2002forecasting and stock2002macroeconomic, which can be considered a well-established benchmark in the literature and can therefore be used as a standard of reference. Letting $\mathbf{S}=T^{-1}\sum_{t}\mathbf{x} _{t}\mathbf{x}_{t}^{\prime}$ and denoting the spectral decomposition of the covariance matrix by
where $\mathbf{V}=(v_{1},\dots,v_{n})$ is the $(n\times n)$ matrix of orthonormal eigenvectors, $\mathbf{V^{\prime}V=I_{N}}$, \ we obtain an estimate of the $r$ common factors as
with $\mathbf{D}$ being the matrix containing the square root of the ordered eigenvalues. For completeness, we also discuss the relationship of FMMDE with the factor model introduced in lam2011estimation, which can be regarded as an analog for FMMDE dealing with the case of serially linearly dependent factors. Such a model is based on the eigendecomposition of the matrix
where $\mathbf{cov}\left( \mathbf{x}_{t},\mathbf{x}_{t-j}\right) $ is a $\left(n \times n\right)$ matrix with $\left(i,h\right)$-th entry being equal to $\text{cov}\left( x_{i,t},x_{h,t-j}\right)$. In particular, $\mathbf{\mathcal{L}}_{k_{0}}$, again a real, symmetric, and positive semidefinite matrix, can be considered as a linear analog of the matrix $\Gamma_{k_{0}}$ used in FMMDE. Less stringent assumptions for the error term are being made in this case, as this is defined to be a white noise sequence. In particular, the model offers a way to retrieve a decomposition between a set of serially correlated factors and a white noise idiosyncratic component. Estimation of the factor loading space works analogously, being based upon the eigenvectors of $\mathcal{L}_{k_{0}}$. The main difference between the models lies in the serial dependence structure characterizing the factors. In particular, while LYB requires factors to exhibit some degree of autocorrelation, this is not the case for FMMDE, as these can show some form of (nonlinear) dependence in the mean even when they are a white noise sequence. This is reflected in the nature of the matrices employed for estimating the models: the cumulative linear matrix $\mathbf{\mathcal{L}}_{k_{0}}$ encodes the linear dependence relationships, whereas cumulative MDDM more generally characterizes conditional mean independence.
We provide an example showing the forecasting advantage granted by FMMDE when factors are characterized by nonlinear serial dependence. For this reason, we consider again the nonlinear DGP seen in Section (ref). For completeness, we additionally consider the linear DGP defined in Example 5 in the original paper. In particular, the process $w_{t}$ generating the factors is assumed to follow a linear MA(1).
where $\boldsymbol{\zeta}\mathbf{_{t}}$ is a random sample of $N\left( 0,\mathbf{\Sigma}\right) $ and is independent of $\mathbf{F}_{t}$, where $\mathbf{\Sigma=}\left( \sigma_{i,j}\right) _{i,j=1}^{n}$ is defined exactly as in Equation (ref).
Assuming that the exact number of factors is known a priori (so that no factor selection procedure is involved), we evaluate each model's one-step-ahead forecast performance. Predictions are obtained for each $j$-th element in the vector $\mathbf{x_{t}}$ using a simple equation of the form
where the correct number of underlying factors, indexed at $i$ and assumed to be known a priori, is employed. The superscript $m$ stands for the model for which the factor component $\widehat{\mathbf{F}}_{t}$ is estimated. We use the first $(T-1)$ observations for each simulated sample to estimate the underlying factor structure; this is then used to obtain the one-step ahead forecasts. For each model, $m$, forecast errors obtained for the individual series are then averaged as follows
where the subscript $S$ indicates that we are considering $S$-th simulated sample. Forecast results, reported in Table (ref) are expressed in terms of the ratio
where $m \in \left\{SW,LYB\right\}$, and results are averaged over 1000 simulated samples. In particular, it emerges how both the SW and LYB models have a predictive advantage over FMMDE when the factor DGP is assumed to be linear. Varying the number of variables $n$ or sample size $T$ does not have a relevant influence on the results. The table highlights the benefits of utilizing FMMDE in the presence of significant nonlinearities in the underlying factor structure. In such cases, the FMMDE outperforms the SW and LYB models. This advantage becomes more pronounced with larger values of $T$ and $n$.
This section aims to assess the predictive capabilities of FMMDE equipped with the factor selection procedure based on the sequential testing procedure (hereafter denoted as FMMDE$_{ST}$) on the macroeconomic aggregates from themccracken2016fred dataset. The performance of FMMDE$_{ST}$ is compared to that of the same model when factor selection is performed according to the eigenvalue ratio methodology (hereafter FMMDE$_{\lambda}$). The two models, SW and LYB, are also included in the analysis: the optimal number of factors is determined following bai2002determining $IC_{p2}$ criterion in the case of SW model, whereas the eigenvalue ratio methodology is used in the case of LYB model (denoted as LYB$_{\lambda}$). All results are compared to those of a univariate autoregressive model, which is used as a benchmark for all empirical comparisons.
We consider an $n$ dimensional stationary process $\mathbf{x}_{t} = (x_{1t},\dots,x_{nt})$ forming a large dataset of economic indicators. Let $y_{t}$ be the variable to be predicted, part of the same set of variables: we are interested in forecasting $y_{t}$ at horizon $h$ by using all the information made available by $\mathbf{x}_{t}$. The dataset consists of monthly observations on 123 U.S. macroeconomic time series observed from January 1959 to December 2019. The time series are grouped into eight categories listed in the Appendix. As reported in mccracken2016fred, the series are all transformed to be stationary by taking first or second differences, logarithms, or first or second differences of logarithms. After transforming the series, the available data range from February 1962 to December 2019, for a total of $T = 695$ observations. Let $y^{h}_{t+h}$ denote the variables to be predicted in a $h$-period ahead forecast. Following stock2002macroeconomic, when predicting real activity variables, $y^{h}_{t+h}$ is defined as the $h$-period growth at an annual rate. As an example, let $IP_t$ the level of the Industrial Production index we obtain
where $y_t = 1200*\text{ln}(IP_t/IP_{t-1})$. When predicting the nominal price and wage series, $y_{t+h}$ is defined as the $h$-period growth of monthly changes. As an example, let $CPI_t$ be the level of the Consumer Price Index; we obtain
where $y_t = 1200*\Delta \text{ln}(CPI_t/CPI_{t-1})$. The forecasts are obtained through a pseudo-real-time forecasting procedure, adopting a rolling window scheme for estimation and prediction. Specifically, we opt for a rolling 10-year window.\footnote{ At the time we started the exercise, the dataset used referred to the vintage dated March 2021. We decided to limit our analysis to December 2019 because we intended to measure the predictive capabilities of FMMDE without interfering with COVID-SARS19 pandemic effects.} Based on the estimated factors, we implement the following Diffusion Index forecasting equation
where $m$ stands for the for model for which the factor component $\widehat{\mathbf{F}}_{t}$ is estimated, so that $m \in \{\text{SW}, \text{FMMDE}_{ST}, \text{FMMDE}_{\lambda},\text{LYB}_{\lambda}\}$, $\mathbf{z}_{t} = \left(\widehat{F}_{1t}^{m},\dots,\widehat{F}_{rt}^{m}, y_{t}, \dots y_{t-p+1}\right)'$ and $\boldsymbol{\vartheta}$ is a $q$-dimensional vector of real parameters, $\hat{\boldsymbol{\vartheta}} = (\widehat{\alpha }_{h}, \widehat{\beta } _{h1}, \dots, \widehat{\beta } _{hr}, \widehat{\gamma }_{h1}, \dots, \widehat{\gamma }_{hp} )$. We select the order $p$ of the autoregressive component by Bayesian Information Criterion (BIC), with $0<p\leqslant 12$. As mentioned, forecasting results obtained using diffusion indexes are compared to those obtained from a univariate autoregression, the lag order $p$ being analogously estimated by BIC, with $0<p\leqslant 12$.
The number of factors to be employed in the prediction equation is determined at each iteration of the forecasting exercise, contextually to the estimation of the latent factor space. However, since $k_{0}$ can influence the rank of matrix $\mathbf{\Gamma}_{k_0}$ and $\mathbf{\mathcal{L}}_{k_{0}}$, for FMMDE and LYB respectively, the Sequential Testing and Eigenvalue Ratio selection procedures are repeated for different values of $k_{0}$ selecting the optimal value through a cross-validation strategy tailored to the case of time series.
The main idea behind cross-validation is to use a sample of past predicted values, referred to as validation set, to perform the selection of the value $k_{0}$ to be employed in the proper out of sample prediction. In particular, given a discrete set of possible values for $k_{0}$, we select the parameter value that minimizes the mean squared forecast error computed over the set of predictions included in the validation set. This value, denoted as $k^{*}_{0}$, will be used to obtain the first prediction outside the validation set. In particular, let $k_0 \in \{1,\dots,20\}$ be the discrete set over which to choose $k^{*}_0$. We indicate validation and estimation sets as $\mathcal{I}_{\upsilon,k_0}^{\mathcal{V}}$ and $\mathcal{I}_{\omega}^{\mathcal{E}}$, respectively. Moreover, $\omega = 119$ will be the size of the estimation window throughout the analysis, while $\nu=49$ is the size of the validation window. The procedure is based on a rolling estimation and validation window and starts by estimating the model with the first 119 observations, namely $\mathcal{I}_{\omega}^{\mathcal{E}} = \left\{\mathbf{z}_{\tau},y_{\tau}\right\}_{\tau=Feb:62}^{Dec:71}$, which produces the forecast $\hat{y}^{h,m}_{Dec:71+h,k_0}$. Subsequently, the observation $y_{Jan:72}$ is included in the estimation sample, namely $\mathcal{I}_{\omega}^{\mathcal{E}} = \left\{\mathbf{z}_{\tau},y_{\tau}\right\}_{\tau=Mar:62}^{Jan:72}$, and the model is estimated again to get the forecast $\hat{y}_{Jan:72+h,k_0}$. The process is repeated until we have a prediction for all $49-h$ out-of-sample observations that will form the first validation set, $\mathcal{I}_{\upsilon,k_0}^{\mathcal{V}} = \left\{y_{\tau},\hat{y}^{h}_{\tau}\right\}_{\tau=Jan:72+h}^{Dec:75}$. The cross-validation mean squared forecast error (CVMSE) is then computed to select the best value for $k_0$, namely $$k^{*}_{0} = \underset{k_0}{\operatorname{argmin}} \; CVMSE(k_0),$$ where $$CVMSE({k_0}) = (49-h)^{-1} \sum_{\tau = Jan:72+h}^{Dec:75} \left(y_{\tau}- \hat{y}^{h,m}_{\tau,k_0}\right)^{2},$$ then the selected $k_{0}^{*}$ is used to forecast the observation $y_{Dec:75+h}$, namely $\hat{y}_{Dec:75+h,k_{0}^{*}}$. To construct the forecast for observation $y_{Jan:76+h}$, we update the validation set, namely $\mathcal{I}_{\upsilon}^{\mathcal{V}} = \left\{y_{\tau},\hat{y}^{h}_{\tau,k_0}\right\}_{\tau=Feb:72+h}^{Jan:76}$, obtaining a new value for $k^{*}_0$ which is used to calculate the second forecast, namely $\hat{y}_{Jan:76+h,k_{0}^{*}}$. The process is repeated until the end of the sample is reached, producing a total of 516 forecasts, ranging from January 1976 to January 2019 for $h=1$, June 1976 to June 2019 for $h=6$, and December 1976 to December 2019 for $h=12$. The predictive performance of the models is finally evaluated in terms of mean square forecast error, $MSFE_{h}^{m}$, defined as \[ MSFE_{h,k^{\star}_{0}}^{m}=\left(T_{1}-T_{0}\right)^{-1}\sum_{\tau=T_{0}+h}^{T_{1}}\left( y_{\tau}-\hat{y}_{\tau,k^{\star}_{0}}^{h,m}\right) ^{2}, \] where $T_{0}$ is the first point in time for out-of-sample evaluation and $T_{1}$ is the last point in time for which we compute MSFE at horizon $h$.
The results of the empirical exercise are reported in Tables (ref) and (ref). In table (ref), the results pertain to some key variables of the US economy, while table (ref) presents aggregated results considering all the series comprising the 8 sectors of the dataset. The models are compared at horizons $h=1$, $6$, and $12$, in terms of the ratio
where rMSFE is the relative mean square forecast error with respect to AR(BIC) model and $ m \in \{FMMDE_{ST}, FMMDE_{\lambda}, LYB, SW\}$.
This section reports the results of the forecasting analysis for the Industrial Production index (IP), Unemployment Rate (UR), Nonfarm Employment (PAYEMS), Consumer Price index inflation (CPI), and Core CPI inflation (CoCPI).
By analyzing the results obtained for single variables, reported in Table (ref)-Panel A, we observe an improvement of the FMMDE$_{ST}$ procedure compared to FMMDE$_{\lambda}$ for the variables IP, UR, and PAYEMS. This result is confirmed across all forecast horizons. However, if we look at the results for the nominal variables, i.e., CPI and cCPI, in this case, FMMDE$_{\lambda}$ has an advantage, especially for $h=12$. When comparing different factor models (martingale difference and SW) they both show a clear improvement over the forecasts obtained with an AR model, especially for the variables IP, UR, and PAYEMS. However, factor models tend to be less accurate in forecasting the US nominal variables. These results are in line with previous findings by eickmeier2008successful. The previous results are further supported by the diebold1995comparing test, as shown in Table (ref) - Panel B. The test confirms that both FMMDE$_{ST}$ and SW significantly outperform AR(BIC) in terms of predictive ability, though this is less evident for FMMDE$_{\lambda}$. In the case of the CPI, we cannot reject the null hypothesis of equal predictive ability for all diffusion index models considered, except for FMMDE$_{\lambda}$ when $h=12$.
As further evidence, to address potential instabilities that might affect the Diebold-Mariano test, we adopt the fluctuation test procedure introduced by giacomini2010forecast (hereafter GR) to locally evaluate the forecasting performance among the procedures at specific points in time.\footnote{This approach is detailed in giacomini2010forecast as a simple way to test against the null of equal local performance of two forecasting methods.} The evolution of the relationship between the forecast performance of FMMDE$_{ST}$ and that of FMMDE$_{\lambda}$, and SW is analyzed in detail in Figure (ref) and (ref) respectively.\footnote{We have not included the comparison with LYB as it exhibits very similar performance with FMMDE$_{\lambda}$} The two figures display results obtained for the fluctuation test for the time series of IP, UNRATE, and CPI at horizons $h=1,12$. The solid line represents the graph of the difference between the square forecast error obtained for the two models, normalized by its estimated standard deviation and smoothed by a centered moving average of length $L = 61$, with the coefficients equal to $1/L$. The zero horizontal line is indicative of the two models having equal performance, while the dotted lines indicate the 5$\%$ critical values, so that FMMDE$_{ST}$ outperforms (underperforms) the other models locally, at the 5$\%$ significance level when the solid line is below (above) the lower (upper) dashed line.\footnote{Because the moving averages are of length 61 and centered, the first and last 30 values are not computed or graphed.}
As illustrated in Figure (ref), for variables such as IP and UNRATE, FMMDE$_{ST}$ and FMMDE$_{\lambda}$ exhibit comparable performance throughout most of the evaluation period. However, during the years following the Great Recession and subsequent recovery (starting from the second half of 2009 until around the end of 2015), FMMDE$_{ST}$ has the best performance. Regarding CPI, FMMDE$_{\lambda}$ has the best results in the mid-2000s for $h = 1$ and during the years between 2011 and 2015 for $h = 12$. Regarding the comparison with the SW model, as illustrated in Figure (ref), FMMDE$_{ST}$ exhibits a significant advantage over the SW model during the period that coincides with the Great Recession for both forecast horizons, with the exception of IPI when $h=12$.
We report the distribution of MSFE relative to AR for a specific group of variables in Table (ref). In particular, we consider seven distinct groups: Output & Income (19 series in total), Labor Market (31 series), Housing Market (10 series), Consumption, Business Inventories & Orders, Money & Credit (12 series), Interest & Exchange Rates (21 series) and Prices (16 series). The results previously observed for individual series are confirmed. FMMDE$_{ST}$ has an advantage over FMMDE$_{\lambda}$ in the prediction of real economic variables (roughly represented by the first four groups in the dataset); this is partly balanced by the higher accuracy displayed by FMMDE$_{\lambda}$ in forecasting nominal variables (represented by the last three groups).
When instead the different procedures are compared, we notice that FMMDE$_{ST}$ delivers more accurate forecasts than FMMDE$_\lambda$ for more than half of the series for the first four groups at all considered horizons. This result extends to more than three-quarters of the series for $h=6,12$. No significant differences emerge in the performance of FMMDE$_\lambda$ and LYB$_\lambda$. The performance of FMMDE$_{ST}$ appears to align with that of the SW model in the prediction of the first four groups, even though FMMDE$_{ST}$ has an advantage for $h=1$, more prominently in the case of the fourth group. All the diffusion index models provide a better forecasting performance than the AR model for more than half of the series in the first four groups, except for FMMDE$_\lambda$ and LYB$_\lambda$ when it comes to longer horizons ($h=12$). On the other hand, FMMDE$_{ST}$ is less performing compared to FMMDE$_\lambda$ and LYB for the groups: Money & Credit, Interest & Exchange Rates, and Prices.
First, we observe that our selection procedure is much less conservative than the eigenvalue ratio methodology in selecting the number of factors, which explains the results obtained with our forecasting experiment see Figure (ref). Selecting a larger number of factors gives a clear advantage in predicting real variables. This result, however, comes at the cost of a lower accuracy in the prediction of nominal variables.
Second, the cross-validation selection of the lags used by the matrix $\widehat{\boldmath{\Gamma}}_{k0}$ influences the factor selection procedure. We can deduce this from the cross-sectional variability observable on the vertical axis of the plot for any fixed time. However, such an effect is limited to specific time intervals, and the degree of variability is limited, as predicting individual variables leads us, in general, to consider latent factor spaces of similar dimensions.
In this paper, we have examined the predictive performance of factor models with martingale difference error (FMMDE). As a primary contribution, we have introduced a novel procedure for selecting the optimal number of factors for FMMDE. Through simulation experiments, we show the good finite sample properties of our proposed methodology. Furthermore, we have evaluated the predictive performance of FMMDE in comparison to other models, such as SW and LYB, using real data. Specifically, we have considered a large dataset of monthly macroeconomic and financial series for the US economy (FRED dataset). The results indicate that, overall, the competing methods exhibit comparable forecasting performance. However, for real activity variables, such as Industrial Production and Unemployment, our selection procedure offers advantages to the FMMDE method, particularly for medium to long-term forecasting horizons.