Forecasting with a Panel Tobit Model

Forecasting with a Panel Tobit Model

JEL CLASSIFICATION: C11, C14, C23, C53, G21

KEY\ WORDS: Bayesian inference, density forecasts, loan charge-offs, panel data, set forecasts, Tobit model.

\thispagestyle{empty} \setcounter{page}{0}

\color{black} Introduction

This paper considers the problem of forecasting a large collection of short time series with censored observations. In the empirical application we forecast charge-off rates on loans for a panel of small banks. A charge-off occurs if a loan is deemed unlikely to be collected because the borrower has become delinquent. The prediction of charge-off rates is interesting to banks, regulators, and investors because they are losses on loan portfolios. If charge-off rates are large, the bank may be entering a period of distress and require additional capital. Due to mergers and acquisitions, changing business models, and changes in regulatory environments the time series dimension that is useful for forecasting is often short. The general methods developed in this paper are not tied to the charge-off rate application and can be used in any setting in which a researcher would like to analyze a panel of censored data with a large cross-sectional and a short time-series dimension.

In a panel data setting, cross-sectional heterogeneity in the data is modeled through unit-specific parameters. The more precisely they are estimated, the more accurate the forecasts are. The challenge in forecasting panels with a short time dimension is that the data set does not contain a lot of information about the heterogeneous parameters. A natural way of adding information to the estimation of these parameters is the use of prior distributions. The key insight in panel data applications is that one can extract information from the cross section and equate the prior distribution with the cross-sectional distribution of unit-specific coefficients. An empirical Bayes implementation of this idea creates a point estimate of the cross-sectional distribution of the heterogeneous coefficients and then conditions the subsequent posterior calculations on the estimated prior distribution. The classic James-Stein estimator for a vector of means can be interpreted as an empirical Bayes estimator.\footnote{Empirical Bayes methods have a long history in the statistics literature going back to Robbins1955; see Robert1994 for a textbook treatment.}

Rather than pursuing an empirical Bayes approach, we conduct a full Bayesian analysis by specifying a hyperprior for the distribution of heterogeneous coefficients and constructing a joint posterior for the coefficients of this hyperprior as well as the actual unit-specific coefficients. This approach can in principle handle quite general nonlinearities and generate predictions under a wide variety of loss functions. It is preferable for interval and density forecasts, because it captures all sources of uncertainty.

The contributions of our paper are threefold. First, we extend the full Bayesian estimation and prediction with a linear panel data model in liu2018density to a dynamic panel Tobit model with heteroskedastic innovations and correlated random effects. {\color{black}We hereby build on work on the Bayesian estimation of static, dynamic, and panel Tobit models by Chib1992, Wei1999, baranchuk2008assessing, and LiZheng2008.}

Second, we construct interval forecasts that target average posterior coverage probability across all units in our panel instead of pointwise coverage probability for each unit. We show that it is optimal to generate these forecasts as highest posterior density sets that use the same threshold for each unit instead of unit-specific thresholds. Because the predictive distributions associated with the Tobit models are mixtures of discrete and continuous distributions, “interval” forecasts may take the form of the union of one or more intervals and the value zero, and thus we refer to them as set forecasts subsequently. {\color{black}We prove that the empirical coverage frequency converges to the average nominal coverage frequency of the sets as the cross-sectional dimension of the panel tends to infinity.} This result is connected to similar findings in the literature on nonparametric function estimation and dates back to Wahba1983 and Nychka1988. {\color{black}The underlying insights also have been recently used in concurrent research by armstrong2020robust to construct empirical Bayes confidence intervals for vectors of means that are valid for multiple priors.} In the Monte Carlo study and the empirical application the proposed Bayesian set forecasts have good finite sample frequentist coverage properties in the cross-section.

Third, we present a novel application in which we forecast bank-level loan charge-off rates. Our empirical analysis is based on more than 100 short panel data sets with a time dimension of $T=10$. These panel data sets include predominantly credit card (CC) and residential real estate (RRE) loans and cover various (overlapping) time periods. We also include local economic conditions as bank-specific regressors with homogeneous coefficients. For each data set, we document the density forecasting performance of several model specifications. We find that allowing for heteroskedasticity is important for good density and set forecasting performance. Overall, a specification with flexibly modeled correlated random effects and heteroskedasticity performs well in terms of density forecasting and is used in the subsequent analysis. In addition, we generate maps that compare the spatial distribution of predicted loan losses during and after the Great Recession and plot cross-sectional distribution of set forecasts. We document how set forecasts change as we move from targeting pointwise coverage probability to targeting average coverage probability. The latter approach smooths out differences among the lengths of the set forecasts and overall improves the forecasts with respect to both coverage probability and average length.

The heterogeneous intercepts in our model can be interpreted as estimates of the quality of the banks' loan portfolios. Loan quality is potentially determined by many factors: the risk taking behavior of the bank, the potential customer base, and its ability to efficiently screen borrowers. In regressing heterogeneous coefficient estimates on bank characteristics we find that bank size as measured in total assets is positively related to inverse quality of the loan portfolio. A favorable interpretation of this finding is that larger banks are able to take higher risks on loans because they are better diversified or have a higher tolerance for risk. However, overall bank characteristics explain only a very small fraction of the estimated heterogeneity.

Because the Tobit model is nonlinear, the effect of a change in local economic conditions that enter the model with homogeneous coefficients depends on the heterogeneous intercept and is thereby bank specific. We are able to compute a posterior distribution of the “treatment” effect for each bank and decompose it into an extensive-margin effect (a bank switches from no charge-offs to positive charge-offs during an economic downturn) and an intensive-margin effect (a bank increases its positive charge-offs during a downturn). We find that the variation in charge-off rates generated by local economic conditions is very small compared to the variation due to the heterogeneous intercept estimates.

Our paper relates to several branches of the literature. {\color{black}We build on the Bayesian literature on the estimation of censored regression models.\footnote{A general survey of the literature on Bayesian estimation of univariate and multivariate censored regression models can be found, for instance, in the handbook chapter by LiTobias2011.} The approach of using data augmentation for limited-dependent variable models that impute the latent uncensored variables dates back to Chib1992 and albert1993bayesian. To sample the latent observations we rely on an algorithm tailored toward dynamic Tobit models by Wei1999. Sampling from Truncated Normal distributions is implemented with a recent algorithm of Botev2017. Bayesian panel Tobit models have been estimated by baranchuk2008assessing and LiZheng2008. Our flexible benchmark model is most closely related to the semiparametric model of LiZheng2008 which we generalize by introducing heteroskedasticity through a latent unit-specific error variance and allowing for a more flexible form of correlated random effects. As mentioned previously, the former is very important for the density and set forecast performance.\footnote{baranchuk2008assessing report some results on point forecasts of the probability of zeros versus non-zeros, whereas we focus on set and density forecasts.}}

We model the unknown distribution of the heterogeneous coefficients (intercepts and innovation variances) as Dirichlet process mixtures (DPM) of Normals. Even though we do not emphasize the nonparametric aspect of this modeling approach (due to a truncation, our mixtures are strictly speaking finite and in that sense parametric), our paper is related to the literature on nonparametric density modeling using DPM.\footnote{KeaneStavrunova introduce a smooth mixture of Tobits to model a cross-section of healthcare expenditures. Our model is related, but different in that we are using a DPM to average across different intercept values and innovation variances.} Examples of econometrics papers that use DPMs in the panel data context are Hirano2002, burda2013panel, 10.2307/j.ctt5hhrfp, and jensen2015mutual. The implementation of our Gibbs sampler relies on IshwaranJames2001,IshwaranJames2002.

{\color{black}As an alternative to a full Bayesian analysis, recent papers by GuKoenkerJAE2016,GuKoenker2014 and LiuMoonSchorfheide2015 have pursued an empirical Bayes strategy to generate predictions based on linear panel data models with heterogeneous coefficients. Forecasts from empirical Bayes and full Bayesian estimation approaches have desirable optimality properties as the cross-sectional dimension of the data set gets large. LiuMoonSchorfheide2015 generalize optimality results for the estimation of a vector of means in BrownGreenshtein2009 to a linear dynamic panel data forecasting setting. liu2018density shows that the predictive density obtained from the full Bayesian analysis of a linear panel data model converges to the predictive density derived from the true cross-sectional distribution of the heterogeneous coefficients as the cross-section gets large.}

There also exists a literature on estimating the determinants of loan losses. This literature often uses nonperforming loans (loans that have not been serviced for more than 90 days) and tends to ignore the censoring which is reasonable if one uses an average across banks but can be problematic if one uses bank-level data. The two papers most closely related to our work are Ghosh2015,Ghosh2017. We base our choice of bank-characteristic regressors on these papers.

The remainder of our paper is organized as follows. Section (ref) presents the specification of our dynamic panel Tobit model, a characterization of the posterior predictive distribution for future observations, and discusses the construction and evaluation of density and set forecasts. Section (ref) provides details on how we model the correlated random effects distribution and heteroskedasticity. It also presents the prior distributions for the parametric and flexible components of the model, and outlines a posterior sampler. We conduct a Monte Carlo experiment in Section (ref) to examine the performance of the proposed techniques in a controlled environment. The empirical application in which we forecast charge-off rates on various types of loans for a panel of banks is presented in Section (ref). Finally, Section (ref) concludes. Detailed derivations and proofs, a description of the data sets, and additional simulation and empirical results are relegated to the Online Appendix.

Model Specification and Forecast Evaluation

Throughout this paper we consider the following dynamic panel Tobit model with heterogeneous intercepts and innovation variances:

where $i=1,\ldots,N$, $t=1,\ldots,T$, and $\mathbb{I}\{ y \ge a \}$ is the indicator function that is equal to one if $y \ge a$ and equal to zero otherwise. {\color{black} Throughout the paper, we abbreviate sequences of the form $(a_1,\ldots,a_n)$ by $a_{1:n}$. For instance, $Y^*_{1:N,0:t-1} = \big\{ (y_{10}^*,\ldots,y_{N0}^*),\ldots,(y_{1t-1}^*,\ldots, y_{Nt-1}^*) \big\}$, and $\lambda_{1:N} = (\lambda_1,\ldots,\lambda_N)$. The $n_x \times 1$ vector $x_{it}$ comprises a set of sequentially exogenous regressors. $\xi$ is a vector of hyperparameters defined in ((ref)) below that does not affect the conditional distribution of $y_{it}^*$. It is assumed that conditional on the parameters and the regressors $x_{it-1}$, the observations $y_{it}$ are cross-sectionally independent.} {\color{black}The distributional assumption in ((ref)) implies that we can write

which we will use subsequently to simplify formulas.} Our specification uses the lagged latent variable $y_{it-1}^*$ on the right-hand side because it is more plausible for our empirical application. The Bayesian computations described in Section (ref) below can be easily adapted to the alternative model, in which the lagged censored variable $y_{it-1}$ appears on the right-hand side.

