Simulation smoothing for nowcasting with large mixed-frequency VARs

TL;DR

This paper introduces two algorithms to improve the computational efficiency of simulation smoothing in large mixed-frequency VARs, enabling better nowcasting in high-dimensional macroeconomic models.

Contribution

The paper proposes an adaptive algorithm that enhances simulation smoothing efficiency for large mixed-frequency VARs, facilitating high-dimensional nowcasting applications.

Findings

Significant speed improvements with the adaptive algorithm

Enables high-dimensional mixed-frequency VAR estimation

Provides a key tool for macroeconomic nowcasting

Abstract

There is currently an increasing interest in large vector autoregressive (VAR) models. VARs are popular tools for macroeconomic forecasting and use of larger models has been demonstrated to often improve the forecasting ability compared to more traditional small-scale models. Mixed-frequency VARs deal with data sampled at different frequencies while remaining within the realms of VARs. Estimation of mixed-frequency VARs makes use of simulation smoothing, but using the standard procedure these models quickly become prohibitive in nowcasting situations as the size of the model grows. We propose two algorithms that alleviate the computational efficiency of the simulation smoothing algorithm. Our preferred choice is an adaptive algorithm, which augments the state vector as necessary to sample also monthly variables that are missing at the end of the sample. For large VARs, we find considerable improvements in speed using our adaptive algorithm. The algorithm therefore provides a crucial building block for bringing the mixed-frequency VARs to the high-dimensional regime.

Equations

$$x_t={\Pi }_c+{\Pi }_1{x}_{t-1}+\cdots +{\Pi }_p{x}_{t-p}+u_t,u_t\sim N(0,{\Sigma }_t)$$

$$z_t=\left(\begin{array}{cccc}{\Pi }_1& {\Pi }_2& \cdots & {\Pi }_p\\ \text{lx@intercol}\text{hfil}{I}_{n(p-1)}\text{hfil}\text{lx@intercol}& 0_{n(p-1)\times n}& & \end{array}\right)z_{t-1}+\left(\begin{array}{c}{\Pi }_c\\ 0_{n(p-1)\times 1}\end{array}\right)+\left(\begin{array}{c}{u}_t\\ 0_{n(p-1)\times 1}\end{array}\right)$$

$$y_t=S_t\Lambda z_t$$

$$p(\Pi ,\Sigma |X)$$

$$p(X|\Pi ,\Sigma ,Y).$$

$$X=\left(\begin{array}{cc}{X}_{m,1:T_b}& X_{q,1:T_b}\\ X_{m,T_b+1:T}& X_{q,T_b+1:T}\end{array}\right)$$

$$\left(\begin{array}{c}{y}_{m,t}\\ y_{q,t}\end{array}\right)$$

$$\left(\begin{array}{c}{x}_{q,t}\\ z_{q,t-1}\end{array}\right)$$

$${\Sigma }_t=W_t{W}_t^{\prime },W_t=\left(\begin{array}{c}{W}_t^{(1,\bullet )}\\ W_t^{(2,\bullet )}\end{array}\right)$$

$$y_t$$

$${\alpha }_t$$

$$v_t$$

$$K_t$$

$$a_{t|t}$$

$$a_{t+1}$$

$$a_{t|T}$$

$$r_{t-1}$$

$$y_t$$

$${\alpha }_t$$

$$Z_t=S_t\left(\begin{array}{cc}{I}_{n_m}& 0_{n_m\times (np+n_q)}\\ 0_{n_m\times n_m}& {\Lambda }_q\end{array}\right),$$

$$T_t=\left(\begin{array}{cccc}{\Pi }_1& \cdots & {\Pi }_p& 0_{n\times n}\\ \text{lx@intercol}\text{hfil}{I}_{np}\text{hfil}\text{lx@intercol}& 0_{n\times n}& 0_{n\times n}& \end{array}\right),$$

$$F_t$$

$$K_t$$

$$L_t$$

$$M_t$$

$$a_{t+1}$$

$$P_{t+1}$$

$$r_{t-1}$$

$$\left(\begin{array}{c}{y}_{m,t}\\ y_{q,t}\end{array}\right)$$

$$\left(\begin{array}{c}{x}_{{\mathcal{U}}_t,t}\\ x_{q,t}\\ ⋮\\ x_{{\mathcal{U}}_t,t-p}\\ x_{q,t-p}\end{array}\right)$$

$$J_t^{\prime }\left(\begin{array}{c}{x}_{{\mathcal{U}}_t,t}\\ x_{q,t}\end{array}\right)=\left(\begin{array}{c}{x}_{{\mathcal{U}}_{t-1},t}\\ x_{q,t}\end{array}\right).$$

$$T_t=\left(\begin{array}{cccccc}{\Pi }_1^{{\mathcal{U}}_t,{\mathcal{U}}_{t-1}}& {\Pi }_1^{{\mathcal{U}}_t,q}& \cdots & {\Pi }_p^{{\mathcal{U}}_t,{\mathcal{U}}_{t-1}}& {\Pi }_p^{{\mathcal{U}}_t,q}& 0_{|{\mathcal{U}}_t|\times |{\mathcal{U}}_{t-1}|}\\ \text{lx@intercol}\text{hfil}{I}_p\otimes J_t\text{hfil}\text{lx@intercol}& 0_{n_q\times |{\mathcal{U}}_{t-1}|}& & & & \end{array}\right),$$

