Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,897 characters · 19 sections · 109 citation commands
Deep Dynamic Factor Models
\affil[1]{\it Independent Researcher} \affil[2]{\it CREST \`Ecole Polytechnique, University of Warwick, CEPR, OFCE-SciencesPo}
\cleanlookdateon
\begingroup \footnote{ This work reflects the analysis and views of the authors, Paolo Andreini, Cosimo Izzo and Giovanni Ricco. No reader should interpret this work to present the views of any third party. Assumptions, opinions, views and estimates constitute the authors' judgment as of the date given and are subject to change without notice and without duty to update.
The replication code for the simulations of this paper is available on the \href{https://github.com/cosimoizzo/DDFM}{GitHub repository}.
We are grateful to Matteo Barigozzi, Antonello D'Agostino, Aldo Lipani, Massimiliano Marcellino, Jasper McMahon, Francesca Medda, Ramin Okhrati, Filippo Pellegrino, Ivan Petrella, Giuseppe Ragusa and Lucrezia Reichlin and to the conference participants of the the 2nd Vienna Workshop on Economic Forecasting 2020, the 40th International Symposium on Forecasting, the 18th Real-Time Data Analysis Methods and Applications conference, and the CFE-CMStatistics 2022 conference for many helpful suggestions and comments. } \addtocounter{footnote}{-1} \endgroup
\thispagestyle{empty}
Keywords: Machine Learning, Deep Learning, Autoencoders, Real-Time data, Time-Series, Forecasting, Nowcasting, Latent Component Models, Factor Models. \\
JEL classification: C22, C52, C53, C55.
An overarching idea in macroeconomics, already shaping the work of RePEc:nbr:nberbk:burn46-1, is that a few common forces can explain the joint dynamics of many macroeconomics variables. This stylised view of the economic data generating process has long informed the effort of economic modelling -- for example, in the Real Business Cycle (RBC) and Dynamic Stochastic General Equilibrium (DSGE) literature -- and is one of the very few robust facts in empirical macroeconomics, motivating the use of factor models STOCK2016415.
In macroeconometrics, factor models were firstly introduced by geweke1977dfm and sargentsims1977dfm and are a very early instance of `big data' in macroeconomics. Dynamic Factor Models (DFMs) deal with a large cross-section of data (`large N problem') by applying a linear dynamic latent state framework to the analysis of economic time-series. The underlying assumption of these models is that there is a small number of pervasive unobserved common factors that stir the economy and inform the comovements of hundreds of economic variables. Economic times series are also possibly affected by variable-specific (idiosyncratic) disturbances. These idiosyncratic disturbances can be due to either measurement error or variable-specific disturbances. Dynamic factor models are workhorse models in macroeconometrics and a large body of empirical evidence has found that, in many applications, a small number of factors -- as many as two -- can capture a dominant share of the variance of all the key macroeconomic and financial variables.\footnote{This family of models has been applied intensively in econometrics to different problems such as forecasting, structural analysis and the construction of economics activity indicators (see, among many others, stockwatson2002a, stockwatson2002b, forni2001generalized, forni2000generalized, forni2005generalized, forni2015dynamic, forni2018dynamic, altissimo2010new).}
Factor models are robust and flexible models, also able to accommodate for missing observations, jagged patterns of data and mixed frequencies.\footnote{Jagged edges arise when there is a varying number of missing observations at the end of multiple time-series due to non-synchronous release dates.} However, two of their important limitations are (i) the almost always assumed linear structure, and (ii) the limited scalability of these models due to the computational challenges that are encountered when estimating factors models with more than a few dozens of variables.
This paper introduces a generalisation of factor models in a deep learning framework -- which we label Deep Dynamic Factor Models (D$^2$FMs) -- that deals effectively with these challenges, while maintaining the same degree of flexibility and of interpretability of a standard DFM. Indeed, our deep learning model can `encode' the information about the state of the economy, as available in real-time, from hundreds of macroeconomic and financial variables at mixed frequency and with `jagged edges', into a handful of unobserved latent states. While similar in spirit to traditional dynamic factor models, differently from those, our model allows for non-linearities both in the encoding -- from variables to factors -- and in the decoding map -- back to the variables from the factors. We also discuss how to generalise it further to nonlinear factor dynamics.
The backbone of our modelling approach is provided by a dynamic autoencoder structure capturing the common information across variables, and whose parameters are estimated via gradient-based backpropagation. An autoencoder is a type of unsupervised learner that maps a number of variables (`input layer') into themselves (`output layer') by first `encoding' the variables' common information into a lower dimensional `code' (viz. non-linear factors), and then `decoding' it. Autoencoders are formed by a series of internal (hidden) layers each formed by a number of nodes (neurons). The encoding happens in the first half of the model, and it is the results of a series of non-linear transformations of linear combinations of inputs coming from the previous layer to the current one. Each neuron in each layer operates one of these transformations. The sequence of layers provides `depth' to the neural net, while the number of neurons per layer provides `width'. Autoencoders can be thought of as a nonlinear generalisation of principal component analysis hinton2006reducing, which are able to transform very high-dimensional data into low-dimensional factors without having to assume a linear factor structure.\footnote{As discussed in baldi1989neural, affine decoder and encoder without any nonlinearity and squared error loss will recover the same subspace of PCA. Moreover, when nonlinearity is added into the encoding network, PCA appears as one of the many possible representations bourlard1988auto, japkowicz2000nonlinear.}
The central methodological contribution of our paper is to show how to embed autoencoders in a dynamic nonlinear factor model structure with idiosyncratic components to tackle macroeconomic problems, thus providing generalisation to linear factor models. The equivalence between maximum likelihood estimation and minimisation of mean squared error together with the Universal Approximation Theorem (cybenko1989approximation, HORNIK1989359,HORNIK1990551) allow to reinterpret the D$^2$FMs and the procedure adopted in estimating them as an efficient computational method to approximate the maximum likelihood estimates of a general nonlinear factor models. A second contribution of this paper is to show how to incorporate in this framework general patterns of missing data, jagged edges and mixed frequencies, by extending gradient-based backpropagation methods for autoencoders.
Our methodology is computationally efficient and provides large gains in terms of computational speed when compared to maximum likelihood methods for DFMs. At the best of our knowledge, this paper is the first to adopt an autoencoder structure in a dynamic model with both factor dynamics and dynamic idiosyncratic components, in a state-space framework for real-time high dimensional mixed frequencies time-series data with arbitrary patterns of missing observations.
The proposed D$^2$FM framework is very general and can be, in principle, applied to many different problems both in forecasting and in structural analysis, as done with DFMs. Indeed, the model is designed to be in spirit as close as possible to DFMs.
We test the performances of the D$^2$FM both in a controlled environment through Monte Carlo experiments, and empirically on US data in a forecasting/nowcasting exercise in the spirit of giannone2008nowcasting. The Monte Carlo experiments show that the D$^2$FM largely outperforms the state-of-the-art DFM when the true data generating process is nonlinear, and offers similar performances in the linear case.
In the empirical application, we employ the D$^2$FM to encode a full real-time version of the mccracken2016fred's FRED-MD `big data' dataset for the US. The specification of the model and its hyperparameters are not fixed ex-ante but are instead selected in an intensive cross-validation exercise. The model is then evaluated in real-time by comparing its forecasting, nowcasting and backcasting performances against two benchmarks: (i) a univariate AR(1) model; and (ii) a state-of-the-art DFM estimated using a quasi maximum likelihood and an Expectation-Maximisation algorithm (see giannone2008nowcasting and banbura2014maximum). The D$^2$FM outperforms the two benchmarks, in forecasting and in nowcasting, with gains of up to around a 14% improvement when measured in terms of the root mean square forecast errors (RMSFE).\newline
The paper is organised as follows. The remainder of this section discusses the related literature. Section (ref) provides some background and the core intuition guiding our methodology: the idea that autoencoders can be seen as static generalisations of PCA, and hence that dynamic versions of these models should be seen as nonlinear generalisations of linear dynamic factor models. Section (ref) presents the methodology proposed and discuss its estimation. Section (ref) illustrates the Monte Carlo experiments. In section (ref) we discuss how to deal with the specificities of economic data. Section (ref) describes the empirical application which has the aim to track in real-time the US GDP, using real-time data. Section (ref) summarises the main results of the paper and sketches some possible future path of research. Additional technical details and a data description are provided in the Appendix.\newline
{\bf Related Literature} This work connects three distinct areas of research: the econometric literature focused on dynamic factor model estimation, research on non-linearity in empirical macro models, and deep learning research on autoencoders and latent factor models for analysing time series data.
The key problem in the factor model literature is that, due to the latency of the factors, maximum likelihood estimators cannot be derived explicitly. geweke1977dfm and sargentsims1977dfm proposed optimised algorithms for small models in the frequency domain, while engle1981one and stock1989new offered solutions in the time domain.\footnote{Specifically, engle1981one estimate the dynamic factor model using a state-space representation in which they apply the Kalman filter to compute the likelihood used for the full maximum likelihood estimation of the parameters. watson1983alternative and shumway1982approach adapt the Expectation Maximisation (EM) algorithm of dempster1977maximum for state-space representation allowing the presence of missing data, but only in the specific case where the matrix of the measurement equation is known.} The common drawback of all these proposed methods is that, in general, the maximum likelihood approach is unfeasible for datasets where the cross-section size is large. To solve this problem in forni1998let, stockwatson2002a and giannone2008nowcasting have proposed non-parametric methods based on the principal component analysis to estimate the latent components with large cross-section data.\footnote{doz2012quasi and barigozzi2019quasi have shown that when the size of the cross-section tends to infinite the estimates obtained by a quasi-maximum likelihood approach are consistent, also when there is a weak cross-sectional correlation in the idiosyncratic components.} The intuition driving this literature and its applications to macroeconomic problems inform our approach, which can be viewed as providing both a generalisation of standard dynamic factor models, and an efficient numerical approach for handling big data.
Our work also connects to the econometric literature that has explored the extent of nonlinearities in macroeconomic data and proposed univariate and multivariate nonlinear time-series models (kock2011forecasting). An important part of this literature has estimated dynamic latent models with nonlinearities, which are explicitly modelled through structural breaks, Markov switching regression or threshold regression (see barnett2016real, camacho2012markov, marcellino2010factor, korobilis2006forecast and nakajima2013bayesian).\footnote{Nonlinearities have been also modelled in structural factor models using DSGE models (Dynamic Stochastic General Equilibrium models) as in amisano2011exact to detect regime switching in volatility.} In this literature, the approach of BAI2008304 is the closest in spirit to ours and an important early effort at including nonlinearities in factor models. In that paper, either principal components of nonlinear transformation of the data are estimated or nonlinear transformation of the factors are added to a linear factor model. Our methodology is more general. In fact, differently from the procedure of BAI2008304 our D$^2$FMs needn't to assume a specific form of nonlinearity either in the encoding or in the decoding map.\footnote{The deep learning literature refers to this as a shift from `feature engineering' to `architecture engineering' PyTorch2019.}
Finally, we connect with the strand of the deep learning literature that has explored different approaches to introduce dynamics in autoencoders. The early work of pmlr-v32-gregor14 proposed the Deep AutoRegressive networks (DARN) in which hidden layers are equipped with autoregressive connections, allowing for dynamics in an autoencoder setting. Successively, Temporal Difference Variation Autoencoders (TD-VAEs) were introduced by gregor2018temporal to model dynamics in autoencoders via long short-term memory networks (LSTM) connections between belief distributions at two distant time steps. A different line of research in this literature has used Deep Learning for State Space models. For example, Krishnan2017 adopt MLPs to estimate the mean and covariance matrix of a state space with Gaussian transition dynamics. In fraccaro2017disentangled a Kalman Variational Autoencoder (K-VAE) is introduced to estimate (locally) linear Gaussian state space models (LGSSM) by disentangling the observations and the latent dynamics.\footnote{NIPS2016_6379 uses Structured Variational Autoencoders (SVAEs) to provide conjugate graphical models} Closer to our approach are the Deep State Space Models (DSSM) of rangapuram2018deep that directly estimate all the state space parameters using a Recurrent Neural Network (RNN) structure that, given the input features, provides both the latent states and all the time-varying parameters of the state space model. This modelling approach can handle both the presence of noise and missing data. While using deep learning techniques to capture salient features of the data generating process, we tackle those two issues differently. For the former, we apply a denoising approach. For the latter, selection matrices are used to mask missing data when computing the objective to be optimised, while the generative spirit of our model allows to do sampling conditional on the data to fill missing inputs. Importantly, in all the aforementioned approaches, variable specific idiosyncratic components and mixed frequencies are not taken into consideration. We propose a framework able to deal with those two additional data issues. Specifically, we include restriction matrices to aggregate the high frequency latent states to the low frequency as in mariano2003new to deal with mixed-frequency, and we design a sequential and iterative (alternated) optimisation scheme between common factors and idiosyncratic components based on a Markov Chain Monte Carlo algorithm.
More generally, deep learning methods have seen early applications in Finance and Economics. In Finance, they have been employed to predict asset prices, stock returns or commodity prices sezer2019financial. Closer to our approach, a few recent works have applied neural net to macroeconomic questions. cook2017macroeconomic employed a number of neural network architectures, including also autoencoder, to forecast the US unemployment rate. loermann2019nowcasting proposed a neural net model to predict the US GDP. holopainen2017toward proposed an horse race among different machine learning methods and showed that such models are able to outperform conventional statistical approaches in predicting crisis periods. Finally, heaton2016deep, gu2018empirical,gu2019autoencoder employ rich datasets incorporating both stock data and macroeconomic aggregates to predict stock returns.\\
Dynamic factor models for econometric times series are multivariate probabilistic models in which a vector of stochastic disturbances are transmitted via linear dynamic equations to the observed variables. They assume that a small number of stochastic unobserved common factors informs the comovements of hundreds of economic variables. In doing so, they combine two core ideas of macroeconomics: the Frisch-Slutsky paradigm that assumes the economic variables to be generated by the stochastic components (the economic shocks) via usually linear dynamic difference equations; and the idea that has guided macro since RePEc:nbr:nberbk:burn46-1 that a few common disturbances explain most of the dynamics of all the macroeconomics variables, with a residual share due to idiosyncratic components. DFMs are similar in intuition to principal component analysis (PCA) but assume stochastic and dynamic structure that allows their application to econometric time series.
Autoencoders (AE) are neural networks trained to map a set of variables into themselves, by first coding the input into a lower dimensional (or undercomplete) representation) and then decoding it back into itself. The lower dimensional representation forces the autoencoder to capture the most salient features of the data. In constructing a nonlinear reflexive map that links the inputs back to itself via a lower dimensional space, autoencoders can be thought of as a nonlinear generalisation of PCA.\footnote{The connection between PCA and Autoencoders is discussed in goodfellow2016deep and in the references therein.}
In this section, we explore the deep connection between factor models and autoencoders to show that a dynamic formulation of autoencoders can be thought of as a nonlinear generalisations of dynamic factor models, in the same way in which standard autoencoders can be seen as generalisations of principal component analysis.
Let us first introduce a general formulation of latent factor models with idiosyncratic components. We define $\boldsymbol{y}_t = (y_{t,1},...,y_{t,n})$ as the vector collecting the $n$ variables of interest at time $t$, usually assumed to be the realisation of a vector stochastic process. A very general latent factor model can be written as
where $\boldsymbol{f}_t$ is an $r \times 1$ (for $r=dim(\boldsymbol{f})<<n=dim(\boldsymbol{y})$) vector of latent common stochastic components -- i.e. the factors --, $\boldsymbol{\varepsilon}_{t}$ are idiosyncratic and unobserved stochastic error terms, and $F(\cdot)$ is a generic function mapping the unobserved factors into the observed variable. Usual assumptions are that $\boldsymbol{f}_t$ and $\boldsymbol{\varepsilon}_{t}$ are independent, with zero mean and finite variance (the variance of $\boldsymbol{f}_t$ is often assumed to be a diagonal matrix). For later reference, we indicate as $\widetilde{\boldsymbol{y}}_t$ the component of $\boldsymbol{y}_t$ that relates to the common factors. By assuming also a linear function $F(\cdot)$, the model reduces to the standard linear factor model
However, in general, $F(\cdot)$ needn't be linear and we can express the factor component of the model as
where $G(\cdot)$ is the function mapping the observables into the `code' $\boldsymbol{f}_{t}$ (encoding function), and $F(\cdot)$ is the function mapping the factors back into $\boldsymbol{y}_t$ (decoding function). In this form, the connection between factor models and autoencoders is more evident. In fact, the map in Equation (ref) can be seen as a very general autoencoder. Linear factor models can be seen as a special case of factor models assuming both a linear encoding and a linear decoding function. It is worth observing that the model in Equations (ref) and (ref), without specifying dynamic equations for the stochastic components, can be seen as purely static model.
Autoencoders belong to the deep neural net family of models and have been introduced for applications involving dimensionality reduction (lecun1987phd, bourlard1988auto, hinton1994autoencoders, hinton2006reducing). Autoencoders solve the parametric problem of finding a mapping (or `learning a representation' in the DNN jargon) of the form $\widetilde{\boldsymbol{y}}_t = F(G(\boldsymbol{y}_t))$ under the constraint of minimising a loss function of choice
where $ \mathcal{L}(\cdot)$ is the loss function and $ \boldsymbol{\theta}$ is the vector collecting all the parameters in $G(\cdot)$ and $F(\cdot)$.
The Principal Components Analysis (PCA) can be seen as the autoencoder minimising the square loss function
and assuming both a linear coding and a linear decoding function, i.e. $\boldsymbol{f}_{t} = G(\boldsymbol{y}_{t}) = W'\boldsymbol{y}_t $ and $\widetilde{\boldsymbol{y}}_t = F(\boldsymbol{f}_{t}) = \boldsymbol{\Lambda} \boldsymbol{f}_{t}$. Figure (ref) provides a neural net representation of PCA.
In principle, $G(\cdot)$ and $F(\cdot)$ can be any nonlinear function and hence finding the correct functional form capturing a data generating process of interest can be a daunting problem. Autoencoders provide a practical implementation of this problem by expressing the composition of two functions as a chain of two multilayer perceptrons (MLPs): the first chain operates the coding, while the second produces the decoded output (see a graphical representation of a symmetric autoencoder in Figure (ref)). A multilayer perceptron is a type of feedforward artificial neural network composed of a number of `hidden' layers each formed by a number of `nodes' (or `neurons'). Each neuron in each layer receives some inputs from the neurons in the previous layer and outputs to the next layer an activation output, $h_{m_l}^l$. The activation function (or `link function') of each neuron is a nonlinear function parametrised as
where $l$ is the layer (for $l = 1, \dots, L$), $m_l$ is the node, and $\boldsymbol{\theta}_{m_l}^l \equiv \{\boldsymbol{W}_{m_l}^l, b_{m_l}^l\}$ are the parameters of the activation function to be determined, in the form of, respectively, a set of weights and a constant (also called bias). If $l=1$ then, $\boldsymbol{h}_{l-1} = \boldsymbol{h}_0 = \boldsymbol{y}_t$ that is the input data vector. Common choices for the activation function $g_{m_l}^l(\cdot)$ are the sigmoid, the hyperbolic tangent (tanh), softplus, and the rectified linear unit (ReLu) functions. Hence, in other words, the neuron in each upper layer is the product of an element wise (usually monotone) transformation applied on an affine mapping of the neurons in the lower layer.
We can define as $\boldsymbol{g_l}(\cdot)$ the vector containing all the activation functions of a layer $\{g_{1}^l(\cdot) \dots g_{M}^l(\cdot)\}'$ for all the nodes $1, \dots M$. Hence the first MLP will be given by the composition of the activation functions of each node in each layer of the encoding feedforward network, i.e.
In a similar way it is usually also defined the MLP operating as decoding network, i.e. as a sequence of layer each containing neurons operating as activation functions over a weighted sum of the inputs plus a constant, i.e.
where $\tilde g_{L'}$ is the vector of link functions, and $L'$ is the number of hidden layers in the decoding network.
While the functional form adopted for the activation functions may seem arbitrary, yet a network with such a structure can approximate any nonlinear continuous function. In fact, the Universal Approximation Theorem, a key result in the neural net literature, states that a feed-forward network even with a single hidden layer, containing a finite number of neurons, can approximate continuous functions on compact subsets of $\mathbb R^n$, under mild assumptions on the activation function. However, it does not guarantee that the algorithm adopted to estimate the network can learn the correct parameters (see cybenko1989approximation, HORNIK1989359,HORNIK1990551 and lu2017expressive).
The autoencoder is said to be symmetric when the number of hidden layers in the encoding and in the decoding networks are the same (i.e. $L'=L$), otherwise asymmetric. Asymmetric autoencoders usually have several layers of encoding but only a single layer of decoding (i.e. $L'=1$), hence the decoding network is not a MLP but an SLP (single layer perceptron).\footnote{Asymmetric autoencoders were introduced by majumdar2017asymmetric and have been found to be more accurate compared to traditional symmetrically autoencoders for classification accuracy and also to yield slightly better results on compression problems (over the following datasets: MNIST, CIFAR-10, SVHN and CIFAR-100). majumdar2017asymmetric argue that the asymmetric structure helps to reduce the number of parameters to estimate and hence the potential extent of overfitting.}
For a given choice of the link functions, the parameter vector of the autoencoder $ \boldsymbol{\theta}$ contains the full set of weights and constants (biases) that define the affine transformations operated by each neuron before the link function is applied. These parameters are determined by minimising the loss function $\mathcal{L}(\boldsymbol{y}_t,\widetilde{\boldsymbol{y}}_t; \boldsymbol{\theta}) = \mathcal{L}(\boldsymbol{y}_t, F(G(\boldsymbol{y}_t)); \boldsymbol{\theta})$, via back-propagation (rumelhart1986learning and lecun1987phd).\footnote{Stochastic gradient descent (SGD) algorithms, proposed by kiefer1952stochastic, are commonly adopted in the deep learning literature, and update the gradient of each parameter using only randomly selected subsamples of the training dataset. These subsamples are called `minibatches' and they are equal partition of the original training datasets. The computational cost of SGD algorithms is independent with respect to the sample size. All the optimisation algorithms tend to analyse the training dataset multiple times in order to reach a better estimation of the parameters that relies less on the starting point. A run of the algorithm over the entire dataset is called `epoch'. Common optimisation algorithms are momentum algorithms, as AdaGrad by duchi2011adaptive, RMSProp and its variation by tieleman2012rmsprop, and ADAM by Adam, that rely all on the well know gradient descent algorithm by cauchy1847methode.} One way to estimate autoencoders is by corrupting the inputs with noise injection vincent2008extracting during the training process. Those are referred as Denoising Autoencoders. The intuition for this procedure is that, as discussed in vincent2008extracting and bengio2013generalized, it forces the model to learn the data distribution and not only the distribution specific of the sample used, thanks to data augmentation. Also, noise injection can be seen as a procedure to improve the robustness of neural networks.
So far we have discussed the general structure of factor models by abstracting from the dynamics and focusing on the `static' map into lower dimensional factors. Dynamics is usually introduced in DFMs by assuming that both $\boldsymbol{f}_t$ and $\boldsymbol{\varepsilon}_{t}$ are generated by linear stochastic vector difference equations. For example, banbura2010nowcasting and banbura2014maximum consider a system specified as
where $\boldsymbol{B}_1, \dots, \boldsymbol{B}_p$ are the $r \times r$ matrices of autoregressive coefficients for the factor and $\boldsymbol{\Phi}_1, \dots, \boldsymbol{\Phi}_d$ are the $n \times n$ diagonal matrices of autoregressive coefficients for the idiosyncratic component (i.e. $\boldsymbol{\Phi}_1 = diag(\phi_1, \dots, \phi_n)$). Specifically, banbura2010nowcasting and banbura2014maximum assume a VAR process of order two ($p= 2$) for factors, and of order one ($d= 1$) for the idiosyncratic component.\footnote{The zero cross-correlation at all leads and lags of the idiosyncratic components has been shown by doz2012quasi and barigozzi2019quasi to be asymptotically valid even when it is mildly violated in small sample.}
Such a structure can be seen, in our framework, as obtained by assuming that:
Autoencoders provide a practical solution to estimate factor models with a more general structure, potentially relaxing one or all of these assumptions to obtain both nonlinear maps from reduced dimension factors to variables and vice-versa, but also to introduce nonlinear dynamic equations. This approach to the generalisation of dynamic factor models, is what we call Deep Dynamic Factor Models (D$^2$FMs). In the next section, we show how to construct and estimate an autoencoder that relaxes assumptions A.1 and A.2, while maintaining the others.\footnote{See doz2012quasi for the robusteness of their DFM estimation procedure with respect to assumptions A.4 and A.5.}
In principle the parameters of a parametric factor model of the form $\boldsymbol{y}_t = F(\boldsymbol{f}_{t}) + \boldsymbol{\varepsilon}_{t} $ would be estimated via maximum likelihood,
where by $\boldsymbol{Y}$ is the full sample of observation, and $p_{model}(\cdot; \cdot) $ is the conditional probability density function of the model.
However, a direct maximum likelihood is rarely feasible, even for linear models, and iterative methods to find maximum likelihood or maximum a posteriori (MAP) estimates of the parameters of the model are preferred. In fact, maximum likelihood estimators of the parameters $\boldsymbol{\theta} = (\boldsymbol{\Lambda}, \boldsymbol{B}(L), \boldsymbol{U}, \boldsymbol{\Phi}(L), \boldsymbol{Q})$ are in general not available in closed form and a direct numerical maximisation can be too demanding when $n$ is large. Indeed, a proposed solution in the linear factor model literature is to adopt the Expectation Maximisation (EM) algorithm, a maximum a posteriori method, and to initialise the common factors $\boldsymbol{f}_t$ with PCA on the observables.\footnote{Estimation of linear factor models was originally carried out via simple principal component analysis (PCA).} The updates of the latent components are performed using the Kalman filter and smoother.
A similar approach can be adopted from a `deep' point of view on factor models by employing the methodologies developed in the deep learning literature, without losing the dynamic model interpretation. As we discuss in the next section, the model parameters of a D$^2$FM can be estimated via Monte Carlo gradient methods, instead of using the EM algorithm. This has computational advantages -- the methods are fast and reliable -- even when the dataset is big. At the same time, estimation results can be thought of as approximating a maximum a posteriori method.
It is well know, in the linear case, that if the innovation $\boldsymbol{\varepsilon}_{t}$ are assumed to be independent (or uncorrelated) of $\boldsymbol{f}_{t}$ and normally distributed, than the maximisation of the likelihood with respect to the parameters of the model yields the same estimate for the parameters as does minimising the mean squared error.
Importantly, this equivalence between maximum likelihood estimation and minimisation of mean squared error holds regardless of the function used to predict the mean of the Gaussian distributed variable $\boldsymbol{y}_t$ goodfellow2016deep. This result allows for an interpretation of estimation results from autoencoders with mean squared error from a Bayesian perspective, using standard likelihood methods, or from a frequentist one as the (approximated) mean estimator of a Gaussian distributed process.
Furthermore, the equivalence between maximum likelihood estimation and minimisation of mean squared error together with the Universal Approximation Theorem allow to reinterpret the autoencoders and the procedure adopted in estimating them as an efficient computational method to approximate the maximum likelihood estimates of nonlinear factor models. These are dynamic models that are defined by a conditionally Gaussian distribution centred around a mean provided by a nonlinear but continuous function of the inputs. In the next section, we provide an algorithm that implements Deep Dynamic Factor Models.
In its general form, the D$^2$FM can be written as
where the assumptions on the linearity of the dynamic equations are maintained, Equations (ref) and (ref), while the model allows for a nonlinear map between variables and factors.\footnote{Alternatively, the dynamic of the common latent states can be estimated directly in the algorithm (ref). This can be achieved by including an autoregressive layer before the decoding layer ($F_{\boldsymbol{\theta}_2}$ of the algorithm). This additional layer would coincide with the state equation of the common part. This layer could be linear for linear dynamics, but also nonlinear and composed of multiple layers such as multi-layer perceptrons or recurrent layers such as LSTM.} The estimation of linear factor dynamics separately creates what stock2011dynamic call a `state space with static (common) factors' as opposed to a `state space with dynamic (common) factors'.
The D$^2$FM can be implemented using a symmetric autoencoder architecture with a MLP capturing the encoding function, Equation (ref), and another MLP providing the nonlinear decoding, Equation (ref). The assumption of i.i.d. and Gaussian innovations allows for an interpretation of the estimated network as MAP of the likelihood of the model goodfellow2016deep.
Importantly, such a model specification encompasses several simplified models, most notably the standard linear DFMs, and hence the estimation algorithm can be specialised to the scope.
The core of the model is provided by an autoencoder with a nonlinear multilayer encoder and either a symmetric structure in the decoding, for nonlinear decoding, or an asymmetric structure with a linear single layer decoder. Linear stochastic autoregressive equations are adopted to model the dynamics of factors and idiosyncratic components. Figure (ref) shows a diagrammatic representation of the model with a single layer decoder.
The number of hidden layers in the encoding network as well as the number of neurons need not be pre-specified but can be selected via cross-validation. This enables the model to capture different degrees of complexity in the data, which is not know a-priori. With respect to the activation functions, we equip each neuron in the coding layers with a link function in the form of the hyperbolic tangent (tanh) for the real-time macroeconomic dataset, and of the rectified linear unit (ReLU) for the Monte Carlo exercises.\footnote{In general, some tuning might be needed in order to find the correct specification for the dataset at use.} In the encoding multilayer perceptron we also include two batch normalisation layers to induce some regularisation, and control over potential covariate shift pmlr-v37-ioffe15.\footnote{Covariate shift is a phenomenon that occurs in deep learning when the distribution of the input data changes between the training and testing phases. This can lead to a decrease in the accuracy of the model because the model has not been exposed to the new distribution during training.}
In the decoding network, an additional linear layer can be included to introduce constraints needed to account for the mixed frequencies of macroeconomic data. This additional layer does not have any additional parameter, and it only includes aggregation weights to map the high frequency latent states to the low frequency observables.
In our estimation of the D$^2$FM, we propose a two-step procedure to differentiate between on-line (out-of-sample) and off-line (in-sample) learning.\footnote{In the deep learning literature, online learning means that the model estimates parameters using the flow of data as it comes in (or using a simulated flow). Offline means that the model employs a static dataset. This is similar to the standard econometric distinction between out-of-sample and in-sample.}
Algorithm (ref) implements the off-line estimation step ({\bf Step 1}) of our D$^2$FM, assuming an AR(d) for $p_{idio}(\cdot)$, but possibly a general encoding $G_{\theta_G}(\cdot)$ and decoding $F_{\theta_F}(\cdot)$ function. The proposed algorithm for estimating D$^2$FM builds on and extends what has been proposed by bengio2013generalized to estimate Generalised Denoising Autoencoders.
Let us present the estimation algorithm. Parameters are first initialised. Line 1 performs a filtering of the input data ${\boldsymbol{y}}_{t}$ by using the conditional mean of the AR(d) of the idiosyncratic components. From lines 2 to 6, the Monte Carlo step and the gradient updates over each epoch and batch are carried out, employing the filtered data $\widetilde{\boldsymbol{y}}_{t}$ and injecting Gaussian noise from $\boldsymbol{\varepsilon}_{t}$ in a denoising fashion to obtain the noisy observations $\boldsymbol{y}^{(mc)}_{t}$. In line 7, the latent states $\boldsymbol{f}_t$ are extracted from the encoding network via Monte Carlo integration, while from line 8 to 10 the algorithm updates the parameters of the idiosyncratic process $\boldsymbol{\varepsilon}_t$, conditional on the factors and the observables. The adoption of an $L_2$ (MSE) loss function $\mathcal{\hat{L}}(\boldsymbol{y}_t,F_{\boldsymbol{\theta}_F}(\boldsymbol{f}_t | \boldsymbol{y}^{(mc)}_{t}))$ allows for interpretability of the results, as discussed. We specify an estimated loss, as missing data prevents us from deriving the exact loss.\footnote{We give details about the treatment of missing data in Section (ref).} Finally, convergence is checked as the $L_1$ norm of the distance between the loss function at two iterations. It is worth noting that the loss, $\mathcal{\hat{L}}(\cdot)$ includes only the common components, since under our assumptions at convergence we have the following decomposition of the log-likehood:
where $\hat{\boldsymbol{\varepsilon}}_t$ is the estimated idiosyncratic autoregressive component. In running over epochs and batches (lines 2 to 6), the algorithm injects uncorrelated noise into the data (it is a Denoising Autoencoder). Hence it searches for a maximum a posteriori of the parameters for the modified model with log-likelihood
where $p_{noisy}(\boldsymbol{y}^{(mc)}_{t} | \boldsymbol{y}_t, \hat{\boldsymbol{\varepsilon}}_t)$ is the corruption distribution, using a Gaussian autoregressive process. The idea behind this procedure is to filter out the foreseeable idiosyncratic part from the input variables, so that only the common component(s) remain(s). Injecting noise from the unconditional idiosyncratic distribution will generate new samples which are not unreasonably far from the old ones. In doing so, we define an appealing and convenient linkage between the corruption process of the denoising approach and the idiosyncratic component distribution ($p_{idio}$).
In {\bf Step 2}, the output of the algorithm is cast in the state-space of equations (ref)-(ref). Dynamics of the common factors are estimated via OLS or Maximum Likelihood.\footnote{The dynamics of the common latent states can also be estimated directly in Algorithm (ref). This can be achieved by including an additional autoregressive layer before the decoding network ($F_{\boldsymbol{\theta}_F}$ of the algorithm).} State updates can then be carried out via either nonlinear filtering procedures for a nonlinear decoder, or via Kalman filtering in the presence of a linear decoder. This allows for online (i.e., out-of-sample) learning with the flow of data.
The D$^2$FM described in this section is subject to the selection of a number of critical parameters determining its structure, beyond the coefficients $\boldsymbol{\theta}$. These parameters are commonly known as hyperparameters, being them set before the training starts and usually selected over a grid with respect to some validation loss, which is estimated via a process called cross-validation.
The D$^2$FM has hyperparameters typical of both deep learning and time-series models. In particular, the deep learning hyperparameters can be divided into two categories. The first relates to the neural network structure and includes: the type of layers, the number of hidden layers, the number of neurons per each hidden layer, penalisation coefficients, dropout layers and relative dropout rates (if included), batch normalisation layers and the link function used. The second category relates to the optimisation algorithm used and comprehends: size of the mini-batches, number of epochs, the learning rate and the momentum coefficients of the gradient optimisation method, if present. Standard time-series factor models have few additional hyperparameters which include: the number of latent common states, the number of lags of the input variables, the number of lags of the latent common states and of the idiosyncratic states. These hyperparameters, in the time series literature, are either fixed a-priori or estimated using information criteria instead of grid-search algorithms.\footnote{The Akaike information criteria (AIC) and the Bayesian information criteria (BIC) can be used to determine the number of lags; while the number of latent factors can be, in principle, estimated using the method proposed by alessi2010improved which improves bai2002determining's methodology.}
It is important to observe that in time-series we cannot apply the common `K-fold' cross-validation method because usually its technical conditions are not met bergmeir2018note. Therefore, we use a standard out-of-sample approach which consists in splitting the set of observations available up to a certain point in time, $T$, between a training set $[0, T-k*h-1]$, and validation set $[T-k*h, T-(k-1)*h]$, where $h$ determines the length of the set, while $k = K,\dots,1$ with $K \ll \frac{T-1}{h}$. By averaging over the losses computed on the $K$ `validation sets', we get an estimate of the `validation loss'. This means that we need to estimate a given model with fixed hyperparameter $K$ times, and this for each possible hyperparameters combination. Therefore, with deterministic search method the computational cost is exponential in the dimensionality of the hyperparameters.\footnote{Alternative methods based on stochastic search are available (see for example bergstra_tpe), but then the results could be not robust when they are computed on a small number of iterations.}
In this section, we compare the performances of DFMs and of D$^2$FM in a controlled environment, by simulating artificial time series data with Monte Carlo experiments. In doing so, we first consider the linear data generating processes studied in doz2012quasi and banbura2014maximum, and then, adopt its nonlinear generalisation proposed by gu2019autoencoder. In particular, we consider data generated the following process:
for $t=1,\dots,T$. As in gu2019autoencoder, we simulate both a linear, and nonlinear factor model, and allow $F(\cdot)$ to take two forms:
where $sgn(\cdot)$ is the sign function ($1$ if positive, $-1$ if negative) and $poly(\cdot,2)$ is a generator of polynomial functions of order $2$. The parameters of the model are set as follows:
We consider a range of possible specifications of the parameters by setting the number of factors $r=\{1, 3\}$, the number of variables $n=\{10, 100\}$, autoregressive coefficient of the factors $\rho=\{0.5, 0.9\}$, autoregressive coefficient of the idiosyncratic component $\alpha=\{0, 0.5\}$, the number of observations $T=\{50, 200\}$, and the fraction of missings is in $\{0, 0.3\}$.
For each setting, we run 100 Monte Carlo simulations and estimate a DFM, and a four layers D$^2$FM, with a ReLU activation function augmented with three BatchNorm layers. The number of factors is set for both models to the true number of factors (i.e. $r$ when the DGP is linear in the factors and $\tilde{r}$ when it is nonlinear). Starting from the factor layer, hidden neurons increase by a factor of two in each layer up to the input layer.
As in stockwatson2002a, doz2012quasi and banbura2014maximum, we compare the models based on the trace $R^2$ of the regression of the estimated factors on the true ones
where $\hat{F}=\mathbb{E}_{\boldsymbol{\hat{\theta}}}[F|H_{T}]$ and $H_{T}$ is the history of the data.
In the linear case, differences in performances between the D$^2$FM and the DFM (Table (ref)) are very small, and indicate a marginal gain for one model or the other, depending on the specific case. This points to the fact that the two framework are generally equivalent when the DGP is linear. In other words, the D$^2$FM can be seen as a computational efficient way to estimate linear factor models, conditional on a linear DGP.
In the nonlinear simulations (Table (ref)), the D$^2$FM is instead clearly superior. Strikingly, all the differences are statistically significant and the D$^2$FM can explain between $15\%$ and $34\%$ more of the total variance of the simulated common factors. These results validate the ability of the D$^2$FM to better handle several forms of nonlinearities in the DGP, as compared to a standard DFM.
An additional key advantage of the D$^2$FM are the superior performances from the point of view of its computational time. In Table (ref) we compare the computational time required to estimate a DFM versus a D$^2$FM. The table shows clear computational advantages in favour of the gradient based Monte Carlo approach of the D$^2$FM as compared to the OLS EM approach of the DFM, when the dataset features many variables (starting with $150$ observable variables, in the case of our experiments).
To apply our model to macroeconomic data in a nowcasting and forecasting exercise, we need to introduce a version of the D$^2$FM able to track and forecast developments in economic variables in real-time. We modify the model to efficiently encode mixed-frequency data with ragged edges. In fact, economic data in real time are generally not available at the same frequency -- be it weekly, monthly or quarterly --, and missing data are a feature of real-time macroeconomic datasets, due to the non-synchronous and staggered data releases of new datapoints from statistical offices.
As in the previous section, we specify a linear mapping between the factors and the variables (see Figure (ref)), i.e.
In this form, the model can be seen as a very flexible generalisation of the approach of BAI2008304 that propose to extract factors from variables as well as their squared values and their crossproducts.
There are a few advantages to considering this simpler D$^2$FM. First, the model maintains the same level of interpretability of a standard DFM, hence making it easy to compare the two models. Indeed, this simple architecture is motivated by the recent work of rudin2019stop that has encouraged the design of models that are inherently interpretable, as opposed to a purely `black box' approach. Second, while interpretable, the autoencoder structure allows us to introduce deep learning methods in this framework to test its potential, towards the construction of more general models. Third, the linear decoding network and the linear state-space framework allow to update in real-time the latent states in an interpretable and statistically grounded framework, by employ a standard Kalman filter. Finally, the adoption of linear filtering techniques, in turn, allows for an easy interpretation of the model forecast revisions coming from the flow of data onto the performances of the model, as in banbura2010nowcasting.
We deal with missing data in two or three steps depending on the dataset. If in the pre-training when dropping missing values we are left with a few number of observations, then we first initialise missing values with a spline method. Otherwise, this step is omitted and the pre-training is carried out only on non missing data points. Second, we iterate the parameters maximisation by replacing the missing data in the full sample with fitted values obtained by conditioning on the estimated model (both parameters and latent states). Maximisation is carried out only on non-missing points, therefore the number of observations over which the gradients are computed can differ across dimensions. Finally, in the real-time online update phase (i.e. the out-of-sample procedure), we employ the Kalman filter to update the missing data, therefore accommodating for ragged edges (banbura2010nowcasting,banbura2014maximum,camacho2012markov).
In dealing with mixed-frequency data there are several options (marcellino2010factor, foroni2013survey, blasques2016weighted). The most popular one, when the dataset includes monthly and quarterly variables, is the mariano2003new approximation, in which the observed quarterly variable is modelled as a partially observed monthly series. By assuming that the (log-)levels of the quarterly variable, $Y_{t}^q$, at the end of the quarter are the sum of an unobservable monthly counterpart $Y_{t}^m, Y_{t-1}^m, Y_{t-2}^m$ , and defining with the growth variable $y_{t}^q$ the quarter over quarter change, we have
where $y_{t}^m = \Delta Y_{t}^m$ denotes the unobserved month-on-month growth rate of a quarterly variable that admits the same factor representation proposed in equation (ref). In our model, this approximation is implemented by including an additional final layer to the decoding network allowing the monthly factors to be mapped into the quarterly variables. This layer has fixed weights not subject to the optimisation.
The core of the model is provided by an asymmetric autoencoder with a nonlinear multilayer encoder and a linear single layer decoding structure. Table (ref) provides a summary of the network design choices, and reports the choices operated for each hyperparameter of our model, a number of which are selected via cross-validation.
Optimisation is carried out by using ADAM Adam with default hyperparameters and 100 epochs, both during pre-training and training. Before starting the training, ADAM is reinitialised and then is run on batches (i.e. subsamples) with size of at least 100 monthly observations (approximately 8 years, the average duration of a business cycle). In the training phase we set again the number of epochs (runs on the full sample) to 100 for each iteration of the MCMC. These iterations are used also to update the idiosyncratic distribution.
In our empirical model, parameters are initialised in a two stage approach. First, by using a Xavier initialisation -- weights in the link functions are sampled from a Gaussian distribution with zero mean and a variance of $2/(n_{in} + n_{out}$), where $n_{in}$ is the number of input units and $n_{out}$ is the number of output units glorot2010understanding, and then by performing a pre-training exercise using a standard autoencoder on a full dataset where the rows that contains missing data are discarded.\footnote{In particular, in the empirical application we check that at least 50 time observations are present when applying this rule. If this is not the case, then we drop time periods that have more than 20% missing values for the corresponding features, and we fill the rest with splines (see Section (ref)).} This pre-training procedure is needed to `warm up' the chain.\footnote{Warming up a chain in deep learning refers to the process of initialising a neural network model with weights that have been pre-trained on a related task or dataset, before fine-tuning the model on the specific dataset of interest. The goal of this process is to provide the model with a good starting point for the optimisation process, as the pre-trained weights may capture useful information that can accelerate convergence and improve performance on the new task.}
We now test the performances of the model presented in the previous section in forecasting, nowcasting and backcasting using a fully real-time `big' US macro dataset, and against three benchmark models:\footnote{Backcast is the estimate of the previous quarter up to the official release date; nowcast is the estimate of the current quarter up to the official release date, and forecast is the estimate of the next quarter up to the official release date. We are able to produce backcast values because the GDP is released usually 5 weeks after the end of the reference quarter, hence we use the releases of the other variables during these 5 weeks to update the backcast figure.} (i) a univariate AR(1) statistical benchmark; and (ii) a state-of-the-art DFM with two and three latent factors, estimated via quasi maximum likelihood as proposed by giannone2008nowcasting and generalised in banbura2014maximum (we refer to this model as DFM-EM). The model is multitarget, i.e. it is optimised against all of the variables in the dataset and not only one of them, but our discussion of the results mainly focus on US GDP. This exercise can be seen as a validation test to check whether the model is able to correctly capture the relevant features of the data generating process, and to benchmark it against other state-of-the-art models.
To test its capability, we estimate the model by encoding a real-time version of the full mccracken2016fred's FRED-MD dataset, a large macroeconomic database for the US economy, specifically designed for the empirical analysis of `big data'.\footnote{We marginally extend the dataset by including two Purchasing Managers' Indices (PMIs), since they are considered to be important indicators for nowcasting and do not get revised over time.} The cross-section of data is mixed frequency because it includes 128 monthly indicators and Real GDP, that is a quarterly variable. All the data are stationarised and standardised following the specifications in mccracken2016fred.\footnote{In the Appendix Tables (ref)-(ref) provide the complete list of the variables used and their transformation code.} In Tables (ref)-(ref), we report also the respective publication delay (in days) of each series. There are substantial differences in the timeliness of different variables. Some of them are more timely (e.g. `soft' indicators or surveys), while others are released with one-two months of delay (usually `hard' data on real activity).
The vintages in the datasets span the period from January 1980 to May 2020. Figure (ref) reports the number of variables available across time periods. We first estimate the model using the data up to December 2005, and then we perform an expanding window forecasting exercise starting form the 1$^{\text{st}}$ of January 2006. Hence our test sample goes from 1$^{\text{st}}$ of January 2006 to 31$^{\text{st}}$ of May 2020, including the Great Recession in 2007-2009. A data vintage is created every time a new time-series data point is released, and it contains all the data available up to that point in time, including also data revisions. The real-time infrastructure adapts automatically to the expanding number of variables used as input for the model. For each iteration, as new data arrive, the model is re-evaluated and outputs a sequence of backcasts-nowcasts-forecasts for GDP and all the other variables. These forecasts are conditional only to the real-time information set, i.e. only data available up to that specific point in time without taking into consideration further revisions.
Many of the hyperparameters of the model are not fixed ex-ante but are instead selected using an intensive cross-validation exercise, as reported in Table (ref). The real-time cross-validation exercise also provides information on the ability of the model to update its optimal hyperparameter specification over time, as new data comes in (the validation length is set to one year). Figure (ref) shows the evolution of the number of factors and lags that are selected via cross-validation over the sample.
Figure (ref) shows the nowcast reconstruction in real-time for the D$^2$FM, the DFM-EM with 2 and 3 factors, and the AR(1) model against the realised quarterly US GDP. Overall, the D$^2$FM and the DFM-EM models provide a similar assessment of the state of the economy, in the nowcasting horizon, though the D$^2$FM is more accurate.
We formally assess the performances of the model -- and of the AR(1), the DFM-EM with 2 and 3 factors -- by computing root mean square forecast error (RMSFE). This metric is updated every time the data vintage gets updated, due to a new data release. We report both an overall RMSFE (Table (ref)) that gives us a synthetic value about the performance of each model on the entire out-of-sample set, and a dynamic RMSFE (Figure (ref)) that illustrates how the RMSFE evolves from the forecast period to the backcast period, until the day before the release. Results indicate that the D$^2$FM is able to outperform all the competitors during the entire forecast period and for most of the nowcast period. The gain in terms of performance achieved by the D$^2$FM in these two periods is quite considerable and reflects the ability of this model to better compress the relevant information in the data, thus reducing the level of uncertainty.
The model also delivers forecasts for all the variables in the model. Table (ref) reports the average of the RMSFEs of the D$^2$FM over all of the monthly variables, in ratio to the AR(1) RMSFEs. The D$^2$FM beats the AR(1) over all the horizons -- the backcast improves by $10\%$, the nowcast improves by $20\%$ and the forecast improves by $18\%$.\footnote{Overall, the D$^2$FM improves the prediction accuracy for roughly the $80\%$ of the monthly variables included in the dataset with respect to the AR(1).}
As a final exercise, we show how to build a composite indicator of the state of economy in real-time using the decoding map (or loadings). We do this by aggregating the latent states trough a weighting scheme. Specifically, we define the composite indicator as
where $f_k$ for $k=1,..., r$ are the common factors (the code) and $F_{k,i}$ is the matrix of the coefficient for the factor $k$ at variable $i$, as in the Equation (ref), while $||F||_{F}^{2} = \sum_{k=1}^{r}\sum_{i=1}^{n} F_{k,i}^{2}$ is the squared Frobenius norm of the coefficients. The sign of the indicator is fixed to have a positive correlation with GDP. Figure (ref) reports the composite indicator using the real time out of sample exercise, that is shown to track well the developments in the US economy.
The central contribution of our paper is to introduce Deep Dynamic Factor Models (D$^2$FMs) by showing how to embed autoencoders in a dynamic nonlinear factor model structure with idiosyncratic components. The equivalence between maximum likelihood estimation and minimisation of mean squared error, together with the Universal Approximation Theorem allow to conceptualise the D$^2$FMs as computationally efficient approximations of the maximum likelihood estimates of generic nonlinear factor models.
The application of a simple version of our D$^2$FM with linear dynamic equations and a linear decoder, both in a Monte Carlo experiment, and in a big data real-time forecasting exercise shows the potential of the methodology.
The model capability can be further generalised and employed in multiple directions. For example, one could consider empirical applications that consider a nonlinear decoding structure, or nonlinear dynamic equations for the factors. Also, the loss function could be changed to allow for quantile estimates koenker1978regression,chen2018quantile. Finally, it is worth observing that the modularity and flexibility of the D$^2$FM allow to easily integrate alternative data (e.g. text data, satellite images, micro-data, etc) into the model. We leave these promising avenues of research to the future.