EconBase
← Back to paper

Estimating Large Mixed-Frequency Bayesian VAR Models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

425,755 characters · 13 sections · 94 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
This text was truncated for display. The citation measures were computed over the complete text.
titlepage\pagestyle{empty} \vskip 1em \begin{center} \let \footnote \vskip 1.5em { \lineskip .5em \begin{tabular}[t]{c} Sebastian Ankargren\footnotemark[2] and Paulina Jon\'{e}us\\ Department of Statistics, Uppsala University \end{tabular}} \vskip 1em \end{center} \vskip 1.5em \begin{abstract} We discuss the issue of estimating large-scale vector autoregressive (VAR) models with stochastic volatility in real-time situations where data are sampled at different frequencies. In the case of a large VAR with stochastic volatility, the mixed-frequency data warrant an additional step in the already computationally challenging Markov Chain Monte Carlo algorithm used to sample from the posterior distribution of the parameters. We suggest the use of a factor stochastic volatility model to capture a time-varying error covariance structure. Because the factor stochastic volatility model renders the equations of the VAR conditionally independent, settling for this particular stochastic volatility model comes with major computational benefits. First, we are able to improve upon the mixed-frequency simulation smoothing step by leveraging a univariate and adaptive filtering algorithm. Second, the regression parameters can be sampled equation-by-equation in parallel. These computational features of the model alleviate the computational burden and make it possible to move the mixed-frequency VAR to the high-dimensional regime. We illustrate the model by an application to US data using our mixed-frequency VAR with 20, 34 and 119 variables. \end{abstract} \footnotetext[1]{We are thankful for comments by Johan Lyhagen and seminar participants at the Department of Statistics, Uppsala University and Sveriges Riksbank. Part of this research was conducted while Ankargren was visiting Sveriges Riksbank, which is gratefully acknowledged. We also wish to thank the Swedish National Infrastructure for Computing (SNIC) through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under projects SNIC 2015/6-117 and 2018/8-361 for providing the necessary computational resources.} \footnotetext[2]{Corresponding author: S. Ankargren, [email removed]}

Introduction

An almost unavoidable feature of macroeconomic datasets used in real-time forecasting is that data are inherently unbalanced. There are mainly two reasons for unbalancedness of datasets: different sampling frequencies, and publication delays causing missing values at the end of the sample. For example, the gross domestic product (GDP), a key economic variable, is available quarterly, whereas many other central economic variables such as inflation, unemployment and industrial production are predominantly available monthly, usually with a delay of one or two months. Some variables, such as variables related to the stock market and interest rates, are often available daily, if not even more frequently. Typically, VAR models are estimated on a single-frequency basis and, consequently, the inclusion of a quarterly variable like GDP determines whether the model will be estimated on monthly or quarterly data.

Econometric models taking information in unequal frequencies into account, and thereby avoiding a loss of information stemming from aggregation to the lower frequency, have in recent years gained attention under the name of mixed-frequency methods. Multiple approaches are available, see for example Foroni2013 for a review. One of the methods put forward in the literature is the mixed-frequency vector autoregressive (VAR) model.

The use of factor models, as in the key contributions by Mariano2003,Mariano2010 is another way of handling mixed-frequency data that is particularly popular when the set of variables available is large; see also Marcellino2016. Other proposed approaches include the Mixed data sampling (MIDAS) regression and MIDAS-VAR proposed by Ghysels2007 and Ghysels2016 respectively. MIDAS regressions cope with the issue of unequal frequencies by regressing a low-frequency dependent variable on its lags as well as on lags of high-frequency variables. In the observation-driven MIDAS VAR model, the vector of endogenous variables includes both high and low-frequency variables and are formulated in terms of observable data. In contrast to the present paper, they do not involve latent processes. In the same spirit, McCracken2015 proposed a high-dimensional variant of the MIDAS-VAR. They use a mixed-frequency Bayesian VAR estimated at the lowest common data frequency and find that mixed-frequency information is important for nowcasting.

Our focus is a situation in which GDP growth is the sole quarterly variable, with the remaining being monthly variables. What we investigate is how a large set of monthly variables can be used when modeling GDP growth in a mixed-frequency VAR. The method we employ in this paper is to cast the mixed-frequency model in a state-space form. By doing so, we can essentially interpolate the latent monthly values of the quarterly variables concurrently with making inference about the parameters in the model. Eraker2015 formulated a Bayesian VAR based on this idea and Schorfheide2015 proceeded with a Gibbs sampling approach based on simulation smoothing and forward-filtering, backward-smoothing along the lines of Carter1994. Schorfheide2015 found that the model improved forecasting as compared to a quarterly VAR model. In an extension of Schorfheide2015, Ankargren2018 developed a steady-state mixed-frequency VAR model for a real-time US dataset and arrived at the same conclusion. Evidence of improved forecasting performance was also presented by Gotz2018.

In parallel with the development of mixed-frequency VARs, the empirical success of large VARs in forecasting macroeconomic variables has been highlighted in a number of papers, including among others Banbura2010 and Koop2013. There are multiple options for working with VARs in a high-dimensional setting, where the number of parameters is large. Setting the level of shrinkage in relation the dimension of the model was discussed by Banbura2010, who subsequently showed that larger models can improve forecasting ability when doing so. In more recent work, new methods drawing from the high-dimensional literature have been developed, see e.g. Korobilis2019 who developed a highly scalable estimation method for BVARs, Koop2019 who proposed a model based on ideas from the compressed regression literature and KastnerHuber2018,Follett2019,Huber2017 who used global-local shrinkage priors.

In addition to the usefulness of large VARs, seminal work by Cogley2005 and Primiceri2005 highlighted the importance of letting parameters and volatilities in VARs vary over time; see also Clark2011,DAgostino2013,Carriero2015,Carriero2016 for more forecasting-related studies. Different solutions have been presented to circumvent the computational difficulties regarding high dimensions in combination with time-varying volatilities. Carriero2019 extended earlier work on estimation of Bayesian VARs with asymmetric priors and time-varying volatilities to the high-dimensional setting, where the assumption of a specific common structure for the volatilities was relaxed and a simple triangularization of the VAR was used to reduce computational complexity. KastnerHuber2018 proposed an alternative route that instead used a factor stochastic volatility model that allows for co-movements in the error covariance structure when estimating large-dimensional VAR models.

In terms of mixed-frequency VARs with stochastic volatility, Cimadomo2016 extended earlier work on stochastic volatility to allow for data sampled at different frequencies. Due to quickly increasing computational complexity, their proposed model is restricted to include only a small number of variables and lags. Gotz2018, also closely related to our work, proposed a model for a moderately large set of variables where the intercept and the common factor in the error variances are allowed to vary over time. Their application used 11 variables and 6 lags, and the authors stated that this model could be enlarged to up to 20 variables. In the case of more variables or when the length of the sample is increased, they mention a considerable obstacle in the running time. Using a standard normal inverse Wishart prior, Brave2019 demonstrated that a larger, 37-variable mixed-frequency VAR outperformed the smaller model of Schorfheide2015.

Our contribution is that we add to this growing literature on large mixed-frequency VARs. We propose a large mixed-frequency VAR where we model the time-varying error covariance matrix by a factor stochastic volatility model along the lines of KastnerHuber2018. In contrast to KastnerHuber2018, we include the possibility to estimate the model on mixed-frequency data. Apart from being a parsimonious approach for modeling time-varying volatility, use of the factor stochastic volatility model comes with computational advantages that we exploit for estimating our model when dimensions are high. The factor stochastic volatility model renders the equations in the model conditionally independent. Consequently, the adaptive simulation smoother developed by Ankargren2019 can be improved further by coupling it with the univariate filtering methodology suggested by Koopman2000. Moreover, the conditional independence of the equations means that we can sample the rows of the regression parameter matrices independently in parallel. These computational improvements ameliorate the computational efficiency of our Markov Chain Monte Carlo (MCMC) sampler, thereby enabling faster estimation for large-scale mixed-frequency VARs. The proposed model is illustrated using the FRED-MD database constructed by McCracken2016, where we consider three models of increasing size with 20, 34 and 119 variables. In summary, the model and computational strategies we propose substantially reduce the computational burden induced by the model, and, as such, provide an efficient and feasible way of estimating VARs with mixed-frequency data and stochastic volatilities for models of dimensions that have up to now been missing in the literature. A benefit of our approach is that we do not assume a Kronecker structure for the prior variance of the regression parameters in contrast to e.g. Gotz2018,Ankargren2018, who used the methodology presented by Carriero2016. Consequently, the model that we will discuss permits asymmetric prior distributions and can therefore be used in conjunction with high-dimensional hierarchical shrinkage priors.

The rest of the paper is organized as follows. Section (ref) introduces the mixed-frequency VAR, Section (ref) presents our suggested approach for large mixed-frequency VARs, Section (ref) develops a fast simulation smoother exploiting the structure of the model, and Section (ref) provides an empirical illustration. Section (ref) concludes.

Mixed-Frequency VARs

The mixed-frequency VAR($p$) is specified at the monthly frequency as

align[align omitted — 101 chars of source]

where $x_t$, $\Pi_c$ and $u_t$ are $n\times 1$ and $\Pi_j$ are $n\times n$. The vector $x_t$ consists of both monthly and quarterly variables and is partitioned as $x_t=(x_{m,t}', x_{q,t}')'$, where the dimension of $x_{m, t}$ is $n_m\times 1$ and and $x_{q, t}$ is of dimension $n_q\times 1$, where $n=n_m+n_q$.

In a mixed-frequency setup, the dataset that is needed for estimating the VAR is not fully observed. Instead, every month $t$ we observe $n_t\leq n$ variables collected in the vector $y_t=(y_{m, t}', y_{q,t}')'$, where $y_{m, t}$ and $y_{q,t}$ are the observed monthly and quarterly subsets of dimensions $n_{m, t}\times 1$ and $n_{q,t}\times 1$, respectively. The vector of observables, $y_t$, relates to the underlying $x_t$ through

equation[equation omitted — 196 chars of source]

where $S_{m,t}$ and $S_{q,t}$ are deterministic selection matrices of dimensions $n_{m, t}\times n_m$ and $n_{q, t}\times n_q$, respectively. The matrix $\Lambda_{qq}$ is an $n_q\times pn_q$ aggregation matrix that aggregates the latent monthly series into observed quarterly variables according to a selected aggregation scheme. We use data transformed to stationarity and therefore employ the triangular aggregation suggested by Mariano2003, whence $y_{q, t}=\frac{1}{9}(x_{q, t}+2x_{q,t-1}+3x_{q,t-2}+2x_{q,t-3}+x_{q,t-4})$ when $t$ corresponds to an end-of-quarter month (with $y_{q,t}=\varnothing$ otherwise).

The posterior distribution that we are interested in is $p(\Pi, \Sigma, x|y)$.\footnote{We will throughout the paper make use of the notational convention that variables without subscripts refer to the whole set of that variable, i.e. $\Pi=\{\Pi_i : i = 1, \dots, p\}$.} Schorfheide2015 developed a Gibbs sampler for estimating the model using the normal-inverse Wishart prior distribution (see Kadiyala1993,Kadiyala1997) for the regression parameters and the constant error covariance matrix. The Gibbs sampler uses as one of its steps a simulation smoothing algorithm Carter1994,FruhwirthSchnatter1994 for making a draw from the conditional posterior distribution $p(X|Y, \Pi, \Sigma)$. The basic Gibbs sampler for the mixed-frequency VAR is to sample repeatedly from:

align[align omitted — 54 chars of source]

The first step is standard in the literature conditional on $x$, and conditionally independent of $y$. Intuitively, the $p(x|\Pi, \Sigma, y)$ step imputes the underlying monthly series of all variables observed quarterly as well as monthly series with missing values at the end of the sample. Next, conditional on $x$, draws of $\Pi$ and $\Sigma$ can be produced as if $x$ were the real data, leading to standard procedures. Given the modular property of MCMC, it is therefore possible to extend the underlying VAR model using existing methods employed for fully-observed VARs by essentially adding the mixed-frequency sampling step to existing MCMC algorithms; see Cimadomo2016,Gotz2018,Ankargren2018 for examples. It is this modularity that we exploit in this paper.

A Large Mixed-Frequency VAR with Factor Stochastic Volatility

In this paper, we estimate three mixed-frequency models of increasing size. We attempt to stay relatively close to models used in the previous literature and use pre-existing single-frequency models now augmented with mixed-frequency data. Our purpose for doing so is to model GDP growth using VARs and a large set of monthly data.

The three models that we use are based on models used by Koop2013 and Carriero2019. Koop2013 used a medium-size quarterly VAR with 40 variables, all transformed to the quarterly frequency, and Carriero2019 used two monthly models with $n=20$ and $n=125$, respectively. After excluding discontinued and short series, we end up with three models of dimensions: $n=20$, $n=34$ and $n=119$, where each contains $n-1$ monthly series and quarterly GDP growth. For ease of presentation, the models are referred to as CCM-20, Koop-34 and CCM-119, respectively. The current section describes the data, some of their important characteristics, and some key aspects of the models. The more technical details, including the complete set of prior specifications, are relegated to Appendix (ref).

Data Structure and Publication Schemes

The data we use are retrieved from FRED-MD McCracken2016. We use the January 2019 vintage and collect a typical observational pattern of the included variables from the St. Louis Federal Reserve's database ALFRED. Because the Koop-34 model includes a number of series no longer published in the FRED-MD vintages, the Koop-34 model uses data through December 2014. The data are transformed according to the transformations suggested by McCracken2016 and thereafter standardized. See Appendix (ref) for a full list of the variables used.

figure[figure omitted — 57,134 chars of source]

A crucial ingredient in real-time forecasting is the observational pattern in the data. Figure (ref) shows the number of monthly variables with zero, one or two missing observations by day-of-month and by model. At the beginning of the month, a number of variables contain zero missing observations indicating that the outcome of the previous month is already known. These variables include interest rates, spreads and stock market indexes. A few days into the month, new outcomes are obtained. These are predominantly labor market variables. The next major change in observational pattern is around the middle of the month. By this time, most real variables have been observed. However, some of the slow real variables are not observed until the end of the month. As indicated by the bottom blue bar in each panel, some of these are also not observed with a single-month lag, but a two-month lag for the most part.

The quarterly variable GDP growth is assumed to be known in the beginning of the last month in the following quarter; that is, GDP growth for the third quarter is available in the beginning of December.

Modeling Stochastic Volatility

The stochastic volatility model used in this paper is based on the idea of parsimoniously modeling the time-varying volatility. To this end, we use a factor stochastic volatility model and decompose the error term in (ref) as

align[align omitted — 37 chars of source]

where $f_t$ is an $r\times 1$ vector of factors, $\nu_t\sim N(0, \Omega_{t}^\nu)$ and $r\ll n$. Furthermore, $f_t\sim N(0, \Omega_t^f)$ and both $\Omega_t^\nu$ and $\Omega_t^f$ are diagonal with diagonal elements $\omega_{i, t}^\nu$ and $\omega_{j, t}^f$ where the law of motion is

alignat{3} &\log \omega_{i, t}^\nu&& = \mu_i+\phi_i^\nu(\log \omega_{i, t-1}^\nu-\mu_i)+\sigma_i^\nu e^\nu_{i, t}, \quad &&e^\nu_{i, t} \sim N(0, 1)\quad i=1,...n\\ &\log \omega_{j, t}^f &&= \phi_j^f\log \omega_{j, t-1}^f+\sigma_j^fe^f_{j, t}, \quad &&e^f_{j, t} \sim N(0, 1) \quad j=1,..., r.

The covariance matrix of the VAR innovation is therefore $\Sigma_t=\Lambda_f\Omega_{t}^f\Lambda_f'+\Omega_{t}^\nu$ and thus allows for idiosyncratic stochastic volatilities as well as time-varying covariances by means of the common component. An efficient way of estimating the factor stochastic volatility model was proposed by Kastner2017, where the authors made use of the interweaving strategy developed by Yu2011 and its application to univariate stochastic volatility models by Kastner2014. Moreover, KastnerHuber2018 used the approach in a VAR framework demonstrating the merits of the method.

Sampling Regression Parameters

Conditional on everything but the regression parameters, the original model in (ref) can be written as

align[align omitted — 83 chars of source]

where $X_t=(1, x_{t-1}', \dots, x_{t-p}')'$. It is evident from this formulation that the multivariate model consists of $n$ heteroskedastic but independent regressions of $\tilde{x}_{i, t}$ on $X_t$. By letting $\ddot{x}_{i, t}=\frac{1}{\sqrt{\omega^\nu_{i,t}}}\tilde{x}_{i,t}$ and $\ddot{X}_{i, t}=\frac{1}{\sqrt{\omega^\nu_{i,t}}}X_t$ we obtain the factorization

align[align omitted — 152 chars of source]

In the new parametrization of the model, draws of each row of $\Pi$ can be made independently of each other. This attractive side-effect of the factor stochastic volatility model for sampling the regression parameters was already noted by KastnerHuber2018, but we want to stress the benefit thereof. Because of the independence, a parallel implementation of this sampling step is easily employed. Carriero2019 present a sequential sampling algorithm to sample $\Pi$ row-by-row that dramatically alleviates the issue of sampling from the posterior of $\Pi$, but because it is sequential it is not parallelizable.

Nevertheless, the task of sampling $\Pi$ is---even using parallelization---challenging when dimensions are high. Standard generators of pseudo-random numbers from multivariate normal distributions involve a Cholesky decomposition of the covariance matrix, an operation of order $O(n^3p^3)$ that is to be carried out once per equation. However, because of the special structure of the posterior distribution some improvements can be achieved by using the sampling methods proposed by Rue2001 or Bhattacharya2016, which utilize the specific structure of the posterior. The complexity of the Bhattacharya2016 sampler is $O(T^2np)$ and tends to be faster than the Rue2001 algorithm when $np>T$ as the latter has a complexity of $O(n^3p^3)$.

Simulation Smoothing Using Univariate Adaptive Filtering

In order to sample from the posterior of the latent variables, Schorfheide2015 proposed an alternative representation---called the compact form---of the mixed-frequency VAR in which the monthly variables are treated as exogenous. The benefit of leaving the monthly variables out of the state equation is that the computational burden is greatly reduced. However, in real-time forecasting situations when data contain ragged edges, some monthly variables are missing at the end of the sample. This characteristic of the data means that these missing monthly variables will need to be included back in to the state vector. Schorfheide2015 solved this issue by simply moving from the efficient compact form to the full companion form VAR when the first monthly variable goes missing. As Ankargren2019 demonstrated, this leads to a large bottleneck when the dimension of the model is high. Consequently, Ankargren2019 developed an adaptive simulation smoothing algorithm based on the compact form suggested by Schorfheide2015. The adaptive part of the algorithm implies that, as some monthly variables go missing at the end of the sample, these are added to the state vector, but not the entire set of monthly variables. By adaptively augmenting the state vector, large efficiency gains can be obtained, particularly when the number of variables and the number of lags are high. For more details, we refer to Ankargren2019; the compact and companion forms are also presented in full in Appendix (ref) for self-containment.

Use of the adaptive procedure improves upon filtering of the unbalanced part to the extent that it no longer constitutes the bottleneck of the algorithm when dimensions get higher. It is now instead the filtering step in the compact form that obstructs an efficient procedure. Fortunately, it is possible to obtain large computational improvements for this chunk of the algorithm by a judicious choice of stochastic volatility model. The method that can be used to further alleviate the procedure of its computational burden is the univariate filtering approach for multivariate time series by Koopman2000. Univariate filtering primarily excels when the dimension of the observation vector is large relative to the state vector, which will be the case in the compact form for a large-dimensional mixed-frequency VAR if the number of quarterly variables remains low. The approach is applicable when the state and observation errors are independent and the observation errors feature a diagonal covariance matrix. In what follows, we will discuss how the factor stochastic volatility model allows us to reap the computational benefits of the univariate filtering approach. Since the unbalanced part of the data in the adaptive procedure is responsible for only a negligible part of the computational burden, the discussion of the implementation of the univariate filtering method focuses on the balanced part of the data. There is no obstacle in applying it also on the unbalanced part, but for ease of exposition we do not discuss that any further. It should also be noted that simulation smoothing consists of three main steps: data generation, filtering, and smoothing. We will focus on filtering and smoothing exclusively in the following and refer to Ankargren2019 for more details on the complete simulation smoothing algorithm.

The Schorfheide2015 compact formulation of the state-space model for the mixed-frequency VAR is

equation[equation omitted — 451 chars of source]

where the system matrices $\{Z_t, C_t, G_t, T_t, D_t, H_t\}$ are functions of the original model parameters and $y_{m, t-1:t-p}=(y_{m,t-1}', \dots, y_{m, t-p}')'$; see Appendix (ref) for more details.

Since the simulation smoothing step discussed in this section is performed as a Gibbs sampling step, the procedure is carried out conditional on all other parameters in the model. Hence, at this stage $\Lambda_f f_t$ and $\Omega^{\nu}_t$ are known. Partition the factor loadings into monthly and quarterly blocks as $\Lambda_f=(\Lambda_{m, f}', \Lambda_{q, f}')'$ and let

align[align omitted — 254 chars of source]

It is now possible to formulate the compact version of the model as a state-space model with time-varying intercepts and diagonal error covariances:

align[align omitted — 373 chars of source]

The univariate filtering method treats the multivariate time series $\{y_t\}$ as a univariate time series. The univariate representation of the observation equation is

align[align omitted — 275 chars of source]

For simplicity, the observation equation does not reflect that $y_{t,i}$ is missing every third period for $i=n_m+1, \dots, n$. Strictly speaking, the observation is empty in these cases. The univariate representation of the state equation is

align[align omitted — 205 chars of source]

The error covariance matrix in the observation equation is diagonal. Univariate filtering exploits the implication of the diagonality in that it brings observations in one at a time instead of bringing the multivariate time $t$ observation in in one go. The computational justification for doing so is that matrix-to-matrix multiplications are avoided in favor of a larger number of scalar or matrix-to-vector multiplications. The latter approach is generally faster if the dimension of the state vector is low in relation to the size of the observation vector.

The standard Kalman filter iterates over the filtering equations to compute $a_t=E(\alpha_t|Y_{t-1})$ and $a_{t|t}=E(\alpha_t|Y_t)$, whereas the univariate approach also computes the intermediate sequence

align[align omitted — 72 chars of source]

for $i=2, \dots, n_t+1$; the connection to the standard filter is $a_{t,1}=a_t$ and $a_{t,n_t+1}=a_{t|t}$. The univariate filtering recursions are, for $t=1, \dots, T_b$ and $i=1, \dots, n_t$:

alignat{3} a_{t,i+1}&=a_{t,i}+K_{t,i}F_{t,i}^{-1}v_{t,i},&\quad P_{t,i+1}&=P_{t,i}-K_{t,i}F_{t,i}^{-1}K_{t,i}'\\ v_{t,i}&=y_{t,i}-Z_{t,i}a_{t,i}-c_{t,i}, &\quad F_{t,i}&=Z_{t,i}P_{t,i}Z_{t,i}'+\omega^\nu_{t, i}, \quad K_{t,i}=P_{t,i}Z_{t,i}',

where $Z_{t,i}$ is row $i$ of $Z_t$ and $c_{t,i}$ is the $i$th element of $c_t$. To transition from time $t$ to time $t+1$, the following relations are used:

align[align omitted — 115 chars of source]

Several of the Kalman filter's familiar terms---such as $F_t$---now appear as scalars instead of matrices. The consequence is that there is no need to invert the full $F_t$ matrix; instead, the only inversion required is that of the scalar $F_{t,i}$.

Finally, to conclude the description of the univariate filtering approach, the recursions used for conducting the smoothing step of the algorithm are:

equation[equation omitted — 198 chars of source]

for $i=n_t, \dots, 1$ and $t=T_b, \dots, 1$. The smoothed estimate of primary interest is computed as $\hat{\alpha}_t=a_t+P_tr_{t,0}$.

Quantifying the Computational Improvement

While the number of variables in the model determines the computational burden to a large extent, the choice of the lag length $p$ naturally plays a crucial role. In the large VAR literature, the choice of lag length is subject to a notable degree of variation. For example, KastnerHuber2018 used $p=1$, whereas Schorfheide2015 and Gotz2018 for mixed-frequency VARs set $p=6$, and Banbura2010,Carriero2019 let $p=13$. Studying VAR specifications more systematically, Carriero2015b found that optimizing the lag length improved the forecasting ability as compared to a baseline specification using $p=13$, although only by a few percentage points. Therefore, the conclusion one can draw is that lag lengths used in practice are relatively idiosyncratic, and that choosing a large lag length is a way to play it safe.

Choosing a large lag length in the mixed-frequency VAR context is more expensive than if data are all on the same frequency, since the dimension of the state vector is $n_q(p+1)$ and thereby enlarged when $p$ increases. Moreover, the more quarterly variables that are included (for a fixed $n$), the heavier the computational burden. Figure (ref) shows the cost of increasing the lag length for the mixed-frequency block of the algorithm for the three models discussed in Section (ref) with $n=20,34$ and 119, respectively, and a single quarterly variable. In addition, the figure displays the cost of increasing $p$ for the Schorfheide2015 model, which uses eight monthly and three quarterly variables. Of the eight monthly variables, three are available with no publication delay and five with a one-month delay. For the CCM-20, Koop-34 and CCM-119 models we use the observational pattern from day 15, see Figure (ref).

figure[figure omitted — 31,970 chars of source]

Figure (ref) shows that our adaptive algorithm with univariate filtering works well for the empirical models when there is a single quarterly variable. The increased computational effort that is needed for a higher lag length is modest, whereas the Schorfheide2015 algorithm does not scale well. The first panel of the figure shows that when there are three quarterly variables, the adaptive algorithm offers an improvement compared with the Schorfheide2015 procedure, while combining it with univariate filtering performs worse. The inferior performance of univariate filtering is not very surprising---the method primarily excels when the dimension of the state vector is small relative to the dimension of the observation vector. For $p>3$, the state vector is larger than the observation vector. In contrast, the state vector is always smaller than the observation vector for the other three models. Comparing the relative performances across model sizes also reveals that when $n$ increases, the relative gains of univariate filtering increase. The reason is that $n_q(p+1)$ as a share of $n$ decreases, which further pushes univariate filtering into its domain of superiority. In particular, when the model includes 119 variables (rightmost panel), the difference in running time between the adaptive algorithm with univaraite filtering and the Schorfheide2015 method is substantial---almost 16 seconds per iteration is required by the latter, as opposed to 0.4 seconds.

To summarize, there are two key points made by Figure (ref). First, the gains from using univariate filtering from a computational perspective become larger as the dimension of the model increases. Second, varying the lag length is relatively cheap using either the adaptive or the adaptive with univariate filtering algorithms, which thereby allow for larger freedom in selecting the number of lags.

Other Stochastic Volatility Models

The preceding discussion of how the univariate filtering approach can be used for the balanced part of the sample focused on the case when a factor stochastic volatility model is used to model time-varying variances in the model. This is by no means the only way to model the stochastic volatility component. A popular alternative is to decompose $\Sigma_t=A_t^{-1}\Lambda_t\left(A_t^{-1}\right)'$ where $A_t$ is a lower triangular matrix and $\Lambda_t$ is a diagonal matrix. In a seminal paper, Primiceri2005 let the non-zero elements of $A_t$ as well as the log of the diagonal of $\Lambda_t$ evolve as driftless random walks, whereas Cogley2005 let $A_t$ be constant over time. Among the literature close to ours, Clark2011 followed Cogley2005 using a constant $A_t$ while Cimadomo2016 let it too vary over time. In order to model larger systems, Carriero2016 proposed the use of a common drifting volatility where $\Sigma_t=f_t\Sigma$ and $f_t$ is a scalar geometric random walk. Thus, a change in volatility for the entire $n$-dimensional system is achieved by scaling the constant covariance matrix $\Sigma$ by $f_t$. This line of modeling was later also pursued in the context of mixed-frequency VARs by Gotz2018,Ankargren2018. In related work, Carriero2019 propose a sampling procedure which greatly mitigates the computational burden in large VARs with $\Sigma_t=A^{-1}\Lambda_t\left(A^{-1}\right)'$ by exploiting the triangular structure of $A$.

We will focus next on two of the aforementioned stochastic volatility models: the common drifting volatility model used by Carriero2016,Gotz2018 and the model used by Cogley2005. The latter model uses $\Sigma_t=A^{-1}\Lambda_t\left(A^{-1}\right)'$ where the Cholesky factor $L_tL_t'=\Sigma_t$ is $L_t=A^{-1}\Lambda_t^{0.5}$. For the common drifting volatility model, we can let $\Sigma=P^{-1}D(P^{-1})'$ where $P$ is lower diagonal with ones on the main diagonal and $D$ is a diagonal matrix. If we let $D_t=f_tD$, the common drifting volatility model can be written as $\Sigma_t=P^{-1}D_tP^{-'}$, i.e. the same form as the Cogley2005 model. In the following, it is therefore, for our purposes, enough to deal with the $\Sigma_t=A^{-1}\Lambda_t \left(A^{-1}\right)'$ case as this encompasses the common drifting stochastic volatility model.