$$D_t=\left(\begin{array}{ccc}{\Pi }_1^{{\mathcal{U}}_t,{\mathcal{O}}_{t-1}}& \cdots & {\Pi }_p^{{\mathcal{U}}_t,{\mathcal{O}}_{t-1}}\\ {\Pi }_1^{q,{\mathcal{O}}_{t-1}}& \cdots & {\Pi }_1^{q,{\mathcal{O}}_{t-1}}\\ \text{lx@intercol}\text{hfil}{I}_p\otimes J_{\perp ,t}\text{hfil}\text{lx@intercol}& & \end{array}\right).$$

$$Z_t=\left(\begin{array}{cccccc}{0}_{n_m\times (|{\mathcal{U}}_t|+n_q)}& {\Pi }_1^{{\mathcal{O}}_t,{\mathcal{U}}_t}& {\Pi }_1^{{\mathcal{O}}_t,q}& \cdots & {\Pi }_p^{{\mathcal{O}}_t,{\mathcal{U}}_t}& {\Pi }_p^{{\mathcal{O}}_t,q}\\ \text{lx@intercol}\text{hfil}{S}_{q,t}{\Lambda }_{qq}[I_{p+1}\otimes I_{-m,|{\mathcal{U}}_t|+n_q}]\text{hfil}\text{lx@intercol}& & & & & \end{array}\right)$$

$$C_t=\left(\begin{array}{ccc}{\Pi }_1^{{\mathcal{O}}_t,{\mathcal{O}}_{t-1}}& \cdots & {\Pi }_p^{{\mathcal{O}}_t,{\mathcal{O}}_{t-1}}\\ \text{lx@intercol}\text{hfil}{S}_{q,t}{0}_{n_q\times p|{\mathcal{O}}_{t-1}|}\text{hfil}\text{lx@intercol}& & \end{array}\right)$$

Full text

1 Introduction

Macroeconomic forecasting is a central activity for many organizations. Central banks, treasuries, and other government institutions as well as non-governmental organizations such as banks and corporations rely on forecasts in order to assess the future development of the local and global economy. Accurate forecasts is therefore a necessity for making informed decisions that minimize adversarial effects and maximize the utility of the decisions that are to be made. There is a vast array of forecasting methods that can be employed for macroeconomic data, but one of the more prominent and popular tools used in macroeconomic forecasting is the vector autoregressive (VAR) model. It is a multivariate time series model relating current values of a vector of variables to lags of all of the variables in the model. In contrast to structural models relying heavily on economic theory, the VAR is a relatively agnostic model that is easy to use.

Traditional vector autoregressive models make use of data that is sampled at the same frequency. Most often, VAR models used for forecasting or structural analysis are estimated on data with exclusively monthly or quarterly variables. In many cases, the inclusion of gross domestic product (GDP) determines whether or not the model is estimated on a quarterly basis. Mixed-frequency methods, in general, are attempts that aim to avoid having to aggregate, say, monthly data to the quarterly frequency when at least one of the included variables is sampled quarterly. One advantage is that the loss of information induced by the aggregation can be avoided, and information can be incorporated into the forecasting model in a more timely manner. Among the more prominent methods are dynamic factor models as developed by Mariano and Murasawa, (2003, 2010); Marcellino et al., (2016) and mixed-data sampling (MIDAS) as introduced by Ghysels et al., (2006) and Ghysels et al., (2007). Related literature has also considered vector autoregressive moving average (VARMA) models in the mixed-frequency setting, see Zadrozny, (2016), Anderson et al., (2016) and Deistler et al., (2017). The relevance of the inclusion of an MA component has also been studied for MIDAS regressions by Foroni et al., (2019).

The MIDAS regression approach is generalized to a VAR setting with a Bayesian approach by Ghysels, (2016), and McCracken et al., (2015) propose a high-dimensional variant of the MIDAS-VAR. Götz et al., (2016) exploit the model by Ghysels, (2016) for testing Granger causality using mixed-frequency data. They consider both reduced rank regressions and a Bayesian approach in situations with a large number of high-frequency observations per low-frequency period.

The approach considered in this paper is the mixed-frequency VAR using a state-space formulation. The primary empirical purpose of the model that motivates our work is the use of a large monthly data set for forecasting GDP growth. We focus on the methodology developed by Schorfheide and Song, (2015), who use Gibbs sampling and simulation smoothing to estimate the model in the spirit of Carter and Kohn, (1994). Other developments of the mixed-frequency VAR include Eraker et al., (2015), who instead of simulation smoothing employ a single-move sampler. Qian, (2016) develop an alternative sampling strategy based on linear transformations of the data. Cimadomo and D’Agostino, (2016) employ a mixed-frequency VAR with time-varying parameters and stochastic volatility to cope with a change in frequency, and Götz and Hauzenberger, (2018) extend the Schorfheide and Song, (2015) model to include stochastic volatility and time-varying intercepts.

The focus of this paper is on developing the simulation smoother proposed by Schorfheide and Song, (2015) in a way that enhances its performance in nowcasting situations. The merits of the mixed-frequency VARs in nowcasting and forecasting have been demonstrated by e.g. Schorfheide and Song, (2015); Eraker et al., (2015); Ankargren et al., (2018); Götz and Hauzenberger, (2018). These papers frame the problem as a latent variable problem, where the low-frequency series have an underlying, high-frequency series which we do not observe. The model thus relates the observed outcome to the unobserved latent variable in a measurement equation, while the VAR itself is specified at the high frequency.

