# sparseDFM: An R Package to Estimate Dynamic Factor Models with Sparse
Loadings

Luke Mosley  Department of Mathematics and Statistics, Lancaster University,
Lancaster, LA1 4YW, United Kingdom  Tak-Shing T. Chan  Department of
Mathematics and Statistics, Lancaster University, Lancaster, LA1 4YW, United
Kingdom  Alex Gibberd111Corresponding author (e-mail:
<EMAIL_ADDRESS>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
Mosley et al. (2023). 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). Mosley et al. (2023) 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.

##  1 Introduction

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 Geweke (1977), it was the application of
Sargent et al. (1977) 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 (Giannone et al., 2008; Banbura et al., 2010; Foroni and
Marcellino, 2014), building indicators of economic activity (Mariano and
Murasawa, 2010; Grassi et al., 2015), counterfactual analysis (Harvey, 1996;
Luciani, 2015), and areas away from economics such as psychology (Molenaar,
1985; Fisher, 2015), the energy sector (Wu et al., 2013; Lee and Baldick,
2016) and many more. See Stock and Watson (2011) and Poncela et al. (2021) for
detailed surveys of the literature.

The state-space framework of a DFM consists of a measurement equation:

$\bm{X}_{t}=\bm{\Lambda}\bm{F}_{t}+\bm{\epsilon}_{t}\,,$ (1)

along with a state equation:

$\bm{F}_{t}=\bm{A}\bm{F}_{t-1}+\bm{u}_{t}\,,$ (2)

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 (1) 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 (Kose et al., 2008; Forni and Gambetti, 2021).

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 Kaiser (1958) for the
well-established varimax rotation and see Carroll (1953) and Jennrich and
Sampson (1966) for the well-established quartmin rotation. For a recent paper
on varimax rotation see Rohe and Zeng (2020). The well-established nowcasting
paper of Banbura et al. (2010) 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. Andreou et al. (2019) 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 Zou et al. (2006) 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 Croux and Exterkate (2011) for a typical
macroeconomic forecasting setting where they consider a robustified version.
Kristensen (2017) use SPCA to estimate diffusion indexes with sparse loadings.
Despois and Doz (2022) 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 Mosley et al. (2023) 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 Mosley et al. (2023). We call the
estimation options:

  * •

`PCA` \- principal components analysis (PCA) for static factors seen in Stock
and Watson (2002).

  * •

`2Stage` \- the two-stage framework of PCA plus Kalman filter & smoother seen
in Giannone et al. (2008) and Doz et al. (2011).

  * •

`EM` \- the quasi-maximum likelihood approach using the EM algorithm to handle
arbitrary patterns of missing data seen in Bańbura and Modugno (2014).

  * •

`EM-sparse` \- the novel sparse EM approach allowing for sparse factor
loadings seen in Mosley et al. (2023).

We allow the user the option of two different idiosyncratic error processes:

  * •

`IID` \- errors are IID white noise: $\epsilon_{i,t}\sim N(0,\sigma_{i}^{2})$.

  * •

`AR1` \- errors follow an AR(1) process:
$\epsilon_{i,t}=\phi_{i}\epsilon_{i,t-1}+e_{i,t}$ with $e_{i,t}\sim
N(0,\sigma_{i}^{2})$.

We also allow two different options for estimating the Kalman filter and
smoother equations:

  * •

`multivariate` \- classic Kalman filter and smoother equations seen in Shumway
and Stoffer (1982).

  * •

`univariate` \- univariate treatment (sequential processing) of the
multivariate equations for fast Kalman filter and smoother seen in Koopman and
Durbin (2000).

Alternative software that implement classic DFMs in R include the MARSS
package of Holmes et al. (2012) which allow for more general state space
structures. The nowcasting packages include nowcasting by De Valk et al.
(2019) and nowcastDFM by Hopp (2021), that allow mixed-frequency time series
nowcasting222These nowcasting packages have now been removed from the CRAN
repository. However, they are still accessible via GitHub.. Recently, the
package dfms by Krantz and Bagdziunas (2022) 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 Koopman and Durbin (2000) 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 2 provides information on the 4
estimation methods we implement in this package: PCA in Section 2.1, 2Stage in
Section 2.2, EM in Section 2.3 and EM-sparse in 2.4. 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 2 with
information on an extension for AR(1) errors in Section 2.5 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 2.6. In Section 3 we outline the structure
of the package and functions it includes. Section 4 provides two real-data
case studies on UK inflation and trade in goods to display the package
capabilities and visualisation features. Section 5 concludes the article.

##  2 Estimation Methods and Tuning

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

$\displaystyle\bm{X}_{t}$
$\displaystyle=\bm{\Lambda}\bm{F}_{t}+\bm{\epsilon}_{t}\,,$ (3)
$\displaystyle\bm{F}_{t}$ $\displaystyle=\bm{A}\bm{F}_{t-1}+\bm{u}_{t}\,.$ (4)

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 2.5 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 `sparseDFM()` function the choice of estimation algorithm can be set
using either `alg = "PCA"`, `alg = "2Stage"`, `alg = "EM"` or `alg = "EM-
sparse"`.

###  2.1 PCA - Stock and Watson (2002)

The classic way to estimate latent factors is through principal components
analysis (PCA), which has been extensively studied in the literature (Stock
and Watson, 2002; Bai, 2003; Doz and Fuleky, 2020). For centered and
complete333$\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 known444Bai and Ng (2002), Onatski (2009) and Ahn and Horenstein
(2013), among others, have developed various approaches to consistently
estimate the number of factors. In the sparseDFM package we favour the
information criteria approach of Bai and Ng (2002). For more discussion see
Section 2.5., PCA allows us to simultaneously estimate the factors and their
loadings by solving the least squares problem:

$\min_{\bm{\Lambda},\bm{F}}\frac{1}{np}\sum_{t=1}^{n}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})=\min_{\bm{\Lambda},\bm{F}}\frac{1}{np}\sum_{t=1}^{n}\bm{\epsilon}_{t}^{\top}\bm{\epsilon}_{t}\,,$
(5)

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 (Doz and Fuleky, 2020). 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 (4) to the
observable equation (3). For this reason, it is preferable to use parametric
methods that are able to model the dynamics in the system.

###  2.2 2Stage - Giannone et al. (2008)

Giannone et al. (2008) proposed a two-stage (2Stage) framework for estimating
DFMs which has since been theoretically studied by Doz et al. (2011) and
successfully applied to the field of nowcasting. This method involves the
following two steps:

  * Step 1:

Preliminary estimates of the loadings $\bm{\Lambda}$ and the factors
$\bm{F}_{t}$ are found via PCA as in (5). $\bm{\Sigma}_{e}$ is estimated as
the empirical covariance of
$\hat{\bm{\epsilon}}_{t}=\bm{X}_{t}-\hat{\bm{\Lambda}}^{\top}\hat{\bm{F}}_{t}$.
The remaining parameters are estimated by fitting a VAR(1) model on
$\hat{\bm{F}}_{t}$.

  * Step 2:

The model is cast into state space form (3)-(4) and the factors are re-
estimated using the Kalman filter and smoother (Shumway and Stoffer, 1982;
Koopman and Durbin, 2000).

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 A provides the classic multivariate Kalman filter and smoother (KFS)
equations of Shumway and Stoffer (1982), 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 Koopman and Durbin
(2000) 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 `sparseDFM()` function it is possible to set `kalman =
"multivariate"` or `kalman = "univariate"` for the respective KFS methods. The
default is `kalman = "univariate"` for fast solutions when the errors are
assumed i.i.d..

###  2.3 EM - Bańbura and Modugno (2014)

Bańbura and Modugno (2014) build on the DFM representation of Watson and Engle
(1983) and adopt an expectation-maximisation (EM) approach to estimate the
system (3)-(4) by quasi-maximum likelihood estimation (QMLE)555It 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 (3)-(4) 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:

$\displaystyle\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})=$
$\displaystyle-\frac{1}{2}\log|\mathcal{\bm{P}}_{0}|-\frac{1}{2}(\bm{F}_{0}-\bm{\alpha}_{0})^{\top}\mathcal{\bm{P}}_{0}^{-1}(\bm{F}_{0}-\bm{\alpha}_{0})$
$\displaystyle-\frac{n}{2}\log|\bm{\Sigma}_{\bm{u}}|-\frac{1}{2}\sum_{t=1}^{n}\bm{u}_{t}^{\top}\bm{\Sigma}_{\bm{u}}^{-1}\bm{u}_{t}$
$\displaystyle-\frac{n}{2}\log|\bm{\Sigma}_{\bm{\epsilon}}|-\frac{1}{2}\sum_{t=1}^{n}\bm{\epsilon}_{t}^{\top}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{\epsilon}_{t}\,,$
(6)

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:

  * E Step:

Compute the expectation of the joint log-likelihood (2.3) conditional on the
available information up to $n$, which we denote by $\bm{\Omega}_{n}$, using
the parameters estimated in the previous iteration, $j$:

$\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta}^{(j)})|\bm{\Omega}_{n}\right]\,;$

  * M Step:

Re-estimate the parameters through the maximisation of the conditional
expected joint log-likelihood:

$\bm{\theta}^{(j+1)}=\operatorname*{arg\,max}_{\bm{\theta}}\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta}^{(j)})|\bm{\Omega}_{n}\right]\,,$

until convergence666We choose to use the same EM convergence criteria as seen
in Doz et al. (2012), which is determined when:
$M_{j}=\frac{\log\mathcal{L}(\bm{X};\hat{\bm{\theta}}^{(j)})-\log\mathcal{L}(\bm{X};\hat{\bm{\theta}}^{(j-1)})}{\left(\log\mathcal{L}(\bm{X};\hat{\bm{\theta}}^{(j)})+\log\mathcal{L}(\bm{X};\hat{\bm{\theta}}^{(j-1)})\right)\Big{/}2}<\text{threshold}\,,$
(7) 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 Giannone et al. (2008) as in Section 2.2. The initial factor mean
$\bm{\alpha}_{0}$ is set to 0 as in Bańbura and Modugno (2014). 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 (Bai, 2003; Doz et al., 2011).

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
(3)-(4) as given in Appendix A. From this we obtain the conditional mean and
covariances of the state which we denote by:

$\displaystyle\bm{a}_{t|n}$
$\displaystyle=\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\,,$
$\displaystyle\bm{P}_{t|n}$
$\displaystyle=\mathrm{Var}[\bm{F}_{t}|\bm{\Omega}_{n}]\,,$
$\displaystyle\bm{P}_{t,t-1|n}$
$\displaystyle=\mathrm{Cov}[\bm{F}_{t},\bm{F}_{t-1}|\bm{\Omega}_{n}]\,,$

conditional on all information we have observed up to $n$, denoted by
$\bm{\Omega}_{n}$.

As shown in Bańbura and Modugno (2014) and fully derived in Appendix B, the
maximisation of the conditional expected log-likelihood (i.e. the M step)
results in the following expressions for the parameter estimates:

$\hat{\bm{\alpha}}_{0}=\bm{a}_{0|n}\quad;\quad\hat{\mathcal{\bm{P}}}_{0}=\bm{P}_{0|n}$

and letting $\bm{S}_{t|n}=\bm{a}_{t|n}\bm{a}_{t|n}^{\top}+\bm{P}_{t|n}$, and
$\bm{S}_{t,t-1|n}=\bm{a}_{t|n}\bm{a}_{t-1|n}^{\top}+\bm{P}_{t,t-1|n}$ we have

$\displaystyle\hat{\bm{A}}$
$\displaystyle=\left(\sum_{t=1}^{n}\bm{S}_{t-1|n}\right)^{-1}\left(\sum_{t=1}^{n}\bm{S}_{t,t-1|n}\right)\,,$
$\displaystyle\hat{\bm{\Sigma}}_{\bm{u}}$
$\displaystyle=\frac{1}{n}\sum_{t=1}^{n}\left[\bm{S}_{t|n}-\hat{\bm{A}}\left(\bm{S}_{t-1,t|n}\right)\right]\,.$

For parameters $\bm{\Lambda}$ and $\bm{\Sigma}_{\bm{\epsilon}}$, we should
also consider there might be missing data in $\bm{X}_{t}$. Bańbura and Modugno
(2014) define a selection matrix $\bm{W}_{t}$ to be a diagonal matrix such
that

$W_{t,ii}=\begin{cases}1&\mathrm{\mathrm{if}\;}X_{i,t}\;\mathrm{observed}\\\
0&\mathrm{\mathrm{if}\;}X_{i,t}\;\mathrm{missing}\end{cases}$

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