We model the heterogeneous parameters as correlated random effects (CRE) with density

assuming cross-sectional independence of the heterogeneous coefficients.\footnote{We consider period $t=-1$ for $x$ in the conditioning set because of the timing assumption that charge-off rates can only respond with a one-period lag to changes in local economic conditions so as to accommodate possible sequentially exogenous regressors. See Section (ref) for more details.} Here $\xi$ is a hyperparameter vector that indexes a family of CRE distributions. For instance, the candidate distribution of $( \lambda_i,y_{i0}^*,\ln\sigma_i^2 )$ could be jointly Normal with a mean that is a linear function of $x_{i,-1}$. In this case $\xi$ would include the parameters of the conditional mean function and the non-redundant parameters of the covariance matrix. To achieve a flexible representation of the distribution of $( \lambda_i,y_{i0}^*,\sigma_i^2 )$ we consider a family of mixtures of Normal distributions in Section (ref). {\color{black}We define the homogeneous parameter $\theta = [\rho,\beta']'$ and complete the model with the specification of a prior distribution for $\big( \theta,\xi \big)$.}

{\color{black}Our model is closely related to the panel Tobit models of baranchuk2008assessing and LiZheng2008, henceforth BC and LZ, respectively. However, the modeling approaches differ with respect to the treatment of coefficient heterogeneity and heteroskedasticity.\footnote{\color{black} As in the panel Probit model of chib2006inference, one could allow for additional lags of $y_{it}^*$.} As in LZ, we restrict regression coefficient heterogeneity to the intercept. We also follow LZ in modeling the CRE distribution in ((ref)) nonparametrically, albeit the details are slightly different. Because the regressors $x_{it}$ in our application are not assumed to be strictly exogenous, we condition the distribution of $(\lambda_i,y_{i0}^*)$ only on the initial values $x_{i,-1}$ and not on other $x_{it}$s. The most important difference between our specification and that of LZ is that we allow for heterogeneous innovation variances $\sigma^2_i$, whereas LZ set $\sigma^2_i = \sigma^2$ for all $i$. As documented in Section (ref), $\sigma^2_i$ heterogeneity is very important for the construction of accurate set and density forecasts in our empirical application.

BC restrict the distribution of the heterogeneous coefficients to be Normal, but they do allow regression coefficients other than the intercept to be heterogeneous.\footnote{\color{black} Our framework can be easily extended to accommodate heterogeneous slope coefficients (see LiuMoonSchorfheide2015 and liu2018density).} Rather than linking the heterogeneity to the regressors $x_{it}$, they let the mean of the distribution depend on additional unit-specific covariates. Instead of embedding additional covariates (such as bank characteristics) {\em ex ante} into ((ref)), we run {\em ex post} regressions of estimates of the ratio $\widehat{\lambda_i/\sigma_i}$ on additional unit-specific covariates to explore potential relationships. The reasons for conducting an {\em ex post} analysis in our application are threefold: (i) it is not clear {\em ex ante} which bank characteristics are relevant, (ii) the relationship between bank characteristics and cross-sectional heterogeneity could be nonlinear, and (iii) bank characteristics may only explain a small fraction of the cross-sectional heterogeneity.

BC's interaction between regressors and the Normal CRE distribution generates heteroskedasticity in what could be interpreted as composite error term that consists of a homoskedastic innovation in the regression equation for $y_{it}^*$ and the randomness in the heterogeneous coefficients scaled by the regressors. In our model specification, the heteroskedasticity is unrelated to the regressors $x_{it}$ because we are treating the $\sigma^2_i$ as random effects. A relationship to the regressors could be generated through a CRE specification for $\sigma^2_i$, but we did not pursue this extension because in our application the regressors, local unemployment and house price growth, cannot explain the dispersion in $\sigma^2_i$. }

In the remainder of this section, we {\color{black}discuss our assumptions about the simultaneous determination of outcomes $y_{it}$ and regressors $x_{it}$ in Section (ref),} the derivation of the posterior predictive density in Section (ref), the density forecast evaluation criteria in Section (ref), and the construction and evaluation of set forecasts in Section (ref).

Simultaneity and Timing Assumptions

{\color{black}In our application $y_{it}$ corresponds to bank-level loan charge-off rates and the regressors $x_{it}$ measure local economic conditions, such as unemployment and house prices, in the state in which the bank operates.\footnote{We consider a sample of small banks that conduct most of their business locally.} In this context it is plausible to assume that there is feedback from the bank charge-offs, which affect profitability and overall health of the banking sector, to the local economic conditions.

The key assumption that we are making throughout the paper is that charge-off rates are only affected by lagged economic conditions and not by contemporaneous economic conditions. For concreteness, suppose that $x_{it}$ corresponds to economic conditions in the state in which bank $i$ operates. We assume that the state-level conditions in period $t=0,\cdots,T$ are described by the conditional density

Thus, we allow current charge-offs to affect current state-level conditions. However, we assume that $X_{1:N,t}$ does not separately depend on the latent variables $Y_{1:N,0:t}^*$ and the heterogeneous coefficients $(\lambda_i,\sigma_i^2)$. In our application only actual charge-off rates are assumed to matter for economic outcomes. $\theta_x$ is a vector of parameters determining the law of motion for the state-level conditions.

Timing restrictions such as the one above have traditionally been widely used in the macroeconometric literature on structural vector autoregressions; see, for instance, the survey by Ramey2016. Here we are assuming that a deterioration of macroeconomic conditions affects banks' decisions to write off loans with a one period delay, where the length of a period is a quarter in our application.\footnote{\color{black} Relaxing this assumption is beyond the scope of this paper.} Combining ((ref)), ((ref)), and ((ref)), we can write

In slight abuse of notation $p(y_{i0}|y_{i0}^*)$ represents the censoring. The distribution of $y_{it} | y_{it}^*$ is a unit point mass that is located at 0 if $y_{it}^* \le 0$ or at $y_{it}^*$ if $y_{it}^* > 0$. Because the system is triangular, the panel Tobit component in ((ref)) can be estimated independently of ((ref)) and without the use of instrumental variables.}

Posterior Predictive Densities

Our goal is to generate forecasts of $Y_{1:N,T+h}$ conditional on the observations $(Y_{1:N,0:T},X_{1:N,-1:T})$. In the empirical analysis in Section (ref) we focus on $h=1$-step-ahead forecasts which require the predictor $x_{iT}$, which is known at the forecast origin $t=T$. The extension to multi-step forecasts is discussed in Section (ref). {\color{black}Because in a Bayesian framework uncertainty with respect to parameters, latent variables, and future shocks is treated identically through the use of random variables, it is conceptually straightforward to construct a predictive distribution of $Y_{1:N,T+1}$ conditional on $(Y_{1:N,0:T},X_{1:N,-1:T})$ by integrating out all sources of uncertainty. The general approach is summarized, for instance, in GewekeWhiteman2006. We subsequently describe the integration steps required for our panel Tobit model.}

According to ((ref)) the distribution of $(Y_{1:N,0},Y_{1:N,0}^*)$ conditional on $X_{1:N,-1}$ does not depend on $\theta_x$. Using the factorization in ((ref)), the CRE density ((ref)), and the prior $p(\theta,\xi) = p(\theta)p(\xi)$, we can write the posterior distribution of the parameters and time-$T$ latent variables as

where $\propto$ denotes proportionality. The posterior predictive distribution for units $i=1,\ldots,N$ is given by

Draws from $p(Y_{1:N,T+1}|Y_{1:N,0:T},X_{1:N,-1:T})$ can be generated by sampling $(Y_{1:N,T}^*, \lambda_{1:N},\sigma^2_{1:N},\theta,\xi)$ from the posterior ((ref)) and then evaluating the autoregressive law of motion for $y_{it}^*$ in ((ref)) for $t=T+1$.

To simplify the notation, we drop $X_{1:N,-1:T}$ from the conditioning set in the remainder of this section. Moreover, we denote the forecast horizon by $h$ again with the understanding that the discussion of multi-step forecasts is deferred to Section (ref). We denote expectations and probabilities under the posterior predictive distribution by $\mathbb{E}_{Y_{1:N,0:T}}^{y_{iT+h}}[\cdot]$ and $\mathbb{P}_{Y_{1:N,0:T}}^{y_{iT+h}}\{\cdot\}$, respectively. More generally, we use subscripts to indicate the conditioning set and superscripts to denote the random variables over which the operators integrate. The predictive distribution is a mixture of a point mass at zero and a continuous distribution for realizations of $y_{iT+h}$ that are greater than zero:

Here $\delta_0(y)$ is the Dirac function with the property $\delta_0(y)=0$ for $y \not=0$ and $\int \delta_0(y)dy = 1$. The density $p_c(y_{iT+h}|Y_{1:N,0:T})$ represents the continuous part of the predictive distribution.

Evaluating Density Forecasts

To compare the density forecast performance of various model specifications $M$ we report the average log predictive scores

and continuous ranked probability scores (CRPSs). The CRPS measures the $L_2$ distance between the cumulative distribution function $F_{Y_{1:N,0:T}}^{y_{iT+h}}(y|M)$ associated with $p(y_{iT+1}|Y_{1:N,0:T})$ and a “perfect” density forecasts which assigns probability one to the realized $y_{iT+h}$. Then,

Both LPS and CRPS are proper scoring rules, meaning that it is optimal for the forecaster to truthfully reveal her predictive density GneitingRaftery.

Constructing and Evaluating Set Forecasts

We construct set forecasts from the posterior predictive distribution $p(y_{iT+h}|Y_{1:N,0:T})$ in ((ref)) of the form:

with the understanding that (i) $C_i = \{0\}$ if $K_i=0$, (ii) $a_{i1}$ may be equal to zero, and (iii) \[ a_{i1} < b_{i1} < a_{i2} < b_{i2} < \ldots < a_{iK_i} < b_{iK_i}. \] The $\{0\}$ value arises from the discrete portion of the predictive density, whereas the interval components are obtained from the continuous portion of the predictive density; see the decomposition in ((ref)).\footnote{Because in our model the support of the posterior predictive distribution of $y_{iT+h}^*$ includes $y < 0$, the probability of censoring is strictly positive and the set that includes $\{0\}$ is strictly shorter than the one without zero.} The disjoint interval segments may arise if the continuous part of the predictive density is multimodal. If we target an average coverage probability in the cross section, then for some units $i$ we might obtain the empty set, i.e., $C_{iT+h|T}(Y_{1:N,0:T}) = \emptyset$.

{\bf Constructing Set Forecasts.} To generate the set forecasts, we adopt a Bayesian approach and require that the probability of $\{y_{iT+h} \in C_{iT+h|T}(Y_{1:N,0:T})\}$ conditional on having observed $Y_{1:N,0:T}$ reaches a pre-specified level. Given that the estimation of the Tobit model is executed with Bayesian techniques, the use of posterior predictive credible sets is natural. We distinguish between forecasts that are constructed to satisfy the coverage probability constraint pointwise, that is,

and sets that are constructed to satisfy the constraint on average:

The latter approach allows the sets $C_{iT+h|T}(Y_{1:N,0:T})$ for some units $i$ to be “shortened” in the sense that their posterior credible level drops below $1-\alpha$, whereas sets for other units are “lengthened.”

It is well known that the shortest credible sets take the form of highest posterior density sets. Suppose that we require to satisfy the coverage constraint for each $i$ individually. If $\mathbb{P}_{Y_{1:N,0:T}}^{y_{iT+h}} \{ y_{iT+h} = 0\} \ge 1-\alpha$, then $C_{iT+h|T}(Y_{1:N,0:T}) = \{0\}$. Otherwise, the set takes the form

where the threshold $\kappa_i$ is chosen such that \[ \int_{y_{iT+h} \in C} p_c(y_{iT+h}|Y_{1:N,0:T})\mathbb{I}\{ y_{iT+h} \ge 0\}dy_{iT+h} = 1 - \alpha - \mathbb{P}_{Y_{1:N,0:T}}^{y_{iT+h}} \{ y_{iT+h} = 0\}. \] Because $p_c(y|\cdot)$ is a continuous density, the HPD set can be represented as a collection of disjoint intervals as in ((ref)).

If the objective is to minimize average length across $i$ conditional on the constraint on coverage probability holding only on average, then the unit-specific thresholds $\kappa_i$ in ((ref)) are replaced by a common threshold $\kappa$ that applies to all units $i$. One can establish the optimality of the common threshold as follows. Suppose that one lowers the threshold for unit $i$ ($\kappa_i < \kappa$) and raises it for unit $j$ ($\kappa_j > \kappa)$. This lengthens the set for unit $i$ by $\delta_i > 0$ and shortens the set for unit $j$ by $\delta_j < 0$. The increase in coverage probability for unit $i$, $\Delta \pi_i > 0$, is less than $\delta_i \kappa$, whereas the decrease in coverage probability for unit $j$, $\Delta \pi_j < 0$, is less than $\delta_j \kappa$. Because we are holding the overall coverage probability constant, we obtain: \[ \delta_i \kappa > \Delta \pi_i = - \Delta \pi_j > - \delta_j \kappa. \] Thus, $\delta_i > - \delta_j$, which means that the overall average length increases and the uniform threshold of $\kappa$ dominates.

{\bf Evaluation of Set Forecasts.} The assessment of the set forecasts in our simulation study and the empirical application is based on the cross-sectional coverage frequency

and the average length of the sets $C_{iT+h|T}(Y_{1:N,0:T})$

Rather than trading off average length against deviations of average coverage frequency from the nominal coverage probability in a single criterion, we simply report both.\footnote{For various approaches to rank set forecasts see AskanaziEtAl2018.}

The relationship between the nominal credible level of the set forecasts and the empirical coverage frequency is delicate. In {\color{black}Theorem (ref)} below we provide high-level regularity conditions under which

in $\mathbb{P}^{Y_{1:N,0:T},Y_{1:N,T+h}}$ probability as $N \longrightarrow \infty$. {\color{black}Underlying this results is the well-known insight -- see, for instance, the textbook by Robert1994 -- that, for a generic parameter $\varsigma$ and data set $Y$, the following relationship between credible sets and confidence sets holds: \[ \mathbb{P}^{\varsigma,Y} \{ \varsigma \in C(Y)\} = \int_Y \mathbb{P}_Y^\varsigma \{ \varsigma \in C(Y) \} d\mathbb{P}^Y= \int_\varsigma \mathbb{P}_\varsigma^Y \{ \varsigma \in C(Y)\} d \mathbb{P}^\varsigma. \] Thus, $1-\alpha$ Bayesian credible sets have on average $1-\alpha$ frequentist coverage probability, but not pointwise for each $\varsigma$. In our framework the cross-sectional averaging across $i$ approximates the integration under the prior distribution.} The basic insight has previously been used in the literature on nonparametric function estimation, dating back to Wahba1983 and Nychka1988, to obtain results that link average coverage probabilities to Bayesian credible levels. {\color{black}More recently, armstrong2020robust constructed empirical Bayes confidence intervals for vectors of means that are valid for multiple priors.}

{\color{black} Let $\vartheta = (\theta,\xi)$.} To state the theorem we define the following probability associated with the interval $[a_{ik,N},\, b_{ik,N}]$ conditional on $(Y_{i,0:T},\vartheta)$:

A proof of this theorem is provided in the Online Appendix. {\color{black}Assumption (i) states that the future observations are generated from the “true” predictive density $p(Y_{1:N,T+h}|Y_{1:N,0:T})$.} In Assumption (ii) we require the posterior distribution of $\vartheta$ to concentrate. Throughout the paper, we represent the CRE distribution through finite-dimensional mixtures; see Section (ref). Thus, $\vartheta$ is finite-dimensional and the concentration results can be obtained from the literature on the consistency and asymptotic Normality of posterior distributions; see Hartigan1983, vanderVaart1998, GoshRamamoorthi2003, or GhosalVaart2017 for textbook treatments. The only difference to many of the results stated in the literature is that we assume that the convergence in probability to occur under the marginal distribution of $Y_{1:N,0:T}$ rather than its distribution conditional on a “true” parameter which imposes some restrictions on the prior for $\vartheta$. Assumptions (iii) and (iv) require the probabilities $F_{ik,N}$ to be smooth functions of $\vartheta$. In our model the probabilities are computed from finite-dimensional mixtures of Normal distributions, which are smooth functions of the underlying parameters. However, the Lipschitz constants are generally sample dependent and one needs to require that their average across $i$ and $k$ is stochastically bounded. In the Online Appendix we verify the conditions for a simple model without censoring.

Correlated Random Effects, Priors, and Posteriors

We provide a characterization of the CRE distribution $p(\lambda_i,y^*_{i0},\sigma^2_i|x_{i,-1},\xi)$ and a specification of the prior distribution for $(\theta,\xi)$ in Section (ref). Section (ref) contains a description of the posterior sampler, and Section (ref) outlines multi-step forecasting approaches.

(Correlated) Random Effects and Prior Distributions

We now describe the prior distribution for $\theta$, the parametrization of the distribution of $(\lambda_i,y_{i0}^*)$, and the prior distribution for the hyperparameter vector $\xi$. We begin with a homoskedastic random effects (RE) setup in which $\lambda_i$ and $y_{i0}^*$ are independent of each other and of $x_{i,-1}$. We then introduce heteroskedasticity and finally extend the model specification to CRE. The prior distribution involves a small number of tuning constants, denoted by $\tau$, that allow the researcher to scale the prior in various dimensions.

The subsequent exposition involves various parametric probability distributions in addition to the Normal distribution that appeared in ((ref)). We use $B(a,b)$, $G(a,b)$, and $IG(a,b)$ to denote the Beta, Gamma, and Inverse Gamma distributions, respectively. The pair $(\theta,\sigma^2)$ has a Normal-Inverse-Gamma distribution $NIG(m,v,a,b)$ if $\sigma^2 \sim IG(a,b)$ and $\theta|\sigma^2 \sim N(m,\sigma^2 v)$. Finally, the pair $(\Phi,\Sigma)$ has a matricvariate Normal-Inverse-Wishart distribution $MNIW(M,V,\nu,S)$ if $\Sigma \sim IW(\nu,S)$ has an inverse Wishart distribution and $\mbox{vec}(\Phi) |\Sigma \sim N(\mbox{vec}(M), \Sigma \otimes V )$.

{\bf Prior for $\theta$.} We standardize the regressors $x_{it}$ to have zero mean and unit variance and use the following Normal prior for the regression coefficients $\theta$:

where $\tau_\theta$ is a tuning constant that controls the prior variance.

{\bf Flexible RE with homoskedasticity.} Under RE, the distribution of $\lambda_i$ and $y_{i0}^*$ does not depend on $x_{i,-1}$. Moreover, we assume that $\lambda_i$ and $y_{i0}^*$ are independent. Thus, \[ p( \lambda_i,y^*_{i0} | x_{i,-1},\xi ) = p(\lambda_i|\xi)p(y^*_{i0}|\xi). \] We consider a mixture representation for $p(\lambda_i|\xi)$ while assuming that the initial values $y_{i0}^*$ are normally distributed:

The maximum number of mixture components $K$ is assumed to be pre-specified.\footnote{We use $K=20$ in the simulation exercise and the empirical analysis. This leads to the following uniform bound on the approximation error (see Theorem 2 of IshwaranJames2001): $ \Vert f^{\lambda,K}-f^{\lambda}\Vert \sim4N\exp[-(K-1)/\alpha]\le2.24\times10^{-5}, $ at the prior mean of $\alpha$ ($=1$) and a cross-sectional sample size $N=1000$.}

A prior over the RE distributions is induced through a prior $p(\xi)$ for the hyperparameter vector \[ \xi = \big[ \phi_{\lambda,1},\Sigma_{\lambda,1},\pi_{\lambda,1}, \ldots, \phi_{\lambda,K},\Sigma_{\lambda,K},\pi_{\lambda,K}, \phi_y,\Sigma_y \big]'. \] During the Bayesian inference stage, the prior is updated in view of the data and we obtain a posterior distribution for $\xi$ and hence a posterior distribution for the RE distribution. The priors for the coefficients of the Normal distributions are

We parameterized the IG distribution such that the variances $\Sigma_{\lambda,k}$ and $\Sigma_{y}$ have a prior distribution with mean $\tau_\sigma$ and variance $\tau_\sigma^2$ (omitting the superscripts).\footnote{Under our parametrization of the $X \sim IG(a,b)$ distribution, $\mathbb{E}[X] = b/(a-1)$ for $a>1$, and $\mathbb{V}[X] = (\mathbb{E}[X])^2/(a-2)$ for $a>2$.} Conditional on $\Sigma$, the mean parameter $\phi$ has a $N(0,\tau_\phi \Sigma)$ distribution (omitting the subscripts). The marginal distribution of $y_{i0}^*$ implied by ((ref)) and ((ref)) is a Student-$t$ distribution, whereas the distribution of $\lambda_i$ is a mixture of Student-$t$ distributions. The tuning constants can be used to control the spread of the means of the mixture components as well as the magnitude and variation of the variances of the mixture components.

The prior for the probabilities $\pi_{\lambda,1:K}$ is generated by a mixture of truncated stick breaking processes $TSB(1,\alpha_\lambda,K)$ of the form

Note that the $B(1,\alpha_\lambda)$ prior has a density $p(\zeta_k) \propto (1-\zeta_k)^{(\alpha_\lambda-1)}$. If $\alpha_\lambda$ is close to zero, then a lot of the mass of the distribution is concentrated near $\zeta_k=1$. This means that the first mixture component has a probability that is close to one, whereas the remaining mixture components have very small probabilities. If $\alpha_\lambda$ is close to two, then most of the mass of the distribution of $\zeta_k$ is concentrated on values of $\zeta_k$ that are close to zero. In turn, a larger number of mixture components receive non-trivial probabilities. The $G(2,2)$ distribution is recommended by Ishwaran and James (2002). It has a mean of one and draws fall with 95% probability into the interval $[0.12,\, 2.75]$ which means that the prior covers both mixtures dominated by few components and mixtures with many non-trivial components.

In the homoskedastic specification, we use the conjugate prior for $\sigma^2$ that arises in the context of a linear regression model:

The IG distribution is parameterized in a similar way as the IG distributions in ((ref)). $V^*=\frac{1}{N} \sum_{i=1}^N \widehat{\mathbb{V}}_{i}(y_{it})$ is the cross-sectional average of the time-series variances of $y_{it}$ and the tuning constant $\tau_v$ provides additional flexibility to scale the prior for $\sigma^2$.

{\bf Heteroskedasticity.} To generate heteroskedasticity one could simply replace ((ref)) by $\sigma^2_i \sim IG\big(3,2 \tau_v V^* \big)$. However, to make the distribution a bit more flexible, we augment the hyperparameter vector $\xi$ and also represent the distribution of $\ln \sigma_i^2$ as a mixture of Normals:\footnote{\color{black}In an earlier version of the paper we used a mixture of IG distributions. We switched to a mixture of Normals for $\ln \sigma^2_i$ for a more symmetric treatment of $\lambda_i$ and $\sigma^2_i$. Alternatively, chib2002semiparametric used Dirichlet process prior with an IG base measure to generate scale mixtures of Normals.}

A straightforward change-of-variables yields the distribution $p(\sigma^2_i|\xi)$. As for the RE distribution, the coefficients $\psi_k$ and $\omega_k$ have $NIG$ priors:

The parametrization is chosen so that the implied prior mean $\mathbb{E}[\sigma_i^2]$ and prior variance $\mathbb{V}[\sigma_i^2]$ for each mixture component $k$ matches the one implied by the prior used in the homoskedastic version of the Tobit model; see ((ref)).\footnote{The marginal IG distribution implies $\mathbb{E}[\omega_k^2] =\ln 2$. Conditional on $\omega_k^2 = \ln 2$, the transformed parameter $\exp(\psi_k)$ has a Lognormal distribution with mean $ \tau_v V_*$ and variance $( \tau_v V_*)^2$.} Moreover, we verified by simulation that the marginal density of $\sigma_i^2$ under this prior is very similar to the $IG(3,(3-1) \tau_v V^*)$ distribution used for the homoskedastic specification. It does, however, have fatter tails as it is a mixture of log $t$ distributions.

{\bf Flexible CRE with heteroskedasticity.} We extend the RE specification in two directions: first, we allow for correlation of $\lambda_i$ and $y_{i0}^*$ with $x_{i,-1}$. Second, we allow $\lambda_i$ and $y_{i0}^*$ to be correlated with each other conditional on $x_{i,-1}$. The CRE distribution is given by the following location and scale mixture of Normal distributions:

where $\Phi_k$ is an $(n_x+1) \times 2 $ matrix and $\Sigma_k$ is a $2\times 2$ matrix. The hyperparameter vector $\xi$ is now defined to include the non-redundant elements of $(\Phi_k,\Sigma_k,\pi_{\lambda,k})$.

For the mixture probabilities $\pi_{\lambda,1:K}$ we use the same prior distribution as in ((ref)). The prior distribution for the coefficient matrices $\Phi_k$ and $\Sigma_k$ is a multivariate generalization of the RE distribution. We assume:

Under this parametrization the marginal IW distribution of the $2 \times 2$ matrix $\Sigma_k$ has mean $D(\tau_\sigma)$. The conditional distribution of $\Phi_k|\Sigma_k$ is $ MN(0,\tau_\phi \Sigma_k \otimes I_{n_z+1})$, where $\tau_\phi$ scales the variance of the Normal distribution. The dimension of $\Sigma_k$ is $2 \times 2$ and, hence, the marginal distribution of $\lambda_i$ is identical to the RE case.\footnote{The marginal distribution of the (1,1) element of the $IW(7,4 D(\tau_\Sigma))$ distribution is $IW(6,4 D_{11}(\tau_\Sigma))$. Converted into the parametrization of the Gamma distribution, this corresponds to an $IG(3,2 D_{11}(\tau_\Sigma))=IG(3,2\tau_\sigma^\lambda)$ distribution.}

{\bf Tuning of the Prior.} The scale of the prior distribution is controlled by a vector of tuning constants: \[ \tau = \big[ \tau_\theta, \tau_\phi, \tau_\sigma^\lambda, \tau_\sigma^y, \tau_v \big]'. \] While these tuning constants could in principle be determined in a data-driven way, using a marginal data density criterion (see the approach used in the Bayesian vector autoregression (VAR) literature, for instance, DelNegro2007b and GiannoneLenzaPrimiceri2015), we do not pursue that route in this paper. Instead we choose $\tau$ ex ante in an informal calibration step. While $\tau_\theta$ has a straightforward interpretation after the regressors have been normalized, the implications of the remaining constants are less transparent because they control priors that are specified over a set of distributions. We recommend the researcher makes an initial choice and then samples from the prior. We found it useful to examine plots of moments or number of modes associated with the distributions. Similar plots can be generated based on the posterior. If a researcher finds that the posterior is located in an area that has essentially no prior mass, then the scaling of the prior can be adjusted to examine whether the initial prior unduly biases the posterior estimates. An example in the context of our empirical application is provided in the Online Appendix.

Posterior Sampling

Draws from the posterior distribution can be obtained with a Gibbs sampling algorithm. We subsequently describe the conditional distributions over which the Gibbs sampler iterates. We focus on the flexible CRE specification with heteroskedasticity, which is the most complicated specification. A key feature of the Gibbs sampler is that it uses data augmentation by sampling the sequences of latent variables $Y_{i,0:T}^*$, $i=1,\ldots,N$. {\color{black} In this regard we are building on TannerWong1987 (data augmentation for a general state-space model), Chib1992 (static Tobit model), albert1993bayesian (Probit model), CarterKohn1994 (linear state space model), and Wei1999 (dynamic Probit model). The general blocking of parameters in the Gibbs sampler is related to baranchuk2008assessing and LiZheng2008. The sampler for the flexible mixture representation of the CRE distribution is based on IshwaranJames2001,IshwaranJames2002.} In terms of the actual implementation, the computations for the Tobit model are very similar to the ones for the linear model studied in liu2018density. The only exception is the treatment of the latent variables $Y_{i,0:T}^*$ which closely follows Wei1999.

In order to characterize the conditional posterior distributions for the Gibbs sampler, we introduce some additional notation. Because $p(\lambda_i,y^*_{i0}|x_{i,-1},\xi)$ and $p(\sigma_i^2|\xi)$ are mixture distributions, {\em ex post} each $(\lambda_i,y_{i0}^*)$ and $\sigma_i^2$ is associated with one of the $K$ mixture components, respectively. We denote the component membership indicators by $\gamma_{i,\lambda}$ and $\gamma_{i,\sigma} \in \{1,\ldots,K\}$, respectively.

{\bf Step 1: Drawing from $Y^*_{i,0:T} | ( Y_{1:N,0:T}, X_{1:N,-1:T}, \lambda_{1:N},\sigma_{1:N}^2,\gamma_{1:N,y},\gamma_{1:N,\sigma},\theta,\xi)$.} To fix ideas, consider the following sequence of observations $y_{i0},\ldots,y_{iT}$: \[ y_{i0}^*, \; y_{i1}^*, \; 0, \; 0, \; 0, \; y_{i5}^*, \; y_{i6}^*, \; 0, \; 0, \; 0, \; y_{i10}^*. \] Our model implies that whenever $y_{it} > 0$ we can deduce that $y_{it}^*=y_{it}$. Thus, we can focus our attention on periods in which $y_{it}=0$. In the hypothetical sample we observe two strings of censored observations: $(y_{i2}, y_{i3}, y_{i4})$ and $(y_{i7}, y_{i8}, y_{i9})$. We use $t_1$ for the start date of a string of censored observations and $t_2$ for the end date. In the example we have two such strings, we write $t_1^{(1)}=2$, $t_2^{(1)}=4$, $t_1^{(2)}=7$, $t_2^{(2)}=9$. The goal is to characterize $p(Y^*_{i,t_1^{(1)}:t_2^{(1)}},Y^*_{i,t_1^{(2)}:t_2^{(2)}} |Y_{i,0:T},\ldots)$. Because of the AR(1) structure, observations in periods $t<t_1-1$ and $t>t_2+1$ contain no additional information about $y^*_{it_1},\ldots,y^*_{it_2}$. Thus, we obtain

which implies that we can sample each string of latent observations independently.

Let $s=t_2-t_1+2$ be the length of the segment that includes the string of censored observations as well as the adjacent uncensored observations. Iterating the AR(1) law of motion for $y_{it}$ forward from period $t_1-1$ we deduce that the vector of random variables $[Y_{i,t_1:t_2}^*,y_{it_2+1}]'$ conditional on $y_{it_1-1}$ is multivariate Normal with mean

The covariance matrix takes the form

We can now use the formula for the conditional mean and variance of a multivariate Normal distribution

to deduce that

Here we use $TN_-(\mu,\Sigma)$ to denote a Normal distribution that is truncated to satisfy $y \le 0$. Draws from this Truncated Normal distribution can be efficiently generated using the algorithm recently proposed by Botev2017.

There are two important special cases. First, suppose that $t_2=T$, meaning that the last observation in the sample is censored. Then the mean vector and the covariance matrix of the Truncated Normal distribution are given by ((ref)) and ((ref)) with the understanding that $s=t_2-t_1+1$. Second, suppose that $t_1=0$, meaning that the initial observation in the sample $y_{i0}=0$. Because in this case the observation $y_{it_1-1} = y_{i,-1}$ is missing, we need to modify the expressions in ((ref)) and ((ref)). According to ((ref)), the joint distribution of $(\lambda_i,y_{i0}^*)$ is a mixture of Normals. Using the mixture component membership indicator $\gamma_{i,\lambda}$, we can express $y^*_{i0}|(\lambda_i ,x_{i,-1}) \sim N(\mu_*(\lambda_i,x_{i,-1}),\sigma^2_*)$. This leads to the mean vector

and the covariance matrix

where the definition of $\rho_{i,j}$ is identical to the definition of $\rho_{i,j|0}$ in ((ref)). One can then use the formulas in ((ref)) to obtain the mean and covariance parameters of the Truncated Normal distribution.

{\bf Step 2: Drawing from $\lambda_i| (Y_{1:N,0:T}, Y^*_{1:N,0:T},X_{1:N,-1:T},\sigma_{1:N}^2,\gamma_{1:N,y},\gamma_{1:N,\sigma},\theta,\xi)$}. Posterior inference with respect to $\lambda_i$ becomes “standard” once we condition on the latent variables $Y^*_{i,0:T}$ and the component membership $\gamma_{i,\lambda}$. It is based on the Normal location-shift model

Because the conditional prior distribution $\lambda_i|(y_{i0}^*,x_{i,-1},\gamma_{i,\lambda})$ is Normal, the posterior of $\lambda_i$ is also Normal and direct sampling is possible.

{\bf Step 3: Drawing from $\sigma_i^2| (Y_{1:N,0:T}, Y^*_{1:N,0:T},,X_{1:N,-1:T},\lambda_{1:N},\gamma_{1:N,y},\gamma_{1:N,\sigma},\theta,\xi)$}. Posterior inference with respect to $\sigma_i^2$ is based on the Normal scale model

However, even conditional on the mixture component membership indicator $\gamma_{i,\sigma}$, the prior for $\sigma_i^2$ in ((ref)) is not conjugate and direct sampling is not possible. Instead, we sample from this non-standard posterior via an adaptive random walk Metropolis-Hastings (RWMH) step.\footnote{We use an adaptive procedure based on AtchadeRosenthalothers2005, which adaptively adjusts the random walk step size to keep acceptance rates around 30%.}

{\bf Step 4: Drawing from $\theta| (Y_{1:N,0:T}, Y^*_{1:N,0:T},,X_{1:N,-1:T},\lambda_{1:N},\sigma_{1:N}^2,\gamma_{1:N,\lambda},\gamma_{1:N,\sigma},\xi)$}. Conditional on the latent variables $Y^*_{i,0:T}$ and the heterogeneous coefficients $\lambda_i,\sigma_i^2$, we can express our model as

The temporal and spatial independence of the $u_{it}$'s allows us to pool observations across $i$ and $t$. Under the Normal prior in ((ref)), the posterior distribution of $\theta = [\rho,\beta']'$ is also Normal and we can obtain draws by direct sampling.

{\bf Step 5: Drawing from $(\gamma_{i,\lambda},\gamma_{i,\sigma})| (Y_{1:N,0:T}, Y^*_{1:N,0:T},,X_{1:N,-1:T},\lambda_{1:N},\sigma_{1:N}^2,\theta,\xi)$.} We describe how to draw the component membership indicator $\gamma_{i,\lambda}$. Straightforward modifications lead to a sampler for $\gamma_{i,\sigma}$. Note that $\xi$ contains the elements of $\Phi_{1:K}$, $\Sigma_{1:K}$, and $\pi_{\lambda,1:K}$. The prior probability that unit $i$ is a member of component $k$ is given by $\pi_{\lambda,k}$. Let $\bar{\pi}_{i,\lambda,k}$ denote the posterior probability of unit $i$ belonging to component $k$ conditional on the set of means $\Phi_{1:K}$ and variances $\Sigma_{1:K}$ as well as $\lambda_i$. The $\bar{\pi}_{i,\lambda,k}$'s are given by

Note that the conditional distribution $\lambda_i| (y_{i0}^*, x_{i,-1}, \Phi_k,\Sigma_k )$ is Normal, indicated by the notation $p_N(\cdot)$, and can be derived from the joint Normal distributions of the mixture components in ((ref)). Thus,

{\bf Step 6: Drawing from $\xi |(Y_{1:N,0:T}, Y^*_{1:N,0:T},X_{1:N,-1:T},\lambda_{1:N},\sigma_{1:N}^2,\gamma_{1:N,\lambda},\gamma_{1:N,\sigma},\theta)$.} Sampling from the conditional posterior of $\Phi_{1:K}$, $\Sigma_{1:K}$, and $\pi_{\lambda,1:K}$ can be implemented as follows. Let $n_{\lambda,k}$ be the number of units and $J_{\lambda,k}$ the set of units that are members of component $k$. Both $n_{\lambda,k}$ and $J_{\lambda,k}$ can be determined based on $\gamma_{1:N,\lambda}$. The conditional posterior of the component probabilities takes the form of a generalized truncated stick breaking process

meaning that the $\zeta_k$'s in ((ref)) have a $B\big(1+n_{\lambda,k},\alpha_\lambda+\sum_{j=k+1}^K n_{\lambda,j}\big)$ distribution. Conditional on $\pi_{\lambda,1:K}$ the hyperparameter $\alpha_\lambda$ has a Gamma posterior distribution of the form

The conditional posterior for $(\Phi_k,\Sigma_k)$ takes the form

Because here the prior $p(\Phi_k,\Sigma_k )$ is MNIW and the likelihood $\prod_{i \in J_{\lambda,k}} p(\lambda_i,y_{i0}^*|x_{i,-1},\Phi_k,\Sigma_k)$ is derived from a multivariate Normal linear regression model, the conditional posterior of $(\Phi_k,\Sigma_k)$ is also MNIW. All three conditional posteriors allow direct sampling. The derivations can be modified to obtain the conditional posterior of $\psi_{1:K}$, $\omega_{1:K}$, and $\pi_{\sigma,1:K}$.

{\bf Step 7: Drawing from the predictive density.} Conditional on $(y_{iT}^*,\lambda_i,\sigma_i^2,\theta)$ and $x_{i,T:T+h-1}$, paths from the predictive distribution for $y_{i,T+1:T+h}$ can be easily generated by simulating ((ref)) forward; see Section (ref) for further details.

{\bf Modifications for the simplified model specifications.} If the CRE distribution is modeled parametrically instead of flexibly, then the drawing of the component membership indicators $(\gamma_{i,\lambda},\gamma_{i,\sigma})$ in Step 5 and the drawing of $\pi_{\cdot,1:K}$ and $\alpha$ in Step 6 are unnecessary. One only has to sample from the MNIW posterior of $(\Phi_1,\Sigma_1)$ and the NIG posterior of $(\psi_1,\omega_1)$. Under homoskedasticity, i.e., $\sigma_i^2=\sigma^2$ for all $i$, we can pool ((ref)) in Step 3 across $t$ and $i$. In combination with the prior in ((ref)) this leads to an IG posterior for $\sigma^2$ from which one can sample directly. The RE specification requires modifications to Step 1, because the distribution of $y_{i0}$ is now simplified to $y_{i0}^* \sim N(\phi_y,\Sigma_y)$, to Step 2 because the prior distribution of $\lambda_i$ is different, and to Step 6 because the pairs of VAR coefficients $(\Phi_k,\Sigma_k)$ are replaced by $(\phi_{\lambda,k},\Sigma_{\lambda,k})$ and $(\phi_y,\Sigma_{y})$, which leads to NIG posteriors.

Multi-Step Forecasting

In general, there are two ways of extending one-step-ahead to multi-step-ahead forecasting: an iterated approach and a direct approach.

First, iterating the law of motion of $y_{it}^*$ in ((ref)) forward by $h$ periods, starting from period $t=T$, yields

{\color{black}Thus, forecasting $y_{iT+h}^*$ iteratively requires the path $x_{i,T:T+h-1}$. We can distinguish the following scenarios: (i) the path is given at time $T$. For instance, in a stress-testing application of our framework the path of the exogenous variables would be specified by the regulator as part of the stressed macroeconomic scenario. (ii) $x_{it}$ is strictly exogenous. In this case the user has to specify a separate model for $x_{it}$ to simulate future trajectories along which ((ref)) is evaluated. Because of the exogeneity, this simulation can be conducted independently of the simulation of ((ref)). Suppose one has draws $(\lambda_i^{(j)},\rho^{(j)}, \beta^{(j)},\sigma^{2(j)}_i)$ and draws $x_{i,T:T+h-1}^{(j)}$ from the posterior predictive distribution of the exogenous regressors, then one can define

One can sample $y_{iT+h}^{*(j)}$ from a $N(\mu_{iT+h|T}^{(j)},\sigma^{2(j)}_{iT+h|T})$ and apply the censoring to obtain a draw $y_{iT+h}^{(j)}$. (iii) The $x_{it}$s are endogenous and interact with the $y_{it}$s, which is the case in our application. To capture the feedback from the dependent variables to the regressors, one has to simulate $(Y_{1:N,T+1:T+h},Y^*_{1:N,T+1:T+h},X_{1:N,T+1:T+h-1})$ jointly; see ((ref))}.

Second, rather than generating $h$-step ahead forecasts iteratively, in practice forecasters often engage in direct estimation of an $h$-step-ahead prediction function. In our framework, this approach amounts to estimating a model of the form \[ y_{it}^* = \lambda_i + \rho y_{it-h}^* + \beta'x_{it-h} + u_{it} \] with the understanding that the serial correlation in $u_{it}$ implied by our original model ((ref)) is ignored. A discussion of the disadvantages and advantages of multi-step estimation in the context of VARs can be found in Schorfheide2005.

Monte Carlo Experiment

{\color{black} We conduct a Monte Carlo experiment to illustrate the performance of the set and density forecasts from the dynamic panel Tobit model in ((ref)) under ideal conditions. We also discuss the estimation of the heterogeneous coefficients. We simplify the model} by omitting the additional predictors $x_{it}$ and using the RE specification. We endow the forecaster with knowledge of the true $p(y^*_{i0})$ and factorize $p(\lambda_i,y^*_{i0},\ln \sigma_i^2|\xi)$ as $p(\lambda_i|\xi) p(y^*_{i0}) p(\ln \sigma_i|\xi)$. The data generating process (DGP) is summarized in Table (ref). We set the autocorrelation parameter to $\rho=0.8$ and consider skewed random effects distributions for $\lambda_i$ and $\ln \sigma_i^2$ that are generated as mixtures of Normals.

The simulated panel data sets consist of $N=1,000$ cross-sectional units and the number of time periods in the estimation sample is $T=10$. We generate one-step-ahead forecasts for period $t=T+1$. The fraction of zeros across all samples is 45% and for roughly 15% of the cross-sectional units the sample consists of $T=10$ zeros (“all zeros”).\footnote{In the Online Appendix we report additional results for Monte Carlo designs with 60% and 75% zeros, respectively. The overall message from the baseline Monte Carlo design is preserved under the alternative specifications.} The measures of forecast accuracy discussed in Sections (ref) and (ref) are first computed for the cross section $i=1,\ldots,N=1,000$ and we then average the performance statistics over the $n_{sim}=100$ Monte Carlo repetitions.

{\bf Model Specifications and Predictors.} We compare the performance of six predictors described below: four Bayes predictors derived from different versions of the dynamic panel Tobit model, a predictor derived from a Tobit model with homogeneous coefficients, and a predictor from a linear model with homogeneous coefficients that ignores the censoring. The prior distributions used for the estimation of the various models were described in Section (ref) and are summarized in Table (ref). Further implementation details are provided in the Online Appendix.

We consider four versions of the dynamic panel Tobit model with random effects (see Section (ref) for details): (i) flexible RE and heteroskedasticity; (ii) Normal RE and heteroskedasticity; (iii) flexible RE and homoskedasticity; and (iv) Normal RE and homoskedasticity. Versions (ii)-(iv) are misspecified in light of the DGP. The pooled Tobit specification ignores the heterogeneity in $\lambda_i$, setting $\lambda_i=\lambda$ for all $i$, and imposes homoskedasticity. Finally, the pooled linear specification imposes $\lambda_i=\lambda$, $\sigma_i=\sigma^2$ for all $i$, and, in addition, ignores the censoring of the observations during the estimation stage (and finally censors the forecasts at 0).

{\bf Density and Set Forecasts.} To assess the density forecasts we compute LPS and CRPS; see Section (ref). The larger LPS and the smaller CRPS the better the forecast. The accuracy statistics are reported in columns 2 and 3 of Table (ref). As expected, the flexible specification with heteroskedasticity that nests the DGP delivers the most accurate density forecasts. While replacing the flexible representations of the RE distributions with Normal distributions only leads to a marginal deterioration of forecast performance, imposing homoskedasticity generates a substantial drop in accuracy. The two “pooled” models that ignore the intercept heterogeneity perform the worst.

We consider two types of set forecasts; see Section (ref). The first type targets the average coverage probability in the cross-section (“average”), whereas the other type targets the correct coverage probability for each unit $i$ (“pointwise”). To assess the set forecasts we compute the coverage frequency and the average length of 90% predictive sets. Results are presented in columns 4 to 7 of Table (ref). The “average” sets constructed from the heteroskedastic specification have good frequentist coverage properties. They attain coverage frequencies of 91.0% and 90.8%, respectively. A comparison between the “average” and the “pointwise” set forecasts from the heteroskedastic models highlights that the average length of the “average” sets is indeed smaller. Moreover, the coverage frequency of the “pointwise” sets exceeds the nominal coverage level of 90% by a larger amount. We observe a similar pattern also for the set forecasts from the homoskedastic model specifications. Overall, the homoskedastic specifications generate worse set forecasts, in terms of coverage frequency {\em and} average length, than the heteroskedastic specifications.

{\bf Parameter Estimates.} The last two columns of Table (ref) summarize the bias and standard deviation of the posterior mean estimator of the homogeneous parameter $\rho$. Under the correctly specified “Flexible & Heterosk.” model the bias is close to zero and the standard deviation is small. Replacing the flexible RE specification by a Normal specification raises the bias by a factor of three. Replacing heteroskedasticity by homoskedasticity approximately increases the standard deviation by 50% because of a loss of efficiency. Imposing intercept homogeneity (pooled Tobit and pooled linear specification) leads to a substantial increase in the bias.

The panels of Figure (ref) show the true RE density $p(\lambda)$, hairlines that represent $p(\lambda|\xi)$ generated from posterior draws of $\xi$, and histograms of the point estimates $\mathbb{E}[\lambda_i|Y_{1:N,0:T}]$. The left panel corresponds to the flexible specification, whereas the panel on the right displays results for the Normal specification. In both cases we allow for heteroskedasticity. The posterior distribution of $p(\lambda|\xi)$ under the flexible specification concentrates near the true density, whereas, not surprisingly, the parametric specification yields larger discrepancies between the true RE density and the draws from the posterior distribution. Because of the shrinkage effect of the prior distribution, we generally expect the cross-sectional distribution of $\mathbb{E}[\lambda_i|Y_{1:N,0:T}]$ (histograms) to be less dispersed than the distribution of $\lambda_i$ (density plots). Moreover, if we observe sequences of all zeros for multiple units $i$, posterior inference of the corresponding $\lambda_i$s should be the same. This will create a spike in the left tail of the $\mathbb{E}[\lambda_i|Y_{1:N,0:T}]$ distribution. Both features are present in the figure.\footnote{\color{black} We provide illustrative analytical examples of these effects in the Online Appendix.}

Empirical Analysis

We now use different versions of the dynamic panel Tobit model to forecast loan charge-off rates (charge-offs divided by the stock of loans in the previous period, multiplied by 400). As mentioned in the introduction, a charge-off occurs if a loan is deemed unlikely to be collected because the borrower has become substantially delinquent after a period of time. The prediction of charge-off rates is interesting from the perspectives of banks, regulators, and investors, because charge-offs generate losses on loan portfolios and are, in fact, a large contributor to bank losses. If these charge-off rates are large, the bank may be entering a period of distress and require additional capital.\footnote{The accounting details are more complicated: bank balance sheets contain a contra asset account called “Allowance for Loan and Lease Losses” (ALLL). Provisions for LLL are created based on estimated credit losses and reduce the income of the bank. Charge-offs reduce the ALLL and the gross loans on the balance sheet, leaving the net amount unchanged. At this stage, the charge-offs do not lead to a further reduction of income. Whether or not a bank takes a loss provision or a charge-off is to some extent a managerial/accounting decision, although regulators require loans they classify as losses to be charged off. We abstract from strategic accounting aspects; see Moyer1990 for a seminal paper.}

We consider a panel of “small” banks, which we define to be banks with total assets of less than one billion dollars.\footnote{Monitoring potential loan losses in small banks is useful by itself. Moreover, the delinquency rates of small banks could foreshadow those rates of large banks since the small banks tend to have more subprime borrowers who are more vulnerable to minor deterioration in economic condition.} For these banks it is reasonable to assume that they operate in local markets. The forecasts are generated from model ((ref)) where $y_{it}$ are charge-off rates. As potential explanatory variables we consider the quarter-on-quarter inflation in the house price index $\Delta \ln \mbox{HPI}_{it-1}$, the change in the unemployment rate $\Delta \mbox{UR}_{it-1}$, and the growth rate in personal income $\Delta \ln \mbox{INC}_{it-1}$. Here $\Delta$ is the temporal difference operator. The term $\beta'x_{it-1}$ therefore captures variation in regional economic conditions which we measure at the state level. Banks located in regions with poor economic conditions may be more likely to encounter loan losses because of a higher fraction of borrowers that are unable to repay their loans. Our baseline model is based on $x_{it} = [\Delta \ln \mbox{HPI}_{it}, \Delta \mbox{UR}_{it}]'$, but we also consider a specification that includes personal income as a third explanatory variable and a specification without any explanatory variables.

The heterogeneous intercept $\lambda_i$ can be interpreted as a bank-specific measure of the quality of the loan portfolio: the smaller $\lambda_i$, the higher the quality of the loan portfolio and the less likely a charge-off is to occur. The autoregressive component in the model captures the persistence of the composition of the loan portfolio over time, and the covariates shift the density of repayment probabilities. We consider various choices of $p \big(\lambda_i,y_{i0}^*,\sigma_i| x_{i,-1},\xi \big)$; see Section (ref). The data set is described in Section (ref). Section (ref) presents density forecast comparisons for various model specifications. Estimates of the heterogeneous and homogeneous parameters are reported in Section (ref). Posterior predictive checks are conducted in Section (ref). Finally, Section (ref) contains the set forecast results.

Data

The raw data are obtained from “call reports” (FFIEC 031 and 041) that the banks have to file with their regulator and are available through the website of the Federal Reserve Bank of Chicago. Due to missing observations and outliers we restrict our attention to four loan categories: credit card (CC) loans, other consumer credit (CON), construction and land development (CLD), and residential real estate (RRE). We construct rolling panel data sets for each loan category that have a time dimension of twelve quarterly observations: one observation $y_0$ to initialize the estimation, $T=10$ observations for estimation, and one observation to evaluate the one-step-ahead forecast. The number of banks $N$ in the cross section varies depending on market size and date availability. The earliest sample considered in the estimation starts ($t=0$) in 2001Q2 and the most recent sample starts in 2016Q1. A detailed description of the construction of the data set is provided in the Online Appendix.

In the remainder of this section, we present two types of results: (i) forecast evaluation statistics and parameter estimates for RRE and CC charge-off rates based on samples that cover the Great Recession and range from 2007Q2 ($t=0$) to 2009Q4 ($t=T$);\footnote{There are, in general, large uncertainties during the Great Recession. Thus, accurate density and set forecasts are important.} (ii) scatter plots summarizing forecast evaluation statistics for the 111 rolling samples that we constructed (based on data availability) for the above-mentioned four loan categories.

Table (ref) contains some summary statistics for the two baseline samples. For the small banks in our sample, RRE loans are an important part of their loan portfolio. For approximately 45% (25%) of the banks RREs account for 20% to 50% (more than 50%) of their loan portfolio. CC loans, on the other hand, make up less than 2% of the loans held by the banks in our sample. Both baseline samples contain a substantial fraction of zero charge-off observations: 76% for RREs and 43% for CC, which makes it challenging to estimate the coefficients of our panel data models. Moreover, 61% of the banks in the RRE sample never write off any loans between 2007 and 2009. The distribution of charge-off rates, across banks and time, is severely skewed. For RREs the 75th percentile is 0 and the maximum is 33.1% annualized. For CCs the corresponding figures are 4.07% and 260%, respectively. A table with summary statistics for the remaining samples is provided in the Online Appendix.

Density Forecasts and Model Selection

{\bf Selected Samples.} We begin the empirical analysis by comparing the density forecast performance of several variants of ((ref)) for the two baseline samples using $x_{it} = [\Delta \ln \mbox{HPI}_{it}, \Delta \mbox{UR}_{it}]'$. This comparison includes forecasts from a Tobit model and a linear model with homogeneous intercepts and homoskedastic innovation variances. Table (ref) reports LPS (the larger the better) and CRPS (the smaller the better). Several observations stand out. First, allowing for heteroskedasticity improves the density forecasts unambiguously. Second, in both RRE and CC samples, all four heteroskedastic specifications lead to very similar density forecasting performance.

{\bf All Samples.} Figure (ref) summarizes the LPS comparisons for all 111 samples. We focus on the comparison of predictive scores from the heteroskedastic specifications versus homoskedastic specifications using flexibly modeled correlated random effects. The solid line is the 45-degree line and the blue and red circles correspond to the scores associated with the baseline RRE and CC samples reported in Table (ref). The figure shows that the results for the baseline samples are qualitatively representative: incorporating heteroskedasticity is important for density forecasting. We provide a figure in the Online Appendix that illustrates that LPS differentials between Normal versus flexible CREs and CREs versus REs are small. In view of these results, we subsequently focus on the flexible CRE specification with heteroskedasticity.

{\bf Tail Probabilities for Selected Samples.} From the density forecasts we can compute probability forecasts for particular events. We consider the tail event $\mathbb{I}\{ y_{iT+1} \ge c\}$ for $c=1$% for now. Figure (ref) visualizes the probabilities of the tail event for RRE charge-off rates for 2010Q1 and 2018Q1, emphasizing the spatial dimension.\footnote{Similar maps for CC charge-off rates are available in the Online Appendix.} We associate each bank $i$ with a particular county. If there are multiple banks in one county, we average the predicted probabilities. 2010Q1 is the immediate aftermath of the Great Recession and the counties that are covered by our sample appear predominantly in dark blue, indicating that predicted probabilities of the event exceed 9.1%. Banks in California, Florida, and the Midwest from Minnesota, Wisconsin, and Michigan down to Arkansas, Mississippi, and Alabama are predicted to write off a considerable fraction of their RRE loans. Eight years later, the situation has improved considerably, as the map now appears in light blue instead of dark blue, in particular in hard hit states such as California and Florida.

While this paper focuses on forecasting problems, the predictive densities derived from our empirical model can be embedded into more complex decision problems that more closely capture the objectives of policy makers or regulators. In this case, the predictive density is used to compute posterior expected losses associated with policy decisions. The accuracy of the loss calculation is tied to the empirical adequacy of the predictive density, which is what we are evaluating in this section.

Parameter Estimates for Selected Samples

{\bf Heterogeneous Parameters.} The distributions of posterior mean estimates of the heterogeneous coefficients for the 2007Q2 sample of RRE charge-off rates are depicted in Figure (ref).\footnote{Similar plots for CC charge-off rates are available in the Online Appendix.} We use the AR coefficient $\rho$ to rescale $\lambda_i$ and $\sigma_i$. The panels on the left and in the center of the figure show histograms for the posterior means of $\lambda_i$ and $\sigma_i$, respectively, whereas the right panel contains a scatter plot that illustrates the correlation between the posterior means of intercepts and shock standard deviations.

A notable feature of the histogram for the posterior means of $\lambda_i/(1-\rho)$ is the spike in the left tail of the distribution. Such spikes were also present in the Monte Carlo simulation; see Figure (ref). The spike corresponds to banks with predominantly zero charge-off rates. For these banks, the sample contains very little information about $\lambda_i$ other than that it has to be sufficiently small to explain the zero charge-off rates. In turn, the posterior mean estimate is predominantly driven by the prior. Similar spikes are visible in the histogram for the posterior means of the re-scaled log standard deviations and the right panel shows that the $\sigma_i$ spike and the $\lambda_i$ spike are associated with the same banks. Small estimates of $\sigma_i$ are associated with near zero estimates of $\lambda_i$, whereas large estimates of $\sigma_i$ are associated with a broad range of $\lambda_i$ estimates. The large dispersion of $\sigma_i$ estimates is consistent with the substantially better density forecast performance of the heteroskedastic models.

Because regional economic conditions have already been controlled for by $\beta'x_{it-1}$, the estimates of $\lambda_i$ are more likely to be related to bank characteristics. Popular explanations for the heterogeneity in loan losses across banks, here captured by the heterogeneity of $\hat{\lambda}_i$, are attitude toward risk, i.e., some banks might have a greater propensity to take risk or have better opportunities to diversify returns on their loan portfolio, and quality of credit management; see KeetonMorris1987 for an early contribution and Ghosh2015,Ghosh2017 more recently.

In Figure (ref) we illustrate the relationship between the posterior mean estimate of $\widehat{\lambda_i/\sigma_i}$, which for $\rho=0$ and $\beta=0$ determines the probability of non-zero charge-offs, and bank size measured by the log of total assets. The top segments of the two panels contain scatter plots with group-wise least-absolute-deviations (LAD) regression lines (left scale). As we have seen previously in Figure (ref) there are two groups of $\hat{\lambda}_i$ estimates. For simplicity, we refer to these groups as low-$\lambda$ and high-$\lambda$ groups respectively. For both the RRE and CC samples the positive relationship between bank size and riskiness of the loan portfolio $\widehat{\lambda_i/\sigma_i}$ is more pronounced for banks in the high-$\lambda$ group. The slope coefficients are 0.18 and 0.19, respectively. The shaded areas at the bottom of the panels contain fitted probabilities (right scale) from a logit model that uses log assets as right-hand-side variable. The larger the assets, the higher the probability that it belongs to the high-$\lambda$ group. These results suggest that larger banks in our sample tend to hold riskier loan portfolios.

In Table (ref) we present estimates from LAD regressions of $\widehat{\lambda_i/\sigma_i}$ on multiple bank characteristics (measured in period $t=0$), separately for the low-$\lambda$ and the high-$\lambda$ group of banks.\footnote{Data definitions and summary statistics for the bank characteristics are provided in the Online Appendix.} We also report estimates for a logit model for $\mathbb{I}\{ i \in \mbox{High} \}$. According to the logit estimates bank size (log assets), the ratio of RRE or CC loans to all loans, lending specialization (ratio of total loans to total assets), and lack of credit quality (ratio of ALLL to total loans) increase the probability that a bank belongs to the high-$\lambda$ group. Capitalization (capital-to-asset ratio) and profitability (return on assets) lower the probability that a bank belongs to the high-$\lambda$ group. For the group-specific regressions only a few bank variables appear to be significant. Foremost, it is bank size measured by log assets. For the RRE high-$\lambda$ group it also includes lending specialization, and for the CC high-$\lambda$ group it includes diversification (share of non-interest income to total income). Operational efficiency, measured by the ratio of overhead costs to assets (OCA) is predominantly insignificant.

Ghosh2017 studies macroeconomic and bank-level determinants of non-performing loans, i.e., loans past due 90 days or more, for the 100 largest commercial banks over the period 1992Q4 to 2016Q1. With the exception of log assets and loan fractions, we followed his study in constructing our bank-level regressors. Although our sample differs from his in several dimensions (selection of banks, measure of loan performance, and time period), we provide a brief comparison of the results for real estate loans as follows.

Ghosh2017 finds the following significant relationships for real estate loans: log capital-to-assets (positive), log loans-to-assets (negative), log inverse credit quality (positive), log return on assets (negative). In our logit regression the same bank characteristics have significant coefficients, but the signs of the estimates for the capital-to-asset and the loan-to-asset ratio differ. As Ghosh2017 points out, the effect of bank capitalization on loan quality is theoretically ambiguous. On the one hand, managers in banks with low capital bases have a moral hazard incentive to engage in risky lending practices (negative relationship). On the other hand, managers in highly capitalized banks may feel confident to engage in risky lending (positive relationship). With respect to the loan-to-asset ratio, our positive estimate for RRE contradicts the notion that banks that are specialized in lending do a better job in selecting high-quality loans, and the positive relationship may reflect that these banks could have more liberal lending policies.

We also report goodness-of-fit ($R^2$) measures in Table (ref). For the LAD regressions we report KoenkerMachado1999's quantile regression $R^2$. For the logit regressions we compute McFadden1973's pseudo $R^2$. For the RRE loans the variation in loan quality ($\widehat{\lambda_i/\sigma_i}$) explained by bank characteristics is low. The $R^2$s for the group-specific LAD regressions are only 0.03 and 0.06, respectively. For the CC sample, bank characteristics are more successful in explaining variations in loan quality. The $R^2$ values are 0.18 and 0.11, respectively. The logit regressions attain pseudo $R^2$ values of 0.32 and 0.47 which indicate that the bank characteristics considered here are partly successful in determining whether a bank belongs to the low-$\lambda$ or high-$\lambda$ group.

{\bf Common Parameters.} Parameter estimates of the common coefficients for the flexible CRE specification with heteroskedasticity are reported in Table (ref) for the 2007Q2 samples. We report posterior means and 90% credible intervals. For each sample we consider three specifications: (i) the baseline specification with $\Delta \ln \mbox{HPI}_{it-1}$ and $\Delta \mbox{UR}_{it-1}$; (ii) an extended version that also includes $\Delta \ln \mbox{INC}_{it-1}$; (iii) and a version without regressors.

Both samples exhibit mild autocorrelation. The point estimate of $\rho$ is 0.21 for RRE and 0.41 for CC. To report the estimates of $\beta$ we undo the standardization of the regressors. The numerical values can be interpreted as follows. For the RRE sample, under the extended specification that includes personal income growth a 1% quarter-on-quarter fall of house prices leads to an increase in charge-off rates by 0.03 percentage points. A 1% increase in the unemployment rate raises the charge-off rates by 0.15 percentage points. Finally, a 1% growth of personal income increases the charge-off rates by 0.001 percentage points. For both samples, the coefficients on persistence, house-price inflation, and unemployment rate changes are “significant," whereas the coefficient on the income growth regressor is “insignificant” in that it is small and its sign is ambiguous. Adding income growth hardly alters the coefficient estimates for house-price inflation and unemployment rate changes. The estimates for the CC sample are qualitatively similar to RRE but about three times larger in magnitude.

In the last column of Table (ref) we report the LPS, now up to four decimal places, that were previously used for the comparison of density forecasts in Table (ref). The values for the three configurations of $x_{it}$ are very close. For the CC sample the LPS criterion favors our baseline specification with $x_{it} = [\Delta \ln \mbox{HPI}_{it}, \Delta \mbox{UR}_{it}]'$, whereas for the RRE sample strictly speaking the model without regressors is preferred. In the Online Appendix we show scatter plots of $\hat{\lambda}_i + \beta'x_{it-1}$ versus $\hat{\lambda}_i$ which indicate that only a very small fraction is explained by local economic conditions. Despite the quantitatively small effect of local economic conditions on charge offs we proceed with $x_{it} = [\Delta \ln \mbox{HPI}_{it}, \Delta \mbox{UR}_{it}]'$, whose coefficients are “significant," and examine the effects of changes in house prices and unemployment more carefully.

Because the Tobit model is nonlinear, the average effect of a change in the regressors (“treatment effect”) depends on $\lambda_i$. We consider a change of the regressor from its sample value $x_{iT}$ to $\tilde{x}_{iT} = x_{iT} + \iota' \Delta x$, where the unit-length vector $\iota$ determines the direction of the perturbation of $x_{iT}$ and $\Delta x>0$ the magnitude. Accounting for censoring, we decompose the treatment effect on $y_{iT+1}$ as follows:

Term $I_i$ captures the intensive margin, i.e., a bank that has non-zero charge-offs conditional on $x_{iT}$ and $\tilde{x}_{iT}$. In this region the Tobit model is linear and the effect is $\beta'\iota$. The second term, $II_i$, captures the extensive margin of banks switching between zero and positive charge-offs.

Figure (ref) depicts the posterior mean and the 90% credible band of the two components of the treatment effect for the banks in the 2007Q2 CC sample.\footnote{A similar figure for the RRE sample is available in the Online Appendix.} We sort the banks based on the posterior means $\widehat{\lambda_i/\sigma_i}$, which for $\rho=0$ and $\beta=0$ would determine the probability of a positive charge-off. We consider two choices for $\iota'\Delta x$: a 5% drop in house prices (left panels) and a 5% rise in the unemployment rate within one quarter (right panels). These are severe shocks to the local economies. For the first approximately 120 banks the posterior mean of $I_i$ (black/grey) is close to zero. These are the banks with low values of $\hat{\lambda}_i$ that appear as a mass in the left tail of the density plot in the left panel of Figure (ref). Under the baseline conditions $x_{iT}$ they are unlikely to have non-zero charge-offs. For the remaining banks the posterior mean of the term $I$ treatment effect rises under the HPI fall scenario from 0.03% to 0.1%, where the latter value is the coefficient estimate reported in Table (ref). The credible intervals are fairly wide, ranging from 0% to 0.15%.

The posterior mean for component $II$ (dark/light blue) of the treatment effect is qualitatively similar under the two economic scenarios. For the first 120 banks term $II$ is small because much of $\beta'(\tilde{x}_{iT}-x_{iT})$ has to compensate for the low estimate of $\lambda_i$ before the latent variable $y_{iT+1}^*$ becomes positive. For the remaining banks the term is also small, but for a different reason: with high probability these banks already have positive charge-offs under the baseline economic conditions. Quantitatively, the effects are larger under the very severe unemployment scenario. The switch of low $\lambda_i$ banks from zero to positive charge-offs leads to a posterior mean of the average treatment effect of 0.04%. As $\widehat{\lambda_i/\sigma_i}$ increases, the expected value of term $II$ decreases because it becomes more likely that the bank has positive charge-offs even under the baseline scenario.

Posterior Predictive Checks for Selected Samples

In order to assess the fit of the estimated panel Tobit model, we report posterior predictive checks in Figure (ref). A posterior predictive check examines the extent to which the estimated model can generate artificial data with sample characteristics that are similar to the characteristics of the actual data that have been used for estimation.\footnote{Textbook treatments of posterior predictive checks can be found, for instance, in lancaster2004 and Geweke2005.} Consider the top left panel of the figure. Here, the particular characteristic, or sample statistic, under consideration is the cross-sectional density of $y_{iT+1}$ conditional on $y_{iT+1}>0$. The black line is computed from the actual RRE loan sample. Each blue hairline is generated as follows: (i) take a draw of $(\rho,\beta,\xi)$ from the posterior distribution; (ii) conditional on these draws generate $\lambda_{1:N}$, $Y_{1:N,0}^*$, and $\sigma^2_{1:N}$; (iii) simulate a panel of observations $\tilde{Y}_{1:N,0:T+1}$; (iv) compute a kernel density estimate based on $\tilde{Y}_{1:N,T+1}$. The swarm of hairlines visualizes the posterior predictive distribution. A model passes a posterior predictive check if the observed value of the sample statistic does not fall too far into the tails of the posterior predictive distribution. Rather than formally computing $p$-values, we focus on a qualitative assessment of the model fit.

By and large, the estimated models for RRE and CC charge-off rates do a fairly good job in reproducing the cross-sectional densities of $y_{iT+1}$ in that some of the hairlines generated from the posterior cover the observed densities. The only discrepancies arise for charge-off values close to zero. With high probability, the densities computed from simulated data have less mass than the observed RRE and CC densities. Moreover, the modes of the simulated densities are slightly to the right and lower than the modes in the two actual densities. The hairlines depict the densities conditional on $y_{iT+1}>0$. In the observed RRE sample the fraction of $y_{iT+1}=0$ is 0.71. The corresponding 90% interval obtained from the estimated model is $[0.73, \, 0.79]$. For CC charge-off rates, the fraction in the data is 0.43 and the corresponding 90% interval obtained from the estimated model is $[0.37, \, 0.47]$.

The center panels of Figure (ref) focus on the estimated models' ability to reproduce the number of zero charge-off observations. For each unit $i$ we compute the number of periods in which $y_{it}=0$. Because $T=10$ the maximum number of zeros between $t=0$ and $t=T+1$ is 12. The histogram is generated from the actual data, whereas the hairlines are computed from the simulated data. For instance, 61% of the banks do not write off any RRE loans in the twelve quarters of the sample and roughly 5% of the banks write off RRE loans in every period. Overall, the estimated models do remarkably well in reproducing the patterns in the data. For RRE loans, the model captures the large number of all-zero samples and the fairly uniform distribution of the number of samples with zero to nine instances of $y_{it}=0$. The only deficiency is that the model cannot explain the absence of samples with ten or eleven instances of zero charge-off rates. In the case of CC loans, the estimated model underpredicts the number of all-zero samples but generally is able to match the rest of the distribution.

The last column of Figure (ref) provides information about the models' ability to capture some of the dynamics of the charge-off data. Here the test statistic is the first-order sample autocorrelation of the $y_{i,0:T+1}$ sequence, conditional on both $y_{it}$ and $y_{it-1}$ being greater than zero. The panels in the figure depict the cross-sectional density of these sample autocorrelations. For the RRE loans the density computed from the actual data is covered by the hairlines generated from the posterior predictive distribution. For the CC loans the estimated model generates somewhat higher sample autocorrelations than what is present in the data.

In the Online Appendix (see Figure (ref)) we consider three additional predictive checks based on (i) the time series mean of $y_{it}$ after observing a zero (and, if applicable, before observing the next zero), (ii) the time series mean of $y_{it}$ before observing a zero (and, if applicable, after observing the previous zero), (iii) a robust estimate of the first-order autocorrelation of $y_{i,0:T+1}$ provided there are sufficiently many non-zero observations. With the exception of the autocorrelations in the CC sample, the two estimated models are able to reproduce the cross-sectional densities of the sample statistics.

Set Forecasts

{\bf Selected Samples.} Set forecasts for 2010Q1, constructed as HPD sets from the posterior predictive distribution, are visualized in Figure (ref). The nominal credible level is 90%. We distinguish forecasts targeting pointwise coverage probability (grey) from forecasts targeting average coverage probability (pink). For each bank $i$ we plot the set forecast, the posterior mean forecast and the actual realization of the charge-off rate. The banks are sorted according to $\mathbb{E}[y_{iT+1}|Y_{1:N,0:T},X_{1:N,-1:T}]$. We don't show forecasts for the first 1,400 (100) banks for the RRE (CC) sample because they are essentially zero.

A comparison of the grey and the pink sets in Figure (ref) shows the effect of targeting average versus pointwise coverage. The upper bound as a function of $i$ increases less under targeting average coverage probability, because the criterion allows us to shorten very wide predictive sets and lengthen narrow sets, while reducing the average length. For the RRE sample set forecasts for banks with large expected charge-off rates $\mathbb{E}[y_{iT+1}|Y_{1:N,0:T},X_{1:N,-1:T}]$ become considerably shorter. In fact for $i>2,500$ many of them become $\{0\}$. Although we plot the actual values of the charge-off rates in Figure (ref), it is not possible to glean how close the empirical coverage frequency is to the nominal coverage probability. Thus, in Table (ref) we report both the average length of the sets and the empirical coverage frequency. For both samples the set forecasts that are constructed by targeting the average coverage probability have a cross-sectional coverage frequency that is close to the nominal coverage probability of 90% and they tend to be shorter than the ones obtained by targeting pointwise coverage probability.\footnote{We also computed evaluation statistics for the homoskedastic specification. It turns out that the set forecasts generated by the homoskedastic specifications are substantially larger than the sets obtained from the models with heteroskedasticity, without improving the coverage probability. This finding is consistent with the density forecast results in Table (ref).}

We also report the frequency of the three types of set forecasts. Due to the large number of zero observations in the RRE sample, there is a large fraction of banks, between 60% and 68%, for which the posterior predictive probability of observing $y_{iT+1}=0$ exceeds 90%. This leads to a forecast of $\{0\}$. For the CC sample the fraction of $\{0\}$ forecasts is considerably smaller.

As one switches from targeting pointwise coverage probability to average coverage probability the composition of the set types changes. Roughly speaking, the forecaster should widen the “narrow” sets (small $\sigma_i$) by lowering their HPD threshold, and tighten the wide sets (large $\sigma_i$) by raising their HPD threshold. For the RRE sample with a relatively high fraction of zeros, when targeting pointwise coverage, the average coverage probability is largely above 90%, so this mechanism manifests itself as reducing wider pointwise sets to $\{0\}$, which decreases the average coverage probability and average length at the same time. Thus, there is an increase in the fraction of $\{0\}$ forecasts; also see the right tail in the left panel of Figure (ref).

For the CC sample with a relatively low fraction of zeros, when targeting pointwise coverage, the average coverage probability is already close to 90%. Switching from targeting pointwise to targeting average coverage, the majority of $\{0\}$ forecasts are converted into $[0,b]$ forecasts by adding a small continuous portion and thereby increasing the pointwise coverage of these units to more than 90%; see the left tail in the right panel of Figure (ref). Moreover, about one third of the disconnected forecasts are converted into connected forecasts, which is due to a lengthening of the sets for small $\sigma_i$ units. In the end, the fraction of $[0,b]$ forecasts increases substantially in this case.

{\bf All Samples.} In Figure (ref) we provide information about the coverage frequency and average length size of the set forecast for all samples. We focus on a comparison between targeting pointwise versus average coverage probability in the flexible CRE specification with heteroskedasticity. Each hairline corresponds to one of the 111 different samples and the two endpoints of the hairlines indicate average length and deviation of the empirical coverage frequency from the 90% nominal credible level. The circled endpoints correspond to targeting average coverage, and unmarked endpoints (or crosses for the baseline samples) represent pointwise coverage targeting. The left panel comprises samples for which targeting average coverage brings the empirical coverage frequency closer to 90% {\em and} reduces the average length. Here the hairlines point into the lower left corner of the graph. The remaining samples are represented by the hairlines in the right panel. Targeting the average coverage unambiguously reduces the average length. For 52% of the samples it also improves the empirical coverage frequency (left panel). For the remaining 48% of the samples the deterioration of the coverage frequency is relatively small. The median improvement in coverage probability in the left panel is 0.022, whereas the median deterioration in the right panel is only 0.007. We conclude that, by and large, directly targeting the average posterior coverage probability improves the empirical coverage frequency in the cross section and produces shorter set forecasts.

Conclusion

The limited dependent variable panel with unobserved individual effects is a common data structure but not extensively studied in the forecasting literature. This paper constructs forecasts based on a flexible dynamic panel Tobit model to forecast individual future outcomes based on a panel of censored data with large $N$ and small $T$ dimensions. Our empirical application to loan charge-off rates of small banks shows that the estimation of heterogeneous intercepts and conditional variances improves density and set forecasting performance in the more than 100 samples considered. Posterior predictive checks conducted for two particular samples indicate that the Tobit model is able to capture salient features of the charge-off panel data sets. Our framework can be extended {\color{black}to allow for stronger forms of simultaneity between the dependent variable and regressors} and to account for dynamic panel versions of more general multivariate censored regression models. We can also allow for missing observations in our panel data set. Finally, even though we focused on the analysis of charge-off data, there are many other potential applications for our methods.

\setstretch{1}

\setstretch{1.3}