One of the key advantages of the mixed-frequency VAR is that it is, in essence, a standard VAR with an additional block in the Markov Chain Monte Carlo (MCMC) sampler. This additional block samples the unobserved latent variable from the conditional posterior, given the observed data and the parameters, by means of simulation smoothing. Next, given the unobserved latent variable, the remaining steps needed to estimate the model are the usual steps drawing from the posterior of the parameters. The important difference is that the conditional posteriors of the parameters here condition upon the draw of the latent variable, which makes the problem standard. Because our focus in this paper is on the simulation smoother used for sampling the latent variables, we do not dwell any further on the other blocks of the MCMC sampler.

Our contribution is that we develop a simulation smoother based on the proposal by Schorfheide and Song, (2015) that better deals with ragged edges in the (mixed-frequency) data, i.e. when the end of the sample is unbalanced. The motivation is that in real-time forecasting situations, the information set that the forecaster possesses can vary from day to day: new data points for the different series are not obtained at the same time, but on different days. In a stylized data set, the forecaster knows at the beginning of April the outcomes of financial variables for March, but the numbers for other variables have not yet been released. In addition, some variables are subject to larger publication delays. At the beginning of April, the forecaster may therefore not know the February outcomes for slow variables. We demonstrate that in such situations, use of our algorithm is a necessity when the dimension of the model is large. On the other hand, if the data is not unbalanced, our algorithm reduces to the Schorfheide and Song, (2015) procedure. There is, in other words, no loss associated with using our procedure. Moreover, our algorithm is purely a computational improvement. The results from the different algorithms are therefore identical given the same random number generator.

For small-scale models, the simulation smoothing algorithm that Schorfheide and Song, (2015) proposed works well. However, there is a growing literature on large VAR models and their usefulness for forecasting, see for example Bańbura et al., (2010); Koop, (2013); Carriero et al., (2019). When the dimension of the model increases—either by the number of variables or the lag length—the simulation smoother becomes prohibitive in estimating the models. The method we propose for dealing with the ragged edge nature of data is intuitive: we sample only those variables which we must sample. The issue with the Schorfheide and Song, (2015) sampler is that, for the ragged edge part of the data, it wastes a substantial amount of computational resources on sampling what is already known, i.e. fully observed monthly variables. We circumvent this problem by developing an adaptive simulation smoother, which augments the state vector with the missing variables only. This allows us to improve the computational efficiency for large-dimensional models. Our method applies to any mixed-frequency VAR estimated using data with ragged edges that employs simulation smoothing.

The rest of the paper is organized as follows. Section 2 describes the mixed-frequency VAR and the approach suggested by Schorfheide and Song, (2015) in more detail, Section 3 describes our two new algorithms amd Section 4 provides some computational experiments illustrating the gains of our algorithms. Section 5 concludes.

2 Mixed-frequency VARs

To model data observed at mixed frequencies we follow Schorfheide and Song, (2015); Ankargren et al., (2018); Götz and Hauzenberger, (2018) and deal with the problem by postulating a VAR model at the high frequency. That is, we start with the VAR($p$) given by

[TABLE]

where $x_{t}$, $\Pi_{c}$ and $u_{t}$ are $n\times 1$ and $\Pi_{j}$ are $n\times n$. The vector $x_{t}$ can be decomposed into monthly and quarterly variables as $x_{t}=(x_{m,t}^{\prime},x_{q,t}^{\prime})^{\prime}$ of dimensions $n_{m}\times 1$ and $n_{q}\times 1$, respectively, where $n=n_{m}+n_{q}$. The model can be formulated in terms of the companion form as

[TABLE]

where $z_{t}=(x_{t}^{\prime},\dots,x_{t-p+1}^{\prime})^{\prime}$. For future use, let also $\Pi=(\Pi_{c},\Pi_{1},\dots,\Pi_{p})$.

The empirical problem motivating our work is that in reality we only observe

[TABLE]

where $S_{t}$ is a deterministic $n_{t}\times n$ selection matrix and $\Lambda$ an $n\times pn$ aggregation matrix. For ease of exposition, we assume that the aggregation employed includes at most $p-1$ previous values.

The posterior distribution of interest is $p(\Pi,\Sigma,X|Y)$. The key feature of Schorfheide and Song, (2015) is the use of the forward-filtering and backwards-smoothing algorithm (Carter and Kohn, , 1994; Frühwirth-Schnatter, , 1994) in a Gibbs sampler. Such an approach utilizes the fact that $(\Pi,\Sigma)$ is conditionally independent of $Y$ given $X$. Hence, the basic Gibbs sampling strategy for the mixed-frequency VAR is to sample iteratively from:

[TABLE]

Note that conditional on $X$, the first step is standard in the literature. Intuitively, the full method can be understood from the perspective of missing values. The $p(X|\Pi,\Sigma,Y)$ step produces a draw of the underlying monthly series of all variables observed quarterly. This draw is produced under the restriction that when aggregated according to the linear aggregation defined by $\Lambda$, the aggregate is equal to the observed outcomes. However, conditional on a particular draw of $X$, draws of $\Pi$ and $\Sigma$ can be produced as if this were the real data, leading to standard procedures.

The main computational obstacle in an off-the-shelf implementation of a mixed-frequency model employing the companion form formulation is to draw from the posterior of the underlying variable, $p(X|\Pi,\Sigma,Y)$. Two ways of circumventing the high dimensionality are to use the computationally convenient reformulation of the model suggested by Schorfheide and Song, (2015) or the method proposed by Qian, (2016). The companion form presented in the beginning of the section should generally be avoided in computations because of the extraneous computational burden it induces. We review next the procedure developed by Schorfheide and Song, (2015), which largely avoids the companion form in the computations. Because of our frequent referral to their procedure, we will refer to it as SS15.

