Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,423 characters · 14 sections · 21 citation commands
A Bayesian Gaussian Process Dynamic Factor Model
\doublespacing
\begingroup \footnotetext[1]{The views expressed in this paper are solely those of the authors and may differ from the views of the Bermuda Monetary Authority. No responsibility for them should be attributed to the Bermuda Monetary Authority.} \endgroup \singlespacing{\footnotesizeCorresponding author: Haroon Mumtaz ([email removed]). We thank Annika Camehl, Yizhou Kuang, Michele Lenza, and Aubrey Poon, as well as the participants of the Advances in Macroeconometrics Workshop in Manchester and the IAAE 2025 in Turin. Hauzenberger acknowledges funding by the Jubil\"aumsfonds of the OeNB, grant no. 18763.}
\thispagestyle{empty} \doublespacing
Nonlinearities are an important feature of macroeconomic and financial data. In just the last decade, the Global Financial Crisis, COVID-19 pandemic, and central banks reaching their effective lower bound provide examples. Capturing these nonlinearities has become an increasingly important issue---there is widening recognition that they are important for understanding and predicting macroeconomic dynamics.\footnote{For example, a broad strand of literature on macro-financial linkages \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[e.g.,][]{Brunnermeier2014, ABG2019} or the Phillips curve \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[e.g.,][]{hamilton2001parametric, BeaudryNBER24, Ball2022inflation, Benigno2023infla} documents instability in the underlying (structural) relationships. This is consequential for effective policy analysis and policy making. Modeling nonlinearities generally has a long tradition in macroeconometrics, with competing approaches ranging from regime-switching \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[e.g.,][]{hamilton1989new,terasvirta1994specification,sims2006were} to time-varying parameter \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[TVP, e.g.,][]{primiceri2005time,cogley2005drifts,koop2013large} models.}
In addition to the focus on modeling various forms of nonlinearities, a key feature of modern macroeconometric approaches is their scalability to high-dimensional datasets. These days it is straightforward to download hundreds of macroeconomic time series with a single click; examples for datasets include the US-based FRED-MD/QD \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[][]{FRED-MD,FRED-QD}, or other well-maintained multi-country databases. Exploiting macroeconomic “Big Data” (which typically involves many variables and few observations) to improve structural analysis and forecasting, and, consequently, policymaking, is a focus of researchers and practitioners alike; see \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{bok2018macroeconomic} for a recent review. And while there is no obvious best way to econometrically model these data \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[][]{giannone2021economic}, there is a popular front-runner---the dynamic factor model (DFM).
DFMs are a workhorse model of empirical macroeconomics \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see,e.g.,][]{aguilar2000bayesian, forni2000generalized, stock2002macroeconomic, kose2003international, GIANNONE2008, ChernisSekkelDFM, kaufmann2019bayesian}.\footnote{sw2016 and doz2020dynamic provide excellent surveys of the earlier literature.} A DFM assumes that there are a few fundamental forces in the economy which explain common dynamics of many time series. These fundamental forces (or latent factors) are usually modeled using a vector autoregression (VAR). A DFM thus compresses the data to work with a more parsimonious VAR, rather than shrinking a large VAR featuring many variables towards a simpler specification, another popular option \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{banbura2010large}.
Despite the overall popularity of (linear) DFMs, few papers model any nonlinear relationships between factors and observable variables. As we pointed out above, these may indeed be crucial to understand major macroeconomic fluctuations. Exceptions include DFMs with time-varying loadings \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{del2008dynamic, mumtaz2012evolving, korobilis2013assessing, zhou2014bayesian}, or DFMs with Markov-switching dynamics \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{Chauvet1998, chauvet2016dynamic, camacho2018markov}. More recently nonlinear DFMs with a squared/quadratic dynamics in the measurement or state equation \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{guerron2023financial} have been proposed. However, these studies impose specific functional forms in the context of inferring the latent factors.\footnote{An exception is Velasco:2024:Chapter3 who uses Bayesian Additive Regression Trees (BART) to model nonlinearities in a DFM. Unlike our approach, she employs a linear approximation to the relationship when estimating the factors. Another related approach is \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2025nonparametric}, who by contrast use a VAR augmented with static nonlinear factors modeled via BART.} In other words, they impose explicit restrictions on the link between the factors and observed variables.
Machine learning techniques have proven useful for modeling nonlinearities of unknown form in big macroeconomic datasets.\footnote{An incomplete list of examples includes \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{bassetti2014beta, kalli2018bayesian, farrell2021deep, medeiros2021forecasting, babii2022machine, goulet2022machine, jin2022infinite, huber2023nowcasting, clark2024investigating, chronopoulos2024forecasting, goulet2024macroeconomy, hauzenberger2024bayesian, hauzenberger2024nowcasting, hauzenberger2025gaussian}.} These techniques are appealing for several reasons. First, machine learning approaches are flexible and only require mild assumptions about the form of nonlinearities. Second, these methods are designed to avoid oversimplification, misspecification, and overfitting. Third, these sophisticated methods are well-suited for learning and identifying common patterns in large datasets, enabling efficient information extraction.
This paper introduces a general nonlinear DFM and develops a computationally feasible and fully Bayesian estimation algorithm. As mentioned above, imposing linearity in a typical DFM may be too restrictive, and we thus relax this assumption. Specifically, we propose to use Gaussian processes \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[GPs,][]{williams2006gaussian}, to obtain a nonparametric Gaussian Process Dynamic Factor Model (GP-DFM). GPs can capture a wide range of possible nonlinear relationships between latent (or observed) factors and high-dimensional data, and they have a successful track record in macroeconomic modeling \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{clark2024forecasting, hauzenberger2025gaussian}. Compared with other recent nonlinear DFM approaches, our approach is more flexible as this novel framework does not impose any specific type of nonlinearity. Instead, we place a prior (which is compatible with a rich menu of functions, and governed by a distance-based kernel function subject to only a few tuning parameters) directly on the functional relationship between common latent factors and the observed series.
Part of our contribution is bridging the machine learning literature on Gaussian Processes with macroeconometrics in a Big data context (where sample sizes are typically rather small). Specifically, we use GPs in a state space framework \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{turnerSSM,Frigola2013}. The main novelty of our framework lies in the nonlinear and flexible treatment of the measurement equation, while the state equation is assumed to follow a standard linear VAR. From both a practical and forecasting perspective, the use of a linear VAR in the state equation offers several appealing features.
First, interpretation of the latent common factors is analogous to a standard linear DFM. So, conventional tools from structural and reduced form VAR analysis can be readily used to the linear VAR in the state equation---for example, to produce forecasts as well as impulse response functions (IRFs) for the latent factors quickly. Once the full paths of these forecasts and the IRFs of the factors are known, the nonlinear mapping between observables and latent factors in the measurement equation can be exploited to generate forecasts or IRFs for the observed time series. This treatment results in the ability to calculate structural and reduced form quantities that can be time-varying or state dependent.
Second, this modeling strategy can also be viewed as a form of nonlinear dimension reduction, where high-dimensional macroeconomic data are assumed to lie on a lower-dimensional manifold or space. Any nonlinearity in the model arises exclusively from the relationship between the latent factors and the observables, which captures the notion of providing a more accurate/precise description of the variation in a large panel of time series. There is a relationship to other dimension reduction techniques such as local linear manifold regression \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{cheng2013local}, nonlinear principal components \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[PCs,][]{BAI2008}, deep/nonlinear dynamic factor models and autoencoders, which use neural networks to uncover complex patterns in high-dimensional data \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{andreini2023deepdynamicfactormodels, hauzenberger2023real, guerron2023financial, klieber2024non, snellman2024nonlinear, luo2025time}. However, unlike existing deep dynamic factor model approaches and popular two-stage procedures, our GP-DFM offers a fully consistent modeling framework. At its core, it represents a Bayesian state space model in which the functional relationship between low-dimensional latent factors and high-dimensional observables is explicitly modeled, with parameters and latent processes with precisely defined priors and posteriors. These posteriors are jointly estimated using a Markov Chain Monte Carlo (MCMC) algorithm (which is part of the contribution of this paper). This framework ensures both model consistency and proper Bayesian uncertainty quantification for all parameters and latent processes in the model, the latter being important for obtaining accurate predictive densities \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{geweke2010comparing}.
A fully nonlinear measurement equation is challenging to estimate, as it requires computationally intensive particle filtering techniques. To facilitate Bayesian estimation and provide a computationally feasible option for inference, we use the algorithm proposed in svensson2016computationally. This framework involves a linear approximation of the GP using basis functions derived from the spectral densities of the GP kernels \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{solin2020hilbert}. In addition, we rely on Particle Gibbs with Ancestor Sampling (PGAS) inspired by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{lindsten14a}, which greatly improves the efficiency of our algorithm. We further improve computational scalability by proposing an alternative parameterization of the GP kernel that uses an additive structure instead of a multiplicative one. This approach sacrifices some flexibility because it restricts interactions between the factors in the measurement equations, but there are significant computational advantages to this specification. Finally, since our approach uses a fully Bayesian MCMC algorithm, we can straightforwardly implement key model features such as stochastic volatility and shrinkage priors for the VAR process governing the state dynamics.
We provide two different empirical applications. In our first application, we forecast key macroeconomic variables from the FRED-QD dataset. The GP-DFM outperforms linear specifications of the DFM---a workhorse model at many policy institutions. Some of the superior performance is due to the GP-DFM forecasting real activity variables well during the COVID-19 period and the Federal funds rate when it is close to the effective lower bound. Interestingly, we find that with the nonlinear GP-DFM, a smaller number of factors is sufficient to extract the high-dimensional information than when imposing linearity. This supports the notion that in large macroeconomic datasets nonlinear factors can capture richer dynamics than their linear counterparts, which require more components to achieve comparable accuracy. In our second application, we measure international inflation dynamics with GP-DFM, and study nonlinearities and asymmetries arising from shocks to the global component of inflation in a large cross-section of economies.
The rest of the paper is structured as follows. Section (ref) lays out the econometric framework. This includes the DFM, GP priors subject to computationally favorable kernel specifications, and the posterior and predictive sampling algorithm. Section (ref) applies the GP-DFM in a forecast horserace with US data, while in Section (ref) we use the model to capture nonlinearities in international inflation dynamics. Section (ref) concludes.
A general nonlinear DFM relates $N$ observed macroeconomic time series (measurements), $\bm{y}_{t} = (y_{1t},\hdots,y_{Nt})'$, to a set of $D$ common latent factors $\bm f_t = (f_{1t}, \dots, f_{Dt})'$:
Equation ((ref)) is a measurement equation, where an unknown function $\bm G (\bullet):\mathbb{R}^D\rightarrow \mathbb{R}^N$ maps latent inputs $\bm f_t$ to observed outputs $\bm y_t$. In the following, we assume $\bm G(\bm{f}_t) = (g_1(\bm{f}_t), \hdots, g_N(\bm{f}_t))'$ collects equation-specific functions, with $g_i(\bullet): \mathbb{R}^D \rightarrow \mathbb{R}$, for $i = 1, \dots, N$. The idiosyncratic component $\bm v_t$ is Gaussian distributed with zero mean and variance-covariance matrix $\bm R = \text{diag}(r_{1},\dots,r_{N})$. While $\bm{G}(\bm{f}_t)$ captures variation common to all time series in $\bm y_t$, $\bm v_t$ is specific to each time series. Moreover, flexible equation-specific functions allow the common factors $\bm f_t$ to affect each variable differently, adequately reflecting the idea that different time series may feature different types of nonlinearities, captured by possibly different functional forms of $g_i$.
The nonlinear measurement equation is combined with a linear state equation for the factors, which follow a joint VAR($P$) process:\footnote{Note that more general multivariate processes, e.g., $\bm{f}_t = \bm{H}(\bm{f}_{t-1},\hdots,\bm{f}_{t-p}) + \bm{\varepsilon}_{t}$, can and have been used as state equations in nonlinear state space models \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[e.g.,][]{Frigola2013,guerron2023financial}; see \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{marcellino2024bookchapter} for a recent review of the macroeconometric literature using such models with (observed) macro-data. For the reasons we outlined in the Introduction (i.e., striking a balance between flexbility, computational efficiency, general ease of use and interpretability), we impose that $\bm{H}(\bullet)$ is linear a priori and leave nonlinear extensions for future research.}
with coefficient matrices $\bm A = (\bm A_1, \dots, \bm A_P)$ that relate $\bm f_t$ to their $P$ lags. In addition, $\bm{\varepsilon}_t = (\varepsilon_{1t}, \dots, \varepsilon_{Dt})'$ are zero mean state innovations with variance-covariance matrix $\bm{Q}$.\footnote{In the empirical application, we also investigate whether allowing for stochastic volatility (SV), with a time-varying $\bm Q_t$ rather than a constant $\bm Q$, improves forecast accuracy \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{C2011JBES, CR2015JAE}.}
In terms of features of this general framework, a few considerations are noteworthy. Eqs. ((ref)) and ((ref)) define a state space model. Specifically, $\bm f_t$ can be viewed as a lower dimensional representation of the high-dimensional data $\bm y_t$, which strikes a balance between explaining most of the variation in the data and parsimony. In this regard, a critical assumption is $D \ll N$.
A standard DFM assumes a linear mapping, $\bm{G}(\bm{f}_t) = \bm{\Lambda} \bm{f}_t$, with $\bm{\Lambda}$ being an $N \times D$ loadings matrix. In such a case, Eqs. ((ref)) and ((ref)) resemble a Gaussian linear state space model and for estimation standard algorithms such as forward sampling backward smoothing \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[FFBS,][]{carter1994gibbs, fruhwirth1994data} or the precision sampler \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{chan2009efficient} can be used. However, this linearity assumption may be too restrictive and is an assumption we wish to relax. To do so, one might consider incorporating nonlinear features, e.g., of the form $\phi(\bm f_t) = (\bm f_t', (\bm f_t \odot \bm f_t)')'$; a related approach is discussed in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{guerron2023financial}.
Related to this so-called “feature space” \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][chapter 2]{williams2006gaussian} two points are worth highlighting. First, $\phi(\bm f_t) = (\bm f_t', (\bm f_t \odot \bm f_t)')'$, or other simple transformations, would denote a deterministic projection of $\bm f_t$ onto $\phi(\bm f_t)$. However, we use a general form for $\phi(\bm f_t)$ which, by virtue of being more flexible, can approximate deterministic transforms but does not require knowing the appropriate transform ahead of time. Second, recall $\bm f_t$ is latent. Hence, any $\phi(\bm f_t)$ that includes nonlinear transforms of $\bm f_t$ leads to a nonlinear state space model, for which standard sampling algorithms can no longer be used. In short, our approach does not require pre-specifying the form of nonlinearity and we provide an estimation algorithm which is feasible for nonlinear measurement equations. The remaining sub-sections of the econometric framework will focus on these two considerations and carefully outline the contribution of this paper.
The key innovation of our model is the measurement equation features unknown nonlinearities. We choose GPs since they have a proven track record in macroeconomic (and many other) applications, and because of several properties that will result in a computationally feasible estimation algorithm. The key feature of GPs is that they model the covariance between observations through a simple kernel function. Kernel functions usually only have few parameters, but they are very expressive. Despite their sparse parametrization, GPs can indeed approximate a wide variety of functions.\footnote{There are theoretical equivalencies between neural nets and GPs \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{Neal1996, lee2018deepneuralnetworksgaussian}.}
We begin describing the model, starting with the $i^\text{th}$ measurement equation:
We assume a GP prior on each observable-specific function:
where $\mathcal{K}_{i}(\bullet)$ is a suitable covariance function (called the kernel) which we define below. Eq. ((ref)) denotes a prior over all possible functions $g_{i}$ that fit the $i^\text{th}$ equation well, without explicitly taking a stance of the functional form of $g_i$ and thus the feature space of $\bm f_t$. Through an appropriate choice of kernel we can allow for an infinite set of basis functions and feature space $\phi(\bm f_t)$. This means a GP can approximate any arbitrary continuous function, see also \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{williams2006gaussian}.
This paper focuses on two distinct versions of the commonly used squared exponential kernel. One variant is multiplicative across covariates and the other additive:
with hyperparameters $\bm{\theta}_{i} = (\xi_i, \{\ell_{ij}\}_{j = 1}^D)'$ defining covariances between two periods $t$ and $\tau$. Here, $\xi_i$ denotes an unconditional variance parameter and $\{\ell_{ij}\}_{j = 1}^D$ refer to factor-specific length scales. We label the former variant in Eq. ((ref)) the {multiplicative} (multi) GP-DFM, and the latter in Eq. ((ref)) the {additive} (hybrid) GP-DFM. Both variants define a stationary covariance function appropriate for macroeconomic data. However, they differ in the interactions between inputs (factors). The additive kernel sums the squared exponential terms instead of multiplying them. This difference has implications for computation and modeling flexibility, which we discuss in more detail below.
In the following, we first discuss the implications of GPs with a generic stationary kernel in the measurement Eq. ((ref)). If $\bm f_t$ were known, this would denote an otherwise standard GP regression, with inference being straightforward. However, $\bm f_t$ is latent which substantially complicates estimation and inference. To see this, define $\bm{g}_i = (g_{i}(\bm{f}_1),\hdots,g_{i}(\bm{f}_T))'$, then the GP takes the form:
with $\mathcal{K}_i(\bm{f}_t,\bm{f}_\tau)$ being the $(t,\tau)^\text{th}$ entry in $\mathcal{K}_i(\bm{F},\bm{F}')$. Considering the weight-space view of our GP model in the $i^{\text{th}}$ equation, we obtain:
Here, $\bm W_i$ denotes the lower Cholesky factor of $\mathcal{K}_i(\bm{F},\bm{F}') = \bm W_i \bm W_i'$, with $w_{i, t \tau}$ referring to the $(t, \tau)^\text{th}$ element of $\bm W_i$. Note that $\mathcal{K}_i(\bm{F},\bm{F}')$ denotes a full (symmetric) $T \times T$ variance-covariance matrix, implying that the measurement equation relates $y_{it}$ to the full history of latent states $\bm f_t$:
As alluded to by the weight space view, the measurement equation is rather involved, rendering a fully Bayesian approach feasible but computationally cumbersome \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{Frigola2013}.
In our paper, we approximate the GP using a set of nonlinear basis functions. These basis functions are based on the spectral density of the kernels in Eq. ((ref)), see \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{svensson2016computationally, solin2020hilbert, riutort2023practical}. Loosely speaking, this facilitates estimation by breaking the dependence on the full history of latent factors $\{\bm f_\tau\}_{\tau = 1}^{t}$, relying instead only on transformations of $\bm f_t$. In general, this approximation is enabled by the fact that any stationary covariance function can be represented in terms of their spectral densities $\mathcal{S}_i(\bm{\omega})$. Using eigenfunction expansion, any stationary covariance function with inputs $\bm{f}_t$ can then be expressed as:
with $\mathcal{S}_i$ denoting the spectral density, $\bm{\lambda}_m$ are vectors of eigenvalues and $\phi_m(\bm{f}_t)$ are eigenfunctions of the Laplacian operator that is used to establish the approximation in domain $\Omega$. More formally, $\Omega \in [-L_1,L_1] \times \hdots \times [-L_D,L_D]$ (i.e., $\Omega \subset \mathbb{R}^D$) describes the support on which the GP approximations are valid \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see also][for details]{solin2020hilbert}. Crucially, these eigenvalues and eigenfunctions do not depend on the specific choice of the kernel or its hyperparameters---only the spectral density $\mathcal{S}_i$ does.
Following \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{riutort2023practical}, the kernels can then be approximated using a set of $M$ basis functions:\footnote{The covariance function within $\Omega$ with inputs $\bm{f}_t,\bm{f}_\tau \in \Omega$ may be written as $\mathcal{K}_{i}(\bm{f}_t,\bm{f}_\tau) = \sum_{m=1}^{\infty} \mathcal{S}_{i}(\sqrt{\bm{\lambda}_m}) \phi_m(\bm{f}_t)\phi_m(\bm{f}_\tau)$. We may truncate this sum so that it can be used to obtain the weight space representation of the GP in Eq. ((ref)).}
These approximations, for $i = 1, \dots, N$, allow us to rewrite Eq. ((ref)) as:
with the $N\times M$-matrix $\bm{C}$ comprising the weights $c_{im}$ across equations $i = 1,\hdots,N,$ and basis functions $m = 1,\hdots,M,$ and the basis functions are stored in an $M \times 1$-vector $\bm{\Phi}(\bm{f}_t) = (\phi_1(\bm{f}_t),\hdots,\phi_{M}(\bm{f}_t))'$. Note that the derivations above are applicable to any generic stationary kernel. Next, we focus on distinct aspects that arise for the multiplicative and additive variants that we mentioned earlier. Indeed, the two squared exponential covariance functions in Eq. ((ref)) are both stationary and therefore the associated spectral densities are Fourier duals:
with a vector $\bm{\omega} = (\omega_1,\hdots,\omega_D)'$ in the frequency domain.
{\sffamilyMultiplicative squared exponential kernel}. Let $\mathbb{S}\in\mathbb{R}^{M \times D}$ be the matrix of all possible combinations of univariate eigenfunctions across dimensions (i.e., all $D$-tuples), and we obtain for $m = 1,\hdots,M$:
With a $D$-dimensional input space, the total number of eigenvalues and eigenfunctions used for the approximation is equal to the number of possible combinations of univariate eigenfunctions across all dimensions, i.e., $M = \tilde{M}^D$ where $\tilde{M}$ is the number of basis functions for each dimension. In case of the multiplicative squared exponential kernel, the approximation changes the computational complexity to $\mathcal{O}((T + 1) M^D)$ from $\mathcal{O}(T^3)$, a significant saving as long as the number of factors remains low.
{\sffamilyAdditive squared exponential kernel}. Large macroeconomic datasets often require a moderate number of factors to fit the data. For example, FRED-QD find seven or eight factors in the full FRED-QD dataset, which would be rather computationally demanding when used with the multiplicative squared exponential kernel. For this reason, we introduce introduce the additive squared exponential kernel as a more computationally efficient alternative.
Start by noting that the GP with this kernel is $g_{i}(\bm{f}_t) \sim \mathcal{GP}(0,\sum_{j = 1}^D \mathcal{K}_{i}(f_{jt},f_{jt}))$, which can be split into the sum of $D$ univariate GPs, $y_{it} = g_{i}(f_{1t}) + \dots + g_{i}(f_{Dt}) + v_{it}$. Let $\tilde{M}$ now denote the number of basis functions for each dimension $(j = 1, \dots, D)$. Then the total number of eigenvalues/eigenfunctions used for the approximation is equal to $M = \tilde{M}^D$ for the multiplicative kernel, while it is $M = D \times \tilde{M}$ for the additive kernel. This is because we can consider $\mathbb{S}\in\mathbb{R}^{\tilde{M}}$ for each dimension separately and simply sum over the factor-specific basis functions given the additive structure:
This seems to be a classic trade-off between model flexibility and computational feasibility. The main limitations of the additive kernel variant is that it ignores interactions between latent factors, potentially underfitting highly complex relationships. On the other hand, it comes with the benefit of substantially improved scalability in the number of latent factors ($D$). One of the empirical contributions of our paper is investigating this trade-off in macroeconomic data.
As we adopt a fully Bayesian approach, we need to assign prior distributions to all unknown parameters in the model. For the measurement equation, we closely follow recent developments in the Gaussian process literature \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{vandvaart2008bayesian, bhattacharya2014anisotropic, lindsten14a, hauzenberger2025gaussian}. Specifically, for the equation-specific (homoskedastic) measurement errors in Eq. ((ref)), we assume independent inverse Gamma priors, $r_{i} \sim \mathcal{G}^{-1}\left(\nu_r, S_r\right)$, with $\nu_r = 3$ and $S_r = 0.3$ for all $i = 1, \dots, N$ equations. For the hyperparameters of the Gaussian processes kernels, we use Gamma priors on the unconditional variances, $\xi_i \sim \mathcal{G}\left(\nu_\xi, S_\xi \right)$, and Gamma priors on the inverse length-scale parameters, $\ell^{-2}_{ij} \sim \mathcal{G}\left(\nu_\ell, S_\ell\right)$. In the empirical application, We set $\nu_\xi = \nu_\ell = 0.5$ and $S_\xi = S_\ell = 0.5$, which is a weakly informative choice.
To specify priors for the parameters in the state equation, Eq. ((ref)), we follow the recent macroeconomic DFM/VAR literature and adopt global--local shrinkage priors where appropriate \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{huber2019adaptive, kaufmann2019bayesian, fruhwirth2025sparse}. Specifically, we use horseshoe \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[HS,][]{carvalho2010horseshoe} priors, as they have been shown to work well in diverse contexts virtually free of tuning. Any hyperparameters associated with these priors are updated in a data-driven manner, making the approach suitable for VARs of any size and providing a highly adaptive prior specification in the state equation. We defer technical details on the prior setup for the state equation to Section (ref) of the Appendix.
To sample from the joint posterior distribution of our nonlinear state space model we rely on a Particle Gibbs approach. This section sketches the main steps of the algorithm, and additional details are provided in Section (ref) of the Appendix. Our Markov Chain Monte Carlo (MCMC) involves three main steps/blocks:
Typically none of the individual parameters of multivariate time series models are of primary interest, but functions of them are. In terms of computing forecasts and IRFs with the GP-DFM model, it suffices to note that the procedure is relatively straightforward. For the latent factors, we can iteratively project them $\tau$ steps ahead, $\bm f_{t+1}, \dots, \bm f_{t+\tau}$, using conventional tools for the linear VAR in the state equation ((ref)). Once the path of the latent factors is known, we can use the approximated GP measurement equation ((ref)) to obtain forecasts/IRFs for the endogenous variables.
In a first application, we assess the overall performance of the GP-DFM in an extensive forecasting exercise using data for the US macroeconomy. We use the popular FRED-QD dataset \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{FRED-QD} which is a large information set spanning a broad array of quarterly macroeconomic data. Our evaluation sample is long enough to include several US macroeconomic events of importance: the dot-com bubble, the great financial crisis (GFC), and the COVID pandemic. We benchmark our main model specifications against a standard linear DFM which is a workhorse model in policy institutions for forecasting \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{GIANNONE2008, ChernisSekkelDFM}. Section (ref) provides details on the forecasting setup as well as the competing specifications and Section (ref) summarizes the main forecasting results.
The FRED-QD dataset covers 248 quarterly time series of which we select $N = 105$ to replicate the information content of stock2012disentangling. This rich macroeconomic dataset covers real activity, price, and financial variables. All variables are transformed to stationarity as suggested in FRED-QD. For a list of variables used and transformations applied, see Appendix (ref). We define real output (GDPC1), non-farm employment (PAYEMS), consumer prices (CPIAUCSL) and the Federal funds rate (FEDFUNDS) as target variables for the forecasting exercise. We use an expanding window estimation scheme. The estimation sample starts in 1965Q1 and is recursively updated each quarter over the hold-out sample which spans 1992Q1 to 2023Q4. We consider forecast horizons of one-quarter, one-year and two-years ahead. For our target variables, we assess both joint and marginal density forecast accuracy. We use the energy score (ES) to evaluate the joint predictive densities and the continuous rank probability score (CRPS) to evaluate the marginal predictive densities \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep{GR2007JASA}.
In the forecasting exercise, we consider 16 different model specifications, which vary across key modeling choices. For all models we use $P = 4$ lags in the state equation. The margins we explore are linearity versus nonlinearity (or kernel choice) of the DFM: linear DFM for the linear specification, GPDFM-A for the GP-DFM with an additive kernel, and GPDFM-M for the GP-DFM with a multiplicative kernel. Moreover, we vary the number of factors ($D = \{2, 4, 8\}$) and consider either homoskedastic variances or SV of the state innovations. The out-of-sample forecasting exercise is used to assess all of these choices because it is robust and we do not have clear theoretical guidance for them. Finally, we need to decide on a few GP-DFM specific choices. For the multiplicative kernel, we set $\tilde{M} = 8$ for $D = 2$ and $\tilde{M} = 4$ for $D = 4$, and we do not consider $D = 8$.\footnote{For the multiplicative kernel, $D = 8$ would result in a total of $M = \tilde{M}^8$ basis functions. Thus, setting $\tilde{M} = 4$ yields $M = 65536$, which is computationally burdensome.} For the additive kernel, we set $\tilde{M} = 8$ for all values of $D$. The boundary condition $L$ is specified in a semi-automatic manner as $1.2$ times the maximum absolute value of the first $D$ PCs.\footnote{We have empirically tested different $\tilde{M}$ values and our choices seem to be sensible ones in terms of the accuracy of the GP approximation.}
{\sffamilySummary of findings}. Before zooming into the details we provide a broad overview of the forecasting results. There are four main take aways. First, overall, the nonlinear DFM specifications consistently perform well and outperform the linear DFM, providing some improvements over the linear benchmark when considering both joint and marginal forecast losses. Second, the computationally efficient additive kernel method performs similarly to the more flexible multiplicative kernel, indicating that there is little trade-off between computation speed and forecast accuracy. Third, using only a relatively small number of factors is sufficient to forecast our target variables accurately. Forth, as is common in the macroeconomic literature, allowing for SV improves forecast accuracy.
In the following, Figs. (ref) to (ref) in this section report forecast evaluation metrics for both the joint forecast performance and the marginal forecast performances. For ease of readability, all figures follow a common convention. The benchmark model is a linear DFM with $D = 4$ and SV, and its absolute forecast score is reported directly. For other model specifications, we show the ratio relative to the benchmark. Ratios below one indicate superior performance relative to the benchmark (indicated by green-shaded cells), while ratios above one indicate inferior performance (indicated by red-shaded cells). Statistical significance is assessed using a one-sided DM1995JBES test, with asterisks denoting the significance of gains relative to the benchmark at the $1\%$ ($^{***}$), $5\%$ ($^{**}$), and $10\%$ ($^{*}$) level. For each horizon, the best performing specification is highlighted in bold.
{\sffamilyJoint forecast performance over different evaluation periods}. Figure (ref) shows the joint forecast performance across horizons for three different evaluation samples. In panel (a) we consider the full evaluation sample, in panel (b) we stop prior to the COVID period in $2019$Q$4$, and in panel (c) we stop prior to the GFC in $2007$Q$4$. This strategy should help identify forecast gains over time. Recall our target variables are real output (GDPC1), non-farm payroll employment (PAYEMS), consumer prices (CPIAUCSL) and the Federal funds rate (FEDFUNDS); the losses in this context are based on the joint predictive distribution of these variables.
We first focus on panel (a) in Fig. (ref), which shows the joint forecast performance for the full evaluation sample and serves as the main measure of model performance. As suggested in the key summarized findings mentioned above, we find that a GP-DFM variant with two factors and SV in the state equation consistently outperforms all other specifications across horizons. This model variant exhibits significant forecast improvements around $6$ to $14\%$ lower losses (as indicated by the dark green shaded cells).
When we vary the evaluation sample, we observe significant differences in the relative forecast performance of the GP-DFMs. Comparing panels (a) and (b) of Fig. (ref) shows that including the COVID period systematically lowers ES ratios for our GP variants relative to the benchmark (but also relative to other linear competitors). This suggests that the impressive performance is partly driven by the period of the pandemic, and indicates that a flexible, nonlinear GP-DFM improves forecast accuracy over a more restrictive linear DFM during this exceptional period. Moreover, the improvements of the GP-DFM variants are largest at shorter horizons, outperforming linear competitors by larger margins. The best performing model improves upon the benchmark by nearly $14\%$ for the full evaluation sample (and around $9\%$ for the pre-COVID sample) at the one-quarter ahead horizon, while the improvement for the two-year ahead horizon is about $6\%$, irrespective of whether the evaluation stops before COVID.
Comparing panels (b) and (c) of Fig. (ref), it is evident that the gains arise at least partially in the context of the GFC and the zero lower bound (ZLB). Indeed, the shorter (pre-GFC) period can be considered a relatively calmer economic period, and while there are some improvements for GP-DFM variants for the one-quarter and one-year horizon, the linear DFM benchmark is only occasionally outperformed, and typically not significantly so. For the two-year ahead horizon, however, the nonlinear model is more accurate, with significantly lower losses.
{\sffamilyMarginal forecast performance}. Figures (ref) and (ref) show the marginal forecast performance for the four target variables for different subsamples and across several forecast horizons. For one-quarter ahead forecasts we specifically look at the performance pre-GFC and the full sample. We do so because particularly at this horizon, GP-DFM shows sizable differences in joint relative performance, and the marginal perspective allows to investigate the sources of this accuracy premium also in terms of the cross-section of variables. For higher-order forecast horizons, we report the marginal scores only for the full evaluation sample in the interest of saving space.
Focusing first on the results at the one-quarter ahead horizon in Fig. (ref), the CRPS ratios echo joint forecast performance patterns, and show that accuracy varies substantially across evaluation samples. Models with SV and a low number of factors perform particularly strongly for the full sample. All GP-DFM variants significantly outperform the linear benchmark at a $1\%$ significance level (as well as the other linear competitors) for the Federal funds rate (FEDFUNDS). Here, the best performing specification, a GP-DFM with two factors, a multiplicative kernel, and SV, outperforms the benchmark by almost $14\%$. This is likely responsible for a sizable portion of the joint forecast performance. But this is not the sole source of the larger forecast accuracy gains for the full evaluation sample, as the GP-DFM also forecasts real output and inflation well. For real output (GDPC1), the best performing model is a GP-DFM equipped with two factors, an additive kernel, and homoskedastic errors, improving about $7\%$ upon the benchmark (but these improvements are not statistically significant). For inflation (CPIAUCSL), most GP-DFM variants significantly outperform the benchmark at the $5\%$ significance level, showing forecast performance similar to the best performing specification (in this case, a linear DFM with eight factors and SV). By contrast, for payroll employment (PAYEMS), there are no improvements.
Turning to the pre-GFC evaluation sample in panel (b), a linear DFM equipped with SV appears highly competitive and is difficult to outperform using GPs during this period (which largely coincides with the comparatively calm Great Moderation). Indeed, for the pre-GFC evaluation sample, virtually all GP-DFM variants perform similarly as the linear benchmark for GDPC1, CPIACUSL, and FEDFUNDS. However, the proposed nonlinear models struggle with forecasting PAYEMS at short horizons, even more so in “normal” economic times.
Moving to the one-year and two-year ahead horizons in Fig. (ref), we see that a GP-DFM variant is consistently the best performing model across all variables. Specifically, the GP-DFM with SV and two factors performs very well. While we observe strong performance for the Federal funds rate similar to the short horizon, forecast performance for real activity variables generally improves at the longer horizons. In addition, inflation forecasts are better by a small but statistically significant amount (at the $5\%$ significance level). The comparatively poor performance for short horizon forecasts of payroll employment also vanishes, with GP-DFM exhibiting more accurate forecasts than the benchmark (although these are not statistically significantly better).
{\sffamilyProperties of real output growth forecast densities}. In the final part of this section, we showcase the features of the predictive densities of a GP-DFM variant using real output growth around the COVID pandemic as a means to detect outliers and nonlinearities. We can highlight what drives differences in forecast performance by comparing the predictive densities of the benchmark model against a GP-DFM variant.
Figure (ref) shows the difference between the predictive densities of a GP-DFM variant with four factors, an additive kernel, and SV versus the linear DFM with four factors and SV.\footnote{Such a comparison between predictive densities has been proposed in DSZ2022JE.} To ensure a fair comparison, we focus on the same number of factors and allow both models to feature SV. By inspecting this figure, we can see that the GP-DFM predictive densities place more mass around the realized outcomes, while the linear DFM is more dispersed. This cannot be attributed to SV, which features in both models. Instead, the GP-DFM's tighter forecast densities are likely a key factor behind the improved forecast performance. Specifically, this is evident during the calmer post-COVID period (from 2021 to 2023) following the sharp downturn and rebound in output growth at the onset of the pandemic.
Our second application demonstrates how to use the model for a semi-structural analysis, by investigating the importance of nonlinearities in macroeconomic data. Specifically, we explore the co-movement in cross-country consumer price index (CPI) inflation rates in the spirit of mumtaz2012evolving and Kose2025. We estimate the following model:
The panel of (quarterly) CPI inflation rates, $\pi_{it}$, is taken from the World Bank inflation database compiled by Kose2025, which runs from $1970$Q$1$ to $2023$Q$4$ and includes $N_{\text{\texttt{D}}} = 27$ developed countries as well as $N_{\text{\texttt{E}}} = 34$ emerging economies (EMDEs), such that $i=1,2,\hdots,N_{\texttt{D}},N_{\texttt{D}}+1,\hdots, N,$ with $N = N_{\text{\texttt{D}}} + N_{\text{\texttt{E}}}$. Let $\mathbb{I}(\bullet)$ refer to an indicator function, and define $s_i = \mathbb{I}(i \leq N_{\text{\texttt{D}}})$. That is, $s_i = 1$ if the respective country is among the developed countries and $0$ when it is an emerging economy. Following Kose2025, inflation is decomposed into a global “world” factor (indicated with $\text{\texttt{W}}$) that is common to all countries ($f_{\text{\texttt{W}}, t}$) and regional factors that are specific to country groups, namely developed countries ($f_{\text{\texttt{D}}, t}$) and EMDEs ($f_{\text{\texttt{E}}, t}$). As each inflation series might show a certain of idiosyncratic persistence, we assume that the idiosyncratic errors in the observation equation follow an AR($2$) process: $e_{it}= \sum_{q=1}^{2} \rho_{iq}e_{it-q}+ v_{it}$, $v_{it} \sim N(0,r_{i})$, which is a straightforward extension.
Unlike the previous literature, we allow the relationship between factors and observables to be nonlinear via the unknown functions $g_i(\bullet)$ and $h_i(\bullet)$. Our measurement equation allows for state-dependence and time-varying contributions from factors to observables. Note Eq. ((ref)) corresponds to using an additive kernel similar to the previous section, with $g_{i}(\bullet)$ specific to $f_{\text{\texttt{W}}, t}$ and $h_{i}(\bullet)$ specific to the relevant regional factor. We assume a GP with a squared exponential kernel for each function that are approximated using $\tilde{M} = 8$ basis functions. The regional factors are disentangled by imposing exclusion restrictions as in Eq. ((ref)). That is, the regional factor for developed economies does not affect EMDEs and vice versa. We normalize the sign and scale of the factor by using a tight prior on the first row of loadings in each regional group.\footnote{The prior is obtained using an initial estimate of the weights obtained conditional on principal component estimate of the factors.} The latent factors $\bm f_{t}= (f_{\text{\texttt{W}}, t}, f_{\text{\texttt{D}}, t}, f_{\text{\texttt{E}}, t})'$ are assumed to evolve according to a VAR(4) process, see Eq. ((ref)), with $P = 4$.
Figure (ref) plots the estimated posterior distribution of the three factors. The global factor has a distinct peaks in the earlier part of the sample in 1974Q1 correspond to the aftermath of the oil shock. After 1990, the global factor displays a declining trend consistent with the Great Moderation. This came to an end with a sharp rise in $2022$Q$1$, possibly associated with supply disruptions in the post-COVID period.
The factor specific to developing economies is persistently high during the mid-1970s and then during the early 1980s. The decline in this factor in 1985 precedes the fall in global factor and is much sharper. In the subsequent period, inflation in developing countries appears to have been stable with the abrupt disruption after 2021. The variation in inflation specific to EMDEs appears radically different. The main peaks in the EMDE factor occurs in 1989 and 1990 capturing the aftermath of debt crises and hyperinflation in several countries.
It is interesting to the note that while the estimated factors from the benchmark model are broadly similar to that obtained from the linear DFM, there are noticeable differences. For example, the global factor from the benchmark model has a smaller peak in 1980 when compared to the estimate from the linear DFM and displays sharper movements in 1994 and 2008. The developed country factor in the benchmark case displays less variation than the linear factor during the 1990s and 2000s while the EMDE factor has a smaller peak during the mid-1980s than the linear counterpart.
This section investigates if there are meaningful nonlinearities between factors and observables. We study this by varying the size and sign of shock applied to the factors and compute generalized forecast error variance decompositions (GFEVD). A linear model would show no difference in the GFEVD for different size/sign shocks whereas if nonlinearities are important there would be substantial differences. For example, a large negative global inflation shock would have a proportionally different impact than a small positive shock in a nonlinear model.
More specifically, we compute generalized impulse responses \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[GIRFs,][]{koop1996impulse}, using a Cholesky decomposition of the reduced form covariance matrix to orthogonalize the structural disturbances \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see also][for related discussions]{pfarrhofer2025scenario}. In the state equation VAR, we order the global factor first, followed by the developed economy and the EMDE factor, respectively. The GIRFs are calculated using the estimated factors in every fourth quarter of the sample as initial conditions. We use the generalised forecast error variance decomposition (GFEVD) to estimate the contribution of the global and group specific factors to country-specific inflation rates. The GFEVD is calculated using the method proposed by lanne2016generalized whereby the GIRFs are used in the standard formula for the FEVD of a linear VAR model. This ensures that the decomposition adds up to one, while retaining properties such as dependence on initial conditions and on the size and sign of the shock being considered.
Panel (a) in Fig. (ref) shows that the contribution of the global and regional factors to the forecast error variance of inflation in developed countries. There is a clear dependence of the contributions on the sign and size of the shock. For example, the contribution of the world (regional) factor is lower (higher) for large positive shocks in Austria, Belgium, Cyprus, Denmark, Italy, Netherlands and Japan. The opposite result appears to hold for countries such as the UK, Greece, New Zealand, and Portugal.
Panel (b) Figure (ref) shows that for many EMDEs the importance of the world factor also depends on the size and sign of the shock. This is especially apparent for countries such as Burkina Faso, Colombia, Morocco and Ghana. For these countries, the regional factor contributes more to inflation fluctuations when shocks to this factor are negative.
In this paper we propose a class of nonlinear dynamic factor models using Gaussian Processes. Our first contribution is developing an estimation approach. Furthermore, we show that this model performs well and is easy to use for structural applications. The estimation algorithm we develop makes estimation of the GP-DFM feasible by addressing many of the drawbacks of Gaussian process models. For example, a GP models have difficulties scaling and filtering unobserved states in nonlinear models is difficult. We show how to solve these problems with an efficient estimation approach using Hilbert space approximations and particle sampling with ancestor sampling. Additionally, we show how to further increase computational speed using an additive kernel so that the model can handle moderate numbers of factors. Succinctly, we make using a nonlinear dynamic factor model feasible in macroeconomic applications.
We assess the GP-DFM in a forecasting exercise using US macroeconomic data. We show that the GP-DFM with SV outperforms linear versions of the model. Specifically, it performs well when there are nonlinearities in the data such as at the effective lower bound and during COVID. On average, its performance is driven by its more precise predictive densities. Finally, there is little cost for using an additive kernel, which excludes interactions between factors, suggesting that interaction terms are not an important type of nonlinearity in this dataset.
We also demonstrate how such a model could be used for a semi-structural analysis. We show that our specification with a nonlinear measurement equation and linear state equation makes for easy inference. The linear state equation allows us to apply structural VAR techniques to the factors while the nonlinear measurement equation allows for flexibility in how the factor maps to observables (state-dependence and time-variation). Our example investigates the co-movement in cross-country CPU inflation rates. First we find differences between the factors estimated in a linear DFM and a GP-DFM. Second, we investigate the importance of nonlinearities by examining how the size and sign of shocks affect contributions to the forecast error variance decomposition. This shows the importance of allowing for nonlinearities.
{\setstretch{0.9} \addcontentsline{toc}{section}{References} }\doublespacing