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.
103,254 characters · 14 sections · 90 citation commands
Bayesian dynamic variable selection in high dimensions
\onehalfspacing \setcounter{page}{1}
Regression models that incorporate stochastic variation in parameters have been used by economists at least since the work CooleyPrescott1976. Thirty years later, Granger2008 argued that time-varying parameter models might become the norm in econometric inference since, as he illustrated via White's theorem, time variation is able to approximate generic forms of nonlinearity in parameters. Indeed, initiated by the unprecedented shocks observed during and after the Global Recession of 2007-9, a large recent literature has established the importance of modeling time variation in the intercept, slopes and variance of regressions for forecasting economic time series; see StockWatson2007 for a representative example of a model using only a stochastic intercept and volatilities. At the same time, the stylized fact that economic predictors are short-lived -- that is, relevant for the dependent variable only in short periods\footnote{An alternative terminology for such periods, which is due to Farmeretal2018, is “pockets of predictability”.} -- has emerged in various forecasting problems such as inflation KoopKorobilis2012, stock returns DanglHalling2012 and exchange rates Byrneetal2018. Following these observations, there is no shortage of recent econometric work on methods for penalized estimation of time-varying parameter models via classical or Bayesian shrinkage, as well as variable selection methods; see for example Belmonteetal2014, Bitto2019, KalliGriffin2014, CallotKristensen2014, Korobilis2019, Kowaletal2019, NakajimaWest2013, RockovaMcAlinn2017, UribeLopes2017 and YousufNg2019.
In this paper we add to this literature by proposing a new dynamic variable selection prior and a novel, for the field of economics, Bayesian estimation methodology. In particular, we propose to use variational Bayes (VB) inference to estimate time-varying parameter regressions using state-space methods. Variational inference has long been used in data science problems such as large-scale document analysis, computational neuroscience, and computer vision Bleietal2017. Nevertheless, it is only relatively recently that posterior consistency and other theoretical properties of these methods have been explored by mainstream statisticians WangBlei2019. Variational inference is a unified estimation methodology which shares similarities with the Gibbs sampler that many economists traditionally use to estimate time-varying parameter models StockWatson2007. Like the Gibbs sampler, parameter updates are derived for one parameter at a time conditional on all other parameters using an iterative scheme. Unlike the Gibbs sampler, there is no repeated sampling involved and the output of VB is typically the first two moments of the posterior distribution of parameters. Our first task is to introduce this estimation scheme in the context of TVP regressions, and contrast it to existing estimation algorithms used in economics for capturing structural change.
Our second contribution lies on the development of a dynamic variable selection prior that is a conceptually straightforward extension of the static variable selection prior of GeorgeMcCulloch1993. The dynamic extension of this prior allows to tackle the non-trivial econometric problem of allowing some predictor variables to enter the TVP regression, model only in some periods of the full estimation sample. With $p$ predictors and $T$ time periods, dynamic variable selection involves choosing the “best” among $2^{p}$ models at each point in time $t$, for $t=1,...,T$. Such procedure is in line with strong, recent empirical evidence that different factors might be driving predictability of economic variables over time; see Rossi2013 for a thorough review of this idea. By specifying our new prior within a variational Bayes framework, we are able to derive an algorithm that is numerically stable and can be extended to much larger $p$ and $T$ than was possible before.\footnote{In particular, many of the algorithms cited above, such as KoopKorobilis2012, KalliGriffin2014, or NakajimaWest2013, are unable to scale up to regressions with hundreds of predictors.}
We show, via a Monte Carlo exercise and an empirical application, that our proposed algorithm works well in high-dimensional, sparse, time-varying parameter settings. Using artificial data we establish that the new algorithm is precise in estimation and in dynamic variable selection, even in settings with more predictors than time-series observations. In a forecasting exercise of various measures of price inflation, we illustrate that our methodology applied to a time-varying parameter regression with 400+ predictors is able to beat a wide range of linear and nonlinear forecasting regressions. The empirical results provide strong evidence that the new algorithm can achieve estimation accuracy comparable to Markov chain Monte Carlo algorithms, while being much faster to run. The additional feature of dynamic variable selection successfully prevents overparametrization, since our high-dimensional TVP specification is able to beat both parsimonious time series models with no predictors as well as factor models and penalized likelihood estimators.
The remainder of the paper proceeds as follows. Section 2 introduces the basic principles of VB inference for approximating intractable posteriors, and applies these principles to the problem of estimating a simplified time-varying parameter regression model. Section 3 introduces the the novel modelling assumptions, namely dynamic variable selection and stochastic volatility, and derives an estimation algorithm within the VB framework. Section 4 assesses the new algorithm on simulated data. In Section 5 we apply the new methodology to the problem of forecasting US inflation using time-varying parameter regressions with many predictors.
as variational Bayes (VB) is not an established estimation methodology in econometrics, we first provide a generic discussion of VB methods in approximating intractable posterior distributions. We then apply the generic concepts and formulas to the specific problem of estimating a simplified time-varying parameter regression model with known measurement error variance.\footnote{Readers already familiar with these concepts can skim through this section, and focus on our novel methodology that is described in the following section} Detailed reviews of variational Bayes can be found in Bleietal2017 and OrmerodWand2010, among several others. Variational Bayes estimation of state-space models is described in detail in the monograph of smidl2006variational, as well as research papers such as BealGhahramani2003, Tranetal2017, and Wangetal2016.
Consider data $y$, latent variables $s$ and (latent) parameters $\theta$. Our interest lies in time-varying parameter models that admit a state-space form. Hence, $s$ represents unobserved state variables, such as time-varying regression coefficients and time-varying measurement error variances, and $\theta$ represents all other parameters, such as the error covariances in the state equation. The joint posterior of interest is $p\left(s, \theta |y\right)$ with associated marginal likelihood $p\left( y\right) $ and joint density of data and parameters $p\left(y,s,\theta \right)$. When the joint posterior is complex and computationally intractable, we can define an approximating density $q\left( s,\theta \vert y \right)$ that belongs to a family $\mathcal{F}$ of simpler distributions defined over the parameter space spanned by $s,\theta$. The main idea behind variational Bayes inference is to make this approximating posterior distribution $q\left( s,\theta \vert y \right)$ as close as possible to $p\left( s,\theta |y\right) $, where “distance” is measured with the Kullback-Leibler divergence\footnote{For notational simplicity we henceforth abbreviate multiple integrals using a single integration symbol.}
That is, the aim is to find the optimal $q^{\star}\left( s,\theta \vert y \right)$ that solves
Insight for why $KL\left(q \vert \vert p \right)$ is a desirable distance metric arises from a simple re-arrangement involving the log of the marginal likelihood OrmerodWand2010 where it can be shown that
Because $KL\left(q \vert \vert p \right)$ is non-negative (it is exactly zero when $q\left( s,\theta \vert y \right) = p\left( s,\theta |y\right) $), the quantity {
}becomes a lower bound for the marginal likelihood $ p(y)$.\footnote{In the following we denote as $\mathbb{E}_{q(\bullet)}$ the expectation w.r.t to a function $q(\bullet)$.} The function $ \mathcal{G} \left(q(s,\theta \vert y )\right)$ is known as the Evidence Lower Bound (ELBO). Therefore, instead of minimizing the objective function $KL\left(q \vert \vert p \right)$ (which cannot be evaluated) we can find an approximating density $q^{\star}\left( s,\theta \vert y \right)$ that maximizes the marginal data density $p\left( y \right)$ by maximizing the ELBO. We emphasize that $\mathcal{G}$ is a functional on the distribution $q(s,\theta \vert y)$. As a result, the ELBO can be maximized iteratively using calculus of variations.
If we assume for simplicity the so-called (in Physics) mean field factorization of the form $q\left( s, \theta \vert y \right)=q\left( \theta \vert y \right) q\left( s \vert y \right)$, it can be shown\footnote{A formal and thorough derivation of these ideas is given in the excellent monograph of smidl2006variational; see Theorem 3.1 and subsequent results.} that the optimal choices for $q\left( s \vert y\right)$ and $q\left( \theta \vert y \right)$ are
The first expression denotes the expectation over $q(\theta \vert y)$ of the conditional posterior for $s$, and the second expression denotes the expectation over $q(s \vert y)$ of the conditional posterior for $\theta$. Because $q(\theta \vert y)$ is a function of $q(s \vert y)$, and vice-versa, the above quantities can be approximated iteratively instead of relying on more computationally expensive numerical optimization techniques. Given an initial guess regarding the values of $(\theta, s)$, VB algorithms iterate over these two quantities until $\mathcal{G} \left(q(s,\theta \vert y)\right)$ has reached a maximum. Due to similarities with the Expectation-Maximization (EM) algorithm of Dempsteretal1977, this iterative procedure in its general form is sometimes referred to as the Variational Bayesian EM (VB-EM) algorithm; see BealGhahramani2003. It is also worth noting the relationship with Gibbs sampling. Like Gibbs sampling, equations (ref) and (ref) involve the full conditional posterior distributions. But unlike Gibbs sampling, the VB-EM algorithm does not repeatedly simulate from them and is computationally much faster.
Before collecting all building blocks of our proposed methodology, we outline a VB algorithm for the univariate TVP regression with known measurement error variance $\underline{\sigma}^{2}$. This simplified model is of the form
where $y_{t}$ is the time $t$ scalar value of the dependent variable, $t=1,..,T$, $\mathbf{x}_{t}$ is a $1 \times p$ vector of exogenous predictors and lagged dependent variables, $\varepsilon _{t} \sim N\left( 0,\underline{\sigma}^{2} \right) $, $\bm \eta_{t} \sim N\left( 0,\bm W_{t} \right)$ with $ \boldsymbol{W}_{t}=diag\left(w_{1,t},...,w_{p,t} \right)$ a $p \times p$ diagonal matrix\footnote{By restricting $\bm W_t$ not to be a full covariance matrix, coefficients $\beta_{it}$ and $\beta_{jt}$ are uncorrelated a-posteriori for $i \neq j$, which might not seem like an empirically relevant assumption. However, allowing for cross-correlation in the state vector $\bm \beta_{t}$ can result in counterproductive increases in estimation uncertainty, with this problem being significantly more pronounced in higher dimensions. A diagonal $\bm W_{t}$ allows for a more parsimonious econometric specification, less cumbersome derivations of posterior distributions, and faster and numerically stable computation; see also Belmonteetal2014, Bitto2019 and RockovaMcAlinn2017 who adopt a similar assumption.} and $\bm w_{t} = [w_{1,t},...,w_{p,t}]^{\prime}$ a $p \times 1$ vector. In likelihood-based analysis of state-space models it simplifies inference if it is assumed that $\varepsilon _{t}$ and $\bm \eta_{t}$ are independent of one another and we do adopt this assumption here. Finally, we use a notational convention where $j,t$ subscripts denote the $j^{th}$ element of a time varying state variable, or parameter, observed only at time $t$, while $1:t$ subscripts denote all the observations of a state variable from period $1$ up to period $t$.
The model in equations (ref) and (ref) has unknown parameters $\left( \bm \beta_{1:T},\bm w_{1:T} \right)$. Following the analysis of the previous subsection we first consider the independent prior on the initial conditions $\bm\beta_{0},\bm w_{0}$ of the form
where $w_{j,0}$ is the $j^{th}$ element of $w_{0}$, and $Gamma(a,b)$ denotes the Gamma distribution with shape parameter $a$ and rate parameter $b$, that is, the definition of the Gamma distribution that has mean $a/b$ and variance $a/b^{2}$. The time $t$ prior, conditional on observing information up to time $t-1$, is given by the Chapman-Kolmogorov equation
where $\mathcal{B}$ is the support of $\bm \beta_{t}$ and $\mathcal{W}$ the support of all $w_{j,t}$. Finally, once the measurement $y_{t}$ is observed, we obtain from Bayes theorem the following time $t$ posterior distribution
This Bayesian joint posterior distribution is rarely analytically tractable, even if conjugate prior densities have been specified. However, posterior conditionals can be tractable, and this is why in macroeconomics TVP models are predominantly estimated using the Gibbs sampler; see StockWatson2007 for an example. Nevertheless, sampling repeatedly using (Markov chain) Monte Carlo methods is computationally prohibitive in high-dimensional settings or in settings with more flexible likelihood and prior distributions. For that reason we define a tractable variational density as an approximation to the exact intractable time $t$ posterior, that is, we define $p\left(\bm \beta_{t},\bm w_{t} \vert \bm y_{1:t} \right) \approx q\left( \bm \beta_{t}, \bm w_{t} \vert \bm y_{1:t} \right)$. Among all possible functions $q\left( \bm \beta_{t}, \bm w_{t} \vert \bm y_{1:t} \right)$ we want to obtain the one that has hyperparameters that minimize the relative entropy with the true posterior. Following the discussion earlier in this section, this problem is equivalent to maximizing the evidence lower bound (ELBO) of the log-marginal likelihood, that is, it is the solution to {
}
This maximization problem is simplified once we assume the mean field factorization of the form $q \left( \bm \beta_{t} , \bm w_{t} \vert \bm y_{1:t} \right) = q\left( \bm \beta_{t} \vert \bm y_{1:t} \right) \prod q\left(w_{j,t} \vert \bm y_{1:t} \right)$ so we can optimize $\bm \beta_{t}$ and $\bm w_{t}$ sequentially. As a result, using variational calculus smidl2006variational we can show that the ELBO is maximized by iterating through the following recursions
Both formulas above become equalities after the addition of a normalizing constant. The first expression is an expectation with respect to the probability density $\prod_{j} q\left(w_{j,t} \vert \bm y_{1:t-1} \right)$, that is, we can write equation (ref) using the following form
where $q\left(\bm \beta_{t} \vert y_{1:t-1}\right)$ is the time $t$ prior of $\bm \beta_{t}$ obtained from the time $t-1$ posterior $q\left( \bm \beta_{t-1} \vert \bm y_{1:t-1}\right)$ using the Kalman filter recursions, and the term $p \left( \bm w_{t} \vert \bm y_{1:t-1} \right)$ in (ref) disappears because the expectation is w.r.t the variational posterior of $w_{t}$. This latter representation of $q\left( \bm \beta_{t} \vert \bm y_{1:t} \right)$ can be trivially updated by a Normal distribution, with moments given by the Kalman filter and smoother; see smidl2006variational for detailed derivations. We can use similar arguments in order to show that equation (ref) is an expectation that leads to a $q\left( \bm w_{t}^{-1} \vert \bm y_{1:t} \right)$ of the form $G \left(c_{j,t},d_{j,t}\right)$ (or equivalently to $q\left( \bm w_{t} \vert \bm y_{1:t} \right)$ that is inverse Gamma).
Algorithm (ref) provides pseudocode for the basic VB estimation problem described in this section, without assuming either a (dynamic) variable selection prior or stochastic volatility in the measurement equation. In the following section we drop these two unrealistic assumptions.
We rewrite for convenience the univariate time-varying parameter model
where we define now $\varepsilon_{t} \sim N(0,\sigma_{t}^{2})$ with $\sigma_{t}^{2}$ a stochastic (time-varying) variance parameter, and we assume that the dimension $p$ of $\bm \beta_{t} = \left(\beta_{1,t},...,\beta_{p,t} \right)^{\prime}$ is large and possibly $p \gg T$.
The core ingredient of our modeling approach is a dynamic variable/model selection strategy. We specify a dynamic variable selection (DVS) prior that extends the “static” variable selection prior of GeorgeMcCulloch1993 that was originally developed for the constant parameter regression using MCMC and is of the form
for $j=1,...,p$, where $\underline{c}$, $g_0$ and $h_{0}$ are fixed prior hyperparameters. Variable selection principles require us to set $\underline{c} \rightarrow 0$, such that the first component in the prior for $\beta_{j,t}$ shrinks the posterior towards zero, while the second component has variance $\tau_{j,t}^{2}$ which is “large enough” in order to allow for unrestricted estimation. The choice between the two components in the prior for $\beta_{j,t}$ is governed by the random variable $\gamma_{j,t}$ which is distributed Bernoulli and takes values either zero or one. If $\gamma_{j,t}=1$ the prior for $\beta_{j,t}$ has a Normal prior with zero mean and variance $\tau_{j,t}^{2}$, while if $\gamma_{j,t}=0$ the prior variance becomes $\underline{c}\tau_{j,t}^{2}$.
Early papers such as GeorgeMcCulloch1993 give very broad guidelines on choosing values for $\underline{c}$ and $\tau_{j,t}^{2}$ such that the first component in equation (ref) has small enough variance (to force shrinkage) and the second component has large enough variance (to allow unrestricted estimation). More recently, NarisettyHe2014 show that selecting and fixing the prior variances of such mixture priors could, as $T$ and $p$ grow, lead to model selection inconsistency. The authors suggest to specify these parameters to be certain deterministic functions of the data dimensions $T$ and $p$. In our case, we do fix $\underline{c}=10^{-4}$ such that the first component has always smaller variance, but we assume $\left(\tau_{j,t}^{2}\right)^{-1}$ is a random variable that has a Gamma prior. That way this parameter is always updated by the information in the data likelihood. The choice of a Gamma prior for $\left(\tau_{j,t}^{2}\right)^{-1}$ implies that the marginal prior for $\beta_{j,t}$ is a mixture of leptokurtic Student's T distributions whose components could tend to shrink $\beta_{j,t}$ towards zero, regardless of whether $\gamma_{j,t}$ is zero or one. Therefore, the proposed prior is able to find patterns of dynamic sparsity as well as impose dynamic shrinkage in time-varying parameters, a property that is very desirable in high-dimensional settings.\footnote{In signal processing a signal (regression coefficient vector) is typically sparse by default, that is, the researcher knows a-priori to expect that estimates of several coefficients will tend to be exactly zero. In economics, the sparsity assumption might not be empirically founded in certain settings; see the discussion in Giannoneetal2017. In such cases, a dense model may be preferred, that is, a model where all predictors are relevant with varying weights. While factor models and principal components have been used widely to model dense models in macroeconomics, shrinkage methods are also quite reliable for this task. In particular, we note the result in DeMoletal2008 that forecasts from Bayesian shrinkage are highly correlated to forecasts from principal components.}
Finally, it becomes apparent that under this variable selection prior setting, $\widehat{\pi}_{0,t} = \mathbb{E} \left( p\left(\pi_{0,t}\right) \right)=\frac{1}{2}$ is the time $t$ prior mean probability of inclusion of all predictors in the TVP regression, while the quantity $\widetilde{\pi}_{j,t} = \mathbb{E} \left( p\left( \gamma_{j,t} \vert \bm y_{1:T} \right)\right)$ is the posterior mean probability of inclusion in the regression of predictor $j$ at time period $t$, simply referred to as the posterior inclusion probability (PIP). Due to the fact that all of the hyperparameters $\bm \gamma,\pi$ and $\bm \tau^{2}$ are time-varying, our prior allows to obtain time-varying PIPs whose interpretation extends this of PIPs in constant parameter settings, such as the one in GeorgeMcCulloch1993, in a straightforward way.
In terms of tackling estimation using this prior we note that adding the prior (ref) to our benchmark TVP specification introduces some peculiarity: by combining equations (ref) and (ref) we end up having two conditional prior structures for $\beta_{j,t}$, namely
where we define $ v_{j,t} = \left( 1 - \gamma_{j,t} \right)^{2}\underline{c} \times \tau_{j,t}^{2} + \gamma_{j,t}^{2} \tau_{j,t}^{2} $ and $V_{t}$ is the $p \times p$ diagonal matrix comprising the elements $v_{j,t}$. Following ideas in Wangetal2016 we combine the two priors for $\beta_{t}$ described above by rewriting the state equation as\footnote{The derivation is straightforward using arguments in the previous subsection, see equation (ref). Define $q \left( \beta_{t} \vert y_{1:t-1} \right)$ to be the time $t$ variational Bayes prior of $\beta_{t}$ given information at time $t-1$. Then we have
where the simplification occurs due to the fact that $\bm \beta_{t-1}$ is known and fixed (i.e. not a random variable) given information at time $t-1$. Therefore, the formula above specifies the new, joint time $t$ prior of $\bm \beta_{t}$ given the two priors in equations (ref)-(ref).}
where $\widetilde{\bm \eta}_{t} \sim N \left( \bm 0, \widetilde{\bm W}_{t} \right)$, with parameter matrices $\widetilde{\bm W}_{t} = \left( \mathbb{E} \left(\bm W_{t}\right)^{-1} + \mathbb{E} \left(\bm V_{t}\right)^{-1} \right)^{-1}$ and $\widetilde{\bm F}_{t} = \widetilde{\bm W}_{t} \times \mathbb{E} \left(\bm W_{t}\right)^{-1}$, where $\bm W_{t} = diag \left(w_{1,t},...,w_{p,t} \right)$ and $\bm V_{t} = diag \left(v_{1,t},...,v_{p,t} \right)$, and all expectation operators are with respect to $q\left( \bm \beta_{t} \vert \bm y_{1:t}\right)$. Under this formulation we can observe that the joint prior variance for $\beta_{j,t}$ is a function of both $w_{j,t}$ and $v_{j,t}$, $\forall j=1,...,p$. Therefore, the TVP regression model with dynamic variable selection prior can be written using a new state-space form, with measurement equation given by (ref) and state equation given by (ref).
Application of algorithm (ref) to the transformed state-space model consisting of equations (ref) and (ref) provides as output estimates $\bm m_{t\vert T}$ $\forall t$, that is, the smoothed posterior mean of $q \left( \bm \beta_{t} \vert y_{1:T} \right)$. Conditional on these estimates, derivation of the update steps for $\gamma_{j,t}$, $\tau^{2}_{j,t}$ and $\pi_{0,t}$ relies also on deriving the expectations of these variables with respect to $q \left( \bm \beta_{t} \vert y_{1:T} \right)$. Therefore, extending the analysis of the previous section to accommodate these new parameters, and similar to derivations found in Gibbs sampling approaches to variable selection GeorgeMcCulloch1993, the updating steps for the parameters in the dynamic variable selection prior are the following
for each $t=1,..,T$ and $j=1,...,p$, where again expectations $\mathbb{E}$ are with respect to the VB posteriors of each of the parameters showing up on the right-hand side of the equations above.
A known regression variance is far from a realistic assumption for most datasets. When forecasting macroeconomic data, so is the assumption of an unknown variance that is constant over time. A vast recent literature highlights the importance of time-varying volatility in improving point and density forecasts ClarkRavazzolo2015, and the purpose of this subsection is to accommodate estimation of the parameter $var(\varepsilon_{t}) = \sigma_{t}^2$ in the VB setting. Several elegant algorithms for VB inference in stochastic volatility models exist in the literature. For example, Naessethetal2017 introduce a variational Bayes Sequential Monte Carlo (SMC) algorithm for stochastic volatility models. Tranetal2017 propose a variational Bayes method for intractable likelihoods that does not rely on the mean field approximation, and apply their algorithm to the estimation of a stochastic volatility model.
Nevertheless, such algorithms assume an explicit time-series model for the stochastic volatility parameter, an assumption that is only useful in a setting where one is interested in forecasting volatility. In a macroeoconomic setting we are interested in forecasting $y_{t}$ and not its volatility (as it would be the case in empirical asset pricing). At the same time, previous empirical work shows that there are no statistically important differences when forecasting with alternative specifications of macroeconomic volatility.\footnote{For example, ClarkRavazzolo2015 compare a range of specifications for time-varying variance parameters in univariate and multivariate autoregressive models, and any differences among such specifications are not statistically important (while all volatility specifications are always better relative to constant variance specifications).} For that reason, our aim here is not only to render estimation of stochastic volatility precise, but at the same time numerically reliable and computationally efficient. In order to achieve this, we build on variance discounting ideas for dynamic linear methods as described in WestHarrison1997; see also RockovaMcAlinn2017.
Define $\phi_{t} = \frac{1}{\sigma_{t}^{2}}$ to be the precision (inverse variance). Following WestHarrison1997 we assume that the time $t-1$ posterior of $\phi$ has the following conjugate form
We do not specify an explicit time series model for the dynamics of $\phi$ (e.g. stochastic volatility or GARCH) because the posterior for $\phi_{t}$ wouldn't be conjugate to the likelihood and we would fail to obtain fast updates. In order to maintain this conjugacy we specify instead the time $t$ prior of the form
for a variance discounting factor $0<\delta<1$, subject to a choice of hyperparameters $a_{0}$ and $b_{0}$. By doing so, we assume that $\phi_{t}$ is centered around $\phi_{t-1}$ as if this parameter had random walk dynamics,\footnote{Even though we haven't specified an explicit time series evolution for $\phi_{t}$, by using results in Uhlig1994 we can show that the proposed variance discounting methodology is equivalent to assuming the following specification:
for a parameter $\gamma_{t} \vert y_{1:t-1} \sim Beta \left( \delta a_{t-1}/2, (1-\delta) a_{t-1} /2 \right)$.} since it holds that $ \mathbb{E} \left( \phi_{t} \vert y_{1:t-1} \right)= \mathbb{E} \left( \phi_{t-1} \vert y_{1:t-1} \right)$. However, based on the properties of the Gamma distribution, the dispersion of $\phi_{t}$ is larger to that of $\phi_{t-1}$.
Under this scheme the variational Bayes update of $\phi_{t}$, that is, its time $t$ posterior mean has the form
where $a_{t} = 1/2 + \delta a_{t-1}$ and $d_{t} = \frac{1}{2} \left[ \left(y_{t} - \bm x_{t} \bm m_{t \vert T}) \right)^{2} + \bm x_{t} \bm P_{t \vert T} \bm x_{t}^{\prime} \right] + \delta b_{t-1}$, where $\bm m_{t \vert T},\bm P_{t \vert T}$ are the smoothed mean and variance of $\beta_{t}$. Using this scheme, past information in the data is discounted exponentially by the factor $\delta$. The scalar $\delta$ can be seen as a prior hyperparameter whose choice determines how much relative weight we give to recent versus older observations, that is, it determines how fast we expect the precision parameter to change over time. For $\delta=1$ we obtain the posterior under a standard recursive update scheme (similar to recursive OLS), while typical values that would allow for faster time-variation in the precision/variance would be between 0.8 and 0.99. Values lower than 0.8 are not empirically advised, since they allow for a large amount of time-variation and stochastic variance estimates become very noisy. In the empirical exercise we set $\delta=0.8$, a choice that reflects our prior expectation that macroeconomic data have many abrupt breaks in their second moments and excess kurtosis during recessions (implying variances that can move very fast over time).
The previous formulas pertain to the iterative updating of $\phi_{t}$ given $\phi_{t-1}$. Estimates of $\phi_{t}$ can be smoothed using subsequent observations $t+1,...,T$. Following WestHarrison1997 we can run a backward recursive filter of the form
for $t=T-1,...,1$, where $\widetilde{\phi}_{t} = \mathbb{E}_{q(\beta_{t} \vert y_{1:T})}\left(\phi_{t} \vert y_{t+1} \right)$ and $\widetilde{\phi}_{T} = \widehat{\phi}_{T}$. Once we obtain this update for the precision $\phi_{t}$, a posterior mean estimate of the volatility $\sigma_{t}^{2}$ can be obtained simply as the inverse of $\widetilde{\phi}_{t}$.
Here we provide details of the exact parameter updates that result from VB inference in our proposed specification. Algorithm (ref) outlines our proposed Variational Bayes Dynamic Variable Selection (henceforth, VBDVS) algorithm. This Algorithm shows an accurate picture of how this would look like when programmed using a language like MATLAB or R: while there are many parameters involved in our specification, the code is short and it involves simple scalar operations (meaning it is very fast). The only cumbersome operation is the inversion of the $p \times p$ matrix $\bm P_{t+1|t}$ in line 14 which has worst case complexity $\mathcal{O}\left(p^3\right)$ for each $t$. There are four main blocks in this algorithm. Lines 4-12 are a result of straightforward application of the Kalman filter on the state-space model of equations (ref) and (ref), and lines 13-17 show the backwards (smoothing) recursions. Lines 18-27 update the prior hyperparameters of the DVS prior for $\bm \beta_{t}$. Finally, lines 28-33 provide updates for the stochastic volatility parameter, as discussed in the previous subsection.
In this section we evaluate the performance of the new estimator using artificial data. Although we view the algorithm as primarily a forecasting algorithm, it is also important to investigate its estimation accuracy in an environment where we know the true data generating process (DGP). Thus, we wish to to establish that the VBDVS is able to track time-varying parameters satisfactorily and establish that the dynamic variable selection prior is able to perform shrinkage and selection with high accuracy (at least in cases where we know that the DGP is that of a sparse TVP regression model). We also wish to investigate the computational gains that arise from application of variational Bayes methods on the complex dynamic variable selection prior structure.
In all our experiments we use the following DGP:
Our benchmark specification sets $\underline{\bm \theta} = \left(-1.7, 2.9, 1.4, -2.3, \bm 0 \right)$, $\underline{\sigma}^{2} = 0.1$, $\underline{\rho} = \underline{\phi} = 0.99$, $\underline{\delta} = \underline{\xi}= T^{-1/2}$. In the specification above $\bm s_{j}$ is $T \times 1$ vector of either zeros or ones, such that $\beta_{j,t}=\theta_{j,t}$ when $s_{j,t}=1$, and zero otherwise. We set $s_{1,t}=1$ for $t=1,...,\lfloor T/3 \rceil-1$ and zero otherwise, $s_{2,t}=1$ $\forall t=1,...,T$, $s_{3,t}=1$ for $t=1,...,\lfloor T/2 \rceil-1$ and zero otherwise, $s_{4,t} = 0$ for $t=1,...,\lfloor T/2 \rceil-1$ and zero otherwise. These choices mean that $\beta_{1,t}$ is zero during the last third of the sample, $\beta_{2,t}$ is a relevant predictor in all periods, $\beta_{3,t}$ is zero during the last half of the sample, and $\beta_{4,t}$ is zero during the first half of the sample. Any other coefficient for $j=5,...,p$ is zero at all periods, i.e. $s_{j,t}=0$ $\forall$ $j>4$, $t=1,...,T$. By doing so, we simulate a situation where only one predictor is relevant in all time periods, three predictors are relevant only in certain subsamples of the data, and all remaining $p-4$ predictors are irrelevant for $y$ at all time periods.
After we generate artificial data, we compare three competing estimation algorithms for TVP models: i) our variational Bayes dynamic variable selection (VBDVS) algorithm, ii) the EM algorithm implementation of the dynamic spike and slab (DSS) of RockovaMcAlinn2017, and iii) Gibbs sampling (MCMC) estimation of the TVP model using the fast algorithm of ChanJeliazkov2009. While there are numerous other algorithms available for estimating TVP models, our limited choice of algorithms reflects our desire to simulate exclusively high-dimensional models. By doing so, we exclude most of the recently proposed Bayesian methodologies cited in the Introduction. These methodologies introduce various flexible parametrizations (like we do) that result, however, in the need for many tuning parameters and estimation via MCMC, such that they become unreasonably cumbersome for $p>50$. Our model instead, as we demonstrate in detail later, requires very straightforward tuning. The default prior setting we use for the VBDVS algorithm is based on the case Prior 3 presented in (ref) in the next section. The settings used in the DSS and MCM algorithms are discussed in the Online Supplement to this paper. In order to compare numerically these algorithms we generate $R=100$ datasets from the above DGP for various choices of sample size and total number of predictors, namely $T=100,200,500$ and $p=50,100,200$. Subsequently squared deviations between true and estimated parameters are calculated, and then averaged over the $T$ time periods, and $p$ predictors. To be precise, if we let $\left( \beta_{t}^{true} \right)$ denote the true artificially generated coefficients and $\left( \beta_{t}^{j}, \sigma_{t}^{j} \right)$, for $j=DVS,DSS,MCMC$, the estimates of these coefficients, we calculate the sum of mean squared deviations (MSD) statistic as
where $r=1,...,500$ denotes the number of Monte Carlo iterations.
(ref) shows the coefficient estimates from VBDVS for the case $T=p=200$. This plot compares the posterior median (green solid lines) versus the true generated coefficients (black dashed lines). The 16$^{th}$ and 84$^{th}$ percentiles over the 100 Monte Carlo iterations are also shown as a shaded area around the posterior median. Only the first 20 coefficients, out of the possible 200, are plotted. The first row shows the four coefficients that, at least in some periods, are non-zero, followed by 16 coefficients that are exactly zero. It is impossible to plot the remaining 180 coefficients in the DGP that are exactly zero, but their estimates are represented fairly well by the estimates of coefficients $\beta_{5,t}$ - $\beta_{20,t}$ shown in (ref). Under the assumption of sparsity in the DGP, the VBDVS algorithm is able to recover the true coefficients with accuracy. Not only the coefficients that are zero in the DGP in all periods are correctly estimated to be zero, but also the three coefficients that are zero only in certain subsamples are estimated precisely. When a coefficient is initially zero and later in the sample becomes important (see coefficient $\beta_{4,t}$), and vice-versa (see coefficients $\beta_{1,t}$ and $\beta_{3,t}$), the dynamic variable selection algorithm is able to identify and jump quickly to the new state. (ref) shows that the true reason why estimation is so precise -- even in such a demanding case with 200 time-varying coefficients for only 200 observations -- is because the estimates of the time-varying posterior inclusion probabilities (PIPs) of each predictor are recovered with precision in the first instance. By identifying correctly which variables should be excluded from the regression model in each period results in shrinking many coefficients to zero and allowing to preserve enough degrees of freedom for estimation of non-zero coefficients.
(ref) shows the values of the MSD statistics for the three algorithms under the different combinations of $T$ and $p$. Given that the MSD statistics measure deviation from the true coefficient, lower values imply that a certain estimation algorithm has done better recovery of the coefficients generated by the DGP. In all cases VBDVS has the best performance among all competing algorithms. The estimation error of the MCMC algorithm is quite large mainly because the algorithm is unable to shrink all $p-4$ coefficients in the DGP that are exactly zero. The DSS algorithm provides a better fit since it is also an algorithm that does dynamic variable selection and shrinkage. Its performance is slightly inferior to VBDVS, but the results should not be taken as final evidence. While we have done all effort to follow the settings suggested by RockovaMcAlinn2017, there might be other priors that could improve the performance of this algorithm.
Another important feature of the VBDVS algorithm is its fast computing time. While it is not surprising that our algorithm is faster compared to MCMC, our algorithm can provide substantial savings in high-dimensional settings compared to the DSS that relies on the EM algorithm. Columns 6-8 in (ref) reveals that VBDVS can be multiple times faster than both DSS and MCMC algorithms.
Following a large literature on time-varying parameter models in macroeconomics, our primary target is to forecast quarterly US inflation. While there exists mixed empirical evidence about the potential of very large datasets to improve forecasts of inflation, our aim is to demonstrate here that the new dynamic variable selection methodology can successfully extract, period-by-period, predictive information from a large number of predictors. For that reason we build a novel, high-dimensional dataset that brings together predictors from several mainstream aggregate macroeconomic and financial datasets.\footnote{While one could also think of potential predictors in disaggregated panels obtained in surveys, internet, or documents (text data), such novel sources are typically proprietary and would make our results hard to replicate.} Our building block is the FRED-QD dataset of McCrackenNg2020, which we augment with portfolio data used in Juradoetal2015, stock market predictors from WelchGoyal2008, survey data from University of Michigan consumer surveys, commodity prices from the World Bank's Pink Sheet database, and key macroeconomic indicators from the Federal Reserve Economic Data for four economies (Canada, Germany, Japan, UK). All data are quarterly, and span the period 1960Q1-2018Q4. All variables are adjusted from their respective sources for seasonality (where relevant), and we additionally remove extreme outliers.\footnote{Following StockWatson2016, we replace outliers using the median of the preceding five observations. An outlier is defined to be any observation that satisfies $\vert y_{t} - m \vert / iqr >\kappa$, where $m$ is the median of $y$, $iqr $ is the interquantile range, and $\kappa = 4.5$.}
The dataset has in total 444. Out of these we forecast the series (FRED-QD mnemonics in parentheses): GDP deflator (GDPCTPI), total CPI (CPIAUCSL), core CPI (CPILFESL), and PCE deflator (PCECTPI). When each of these price series, $P_{t}$, is used as the dependent variable to be forecasted $h$-quarters ahead we transform it according to the formula $y_{t+h} = (400/h) \ln \left( P_{t+h} - P_{t} \right)$. We forecast these transformed series one at a time, and the remaining three price series are included in the list of exogenous predictor variables (443 in total). The predictor variables are transformed using standard norms in the literature McCrackenNg2020: i) levels for variables that are already expressed in rates (e.g. unemployment, interest); ii) first differences of logarithm for variables measuring population (e.g. employment), variables expressed in dollars (e.g. GDP), commodity prices, and some indexes (e.g. Industrial production); and iii) second differences of logarithm for price and consumption indexes, as well as deflator series. The online supplement describes in detail all variables and transformations, and provides links to all sources.
Before we set up a comprehensive out-of-sample forecasting exercise, we first assess in-sample estimates from the VBDVS by doing small sensitivity analysis to various prior choices. This exercise is intended to demonstrate that the new algorithm provides reasonable estimates of trends, volatilities and other parameters. Most importantly it serves as a way to clarify that, despite the fact that our prior is heavily parametrized, prior elicitation in the VBDVS algorithm becomes a reasonably straightforward task. As it is impossible to present estimates of the TVP model using all variables in our dataset as predictors, we focus on a small TVP model where GDP deflator regressed on an intercept, two own lags, and the first five principal components from the 443 exogenous predictors (eight predictors in total).
Out of all parameters and hyperparameters defined in our algorithm it is only a handful that are crucial for inference and forecasting, while others can be fixed to reasonable or uninformative values and possibly have little effect on forecasting. (ref) lists all hyperparameters one need to choose in the VBDVS algorithm, and does an explicit separation into “Important” and “Fixed” hyperparameters. Starting from the latter, $a_{0}$ and $b_{0}$ are the initial scale and rate parameters of the initial condition of the precision parameter in equation (ref). Setting $a_{0}=b_{0}=0.01$ implies that the precision has prior mean one and variance 10, which is a reasonable uninformative choice for an inverse variance parameter. Next, we set $\delta=0.8$ for reasons explained in (ref). Given that $p$ is very large to allow us to obtain meaningful prior information about the regression coefficients $\bm \beta_{t}$ (e.g. using a training sample), we allow their initial condition $\bm \beta_{0}$ to be fairly uninformative by setting $\bm m_{0}= \bm 0$ and $\bm P_{0}=4 \bm I_{p}$. The parameter $\underline{c}$ in the dynamic variable selection prior has to be small (see discussion in (ref)) and how small it exactly is, affects the way the algorithm selects each of the two Normal components in the spike and slab prior -- that is, it affects the choice between a certain $\beta_{j,t}$ being restricted or not. We prefer to fix this parameter to $\underline{c}=0.0001$ and allow only $\tau_{j,t}^{2}$ and its prior to determine the ratio of the prior variances of the two Normal components in the mixture prior.
The parameters that are important in our high-dimensional setting are the ones affecting the two prior variances of the time-varying coefficients $\bm \beta_{t}$, namely the hyperparameters of $\tau_{j,t}^{2}$ and $w_{j,t}$. Our first prior choice, denoted as “Prior 1” in (ref), selects $c_{j,t}=100, d_{j,t}=1$ such that $w_{j,t}$ has a prior mean of 0.01 and prior variance 0.0001. This conservative choice restricts movements $\beta_{j,t}$ to be very persistent and excludes the case of frequent, noisy jumps. Such prior is used widely in empirical macroeconomic applications, see for example the “business as usual” prior motivated in CogleySargent2005 for the case of a vector autoregression with time-varying parameters. We subsequently set an uninformative prior on $\tau_{j,t}^{2}$ by setting $g_{0}=h_{0}=0.01$. The dashed lines in (ref) represent (posterior mean) coefficient estimates from our eight-predictor model: coefficient $\beta_{1,t}$ is the time-varying intercept (trend inflation), coefficients $\beta_{2,t},\beta_{3,t}$ correspond to the first two lagged values of GDP deflator, and coefficients $\beta_{4,t}$ to $\beta_{8,t}$ correspond to the five principal components. As a comparison, we plot posterior mean estimates from the same time-varying parameter regression estimated with MCMC (using identical settings as in the Monte Carlo comparison). The MCMC-based estimates can be broadly thought of as the unrestricted equivalents of the VBDVS algorithm, since they are not based on any form of dynamic variable selection or hierarchical shrinkage. The intercept and first lag coefficients are virtually identical using the two algorithms. However, all remaining coefficients are penalized heavily by the VBDVS algorithm. Variation over time of these coefficients is very moderate and restricted to be close to zero for many time periods.
In order to examine the effect that the prior has on the time evolution of the coefficients, we change the initial condition for $w_{j,t}$ to have hyperparameters $c_{j,0}=d_{j,0}$ and we leave the same uninformative prior for $\tau_{j,t}^{2}$. The posterior mean coefficient estimates in (ref) exhibit an interesting pattern. By allowing a looser prior on $w_{t}$ the parameters that are unrestricted (intercept and first lag), do exhibit larger amount of time-variation compared to the MCMC estimates. However, the remaining coefficients that were previously restricted to be close to zero, are now forced more aggressively towards zero in all time periods. This demonstrates the fact that our algorithm imposes the state-space model in equation (ref), where the variance of $\beta_{j,t}$ is a function of both $w_{j,t}$ and $v_{j,t}$ (where the latter, is in turn a linear function of $\tau_{j,t}^{2}$). Therefore, allowing for a looser $w_{j,t}$ tends to introduce more noise in the state-space model, and for that reason the dynamic variable selection prior compensates for this increased noise by shrinking more aggressively. While there is this compensation effect and coefficient estimates won't explode as quickly as the model without the dynamic variable selection prior (recall that $\beta_{j,t}$ evolves as a non-stationary random walk), it is not advisable to use such a lose prior on $w_{j,t}$.
For that reason, our final prior (called Prior 3 in (ref)) returns to the conservative choice $c_{j,0}=100$ and $d_{j,0}=1$, and sets instead $g_{0}=1$ and $h_{0}=12$. (ref) shows the estimates from this prior. Once again the VBDVS estimates of the intercept and first lag coefficients are identical to the estimates from the MCMC algorithm. The remaining coefficients are again heavily penalized but there are also many time periods where these evolve unrestrictedly. As a matter of fact, this prior allows the time-varying coefficients to exhibit distinct and abrupt jumps between periods where they are zero and periods where they are unrestricted. This pattern of time-variation is more in line with the findings of the previous literature that there are pockets of predictability or, put differently, that economic predictors are short-lived (see discussion in the Introduction).
In order to have a visual assessment of the time pattern of dynamic variable selection and shrinkage, panel (a) of (ref) plots the posterior inclusion probabilities of each regressor associated with the time-varying coefficient estimates presented in (ref). These seem to show the exact periods where each coefficient moves from a state of being restricted to zero to a state where it is not zero. Panel (b) of the same figure shows the posterior mean of the stochastic volatility estimate from VBDVS versus the estimate from MCMC. These two estimates are fairly similar, showing that the specification of time-varying variances in the VBDVS does a good job at capturing known peaks in GDP deflator inflation volatility. Any differences in volatility estimates reflect the fact that the two algorithms assume different specification of $\sigma_{t}^{2}$ and also use different priors in the estimation of $\bm \beta_{t}$.
For all these reason, we build all of our forecasting models in the next subsection based on this last prior.\footnote{Due to the fact that the choice $h_{0}=12$ looks in (ref) to penalize possibly excessively the small model with just eight coefficients, in the next subsection we adapt only this hyperparameter depending on the number of predictors we have available. Otherwise, all other hyperparameters are identical to the ones in the column labelled Prior 3 in (ref).}
We forecast inflation using models of the form
where $y_{t+h} $ is $h$-step ahead inflation (see (ref) for a definition) regressed on an intercept, two own lags and exogenous predictors. We use a variety of forecasting models. Some benchmark models are based on equation (ref) but assume constant coefficients (i.e. $\alpha_{t}=\alpha$, $\phi_{1,t} = \phi_{1}$ and so on), while others assume different sets of exogenous predictors. However, what all models have in common is that they always include an intercept and two own lags of inflation. Given that our dataset is much larger than datasets used before for forecasting inflation, in order to avoid confusion by specifying different combinations or subsets of predictors, we only distinguish four simple categories of models: i) models with no predictors (i.e. only intercept and autoregressive terms); ii) models with first five principal components as predictors; iii) models with sixty principal components as predictors; and iv) models with all 443 predictors. Our list of models representing each category is the following
The choice of models is based on their simplicity and replicability. In particular, the Gaussian Process Regression, Partial Least Squares, and Elastic Net algorithms are based on built-in functions in MATLAB's Statistics and Machine Learning Toolbox MatlabSTB, and are fairly easy to set up. Estimation of these models is done using default settings in MATLAB or default choices proposed by their respective creators.\footnote{As an example, the penalty parameter in the Elastic Net is estimated using 10-fold cross-validation.} Exact details of these algorithms and their default settings is provided in the Online Supplement.
In terms of statistical properties, all these models cover a wide spectrum of forecasting specifications. The AR(2) is a standard benchmark in economic time series forecasting, and typically performs better than a random walk (which is the benchmark for financial data). Its time-varying parameter counterpart, our second model on the list, allows for proxying for similar specifications that have been shown to forecast inflation well, see StockWatson2007 and Bauwensetal2015. Extracting the first few principal components (factors) is possibly the most popular way of representing parsimoniously the information in a large dataset, see StockWatson2016. A naive factor model uses least squares estimation on a model that has the first five principal components as exogenous predictors, while a second factor model replaces OLS with the Bagging algorithm of Breiman1996 that allows to select the “best” factors in a static way. Next the Dynamic Model Averaging (DMA) algorithm described in KoopKorobilis2012 as well as our VBDVS algorithm allow to implement dynamic variable selection in a TVP setting using the same first five principal components. The Gaussian Process Regression is a very flexible nonparametric method that allows us to understand whether inflation is better described by time-varying parameters or some more complex form of nonlinearity. Moving on to models with 60 factors, we have to drop many previous specifications for computational reasons.\footnote{For example, DMA cannot scale up to these large dimensions, Gaussian Process Regression becomes overparametrized, and Bagging becomes numerically unstable in some periods of the forecasting exercise.} For that reason we use the SSVS algorithm of GeorgeMcCulloch1993, which can be thought of as the static equivalent of our VBDVS algorithm. The Elastic Net of ZouHastie2005 is a popular penalized likelihood estimator for high-dimensional data. Finally, our VBDVS algorithm is also estimated with a larger number of factors to find out whether its dynamic shrinkage properties are useful relative to the naive selection of the first five factors. Finally, we estimate models using all 443 exogenous predictors. The Elastic Net is again on the list, and we also include Partial Least Squares (PLS) regression. PLS is similar to principal component analysis, with the main difference being that factors are extracted with reference to the variable to be predicted. Principal components instead only explain the variability in the exogenous predictors, and it may be the case that they do not carry predictive information for the predicted variable. Finally, our VBDVS algorithm is applied to this full model with all predictors.
In terms of the prior choices used when forecasting with our VBDVS algorithm, these are based on Prior 3 described in the previous subsection, see (ref). We only adapt how “aggressively” we shrink based on the total number of predictors in each model. For model VBDVS/FAC5 we set $h_{0}=1$, for VBDVS/FAC60 we set $h_{0}=12$ and for VBDVS/X we set $h_{0}=100$.
We forecast $h=1,2,4$ and $8$ quarters ahead. We use 50% of the sample as our initial estimation period which, for example, for $h=1$ translates to using data for the period 1960Q4-1989Q2 in order to forecast 1989Q3. We then add one new observation to the estimation sample and forecast $h$-step ahead, until the full sample is exhausted. Since all models that have predictors rely on the direct forecasting regression (ref), for comparability we produce direct AR(2) forecasts as a special case of this equation with no predictors.\footnote{The alternative would be to specify an AR(2) model linking $y_{t}$ with $y_{t-1}$ and $y_{t-2}$ and then iterate the process $h$ periods ahead, a procedure also known as iterative forecasting. By using direct AR(2) forecasts as the benchmark we can explicitly assess the exact contribution of various models that introduce exogenous predictors.} We measure forecast accuracy using the mean squared forecast error (MSFE) and the average log-predictive likelihood (ALPL). The first measure is the square of the forecast error (difference between forecast and real value of $y_{t+h}$) averaged over the out-of-sample evaluation period, while the second measure is calculated as the logarithm of the predictive distribution evaluated at the observation $y_{t+h}$ and also averaged over the out-of-sample evaluation period; see Bauwensetal2015 for more details on these two metrics.
(ref) present the MSFEs and ALPLs for GDP deflator, PCE deflator, CPI and Core CPI, for all competing models and all considered forecast horizons. To be precise results for the benchmark AR(2) are the values of the MSFE and ALPL statistics, while results for all other models are relative to those for the AR(2). For the MSFE this means calculating the ratio such that a number lower than one means that a certain model performs better than the AR(2). For the ALPL relative quantities are obtained as the spread from the ALPL of the AR(2) (i.e. the logarithm of the ratio) such that positive numbers indicate that a certain model performs better than the AR(2).
The immediate message from these tables is that the VBDVS/X is the model that performs best, especially when looking at point forecast evaluation (MSFEs) for $h=2,4,8$. In terms of density forecasts, VBDVS/X and VBDVS/FAC60 are jointly the best performing specifications. While VBDVS/FAC5 is also doing well in longer horizons, this model is always underperforming the TVPAR, that is, the TVP model that doesn't consider any predictors.
How can we explain these results? There are various stylized facts we can derive from the information in these tables. Our discussion here focuses on point forecasts (MSFE criterion), due to the fact for that metric the picture is much clearer. First, time variation seems to matter a lot, especially in the long-run. TVPAR, DMA/FAC5, and the three VBDVS specifications can improve dramatically over their constant parameter counterparts, regardless of whether these consider exogenous predictors or not. Are exogenous predictors important for forecasting? The answer depends on the variable to be forecast, the horizon considered, as well as the way each model specification utilizes the predictors. For example, for GDP deflator for $h=8$ the differences in MSFE between VBDVS/X (TVP model with all available predictors) and TVPAR (TVP model with no predictors) is vast, suggesting that not only time-variation is important but also the information in exogenous predictors. However, looking at all constant parameter models with exogenous predictors, whether these predictors are observed or enter each regression via factor methods, all these methods struggle to beat the simple AR(2). This suggests the argument in the Introduction about pockets of predictability. For that reason, DMA (which is the best performing model for $h=1$ and $h=2$) and the three VBDVS specifications perform very well, with VBDVS/X providing the most dramatic improvements for $h=8$ when at the same time ELN/X and PLS/X perform 24% and 39% worse than the benchmark AR(2).
For the next two inflation variables (PCE deflator and total CPI) a large number of predictors does seem to be important in the short-run, but in the long-run it looks like the largest contribution in forecasting accuracy is due to time-variation in parameters. For example, for PCE deflator and total CPI, for horizons $h=1,2$, ELN/X seems to be performing much better than the AR and the TVPAR specifications. However, for $h=4,8$ the TVPAR overtakes substantially both the AR and ELN/X specifications. While the VBDVS/X is still the best performing model for $h=4,8$, its differences to the TVPAR are statistically much smaller compared to the differences of these two models when forecasting GDP deflator. In any case, whether predictors are important or not, the VBDVS algorithm seems to be doing a very good job in shrinking irrelevant coefficients and making sure that there is not overfitting -- if there was, the VBDVS/X forecasts would be inferior to those from the TVPAR.
Finally, for core CPI all methods struggle to beat the simple AR for very short-run forecasts. The VBDVS/FAC60 and VBDVS/X models do so marginally, while many others perform as much as 150% worse than the benchmark. For longer horizons all constant parameter models continue to underperform, however, the TVP models seem to provide the most dramatic improvements, with the VBDVS/X improving almost 60% over the benchmark. Combined with the observation that the differences between the three VBDVS specifications and the TVPAR are minimal, it looks like that exogenous predictors are not relevant for core CPI. Since core CPI is based on the total CPI by removing its most volatiles components (food and energy), it might be the case that this variable is basically a random walk and even a simple time-varying intercept model (that is, a local level model as in StockWatson2007) would forecast this variable well.
It is harder to extract stylized facts for inflation forecasting based on ALPLs. This is because this metric is based on all the features of the predictive density, that is, all its moments and not just the mean. Given that predictive densities can differ a lot between specifications (e.g. they can multimodal in time-varying parameter models), it is not possible to attribute differences in ALPLs to specific modeling assumptions. However, a clear pattern that emerges is that predictors do help to improve predictive density forecasting relative to the simple AR benchmark, but the largest gains overall are achieved by time-varying parameter models. In all these comparisons the VBDVS/X is the clear winner showing that, even though this is a heavily parametrized model and could easily produce erroneous forecasts, our algorithm ensures sufficient penalization and impressive forecasting gains.
\addcontentsline{toc}{section}{\refname}
\singlespace \setcounter{page}{1}
While all technical details regarding our methodology are provided in detail in the paper, we have skipped details for the numerous competing algorithms used in the Monte Carlo and empirical exercises.