The companion form expresses the model in such a way that all of the endogenous variables in the VAR enter the state equation of the state-space model. Furthermore, lagged values of the variables are also included, which yields a state vector of dimension $np\times 1$. In empirical applications dealing with monthly and quarterly data, it is usually the case that the set of monthly variables is fully observed up through a time point $T_{b}$.

The SS15 method exploits the balanced nature of the monthly part of the sample for $t=1,\dots,T_{b}$. Intuitively, the primary purpose of employing a state-space model formulation and its associated simulation smoother is that it provides a tractable way of imputing the missing observations of the quarterly variables through a data augmentation approach common in Bayesian statistics (see Tanner and Wong, , 1987). The simulation smoothing sampling step consists of sampling from $p(X|\Pi,\Sigma,Y)$, but since

[TABLE]

it is evident that including the observed monthly variables in the state vector means that we are also drawing from $p(X_{m,1:T_{b}}|\Pi,\Sigma,Y)$ — a degenerate distribution with all of its mass on $Y_{m,1:T_{b}}$. To alleviate the computational cost, SS15 therefore formulate the model differently for $t=1,\dots,T_{b}$, where the $X_{m,1:T_{b}}$ is moved out of the state equation and treated as exogenous and the dimension of the state vector is thereby reduced from $np\times 1$ to $n_{q}(p+1)\times 1$. The new formulation of the model does not alter the underlying statistical model, but is purely computational. Throughout, we will refer to the following as the compact form:

[TABLE]

where $e_{t}=\begin{pmatrix}e_{m,t}^{\prime}&e_{q,t}^{\prime}\end{pmatrix}^{\prime}$ is independently and identically distributed as $N(0,I_{n})$ and

[TABLE]

and the partitioning in $W_{t}$ consists of blocks of size $n_{m}\times n$ and $n_{q}\times n$. Furthermore, the matrices $\Pi_{ij}$ for $i,j\in\{m,q\}$ represent submatrices of $\Pi$ consisting of parameters relating the set of $j$ variables to the $i$ variable; i.e. $\Pi_{mq}$ contains the parameters pertaining to the lagged effects of the quarterly variables in the equations for the monthly variables. In the preceding description of the model we explicitly allow for a time-varying error covariance matrix, but assume constant regression parameters. Relaxing the assumption of constant regression parameters, as Cimadomo and D’Agostino, (2016); Götz and Hauzenberger, (2018) have done for mixed-frequency VARs, can easily be done within the algorithms we develop.

To simplify the description of the filtering and smoothing algorithm, we write the model in more general notation as

[TABLE]

The full SS15 procedure for producing a draw from the conditional posterior uses the Durbin and Koopman, (2002) simulation smoothing technique, which consists of the steps:222The constants in the model need to be excluded either in the simulation step, or in the filtering/smoothing steps, see Jarociński, (2015).

1. Generate a pseudo sample $y_{t}^{+}$ and $\alpha_{t}^{+}$ based on the recursions in (1). For the difference $y_{t}^{*}=y_{t}-y_{t}^{+}$:

Filter for $t=1,\dots,T_{b}$ using the compact form 3. 2.

Move to the companion form and filter for $t=T_{b}+1,\dots,T$ 4. 3.

Smooth for $t=T,\dots,T_{b}+1$ using the companion form 5. 4.

Move back to the compact form and smooth for $t=T_{b},\dots,1$ 6. Compute the draw as $x_{t}=\hat{x}_{t}^{*}+x_{t}^{+}$ where $\hat{x}_{t}^{*}$ is the smoothed series from step 4.

Because $n_{q}(p+1)\ll np$ in most applications, the compact procedure typically offers a substantial improvement in speed.

The key parts of the simulation smoothing algorithm are the filtering and smoothing steps. In the following when discussing the simulation smoothing step, we will for simplicity discuss filtering and smoothing on their own, but note that they should be interpreted as the main building blocks for the full simulation smoothing algorithm. The filtering recursions are, for $t=1,\dots,T_{b}$:

[TABLE]

Having transitioned back to the compact form, the algorithm is completed by the following smoothing recursions for $t=T_{b},\dots,1$:

[TABLE]

The details for the transitions between the compact and companion forms can be found in A.

3 Handling mixed frequencies in large-dimensional VARs

Unsurprisingly, if $n_{q}$ is fixed and $n$ grows an increasingly large share of the total computational time of the SS15 procedure can be attributed to the companion form, even if the companion form is employed for only one or two time points. The main bottleneck is the computation of the variance of the one-step ahead prediction. The variance is computed by multiplying three square matrices of dimensions $n(p+1)\times n(p+1)$, an operation which is prohibitive even if it needs to performed only once or twice if, say, $n=100$ and $p=4$.

In the next two subsections, we describe two ways to increase the computational efficiency of filtering and smoothing for the unbalanced part of the sample. The first is blocked filtering, in which the structure of the system matrices in the state-space model is exploited to reduce large matrix multiplications to a larger number of smaller block operations. The second approach is adaptive filtering, which adaptively augments the state vector with the missing monthly variables and by doing so avoids the use of a full-blown companion form altogether. It should be emphasized that the two alternative procedures that we propose produce draws which, given the same seed, are numerically identical to draws made by the algorithm described in the previous section.

For use in the following, we formulate the state-space model for the companion form as:333When the companion form is used for the unbalanced part of the data following use of the compact form, it is formulated with $p+1$ lags to aid in the transition between the model formulations.

[TABLE]

where

[TABLE]

For simplicity, we abstract here from including the intercept in the following, but inclusion thereof is straightforward.

3.1 Blocked filtering