To transform the model into a form that enables the use of the univariate filtering approach, it is necessary to transform the monthly VAR model such that the errors are independent:

align[align omitted — 92 chars of source]

where $x_t^*=Ax_t$ and $\Pi_j^*=A\Pi_jA^{-1}$. The observation equation is similarly adjusted to account for the transformation:

align[align omitted — 86 chars of source]

where $y^*_t=Ay_t$ and $\Lambda^*=A\Lambda (I_p\otimes A^{-1})$. However, this seemingly innocuous transformation comes with one drawback: for most choices of aggregation matrices $\Lambda$, the rows of $\Lambda^*$ corresponding to quarterly variables will consist of only non-zero elements. The consequence of this is that $x_{q,t}^*$ depends concurrently on $x_{m, t}^* $---a relation that is not predetermined at time $t$ and hence cannot be moved to the exogenous component of the model. It would therefore be necessary to keep $x_{m,t}^*$---although not its lags---in the state vector of the compact form, thereby creating a state vector of dimension $(n+n_qp)\times1$. Differently put, for these stochastic volatility models, enabling the use of the univariate filtering procedure comes at the expense of requiring an additional $n_m$ terms to be put into the state vector. However, the adaptive algorithm proposed by Ankargren2019 is still applicable and may provide considerable computational improvements relative to the original Schorfheide2015 procedure as demonstrated by Figure (ref).

