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.
88,793 characters · 21 sections · 5 citation commands
Nowcasting with Mixed Frequency Data Using Gaussian Processes
\doublespacing
\singlespacing{\footnotesizeContact: Massimiliano Marcellino ([email removed]), Department of Economics, Bocconi University. The views expressed in this paper are those of the authors and do not necessarily reflect the views of the Oesterreichische Nationalbank (OeNB) or the Eurosystem. We thank Eric Ghysels and participants at research seminars at the IHS, IIASA, OeNB, CIREQ-CMP Econometrics Conference and the IAAE 2024 for insightful comments and discussions. Hauzenberger and Pfarrhofer acknowledge funding by the Jubil\"aumsfonds of the OeNB, grants 18763, 18765. Codes and replication files are available at \href{https://github.com/mpfarrho/gp-midas}{github.com/mpfarrho/gp-midas}.}
\thispagestyle{empty} \doublespacing
This paper develops flexible nowcasting and short-horizon forecasting methods by combining elements from three strands of econometric literature. {First}, drawing from the mixed data sampling (MIDAS) framework proposed by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{ghysels2007midas} and \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{andreou2010regression}, we leverage techniques that permit the efficient use of predictors sampled at a higher frequency than the target variable \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][for a recent review]{ghysels2024econometric}. {Second}, the Big Data literature, which is based on the notion that exploiting a large set of predictors can improve predictive accuracy \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][in the context of mixed frequency models]{babii2022machine,babii2024nowcasting,mogliani2021bayesian,mogliani2024bayesian,beyhum2023factor,borup2023mixed}. {Third}, the machine learning literature, which suggests that flexible models can improve predictive accuracy by uncovering complex relationships among variables (see, e.g., \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{hasti2009elements} for a textbook, or \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{goulet2022machine} in a macroeconomic context).
We use Gaussian Processes (GPs) to estimate the unknown and potentially nonlinear relationships between a target variable and a large set of mixed frequency predictors. When modeling the conditional mean of the MIDAS framework with a GP, we refer to this model as GP-MIDAS. GPs have been used previously in single frequency economic applications, see e.g., \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2024forecasting} in the context of inflation forecasting, or \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{hauzenberger2024gaussian} who estimate effects of uncertainty shocks. By contrast, to handle mixed frequencies in this nonlinear context, we use several variants of restricted MIDAS polynomials in the spirit of \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{ghysels2007midas}, or the unrestricted MIDAS (UMIDAS) approach by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{foroni2015unrestricted}, where each single high-frequency predictor is split into many low-frequency ones.
The restricted GP-MIDAS can be interpreted as a structured (i.e., not randomly) compressed GP \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][for a discussion of the latter]{guhaniyogi2016compressed}, where the number of relevant predictors in our nonparametric model lies on a lower dimensional space. The MIDAS framework offers a natural angle to reduce the dimensionality of covariates. Therefore, it can improve the quality of feature extraction in a high-dimensional Big Data setting but also sharpen predictive inference \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{snelson2012variable, paananen2019variable, liu2020gaussian}. Such a modeling strategy is therefore expected to be effective in disentangling signal from noise in a MIDAS regression, particularly when dealing with either noisy measurements of high-frequency predictors and/or substantial correlations among predictors (which are inherent due to the design of the MIDAS setup).
To assess its merits, we benchmark GP-MIDAS against other machine learning techniques that have recently gained popularity in the macroeconomic (mixed frequency) forecasting literature. As a nonparametric alternative to GPs for modeling the conditional expectation of the dependent variable in MIDAS models, we consider Bayesian Additive Regression Tree \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[BART, originally proposed by][]{chipman2010bart} models. While GPs estimate nonlinear functions via a (theoretically infinite) mixture of Gaussian distributions, BART does so through a sum of regression trees. BART has been shown to perform quite well for nowcasting and forecasting with time series data in macroeconomic applications \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{clark2023tail,clark2024investigating,huber2023nowcasting}. A recent review of GPs and BART in a (multivariate) time series context is provided in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{marcellino2024chapter}. Parametric machine learning regressions equipped with global-local shrinkage priors reflect the recent (Bayesian) MIDAS literature \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{rodriguez2010mixed,mogliani2021bayesian,babii2022machine, kohns2023flexible}.\footnote{\bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{breitung2015forecasting} use small-scale homoskedastic nonparametric MIDAS models to predict inflation using daily data; \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{chan2024tvpmidas} discuss how to efficiently introduce time-varying parameters in Bayesian MIDAS regressions. Other recent examples of how to use machine learning approaches in a MIDAS context are \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{schnorrenberger2024harnessing} or \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{arro2024nowcasting}.} This competitor is conditionally linear in parameters and thus allows us to gauge the role of considering nonlinearities of an unknown form. Besides, we implement all competing model specifications with and without stochastic volatility (SV) in the error terms, as SV has been useful for improving short-horizon density forecasts in previous work based on linear mixed frequency models \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{carriero2015realtime,pettenuzzo2016midas,carriero2022nowcasting}.
We develop efficient Bayesian estimation algorithms based on Markov chain Monte Carlo (MCMC) sampling, which permits the use of very large sets of predictors to produce point, density, and tail forecasts in a computationally efficient manner. After a careful examination of the relative empirical performance of the various methods with simulated data, we apply these models to now- and forecast quarterly annualized real GDP growth and inflation in the GDP deflator with monthly predictors in a large-scale pseudo real-time evaluation scheme. The competing model combinations differ in the size of the information set (with up to about $120$ predictors), their conditional mean and variance parameterizations, and several MIDAS weighting schemes. We consider a long holdout period starting in the early 1990s.
There are several findings to report. First, nonlinear means are relatively more important when trying to improve forecasts than nowcasts, although they remain competitive for nowcasting as well. GP outperforms BART among the nonparametric methods, and SV is an important model feature irrespective of other specification details. Second, the size of the information set matters slightly more for linear competitors; the nonparametric versions typically perform well already with small data input. Relatedly, unrestricted MIDAS versions are rarely among the best-performing specifications, and it is beneficial in most cases to structure the predictors with specific lag polynomials. Third, the nonparametric models are strongest when the least information is available. For instance, on average, they perform relatively better for somewhat longer-horizon predictions.
The rest of the paper is structured as follows. Section (ref) presents the econometric framework. Section (ref) assesses the performance of the methods with simulated data. Section (ref) contains the empirical application on nowcasting and short-term forecasting of US GDP growth and inflation. Section (ref) summarizes the main findings and concludes. Additional details on the models, estimation and empirical results are provided in an Appendix.
Let $\{y_{t}\}_{t=1}^{T_L}$ denote a scalar target variable on the lower of two or more frequencies, observed $T_L$ times; and $\{z_{t}\}_{t=1}^{T_H}$ is a single high-frequency predictor observed $T_H$ times. Assuming an evenly spaced frequency mismatch between the target and the predictor yields an integer ratio $m=T_H/T_L$; $m$ thus indicates how often the high-frequency variables are observed in terms of one observation of the low-frequency variable.\footnote{For instance, when linking a quarterly dependent variable to monthly predictors there are three months per quarter and we have $m=3$. In this case, the fractional time indexes $t-2/3$ and $t-1/3$ refer to the first and second month of quarter $t$, respectively. It is worth noting that our approach can easily be extended to feature uneven frequency mismatches and multiple frequencies. Each frequency could then be assigned its own lag order and mismatch ratio, which we avoid here to keep the notation simple.}
Our more general framework is presented in the next section. For illustration, we begin with a simple MIDAS regression model $y_t = \mathfrak{B}(\mathfrak{L}^{1/m},\tilde{\bm{b}}) z_t + \epsilon_t$ where $\epsilon_t$ is a zero-mean error term, $\mathfrak{B}(\mathfrak{L}^{1/m},\tilde{\bm{b}})$ is a weighting function indexed by the high-frequency lag operator, $\mathfrak{L}^{p/m}z_t = z_{t-p/m}$, and $\tilde{\bm{b}} = (\tilde{b}_0,\hdots,\tilde{b}_{\mathbb{L}})'$ is an $(\mathbb{L}+1)\times1$-vector of parameters. In the spirit of distributed lag models, the weighting function is parameterized, as a finite-dimensional approximation, with $\mathfrak{B}(\mathfrak{L}^{1/m},\tilde{\bm{b}}) = \sum_{p=0}^{P_H - 1} \mathrm{B}(p,\tilde{\bm{b}}) \mathfrak{L}^{p/m}$, where $\mathrm{B}(p,\tilde{\bm{b}}) = \sum_{l = 0}^{\mathbb{L}} \tilde{b}_{l}\varphi_{l}(p)$. The $\varphi_{l}(p)$'s are basis functions which we stack in $(\mathbb{L} + 1)$-vectors $\bm{\mathrm{w}}_p = (\varphi_0(p),\hdots,\varphi_{\mathbb{L}}(p))'$, with associated weights in $\tilde{\bm{b}}$. We may collect these in a $P_H \times (\mathbb{L}+1)$-matrix $\bm{W} = (\bm{\mathrm{w}}_0,\bm{\mathrm{w}}_1,\hdots,\bm{\mathrm{w}}_{P_H - 1})'$, such that the linear MIDAS regression can be written as:
where $\tilde{\bm{z}}_t = (z_t, z_{t-1/m},\hdots,z_{t-(P_H-1)/m})'$ collects high-frequency lags of the predictor and $\bm{x}_t = \bm{W}'\tilde{\bm{z}}_t$. Define the $P_H$-vector $\bm{b} = \bm{W}\tilde{\bm{b}}$, then it is easy to see that this imposes specific shapes on the implied lag-specific parameters $\bm{b} = (b_{0},\hdots,b_{P_H-1})'$, that are weighted averages of the basis functions: $b_p = \sum_{l=0}^{\mathbb{L}} \tilde{b}_l \varphi_{l}(p)$ for $p = 0,\hdots,P_H-1$.
Loosely speaking, the matrix of weights maps the $P_H$-sized predictor vector to a (usually) lower dimension $\mathbb{L} + 1$, thereby reducing the number of parameters that need to be estimated at the cost of some lost flexibility. We break the linearity in the relationship between ${y}_t$ and $\bm{x}_t$ imposed in Eq. ((ref)) by assuming:
where $f(\bullet)$ is an unknown function that links a low-dimensional aggregate of the high-frequency lags nonlinearly to the target variable. On the one hand, we thus exploit the MIDAS setup with distinct basis functions to reduce the dimensionality of the problem in a structured way (i.e., subject to sensible parametric restrictions). On the other hand, we regain flexibility (which was lost due to finite-dimensional approximations and choices about basis functions) by estimating a possibly nonlinear function instead of a deterministic functional form. Rather than a linear combination of basis functions we obtain a nonlinear one. In our favored implementation, we rely on a GP prior to infer this function from the data.
The MIDAS structure (and underlying weighing approach) has several interesting and potentially favorable implications for common kernels used with GPs. It also offers computational advantages. Before we discuss these in detail, we lay out our general framework, which allows for many high-frequency predictors.
Let $\{y_{t}\}_{t=1}^{T_L}$ again be the scalar target variable and $\{\bm{x}_t\}_{t=1}^{T_L}$ now denotes an $M\times1$-vector comprised of predictors on the lowest of two or more frequencies, observed $T_L$ times. We consider regressions of the form:
The (potentially unknown and nonlinear) conditional mean function $f:\mathbb{R}^M\rightarrow\mathbb{R}$ will be specified below, and $\epsilon_{t}$ is a zero mean conditionally Gaussian error term with variance $\sigma^2_{t}$. As discussed in the single-predictor case above, the mixed frequency aspects of our work are due to the way the vector $\bm{x}_t$ is constructed. The goal is to model the low-frequency variable as a function of high-frequency variables. Let $\{\bm{z}_t\}_{t=1}^{T_H}$ with $\bm{z}_t = (z_{1t},\hdots,z_{Kt})'$ denote a $K \times 1$-vector of variables observed $T_H$ times on a higher frequency. Further, $\tilde{\bm{z}}_{kt} = (z_{kt},z_{kt-1/m},z_{kt-2/m},\hdots,z_{kt-(P_H-1)/m})'$ is a $P_H\times1$-vector of high-frequency lags of the $k$th predictor.
In our empirical work, the vector $\bm{x}_t$ contains $P_L$ lags of the dependent variable and $P_H$ lags of the $K$ predictors:
where $\bm{W}$ is the $P_H \times (\mathbb{L}+1)$ matrix of weights determined by the basis functions used to compress the high-frequency lags. Effectively, we thus have $M = P_L + K(\mathbb{L} + 1)$ predictors. We note that $\bm{x}_t$ could be augmented to include deterministic terms, lags of other low-frequency variables and latent or observed factors.
{\noindentRestricted and Unrestricted MIDAS}. When using UMIDAS (abbreviated u later) we have $\mathbb{L} + 1 = P_H$ and $\bm{W} = \bm{I}_{P_H}$ where $\bm{I}_N$ refers to an $N$-dimensional identity matrix. Hence, for each additional high-frequency indicator, UMIDAS results in $P_H$ additional covariates in Eq. ((ref)), leading to a rapid increase in the number of parameters, since $M = P_L + K P_H$. For an overview and a more detailed discussion, see \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{foroni2015unrestricted}. It suffices to note here that UMIDAS neither imposes restrictions nor explicitly models the distributed lag structure of the predictor-specific high-frequency lags.
Restricted MIDAS specifications alleviate overparameterization concerns via setting $\mathbb{L} \leq P_H$. Several parsimonious parameterizations of the weight function are available. We briefly summarize the main approaches that will serve as competing variants in our empirical work. In many cases, it is convenient to impose specific shapes on the weights and thus emphasize distinct high-frequency lags explicitly --- for instance, to reflect that more recent lags are usually more important than more distant ones. This is the first option we consider, following \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{ghysels2007midas}, referred to as the exponential Almon lag (xalm). Here, $\bm{W}$ collapses to a $P_H\times1$-vector and an implementation with two parameters sets its $(r+1)$th element to:
for $r = 0,\hdots,P_H-1$. While allowing for a nonlinear weighing function, we consider this implementation as an $\mathbb{L} = 1$ case due to $\bm{W}$ being a vector. Estimating $\theta_1,\theta_2$ allows for straightforward high-frequency lag selection in a data-driven manner. In particular, the shape of the weight function translates immediately into the implied lag polynomial on the coefficients (when assuming linearity). An important special case arises for $\theta_1 = \theta_2 = 0$ which yields equal weights. This is commonly referred to as a bridge model (br later on), and represents the simplest case we consider in our empirical work.
More generally, “indirect” parameterizations of the weights may be superior because they offer additional flexibility. We refer to these cases as $\mathbb{L} > 1$. Among these are non-orthogonalized power polynomials, usually referred to as \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{almon1965distributed} polynomials of degree $\mathbb{L}$, abbreviated as alm. These can be implemented by setting $\bm{\mathrm{w}}_p = (1,p,p^2,\hdots,p^{\mathbb{L}})'$, i.e., the associated basis functions are $\varphi_l(p) = p^l$. However, orthogonal polynomials are often preferable due to the reduced multicollinearity of compressed predictors. \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{babii2022machine}, for instance, suggest to use Legendre polynomials (leg) shifted to the interval $[0,1]$. \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{mogliani2024bayesian} propose to use Bernstein polynomials (ber), while \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{chan2024tvpmidas} discuss using Fourier basis functions (fou). We implement several of these options and provide additional details in Appendix (ref).
The choices about $\bm{W}$ and $\mathbb{L}$ (i.e., how to compress the predictors in an explicitly structured way) have several compelling implications when used in conjunction with GPs. We introduce GPs and discuss what these implications are next.
Rather than assuming a specific functional form we treat the conditional mean function $f(\bm{x}_{t})$ as unknown in our most general implementations and impose a GP prior directly on the functional relationship \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{williams2006gaussian}:
where $\mathcal{K}_{\bm{\kappa}}(\bullet)$ denotes a suitable covariance function (the kernel) which we define below. The properties of the conditional mean function, such as stationarity or smoothness, can be defined with distinct functional forms subject to a moderate number of tuning parameters. These are encoded in the vector $\bm{\kappa}$. Stacking the function values and data for a finite sample of $T_L$ observations, we obtain $\bm{f} = (f(\bm{x}_{1}),\hdots,f(\bm{x}_{T_L}))'$ and $\bm{X} = (\bm{x}_{1},\hdots,\bm{x}_{T_L})'$. The Gaussian process prior then takes the form of a $T_L$-dimensional multivariate Gaussian distribution:
where we assume zero means, and the $(t,\tilde{t})$th element of $\mathcal{K}_{\bm{\kappa}}(\bm{X},\bm{X})$ is given by $\mathcal{K}_{\bm{\kappa}}(\bm{x}_t,\bm{x}_{\tilde{t}})$. Put simply, the value of the covariance function governs the prior association of the functional values between time $t$ and $\tilde{t}$ conditional on the input vectors $\bm{x}_t$ and $\bm{x}_{\tilde{t}}$.
{\noindentKernel}. It is common to define kernels based on a distance metric and a few hyperparameters. We rely on the Euclidean distance and use a squared exponential kernel:
where $\bm{\Lambda}=\text{diag}(\lambda_1,\hdots,\lambda_M)$. A few general comments about the hyperparameters are in order. The unconditional variance of the prior, $\xi$, is sometimes also referred to as the signal variance (as opposed to the noise variance $\sigma_t^2$). This labeling can best be justified by noting that the marginal variance of the observations is Var$(y_t) = \xi + \sigma_t^2$. Notice also that $\lim_{(\bm{x}_t - \bm{x}_{\tilde{t}}) \rightarrow \infty}\mathcal{K}_{\bm{\kappa}}(\bm{x}_t,\bm{x}_{\tilde{t}}) = 0$; that is, the more the values of the predictors differ between $t$ and $\tilde{t}$, the less likely there is any noteworthy association in the function values a priori.
The so-called inverse length scales $\lambda_i$ in $\bm{\Lambda}$, by contrast, govern how quickly the conditional mean function varies with respect to changes in the $i$th input. As $\lambda_i\rightarrow0$, the impact of the $i$th predictor vanishes. Mechanically, this is due to the exclusion of $x_{it}$ from the distance metric, and $\lambda_i$ thus provides a rough gauge of variable importance. In practice, having predictor-specific inverse length scales is computationally costly and unreliable \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][for a recent discussion]{dance2022fast}. Usually one thus assumes a common inverse length scale, such that $\bm{\Lambda} = \lambda\bm{I}_M$, which reduces the number of hyperparameters in $\bm{\kappa} = (\xi,\lambda)'$, see also \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{guhaniyogi2016compressed} for a recent example.\footnote{We have experimented with alternative kernel parameterizations, such as using a Mat\'ern covariance function, but found the squared exponential kernel to perform very similarly, albeit with less hyperparameters to estimate.}
{\noindentInference}. Bayesian inference is based on manipulating the Gaussian distribution given by Eq. ((ref)). Using stacked notation such that $\bm{y} = (y_1,\hdots,y_{T_L})'$ and $\bm{\epsilon} = (\epsilon_1,\hdots,\epsilon_{T_L})'$, we obtain $\bm{y}\sim\mathcal{N}(\bm{0},\mathcal{K}_{\bm{\kappa}}(\bm{X},\bm{X}) + \bm{\Sigma})$ and $\bm{\Sigma} = \text{diag}(\sigma_{1}^2,\hdots,\sigma_{T_L}^2)$. In the machine learning jargon, $\bm{X} \in \mathbb{R}^{T_L\times M}$ is usually referred to as “training” inputs. For a general discussion, let $\tilde{\bm{X}} \in \mathbb{R}^{\tilde{T}_L\times M}$ denote the corresponding “test” values for which we would like to obtain inference about the unknown function values \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][for details]{williams2006gaussian}. This can be achieved by considering their joint distribution:
where $\tilde{\bm{f}}$ denotes the conditional mean function at the test inputs. Exploiting the properties of the multivariate Gaussian, the posterior is yet another multivariate Gaussian, with moments:
Computationally, this implies operations with matrices that are at most $\overline{T}_L = \max(T_L,\tilde{T}_L)$ dimensional. Specifically, the computationally costly operations are a matrix inversion, and, when it is required to sample from this distribution, a Cholesky decomposition. The computational complexities in general are thus cubic in $\overline{T}_L$ \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][for a discussion of the scalability of GPs]{banerjee2013efficient}. The number of predictors $M$, however, is by and large unimportant from a computational viewpoint.
This makes the MIDAS case particularly attractive for applying GPs because, on the one hand, we work with the lowest of all available frequencies, which means the least amount of observations. The side effect of increasing the number of predictors (particularly in the case of UMIDAS), on the other hand, and different from the case when assuming linearity, does not affect our framework other than implying an alternative covariance structure as captured with the kernel.
Another way of writing the GP is to consider its so-called weight-space representation. The prior in Eq. ((ref)) may be written as a conditionally linear regression:
It can be shown that this corresponds to a Bayesian linear regression with infinitely many basis functions \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see chapter 4.3 of][for details]{williams2006gaussian}. For an observation at time $t$ we have:
where $\bm{\psi}_t(\bm{X})'$ represents the $t$th row of $\bm{\Psi}(\bm{X})$ and $\psi_{t\mathrm{t}}(\bm{X})$ is the $\mathrm{t}$th element of this vector. This representation illustrates that our approach shares similarities with \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{mogliani2024bayesian}; in our case, the conditional distribution of $y_t$ is defined in terms of all observations --- due to our approach of modeling nonlinearities in the conditional mean --- rather than groups of predictors that arise from MIDAS-related schemes.
It is also worthwhile to assess a closely related alternative way of writing the posterior in Eq. ((ref)); using $\bm{\mathfrak{w}} = (\mathfrak{w}_1,\hdots,\mathfrak{w}_{T_L})' = \left[\mathcal{K}_{\bm{\kappa}}(\bm{X},\bm{X}) + \bm{\Sigma}\right]^{-1}{\bm{y}}$, we obtain $\mathbb{E}(\tilde{\bm{f}}|\bm{y}) = \sum_{t=1}^{T_L} \mathfrak{w}_t \mathcal{K}_{\bm{\kappa}}(\tilde{\bm{X}},\bm{x}_t)$. Let us revisit the single high-frequency predictor case of Section (ref) and consider the in-sample functional values where $\tilde{\bm{Z}} = (\tilde{\bm{z}}_1,\hdots,\tilde{\bm{z}}_{T_L})'$ such that $\bm{X} = \tilde{\bm{Z}}\bm{W}$:
These alternative expressions illustrate two main aspects. First, Eq. ((ref)) defines a theoretically infinite dimensional set of basis functions; conditioning on the actual sample yields a finite-dimensional representation. Second, the functional values are a linear combination of the kernel functions (which act as basis functions). These are in turn shaped by the presence of the MIDAS-related basis functions encoded in $\bm{W}$, i.e., how the high-frequency predictors are compressed.
This structure relates to the literature on compressed GP regression \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{snelson2012variable, guhaniyogi2016compressed, paananen2019variable, liu2020gaussian}. The key assumption for compressed nonparametric regressions is that the relevant predictors can be projected into a lower dimensional feature space. In a MIDAS framework, the compression of the feature and/or predictor space arises from taking into account the mixed frequency nature of the target variable $y_t$ and the raw predictors $\tilde{\bm{z}}_t$. Specifically, the weighting scheme implies variation in length scales across high-frequency predictors.
To investigate these claims more formally, it is instructive to consider the kernel written explicitly in terms of the high-frequency lags. To simplify notation, we ignore the lags of the dependent variable for the moment and let $\tilde{\bm{z}}_t = (\tilde{\bm{z}}_{1t}',\hdots,\tilde{\bm{z}}_{Kt}')'$, such that $\bm{x}_t = (\bm{I}_K \otimes \bm{W}')\tilde{\bm{z}}_t$. We may then write the implied kernel as:
Notice that this corresponds to Eq. ((ref)), which results from $\bm{\Lambda} = \lambda \tilde{\bm{\Lambda}}$, and thus the weights can be interpreted as emphasizing or de-emphasizing which lags enter the distance metric in the covariance function.
The implied inverse length scale matrix $\bm{\Lambda} = \bm{I}_K \otimes (\lambda\bm{W}\bm{W}')$ is block diagonal, and each of the blocks is associated with the high-frequency lags grouped by individual predictor. The common inverse length scale, $\lambda$, governs the “global” level of variability of the conditional mean function, and $\bm{W}$ provides distinct “local” adjustments through the MIDAS lags. This yields a family of “MIDAS kernels” among the class of squared exponential kernels, each with different implications for the covariance structure of the low-frequency functional values.
Figure (ref) provides a visual representation of the weights associated with each of the high-frequency lags, and the resulting implied matrix of inverse length scales per predictor, $\bm{\Lambda} = \lambda \bm{W}\bm{W}'$. We show only the $\mathbb{L}>1$ versions except UMIDAS which does not impose any structure and the resulting $\bm{\Lambda}$ is an identity matrix; the br implementation yields a matrix whose values are all equal to $P_H^{-2}$ and xalm specifies the shape of the weights explicitly. The ber and fou polynomials result in rather similar patterns, the former being more sparse; alm by comparison has more one-sided emphasis, whereas leg exhibits the most complex pattern across HF-lags.
We argued above that the GP prior is well-suited for the MIDAS framework, because of the way it handles variables and the rather general form of nonlinearities it can capture. There are, however, alternative ways of how to estimate the function $f(\bm{x}_{t})$. We briefly discuss the main features of such alternatives below and provide further details in Appendix (ref).
The first option we consider are regression trees. In our empirical work, we consider a Bayesian implementation as a competing specification and rely on BART:
where $\mathcal{T}_s$ denotes the tree structure, $\bm{\mu}_s$ is a vector of terminal node parameters associated with each tree, and we have $s=1,\hdots,S,$ regression trees. For these, we use the default prior setup and tuning parameters suggested in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{chipman2010bart}, which have been shown to work well also in a time series context \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{clark2023tail}.
Another, and arguably our simplest variation (see also Section (ref)), is the linear MIDAS regression:
where $\bm{\beta}$ is an $M \times 1$ vector of regression coefficients. This represents the basic framework used in recent related papers that are inspired by the machine learning literature. \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{babii2022machine} serves as an example of a classical econometric implementation, while Bayesian versions are developed in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{mogliani2021bayesian,mogliani2024bayesian}. Methodologically, these papers use variants of grouped shrinkage on the regression parameters, where groups are comprised of the compressed high-frequency lags for each of the $K$ predictors. Another noteworthy example is \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{kohns2023flexible}, who add a time-varying intercept term and $t$-distributed errors with stochastic volatility.
Our approach to inference is closely related to these implementations of (penalized or regularized) linear MIDAS regressions. We impose global-local shrinkage on the parameters $\bm{\beta}$ with a horseshoe prior \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[HS,][]{carvalho2010horseshoe}. Below we refer to this specification as Bayesian linear regression (BLR).
{\noindentError Variance}. For the case of homoskedastic errors (labeled hom), we assume $\sigma_{t}^2 = \sigma^2$ and $\sigma^2\sim\mathcal{G}^{-1}(a_0,b_0)$. To introduce stochastic volatility (SV), similar to \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{pettenuzzo2016midas} but without predictors in the state equation, let $\varsigma_{t} = \log(\sigma_{t}^2)$ and assume:
We rely on the prior setup of kastner2014ancillarity, and assume a Gaussian prior for the unconditional mean $\mu_\varsigma\sim\mathcal{N}(0,10)$, a transformed Beta prior for the autoregressive parameter $(\phi_\varsigma+1)/2\sim\mathcal{B}(5,1.5)$, and a Gamma prior on the state innovation variances $\sigma_\varsigma^2\sim\mathcal{G}\left(1/2,1/2\right)$. The prior on the initial state is $\varsigma_{h}\sim\mathcal{N}(\mu_\varsigma,\sigma_\varsigma^2/(1-\phi_\varsigma^2))$.\footnote{{We have experimented with heavy tailed error distributions such as $t$-distributed errors \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][in this context]{carriero2022addressing}. There were no noteworthy differences in predictive accuracy relative to our SV implementation with conditionally Gaussian errors.}}
{\noindentTuning Parameters}. We also need to choose priors that govern xalm. Here we use independent Gaussian priors: $\theta_1,\theta_2\sim\mathcal{N}(0,0.1^2)$. The final set of parameters that we need to equip with a prior are the hyperparameters that govern the GP kernel (the unconditional variance $\xi$, and the inverse length scale $\lambda$). We use Gamma priors on both:
Specifically, we set $a_\xi = 0.5$ and $b_\xi = 1$ to establish a rather vague prior on the unconditional variance of the kernel. It is only weakly informative because we normalize the target and predictor variables before estimation, to have unconditional means and variances equal to zero and one, respectively. {We apply this normalization to avoid unintended influences of any priors due to different scalings of the predictors. This aspect is important both in the context of the global-local prior in BLR, but also considering potential sensitivities of the kernel distance metric.} An inverse transformation is then used to recover predictions in the original scale of the data after estimation.
Prior information on $\lambda$ is imposed to reduce the risk of overfitting by pushing the inverse length scale to smaller values, in line with \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{hauzenberger2024gaussian} and \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2024forecasting}. We achieve this by setting $a_\lambda = 0.5$ and $b_\lambda = 0.1 s_y^2$ where $s_y^2$ is the residual variance of running an auxiliary linear AR$(1)$ model for the normalized dependent variable. Given our normalization, $s_y^2\in(0,1]$, this tightens the prior and shifts its expectation to smaller values when the target variable is more persistent (and vice versa).\footnote{In our empirical application, $s_y^2$ usually takes values between $0.7$ and $0.9$ for output growth, and smaller values between $0.3$ and $0.5$ for inflation, depending on the respective training sample. These numbers imply that our proposed prior scaling puts emphasis on smaller inverse length scales, but is not overly tight. Recall also that the MIDAS weights in $\bm{W}$ implicitly interact with $\lambda$ to determine $\bm{\Lambda}$ in Eq. ((ref)).} We have found this to be a useful semi-automatic scaling choice in this (time series) context.
The joint posterior distribution is not of a well-known form, and we therefore use several Metropolis-Hastings steps within a Gibbs sampling algorithm. Our algorithm cycles through the following steps:
In our empirical work, we iterate the sampling algorithm 12,000 times; then we discard the initial 3,000 as burn-in and use each third of the remaining draws for inference. This yields a set of 3,000 draws to be used for computing any object or moment of interest. {Our algorithm exhibits good convergence properties, as measured with standard MCMC diagnostics such as inefficiency factors.} Further details about the sampling algorithm are provided in Appendix (ref).
In this section, we assess the merits of our proposed model using synthetic data in a controlled environment. Our specified DGPs are inspired by recent literature \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{ghysels2007midas, andreou2010regression, ghysels2019estimating, mogliani2021bayesian}. However, we also wish to investigate the performance of the proposed MIDAS specifications (along with a large set of competitors) by assuming not only different variable importance profiles of high-frequency covariates (i.e., whether a DGP is relatively sparse or dense), but also different functional forms for the conditional mean. Consequently, to obtain a rich picture through the Monte Carlo exercise, our DGPs vary along these dimensions.
Across all DGPs, we assume $m = T_H/T_L = 3$, capturing the notion of a quarterly/monthly frequency mismatch we also encounter in our empirical application. Similar to mogliani2021bayesian and babii2022machine, we assume an independent AR(1) process for each high-frequency indicator, $z_{kt} = \rho_z z_{kt-1/3} + \varepsilon_{z,kt}$, where $\varepsilon_{z,kt} \sim \mathcal{N}(0, 1)$, and we specify a moderate degree of serial correlation with $\rho_z = 0.3$. This degree of persistence is reasonable for a typical monthly macroeconomic variable once it is transformed to stationarity (e.g., the monthly growth rate of industrial production or monthly inflation of an aggregate price series).
To specify $\bm W$, we closely follow the MIDAS literature and use three typical/characteristic weighting schemes --- fast-decaying, hump-shaped, and equal weights for $P_H = 12$ high-frequency lags of $z_{kt}$. For $k = 1, \dots, K$, these weighting schemes are defined by an exponential Almon lag (xalm), where we use $\theta_1 = 0$ and $\theta_2 = -1 \times 10^{-1}$ for fast-decaying weights, $\theta_1 = 5 \times 10^{-1}$ and $\theta_2 = -5 \times 10^{-2}$ for hump-shaped weights, and $\theta_1 = \theta_2 = 0$ for equal weights. We map the high-frequency covariates to the lower frequency as specified in Section (ref): $\tilde{\bm {x}}_t = (\bm I_K \otimes \bm W')\tilde{\bm {z}}_t$, with $\tilde{\bm{z}}_t = (\tilde{\bm{z}}_{1t}',\hdots,\tilde{\bm{z}}_{Kt}')'$ and $\tilde{\bm{z}}_{kt} = (z_{kt},z_{kt-1/3},z_{kt-2/3},\hdots,z_{kt-11/3})'$.
To relate our low-frequency endogenous variable $y_t$ to the low-frequency predictors $\tilde{\bm{x}}_t$, we assume the following autoregressive distributed lag (ADL) process:
with $\epsilon_{t}$ denoting homoskedastic errors with $\sigma^2 = 0.5$. The parameter $\rho_y$ introduces persistence in $y_t$. Similar to our high-frequency indicators, we assume $\rho_y = 0.3$. For $f(\tilde{\bm {x}}_t)$, we vary the functional form from highly nonlinear (NL) to linear (L):
Moreover, note that for both cases (NL and L), we assume that only the first five (compressed) predictors are important in defining movements in $y_t$, while all other covariates are considered to be irrelevant. In this sense, by varying $K$, we also vary the degree of sparsity. We use two different settings in terms of model complexity: one small but relatively dense regression with $K = 10$ --- resulting in a degree of sparsity of $50\%$ (five out of ten covariates are important) --- and one medium-sized but relatively sparse regression with $K = 25$ --- resulting in a degree of sparsity of $80\%$ (five out of $25$ covariates are important). In total, we have twelve DGPs with distinct features.
To evaluate the different model specifications, we simulate $R = 50$ replications for each DGP. We obtain (in-sample) $T_H = 750$ (and thus $T_L = 250$) observations. In addition, we simulate one out-of-sample observation at a lower frequency, used for model evaluation. For each DGP simulation, we predict this out-of-sample observation. As loss measures, we use the continuous ranked probability score \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[CRPS, see][and Section (ref)]{gneiting2007strictly}, which we average across the $R$ replications for each DGP.
Table (ref) shows CRPSs for the various MIDAS specifications (the results for mean absolute errors, which only measure point forecast accuracy, are very similar and shown in Appendix (ref)). MIDAS variants vary according to the specification of the conditional mean and weighting scheme. The left-hand panel displays results for our proposed GP-MIDAS variants, while the middle and right-hand panels show results with alternative conditional mean specifications (BART and linear). For the exercise using synthetic data, we consider a linear bridge regression (BLR-br) as a benchmark.
Overall, Table (ref) demonstrates that our proposed GP specifications perform well and are competitive across all DGPs. For five out of twelve DGP specifications, a GP variant is the overall winner. Particularly for highly nonlinear DGPs, the GP-MIDAS yields large margins compared to the BLR-br benchmark and the competitors. Interestingly, BART shows a somewhat more consistent performance across weighting schemes and the $\mathbb{L} > 1$ cases are useful when using a tree-based approach.
Taking a closer look, we observe that GP-MIDAS specifications, particularly those equipped with an exponential Almon lag (xalm), show the strongest performance for a nonlinear DGP with fast-decaying or hump-shaped weights. It is also noteworthy that since DGPs are simulated with $\mathbb{L} = 1$, for GP-MIDAS it pays off to compress the raw predictors to a parsimonious specification, while the unrestricted GP-MIDAS or the richer specifications with $\mathbb{L} > 1$ are in a sense overparameterized. By contrast, for BLR, the added flexibility when $\mathbb{L} > 1$ can be marginally beneficial in terms of predictive accuracy. This finding captures the notion that a misspecified linear model can partially approximate nonlinearities through a larger information set.
A GP-MIDAS model with the xalm weighting scheme is capable of recovering the truly underlying weight structure, thus improving upon linear benchmarks. For example, the NL-fast DGPs can be considered as the DGP specification featuring the most nonlinearities, as nonlinearities do not only arise from the conditional mean but also arise from a highly nonlinear weighting scheme implied by the fast-decaying weights. In such a case, a linear bridge regression is extremely misspecified, while a GP-xalm detects those nonlinearities, resulting in large predictive gains. Conversely, when DGPs are linear (which is when the nonparametric versions will be less efficient than BLR by construction), both GP and BART still show very competitive predictive metrics. This implies that our setup achieves sufficient regularization, and alleviates concerns about overfitting.
The quarterly target series are sourced from FRED-QD and the monthly predictors are from the FRED-MD database \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see][]{mccracken2016fred,mccracken2020fred}. Our target variables are real GDP growth (GDPC1) and inflation in the GDP deflator (GDPCTPI). Both are used in the form of annualized growth rates. Our full sampling period runs from 1963Q1 to 2023Q2, and in our out-of-sample evaluation scheme, the initial estimation (training) sample ends in 1989Q4. For nowcasts and forecasts, the first targeted quarter in our holdout sample is 1990Q1. Charts of the target variables and subsamples we refer to below are provided in Appendix (ref).\footnote{Full denotes the whole holdout; sample splits are Pre Covid and Post Covid, where the full sample is split at 2019Q4. We further differentiate Recession versus \texttt{Expansion} by considering all dates qualified as falling into a recessionary period by NBER as the former, and the rest of the sample as the latter.}
In the out-of-sample evaluation scheme, we rely on final vintage data, but truncate each of the training samples according to the release calendar as we simulate going forward in time. That is, in our “pseudo real-time” evaluation setup, we abstract from data revisions but respect the publication schedule of the variables. We focus on the current quarter nowcast indexed $h \in \{0,1/3,2/3\}$ and the one-quarter-ahead forecast $h \in \{1,4/3,5/3\}$ --- both are updated at the end of each month during a quarter, and $h$ encodes the distance to the target quarter in months. To give an example, a nowcast for 2023Q1 made in 2023M2 (which is referred to as $h = 1/3$), due to publication delay, uses monthly information up to 2023M1 \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see also][]{babii2022machine}. A typical monthly predictor would thus take the form $\tilde{\bm{z}}_{kt} = (z_{kt-2/3},z_{kt-1},\hdots,z_{kt-(P_H+1)/3})'$. More generally, for any $h$, we have $\tilde{\bm{z}}_{kt} = (z_{kt-h-1/3},z_{kt-h-2/3},\hdots,z_{kt-h-(P_H)/3})'$.
Potential predictors are defined based on information sets of different sizes: small (s, $K=12$), medium (m, $K=23$) and Big Data (b, $K=116$). These sets are constructed so that the respective larger one nests the smaller ones. A detailed description of which predictors are used in each of them is provided in Appendix (ref). All these predictor variables are processed using the suggested transformation codes to achieve approximate stationarity. We use $P_L = 4$ lags of the target variable as is common with quarterly data; and, reflecting the month/quarter frequency mismatch, we set $P_H = 12$ based on $P_H = m P_L$ for consistency. In line with \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{mogliani2021bayesian,mogliani2024bayesian} and \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{babii2022machine} we choose $\mathbb{L} = 3$ if applicable, and also consider $\mathbb{L} = 5$ for Bernstein and Legendre polynomials.
Besides differently sized information sets, our set of competing models is defined as follows. We vary the respective approach to estimating the conditional means $f(\bullet)$ by using linear (BLR), GP and BART versions. In addition, we consider homoskedastic (hom) and heteroskedastic (sv) implementations. The final variation is due to how we construct the design matrices for distinct MIDAS types. Here, we use both UMIDAS and the restricted approaches we discussed in Section (ref). As a simplistic benchmark, we also estimate univariate $AR(P_L)$ models which are updated each quarter. An overview is provided in Table (ref). When we explicitly refer to a model in the text, we structure their IDs as mean-variance-midas-size.
{\noindentLoss Functions}. We evaluate our predictions using several distinct loss functions that measure different aspects of predictive accuracy. Specifically, we use the quantile score (QS) and the (quantile-weighted) continuous ranked probability score (CRPS), see \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{giacomini2005evaluation,gneiting2007strictly}.
The QS for quantile $\tau\in(0,1)$ of the forecast of target variable $y_t$ is defined as $\text{QS}_{\tau,t} = 2\left(y_{t}-\hat{y}_{\tau,t}\right)\left(\tau - \mathbb{I}\left\{y_{t}\leq \hat{y}_{\tau,t}\right\}\right)$, where $\hat{y}_{\tau,t}$ indicates the $\tau$th quantile of the predictive distribution. The indicator function $\mathbb{I}\left\{y_{t}\leq \hat{y}_{\tau,t}\right\}$ has a value of $1$ if the realized value is at or below the quantile forecast, and $0$ otherwise. Note that $\text{QS}_{0.5,t}$ is the absolute error (which we use in the form of the mean absolute error, MAE, as our point forecast loss). Using the QS, we can then define the (quantile-weighted) CRPS following \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{gneiting2011comparing}:
where $\mathfrak{w}_{\text{V},\tau}$ indicate weights that emphasize different parts of the distribution. The default CRPS results when we use equal weights, but we also consider different weighting schemes that target the left tail $\mathfrak{w}_{\text{\texttt{L}},\tau} = (1-\tau)^2$ and the right tail $\mathfrak{w}_{\text{\texttt{R}},\tau} = \tau^2$. We approximate the integral above with a sum based on a grid of quantiles, $\tau\in\{0.05,0.06,\hdots,0.94,0.95\}$.
The set of competing models is vast. We therefore slice our empirical results along several dimensions to begin with a concise overview before drilling deeper into details.
First, we take a bird's-eye view of which broad model features --- on average --- produce better forecasts. We choose the CRPS averaged for various subsamples over the holdout as our metric of choice for measuring overall forecast performance by target variable. The logarithms of these losses across specifications are then regressed on dummies comprised of the categories determined by the underlying means, variances, sizes of the information set, and MIDAS type (i.e., our broad model features).
Our “baseline categories” reflect the linear homoskedastic bridge model and the small information set --- which are the categories left out by the groups of dummies. The parameters are estimated with OLS and multiplied by $100$. These are the numbers shown in Figure (ref), which can be interpreted as percentage gains/losses due to the extensions indicated in the rows of each panel relative to the losses of the benchmark. {Panel (a) shows the results including all dummies, panel (b) subsets to only GP implementations and thus omits the dummy for the mean component.} The columns refer to the nowcast/forecast horizon $h$.\footnote{When using the mean absolute error (MAE) as explained loss instead, we obtain qualitatively very similar patterns.}
Starting with panel (a), there is heterogeneity concerning the predicted variable. Addressing choices about the conditional mean first, we find that the nonparametric models pay off particularly as $h$ becomes larger. By contrast, for the nowcasts, we find that GP and BART are often performing worse than when assuming linearity a priori (sometimes even significantly so). It is worth noting, however, that this overview abstracts from interactions of model features. So even when on average losses are larger when moving from BLR to GP or BART, keeping everything else fixed, this is not necessarily the case for distinct model specifications. Indeed, we zoom into model-specific performance below and find that several versions of the GP implementation are the overall best-performing models.
Adjusting the MIDAS weights provides the greatest relative gains, at least for GDPC1. The two-parameter specification xalm seems to succeed in striking a balance between flexibility and simplicity. For GDPCTPI, this aspect is somewhat less important. The corresponding accuracy premium appears consistently across horizons and target variables, although the magnitudes differ somewhat. Another specification detail that consistently improves scores is sv. Usually, these improvements are larger as the forecast horizon increases. The size of the information set does not matter much, although it must be said that using b sometimes significantly hurts predictive accuracy by a few percentage points. These general patterns also apply when subsetting to GP implementations in panel (b), with differentials due to varying these model features being somewhat more muted.
Our second bird's-eye view set of results aims at identifying those specifications that perform best consistently (meaning across all now-/forecast horizons and all subsample splits). To answer the question posed in the title of this subsection, we compute the model confidence set (MCS) of \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{hansen2011model} for all possible combinations of losses, horizons, and subsamples.
{The output of the MCS procedure is a subset of models which includes the best-performing one with a specified level of confidence. The specifications included in the MCS are statistically indistinguishable at this pre-defined level. Our results are based on the T$_{R}$ statistic defined in \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{hansen2011model} and we use a $10$ percent level of significance. This procedure provides us with a framework to distinguish models that consistently perform well, as opposed to ones that are strictly dominated by others. Having computed the MCS across all subsets of results, we count how often each model was among the superior set of models, and present inclusion percentages alongside the minimal average loss (across horizons and subsamples) for the best three specifications in Table (ref).}
Again there are differences between target variables, but also, with respect to predictive loss metrics. Consistent with our exploratory regressions, xalm features very prominently, irrespective of how the conditional mean is specified or how large the information set is. Interestingly, short-horizon predictions of inflation as measured by the GDP deflator tend to require larger models (the Big Data information set appears more often) than for real GDP growth, where small and medium models are more common.
One clear lesson of this exercise is that BART is dominated (both in terms of averaged losses and selection frequency) by GP and BLR, and the latter two often exhibit similar performance. For GDPC1, we find GP-sv-xalm-s to perform well very consistently. This model has the lowest average losses for both point forecasts (as measured by the MAE) and the density forecasts captured by CRPSs and is included in the MCS in all cases for the former metric, and 93% of all cases for the latter. By contrast, for GDPCTPI, while the average loss metrics for the best-performing model specifications are very close in many cases, the MCS favors BLR. As we will see below, this is due to heterogeneities in the performance of GP both across horizons and across subsamples.
Before proceeding with a more granular discussion of our nowcast and forecast results, it is also worth briefly discussing estimation times. Theoretical considerations about the computational complexity were alluded to in Section (ref); below we provide a summary based on estimation times for our full out-of-sample forecast simulation.\footnote{{This simulation was run in parallel on a high-performance computing cluster with Intel Xeon Platinum 8358 32C 2.6GHz CPUs, a total number of 1,792 cores and 16,896GB of RAM.}} It is worth reiterating that the bottleneck for BLR in its baseline implementation is the number of predictors $M$, while for GP the relevant dimension is the number of observations $T_L$. When we employ a fast-sampling algorithm for BLR (we use this approach when $T_L < M$ because the computational burden for some of our Big Data models is otherwise insurmountable, especially in a recursive prediction exercise), the limiting factor for BLR is also $T_L$. Due to how BART models the conditional mean, neither of the two dimensions is particularly important.
{Figure (ref) shows the average duration of producing 1,000 draws (in seconds). Since we use an expanding window, the number of observations differs for each training and holdout sample, which is why we indicate the median with colored bars. The black error bars (lines) mark the $25$th and $75$th percentiles of estimation times across all training samples. These results are based on models featuring SV (they are roughly the same for the homoskedastic versions).}
A few aspects are worth noting. Disregarding xalm, it always takes virtually the same amount of time to estimate the GP specification regarding the size of the information set and polynomial degree. This is due to the number of observations being the constraint, which is reflected in rather wide bars across differently-sized training samples. The same is true when $M > T$ for BLR due to the properties of the fast-sampling algorithm we use in this case. It is worth mentioning that when using a default algorithm for BLR, the scales of this chart would be obscured because BLR turns out to be computationally infeasible as the number of variables becomes huge. At least for the Big Data information set, it roughly takes the same amount of time to estimate GP and BLR, and GP is computationally more efficient when using UMIDAS (i.e., when the number of predictors is huge). For BART, the number of predictors solely affects the number of potential splitting variables, which is why it is on average the fastest of our competitors across model variations.
{It is also worth noting why the xalm implementation is somewhat slower than the others. This is because it requires several comparatively costly matrix operations in the context of the MH updates. If one were to calibrate the parameters $\theta_1,\theta_2$ beforehand, the computational burden would be roughly the same as in the other cases.}
To drill deeper into model-specific predictive accuracy, we select a subset of models that perform best on average as we did for the MCS. We now benchmark these results relative to an AR$(P_L)$-model which is a useful baseline for short-horizon forecasts of output growth and inflation. It reflects the case where no within-quarter information is used. Relative losses are shown as ratios in Table (ref), and the models are selected based on the sum of the CRPS across all horizons being minimal. We then display the best three GP models and the respective best BLR and BART one, for each target variable. Below we mostly focus on discussing the CRPS, since qualitative differences between point and density forecasts are often muted. Additional metrics for other specifications are collected in Appendix (ref).
An obvious starting point when discussing these results is to note that virtually all cells of Table (ref) are colored blue, which is used to indicate that the respective model outperforms the benchmark. In other words, adding within-quarter information for short-horizon predictions --- either when forecasting the next quarter, or when nowcasting the current quarter --- pays off irrespective of whether one addresses any nonlinearities or not. Turning to a discussion of which models end up specifically in this table, we again find MIDAS-type xalm to perform well across all conditional mean specifications, and only in the case of BART, we find a different implementation being included in this set of models. Interestingly, and different from BLR and GP, BART is capable of better extracting information without introducing prior weighting via lag polynomials, though it must be said that this specification is always beaten by another conditional mean implementation. Second, controlling for heteroskedasticity (which we do with sv) pays off. The only exception, again, is BART. This corroborates the findings of \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2023tail}, who provide a discussion and some empirical evidence of how the regression trees can capture specific forms of heteroskedasticity even when the error terms are homoskedastic.
Focusing on the size of the information set, it is worth mentioning that for the GP it appears that s is sufficient to attain superior relative predictive losses, and the b-sized version only occasionally improves upon the small dataset. By contrast, for BLR, assuming linearity a priori appears to be slightly restrictive, but more information in the form of either the medium or Big Data input appears to help. Another interpretation of this pattern is that the nonparametric features offset omitted variables partly, and the bigger the information set becomes, the less there is a need for modeling nonlinearities explicitly.
A related important observation is that BLR is ranked first (indicated by the bold numbers) in many cases for the nowcasts when considering the full evaluation period. For forecasts, the GP improves upon the linear version for all metrics and both target variables. We conjecture that this is due to the notion that as more information becomes available, the nonlinearities become less emphasized and there is less need for the flexibility in conditional means that our proposed specifications provide. Indeed, this finding is in line with the literature; e.g., \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2023tail} find that nonparametric models offer larger improvements at longer horizons (where nonlinearities are typically stronger).
When comparing the full holdout sample to the pre-Covid period, we find that the nonparametric models perform better, relatively speaking, when excluding the pandemic observations from the evaluation period. This finding at first seems puzzling, and, in our shorter-horizon predictive context, at odds with the longer-horizon, forecasting-oriented literature. Specifically, earlier related papers have often found more flexible models to offer gains in predictive accuracy in the presence of such outlying observations, due to better robustness properties. Indeed, it is a lack of responsiveness of some of our nonparametric models to huge-variance shocks (reflected timely in the monthly series) that can explain this puzzle. We investigate this claim in the next section, after a more thorough analysis of the resulting predictive distributions.
For the following analysis, we pick the two best-performing models, GP-sv-xalm-s (in blue) and BLR-sv-xalm-m (in green). We refrain from showing BART because it is dominated by the others for the most part. Figure (ref) shows the predictive distributions (68 percent credible set) over time for the end-of-quarter forecast ($h = 1$) and nowcast ($h = 0$) alongside realized values (black crosses).
As one would expect, the credible sets narrow as additional information becomes available to produce the predictions, for both specifications. Deviations between the predicted quantiles, comparing the GP and BLR, are strongest in and around recessionary episodes of the US economy. But we do not detect any noteworthy asymmetries or multimodalities and the like --- features that often arise when related econometric methods are used for forecasting. To some extent, this is by design: we use a direct predictive equation and are interested in very short-horizon predictions. So any nonlinearities that might arise from iterative propagation mechanisms and feedback effects \bibpunct{(}{)}{;}{a}{,}{;}\natbibcitep[see, e.g.,][]{huber2023nowcasting}, or generally when focusing on higher-order forecasts, do not appear in our out-of-sample exercise.
Due to this apparent symmetry of the predictive distributions, the tail forecast metrics we consider (CRPS-L for downside risk and CRPS-R for upside risk) are mostly in line with the unweighted CRPS density metric regarding model performance rankings. This is why we report them in Appendix (ref). One noteworthy pattern is that upside risks in both target variables seem to be captured relatively better by the GP-versions, in line with the findings of \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2024investigating}. In terms of the magnitudes of improvements over the autoregressive benchmark, the monthly predictors help even more when the interest is on predicting tail risk. Which of the predictors are the most important is what we aim to investigate next.
{\noindentMeasuring Variable Importance}. One shortcoming of machine learning approaches, while usually offering good forecast performance, is that in most cases they lack interpretability. While in a linear regression, the importance of any variable is naturally measured with the size of its regression coefficient, this is not as easily assessed in nonlinear frameworks. To identify which predictors drive the respective predictive distributions we thus resort to auxiliary tools. Specifically, we rely on linearly approximating the relationship between the predictive distribution and the corresponding predictors, inspired by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{woody2021model} and applied recently in a similar context to ours by \bibpunct{(}{)}{,}{a}{,}{,}\natbibcitet{clark2024forecasting}.
Let $\hat{y}_{t+h}$ denote the median of the predictive distribution for horizon $h$;\footnote{This procedure may also be used for different quantiles along the predictive distribution. Due to the predictive distributions being mostly symmetric in our empirical exercise, we use the median.} we define the following Lasso problems for each model specification:
where $T_0$ indexes the first quarter (i.e., 1990Q1) of our holdout sample and $T$ marks the last quarter (2023Q2), $\bm{\mathfrak{b}}_{h} = (\mathfrak{b}_{h1},\hdots,\mathfrak{b}_{hM})'$ are regression coefficients and $\varrho \geq 0$ is a tuning parameter (determined through cross-validation) that governs the weight of the penalty term and degree of sparsity in $\bm{\mathfrak{b}}_{h}$. This procedure recovers a sparse vector $\hat{\bm{\mathfrak{b}}}_{h}$ that measures which predictors in $\bm{x}_t$ shift the location of the predictive distributions (specific to each model variant) throughout the holdout sample.
{We summarize our results for the respective end-of-quarter forecast and nowcast in Table (ref). Additional results for other horizons are provided in Appendix (ref). For this exercise, we pick the Big Data variants of BLR and GP, with stochastic volatility and an xalm MIDAS piece. We pick the largest information set for illustration, because these implementations (while not being the overall outright best-performing models at least in the GP case) are competitive, and none of the predictors are excluded a priori. When there are multiple lags of the same variable, we take the sum of the respective coefficients.}
A general insight from this approach to measuring variable importance is that models for inflation are denser, in the sense that more predictors shape the forecast distributions, than those for real output growth. This is true for both the linear and nonlinear models. Interestingly, BLR and GP by and large tend to select the same variables, although there are some modest differences in magnitudes of importance. A noteworthy deviation from this pattern occurs for forecasts of output, where GP-MIDAS only has two non-zero predictors while many more are present for BLR. And, while persistence via the autoregressive low-frequency lags is important for inflation, this is not necessarily the case for output.
{Turning to individual predictors, it is noteworthy that components of industrial production (most notably, durable materials and nondurable consumer goods) shift predictions of output growth, as do personal consumption expenditures (DPCERA3M086SBEA) when using BLR. By contrast, GP puts emphasis on manufacturing and trade industry sales (CMRMTSPLx) and industrial production of the manufacturing sectors. For inflation, selected variables are more homogeneous across the two competing specifications, and the GP model relies on a few more predictors. Besides the autoregressive lags of the target variable, labor market variables (e.g., average weekly hours in manufacturing, AWHMAN), business inventories, and variables related to housing (HOUST and PERMIT) seem to matter, as does the real M2 money stock.}
{\noindentA Case Study of Two Recessions}. In terms of our earlier note about the role and importance of outlying observations, it is particularly illustrative to focus on two distinct economic episodes: the Great Recession and the Covid pandemic. We do so for real GDP growth in Figure (ref), and note that similar arguments can be made in the context of inflation as measured by the GDP deflator. This chart is to be understood as follows. Each of the colored lines and shades refer to the predictive median and 68 percent credible set for a single target observation. Since we update both our forecasts ($h = 5/3, 4/3, 1$) and nowcasts ($h = 2/3, 1/3, 0$) every month, we have six predictions for each quarterly reading of the data (which is marked with a filled black square).
Starting with a review of the Great Recession through the lens of both of these two models, it is noteworthy that during the initial phase of the recession that started in 2008Q1, the prediction paths across horizons are almost identical for both specifications. A key difference then occurs when the trough is reached in 2008Q4: BLR nowcasts this low to persist and overpredicts the severity of the observation in 2009Q1 at $h=0$. By contrast, the GP version is pretty much spot on at the end-of-month horizon for this quarter. These dynamics are due to the nowcast of the nonparametric model responding less strongly to new information in the monthly predictors.
Indeed, an even more prominent example of this phenomenon is visible in the right panel of Figure (ref), which shows the quarters most strongly affected by the Covid pandemic. The prediction paths targeting 2020Q2 are especially important for our argument. The difference between the full evaluation and the pre-Covid subsample is essentially due to this period (which is also why the differences are often statistically insignificant): the nonparametric predictions are only very modestly affected by the updated monthly information, while the BLR predictions shift strongly because it linearly extrapolates past patterns.
For instance, the huge downward spike at $h = 1/3$ (which is the nowcast made in May) reflects the data release for the April reading of the unemployment rate at a historic all-time high and industrial production growth at an all-time low. This lazy reaction of the nonparametric predictions is indeed what provides robustness in a longer-horizon forecasting context, where the goal is to predict future dynamics in periods strictly after a shock materializes. By contrast, when nowcasting, we partly observe the shocks through the monthly variables as they happen. And these more sluggish updates of the GP prediction, as additional information becomes available, can also hurt predictive accuracy in this context.
We develop Bayesian nonparametric methods to be used in combination with MIDAS regressions and discuss their implications. Specifically, we consider GP and BART as flexible alternatives to the penalized linear framework. Our models are also equipped with stochastic volatility, and they are further differentiated with respect to specifics about MIDAS weighting schemes. We show how different MIDAS weights used in GP-MIDAS yield a framework related to compressed GP regression.
After showing their good performance in a set of simulation experiments, the various nonparametric MIDAS specifications are applied to nowcasting and short-term forecasting US real output growth and inflation in the GDP deflator in a large-scale out-of-sample exercise. The set of potential predictors features up to almost $120$ variables for our evaluation period, which starts in the early 1990s. The new models are computationally efficient, and they are competitive while offering gains in predictive accuracy for point, density, and tail forecasts.
{\setstretch{0.9} \addcontentsline{toc}{section}{References} }\doublespacing