While the bottleneck of the Kalman filter lies in the multiplication of three potentially large matrices, improvements can be made throughout the entire filtering algorithm by exploiting the fact that many of the system matrices in general are not dense, but have inherent structures—including blocks of zeros and identity matrices. Strid and Walentin, (2009) successfully used the same idea to reduce the computational burden associated with a large-scale dynamic stochastic general equilibrium (DSGE) model. Inspired by their work for DSGE models, it is clear from Equations (15)–(16) that there is a large degree of structure in the recursions for the companion form. Hence, it is likely subject to a speed improvement if the structure is efficiently utilized. In addition to what (15)–(16) display, the matrix $\Lambda_{q}$ has a sparse structure where at most $n_{q}p$ columns (out of $np+n_{q}$) contain non-zero elements, the precise number depending on the specific aggregation scheme used. Let $p_{q}$ denote the number of lags included in the aggregation scheme. If the intra-quarterly average is used for aggregation, then $p_{q}=3$ and $3n_{q}$ columns contain non-zero values with a maximum of one non-zero value per column. Thus, when the system matrices become larger, substantial gains in computational efficiency can be obtained by simply exploiting these structures.

Let $\Lambda_{qq}$ be the $n_{q}\times n_{q}p_{q}$ submatrix of $\Lambda_{q}$ obtained by deleting zero-columns in $\Lambda_{q}$, i.e. columns corresponding to monthly variables in $\alpha_{t}$ and additional lags not included in the aggregation. Let also $a_{t}$ (and $a_{t|t}$) be partitioned into $p+1$ blocks of $n$ values, i.e. $a_{t}=(a_{t}^{1^{\prime}},a_{t}^{2^{\prime}},\dots,a_{t}^{p+1^{\prime}})^{\prime}$. We present next the filtering algorithm for the companion form where the use of the structure is more explicit. In doing so, we let square brackets indicate that the enclosed quantity is computed before it is used to facilitate multiple uses efficiently.

First, we can utilize the sparsity of $Z_{t}$ when computing $F_{t}$:

[TABLE]

Here, $P_{t}^{mm}$ refers to the upper left $n_{m}\times n_{m}$ block of $P_{t}$ containing the variances of the one-step ahead predictions for the monthly variables. Similarly, $P_{t}^{qq}$ denotes the $p_{q}n_{q}\times p_{q}n_{q}$ submatrix of $P_{t}$ formed by extracting the elements corresponding to the first $p_{q}$ lags of all quarterly variables, and $P_{t}^{mq}$ is the $n_{m}\times p_{q}n_{q}$ submatrix of covariances between $x_{m,t}$ and $x_{q,t-i}$ for $i=0,\dots,p_{q}-1$.

Next, compute $P_{t}^{1:p,m}(F_{t}^{-1})_{m,\bullet}$, where $P_{t}^{1:p,m}$ is the $np\times n_{m}$ submatrix retrieved from the upper left part of $P_{t}$, and $(F_{t}^{-1})_{m,\bullet}$ is the first $n_{m}$ rows of $F_{t}^{-1}$. Analogously, we compute also $P_{t}^{1:p,q}(F_{t}^{-1})_{q,\bullet}$ where $P_{t}^{1:p,q}$ is the $np\times n_{q}$ matrix obtained from taking the first $np$ rows of $P_{t}$ and the $n_{q}$ columns corresponding to $x_{q,t}$, and $(F_{t}^{-1})_{q,\bullet}$ is the bottom $n_{q}$ rows of $F_{t}^{-1}$. We can then construct:

[TABLE]

Let the left block of columns in $K_{t}$ be denoted by $K_{t}^{m}$ and the right by $K_{t}^{q}$ and compute

[TABLE]

and

[TABLE]

The state prediction is

[TABLE]

Finally, a sizable improvement can be obtained for the variance of the state prediction by first computing $\Pi P_{t|t}^{1:p,1:p}$ and then constructing

[TABLE]

There are two main points to make from the recursions for blocked filtering. First, there are several instances where the same quantity appears in multiple places. It is therefore unnecessary to compute it over and over. Instead, efficiency can be gained by computing it once and copying the object to the appropriate locations. Second, because of the large number of zeros as well as the identity matrices in the system matrices, several computations can be avoided altogether in favor of subsetting the appropriate matrices. The two approaches can be used in conjunction in some cases. For example, to compute $K_{t}Z_{t}$ one would first select the first $n_{m}$ columns of $K_{t}$. To compute the second part of $K_{t}Z_{t}$, i.e. $K_{t}^{q}\Lambda_{q}$, one can compute $K_{t}^{q}\Lambda_{qq}$ and subtract the column vectors of the obtained matrix from the corresponding columns of $T_{t}$; because $\Lambda_{q}$ typically contains a large number of zero columns, so will $K_{t}^{q}\Lambda_{q}$ and these need not be computed. The same type of computation can also be carried out for $v_{t}$, which involves $Z_{t}a_{t}$ for which copying of pre-computed data and multiplication for a subset improves efficiency.

The smoothing step is cheap in comparison with the filtering part, but is nevertheless subject to a modest improvement through:

[TABLE]

where $[F_{t}^{-1}v_{t}]_{m,\bullet}$ and $[F_{t}^{-1}v_{t}]_{q,\bullet}$ are the monthly and quarterly, respectively, blocks of rows of $F_{t}^{-1}v_{t}$, which is computed during filtering.

3.2 Adaptive filtering

The advantage of blocked filtering is its ability to improve upon the existing SS15 compact filtering procedure with relatively little effort. However, the downside is that it improves on the problem without attempting to fix it.

