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.
141,373 characters · 24 sections · 39 citation commands
Forecasting with Dynamic Panel Data Models
\thispagestyle{empty}
JEL CLASSIFICATION: C11, C14, C23, C53, G21
KEY\ WORDS: Bank Stress Tests, Empirical Bayes, Forecasting, Panel Data, Ratio Optimality, Tweedies Formula
\thispagestyle{empty} \setcounter{page}{0}
The main goal of this paper is to forecast a collection of short time series. Examples are the performance of start-up companies, developmental skills of small children, and revenues and leverage of banks after significant regulatory changes. In these applications the key difficulty lies in the efficient implementation of the forecast. Due to the short time span, each time series taken by itself provides insufficient sample information to precisely estimate unit-specific parameters. We will use the cross-sectional information in the sample to make inference about the distribution of heterogeneous parameters. This distribution can then serve as a prior for the unit-specific coefficients to sharpen posterior inference based on the short time series.
More specifically, we consider a linear dynamic panel model in which the unobserved individual heterogeneity, which we denote by the vector $\lambda_i$, interacts with some observed predictors:
Here, $(W_{it-1},X_{it-1},Z_{it-1})$ are predictors and $U_{it}$ is an unpredictable shock. Throughout this paper we adopt a correlated random effects approach in which the $\lambda_i$s are treated as random variables that are possibly correlated with some of the predictors. An important special case is the linear dynamic panel data model in which $W_{it-1}=1$, $\lambda_i$ is a heterogeneous intercept, and the sole predictor is the lagged dependent variable: $X_{it-1}=Y_{it-1}$.
We develop methods to generate point forecasts of $Y_{iT+1}$, assuming that the time dimension $T$ is short relative to the number of predictors $(W_{iT},X_{iT},Z_{iT})$. The forecasts are evaluated under a quadratic loss function. In this setting an accurate forecasts not only requires a precise estimate of the common parameters $(\alpha,\rho)$, but also of the parameters $\lambda_i$ that are specific to the cross-sectional units $i$. The existing literature on dynamic panel data models almost exclusively studied the estimation of the common parameters, treating the unit-specific parameters as a nuisance. Our paper builds on the insights of the dynamic panel literature and focuses on the estimation of $\lambda_i$, which is essential for the prediction of $Y_{it}$.
The benchmark for our prediction methods is the so-called oracle forecast. The oracle is assumed to know the common coefficients $(\alpha,\rho)$ as well as the distribution of the heterogeneous coefficients $\lambda_i$, denoted by $\pi(\lambda_i|\cdot)$. Note that this distribution could be conditional on some observable characteristics of unit $i$. Because we are interested in forecasts for the entire cross section of $N$ units, a natural notion of risk is that of compound risk, which is a (possibly weighted) cross-sectional average of expected losses. In a correlated random-effects setting, this averaging is done under the distribution $\pi(\lambda_i|\cdot)$, which means that the compound risk associated with the forecasts of the $N$ units is the same as the integrated risk for the forecast of a particular unit $i$. It is well known, that the integrated risk is minimized by the Bayes predictor that minimizes the posterior expected loss conditional on time $T$ information for unit $i$. Thus, the oracle replaces $\lambda_i$ by its posterior mean.
The implementation of the oracle forecast is infeasible because in practice neither the common coefficients $(\rho,\alpha)$ nor the distribution of the unit-specific coefficients $\pi(\lambda_i|\cdot)$ is known. To obtain a feasible predictor, we extend the classical posterior mean formula attributed to separate works of Arthur Eddington and Maurice Tweedie to our dynamic panel data setup. According to this formula, the posterior mean of $\lambda_i$ can be expressed as a function of the cross-sectional density of certain sufficient statistics. Conditional on the common parameters, this distribution can then be estimated either parametrically or non-parametrically from the panel data set. The unknown common parameters can be replaced by a generalized method of moments (GMM) estimator, a likelihood-based correlated random effects estimator, or a Bayes estimator.
Our paper makes three contributions. First, we show in the context of the linear dynamic panel data model that a feasible predictor based on a consistent estimator of $(\rho,\alpha)$ and a non-parametric estimator of the cross-sectional density of the relevant sufficient statistics can achieve the same compound risk as the oracle predictor asymptotically. Our main theorem extends a result from BrownGreenshtein2009 for a vector of means to a panel data model with estimated common coefficients. Importantly, this result also covers the case in which the distribution $\pi(\lambda_i|\cdot)$ degenerates to a point mass. As in BrownGreenshtein2009, we are able to show that the rate of convergence to the oracle risk accelerates in the case of homogeneous $\lambda$ coefficients. Second, we provide a detailed Monte Carlo study that compares the performance of various implementations, both non-parametric and parametric, of our predictor. Third, we use our techniques to forecast pre-provision net-revenues of a panel of banks.
If the time series dimension is small, our feasible predictor performs much better than a naive predictor of $Y_{iT+1}$ that is based on within-group estimates of $\lambda_i$. A small $T$ leads to a noisy estimate of $\lambda_i$. Moreover, from a compound risk perspective, there will be a selection bias. Consider the special case of $\alpha=\rho=0$ and $W_{it}=1$. Here, $\lambda_i$ is simply a heterogeneous intercept. Very large (small) realizations of $Y_{it}$ will be attributed to large (small) values of $\lambda_i$, which means that the within-group mean will be upward (downward) biased for those units. The use of a prior distribution estimated from the cross-sectional information essentially corrects this bias, which facilitates the reduction of the prediction risk if it is averaged over the entire cross section. Alternatively, one could ignore the cross-sectional heterogeneity and estimate a (misspecified) model with a homogeneous coefficient $\lambda$. If the heterogeneity is small, this procedure is likely to perform well in a mean-squared-error sense. However, as the heterogeneity increases, the performance of a predictor that is based on a pooled estimation quickly deteriorates. We illustrate the performance of various implementations of the feasible predictor in a Monte Carlo study and provide comparisons with other predictors, including one that is based on quasi maximum likelihood estimation of the unit-specific coefficients and one that is constructed from a pooled OLS estimator that ignores parameter heterogeneity.
In an empirical application we forecast pre-provision net revenues of bank holding companies. The stress tests that have become mandatory under the Dodd-Frank Act require banks to establish how revenues vary in stressed macroeconomic and financial scenarios. We capture the effect of macroeconomic conditions on bank performance by including the unemployment rate, an interest rate, and an interest rate spread in the vector $W_{it-1}$ in ((ref)). Our analysis consists of two steps. We first document the one-year-ahead forecast accuracy of the posterior mean predictor developed in this paper under the actual economic conditions, meaning that we set the aggregate covariates to their observed values. In a second step, we replace the observed values of the macroeconomic covariates by counterfactual values that reflect severely adverse macroeconomic conditions. We find that our proposed posterior mean predictor is considerably more accurate than a predictor that does not utilize any prior distribution. The posterior mean predictor shrinks the estimates of the unit-specific coefficients toward a common prior mean, which reduces its sampling variability. According to our estimates, the effect of stressed macroeconomic conditions on bank revenues is very small relative to the cross-sectional dispersion of revenues across holding companies.
Our paper is related to several strands of the literature. For $\alpha=\rho=0$ and $W_{it}=1$ the problem analyzed in this paper reduces to the problem of estimating a vector of means, which is a classic problem in the statistic literature. In this context, Tweedie's formula has been used, for instance, by Robbins1951 and more recently by BrownGreenshtein2009 and Efron2011 in a “big data” application. Throughout this paper we are adopting an empirical Bayes approach, that uses cross-sectional information to estimate aspects of the prior distribution of the correlated random effects and then conditions on these estimates. Empirical Bayes methods also have a long history in the statistics literature going back to Robbins1955 (see Robert1994 for a textbook treatment).
We use compound decision theory as in Robbins1964, BrownGreenshtein2009, JiangZhang2009 to state our optimality result. Because our setup nests the linear dynamic panel data model, we utilize results on the consistent estimation of $\rho$ in dynamic panel data models with fixed effects when $T$ is small, e.g., AndersonHsiao1981, ArellanoBond1991, ArellanoBover1995, BlundellBond1998, AlvarezArellano2003. Fully Bayesian approaches to the analysis of dynamic panel data models have been developed in ChamberlainHirano1999, Hirano2002, Lancaster2002.
The papers that are most closely related to ours are GuKoenkerJAE2016,GuKoenker2014. They also consider a linear panel data model and use Tweedie's formula to construct an approximation to the posterior mean of the heterogeneous regression coefficients. However, their papers focus on the use of the Kiefer-Wolfowitz estimator for the cross-sectional distribution of the sufficient statistics, whereas our paper explores various plug-in estimators for the homogeneous coefficients in combination with both parametric and nonparametric estimates of the cross-sectional distribution. Moreover, our paper establishes the ratio-optimality of the forecast and presents a different application. Finally, Liu2016 develops a fully Bayesian (as opposed to empirical Bayes) approach to construct density forecast. She uses a Dirichlet process mixture to construct a prior for the distribution of the heterogeneous coefficients, which then is updated in view of the observed panel data.
There is an earlier panel forecast literature (e.g., see the survey article by Baltagi2008 and its references) that is based on the best linear unbiased prediction (BLUP) proposed by Goldberger1962. Compared to the BLUP-based forecasts, our forecasts based on Tweedie's formula have several advantages. First, it is known that the estimator of the unobserved individual heterogeneity parameter based on the BLUP method corresponds to the Bayes estimator based on a Gaussian prior (see, for example, Robinson1991), while our estimator based on Tweedie's formula is consistent with much more general prior distributions. Second, the BLUP method finds the forecast that minimizes the expected quadratic loss in the class of linear (in $(Y_{i0},...,Y_{iT})'$) and unbiased forecasts. Therefore, it is not necessarily optimal in our framework that constructs the optimal forecast without restricting the class of forecasts. Third, the existing panel forecasts based on the BLUP were developed for panel regressions with random effects and do not apply to correlated random effects settings.
There is a small academic literature on econometric techniques for stress test. Most papers analyze revenue and balance sheet data for the relatively small set of bank holding companies with consolidated assets of more than 50 billion dollars. There are slightly more than 30 of these companies and they are subject to the Comprehensive Capital Analysis and Review conducted by the Federal Reserve Board of Governors. An important paper in this literature is CovasRumpZakrajsek2014, which uses quantile autoregressive models to forecast bank balance sheet and revenue components. We work with a much larger panel of bank holding companies that comprises, depending on the sample period, between 460 and 725 institutions.
The remainder of the paper is organized as follows. Section (ref) introduces the panel data model considered in this paper, derives the likelihood function, and provides an important identification result. Decision theoretic foundations for the proposed predictor and a derivation of the oracle forecast are provided in Section (ref). Section (ref) discusses feasible implementation strategies for the predictor and we show in Section (ref) in the context of a basic dynamic panel data model that our proposed predictor asymptotically has the same risk as the oracle forecast. A simulation study is provided in Section (ref). The empirical application is presented in Section (ref) and Section (ref) concludes. Technical derivations, proofs, the description of the data set used in the empirical analysis, and further empirical results are relegated to the Appendix.
We consider a panel with observations for cross-sectional units $i=1,\ldots,N$ in periods $t=1,\ldots,T$. Observation $Y_{it}$ is assumed to be generated by ((ref)). We distinguish three types of regressors. First, the $k_w\times 1$ vector $W_{it}$ interacts with the heterogeneous coefficients $\lambda_i$. In many panel data applications $W_{it} = 1$, meaning that $\lambda_i$ is simply a heterogenous intercept. We allow $W_{it}$ to also include deterministic time effects such as seasonality, time trends and/or strictly exogenous variables observed at time $t$. To distinguish deterministic time effects $w_{1,t+1}$ from cross-sectionally varying and strictly exogenous variables $W_{2,it}$, we partition the vector into $W_{it} = (w_{1,t+1},W_{2,it})$.\footnote{Because $W_{it}$ is a predictor for $Y_{it+1}$ we use a $t+1$ subscript for the deterministic trend component $w_{1}$.} The dimensions of the two components are $k_{w_1}$ and $k_{w_2}$, respectively. Second, $X_{it}$ is a $k_x \times 1$ vector of sequentially exogenous predictors with homogeneous coefficients. The predictors $X_{it}$ may include lags of $Y_{it+1}$ and we collect all the predetermined variables other than the lagged dependent variable into the subvector $X_{2,it}$. Third, $Z_{it}$ is a $k_z$-vector of strictly exogenous regressors, also with common coefficients.
Our main goal is to construct optimal forecasts of $(Y_{1T+1},...,Y_{NT+1})$ conditional on the entire panel observations $\{(Y_{it},W_{it-1},X_{it-1},Z_{it-1})$, $i = 1,\ldots,N$ and $t=1,...,T$ using the forecasting model ((ref)). An important special case of model ((ref)) is the basic dynamic panel data model
which is obtained by setting $W_{it} = 1$, $X_{it} = Y_{it}$ and $\alpha=0$. The restricted model ((ref)) has been widely studied in the literature. However, most studies focus on consistently estimating the common parameter $\rho$ in the presence of an increasing (with the cross-sectional dimension $N$) number of $\lambda_i$s. In forecasting applications, we also need to estimate the $\lambda_i$s. In Section (ref) we specify the likelihood function for model ((ref)) and in Section (ref) we establish the identifiability of the model parameters, including the distribution of the heterogeneous coefficients $\lambda_i$.
Let $Y_i^{t_1:t_2} =(Y_{it_1},...,Y_{it_2})$ and use a similar notation to collect $W_{it}$s, $X_{it}$s, and $Z_{it}s$. We begin by making some assumptions on the joint distribution of $\{Y_i^{1:T+1},X_i^{0:T},W_{2,i}^{0:T},Z_i^{0:T},\lambda_i\}_{i=1}^N$ conditional on the regression coefficients $\rho$ and $\alpha$ and the vector of volatility parameters $\gamma$ (to be introduced below). We drop the deterministic trend regressors $w_{1,t}$ from the notation for now. We use $\mathbb{E}[\cdot]$ to denote expectations and $\mathbb{V}[\cdot]$ to denote variances.
Assumption (ref)(i) states that conditionally on the predictors, the $Y_{it+1}$s are cross-sectionally independent. Thus, we assume that all the spatial correlation in the dependent variables is due to the observed predictors. Assumption (ref)(ii) formalizes the correlated random effects assumption. The subsequent Assumptions (ref)(iii) and (iv) imply that $\lambda_i$ may affect $X_{it}$ only indirectly through $Y_i^{1:t}$ -- an assumption that is clearly satisfied in the dynamic panel data model ((ref)) -- and that the strictly exogenous predictors do not depend on $\lambda_i$. In Assumption (ref)(v), we allow the unpredictable shocks $U_{it}$ to be conditionally heteroskedastic in both the cross section and over time. We allow $\sigma_t(\cdot)$ to be dependent on the initial condition of the sequentially exogenous predictors, $X_{i0}$, and other exogenous variables. Because throughout the paper we assume that the time dimension $T$ is small, the dependence through $X_{i0}$ can generate a persistent ARCH effect.
We now turn to the likelihood function. We use lower case $(y_{it}, w_{it}, x_{it}, z_{it})$ to denote the realizations of the random variables $(Y_{it},X_{it},W_{it},Z_{it})$. The parameters that control the volatilities $\sigma_t(\cdot)$ are stacked into the vector $\gamma=[\gamma_1',...,\gamma_T']'$ and we collect the homogeneous parameters into the vector $\theta = [\alpha',\rho',\gamma']'$. We use $H_i = (X_{i0},W_{2,i}^{0:T},Z_i^{0:T})$ for the exogenous conditioning variables and $h_i = (x_{i0},w_{2,i}^{0:T},z_i^{0:T})$ for their realization. Finally, we denote the density of $V_i$ by $\varphi(v)$. Recall that we used $x_{2,it}$ to denote predetermined predictors other than the lagged dependent variable. According to Assumption (ref)(iii) the density $q_{t}(x_{2,it}|y_i^{1:t}, x_i^{0:t-1},w_{2i},z_i)$ does not provide any information about $\lambda_i$ and will subsequently be absorbed into a constant of proportionality. Combining the likelihood function for the observables with the conditional distribution of the heterogeneous coefficients leads to
Because conditional on the predictors the observations are cross-sectionally independent, the joint densities for observations $i=1,\ldots,N$ can be obtained by taking the product across $i$ of ((ref)).
We now provide conditions under which the forecasting model ((ref)) is identifiable. While the identification of the finite-dimensional parameter vector $\theta$ is fairly straightforward, the empirical Bayes approach pursued in this paper also requires the identification of the correlated random effects distribution $\pi(\lambda_i|h_i)$ from the cross-sectional information in the panel. Before presenting a general result which is formally proved in the Online Appendix, we sketch the identification argument in the context of the restricted dynamic model ((ref)) with heterogeneous intercept and heteroskedastic innovations.
The identification can be established in three steps. First, the identification of the homogeneous regression coefficient $\rho$ follows from a standard argument used in the instrumental variable (IV) estimation of dynamic panel data models. To eliminate the dependence on $\lambda_i$ define $Y_{it}^* = Y_{it} - \frac{1}{T-t} \sum_{s=t+1}^T Y_{is}$ and $X_{it-1}^* = Y_{it-1} - \frac{1}{T-t} \sum_{s=t+1}^T Y_{is-1}$. Then, because $\mathbb{E}[U_{it}|Y_i^{0:t-1},\lambda_i] = 0$, the orthogonality conditions $\mathbb{E}\big[ (Y_{it}^* - \rho X_{it-1}^*) Y_{it-1} \big] = 0$ for $t=1,\ldots,T-1$ in combination with a relevant rank condition can be used to identify $\rho$ (see, e.g., ArellanoBover1995). Second, to identify the variance parameters $\gamma$, let $Y_i$, $X_{i}$, and $U_i$ denote the $T \times 1$ vectors that stack $Y_{it}$, $Y_{it-1}$, and $U_{it}$, respectively, for $t=1,\ldots,T$. Moreover, let $\iota$ be a $T\times 1$ vector of ones and define $\Sigma_i^{1/2}(\tilde{\gamma}) = \mbox{diag}\big( \sigma_1(h_i, \tilde{\gamma}_1), \ldots, \sigma_T(h_i, \tilde{\gamma}_T) \big)$, $S_i(\tilde{\gamma}) = \Sigma_i^{-1/2}(\tilde{\gamma}) \iota$, and $M_i(\tilde{\gamma}) = I - S_i(S_i'S_i)^{-1} S_i'$. Using this notation, we obtain
This leads to the conditional moment condition
if and only if $\tilde{\gamma} = \gamma$, which identifies $\gamma$. Third, let
The identification of $\pi(\lambda_i|h_i)$ can be established using a characteristic function argument similar to that in ArellanoBonhomme2012. For the general model ((ref)) we make the following assumptions:
Because the identification of $\alpha$ and $\rho$ in panel data models with fixed or random effects is well established, we make the high-level Assumption (ref)(i) that the homogeneous parameters are identifiable.\footnote{Textbook / handbook chapter treatments can be found in, for instance, Baltagi1995, ArellanoHonore2001, Arellano2003 and Hsiao2014.} We discuss in the appendix how the identification argument for $\rho$ in the basic dynamic panel data model can be extended to a more general specification as in ((ref)). Assumption (ref)(ii) enables us to identify the volatility parameters $\gamma$, and (iii) and (iv) deliver the identifiability of the distribution of heterogeneous coefficients. The following theorem summarizes the identification result and is proved in the Appendix.
We adopt a decision-theoretic framework in which forecasts are evaluated based on cross-sectional sums of mean-squared error losses. Such losses are called compound loss functions. Section (ref) provides a formal definition of the compound risk (expected loss). In Section (ref) we derive the optimal forecasts under the assumption that the cross-sectional distribution of the $\lambda_i$s is known (oracle forecast). While it is infeasible to implement this forecast in practice, the oracle forecast provides a natural benchmark for the evaluation of feasible predictors. Finally, in Section (ref) we introduce the concept of ratio optimality, which describes forecasts that asymptotically (as $N \longrightarrow \infty$) attain the same risk as the oracle forecast.
Let $L( \widehat{Y}_{iT+1},Y_{iT+1}) $ denote the loss associated with forecast $\hat{Y}_{i,T+1}$ of individual $i^{\prime }s$ time $T+1$ observation, $Y_{iT+1}$. In this paper we consider the conventional quadratic loss function,
The main goal of the paper is to construct optimal forecasts for groups of individuals selected by a known selection rule in terms of observed data. We express the selection rule as
where $D_i({\cal Y}^N)$ is a measurable function of the observations ${\cal Y}^N$, ${\cal Y}^N = ({\cal Y}_1,\ldots,{\cal Y}_N)$, and ${\cal Y}_i = (Y_i^{0:T},X_i^{1:T},H_i )$. For instance, suppose that $D_i({\cal Y}^N) = \mathbb{I}\{Y_{iT} \in A \}$ for $A \subset \mathbb{R}$. In this case, the selection is homogeneous across $i$ and, for individual $i$, depends only on its own sample. Alternatively, suppose that units are selected based on the ranking of an index, e.g., the empirical quantile of $Y_{iT}$. In this case, the selection dummy $D_i$ depends on $(Y_{1T},...,Y_{NT})$ and thereby also on the data for the other $N-1$ individuals.
The compound loss of interest is the average of the individual losses weighted by the selection dummies:
where $Y^N_{T+1} = (Y_{1T+1},\ldots,Y_{NT+1})$. The compound risk is the expected compound loss
We use the $\theta$ subscript for the expectation operator to indicate that the expectation is conditional on $\theta$.\footnote{Strictly speaking, the expectation also conditions on the deterministic trend terms $W_1$}. The superscript $({\cal Y}^N,\lambda^N,U^N_{T+1})$ indicates that we are integrating with respect to the observed data ${\cal Y}^N$ and the unobserved heterogeneous coefficients $\lambda^N = (\lambda_1,\ldots,\lambda_N)$ and $U^N_{T+1}=(U_{1T+1},\ldots,U_{NT+1})$.
We now derive the optimal forecast that minimizes the compound risk. The risk achieved by the optimal forecast will be called the oracle risk, which is the target risk to achieve. In the compound decision theory it is assumed that the oracle knows the vector $\theta$ as well as the distribution of the heterogeneous coefficients $\pi(\lambda_i,h_i)$ and observes ${\cal Y}^N$. However, the oracle does not know the specific $\lambda_i$ for unit $i$. In order to find the optimal forecast, note that conditional on $\theta$ the compound risk takes the form of an integrated risk that can be expressed as
The inner expectation can be interpreted as posterior risk, which is obtained by conditioning on the observations ${\cal Y}^N$ and integrating over the heterogeneous parameter $\lambda^N$ and the shocks $U^N_{T+1}$. The outer expectation averages over the possible trajectories ${\cal Y}^N$.
It is well known that the integrated risk is minimized by choosing the forecast that minimizes the posterior risk for each realization ${\cal Y}^N$. Using the independence across $i$, the posterior risk can be written as follows:
where $\mathbb{V}_{\theta,{\cal Y}_i}^{\lambda_i,U_{iT+1}}[\cdot]$ is the posterior variance. The decomposition of the risk into a squared bias term and the posterior variance of $Y_{iT+1}$ implies that $\mathbb{E}_{\theta,{\cal Y}_i}^{\lambda_i,U_{iT+1}} [Y_{iT+1}]$ is the optimal predictor. Because $U_{iT+1}$ is mean-independent of $\lambda_i$ and ${\cal Y}_i$, we obtain
Note that the posterior expectation of $\lambda_i$ only depends on observations for unit $i$, even if the selection rule $D_i({\cal Y}^N)$ also depends on the data from other units $j \not=i$. The result is summarized in the following theorem:
According to ((ref)), the compound oracle risk has two components. The first component reflects uncertainty with respect to the heterogeneous coefficient $\lambda_i$ and the second component captures uncertainty about the error term $U_{iT+1}$. Unfortunately, the direct implementation of the optimal forecast is infeasible because neither the parameter vector $\theta$ nor the correlated random effect distribution (or prior) $\pi(\cdot)$ are known. Thus, the oracle risk $R_{N}^{\text{opt}}$ provides a lower bound for the risk that is attainable in practice.
The identification result presented in Section (ref) implies that as the cross-sectional dimension $N \longrightarrow \infty$, it might be possible to learn the unknown parameter $\theta$ and random-effects distribution $\pi(\cdot)$ and construct a feasible estimator that asymptotically attains the oracle risk. Following BrownGreenshtein2009, we say that a predictor achieves ratio optimality if the regret $R_{N}(\widehat{Y}^N_{T+1}) - R_N^{\text{opt}}$ of the forecast $\widehat{Y}^N_{T+1}$ is negligible relative to the part of the optimal risk that is due to uncertainty about $\lambda_i$:
Using ((ref)), the risk differential in the numerator (called regret) can be written as
For illustrative purposes, Consider the basic dynamic panel data model ((ref)). For this model $\mathbb{E}_{\theta,{\cal Y}_i}^{\lambda_i,U_{iT+1}}[Y_{iT+1}] = \mathbb{E}_{{\cal Y}_i}^{\lambda_i}[\lambda_i] + \rho Y_{iT}$. A natural class of predictors is given by $\widehat{Y}_{iT+1} = \widehat{\mathbb{E}}_{{\cal Y}_i}^{\lambda_i}[\lambda_i] + \hat{\rho} Y_{iT}$, where $\widehat{\mathbb{E}}_{{\cal Y}_i}^{\lambda_i}[\lambda_i]$ is an approximation of the posterior mean of $\lambda_i$ that replaces the unknown $\rho$ and distribution $\pi(\cdot)$ by suitable estimates. The autoregressive coefficient in this model can be $\sqrt{N}$-consistently estimated, which suggests that $\sum_{i=1}^N (\hat{\rho}-\rho)^2Y_{iT}^2 = O_p(1)$. Thus, whether a predictor attains ratio optimality crucially depends on the rate at which the discrepancy between $\mathbb{E}_{{\cal Y}_i}^{\lambda_i}[\lambda_i]$ and $\widehat{\mathbb{E}}_{{\cal Y}_i}^{\lambda_i}[\lambda_i]$ vanishes.
The denominator of the ratio in Definition (ref) is divergent. The rate of divergence depends on the posterior variance of $\lambda_i$. If the posterior variance is strictly greater than zero, then the denominator is of order $O(N)$. Note that for each unit $i$, the posterior variance is based on a finite number of observations $T$. Thus, for the posterior variance to be equal to zero, it must be the case that the prior density $\pi(\lambda)$ is a pointmass, meaning that there is a homogeneous intercept $\lambda$. In this case the definition of ratio optimality requires that the regret vanishes at a faster rate, because the rate of the numerator drops from $O(N)$ to $N^{\epsilon_0}$. Subsequently, we will pursue an empirical Bayes strategy to construct an approximation $\widehat{\mathbb{E}}_{{\cal Y}_i}^{\lambda_i}[\lambda_i]$ based on the cross-sectional information and show that it attains ratio-optimality.
In the linear panel literature, researchers often use the first difference to eliminate $\lambda_i$. In this case, the natural forecast of $Y_{iT+1}$ in the basic dynamic panel data model ((ref)) would be $\widehat{Y}_{iT+1}^{FD}(\rho) = Y_{iT} + \rho (Y_{iT} - Y_{iT-1})$, which is different from $\widehat{Y}_{iT+1}^{opt}$ in ((ref)). Thus, we can immediately deduce from Theorem (ref) that $\widehat{Y}_{iT+1}^{FD}(\rho)$ is not an optimal forecast. The quasi-differencing of $Y_{it}$ introduces a predictable moving-average error term that is ignored by the predictor $\widehat{Y}_{iT+1}^{FD}(\rho)$.
We will construct a consistent approximation of the posterior mean $\mathbb{E}_{\theta,{\cal Y}^i}^{\lambda_i,U_{iT+1}}[\lambda_i]$ using a convenient formula which is named after the statistician Maurice Tweedie (though it had been previously derived by the astronomer Arthur Eddington). This formula is presented in Section (ref). In Section (ref) we discuss the parametric estimation of the correction term and in Section (ref) we consider a nonparametric kernel-based estimation. The QMLE and Generalized Method-of-Moments (GMM) estimation of the parameter $\theta$ are discussed in Sections (ref) and (ref).
When the innovations $U_{it}$ are conditionally normally distributed, we can derive a convenient formula for the posterior expectation $\mathbb{E}_{\theta,{\cal Y}_i}^{\lambda_i} [\lambda_i]$ of the individual heterogeneous parameter $\lambda_i$.
The assumption of normally distributed $V_{it}$'s is not as restrictive as it may seem. Recall that the shocks $U_{it}$ are defined as $V_{it} \sigma_t(X_{i0},W_{2,i}^{0:T},Z_i^{0:T},\gamma_t)$. Thus, due to the potential heteroskedasticity, the distribution of shocks is a mixture of normals. The only restriction is that the random variables characterizing the scale of the mixture component are observed. Moreover, even in the homoskedastic case $\sigma_t = \sigma$, the distribution of $Y_{it}$ given the regressors is non-normal because the distribution of the $\lambda_i$ parameters is fully flexible. Using Assumption (ref) we will now further manipulate the density $p(y_i,x_{2,i},\lambda_i|h_i,\theta)$ in ((ref)).\footnote{In principle, the normality assumption could be generalized to the assumption that the distribution of $V_{it}$ belongs to the exponential family.} To simplify the notation we will drop the $i$ subscript. Define
and let $\tilde{y}(\theta)$ and $w$ be matrices with rows $\tilde{y}_t(\theta)$ and $w_{t-1}'$, $t=1,...,T$. Because the subsequent calculations condition on $\theta$ we will omit the $\theta$-argument from $\tilde{y}$, $\Sigma$, and functions thereof. Replacing $\varphi(v)$ in ((ref)) with a Gaussian density function we obtain:
The factorization of $p(y,x_2,\lambda|h,\theta)$ implies that
is a sufficient statistic and that we can express the posterior distribution of $\lambda$ as \[ p(\lambda|y,x_2,h,\theta) = p(\lambda|\hat{\lambda},h,\theta) = \frac{ p(\hat{\lambda}|\lambda,h,\theta) \pi(\lambda|h) } { p(\hat{\lambda}| h,\theta) }, \] where
To obtain a representation for the posterior mean, we now differentiate the equation $\int p(\lambda|\hat{\lambda}, h,\theta) d\lambda = 1$ with respect to $\hat{\lambda}$. Exchanging the order of integration and differentiation and using the properties of the exponential function, we obtain
Solving this equation for the posterior mean yields Tweedie's formula, which is summarized in the following theorem.
Tweedie's formula was used by Robbins1951 to estimate a vector of means $\lambda^N$ for the model $Y_i|\lambda_i \sim N(\lambda_i,1)$, $\lambda_i \sim \pi(\cdot)$, $i=1,\ldots,N$. Recently, it was extended by Efron2011 to the family of exponential distribution, allowing for a unknown finite-dimensional parameter $\theta$. Theorem (ref) extends Tweedie's formula to the estimation of correlated random effect parameters in a dynamic panel regression setup.
The posterior mean takes the form of the sum of the sufficient statistic $\hat{\lambda}_i(\theta)$ and a correction term that reflects the prior distribution of $\lambda_i$. The correction term is expresses as a function of the marginal density of the sufficient statistic $\hat{\lambda}_i(\theta)$ conditional on $H_i$ and $\theta$. Thus, it is not necessary to solve a deconvolution problem that separates the prior density $\pi(\lambda_i|h_i)$ from the distribution of the error terms $V_{it}$. We expressed Tweedie's formula in ((ref)) in terms of the conditional density $p(\hat{\lambda}_i(\theta)|H_i,\theta)$. However, because the posterior mean is a function of the log density differentiated with respect to $\hat{\lambda}_i(\theta)$, the conditional density can be replaced by a joint density: \[ \frac{\partial}{\partial \hat{\lambda}_i(\theta)} \ln p(\hat{\lambda}_i(\theta)|H_i,\theta) = \frac{\partial}{\partial \hat{\lambda}_i(\theta)} \ln p(\hat{\lambda}_i(\theta),H_i|\theta). \] The construction of ratio-optimal forecasts relies on replacing the density $p(\hat{\lambda}_i(\theta),H_i|\theta)$ and the common parameter $\theta$ by consistent estimates.
If the random-effects distribution $\pi(\lambda|h_i)$ is Gaussian, then it is possible to derive the marginal density of the sufficient statistic $p(\hat{\lambda}_i(\theta)|h_i,\theta)$ analytically. Let
Moreover, define $\xi = \big( \mbox{vec}(\Phi), \, \mbox{vech}(\underline{\Omega}) \big)'$. To highlight the dependence of the correlated random-effects distribution on the hyperparameter $\xi$ we will write $\pi(\lambda_i|h_i,\xi)$. The marginal density (omitting the $i$ subscripts and the $\theta$-argument of $\hat{\lambda}$) is given by
Here, we used the likelihood of $\hat{\lambda}$ in ((ref)), the density associated with the Gaussian prior in ((ref)), and then the properties of a multivariate Gaussian density to integrate out $\lambda$. The terms $\bar{\lambda}$ and $\bar{\Omega}$ are the posterior mean and variance of $\lambda$, respectively: \[ \bar{\Omega}^{-1} = \underline{\Omega}^{-1} + w'\Sigma^{-1} w, \quad \bar{\lambda} = \bar{\Omega} \big( \underline{\Omega}^{-1} \Phi h + w'\Sigma^{-1} w \hat{\lambda} \big). \]
Conditional on $\theta$ the vector of hyperparameters $\xi$ can be estimated by maximizing the marginal likelihood
using the cross-sectional distribution of the sufficient statistic. Tweedie's formula can then be evaluated based on $p \big(\hat{\lambda}_i(\theta)|h_i,\theta,\hat{\xi}(\theta) \big)$. In principle it is possible to replace the Gaussian prior distribution with a more general parametric distribution. However, in general it will not be possible to derive an analytical formula for the marginal likelihood.
A nonparametric implementation of the Tweedie correction can be obtained by replacing $p(\hat{\lambda}_i(\theta),h_i|\theta)$ and its derivative with respect to $\hat{\lambda}_i(\theta)$ with a Kernel density estimate, e.g.,
where $B_N$ is the bandwidth and $V_{\hat{\lambda}}$ and $V_h$ are tuning matrices. Note that even if the prior distribution $\pi(\lambda)$ is a pointmass, the sufficient statistic $\hat{\lambda}$ in ((ref)) has a continuous distribution and one can use a kernel density estimator to construct the Tweedie correction.
If the dimension of the conditioning variables $H_i$ is large, the nonparametric estimation suffers from the curse of dimensionality. In this case, one may reduce the dimension of the conditioning set with some smaller dimensional indices, e.g., by assuming that $\lambda_i$ and $H_i$ dependent only through $\bar{H}_{i} = \frac{1}{T} \sum_{t=1}^T H_{it}$, that is, $\pi(\lambda|h) = \pi(\lambda|\bar{h})$. In Section (ref) we provide a detailed analysis of the Gaussian kernel estimator in the context of the basic dynamic panel data model in ((ref)) with time-homoskedastic innovations.
Notice that under Assumption (ref), $\hat{\lambda}_i(\theta)$ in ((ref)) is a sufficient statistic of $\lambda_i$ conditional on $\theta,h_i$, and $ \pi_{\lambda}(\lambda_i|h_i,\xi)$ is the parametric version of the correlated random effect density. Integrating out $\lambda$ under a parametric correlated random effect (or prior) distribution $\pi_{\lambda}(\lambda|x_{0},w_2,z,\xi)$, we have (omitting the $i$ subscripts)
Here, we used the definition of $\tilde{y}(\theta)$ in ((ref)) and the product of Gaussian likelihood and prior in ((ref)). Note that the term $p(\hat{\lambda}(\theta)|h,\theta,\xi)$ in the last line of ((ref)) is identical to the objective function for $\xi$ used in ((ref)). Thus, we can now jointly determine $\theta$ and $\xi$ by maximizing the integrated likelihood as a function:
We refer to this estimator as {\em quasi} (Q) maximum likelihood estimator (MLE), because the correlated random effects distribution could be misspecified.
Without a convenient assumption about the random effects distribution, one can estimate the parameter $\theta$ using a sample analogue of the moment conditions that were used in the identification analysis in Section (ref). For $t=1,\ldots,T-k_w$, define
Moreover, define $X_{it-1}^{*}$ and $Z_{it-1}^{*}$ by replacing $Y_{i\cdot}$ in ((ref)) with $X_{i\cdot}$ and $Z_{i\cdot}$, respectively, and let \[ g_{it}(\rho,\alpha) = (Y_{it}^{*} - \rho^{\prime}X_{it-1}^{*} - \alpha^{\prime}Z_{it-1}^{*}) \left[
\right], \quad g_i(\rho,\alpha) = \big[ g_{i1}(\rho,\alpha)', \ldots ,g_{iT-k_{w}}(\rho,\alpha)' \big]'. \] The continuous-updating GMM estimator of $\rho$ and $\alpha$ solves
This estimator was proposed by ArellanoBover1995 and we will refer to it as GMM(AB) estimator in the Monte Carlo simulations (Section (ref)) and the empirical application (Section (ref)).\footnote{There exists a large literature on the estimation of dynamic panel data models. Alternative estimators include ArellanoBond1991 and BlundellBond1998.}
To estimate the heteroskedasticity parameter $\gamma = [\gamma_1,...,\gamma_T]'$ in $\sigma^2_t(H_i,\gamma_t)$, define:
where $\hat{\rho}$ and $\hat{\alpha}$ could be the estimators in ((ref)). We use the sample analogue to a set of moment condition implied by a generalization of ((ref)):
where $B$ is a selection matrix that can be used to eliminate off-diagonal elements of the covariance matrix. In population, these off-diagonal elements should be zero, because the $U_{it}$'s are assumed to be uncorrelated across time.
While this paper focuses on single-step forecasting, we briefly discuss in the context of the basic dynamic panel data model how the framework can be extended to multi-step forecasts. We can express \[ Y_{iT+h} = \left(\sum_{s=0}^{h-1} \rho^s\right) \lambda_i + \rho^h Y_{iT} + \sum_{s=0}^{h-1} \rho^2 U_{iT+h-s}. \] Under the assumption that the oracle knows $\rho$ and $\pi(\lambda_i,Y_{i0})$ we can express the oracle forecast as \[ \widehat{Y}^{opt}_{iT+h} = \left(\sum_{s=0}^{h-1} \rho^s\right) \mathbb{E}_{\theta,{\cal Y}_i}^{\lambda_i}[\lambda_i] + \rho^h Y_{iT}. \] As in the case of the one-step-ahead forecasts, the posterior mean $\mathbb{E}_{\theta,{\cal Y}_i}^{\lambda_i}[\lambda_i]$ can be replaced by an approximation based on Tweedie's formula and the $\rho$'s can be replaced by consistent estimates. A model with additional covariates would require external multi-step forecasts of the covariates, or the specification in ((ref)) would have to be modified such that all exogenous regressors appear with an $h$-period lag.
Throughout this section we will consider the basic dynamic panel data model with homoskedastic Gaussian innovations:
We will prove that ratio optimality for a general prior density $\pi(\lambda_i|h_i)$ can be achieved with a Kernel estimator of the joint density of the sufficient statistic and initial condition: $p(\hat{\lambda}_i(\theta),H_i|\theta)$. The proof of the main result is a significant generalization of the proof in BrownGreenshtein2009 for a vector of means to the dynamic panel data model with estimated common coefficients.
For the model in ((ref)), the sufficient statistic is given by
and the posterior mean of $\lambda_i$ simplifies to
The formula recognizes that the heterogeneous coefficient is a scalar intercept and that the errors are homoskedastic. We simplified the notation by writing $p(\hat{\lambda}_i(\rho),Y_{i0})$ instead of $p(\hat{\lambda}_i(\rho),Y_{i0}|\theta)$. This simplification is justified because we will estimate the density of $(\hat{\lambda}_i(\rho),Y_{i0})$ directly from the data; see ((ref)) below. We will use the notation $\mu(\cdot)$ to refer to the conditional mean as function of the sufficient statistic $\hat{\lambda}$, the scale factor $\sigma^2/T$, and the density $p(\hat{\lambda}_i,Y_{i0})$.
To facilitate the theoretical analysis, we make two adjustments to the posterior mean predictor of $Y_{iT+1}$. First, we replace the kernel density estimator of $(\hat{\lambda}_i(\rho),Y_{i0})$ given in ((ref)) by a leave-one-out estimator of the form:
where $\phi(\cdot)$ is the pdf of a $N(0,1)$. Using the fact that the observations are cross-sectionally independent and conditionally normally distributed one can directly compute the expected value of the leave-one-out estimator:
Taking expectations of the kernel estimator leads to a variance adjustment for conditional distribution of $\hat{\lambda}_i|\lambda_i$ ($\sigma^2/T+B_N^2$ instead of $\sigma^2/T$) and the density of $y_{i0}|\lambda_i$ is replaced by a convolution.
Second, we replace the scale factor $\hat{\sigma}^2/T$ in the posterior mean function $\mu(\cdot)$ by $\hat{\sigma}^2/T+B_N^2$, which is the term that appears in ((ref)). Moreover, we truncate the absolute value of the posterior mean function from above. For $C > 0$ and for any $x \in \mathbb{R}$, define $\left[x\right]^C := \operatorname*{sgn}(x) \min \{ |x|,C\}$. Then
where $C_N \longrightarrow \infty$ slowly. Formally, we make the following technical assumptions.
We factorize the correlated random effects distribution as $\pi(\lambda_i,y_{i0}) = \pi(\lambda_i)\pi(y_{i0}|\lambda_i)$ and impose regularity conditions on the marginal distribution of the heterogeneous coefficient and the conditional distribution of the initial condition. In Assumption (ref) we let the support of $\pi(\lambda_i)$ slowly expand with the sample size by assuming that $C_N$ grows at a subpolynomial rate. Assumption (ref) provides an upper and a lower bound for the rate at which the bandwidth of the kernel estimator shrinks to zero. Note that for technical reasons the assumed rate is much slower than in typical density estimation problems.\footnote{In a nutshell, we need to control the behavior of $\hat{p}(\hat{\lambda}_i,Y_{i0})$ and its derivative uniformly, which, in certain steps of the proof, requires us to consider bounds of the form $M/B_N^2$, where $M$ is a generic constant. If the bandwidth shrinks too fast, the bounds diverge too quickly to ensure that it suffices to standardize the regret in Definition (ref) by $N^{\epsilon_0}$ if the $\lambda_i$ coefficients are identical for each cross-sectional unit.}
Assumption (ref) imposes regularity conditions on the conditional density of the initial observation. In (i) we assume that $\pi(y_{i0}|\lambda_i)$ is bounded. In (ii) we control the tails of the distribution. In the first constraint on $m(C,\lambda)$ we essentially assume that the density of $y_{i0}$ has exponential tails. This also guarantees that the fourth moment of $Y_{i0}$ exists. In part (iii) we assume that $\pi(y|\lambda)$ is sufficiently smooth with respect to $y$ such that the convolution on the left-hand side uniformly converges to $\pi(y|\lambda)$ as the bandwidth $B_N$ tends to zero. We verify in the Appendix that a $\pi(y|\lambda)$ that satisfies Assumption (ref) is $\pi(y|\lambda) = \phi( y - \lambda)$, where $\phi(x) = \exp(-\frac{1}{2}x^2)/\sqrt{2\pi}$. Finally, Assumption (ref) postulates the existence of finite sample moments of the estimators of the common parameter. The main result is stated in the following theorem:
The result in Theorem (ref) is pointwise with respect to $\theta$. However, the convergence of the predictor $\widehat{Y}_{iT+1}$ to the oracle predictor is uniform with respect to the unobserved heterogeneity and the observed trajectory ${\cal Y}_i$ in the sense that the integrated risk (conditional on $\theta$) of the feasible predictor converges to the integrated risk of the oracle predictor. The proof of the theorem is a generalization of the proof in BrownGreenshtein2009, allowing for the presence of estimated parameters in the sufficient statistic $\hat{\lambda}(\cdot)$. The remarkable aspect of the results is the acceleration of the convergence ($N^\epsilon_0$ instead of $N$ in the denominator of the standardized regret in Definition (ref)) in cases in which the intercepts are identical across units and $\pi(\lambda)$ is a pointmass.
We will now conduct several Monte Carlo experiments to illustrate the performance of the empirical Bayes predictor.
The first Monte Carlo experiment is based on the basic dynamic panel data model in ((ref)). The design of the experiment is summarized in Table (ref). We assume that the $\lambda_i$'s are normally distributed and uncorrelated with the initial condition $Y_{i0}$. The innovations $U_{it}$ and the heterogeneous intercepts $\lambda_i$ have unit variances. We consider two values for the autocorrelation parameter: $\rho \in \{0.5, 0.95\}$. The panel consists of $N=1,000$ cross-sectional units and the number of time periods is $T=3$. Generally, the smaller $T$ relative to number of right-hand-side variables with heterogeneous coefficients, the larger the gain from using a prior distribution to compute posterior mean estimates of the $\lambda_i$'s. We will compare the performance of the following predictors:
{\bf Oracle Forecast.} The oracle knows the parameters $\theta = (\rho,\gamma)$ as well as the random effects distribution $\pi(\lambda_i|Y_{i0},\xi)$, where $\xi = (\phi_0,\phi_1,\underline{\Omega})$. However, the oracle does not know the specific $\lambda_i$ values. Its forecast is given by ((ref)).
{\bf Posterior Predictive Mean Approximation Based on QMLE}. The random effects distribution is correctly modeled as belonging to the family $\lambda_i|(Y_{i0},\xi) \sim N(\phi_0+\phi_1Y_{i0},\underline{\Omega})$. The estimators $\hat{\theta}_{QMLE}$ and $\hat{\xi}_{QMLE}$ are defined in ((ref)). Tweedie's formula (see ((ref)) for the simplified version) is evaluated based on $p\big(\hat{\lambda}_i(\hat\theta_{QMLE})|y_{i0},\hat\theta_{QMLE},\hat{\xi}_{QMLE} \big)$.
{\bf Posterior Predictive Mean Approximation Based on GMM Estimator}. We use the Arellano-Bover estimator described in Section (ref). The estimator for $\rho$ is given by ((ref)) and the estimator for $\gamma$ by ((ref)). The formulas simplify considerably. We have $W_{it}=1$, $X_{it-1}=Y_{it-1}$, $Z_{it-1}=\emptyset$ and $\alpha = \emptyset$. Moreover, $\Sigma_i^{1/2} = \gamma I$, $M_i(\gamma) = I - \iota\iota'/T$, where $\iota$ is a $T\times 1$ vector of ones. Let $\bar{\tilde{Y}}_i(\hat{\rho})$ be the temporal average of $\tilde{Y}_i(\hat{\rho})$. Then \[ \hat{\gamma}_{GMM}^2 = \frac{1}{NT} \frac{T}{T-1} \sum_{i=1} \mbox{tr}\big[ (\tilde{Y}_i(\hat{\rho}) - \iota \bar{\tilde{Y}}_i(\hat{\rho}) )(\tilde{Y}_i(\hat{\rho}) - \iota \bar{\tilde{Y}}_i(\hat{\rho}) )' \big]. \] The estimator $\hat{\xi}(\hat{\theta}_{GMM})$ is obtained from ((ref)). Finally, Tweedie's formula is evaluated based on $p\big(\hat{\lambda}_i(\hat\theta_{GMM})|y_{i0},\hat\theta_{GMM},\hat{\xi}(\hat{\theta}_{GMM}) \big)$.
{\bf GMM Plug-In Predictor.} We use the Arellano-Bover estimator to obtain $\hat{\rho}_{GMM}$. Instead of using the posterior mean for $\lambda_i$, the plug-in predictor is based on the MLE $\hat{\lambda}_i(\hat{\rho}_{GMM})$. The resulting predictor is $\widehat{Y}_{iT+1} = \hat{\lambda}_i(\hat{\rho}_{GMM}) + \hat{\rho}_{GMM} Y_{iT}$.
{\bf Loss-Function-Based Predictor.} We construct an estimator of $(\rho,\lambda^N)$ based on the objective function:
This estimator minimizes the loss function under which the forecasts are evaluated in sample. It is well-known that due to the incidental parameter problem, the estimator $\hat{\rho}_L$ is inconsistent under fixed-$N$ asymptotics. The resulting predictor is $\widehat{Y}_{iT+1} = \hat{\lambda}_i(\hat{\rho}_L) + \hat{\rho}_L Y_{iT}$.
{\bf Pooled-OLS Predictor.} Ignoring the heterogeneity in the $\lambda_i$'s and imposing that $\lambda_i=\lambda$ for all $i$, we can define
The resulting predictor is $\widehat{Y}_{iT+1} = \hat{\lambda}_P + \hat{\rho}_P Y_{iT}$.
{\bf First-Difference Predictor.} In the panel data literature it is common to difference-out idiosyncratic intercepts, which suggests to predict $\Delta Y_{iT+1}$ based on $\Delta Y_{iT}$. We evaluate the first-difference predictor at the Arellano-Bover GMM estimator of $\rho$ to obtain $\widehat{Y}_{iT+1}^{FD}(\hat{\rho}_{GMM})$.
In Table (ref) we report the regret associated with each predictor relative to the posterior variance of $\lambda_i$, averaged over all trajectories ${\cal Y}^N$, as specified in Definition (ref) (setting $N^\epsilon=1$). For the oracle predictor the regret is by definition zero and we tabulate the risk $R_N^{opt}$ instead (in parentheses). We also report the median forecast error $\widehat{e}_{iT+1|T} = Y_{iT+1}-\widehat{Y}_{iT+1}$ to highlight biases in the forecasts.
\afterpage
The columns titled “All Units” correspond to $D_i({\cal Y}^N)=1$. As expected from the theoretical analysis, the posterior mean predictors have the lowest regret among the feasible predictors. The density of $\hat{\lambda}_i$ is estimated parametrically, using a family of distributions that nests the true random effects distribution. Because it is based on a correctly specified likelihood function, the predictor based on $\hat{\theta}_{QMLE}$ performs slightly better than the predictor based on $\hat{\theta}_{GMM}$. Consider $\rho=0.5$: for the QMLE-based predictor the regret is 0.5% of the average posterior variance, whereas it is 3% for the GMM-based predictor. The plug-in predictor that replaces the unknown $\lambda_i$'s by the sufficient statistic $\hat{\lambda}_i$ (which is also the maximum likelihood estimator) instead of the posterior mean is associated with a much larger relative regret, which is about 37%.
The remaining three predictors are also strictly dominated by the posterior mean predictors. Ignoring the serial correlation in $\Delta Y_{it}$, the first-difference predictor performs the worst for both choices of $\rho$. The second-to-worst predictor is the pooled-OLS predictor which ignores the cross-sectional heterogeneity in the $\lambda_i$'s. A reduction of the variance $\underline{\Omega}$ of the heterogeneous intercepts would improve the relative performance of the pooled-OLS predictor. Finally, the loss-function-based predictor dominates the pooled-OLS and the first difference predictor. As mentioned above, while conceptually appealing, the loss-function-based predictor relies on an inconsistent estimate of $\rho$, which in comparison to the GMM plug-in predictor is unappealing if the cross-sectional dimension $N$ is very large.
Across all units, the predictions under the loss-function-based estimator and the pooled-OLS estimator appear to be biased. To study this bias further we now consider level-based selection rules $D_i({\cal Y}^i)$. Using the 5%, 47.5%, 52.5%, and 95% quantiles of the population distribution of $Y_{iT}$, we define cut-offs for a bottom 5% group, a middle 5% group, and a top 5% group. Because the cut-offs are computed from the population distribution of $Y_{iT}$, for unit $i$ the selection rules only depends on ${\cal Y}_{iT}$ and not on $Y_{jT} $ with $j\not=i$.
For the top and bottom groups only the posterior mean predictors lead to unbiased forecast errors. The sufficient statistic $\hat{\lambda}_i$ tends to overestimate (underestimate) $\lambda_i$ for the top (bottom) group, because it interprets a sequence of above-average (below-average) $U_{iT}$'s as evidence for a high (low) $\lambda_i$. This is reflected in the bias: the plug-in predictors' forecast errors for the top group are on average positive, whereas the forecast errors for the bottom group tend to be negative. The posterior mean tends to correct these biases because it shrinks toward the mean of the prior distribution of the $\lambda_i$'s. This reduces the regrets for the top and bottom groups, and is also reflected in the risk calculated across all units. The bias correction is illustrated in Figure (ref), which compares the cross-sectional distribution of the sufficient statistics $\hat{\lambda}_i(\hat\theta)$ to the distribution of the posterior mean estimates $\widehat{\mathbb{E}}_{\hat{\theta},{\cal Y}_i}^{\lambda_i}[\lambda_i]$ obtained with Tweedie's formula. Due to the shrinkage effect of the prior, the distribution of the posterior means, in particular for the top and bottom groups, is more compressed.
We now change the Monte Carlo design in two dimensions. First, we replace the Gaussian random effects specification with a non-Gaussian specification in which the heterogeneous coefficient $\lambda_i$ is correlated with the initial condition $Y_{i0}$. Second, we consider a Tweedie correction based on a kernel density estimate of $p(\hat{\lambda}_i|Y_{i0})$ as discussed in Section (ref).
The Monte Carlo design is summarized in Table (ref). Starting point is a joint normal distribution for $(\lambda_i,Y_{i0})$, factorized into a marginal distribution $\pi_*(\lambda_i)$ and a conditional distribution $\pi_*(Y_{i0}|\lambda_i)$. We assumed $\lambda_i \sim N(\underline{\mu}_\lambda,\underline{V}_\lambda)$ and that $Y_{i0}|\lambda_i$ corresponds to the stationary distribution of $Y_{it}$ associated with its autoregressive law of motion. The implied marginal distribution for $Y_{i0}$ is used as $\pi(Y_{i0})$ in the Monte Carlo design. To obtain $\pi(\lambda_i|Y_{i0})$ we took $\pi_*(\lambda_i|Y_{i0})$ from the Gaussian model and replaced it with a mixture of normals described in Table (ref). For $\delta=0$ the mixture reduces to $\pi_*(\lambda_i|Y_{i0})$, whereas for large values of $\delta$ it becomes bimodal. This bimodality also translates into the distribution of $\hat{\lambda}|Y_{i0}$, which is depicted in Figure (ref) for $\delta=1/10$ (almost Gaussian) and $\delta=1$ (bimodal).
In this experiment we consider a parametric Tweedie correction (same as in Experiment 1, but now misspecified in view of the DGP) and two nonparametric Tweedie corrections. First, we compute the correction based on the simple Gaussian kernel in ((ref)). The bandwidth is chosen in accordance with the theory in Section (ref). We set $B_N = c/(\ln N)^{0.55}$, which would be consistent with a truncation of the form $C_N = c \sqrt{\ln N}$, and let $c \in \{1/2, 1, 2\}$.\footnote{The tuning matrices $V_{\hat{\lambda}}$ and $V_h$ are set equal to the sample variances of $\hat{\lambda}_i$ and $y_{i0}$, respectively.} Second, we use the adaptive estimator proposed by BotevGrotowskiKroese2010, henceforth BGK estimator, which is based on the solution of a diffusion partial differential equation. This estimator is associated with a plug-in bandwidth selection rule that requires no further tuning.\footnote{Our estimates are based on Algorithms 1 and 2 in BGK. We use the authors' MATLAB code to implement the density estimator.} Unless otherwise noted, the subsequent results are based on the BGK estimator.
Figure (ref) shows the “true” density $p(\hat{\lambda}_i|y_{i0},\theta)$ as well as Gaussian and nonparametric approximations. Under the Gaussian correlated random effects distribution we can directly calculate the conditional distribution of $\hat{\lambda}_i$ given $y_{i0}$. The nonparametric approximation is obtained by dividing an estimate of the joint density of $(\hat{\lambda}_i,y_{i0})$ by an estimate of the marginal density of $y_{i0}$ (this normalization is not required for the Tweedie correction). Each hairline in Figure (ref) corresponds to a density estimate from a different Monte Carlo run. For $\delta=1/10$ the Gaussian approximation is accurate and the variability of the estimates is much smaller than that of the kernel estimates. For $\delta=1$ the Gaussian density is unable to approximate the bimodal $p(\hat{\lambda}_i,y_{i0}|\theta)$, whereas the non-parametric approximation, at least for $y_{i0}=2.0$ captures the key features of the density of $\hat{\lambda}_i$.
For the prediction, the relevant object is the correction $(\sigma^2/T)\partial \ln p(\hat{\lambda}_i,y_{i0}|\theta) / \partial \hat{\lambda}_i$, which is depicted in Figure (ref). Under a Gaussian correlated random effects distribution, the Tweedie correction is linear in $\hat{\lambda}_i$ because the posterior mean is a linear combination of the prior mean and the maximum of the likelihood function. Thus, the corrections based on the Gaussian density estimate are linear regardless of $\delta$. For $\delta=1/10$ the correction under the “true” random effects distribution is nearly linear, and thus well approximated by the Gaussian correction. The nonparametric correction is fairly accurate for values of $\hat{\lambda}$ in the center of the conditional distribution $\hat{\lambda}_i|(y_{i0},\theta)$, but it becomes less accurate in the tails. For $\delta=1$, on the other hand, the kernel-based correction provides a much better approximation of the optimal correction than the Gaussian correction.
Table (ref) compares the performance of twelve predictors; half of them based on QMLE and the other half based on GMM. It is well-known that the GMM estimator of $\theta$ is consistent under the DGP described in Table (ref). We show in the Appendix that the QMLE estimator is also consistent for $\theta$ under this DGP, despite the fact that the correlated random effects distribution is misspecified. For each of the two $\theta$ estimators we construct posterior mean predictors using four different nonparametric Tweedie corrections as well as the Gaussian Tweedie correction. Moreover, we compute the plug-in predictor based on $\hat{\lambda}_i(\hat{\theta})$.
Among the nonparametric predictors, the one based on the BGK density estimator clearly dominates the ones derived from the simple kernel density estimator. If the random effects distribution is almost normal, i.e., $\delta=1/10$, setting $c=2$ is preferable to the other choices of $c$. For the bimodal random effects distribution, i.e., $\delta=1$, the best performance of the simple kernel estimator is attained for $c=1/2$. The predictors that rely on posterior mean approximations generally outperform the naive predictors based on $\hat{\lambda}_i(\hat{\theta})$. The benefits from shrinkage are most pronounced for the bottom and top groups. If the misspecification is small $(\delta=1/10)$, the parametric correction leads to more precise forecasts than the nonparametric correction because it is based on a more efficient density estimator. As the degree of misspecification increases, the nonparametric correction starts to perform better and for $\delta=1$ it clearly dominates the parametric competitor. This is consistent with the accuracy of the underlying density estimators shown in Figures (ref) and (ref).
In the third experiment, summarized in Table (ref), we consider a misspecification of the Gaussian likelihood function by replacing the Normal distribution in the DGP with two mixtures. We consider a scale mixture that generates excess kurtosis and a location mixture that generates skewness. The innovation distributions are normalized such that $\mathbb{E}[U_{it}]=0$ and $\mathbb{V}[U_{it}]=1$. For the heterogeneous intercepts $\lambda_i$ we adopt the Gaussian random effects specification of Experiment 1. In this experiment we compute the relative regret for five predictors:\footnote{The computation of the oracle predictor and the normalization of the regret by the posterior variance of $\lambda$ require a Gibbs sampler which is described in the Appendix.} the posterior mean predictor based on the non-parametric Tweedie correction and the plug-in predictor based on $\hat{\theta}_{QMLE}$ and $\hat{\theta}_{MLE}$, respectively. Note that both the QMLE and the GMM estimator of $\theta$ remain consistent under the likelihood misspecification. However, the (non-parametric) Tweedie correction no longer delivers a valid approximation of the posterior mean.
The results are summarized in Table (ref). The risk of the oracle predictors can be compared to that reported in Table (ref). The excess kurtosis of the scale mixture and the skewness of the location mixture slightly reduce the posterior variance of $\lambda$ compared to the standard normal benchmark in Experiment 1. Due to the misspecification of the likelihood function, the relative regret of the various predictors increases considerably, but the relative ranking is essentially unchanged. The posterior mean predictors based on the nonparametric Tweedie correction dominate all the other predictor, attaining a relative regrets of about 1 and 0.4, respectively. Compared to the plug-in and loss-function based predictors, the Tweedie correction still reduces the regret 40% to 50%. The predictor based on the pooled OLS estimation performs the worst among the five predictors in this experiment.
We will now use the previously-developed predictors to forecast pre-provision net revenues (PPNR) of bank holding companies (BHC). The stress tests that have become mandatory under the 2010 Dodd-Frank Act require banks to establish how PPNR varies in stressed macroeconomic and financial scenarios. A first step toward building and estimating models that provide trustworthy projections of PPNR and other bank-balance-sheet variables under hypothetical stress scenarios, is to develop models that generate reliable forecasts under the observed macroeconomic and financial conditions. Because of changes in the regulatory environment in the aftermath of the financial crisis as well as frequent mergers in the banking industry our large $N$ small $T$ panel-data-forecasting framework seems particularly attractive for stress-test applications.
We generate a collection of panel data sets in which pre-provision net revenue as a fraction of consolidated assets (the ratio is scaled by 400 to obtain annualized percentages) is the key dependent variable. The data sets are based on the FR Y-9C consolidated financial statements for bank holding companies for the years 2002 to 2014, which are available through the website of the Federal Reserve Bank of Chicago. Because the balance sheet data exhibit strong seasonal features, we time-aggregate the quarterly observations into annual observations and take the time period $t$ to be one year.
We construct rolling samples that consist of $T+2$ observations, where $T$ is the size of the estimation sample and varies between $T=3$ and $T=11$ years. The additional two observations in each rolling sample are used, respectively, to initialize the lag in the first period of the estimation sample and to compute the error of the one-step-ahead forecast. For instance, with data from 2002 to 2014 we can construct $M=9$ samples of size $T=3$ with forecast origins running from $\tau = 2005$ to $\tau = 2013$. Each rolling sample is indexed by the pair $(\tau,T)$. The cross-sectional dimension $N$ varies from sample to sample and ranges from approximately $=460$ to 725. Further details about the data as well as a description of our procedure to create balanced panels and eliminate outliers are provided in the Appendix.
In Section (ref) we use the basic dynamic panel data model to generate PPNR forecasts. In Section (ref) we extend the model to include covariates and compare forecasts under the actual realization of the covariates and stressed scenarios in which we set the covariantes to counterfactual levels.
We begin by evaluating forecasts from the basic dynamic panel model in ((ref)). The parametric Tweedie correction is based on $\lambda_i | (H_i,\theta) \sim N(\phi_0 + \phi_1 Y_{i0}, \underline{\omega}^2)$. The forecast evaluation criterion is the mean-squared error (MSE) computed across institutions and across time:
where $M$ is the number of rolling samples. Table (ref) summarizes the MSEs for different estimators and different sizes $T$ of the estimation samples. Recall that the unit of $\widehat{Y}_{i\tau}$ is annual revenue as fraction of total assets converted into annualized percentages.
For the short samples, i.e., $T=3$ and $T=5$, the QMLE-based predictors are more accurate than the GMM-based predictors. This discrepancy vanishes as the sample size is increased to $T=11$. The posterior mean predictors computed with the Tweedie correction are more accurate than the plug-in predictors. As expected, the MSE differential is largest in the small $T$ samples, because the unit-specific likelihood function contains fairly little information and the prior strongly influences the posterior. The parametric Tweedie correction delivers more accurate predictions than the non-parametric Tweedie correction, in particular for small $T$. In Figure (ref) we compare the Tweedie corrections for $T=5$ and $\tau=2012$. While the corrections are quite similar for values of the sufficient statistic $\hat{\lambda}_i(\rho) = \frac{1}{T} \sum_{t=1}^T (Y_{it}-\rho Y_{it-1})$ between -1% and 1%, the non-parametric correction behaves somewhat erratic outside of this interval which hurts the predictive performance.
Returning to the MSE results in Table (ref), the posterior mean predictor yields roughly the same MSE as pooled OLS. This suggests that {\em a posteriori} the data sets contain only weak evidence for heterogeneous intercepts. In this regard, the parametric specification is more efficient in shrinking the intercept estimates toward a common value. Finally, for all sample sizes except $T=11$, the posterior-mean predictor based on $\hat{\theta}_{QMLE}$ and the parametric Tweedie correction is more accurate than the loss-function-based predictor.
In Table (ref) we focus on the sample size $T=5$. In addition to averaging forecast errors across all $T=5$ samples, we also report results for specific forecast origins, namely choices of $\tau$ that correspond to the years 2007, the onset of the Great Recession, and 2012, which is during the recovery period. Moreover, we compute MSEs based on cross-sectional selection rules that depend on the level of PPNR at the forecast origin $\tau$. We focus on institutions with PPNR less than 0%, -1%, -2%, and -3%, respectively. Because the QMLE predictors dominate the GMM predictors and the parametric Tweedie correction was preferable to the nonparametric correction, we now restrict our attention to the posterior-mean predictor based on $\hat{\theta}_{QMLE}$ and the parametric Tweedie correction, the $\hat{\theta}_{QMLE}$ plug-in predictor, and predictors constructed from loss-function-based estimates and pooled OLS, respectively.
For the 2007 sample, the plug-in and the loss-function-based predictor are dominated by the other two predictors. The performance of the posterior-mean and the pooled-OLS predictor are essentially identical. For the 2012 sample, the posterior-mean predictor performs better than the plug-in predictor if we average across all institutions or if we condition on BCHs with PPNR of less than -3%. In the other cases the ranking is reversed. Across all rolling samples, the posterior mean predictor dominates. Across all institutions its performance is only slightly better than pooled OLS, but if we condition on BCHs with PPNR of less than -1%, -2%, or -3% then the accuracy relative to pooled OLS is more pronounced.
Table (ref) in the Appendix provides point estimates of the parameters of the basic dynamic panel model and the parametric correlated random effects distribution for $T=5$ and $\tau = 2007, \ldots, 2013$. Until 2010 the estimated variance of the correlated random effects distribution is essentially zero, which implies that $\lambda_i \approx \phi_0+\phi_1 Y_{i0}$. Because of a non-zero $\hat{\phi}_1$ the resulting predictor is not exactly pooled OLS but it is very similar as we have seen from the results in Table (ref). Starting in 2011, we obtain non-trivial estimates of $\hat{\underline{\omega}}^2$ which imply non-trival {\em a priori} dispersion of the intercepts (that is not due to the dispersion in initial conditions). Overall, the estimates $\hat{\underline{\omega}}^2$ imply a large degree of shrinkage. The positive estimate $\hat{\phi}_1$ generates positive correlation between $\lambda_i$ and $Y_{i0}$. The intercept of the correlated random effects distribution drops during the Great Recession\footnote{Recall that the $\tau=2010$ estimation sample comprises the observations for 2006-2010.}, which is consistent with the fact that bank revenues eroded during the financial crisis. The estimated common autoregressive coefficients range from 0.7 to 0.9.
To analyze the performance of the banking sector under stress scenarios it is necessary to add predictors to the dynamic panel data model that reflect macroeconomic and financial conditions. We consider three aggregate variables: the unemployment rate, the federal funds rate, and the spread between the federal funds rate and the 10-year treasury bill. Because these predictors are not bank-specific, the effect of the predictors on PPNR has to be identified from time-series variation, which is challenging given the short time-dimension of our panels. We consider two specifications: the first model only includes the unemployment rate as additional predictor and we focus on the $T=5$ data sets. The second model includes all three aggregate predictors and we estimated it based on the $T=11$ sample.
We generate forecasts using the actual values of the aggregate predictors (which we can evaluate based on the actual PPNR realizations for the forecast perior) and compare these forecasts to predictions under a stressed scenario, in which we use hypothetical values for the predictors. When analyzing stress scenarios, one is typically interested in the effect of stressed economic conditions on the current performance of the banking sector. For this reason, we are changing the timing convention slightly and include the time $t$ macroeconomic and financial variables into the vector $W_{it-1}$. We are implicitly assuming that there is no feedback from disaggregate BCH revenues to aggregate conditions. While this assumption is inconsistent with the notion that the performance of the banking sector affects macroeconomic outcomes, elements of the Comprehensive Capital Analysis and Review (CCAR) conducted by the Federal Reserve Board of Governors have this partial equilibrium flavor.
{\bf Results From a Model with Unemployment.} We use the unemployment rate (UNRATE) from the FRED database maintained by the Federal Reserve Bank of St. Louis and convert it to annual frequency by temporal averaging. We begin by computing MSEs, which are reported in Table (ref). This table has the same format as Table (ref): we consider MSEs for 2007, 2012, and averaged across all rolling samples. Moreover, we compute MSEs conditional on the level of PPNR at the forecast origin. A few observations stand out. First, the MSE for the posterior mean predictor is slightly reduced by including unemployment for the 2007 and 2012 samples, but across all of the rolling samples it slightly increases. Second, the gain of using the Tweedie correction, that is, the MSE differential between the plug-in predictor and the posterior mean predictor, becomes larger as we include unemployment. This is very intuitive: the more coefficients need to be estimated based on a given time-series dimension, the more important the shrinkage induced from the prior distribution. Third, the performance of the posterior-mean predictor and the pooled-OLS predictors remain very similar, meaning that the Tweedie correction shrinks toward pooled OLS.\footnote{This is supported by the estimates of $\hat{\underline{\omega}}_1^2$ and $\hat{\underline{\omega}}_2^2$ reported in the Online Appendix.}
We now impose stress by increasing the unemployment rate by 5%. This corresponds to the unemployment movement in the {\em severely adverse} macroeconomic scenario in the Federal Reserve's CCAR 2016. In Figure (ref) we are comparing one-year-ahead predictions for forecast origins $\tau=2007$ and $\tau=2012$ under the actual period $\tau+1$ unemployment rate and the stressed unemployment rate. Each circle in the graphs corresponds to a particular BHC. We indicate institutions with assets greater than 50 billion dollars\footnote{These are the BHCs that are subject to the CCAR requirements.} by red circles, while the other BHCs appear as blue circles. The large institutions have in general smaller revenues than the smaller BHCs. According to the plug-in predictor (the two right panels), the response to the unemployment shock is very heterogeneous. For about half of the intitutions a rise in unemployment leads to a drop in revenues, whereas for the other half higher unemployment is associated with larger revenues. However, we know from Table (ref) that forecasts from the plug-in predictor are fairly inaccurate. The stress-test implications of the posterior mean predictor are markedly different. Due to the strong shrinkage the effect is more homogeneous across institutions and appears to be slightly positive.
{\bf A Model with Unemployment, Federal Funds Rate, and Spread.} We now expand the list of covariates and in addition to the unemployment rate include the federal funds rate and the spread between the federal funds rate and the 10-year treasury bill. Both series are obtained from the FRED database (FEDFUNDS and DGS10). We convert the series into annual frequency by temporal averaging. Because we now have three regressors that do not vary across units (meaning all BHCs are operating within the same macroeconomic conditions, but may have hetereogeneous responses to these conditions), we focus on the data set with the largest time series dimension, namely $T=11$. MSEs are presented in Table (ref). The forecast origin is $\tau=2013$. As before, the posterior mean predictor with the Tweedie correction strongly dominates the plug-in predictor. Moreover, the posterior mean predictor is also slightly more accurate than the predictor based on pooled OLS.\footnote{While the estimates of the conditional variances of the $\lambda_{ij}$ coefficients are close to zero, the estimated conditional means of $\lambda_{ij}$ vary with $Y_{i0}$. This explains the difference between the posterior mean and the pooled-OLS predictor.} Unlike in the previous cases, the predictor constructed from the loss-function-based estimate of the model coefficients now performs slightly better than the posterior mean predictor.
Figure (ref) compares PPNR predictions under the actual macroeconomic conditions and a stressed macroeconomic scenario. The stressed scenario comprises an increase in the unemployment rate by 5% (as before) and an increase in nominal interest rates and spreads by 5%. This scenario could be interpreted as an aggressive monetary tightening that induced a sharp drop in macroeconomic activity. The plug-in predictor generates very heterogeneous responses to the macroeconomic stress scenario. Some banks benefit from the monetary tightening and others experience a substantial fall in revenues. The posterior mean predictor implies a much more homogeneous response of the banking sector under which there is a very small (relative to the cross-sectional dispersion) increase in predicted revenues.
{\bf Discussion.} We view this analysis as a first-step toward applying state-of-the-art panel data forecasting techniques to stress tests. First, it is important to ensure that the empirical model is able to accurately predict bank revenues and balance sheet characteristics under observed macroeconomic conditions. Our analysis suggests that there are substantial performance differences among various plausible estimators and predictors. Second, a key challenge is to cope with model complexity in view of the limited information in the sample. There is a strong temptation to over-parameterize models that are used for stress tests. We decided to time-aggregate the revenue data to smooth out irregular and non-Gaussian features of the accounting data at the quarterly frequency. This limits the ability to precisely measure the potentially heterogeneous effects of macroeconomic conditions on bank performance. Prior information is used to discipline the inference. In our empirical Bayes procedure, this prior information is essentially extracted from the cross-sectional variation in the data set. While we {\em a priori} allowed for heterogeneous responses, it turned out {\em a posteriori}, trading-off model complexity and fit, that the estimated coefficients exhibited very little heterogeneity. Third, our empirical results indicate that relative to the cross-sectional dispersion of PPNR, the effect of severely adverse scenarios on revenue point predictions are very small. We leave it future research to explore richer empirical models that focus on specific revenue and accounting components and consider a broader set of covariates. Finally, it would be desirable to allow for a feedback from the performance of the banking sector into the aggregate conditions.
The literature on panel data forecasting in settings in which the cross-sectional dimension is large and the time-series dimension is small is very sparse. Our paper contributes to this literature by developing an empirical Bayes predictor that uses the cross-sectional information in the panel to construct a prior distribution that can be used to form a posterior mean predictor for each cross-sectional unit. The shorter the time-series dimension, the more important this prior becomes for forecasting and the larger the gains from using the posterior mean predictor instead of a plug-in predictor. We consider a particular implementation of this idea for linear models with Gaussian innovations that is based on Tweedie's posterior mean formula. It can be implemented by estimating the cross-sectional distribution of sufficient statistics for the heterogeneous coefficients in the forecast model. We consider both parametric and nonparametric techniques to estimate this distribution. We provide a theorem that establishes a ratio-optimality property for the nonparametric estimator of the Tweedie correction. The nonparametric estimation works well in environments in which the cross-sectional distribution of heterogeneous coefficients is irregular. If it is well approximated by a Gaussian distribution, then a parametric implementation of the Tweedie correction is preferable. We illustrate in an application that our forecasting techniques may be useful to execute bank stress tests. Our paper focuses on one-step-ahead point forecasts. We leave extensions to multi-step forecasting and density forecasting for future work.
\setstretch{1}
\ifx\undefined\leavevmode\rule[.5ex]{3em}{.5pt}\ \fi \ifx\undefined\textsc \let\tmpsmall\tmpsmall\sc \fi
\setstretch{1.3}