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.
74,009 characters · 16 sections · 39 citation commands
Bayesian Nonlinear Regression using Sums of Simple Functions
\doublespacing
Nonlinear modeling of large datasets has received increasing attention in recent years. Extreme events such as the Covid-19 pandemic and the surge in inflation in the aftermath of the pandemic have raised the interest in more flexible econometric models Coulombe2020randomforest, Coulombe2021covid, carriero2022addressing, hauzenberger2022fast, clark2023tail, huber2023nowcasting,koop2023bayesian.
Capturing nonlinearities in economic time series is predominantly achieved through estimating models with particular assumptions on the form of nonlinearities. For instance, Markov switching or structural break regressions and vector autoregressions (VARs) assume that the parameters in the conditional mean change abruptly and there are only few but large breaks sims2006were, koop2007estimation, bauwens2015contribution. By contrast, time-varying parameter (TVP) models primiceri2005time, COGLEY2005262, koop2013large, bitto2019achieving assume that the parameters evolve smoothly over time and thus feature a large number of small breaks in the regression coefficients.
All these methods have in common that they postulate a linear relationship between the endogenous variables and the regressors at particular points in time. By contrast, nonlinear regression assumes a nonlinear relationship between the endogenous variables and the predictors. This relationship remains constant over time. Some examples are white1984nonlinear, hamilton2001parametric, lubrano2001smooth, hamilton2003oil, gerlach2008bayesian, gefang2009nonlinear, bruns2023tractable. However, assuming a particular form of nonlinearity might give rise to model mis-specification and can be interpreted as a dogmatic Bayesian prior on the space of (nonlinear) conditional mean functions.
Another strand of the literature does not take a strong stance on the precise form of the conditional mean and uses nonparametric techniques to infer probable functional forms or detect structural breaks in the conditional mean. These methods remain agnostic on nonlinearities in the conditional mean and variances and try to infer them from the data. In recent years, nonparametric techniques have been increasingly used to forecast macro and financial aggregates clark2023tail, huber2023nowcasting, estimate nonlinearities in key macroeconomic relations such as the Phillips curve Coulombe2020randomforest, flexibly combine forecasts bassetti2018bayesian, to construct shrinkage priors for vector autoregressions billio2019bayesian and for pooling coefficients casarin2023bayesian. The key shortcomings of these methods is that they are difficult to implement, customize and to tune.
These techniques all have their own pros and cons. However, what they share is the lack of scalability to very high dimensions. While there has been much progress in recent years chan2023large the largest nonlinear models often feature less than 20 endogenous variables. These 20 indicators often represent only a small fraction of the series available in different macroeconomic databases provided by major central banks such as the US Federal Reserve or the Bank of England. For forecasting and structural analysis, exploiting as much information as possible can be important, increasing the demand for flexible models that can handle large datasets.
The last two paragraphs provide the main motivation for the current paper. We wish to develop techniques that are relatively simple to implement, modify and have the ability to handle large datasets commonly used in macroeconomics. These characteristics, however, should not come at the cost of reduced flexibility. We achieve this through a new parametric Bayesian nonlinear regression model that can be applied to univariate and multivariate time series and is inherently related to popular methods such as Bayesian additive regression trees chipman2010bart and shallow neural networks. Our main assumption is that the conditional mean is modeled through a sum of simple functions. These functions are two-component location mixtures with transition between regimes driven by a logistic transition function. The logistic function is parameterized by a speed of adjustment coefficient, a threshold variable and a threshold parameter. These are all estimated through Bayesian techniques. When viewed individually, each of these simple models explains only a small fraction of the variation in the response (i.e., it acts as a 'weak learner'). However, when we sum over a moderate to large number of logistic functions we obtain a great deal of representation flexibility and end up with a model that is straightforward to estimate and to implement.
The logistic function, while being tightly parameterized, is also flexible. For instance, if the speed of adjustment parameter becomes large, the transition function reduces to the indicator function that equals one if the threshold variable exceeds a threshold. If this applies for each of the individual functions we end up with an extreme version of BART with very simple trees.
Computation is carried out under conjugate priors. These provide further regularization but, more importantly, give rise to substantial computational gains. In particular, the algorithms we develop are highly scalable and can handle systems with hundreds of endogenous variables, leading to a huge dimensional nonlinear VAR model.
We start by illustrating our techniques by means of simulated data. Using a highly nonlinear DGP, we show that our parametric Bayesian model produces point and density forecasts that are often better than the ones produced by BART. We find that, as opposed to BART, the optimal number of functions to sum over is between 5 and 15 and thus much smaller. Moreover, we also find that fixing the speed of adjustment parameter so that the transition between regimes is instantaneous yields results that are only slightly worse than the ones from the model that estimates all parameters of the transition function.
We then move on to the real data analysis and estimate a large nonlinear VAR of the US economy. This analysis consists of two parts. In the first, we show that our approach yields highly competitive density forecasts relative to the BART-VAR of clark2023tail. In the second, we illustrate how our model can be used to analyze the nonlinear effects of financial shocks on the US economy. This exercise shows that for a dataset comprising of 80 endogenous variables, substantial asymmetries arise between benign and adverse shock. But this only holds if the shock is sufficiently large.
Our plan for the remainder of the paper is the following. We will introduce our main techniques in the next section. In this section, our focus is on approximating a nonlinear univariate regression model using sums of simple logistic functions. We provide an illustrating example, derive the likelihood, specify conjugate Bayesian priors for the parameters of the model and discuss posterior simulation. The next section, (ref), provides simulation evidence for this model. Then, in (ref), we generalize the model to the multivariate case. (ref) applies the model to a large US dataset and includes a forecasting exercise and the structural application. The final section summarizes and concludes the paper.
We start our discussion by focusing on the univariate case. Suppose that we have a time series $\{ y_t \}_{t=1}^T$ and model it as a nonlinear function of a large panel of $K$ predictors $\bm x_t = (x_{1,t}, \dots, x_{K,t})' \in \mathbb{R}^K$. We approximate this nonlinear function using a sum of $J$ simpler functions (also called base learners):
where $g: \mathbb{R} \to \mathbb{R}$ is a simple function that is fully parameterized by a low-dimensional vector $\bm \theta_j$. We will assume that $g$ is given by:
with $\beta_{i,j}~(i=0, 1)$ denoting a switching intercept term and $S_{j,t} \in [0, 1]$ is a transition function. We let $\Tilde{x}_{j,t} = \bm \delta'_j \bm x_t$ denote an element of $\bm x_t$ and $\bm \delta_j$ is a $K-$dimensional selection vector. If the $s^{th}$ element of $\bm \delta_j$ equals $1$, the $s^{th}$ variable in $\bm x_t$ is selected and hence $\Tilde{x}_{j,t} = x_{s,t}$. In what follows, we suppress the dependence of $S_{j, t}$ on $\Tilde{x}_{j,t}$.
It is worth stressing that, as opposed to other algorithms, we only assume that a single variable informs the transition between two regimes. Using the jargon of the boosting literature schapire2003boosting, the function in (ref) acts as a weak learner bai2009boosting and is expected to explain only a small amount of the variation in $y_t$. However, summing over multiple functions will provide sufficient representation flexibility to approximate any conditional mean function. This finding builds on theoretical results in cybenko1989approximation and is closely related to the universal approximation theorem in the literature on machine learning.
For Theorem 1 in cybenko1989approximation to work we need to make a few additional assumptions on the transition function $S_{j,t}$. In particular, we need to assume that $S_{j,t} = 0$ if $(\Tilde{x}_{j,t} \to -\infty$ and $S_{j,t} = 1$ if $(\Tilde{x}_{j,t} \to \infty$. A general function that fulfills this is the logistic function:
whereby $\nu_j \in \mathbb{R}^+$ is a speed of adjustment parameter and $\mu_j \in \mathbb{R}$ is a threshold parameter. The parameter $\nu_j$ controls the smoothness of the transition function. If it equals $0$, $S_{j,t}$ equals 1/2 and $\Tilde{x}_{j,t}$ does not enter $S_{j,t}$. If it is greater than zero but not too large we have a smooth transition between regimes with the transition being driven by the movements in $\Tilde{x}_{j,t}$. In this case, we would end up observing an S-shaped function. By contrast, if $\nu_j$ becomes large, we end up with an indicator function that equals zero if $\Tilde{x}_{j,t} > \mu_j$ and one otherwise. We call this model additive smooth transition (AST) model and, for later convenience, we let $\bm \theta_j = (\mu_j, \nu_j, \bm \delta_j, \beta_{0, j}, \beta_{1, j})'$ denote the vector of component-specific parameters.
The main advantage of (ref) is that if $\Tilde{x}_{j,t}$ exerts a smooth effect (implying a gradual transition between regimes), the logistic function captures this through estimates of $\nu_j$ closer to zero. By contrast, if $\Tilde{x}_{j,t}$ might only have a threshold effect, the model would estimate $\nu_j$ to be large and thus lead to a heavy side function. By summing over many of these functions and allowing for the different parameters (thresholds and speed of adjustment coefficients) to vary our model provides a great deal of flexibility.
Our model is related to, at least, two popular models in the literature: BART and neural networks (NNs). BART is obtained if $g$ is replaced with a considerably more complex tree function. In this case, the dimension of the parameter vector $\bm \theta_j$ is not known a priori, rendering the model nonparametric. In many applications in a vast range of different fields, BART has been among the best performing specifications in terms of achieving low out-of-sample forecast errors chipman2010bart. However, as opposed to our approach, if one wishes to apply BART to customized models (such as VARs) substantial coding efforts are required and while estimation of larger models is possible,\footnote{The largest BART model clark2023tail consider features around 20 endogenous variables.} scalability to large simultaneous equation models such as the one we consider in our applied work, is currently unfeasible.
Another model closely related to the one presented in this section is the (shallow) NN. A shallow NN sets $\bm \delta_j \in \mathbb{R}^K$ equal to a weight vector. By doing so, every element in $\bm x_t$ informs the corresponding component-specific function. In addition, the transition functions often take different forms and enter the conditional mean equation as transformed regressors with separate coefficients. The key disadvantage relative to our approach is that it requires estimating a (possibly huge dimensional) coefficient vector per component function $J$. If $J$ becomes large, this becomes computational intensive and fully Bayesian inference is difficult to carry out in large models.
Our model is best understood by considering a simple illustrative example where we fix $\nu_j$ and $\mu_j$. In this case, the intercept parameters $\beta_{0,j}$ and $\beta_{1,j}$ can be obtained through OLS. We consider the quarterly growth rate of US industrial production (IP) from $1990$:Q$1$ to $2019$:Q$4$. Our goal is to model IP growth as an unknown function of lagged IP growth and the excess bond premium (EBP) of gilchrist2012credit.
Consider the case of $J=1$ first and, let us assume, that $\Tilde{x}_{1,t}$ is the EBP, the threshold $\mu_1$ is the mean of the EBP and the speed of adjustment parameter is $\nu_1=0.3$, implying a smooth transition between regimes. The parameters $\beta_{0,1}$ and $\beta_{1,1}$ are estimated to be $-0.9$ and $1.4$, respectively.
The resulting transition function $S_{1,t}$ and fitted values are depicted in (ref). Starting from top of the figure shows, in the left panel, the transition function $S_{1,t}$. Comparing the transition function with the outcome (right panel, black line) reveals that $S_{1,t}$ becomes large (approaches $1$) if IP growth is (strongly) negative. When we consider the fitted values, defined as:
we find that the first function already captures a considerable amount of variation in $y_t$. In particular, it succeeds in matching the slowly evolving local trends in IP growth. But it fails to capture much of the idiosyncratic behavior and, in particular, the substantial decline in IP growth during the 2008/2009 global financial crisis (GFC).
Consider adding a second component function. In this case, we fix $\nu_2=10$ and let $\mu_2$ be equal to the $0.99$ quantile of the EBP. With these parameter values, the transition function reduces to the indicator function that equals one if the EBP exceeds its $99$ percent quantile. Considering the transition function reveals that this is only the case during the GFC. In all other periods, the corresponding transition function is (almost) equal to zero with estimated parameters $\beta_{0, 2}= - 4.2$ and $\beta_{1, 2}=0.1$. Since our model is additive this implies that during the GFC, the growth rate is shifted downwards to reach approximately $-5.1$ percent. Notice that a model with $J=2$ component functions is already capable of learning a great deal of variation in IP growth.
Increasing $J$ beyond two further improves the fit, but only slightly so. Using $J=3$ or $J=4$ (with transition functions being informed by the EBP in the case of $J=3$ and lagged IP growth for $J=4$) indicates that the estimated model fit displays more high frequency variation (in consistence with the actual time series). The key question, empirically, however is whether capturing more high frequency noise pays off for predictive performance. In our simulation study, we will return to this question and analyze the relationship between $J$ and predictive performance in more detail.
To sum up, in this simple toy example we find that summing over two logistic functions already provides a decent model fit. The first function, which is a smooth logistic function, explains low frequency trends whereas the second function captures the abrupt downturn during the GFC.
Next we define the likelihood function of our model. To simplify the exposition, we let $\bm Z_t$ denote a $2J$-dimensional vector of (generated) regressors so that:
In this case, we can rewrite (ref) as follows:
with $\bm \beta = (\beta_{0, 1}, \beta_{1,1}, \dots, \beta_{0, J}, \beta_{1, J})'$ being a vector of stacked coefficients. Stacking over $t$ gives rise to the full-data representation of the model:
where $\bm y =(y_1, \dots, y_T)', \bm Z = (\bm Z_1, \dots, \bm Z_T)'$ and $\bm \varepsilon = (\varepsilon_1, \dots, \varepsilon_T)'$ are $T \times 1$, $T \times 2J$ and $T \times 1$ matrices, respectively.
Standard textbook results chan2019bayesian show that the likelihood function can be rewritten as:
Here, we let $\hat{\bm \beta} = (\bm Z' \bm Z)' \bm Z' \bm y$ denote the OLS/maximum likelihood estimator of $\bm \beta$, $w = T - 2J$ the degrees of freedom, and $s^2 = \frac{(\bm y - \bm Z \hat{\bm \beta})' (\bm y - \bm Z \hat{\bm \beta})}{w}$ is the OLS estimator of the error variance. Notice that since $\bm Z$ depends on the speed of adjustment parameters, thresholds and threshold variables, we do not condition on $\{\nu_j\}, \{\mu_j\}$ and $\{ \bm \delta_j \}$. The likelihood function consists of two terms. The first term implies a dependence between $\bm \beta$ and $\sigma^2$ whereas the second term is independent of $\bm \beta$ and looks like the kernel of an inverse Gamma distribution. We will use these observations to construct a standard conjugate prior in the next section.
The model in (ref) might be subject to overfitting if $J$ is set too large. Hence, we need to regularize the estimates of $\bm \beta$. This is achieved through shrinkage priors that are inspired by the priors stipulated in chipman2010bart. Our joint prior on the parameters of the model can be factorized as follows:
Note that the prior on $\bm \beta$ depends on $\sigma^2$ while the priors on the other parameters are independent of each other. We assume that $p(\bm \beta|\sigma^2)$ is Gaussian:
where $\underline{\bm V} = \phi J^{-1} \times \bm I_{2J}$ and $\phi$ is a positive prior scaling parameter. The prior variance decreases in $J$ and hence, for a large number of component functions, we shrink the parameters stronger to zero so that each function is expected to contribute less to explain the variation in $y_t$. Give that $\bm Z_t$ is bounded between $0$ and $1$, we set $\phi=1$.
On the error variances we use the usual inverse Gamma prior $p(\sigma^2) = \mathcal{G}^{-1}(\underline{a}_\sigma, \underline{b}_\sigma)$. We let $\underline{a}_\sigma$ and $\underline{b}_\sigma$ denote the prior degree of freedom and a prior scaling parameter, respectively. To render this prior effectively uninformative, we set $\underline{a}_\sigma = \underline{b}_\sigma = 0.01$. Notice that one could also use a data-driven prior that is scaled by, e.g., the OLS standard deviation or other estimates of the error variances. If this estimate implies under-dispersion one could then place more weight on the prior to shrink the error variances towards zero and thus force the conditional mean to soak up more variation in $y_t$.
For the thresholds we use weakly informative Gaussian priors $p(\mu_j) = \mathcal{N}(0, \underline{\sigma}^2_j)$ where $\underline{\sigma}^2_j$ is a hyperparameter which we set to a large value. In our case, we standardize the input data by subtracting the mean and normalizing by the standard deviation. Hence, the prior centers the threshold over the mean of the non-normalized version of $\Tilde{x}_{j, t}$. Using more informative priors on the thresholds is difficult, in particular given the fact that we do not consider a standard smooth transition model where prior information about possible thresholds could exist.
On $\nu_j$ we use an weakly informative inverse Gamma prior $p(\nu_j) = \mathcal{G}^{-1}(\underline{a}_\nu, \underline{b}_\nu)$ where the hyperparameters $\underline{a}_\nu =\underline{b}_\nu =0.01$ are set close to zero throughout the paper. This choice has been used in, e.g., lopes2006bayesian.\footnote{lubrano2001smooth discusses an alternative based on a truncated Cauchy prior on $\nu_j$. This choice would be straightforward to adopt in our setting.} Finally, we use a discrete uniform prior on $\bm \delta_j$ so that each element is equally likely to be equal to zero.
The prior in (ref) can be combined with the likelihood in (ref) to derive the joint posterior:
This joint posterior takes no well known form. However, given that the prior on $\bm \beta$ and $\sigma^2$ are conditionally (on $\bm Z$) conjugate, we can make use of the fact that $p(\bm \beta, \sigma^2|\bm y, \bm Z)$ takes a well known form koop2003bayesian:
where NIG denotes the Normal-Inverse Gamma distribution with four parameters:
The parameters associated with the transition functions are then sampled using the Bayesian backfitting strategy outlined in hastie2000bayesian and used in chipman2010bart. We let:
denote the partial residual vector and $\bm R_j = (R_{j1}, \dots, R_{jT})'$. Moreover, let $\bm Z_j$ denote a $T \times 2$ matrix with $t^{th}$ row equal to $(S_{jt}, (1-S_{jt}))$ and $\bm \beta_j = (\beta_{0, j}, \beta_{1, j})'$.
The conjugacy of our prior setup implies that we can integrate out $\bm \beta_j$ and $\sigma^2$ to obtain:
Here, we let $\overline{\bm \beta}'_j =\overline{\bm V}_j \bm Z_j \bm R_j$ and $\overline{\bm V}_j = (\bm Z'_j\bm Z_j + J \phi^{-1} \bm I_2)^{-1}$ denote the posterior mean and variance of $\bm \beta_j$, respectively. Notice that we implicitly condition on the other $\bm \beta_s$ and $\bm Z_s$ for $s \neq j$.
We then sample $p(\nu_j, \mu_j, \bm \delta_j|\bm R_j) \propto p(\bm R_j | \nu_j, \mu_j, \bm \delta_j) \times p(\nu_j, \mu_j, \bm \delta_j)$ in two blocks. First, we let $\Tilde{\delta}_j \in \{1, \dots, K\}$ denote a categorical auxiliary variable that indicates the element in $\bm \delta_j$ which equals one. The posterior probability that $\Tilde{\delta}_j = i$ is then given by:
and we can easily compute (ref) for all $j=1, \dots, K$.
Conditional on $\bm \delta_j$ we sample $\nu_j$ and $\mu_j$ jointly using a single random walk Metropolis Hastings step where we propose $(\nu^*_j, \mu^*_j)' \sim \mathcal{N}((\nu^a_j, \mu^a_j)', \text{diag}(s_\nu, s_\mu))$, with the superscript $a$ indicating the previous accepted draw. The scaling parameters of the proposal distribution are tuned during the first half of the burn-in stage of our algorithm so that the acceptance probability is between 30 and 60 percent. After proposing $\nu^*_j, \mu^*_j$, we accept the proposed values with probability equal to:
This completes the different steps to sample from the relevant full conditional posterior distributions. Since we sample some parameters marginal of the others, the ordering of the steps of the sampler play an important role van2008partially. Taking this into account, our algorithm cycles between the following steps:
The first two steps are marginal of $\bm \beta$ and $\sigma^2$ while Step 3 is conditional on $\bm Z$ (i.e. $\{\nu_j, \mu_j, \bm \delta_j\}$) but marginal of $\bm \beta$. The final step is conditional on the error variances and the different component functions. The key property of this algorithm is that we exploit conjugacy to sample the parameters of the component functions independently from the error variances and the regression coefficients. This improves mixing substantially and we found that in our empirical work and the simulations that our algorithm converges rapidly towards the desired stationary distribution.
In this section we put our proposed model to a test within a controlled environment. In particular, we show that under a nonlinear data generating process (DGP), our model yields predictions which are accurate and can compete with the ones of BART. The reason why we benchmark the results to BART is due to the empirical success of BART across many fields. We use a standard BART implementation with precisely the same set of hyperparameters on the trees and the error variances as in chipman2010bart.
We assume that $\{y_t\}_{t=1}^{T=300}$ is generated as follows:
where $\bm x_t \sim \mathcal{N}(\bm 0, \bm I)$ is a $K=25$-dimensional vector, $ \bm \beta_{\text{true}} \sim \mathcal{N}(3, 9)$ is a $25$-dimensional vector of true linear coefficients and $\bm \kappa_{\text{true}} \sim \mathcal{N}(2, 9)$ is a $K-$dimensional vector of nonlinear coefficients. To have a sparse model we zero out 60% of the elements in both $\bm \beta_{\text{true}}$ and $\bm \kappa_{\text{true}}$. Finally, we initialize $y_0 = 0$. This DGP produces time series that match patterns commonly observed in macroeconomics and finance.
We estimate four variants of the AST model for different values of $J$. The first estimates $\nu_j$ and $\mu_j$ using the prior setup discussed in the previous section. The second fixes $\nu_j = 10$, leading to a model that sums over mixtures connected by a threshold function. The third assumes $\mu_j = \hat{\mu}_j = \sum_{t=1}^T \bm \delta_j \bm x_t/T$, implying that the threshold is the empirical mean of the corresponding variable selected by $\bm \delta_j$. Finally, the last specification fixes $\nu_j = 1$ and $\mu_j = \hat{\mu}_j$, leading to a model which sets $S_{jt} = 1$ if $\tilde{x}_{j, t}$ exceeds its mean. The first model is the most flexible one and allows for different threshold values and different speed of adjustments of the transition functions. The last one introduces strong restrictions. The intermediate specifications provide slightly more flexibility and by doing so reduce the number of free parameters. In our simulation we investigate how these choices impact the point and density forecasting performance. To simulate a high dimensional setting, we include four lags of the regressors.
We carry out our forecasting exercise by taking each generated series $\{y_t\}$ and splitting it into two halves of equal size. The first half $t=1, \dots, T_0 ( = T/2 = 150)$ is used to train each model whereas we predict the second half $T_0+1 (=151), \dots, T (=300)$. To speed up computation and due to the fact that our DGP features no structural breaks, we only estimate the model once and then compute one-step-ahead predictions. To control for sampling uncertainty with respect to the DGP we repeat these experiments $50$ times.
To analyze forecast accuracy we compute the root mean squared error (RMSE) as follows:
where $\overline{y}_{t|t-1}$ denotes the median of the one-step-ahead predictive density. To measure the accuracy of density forecasts we compute the log predictive likelihood (LPL) using a Gaussian approximation:
with $p(y_t |\overline{y}_{t|t-1}, \sigma^2_{t|t-1})$ being the predictive distribution evaluated at the actual outcome and $\overline{\sigma}^2_{t|t-1}$ denoting the predictive variance.
Table (ref) shows the results of this simulation exercise. The upper panel of the table shows the RMSEs relative to BART so that numbers greater than one suggest a weaker point forecasting accuracy whereas number smaller than one point towards outperformance of a corresponding AST model. The lower panel shows differences in LPLs between a given AST specification and the BART benchmark, with numbers greater than zero suggesting more precise density forecasts and negative numbers point towards a weaker average density forecast performance.
Considering RMSE results reveals that our baseline specification that estimates $\nu_j$ and $\mu_j$ yields forecasts that improve upon the BART forecasts for $J$ between five and $25$. The improvements in relative RMSEs are U-shaped and first increase until $J=10$, becoming smaller afterwards. At a first glance, this suggests that careful selection of $J$ is necessary to produce accurate forecasts. However, it is worth stressing that BART-based forecasts are typically very precise and our approach, being simpler to implement and, as we will see in the next sections, more scalable, never loses against BART as long as $J > 1$.
If we consider the specification that fixes the threshold parameters, we find a weaker overall performance but RMSE ratios are still $\pm 10$ percent within the absolute RMSEs of BART for all values of $J$. Similar results arise if we fix the speed of adjustment parameters but estimate thresholds. In this case, the $J=1$ case performs poorly. This is expected given that this model is a simple switching model with endogenous selected threshold variable and estimated threshold. If we increase the number of component functions the performance increases until $J=25$. Finally, not estimating $\nu_j$ nor $\mu_j$ is not a good idea. In this case, we lose against the BART benchmark by large margins.
Next, we consider the density forecasting performance in the lower panel of (ref). Recall that numbers greater than zero indicate outperformance of AST whereas negative numbers suggest the opposite. In principle, the density forecasting results tell a story similar (but slightly more pronounced) to the RMSE results. Depending on the choice of $J$, AST improves upon BART and the most flexible version does best on average. The key difference, however, is that for the model that estimates $\nu_j$ and $\mu_j$ we find improvements in accuracy for all values of $J$. But, similar to the findings for point forecasts, these improvements first increase with $J$ and then slowly decay. Out of the restricted versions we find that the model which estimates the speed of transition parameter performs best, yielding gains for different values of $J$. The one that performs worst is, again, the model that fixes both $\nu_j$ and $\mu_j$.
Our previous discussion has established that AST yields forecasts which are often better than the ones produced by BART. A general conclusion stemming from our synthetic data exercise is that for the most flexible version, setting $J>1$ yields point forecasts which are, in the worst case, very close to the ones produced by BART and always produces slightly more accurate density predictions. Another question relevant for practitioners, however, is whether the model performs well in selecting the correct covariates. This is what we investigate in (ref) by looking at a particular realization from the DGP.
The first two columns of the table show the actual values of $\bm \beta_{\text{true}}$ and $\bm \kappa_{\text{true}}$ and the remaining columns show variable relevance scores for the different lags of $\bm x_t$. These are computed by taking the posterior mean of $\bm \delta_i, \overline{\bm \delta}_i$, and then summing over all $i$. A given number hence indicates how often a variable shows up in all component functions and greater numbers thus imply a higher variable relevance. If a score is close to one it implies that only one of the base functions includes a given variable in $\bm x_t$.
At a very general level, we find a close association between regressors that feature large values of $\bm \beta_{true}$ and/or $\bm \kappa_{true}$. If this is the case, variable relevance scores are often above one. There are some rare cases where this does not hold (such as $x_{13,t}$ and $x_{14,t}$), but for the vast majority our model attributes appreciable relevance to covariates that feature large coefficients (in absolute terms). Variables that do not enter the DGP are, without any exception, never included in the corresponding base learners and thus do not impact our model. The single most variable, as expected, is the first lag of the endogenous variable, which shows up almost four out of $J$ times in the corresponding functions.
To shed light on the differences in the predictive densities across different values of $J$, (ref) plots the in-sample and out-of-sample predictive densities for different values of $J$ for the baseline model that estimates both $\nu_j$ and $\mu_j$. In both cases, the results reveal that the model does a good job in fitting the data, irrespective of the choice of $J$. In principle, there are no discernible differences in terms of the posterior medians. The only feature that stands out is that credible sets become smoother and slightly tighter for larger values of $J$ in-sample and, to a somewhat lesser degree, out-of-sample.
Next, we investigate the relationship between $J$ and the shape of the transition functions. To this end, we compute the average transition function based on the (normalized) predictors. This is achieved as follows. For each $j=1, \dots, J$, we take the posterior median of $\nu_j, \mu_j$ and $\bm \delta_j$ and plot the transition function for $\tilde{x}_{j, t}$ ranging from $-10$ to $10$. This gives rise to $J$ different transition functions. When then simply compute the average across these $J$ transition functions and plot these. Hence, the resulting average transition function reflects how the average base learner moves from $S_{j,t} = 0$ to $S_{j, t}=1$.
(ref) shows the shape of these average functions. The single most striking observation is that the speed of adjustment parameter seems to increase with $J$. Whereas we find a rather smooth transition for $J=1$ and $J=5$, going from $J=10$ to $J=15$ implies a transition that is very close to using an indicator function. This indicates that if we use only few base functions to learn the conditional mean relations, our algorithm places substantial posterior mass on transition functions that feature more complex patterns. But for larger $J$, the individual functions become simpler. It is, however, worth stressing that the average transition function for $J=50$ is still not exactly equal to an indicator function and still implies a somewhat gradual transition between regimes for values of $\tilde{x}_{j,t}$ close to the mean.
In the previous sections we have developed the AST model and illustrated its usefulness in Monte Carlo simulations. In this section, we generalize the model to the multivariate case and develop a scalable, conjugate version of it to model a possible large panel of $M$ macroeconomic time series which we store in $\bm y_t$. This model is henceforth labeled the vector additive smooth transition (VAST) model.
We assume that $\bm y_t$ depends nonlinearily on its $P$ lags. These are stored in a $K$-dimensional vector $\bm x_t = (\bm y'_{t-1}, \dots, \bm y'_{t-P})'$ with $K = MP$. The vector additive smooth transition (VAST) model is then given by:
where $g: \mathbb{R} \to \mathbb{R}^M$ is a function that maps a scalar input $\tilde{x}_{j, t}$ into an $M$-dimensional output and $\bm \varepsilon_t$ is a Gaussian white noise process with zero mean and covariance matrix $\bm \Sigma$.
The component function $g$ takes the following form:
This transition function looks similar to (ref) but the location parameters $\bm \beta_{i, j} = (\beta_{ij, 1}, \dots, \beta_{ij, M})'$ are now $M$-dimensional vectors. We will again assume that $S_{j, t}$ takes precisely the same form as (ref).
Under Eqs. ((ref)) to ((ref)), the model can be written as:
where $\bm Z_t$ is given by (ref) and $\bm \beta = (\bm \beta'_{0, 1}, \bm \beta'_{1, 1}, \dots, \bm \beta'_{0, J}, \bm \beta'_{1, J'})'$ is a $N-$dimensional vector with $N=2JM$.
The number of free parameters in the model is $v^{\text{VAST}} = J(3 + 2M) + M (M+1)/2$. This number, for moderate values of $J$, is much smaller than $v^{\text{VAR}} = M^2P + M (M+1)/2$, the number of parameters of an unrestricted (but linear) VAR. Notice that the Kronecker structure implies that each equation in the model features the same set of nonlinear transformations of selected covariates. At a first glance, this assumption might be restrictive but if $J$ is set to be large, the model is still flexible enough to capture arbitrary nonlinearities across equations and, specifically, equation-specific idiosyncrasies in terms of nonlinear behavior of the time series. This is because the corresponding equation-specific parameters for each $\bm Z_t$ can differ. So in case that there is strong evidence that one (or more) variable(s) in the system evolve according to, e.g., a threshold process, our algorithm would add appropriate base learners to the conditional mean model. In this case, the corresponding coefficients would be non-zero whereas the coefficients associated with other transformations would then be close to zero.
In terms of computation, the Kronecker structure in the likelihood gives rise to substantial computational advantages. This not only relates to posterior sampling (see (ref)) but also to the computation of generalized impulse responses (GIRFs), see (ref) of the Online Appendix. GIRF computation in models such as the BART-VAR of huber2020inference require computing forecast distributions (both unconditional and conditional on a shock of interest). Since nonlinear models imply that GIRFs are state-dependent, one needs to integrate over the economic conditions. If each equation is determined by its own equation-specific function, this becomes excessively slow and turns out to be the computational bottleneck in these models. The reason is that each equation-specific function, $f_j(\bm x_t)$, needs to be approximated and, for large $M$, the computational burden becomes large. By contrast, our approach only requires us to compute $g(\Tilde{x}_{j, t}|\bm \theta_j)$ for all $j$ and then use the location coefficients to obtain a draw from the conditional mean for $\bm y_t$. Hence, we do not need to evaluate $g$ for each equation due to the Kronecker structure. And this translates into substantial speed improvements when it comes to computing nonlinear functions such as GIRFs.
Most priors and steps in the posterior simulator remain untouched by moving from the univariate to the multivariate model. Hence, we briefly summarize differences in priors first and then discuss differences to the MCMC algorithm sketched in (ref)
The priors of the model exactly resemble the ones used for the univariate AST model with two exceptions. We use a Gaussian prior on $\bm \beta$ that conditions on $\bm \Sigma$:
The prior covariance matrix thus features a Kronecker structure similar to the one in the likelihood. Again, we set $\underline{\bm V} = \phi J^{-1} \times \bm I_{2J}$ and set $\phi = 1$.
The prior on $\bm \Sigma$ is inverse Wishart:
with $\underline{a}_\Sigma$ denoting prior degrees of freedom and $ \underline{\bm S}_\Sigma$ is a prior scaling matrix. We set $\underline{a}_\Sigma = M$ and $\underline{S}_\Sigma = 1/100 \times \bm I$. This choice yields a proper prior that is relatively uninformative. If one wishes to force a more aggressive model fit, one could set $\underline{\bm S}_\Sigma$ equal to a variance estimator that would imply overfitting and place more weight on the prior by increasing the prior degrees of freedom.
The posterior simulator differs in three respects from the one associated with the univariate model. First, the particular form of $p(\bm R_j|\nu_j, \mu_j, \bm \delta_j)$ differs, where $\bm R_j$ is now $T \times M$ matrix defined with the $t^{th}$ row given by $(\bm y_t - \sum_{s \neq j} g(\tilde{x}_{s, t}|\bm \theta_s))'$. When we integrate over $\bm \beta$ and $\bm \Sigma$, we end up with the following standard expression for the marginal likelihood of the Bayesian seemingly unrelated regression (SUR) model:
where $\overline{\bm \beta}_j$ and $\overline{\bm V}_j$ is defined below (ref). This expression is used to set up the Metropolis Hastings updates or inverse transform steps employed to sample the thresholds, threshold variables and speed of adjustment parameters.
The next difference relates to how we sample the regression coefficients $\bm \beta$. The full conditional posterior of the multivariate model takes the following form:
with $\overline{\bm V}$ being defined as before, $\bm Y = (\bm y_1, \dots, \bm y_T)'$ and $\overline{\bm \beta} = \text{vec}(\overline{\bm V}_\beta \bm Z' \bm Y)$.
Finally, the posterior of $\bm \Sigma$ is inverse Wishart:
The resulting MCMC algorithm closely resembles the one discussed in (ref) with the sampling steps for $\bm \beta$, $\bm \Sigma$ and the acceptance/posterior probabilities adjusted accordingly.
It is worth stressing that this algorithm is only slightly more costly than the one for the univariate model. In particular, sampling from the posterior of $\bm \beta$ is more expensive but the Kronecker structure implies that high dimensional matrix operations can be avoided. Hence, sampling from $p(\bm \beta|\bm \Sigma, \bm Y, \bm Z)$ is fast and one can easily estimate nonlinear VARs with more than 100 equations.
In this section we apply the VAST model to US macroeconomic data. We start by providing a brief overview on the dataset and then move on to provide some evidence on the predictive performance of our model. Finally, we discuss how the US economy reacts to financial shocks.
We apply the VAST model to the FRED-QD dataset mccracken2020fred. Our sample runs from $1973$Q1 to $2019$Q4. In $\bm y_t$, we include $M=80$ variables. These are given in Table (ref). Notice that this set of variables implies that we include a large number of quantities that measure the real side of the economy as well as several factors that capture movements in financial markets. When we consider the effects of financial shocks on the US economy, we also add the EBP stipulated in gilchrist2012credit as a measure of financial conditions. Since this series is only available up to $2016$Q4, we use a slightly shorter sample for the structural analysis.
In our forecasting exercise we also consider two smaller-sized datasets. These are formed as sub-groups out of this large-scale dataset and defined in Footnote (ref). For the predictive exercise, we drop the EBP to use data through $2019$Q4.
All the models we consider include $p=5$ lags of $\bm y_t$. The number of base learners $J$ is set equal to $50$ when we discuss full sample results (such as the ones in the next sub-section and Sub-section (ref)). We analyze predictive performance over $J$ and find that setting it equal to $40$ or $50$ generally yields the best density forecasting performance. In our structural analysis, we find that changing $J$ leads to impulse responses which are similar in qualitative terms.
In this section, we analyze whether our VAST model is capable of outperforming the BART-VAR proposed in clark2023tail.\footnote{The setup is precisely the same as the one used in clark2023tail.} We include this model because, on a very similar dataset, we have shown that it works well for density and tail forecasts, often improving upon a BVAR with SV and never being substantially outperformed by the BVAR-SV. To analyze the relationship between model size, density forecasting performance and $J$, we consider three different model sizes and set $J \in \{10, 15, 20, 25, 30, 40, 50\}$. The model sizes we consider are a small-scale (S) model that includes $M=3$ variables. These are the unemployment rate (UNRATE), CPI inflation (CPIAUCSL) and the Federal Funds Rate (FEDFUNDS). The next larger model is a medium-sized (M) one that includes $M=23$ variables. This model uses the small dataset and adds additional real quantities and financial market variables.\footnote{More precisely, we include the following series from the FRED-QD database: GDPC1, PCECC96, FPIx, GCEC1, INDPRO, CE16OV, UNRATE, CES0600000007, HOUST, PERMIT, PCECTPI, PCEPILFE, GDPCTPI, CPIAUCSL, CPILFESL, CES0600000008, FEDFUNDS, GS1, GS10, M2REAL, TOTRESNS, NONBORRES, S.P.500. The definition of the different abbreviations is given in (ref). } The large-scale (L) dataset is the one described in (ref) bar the EBP and thus $M=79$.
We use a recursive forecasting design that starts with using data through $1989$Q4 to initially train the models. We then compute one-quarter-ahead forecast distributions for $1990$Q1 and evaluate these at the actual outcome using log predictive likelihoods (LPLs). After obtaining the LPLs for $1990$Q1, we add this data point to the training sample and estimate the one-quarter-ahead predictive density for $1990$Q2 and compute the corresponding LPLs. We repeat this procedure until we reach $2019$Q3 and thus compute the forecasts for $2019$Q4 (the end of the sample). This yields a sequence of time-specific LPLs which we average to end up with the average LPLs we report in (ref). This table includes differences between the average LPLs of a particular model and the BART-VAR of a given size. There are two types of LPLs in the table. One is the marginal LPL for a particular focus variable (UNRATE, CPIAUCSL or FEDFUNDS) whereas the second is the joint LPL for the three focus variables.
The table indicates that, for specific values of $J$, VAST is capable of improving upon the BART-VAR for almost all target variables and most model sizes. In particular, we find gains that range from being almost zero (such as for small datasets and unemployment forecasts) to modest (such as for inflation arising from small datasets or interest rate forecasts and large datasets).
For inflation forecasts using medium and large-sized datasets, we find that VAST does not outperform the BART-VAR. But in these cases, setting $J$ either to $40$ or $50$ yields LPLs that are almost identical to the one of the benchmark model. Another case where the BART-VAR produces more accurate density forecasts is the FEDFUNDS rate when the small dataset is adopted.
When we focus on the joint forecasting performance a similar picture arises. We find that VAST yields improvements for small and large models and similar predictive likelihoods for the medium-sized models if $J$ exceeds $10$. For small models, these improvements are muted but consistent across different values of $J$. For the large model, we find that forecast performance varies with $J$ and larger values of $J$ translate into the most precise density forecasts. These joint density forecasts are obtained when we set $J=40$ or $50$. This finding is not surprising given that a larger number of base learners improves flexibility to capture equation-specific nonlinear patterns.
Next we ask whether the forecasting performance in terms of one-quarter-ahead joint LPLs is heterogenous over time. To do so, we consider the model with $J=40$ as this specification performs well across all focus variables and for joint LPLs and benchmark it against the large BART-VAR. To understand how performance changes over time, we compute cumulative relative LPLs to the large BART-VAR for VAST across the three different model sizes.
The results are depicted in (ref). The figure shows that VAST ($J=40$) outperforms the large BART-VAR consistently across all model sizes. Only during the GFC in 2008/2009 we find a slight decline in relative forecasting accuracy for the small and large VAST specifications. In this period, interestingly the medium-scale model outperforms the BART-VAR and relative model performance increases. Apart from the GFC, the consistent outperformance of VAST remains visible throughout the hold-out period.
This discussion has shown that VAST can improve upon BART, a very competitive benchmark model that has a proven track record in density forecasting. If it is outperformed by the BART-based VAR, the losses in predictive accuracy are typically quite limited. In light of this, it is worth stressing that obtaining the predictive densities of VAST is quick relative to the benchmark. Producing the one-step-ahead density for a particular point in the hold-out takes around five minutes for the large model whereas for the BART-VAR it takes over 1.5 hours on a state-of-the-art Macbook pro.
Next, we turn to the analysis of nonlinearities in the transmission of financial shocks to the US economy. This issue has gained increasing attention in the recent literature barnichon2022effects, mumtaz2022impulse, forniJMCB. Most of these studies find that uncertainty shocks trigger important effects only when they are contractionary and sizable. In all other cases, the effects appear to be muted. This is in stark contrast to the literature utilizing linear VARs gilchrist2012credit which find sizable reactions to financial shocks. This is because linear models mix over positive and negative shocks and thus over-exaggerate the effect of a benign shock while underestimating the effect of contractionary shock.
To identify the shock we use zero restrictions that imply that real variables and the Federal Funds rate react sluggishly with respect to a financial shock while financial markets react immediately. This choice is consistent with other papers del2020sthat deal with estimating the effects of financial shocks in multivariate time series models. We use the large dataset with $M=80$ (the $79$ macro series and the EBP) and set $J=50$.
Our empirical focus will be on two forms of asymmetries. The first is whether benign and adverse shocks trigger different reactions of $\bm y_t$. The second form is whether small shocks trigger different reactions from large shocks. Since our model is nonlinear, asymmetries here could imply that shocks are disproportionally stronger for larger shock sizes or that the shape of the impulse responses differ for small versus large shocks. We consider a one standard deviation (S.D.) shock to be a small shock whereas a five S.D. shock is perceived as a large shock.
We first discuss the endogenous reaction of the EBP to financial shocks. gilchrist2012credit argue that fluctuations in the EBP represent movements in investor sentiment or changes in risk preferences in the corporate bond market. (ref), panels (a) and (b), shows the reaction of the EBP to a small (panel (a)) and a large (panel (b)) financial shock.
Starting with panel (a) reveals that, if the shock is small, benign and adverse financial shocks trigger a symmetric increase in the EBP which slowly fades out, turning insignificant after around 8 to 10 quarters. By contrast, if the shock size becomes large, we find differences in the shapes of of the EBP reactions. A large and benign financial shock induces a strong immediate reaction that abruptly dies out. An adverse financial shock translates into a strong but more persistent increase in the EBP. Notice that posterior uncertainty is slightly smaller in the case of a benign shock.
The reaction of the EBP indicates symmetric responses to small shocks (irrespective of sign) but increasing asymmetries if the shock is sizable. This finding carries over to the impulse responses of GDP growth and the unemployment rate in (ref). In both cases, small financial shocks of either sign trigger reactions of GDP growth and the unemployment rate with the correct sign (i.e. decreasing levels of real activity if the shock is adverse and increasing levels of real activity if the shock is benign). When we consider the effects of large shocks we find substantial asymmetries. In particular, both GDP growth and the unemployment rate display a much stronger reaction to an adverse shock. This is consistent with, e.g., barnichon2022effects and mumtaz2022impulse, who also document stronger reactions of output growth to contractionary financial shocks. Apart from the stronger short-run effects, we also find that reactions to an adverse shock are slightly more persistent and turn insignificant around two years after the shock hit the system.
(ref) shows the responses of short-term interest rates and 10-year US treasury yields. As opposed to output growth and unemployment reactions, the Federal Funds Rate and the 10-year yield reactions feature some asymmetries. For the Federal Funds rate, these asymmetries relate to short-run responses (between one to 1.5 years), with the central bank displaying a stronger reaction in response to a contractionary financial shock. In this case, short-rates are decreased by around 15 basis points (bps) whereas in the benign case, the central bank lowers short-term interest rates by around 10 bps. Treasury yields, by contrast, display a slightly stronger impact reaction to a benign shock but, in the short-run, the decline in response to adverse financial shocks is stronger (in absolute terms) than the increase in yields in response to benign shocks.
When we consider large shocks, we find substantial evidence for asymmetries. The Federal Funds rate declines by around 45 bps after one year in response to an adverse shock. The reaction to a benign shock is much more muted, with median increases in short-term interest rates of around 22 bps. For treasury yields, we find similar impact reactions to a large shock (with flipped sign). This is in contrast to a small shock. However, treasury markets exhibit a much stronger reaction to adverse shocks than to benign shocks during the first two years.
Finally, we consider how other financial markets react to financial shocks. In particular, we focus on the reactions of BAA-rated corporate bond yields and the S&P 500. Responses of BAA-rated bond yields to small financial shocks point towards no asymmetries for immediate reactions but slightly stronger reactions after around four quarters. Stock markets, by contrast, react symmetrically to small shocks of either sign. When we consider larger shocks we, again, find substantial asymmetries. Reactions to large adverse financial shocks are stronger, in particular between two and 8 quarters. These results, in combination with the reaction of the 10-year treasury yields, can be interpreted in the sense that financial markets are much more reactive to contractionary financial shocks, triggering a 'risk-off' mood of investors. This is reflected by declining stock market returns, increases in BAA-rated bond yields and declines in treasury yields which act as a safe asset.
To sum up, our findings confirm findings in the literature that reactions to adverse shocks are stronger than the ones to benign shocks. However, we also find that this asymmetry result only arises if the shock becomes large.
The main goal of this paper is to develop a highly scalable yet flexible econometric model that is capable of capturing asymmetries in possible macroeconomic relations. We achieve this by summing over $J$ simple functions which are location mixtures with transition between regimes driven by a logistic function. Monte Carlo evidence suggests that, in the univariate case, our approach produces predictions which are close to the one of BART, a popular machine learning tool, and in some cases slightly more precise. We generalize our approach to the multivariate case, leading to the VAST. This model is highly scalable and can be applied to large datasets.
After showing that our approach works well in predictive terms, we apply the VAST to a dataset with $80$ endogenous variables and consider the asymmetries in the responses to financial shocks. Our results indicate that macro reactions to adverse financial shocks are asymmetric with respect to sign and shock. In particular, we find that positive and negative financial shocks trigger similar reactions if the shock is small. But if the shock becomes large, adverse shocks lead to much stronger reactions.
{\setstretch{1.1} \addcontentsline{toc}{section}{References} } \setcounter{page}{1} \counterwithin{figure}{section} \counterwithin{table}{section}