Figure 1 offers an illustrative visual understanding of the compact filtering procedure: for $t\leq T_{b}$, the monthly block is fully observed, but for $t>T_{b}$ some of the monthly variables are missing. In this particular data set, based on the January vintage of the FRED database (McCracken and Ng, , 2016), at $t=T_{b}+1$ only a single monthly indicator is missing. Nevertheless, in the compact procedure a move to the companion form is still triggered. Using the companion form, the 18 monthly variables which are observed at $t=T_{b}+1$ will have their true values imputed with themselves. For $t=T_{b}+2$, a larger number of monthly variables are unobserved. Yet, use of the full companion form would still need to carry out additional computations for the sole purpose of imputing six variables with their already-known outcomes.

Motivated by the large degree of unnecessary computations performed in the companion form for larger models, we therefore develop an adaptive procedure. Intuitively, the adaptive filtering procedure is simple and appealing: for every $t=T_{b}+1,\dots,T$, augment the state vector with the monthly variables missing at that time period. In relation to Figure 1, the procedure consists of adding the single missing monthly variable to the state vector at $t=T_{b}+1$ and at $t=T_{b}+2$ include all but the first six monthly variables in the state vector.

While the adaptive procedure is intuitive and simple in principle, it warrants some additional notation to be introduced. Let $\mathcal{U}_{t}$ denote the set of indexes of the monthly set of variables which are unobserved at time $t$. Correspondingly, let $\mathcal{O}_{t}$ denote the set of indexes of the observed set of monthly variables at time $t$. For $t=1,\dots,T_{b}$, we have that $y_{m,t}=y_{\mathcal{O}_{t},t}$ is of dimension $n_{m}\times 1$. For $t=T_{b}+1,\dots,T$, it still holds that $y_{m,t}=y_{\mathcal{O}_{t},t}$, but now the dimension is $|\mathcal{O}_{t}|<n_{m}$. The notation $y_{\mathcal{O}_{t},t-\tau}$ is taken to mean the observed values at time $t-\tau$ of the monthly variables observed at time $t$. Similarly, $x_{\mathcal{U}_{t},t-\tau}$ refers to the latent process at time $t-\tau$ for the set of unobserved monthly variables that are unobserved at time $t$.

Equipped with this new set of notation, we can formulate the state-space model as

[TABLE]

The preceding model formulation is valid for any $t=1,\dots,T$, with no particular distinction between the balanced and unbalanced part. However, because $\mathcal{U}_{t}=\varnothing$ for $t=1,\dots,T_{b}$, the model is mostly constant and equal to the compact formulation. The SS15 compact-companion procedure fits into the above formulation by specifying $\mathcal{U}_{t}=\varnothing$ and $\mathcal{O}_{t-1}=\{1,\dots,n_{m}\}$ for $t=1,\dots,T_{b}$, and $\mathcal{U}_{t}=\{1,\dots,n_{m}\}$ and $\mathcal{O}_{t-1}=\varnothing$ for $t=T_{b}+1,\dots,T$. One particular feature is that if $\mathcal{U}_{t}\neq\mathcal{U}_{t-1}$, $T_{t}$ will not be a square matrix. It is this time-varying size and shape that automatically lets the model augment new variables with missing observations to the state vector.

Let $J_{t}$ be a $(|\mathcal{U}_{t}|+n_{q})\times(|\mathcal{U}_{t-1}|+n_{q})$ selection matrix with the property that

[TABLE]

We can then construct the $T_{t}$ matrix as:

[TABLE]

where $\Pi_{j}^{\mathcal{U}_{t},\mathcal{U}_{t-1}}$ refers to the elements in the $\mathcal{U}_{t-1}$ columns and $\mathcal{U}_{t}$ rows of $\Pi_{j}$. In practice, there is no need to compute the Kronecker product as the bottom-left block of $T_{t}$ can more efficiently be filled programmatically. Depending on the observational pattern of the data, it may be worthwhile to use a block approach when dealing with multiplications involving $T_{t}$ because of the sparse structure of $I_{p}\otimes J_{t}$. However, if $n_{q}$ is low—our interest lies primarily in $n_{q}=1$—and if the unbalanced part contains many observations, there is likely little to gain.

The selection matrix $J_{t}$ is the $|\mathcal{U}_{t}|+n_{q}$ identity matrix with columns corresponding to $\mathcal{U}_{t}\bigcap\mathcal{O}_{t-1}$ deleted, i.e. columns for variables which were observed at $t-1$ but are unobserved at $t$. The matrix formed by collecting the deleted columns constitutes the orthogonal complement of $J_{t}$, denoted by $J_{\perp,t}$. The $J_{\perp,t}$ matrix is used in constructing $D_{t}$ to ensure that variables observed at $t-1$ but unobserved at $t$ get their lag identities—i.e. the bottom equations of the state equation—correct:

[TABLE]

For the observation equation, let $I_{-m,|\mathcal{U}_{t}|+n_{q}}$ be the $|\mathcal{U}_{t}|+n_{q}$ identity matrix with the first $|\mathcal{U}_{t}|$ rows deleted, so that $I_{-m,|\mathcal{U}_{t}|+n_{q}}\begin{pmatrix}x_{\mathcal{U}_{t},t}\\ x_{q,t}\end{pmatrix}=x_{q,t}$. Then we can construct the $Z_{t}$ system matrix as:

[TABLE]

[TABLE]

The only pieces missing for concluding the description of the adaptive procedure are the $H_{t}$ and $G_{t}$ matrices, which consist of:

[TABLE]

where the superscripts refer to all columns and the subsets of rows corresponding to: the observed monthly variables ($\mathcal{O}_{t}$), the unobserved monthly variables ($\mathcal{U}_{t}$) or the quarterly variables ($q$).