$\displaystyle\mathrm{vec}(\hat{\bm{\Lambda}})$
$\displaystyle=\left(\sum_{t=1}^{n}\bm{S}_{t|n}\otimes\bm{W}_{t}\right)^{-1}\mathrm{vec}\left(\sum_{t=1}^{n}\bm{W}_{t}\bm{X}_{t}\bm{a}_{t|n}^{\top}\right)\,,$
$\displaystyle\hat{\bm{\Sigma}}_{\bm{\epsilon}}$
$\displaystyle=\frac{1}{n}\sum_{t=1}^{n}\mathrm{diag}\Bigg{[}\bm{W}_{t}\bigg{(}\bm{X}_{t}\bm{X}_{t}^{\top}\
-2\bm{X}_{t}\bm{a}_{t|n}^{\top}\hat{\bm{\Lambda}}^{\top}+\hat{\bm{\Lambda}}\bm{S}_{t|n}\hat{\bm{\Lambda}}^{\top}\bigg{)}\bm{W}_{t}+(\bm{I}-\bm{W}_{t})\hat{\bm{\Sigma}}_{\bm{\epsilon}}^{*}(\bm{I}-\bm{W}_{t})\Bigg{]}\,,$

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 B
for full detail of the EM algorithm.

Doz et al. (2012) and Barigozzi and Luciani (2019) 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 Doz et al. (2012). It is a very flexible framework able to deal
with arbitrary patterns of missing data in $\bm{X}_{t}$, as shown in Bańbura
and Modugno (2014). It enables one to impose restrictions on parameters when
updating them in the EM algorithm. For example, Banbura et al. (2010) 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.

###  2.4 EM-sparse - Mosley et al. (2023)

Mosley et al. (2023) extend the Bańbura and Modugno (2014) 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

$\hat{\bm{\theta}}=\arg\min_{\bm{\theta}}\left[-\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]+\alpha\|\bm{\Lambda}\|_{1}\right]\;,$
(8)

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 (8), Mosley et al. (2023) propose to solve the LASSO problem
using the Alternative Directed Method of Moments (ADMM) by Boyd et al. (2011)
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

$\mathcal{C}(\bm{\Lambda},\bm{Z},\bm{U}):=-\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]+\alpha\|\bm{Z}\|_{1}+\frac{\nu}{2}\|\bm{\Lambda}-\bm{Z}+\bm{U}\|_{F}^{2}\,,$
(9)

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 (8), e.g.

$\displaystyle\arg\min_{\bm{Z}=\bm{\Lambda}}\max_{\bm{U}}\mathcal{C}(\bm{\Lambda},\bm{Z},\bm{U})$
$\displaystyle=\arg\min_{\bm{\Lambda}}\left[-\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]+\alpha\|\bm{\Lambda}\|_{1}\right]\,,$

as (8) is convex in the argument $\bm{\Lambda}$ with all other parameters
fixed, this argument holds for any $\nu>0$.

The augmented Lagrangian (9) can be sequentially minimised via the following
updates777The full derivation of these equations are given in Appendix C.

$\displaystyle\bm{\Lambda}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{\Lambda}}\mathcal{C}(\bm{\Lambda},\bm{Z}^{(k)},\bm{U}^{(k)})$
$\displaystyle\bm{Z}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{Z}}\mathcal{C}(\bm{\Lambda}^{(k+1)},\bm{Z},\bm{U}^{(k)})$
$\displaystyle=\mathrm{soft}(\bm{\Lambda}^{(k+1)}+\bm{U}^{(k)};\alpha/\nu)$
$\displaystyle\bm{U}^{(k+1)}$
$\displaystyle=\bm{U}^{(k)}+\bm{\Lambda}^{(k+1)}-\bm{Z}^{(k+1)}\,.$

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

$\mathrm{vec}(\bm{\Lambda}^{(k+1)})=\left(\sum_{t=1}^{n}\bm{S}_{t|n}\otimes\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{W}_{t}+\nu\bm{I}_{pr}\right)^{-1}\mathrm{vec}\bigg{[}\sum_{t=1}^{n}\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{W}_{t}\bm{X}_{t}\bm{a}_{t|n}^{\top}+\nu(\bm{Z}^{(k)}-\bm{U}^{(k)})\bigg{]}\,.$
(10)

Mosley et al. (2023) further detail how the computation of (10) can be made
significantly faster by exploiting dimensionality reduction, going from a per-
iteration cost of the ADMM of $\mathcal{O}(r^{3}p^{3})$ to
$\mathcal{O}(r^{3}p)$. For details of this speed up refer to Mosley et al.
(2023). We implement this speed up in the package.

###  2.5 An Extension for AR(1) Idiosyncratic Errors

To extend the state-space framework (3)-(4) to AR(1) idiosyncratic errors, as
opposed to i.i.d. white noise, we just follow the methodology laid out in
Bańbura and Modugno (2014) 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

$\displaystyle\bm{\epsilon}_{i,t}$
$\displaystyle=\tilde{\bm{\epsilon}}_{i,t}+\bm{\eta}_{i,t}\,,\quad\bm{\eta}_{i,t}\sim
N(0,\kappa)\,,$ $\displaystyle\tilde{\bm{\epsilon}}_{i,t}$
$\displaystyle=\phi_{i}\tilde{\bm{\epsilon}}_{i,t-1}+\bm{e}_{i,t}\,,\quad\bm{e}_{i,t}\sim
N(0,\sigma_{i}^{2})\,,$

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

$\displaystyle\bm{X}_{t}$
$\displaystyle=\tilde{\bm{\Lambda}}\tilde{\bm{F}}_{t}+\bm{\eta}_{t}\,,\quad\bm{\eta}_{t}\sim
N(0,\bm{\Sigma}_{\bm{\eta}})\,,$ $\displaystyle\tilde{\bm{F}}_{t}$
$\displaystyle=\tilde{\bm{A}}\tilde{\bm{F}}_{t-1}+\tilde{\bm{u}}_{t}\,,\quad\tilde{\bm{u}}_{t}\sim
N(0,\bm{\Sigma}_{\tilde{\bm{u}}})\,,$

where

$\tilde{\bm{F}}_{t}=\left[\begin{array}[]{c}\bm{F}_{t}\\\
\tilde{\bm{\epsilon}}_{t}\end{array}\right]\,,\;\tilde{\bm{u}}_{t}=\left[\begin{array}[]{l}\bm{u}_{t}\\\
\bm{e}_{t}\end{array}\right]\,,\;\tilde{\bm{\Lambda}}=\left[\begin{array}[]{cc}\bm{\Lambda}&\bm{I}\end{array}\right]\,,\;\tilde{\bm{A}}=\left[\begin{array}[]{cc}\bm{A}&\bm{0}\\\
\bm{0}&\bm{\phi}\end{array}\right]\,,$
$\bm{\Sigma}_{\tilde{\bm{u}}}=\left[\begin{array}[]{cc}\bm{\Sigma}_{\bm{u}}&\bm{0}\\\
\bm{0}&\bm{\Sigma}_{\bm{e}}\end{array}\right]\;\text{and}\quad\bm{\Sigma}_{\bm{\eta}}=\kappa\bm{I}_{r+p}\,,$

for $\bm{\phi}=\text{diag}(\phi_{1},\dots,\phi_{p})$ and
$\bm{\Sigma}_{\bm{e}}=\mathbb{E}[\bm{e}_{t}\bm{e}_{t}^{\top}]=\text{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 Bańbura and Modugno (2014) 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}})$.

###  2.6 Tuning the Model

####  2.6.1 For the number of factors

To calculate the number of factors to use in the model we opt to take the
information criteria approach of Bai and Ng (2002). This can be done before
the `SparseDFM()` model is fitted. Bai and Ng (2002) consider 3 information
criteria with different penalties of the form:

$\displaystyle IC_{1}(r)$
$\displaystyle=\log\left(V_{r}(\hat{\bm{F}},\hat{\bm{\Lambda}})\right)+r\left(\frac{n+p}{np}\right)\log\left(\frac{np}{n+p}\right)\,,$
$\displaystyle IC_{2}(r)$
$\displaystyle=\log\left(V_{r}(\hat{\bm{F}},\hat{\bm{\Lambda}})\right)+r\left(\frac{n+p}{np}\right)\log\left(\min\\{n,p\\}\right)\,,$
$\displaystyle IC_{3}(r)$
$\displaystyle=\log\left(V_{r}(\hat{\bm{F}},\hat{\bm{\Lambda}})\right)+r\frac{\log\left(\min\\{n,p\\}\right)}{\min\\{n,p\\}}\,.$

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 2.3.

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 `tuneFactors()` function.

####  2.6.2 For the LASSO sparsity parameter

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

$BIC_{\alpha^{*}}=\log\left(V_{\alpha^{*}}(\hat{\bm{F}},\hat{\bm{\Lambda}})\right)+m\left(\frac{\log(np)}{np}\right)\,,$

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.

##  3 Package Structure

###  3.1 Installation

The sparseDFM package can be installed from CRAN as follows:

        install.packages("sparseDFM")

The development version can be installed from GitHub using the devtools
package (Wickham et al., 2022):

        devtools::install_github("mosleyl/sparseDFM")

To load the package into the R environment, type:

        library(sparseDFM)

###  3.2 Core Function

The core function of the package is `sparseDFM`. This function allows
different estimation methods for dynamic factor models including: principal
components analysis (PCA) for static factors (Stock and Watson, 2002), the
two-stage (2Stage) framework of PCA plus Kalman filter/smoother (Giannone et
al., 2008), the quasi-maximum likelihood approach using the expectation-
maximisation (EM) algorithm to handle arbitrary patterns of missing data
(Bańbura and Modugno, 2014), and the novel sparsified EM approach (EM-sparse)
allowing LASSO regularisation on factor loadings (Mosley et al., 2023). The
function can be implemented using:

        sparseDFM(X, r, q = 0, alphas = logspace(-2,3,100), alg = "EM-sparse",
                  err = "IID",  kalman = "univariate", store.parameters = FALSE,
                  standardize = TRUE, max_iter = 100, threshold = 1e-4)

It takes the following arguments:

  * •

`X`: An $n\times p$ numeric data matrix or data frame representing $p$
stationary time series variables with $n$ observations. Missing data should be
represented by `NA`.

  * •

`r`: The number of factors the model should estimate. This can be selected
using the function `tuneFactors()`. Usually less than 10 factors should be
sufficient.

  * •

`q`: The first `q` variables in `X` will not be regularised. This parameter is
only relevant when `alg = "EM-sparse"`. It is useful when you want to ensure
variables are always loaded onto all factors and hence regularisation is not
applied to these variables. Set the first `q` series in `X` to be these
variables. The default is `q = 0` where all variables are regularised.

  * •

`alphas`: Numeric vector of $\ell_{1}$-norm penalty parameters. This parameter
is only relevant when `alg = "EM-sparse"`. The default grid is set to be
`alphas = logspace(-2, 3, 100)` to ensure a wide range of $\alpha$ values.
Each iteration of the EM algorithm uses the previous $\alpha$’s parameter
estimates as a warm-start to the EM algorithm. The grid search will stop when
all of the values in `alphas` are used for estimation or when a column of the
estimated loadings $\hat{\bm{\Lambda}}$ becomes entirely 0. This will occur if
$\alpha$ is too high. Single $\ell_{1}$-norm penalty parameters can be used as
well as vectors.

  * •

`alg`: The choice of estimation algorithm. Options include:

    * -

`"PCA"`: principal components analysis (PCA) for static factors (Stock and
Watson, 2002);

    * -

`"2Stage"`: the two-stage framework of PCA plus Kalman filter/smoother
(Giannone et al., 2008);

    * -

`"EM"`: The quasi-maximum likelihood approach using the EM algorithm to handle
arbitrary patterns of missing data (Bańbura and Modugno, 2014);

    * -

`"EM-sparse"`: The novel sparsified EM approach allowing LASSO regularisation
on factor loadings (Mosley et al., 2023). This is the default.

  * •

`err`: The choice of distribution for idiosyncratic residuals. Options are
`err = "IID"` for IID white noise residuals and this is the default. Or `err =
"AR1"` for AR(1) residuals. It is recommended to use IID errors for
computational efficiency. See Appendix H for details.

  * •

`kalman`: The choice of Kalman filter and smoother equations to use. Options
include:

    * -

`"multivariate"`: classic Kalman filter and smoother equations seen in Shumway
and Stoffer (1982).

    * -