Empirical Application: Estimating Large Mixed-Frequency VARs for the US

In this section, we estimate the mixed-frequency VARs with factor stochastic volatility using the three sets of data discussed in Section (ref). We estimate the three models on the full data and discuss some of the key features.

Implementation Details

The sample we use for estimation starts in January 1980. We follow previous work that has used mixed-frequency VARs for US data, see e.g. Schorfheide2015,Gotz2018, and use $p=6$ lags. For the number of factors, KastnerHuber2018 found that a factor structure with two factors was preferable in terms of joint log predictive score for a selected subset of variables using their quarterly model with 215 variables, while one or no factor was preferable with respect to univariate density forecasting of GDP. Thus, since their results suggest zero to two factors for a substantially larger model than ours, we settle on a single $r=1$ factor in all three of our models. Each model is estimated using 30,000 MCMC draws, where the initial 10,000 are discarded for burn-in. Of the remaining 20,000, we save every 20$^{\text{th}}$ draw to improve mixing and reduce storage costs. It can be noted that the CCM-119 model draws $n(np+1)=1$85,085 regression parameters and $T(n+r)=55,320$ volatilities at every iteration and therefore requires a substantial amount of storage, which makes thinning a necessity. We leverage multiple cores to speed up our computations and draw regression parameters in parallel. Using an Intel Xeon E5-2630 v4 processor with 10 cores, roughly half a second per draw is needed for CCM-119 when estimated on the full sample. While computational resources, models and implementations are inevitably different, half a second per draw can be contrasted to the work by Carriero2019 who document roughly five seconds per draw. Therefore, by current standards the computational burden of our large mixed-frequency VAR is undoubtedly within the realms of acceptability.