It should be noted that while the system matrices can be characterized using zero, identity and selection matrices, it is usually inefficient to use these forms explicitly. Faster handling of the matrices and the operations involved can be achieved by subsetting and copying previously computed objects appropriately.

The procedure itself is simple and straightforward, but the description of the system matrices is somewhat involved notationally. To illustrate how the model adapts to the structure of observed data, let us demonstrate it using a simple example. Figure 2 displays the observational pattern for our hypothetical data. The system matrices are not displayed here for space considerations, but can be found in full in Appendix B.

As before, the compact form is used for the balanced part of the sample as $\mathcal{U}_{t}=\varnothing$ and $\mathcal{O}_{t-1}=\{1,2,3\}$ for $t=1,\dots,T_{b}$. We use the short form $y_{1:j,t-1:t-3}$ to denote the vector $y_{1:j,t-1:t-3}=(y_{1,t-1},\dots,y_{j,t-1},y_{1,t-2},\dots,y_{j,t-3})^{\prime}$. For the current example, the observation and state equations for a model with three lags are:

[TABLE]

At $t=T_{b}+1$, the third monthly variable is missing so that $\mathcal{U}_{T_{b}+1}=\{3\}$ and $\mathcal{O}_{T_{b}}=\{1,2,3\}$. The third monthly variable is thus added to the state equation:

[TABLE]

At the end of the sample ($t=T_{b}+2=T)$, we obtain $\mathcal{U}_{T_{b}+2}=\{2,3\}$ and $\mathcal{O}_{T_{b}+1}=\{1,2\}$. The second monthly variable is now added to the state equation:

[TABLE]

In principle, the procedure thus simply amounts to continuing to perform the compact filtering procedure, but adjusting which variables are included as exogenous based on the observational pattern. The filtering and smoothing recursions are therefore the same as for the compact procedure in (11)–(12); what changes is instead the composition in the model and the construction of the system matrices.

By studying the recursions in more detail, the benefit of the adaptive procedure becomes clear. Take the main bottleneck—the computation of the variance of the one-step ahead prediction—as an example. Using the companion form with blocked filtering for the observational pattern in Figure 1 with $p=3$ requires multiplying $\Pi\,(20\times 80)$, $P_{t|t}^{1:p,1:p}\,(80\times 80)$ and the transpose $\Pi^{\prime}$ together. The number of scalar multiplications required to compute the full product is $20^{3}\times 80^{3}=4,096,000,000$.444The calculated number of scalar multiplications does not exploit the symmetric nature of the resulting matrix. Using the adaptive procedure, it is (for $P_{T_{b}+2}$) instead necessary to multiply $T_{T_{b}+2}\,(14\times 8)$, $P_{T_{b}+1|T_{b}+1}\,(8\times 8)$ and the transpose $T_{T_{b}+2}^{\prime}$ together. The number of scalar multiplications for this latter operation to be executed is $14^{3}\times 8^{3}=1,404,928$—almost four orders of magnitude smaller. There is, in other words, large possible gains to be achieved from adaptively adjusting the state vector to include only the bare minimum. The relative gain depends on the number of observed variables in the unbalanced part and the length of the unbalanced part.

4 Computational experiments

In this section, we demonstrate the computational benefits of our proposed algorithms. In designing the setting for which we conduct the computational experiments, we aim to stay close to a realistic scenario. Our primary empirical interest lies in using variables sampled monthly in forecasting a single quarterly variable, the canonical example in mind being use of a large set of monthly indicators for forecasting GDP growth. In order to see also the effects of including more quarterly variables (cf. Schorfheide and Song, , 2015; Götz and Hauzenberger, , 2018), we also consider $n_{q}=3$. For the total number of variables, we consider a grid of $n=10,20,\dots,120$ to capture the computational gains in both small and large models. The sample size is $T=500$ and we let the number of monthly variables with missing values at $t=499,500$ be: one for $n=10,\dots,40$, two for $n=50,\dots,80$ and three for $n=90,\dots,120$. The number of fully observed monthly series is set to $0.3n$ (rounded up), with the remaining containing missingness at $t=500$ only. The composition of the ragged edge is roughly in line with the January 2019 vintage of the FRED-MD database (McCracken and Ng, , 2016). We first let $p=6$ and vary $n$, and then consider the reverse situation. We implement the simulation smoothing routines in C++ using the Armadillo library for linear algebra (Sanderson and Curtin, , 2016) via R and RcppArmadillo (R Core Team, , 2019; Eddelbuettel and Sanderson, , 2014).

Figure 3 displays the computational cost per iteration. At the top in the figure, we find the standard SS15 compact procedure followed by, in order, the blocked and adaptive approaches. It is evident from the figure that for low-dimensional models, there is no need to care about the computational efficiency as the cost per iteration is low to begin with. Nevertheless, when the dimension of the model is large the choice of algorithm starts to matter. There is also little relative difference between using one or three quarterly variables in the model when the model size is large, which means that in large models there is no need to restrict oneself from adding more quarterly variables to the model for computational reasons. For a model using a data set similar in size to Bańbura et al., (2010) and Carriero et al., (2019), use of the adaptive method yields 5,000 draws in approximately one hour, whereas the SS15 procedure requires over 3 hours. While an hour is still needed, avoiding every bottleneck possible is highly desirable in a real-time forecasting situation where the model may be subject to frequent re-estimation.