`"univariate"`: univariate treatment (sequential processing) of the
multivariate equations for fast Kalman filter and smoother seen in Koopman and
Durbin (2000). This is the default.

It is recommended to use `kalman = "univariate"` when `err = "IID"` and
`kalman = "multivariate` when `err = "AR1"`. See Appendix H for details.

  * •

`store.parameters`: A logical parameter indicating whether the estimation
output of parameters and factors should be stored for every value of $\alpha$
considered in `alphas`. This parameter is only relevant when `alg = "EM-
sparse"`. This may be useful when checking alternative outputs other than the
one for the optimal $\alpha$. The default is `store.parameters = FALSE`.

  * •

`standardize`: A logical parameter indicating whether the data should be
standardised or not before estimating the model. The default is `standardized
= TRUE`.

  * •

`max_iter`: An integer representing the maximum number of EM iterations. The
default is `max_iter = 100`.

  * •

`threshold`: Tolerance on EM iterates. See equation (7). Default is `threshold
= 1e-4`.

The function `sparseDFM` returns an S3 object of class `‘sparseDFM’`. This
allows the user to pass the generic functions: `print()`, `summary()`,
`plot()`, `residuals()`, `resids()`, `fitted()`, and `predict()` to the
returned `sparseDFM` object to display generic outputs. The `predict()`
function will itself return an S3 object of class `‘sparseDFM_forecast’` to
allow the user to `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 D). We go into detail on the specific usage
of these options in the empirical applications of Section 4.

When called, the S3 object `sparseDFM` returns a list-of-lists containing the
following:

  * •

`data`: A list containing information on the data used for the model fit.
Including the original data matrix, the initial balanced data matrix from
`fillNA()`, the fitted values $\hat{\bm{\Lambda}}\hat{\bm{F}}$ and more.

  * •

`params`: A list containing the estimated model parameters of $\hat{\bm{A}}$,
$\hat{\bm{\Lambda}}$, $\hat{\bm{\Sigma}}_{u}$ and
$\hat{\bm{\Sigma}}_{\epsilon}$.

  * •

`state`: A list containing the estimated states and state covariances for each
time point considered. I.e. the factors and factor covariances, along with
AR(1) errors and covarainces if `err = AR1`.

  * •

`em`: A list containing information about the convergence of the EM algorithm
and $\ell_{1}$-norm penalty parameter tuning.

  * •

`alpha.output`: Parameter and state outputs for each $\ell_{1}$-norm penalty
parameter in `alphas`. Only stored if `store.parameters = TRUE`.

###  3.3 Other Functions

This package also includes eight other exported functions:

Function |  Brief Description
---|---
`kalmanMultivariate()` |  Function to run the classic multivariate Kalman filter and smoother equations of Shumway and Stoffer (1982).
`kalmanUnivariate()` |  Function to run the fast univariate Kalman filter and smoother equations of Koopman and Durbin (2000).
`fillNA()` |  Function to interpolate missing data in $\bm{X}$. Used to obtain a balanced dataset for the PCA initialisation fit.
`transformData()` |  Function to transform a non-stationary data matrix to become stationary by taking first/second (log) differences/growth rates.
`tuneFactors()` |  Function to determine the number of factors to use based on the Bai and Ng (2002) information criteria and variance explained from PCA.
`missing_data_plot()` |  Visualise the amount of missing data in a data matrix or data frame.
`raggedEdge()` |  Function to generate a ragged edge structure for the end of sample of a data matrix/frame.
`logspace()` |  Function to produce a vector of log10 space values.

##  4 Using the sparseDFM Package

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) website888We use the Q4 2022 release (benchmarked to
2015=100) from
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 D that summarises the various types of plot and other S3 generic
functions that the package contains.

###  4.1 Case Study 1: Exploring UK Inflation

The inflation data is included in the sparseDFM package and we can load it
into the R environment as a data frame called `data` using:

      R> data <- inflation

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 `missing_data_plot()` allows the user to
visualise where missing data is present:

        R> missing_data_plot(data)

Figure 1: Output of _missing_data_plot(data)_ showing the amount of missing
data in the inflation example.

Figure 1 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 `transformData()` and set the `stationary_transform`
parameter to represent a first-difference transform for all variables:

        R> new_data <- transformData(data, stationary_transform = rep(2, ncol(data))))

We store the stationary first-difference data as the variable `new_data`. The
`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: `1` (no change), `2` (first-difference), `3`
(second-difference), `4` (log-first-difference), `5` (log-second-difference),
`6` (growth rate) or `7` (log growth rate). The returned transformed data set
is of the same dimension as the original data with `NA` where missing data is,
for example, the first row of `new_data` will all be `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
`tuneFactors()` that selects the best number of factors to use based on Bai
and Ng (2002) information criteria:

        R> tuneFactors(new_data, type = 2, standardize = TRUE,
                      r.max = min(15,ncol(X)-1), plot = TRUE)

The default information criteria (IC) is `type = 2` (see Section 2.6 for
details). We search for a maximum number of 15 factors and plot the results.
Figure 8 (in Appendix E) shows the output of `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 `summary()` to
summarise the call:

    fit.sdfm <- sparseDFM(new_data, r = 3, alphas = logspace(-2,3,100),
                          alg = ’EM-sparse’, err = ’IID’, kalman = ’univariate’)
    summary(fit.sdfm)
    #>
    #> Call:
    #>
    #>  sparseDFM(X = new_data, r = 3, alphas = logspace(-2, 3, 100),
                  alg = "EM-sparse", err = "IID", kalman = "univariate")
    #>
    #>  Sparse Dynamic Factor Model using EM-sparse with:
    #>
    #>  n = 135 observations,
    #>  p = 36 variables,
    #>  r = 3 factors,
    #>  err = IID
    #>
    #>  The r x r factor transition matrix A
    #>            F1         F2          F3
    #> F1  0.4654857  0.9849589 -0.96021652
    #> F2  0.6433499 -0.3851340 -0.55318409
    #> F3 -0.2577013  0.9030837 -0.02329362
    #>
    #>
    #>  The r x r factor transition error covariance matrix Sigma_u
    #>            u1         u2         u3
    #> u1  3.5334868 -0.2695235 -0.5345583
    #> u2 -0.2695235  2.4666091  0.7245212
    #> u3 -0.5345583  0.7245212  2.5614458

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 8 in Appendix E) using:

    # The best alpha chosen
    fit.sdfm$em$alpha_opt
    #> [1] 0.4641589
    # Plot the BIC values for each alpha
    plot(fit.sdfm, type = ’lasso.bic’)

We are able to display the estimation results using various types of plots -
see Appendix D for options. We can plot the three estimated factors overlaid
on top of the standardised inflation data using

    # Plot all of the estimated factors
    plot(fit.sdfm, type = ’factor’)

and make a boxplot of the residuals for each inflation series using

    # Plot boxplots for the residuals of each variable
    plot(fit.sdfm, type = ’residual’, use.series.names = TRUE)

See Figure 8 and 8 in Appendix E for these plots, respectively. Note, we can
set `use.series.names = FALSE` to avoid plotting the series names. A very
useful plot to understand the factor loading structure is to use:

    # Plot a heatmap of the estimated loadings
    plot(fit.sdfm, type = ’loading.heatmap’, use.series.names = TRUE)

This produces a heatmap (shown in Figure 2) of the estimated factor loadings
with a clear display of sparse loadings displayed by blank squares.

Figure 2: Output of _plot()_ with _type = “loading.heatmap”_ showing a heatmap
of the estimated loadings from Sparse DFM fit on the inflation data.

###  4.2 Case Study 2 - Nowcasting UK Trade in Goods (Exports)

Nowcasting999For a detailed survey on the nowcasting literature see Bańbura et
al. (2013) 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,
Chemicals, Material manufactures, Machinery and transport and Miscellaneous
manufactures.

To try and produce accurate nowcasts for the targets we use a large collection
of 434 monthly indicator series including:

  * •

Index of Production (IoP) \- Movements in the volume of production for the UK
production industries - _2 month lag_ \- _89 series_.

  * •

Consumer Price Inflation (CPI) \- The rate at which the prices of goods and
services bought by households rise or fall - _1 month lag_ \- _166 series_.

  * •

Producer Price Inflation (PPI) \- Changes in the prices of goods bought and
sold by UK manufacturers - _1 month lag_ \- _153 series_.

  * •

Exchange rates \- Sterling exchange rates with 12 popular currencies - _1
month lag_ \- _12 series_.

  * •

Business confidence Index (BCI) \- Opinion surveys on developments in
production, orders and stocks of finished goods in the industry sector - _1
month lag_ \- _1 series_.

  * •

Consumer Confidence Index (CCI) \- Opinion surveys on future developments of
households’ consumption and saving - _1 month lag_ \- _1 series_.

  * •

Google Trends (GT) \- Popularity scores of 14 google search queries related to
trade in goods - _real-time_ \- _14 series_.

The sparseDFM package contains a data frame `exports` containing the 9 monthly
target series of UK exports and the 434 monthly indicator series observed from
January 2004 to October 2022101010All data is openly available from the Office
for National Statistics’ website: https://www.ons.gov.uk/. and can be loaded
into the R environment using

      R> data <- exports

Figure 3 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:

    # Plot the 9 target series using ts.plot with a legend on the right
    def.par <- par(no.readonly = TRUE) # initial graphic parameters
    goods <- data[,1:9]
    layout(matrix(c(1,2),nrow=1), width=c(4,3))
    par(mar=c(5,4,4,0))
    ts.plot(goods, gpars= list(col=10:1,lty=1:10))
    par(mar=c(5,0,4,2))
    plot(c(0,1),type="n", axes=F, xlab="", ylab="")
    legend("center", legend = colnames(goods), col = 10:1, lty = 1:10, cex = 0.7)
    par(def.par) # reset graphic parameters to initial

Figure 3: The target series of the 9 main commodities of UK exports worldwide.

To make the data stationary, we apply first-difference transformations to each
series using:

        R> new_data <- transformData(data, stationary_transform = rep(2, ncol(data))))

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 10 in Appendix F which plots a `missing_data_plot()`
for the last 12 months of the sample using:

    # last 12 months
    data_last12 = tail(data, 12)
    # Missing data plot. Too many variable names so use.names is set to FALSE for
    # clearer output.
    missing_data_plot(data_last12, use.names = FALSE)

We are interested in nowcasting both September and October for the export
targets. We tune for the number of factors using (see Figure 10 in Appendix F
for the output)

        R> tuneFactors(new_data, type = 2, standardize = TRUE,
                      r.max = min(15,ncol(X)-1), plot = TRUE)

According to the Bai and Ng (2002) 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:

    ## Regular DFM fit - takes around 18 seconds
    fit.dfm <- sparseDFM(new_data, r = 4, alg = ’EM’)

    ## Sparse DFM fit - takes around 2 mins to tune
    # set q = 9 as the first 9 variables (targets) should not be regularised
    # L1 penalty grid set to logspace(0.4,1,15) after exploration
    fit.sdfm <- sparseDFM(new_data, r = 4, q = 9, alg = ’EM-sparse’,
                          alphas = logspace(0.4,1,15))

Note, we set the parameter `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:

    # Number of iterations the DFM took to converge
    fit.dfm$em$num_iter
    #> [1] 14

    # Number of iterations the Sparse DFM took to converge at each L1 norm penalty
    fit.sdfm$em$num_iter
    #>  [1] 17  6 14  2  2  2  3  3  3  3 18  3 12  5  5

    # Optimal L1 norm penalty chosen
    fit.sdfm$em$alpha_opt
    #> [1] 4.54091

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 `type = "loading.grouplineplot"`
setting in `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:

    ## Plot the estimated factors for the DFM
    plot(fit.dfm, type = ’factor’)

    ## Plot the estimated factors for the Sparse DFM
    plot(fit.sdfm, type = ’factor’)

    ## Specify the name of the group each indicator belongs too
    # this is a vector of length matching the number of columns of the data
    groups = c(rep(’TiG’,9), rep(’IoP’,89), rep(’CPI’,166), rep(’PPI’,153),
              rep(’Exch’,12), rep(’Conf’,2), rep(’GT’,14))

    ## Specify the colours for each of the groups
    # this is a vector of length matching the number of different groups
    group_cols = c(’black’,’blue’,’red’,’pink’,’green’,’navy’,’brown’)

    ## Plot the DFM group lineplot in a 2 x 2 grid - uses the gridExtra package from CRAN.
    p1 = plot(fit.dfm, type = ’loading.grouplineplot’, loading.factor = 1,
              group.names = groups, group.cols = group_cols)
    p2 = plot(fit.dfm, type = ’loading.grouplineplot’, loading.factor = 2,
              group.names = groups, group.cols = group_cols)
    p3 = plot(fit.dfm, type = ’loading.grouplineplot’, loading.factor = 3,
              group.names = groups, group.cols = group_cols)
    p4 = plot(fit.dfm, type = ’loading.grouplineplot’, loading.factor = 4,
              group.names = groups, group.cols = group_cols)

    grid.arrange(p1, p2, p3, p4, nrow = 2)

    ## Do the same for the Sparse DFM fit
    p5 = plot(fit.sdfm, type = ’loading.grouplineplot’, loading.factor = 1,
              group.names = groups, group.cols = group_cols)
    p6 = plot(fit.sdfm, type = ’loading.grouplineplot’, loading.factor = 2,
              group.names = groups, group.cols = group_cols)
    p7 = plot(fit.sdfm, type = ’loading.grouplineplot’, loading.factor = 3,
              group.names = groups, group.cols = group_cols)
    p8 = plot(fit.sdfm, type = ’loading.grouplineplot’, loading.factor = 4,
              group.names = groups, group.cols = group_cols)

    grid.arrange(p5, p6, p7, p8, nrow = 2)