Mixing of MCMC Samplers

Thinning the draws has, in addition to mitigating storage costs, the advantage of reducing the correlation between the draws that are stored, thereby improving the share of information contained in the draws saved. To evaluate the mixing of an MCMC sampler, a common metric is the inefficiency factor (see e.g. Chib2001). The inefficiency factor for the $i$th parameter is

align[align omitted — 51 chars of source]

where $\rho_{i,j}$ is the $j$th autocorrelation of the $i$th parameter. A value of one is obtained when there is no autocorrelation. Consequently, values close to one indicate low autocorrelation in the MCMC sampler and a large degree of information in the samples. Large values, on the other hand, imply that the chain contains little information due to high autocorrelation. A value below 20 is often used as a rule of thumb for good mixing Primiceri2005. We estimate the inefficiency factor using the coda package for R Plummer2006. Table (ref) provides a summary of the distributions of inefficiency factors by group of parameter.

table[table omitted — 2,774 chars of source]

Table (ref) shows that the inefficiency factors for CCM-20 and Koop-34 are generally well below or around 20, with the exception of a handful of outliers. The vast majority of parameters in the CCM-119 model also display satisfactory mixing, but the size of the outliers is somewhat larger. The most problematic aspect of the model is the mixing of the latent factor. It is encouraging, however, that mixing is excellent for the other parameters in the model. The approach we have discussed in this paper enables us to speed up computations and, for a fixed amount of time, produce more draws from the posterior distribution. We have noticed that mixing of the the latent high-frequency series can be troublesome if the number of draws is limited and so while other methods for large VARs can be extended to handle mixed-frequency data, the benefit of the route taken in this paper is that we can obtain draws of higher quality in the same amount of time.