Because the three methods all make use of the compact filtering procedure for the balanced part, the difference in computational cost is due to the different treatments of the unbalanced part. Thus, the figure clearly shows that using the companion form—even for two time periods only—can be costly, as substantial reductions can be obtained by handling it more carefully. In fact, using the adaptive procedure, dealing with the unbalanced part constitutes only a small part of the overall computational cost. Consequently, for further speed improvements one needs to focus on filtering and smoothing in the compact form. Further enhancements of the compact form may be possible depending on the remaining parts of the model, if e.g. there are structures imposed on the parameters that can be exploited.

Figure 4 displays the cost per iteration as a function of the number of lags $p$. The number of variables is $n=20$ (upper panel) and $n=120$ (bottom panel). Note that these two model sizes correspond to the models used by Bańbura et al., (2010); Carriero et al., (2019). In both of the aforementioned papers, the authors use $p=13$. The figure shows that as the number of lags increases, the added computational cost is lower using the adaptive procedure than if the blocked approach is used. Going from 3 to 12 lags roughly means an increase from 0.8 to 1.1 seconds using the adaptive procedure. In contrast, the cost of the SS15 procedure increases from 1.1 seconds to 14.1 seconds. Carriero et al., (2019) report that producing 5,000 posterior draws in their large model takes about 5.5 hours; adding a quarterly variable to their model would increase the overall computational time to approximately 7 hours using the adaptive procedure. If instead SS15 is used, estimating the model would take approximately 24 hours.

To more clearly see where gains of the adaptive algorithm are larger, Figure 5 presents the relative cost of the adaptive algorithm relative to the SS15 algorithm. It more explicitly details the pattern that previous figures have documented: the adaptive algorithm scales better with $p$, and it is primarily for large-dimensional models that the adaptive algorithm excels and provides substantial improvements. The precise level of the gains also depends on the observational pattern in the data. While the SS15 and the blocked procedures are unaffected by the specific observational pattern (given $T_{b}$ and $T$), the adaptive procedure is more apt the more variables are observed. Consequently, situations in which the data set includes a large number of financial variables that are observed immediately are especially suitable for the adaptive algorithm. The adaptive algorithm therefore makes use of large-dimensional VARs for nowcasting in the presence of mixed-frequency data feasible.

5 Conclusion

The increasing interest in large vector autoregressive models means that existing methods and algorithms are put to the test as large models are substantially more computationally demanding. In order to facilitate estimation of large-dimensional VARs with mixed-frequency data for real-time nowcasting we have developed methods that cope better with the large dimensionality. Our preferred algorithm is adaptive and augments the state vector in the simulation smoother at each time point as necessary, thereby avoiding to sample series we already know. The adaptive algorithm exploits the ragged edges of the data to ameliorate estimation when the data set is unbalanced. The algorithm consistently improves upon the Schorfheide and Song, (2015) simulation smoother, with the largest gain obtained when the number of variables and the number of lags are high. Given that recent influential papers such as Bańbura et al., (2010) and Carriero et al., (2019) estimate models with the number of variables being in the range 125–131 while the number of lags is set to 13, the situation in which we find the largest improvements is clearly of high interest in the literature. Our algorithm can therefore be used as an efficient building block in incorporating additional variables sampled at lower frequencies in existing large-dimensional VARs.

Appendix A Transitions in the SS15 simulation smoother

From compact to companion

At $t=T_{b}$, we obtain $a_{T_{b}+1}$ and $P_{T_{b}+1}$ for $\alpha_{t}=(x_{q,t}^{\prime},z_{q,t-1}^{\prime})^{\prime}$. If we instead use $\tilde{\alpha}_{t}=(z_{t}^{\prime},x_{t-p}^{\prime})^{\prime}$, we see that these now represent

[TABLE]

where $a_{t|t}=(a_{t|t}^{(1)^{\prime}},\dots,a_{t|t}^{(p)^{\prime}})^{\prime}$. Moreover,

[TABLE]

where $P_{T_{b}|T_{b}}^{(i,j)}$ is the $n_{q}\times n_{q}$ $(i,j)$ block of $P_{T_{b}|T_{b}}$.

Now, we can see that

[TABLE]

The companion-based filter is obtained by using $\tilde{a}_{T_{b}+1}$ and $\tilde{P}_{T_{b}+1}$ as starting values and then proceeding as usual. By smoothing, we obtain $\hat{\tilde{\alpha}}_{s}=E({\tilde{\alpha}}_{s}|y_{T},y_{T-1},\dots)$ for $s=T,\dots,T_{b}+1$.

From companion to compact

The final step is to move from the companion form to the compact representation. Doing so is accomplished by

[TABLE]

Hence, to move from the companion form back to the compact representation, we extract $\hat{\alpha}_{T_{b}+1}=E(\alpha_{T_{b}+1}|y_{T},y_{T-1},\dots)$ from $\hat{\tilde{\alpha}}_{T_{b}+1}$. After that, we can smooth using (12) initialized with $r_{T_{b}}$.

Appendix B Adaptive filtering: Example

System matrices for $t=1,\dots,T_{b}$

For $t=1,\dots,T_{b}$, monthly variables are fully observed and the quarterly variable is observed every third month. When $t$ corresponds to a month in which the quarterly variable is not observed, the fourth row of $Z_{t}$, $C_{t}$ and $G_{t}$ is deleted. We let $\Pi_{l}^{i,j}$ denote element $(i,j)$ of $\Pi_{l}$, and $w_{t}^{i,j}$ element $(i,j)$ of the lower-triangular Cholesky factor $W_{t}$.

[TABLE]

System matrices for $t=T_{b}+1$

[TABLE]

System matrices for $t=T_{b}+2$

[TABLE]