Time-Varying Parameters as Ridge Regressions

-0.54cm Time-Varying Parameters as Ridge Regressions

}

\thispagestyle{empty}

\pagenumbering{arabic}

Introduction

Economies change. Intuitively, this should percolate to the parameters of models characterizing them. To econometrically achieve that, a popular approach is Time-Varying Parameters (TVPs), where a linear equation’s coefficients follow stochastic processes --- usually random walks. {\color{black} There is a wide body of work utilizing TVPs to study structural changes in key macroeconomic equations, such as Phillips curves and Taylor rules stock1996evidence,boivin2005has,blanchard2015inflation. In the multivariate setting, } classic papers consider TVP Vector Autoregressions (TVP-VARs) to examine changing monetary policy primiceri2005 and evolving inflation dynamics cogley2001evolving,cogley2005drifting. Recently, these ideas have been extended to jorda2005's local projections (LPs) to obtain directly time-varying impulse response functions ruisi2019time,lusompa2020local.

In both the VAR and LP cases, important practical obstacles reduce the scope and applicability of TVPs. One is prohibitive computations limiting the size of the model. Another is the difficulty of tuning the crucial amount of time variation. To address those and other problems, I leverage the underutilized fact that TVP models are ridge regressions (RR). The connection is useful: 50 years hoerl1970ridge of RR widespread use, research and wisdom is readily transferable to TVPs. {\color{black}Among other things, this provides fast computations via a closed-form dual solution only using matrix operations. Therefore, it avoids potential initialization and convergence issues because it does not rely on MCMC simulations and filtering. Next, we have that the amount of time variation is automatically tuned by cross-validation (CV). Adjustments for evolving residuals' volatility and heterogeneous parameter drifting speeds (random walk variances) are implemented via a 2-step ridge regression (henceforth 2SRR), the flagship model of this paper. As a result, it is easy to implement and operate, and it integrates seamlessly into existing machine learning-based macroeconomic forecasting pipelines. } Moreover, the setup is directly extendable to deploy additional shrinkage schemes (sparse TVP, reduced rank TVPs) recently proposed in the literature stevanovic2016common,bitto2018achieving.

{\color{black} In sum, by introducing 2SRR and related methods, the paper's key contributions are threefold: (i) easier computations, (ii) easier tuning, and (iii) a deeper understanding of standard TVP models used in empirical macroeconomics through the lenses of ridge regression. Since steps have already been taken in the direction of each of these contributions, I will now review the relevant literature and discuss where the ridge approach sits.}

\vskip 0.2cm

{\sc Computations.} Using standard implementations allowing for stochastic volatility (SV) in TVP-VARs, researchers are limited to a few lags and a small number of variables (not more than 4 or 5) kilian2017SVAR. This constraint leaves the reader ever-wondering whether time variation is not merely the byproduct of omitted variables. A similar problem occurs if one wants to study time-varying local projections of a certain size, which imply repeated estimation over many horizons. Consequently, a growing number of contributions attempt to deal with the computational problem. Within the state-space paradigm, koop2013large and huber2020bayesian proposed approximations to speed up MCMC simulations. giraitis2014inference and kapetanios2017large drop the state-space paradigm altogether in favor of a nonparametric kernel-based estimator. chen2012testing consider a similar approach to develop a test for smooth structural change while petrova2016quasi develops a Bayesian version particularly apt with large multivariate systems. Also in the Bayesian arena, hauzenberger2022fast drop the random walk in favor of a hierarchical prior that lends itself well to computations through matrix operations using the singular value decomposition.

This paper's contribution along the computational angle is to provide quick and off-the-shelf point estimates for the original random-walk based model. It demonstrates that this can be achieved using reasonably elementary methods without compromising reliability. However, unlike Bayesian methods, it does not inherently provide inference, necessitating an additional procedure to obtain it.

\vskip 0.2cm

{\sc Tuning.} AAG2013, BaumeisterKilian2014, pettenuzzo2017forecasting, and several others have investigated, with varying angles, the usefulness of time variation to increase prediction accuracy. A critical choice underlying forecasting successes and failures is the amount of time variation. Notoriously, tuning parameter(s) determining it can largely influence prediction results and estimated coefficients, accounting for much of the transparency and reliability concerns regarding TVP models. amir2018choosing and cadonna2020triple propose to treat those pivotal hyperparameters as another layer of parameters to be estimated within the Bayesian procedure --- and find this indeed helps.

By demonstrating the TVP-RR equivalence, this paper contributes to the literature by clarifying further the fundamental tuning problem for this class of models. TVP models are standard high-dimensional regressions that require regularization. When rewritten as a ridge regression, the overall level of regularization is determined by the well-known \(\lambda\). Thus, the proposed solution to the tuning problem is intuitive: \(\lambda\) is determined using standard cross-validation techniques. This measure is crucial for both predictive accuracy and economic analysis, making it reassuring that it can be tuned using the same methods that have been applied to ridge \(\lambda\) for decades. Additionally, since ridge regression can equivalently be rewritten as a constrained optimization problem where the constraint is a "parameter budget constraint"—here translating to a "time-variation budget constraint"—we gain a new perspective on how the amount of time variation should be set. That is, it should balance accommodating relevant time-variation in the data while addressing the reality of short estimation samples.

\vskip 0.2cm

{\sc Interpretation Trinity: Ridge Regressions, Splines, and State Space Models.} The TVP regression is closely related to the Hodrick-Prescott (HP) filter. HP filtering can be seen as a special case of the ridge approach to TVP with only an intercept as the regressor, and it penalizes double rather than first differences. In fact, schlicht2005estimating, paige2010ridge, and yamada2016hodrick rewrite the HP filter as a ridge regression. This connection allows for more principled methods to determine \(\lambda\), which dictates the allocation between the cyclical and trend components. From this perspective, my work extends their contributions to a multiple time-series model with evolving volatility.

Furthermore, there appears to be an artificial division in the recent TVP literature between nonparametric and law-of-motion approaches. By demonstrating that random walk TVPs can be represented as a ridge regression, which is a special case of smoothing splines in this context, it becomes evident that random walk TVPs are just as nonparametric as those derived from the "nonparametric approach." Such equivalences are well-known in the context of local-level models kimeldorf1970, koop2003bayesian, durbin2012time, which are TVP models with only an intercept.

Basis expansion and reparametrization arguments are crucial for establishing the "ridge link", and have seldom been used in related applications. For example, tibshirani2015statistical evokes reparametrization as a possibility to estimate "fused" Lasso via plain Lasso, and korobilis2019high incorporates related ideas as a core component of his "message-passing" algorithm. Additionally, the structure of the expanded matrix in 2SRR mirrors the "indicator-saturation" approach by castle2015detecting for detecting location shifts in the intercept using model selection tools. With the growing use of machine learning methods in time series analysis, the basis expansion and reparametrization ideas developed in this paper are particularly relevant. These concepts can be applied to easily induce random walk parameters in any regularized models, thereby achieving similar benefits beyond the realm of penalized regressions and those with explicitly specified regularization MDTM,HNN.

\vskip 0.2cm

{\sc Extensions and Additional Features.} Extending the two-step RR to a multistep one brings forth two interesting refinements that have been of interest in the literature. Firstly, I consider Sparse TVPs. This means that some parameters will vary and some will not, which can provide efficiency and interpretability gains. belmonte2014hierarchical, korobilis2014data, bitto2018achieving, cadonna2020triple, huber2021inducing, and hauzenberger2024dynamic have proposed such extensions to MCMC-based procedures. In the Ridge paradigm, similar shrinkage is obtainable by iterating 2SRR (thus, a multistep RR) and continuously updating heterogeneous parameter variance estimates.

Another natural way to discipline TVPs is to impose a factor structure. This means that instead of trying to filter, say, 20 independent states, we can span these with a parsimonious set of latent factors. This extension is considered in stevanovic2016common, dewindgambetti2014 and chan-ftvp2018. Such reduced rank restrictions are brought to this paper’s arsenal using a different multistep RR that alternates between loading and a factor estimation steps.

The paper also discusses how to tackle efficiently multivariate models, and in particular, how the same ridge apparatus can be used to obtain a time-varying covariance matrix, a necessary input for structural TVP-VAR analysis. Lastly, a scalable Weighted Bayesian Bootstrap procedure (inspired from newton2018weighted and ng2022random) is proposed to quantify parameter uncertainty taking into account that of $\lambda$.

\vskip 0.2cm

{\sc Simulations and Empirical Results.} I first evaluate the method with simulations. For models of smaller size, where traditional Bayesian procedures can also be used, 2SRR does as well and sometimes better at recovering the true parameter path than the (Bayesian) TVP-VAR. This is true whether SV is involved or not. This is practical given that running and tuning 2SRR for such small models (300 observations, 6 TVPs per equation) takes less than 5 seconds to compute. Since large single equations will be at the heart of the forecasting and empirical application, I also benchmark with a state-of-the-art single equation Bayesian method cadonna2020triple that (i) optimizes the amount of shrinkage and (ii) can work in modestly higher dimensions with an acceptable computing time. Again, 2SRR provides estimates of comparable quality at a fraction of the computational cost. Additionally, I evaluate the performance of 3 variants of the RR approach in a substantive forecasting experiment. 2SRR and its iterated extension provide sizable gains for interest rates and inflation, two variables traditionally associated with the need for time variation.

I complete with an application to estimating large time-varying LPs (more than 4\,600 TVPs) and equally demanding TVP-VARs in a Canadian context using champagnesekkel2018's narrative monetary policy (MP) shocks. It is found that MP shocks' long-run impact on inflation increased substantially starting from the 1990s (onset of inflation targeting), whereas the effects on real activity indicators (GDP, unemployment) became milder. To give a sense of the practicality of the method, estimating the necessary elements to generate a 3D plot of all impulse response functions (IRFs) from a TVP-VAR with 8 variables and 24 (monthly) lags takes 1:34 minutes and requires very limited user input beyond specifying $Y$.

\vskip 0.2cm

{\sc Outline.} Section (ref) presents the ridge approach, its extensions, and related practical issues. Sections (ref) and (ref) report simulations and forecasting results, respectively. Section (ref) applies 2SRR to (large) time-varying LPs. Tables, additional graphs and technical details are in the Appendix.

\vskip 0.2cm

{\sc Notation.} $\beta_{t,k}$ refers to the coefficient on regressor $X_k$ at time $t$. To make things lighter, $\beta_t \in {\rm I\!R}^K$ or $\beta_0 \in {\rm I\!R}^K$ always refers to {all} coefficients at time $t$ or time zero, respectively. Analogously, $\beta_k$ represents the whole time path for the coefficient on $X_k$. $\boldsymbol{\beta} \in {\rm I\!R}^{KT}$ is stacking all $\beta_k$'s one after the other, for $k=1,...,K$. All this also applies to $u$ and $\theta$.

Time-Varying Parameters and Ridge Regressions

This covers the ridge approach, its extensions, and associated practical considerations.

A Useful Observation

Consider a generic linear model with random walk time-varying parameters

where $\beta_t \in {\rm I\!R}^K$, $X_t' \in {\rm I\!R}^K$, $u_t \in {\rm I\!R}^p$ and both $y_t$ and $\epsilon_t$ are scalars. This paper first considers a general single equation time series model and then discuss its generalization to the multivariate case in Section (ref). For clarity, a single equation in a VAR with $M$ variables and $P$ lags has $K=PM+1$ parameters for each equation. For simplicity of exposition, I first impose $\Omega_u = \sigma_{u}^{2} I_K$ and $\sigma_{\epsilon_t}^{2}=\sigma_{\epsilon}^{2} \quad \forall t$. This means all parameters are assumed to vary equally a priori and constant variance of residuals. These assumptions will be relaxed in Section ((ref)). The textbook way of estimating (ref) is maximum likelihood using the Kalman filter for linear Gaussian model hamilton1994. The advantages of the newly proposed methods will be more apparent when considering complications typically encountered in macroeconometric modeling (e.g. evolving volatility, heterogeneous variances for the coefficients paths, and a large $X_t$).

A useful observation is that (ref)'s parameters can be equivalently obtained from solving the penalized regression problem

This is merely an implication of the well-known fact that $l_2$ regularization is equivalent to opting for a standard normal prior on the penalized quantity (see, for instance Sections 7.5-7.6 in murphy2012). Hence, model ((ref)) implicitly poses ${\beta_{t}}-{\beta_{t-1}}\sim{N(0,\sigma_{u}^{2}})$, which is exactly what model ((ref)) also does. Defining $\lambda \equiv \frac{\sigma^2_{\epsilon}}{\sigma^2_{u}}\frac{1}{K}$, the problem has the more familiar look of

The sole hyperparameter of the model is $\lambda$ and it can be tuned by cross-validation.\footnote{This definition of $\lambda$ does not imply it decreases in $K$ since $\sigma^2_{u}$ will typically decrease with $K$ to avoid overfitting.} This model has a closed-form solution as an application of generalized ridge regression hastie2015statistical. In particular, it can be seen as the $l_2$ norm version of the "Fused" Lasso of fusedlasso2005 and embeds the economic assumption that coefficients evolve slowly. However, as currently stated, solving directly (ref) may prove unfeasible even for models of medium size. Note that the layout of the regularization apparatus implies that a priori, \(u_{t}\) and \(\epsilon_{t}\) are serially and mutually uncorrelated. However, estimation results may, for instance, exhibit some persistence in \(u_{t}\), as regularization is a soft rather than a hard constraint.

Getting a Ridge Regression by Reparametrization

The goal of this subsection is to rewrite the problem (ref) as a ridge regression so it is easy and fast to implement in any software. Doing so will prove useful at the conceptual level, but also to dramatically alleviate the computational burden. From now on, it is less tedious to use matrix notation. The fused ridge problem reads as

where $\boldsymbol{D}$ is the first difference operator. $\boldsymbol{W}=[diag(X_1) \quad \dots \quad diag(X_K)]$ is a $T \times KT$ matrix. To make matters more visual, the simple case of $K=2$ and $T=4$ gives rise to $$ \boldsymbol{W}=

. $$ The first step is to reparametrize the problem by using the relationship $\beta_k = C \theta_k$ that we have for all $k$ regressors. $C$ is a lower triangular matrix of ones (for the random walk case) and I define ${\theta_k} = [{u_k} \quad {\beta_{0,k}}]$. For the simple case of one parameter and $T=4$: \[

=

\] For the general case of $K$ parameters, we have $$ \boldsymbol{\beta}= \boldsymbol{C} \boldsymbol{\theta}, \quad \boldsymbol{C} \equiv I_K \otimes C $$ and $\boldsymbol{\theta}$ is just stacking all the $\theta_k$ into one long vector of length $KT$. Note that the summation matrix $C$ could accommodate easily for a wide range of law of motions just by changing summation weights. Actually, any process that can be rewritten (a priori) in terms of uncorrelated $u$'s could be used. For instance, AR models of arbitrary order and RW with drifts would be straightforward to implement.\footnote{In the latter case, it can be shown that one simply needs to add regressors $t*X_t$ to those implied by the RW without drift, that is the $Z_t$'s to be detailed later.} Furthermore, one could use $C^2$ in the random walk setup and obtain smooth second derivatives, e.i. a local-level model. While it is clear that many more exotic configurations are only a $C$ choice away, there is a clear advantage to random walks-based processes: the corresponding $C$ has no parameter to estimate. If we wanted to consider an AR(1) process with a coefficient $\phi \in (0,1]$, either a 2-step estimation procedure or cross-validating $\phi$ would be necessary. Thus, the ridge approach is possible, whether $\beta_t$'s are random walks or not.

Using the reparametrization $\boldsymbol{\beta}= \boldsymbol{C} \boldsymbol{\theta}$, the fused ridge problem becomes

and it is now clear what should be done. Let $\boldsymbol{Z} \equiv \boldsymbol{WC}$ and use the fact that $\boldsymbol{D} = \boldsymbol{C}^{-1}$ to obtain the desired ridge regression problem

Again, for concreteness, the matrix $\boldsymbol{Z} = \boldsymbol{WC}$ looks like $$ \boldsymbol{Z}=

$$ in the $K=2$ and $T=4$ case. The solution to the original problem is thus

This is really just a standard (very) high-dimensional Ridge regression.\footnote{In this Section, I assumed for simplicity that we wish to penalize equally each member of $\theta$ which is not the case in practice. It is easy to see why starting values $\beta_0$ should have different (smaller) penalty weights and this will be relaxed as a special case of the general solution presented in Section (ref).} These derivations are helpful to understand TVPs, which is arguably one of the most popular nonlinearity in modern macroeconometrics. (ref) is equivalent to that of a first-order (multivariate) smoothing splines estimator. More generally, the equivalence between Bayesian stochastic process estimation and splines has been known since kimeldorf1970. Following along, considering a local-level model for $\beta_t$ would yield second order smoothing splines. Clearly, random walk TVPs and their derivatives can hardly be described as "more parametric" than kernel-based approaches: splines methods are prominent within the nonparametric canon. Furthermore, the observation by itself can be leveraged to think about how inducing RW-like structural change in more modern ML tools, like Random Forests MRFjae.

However, we are not there yet, since (ref) implies inverting a $KT \times KT$ matrix, a major computational bottleneck that would limit the applicability of the technique to models of similar size to primiceri2005. Fortunately, there is no need to invert that matrix.

Solving the Dual Problem

The goal of this subsection is to introduce a computationally tractable way of obtaining the ridge estimator $\boldsymbol{\hat{\beta}}$ in (ref). It is well known from the splines literature wahba1990spline and later generalized by scholkopf2001generalized that for a $\boldsymbol{\hat{\theta}}$ that solves problem ((ref)), there exist a $\boldsymbol{\hat{\alpha}} \in {\rm I\!R}^T$ such that $\boldsymbol{\hat{\theta}}=\boldsymbol{Z'\hat{\alpha}}$. While this may appear to be a constraint, it actually arises directly from the first-order conditions of the original problem, just expressed differently. To see this, note that the primal problem (ref) can be converted into the constrained minimization problem

as observed in saunders1998ridge and others. Its Lagrangian is

where $\boldsymbol{a} \in {\rm I\!R}^T$ is a vector of Lagrange multipliers. Setting derivatives with respect to the primal variables $(\boldsymbol{\theta}, \boldsymbol{r})$ to zero, we obtain from first order conditions that the solution should satisfy $\boldsymbol{\theta} = \frac{1}{\lambda} \boldsymbol{Z}' \boldsymbol{a} \quad \text{and} \quad \boldsymbol{r} = -\boldsymbol{a}$. Making these substitutions to eliminate \(\boldsymbol{r}\) and \(\boldsymbol{\theta}\) gives the dual problem

where everything is expressed in terms of \(\boldsymbol{a}\) rather than \(\boldsymbol{\theta}\). We can also reparametrize, defining \(\boldsymbol{\alpha} = \frac{1}{\lambda} \boldsymbol{a}\), and use directly the knowledge about the dual solution (\( \boldsymbol{\theta}= \boldsymbol{Z}' \boldsymbol{\alpha}\)) in the primal problem (ref) to obtain

The solution to the original problem becomes

When the number of observations is smaller than the number of regressors, the dual problem allows to obtain numerically identical estimates by inverting a smaller matrix of size $T$. Since sample sizes for macroeconomic applications quite rarely exceed 700 observations (US monthly data from the 1960s), the need to invert that matrix is not prohibitive. While computational burden does still increase with $K$, it increases much slower since the complexity of matrix multiplication is now $O(KT^3)$ and $O(T^3)$ for matrix inversion. Solving the primal problem, one would be facing $O(K^2T^3)$ and $O(K^3T^3)$ complexities respectively. Concretely, solving the dual problem makes the estimation of high-dimensional TVP models significantly more practical because the computations scale linearly with \(K\) for a fixed \(T\).

Heterogenous Parameter and Residual Variances

For pedagogical purposes, previous sections considered the simpler case of $\Omega_u = \sigma_u^2 I_K$ and no evolving volatility of residuals. I now generalize the solution (ref) to allow for heterogeneous $\sigma_{u_k}^2$ (a diagonal $\Omega_u \neq \sigma_u^2 I_K$) and $\sigma_{\epsilon,t}^2$. The end product is 2SRR, this paper's flagship model.

New matrices must be introduced. First, we have the standard matrix of time-varying residuals variance $\Omega_{\epsilon}=diag( [\sigma_{\epsilon_1}^2 \quad \sigma_{\epsilon_2}^2 \quad ... \quad \sigma_{\epsilon_T}^2]).$ I assume in this section that both $\Omega_{\epsilon}$ and $\Omega_{u}$ are given and will provide a data-driven way to obtain them later. Departing from the homogeneous parameter variances assumption implies that the sole hyperparameter $\lambda$ must now be replaced by an enormous $KT \times KT$ diagonal matrix $\boldsymbol{\Omega_u}=\Omega_u \otimes I_T$ which is fortunately only used for mathematical derivations. For convenience, I split $\boldsymbol{Z}$ in two parts so they can be penalized differently. Hence, the original $\boldsymbol{Z} \equiv [X \quad \boldsymbol{Z_{-0}}]$. The new primal problem is

For convenience, let the $\boldsymbol{\Omega_\theta}$ matrix that stacks on the diagonal all the parameters prior variances, which allow rewriting the problem in a more compact form

Using a GLS re-weighting scheme on observations and regressors, we get a "new" standard primal ridge problem

where $\boldsymbol{\tilde{\theta}}= \boldsymbol{\Omega^{-\frac{1}{2}}_{\theta}}\boldsymbol{{u}}$, $\boldsymbol{\tilde{Z}}=\Omega^{-\frac{1}{2}}_{\epsilon}\boldsymbol{Z} \boldsymbol{\Omega^{\frac{1}{2}}_{\theta}}$ and $\boldsymbol{\tilde{y}}=\Omega^{-\frac{1}{2}}_{\epsilon} \boldsymbol{y}$. Solving this problem by the "dual path" and rewriting it in terms of original matrices gives the general formula

Equation (ref) contains all the relevant information to back out the parameters paths, provided some estimates of matrices $ \boldsymbol{\Omega_{\theta}}$ and $ \Omega_{\epsilon}$.\footnote{This two-step procedure is partly reminiscent of ito2014GLS and ito2017alternative's non-Bayesian Generalized Least Squares (GLS) estimator, where a two-step strategy is also proposed for reasons similar to the above. Their approach could have a ridge regression interpretation with certain tuning parameters fixed. However, the absence of tuning leads to overfitting and the GLS view cannot handle bigger models because the implied matrices sizes are even worse than that of the primal ridge problem discussed earlier.}

Implementation

The solution (ref) takes $ \boldsymbol{\Omega_{\theta}}$ and $ \Omega_{\epsilon}$ as given. In this section, I provide a simple adaptive algorithm to get the heterogeneous variances model estimates empirically. Multi-step approaches to obtain the obtain analogs of $ \boldsymbol{\Omega_{\theta}}$ and $ \Omega_{\epsilon}$ have been proposed in ito2017alternative and giraitis2014inference. Algorithm (ref) follows along and describes a two-step ridge regression (2SRR) which uses a first stage plain RR to gather the necessary hyperparameters.

Here are some necessary details. In step 2, the prior variances for "starting values" $\beta_0$ can either be fixed of cross-validated -- which will be important for bigger models. The prior mean can be zero or the OLS solution as in primiceri2005. Regarding step 3, I use GARCH(1,1) as the volatility model. Better (and faster) alternatives are available for multivariate cases where we may also want to model time-varying covariances. Those are discussed in Section (ref). Lastly, the optimal $\lambda$ may change when $\Omega_{u}$'s entries are heterogeneous, justifying a second CV run in step 4.

2SRR (and eventually MSRR$_{\text{S}}$ in Section (ref)) makes use of adaptive (or data-driven) shrinkage. Adaptive prior tuning has a long tradition in Bayesian hierarchical modeling murphy2012 but the term itself came to be associated with the Adaptive Lasso of zou2006. To modulate the penalty's strength in Lasso, the latter suggest weights based on preliminary OLS (or Ridge) estimates. Those, taken as given, may be contaminated with a considerable amount of noise, especially when the regression problem is high-dimensional (like the one considered here). Adaptive weights in 2SRR have a natural group structure which substantially improves their accuracy through averaging.

Choosing $\lambda$

Derivations from previous sections rely on a given $\lambda$. This section explains how to obtain the amount of time variation by CV (as alluded to in Algorithm (ref)), and how that new strategy compares to more traditional approaches to the problem.

\vskip 0.2cm

{\sc Cross-Validating $\lambda$.} I use k-fold CV for convenience, but anything could be used -- conditional on some amount of thinking about how to make it computationally tractable. This is also what standard RR implementations use, like glmnet in R. A concern is that k-fold CV might be overoptimistic with time series data. Fortunately, bergmeir2018note show that without residual autocorrelation, k-fold CV is consistent. Assuming models under consideration include enough lags of $y_t$, this condition will be satisfied for one-step ahead forecasts. Moreover, GCLSS2018 report that macroeconomic forecasting performance can often be improved by using k-fold CV rather than a CV procedure that mimics the recursive pseudo-out-of-sample experiment.

\vskip 0.2cm

{\sc Time-Variation Budget Interpretation.} Like any ridge regression problem, (ref) can be written as a minimization problem with a "parameter budget" constraint hastie2015statistical. In our context, the parameter budget has a specific meaning as "Time-Variation Budget" given \(\boldsymbol{\theta}\)'s content. Precisely, we have

where TVB stands for the "Time-Variation Budget", a constant corresponding to the maximal total magnitude of \(\boldsymbol{\theta}\). Technically (but less so practically), one could tune TVB instead of $\lambda$ and achieve identical results. Still, this representation is useful to understand how the amount of time-variation is set in the 2SRR context using cross-validation techniques. It asks the question: given a specific dataset finite sample limitations, how much time-variation can we afford? The CV answer is: as much as out-of-sample performance in predicting $y_t$ does not degrade from overly wiggly parameter path estimates.

This representation also helps in understanding the hypothesized behavior of \(\lambda\) as a function of \(T\) and \(K\). In a standard ridge context, where the number of regressors is fixed as the sample size increases, \(\lambda \rightarrow 0\) as \(T\) grows, and we return to the OLS solution. This will not happen in the TVP setup since the effective number of regressors is \(KT\), which grows as fast as \(T\). Regarding \(K\), \(\lambda\) will tend to increase with it, simply because a larger model requires more regularization. This is a feature of the model, not the ridge estimation strategy, and the representation in (ref) clarifies that, for a fixed TVB, the total amount of time-variation must now be shared across more predictors as \(K\) increases. This results in less time-variation allocated to each coefficient individually.

\vskip 0.2cm

{\sc Comparison To Bayesian Approaches.} In the Bayesian TVP-VAR literature, it is common to use a sparse grid of values for the hyperparameters of the smoothing priors, inspired by primiceri2005. Primiceri employed a reversible jump MCMC algorithm where, in an initial step, the model must be estimated conditional on every possible combination of hyperparameters from a grid. However, as noted in amir2018choosing, it is unclear whether this grid is appropriate for applications beyond those originally considered by Primiceri, and simply using wider grids introduces significant computational challenges. Consequently, amir2018choosing and cadonna2020triple propose estimating the hyperparameters in the priors of smoothing parameters through a hierarchical structure within the entire Bayesian procedure, finding that this approach can significantly improve forecasting results.

At a conceptual level, cross-validating \(\lambda\) in the ridge approach aligns with recent Bayesian advances. That is, hyperparameters ultimately influencing the smoothness of coefficients should ideally be obtained in a data-driven manner, and allowing for a reasonably wide range of possible values. However, my approach estimates \(\lambda\) directly (resulting in a single value), whereas amir2018choosing estimate the hyperparameters of the prior for parameters that play an equivalent role to that of \(\lambda\) in their framework.

\vskip 0.2cm

{\sc Comparison to Classical Approaches.} Within the older literature where TVPs were obtained via classical methods, estimating the parameters variances ($\sigma_u^2$ and $\sigma_{\epsilon}^2$ in my notation) made plain MLE's life particularly difficult stock1998median,boivin2005has. One of such issues is the "pile-up" problem, which can occur when a variance parameter is estimated to be 0. In the Ridge Regression (RR) paradigm, where parameter variances are expressed through \(\lambda\), it becomes clear why some of these issues emerged: directly maximizing the in-sample likelihood to determine ridge's \(\lambda\) is not a common practice. Seeing the TVP problem as a high-dimensional regression, one avoids the "pile-up" problem inherent to in-sample maximum likelihood estimation grant2017bayesian by estimating it separately on a hold-out sample via cross-validation.

It is natural to ask how we should think of \(\lambda\) obtained from cross-validation compared to estimating \(\sigma_u\) and \(\sigma_{\epsilon}\) jointly with \(\boldsymbol{\theta}\) using (well-behaved) maximum likelihood estimation. This approach is standard for classical state-space estimation of more compact models, such as local-level models. In that context, for the same ratio of variances, the smoothing splines (or ridge) version of the problem yields the same result as the filtering approach durbin2012time. In the context of HP filtering, which can be seen as a special case of the ridge approach with only an intercept and penalizing double rather than first differences, paige2010ridge also discusses that setting \(\lambda\) to the true ratio of variances \(\frac{\sigma_{\epsilon}^2}{\sigma_{u}^2}\) delivers the best linear unbiased predictor. Possible ways to estimate \(\lambda\) in this context are Generalized CV golub1979gcv or Restricted MLE. However, for TVP models, where \(K\) is often much larger than 1, it is evident that CV may select \(\lambda\) 's yielding less time-variation than MLE would command (thus, not unbiased). This occurs to maximize out-of-sample prediction accuracy and win at the bias-variance trade-off, simply because \(K\) may be large and \(T\) is small.

From Univariate to Multivariate

So far, derivations have focused on the univariate case. This section details the modifications necessary for a multivariate 2SRR and covers de facto the canonical case of a TVP-VAR and also that of estimating local projections over a large horizons range. These will be the two applications considered in Section (ref).

\vskip 0.2cm

{\sc Tackling Many Conditional Means Efficiently .} Since both $\Omega_{u}$ and $\Omega_{\epsilon}$ are equation-specific, we must use ((ref)) for each $\boldsymbol{y}$. However, all estimation procedures proposed in this paper have the homogeneous case of Section (ref) as a first step. This is usually the longer step since it is where cross-validation is done. Hence, it is particularly desirable not to have computations of the first step scaling up in $M$, the number of variables in the multivariate system. Thankfully, in the plain ridge case, we can obtain all parameters of the system in one swift blow, by stacking all $\boldsymbol{y}$'s into $\boldsymbol{Y}$ (a $T \times M$ matrix) and computing

This is precisely the approach that will be used as a first step for any multivariate extension. Of course, this works because the multivariate model has the same regressor matrix for each equation (like VARs and LPs). In this homogeneous variance case, $\boldsymbol{Z'}(\boldsymbol{ZZ'}+\lambda\boldsymbol{I}_T)^{-1}$ is the same for all equations and cross-validation still implies inverting $(\boldsymbol{ZZ'}+\lambda\boldsymbol{I}_T)$ as many times as we have candidates for $\lambda$. That is, even if we wish to have a different $\lambda$ for each equation in the first step, computing time does not increase in $M$, except for matrix multiplication operations which are much less demanding.\footnote{Precisely, cross-validation implies calculating \# of folds $\times$ \# of $\lambda$'s the $\boldsymbol{P}^\lambda_Z$. Then, these matrices can be used for the tuning of every $m$ equation, which is precisely why the computational burden only very mildly increases in $M$.}

\vskip 0.2cm

{\sc Time-Varying Covariance Matrix.} When entering multivariate territory, modeling the residuals covariances -- a necessary input to structural TVP-VAR analysis (but not point forecasting) -- arises as an additional task. The number of TVPs entering the $\boldsymbol{\Omega_\epsilon}$ matrix grows quickly with $M$ as there are $\frac{M(1+M)}{2}$ of them. This is the other key part of the challenges facing traditional Bayesian TVP-VARs computations. So far, we have seen how the ridge approach can help in computing the various conditional means efficiently. The following discusses how it can also be applied to the time-varying covariance matrix.

I now describe how similar ideas as above -- applying the same smoothing matrix to many series at once -- can also be applied so that each parameter of the covariance matrix follows a random walk. First, let us denote $\boldsymbol{\epsilon_t} \in {\rm I\!R}^M$ as the residuals obtained from running multivariate 2SRR with $\boldsymbol{Y}$ as dependent variables. Let $\tilde{\eta}_t=vech(\boldsymbol{\epsilon_t}\boldsymbol{\epsilon_t}')$ where the $vech$ operator means that we vectorize but only keep lower-diagonal elements (i.e., only the non-redundant entries). Also, denote $\bar{\eta}_t$ to be a vector corresponding to the full-sample average of the columns of $\tilde{\eta}_t$, or, in other words, a "time series" of the lower-diagonal elements of what would be the time-invariant covariance matrix for those residuals. To be more concrete, in the $M=2$ case, we have $\tilde{\eta}_t = [ \epsilon_{1,t}^2 \enskip \epsilon_{2,t}^2 \enskip \epsilon_{1,t}\epsilon_{2,t} ]$ and $\bar{\eta}_t = [ \frac{1}{T}\sum_{\tau=1}^{T}\epsilon_{1,\tau}^2 \enskip \frac{1}{T}\sum_{\tau=1}^{T} \epsilon_{2,\tau}^2 \enskip \frac{1}{T}\sum_{\tau=1}^{T} \epsilon_{1,\tau}\epsilon_{2,\tau} ]$. Finally, let $\tilde{\boldsymbol{\eta}}$ be the $T \times \frac{M(1+M)}{2}$ matrix stacking those $\frac{M(1+M)}{2}$ time series of residuals products together and $\bar{\boldsymbol{\eta}}$ doing an analogous stacking for $\bar{\eta}_t $. Treating those as observed data, we can get the whole set of paths $\hat{\boldsymbol{\eta}}$ by writing a multivariate TVP-ridge problem as the above, but with $\frac{M(1+M)}{2}$ targets and a conditional mean model that is only a time-varying intercept. Precisely, we can compute

where $\varphi$ is a smoothness hyperparameter (just like $\lambda$) and $\boldsymbol{D}$ is the matrix difference operator described in Section (ref). This implies, again, that each column of the $\tilde{\boldsymbol{\eta}}$ matrix is assumed to follow a random walk around its mean value with a velocity inversely linked to $\varphi$. Clearly, large $\varphi$'s push the solution towards a time-invariant covariance matrix of residuals ($\bar{\boldsymbol{\eta}}$) and extremely low values will make it overly wiggly.\footnote{I also note that overly small values of $\varphi$ could threaten the positive-definiteness (and thus validity) of the resulting matrix because, intuitively, some minimal form of averaging is necessary to obtain covariances. This problem has not been encountered in practice because results become implausibly wiggly well before positive-definiteness breaks down. This also highlights a homogeneous $\varphi$ is necessary here, elements of the same matrix averaged differently at different points in time could create momentousness incoherency problems when putting things back together. } Note that without the pre-demeaning of the $\tilde{\eta}_t$ series (which the use of $\bar{\boldsymbol{\eta}}$ is the formal mathematical characterization of) and adding back the mean is only so that shrinkage pushes towards the solution of time-invariant covariance matrix rather than 0. Computations-wise, for a fixed $\varphi$, this amounts to one matrix $T \times T$ inversion and then multiplying it with $\tilde{\boldsymbol{\eta}}$, with both operations fast for $T$'s encountered in macroeconomics.

In applications, I opt for a user-specified $\varphi$. Unlike the conditional mean targets, covariances in such applications are fundamentally unobserved, and the elements of $\tilde{\boldsymbol{\eta}}$ in themselves cannot be regarded as realized variances or covariances. Hence, it is not entirely clear that the target one would be cross-validating for is a proper target at all. With some ingenuity (left for future work\footnote{One possibility would be to cross-validate $\lambda$ and $\varphi$ together and use multivariate out-of-fold log scores as metric.}) on how to create a proper scoring metric in this case, CV could be conducted very fast because it only implies inverting as many $T \times T$ matrices as one has $\varphi$'s, which by all reasonable means, should be below 50.

The algorithm proposed in this subsection provides a covariance matrix at each $t$. This sequence of matrices can then be used to conduct structural TVP-VAR analysis like in Section (ref) where a Cholesky decomposition is applied at each $t$. Evidently, this is simplest thing one can do, but as long as we have the sequence of covariance matrices, there are quite a few more sophisticated things one could opt for.

Finally, two important notes for the practical use of the tools above. First, it is best practice to add a sufficient number dummy observations duplicating the boundary value at the beginning and end of $\tilde{\boldsymbol{\eta}}$.\footnote{I include 30 such duplicates in the applications.} This is to avoid that the $\boldsymbol{D'D}$ matrix, in absentia of observations before 1 and after $T$, overly shrinks $(\tilde{\boldsymbol{\eta}}-\bar{\boldsymbol{\eta}})$ to 0, and thus in our case, mechanically forcing the top and bottom rows of $\hat{{\boldsymbol{\boldsymbol{\eta}}}}$ towards $\bar{\boldsymbol{\eta}}$. This remedy is obviously very much related to a similar tactic employed to deal with boundary problems in kernel regressions.

Second, and far more subtle, is the case of the "ridgeless" solution. For very large $\frac{K}{T}$ ratios, it has been observed that cross-validation will sometimes pick $\lambda \approx 0$, which is the so-called ridgeless regression of hastie2019surprises. This is a reflection of the double descent phenomenon in extremely high-dimensional models, which is ubiquitous in modern deep learning belkin2019reconciling,bartlett2020benign. In a nutshell, the unregularized model might provide a more regularized (and better) solution out-of-sample since, with $\lambda \rightarrow 0$, the solution converges to the minimum-norm solution because the inverse of the covariance matrix is Moore-Penrose pseudo-inverse (see kellyvirtue for a more detailed discussion applied to finance). However, this "benign overfitting" solution (because forecasting is not impacted) provides perfect fit in-sample. In that scenario, there are no residuals, no covariance matrix, and $\beta_t$ certainly overfit, rendering historical analysis through IRFs unfeasible. Fortunately, there is an easy mechanical fix, which is limiting the range of $\lambda$ to solutions that exists the classical bias-variance regime, rather than the so-called double descent regime. In many applications, it is found that there is a solution in the "usable" regime that is almost as good (in terms of cross-validated MSE) as the ridgeless solution hastie2019surprises.

Quantifying Uncertainty

Clearly, 2SRR advantage is in getting point estimates rapidly, which already goes a long way for many of TVP models applications. Yet, there are many others where quantifying $\boldsymbol{\beta}$'s uncertainty is instrumental. In principle, this would be possible for 2SRR by leveraging the link between ridge and a plain Bayesian regression murphy2012. In the general case, one would need to obtain

which is precisely the large matrix we were avoiding to invert earlier. While this is not an issue for smaller models, it can become one in medium-sized models because it relies on the primal solution. Also, note that in the simple case where $ \Omega_{u}=\sigma_u^2 I_K$ and $ \Omega_{\epsilon}=\sigma_\epsilon^2 \ I_T$, there is a clear Bayesian interpretation allowing the use of the posterior variance formula for linear Bayesian regression. However, it treats the cross-validated $\lambda$ as known. In a similar line of thought, it treats the hyperparameters inherent to 2SRR as given when computing the bands. While in traditional (modestly dimensional) applications of ridge regression, this could be regarded as a minor impediment, it is not the case here because, in very high-dimensional models, a small change in regularization structure can lead to rather different coefficients estimates. Therefore, it is important to account not only for uncertainty conditional on hyperparameters, but also that which emanate directly from them.

\vskip 0.2cm

{\sc Weighted Bayesian Bootstrap.} I provide a way for users to quantify TVP uncertainty that circumvents both of the aforementioned issues, namely, it will not need to invert a $KT \times KT$ matrix and will account for regularization uncertainty. I adapt the Weighted Bayesian Bootstrap (WBB) of newton2018weighted, which is an extension of newton1994approximate's ideas to general machine learning models with regularization.\footnote{More closely related to this paper's framework, newton2018weighted consider a "trend filtering" example, which is a TVP model with only a intercept where the second difference is penalized with a $l_1$ norm. } In short, randomization via exponential draws is repeatedly applied to observations and regularization weights, thereby accounting for $\lambda$ and overall regularization uncertainty. This approximate Bayesian inference technique is increasingly popular for scalable inference with machine learning methods or Bayesian models with demanding posterior computations nie2022bayesian. Clearly, there are many values of $K$ and $T$ for which (ref) qualifies as such.

Implementation is convenient, quick, and is easily parallelizable if need be. It implies looping over ridge estimation (without CV and only using the dual solution (ref)), which is rapid and amenable also to extensions discussed in the following section. As an introductory example, applying WBB in its plainest form to the simplest TVP formulation found in equation (ref) implies repeatedly estimating

with $\omega_t \sim \text{Exp}(1)$, $\omega_{\lambda} \sim \text{Exp}(1)$, and where $\lambda_{\text{CV}}$ stands for $\lambda$ as estimated from running cross-validation on the unweighted problem. As discussed in ng2022random, it may be worthwhile in high-dimensional models to have a distinct $\omega_{\lambda}$ for the regularization strength of each predictor. This therefore account not only for overall regularization strength uncertainty, but that of relative regularization, i.e., how heavily certain regressors' coefficients are shrunk versus others. Clearly, this is relevant in the 2SRR context because it explicitly contains regularization heterogeneity via $\Omega_{u}$ and $\Omega_{\epsilon}$ estimated with from a homogeneous RR first step. Therefore, what we want is something along the lines of a $\omega_{\lambda,k}$. Accordingly, the WBB procedure proposed to accompany 2SRR (and tested via simulations in Section (ref)) is to divide each diagonal element of the two variance matrices (${\Omega_{\theta}}$ and $\Omega_{\epsilon}$ in equation (ref)) by a $\text{Exp}(1)$ draw, and give $\lambda_0$ its own $\omega_0$.\footnote{In practice, some trimming of the draws is necessary to avoid anything close to dividing by 0, and thus multiplying variances by very large values.} Then, one can use the draws

where $\boldsymbol{\Omega_{\theta}^b}$ and $\Omega_{\epsilon}^{\boldsymbol{b}}$ denote the modulated matrices to calculate percentiles and other relevant metrics. Note that by the virtue of this being an approximate Bayesian method, it provides no frequentist coverage guarantees (along with most Bayesian methods). Nonetheless, to get a sense how credible regions constructed from a WBB approach behave in comparison a full-fleshed Bayesian procedure, I investigate nominal coverage in Section (ref) and find they are comparable.\footnote{One could further refine the above by considering a blocked Bayesian Bootstrap approach for the elements of $\Omega_{\epsilon}^{\boldsymbol{b}}$ (as considered for other kinds of ML algorithms in MRFjae and HNN), but in the conducted experiments (which includes persistent targets and regressors), this did not seem to play an appreciable role. Such a modification is easily implementable by drawing exponentials in groups. } {Note that in the context of VAR modeling, we are typically interested in making inference on impulse responses, which are nonlinear transformation of $\boldsymbol{\beta}$. In that scenario, $\hat{\boldsymbol{\theta}}_{\boldsymbol{b}}$ draws can be transformed in IRFs draws for each ${\boldsymbol{b}}$.

\vskip 0.2cm

{\sc Other Avenues for Density Forecasting.} Despite not being exploited in this paper, the ridge and TVP connection also opens the door for the usage of conformal prediction methods that are increasingly popular in the ML literature to obtain model free prediction intervals with coverage guarantees. WBB draws could be used as original inputs and time series adaptation of these methods could be utilized chernozhukov2018exact.

Refining the Conditional Mean: Two Extensions

This subsection extends 2-step ridge regression to a multistep ridge regression (MSRR) which, as the name suggests, by iterating further can deliver two refinements of the conditional mean function. The first will deliver Sparse TVPs (MSRR$_{\text{S}}$) and the second Dense TVPs (MSRR$_{\text{D}}$) through a factor structure.

Sparse TVPs via Adaptive Ridge

I provide a way to iterate it so that not only it fine tunes $\Omega_u$ but also set some of its elements to zero, obtaining Sparse TVPs. That is, some parameters will vary and some will not, which can provide efficiency gains. Those have already been proposed in the standard Bayesian MCMC paradigm most notably by bitto2018achieving and belmonte2014hierarchical, so it worthwhile to demonstrate how a similarly sparse behavior can be obtained from a relatively minor twist on 2SRR. The new primal problem is

which is just adding the penalty $\xi tr(\Omega_{u}^{\frac{1}{2}})$ to ((ref)). In the RR paradigm, it is quite straightforward to implement since making a parameter constant amounts to dropping the group of regressors $\boldsymbol{Z}_k$ corresponding to the basis expansion of $X_k$ making it time-varying. Thereby, it corresponds to some blend of ridge regression and Group Lasso. This is done by setting this group's $\sigma^2_{u,k}$ to 0. The proposed implementation, formalized by Algorithm (ref) (Appendix (ref)), amounts to iterating ridge regressions (hence the name) and constantly updating the elements of $\Omega_u$ in a particular way. This estimation approach is inspired by grandvalet1998's proposition of using Adaptive Ridge to compute the Lasso solution. The insight has since been recuperated by frommletnuel2016 and liuli2014 to implement $l_0$ regularization in a way that makes computations tractable. MSRR$_{\text{S}}$ goes back to implementing the $l_1$ norm by Adaptive Ridge as in grandvalet1998 but extend it to do Group Lasso and add a Ridge penalty within selected groups. Derivation details are omitted from the main text and can be found in Appendix (ref).

Dense TVPs with Reduced Rank Restrictions

A frequent empirical observation, dating back to cogley2005drifting, is that $\beta_t$'s can be spanned very well by a handful of latent factors. dewindgambetti2014, stevanovic2016common and chaneisenstat2018 exploit this that directly by implementing directly a factor structure within the model. It is clear that dimensionality can be greatly reduced if we only track a few latent states and impose that evolving parameters are linear combinations of those -- Dense TVPs. Precisely, the model under consideration here is univariate but reduced rank restrictions will be applied to a matrix $\boldsymbol{U}=vec^{-1}_{K \times T}(\boldsymbol{u})$ where $\boldsymbol{u}$ are the coefficients from an univariate ridge regression the likes of which we have been running all along.\footnote{This can be done because $\boldsymbol{u}$ has an obvious block structure. It has two dimensions, $K$ and $T$, that we can use to create a matrix. } The primal problem from TVP model with a dense structure for TVPs and sparse loadings can be written as

where $\boldsymbol{f}=vec(\boldsymbol{F})$, $l=vec(\Lambda)$, $\Lambda$ being $K \times r$, and $F$ being $r \times T$.\footnote{To make the exposition less heavy, I assume throughout this section that $\Omega_{\epsilon}$ is given and that $\beta_0$ are not penalized in any way.} Appendix (ref) goes through the details of how to estimate this model with MSRR$_{\text{D}}$, which, by using two vectorization identities, splits the optimization into a ridge step that estimates "factors" (the same way we typically estimate TVPs in this paper, but now with a smaller number of them ($r<K$)) and a Lasso step that updates loadings. MSRR$_{\text{D}}$ estimates are obtained from the alternation of those two steps.

Simulations

The simulation study investigates how accurately the different estimators proposed in this paper can recover the true parameters path. I consider three numbers of observations $T\in\{150,300,600\}$. Most of the attention will be dedicated to $T=300$ since it is roughly the number of US quarterly observations we will have 15 years from now. The size of the original regressor matrix $X$ is $K \in \{6,20,100\}$ and the first regressor in each is the first lag of $y$. Figure (ref) display the 5 types of parameters path $f_i$ that will serve as basic material: cosine, quadratic trend, discrete break, a pure random walk and a linear trend with a break. $f_1$, $f_2$ and $f_4$ "fit" relatively well with the prior that coefficients evolve smoothly whereas $f_3$ and $f_5$ can pose more difficulties. In those latter situations, TVP models are expected to underperform.\footnote{This partly motivates the creation of Generalized TVPs via Random Forest in MRFjae.} The design for simulations $S_1$, $S_2$, $S_3$ and $S_4$ can be summarized in a less cryptic fashion as

The considered proportions of TVPs within the $K$ parameters are $\sfrac{K^*}{K} \in \{0.2,0.5,1\}$. Formally, we have

The scale of coefficients is manually adjusted to prevent explosive behavior and/or overwhelmingly high $R^2$'s. The most important transformation in that regard is a min-max normalization on the coefficient of $y_{t-1}$ to prevent unit/explosive roots or simply persistence levels that would drive the true $R^2$ above its targeted range. Regarding the latter, I consider four different types of noise process. Three of them are homoscedastic and have a \{Low, Medium, High\} noise level. Those are calibrated so that $R^2$'s are around 0.8, 0.5 and 0.3 for low, medium and high respectively. The last two noise processes are SV, which is the predominant departure from the normality of $\epsilon_t$ in applied macroeconomics. For better comparison with time-invariant volatility cases, those are "manually" forced (by a min-max normalization) to oscillate between a predetermined minimum and maximum. The first SV process is constrained within the Low and Medium noise level bounds. For the second, it is Low and High, making the volatility spread much higher than in the first SV process case.

Four estimators are considered: the standard TVP-BVAR with SV\footnote{For the TVP-BVAR, I use the R package by Fabian Krueger that implements primiceri2005's procedure (with the delnegroprimiceri2015 correction), available \href{https://sites.google.com/site/fk83research/code}{here}. Default parameters are used. The total number of MCMC iterations is 15\,000 with burn-in of 5\,000. }, the two-step Ridge Regression (2SRR), the sparse version (MSRR$_{\text{S}}$) and the dense one (MSRR$_{\text{D}}$).\footnote{The maximal number of factors for MSRR$_{\text{D}}$ is set to 5 and the chosen number of factors is updated adaptively in the EM procedure according to a share of variance criteria.} TVP-BVAR results are only obtained for $K=6$ for obvious computational reasons. Performance is assessed using the mean absolute error (MAE) with respect to the true path. I then take the mean across 100 simulations for each setup. To make this multidimensional notation more compact, let us define the permutation $\mathcal{J}=\{K,\sfrac{K^*}{K},\sigma_{\epsilon},S_i \}$. I consider simulations $s=1,...,50$ for all $\mathcal{J}$'s. Formally, for model $m$ and setup $\mathcal{J}$, the reported performance metric is $\frac{1}{50}\sum_{s=1}^{50}MAE_\mathcal{J}^{s,m}$ where

Results

The results for $T=300$ are in (ref). With these simulations, I am interested in verifying many things, which are covered aspect by aspect in paragraphs that follow. First, I want to verify that 2SRR's performance is comparable to that of the BVAR for models' size that can be handled by the latter. Second, I want to demonstrate that additional shrinkage embedded in MSRR$_{\text{S}}$ and MSRR$_{\text{D}}$ can help under DGPs that more or less fit the prior of reduced-rank and/or sparsity. To make the investigation of these two points visually easier by looking at the tables, the lowest MAE out of BVAR/2SRR for each setup is in blue while that of the best one out of all algorithms is in bold. Then, follows a series of investigations, like the effects of changing sample size, time-varying volatility in DGP, benchmarking in higher dimensions, computational time, and a look at the empirical properties of the WBB procedure for uncertainty quantification.

\vskip 0.2cm

{\sc 2SRR vs. BVAR.} Overall, results for 2SRR and the BVAR are very similar and their relative performance depends on the specific setup. These two models are interesting to compare because they share the same prior for TVPs (no additional shrinkage) but address evolving residuals volatility differently. Namely, the BVAR models SV directly within the MCMC procedure whereas 2SRR is a two-step GLS-like approach using a GARCH(1,1) model of the first step's residuals. Figure (ref) summarizes results of the 2SRR/BVAR comparison by reporting boxplots showcasing the distribution of relative MAEs (2SRR/BVAR) for different subsets. Overall, 2SRR does marginally better in almost all cases. A lower noise level seems to help its cause. It is plausible that cross-validating $\lambda$ as implemented by 2SRR plays a role (BVAR uses default values).\footnote{Replacing the absolute distance by the squared distance in (ref) produces similar looking boxplots as in Figure (ref), with wider dispersion but a near-identical ranking of methods by simulations and noise processes.}

I now dive in specific DGPs. In Table (ref), where the DGP is the rather friendly cosine-based TVPs, it is observed that the BVAR will usually outperform 2SRR by a thin margin when the level of noise is high. The reverse is observed for low noise environment and results are mixed for the medium one. Results for the SV cases will the subject of its own discussion later. In Table (ref), where the DGP is a structural break -- at odds with the prior of slow change -- 2SRR is clearly performing better than the BVAR, providing a smaller average MAE in 12 out of 15 cases for $K=6$. In Table (ref), where the DGP is a mix of trending coefficients and cosine ones, Table (ref) reports very similar results to that of Simulation 1. 2SRR is better than BVAR except in the high noise setups, where the latter has a minor advantage. Lastly, for the simulation in Table (ref), a sophisticated mixture of TVP-friendly and -unfriendly parameters paths, BVAR does a better job than 2SRR for 8 out of 15 cases. The gains are, as before, quantitatively small. When 2SRR does better, gains are also negligible, suggesting that BVAR and 2SRR provide very similar results in this environment.

\vskip 0.2cm

{\sc MSRR$_{\text{S}}$ and MSRR$_{\text{D}}$.} Overall, the two refinements can help further improving results, and unsurprisingly, they tend to do so where the underlying DGP best fits their prior. In Table (ref), for $K=6$, MSRR$_{\text{S}}$ will marginally improve on 2SRR for most setups, especially those where 2SRR is already better than BVAR. In higher dimensions ($K=20$ or $K=100$), MSRR$_{\text{S}}$ constantly improves on 2SRR (albeit minimally) whereas MSRR$_{\text{D}}$ can provide important gains (see the $\sfrac{K^*}{K}=1$ block for instance) but is more vulnerable in the high noise environment. Moving to Table (ref), for the small dimensional case, it is observed that MSRR$_{\text{S}}$ can further reduce the MAE --- albeit by a very small amount --- in many instances. The same is true for MSRR$_{\text{D}}$ when all parameters vary ($\sfrac{K^*}{K}=1$). For MSRR$_{\text{D}}$, this observation additionally extends to $K=20$, an environment where it is expected to thrive. In Table (ref), MSRR$_{\text{S}}$ often marginally improves upon 2SRR whereas MSRR$_{\text{D}}$'s edge is more visible in low-noise and high-dimensional environments --- factors being more precisely estimated with a large cross-section. Lastly, in Table (ref), when it comes to higher-dimensional setups ($K=20$ or $K=100$), MSRR$_{\text{S}}$ emerges as the better option with (now familiar) marginal improvements with respect to 2SRR. This recurrent observation is potentially due to the iterative process producing a more precise $\hat{\Omega}_u$ when $\sigma_{u_k}^2$'s are heterogeneous whether sparsity is involved or not.

\vskip 0.2cm

{\sc Changing Sample Size.} The results for $T=150$ and $T=600$ are in (ref). When $T$ is reduced from 300 to 150, the performance of 2SRR relative to that of BVAR remains largely unchanged: both report very similar results. When bumping $T$ to 600, overall performance of all estimators improves, but not by a gigantic leap. This is, of course, due to the fact that increasing $T$ also brings up the number of effective regressors. BVAR has a small edge on Simulation 1 in Table (ref) whereas 2SRR wins marginally for the more complicated Simulation 4 (Table (ref)).

What is perhaps most noticeable from simulations with a larger $T$ is how much more frequently MSRR$_{\text{S}}$ and especially MSRR$_{\text{D}}$ are preferred, especially in the medium- and high-dimensional cases. For instance, for the Cosine DGP ($S_1$) with $K\in\{20,100\}$, MSRR$_{\text{D}}$ almost always deliver the lowest MAE, and sometimes by good margins (e.g., $\{S_1,\sfrac{K^*}{K}=1,\sigma_{\epsilon}=\text{Low} \}$ for both $K$'s). Similar behavior is observed for $S_3$ in almost all cases of $K=20$. This noteworthy amelioration of MSRR$_{\text{D}}$ is intuitively attributable to factor loadings being more precisely estimated with a growing $T$. Thus, unlike 2SRR whose performance is largely invariant to $T$ by model design, algorithms incorporating more sophisticated shrinkage schemes may benefit from larger samples.

\vskip 0.2cm

{\sc Time-Varying Volatility.} A pattern emerges across the four simulations: when SV is built in the DGP ($\sigma_{\epsilon,t}$ in tables), 2SRR either performs better or deliver roughly equivalent results to that of the BVAR. Indeed, with $T=300$, for 17 out of 24 setups with SV-infused DGPs, 2SRR supplants BVAR. The wedge is sometimes small ($\{S_1,\sfrac{K^*}{K}=0.2,\text{SV Low-Med} \}$,$\{S_4,\sfrac{K^*}{K}=1,\text{both SV} \}$), sometimes large ($\{S_1,\sfrac{K^*}{K}=1, \text{SV Low-High} \}$,$\{S_2,\sfrac{K^*}{K}=0.5, \text{both SV}\}$). However, it is fair to say that small gaps between 2SRR and BVAR performances are the norm rather than the exception. Nonetheless, these results suggest that 2SRR is not merely a suboptimal approximation to the BVAR in the wake of computational adversity. It is a viable statistical alternative with the additional benefit of being easy to compute and to tune.

\vskip 0.2cm

{\sc Benchmarking in Higher Dimensions.} In the years following the original draft of this paper, new Bayesian procedures have been proposed and packaged in statistical software. This allows for an easier evaluation of 2SRR in higher dimensions, both statistically and computationally. One such package is shrinkTVP from knaus2021shrinkage. It conveniently features C++ based computations, a very wide set of potential prior specifications, focus on single equation modeling, and features a data-driven way of choosing shrinkage intensity. Two default specifications are used, one with a ridge prior (ShTVP-R) on coefficients and another with the triple gamma prior (ShTVP-3G) of cadonna2020triple). The former is the closest relative to 2SRR within the suite of available configurations. The latter should perform well in the sparse TVP environments for which it has been designed. The high-dimensional case of $K=100$ has been reduced to $K=50$ since it took shrinkTVP 38 minutes to estimate a single model in the original large $K$ case (whereas 2SRR delivers point estimates in under 30 seconds).

Results in Table (ref) further establish the relevance of 2SRR for the estimation of models with many TVPs. Its statistical performance is mostly on par with ShTVP-R, with differences always being marginal in favor of one or another. ShTVP-3G can provide some less trivial improvements over both 2SRR andShTVP-R under a sparse DGP or in the hostile high noise and $K=50$ environment. Otherwise, it does not meaningfully outperform 2SRR. This is evidently the prior's doing, and it is an empirical question as to which setup in Table (ref) most closely resembles real data.

\vskip 0.2cm

{\sc Computational Time.} In Table (ref), we see that the Bayesian computations are manageable for $K=20$, and becomes rather burdensome for $K=50$ (about 10 minutes vs. less than 15 seconds). Even at $K=20$, the fact that 2SRR gets similarly valid estimates in 5 seconds rather than in over a minute makes it more convenient in applied work where the trial and error of specifications (and sometimes even grid search of them) is widespread. Of course, full-blown Bayesian estimation directly provides many additional quantities that 2SRR does not procure readily, and its computational time could be brought down to a certain extent by cutting MCMC iterations. Nonetheless, the point is that, depending on the intended use, the ridge-based methods have non-trivial comparative advantages in many ways that are relevant for applied work.

It is apparent that for large models, such as those used in rolling/expanding window forecasting evaluations (Section (ref)) and large local projections (Section (ref)), 2SRR remains the only practical approach. As discussed in Section (ref), we expect computations for 2SRR (for a single equation) to scale (approximately) linearly with \( K \) for a fixed \( T \). Table (ref) (Appendix) provides additional results, verifying this by showing how 2SRR computational time varies with \( K \) and \( T \) for the Mixture DGP. I also find, as expected, that an increasing sample size has a nonlinear effect on computing time. For instance, the most demanding combination (\( K = 100 \), \( T = 600 \)) takes about 153 seconds to compute, whereas the same high-dimensional \( K = 100 \) results in 3 seconds when \( T = 150 \) and 22 seconds when \( T = 300 \).

The table also reports running times for estimating \(K\) equations, rather than just one, to illustrate how the techniques discussed in Section (ref) help in preventing the running time from being multiplied by \(K\) for large systems. Estimating 50 equations, each with 50 time-varying parameters (TVPs), resulting in a total of 2,500 TVPs, takes one minute with \(T = 300\). This is about six times longer than the estimation time for a single equation (\(K = 50\), \(T = 300\)). The most extreme case tested involves estimating 100 equations, each with 100 parameters (totaling 10,000 TVPs) over 600 time periods, which takes 30 minutes to run.

\vskip 0.2cm

{\sc Some Evaluation of Uncertainty Quantification for 2SRR.} How reliable are credible regions obtained from the WBB approach proposed in Section (ref) versus that of fully Bayesian methods deployed above? I provide some answers in Table (ref) (Appendix) where I report nominal coverage for 2SRR (obtained from 250 WBB draws) and those from ShTVP-R and ShTVP-3G. As mentioned earlier, Bayesian techniques (approximate or not) are not formally expected to hit nominal coverage targets, an inherently frequentist object. Nonetheless, looking at such metrics gives an idea of how those methods will behave when using them on real data, e.g., how conservative or not they may be. For instance, it could be worrisome if 2SRR is always much higher or lower than either ShTVP-R or ShTVP-3G.

In all cases, I use the relevant percentiles from the draws coming out of the procedure and calculate the frequency (over $t$ and 20 simulations) at which the true $\beta_{t,k}$ falls in the interval. Also, to give a sense how methods adapt the local width of bands along with heteroscedasticity, in the two DGPs where errors follow a SV process, I calculate the average correlation between the empirical 12% and 84% interquantile range (IQR) of $\beta_{t,k}$ (averaging over $k$) and the true SV process. Since we expect heteroscedasticity to drive time-varying width of bands, those two quantities should be correlated in a non-trivial fashion.

Again, in Table (ref), I focus on the $T=300$ case, and look at a subset of simulations for this particular aspect. Namely, I center the attention on the mixture DGP of Table (ref) above and look at $K=6$ and $K=20$ in the middle ground case of ${\sfrac{K^*}{K}=0.5}$. Except in a few instances, 2SRR provides coverage that sits somewhere between ShTVP-R and ShTVP-3G. ShTVP-R is conservative, and so is for the most part 2SRR. ShTVP-3G is almost always under the other two, and sometimes by a lot (like for the 68% high noise case). 2SRR often remains in the middle, but more on the side of ShTVP-R, especially for $K=20$. The discrepancies are less substantial for the 95% target at $K=6$, with 2SRR oscillating around 95%, ShTVP-R being above by 2-3%, and \texttt{ShTVP}-3G remaining below. While 2SRR and \texttt{ShTVP}-R coverage is mostly unchanged moving from $K=6$ to $K=20$ and remains above nominal targets, \texttt{ShTVP}-3G often provides sharper posteriors than what frequentist coverage would commend.

Perhaps the most intriguing observation from Table (ref) is how the relative width of credible regions around TVPs correlates with the true SV process. In all cases, 2SRR provides bands which IQR is most correlated with the true time-varying volatility. Its correspondence with it is up to twice as high vis-à-vis ShTVP-R and even more so when compared to ShTVP-3G. Thus, their ability to reflect changing levels of uncertainty in the width of bands, is at the very least, within the simulation setup under study, on par with two Bayesian methods that explicitly feature a SV component. Some additional intuition on this -- beyond the fact that 2SRR models volatility in its own way between the first and the second step -- can be obtained from noting that in the frequentist domain, the analogous nonparametric "pairs" bootstrap is valid under heteroscedasticity mackinnon2006bootstrap. In fact, from a Bayesian point of view, lancaster2003 show that the obtained variance for a linear regression model from using such a bootstrap is asymptotically equivalent to what one would get from White's sandwich formula.

Forecasting

Machine Learning models' merits are typically evaluated out-of-sample before taking to interpret them. 2SRR lends itself nicely to forecasting because everything is automatically tuned and more importantly, its quick computations makes it amenable to recursive forecasting backtests where the models needs to be re-estimated at every step. This section describes forecasting results and discusses several practical aspects (like variations on CV) that may help improve forecasting performance on real data and further inform choices for macroeconomic applications in Section (ref).

I present results for a pseudo-out-of-sample forecasting experiment at the quarterly frequency using the dataset FRED-QD mccracken2020fred. The latter is publicly available at the Federal Reserve of St-Louis's web site and contains 248 US macroeconomic and financial aggregates observed from 1960Q1. The forecasting targets are real GDP, Unemployment Rate (UR), CPI Inflation (INF), 1-Year Treasury Constant Maturity Rate (IR) and the difference between 10-year Treasury Constant Maturity rate and Federal funds rate (SPREAD). These series are representative macroeconomic indicators of the US economy which is based on GCLSS2018's exercise for many ML models, itself based on KLS2019 and a whole literature of extensive horse races in the spirit of SW1998comparison. The series transformations to induce stationarity for predictors are indicated in mccracken2020fred. For forecasting targets, GDP, CPI and UR are considered $I(1)$ and are first-differenced. For the first two, the natural logarithm is applied before differencing. IR and SPREAD are kept in "levels". Forecasting horizons are 1, 2, and 4 quarters. For variables in first differences (GDP, UR and CPI), average growth rates are targeted for horizons 2 and 4.

The pseudo-out-of-sample period starts in 2003Q1 and ends 2014Q4. I use expanding window estimation from 1961Q3. Models are estimated and tuned at each step. I use direct forecasts, meaning that $\hat{y}_{t+h}$ is obtained by fitting the model directly to $y_{t+h}$ rather than iterating one-step ahead forecasts. Following standard practice in the literature, I evaluate the quality of point forecasts using the root Mean Square Prediction Error (MSPE). For the out-of-sample (OOS) forecasted values at time $t$ of variable $v$ made $h$ steps ahead, I compute $$RMSPE_{v,h,m}= \sqrt{ \frac{1}{\#\text{OOS}}\sum_{t \in \text{OOS}} (y_{t}^v-\hat{y}_{t-h}^{v,h,m})^2}.$$ The standard dieboldmariano (DM) test procedure is used to compare the predictive accuracy of each model against the reference AR(2) model. $RMSPE$ is the most natural loss function given that all models are trained to minimize the squared loss in-sample.

Three types of TVPs will be implemented: 2SRR (Section (ref)), MSRR$_{\text{S}}$ (Section (ref)), MSRR$_{\text{D}}$ (Section (ref)). I consider augmenting four standard models with different methodologies proposed in this paper. The first will be an AR with 2 lags. The second is the well-known sw2002 ARDI (Autoregressive Diffusion Index) with 2 factors and 2 lags for both the dependent variable and the factors. The third is a VAR(5) with 2 lags and the system is composed of the 5 forecasted series. Finally, I consider as a fourth model a VAR(20) with 2 lags in the spirit of banbura2010large's medium VAR. Thus, there is a total of $4 \times 4 = 16$ models considered in the exercise. The BVAR used in Section (ref) is left out for computational reasons --- models must be re-estimated every quarter for each target. Moreover, the focus of this section is rather single equation direct (as opposed to iterated) forecasting.

The first three constant coefficients models are estimated by OLS, which is standard practice. Since the constant parameters VAR(20) has 41 coefficients and around 200 observations, it is estimated with a ridge regression. Potential outliers are dealt with as in GCLSS2018 for Machine Learning models. If the forecasted values are outside of $[\boldsymbol{\bar{y}}+2*min(\boldsymbol{y}-\boldsymbol{\bar{y}}),\boldsymbol{\bar{y}}+2*max(\boldsymbol{y}-\boldsymbol{\bar{y}})]$, the forecast is discarded in favor of the constant parameters forecast.

Results

I report three sets of results. Table (ref) corresponds exactly to what has been described beforehand. Table (ref) gathers results where TVPs have been additionally shrunk to their constant parameters counterparts by means of model averaging with equal weights. The virtues of this Half & Half strategy are two-fold. First, k-fold CV can be over-optimistic for horizons $h$>1 because of imminent serial correlation. Second, k-fold CV ranks potential $\lambda$'s using the whole sample, whereas in the case of "forecasting", prediction always occurs at the boundary of the implicit kernel. In that region, the variance is mechanically higher and ensuing predictions could benefit from extra shrinkage. Shrinking to OLS in this crude and transparent fashion is a natural way to attempt getting even better forecasts. Finally, Table (ref) reports results where k-fold CV has been replaced by Blocked k-fold CV (BCV) with blocks of 8 quarters as in, e.g., MRFjae. This is a more principled strategy to curb the downward bias of $\lambda$ induced by serial dependence. Among other things, it does not rely on blending two models using intuitive but admittedly arbitrary weights.

\vskip 0.2cm

{\sc Results with Basic CV and Half & Half.} Overall, results are in line with evidence previously reported in the TVP literature: very limited improvements are observed for real activity variables (GDP, UR) whereas substantial gains are reported for INF and IR. For the latter, allowing for time variation in either AR or a compact factor model (ARDI) generate very competitive forecasts. For instance, ARDI-2SRR is the best model for IR with a reduction of 36% in RMPSE over the AR(2) benchmark which is strongly statistically significant. Still for IR, at horizon 2 quarters, iterating 2SRR to obtain MSRR$_{\text{S}}$ generate sizable improvements for both AR and ARDIs. VAR(20) is largely inferior to alternatives in any of its forms. Two exceptions are IR forecasts at a one-year horizon where combining VAR(20) with MSRR$_{\text{D}}$ yields the best forecast by a wide margin with improvement of 19% in RMSPE. VAR(20)-MSRR$_{\text{D}}$ also provide a very competitive forecast for IR at an horizon of one quarter. Finally, at horizon 1 quarter, any form of time variation (2SRR, MSRR$_{\text{S}}$, MSRR$_{\text{D}}$) at least increases SPREAD's forecasting accuracy for for all models but the VAR(20). Precisely, it is a 16% reduction in RMSPE for AR, about 5% for ARDI and up to 14% in the VAR(5) case. For the latter, its combination with 2SRR provides the best forecast with a statistically significant improvement of 18% with respect to the AR(2) benchmark.

A notable absence from the relatively cheerful discussion above is inflation, which is the first (or second) variable one would think should benefit from time variation. It is clear that, in Table (ref), any AR at horizon 1 profits rather timidly from it. A similar finding for Half & Half is reported in Table (ref). What differs, however, are longer horizons results for INF. Indeed, mixing in additional shrinkage to OLS strongly helps results for those targets: every form of time variation now improves performance by a good margin. For instance, any time-varying ARs improves upon the constant benchmark by around 15%. It is now widely documented that inflation prior to the Pandemic was better predicted by past values of itself and not much else -- besides maybe for recessionary episodes KLS2019. Results of Table (ref) comfortably stand within this paradigm except for the noticeable efforts from MSRR$_{\text{S}}$ versions of both ARDI and VAR(20). While those are the best models, they are closely matched in performance by their AR counterparts. Nevertheless, it is noteworthy that this surge in performance mostly occurs for their sparse TVP versions, suggesting time variation is likely crucial for more sophisticated inflation forecasts not to be off the charts. Finally, additional shrinkage marginally improves GDP forecasting at the two longer horizons, with the Half & Half ARDI-MSRR$_{\text{S}}$ providing the best forecasts.

\vskip 0.2cm

{\sc Blocked CV and an Overall Assessment.} Coming to Table (ref), we see that switching from CV to Blocked CV offers comprehensive gains. Large losses to constant parameters models for real activity variables have been substantially mitigated and sometimes turned into marginal gains (for instance, ARDI-2SRR, unemployment at $h=1$). Overall, at $h=1$, gains from BCV are located where those of Half & Half are, which is intuitive. Things are more heterogeneous at longer horizons. The most notable amelioration is that of inflation, where 2SRR (and variants) now brings important gains across the board. BCV being such a game changer for INF requires little creativity in explaining it. Inflation is highly persistent for a great part of the training sample, which may have misled plain CV in choosing too low of a $\lambda$. Note that, interestingly, BCV does not fare as well for IR as other strategies did. One way to rationalize this is that "arbitrarily" opting for a lower $\lambda$ through non-blocked CV may have helped with the massive and rapid structural change that occurred when the zero-lower bound was hit in the test sample. We can posit that this event with no historical precedent likely induced desirable time variation in $\lambda$ itself. Thus, for a panoply of unforeseen reasons, it is absolutely possible that $\lambda^\text{CV}$ is closer to the ex-post optimal $\lambda$ than $\lambda^\text{BCV}$, even if BCV dominates CV in principle for dependent data.

To a large extent, forecasting results suggest that the three main algorithms presented in Section (ref) can procure important gains for forecasting targets that are frequently associated with the need for time variation. This subset is put on the spotlight by Figure (ref). The gains for IR at $h=1$ and INF at the two quarters and one-year horizon are particularly visible. For those kinds of targets, it is observed that any form of time-variation will usually ameliorate the constant parameters benchmark, especially in the Half & Half and BCV cases. We also see that those perform comparably well or better than a state-of-the-art Bayesian approach, precisely the Triple Gamma prior model from cadonna2020triple using default hyperparameters.\footnote{Marginally less competitive results were obtained using the ShTVP-R configuration also used in simulations, so I only reportShTVP-3G for parsimony. } Of course, both the ridge and Bayesian approach have distinct merits, but it is interesting to note that ShrinkTVP marginal advantage in simulations seems to be reversed in this real data application.

Time-Varying Effects of Monetary Policy in Canada

This section shows how to conduct both a TVP-VAR analysis and TVP-local projections (henceforth TVP-LP) on a data set where currently available methods would struggle. In particular, I study the changing effects of monetary policy (MP) in Canada using the MP shocks series of champagnesekkel2018. The computational difficulties come from mainly two reasons. First, it is monthly data for 40 years, which amounts to 480 (generally noisy) observations. Many traditional methods have computations which scale badly in $T$, and as such, fare better on quarterly data. Aggregating to quarterly is not an option because it would complicate identification, by needlessly aggregating into a simultaneous unit things that happened in different months. A far more challenging obstacle is that 24 lags need to be included in the model for VAR-based impulse response functions to have the "right" sign. This is true for OLS but also for TVP specifications tried here. Accordingly, the TVP-VAR generalization of the "correct" OLS specification will have 97 parameters in each equation for the 4-variate case, and 193 in the 8-variate case.\footnote{Thus, computational struggle mainly comes from conditional means rather the covariance matrix, but as we have seen computations for the latter scales rather gently in $M$ when using the tools provided in Section (ref).} Similarly daunting number of lags are utilized in their LP specifications.

I use the same monthly Canada data set as in champagnesekkel2018 and the analysis spans from 1976 to 2015. The target variables are unemployment, CPI Inflation and GDP.\footnote{For reference, the three time series being modeled and the shock series can be visualized in Figure (ref). Notably, we can see that the conquest of Canadian inflation (prior to Covid) was done in two steps: reducing the mean from roughly 8% to 5% in the 1980s and from 5% to 2% in the early 1990s.} The small open economy went through important structural change over the last 30-40 years. Most importantly, from a monetary policy standpoint, it became increasingly open (especially following NAFTA) and an inflation targeting regime (IT) was implemented in 1991 -- a specific and publicly known date. Both are credible sources of structural change in the transmission of monetary policy.

\vskip 0.2cm

{\sc Specification & Implementation : TVP-VAR.} champagnesekkel2018 estimate a parsimonious VARs (4 variables) over two non-overlapping subsamples to check visually whether a break occurred in 1992 following the onset of IT. The reported evidence for a break is rather weak with GDP's response increasing slightly while that of inflation decreasing marginally. While the sample-splitting approach has many merits such as transparency and simplicity, there is arguably a lot it can miss.

I consider a TVP-VAR(4) with their four original predictors (GDP, inflation, commodity prices, MP shocks), and a TVP-VAR(8) including additionally unemployment, USD-CAD exchange rate, exports, and imports. In such dimensions where the time-invariant model itself can overfit, it is preferable to set $\lambda_0$ in equation (ref) (regularization for $\beta_0$) to a value greater than 0. I set it to the value obtained by cross-validating a plain (thus time-invariant) ridge regression with $\boldsymbol{X}$ as predictors. Evidently, this is an educated guess for what one could obtain cross-validating the whole thing over a grid of $\lambda$ and $\lambda_0$. Also, the "prior mean" for $\beta_0$ is re-centered at the values corresponding to the coefficients of the time-invariant ridge regression. As discussed in the forecasting experiments, opting for a Blocked CV that takes into account the dependence between data points can help in avoiding overfitting. The block size is set to 24 months. $\varphi$, the sole tuning parameter of the covariance matrix, is set to 1000, which made the variances processes move at about the same speed as we typically see in stochastic volatility models. Increasing it in 2000-3000 range was found to eventually eradicate time variation and decreasing it below 500 made appear unrealistic higher frequency movements (especially in the covariances) that random walks are not built to capture efficiently. The ordering in the TVP-VAR(4) is the same as in champagnesekkel2018 (with the MP shocks entering last) and the additional 4 variables in TVP-VAR(8) are also added before the MP shock.\footnote{The exact ordering for the TVP-VAR(4) is GDP, inflation, commodity price index, MP shock. In TVP-VAR(8), we have GDP, inflation, commodity price index, unemployment, USD-CAD exchange rate, exports, imports, MP shock. } Lastly, relative IRFs (dividing by the contemporaneous effect of MP shock on itself) are reported as they are directly comparable to TVP-LPs plagborg2021local and that each IRFs desirably corresponds to a unit shock at each $t$.

In terms of computational results, the TVP-VAR(4) with a total $4 \times 97$ conditional mean TVPs took 25 seconds to estimate and tune. It took an additional 6 seconds to calculate 10 (co)variances TVPs and calculate/generate IRFs for 48 horizons and 4 variables. The TVP-VAR(8), with $8 \times 193$ conditional mean TVPs and 36 additional TVPs in the covariance matrix, took 1:17 minute for the first step, and 17 seconds for the second. Hence, 2SRR lives up to the claim that it can estimate demanding models very quickly with little user interaction.

\vskip 0.2cm

{\sc Specification & Implementation : TVP-Local Projections.} VARs do not have a monopoly on the proliferation of parameters. jorda2005 local projections' -- by running a separate regression for each horizon -- are also densely parametrized. For that reason, constructing a large LP-based time-varying IRF via a MCMC procedure would either be burdensome or unfeasible. The use of 2SRR for estimation of LPs constitutes a very useful methodological development given how popular local projections have become over recent years. It is now known that IRFs from LPs and VARs (coming from linear models) are intimately linked plagborg2019local. However, here, since we are modeling time variation in a different space than the TVP-VAR (i.e., directly in IRF-space rather than in reduced-form VAR coefficients and then the covariance matrix), it is natural to expect that results can differ for reasons other than finite sample estimation variability. Moreover, cross-validation is done for direct forecasts at each horizon in LPs whereas it is done for each variable and one-step ahead forecasts in the TVP-VAR case. Hence, from a methodological standpoint, it is interesting to compare TVP-VAR and TVP-LP results on a terrain where full agreement is not necessarily the expected outcome.

champagnesekkel2018 original LP specification includes 48 lags of the narrative monetary policy (MP) shock series which is constructed in the spirit of romerromer2004 and carefully adapted to the Canadian context.\footnote{For details regarding the construction of the crucial series --- especially on how to account for the 1991 shift to IT, see champagnesekkel2018. Note that a positive shock means (unexpected) MP tightening.} Furthermore, their regression comprises 4 lags for the controls which are first differences of the log GDP, log inflation and log commodity prices. To certify that time-variation will not be found as a result of omitted variables, I increase the lag order from 4 to 6 months and augment the model with the USD/CAD exchange rate, exports, imports and CPI excluding Mortgage Interest Cost (MIC). In terms of TVP accounting, $\boldsymbol{X}$ contains 97 regressors (including the constant) and $\boldsymbol{Y}$ is $48 \time 419$. Thus, a single TV-LP is assembled from a staggering total of $97 \times 48=4\,656$ TVPs. The computing time one IRF is 26 minutes. Thus, an evident advantage of the VAR approach is that by estimating a single generative model, one can get results considerably faster.

\vskip 0.2cm

{\sc Results.} Current evidence on the effects of IT is mixed. ramey2016hb and barakchian2013monetary obtain counterintuitive results (a price puzzle and MP tightening increases GDP) for post-1988 US data. cloyne2016macroeconomic and champagnesekkel2018 report more intuitive findings for the UK and Canada: GDP response increases marginally after IT and that of inflation shrinks. Figure (ref) suggests new answers to the evolving effect of monetary policy on the Canadian economy based on 3 TVP models. Results mostly follow a qualitative common thread. There are quantitative differences, but this is expected that as even OLS versions of LPs vs. VARs give quantitatively different answers in champagnesekkel2018. Yet, in a pairwise fashion, TVP-LPs signs and magnitudes are consistent with those reported in champagnesekkel2018, and TVP-VARs similarly concur with time-invariant VARs. Thus, when comparing IRFs, it is good to keep in mind that there are two potential sources of dissimilarity: those that emanate from the linear specification itself, and those that will be attributable to how time variation is captured.

Generally, short- and medium-run effects ($h<24$ months) have been much more stable than longer-run ones. This is especially true of unemployment which exhibits a very homogeneous response (through time) for the first year and half after the shock. GDP's response roughly follows a similar pattern. When it comes to inflation, its usual long response lag has mildly shortened in the 2000s, a phenomenon that is visible in all 3 specifications. For instance, in the TVP-LP case, at a horizon of 24 months, the effect of a positive one standard deviation shock was -0.3% from 1976 to the late 1990s, then slowly increased (in absolute terms) to nearly double at -0.6%. A similar relative decrease is observed for TVP-VARs, but their overall absolute size throughout the sample is smaller. This is in line with champagnesekkel2018's time-invariant LPs vs. VAR where the overall IRF magnitudes of the latter are always smaller than the former. Nonetheless, all results for horizons up to 18 months suggest that the ability of the central bank to "rapidly" impact inflation has increased. Despite these insights coming from a sample ending in the mid-2010s, they appear roughly in accord with the post-Pandemic Canadian inflation experience, where the Bank of Canada's rate hikes arguably had a material impact (among a myriad of confounding factors) within 1.5 years.

Given the long lags of monetary policy, most of the relevant action from an economic standpoint is also where most time variation is found: from 1.5 to 4 years after the shock. For the two models that include it, the cumulative long-run effects of MP shocks on unemployment have substantially shrunk over the sample period. For TVP-LP , the decrease from a 0.6 to 0.4 unemployment percentage points mostly occurred throughout the 1980s, and stabilized at 0.4 thereafter. In TVP-VAR(8), the overall effect (even in the OLS case) is more muted, but the reduction of long-term effects is also visible and of similar relative size. Additionally, both models showcase the same attenuation pattern for mid-range horizons, that is, the peak effect shrinking to about $\sfrac{2}{3}$ of its early 1980s value. All inflation IRFs show a declining pattern starting from the early 1990s, but that of TVP-LP is more visible given its overall smoother behavior. Pre-1990, we see TVP-VAR based IRFs displaying mild business-cycle frequency movement, but those disappear post-1990.

Finally, GDP IRFs all display a U-shape pattern, with a trend reversal in the strength of long-run effects happening in the 1990s. However, there is discord between TVP-VAR-based estimates as to whether the effects were milder at the beginning or the end of the sample. TVP-VAR(4) and TVP-LPs agree on it being at the end, in line with results for unemployment, but TVP-VAR(8) sees this happening in the late 1970s instead. Obviously, decisively choosing between one model would require further work, and in the machine learning spirit, the model (here, direct vs. iterative TVP-VAR forecasts) that predicts best for the contentious horizons should be preferred. Nonetheless, key findings are homogeneous across specifications, that is, there is (i) a trend reversal in the strength long-run effects of MP shocks on inflation and GDP (getting stronger for the former and weaker for the latter) and (ii) the effects on unemployment became milder starting from the 1990s for all $h$'s greater than 18 months.

\vskip 0.2cm

{\sc Zooming on 1992.} An important question is what happens to $\beta_t$ around 1992, after the onset of IT. Neither TVP-VAR nor TVP-LPs suggest the occurrence of a structural break in 1992, which is in line with most of the international evidence on IT implementation. For inflation, TVP-LPs' results point to a change in coefficients' trending behavior, a subtle phenomenon which would effectively stay under the radar of simpler approaches. For illustrative purposes, Figure (ref) (Appendix) reports $\beta_{t}^{\text{TVP-LP}}-\beta^{\text{OLS}}$ for the dynamic effect of MP shocks on the three variables. The response of inflation (in absolute terms) is much larger at the end of the sample than what constant coefficients would suggest. This is especially true at the 24 months horizon. It is quite clear that, for all horizons, the effect of MP shocks on inflation starts increasing in the years following the implementation of IT. The downward trend in MP shocks' impact on unemployment seems to have started at least since the 1970s and slowed down in recent years.

More generally, it is interesting to note that those results are consistent with some flattening of the unemployment- or GDP-based Phillips' curve (but maybe not a more general one, HNN) which was a well-established observation, at least, before the COVID-19 pandemic blanchard2015inflation,del2020s. The combined evidence displayed in Figure (ref) is that, while the cumulative effect of MP shocks either stood still or became more muted starting from the 1990s for real activity variables, it has increased for inflation. This suggests that stabilizing Canadian inflation is less costly (in terms of unemployment/GDP variability) within the IT framework. Current evidence from the post-Pandemic rates hikes cycle seems to partly validate this view.

\vskip 0.2cm

{\sc Testing Limits.} It is not because we can compute something that we necessarily want to. To test the limits of the ridge approach, I also computed results from a TVP-VAR(23), again including the necessary 24 lags, which makes it for a model that, even with time-invariant coefficients, has more parameters than observations ($23 \times 553 + 276$ in total). The first 8 variables and that of TVP-VAR(8) and the remaining 15 are the US variables included in UK database constructed in GCMS and can be seen as a reasonably sparse characterization of economic conditions facing Canada's main commercial partner. In such an environment, it is difficult to avoid the ridgeless solution mentioned earlier, and all IRFs are much smaller in magnitude falling under the weight of heavy shrinkage. Nonetheless, it can be computed. It takes 10 minutes to estimate and less than 3 minutes to run the code evaluating the covariance matrix, calculating IRFs, and generating 3D plots. Hence, it takes less than 15 minutes to obtain results. However, in such dimensions, things are heavily compressed and time-variation patterns obtained do not seem as economically plausible. Reducing the number of lags to 12 and 6, IRFs become less compressed and ridgeless issues go away, but as discussed in the intro of this section, puzzles emerge. Computations are also much less demanding, with the whole operation taking 3:30 minutes. Thus, in this environment where an abundant number of lags appear to be necessary, the number of variables in the system cannot be too large.

Conclusion

I provide a new framework to estimate TVP models with potentially evolving volatility of shocks. It is conceptually enlightening and computationally very fast. Moreover, seeing such models as ridge regressions suggest a simple way to tune the amount of time variation, a consequential quantity. The approach is easily extendable to have additional shrinkage schemes like sparse TVPs or reduced-rank restrictions. The proposed variants of the methodology are very competitive against the standard Bayesian TVP-VAR in simulations. Furthermore, they improve forecasts against standard forecasting benchmarks for variables usually associated with the need for time variation (US inflation and interest rates). Finally, I apply the tool to estimate time-varying IRFs via VAR and local projections. The large specifications necessary to adequately characterize the evolution of monetary policy in Canada rendered this application nearly unfeasible without the newly developed tools. I report that monetary policy shocks long-run impact on the price level increased substantially starting from the early 1990s (onset of inflation targeting), whereas the effects on unemployment became milder.

\setlength\bibsep{5pt}

\setstretch{0.3}