Results

We now turn to the estimated models and illustrate some of their core features. To begin, we study the volatility of the latent factor. Figure (ref) plots the volatilities estimated by each model. To simplify the interpretation, the figure includes periods of contraction as defined by the National Bureau of Economic Research.

figure[figure omitted — 98,610 chars of source]

Figure (ref) shows a remarkable similarity between the factor volatilities estimated by Koop-34 and CCM-119, with major peaks around all recessions in the sample. The volatility estimated by CCM-20, however, is quite different as it estimates the volatility of the factor to be substantially higher in the beginning of the sample during the 1981--1982 recession. Nevertheless, the volatility in the CCM-20 model did show signals of increased volatility during the most recent two contractions. That volatility is heightened is more evident in the volatilities estimated by the Koop-34 and CCM-119 models, where all recessions coincide with peaks in volatility.

In order to better understand the meaning of the factor and its volatility, Figure (ref) displays boxplots of the posterior distribution of the factor loadings in the CCM-119 model.\footnote{Without further restrictions, the sign of the vector of loadings is not identified. We identify the sign of the vector of loadings a posteriori using the “maximin” method suggested by Kastner2017. The sign at every iteration is selected such that the loading that has the largest minimum of absolute values of draws is positive at every iteration.} From the figure, it is evident that the factor largely captures comovements in prices as most of the variables in this group are associated with large posterior loadings. Moreover, most of the variables in the output and income category load on the factor while the volatility of, e.g., labor market and stock market variables are mainly driven by their respective idiosyncratic components. We can also see that the posterior distribution of the GDP loading is wider than for other loadings. It is hardly surprising that we find a wider posterior distribution for the GDP loading as the factor, and hence the associated loadings, relate to the monthly frequency, at which the GDP series changes from iteration to iteration.

figure[figure omitted — 88,385 chars of source]

As a final illustration of the model's features, we next study the implied volatility of GDP. It follows from the construction of the model that the conditional variance of the error term in the equation for monthly GDP is

align[align omitted — 111 chars of source]

where $\lambda_{GDP, f}$ is the GDP loading. In order to improve interpretability and make the volatility more comparable to the observed outcomes of GDP, we compute and plot “aggregated” volatility where $\sigma^2_{t, GDP}$ is aggregated using the triangular aggregation scheme discussed in Section (ref). Figure (ref) shows the volatilities over time.

figure[figure omitted — 75,349 chars of source]