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.
72,895 characters · 19 sections · 98 citation commands
sparseDFM: An R Package to Estimate Dynamic Factor Models with Sparse Loadings
\affil{Department of Mathematics and Statistics, Lancaster University, Lancaster, LA1 4YW, United Kingdom}
\abstract{ sparseDFM is an R package for the implementation of popular estimation methods for dynamic factor models (DFMs) including the novel Sparse DFM approach of mosleySAC. The Sparse DFM ameliorates interpretability issues of factor structure in classic DFMs by constraining the loading matrices to have few non-zero entries (i.e. are sparse). mosleySAC construct an efficient expectation maximisation (EM) algorithm to enable estimation of model parameters using a regularised quasi-maximum likelihood. We provide detail on the estimation strategy in this paper and show how we implement this in a computationally efficient way. We then provide two real-data case studies to act as tutorials on how one may use the sparseDFM package. The first case study focuses on summarising the structure of a small subset of quarterly CPI (consumer price inflation) index data for the UK, while the second applies the package onto a large-scale set of monthly time series for the purpose of nowcasting nine of the main trade commodities the UK exports worldwide.}
Dynamic factor models (DFMs) are one of, if not the most, popular tool for summarising the sources of variation among large collections of time series variables. Originally formalised by geweke1977dynamic, it was the application of sargent1977business showing how just two dynamic factors were able to explain the majority of variance in headline US macroeconomic variables that has led to DFMs popularity in modern day forecasting and structural analysis with `big data'. Their natural state-space representation makes it possible to conduct inference based on Kalman filtering and smoothing techniques and so we obtain a compact framework to handle mixed-frequency and missing data issues in economic data sets and the ability to update forecasts in real-time. DFMs have been intensively used in many contexts, including macroeconomic nowcasting giannone2008nowcasting,banbura2010nowcasting,foroni2014comparison, building indicators of economic activity mariano2010coincident, grassi2015euromind, counterfactual analysis harvey1996intervention, luciani2015monetary, and areas away from economics such as psychology molenaar1985dynamic, fisher2015toward, the energy sector wu2013new, lee2016load and many more. See stock2011dynamic and poncela2021factor for detailed surveys of the literature.
The state-space framework of a DFM consists of a measurement equation:
along with a state equation:
for observations $t=1,\dots,n$, which together assume that the common dynamics of a large number of stationary time series variables $\bm{X}_t = (X_{1,t}, \dots, X_{p,t})^\top$ are driven by a small number $r \ll p$ of latent factors $\bm{F}_t = (F_{1,t},\dots,F_{r,t})^\top$ which evolve dynamically as a VAR(1) process over time. The $p \times r$ matrix $\bm{\Lambda}$ contains factor loadings representing the relationship between the measurements and the underlying latent factors that drive their behaviour. The factor loadings play a critical role in the analysis of DFMs, providing an insight into the underlying structure of the data and the strength of correlation between factors and variables. As we will see, factor loadings are the key parameters in the Sparse DFM extension.
The measurement equation ((ref)) can be decomposed into a common component $\bm{\chi}_t = \bm{\Lambda}\bm{F}_t$ and an idiosyncratic component $\bm{\epsilon}_t = (\epsilon_{1,t},\dots,\epsilon_{p,t})^\top$. The common component captures the variation in the time series variables that can be attributed to the underlying latent factors. In a classic DFM, the factors are linear combinations of all variables and therefore $\bm{\chi}_t$ reflects the common sources of variation amongst all variables. The idiosyncratic errors capture the sources of variation that cannot be explained by the factors and are specific to individual series, such as measurement errors. Interpreting the common and idiosyncratic components is important when trying to understand data. For example, in an econometric application, it might be that the common component represents an underlying business cycle, while the idiosyncratic component may represent shocks specific to individual industries or regions kose2008understanding,forni2021policy.
To avoid the proliferation of parameters, it is generally assumed that the idiosyncratic errors are zero-mean and cross-sectionally uncorrelated, meaning their covariance matrix, which we denote $\bm{\Sigma}_{\bm{\epsilon}}$, is diagonal - this is termed an exact DFM. Under this assumption, the correlation of one series with another occurs only through the latent factors $\bm{F}_t$. When this assumption is relaxed, and the idiosyncratic errors are allowed to be (weakly) cross-correlated, termed an approximate DFM, consistent estimation of the factors is possible only in a high-dimensional setting, i.e. when $p \rightarrow \infty$. Therefore, the `curse of dimensionality' problem, often a burden for analysing time series models, is beneficial when estimating approximate DFMs. In the sparseDFM package we consider two distributions for the idiosyncratic errors: $\bm{\epsilon}_t$ is serially uncorrelated IID white noise, or $\bm{\epsilon}_t$ follow an AR(1) process.
Despite the many qualities DFMs possess in handling high-dimensional data, dealing with arbitrary patterns of missing data and being able to update predictions in real time, one may argue they suffer from a weakness - a lack of interpretability. In the classic DFM, every variable in $\bm{X}_t$ is loaded onto the factors $\bm{F}_t$ via the factor loadings, i.e. $\bm{\Lambda}$ is a dense matrix. Often these loading weights are quite large and it is very difficult, if not impossible, to understand how individual series, or groups of series, are influencing the dynamics of the estimated latent factors. If instead, we could obtain factor estimates that are linear combinations of small and distinct groups of variables, where we could associate a given factor estimate with an underlying feature of the observed data, then this would allow us to make statements like: ‘this factor is largely associated with data coming from this group of variables’. This would be highly informative for users of such data.
A natural way we could impose this kind of interpretability into the factor structure would be through the factor loadings $\bm{\Lambda}$. A classic approach to accomplish this is based on factor rotations, which has been widely used to minimise the complexity in the factor loadings to make the structure simpler to interpret. This is useful when you do not have, or would not like to use prior knowledge on factor structure. See kaiser1958varimax for the well-established varimax rotation and see carroll1953analytical and jennrich1966rotation for the well-established quartmin rotation. For a recent paper on varimax rotation see rohe2020vintage. The well-established nowcasting paper of banbura2010nowcasting use prior knowledge of the variables to partition the factors into mutually independent global, real and nominal factors, assuming that the global factor is loaded by all the variables, while real and nominal factors are specific to real and nominal variables. They show this is helpful when calculating a global factor. andreou2019inference is a similar more recent application for grouped factor structure.
An alternative approach is to utilise modern tools from sparse modelling and regularised estimation to constrain the loading matrix to only have a few non-zero loadings for each factor, i.e. is sparse. Sparse principal components analysis (SPCA) from zou2006sparse is a variation of traditional PCA that adds an elastic net penalty (i.e. a mix of $\ell_1$ and $\ell_2$ norms) to induce loadings that are sparse. The interpretable principal components can be viewed as static factors in a factor model and used alongside the Kalman smoother for prediction purposes. See croux2011sparse for a typical macroeconomic forecasting setting where they consider a robustified version. kristensen2017diffusion use SPCA to estimate diffusion indexes with sparse loadings. despois2022identifying prove that SPCA consistently estimates the factors in an approximate factor model if the $\ell_1$ penalty is of $\mathcal{O}(1/\sqrt{p})$. They compare SPCA with factor rotation methods and show an improved performance when the true loadings structure is sparse.
A new paper by mosleySAC implements regularisation (via the $\ell_1$ norm) of factor loadings within an expectation maximisation (EM) algorithm framework allowing the ability to robustly handle arbitrary patterns of missing data, model temporal dependence in the process (i.e. obtain dynamic factors rather than static), and impose weakly informative (sparse) prior knowledge on the factor loadings. Numeric results convey the importance of an EM framework when recovering the true factor loadings and producing accurate forecasts. When the autocorrelation between factors is moderately persistent, they show an improved performance over SPCA for loading recovery and improved performance over a sparse VAR model in forecasting.
The aim of the package SparseDFM is to offer the tools for the R user to implement dynamic factor models with the option to induce sparse loadings. In summary we offer three different popular DFM estimation techniques plus the novel Sparse DFM estimation technique of mosleySAC. We call the estimation options:
We allow the user the option of two different idiosyncratic error processes:
We also allow two different options for estimating the Kalman filter and smoother equations:
Alternative software that implement classic DFMs in R include the MARSS package of holmes2012marss which allow for more general state space structures. The nowcasting packages include nowcasting by de2019nowcasting and nowcastDFM by dhopp2021, that allow mixed-frequency time series nowcasting\footnote{These nowcasting packages have now been removed from the CRAN repository. However, they are still accessible via GitHub.}. Recently, the package dfms by krantz2022 has been uploaded to CRAN which implement the regular DFM methods with IID errors. They make use of C++ code and therefore is computationally faster than the previous three packages whose implementations are solely in R. Our package SparseDFM is novel in multiple ways. We implement the general DFM estimation methods found in the above packages also using C++ code for speed. We implement the univariate treatment of the Kalman filter/smoother equations of koopman2000fast which is computationally much faster than the classic multivariate approach when $p$ is large. We consider AR(1) idiosyncratic errors as well as IID. Finally, we allow the ability to estimate a sparse DFM using an efficient sparsified EM algorithm framework.
This article is structured as follows. Section (ref) provides information on the 4 estimation methods we implement in this package: PCA in Section (ref), 2Stage in Section (ref), EM in Section (ref) and EM-sparse in (ref). We provide derivations of the Kalman filter and smoother equations and parameters estimated in the EM algorithm. In the EM-sparse section we go into detail on methods used to to make the computation of sparse loadings faster. We conclude Section (ref) with information on an extension for AR(1) errors in Section (ref) and how we tune for the number of factors to use in general and the $\ell_1$-norm penalty parameter of EM-sparse in Section (ref). In Section (ref) we outline the structure of the package and functions it includes. Section (ref) provides two real-data case studies on UK inflation and trade in goods to display the package capabilities and visualisation features. Section (ref) concludes the article.
For a $p$-dimensional collection of stationary time series variables $\bm{X}_t$ over $t=1,\dots,n$, we can assume the dynamic factor model for $r$ factors $\bm{F}_t$ to be given by
The key assumption for this to be a sparse DFM is that the $p \times r$ loading parameter $\bm{\Lambda}$ is sparse. The idiosyncratic errors $\{\bm{\epsilon}_t\}$ are assumed to be multivariate white noise and for simplicity we assume $\mathbb{E}[\bm{\epsilon}_t\bm{\epsilon}_t^\top] = \bm{\Sigma}_{\bm{\epsilon}} = \text{diag}(\bm{\sigma}^2_{\bm{\epsilon}})$ and $\bm{\sigma}^2_{\bm{\epsilon}} \in \mathbb{R}^p_+$ is a vector of idiosyncratic variances. In Section (ref) we give an extension to AR(1) idiosyncratic errors. We let $\mathbb{E}[\bm{u}_t\bm{u}_t^\top] = \bm{\Sigma}_{\bm{u}}$ and assume $\|\bm{A}\|_2 <1$, thus the latent VAR(1) model is stationary. This model corresponds to an exact DFM, where all the temporal dependence is modelled via the latent factors.
We now provide the four estimation algorithms the sparseDFM package can implement to obtain estimates of latent factors and model parameters. In the main \verb+sparseDFM()+ function the choice of estimation algorithm can be set using either \verb+alg = "PCA"+, \verb+alg = "2Stage"+, \verb+alg = "EM"+ or \verb+alg = "EM-sparse"+.
The classic way to estimate latent factors is through principal components analysis (PCA), which has been extensively studied in the literature stock2002forecasting, bai2003inferential, doz2020dynamic. For centered and complete\footnote{$\bm{X}_t$ can be made complete by simply removing rows with missing observations or by interpolation. In the package, we interpolate internal missing data using cubic splines and fill in missing data at the start and end of the sample using the median of the series which is then smoothed with an MA(3) process.} $\bm{X}_t$ and assuming the number of $r$ factors is known\footnote{bai2002determining, onatski2009testing and ahn2013eigenvalue, among others, have developed various approaches to consistently estimate the number of factors. In the sparseDFM package we favour the information criteria approach of bai2002determining. For more discussion see Section (ref).}, PCA allows us to simultaneously estimate the factors and their loadings by solving the least squares problem:
subject to the normalization $\bm{\Lambda}^\top \bm{\Lambda}/p = \bm{I}_r$ and $\bm{\Sigma}_F$ is diagonal (where $\bm{\Sigma}_F = \mathbb{E}(\bm{F}_t\bm{F}_t^\top)$). The solution to this least squares problem is the PC estimator of the factors $\hat{\bm{F}}_t = p^{-1}\hat{\bm{\Lambda}}^\top \bm{X}_t$ where the estimated loadings $\hat{\bm{\Lambda}}$ are the matrix of eigenvectors of the sample covariance matrix of $\bm{X}_t$ associated with the largest $r$ eigenvalues.
When mild conditions are placed on the correlation structure of idiosyncratic disturbances, the PCA estimator is the optimal non-parametric estimator for a large approximate DFM. With even tighter conditions of spherical idiosyncratic components (i.i.d. Gaussian), the PCA estimator is equivalent to the maximum likelihood estimator doz2020dynamic. The problem with using non-parametric PCA methods to estimate DFMs is that there is no consideration of the dynamics of the factors or idiosyncratic components. In particular, there is no feedback from the estimation of the factor equation ((ref)) to the observable equation ((ref)). For this reason, it is preferable to use parametric methods that are able to model the dynamics in the system.
giannone2008nowcasting proposed a two-stage (2Stage) framework for estimating DFMs which has since been theoretically studied by doz2011two and successfully applied to the field of nowcasting. This method involves the following two steps:
The Kalman filter (KF) enables the factor estimates to be continually updated as new observations become available. This way, we are able to assess the impact new data releases of certain variables have in forecasting models. In the nowcasting literature this is referred to as the `news' of data releases. The KF provides the basis for smoothing, where the factor at a time $t$ can be estimated based on all data in the sample. This is very helpful for missing data interpolation and backcasting.
Appendix (ref) provides the classic multivariate Kalman filter and smoother (KFS) equations of shumway1982approach, often used in DFM estimation. We also provide the univariate treatment (sequential processing) of the multivariate equations for fast Kalman filter and smoother seen in koopman2000fast that can lead to substantial computational gains for exact DFMs. As it is assumed in exact DFMs that $\bm{\Sigma}_{\bm{\epsilon}}$ is diagonal, it becomes possible to filter the observations $\bm{X}_t$ one element at a time, as opposed to altogether as in the classic multivariate approach. This way, matrix inversions become scalar divisions and thus huge speedups are possible. In the main \verb+sparseDFM()+ function it is possible to set \verb+kalman = "multivariate"+ or \verb+kalman = "univariate"+ for the respective KFS methods. The default is \verb+kalman = "univariate"+ for fast solutions when the errors are assumed i.i.d..
banbura2014maximum build on the DFM representation of watson1983alternative and adopt an expectation-maximisation (EM) approach to estimate the system ((ref))-((ref)) by quasi-maximum likelihood estimation (QMLE)\footnote{It is described as quasi-maximum likelihood as consistency results are robust to cross-sectional misspecification, serial correlation of idiosyncratic components, and non-Gaussianity.}. Assuming residuals in the system ((ref))-((ref)) are Gaussian, and collecting all parameters in $\bm{\theta} = (\bm{\Lambda}, \bm{A}, \bm{\Sigma}_{\bm{\epsilon}}, \bm{\Sigma}_{\bm{u}})$, the idea is to write the joint log-likelihood of $\bm{X}_t$ and $\bm{F}_t$ for all $t=1,\dots,n$ in terms of both the observed and unobserved data as:
where $\bm{\epsilon}_t = \bm{X}_t-\bm{\Lambda} \bm{F}_t$, $\bm{u}_t=\bm{F}_t-\bm{A}\bm{F}_{t-1}$, and we have assumed an initial distribution at $t=0$ of the factors as $\bm{F}_0 \sim N(\bm{\alpha}_0, \mathcal{\bm{P}}_0)$. We then iterate between the two steps:
until convergence\footnote{We choose to use the same EM convergence criteria as seen in doz2012quasi, which is determined when:
where the user has on option in the sparseDFM() function to set the threshold value. The default is $10^{-4}$.}.
The first iteration of the EM algorithm requires initial parameter values $\bm{\theta}^{(0)}$. These are obtained by doing stage one of the 2Stage framework of giannone2008nowcasting as in Section (ref). The initial factor mean $\bm{\alpha}_0$ is set to 0 as in banbura2014maximum. The initial factor covariance $\mathcal{\bm{P}}_0$ is set to the matrix of $\mathrm{vec}(\mathcal{\bm{P}}_0) = (\bm{I}_{rp} - \hat{\bm{A}} \otimes \hat{\bm{A}})^{-1}\mathrm{vec}(\hat{\bm{\Sigma}}_{\bm{u}})$. This initialisation will likely to be good when $p$ is large, as PCA estimation of factor loadings have been shown to be consistent for large cross-sections bai2003inferential, doz2011two.
The EM algorithm then proceeds iterating between the E step and M step. Assuming Gaussian errors, the conditional expectation in the E step can be found using the Kalman filter and smoother equations of the state-space ((ref))-((ref)) as given in Appendix (ref). From this we obtain the conditional mean and covariances of the state which we denote by:
conditional on all information we have observed up to $n$, denoted by $\bm{\Omega}_n$.
As shown in banbura2014maximum and fully derived in Appendix (ref), the maximisation of the conditional expected log-likelihood (i.e. the M step) results in the following expressions for the parameter estimates:
and letting $\bm{S}_{t\vert n} = \bm{a}_{t\vert n}\bm{a}_{t\vert n}^{\top} + \bm{P}_{t\vert n}$, and $\bm{S}_{t,t-1\vert n} = \bm{a}_{t\vert n}\bm{a}_{t-1\vert n}^{\top} + \bm{P}_{t,t-1\vert n}$ we have
For parameters $\bm{\Lambda}$ and $\bm{\Sigma}_{\bm{\epsilon}}$, we should also consider there might be missing data in $\bm{X}_t$. banbura2014maximum define a selection matrix $\bm{W}_t$ to be a diagonal matrix such that
and note that $\bm{X}_t = \bm{W}_t\bm{X}_t + (\bm{I}-\bm{W}_t)\bm{X}_t $. The update for the loadings and the idiosyncratic error covariance is then given by
where $\mathrm{vec}(\cdot)$ represents vectorisation of a matrix and $\hat{\bm{\Sigma}}_{\bm{\epsilon}}^*$ is obtained from the previous EM iteration. In practice we update $\hat{\bm{\Sigma}}_{\bm{\epsilon}}$ after estimating $\hat{\bm{\Lambda}}$, as the former is based on the difference between the observations and the estimated common component. See Appendix (ref) for full detail of the EM algorithm.
doz2012quasi and barigozzi2019quasi provide consistency results of the EM approach as $(n,p) \rightarrow \infty$. In terms of macroeconomic forecasting, the EM algorithm approach has many gains over the PCA and 2Stage procedures. It is much more efficient in small samples, as shown by simulation studies of doz2012quasi. It is a very flexible framework able to deal with arbitrary patterns of missing data in $\bm{X}_t$, as shown in banbura2014maximum. It enables one to impose restrictions on parameters when updating them in the EM algorithm. For example, banbura2010nowcasting impose restrictions on the loadings to reflect temporal aggregation and introduce factor specific groups of variables in their nowcasting application. Furthermore, it is a computationally efficient algorithm allowing fast forecasts for policy makers/analysts, as opposed to computationally demanding MCMC based alternatives.
mosleySAC extend the banbura2014maximum EM algorithm framework by imposing LASSO regularisation on the conditional expected log-likelihood for the $\bm{\Lambda}$ update with the endeavour of obtaining sparse factor loadings. In particular, they find estimators of the form
where $\alpha$ is the LASSO tuning parameter controlling the amount of sparsity in $\bm{\Lambda}$. Due to the complexity of the objective being maximised in ((ref)), mosleySAC propose to solve the LASSO problem using the Alternative Directed Method of Moments (ADMM) by boyd2011distributed as this approach provides a way to split up the complex objective function into separate, simpler optimisation problems. This involves considering the penalised and augmented Lagrangian
where $\bm{Z}\in\mathbb{R}^{p\times r}$ is an auxiliary variable, $\bm{U}\in\mathbb{R}^{p\times r}$ are the (scaled) Lagrange multipliers and $\nu$ is the scaling term. Under equality conditions relating the auxilary ($\bm{Z}$) to the primal ($\bm{\Lambda}$) variables, this is equivalent to minimising ((ref)), e.g.
as ((ref)) is convex in the argument $\bm{\Lambda}$ with all other parameters fixed, this argument holds for any $\nu>0$.
The augmented Lagrangian ((ref)) can be sequentially minimised via the following updates\footnote{The full derivation of these equations are given in Appendix (ref). }
for $k=0,1,2,\ldots,$ until convergence. The first (primal) update is simply a least-squares type problem, whereby on vectorising $\bm{\Lambda}$ one finds
mosleySAC further detail how the computation of ((ref)) can be made significantly faster by exploiting dimensionality reduction, going from a per-iteration cost of the ADMM of $\mathcal{O}(r^3p^3)$ to $\mathcal{O}(r^3p)$. For details of this speed up refer to mosleySAC. We implement this speed up in the package.
To extend the state-space framework ((ref))-((ref)) to AR(1) idiosyncratic errors, as opposed to i.i.d. white noise, we just follow the methodology laid out in banbura2014maximum and augment the state vector so it contains the latent factors and the latent AR(1) errors. Formally, we set up the AR(1) process like so
for $i=1,\dots,p$ and $t=1,\dots,n$, where $\kappa$ is a very small number which allows us to easily set up the new state-space framework.
This new state-space framework is given by
where $$ \tilde{\bm{F}}_t=\left[
\right] \, , \; \tilde{\bm{u}}_t=\left[
\right] \, , \; \tilde{\bm{\Lambda}}=\left[
\right] \, , \; \tilde{\bm{A}}=\left[
\right] \, , $$ $$ \bm{\Sigma}_{\tilde{\bm{u}}}=\left[
\right] \; and \quad \bm{\Sigma}_{\bm{\eta}} = \kappa \bm{I}_{r+p} \, , $$ for $\bm{\phi} = diag(\phi_1, \dots, \phi_p)$ and $\bm{\Sigma}_{\bm{e}} = \mathbb{E}[\bm{e}_t\bm{e}_t^\top] = diag(\bm{\sigma}_{\bm{e}}^2)$ and $\bm{\sigma}_{\bm{e}}^2 \in \mathbb{R}^p_{+}$ is a vector of error variances.
The parameter estimates of this new augmented state-space model can be obtained very similarly to before. See banbura2014maximum for the estimates of the new set of parameters $\tilde{\bm{\theta}}=(\bm{A}, \bm{\phi}, \bm{\Lambda}, \bm{\Sigma}_{\bm{u}}, \bm{\Sigma}_{\bm{e}})$.
To calculate the number of factors to use in the model we opt to take the information criteria approach of bai2002determining. This can be done before the \verb+SparseDFM()+ model is fitted. bai2002determining consider 3 information criteria with different penalties of the form:
The sum of squared residuals for $r$ factors $V_r(\hat{\bm{F}},\hat{\bm{\Lambda}}) = \sum_{i=1}^p\sum_{t=1}^n\mathbb{E}[\hat{\epsilon}_{i,t}^2]\big/{np}$ with $\hat{\epsilon}_{i,t} = X_{t,i}-\hat{\bm{F}}_t\hat{\bm{\Lambda}}_i$ is found using PCA on the standardised dataset $\bm{X}$. The estimated factors $\hat{\bm{F}}$ corresponding to the principal components and the estimated loadings $\hat{\bm{\Lambda}}$ corresponding to the eigenvectors. Should the data contain missing values, we first interpolate the missing values using the median of the series that are then smoothed by a simple moving window as discussed in Section (ref).
The number of factors to use will correspond to $\operatorname*{arg\,min}_r IC_i(r)$ for $i=1,2$ or $3$. The second penalty is the highest when working in finite samples and therefore is set to default in our \verb+tuneFactors()+ function.
We tune the model for the best $\alpha$ by doing a grid search over a range of $\alpha$'s using BIC. For factors $\hat{\bm{F}}$ and $\hat{\bm{\Lambda}}$ estimated by a run of the EM algorithm for a particular $\alpha$ value, $\alpha^*$, the BIC is
where $V_{\alpha^*}(\hat{\bm{F}},\hat{\bm{\Lambda}})$ is the sum of squared residuals using $\alpha^*$ and $m$ is the number of non-zero entries in $\hat{\bm{\Lambda}}$.
The default search grid for $\alpha$ is a log-space from $10^{-2}$ to $10^3$. Any grid can be inputted by the user and we make sure the grid is ordered from low to high. Each time the EM algorithm is run with a new $\alpha$ value, the estimated parameters from the previous $\alpha$ value are used as a warm start to the EM algorithm. This should reduce the number of iterations needed for the EM algorithm to converge as we search through the $\alpha$ grid. We implement a stopping rule in the grid search if a column of the loadings parameter $\bm{\Lambda}$ gets set entirely to zero.
The sparseDFM package can be installed from CRAN as follows:
The development version can be installed from GitHub using the devtools package devtools:
To load the package into the R environment, type:
The core function of the package is \verb+sparseDFM+. This function allows different estimation methods for dynamic factor models including: principal components analysis (PCA) for static factors stock2002forecasting, the two-stage (2Stage) framework of PCA plus Kalman filter/smoother giannone2008nowcasting, the quasi-maximum likelihood approach using the expectation-maximisation (EM) algorithm to handle arbitrary patterns of missing data banbura2014maximum, and the novel sparsified EM approach (EM-sparse) allowing LASSO regularisation on factor loadings mosleySAC. The function can be implemented using:
It takes the following arguments:
The function \verb+sparseDFM+ returns an S3 object of class \verb+`sparseDFM'+. This allows the user to pass the generic functions: \verb+print()+, \verb+summary()+, \verb+plot()+, \verb+residuals()+, \verb+resids()+, \verb+fitted()+, and \verb+predict()+ to the returned \verb+sparseDFM+ object to display generic outputs. The \verb+predict()+ function will itself return an S3 object of class \verb+`sparseDFM_forecast'+ to allow the user to \verb+print()+ forecast results. There are plenty of helpful types of plot included in the package to allow the user to visualise estimation results (see Appendix (ref)). We go into detail on the specific usage of these options in the empirical applications of Section (ref).
When called, the S3 object \verb+sparseDFM+ returns a list-of-lists containing the following:
This package also includes eight other exported functions:
In this section we provide a tutorial of how one might use the package for empirical applications to real data. We provide two case studies both with the intention to display the package capabilities and plotting features. The first case study focusses on a small subset of quarterly CPI (consumer price inflation) index data for the UK taken from the Office from National Statistics' (ONS) website\footnote{We use the Q4 2022 release (benchmarked to 2015=100) from \url{https://www.ons.gov.uk/economy/inflationandpriceindices/datasets/consumerpriceindices}}. The data contains 36 variables of different classes of the inflation index and 135 observations from 1989 Q1 to 2022 Q3. The purpose of this small, 36 variable example is to demonstrate the core functionality and the ability to graphically display results with the sparseDFM package. The second case study applies the package onto a large-scale data set for the purpose of nowcasting UK trade in goods. The data contains 445 columns, including 9 target series (UK exports of the 9 main commodities worldwide) and 434 monthly indicator series, and 226 rows representing monthly values from January 2004 to October 2022. In this large scale example we make comparisons between a regular DFM and a Sparse DFM in terms of interpreting factor structure and accuracy of predictions.
Before going through the case studies, we refer the reader to the table in Appendix (ref) that summarises the various types of plot and other S3 generic functions that the package contains.
The inflation data is included in the sparseDFM package and we can load it into the R environment as a data frame called \verb+data+ using:
Before we fit a model it is worthwhile to first perform some exploratory data analysis. Two main properties to look out for when working with dynamic factor models are the amount of missing data present and if the data series are stationary or not. The function \verb+missing_data_plot()+ allows the user to visualise where missing data is present:
Figure (ref) shows the output of this missing data plot. We can see that $6$ of the $36$ variables have missing data towards the start of the sample.
In a DFM, we assume the latent factors to follow a stationary VAR(1) process and hence we require the input data to be stationary. After exploration it appears the inflation data is stationary in first-differences. Hence, we can apply the function \verb+transformData()+ and set the \verb+stationary_transform+ parameter to represent a first-difference transform for all variables:
We store the stationary first-difference data as the variable \verb+new_data+. The \verb+stationary_transform+ parameter is a vector of length matching the number of variables in the data that contains a stationary transform index for each variable with the options: \verb+1+ (no change), \verb+2+ (first-difference), \verb+3+ (second-difference), \verb+4+ (log-first-difference), \verb+5+ (log-second-difference), \verb+6+ (growth rate) or \verb+7+ (log growth rate). The returned transformed data set is of the same dimension as the original data with \verb+NA+ where missing data is, for example, the first row of \verb+new_data+ will all be \verb+NA+ as we do not have an observed 1988 Q4 value to difference with.
Now we have a stationary data set, the next step would be to determine the number of factors to use with our DFM. To do this we can apply the function \verb+tuneFactors()+ that selects the best number of factors to use based on bai2002determining information criteria:
The default information criteria (IC) is \verb+type = 2+ (see Section (ref) for details). We search for a maximum number of 15 factors and plot the results. Figure (ref) (in Appendix (ref)) shows the output of \verb+tuneFactors()+: the IC values for different number of factors along with a screeplot showing the percentage variance explained obtained from the eigenvalues of the covariance of the data.
It appears 3 factors is sufficient for the inflation exercise. We can now fit a Sparse DFM using 3 factors and call the generic S3 function \verb+summary()+ to summarise the call:
We choose to tune the $\ell_1$-norm penalty parameter using BIC over the default logspace grid of $100$ values between $10^{-2}$ and $10^3$. The algorithm will stop searching over this grid as soon as an entire column of $\bm{\Lambda}$ becomes $0$. We can find the optimal $\ell_1$-norm parameter $\hat{\alpha}$ and plot the BIC curve (shown in Figure (ref) in Appendix (ref)) using:
We are able to display the estimation results using various types of plots - see Appendix (ref) for options. We can plot the three estimated factors overlaid on top of the standardised inflation data using
and make a boxplot of the residuals for each inflation series using
See Figure (ref) and (ref) in Appendix (ref) for these plots, respectively. Note, we can set \verb+use.series.names = FALSE+ to avoid plotting the series names. A very useful plot to understand the factor loading structure is to use:
This produces a heatmap (shown in Figure (ref)) of the estimated factor loadings with a clear display of sparse loadings displayed by blank squares.
Nowcasting\footnote{For a detailed survey on the nowcasting literature see banbura2013now} is a method used in econometrics that involves estimating the current state of the economy based on the most recent data available. It is an important tool because it allows policy makers and businesses to make more informed decisions in real-time, rather than relying on outdated information due to publication delays that may no longer be accurate. Trade in Goods (imports and exports) is currently published with a 2 month lag by the UK's Office for National Statistics (ONS), which is quite a long time to wait for current assessments of trade, especially during times of economic uncertainty or instability. Nowcasting UK trade information has become particularly important in recent years due to the key events of the Brexit referendum, held in 2016, and the coronavirus pandemic, reaching UK shores early 2020. While the cause of these shocks are drastically different, both have imposed restrictions on trade in goods.
Using sparseDFM, we consider the task of understanding and nowcasting the movements of 9 monthly target series representing 9 of the main commodities the UK exports worldwide. These include: Food & live animals, Beverages and tobacco, Crude materials, Fuels, Animal and vegetable oils and fats, \textbf{Chemicals}, \textbf{Material manufactures}, \textbf{Machinery and transport} and \textbf{Miscellaneous manufactures}.
To try and produce accurate nowcasts for the targets we use a large collection of 434 monthly indicator series including:
The sparseDFM package contains a data frame \verb+exports+ containing the 9 monthly target series of UK exports and the 434 monthly indicator series observed from January 2004 to October 2022\footnote{All data is openly available from the Office for National Statistics' website: \url{https://www.ons.gov.uk/}.} and can be loaded into the R environment using
Figure (ref) provides a plot of the 9 target series we are interested in nowcasting. We see exports of machinery and transport being the largest. We also see two main drops during the 2009 and 2020 recessions and an upwards trend in the past year or so. This figure can be generated using:
To make the data stationary, we apply first-difference transformations to each series using:
The only missing data present in the data is at the end of the sample during the months of September and October 2022 depending on publication delays of the variables. This forms a ragged edge at the end of the sample with both months missing for the targets and IoP, just October missing for CPI, PPI, exchange rates, BCI and CCI, and nothing missing for google trends. This can be visualised with Figure (ref) in Appendix (ref) which plots a \verb+missing_data_plot()+ for the last 12 months of the sample using:
We are interested in nowcasting both September and October for the export targets. We tune for the number of factors using (see Figure (ref) in Appendix (ref) for the output)
According to the bai2002determining information criteria, the best number of factors to use is 7. However, the screeplot seems to suggest that after 4 factors, the addition of more factors does not add that much in terms of explaining the variance of the data. For this reason, we choose to use 4 factors when modelling.
We now fit a regular DFM and a Sparse DFM to the data with 4 factors:
Note, we set the parameter \verb+q = 9+ in the Sparse DFM fit as we do not want to regularise the first 9 variables of the data set corresponding to the 9 target series.
We can explore the convergence and tuning of each algorithm like so:
We now explore the estimated factors and loadings for the regular and sparse DFM fits. We are able to group the indicator series into colours depending on the source of the indicator and use the \verb+type = "loading.grouplineplot"+ setting in \verb+plot()+. We set the trade in goods (TiG) target black, IoP blue, CPI red, PPI pink, exchange rate (Exch) green, BCI & CCI (Conf) navy and google trends (GT) brown. This will make it easier to visualise which indicators are loading onto specific factors. Factor and loading plots can be generated by:
Figure (ref) displays the generated factor and loadings plots for the regular DFM fit on the top row and the Sparse DFM fit on the bottom row. We have used the R package gridExtra to neatly display the loading plots in a $2\times 2$ grid. As all variables are loaded onto all the factors in a regular DFM it is very difficult to interpret the factor structure from these loading plots. The loadings from all groups in every factor are quite large and it is impossible to make conclusions on which data groups are related to specific factors. On the other hand, with sparse factor loadings, we are able to make a lot clearer conclusions on factor structure. We can clearly visualise which indicator series are the driving force between each factor. Factor 2 for example, with the obvious drop in early 2020 due to the COVID pandemic, is heavily loaded with indicators coming from the Index of Production, confidence indices and google trends. Index of Production does not actually appear in any other factors and we can view this as a clear indicator of the COVID drop. It is interesting that google search words to do with trade in goods are present in factor 2 as well. With lots of economic volatility in recent years, using google trends search words may be a very useful indicator of economic activity. Factor 1 seems to be mainly loaded with PPI data, while factor 3 is heavily loaded with CPI data. Some inflation data and exchange rate indicators are present in factor 4, which seems to have shocks during the 2009 and 2020 recessions. Note that in all 4 factor loading plots the trade in goods target series are loaded as we specified q = 9 in the sparseDFM fit to ensure these variables are not regularised.
We now turn our attention to nowcasts of the missing months of September and October 2022. It is very easy to extract nowcasts from the \verb+sparseDFM+ fit. As the data was inputted with the ragged edge structure, with NAs coded in for September and October 2022 for the target series, the sparseDFM output will provide us with estimates for these missing months. There are two ways these values can be extracted:
To transform these first-differenced-nowcasts into nowcasts on the original level data, we need to undifference. To do this we can take the most recent observed value (August 2022) and add the September first-difference-nowcast for the September level nowcast, and then add the October first-difference-nowcast to this value to get the October level nowcast:
To determine which model is producing the most accurate nowcasts we can set up a pseudo real-time nowcasting exercise where we try and predict the missing ragged edge in an expanding window from August 2018 to July 2022. As it would be difficult to store all previous vintages of data back to 2018, we just use the complete data set (up to October 2022) and at each nowcasting window create a ragged edge following the 2 month lag for the targets and IoP and 1 month lag for everything else except google trends. We begin with the window Jan 2004 - Aug 2018, where the current available target exports data is June 2018. We produce nowcasts for July 2018 (representing a nowcast 1 month prior to release) and for August 2018 (representing a nowcast 2 months prior to release). We produce these nowcasts for horizon 1 and 2 for every month up to July 2022 - a total of 48 months. The code for the pseudo real-time nowcasting exercise is presented in Appendix (ref). Note, we fit a new model at each nowcasting window and re-tune for the best $\ell_1$-norm penalty parameter. An alternative could be to store the model parameters estimated from one fit of the SDFM to a particular time window and then just run the \verb+kalmanUnivariate()+ function on these stored parameters at each time window.
Table (ref) displays the average, 25th, 50th and 75th percentiles of the mean absolute error across the 48 month expanding window for horizon 1 and 2 in both models. It appears that Sparse DFM is performing better than a regular DFM. It has a lower average mean absolute error and tighter bands around the median. As expected, the error for horizon 1 is slightly lower than horizon 2 as it is able to exploit all indicators with a 1 month lag in its estimation.
The R package sparseDFM has been developed to allow the fitting of various estimation methods for dynamic factor models (DFMs) including PCA stock2002forecasting, 2Stage giannone2008nowcasting, EM banbura2014maximum, and the novel EM-sparse approach for sparse DFMs mosleySAC. We implement both multivariate shumway1982approach, and univariate koopman2000fast Kalman filter and smoother equations in C++ using the Armadillo library rcpparmadillo, as well as the ADMM updates for the EM-sparse algorithm for computational speed.
The implementation of an efficient and interpretable EM framework for sparse DFM estimation is new to the R programming language and is available from the Comprehensive R Archive Network (CRAN) at \url{https://cran.r-project.org/}. The real-data case studies on UK inflation and nowcasting Trade in Goods displayed the benefits of a sparse DFM; from understanding factor structure to nowcasting performance. The sparseDFM package will be maintained at \url{https://github.com/mosleyl/sparseDFM} and we hope add new features in the future including new plotting capabilities and extended assumptions on the underlying VAR(1) latent factor process.
AG and TTC acknowledge funding from the EPSRC grant EP/T025964/1. AG and LM were supported by the ESRC grant ES/V006339/1. LM acknowledges support from the STOR-i Centre for Doctoral Training and the Office for National Statistics.