Figure 4: Estimated factor estimates and corresponding loadings for the
exports application with a regular DFM fit on the top row and a Sparse DFM fit
on the bottom row.

Figure 4 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 `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:

    ## DFM nowcasts (on the differenced data)

    # (1) directly from fit.dfm
    dfm.nowcasts = tail(fit.dfm$data$fitted.unscaled[,1:9],2)

    # which is the same as

    # (2) from fitted()
    dfm.nowcasts = tail(fitted(fit.dfm)[,1:9],2)

    ## Sparse DFM nowcasts (on the differenced data)

    sdfm.nowcasts = tail(fit.sdfm$data$fitted.unscaled[,1:9],2)

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:

    ## August 2022 figures for targets

    obs_aug22 = tail(data,3)[1,1:9]

    ## DFM nowcast for original level

    dfm_sept_nowcast = obs_aug22 + dfm.nowcasts[1,]
    dfm_oct_nowcast = dfm_sept_nowcast + dfm.nowcasts[2,]

    ## Sparse DFM nowcast for original level

    sdfm_sept_nowcast = obs_aug22 + sdfm.nowcasts[1,]
    sdfm_oct_nowcast = sdfm_sept_nowcast + sdfm.nowcasts[2,]

    # Print
    cbind(dfm_sept_nowcast,
    dfm_oct_nowcast,
    sdfm_sept_nowcast,
    sdfm_oct_nowcast)
    #>                                           dfm_sept_nowcast dfm_oct_nowcast
    #> target.Food & live animals                      1466.20494      1487.89722
    #> target.Beverages and tobacco                     892.52105       910.32919
    #> target.Crude materials                           931.82943       935.86939
    #> target.Fuels                                    5434.76988      5281.31709
    #> target.Animal and vegetable oils and fats         81.11094        82.15963
    #> target.Chemicals                                5450.35939      5496.92373
    #> target.Material manufactures                    4250.23561      4348.36403
    #> target.Machinery and transport                 14966.96739     15394.04752
    #> target.Miscellaneous manufactures               4575.18045      4759.21107
    #>                                           sdfm_sept_nowcast sdfm_oct_nowcast
    #> target.Food & live animals                       1373.09109       1379.51984
    #> target.Beverages and tobacco                      806.32814        810.16403
    #> target.Crude materials                            843.30349        828.98577
    #> target.Fuels                                     5444.14078       5284.05460
    #> target.Animal and vegetable oils and fats          79.86784         80.48402
    #> target.Chemicals                                 5233.46691       5237.28705
    #> target.Material manufactures                     3881.44541       3923.95003
    #> target.Machinery and transport                  12776.89494      12788.07209
    #> target.Miscellaneous manufactures                3896.18607       3951.78687

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 G. 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
`kalmanUnivariate()` function on these stored parameters at each time window.

Table 1 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.

|  | Mean | $\bm{25\%}$ | $\bm{50\%}$ | $\bm{75\%}$
---|---|---|---|---|---
Horizon 1 | DFM | 373.72 | 152.403 | 232.856 | 455.077
SDFM | 297.233 | 141.118 | 212.392 | 306.925
Horizon 2 | DFM | 437.15 | 201.821 | 302.672 | 515.081
SDFM | 357.602 | 161.937 | 233.790 | 366.392
Table 1: Pseudo real-time nowcasting exercise results.

##  5 Conclusions

The R package sparseDFM has been developed to allow the fitting of various
estimation methods for dynamic factor models (DFMs) including PCA (Stock and
Watson, 2002), 2Stage (Giannone et al., 2008), EM (Bańbura and Modugno, 2014),
and the novel EM-sparse approach for sparse DFMs (Mosley et al., 2023). We
implement both multivariate (Shumway and Stoffer, 1982), and univariate
(Koopman and Durbin, 2000) Kalman filter and smoother equations in C++ using
the Armadillo library (Eddelbuettel and Sanderson, 2014), 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 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
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.

### Acknowledgments

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.

## References

  * Ahn and Horenstein (2013) Ahn, S. C. and A. R. Horenstein (2013).  Eigenvalue ratio test for the number of factors.  Econometrica 81(3), 1203–1227.
  * Andreou et al. (2019) Andreou, E., P. Gagliardini, E. Ghysels, and M. Rubin (2019).  Inference in group factor models with an application to mixed-frequency data.  Econometrica 87(4), 1267–1305.
  * Bai (2003) Bai, J. (2003).  Inferential theory for factor models of large dimensions.  Econometrica 71(1), 135–171.
  * Bai and Ng (2002) Bai, J. and S. Ng (2002).  Determining the number of factors in approximate factor models.  Econometrica 70(1), 191–221.
  * Bańbura et al. (2013) Bańbura, M., D. Giannone, M. Modugno, and L. Reichlin (2013).  Now-casting and the real-time data flow.  In Handbook of Economic Forecasting, Volume 2, pp. 195–237. Elsevier.
  * Banbura et al. (2010) Banbura, M., D. Giannone, and L. Reichlin (2010).  Nowcasting.  ECB Working Paper.
  * Bańbura and Modugno (2014) Bańbura, M. and M. Modugno (2014).  Maximum likelihood estimation of factor models on datasets with arbitrary pattern of missing data.  Journal of Applied Econometrics 29(1), 133–160.
  * Barigozzi and Luciani (2019) Barigozzi, M. and M. Luciani (2019).  Quasi maximum likelihood estimation and inference of large approximate dynamic factor models via the em algorithm.  arXiv Preprint arXiv:1910.03821.
  * Boyd et al. (2011) Boyd, S., N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al. (2011).  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine Learning 3(1), 1–122.
  * Carroll (1953) Carroll, J. B. (1953).  An analytical solution for approximating simple structure in factor analysis.  Psychometrika 18(1), 23–38.
  * Croux and Exterkate (2011) Croux, C. and P. Exterkate (2011).  Sparse and robust factor modelling.  Tinbergen Institute Discussion Paper TI 122/4.
  * De Jong and Mackinnon (1988) De Jong, P. and M. J. Mackinnon (1988).  Covariances for smoothed estimates in state space models.  Biometrika 75(3), 601–602.
  * De Valk et al. (2019) De Valk, S., D. de Mattos, and P. Ferreira (2019).  Nowcasting: An r package for predicting economic variables using dynamic factor models.  The R Journal 11(1), 230–244.
  * Despois and Doz (2022) Despois, T. and C. Doz (2022).  Identifying and interpreting the factors in factor models via sparsity: Different approaches.  HAL Id: halshs-02235543v3.
  * Doz and Fuleky (2020) Doz, C. and P. Fuleky (2020).  Dynamic factor models.  In Macroeconomic Forecasting in the Era of Big Data, pp. 27–64. Springer.
  * Doz et al. (2011) Doz, C., D. Giannone, and L. Reichlin (2011).  A two-step estimator for large approximate dynamic factor models based on kalman filtering.  Journal of Econometrics 164(1), 188–205.
  * Doz et al. (2012) Doz, C., D. Giannone, and L. Reichlin (2012).  A quasi-maximum likelihood approach for large, approximate dynamic factor models.  Review of Economics and Statistics 94(4), 1014–1024.
  * Durbin and Koopman (2012) Durbin, J. and S. J. Koopman (2012).  Time Series Analysis by State Space Methods.  Oxford University Press.
  * Eddelbuettel and Sanderson (2014) Eddelbuettel, D. and C. Sanderson (2014).  Rcpparmadillo: Accelerating r with high-performance c++ linear algebra.  Computational Statistics and Data Analysis 71, 1054–1063.
  * Fisher (2015) Fisher, A. J. (2015).  Toward a dynamic model of psychological assessment: Implications for personalized care.  Journal of Consulting and Clinical Psychology 83(4), 825\.
  * Forni and Gambetti (2021) Forni, M. and L. Gambetti (2021).  Policy and business cycle shocks: A structural factor model representation of the us economy.  Journal of Risk and Financial Management 14(8), 371.
  * Foroni and Marcellino (2014) Foroni, C. and M. Marcellino (2014).  A comparison of mixed frequency approaches for nowcasting euro area macroeconomic aggregates.  International Journal of Forecasting 30(3), 554–568.
  * Geweke (1977) Geweke, J. (1977).  The dynamic factor analysis of economic time series.  In Latent Variables in Socio-Economic Models. North-Holland.
  * Giannone et al. (2008) Giannone, D., L. Reichlin, and D. Small (2008).  Nowcasting: The real-time informational content of macroeconomic data.  Journal of Monetary Economics 55(4), 665–676.
  * Grassi et al. (2015) Grassi, S., T. Proietti, C. Frale, M. Marcellino, and G. Mazzi (2015).  Euromind-c: A disaggregate monthly indicator of economic activity for the euro area and member countries.  International Journal of Forecasting 31(3), 712–738.
  * Harvey (1996) Harvey, A. (1996).  Intervention analysis with control groups.  International Statistical Review/Revue Internationale de Statistique 64(3), 313–328.
  * Harvey (1990) Harvey, A. C. (1990).  Forecasting, Structural Time Series Models and the Kalman Filter.  Cambridge University Press.
  * Holmes et al. (2012) Holmes, E. E., E. J. Ward, and W. Kellie (2012).  Marss: Multivariate autoregressive state-space models for analyzing time-series data.  The R Journal 4(1), 11.
  * Hopp (2021) Hopp, D. (2021).  nowcastdfm.  https://github.com/dhopp1/nowcastDFM.
  * Jennrich and Sampson (1966) Jennrich, R. I. and P. Sampson (1966).  Rotation for simple loadings.  Psychometrika 31(3), 313–323.
  * Kaiser (1958) Kaiser, H. F. (1958).  The varimax criterion for analytic rotation in factor analysis.  Psychometrika 23(3), 187–200.
  * Koopman and Durbin (2000) Koopman, S. J. and J. Durbin (2000).  Fast filtering and smoothing for multivariate state space models.  Journal of Time Series Analysis 21(3), 281–296.
  * Kose et al. (2008) Kose, M. A., C. Otrok, and C. H. Whiteman (2008).  Understanding the evolution of world business cycles.  Journal of International Economics 75(1), 110–130.
  * Krantz and Bagdziunas (2022) Krantz, S. and R. Bagdziunas (2022).  dfms: Dynamic Factor Models.  R Package Version 0.1.3.
  * Kristensen (2017) Kristensen, J. T. (2017).  Diffusion indexes with sparse loadings.  Journal of Business & Economic Statistics 35(3), 434–451.
  * Lee and Baldick (2016) Lee, D. and R. Baldick (2016).  Load and wind power scenario generation through the generalized dynamic factor model.  IEEE Transactions on Power Systems 32(1), 400–410.
  * Luciani (2015) Luciani, M. (2015).  Monetary policy and the housing market: A structural factor analysis.  Journal of Applied Econometrics 30(2), 199–218.
  * Mariano and Murasawa (2010) Mariano, R. S. and Y. Murasawa (2010).  A coincident index, common factors, and monthly real gdp.  Oxford Bulletin of Economics and Statistics 72(1), 27–46.
  * Molenaar (1985) Molenaar, P. (1985).  A dynamic factor model for the analysis of multivariate time series.  Psychometrika 50(2), 181–202.
  * Mosley et al. (2023) Mosley, L., T.-S. Chan, and A. Gibberd (2023).  The sparse dynamic factor model: A regularised quasi-maximum likelihood approach.
  * Onatski (2009) Onatski, A. (2009).  Testing hypotheses about the number of factors in large factor models.  Econometrica 77(5), 1447–1479.
  * Poncela et al. (2021) Poncela, P., E. Ruiz, and K. Miranda (2021).  Factor extraction using kalman filter and smoothing: This is not just another survey.  International Journal of Forecasting 37(4), 1399–1425.
  * Rohe and Zeng (2020) Rohe, K. and M. Zeng (2020).  Vintage factor analysis with varimax performs statistical inference.  arXiv Preprint arXiv:2004.05387.
  * Sargent et al. (1977) Sargent, T. J., C. A. Sims, et al. (1977).  Business cycle modeling without pretending to have too much a priori economic theory.  New Methods in Business Cycle Research 1, 145–168.
  * Shumway and Stoffer (1982) Shumway, R. H. and D. S. Stoffer (1982).  An approach to time series smoothing and forecasting using the em algorithm.  Journal of Time Series Analysis 3(4), 253–264.
  * Stock and Watson (2011) Stock, J. H. and M. Watson (2011).  Dynamic factor models.  Oxford Handbooks Online.
  * Stock and Watson (2002) Stock, J. H. and M. W. Watson (2002).  Forecasting using principal components from a large number of predictors.  Journal of the American Statistical Association 97(460), 1167–1179.
  * Watson and Engle (1983) Watson, M. W. and R. F. Engle (1983).  Alternative algorithms for the estimation of dynamic factor, mimic and varying coefficient regression models.  Journal of Econometrics 23(3), 385–400.
  * Wickham et al. (2022) Wickham, H., J. Hester, W. Chang, and J. Bryan (2022).  devtools: Tools to Make Developing R Packages Easier.  R Package Version 2.4.5.
  * Wu et al. (2013) Wu, H., S. Chan, K. Tsui, and Y. Hou (2013).  A new recursive dynamic factor analysis for point and interval forecast of electricity price.  IEEE Transactions on Power Systems 28(3), 2352–2365.
  * Zou et al. (2006) Zou, H., T. Hastie, and R. Tibshirani (2006).  Sparse principal component analysis.  Journal of Computational and Graphical Statistics 15(2), 265–286.

##  Appendix A Multivariate and Univariate Kalman Filter and Smoother
Equations

In this Appendix we provide the classic multivariate Kalman filter and
smoother (KFS) equations of Shumway and Stoffer (1982) and the fast univariate
treatment KFS equations of Koopman and Durbin (2000) that can be used to
estimate the mean and covariance of the latent states in a state space
framework (3)-(4) based on observations of the system. In a DFM, the latent
states are the factors111111Note that in the extension to AR(1) idiosyncratic
errors, as discussed in Section 2.5, the latent state involves the factors
augmented with the AR(1) errors. The KFS equations in this Appendix are
analogous for both cases with the notation $\bm{F}_{t}$ representing the
latent state - either just the factors or the factors augmented with AR(1)
errors. Similarly, the parameters can be swapped, $\bm{\theta}$ with
$\tilde{\bm{\theta}}$. I.e. these equations are for a general state-space
model of the form (3)-(4)., $\bm{F}_{t}$. We use the following notation for
the conditional mean and covariances of the state:

$\displaystyle\bm{a}_{t|s}$
$\displaystyle=\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{s}]\,,$
$\displaystyle\bm{P}_{t|s}$
$\displaystyle=\mathrm{Var}[\bm{F}_{t}|\bm{\Omega}_{s}]\,,$
$\displaystyle\bm{P}_{t,t-1|s}$
$\displaystyle=\mathrm{Cov}[\bm{F}_{t},\bm{F}_{t-1}|\bm{\Omega}_{s}]\,,$

conditional on all information we have observed up to a time $s$, denoted by
$\bm{\Omega}_{s}$.

Filtering aims to make recursive predictions of the mean and covariance of the
state at $t$ given information up to $t-1$, i.e. finds $\bm{a}_{t|t-1}$ and
$\bm{P}_{t|t-1}$ starting from $t=1$ to $t=n$. Smoothing then works backwards
in time to find smoothed estimates of the mean and covariance of the state
given all information, i.e. finds $\bm{a}_{t|n}$ and $\bm{P}_{t|n}$ starting
from $t=n$ to $t=1$.

###  A.1 Multivariate KFS Equations - Shumway and Stoffer (1982)

The classic approach to filtering and smoothing found in Shumway and Stoffer
(1982) and detailed in Harvey (1990) is based on considering the contribution
of every variable in the observational vector at each successive time point.
We therefore describe this as a multivariate approach and set the input in
`SparseDFM()` as `kalman = multivariate`.

Assuming the set of parameters
$\bm{\theta}=(\bm{\Lambda},\bm{A},\bm{\Sigma}_{\bm{\epsilon}},\bm{\Sigma}_{\bm{u}})$
and initial conditions $\bm{a}_{0|0}$ and $\bm{P}_{0|0}$ are known, the
filtering equations for $t=1,\dots,n$ are

$\displaystyle\bm{a}_{t|t-1}$ $\displaystyle=\bm{A}\bm{a}_{t-1|t-1}$
(predicted factor estimate) $\displaystyle\bm{P}_{t|t-1}$
$\displaystyle=\bm{A}\bm{P}_{t-1|t-1}\bm{A}^{\top}+\bm{\Sigma}_{u}$ (predicted
factor covariance) $\displaystyle\bm{v}_{t}$
$\displaystyle=\bm{X}_{t}-\bm{\Lambda}\bm{a}_{t|t-1}$ (innovation error)
$\displaystyle\bm{C}_{t}$
$\displaystyle=\bm{\Lambda}\bm{P}_{t|t-1}\bm{\Lambda}^{\top}+\bm{\Sigma}_{\bm{\epsilon}}$
(innovation covariance) $\displaystyle\bm{K}_{t}$
$\displaystyle=\bm{P}_{t|t-1}\bm{\Lambda}^{\top}\bm{C}_{t}^{-1}$ (Kalman gain)
$\displaystyle\bm{a}_{t|t}$
$\displaystyle=\bm{a}_{t|t-1}+\bm{K}_{t}\bm{v}_{t}$ (updated factor estimate)
$\displaystyle\bm{P}_{t|t}$
$\displaystyle=\bm{P}_{t|t-1}-\bm{K}_{t}\bm{\Lambda}\bm{P}_{t|t-1}$ (updated
factor covariance)

and the smoothing equations for $t=n,\dots,1$ are

$\displaystyle\bm{J}_{t-1}$
$\displaystyle=\bm{P}_{t-1|t-1}\bm{A}^{\top}(\bm{P}_{t|t-1})^{-1}$ (smoother
gain) $\displaystyle\bm{a}_{t-1|n}$
$\displaystyle=\bm{a}_{t-1|t-1}+\bm{J}_{t-1}(\bm{F}_{t|n}-\bm{A}\bm{a}_{t-1|t-1})$
(smoothed factor estimate) $\displaystyle\bm{P}_{t-1|n}$
$\displaystyle=\bm{P}_{t-1|t-1}+\bm{J}_{t-1}(\bm{P}_{t|n}-\bm{P}_{t|t-1})\bm{J}_{t-1}^{\top}$
(smoothed factor covariance)

where we initialise $\bm{a}_{n|n}$ and $\bm{P}_{n|n}$ based on the last filter
estimates. At any stage where an observation $X_{t,i}$ is missing at time $t$,
we omit this variable from the calculation of the Kalman gain $\bm{K}_{t}$.
The lagged-covariance matrix $\bm{P}_{t,t-1|n}$ is found using backwards
recursions on:

$\bm{P}_{t-1,t-2|n}=\bm{P}_{t-1|t-1}\bm{J}_{t-2}^{\top}+\bm{J}_{t-1}(\bm{P}_{t,t-1|n}-\bm{A}\bm{P}_{t-1|t-1})\bm{J}_{t-2}^{\top}\,,$
(11)

for $t=n,\dots,2$, where

$\bm{P}_{n,n-1|n}=(\bm{I}-\bm{K}_{n}\bm{\Lambda})\bm{A}\bm{P}_{n-1|n-1}\,.$

The filtering equations require the covariance of innovation errors
$\bm{C}_{t}$ to be inverted. This can be difficult when $p$ is large. Harvey
(1990) show this inversion can be made easier by making use of the Woodbury
identity and find:

$\bm{C}_{t}^{-1}=\bm{\Sigma}_{e}^{-1}-\bm{\Sigma}_{e}^{-1}\bm{\Lambda}(\bm{P}_{t|t-1}^{-1}+\bm{\Lambda}^{\top}\bm{\Sigma}_{e}^{-1}\bm{\Lambda})^{-1}\bm{\Lambda}^{\top}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\,,$
(12)

which is valid when $\bm{\Sigma}_{\bm{\epsilon}}$ and $\bm{P}_{t|t-1}$ are
non-singular. This inversion is easy to evaluate if $\bm{\Sigma}_{e}$ is
diagonal (as assumed in exact factors models). In cases when $p$ is much
larger than $r$, as is the case in the majority of DFM problems, the inversion
of $\bm{P}_{t|t-1}$ is a much easier operation than the inversion of
$\bm{C}_{t}$ directly.

Harvey (1990) further show that the determinant of $\bm{C}_{t}$ can be written
as:

$|\bm{C}_{t}|=|\bm{\Sigma}_{\bm{\epsilon}}|\times|\bm{P}_{t|t-1}|\times|\bm{P}_{t|t-1}^{-1}+\bm{\Lambda}^{\top}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{\Lambda}|\,.$
(13)

Both (12) and (13) are used in the calculation of the log-likelihood of the
innovations:

$\log\mathcal{L}(\bm{X};\bm{\theta})=-\frac{1}{2}\sum_{t=1}^{n}\left(p_{t}\log(2\pi)+\log|\bm{C}_{t}|+\bm{v}_{t}\bm{C}_{t}^{-1}\bm{v}_{t}\right)\,.$
(14)

where $p_{t}$ is the number of variables observed (not-missing) at $t$ and
$\bm{\theta}$ is the collection of all parameters in (1)-(2). Note, this log-
likelihood is used to determine convergence of the EM algorithm as in (7) of
the EM approaches in Sections 2.3 and 2.4.

###  A.2 Univariate KFS Equations - Koopman and Durbin (2000)

The classic KFS equations can be slow when $p$ is large. Since we assume
$\bm{\Sigma}_{\bm{\epsilon}}$ is diagonal, we can equivalently filter the
observations $\bm{X}_{t}$ one element at a time, as opposed to all together as
in the classic approach (Durbin and Koopman, 2012; Koopman and Durbin, 2000).
As matrix inversion becomes scalar divisions, huge speedups are possible. The
input to the `SparseDFM()` function is `kalman = univariate`.

Let us define the individual elements
$\bm{X}_{t}=(X_{t,1},\dots,X_{t,p})^{\top}$,
$\bm{\Lambda}=(\bm{\Lambda}_{1}^{\top},\dots,\bm{\Lambda}_{p}^{\top})^{\top}$,
$\bm{\Sigma}_{\bm{\epsilon}}=diag(\sigma_{\epsilon
1}^{2},\dots,\sigma_{\epsilon p}^{2})$ and define new notation similar to that
found in Durbin and Koopman (2012):

$\displaystyle\bm{a}_{t,i}$
$\displaystyle=\mathbb{E}[\bm{F}_{t,i}|\bm{\Omega}_{t-1},X_{t,1},\dots,X_{t,i-1}]\,,$
$\displaystyle\bm{a}_{t,1}$
$\displaystyle=\mathbb{E}[\bm{F}_{t,1}|\bm{\Omega}_{t-1}]\,,$
$\displaystyle\bm{P}_{t,i}$
$\displaystyle=\mathrm{Var}[\bm{F}_{t,i}|\bm{\Omega}_{t-1},X_{t,1},\dots,X_{t,i-1}]\,,$
$\displaystyle\bm{P}_{t,1}$
$\displaystyle=\mathrm{Var}[\bm{F}_{t,1}|\bm{\Omega}_{t-1}]\,,$

for $i=1,\dots,p$ and $t=1,\dots,n$.

The univariate treatment (Koopman and Durbin, 2000) iterates the below
equations. These update equations are equivalent in form to the multivariate
ones of the classic Shumway and Stoffer (1982) approach except that the $t$
subscript now becomes a $t,i$ subscript, and the $t|t$ subscript now becomes
$t,i+1$. This means that we are filtering spatiotemporally, processing one $i$
after another in a sequential manner, instead of updating all $p$ measurements
at once:

$\displaystyle v_{t,i}$ $\displaystyle=X_{t,i}-\bm{\Lambda}_{i}\bm{a}_{t,i},$
$\displaystyle C_{t,i}$
$\displaystyle=\bm{\Lambda}_{i}\bm{P}_{t,i}\bm{\Lambda}_{i}^{\top}+\sigma^{2}_{\epsilon,i},$
$\displaystyle\bm{K}_{t,i}$
$\displaystyle=\bm{P}_{t,i}\bm{\Lambda}_{i}^{\top}C_{t,i}^{-1},$
$\displaystyle\bm{a}_{t,i+1}$
$\displaystyle=\bm{a}_{t,i}+\bm{K}_{t,i}v_{t,i},$
$\displaystyle\bm{P}_{t,i+1}$
$\displaystyle=\bm{P}_{t,i}-\bm{K}_{t,i}C_{t,i}\bm{K}_{t,i}^{\top},$

for $i=1,\ldots,p$ and $t=1,\ldots,n$. If $X_{t,i}$ is missing or $C_{t,i}$ is
zero, omit the term containing $\bm{K}_{t,i}$. The transition to $t+1$ is
given by the following prediction equations:

$\displaystyle\bm{a}_{t+1,1}$ $\displaystyle=\bm{A}\bm{a}_{t,p+1},$
$\displaystyle\bm{P}_{t+1,1}$
$\displaystyle=\bm{A}\bm{P}_{t,p+1}\bm{A}^{\top}+\bm{\Sigma}_{\bm{u}}.$

These prediction equations are exactly the same as the multivariate ones
(i.e., predictions are not treated sequentially but all at once). From our
perspective, univariate treatment is better called univariate updates plus
multivariate predictions.

Unlike Shumway and Stoffer (1982), the measurement update comes before the
transition; however, we can revert to doing the transition first if our
initial state means and covariances start from $t=0$ instead of $t=1$.
Likewise, univariate smoothing is defined by:

$\displaystyle\bm{L}_{t,i}$
$\displaystyle=\bm{I}_{m}-\bm{K}_{t,i}\bm{\Lambda}_{i},$
$\displaystyle\bm{r}_{t,i-1}$
$\displaystyle=\bm{\Lambda}_{i}^{\top}C_{t,i}^{-1}v_{t,i}+\bm{L}_{t,i}^{\top}\bm{r}_{t,i},$
$\displaystyle\bm{N}_{t,i-1}$
$\displaystyle=\bm{\Lambda}_{i}^{\top}C_{t,i}^{-1}\bm{\Lambda}_{i}+\bm{L}_{t,i}^{\top}\bm{N}_{t,i}\bm{L}_{t,i},$
$\displaystyle\bm{r}_{t-1,p}$ $\displaystyle=\bm{A}^{\top}\bm{r}_{t,0},$
$\displaystyle\bm{N}_{t-1,p}$ $\displaystyle=\bm{A}^{\top}\bm{N}_{t,0}\bm{A},$

for $i=p,\ldots,1$ and $t=n,\ldots,1$, with $\bm{r}_{n,p}$ and $\bm{N}_{n,p}$
initialised to $0$. Again, if $X_{t,i}$ is missing or $C_{t,i}$ is zero, drop
the terms containing $\bm{K}_{t,i}$. Finally, the equations for $\bm{a}_{t|n}$
and $\bm{P}_{t|n}$ are:

$\displaystyle\bm{a}_{t|n}$
$\displaystyle=\bm{a}_{t,1}+\bm{P}_{t,1}\bm{r}_{t,0},$
$\displaystyle\bm{P}_{t|n}$
$\displaystyle=\bm{P}_{t,1}-\bm{P}_{t,1}\bm{N}_{t,0}\bm{P}_{t,1}.$

These results will be equivalent to $\bm{a}_{t|n}$ and $\bm{P}_{t|n}$ from the
classic multivariate approach, yet obtained with substantial improvement in
computational efficiency.

In order to calculate the cross-covariance matrix $\bm{P}_{t,t-1|n}$, we use
De Jong and Mackinnon (1988)’s theorem:

$\bm{P}_{t-1,t|n}=\bm{P}_{t-1|t-1}\bm{A}^{\top}(\bm{P}_{t|t-1})^{-1}\bm{P}_{t|n},$

which, after transposition, gives

$\bm{P}_{t,t-1|n}=\bm{P}_{t|n}(\bm{P}_{t|t-1})^{-1}\bm{A}\bm{P}_{t-1|t-1}.$

The log-likelihood for the above filter (Durbin and Koopman, 2012; Koopman and
Durbin, 2000) is given below:

$\displaystyle\log\mathcal{L}(\bm{X};\bm{\theta})$
$\displaystyle=-\frac{1}{2}\sum_{t=1}^{n}\sum_{i=1}^{p}\iota_{t,i}\left(\log(2\pi)+\log|C_{t,i}|+v_{t,i}^{\top}C_{t,i}^{-1}v_{t,i}\right)$
(15)
$\displaystyle=-\frac{1}{2}\sum_{t=1}^{n}\sum_{i=1}^{p}\iota_{t,i}\left(\log(2\pi)+\log
C_{t,i}+\frac{v_{t,i}^{2}}{C_{t,i}}\right),$ (16)

as $C_{t,i}$ and $v_{t,i}$ are scalars. Here, $\iota_{t,i}$ is zero if
$X_{t,i}$ is missing or $C_{t,i}$ is zero, or one otherwise.

##  Appendix B EM Algorithm Derivations for Section 2.3

The EM algorithm of Bańbura and Modugno (2014) involves calculating and
maximising the expected log-likelihood of the DFM conditional on the available
information $\bm{\Omega}_{n}$. Given the log-likelihood in (2.3), the
conditional expected log-likelihood is

$\displaystyle\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]=$
$\displaystyle-\frac{1}{2}\log|\mathcal{\bm{P}}_{0}|-\mathrm{tr}\left\\{\mathcal{\bm{P}}_{0}^{-1}\mathbb{E}\left[(\bm{F}_{0}-\bm{\alpha}_{0})(\bm{F}_{0}-\bm{\alpha}_{0})^{\top}|\bm{\Omega}_{n}\right]\right\\}$
(17)
$\displaystyle-\frac{n}{2}\log|\bm{\Sigma}_{\bm{u}}|-\frac{1}{2}\sum_{t=1}^{n}\mathrm{tr}\left\\{\bm{\Sigma}_{\bm{u}}^{-1}\mathbb{E}\left[\bm{u}_{t}\bm{u}_{t}^{\top}|\bm{\Omega}_{n}\right]\right\\}$
(18)
$\displaystyle-\frac{n}{2}\log|\bm{\Sigma}_{\bm{\epsilon}}|-\frac{1}{2}\sum_{t=1}^{n}\mathrm{tr}\left\\{\bm{\Sigma}_{\bm{\epsilon}}^{-1}\mathbb{E}\left[\bm{\epsilon}_{t}\bm{\epsilon}_{t}^{\top}|\bm{\Omega}_{n}\right]\right\\}\;.$
(19)

We refer to
$\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]$
just as $\mathcal{E}$ to make notation easier. In each derivation we use the
linearity of conditional expectation:
$\mathrm{Cov}(X,Y|A)=\mathbb{E}(XY|A)-\mathbb{E}(X|A)\mathbb{E}(Y|A)$ for two
random variables $X,Y$ and non-zero event $A$. We also use the derivative
property $\partial(\log|\bm{Z}|)/\partial(\bm{Z}^{-1})=-\bm{Z}$ for a matrix
$\bm{Z}$. Recall, we define the conditional mean and covariances of the state
as:

$\displaystyle\bm{a}_{t|n}$
$\displaystyle=\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\,,$
$\displaystyle\bm{P}_{t|n}$
$\displaystyle=\mathrm{Var}[\bm{F}_{t}|\bm{\Omega}_{n}]\,,$
$\displaystyle\bm{P}_{t,t-1|n}$
$\displaystyle=\mathrm{Cov}[\bm{F}_{t},\bm{F}_{t-1}|\bm{\Omega}_{n}]\,,$

conditional on all information we have observed up to $n$, denoted by
$\bm{\Omega}_{n}$, and let
$\bm{S}_{t|n}=\bm{a}_{t|n}\bm{a}_{t|n}^{\top}+\bm{P}_{t|n}$, and
$\bm{S}_{t,t-1|n}=\bm{a}_{t|n}\bm{a}_{t-1|n}^{\top}+\bm{P}_{t,t-1|n}$.

Expanding out the conditional expectation in the initial state distribution
part (17) we obtain:

$\displaystyle\mathbb{E}[(\bm{F}_{0}-\bm{\alpha}_{0})(\bm{F}_{0}-\bm{\alpha}_{0})^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\mathbb{E}[\bm{F}_{0}\bm{F}_{0}^{\top}-2\bm{F}_{0}\bm{\alpha}_{0}^{\top}+\bm{\alpha}_{0}\bm{\alpha}_{0}^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\mathbb{E}[\bm{F}_{0}\bm{F}_{0}^{\top}|\bm{\Omega}_{n}]-2\mathbb{E}[\bm{F}_{0}|\bm{\Omega}_{n}]\bm{\alpha}_{0}^{\top}+\bm{\alpha}_{0}\bm{\alpha}_{0}^{\top}$
$\displaystyle=\mathbb{E}[\bm{F}_{0}|\bm{\Omega}_{n}]\mathbb{E}[\bm{F}_{0}^{\top}|\bm{\Omega}_{n}]+\mathrm{Var}[\bm{F}_{0}|\bm{\Omega}_{n}]-2\mathbb{E}[\bm{F}_{0}|\bm{\Omega}_{n}]\bm{\alpha}_{0}^{\top}+\bm{\alpha}_{0}\bm{\alpha}_{0}^{\top}$
$\displaystyle=\bm{a}_{0|n}\bm{a}_{0|n}^{\top}+\bm{P}_{0|n}-2\bm{a}_{0|n}\bm{\alpha}_{0}^{\top}+\bm{\alpha}_{0}\bm{\alpha}_{0}^{\top}=\bm{S}_{0|n}-2\bm{a}_{0|n}\bm{\alpha}_{0}^{\top}+\bm{\alpha}_{0}\bm{\alpha}_{0}^{\top}\,.$

Differentiating $\mathcal{E}$ with respect to $\bm{\alpha}_{0}$ and
$\mathcal{\bm{P}}_{0}^{-1}$ we obtain:

$\displaystyle\frac{\partial\mathcal{E}}{\partial\bm{\alpha}_{0}}$
$\displaystyle=-\frac{1}{2}\mathcal{\bm{P}}_{0}^{-1}\frac{\partial}{\partial\bm{\alpha}_{0}}\mathbb{E}[(\bm{F}_{0}-\bm{\alpha}_{0})(\bm{F}_{0}-\bm{\alpha}_{0})^{\top}|\bm{\Omega}_{n}]=-\frac{1}{2}\mathcal{\bm{P}}_{0}^{-1}\left(-2\bm{a}_{0|n}+2\bm{\alpha}_{0}\right)$
$\displaystyle\frac{\partial\mathcal{E}}{\partial\mathcal{\bm{P}}_{0}^{-1}}$
$\displaystyle=\frac{1}{2}\mathcal{\bm{P}}_{0}-\frac{1}{2}\mathbb{E}[(\bm{F}_{0}-\bm{\alpha}_{0})(\bm{F}_{0}-\bm{\alpha}_{0})^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\frac{1}{2}\mathcal{\bm{P}}_{0}-\frac{1}{2}\left(\bm{S}_{0|n}-2\bm{a}_{0|n}\bm{\alpha}_{0}^{\top}+\bm{\alpha}_{0}\bm{\alpha}_{0}^{\top}\right)\,.$

Setting both derivatives equal to 0, we first get the result

$\hat{\bm{\alpha}}_{0}=\bm{a}_{0|n}\,,$

and use this to obtain

$\hat{\mathcal{\bm{P}}}_{0}=\bm{P}_{0|n}\,.$

Expanding out the expectation in the state transition part (18) we obtain:

$\displaystyle\mathbb{E}\left[\bm{u}_{t}^{\top}\bm{u}_{t}|\bm{\Omega}_{n}\right]$
$\displaystyle=\mathbb{E}[(\bm{F}_{t}-\bm{A}\bm{F}_{t-1})(\bm{F}_{t}-\bm{A}\bm{F}_{t-1})^{\top}|\bm{\Omega}_{n}]=\mathbb{E}[\bm{F}_{t}\bm{F}_{t}^{\top}-2\bm{F}_{t}\bm{F}_{t-1}^{\top}\bm{A}^{\top}+\bm{A}\bm{F}_{t-1}\bm{F}_{t-1}^{\top}\bm{A}^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\mathbb{E}[\bm{F}_{t}\bm{F}_{t}^{\top}|\bm{\Omega}_{n}]-2\mathbb{E}[\bm{F}_{t}\bm{F}_{t-1}^{\top}|\bm{\Omega}_{n}]\bm{A}^{\top}+\bm{A}\mathbb{E}[\bm{F}_{t-1}\bm{F}_{t-1}^{\top}|\bm{\Omega}_{n}]\bm{A}^{\top}$
$\displaystyle=\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\mathbb{E}[\bm{F}_{t}^{\top}|\bm{\Omega}_{n}]+\mathrm{Var}[\bm{F}_{t}|\bm{\Omega}_{n}]-2\left(\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\mathbb{E}[\bm{F}_{t-1}^{\top}|\bm{\Omega}_{n}]+\mathrm{Cov}[\bm{F}_{t},\bm{F}_{t-1}^{\top}|\bm{\Omega}_{n}]\right)\bm{A}^{\top}$
$\displaystyle\quad+\bm{A}\left(\mathbb{E}[\bm{F}_{t-1}|\bm{\Omega}_{n}]\mathbb{E}[\bm{F}_{t-1}^{\top}|\bm{\Omega}_{n}]+\mathrm{Var}[\bm{F}_{t-1}|\bm{\Omega}_{n}]\right)\bm{A}^{\top}$
$\displaystyle=\bm{S}_{t|n}-2\bm{S}_{t,t-1|n}\bm{A}^{\top}+\bm{A}\bm{S}_{t-1|n}\bm{A}^{\top}\,.$

Differentiating $\mathcal{E}$ with respect to $\bm{A}$ and
$\bm{\Sigma}_{\bm{u}}^{-1}$ we obtain:

$\displaystyle\frac{\partial\mathcal{E}}{\partial\bm{A}}$
$\displaystyle=-\frac{1}{2}\bm{\Sigma}_{\bm{u}}^{-1}\frac{\partial}{\partial\bm{A}}\sum_{t=1}^{n}\mathbb{E}\left[\bm{u}_{t}^{\top}\bm{u}_{t}|\bm{\Omega}_{n}\right]=-\frac{1}{2}\bm{\Sigma}_{\bm{u}}^{-1}\sum_{t=1}^{n}\left(-2\bm{S}_{t,t-1|n}+2\bm{A}\bm{S}_{t-1|n}\right)\,,$
$\displaystyle\frac{\partial\mathcal{E}}{\partial\bm{\Sigma}_{\bm{u}}^{-1}}$
$\displaystyle=\frac{n}{2}\bm{\Sigma}_{\bm{u}}-\frac{1}{2}\sum_{t=1}^{n}\mathbb{E}\left[\bm{u}_{t}^{\top}\bm{u}_{t}|\bm{\Omega}_{n}\right]=\frac{n}{2}\bm{\Sigma}_{\bm{u}}-\frac{1}{2}\left(\bm{S}_{t|n}-2\bm{S}_{t,t-1|n}\bm{A}^{\top}+\bm{A}\bm{S}_{t-1|n}\bm{A}^{\top}\right)\,.$

Setting both derivatives equal to 0, we first get the result

$\hat{\bm{A}}=\left(\sum_{t=1}^{n}\bm{S}_{t-1|n}\right)^{-1}\left(\sum_{t=1}^{n}\bm{S}_{t,t-1|n}\right)\,,$

and use this to obtain

$\hat{\bm{\Sigma}}_{\bm{u}}=\frac{1}{n}\sum_{t=1}^{n}\left[\bm{S}_{t|n}-\hat{\bm{A}}\left(\bm{S}_{t-1,t|n}\right)\right]\,.$

Using the selection matrix $\bm{W}_{t}$ defined in Section 2.3, with the
property $\bm{X}_{t}=\bm{W}_{t}\bm{X}_{t}+(\bm{I}-\bm{W}_{t})\bm{X}_{t}$, we
can expand out the expectation in the measurement equation part (19) and
obtain:

$\displaystyle\mathbb{E}\left[\bm{\epsilon}_{t}\bm{\epsilon}_{t}^{\top}|\bm{\Omega}_{n}\right]=\mathbb{E}[(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\mathbb{E}[[\bm{W}_{t}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})+(\bm{I}-\bm{W}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})][\bm{W}_{t}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})+(\bm{I}-\bm{W}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})]^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\mathbb{E}[\bm{W}_{t}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})\bm{W}_{t}+(\bm{I}-\bm{W}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}\bm{W}_{t}$
$\displaystyle\quad+\bm{W}_{t}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}(\bm{I}-\bm{W}_{t})+(\bm{I}-\bm{W}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}(\bm{I}-\bm{W}_{t})|\bm{\Omega}_{n}]\,.$

Using the law of iterated expectations:

$\mathbb{E}[(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}|\bm{\Omega}_{n}]=\mathbb{E}\left[\mathbb{E}[(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}|\bm{F},\bm{\Omega}_{n}]|\bm{\Omega}_{n}\right]\,,$

Bańbura and Modugno (2014) show

$\displaystyle\mathbb{E}[\bm{W}_{t}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}(\bm{I}-\bm{W}_{t})|\bm{F},\bm{\Omega}_{n}]=0\,,$
$\displaystyle\mathbb{E}[(\bm{I}-\bm{W}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}(\bm{I}-\bm{W}_{t})|\bm{F},\bm{\Omega}_{n}]=(\bm{I}-\bm{W}_{t})\bm{\Sigma}_{\bm{\epsilon}}(\bm{I}-\bm{W}_{t})\,,$

and

$\displaystyle\mathbb{E}[\bm{W}_{t}(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}\bm{W}_{t}|\bm{\Omega}_{n}]$
$\displaystyle=\bm{W}_{t}\bm{X}_{t}\bm{X}_{t}^{\top}\bm{W}_{t}-\bm{W}_{t}\bm{X}_{t}\mathbb{E}[\bm{F}_{t}^{\top}|\bm{\Omega}_{n}]\bm{\Lambda}^{\top}\bm{W}_{t}$
$\displaystyle\quad-\bm{W}_{t}\bm{\Lambda}\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\bm{X}_{t}^{\top}\bm{W}_{t}+\bm{W}_{t}\bm{\Lambda}\mathbb{E}[\bm{F}_{t}\bm{F}_{t}^{\top}|\bm{\Omega}_{n}]\bm{\Lambda}^{\top}\bm{W}_{t}\,.$

Hence, we get the result

$\displaystyle\mathbb{E}\left[\bm{\epsilon}_{t}\bm{\epsilon}_{t}^{\top}|\bm{\Omega}_{n}\right]$
$\displaystyle=\mathbb{E}[(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})(\bm{X}_{t}-\bm{\Lambda}\bm{F}_{t})^{\top}|\bm{\Omega}_{n}]$
$\displaystyle=\bm{W}_{t}\bm{X}_{t}\bm{X}_{t}^{\top}\bm{W}_{t}-\bm{W}_{t}\bm{X}_{t}\mathbb{E}[\bm{F}_{t}^{\top}|\bm{\Omega}_{n}]\bm{\Lambda}^{\top}\bm{W}_{t}-\bm{W}_{t}\bm{\Lambda}\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\bm{X}_{t}^{\top}\bm{W}_{t}$
$\displaystyle\quad+\bm{W}_{t}\bm{\Lambda}\left(\mathbb{E}[\bm{F}_{t}|\bm{\Omega}_{n}]\mathbb{E}[\bm{F}_{t}^{\top}|\bm{\Omega}_{n}]+\mathrm{Var}[\bm{F}_{t}|\bm{\Omega}_{n}]\right)\bm{\Lambda}^{\top}\bm{W}_{t}+(\bm{I}-\bm{W}_{t})\bm{\Sigma}_{e}^{*}(\bm{I}-\bm{W}_{t})$
$\displaystyle=\bm{W}_{t}\bm{X}_{t}\bm{X}_{t}^{\top}\bm{W}_{t}-2\bm{W}_{t}\bm{X}_{t}\bm{a}_{t|n}\bm{\Lambda}^{\top}\bm{W}_{t}+\bm{W}_{t}\bm{\Lambda}\bm{S}_{t|n}\bm{\Lambda}^{\top}\bm{W}_{t}+(\bm{I}-\bm{W}_{t})\bm{\Sigma}_{e}(\bm{I}-\bm{W}_{t})$
$\displaystyle=\bm{W}_{t}\left(\bm{X}_{t}\bm{X}_{t}^{\top}-2\bm{X}_{t}\bm{a}_{t|n}\bm{\Lambda}^{\top}+\bm{\Lambda}\bm{S}_{t|n}\bm{\Lambda}^{\top}\right)\bm{W}_{t}+(\bm{I}-\bm{W}_{t})\bm{\Sigma}_{e}^{*}(\bm{I}-\bm{W}_{t})\,,$
(20)

where $\bm{\Sigma}_{\bm{\epsilon}}^{*}$ is obtained from the previous EM
iteration.

Differentiating $\mathcal{E}$ with respect to $\bm{\Lambda}$ we obtain:

$\frac{\partial\mathcal{E}}{\partial\bm{\Lambda}}=-\frac{1}{2}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\frac{\partial}{\partial\bm{\Lambda}}\sum_{t=1}^{n}\mathbb{E}\left[\bm{\epsilon}_{t}\bm{\epsilon}_{t}^{\top}|\bm{\Omega}_{n}\right]=-\frac{1}{2}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\sum_{t=1}^{n}\bm{W}_{t}\left(-2\bm{X}_{t}\bm{a}_{t|n}+2\bm{\Lambda}\bm{S}_{t|n}\right)\bm{W}_{t}\,.$
(21)

Setting this equal to 0, dividing through by
$\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}$, and using the identity
$\mathrm{vec}(ABC)=(C^{\top}\otimes A)\,\mathrm{vec}(B)$, for matrices $A,B$
and $C$, we obtain the result

$\mathrm{vec}(\hat{\bm{\Lambda}})=\left(\sum_{t=1}^{n}\bm{S}_{t|n}\otimes\bm{W}_{t}\right)^{-1}\,\mathrm{vec}\left(\sum_{t=1}^{n}\bm{W}_{t}\bm{X}_{t}\bm{a}_{t|n}^{\top}\right)\,.$

Differentiating $\mathcal{E}$ with respect to
$\bm{\Sigma}_{\bm{\epsilon}}^{-1}$ we obtain:

$\frac{\partial\mathcal{E}}{\partial\bm{\Sigma}_{\bm{\epsilon}}^{-1}}=\frac{n}{2}\bm{\Sigma}_{\bm{\epsilon}}-\frac{1}{2}\sum_{t=1}^{n}\bm{W}_{t}\left(\bm{X}_{t}\bm{X}_{t}^{\top}-2\bm{X}_{t}\bm{a}_{t|n}\bm{\Lambda}^{\top}+\bm{\Lambda}\bm{S}_{t|n}\bm{\Lambda}^{\top}\right)\bm{W}_{t}+(\bm{I}-\bm{W}_{t})\bm{\Sigma}_{e}^{*}(\bm{I}-\bm{W}_{t})\,.$

Setting this equal to 0 and subbing in $\hat{\bm{\Lambda}}$ we obtain the
result

$\hat{\bm{\Sigma}}_{\bm{\epsilon}}=\frac{1}{n}\sum_{t=1}^{n}\mathrm{diag}\Bigg{[}\bm{W}_{t}\bigg{(}\bm{X}_{t}\bm{X}_{t}^{\top}\
-2\bm{X}_{t}\bm{a}_{t|n}^{\top}\hat{\bm{\Lambda}}^{\top}+\hat{\bm{\Lambda}}\bm{S}_{t|n}\hat{\bm{\Lambda}}^{\top}\bigg{)}\bm{W}_{t}+(\bm{I}-\bm{W}_{t})\bm{\Sigma}_{\bm{\epsilon}}^{*}(\bm{I}-\bm{W}_{t})\Bigg{]}\,.$

##  Appendix C ADMM Algorithm Derivations for Section 2.4

For the ADMM algorithm of Section 2.3, used to incorporate $\ell_{1}$
regularisation into the estimation of the loadings parameter $\bm{\Lambda}$,
we consider the penalised and augmented Lagrangian

$\mathcal{C}(\bm{\Lambda},\bm{Z},\bm{U}):=-\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]+\alpha\|\bm{Z}\|_{1}+\frac{\nu}{2}\|\bm{\Lambda}-\bm{Z}+\bm{U}\|_{F}^{2}\,,$
(22)

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. The derivation for minimising with respect to
$\bm{Z}$ and $\bm{U}$ are standard and can be found in Boyd et al. (2011). We
now go into detail on the minimisation with respect to $\bm{\Lambda}$ for
iteration $(k+1)$ of the ADMM algorithm:

$\bm{\Lambda}^{(k+1)}=\operatorname*{arg\,min}_{\bm{\Lambda}}\mathcal{C}(\bm{\Lambda},\bm{Z}^{(k)},\bm{U}^{(k)})\,.$

We can re-write the Lagrangian (22) as

$\mathcal{C}(\bm{\Lambda},\bm{Z},\bm{U})\propto-\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]+\frac{\nu}{2}\,\mathrm{tr}\left(\bm{\Lambda}\bm{\Lambda}^{\top}-2(\bm{Z}^{(k)}-\bm{u}^{(k)})\bm{\Lambda}^{\top}\right)\,,$
(23)

by dropping terms not involving $\bm{\Lambda}$ and using the property
$\mathrm{tr}(AB)=\mathrm{tr}(BA)$ for matrices $A$ and $B$. Taking the
derivative and using (21) we obtain

$\displaystyle\frac{\partial\mathcal{C}(\bm{\Lambda},\bm{Z}^{(k)},\bm{U}^{(k)})}{\partial\bm{\Lambda}}$
$\displaystyle=-\frac{\partial}{\partial\bm{\Lambda}}\mathbb{E}\left[\log\mathcal{L}(\bm{X},\bm{F};\bm{\theta})|\bm{\Omega}_{n}\right]+\frac{\nu}{2}\,\mathrm{tr}\left(2\bm{\Lambda}-2(\bm{Z}^{(k)}-\bm{u}^{(k)})\right)$
$\displaystyle=\frac{1}{2}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\sum_{t=1}^{n}\bm{W}_{t}\left(-2\bm{X}_{t}\bm{a}_{t|n}+2\bm{\Lambda}\bm{S}_{t|n}\right)\bm{W}_{t}+\frac{\nu}{2}\,\mathrm{tr}\left(2\bm{\Lambda}-2(\bm{Z}^{(k)}-\bm{u}^{(k)})\right)\,.$

Setting this equal to 0 and re-arranging for $\bm{\Lambda}$ we have

$\sum_{t=1}^{n}\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{W}_{t}\bm{\Lambda}\bm{S}_{t|n}+\nu\bm{\Lambda}=\sum_{t=1}^{n}\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{W}_{t}\bm{X}_{t}\bm{a}_{t|n}+\nu(\bm{Z}^{(k)}-\bm{U}^{(k)})\,.$

Using the identity $\mathrm{vec}(ABC)=(C^{\top}\otimes A)\,\mathrm{vec}(B)$,
for matrices $A,B$ and $C$, we obtain the result

$\mathrm{vec}(\bm{\Lambda}^{(k+1)})=\left(\sum_{t=1}^{n}\bm{S}_{t|n}\otimes\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{W}_{t}+\nu\bm{I}_{pr}\right)^{-1}\mathrm{vec}\bigg{[}\sum_{t=1}^{n}\bm{W}_{t}\bm{\Sigma}_{\bm{\epsilon}}^{-1}\bm{W}_{t}\bm{X}_{t}\bm{a}_{t|n}^{\top}+\nu(\bm{Z}^{(k)}-\bm{U}^{(k)})\bigg{]}\,.$

##  Appendix D Table to Summarise Generic Functions for S3 Object of Class
“sparseDFM”

Function |  Brief Description
---|---
`print()` |  Print out the inputs to the `sparseDFM()` call.
`summary()` |  Brief summary information on the estimation results for the factor parameters.
`plot()` |  Various types of plots available: • `factor` \- plot factor estimate series on top of the original standardized stationary data. • `loading.heatmap` \- make a heatmap of the loadings matrix. • `loading.lineplot` \- make a lineplot of variable loadings for a given factor. • `loading.grouplineplot` \- separate variable groups into colours for better visualisation. • `residual` \- boxplot or scatterplot of residuals. • `lasso.bic` \- BIC values for the LASSO tuning parameter. • `em.convergence` \- log-likelihood convergence of EM iterations.
`fitted()` |  Fitted values of the `sparseDFM()` fit, i.e. $\hat{\bm{\Lambda}}\hat{\bm{F}}$.
`residuals()/resids()` |  Residual values of the `sparseDFM()` fit, i.e. $\tilde{\bm{X}}-\hat{\bm{\Lambda}}\hat{\bm{F}}$ where $\tilde{\bm{X}}$ is the balanced (nothing-missing) data.
`predict()` |  Predict the next h steps ahead for the factor estimates and the data series. Returns S3 object of class “sparseDFM_forecast”. Use `print()` to print results.

##  Appendix E Plots for Section 4.1 - Exploring UK Inflation

Figure 5: Output of _tuneFactors()_ showing the best number of factors to use
for the inflation data set.

Figure 6: Output of _plot()_ with _type = “lasso.bic”_ showing BIC values for
$\alpha$’s considered.

Figure 7: Output of _plot()_ with _type = “factor”_ showing the 3 estimated
factors on top of the standardised stationary inflation data.

Figure 8: Output of _plot()_ with _type = “residual”_ showing the residuals
for each inflation series.

##  Appendix F Plots for Section 4.2 - Nowcasting UK Trade in Goods (Exports)

Figure 9: Missing data plot for the last 12 months of the exports full data
set.

Figure 10: Output of _tuneFactors()_ showing the best number of factors to use
for the exports data set.

##  Appendix G Code for the Pseudo Real-Time Nowcasting Exercise of Section
4.2

The expanding window nowcasting loop from August 2018 to July 2022 can be
implemented in R using:

    ## Start and end dates
    start_date <- as.Date("2018/08/01")
    end_date = as.Date("2022/07/01")

    ## Generating range of dates
    range <- seq(start_date, end_date,"months")

    ## Ragged edge variable information (match to column variables)
    # 98 series with 2 month lag (targets,IoP)
    # 333 series with 1 month lag (CPI,PPI,exchange,conf)
    # 14 series with 0 month lag (google trends)
    ragged_edge_values = c(rep(2,98), rep(1,333), rep(0,14))

    ## Remove missing jan 2004 and make time series from Feb 2004
    new_data = new_data[-1,]
    new_data = ts(new_data, start = c(2004,2), frequency = 12)

    ## Expanding window nowcasting loop and store the DFM and SDFM fits
    fit.dfm <- fit.sdfm <- list()

    for(t in 1:length(range)){

      date = range[t]
      date.year = as.numeric(format(date, "%Y"))
      date.month = as.numeric(format(date, "%m"))

      Y = window(new_data, end = c(date.year, date.month))
      Y = raggedEdge(Y, ragged_edge_values)

      fit.dfm[[t]] <- sparseDFM(Y, r = 4, alg = ’EM’)
      fit.sdfm[[t]] <- sparseDFM(Y, r = 4, q = 9, alg = ’EM-sparse’)

    }

Using the DFM fits, we can extract nowcasts and then obtain the mean absolute
error values in Table 1 using:

    ## Get the original data in levels for the targets (with Jan 2004 removed)

    data_exports = data[-1,]
    data_exports = data_exports[,1:9]

    ## Compute errors for the fits at horizon 1 and 2 for each nowcasting window

    # n = 175 represents August 2018 (for the 1st nowcasting window)
    n = 175

    err1.dfm = err2.dfm = err1.sdfm = err2.sdfm = c()

    for(t in 1:48){

      # extract nowcasts for horizon 1 and 2 of the targets
      hor1.dfm = fit.dfm[[t]]$data$fitted.unscaled[n-1,1:9]
      hor2.dfm = fit.dfm[[t]]$data$fitted.unscaled[n,1:9]
      hor1.sdfm = fit.sdfm[[t]]$data$fitted.unscaled[n-1,1:9]
      hor2.sdfm = fit.sdfm[[t]]$data$fitted.unscaled[n,1:9]

      # true data
      data.obs = data_exports[n-2,]

      # nowcasts
      for1.dfm = data.obs + hor1.dfm
      for2.dfm = for1.dfm + hor2.dfm
      for1.sdfm = data.obs + hor1.sdfm
      for2.sdfm = for1.sdfm + hor2.sdfm

      # mean absolute error
      err1.dfm[t] = mean(abs(for1.dfm - data_exports[n-1,]))
      err2.dfm[t] = mean(abs(for2.dfm - data_exports[n,]))
      err1.sdfm[t] = mean(abs(for1.sdfm - data_exports[n-1,]))
      err2.sdfm[t] = mean(abs(for2.sdfm - data_exports[n,]))

      # move window along a month
      n = n + 1
    }

    ## Print out the results
    # DFM horizon 1
    quantile(err1.dfm)
    #>        0%       25%       50%       75%      100%
    #>   80.7734  152.4034  232.8558  455.0769 1990.3548

    # Sparse DFM horizon 1
    quantile(err1.sdfm)
    #>         0%        25%        50%        75%       100%
    #>   82.47065  141.11771  212.39231  306.92548 2030.41868

    # DFM horizon 2
    quantile(err2.dfm)
    #>         0%        25%        50%        75%       100%
    #>   82.87274  201.82086  302.67172  515.08069 2162.27714

    # Sparse DFM horizon 2
    quantile(err2.sdfm)
    #>        0%       25%       50%       75%      100%
    #>   84.0062  161.9366  233.7898  366.3923 2129.1993

    ## Average

    # DFM horizon 1
    mean(err1.dfm)
    #> [1] 373.72

    # Sparse DFM horizon 1
    mean(err1.sdfm)
    #> [1] 297.2334

    # DFM horizon 2
    mean(err2.dfm)
    #> [1] 437.1499

    # Sparse DFM horizon 2
    mean(err2.sdfm)
    #> [1] 357.6024

##  Appendix H Timing Comparison Between Univariate and Multivariate KFS
Equations

It is recommended to use `kalman = "univariate"` when setting `err = "IID"`
and `kalman = "multivariate"` when `err = "AR1"` for the fastest results. In
Section 2.2, it was shown how the univariate treatment of the classic
multivariate KFS equations (Koopman and Durbin, 2000) can lead to substantial
computational gains over the classic multivariate approach (Shumway and
Stoffer, 1982) since matrix inversions become scalar divisions. The only
matrix inversion that is required is in the update of the lagged cross-
covariance matrix $\bm{P}_{t-1,t|n}$ that uses De Jong and Mackinnon (1988)
theorem:

$\bm{P}_{t,t-1|n}=\bm{P}_{t|n}(\bm{P}_{t|t-1})^{-1}\bm{A}\bm{P}_{t-1|t-1}.$

The matrix needing to be inverted, $\bm{P}_{t|t-1}$, is of dimension equal to
the number of states. When idiosyncratic errors are assumed to be IID, the
number of states is just equal to the number of factors $r$, a small number,
and hence matrix inversion is rapid. However, in the AR(1) case, the errors
augment the state variable and the number of states becomes $p+r$; leading to
a matrix that is much slower to invert.

To test the timing comparisons we simulate data from a DFM with
$\bm{\Lambda}=N(\bm{0}_{p},\bm{I}_{p})$, $\bm{\Sigma}_{\epsilon}=\bm{I}_{p}$,
$\bm{A}=0.8\times\bm{I}_{r}$ and $\bm{\Sigma}_{u}=(1-0.8^{2})\times\bm{I}_{r}$
using 2 factors and fit a DFM using `alg = "EM"` (Bańbura and Modugno, 2014)
comparing the multivariate and univariate KFS methods when errors are assumed
to be IID and AR(1). Figure 11 shows the average time from 10 experiments with
fixed number of observations, $n=100$, and varying dimension $p$, with the
`err = "IID"` case on the left and the `err = "AR1"` case on the right. The
IID case shows the clear computational gain of the univariate filtering KFS
approach, and this becomes even more prominent for larger $p$. The AR(1)
setting clearly takes a longer amount of time to run and the univariate
filtering is slower than the classic multivariate approach.

Due to the significant time difference between assuming IID errors compared
with AR(1) errors, it is recommended to assume IID errors and then perform
residual analysis afterwards for forecasting purposes.

Figure 11: Average time from 10 experiments the EM algorithm (setting _alg =
”EM”_) takes with the multivariate and univariate KFS methods with _err =
”IID”_ on the left and _err = ”AR1”_ on the right.