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.
83,134 characters · 9 sections · 0 citation commands
High-dimensional macroeconomic forecasting using message passing algorithms
\\ University of Glasgow}
Keywords: high-dimensional inference; factor graph; Belief Propagation; Bayesian shrinkage; time-varying parameter model
JEL Classification:\ C11, C22, C52, C55, C61
\doublespacing
As a response to the increasing linkages between the macroeconomy and the financial sector, as well as the expanding interconnectedness of the global economy, empirical macroeconomic models have increased both in complexity and size. For that reason, estimation of modern models that inform macroeconomic decisions -- such as linear and nonlinear versions of dynamic stochastic general equilibrium (DSGE) and vector autoregressive (VAR) models -- many times relies on Bayesian inference via powerful Markov chain Monte Carlo (MCMC) methods.\footnote{See Herbst and Schorfheide (2015) and Koop and Korobilis (2010) for detailed discussion of Bayesian computation in DSGE and VAR models, respectively.} However, existing posterior simulation algorithms cannot scale up to very high-dimensions due to the computational inefficiency and the larger numerical error associated with repeated sampling via Monte Carlo; see Angelino et al. (2016) for a thorough review of such computational issues from a machine learning and high-dimensional data perspective. In that respect, while Bayesian inference is a natural probabilistic framework for learning about parameters by utilizing all information in the data likelihood and prior, computational restrictions might make it less suitable for supporting real-time decision-making in very high dimensions.
This paper introduces to the econometric literature the framework of factor graphs (Kschischang et al., 2001) for the purpose of designing computationally efficient, and easy to maintain, Bayesian estimation algorithms. The focus is not only on “faster” posterior inference broadly interpreted, but on designing algorithms that have such low complexity that are future-proof and can be used in high-dimensional econometric problems with possibly thousands or millions of coefficients. While a graph, in general, is a structure that allows the representation of objects that are related in some sense\footnote{The most popular use of graphs in economics is to represent networks of agents, banks, social networks etc; see Jackson (2008).}, a factor graph representation of a high-dimensional vector of model parameters, in particular, depicts how each of its scalar elements is connected with each other based on the functional form of their joint posterior distribution. As a result, the factor graph representation provides a visual tool for the decomposition of a high-dimensional joint posterior distribution into smaller, tractable parts. By doing so, factor graphs can be used to design parallel versions of MCMC algorithms, as well as efficient iterative algorithms called message passing algorithms -- the latter being the concept of interest in this paper.\footnote{Message passing algorithms are dynamic programming methods designed for efficiently performing large computations by distributing calculations among a number of simpler processors. Readers working with High-Performance Clusters (HPC) might be familiar with the related concept of message passing interface (MPI) which is a standardized means for exchanging data/commands between multiple processors in a computer cluster.}
Having the factor graph as the starting point, interest lies in an estimation strategy called the sum-product algorithm which is not well known in mainstream statistics, despite the fact that it is computationally powerful (Wand, 2017, p. 137-138). The sum-product algorithm is a general rule in factor graphs that allows to iteratively approximate marginal (posterior) distributions. When applied to a parametric problem with arbitrary likelihood and prior functions, the so-called Generalized Approximate Message Passing (GAMP) algorithm introduces further Gaussian and quadratic approximations to the possibly complicated expressions derived by the sum-product iterative algorithm. Proposed by Rangan (2011), GAMP is an extension of the popular Approximate Message Passing (AMP) algorithm of Donoho et al. (2009). The GAMP algorithm has desirable properties, namely, high-dimensional scalability, parallelizability, and effortless maintenance. Therefore, the first task of this paper is to analyze the concept of message passing algorithms in general; simplify the jargon stemming from signal processing, computing science, and similar literatures that have introduced such algorithms; and show how GAMP, in particular, can lead to efficient posterior inference in very high-dimensions.
At the same time, a second important task is to provide compelling evidence that the proposed algorithm is relevant for modeling macroeconomic variables. For that reason, I utilize a regression model setting with time-varying coefficients, stochastic volatility, and exogenous predictors. Regression models featuring time-varying parameters (TVPs) have been popular in economics at least since the seminal work of Cooley and Prescott (1976). More recently, there has been a systematic effort to introduce efficient MCMC algorithms for flexible estimation and shrinkage in Bayesian TVP models; see Belmonte et al. (2014), Chan et al. (2012), Giordani and Kohn (2008), Groen et al. (2013), Kalli and Griffin (2014), Koop and Potter (2007), Kowal et al. (2018), Nakajima and West (2013), Ro\v{c}kov\'{a} and McAlinn (2018) and Stock and Watson (2007) among others. These are examples of carefully designed MCMC algorithms that result in flexible joint modeling of structural instabilities and parameter shrinkage, but that may not be scalable to very high dimensions due to their reliance on repeated sampling via Monte Carlo.
As a consequence, a novel empirical contribution introduced in this paper is to estimate a time-varying parameter regression model by using an observationally equivalent high-dimensional static regression form, and to address computational concerns by using message passing inference. With $T$ observations and $p$ predictors, the TVP model can be written as a static regression with the same $T$ observations but $(T+1)p$ covariates -- where the product $(T+1)p$ can easily be in the order of tens of thousands in standard macroeconomic applications. This static representation of the time-varying parameter model is anything but new, however, its estimation in the past has been exclusively tackled by specifying an additional hierarchical random walk (or some times stationary autoregressive) model for all time-varying parameters. This hierarchical form allows for inference using state-space methods and at the same time it can be interpreted as an informative shrinkage prior that makes estimation of this high-dimensional problem feasible. Instead I propose to completely drop this “random-walk prior” and the resulting state-space representation, and estimate the time-varying parameter model as a high-dimensional static regression with the assistance of a flexible Bayesian hierarchical shrinkage prior inspired by Tipping (2001). That way, by casting the TVP regression model into equivalent static form, standard shrinkage principles can be used in order to determine by how much coefficients evolve over time, or whether their value is zero and they are completely irrelevant. Most importantly, the use of the low-complexity GAMP algorithm ensures that the static form of the TVP regression with $(T+1)p$ covariates can be estimated quickly. The benefits of this algorithm and modeling strategy are illustrated using a forecasting exercise for monthly U.S. inflation that extends Stock and Watson (1999) to the TVP setting. The static form of the TVP regression estimated with GAMP is contrasted with powerful but slow MCMC algorithms for TVP models, such as Chan et al. (2012) and Kalli and Griffin (2014). The proposed approach, by incorporating a larger number of predictors and by shrinking coefficients flexibly, does perform significantly better compared to competitors in out-of-sample forecasting.
In the next section I introduce the general framework of factor graphs on random variables (parameters) and with the help of a toy example I show how this framework allows for efficient calculation of marginal distributions. Next, in Section 3 I introduce the TVP regression setting, rewrite the likelihood in static regression form and specify a shrinkage “sparse Bayesian learning” (SBL) prior. Under the given functional forms for the likelihood and prior, I proceed to derive a GAMP algorithm for this particular problem. In Section 4 the benefits of the proposed high-dimensional modeling approach are evaluated in a forecasting exercise for U.S. price inflation. Section 5 concludes the paper.
A factor graph represents the way a global function of several variables can be decomposed into a product of simpler functions (“factors”). Consider a generic example with discrete random variables $x=\left(x_{1}, x_{2}, x_{3} \right)$ and a joint mass function $p$ that we can decompose, say, as
where $f_{a},f_{b},f_{c},f_{d}$ are the factors that have known functional forms.\footnote{In the next section, the discrete random variables $x$ are replaced by continuous model parameters, and the factors/functions are conditional or marginal probability distributions over these parameters.} This simple example can be depicted using the factor graph of (ref), where circles denote the place of random variables in the graph and filled boxes denote the factors/functions.\footnote{In graph theory, symbols like the boxes and the circles in this example are called nodes or vertices. Nodes that depend to each other are connected with a solid line, and each connected pair of nodes is called an “edge”. }
Consider now calculation of the marginal distribution of $x_{i}$. This is a computationally demanding task due to the fact that it involves integration (summation, in the discrete variable case) over all variables other than $x_{i}$
where $x \setminus x_{i}$ denotes the set $x$ with the element $x_{i}$ removed. As an example, if the variables in $x$ have two states (e.g. they are binary variables), then the above sum would only require $2^3$ operations. However, for high number of states and/or variables computational requirements proliferate substantially. Nevertheless, if $p(x_{1},x_{2},x_{3})$ is replaced with the expression in (ref) it can be seen that not each variable is coupled to every other one, and this feature can be exploited in order to simplify the summation. For example, in the case of variable $x_{1}$, (ref) depicts that it is directly connected to $x_{2}$ and the factors $f_{a}(x_{1})$ and $f_{d}(x_{1},x_{2})$, but it is only indirectly connected to $x_{3}$ and the remaining factors. Put differently, we can simplify (ref) via identity (ref) as follows
The second line of the equation above implies less algorithmic operations compared to the expression in the first line.
It should be clear at this point that the role of the factor graph representation is to allow to pin down the full path of influence that each variable $x_{i}$ exerts on other variables. As a consequence, by having this path of influence, only the required factors $f_{j}$ can be used when calculating marginal distributions, which increases computational efficiency. This is where the concept of message passing formalizes such an efficient procedure for computing marginals. Each variable node passes messages to the next variable, where these messages are real-valued functions showing the influence that this variable exerts on all other variables. In the remainder of this Section message passing inference is introduced and the sum-product algorithm is derived, such that simplifications similar to the ones in equations (ref)-(ref) are formalized mathematically. Subsequently, in Section 3 the results of this toy example with three discrete random variables (parameters) can be generalized to a high-dimensional regression setting with possibly millions of parameters. More detailed introductions to these concepts can be found in popular machine learning textbooks, such as Barber (2012) and Bishop (2006). A recent introduction of message passing inference in factor graphs from a statistician's perspective is provided in Wand (2017).
Denote with $\mu_{x_{i} \rightarrow f_{j}}$ the message sent from variable $x_{i}$ to function $f_{j}$, and with $\mu_{f_{j} \rightarrow x_{i}}$ the message sent from factor node $f_{j}$ to variable node $x_{i}$, where $i=1,2,3$ and $j=a,b,c,d$ in our simple example with three variables and four factors. The message sent from variable $x_{i}$ to factor node $f_{j}$ is equal to the product of all messages arriving to node $x_{i}$ except from the message coming from the target node $f_{j}$:
where $N(x_{i})$ is the set of neighboring (factor) nodes to $x_{i}$. Similarly, the message sent from factor node $f_{j}$ to variable node $x_{i}$ is given by the sum over the product of the factor function $f_{j}$ itself and all the incoming messages, except the messages from the target variable node $x_{i}$:
Due to the form of the equation above, algorithms that are designed to iterate between (ref) and (ref) are called sum-product algorithms; see also respective equations for the regression model in the next section.
In the special case where $x_{i}$ is an external node (as is the case with $x_{1}$ and $x_{3}$ in this example) it holds that $\mu_{x_{i} \rightarrow f_{j}}=1$. Similarly, if $f_{j}$ is an external factor node (see $f_{a}(x_{1})$ and $f_{d}(x_{3})$ in (ref)) it holds that $\mu_{f_{j} \rightarrow x_{i}}=f_{j}\left(x_{i}\right)$. Equations (ref)-(ref) define the iterations of the so-called sum-product algorithm (also called Belief Propagation; see Pearl, 1982), that allows calculation of marginal distributions (also called “beliefs” in computing science and the Bayesian networks literature). Upon convergence, it can be shown\footnote{It is beyond the scope of this paper to derive and prove the algorithm, and the reader is referred to the excellent machine learning books of Barber (2012) and Bishop (2006).} that
that is, the marginal distribution of variable $x_{i}$ is simply the product of all messages received only from factor nodes that are connected to $x_{i}$.
Consider for example calculation of $p\left( x_{2} \right)$. Starting from the left of the graph, the messages emitted to node $x_{2}$ are:
where the first identity holds because $f_{a} \left(x_{1}\right)$ is an external factor node, the second identity is a result of equation (ref), and the third identity is a result of (ref). Similarly, the messages that arrive to $x_{2}$ stating from the right of the graph are
where again the first identity results from the fact that $f_{d}\left( x_{3},x_{4} \right)$ is an external factor node, the second results from equation (ref) and the third from equation (ref). Therefore, the marginal distribution of $x_{2}$ is now
Using similar arguments we can derive $p\left(x_{1} \right)$ and $p\left( x_{3} \right)$.
In this particular example, the formula derived in (ref) might seem redundant as for a wide class of distributions $p\left( \bullet \right)$, one can simply calculate the marginal distribution of $x_{2}$ using numerical integration. However, in high dimensions with many random variables, the sum-product rule can provide us with scalable and parallel posterior inference algorithms that can be several times faster compared to conventional algorithms that iterate sequentially (e.g. Gibbs sampler). It can be shown that the sum-product (Belief Propagation) algorithm is a special case of the more general expectation propagation algorithms that have been very popular in Bayesian machine learning; see Vehtari et al. (2018). Finally, note at this point that there is no mention about how to approximate the summations in (ref), which will not necessarily be tractable. Given the sum-product formula, there are several algorithms that would allow for the approximation of the required messages which are functions of the factors $f_{j}$. For example, Wand (2017) develops message passing inference inspired by the variational Bayes method. In the next section I adopt a recently developed algorithm (Generalized Approximate Message Passing) that performs Normal approximations to the functions implied by the sum-product iterations.
The starting point is the following time-varying parameter (TVP) regression with stochastic volatility of the form
subject to an initial condition for $\beta_{t}$ at $t=0$ (denoted as $\beta_{0}$), where $y_{t}$ is the $t^{th}$ observation on the variable of interest, $t=1,...,T$, $x_{t}$ is a $1 \times p$ vector of predictors (possibly including lags of $y_{t}$), $\beta_{t}$ is a $p \times 1$ vector of coefficients, and $\varepsilon_{t} \sim N \left( 0 , \sigma_{t}^2 \right)$ with $\sigma_{t}^2$ the time-varying variance parameter. It is desirable to estimate the initial condition in this model, rather than assume it is knonw. For that reason, following Fr\"{u}hwirth-Schnatter and Wagner (2010), this model can be written using an equivalent non-centered parametrization that allows to split the parameter $\beta_t$ into a part that is constant (which is equivalent to its initial condition $\beta_{0}$), and an “add-on” time-varying part with initial condition fixed to zero. The equivalent specification is
where now $\widetilde{\beta}_{t}$ has initial condition zero and it holds that $\beta_{t} = \widetilde{\beta} + \widetilde{\beta}_{t}$. As shown in Belmonte et al. (2014) this parametrization allows to use shrinkage priors to determine whether a variable has constant coefficient (by only shrinking the time-varying part), or it is completely irrelevant for modeling $y$ (by shrinking both the constant and time-varying parts to zero). More details of this approach are provided in the Online Appendix, Section D.1.
The TVP regression can be written in the following equivalent static regression form
where $y = \left[y_{1},...,y_{T} \right]^{\prime}$ and $\varepsilon = \left[\varepsilon_{1},...,\varepsilon_{T} \right]^{\prime}$ are column vectors stacking the observations $y_{t}$ and $\varepsilon_{t}$ respectively, $\beta = \left[\widetilde{\beta}^{\prime},\widetilde{\beta}_{1}^{\prime},...,\widetilde{\beta}_{T}^{\prime} \right]^{\prime}$ is a $(T+1)p \times 1$ vector, and
is a $T \times (T+1)p$ matrix. It is evident that the first $p$ columns of $\EuScript{X}$ specify a constant parameter regression and its remaining columns add “time-dummies” to that regression. The Gram matrix $\left( \EuScript{X}^{\prime} \EuScript{X} \right)$ is of rank $T$ and the $q = (T+1)p$, in total, regression coefficients in (ref) cannot be estimated with OLS. For that reason, following a long-standing tradition in engineering, economists tend to assume that $\beta_{t}$ (similarly for $\widetilde{\beta}_{t}$ in the non-centered parametrization) typically follows a random walk of the form $\beta_{t} = \beta_{t-1} + \eta_{t}$, where $\eta_{t} \sim N\left(0,Q\right)$ for some $p \times p$ symmetric, positive-definite covariance matrix $Q$. This random walk regression for $\beta_{t}$ allows to write the full time-varying parameter regression model in familiar state-space form, and also provides the additional information needed to estimate $\beta_{t}$ using data $y$ and $\EuScript{X}$. By doing so, estimation typically relies on Markov chain Monte Carlo methods by means of a simulation smoother; see Primiceri (2005) for a representative example. From a Bayesian point of view this additional information can be viewed as a conditional hierarchical prior of the form $p\left(\beta_{t} \vert \beta_{t-1} \right) \sim N\left(\beta_{t-1}, Q\right)$ that provides appropriate level of shrinkage. Put differently, equation (ref) alone can be seen as an ill-posed problem where OLS does not have a unique solution and regularization is imperative for estimation.
In this paper I adopt this shrinkage view of the time-varying parameter regression model and propose an alternative inference strategy. That is, inference is done without reference to the useful but rather informative and subjective conditional hierarchical prior for $\beta_{t}$ given $\beta_{t-1}$ outlined above. Instead, the time-varying parameters are recovered by estimating directly equation (ref) using data-based hierarchical shrinkage priors. In particular, I follow Tipping (2001) and define the following independent hierarchical prior for each element $\beta_{i}$ of the vector $\beta$, $i=1,2,...,(T+1)p$,
This conditionally Normal prior for $\beta_{i}$ and Gamma prior for the precision parameter $\alpha_i$ is a scale mixture of Normal representation of a Student-t prior. Tipping (2001) calls this heavy-tailed prior a sparse Bayesian learning (SBL) prior, and I adopt this name henceforth; see also Korobilis (2013) for a detailed explanation why such hierarchical priors have good shrinkage properties. I follow Tipping (2001) and present all empirical results using the uniform hyperpriors (over a logarithmic scale) $\underline{a}=\underline{b}=1\times 10^{-10}$.
Two additional comments are in order regarding this time-varying parameter regression. First, the number of columns of $ \EuScript{X}$ is $q = (T+1)p$, therefore, the number of coefficients grows rapidly. For example, with 700 monthly observations and only 100 predictors, we end up with 70,100 regression coefficients. As a consequence, it is imperative to choose a fast estimation algorithm that approximates the parameter posterior, and this is where the scalability of message passing algorithms comes into play. Second, there is no mention yet of inference on $\sigma_{t}^{2}$, as this issue is covered later in this section after the GAMP inference algorithm is outlined. In a nutshell, estimation of stochastic volatility $\sigma_{t}^{2}$ also follows the same shrinkage principles defined for $\beta_{t}$. That is, it is shown that we can write estimation of $\sigma_{t}^{2}$ as a high-dimensional regression problem, without having to assume any kind of first-order Markov dependence to $\sigma_{t-1}^{2}$.
At this point we have all the necessary ingredients in order to cast the static form of the time-varying parameter regression in equation (ref) into a factor graph form.\footnote{For the sake of brevity, notation for prior, posterior and likelihood distributions is generic, that is, there is no reference to their exact functional forms. Exact details and parametric formulas can be found in the Online Appendix, Section B.} Consider first an independent (but not necessarily i.i.d) prior for $\beta$, denoted $p\left( \beta \right) = \prod _{i=1}^{q} p\left( \beta_{i} \right)$, and the resulting posterior from Bayes Theorem
The exact marginal posterior of $\beta_i$, $i=1,...,q$ is of the form
where $d\beta_{j \neq i}$ denotes integration over the whole set of $q-1$ parameters $\beta_j$ for $j \neq i$. Therefore, the formula above requires integration over a $(q-1)$-dimensional integral, a numerical problem that can become computationally infeasible for a high-dimensional vector $\beta$.
We can now call the framework of factor graphs in order to factorize efficiently the marginal posteriors of $\beta$. The factor graph representation of the regression model is depicted in Figure (ref). Based on this figure, the marginal posterior of $\beta_i$, presented in equation (ref), can be defined as the product of incoming messages at node $\beta_i$ in the graph
Similar to equation (ref) in the example of Section 2, the message $\mu_{ p(\beta_i) \rightarrow \beta_i }$ is an external factor node and for that reason it is equal to the prior $p \left( \beta_i \right)$. Generalizing the example sum-product rule derived in equations (ref) - (ref) of the previous section, we can write the messages from $p(y_{t} | \beta)$ $\forall t$ to $\beta_{i}$ using the following expression
In the decomposition above, the message from node $\beta_j$ to function (factor) $p \left( y_t | \beta \right)$ is the product of all incoming messages to node $\beta_i$, excluding the message coming from $p \left( y_t | \beta \right)$ itself
We can see in equations (ref)-(ref) that in order to obtain the message $\mu_{ p(y_{t} | \beta) \rightarrow \beta_i }$ we need $\mu_{ \beta_j \rightarrow p(y_{t} | \beta) }$ and vice-versa. Therefore, one can simply update both equations iteratively using the following iterative sum-product scheme
where the superscript $(r)$ denotes the $r^{th}$ iteration of the algorithm. In graphs with a tree structure, one iteration of the algorithm above will always recover the exact marginal posteriors for the parameters $\beta_i$. In a factor graph with loops there are no guarantees that the sum-product rule will converge to a good fixed point. However, the sum-product rule can still achieve a good approximation and this is the reason why it is used extensively in applications of coding theory, machine vision, and compressive sensing that have a loopy graph representation (Mooij and Kappen, 2007). Translating these facts into familiar jargon for the static regression in equation (ref), algorithmic convergence is achieved if the correlation of right-hand side predictors is not excessively high. If this is not the case, the joint posterior of the coefficients $\beta$ might also be highly correlated, which would make inference solely based on the marginal posteriors $p(\beta_{i})$ less accurate. In our benchmark time-varying parameter regression in (ref), correlation is by default not excessively high due to the fact that the Gram matrix $ \EuScript{X}^{\prime} \EuScript{X}$ has a certain block-diagonal structure that allows for a general sparse correlation pattern -- even if within a given block correlation may be high. In the empirical application, predictor variables are mainly principal components or lags thereof, such that correlation within each block is also low. Finally, note that the specific time-decomposition of the likelihood function does not accommodate autoregressive and general time-series models, where the likelihood at time $t$ may be written conditional on past observations. In the empirical application it is found that, despite this approximation, autoregressive coefficients are recovered accurately.\footnote{A simulation exercise in the Online Appendix, Section C.3, generating artificial data from an AR(4) model, also verifies that the proposed GAMP algorithm performs well even if the likelihood function is not i.i.d. Another assumption that affects performance of GAMP is that $\EuScript{X}$ is mean-zero Gaussian; see the discussion in Al-Shoukairi et al. (2018) and references therein. In a time series context this means that GAMP will have better convergence when right hand-side predictors are strictly stationary, although the use of weakly stationary predictors is not excluded.}
While the core of any message passing algorithm is fully described by the sum-product iterations, deriving the exact functional form of the messages in equations (ref) and (ref) under the regression likelihood and the Student-t hierarchical prior implies that cumbersome integrations might be necessary. The GAMP algorithm introduces certain Gaussian approximations to the sum-product iterations. Unlike Laplace approximations, that is, Gaussian approximations to parameter posteriors that many times can be poor, the GAMP approximation is fully based on asymptotic results that make it more reliable as the number of predictors grows large. First, when $q \rightarrow \infty$ a central limit theorem (CLT) postulates that the messages $\prod_{j=1,j \neq i} ^{q} \mu_{ \beta_j \rightarrow p(y_{t} | \beta) }$ can be approximated by a Gaussian distribution with respect to the uniform norm.\footnote{This is a result of the Berry-Esseen central limit theorem which states that a sum of random variables converge to a Gaussian density; see a proof of this theorem in Donoho et al. (2011). Given that the sum-product equations involve products of random variables, rather than sums, derivations of GAMP based on this central limit theorem typically proceed by taking logarithms of equations (ref)-(ref). The marginal posterior $p \left( \beta_i | y \right)$ is then recovered by performing an exponential transformation of the log messages, and by normalizing so that the posterior integrates to one; see the Online Appendix, Section A, for details.} This result means that messages in (ref) can be represented to be proportional to a Gaussian distribution. A second approximation involves taking the Taylor-series expansion of terms in the messages, so that the first two moments (mean and variance) of $p \left( \beta_i | y \right)$ can be obtained analytically up to the omission of $O\left(1/q \right)$ terms. Exact derivation of these approximations involves many tedious steps and transformations, and the reader is referred to the Online Appendix for more details. What is important to stress at this point is that both the CLT and Taylor-series approximations vanish as $q \rightarrow \infty$ with $q/T \rightarrow \delta$ for some constant $\delta$; see Rangan (2011) and Rangan et al. (2016) for more details. This is an example of the “blessing of Big Data” -- rather than the “curse of dimensionality” embedded in many traditional estimation algorithms -- as the GAMP algorithm fully facilitates the large $q$ asymptotics.
Deriving the GAMP algorithm involves several steps and lengthy proofs which are left for the Online Appendix. The final product of all the approximations to the two sum-product update equations (ref) - (ref), is a simple iterative algorithm that provides an approximation to the mean and variance of $p \left( \beta_{i} \vert y \right)$. The algorithm iterates through computationally trivial scalar multiplications and additions that result in worst case algorithmic complexity of $\mathcal{O}(Tq)$. That is, estimation of the marginal parameter posterior distribution does not involve costly operations such as high-dimensional integration or inversion of large matrices. This feature implies that the algorithm can handle regressions with an excessively large number of predictors with the same ease it can handle smaller regression models. Convergence is achieved when the difference between estimates of the posterior mean of $\beta$ between two consecutive iterations is below a pre-specified tolerance level. Other parameters can be updated by combining the GAMP algorithm with EM updates.\footnote{See Al-Shoukairi et al. (2018) and Zou et al. (2016) for examples of how to derive EM updates for prior hyperparameters.} This feature is explained in the Online Appendix, where it is shown how to update the hyperparameter $\alpha_{i}$ introduced in the hierarchical prior of equation (ref).
A sketch of the algorithm is provided in Algorithm (ref). This is a simplified version that focuses on estimation of $\beta$ by assuming that the regression variance and prior hyperparameters are all known and fixed. Following the analysis in Section 2, the algorithm can be split into two steps: i) evaluating all messages that leave each variable node $\beta_{j}$ (output), and ii) evaluating all messages that arrive at each variable node $\beta_{j}$ (input). The final product is estimates of the posterior mean and variance of $\beta_{j}$ which are denoted as $\widehat{\beta}_{i}$ and $\widehat{\tau}_{i}^{\beta}$, respectively. At the core of the calculation of the posterior mean and variance are the scalar functions $g_{in}$ and $g_{out}$. Derivation of the exact form of these two functions depends on the form of the prior distribution and the likelihood. Online Appendix, Section B, provides a detailed algorithm in the case of the regression likelihood in equation (ref) and the prior in (ref)-(ref). In any case, Rangan (2011) shows that regardless of the form of the nonlinear scalar functions $g_{in}$ and $g_{out}$, the worst-case complexity of the GAMP algorithm is not affected and is always $\mathcal{O}(Tq)$.
The algorithm above assumes a known regression variance, e.g. normalized to be one. Of empirical interest is the derivation of an update rule for the variance parameter when this is both unknown and time varying. Here I propose a novel, computationally trivial estimator of the variance that builds on approximations used in the Bayesian stochastic volatility estimator of Kim et al. (1998). First, we write the regression model in (ref) in the following form
where $\Sigma$ is a $T \times T$ diagonal matrix with the time-varying standard deviations $\sigma_{t}$ on its main diagonal. Subsequently, conditional on knowing $\beta$ by means of some estimate $\widehat{\beta}$, we can re-write the above model as
where $diag\left(\Sigma\right)^{2}$ is a $T \times 1$ vector with elements $\sigma_{t}^{2}$ $\forall t \in [1,T]$, and variables with a $\widetilde{\bullet}$ denote quantities in log-squares. In particular, the distribution of $\widetilde{v}$ is $\log-\chi^{2}$ with one degree of freedom. Following Kim et al. (1998) we can approximate this with a mixture of seven Normal distributions with means $\mu_{i}$, variances $V_{i}$ and component weights $\pi_{i}$, where $i=1,...,7$ and $\sum_{i} \pi_{i}=1$.\footnote{The exact values of $\mu_{i}$, $V_{i}$, $\pi_{i}$ for all seven components is provided in the Online Appendix, Section B.} Then equation (ref) can be replaced with the following set of seven equations
where $u_i \sim N\left(\mu_{i}, V_{i} \right)$. An estimator of the $T \times 1$ vector of log-volatilities is of the form $E_{i}\left(\widetilde{\sigma}^{2} \right) = \widetilde{y} - \mu_{i}$, and the final volatility estimate at time $t$ is
Similar expressions can also be derived for the posterior variance of $\sigma_{t}^{2}$ if desired, for example, when computing the posterior predictive density via simulation. It turns out that the resulting estimate of volatility is similar to the standard stochastic volatility estimator of Kim et al. (1998), but it is much less persistent due to the lack of dependence of $\sigma_{t}^{2}$ on $\sigma_{t-1}^{2}$. More evidence on the excellent properties of this simple estimator of stochastic volatility is provided in the Online Appendix, Section D.1.
Finally, Online Appendix, Section C, provides detailed Monte Carlo evidence on the usefulness of the proposed econometric specification and algorithm. By simulating artificial data from models with various patterns of time-variation in parameters, it is assessed how good the specification in equation (ref), with the assistance of the sparse Bayesian learning prior, is at recovering the true time-varying parameters. At the same time, a second simulation exercise shows the ability of the GAMP algorithm with shrinkage prior to perform high-dimensional shrinkage even in cases with more predictors than observations. A final simulation exercise discusses the stability of the GAMP algorithm in models with correlated predictors, and assesses numerically the case where the likelihood function is not i.i.d. While the results of the simulated data exercises suggest that the proposed algorithm provides a reasonable balance between computational speed and estimation accuracy, the next section establishes that the proposed algorithm is also very useful in a forecasting application using real macroeconomic data.
This section describes the set-up and results of a comprehensive forecasting exercise that demonstrates the merits of the modeling approach outlined in the previous section. Most applications of time-varying parameter regressions focus in particular on inflation. Of course, this class of models is flexible enough to provide useful forecasts of any other variable of interest; see Bauwens et al. (2015) for assessing structural breaks in several monthly and quarterly macroeconomic time series. Nevertheless, there is ample evidence that structural breaks in inflation are so evident and complex, such that TVP models are particularly useful for forecasting this variable; see Chan et al. (2012), Groen et al. (2013), Koop and Korobilis (2012), Pettenuzzo and Timmermann (2017) and Stock and Watson (2007) among many others.
The data collected for this exercise are 115 macroeconomic variables from Federal Reserve Economic Data (FRED) of St. Louis Federal Reserve Bank website. The data originally span the period 1959M1 to 2016M6, but the effective sample is smaller after taking stationarity transformations and lags. The stationarity transformations follow standard norms in this literature (see Stock and Watson, 1999) and exact details are provided in the Online Appendix, Section A.
The empirical application builds on the seminal work of Stock and Watson (1999) for forecasting inflation. These authors specify the following benchmark forecasting model
where $\pi_{t}^{h} = (1200/h)log\left(P_{t}/P_{t-h}\right)$ is the $h$-period inflation in the price level $P_{t}$. As Stock and Watson (1999; Section 2) explain in detail the assumption here is that inflation is $I(1)$ while the exogenous variables in $z_{t}$ are $I(0)$. Two modifications of this basic forecasting model are in order. First, as Stock and Watson (1999, 2002) also suggest, the high-dimensional variables $z_{t}$ are replaced by factors $f_{t}$ estimated using principal components. Second, the forecasting equation is enhanced with time-varying parameters and stochastic volatility. The final forecasting model used in this paper is of the form
where $e_{t} \sim N \left(0, \sigma^{2}_{t} \right)$ and $f_{t}$ is a lower-dimensional vector of factors.
The forecasting exercise is run for two measures of inflation, namely the consumer price index for all items (CPIAUCSL) and the personal consumption expenditures price index (PCEPI). The forecast horizons evaluated are $h=1,3,6,12$ which correspond to one-month, one-quarter, one-semester and one-year ahead forecasts, respectively. Following Bauwens et al. (2015) evaluation of forecasts is based on the mean square forecast error (MSFE) for point forecasts, and on the logarithm of the average predictive likelihoods (log APL) for comparing whole forecast densities. Exactly $50\%$ of the sample is used for evaluation of out-of-sample forecasts, leading to a period of $343-h$ months where MSFEs and log APLs are calculated. Note that while estimation entails the spread $\pi_{t+h}^{h} - \pi_{t}$, all forecast evaluations in this Section (see also alternative model in equation (ref)) pertain to $\pi_{t+h}^{h}$.
When applying the proposed GAMP estimation methodology, equation (ref) is estimated using two own lags of the dependent variable, the first 20 principal component estimates of the factors $f_{t}$ (updated recursively using only information up to time $t$) and two lags of these factors (that is, their values in periods $t$ and $t-1$). As explained in the main text, this TVP model can be estimated using GAMP by casting it into the form (ref) by setting $y_{t} = \pi_{t+h}^{h} - \pi_{t}$, $x_{t} = \left[1,f_{t},\Delta \pi_{t}\right]$, $\beta_{t} = \left(\phi_{t,0},\theta_{t}(L)^{\prime}, \gamma_{t}(L)^{\prime}\right)^{\prime}$ and $e_{t+h} = \varepsilon_{t}$. Written in this static form and using all available observations, the proposed empirical model has nearly 30000 regression coefficients and another 700 volatility parameters to estimate. The only input that the GAMP algorithm requires is choice of two scalar prior hyperparameters. For the sparse Bayesian learning prior of equations (ref) - (ref) these hyperparameters are set, as explained in Section 3, to the uniform values $\underline{a} = \underline{b} = 1 \times 10^{-10}$. This approach to estimating the TVP regression of (ref) using GAMP is abbreviated as TVP-GAMP in the results presented next.
The benchmark time-varying regression approach estimated with the GAMP algorithm is contrasted against a range of popular algorithms for inference in models with many predictors and/or stochastic variation in coefficients. The list of competing specifications and estimation algorithms is the following:
All these models collapse to being special cases of the benchmark equation (ref), despite the fact that different specifications might imply various additional assumptions about how the coefficients might evolve over time (whereas TVP-GAMP does not rely on such additional assumptions). All models except for the UCSV have in common an intercept and the two own lags of inflation.\footnote{In order to understand better whether forecast gains can be achieved from specifying a model with many predictors, or with flexible time-variation, or both, I only calculate direct multi-step forecasts from all competing models. That way all algorithms are used to estimate different versions of the same regression with $y_{t+h}$ on the left hand side (for each $h$) and information dated $t$ or earlier on the left hand side. However, iterated forecasts can be computed from models with no exogenous predictors (e.g. TVP-AR or UCSV). Direct forecasts are better when the model is misspecified, while iterated forecasting models in general result in more efficient econometric estimates and sharper predictive densities. Examination of $h=12$ month ahead iterated forecasts from the KP-AR, GK-AR, TVP-AR and UCSV models reveals that these are, most times, slightly inferior to respective direct forecasts in terms of MSFE, but they can be in some cases up to $15\%$ better in terms of average log predictive likelihoods. Iterated forecasting results are not presented here, but they are available from the author.} For those algorithms that rely on shrinkage priors (TVP-GAMP, TVD, TVS, TVP-BMA, and BMA) the intercept and the two lags of inflation are never allowed to shrink by using a noninformative prior on them. Therefore, whenever shrinkage (static or dynamic) is implemented this only applies to the exogenous information in the factors. Exact details of the econometric specifications and prior settings associated with the competing models is provided in the Online Appendix, Section E.
A final note is on computation. All of the competing models listed above are based on estimation using MCMC and in particular the Gibbs sampler. Most of these models were originally developed by their respective authors for forecasting inflation. This is due to the fact that time-varying parameter regressions have consistently been found to be superior for this series. However, even though one would normally expect more breaks to be present in higher frequency monthly inflation, all of these papers estimate their models using quarterly data. This is done for computational reasons. Due to the fact that here these models are estimated for monthly data, I follow Bauwens et al. (2015) and base inference only on 5000 samples from the posterior after a burn-in period of 1000 draws, that is, a total of 6000 MCMC iterations. Convergence criteria suggest that such low number of iterations is sufficient for forecasting, even though it might not be satisfactory for other econometric exercises. Despite the low number of MCMC iterations, computation is quite cumbersome taking several hours for some models. In contrast, it takes only minutes to run the full recursive exercise using the TVP-GAMP model that features both time-varying parameters and the full set of available predictors. The GAMP algorithm not only involves simple scalar computations, but also converges fairly quickly after 10 to 100 iterations. Once convergence is achieved, the first two posterior moments are readily available for further inference, rather than having to store thousands of samples from the posterior of a high-dimensional parameter vector.
The results from this forecasting exercise are presented in Tables 1 and 2, and are very encouraging for the proposed TVP-GAMP method. (ref) shows MSFEs relative to an AR(2) benchmark (with an intercept), such that numbers lower than one signify better performance of a competing model relative to that benchmark AR(2) specification. It can be seen that under the specified regression model, point forecasts from TVP-GAMP dominate alternatives by a substantial amount, both for CPI and PCE inflation. The forecast gains are increasing with the horizon. (ref) shows the logarithm of the average predictive likelihood (log APL), and this metric is quoted as a spread from the log APL of the simple AR(2) specification. Positive values signify better performance relative to the benchmark AR(2). Using this metric, TVP-GAMP is either the top performing model or among the top, for the four forecast horizons and the two measures of inflation.
It is notable that these results contradict the previous claims that time-variation in parameters is important for inflation. The three models with the largest number of predictors, namely BMA and TVP-GAMP, and to a lesser degree TVP-BMA, seem to be improving a lot over time-varying parameter models with no predictors. The results seem to suggest that information in predictors is more important than the specification of time variation in regression parameters. This observation is not undermined by the fact that point forecasts from TVP-BMA are not significant, and that density forecasts from BMA are quite poor relative to TVP-BMA and TVP-GAMP. First, TVP-BMA is overparametrized\footnote{Shrinkage in TVP-BMA is only across predictors, but this model does not restrict the amount of time-variation in parameters.} its point forecast performance is not as good as the more conservative (in terms of time-variation in parameters, not available number of predictors) BMA and TVP-GAMP specifications. Second, when considering density forecasts, BMA is definitely misspecified since it does not allow for stochastic volatility, and it naturally doesn't perform as well as TVP-BMA and TVP-GAMP that allow for changing variance. Therefore, these findings suggest that TVP-GAMP is overall the best model and that its specification is flexible enough to capture both structural change and utilize information in a large set of predictors at the same time. Most importantly, the SBL prior allows to strike a good balance between these two modeling characteristics by removing irrelevant predictors as well as regularizing time variation.
These results are in stark contrast to existing results for TVP models presented in the papers cited above (see e.g. footnotes in Table 1). The culprit is simply the assumption that inflation is I(1) that Stock and Watson (1999) introduce in their seminal paper, and that it is adopted in equation (ref). Once the random walk dynamics are removed from inflation (i.e. inflation gap becomes the dependent variable), the role of time-varying parameters in forecasting becomes less important and the most significant feature is the information included in exogenous predictors. It would be interesting then, as a robustness check, to specify the forecasting regression for inflation using the following form
This equation is more in line with the forecasting model estimated in papers such as Chan et al. (2012), Groen et al. (2013), or Pettenuzzo and Timmermann (2017).
(ref) shows results based on this alternative specification of equation (ref) for CPI inflation only. The left part of the table presents MSFE results, while the right panel presents log APLs. In this case it is evident that the various variants of TVP models considered improve tremendously over the benchmark. As a matter of fact, models such as the KP-AR, TVP-AR and UCSV also improve a lot relative to the constant parameter BMA. Looking at point forecasts and the associated MSFE results, we can observe many differences among TVP models, especially as the forecast horizon increases. For example, the structural breaks KP-AR specification has the lowest relative MSFE for $h=12$ among all models, but the also structural breaks GK-AR specification is among the worst performing (but still much better than the simple AR model). TVD and TVS estimated with the monthly data are not only cumbersome, but also do not perform as well as TVP models with no predictors. In contrast, the TVP-BMA algorithm is performing quite well, even though it still doesn't beat TVP models with no predictors. In this alternative forecasting regression, TVP-GAMP is not the top forecasting model but its performance is still quite good. If it wasn't for the exceptional performance of the KP-AR model, TVP-GAMP would have been a top model for $h=1,3,6$.
When looking at density forecast evaluation the results might not comply with the results for the point forecasts. Still good performing models are the KP-AR and the TVP-AR, but now the BMA and TVP-GAMP beat models such as the UCSV. With such diverse set of flexible models it is hard to pin down which exact features help in point and density forecasts. Nevertheless, for the forecasting regression (ref) it seems that the way time variation in parameters is specified is more important than information in exogenous predictors. Further numerical evidence on the relative forecast performance of some of the competing models, is provided in Online Appendix, Section D.2.
This paper evaluates a new methodology for performing Bayesian inference in high-dimensional regression models. The proposed Generalized Approximate Message Passing (GAMP) is a fast algorithm for approximating iteratively the first two moments of the marginal posterior distribution of a high-dimensional vector of coefficients. It is established how effortlessly the GAMP algorithm can be extended with interesting modeling features such as hierarchical shrinkage priors, time-varying coefficients and stochastic volatility, and many predictors. The benefit of the proposed approach is demonstrated using an inflation forecasting exercise that leads to the recursive estimation of regression models with thousands of covariates. Due to the low algorithmic complexity, GAMP could be generalized to much higher dimensions with millions of predictors/covariates, as it is also trivially parallelizable.
The current study opens up new avenues for research. First, the proposed framework for modeling time-vayring parameters using hierarchical shrinkage priors can be extended in interesting ways. For example, shrinkage estimators/priors that apply on group of coefficients (such as the Group Lasso) can be used in this setting so that coefficients are shrunk either in groups of predictors for a given time period or in groups of consecutive time periods for a given predictor. This is because in the TVP setting the vector of regression coefficients $\beta$ has elements that correspond both to predictor $j$, $j=1,...,p$, but also to time period $t$, $t=1,...,T$. One can think of other shrinkage priors in order to perform a more structured approach to uncovering patterns of time-variation in parameters, such as various pooling priors used in the panel data literature. Finally, the paper proposes the framework of factor graphs for designing efficient algorithms. Many macroeconomic problems currently do not typically involve extensive use of Big Data sets, however, they involve multivariate models with possibly thousands of coefficients, such as VAR, factor, and DSGE models. Bayesian estimation of these models is quite cumbersome, many times relying on linear or nonlinear state-space methods. As empirical macroeconomic models become larger and more complex, factor graph inference could help economists come up with novel efficient algorithms and unveil new features in macroeconomic data.