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.
76,892 characters · 17 sections · 79 citation commands
Bayesian Modeling of TVP-VARs Using Regression Trees
\thispagestyle{empty}
\doublespacing
Econometric models used in macroeconomics have traditionally been linear and homoskedastic. Examples include vector autoregressions (VARs), dynamic factor models (DFMs), and linearized dynamic stochastic general equilibrium (DSGE) models, inter alia. However, in recent decades, there has been growing awareness of the empirical need to allow for parameter change and to relax homoskedasticity assumptions.
Most of the models used to capture such parameter change allow for parameters to vary over time, but at a specific point in time the model remains linear. Structural break models, Markov switching models sims2006were, and time-varying parameter (TVP) models are some prominent examples. TVP regressions and TVP-VARs, in particular, have been highly successful for structural macroeconomic analysis and forecasting dangl2012predictive, giannone2013, koop2013large, belmonte2014hierarchical, bitto2019achieving, korobilis2021high, huber2021inducing, hauzenberger2021fast. In the TVP-VAR literature, parameters are assumed to evolve according to random walks or autoregressive processes.\footnote{It is also common to assume these random walks are independent of one another. But this assumption can translate into overfitting, since coefficients often feature substantial co-movement stevanovic2016common, chan2020reducing.} This also holds true for stochastic volatility processes, the addition of which is the most common way of relaxing the homoskedasticity assumption.
A drawback of these parametric models of parameter change is that they risk mis-specification. That is, there is uncertainty about the specific law of motion driving changes in coefficients and volatilities. This not only includes how coefficients (volatilities) evolve over time but also whether changes in the law of motion depend on other observed factors (resulting in nonlinear interactions). Since policy makers often have a keen interest in how structural factors affect not only observed quantities (such as output or inflation) but also latent quantities that are often highly nonlinear functions of the reduced-form parameters of a TVP-VAR (such as the long-run unconditional mean, impulse response functions, or predictive distributions), this ignorance of possible relations between the parameters of a model and additional covariates is not merely a statistical problem but also has substantial practical relevance. For instance, a policy maker might want to know whether the effectiveness of policy measures depend on the state of the business cycle or whether macroeconomic trends (such as trend inflation or the equilibrium interest rate) depend dynamically on other factors.
The model developed in this paper allows for flexibility in how the parameters and the error volatilities in a VAR evolve over time. The model also allows for a (possibly) nonlinear relationship between a set of covariates, labeled effect modifiers, and the reduced-form parameters of the VAR. As will be explained, this simplifies use and interpretation of the model. Our model builds on those in two recent papers, deshpande2020vcbart and coulombe2020macroeconomy, which also model the dynamic evolution of parameters in a nonparametric manner.\footnote{fischer2022TVP propose a parametric alternative that controls for uncertainty with respect to the form of the state evolution equation of a TVP-VAR.} Both of these papers --- that unlike this paper focus on univariate models --- assume that each coefficient has its own nonparametric law of motion. deshpande2020vcbart approximate the law of motion using Bayesian additive regression trees chipman2010bart, while coulombe2020macroeconomy uses random forests. However, this flexibility might translate into overfitting and a lack of scalability to high dimensions. Given our interest in (potentially) large-scale TVP-VARs with many more parameters, we address these issues by introducing further restrictions and Bayesian shrinkage priors so as both to lessen the risk of overfitting and to maintain computational tractability.
In particular, the main contribution of the paper is the development of a flexible TVP-VAR that has several key model features that are important for inference. First, it allows for any time-variation in the parameters to be driven by a low number of latent factors. Second, it allows for factor stochastic volatility in the reduced-form VAR shocks. Third, it is nonparametric along several key dimensions: latent factors which define the law of motion of parameters (including both the TVPs and the error variance-covariance matrix) follow independent BART models. These BART models, in turn, depend on additional covariates that can be exogenous, endogenous, or latent. Novel shrinkage priors control for overfitting and allow for selection of an appropriate number of BART models.
One advantage of our TVP-VAR, also noted by coulombe2020macroeconomy in the univariate case, is that, by giving a nonparametric treatment to the parameters rather than the variables, our model remains conditionally linear in the parameters. This helps avoid the black-box nature of other nonparametric techniques, such as BART, since it facilitates interpretation of the model. Since we consider TVP-VARs and use them to carry out impulse response analysis, we can, for example, use our model to ask questions such as how impulse responses to structural shocks change if the effect modifiers are altered. This permits scenario analysis that is otherwise not possible (even with fully fledged nonparametric models) or strongly depends on the assumed relationship between the parameters and the effect modifiers.\footnote{Particular examples are hubrich2015financial, aastveit2017economic, caggiano2017estimating, alessandri2019financial, and hauzenberger2021effectiveness.} Our model remains agnostic on this relationship and thus can capture nonlinearities of arbitrary form. In addition, the heteroskedastic factor structure facilitates structural identification of the VAR model korobilis2022new, chan2022large.
We also devise a new efficient and scalable Markov chain Monte Carlo (MCMC) algorithm. The key features of our algorithm are that we simulate the trees of the TVPs marginally and thus avoid mixing issues that arise if one is sampled conditionally on the other. The factor model assumed for the shocks not only ensures parsimony but also allows us to exploit computational gains by rendering the different equations of the model conditionally independent. Hence, we can use fast equation-by-equation updating.
To illustrate use of the model, we revisit the debate on the possibly evolving nature of the Phillips curve in the US. Given heteroskedasticity, we first identify the error volatility factors up to sign and scale chan2022large, and then isolate a business cycle shock in a narrative fashion, by assuming that it is the factor that explains the largest share of variation in the reduced-form VAR residuals to output and unemployment variations during recessionary periods. Using this identified shock, we investigate the dynamic reactions of a panel of macroeconomic quantities with a particular focus on prices. Considering how different effect modifiers impact the posterior distributions of the impulse responses reveals that the effects of business cycle shocks vary nonlinearly with uncertainty and according to whether the economy is in a recessionary regime.
The remainder of the paper is structured as follows. Section (ref) introduces the general econometric framework. This section includes information on the likelihood function and our nonparametric treatment of the TVPs. Section (ref) introduces our prior setup and discusses our posterior simulation algorithm, while Section (ref) illustrates our techniques using US data. The final section concludes.
In this section we develop our flexible econometric model. Our goal is to model the evolution of an $M \times 1$-vector of macroeconomic time series, which we denote by $\{\bm y_t\}_{t=1}^T$. The elements in $\bm y_t$ might feature structural breaks, changing cross-variable dependencies, and/or different persistence behavior over time. In addition, it could be that the number of time series $M$ is large. Our framework will be capable of simultaneously handling many time series that might feature complex dynamics.
We assume that $\bm y_t$ evolves according to a VAR model with drifting parameters. Our TVP-VAR with $P$ lags is given by:
Note that we have written the TVP-VAR as involving $\bm A_{p}$ for $p=1,..,P$, which are $M \times M$-matrices of constant coefficients, and $\bm B_{pt}$ for $p=1,..,P$, which are $M \times M$-matrices of TVPs capturing deviations from the constant part of the model. $ \bm \epsilon_t$ is an $M \times 1$-vector of Gaussian shocks with mean zero and time-varying variance-covariance matrix $\bm \Omega_t$.
There is substantial evidence that macroeconomic data are driven by a small set of fundamental shocks bai2007determining. We incorporate this insight into our model by assuming that $\bm \epsilon_t$ features a common factor structure:
where $\bm q_t$ is a $Q_q ( \ll M) \times 1$-vector of Gaussian-distributed factors with mean zero and diagonal time-varying variance-covariance matrix $\bm R_t = \text{diag}(r_{1t}, \dots, r_{Q_q t})$. $\bm \Gamma = (\bm \gamma_1, \dots, \bm\gamma_{Q_q})$ refers to a $M \times Q_q$-matrix of factor loadings, while $\bm \varepsilon_t = (\varepsilon_{1t}, \dots, \varepsilon_{Mt})'$ is an $M \times 1$-vector of Gaussian idiosyncratic shocks with mean zero and diagonal variance-covariance matrix $\bm \Sigma_t = \text{diag}(\sigma_{1t}^2, \dots, \sigma_{Mt}^2)$. One key observation is that, conditional on the factors $\bm q_t$, the shocks $\bm \varepsilon_t$ are independent and equation-by-equation estimation is possible. This factor structure has been used in other papers to facilitate estimation of large VARs kastner2020sparse, chan2021comparing, clark2021tail and, in addition, this structure can also facilitate identification of the factors as structural VAR disturbances korobilis2022new,chan2022large. Moreover, in contrast to VAR-based estimation using a Cholesky decomposition of the error covariances, another convenient feature of our model is that it is invariant to how the variables are ordered in $\bm y_t$ chan2022large.
It is worth stressing that the factor model, without additional restrictions, is not point identified. This is because the factors and the loadings enter the likelihood in product form and are, thus, not invariant to rotation, column and sign switching. Point identification can then be achieved by introducing restrictions on the loadings. A standard restriction assumes that the first $Q_q \times Q_q$ leading matrix of $\bm \Gamma$ is lower uni-triangular. This immediately implies that the resulting estimates will depend on the ordering of the elements in $\bm y_t$, a property that we would like to avoid. Hence, in what follows we do not impose identification restrictions on our factor model during MCMC estimation. Results in chan2022large suggest that the decomposition in (ref) is identified up to column and sign switching. Since our focus is on impulse responses to changes in particular elements in $\bm q_t$, we tackle column and sign switching ex-post. kaufmann2019bayesian follow a similar strategy, post-processing the posterior factor draws to point-identify the factors/shocks. We discuss this detail further in the empirical application below.
Up to this point we have remained silent on how the latent states (which include both the VAR coefficients and the time-varying elements of the error variances) evolve over time. In the next section, we introduce a flexible law of motion for the latent states.
A standard assumption in the TVP-VAR literature is that the elements in $\bm B_{pt}$ $(p = 1, \dots, P)$ and $\bm R_t$ evolve according to simple parametric stochastic processes, most often random walks primiceri2005, cogley2005drifts, belmonte2014hierarchical, bitto2019achieving. Assuming that the states evolve according to random walks introduces parsimony, because it implies a prior on the smoothness of the time variation in the coefficients. However, in turbulent periods, such as during the global financial crisis or the COVID-19 pandemic, it could be that parameters change rapidly and display sharp structural breaks. In such a case, a mixture model that models the evolution of the parameters as characterized by a low number of breaks sims2006were,koop2007estimation, kaufmann2015k would be more appropriate. Another possibility is that parameter change could depend on exogenous effect modifiers, translating into a specification with interaction effects. The nonparametric approach that we develop allows for all these possibilities.
Although most TVP-VARs allow for each coefficient to have its own random walk process, this is probably too flexible. That is, it is an empirical regularity that there is a high degree of co-movement in the parameters. This motivates the inclusion of a factor structure in the TVPs (that is, allowing for the the process innovation variance-covariance matrix to be of reduced-rank). In the parametric TVP-VAR literature, chan2020reducing propose a model that assumes a factor structure on the TVPs and assumes that the factors driving the states evolve according to a random walk. fischer2022TVP modify this approach by allowing for different forms of parameter change. This is achieved through including effect modifiers that can be either observed or latent.
In this paper, we do something similar, but we do it nonparametrically. That is, we assume there are a small number of latent nonlinear factors driving parameter change, which we estimate nonparametrically. In other words, we remain agnostic on the precise law of motion of the latent states, and let the data decide on the appropriate state dynamics, while achieving parsimony by introducing a factor structure to the TVPs.
We begin by writing the TVP-VAR in more compact form. Let $\bm x_t = (\bm y'_{t-1}, \dots, \bm y'_{t-p})'$ denote a $K = (MP) \times 1$-vector of covariates. Moreover, let $\bm A = (\bm A_{1}, \dots, \bm A_{P})$ and $\bm B_t = (\bm B_{1t}, \dots, \bm B_{Pt})$ refer to $M \times K$-matrices that stack the VAR coefficients. Since our model, conditional on the latent factors, is a system of independent regression models we can focus on the $m^{th}$ equation of $\tilde{\bm y}_t = \bm y_t - \bm A \bm x_t$. This regression model can be expressed as:
Here, $\bm \beta_{mt}$ and $\bm \gamma_m$ refer to the $m^{th}$ rows of $\bm B_t$ and $\bm \Gamma$, respectively. The factors $\bm q_t$ arise from a Gaussian distribution with variance $\bm R_t = R(\bm z_t) = \text{diag}(r_{1}(\bm z_t), \dots, r_{Q_q}(\bm z_t))$ where $r_s: \mathbb{R}^N \to \mathbb{R}^+$ is an unknown function. Notice that the error variances depend on a set of effect modifiers in $\bm z_t$.
We assume that $\bm \beta_{mt}$ evolves according to:
where $F_m(\bm z_t) = (f_{m1}(\bm z_t), \dots, f_{mQ_\beta}(\bm z_t))'$ denotes an unknown function with $Q_\beta$ components $f_{mq}: \mathbb{R}^N \to \mathbb{R}$, $\bm \Lambda_m$ is a $K \times Q_\beta$-matrix of factor loadings, and $Q_\beta$ is the number of latent factors that drive the TVPs. In addition, we assume that $\bm \eta_{mt} \sim \mathcal{N}\left(\bm 0_K, \bm V_m\right)$ is a vector of Gaussian shocks with $\bm V_m = \text{diag}(v_{m1}^2, \dots, v_{mK}^2)$ denoting the process innovation variances.
Conditional on choosing an appropriate number of factors $Q_\beta$, this specification is extremely flexible. It allows for (potentially) nonlinear interactions between $\bm z_t$ and $\bm \beta_{mt}$ (and thus implicitly $\bm x_t$). If elements in $\bm \beta_{mt}$ do not depend on $\bm z_t$ the corresponding loadings are zero and time-variation can still be captured through the presence of the idiosyncratic shocks in $\bm \eta_{mt}$.\footnote{Our model can also be related to random coefficient models, see fruhwirth2004bayesian.} Notice that the variances in $\bm V_m$ also control the weight put on the nonlinear factor component. For instance, if the $j^{th}$ coefficient closely co-moves with the other coefficients in a nonlinear manner, $v_{mj}^2$ will be close to zero.
To make this model operational we have to learn the functions $F$ and $R$ and decide on appropriate effect modifiers in $\bm z_t$. Our approach remains agnostic on the specific shape of both $F$ and $R$ and uses BART to estimate them. The next sub-sections show how this is achieved.
The choice of effect modifiers should depend on the application. The modifiers could include exogenous regressors, deterministic functions of time, lagged elements of $\bm y_t$, or latent quantities. If elements in $\bm z_t$ are endogenous and interest centers on higher-order impulse responses or multi-step-ahead predictive densities, one could either set up a separate law of motion for $\bm z_t$ or introduce hard restrictions on how the $\bm z_t's$ are expected to evolve over the forecast/impulse response horizon. In our empirical application below, we follow the latter approach, not only for simplicity, but because we are interested in how the dynamic reactions of $\bm y_t$ to shocks depend on the elements in $\bm z_t$ taking on certain values.\footnote{This would resemble common practice in, e.g., threshold or Markov switching models that condition on the prevailing regime when computing impulse responses.} This enables us to answer what-if questions, such as, “How would inflation react to business cycle shocks if uncertainty is (and remains) high?", or, “How do price reactions to business cycle movements change if the population becomes increasingly over-aged?"
Another interesting possibility would be to set $\bm z_t= \bm x_t$. In this case, however, interpretation becomes more difficult since the model then becomes nonlinear in $\bm y_t$. This would then necessitate the use of generalized impulse responses koop1996impulse to carry out dynamic analysis. With our application seeking to characterize features of the US business cycle, we choose effect modifiers that either slowly evolve independently of the business cycle, like factors related to the age-composition of the population, or binary variables (such as recession indicators), or other variables not included in $\bm y_t$. One such variable we consider is the uncertainty measure proposed in jurado2015measuring which can be interpreted as a proxy of (unobserved) macroeconomic uncertainty.
Another strategy to selecting the elements of $\bm z_t$ would be to entertain a large set of potential effect modifiers, and then use regularization techniques. As we will describe below, our approach is capable of handling all these cases without additional modification.
We approximate each function $f_{mq}$ through a sum-of-trees model chipman2010bart:
with the $T \times N$-matrix $\bm Z$ having a typical $t^{th}$ row $\bm z'_t$ and $u$ being a regression tree function that depends on a tree structure, $\mathcal{T}^{\beta}_{mq} = \{\mathcal{T}^\beta_{mq,1}, \dots, \mathcal{T}^\beta_{mq, S_\beta}\}$, which is a sequence of disjoint sets that partition the input space and a vector of terminal node parameters $\bm \mu_{mq} = \{\bm \mu_{mq,1}, \dots, \bm \mu_{mq,S_\beta}\}$ of dimension $b_{mq}$. These partitions are driven by splitting rules of the form $z_{jt} \le c_j$ or $z_{jt} > c_j$, with $z_{jt}$ denoting the $j^{th}$ element of $\bm z_t$ and $c_j$ being a threshold parameter. Moreover, $S_\beta$ is the number of trees used to approximate each of the functions (factors) $f_{mq}$. (ref) is a standard BART model. To avoid issues associated with overfitting when $S_\beta$ is large, chipman2010bart propose using a regularization prior to force the trees to take a particularly simple form and thus explain only a small fraction of the variation of the response variable. Adding together many simple trees (weak learners) has been found to work better than working with a single more complicated tree. We follow such an approach in this paper.
Plugging ((ref)) into ((ref)) yields our state equation:
with $\bm \lambda_{mq}$ denoting the $q^{th}$ column of $\bm \Lambda_m$. This shows that we combine $Q_\beta$ BART models (each used to approximate one of the $Q_\beta$ functions). Setting $Q_\beta=1$ implies that all coefficients are driven by a single factor (if $\bm \lambda_{mq} \neq \bm 0_K$), while when setting $Q_\beta=K$ we obtain a model closely related to the one proposed in deshpande2020vcbart and coulombe2020macroeconomy. Since the latter specification, in light of large $K$, does not scale well to high dimensions we will focus on the case where $K \gg Q_\beta$, which frequently arises in the analysis of large TVP-VAR models. We will call models that assume this nonparametric factor form for the conditional mean, TVP-FBART. Ahead of our main empirical application and to help the reader further understand our model, Sub-section (ref) in the Online Appendix provides a toy empirical example to illustrate how BART can be used to approximate TVPs.
Recall that our model also assumes that the shocks $\bm \epsilon_t$ feature a factor structure. We will again approximate the factor-specific functions, $r_s(\bm z_t)$, in $\bm R_t$ with BART. More precisely, our approach can be interpreted as a variant of heteroskedastic BART pratola2020heteroscedastic. heteroBART is a multiplicative version of BART and assumes that the trees enter the model in product form. In this paper, we follow clark2021tail and linearize the model so that standard BART techniques can be used.
Let the $s^{th}$ element of $\bm q_t$ be given by:
with $S_q$ being the number of trees used to approximate the variance functions, where $\mathcal{T}^{q}_{sd}$ and $\bm \pi_{sd}$ denote the corresponding tree structures and terminal node parameters, respectively. To render (ref) linear we square it and take logs. This yields a linear equation with shocks that are log-$\chi^2$ distributed with one degree of freedom, a distribution which can be well approximated using a ten-component mixture approximation omori2007stochastic:
Here, $m_i$, $w_i$, and $\mathfrak{r}_i^2$ are fixed numbers defining the mixture components taken from Table 1 in omori2007stochastic.
This specification of heteroBART implies that the factor volatilities are allowed to change rapidly, but can also move more gradually. This feature might pay off during recessions, where large jumps in error volatilities are common. Traditional stochastic volatility models will be unable to match this pattern, since they assume that the log-volatilities evolve according to a stochastic process that translates into a more gradual evolution of the error variances. Such behavior is warranted if the trend movement in volatility is persistent (such as during the Great Moderation). For heteroBART, matching slowly evolving trends is also possible but considerably harder. To allow for smoothly evolving stochastic trends we combine heteroBART with a standard stochastic volatility model in the measurement errors (that is, the elements in $\bm \Sigma_t$). The combination between a parametric law of motion for $\log(\sigma^2_{mt})$ and $q_{st}$ allows for rich dynamics in terms of $\bm \Omega_t$.\footnote{Another option to capture smoothly varying trends with heteroBART would be through the specification of a latent component which enters $\bm z_t$. But this would require nonlinear filtering algorithms or linear approximations (which can fail in certain environments) such as the ones proposed in huber2020nowcasting.} We will use the abbreviation FHB (using a factor structure involving heteroBART) for models which adopt this specification. Thus, our most general model is TVP-FBART-FHB.
The model described in the previous sub-sections is very flexible and nests a wide variety of competing models. In this sub-section, we first summarize key model features and then discuss how our model is related to alternative models commonly used in the literature.
Flexible machine learning techniques such as BART have the shortcoming that interpretability is difficult. As noted, for example in coulombe2020macroeconomy, using regression trees to model the parameters of a TVP regression allows for flexibility but also maintains simplicity of interpretation. In our case, once we have learned the TVPs and the functions driving them using BART, interpretation of the model works analogously to a standard TVP-VAR model. Hence, one can compute functions of the parameters such as impulse responses, forecast error variance or historical decompositions, and conditional forecasts using standard techniques. This constitutes a big advantage of our approach relative to models such as the one proposed in huber2022inference. Computation of (generalized) impulse response function in traditional BART-based VAR models is much more involved, as the model remains nonlinear.\footnote{koop1996impulse discuss how to compute generalized impulse response functions in nonlinear multivariate models.}
The previous paragraph is related to the effect that shocks might have on $\bm y_t$. Since the effect modifiers influence $\bm y_t$ indirectly through the BART modeling of the TVPs, we can also assess how $\bm z_t$ affects $\bm y_t$. This can be easily achieved in our framework since one can compute different realizations of the TVPs for different configurations of $\bm z_t$. Doing so allows us to study how (higher-order) interaction effects, which might take an unknown form, impact quantities such as forecast distributions, impulse responses, or even long-run trends such as the (time-varying) unconditional mean of the TVP-VAR. We will illustrate these features in our empirical work that follows in Section (ref) below.
Apart from the ease of interpretation and the additional inferential possibilities, our model, for appropriately chosen values of $Q_\beta, Q_q, S_\beta$, and $S_q$, provides a great deal of flexibility when it comes to capturing different forms of parameter change. While our aim is to introduce as few restrictions on the state evolution as possible, we can nevertheless control the dynamics of the TVPs by choosing appropriate values of $S_\beta$ and $S_q$. In principle, larger values of $S_\beta$ and $S_q$ are consistent with smooth law of motions of the parameters, whereas smaller values imply parameter dynamics closer to the ones generated by a structural break model. An extreme case of our model would set $Q_\beta = S_\beta = 1$. This specification would imply that parameters follow a single regression tree and are proportional to each other. In our empirical work we will explore the sensitivity of results by varying these parameters.
We start our discussion with the priors relating to the regression trees and the process innovation variances. The remaining priors are relatively standard and a discussion can be found in Section (ref) in the Online Appendix.
chipman1998bayesian and chipman2010bart specify a tree-generating stochastic process on the tree structures, $\mathcal{T}^{\beta}_{mq, s}$ and $\mathcal{T}^{q}_{sd}$. Our approach is similar, but specifies the prior such that the probability of growing more complex trees decreases with the number of factors $\nu =1,\dots, Q_j$ for $j \in \{\beta, q\}$. This process is designed to penalize complex trees and consists of three features:
This prior encourages smaller trees and is thus consistent with the notion that each individual tree is a “weak learner," but the composite model is capable of capturing complex dynamics in the parameters.
The prior on the terminal node parameters is Gaussian. Following chipman2010bart, we scale the data such that the dependent variable is between $-0.5$ and $0.5$ and our prior covers this range. Let $\mu_{mq, ij}$ denote the $j^{th}$ element of $\bm \mu_{mq, i}$ and $\pi_{sd,j}$ the $j^{th}$ element of $\bm \pi_{sd}$. The prior for the respective element is then given by:
Here, $\kappa$ is a parameter that controls the prior variance. Shrinkage is introduced by increasingly forcing $\mu_{mq, ij}$ ($\pi_{sd, j}$) towards zero if $S_\beta$ ($S_q$) is large. Since $S_j$ ($j \in \{\beta, q\}$) is typically between 50 and 200, this prior is the second ingredient of BART used to capture the notion that each tree explains only a small amount of variation in $\bm \beta_t$ (and $q_{st}$).
On the different elements of $\bm V_m$, several priors are possible. The simple conjugate inverse Gamma prior can be used. This prior, however, has implications for our model, since it rules out values of $v^2_{mj}$ very close to zero. Hence, it would artificially push the likelihood away from the factor part in (ref). We follow recommendations in fs_wagner and use a prior that introduces shrinkage on $\bm V_m$. Our prior assumes that $v^2_{mj}$ arises from a Gamma distribution:
with $B_v$ being a scalar hyperparameter that controls the amount of shrinkage towards a factor structure in the TVPs. Since there exists strong evidence that the TVPs feature a factor structure, we set $B_v=0.01$ to have a tight prior on the idiosyncratic deviations of the TVPs from the common factor structure.
We sample from the joint posterior distribution of the model by using an MCMC algorithm that, conditional on the latent factors, simulates the coefficients and latent states for each equation separately. Since for some of the steps in the sampler we integrate out other parameters the precise ordering of the steps of the MCMC algorithm is important to simulate from the correct stationary distribution. Our algorithm cycles between the following steps. For each equation $m=1,\dots, M$:
The following quantities are not estimated in an equation-by-equation manner:
Notice that steps (1) to (3) yield a draw from $p(\{\mathcal{T}^\beta_{mq}, \bm \mu_{mq}, \bm \lambda_{mq}\}_{q=1}^{Q_\beta}|\bm \bullet_{/\bm \beta_{mt}})$ where the notation $\bullet_{/\bm \beta_{mt}}$ indicates the remaining model parameters except the TVPs and the data. The TVPs are then simulated from $p(\{\bm \beta_{mt}\}_{t=1}^T|\bullet)$ where $\bullet$ means all other model parameters, latent quantities and the data. This step differs from the one used in deshpande2020vcbart since we improve mixing by integrating out the TVPs. In principle, the loadings and trees can also be sampled conditionally on the TVPs but in cases where the loadings are very small substantial mixing issues arise.
We repeat this algorithm $15,000$ times and discard the first $5,000$ draws as burn-in.\footnote{To obtain $15,000$ draws, the actual computation time is about $124$ minutes, based on a MacBook Pro with an M1 8-core processor.} From a computational perspective, this algorithm is quite efficient. This is because the sampling step associated with the TVPs can be sped up enormously by exploiting the fact that $\bm W_m' \bm W_m$ is a block-diagonal matrix of rank $T$.
We use the quarterly version of the mccracken2016fred data set and focus on a sample ranging from $1975$:Q$1$ to $2019$:Q$4$. In our empirical work, we aim to investigate how business cycle shocks impact a range of different price measures and whether these dynamic reactions depend on the effect modifiers. To this end, we follow del2020s and estimate medium-sized VAR models that are rich in wage, price, and labor market measures. We consider $M=12$ endogenous variables, where $\bm y_t$ includes output growth, employment, unemployment, average weekly hours worked, personal consumption expenditure (PCE) inflation, PCE inflation excluding food and energy, (core) consumer price inflation, the GDP deflator, wage inflation, the federal funds rate, and ten-year government bond yields to capture movements in treasury markets. But unlike del2020s, we allow for nonlinear relationships between these variables and for these effects to vary over time. del2020s accommodate temporal change only, by simply estimating their linear VAR model over two non-overlapping samples.
As effect modifiers in $\bm z_t$, we consider five indicators that may affect the TVPs, and in turn the impulse response functions, in a nonlinear manner. Specifically, we consider the old-age dependency ratio, a financial globalization indicator, the (lagged) ex-post real rate, a binary recession indicator (taken from the NBER), and the economic uncertainty index proposed in jurado2015measuring. Secular stagnation factors, such as a boost in financial globalization, the rising old-age dependency ratio, and a declining real rate, may affect the dynamics of business cycle phases in a nonlinear manner jones2022aging. These factors have also been identified as one cause of the flattening of the Phillips curve forbes2019inflation, forbes2021low. The last two effect modifiers allow for possible structural breaks in recessionary and high uncertainty periods aastveit2017economic, alessandri2019financial.
Some of the effect modifiers are clearly endogenous and should depend on the other quantities of our model. This does not cause any issues for the validity of our econometric approach. However, when we focus on impulse responses it has the implication that $\bm z_t$ is not allowed to react to changes in $\bm y_t$. As discussed in Section (ref), this is an assumption made for the sake of interpretability. The main implication is that impulse responses can be understood as being conditional on $\bm z_t$ remaining at the current level over the impulse response horizon. Since we are going to construct “scenarios,” based on assumptions about how $\bm z_t$ behaves, this restriction can be interpreted as similar in nature to conditional forecasts when the restricted variables are not located in $\bm y_t$ but in $\bm z_t$. If the researcher wishes to relax these assumptions, they can set up auxiliary models for $\bm z_t$, such that $\bm z_t$ is again a function of $\bm y_t$.
Table (ref) in the Online Appendix provides additional information on the time series and associated data transformations used. All models we consider in this paper feature $p=5$ lags. In Sub-section (ref) we assess how different model features impact model fit and compare our proposed model to standard models in the literature. This analysis evidences that our model generally captures the data well, often improving upon competitors commonly used in the literature. Based on the results in Table (ref), we use the model that sets $Q_\beta=25$, $S_\beta=1$, $Q_q=3$, and $S_q=250$.
In this sub-section we consider what is driving the time variation in the VAR coefficients in our TVP-BART model with FHB. (ref)(a) shows a heatmap of the total share of time-variation of $\bm \beta_{mt}$ explained by the nonlinear factors $F_m(\bm z_t)$ across equations $m=1, \dots, M$. This quantity, closely related to the familiar $R^2$, is computed as follows:
with $\text{Var}(g(\bm z_t|\mathcal{T}^\beta_{mq}, \bm \mu_{mq})$ denoting the empirical variance of the function $g$. Dark red values indicate that a given TVP is driven almost exclusively by $F_m(\bm z_t)$, whereas white values suggest that most of the variation is driven by idiosyncratic movements in the TVPs.
Panel (b) of (ref) displays a heatmap of posterior means of the number of tree splits induced by one of the effect modifiers in $\bm z_t$ across coefficients and equations. This serves as a way to assess the relative importance of different effect modifiers in shaping the coefficient dynamics over time.
Starting with panel (a) of the figure, we see that the explanatory power of the TVP factors varies substantially across equations (and also across variables). While we find that TVPs in the interest rate and CPI core equations are strongly shaped by the effect modifiers, this share is considerably lower for the other equations. With two exceptions (PCETCPI and GDPCTPI), the shares are, however, sizable and often above 50 percent. Turning to PCETCPI and GDPCTPI, the effect modifiers explain a rather small amount of variation. Interestingly, for labor market quantities (EMPL, UNRATE, AWH) and real GDP we also find that the intercept (which determines the unconditional mean of the model) is strongly influenced by different effect modifiers. This indicates that long-run properties of these time series depend on covariates that may be interpreted as capturing structural change in the macroeconomy.
Focusing on panel (b) of the figure provides additional insights. First, the old-age dependency ratio, financial globalization, and the real rate play only a limited role in explaining parameter dynamics. Second, for several variables we find that uncertainty shapes TVP dynamics. Among these are coefficients in the CPI and CPI core equations, the short-term interest rate equation, and the ten-year government bond yield. Third, for other variables such as output, employment, and the unemployment rate, we observe that uncertainty plays a more limited role. However, in these equations we instead find that the NBER's recession indicator is frequently included in the splitting rules.
One of the main advantages of our nonparametric model is that, conditional on knowing the TVPs and error covariances, the model is a standard linear TVP-VAR model. Hence structural analysis, using identified impulse responses, can be readily carried out. In principle, an economist's preferred identification strategy based on, for example, sign restrictions benati2008, zero impact restrictions primiceri2005, koop2009evolution, or long-run restrictions can be implemented within our TVP-BART framework.
In this application, we focus on the question of how adverse business cycle shocks impact a set of inflation measures. To do so, we exploit the factor structure on the reduced-form VAR shocks to identify a business cycle shock korobilis2022new, chan2022large. As emphasized by Gorodnichenko2005, in VAR models like ours where the number of variables is relatively large (we have $M=12$) it can facilitate structural interpretation to have fewer structural shocks than $M$. In the next step, we trace out the dynamic evolution of our inflation measures to such a business cycle shock.
One can decompose, as in (ref), the reduced-form VAR shocks into a factor component and an idiosyncratic measurement-error component (both of which are independent) under standard conditions anderson1956statistical, fruhwirth2018sparse, kaufmann2019bayesian.\footnote{These conditions relate to the number of factors being smaller then the Ledermann bound and the number of non-zero elements in $\bm \Gamma$ being sufficiently large so that the decomposition in (ref) is unique.} Absent heteroskedasticity, the resulting factors still have no economic interpretation and thus additional structure is required to identify the shocks, given that the factors and the factor loadings can be rotated by any random orthogonal matrix. But, given the heteroskedasticity in $\bm R_t$, we can follow chan2022large and identify, up to sign and scale, a business cycle shock as that factor (shock) that explains the largest amount of variation in innovations to output and unemployment variations during recessionary periods (as identified by the NBER). This identification strategy resembles the one proposed in bianchi2023inflation. They identify business cycle shocks by searching for linear combinations of the reduced-form shocks of a trend-cycle VAR so as to maximize the amount of variation in unemployment or cyclical output.\footnote{Alternative approaches to identify business cycle shocks are proposed in del2020s and Angeletos.}
Specifically, our business cycle shock is obtained by computing:
for all $j$ and finding that factor that maximizes the variances explained for real GDP and the unemployment rate during recessionary episodes. This yields, for each MCMC draw, a factor that can be interpreted as a business cycle shock. To point-identify the sign of the factors and the associated loadings, we normalize the factors and loadings to identify the business cycle shock as having a negative impact effect on output growth and a positive impact effect on unemployment.
In summary, we identify the business cycle shock and the associated impulse responses via the following steps:
These steps yield partial identification, implying that the business cycle shock is uniquely identified whereas the remaining $(Q_q - 1)$ factors (and the associated columns in $\bm \Gamma$) are left unrestricted. This identification approach is related to ones developed in recent papers korobilis2022new, chan2022large which advocate using sign restrictions on the factor loadings to pin down a shock of interest. But our approach differs in the sense that we solve the column switching problem (which is required to attach an economic meaning to the different factors) through a narrative approach that builds on the notion that business cycle shocks are the ones that determine the largest amount of variation in real activity quantities during recessions. Our approach could easily be combined with sign-restricted factor stochastic volatility models, by introducing certain restrictions on the prior associated with $\bm \Gamma$.
(ref) plots the posterior mean of the proportion of the variation, $\bm \zeta_{jt}$, in the unemployment rate and in output growth explained by the three factors over time. This figure shows that the second factor explains the largest amount of variation in the early part of the sample (until the twin recession of the early 1980s) and during all recessions in our sample. During recessions, this factor explains close to 70 percent of the variation in the reduced-form shocks to both unemployment and output growth.
In this sub-section, we look at the dynamic effects of our business cycle shock. Impulse response functions are computed by shocking the business cycle (second) factor and tracing out the dynamic reactions of $y_{t+h}$ for $h=1, \dots, 16$.
Since our model features TVPs, the impulse responses can be computed at each point in time. This gives us a posterior distribution over $T$ period-specific IRFs, a statistical object that is difficult to visualize. To aid exposition, we start our analysis by considering average impulse responses. These are obtained by averaging the time-specific impulse responses over time and are depicted in (ref).
(ref) shows that, averaged over time, a contractionary business cycle shock leads to unemployment rising and inflation (including core and wage inflation) falling. The dynamic effects on the different inflation measures are similar, but long-lasting. Like the main business shock of Angeletos, the peak effect of our business cycle shock on the real variables also occurs within a year or two. Specifically, we observe that output and employment decline while the unemployment rate increases. Real GDP growth reacts rapidly by declining by around one percentage point on impact. For employment growth, the peak effect materializes after about three quarters. The unemployment rate quickly increases and displays a peak reaction of around 0.5 percentage points after around one year. These reactions are largely consistent (both in terms of shape and size) with the ones reported in bianchi2023inflation.
When we consider the reactions of our different inflation measures, we find that prices decline on impact. This reduction appears to be quite persistent. As we will show below (see (ref)), this persistent reaction of prices is mainly driven by strong and persistent declines of inflation up to the early 1990s. These results are consistent with the existence of a negatively sloped Phillips curve --- at least on average through the $1975$:Q$1$ to $2019$:Q$4$ period.
To hone in on this relationship between inflation and unemployment, we normalize the IRFs of the different price measures by the IRFs of the unemployment rate. barnichon2021phillips call this quantity the Phillips curve multiplier. The multipliers are shown in (ref). We again see that, on average over time, the Phillips curve multipliers are negative and statistically significant. (ref) also reveals that these negative effects persist for two to three years. And they vary by inflation measure. As we should expect, the Phillips curve is stronger for headline than core measures of inflation. The strongest effects on inflation are typically seen two years after the business cycle shock.
To understand to what degree averaging over time is masking temporal variations in the Phillips curve relationship, (ref) plots, at the one-year-ahead horizon ($h=4$), the impulse responses due to the contractionary business cycle shock at each point in time. The ability to identify and capture structural change of different forms is a key feature of our model. (ref) reveals that there are indeed important temporal variations. The responses of, in particular, the headline inflation measures become more muted over time. Focusing in on the effects on CPI inflation, we see that the business cycle shock lowers inflation significantly through the $1970$s and $1980$s. But the responses thereafter are more muted. They become increasingly muted as we look to the period after the global financial crisis. Interestingly, evidencing a clear nonlinearity, there is a strong negative effect on inflation during the recessionary period associated with the global financial crisis itself. Our findings therefore provide ex-post justification for the decision by del2020s to estimate their VAR model, designed to understand the Phillips curve, on samples before and after $1990$. But our results also reveal important temporal instabilities and changes within these two periods that are lost by simple sample-slit or indeed rolling regressions as also often used in the literature.
Turning to the effects on unemployment, again consistent with del2020s, (ref) shows that the response of unemployment to a business cycle shock becomes more persistent over time. This is consistent with economic expansions lasting longer in more recent decades. But (ref) adds texture to this narrative by revealing that recessionary periods, except for $2008$-$9$, are marked by especially strong responses.
Bringing together the price inflation and unemployment responses, we conclude that the sensitivity of price inflation to unemployment has weakened markedly since $1990$. The response of wage inflation to the the business cycle shock is weaker throughout the sample. This casts doubt on the view Knotek, HOOPER2020 that the Phillips curve is stronger for wage than for price inflation. Since the $1990$s the impulse responses for wage inflation and price (CPI) inflation look broadly similar; see (ref). This includes evidence that wage as well as price inflation did decline in response to a business cycle shock during the Great Recession, with prices declining by more than wages.
One key feature of our model is that it allows us to link the time variation in the parameters (and thus functions thereof such as IRFs) to the effect modifiers in $\bm z_t$. Since $\bm z_t$ influences the TVPs using a nonparametric model, it is difficult to clearly answer how changes in $\bm z_t$ impact the TVPs. Since our interest centers on the implied IRFs, we can, however, carry out simulations that show how the dynamic responses to a business cycle shock change as we vary $\bm z_t$. The results of this exercise are shown in (ref). This figure depicts the price responses in the rows of the panel and in the columns shows different assumptions on $\bm z_t$.
To analyze whether IRFs differ in expansion and recessions, we set the NBER recession indicator to zero (that is, we assume that the economy is in an expansion) or to one (that is, we assume that the economy is in a recession). Based on this, we vary one of the effect modifiers while setting the remaining effect modifiers to some pre-specified value. This pre-specified value is either the average value over the period $1975$ to $1985$ (which are the blue-shaded IRFs in the figure) or the period after $2010$ (which are the red-shaded IRFs). This allows us to capture the general macroeconomic environment in the respective time periods. This gives us four overall combinations for the IRFs. We consider how prices react in expansions and recessions and whether there are discernible differences in the transmission of business cycle shocks in these two regimes. Based on one of these four general scenarios, we set each effect modifier (for example, the dependency ratio, financial globalization, the real rate, and uncertainty) equal to different sample quantiles and then compute the implied price IRFs. This provides a detailed picture on how impulse responses depend on the effect modifiers.
(ref) confirms that the headline (non-core) inflation measures were more strongly affected by business cycle shocks before 1985. The most striking nonlinearity for the effect modifiers is seen with respect to uncertainty. As uncertainty increases beyond its $75^{th}$ percentile, we see much stronger negative effects on all the inflation measures, including with post-2010 data. This effect is especially pronounced during recessionary periods. This all supports a view that the Phillips curve remains alive and well during times of recession and greater-than-average uncertainty, events that empirically tend to co-exist. This is consistent with theories of the financial accelerator, suggesting that shocks have amplified effects in recessions.
In this paper, we have developed a nonparametric model that uses Bayesian additive regression trees (BART) methods to allow for change of an unknown form in both the conditional means and variances of a multivariate time series model (a VAR). Unlike existing nonparametric approaches, interpretation and macroeconomic inference including structural analysis is easier, since, as the model gives a nonparametric treatment to the parameters rather than the variables, it remains conditionally linear in the mean. An additional novel feature is that the new model allows for nonparametric factor structures for parameters in the conditional means and variances, thus reducing the number of nonparametric functions to estimate and ensuring parsimony.
In an empirical exercise we show how the proposed nonparametric VAR model contributes to our understanding of the time-varying nature of the Phillips curve. Inflation has become considerably less sensitive to business cycle shocks, in particular since $1990$. However, the flexible nonlinear features of the model show that the effects on inflation remain strong when uncertainty rises to high levels.
{\setstretch{0.9} \addcontentsline{toc}{section}{References} }
\setcounter{page}{1} \setcounter{footnote}{0}