Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
78,550 characters · 14 sections · 40 citation commands
Imputation of Counterfactual Outcomes when the Errors are Predictable
Keywords: treatment effect, synthetic control, missing values, blup.
JEL classification: C22, C23
\thispagestyle{empty} \setcounter{page}{0} \baselineskip=18.0pt \setcounter{table}{0}
Understanding the effects of policies is an important aspect of economic analysis and many questions of interest involve an individual's or a group's response to multiperiod interventions. Given treatment status, a researcher observes the outcome of unit $i$ after intervention at $t=T_0+1$ (denoted $Y_{it}(1)$) and wants to compare it to the hypothetical outcome without intervention (denoted $Y_{it}(0)$). Since we do not observe $Y_{it}(0)$ for $t>T_0$, these values need to be imputed. Now imputation is concerned with the prediction of values that will never be sampled, and from results in the prediction literature, we know that in-sample estimation uncertainty should diminish with the sample size and what dominates total prediction error asymptotically is the variation not explained by the model. Yet, in applications, we tend to perform robust inference taking the residuals as given, when an improved prediction is possible by removing predictable variations that might still be in the residuals.
To illustrate, consider Figure (ref) which studies the impact of the German reunification in 1990 on $Y_{1t}=\log$ GDP. Because the GDP data are non-stationary, we estimate common factors from a 16 country panel of GDP growth $(\Delta Y_{2,1:T}(0),\ldots, \Delta Y_{17,1:T}(0))$, where $\Delta Y_{i,1:T}(0)\equiv (\Delta Y_{i1}(0),\ldots, \Delta Y_{iT}(0))'$. The top-left panel displays actual (log) GDP over the full sample, along with the in-sample fit $\hat Y_{1,1:T_0}(0)$, and the counterfactual values $\hat Y_{1,T_0+1:T}(0)$. The effect of reunification on GDP is stark, but masks the fact that the (in-sample) residuals $\hat e_{1,1:T_0}=\Delta Y_{1,1:T_0}(0)-\hat{\Delta Y}_{1,1:T_0}(0)$ are persistent. This can be seen from the correlogram in the bottom left panel, or from the plot of the series itself in the top right panel. The series $\hat e_{1,1:T_0}$ is also cross-correlated with other errors, though many of the $\text{corr}(\hat e_{1,1:T_0},\hat e_{j,1:T_0})$ are not statistically significant as shown in the bottom right panel. The in-sample residuals of the log level model are also serially correlated, as shown in Figure (ref). Time series and cross-section correlation of the in-sample residuals is not specific to this example.
This paper considers the implications of non-spherical errors for model-based imputation. Non-spherical errors, which induce predictability, can arise because the model for $Y_{it}(0)$ is mis-specified or because $Y_{it}(0)$ cannot be adequately captured by observed information without further signal extraction. We build on the best-linear-unbiased predictor (hereafter, blup) developed in goldberger:62 for linear models with non-spherical errors. The key to blup\; is not that it is based on GLS estimation, but that it has a correction term that depends on the covariance structure of the errors. We suggest a practical predictor ( plup) that is asymptotically equivalent to the infeasible blup\; to a first order. Furthermore, if predictability is due to serial correlation, a simple AR(1) correction will reduce the mean-squared prediction error for a large class of stationary mixing error processes. This is not to say that an AR(1) correction is best as the desired adjustment will necessarily be data dependent, but to point out that simple modifications can reduce the mean-squared error of the standard prediction.
We adapt Goldberger's result for linear prediction to an imputation setting when the counterfactual outcomes are never observed. We derive infeasible blup\; for linear panel data models and make precise its dependence on the covariance structure of the errors. We show that when $T_0$ is large, a plup\; that controls for time and/or cross-section correlation can be constructed. But the idea of correcting the standard prediction for predictable errors is more general and can be applied to non-linear models when direct modeling of dynamics is not so straightforward. Thus when linearity is not required, we refer to the practical unbiased predictor as pup.\footnote{We thank Bruce Hansen for this suggestion.}
In addition to inefficient point predictions, ignoring correlation in the errors also has implications for inference. Though a standard prediction will yield asymptotically valid unconditional inference, the prediction interval will be wider because correlated residuals inflate error variance. More concerning is that standard prediction is biased conditional on pre-treatment outcomes of both the treated and the untreated, as well as the post-treatment outcomes of the untreated. This bias may distort inference and the precise impact will depend not just on persistence of the residuals, but also on the realized values of the residuals relevant for imputing $Y_{i,T_0+h}(0)$. In the German unification example, $\hat e_{1,T_0}$ is 1.03 in the growth rate model. This non-zero value yields pup\; growth rates that are slightly different from the standard prediction.
pup\; is concerned with reducing the out-of-sample error and does not preclude the use of robust standard errors or resampling schemes to account for correlation in the in-sample residuals. In practice, a pup\; for unit $i$ will use residuals of unit $i$ before treatment, and possibly of the untreated units after $T_0$. For serial correlation type error dependence that should die off as $h$ increases, a pup\; correction is most effective in imputating $Y_{i,T_0+h}(0)$ at small $h$. Recent work by cwz-jasa:21, fmm-jasa:22, and ferman-jae:23 can be seen from a pup\; perspective.
Our central message is that in-sample uncertainty is asymptotically dominated by variability of the out-of-sample prediction error, and more attention should be paid to improving the point-prediction before turning to inference. The rest of the paper proceeds as follows. Section 2 sets up the econometric framework and provides motivating examples for predictable errors. Section 3 summarizes the properties of blup\; and then presents plup. Its mean-squared error is analyzed using asymptotic expansions of the population prediction error under mixing conditions. Section 4 switches focus from predicting future outcomes to imputing missing values. Unconditional and conditional coverage of prediction intervals are analyzed in Section 5. With some abuse of language, we sometimes use 'imputation' and 'prediction' interchangeably. Our discussion will focus on time dependence but the arguments also hold for spatial and cross-section dependence.
We will use the standard potential outcome framework for analysis. Let $Y_{it}(1)$ be the potential response for unit $i$ at time $t$ if it was exposed to treatment (or policy intervention), and $Y_{it}(0)$ be the potential response of a (control) unit $i$ that was not exposed to intervention at $t$. We observe $Y_{it}=Y_{it}(0)(1-D_{it})+Y_{it}(1) D_{it}$ where treatment status $D_{it}=1$ if unit $i$ is exposed in period $t$ and is zero otherwise. Without loss of generality, we order the $N_1$ exposed units before the $N_0=N-N_1$ unexposed units. We observe \[ Y_{it}=
\] and are interested in the effect on unit $i\in[1,N_1]$ in $h>0$ periods after treatment at $T_0$. Different average treatment effects can be derived from the individual treatment effect, defined as \[ \delta_{i,T_0+h}=\underbrace{Y_{i,T_0+h}(1)}_{\substack{\text{observed outcome }\\ \text{h periods after treatment }}}-\underbrace{Y_{i,T_0+h}(0)}_{\substack{\text{ outcome without treatment}\\ \text{in period } T_0+h}}.\] The econometrics challenge is that $Y_{it}(0)$ is not observed for $i\le N_1$ when $t>T_0$.
Following the literature, we assume that $Y_{it}(0)$ has a pseudo-true conditional mean (or mean-unbiased proxy) $m_{it}=\mathcal M(\beta; \mathcal H)$ that is parameterized by a vector $\beta$ given some information set $\mathcal{H}$, and $e_{it}=Y_{it}(0)-m_{it}$ is such that $E(e_{it})=0$. For example, an AR(1) approximation would make $m_{it}=\rho y_{it-1}$ and $\mathcal H$ would be $y_{is}$ for $s\le T_0$. Being a pseudo-true mean, $m_{it}$ may not coincide with the true conditional mean say, $m^*_{it}$, where $e^*_{it}=Y_{it}(0)-m^*_{it}$. For each $i=1,\ldots, N_1$,
Let $\hat m_{it}$ be a consistent estimate of $m_{it}$. Then \[ Y_{it}(0)=\hat m_{it}+\hat e_{it}. \] Since $\hat Y_{it}(0)=\hat m_{it}$, the treatment effect on unit $i$ at a given $t=T_0+h$ is then estimated by
The pointwise imputation/prediction error is
This error has two sources of variation: one from in-sample estimation of $m_{it}$, and one due to the out-of-sample error $e_{i,T_0+h}$ not captured by the model. The first error will be negligible as $T_0$ increases provided that $\hat m_{it}$ is consistent for $m_{it}$ in some well defined sense, but the second error does not vanish with the sample size and thus total prediction error variance is minimized asymptotically if $m_{it}$ is chosen such $e_{it}$ does not contain predictable information. However, theory actually allows $e_{it}$ to be serially and/or mutually correlated, and while the assumption of no correlation is convenient, it is not always appropriate. In the next subsection, we provide some examples for dependence in the errors.
We will first clarify what we mean by in-sample and out-of-sample errors. To fix ideas, suppose that unit 1 is being treated and the model is linear so that $m_{1t}=x_{t}'\beta$. Single equation estimation yields the imputed value $\hat \delta_{1,T_0+1}=x_{T_0+1}'\hat\beta$ and imputation error $ \hat \delta_{1,T_0+1}-\delta_{1,T_0+1}= -x_{T_0+1}'(\hat\beta-\beta) +e_{1,T_0+1}$ whose variance is \[\text{var}(\hat \delta_{1,T_0+1}-\delta_{1,T_0+1})=\sigma_e^2+ x_{T_0+1}' \text{var}(\hat\beta) x_{T_0+1}.\] Correlation in $x_te_t$ may necessitate robust standard errors for $\hat\beta$, but provided that $E[x_te_{1t}]=0$, $\hat\beta$ is consistent in the sense that $\text{var}(\hat\beta)\rightarrow 0$ as $T_0\rightarrow \infty$. Thus, the variance of imputation error is dominated by the out-of-sample error variance $\sigma^2_e\equiv\text{var}(e_{1,T_0+1})$ asymptotically. This variance is minimized when $e_{1,T_0+1}$ is uncorrelated. Serial correlation can arise because of temporal aggregation, but residual correlation (temporally or mutually) is usually a symptom of misspecification of the model or conditioning information. We give some examples below.
\paragraph{Example misspecification 1:} Suppose that $Y_{1,t}(0)=\phi_1 Y_{1,t-1}(0)+\phi_2 Y_{1,t-2}(0)+v_{1t}$ is an AR(2) process with iid innovations $v_{1t}$, but the researcher assumes an AR(1) model. Then $m_{1t}= \beta Y_{1,t-1}(0)$, the pseudo true parameter is $\beta=\frac{\phi_1}{1-\phi_2}\ne \phi_1$, and $e_{1t}=v_{1t}+(\phi_1-\beta) Y_{1,t-1}(0)+\phi_2 Y_{1,t-2}(0)$ is serially correlated when $\phi_2\ne 0$.
\paragraph{Example misspecification 2:} Suppose that the potential outcome has an interactive fixed effect structure: $Y_{1t}(0)=\lambda_1'F_t+\epsilon_{1t}$ but a researcher specifies an additive fixed effect model $Y_{1t}(0)=\lambda_1+F_t+e_{1t}$. Then $e_{1t}$ will be serially correlated if $F_t$ is serially correlated, even if $\epsilon_{1t}$ is white noise.
Factor-based imputation assumes $X_{it}=\lambda_i'F_t+e_{it}$ where $F_t$ is a vector of $r$ latent common factors with $\lambda_i$ as loadings and $e_{it}$ is an idiosyncratic error. An appeal of factor-based imputation is that under some conditions, the space spanned by $F$ can be consistently estimated without modeling the (weak cross-section or time) dependence in the idiosyncratic errors. While xu:17 iteratively estimates $F$ and the missing values jointly by principal components (PCA), baing-jasa:21 impute the missing values (or complete the matrix) using two full-sample applications of PCA. athey-jasa:21, syndid-aer:21 estimate the low rank component using singular value thresholding (SVT).\footnote{Regularization is not necessary to consistently estimate the missing values, but could give a lower rank common component than the one in baing-jasa:21.} Though all consistent estimators of $F$ imply that $\hat F_{T_0+h}$ can be used as though they were observed predictors, the possibility remains that $e_{1,T_0+h}$ can be predicted by information available.
\paragraph{Example Correlated Idiosyncratic Errors} (from fmm-jasa:22) Suppose that $Y_{it}(0)=\lambda_i'F_t+e_{it}$ and the researcher correctly assumes $m_{1t}=\lambda_1'F_t$ but $e_{it}$ is correlated with $e_{jt}$ for $j$ in some index set $C$. Then $e_{1t}=v_{1t}+\sum_{j\in C} \theta_j e_{jt}$ is predictable by those $e_{jt}$ where $j\in C$.
The method of synthetic control (SC) developed in abadie:03 assumes that there exist weights $\beta_j^*$ such that a perfect fit $Y_{1t}(0)=\sum_{j=2}^{N} \beta^*_j Y_{j,t}(0)$ exists for every $t\le T_0$. abadie-jasa:10 make additional use of $K$ economic predictors $X_t=(X_{1t},X_{0t})$, where $X_{0t}=(X_{2t},\ldots, X_{N,t})$ for the unexposed. The Synthetic Difference-in-Difference (SDID) in syndid-aer:21 also reweights the pre-treatment time periods to balance the pre-and post exposure time periods and nests SC and DID as special cases. However, an increasing number of papers suggest that an `imperfect pretreatment fit' may prevent recovery of $\beta^*$.
\paragraph{Example Imperfect Fit 1:} (from bfr:21) Suppose that $Y_{it}(0)=\phi_1 Y_{i,t-1}(0)+v_{it}$ for all $i=1,2,\ldots,N$, and one constructs $m_{1t}=\sum_{j=2}^{N}\beta^*_j Y_{j,t}$. We can show that the error $e_{1t}=Y_{1t}(0)-m_{1t}$ can be decomposed as $e_{1t} =\phi_1 e_{1,t-1} +v_{1t}-\sum_{j=2}^{N}\beta^*_j v_{jt}$. The error $e_{1,T_{0}+1}$ contains an imbalance component $\phi_1 e_{1,T_0}=\phi_1 (Y_{1,T_0}-\sum_{j=2}^N \beta^*_j Y_{j,T_0})$ (which is zero if there is perfect fit but not otherwise) as well as a noise component $v_{1,T_0+1}-\sum_{j=2}^{N} \beta^*_j v_{j,T_0+1}$, and both can contribute to serial correlation.
\paragraph{Example Imperfect Fit 2:} (from ferman-pinto:21) Suppose that $Y_{it}(0)=c_i+\Lambda_i'F_t+\epsilon_{it}$ and one estimates $\hat\beta=\text{argmin}_{b} \|Y_1(0)-X_0b)||_2^2$ where $X_0=(Y_2,\ldots, Y_N)$. With $\beta^*=\text{plim} \hat\beta$, the population imputation error is \[ e_{1t}=Y_{1t}(0)-X_{0t}'\beta^*=\bigg(c_1-\sum_{j=2}^{N}\beta^*_j c_j \bigg)+F_t'\bigg(\lambda_1-\sum_{j=2}^{N} \beta_j^* \lambda_j\bigg)+\bigg(\epsilon_{1t}-\sum_{j=2}^{N} \beta^*_j \epsilon_{jt}\bigg).\] The first two terms vanish only if $\sigma_{\epsilon}\equiv \text{var}(\epsilon_{it})= 0$; otherwise, $(c_1,\lambda_1)\ne (\sum_{j=2}^N \beta_j^*c_j, \sum_{j=2}^N \beta_j^*\lambda_j)$. ferman-pinto:21 suggest to remove the bias with a mean adjustment but this may not remove serial or mutual correlation in the errors.
In the above examples, $e_{it}$ absorbs all sorts of deficiencies in $m_{it}$ and thus contains information about $Y_{it}(0)$. Cross-section, spatial, and time dependence in $e_{it}$ are examples of non-spherical errors.
This section uses classical results in linear prediction to motivate how information in the errors can be used to improve prediction. We will consider optimal linear prediction in a time series setting so that the $i$ subscript can be dropped.
This subsection summarizes results for best linear unbiased prediction, blup. The concept seems to be first considered in henderson:50 in the animal breeding literature to predict the quality of offsprings. It is still widely used in estimation of random effects in linear mixed models for cross-section data.\footnote{robinson:91 provides a survey of its many derivations, including a Kalman filter interpretation, see also spall:91. taub:79 and baltagi-jf:08,baltagi-handbook use it in variance components analysis of panel data.} goldberger:62 formalizes the idea in a setting where the linear model for predicting a scalar variable $y_t$ is given by
where $X_t$ is a $K\times 1 $ vector of completely observed predictors assumed to be fixed in repeated samples, $\beta$ is a vector of time invariant parameters, $e_{t}$ is a zero mean stationary process that is possibly serially correlated, and $ \Omega$ is the $T_0\times T_0$ covariance matrix of the $T_0\times 1$ vector $e$.
goldberger:62 is interested in a linear unbiased prediction of $y_m$ at some $m>T_0$ given information up to $T_0$ when $\Omega$ is positive definite but has non-zero off-diagonal entries. Let $X$ be at $T_0\times K$ matrix of regressors. The assumption of squared loss $E[(y_{m}-y_{m|T_0})^{2}]$ implies a linear predictor of the form $ y_{m|T_0}=A^{\prime }y$ with prediction error $y_{m|T_0}-y_{m}=(A^{ \prime }X-X_{m}^{\prime })\beta +A^{\prime }e-e_{m}$. The unbiasedness constraint $E[y_{m|T_0}]=y_{m}$ requires that $ A^{\prime }X-X_{m}^{\prime }=0$, implying a prediction variance of
where $\omega \equiv E[e_{m}e]$. Let $\lambda$ be the Lagrange multiplier on the unbiasedness constraint. Minimizing $A^{\prime }\Omega A-2A^{\prime }\omega -2\lambda ^{\prime }(X^{\prime }A-X_{m})$ with respect to $A$ gives the best linear unbiased prediction ( blup)
where $\beta _{GLS}=(X^{\prime }\Omega ^{-1}X)^{-1}X^{\prime }\Omega ^{-1}y$ is the infeasible GLS estimator, and ${e}_{GLS}=y-X\beta _{GLS}$ is a $T_0\times 1$ vector of errors. Notably, $ y^*_{m|T_0}$ depends on assumptions about $\Omega $, and in a time series setting, this depends on the dynamics of $e_t$. If $e_{t}=\phi_{1} e_{t-1}+v_{t}$, where $v_t\smash{\mathop{\sim}\limits^d} (0,\sigma^2_v)$, $\Omega$ is $\sigma^{2}_{e}=\sigma^2_{v}/(1-\phi_{1}^2)$ times a $T_0\times T_0$ Toeplitz matrix with $\phi_{1}^{i-1}$ on the $i$-th diagonal. Then $\omega \equiv E[e_{m}e]=\phi_{1} ^{-T_0+m}\Omega _{T_0}$, where $\Omega _{T_0}$ is the last column of $\Omega $. The AR(1) assumption implies a blup\;at $m=T_0+h$ of
with prediction error $e^*_{T_0+h|T_0}=y_{T_0+h}-y_{T_0+h|T_0}^{\ast }=e_{T_0+h}-\phi^h_{1} e_{T_0} +o_p(1)= v_{T_0+h}+o_p(1)$ where the $o_p(1)$ term converges to 0 as $T_0\rightarrow\infty$.
blup\; is infeasible because $\phi_{1} $ is not observed. Feasible blup\; requires iterative Cochrane-Orcutt or Prais-Winsten estimation of $\phi_{1}$, or direct estimation of $\beta$ and $\phi_{1}$ from a Durbin equation.\footnote{Cochrane-Orcutt performs least squares regression of $y_{t}-\phi_{1} y_{t-1}$ on $X_{t}-\phi_{1} X_{t-1}$ for given $\phi_{1} $ using data from $t=2,\ldots ,T_0$, and then estimates $\phi_{1} $ from an autoregression in $y_{t}-X_{t}^{\prime }\hat{\beta}$ till convergence. The Prais-Winsten estimator additionally exploits information in $t=1$. It is also possible to estimate $\beta$ and $\phi_{1}$ directly from the Durbin equation $y_{t}=X_{t}^{\prime }\beta+ y_{t-1}\phi_{1}+X_{t-1}^{\prime }{\gamma}+$ error.} These feasible estimators are all efficient and consistent. cochrane-orcutt:49 suggest to improve the standard prediction by incorporating lags of regressors and the dependent variable. The difference is that $\hat e_{GLS,t}$ summarizes the dynamic relation between $y$ and $X$ into a single signal and can be appealing when $X$ is of high dimension.
A feasible blup\;differs from the OLS prediction in two ways. First, it uses $\hat{\beta} _{GLS}$ instead of $\hat{\beta}_{OLS}$ and thus requires the dynamics of $e_t$ to be specified. Second, blup\; adds to the GLS prediction a term that adjusts for serial correlation in $e$ which in this AR(1) example is $\phi_1 e_{T_0}$ for $h=1$. Since blup\; is an optimal prediction, it is more efficient than an OLS prediction. To make this point precise, consider again the AR(1) case. At $h=1$, feasible blup\;
has prediction error $\hat e_{T_0+1|T_0}^*=y_{T_0+1}-\hat y_{T_0+1|T_0}^*$, or
where the $o_{p}(1)$ term comes from the fact that the jointly estimated $ \hat{\beta}_{GLS}$ and $\hat\phi_{1,GLS}$ are $\sqrt{T_0}$ consistent for $ \beta $ and $\phi_{1} $. Since $\hat e_{T_0+1|T_0}^*$ is asymptotically $v_{T_0+1}$ whose variance is $\sigma^2_v$, feasible blup\; achieves the same asymptotic efficiency as infeasible blup.
In contrast, the OLS prediction error is
where the $o_p(1)$ term comes from $\sqrt{T_0}$ consistency of $\hat{\beta}_{OLS}$ for $\beta$. But $\hat e_{T_0+1|T_0}$ is asymptotically $e_{T_0+1}$ whose variance is $\sigma^2_e\ge \sigma^2_v$. Thus, the MSE improvement of blup\; over OLS is due to the additional term $\phi_{1} e_{T,GLS}$ in the prediction, not because of GLS versus OLS estimation per se. Building on this idea, we will consider a linear prediction that is also asymptotically unbiased but can improve upon the OLS prediction without a priori knowledge of the precise dynamic structure of $e_t$.
This subsection suggests a practical variant ( plup) and studies its mean-squared error (MSE) using asymptotic expansions, first for $h=1$, and then for $h>1$ when direct and iterative forecasts are possible.
Our point of departure is that any predictor that controls for serial correlation will have the same first order effect as feasible blup. Let $ \hat\beta$ denote the least squares estimate of $\beta$. Consider modifying the (standard) least-squares prediction $\hat y_{T_0+1|T_0}=X_{T_0+1}^{\prime}\hat{ \beta}$ as follows:
where $\hat{e}_{T_0}=y_{T_0}-X_{T_0}^{\prime}\hat{\beta}$ is the OLS residual, and
is the least squares estimate of the first order autocorrelation coefficient of $\hat e_{t}$. Note that unlike feasible blup\;which re-estimates $\beta$ after $\hat\rho_{1}$ is available, we simply adjust the OLS prediction $\hat y_{T_0+1|T_0}$ for serial correlation with $\hat\rho_1\hat e_{T_0}$. The prediction error $\hat e_{T_0+1|T_0}^+= y_{T_0+1}-\hat{y}^+_{T_0+1|T_0} = y_{T_0+1}-\hat y_{T_0+1|T_0}-\hat\rho_{1}\hat e_{T_0}$ is
where the last equality follows because $\hat\beta\smash{\mathop{\longrightarrow}\limits^p} \beta$ and $\hat \rho_{1}\smash{\mathop{\longrightarrow}\limits^p} \rho_1$. If $e_t$ is indeed an AR(1) model, then $ \rho_{1}=\phi_{1}$ and
which is asymptotically equal to $v_{T_0+1}$, the prediction error of the infeasible blup.
As it turns out, adding the term $\hat{\rho}_{1}\hat{e}_{T_0}$ to any consistent prediction $\hat{y}_{T_0+1|T_0}$ will yield an efficiency gain even when the true model is not an AR(1). We will refer to the prediction
as practical blup\; (or plup), practical because it does not require GLS estimation and it is asymptotically as efficient as blup. To formalize the properties of the plup\; error $\hat{e}_{T_0+1|T_0}^{+}=y_{T_0+1}-\hat{y} _{T_0+1|T_0}^{+}$, we assume the following.
We define the strong mixing coefficients as $\alpha \left( k\right) =\sup_{A,B}\left\vert P\left( A\cap B\right) -P\left( A\right) P\left( B\right) \right\vert $ where $A$ and $B$ vary over events in the sigma fields generated by $\left\{ e_{s}:s\leq 0\right\} $ and $\left\{ e_{s}:s\geq k\right\} $, respectively. Assumption A1 includes linear processes $ e_{t}=\sum_{j=0}^{\infty }\psi _{j}v_{t-j}$, where $\sum_{j=0}^{\infty }\left\vert \psi _{j}\right\vert <\infty $ and $v_{t}$ is i.i.d.$\left( 0,\sigma _{v}^{2}\right) $ with $E\left\vert v_{t}\right\vert ^{r}<\infty $, which includes stationary invertible ARMA(p,q) processes and nonlinear weakly dependent processes with GARCH and ARCH innovations.
Part (i) shows that $\hat{\rho}_{1}$ converges to the first order autocorrelation coefficient of $e_{t}$. Parts (ii) and (iii) describe the asymptotic expansion of the prediction errors ignoring the estimation error uncertainty. Part (ii) implies that the standard prediction is asymptotically unconditionally unbiased in spite of not accounting for serial correlation because $ e_{T_0+1}$ is mean zero by assumption, and it has asymptotic variance $\gamma _{0}\equiv \text{var}\left( e_{T_0+1}\right) =\sigma _{e}^{2}$. The plup\; error in (iii) also has an unconditional mean of zero, but its variance is $\gamma _{0}(1-\rho _{1}^{2})$. Since $\left\vert \rho _{1}\right\vert \leq 1$,
implying that the plup\; mean-squared dominates the standard prediction. If $e_{t}$ is truly generated as $e_{t}=\phi _{1}e_{t-1}+v_{t}$, the plup\;error $\hat{e}_{T_0+1|T_0}^{+}$ will be asymptotically serially uncorrelated since $\rho _{1}=\phi_1$.
However, an AR(1) correction will improve upon the standard prediction even when $e_t$ is not an AR(1), provided that $e_t$ is a mixing process satisfying A1. For instance, if $e_{t}$ is an AR(2) defined by $e_{t}=\phi _{1}e_{t-1}+\phi _{2}e_{t-2}+v_{t}$, then $\rho _{1}=\frac{\phi _{1}}{1-\phi _{2}}\neq \phi _{1}$. But it will still be the case that $\hat{\rho}_{1}\smash{\mathop{ \longrightarrow}\limits^p}\rho _{1}$ as stated in (i). The plup\; error is now $\hat{e}_{T_0+1|T_0}^{+}=e_{T_0+1}-\rho _{1}e_{T_0}+o_{p}(1)=(\phi _{1}-\rho _{1})e_{T_0-1}+\phi _{2}e_{T_0-2}+o_{p}(1)$, while the standard prediction error is $\hat{e}_{T_0+1|T_0}=\phi _{1}e_{T_0-1}+\phi _{2}e_{T_0-2}+o_p(1)$. Both have a mean of zero, implying that misspecifying the dynamics will not contribute to unconditional bias. Nonetheless, the plup\; variance is $(1-\rho _{1}^{2})\sigma _{e}^{2}$, which is smaller than the standard prediction error variance of $ \sigma _{e}^{2}$ since $\left\vert \rho _{1}\right\vert \leq 1$. Thus the AR(1) correction unambiguously reduces one-step asymptotic mean squared prediction error. In theory, an AR(p) correction with known parameters should improve prediction when $e_t$ is an AR(p). But in practice, sampling uncertainty may offset some gains. Furthermore, when the assumed AR(p) is not the true model, the mean-squared error is no longer tractable as shown in kunitomo-yamamoto:85 even without sampling error. The AR(1) correction is appealing because it is simple to implement, and precise mean-squared error statements can be made when there is no sampling uncertainty as stated in Lemma (ref).
Next, consider cases when $h>1$. The standard prediction $\hat{y}_{T_0+h|T_0}=X_{T_0+h}^{\prime }\hat{ \beta}$ has error $\hat{e}_{T_0+h|T_0}=y_{T_0+h}-\hat{y}_{T_0+h|T_0} =e_{T_0+h}+o_{p}(1)$, where $E[e_{T_0+h}]=0$ and $\text{var}(e_{T_0+h})=\gamma _{0}$. There are two ways to implement plup. The first to use the AR(1) model for $e_t$ to iteratively predict $e_{T_0+h}$. {\em Iterated plup} (or plupi), defined as
has error $\hat{e}_{T_0+h|T_0}^{+I}=e_{T_0+h}-\rho _{1}^{h}e_{T_0}+o_{p}(1)$. As shown in the Appendix, \[e_{T_0+h}-\rho _{1}^{h}e_{T_0}\smash{\mathop{\sim}\limits^d} (0, \gamma _{0}[1+\rho _{1}^{2h}-2\rho _{1}^{h}\rho _{h}]).\] The second approach is to directly predict $e_{T_0+h}$ using information up to $T_0$. Let $\hat{\rho}_{h}=\hat{\gamma}_{h}/\hat{\gamma}_{0}$ be the $h^{th}$ order sample autocorrelation coefficient of $\{\hat{e}_{t}:t=1,\ldots ,T_0\}$. Direct plup\; (or plupd) is defined as
plupd\; has error $\hat{e}_{T_0+h|T_0}^{+d}=e_{T_0+h}-\rho _{h}e_{T_0}+o_{p}\left( 1\right) $, where \[ e_{T_0+h}-\rho _{h}e_{T_0}\smash{\mathop{\sim}\limits^d} (0, \gamma _{0}\left( 1-\rho _{h}^{2}\right)) .\]
The proof is given in the Appendix. The standard predictor, plupi\; and plupd\; are all asymptotically unbiased provided that $E[e_t]=0$. But
with equality when the AR(1) model is correctly specified. Furthermore, \[\text{var}(e_{T_0+h}-\rho _{h}e_{T_0})-\text{var}(e_{T_0+h})= -\rho_h^2 \gamma_0 \le 0\] with equality when $\rho_h=0$. Hence, absent sampling uncertainty, the asymptotic MSE of plupd\; using $e_{T_0}$ to improve the standard prediction can be no larger than plupi\; which uses the same information for correction, or the standard predictor which ignores $e_{T_0}$ for any $h\ge 1$. This is not to say that a richer dynamic model would not produce further improvements. What's noteworthy is that even a simple correction will reduce the prediction MSE at any $h$. We can also expect the plupd\; gains to be largest at $h=1$ and diminish with $h$ because the long horizon forecast of a covariance stationary process is the unconditional mean. Though precise statements can be made for plupd, we can only say that the asymptotic MSE of plupi\; is smaller or equal than that of the standard predictor if $1+\rho _{1}^{2h}-2\rho _{1}^{h}\rho _{h}\leq 1$ (see the Appendix for a proof).
Once a prediction is made, we can construct prediction intervals. We will be studying plup\;based inference under normality in the context of causal inference. As we will see in Section 5, while the standard prediction is unconditionally unbiased, it is conditional biased and conditional inference will, in general, have the wrong coverage.
This subsection evaluates the unconditional and conditional prediction bias, MSE, and coverage with and without plup\;correction in a single equation setting where by unconditional inference, we mean that $e_{i,T_0}$ is random in repeated sampling, and by conditional inference, we mean that $e_{i,T_0}$ is treated as fixed with respect to some conditioning information, as would be the case in practice.
In each of the 5000 replications, we first simulate $K=2$ regressors and $y_t=X_t'\beta+e_t$ where for $T=1,\ldots, 200$, $e_t=\phi_{1} e_{t-1}+\phi_{2} e_{t-2}+v_{t}$. With $v_t\smash{\mathop{\sim}\limits^d} N(0,.05)$, the $R^2$ of the regression is about 2/3. In Case 1, $e_{t}$ is an AR(1) with $\phi_1=.8$, and in Case 2, $e_t$ is an AR(2) with $(\phi_1,\phi_2)=(1.3,-.4)$. Table (ref) reports four sets of errors in predicting $y_{T_0+h}$. The column labeled 'best' is the infeasible prediction when $\beta ,\phi_1,\phi_2$, and $ e_{T_0} $ are known. The column labeled 'noadj' is also infeasible but unlike 'best', it does not take into account information in $e_{T_0}$. The column labeled 'ols' is the standard prediction $ \hat y_{T_0+h}$ using the least squares estimate $\hat\beta$. The columns plupi\;and plupd\;are iterative and direct plup\; respectively. Both are based on a simple AR(1) correction, ie. even when the true DGP is AR(2). Note that they are identical when $h=1$.
The top panel of Table (ref) reports the unconditional bias and MSE for horizons $h=1,2,5,10$. The average over all 10 horizons is reported in the row labeled 'avg'. Since $E[e_{T_0}]=0$, the unconditional prediction bias is close to zero. However, the unconditional prediction MSE is much smaller with plup\;corrections. In the AR(1) case, the MSE for the standard (OLS) prediction is 0.14 at $h=1$, but the plup\;corrections reduce the MSE to 0.05. In the AR(2) case, the OLS prediction has an MSE of 0.43 while the plup\;corrections reduce the MSE to 0.06. The MSE improvements are smaller when $h>1$, but still non-trivial.
The middle panel of Table (ref) shows conditional prediction errors when $(e_{T_0-1},e_{T_0})$ are fixed to $(0.5, 1)$. All predictions are conditionally biased, but the plup\; biases are significantly smaller. When $h=1$ and the errors are AR(1), the standard prediction has a conditional bias of 0.79 while the plup\;corrections reduce it to 0.01. When the errors are an AR(2) process but an AR(1) correction is implemented, the conditional OLS bias at $h=1$ is reduced from 1.09 to 0.18. Correspondingly, the MSE is reduced from 1.26 to 0.09. Note that the biases are largest when $h=1$ because the predictability of a stationary ergodic process decreases with the forecast horizon. The improvements in MSE at $h=1$ translate into improved average predictions over 10 periods. Without the corrections, the average prediction in the AR(2) case has a bias of 0.61 and an MSE of 0.56. The AR(1) plup\; direct correction reduces bias to 0.17 and MSE to 0.22.
Imputation concerns prediction of values that are never observed. The problem is widely studied in a static setup, but there are few results for a dynamic setting. little-rubin:19 consider an AR(1) model where $y_1,y_3,\ldots, y_{T-1}$ are observed but not $y_2,y_4,\ldots, y_T$. The adjustments, shown to require an implicit regression of $y_t$ on $y_{t-1}$ and $ y_{t+1}$, can be seen as smoothed estimates of a suitably defined Kalman filter. chow-lin:71 consider missing values occurring between two releases of low frequency data and show that the best prediction involves a correction term that has a blup\; form. ng-scanlan:24 consider factor-based imputation of weekly missing values of a scalar series occurring throughout the sample.
Causal inference concerns imputation of missing potential outcomes that tend to occur at the end of the sample. The problem is typically studied for an iid setting when it is natural to assume that the errors are uncorrelated\footnote{causalimpact:15 consider state space estimation of the counterfactual outcomes in the presence of trends, but serial correlation in idiosyncratic shocks and/or the factors are not allowed. carvalho-masini-medeiros:18,masini/medeiros-jasa:21,masini-medeiros:22 consider causal inference in a high-dimensional setting when the data are persistent and possibly non-stationary.}. As suggested in Section 2, correlation in the residuals cannot be ruled out. We will consider the imputation problem from the perspective of optimal prediction, with the goal of using the insights of blup\; to improve the imputation of $Y_{i,T_0+h}(0)$. We assume that $e_{it}=Y_{it}(0)-m_{it}$ are strong mixing processes and rule out non-stationary data. In addition, we impose the following high level assumption: \paragraph{Assumptions A2:} For $h\ge 1$, $\hat{m}_{i,T_0+h}-m_{i,T_0+h}=o_{p}\left( 1\right) $ and $T_{0}^{-1}\sum_{t=1}^{T_{0}}\left( \hat{m}_{it}-m_{it}\right) ^{2}=o_{p}\left( 1\right) $.
Assumption A2 is verified in cwz-jasa:21 for estimators including synthetic control, matrix completion, factor-based methods. Given an asymptotically unbiased $\hat m_{i,T_0+h}$ satisfying Assumptions A1 and A2, the estimated treatment effect
has error
This error has two components: an in-sample estimation uncertainty component that depends on the estimator but vanishes as $T_0 \rightarrow \infty$, and an out-of-sample prediction component that depends on the choice of $m_{it}$ and the information $\mathcal H$ used in the imputation.
In order to extend Goldberger's blup\; from a complete data setting to a potential outcomes setting, define the $n\times 1$ vector of (observed) control outcomes $\mathcal Y(0)$ by \[ \underbrace{\mathcal Y(0)}_{n\times 1}=
\equiv
\] where $n=(NT -N_1T_1)$. Note that $\mathcal Y(0)$ includes not only the pre-intervention outcomes on all units $\mathcal Y^{pre}(0)\equiv (Y_{1,1:T_0} \ldots Y_{N,1:T_0})'$ but also the post-intervention outcomes on the control units $\mathcal Y^{post}(0)\equiv (Y_{N_1+1,T_0+1:T}\ldots Y_{N,T_0+1:T})'$.
Let $\mathcal M$ be the pseudo-conditional mean for $\mathcal Y(0)$ and $\mathcal E$ be the corresponding errors. The matrices $\mathcal M$ and $\mathcal E$ contain typical elements $m_{it}$ and $e_{it}$, respectively. With this notation,
The $n\times n$ matrix $\Gamma$ depends on $\Sigma=E[e_te_t']$, which is the $N\times N$ covariance matrix of $e_t=(e_{1:N_1,t}',e_{N_1+1:N,t}')'$. It also depends on $\Omega_i=E[e_ie_i']$, which is the time series covariance of unit $i$, and its dimension can be $T\times T$ or $T_0\times T_0$ depending on whether $i$ is treated.
Consider the linear case $\mathcal M =\mathcal X\beta$, where $\mathcal X$ contains observed predictors. Consider obtaining blup\; given information on $\mathcal Y(0)$ and $\mathcal X$. Re-doing Goldberger's problem gives the following:
Proposition (ref) provides the individual level best linear unbiased prediction in a treatment effects setting. Two cases are of special interest.
\paragraph{Case 1: serial correlation only:} If $E[e_{\ell t}e_{js}]=0$ for $\ell\ne j$ and for all $t,s$, then for $i\in[1,N_1]$,
where $\theta_i =(E[e_{i,1:T_0} e_{i,1:T_0}'])^{-1} E[e_{i,1:T_0}e_{i,1:T_0+h}]$ is the $T_0\times 1$ vector of coefficients from projecting $e_{i,T_0+h} $ on $(e_{i,1},\ldots, e_{i,T_0})$.
The result in ((ref)) follows from the fact that when there is no cross-section dependence in the errors, then the correction for unit $i$ only depends on $\Omega_i$, the $T_0\times T_0$ autocovariance structure of $e_i$. For instance, if $N_1=1$ and $i=1$ is the treated unit, $$\omega_{ih}=
,$$ where $N_0=N-N_1$, and the prediction simplifies to $ m_{i,T_0+h}+\rho^h_{i} e_{GLS,i,T_0}$ if $e_{it}$ is assumed to be an AR(1), where $\rho_i$ is the first order autocorrelation coefficient of $e_{it}$. This coincides with Goldberger's correction reviewed in Section 3.
\paragraph{Case 2: cross-section correlation only:} If $E[e_{\ell t}e_{js}]=0$ for $t\ne s$ and for all $\ell,j$, then for $i \in [1,N_1]$,
where $\theta_{i,N_1+j}$ are the slope coefficients from projecting $e_{it}$ on $e_{N_1+1,t},\ldots,e_{N,t}$ using $t=1,\ldots,T_0$.
blup\; corrections with cross-section dependence have been derived in a static variance components setting but not in our set up. The result in ((ref)), which is new, is based on two features that follow from no serial correlation. First, for any treated unit $i\in [1,N_1]$, $$\omega_{ih}=
=
.$$ Second, the covariance matrix of errors $\Gamma$ has a block diagonal structure \[ \Gamma=
, \quad \Sigma =
,\] where $\Sigma_{11}=E(e_{1:N_1,t}e_{1:N_1,t}')$, $\Sigma_{00}=E(e_{N_1+1:N,t}e_{N_1+1:N,t}')$, and $\Sigma_{10}=E(e_{1:N_1,t}e_{N_1+1:N,t}')$. Let $[\Sigma_{10}]_{i,.}$ be the $i$-th row of the matrix $\Sigma_{10}$. Then, for any $i\in [1,N_1]$ and $h>0$, the non-zero entries of the $n\times 1$ vector $\omega_{ih}$ are given by \[ E(\mathcal E^{post}e_{i,T_0+h})= [\Sigma_{01}]_{i.}'\otimes J_h\] where $J_h$ is the $h$-th column of the identity matrix of dimension $T_1$. As a consequence of the two features, the blup\; correction is \[ \omega_{ih}'\Gamma^{-1} \mathcal E_{GLS}= ([\Sigma_{10}]_{i,.}\Sigma_{00}^{-1} \otimes J_{h}') \mathcal E^{post}_{GLS}. \] It is a linear combination of $\mathcal E_{GLS}^{post}$, the GLS errors of the control units in the post-treatment sample, with weights given by the $N-N_1$ vector $\theta_i'=(\theta_{i,N_1+1},\ldots, \theta_{i,N})=\Sigma_{00}^{-1}[\Sigma_{10}]_{i,.}'$, where $\theta_{i,N_1+j}$ is the population coefficient associated with $e_{N_1+j,t}$ in the regression of $e_{it}$ on $e_{N_1+1:N,t}=(e_{N_1+1,t},\ldots,e_{N,t})'$. Given the definition of $J_h$, we can write the blup\;correction for unit $i$ in period $T_0+h$ as $ \omega_{ih}'\Gamma^{-1} \mathcal E_{GLS}= \theta_i'e_{GLS,N_1+1:N,T_0+h}$, as given in the Proposition. The $\theta_i$ vector can be sparse such as in factor models when the idiosyncratic errors can only be weakly dependent in the sense that if $E(e_{it}e_{jt})=\tau_{ij,t}$, $|\tau_{ij,t}|\le |\bar \tau_{ij}|$ for some $\bar \tau_{ij}$ for all $t$, and $\sum_{j=1}^N |\bar \tau_{ij}|\le M\le \infty$ for all $i$.
In Section 3, we take as a starting point that the first order improvement of blup\;comes from controlling for the predictability in $e_{it}$. While blup\; is developed for linear predictions, linearity is not necessary to obtain an improved predictor. Our {\em practical unbiased predictor} (henceforth, pup) can be used with any choice of $\mathcal M$ that can be consistently estimated so that the residuals $\hat{\mathcal E}=\mathcal Y (0)-\hat{\mathcal M}$ can be used to improve prediction. In practice, implementation still requires parametric assumptions on $\Omega$ and $\Sigma$. For a serially correlated process $e_{it}$ satisfying mixing conditions, we have the following
Lemma (ref) follows immediately from Lemmas (ref) and (ref) but the results are presented in a treatment effect setting where the $o_p(1)$ term vanishes with $T_0$. For $h=1$, the MSE of a pup\;imputation will always be smaller than that of a standard imputation. We focus on $h=1$ since the gain in MSE should be largest, and furthermore, it is also the case when the direct and iterative corrections coincide. In theory, a direct pup\; using $\rho_{i,h}$ has better properties than an iterative pup\; using $\rho_{i,1}^h$ when $h>1$. But in simulations when sampling uncertainty is present, the two behave similarly and both dominate the standard predictor.
It is also possible to entertain both time and cross-section dependence in $e_{it}$. For example,
The first correction captures the time series information from $\hat e_{i}$'s own history, while the second correction captures cross-section dependence using estimates of the current and past idiosyncratic errors of the control units.
Though optimal prediction of counterfactual outcomes has not been studied, pup\;like corrections have recently been considered. cwz-jasa:21 consider time dependence in $e_{1t}$ and suggest adding to the standard imputed value an AR(p) estimate of the residuals, but the idea is not flushed out. fmm-jasa:22 consider a factor model with observables $W_{it}$ and assume contemporaneously correlated idiosyncratic errors $e_{it}=\theta_i' e_{-i,t}+ v_{it}$, where $e_{-i,t}$ is a $N-1$ vector that excludes $e_{it}$. This correlation can be controlled by $\hat e_{jt}$ from estimation of the factor model. To circumvent overfitting the augmented prediction model $ Y_{it}(0) =\gamma_i' W_{it}+\lambda_i'F_t+ \theta_1'\hat e_{-i,t}+v_{it}$ when $N_0$ is large, lasso is used to select which of the $\hat e_{j,t}$ to keep. Our analysis provides a framework for thinking about these pup\;like modifications.
In Table (ref), we saw in a linear prediction setting that plup\; yields significant reductions in bias and mean-squared error. Since imputation is a form of prediction, we can anticipate improvements from using pup\; in imputation settings. To illustrate this, we generate $Y_{it}(0)$ using a factor model:
We assume that unit 1 is treated with $\delta_{1t}=0.1$ and $e_{it}=\phi_i e_{i,t-1}+v_{it}$, $v_{it}\smash{\mathop{\sim}\limits^d} N(0,.25)$, $\phi_1=0.6$, $\phi_i=0$ for $i>1$. We set $c_i=0$ for all $i$, $r=2$, $F_{1t}=.8F_{1t-1}+e^F_{1t}$ and $F_{2t}=.5 F_{2,t-1}+e^F_{2t}$, $\Lambda_{ik}\smash{\mathop{\sim}\limits^d} N(0,1) $, $e^F_{1t}\smash{\mathop{\sim}\limits^d} N(0,.5)$, $e^F_{2t}\smash{\mathop{\sim}\limits^d} N(0,.3)$. The `best' prediction is $Y_{1,T_0+h}=\Lambda_1'F_{T_0+h}+ \phi_1^h e_{1,T_0}$ and the standard prediction is based on the principal components $\hat \Lambda_i'\hat F_{T_0+h}$ of the demeaned data for the control group, which will be denoted pca.
Table (ref) reports the error in imputing $\hat \delta_{1,T_0+h}$ for $(T,N)=(50,20)$. Because $Y_{1t}(1)$ is observed, $\hat\delta_{1t}-\delta_{1t}=Y_{1t}(1)-\hat Y_{1t}(0)$ is due entirely to $\hat Y_{1t}(0)$. The top panel considers $e_{1t}=0.6 e_{1t-1}+v_{1t}$. We see that though pca is unconditionally unbiased, its MSEs are larger compared to the pup\; ones. Furthermore, the conditional biases upon fixing $e_{1,T_0}=1.0$ are smaller with pup, as is MSE. For mutually correlated errors when $e_{1t}=0.5 e_{2t}+v_{1t}$, the conditional case randomly draws $e_{2,T_0+1:T_1}$ once and keeps the vector fixed in the replications. Again, the two imputations based on pup\; reduce MSE unconditionally and conditionally. Unlike in the time series case, the improvements can occur at any horizon $h$. Though direct pup\; has slightly better population properties than iterative pup\;, they have rather similar properties in simulations. To simplify notation, we use the term \textsc{ pup}\; in the next subsection, with the understanding that either direct or iterated correlation can be used.
Let $F\left( x\right) =P\left( e_{it}\leq x\right) $ be the marginal distribution function of $ e_{it}$, and $F^+\left( x\right) $ be the marginal distribution of an adjusted pup\; error. Lemma 3 suggests that
Though both $F$ and $F^+$ are centered at zero, $F$ is more dispersed than $F^+$. There is surprisingly little work on inference based on a feasible blup\; even in a single equation setting presumably because further assumptions on $F$ are needed. We will assume that $\left\{ e_{it}\right\} $ is a Gaussian process with autocovariance at lag $j$ of $\gamma_{j,i}$, and let $z_{1-\alpha }$ be such that $\Phi \left( z_{1-\alpha }\right) =1-\alpha$, where $\Phi$ is the cdf of the standard normal distribution. The Gaussian assumption is only made to illustrate the issues created by omitting predictability. Other distributions can be used in its place so long as $e_{it}$ satisfies our mixing assumption.
We consider intervals of the form $\hat\delta\pm \hat\sigma_\delta z_{1-\alpha/2}$ where $\hat\delta$ is either the standard predictor, or a pup\; predictor with variance defined in Lemma (ref). These intervals, denoted PI$_{ih}$, and PI$^{+}_{ih}$, can be based on asymptotic theory or resampling methods as in li-jasa:20 and cattaneo-jasa:21, among others. The intervals will be used for both unconditional and conditional inference. By unconditional inference, we mean that $e_{i,T_0}$ is random in repeated sampling. By conditional inference, we mean that $e_{i,T_0}$ is treated as fixed with respect to some conditioning information, as would be the case in practice. phillips-joe:79 notes that while unconditional inference is useful for evaluating econometric methods, conditional inference has a role in applications.
We begin with unconditional inference. We say that a prediction interval is unconditionally asymptotically valid if it contains $\delta _{i,T_{0}+h}$ with probability $1-\alpha $ as $ T_{0}\rightarrow \infty $. We will focus on pointwise results for unit $i$ where $ e_{i,T_0+h}\equiv Y_{i,T_0+h}( 0) -m_{i,T_0+h}$. It is easy to show that under the assumptions of Lemma (ref) and assuming that $\left\{ e_{it}\right\} $ is Gaussian, all three intervals are asymptotically valid unconditionally. Consider for instance $\text{\textsc{PI}}_{ih}$, an interval for the standard prediction error for unit $i$. Since $\hat{\delta}_{i,T_{0}+h}-\delta _{i,T_{0}+h}=e_{i,T_{0}+h}+o_{p}(1)$ and $e_{i,T_{0}+h}\smash{\mathop{\sim}\limits^d} N\left( 0,\gamma _{0,i}\right) $, where $\gamma_{0,i}\equiv\sigma^2_{e,i}$, we have
In the above, the third equality follows because $\hat{\sigma}_{e,i}^{-1}(\hat{\delta} _{i,T_{0}+h}-\delta _{i,T_{0}+h})=\sigma_{e,i} ^{-1}e_{i,T_{0}+h}+o_{p}(1)$ and the fourth equality uses the Gaussianity assumption on $e_{it}$. The argument for the pup\;prediction intervals is similar and thus all three intervals are (unconditionally) asymptotically valid. In the absence of sampling uncertainty, PI$^{+}$ is asymptotically narrower than PI for all $h$ whether or not $e_{it}$ is truly an AR(1) process.
While the three intervals provide correct unconditional coverage, will they all have correct conditional coverage? In particular, will they cover $\delta_{i,T_0+h}$ with the nominal coverage probability of $1-\alpha$, conditionally on some information set $\mathcal H$? To answer this question, consider again the AR(1) case when $e_{it}=\phi_i e_{i,t-1}+v_{it}$, $v_{it}\smash{\mathop{\sim}\limits^d} N(0,\sigma^2_{v,i})$, and $\mathcal H$ is an information set containing $e_{i,T_0}$. Under the assumption of normality, $e_{i,T_0+1} \smash{\mathop{\sim}\limits^d} N(\phi_i e_{i,T_0},\sigma^2_{v,i})$ conditionally on $e_{i,T_0}$ and hence $\frac{e_{i,T_0+1}-\phi_i e_{i,T_0}}{\sigma_{v,i}}\smash{\mathop{\sim}\limits^d} N(0,1)$. But as $ \frac{e_{i,T_0+1}}{\sigma_{e,i}} \smash{\mathop{\sim}\limits^d} N(\phi_i \frac{ e_{i,T_0}}{\sigma_{e,i}},\frac{\sigma^2_{v,i}} {\sigma^2_{e,i}}) \ne N(0,1)$, the standard prediction will not usually have the correct coverage unless $\phi_i$ or $e_{i,T_0}$ are zero. Indeed, PI$_{i1}$ will not have the correct conditional coverage probability even asymptotically because $P\left( \delta _{i,T_{0}+1}\in \text{\textsc{PI}}_{i1}\vert e_{i,T_0}\right)$ is asymptotically equal to
When $\phi _{i}\neq 0$, the first term will not return the normal cdf at level $1-\alpha/2$ unless ${e}_{i,T_0}$ is zero.
The problem of distorted inference extends to a multi-period ahead conditional inference. For any $h>1$, we see that
where $u_{i,T_{0}+h}\smash{\mathop{\sim}\limits^d} (0,\omega^2_{h,i})$, $\omega^2_{h,i}=\gamma _{0,i}\left( 1-\phi _{i}^{2h}\right)$,\footnote{This follows because $u_{i,T_{0}+h}=v_{i,T_{0}+h}+\ldots +\phi _{i}^{h-1}v_{i,T_{0}}$. Thus, $\omega^2_{h,i}=\sigma _{v,i}^{2}\sum_{j=0}^{h-1} \phi_i^{2h} =\sigma _{v,i}^{2}\frac{1-\phi _{i}^{2h}}{ 1-\phi _{i}^{2}}\equiv \gamma _{0,i}\left( 1-\phi _{i}^{2h}\right)$, $\gamma _{0,i} =\sigma _{v,i}^{2}\left( 1-\phi _{i}^{2}\right) ^{-1}$.} with $\gamma _{0,i}\equiv\sigma^2_{e,i} =\sigma _{v,i}^{2}\left( 1-\phi _{i}^{2}\right) ^{-1}$. The problem arises because the standard prediction error is not centered at zero when we condition on $e_{i,T_0}$. Thus for any $h\ge 1$, $P(\delta _{i,T_{0}+h}\in \text{\textsc{PI}}_{ih}|e_{i,T_0})$ is asymptotically equal to
where $\omega^2_{h,i}$ is the h-period forecast error variance defined above. For fixed $h$, conditional inference is distorted unless $e_{i,T_0}= 0$, though the distortion decreases with $h$ because $\phi_i^h$ tends to zero.
To illustrate the extent of conditional bias, consider two models for $e_{it}$: one when $e_{it}= \phi_i e_{i,t-1}+v_{it}$ is an AR(1), and one when $e_{it}=v_{it}+\theta_i v_{i,t-1}$ is an MA(1). In both cases, the conditional forecast is biased. Analytically evaluating actual coverage for any $i$ for a nominal $95\%$ interval, we have
The coverage probability is 0.95 at $\phi_1=\theta_1=0$ for all values of $h$ because there is no conditional bias. When $\phi_1\ne 0$ in the AR(1) case, bias decreases with $\sigma^2_{v,1}$ and $h$. When $\theta_1\ne 0$ in the MA(1) case, the bias is limited to $h=1$ because the process has a memory of one period. Coverage is distorted by bias, as suggested by theory.
In contrast, pup\; has error $\hat{\delta}^{+}_{i,T_{0}+h}-\delta _{i,T_{0}+h}=v_{i,T_{0}+h}+o_{p}(1)$. It is conditionally centered at zero with conditional variance $\sigma^2_{v,i}$. Let $\hat {\sigma} _{i,h}^{+}$ denote a consistent estimator of this variance (which will depend on whether an iterative or a direct estimator is used). Under normality, the pup\; prediction interval $\text{PI}_{ih}^{+}=\hat{\delta}_{i,T_{0}+h}^{+}\pm \hat{\sigma} _{i,h}^{+}z_{1-\alpha /2}$ has the correct coverage asymptotically. At $h=1$ when iterative and direct pup\; coincide, $P\left( \delta _{i,T_{0}+1}\in \text{\textsc{PI}}_{i1}^{+}|e_{i,T_0}\right)$ is asymptotically
where the last equality uses the fact that $\hat{\sigma}_{i,1}^{+}\smash{\mathop{\longrightarrow}\limits^p} \sigma _{v,i}$.
Analogous to conditional bias due to time dependence, a similar bias occurs when the errors are cross-sectionally correlated. In particular, suppose that $e_{it}$ is serially uncorrelated for all $i$, but $E[e_{it}e_{jt}]\ne 0$ for at least one $j\ne i$. If $i=1 $ is the only treated unit, the model of interest is $Y_{1t}(0)=m_{1t}+e_{1t}$ with $e_{1t}=\theta_1'e_{2:N,t}+v_{1t}$. Consider two treatment effects estimators, one based on the standard prediction of $Y_{1t}(0)$ and another based on pup\;. The standard prediction yields $\hat{\delta}_{1,T_0+1}-\delta_{1,T_0+1}=e_{1,T_0+1}+o_{p}(1)$ which has asymptotic mean zero and asymptotic variance $\sigma^2_{e,1}\equiv E(e^2_{1t})$. Instead, the pup\;prediction yields a prediction error $\hat{\delta}^{+}_{1,T_0+1}-\delta_{1,T_0+1}=e_{1,T_0+1}-\theta_1'e_{2:N,T_0+1}+o_{p}(1)$, whose asymptotic variance is $\sigma^2_{v,1}=\sigma^2_{e,1}-\Sigma_{01}'\Sigma_{00}^{-1}\Sigma_{01}<\sigma^2_{e,1}$. With cross-sectionally correlated errors, the standard prediction is conditionally biased because $e_{1,T_0+1}\vert e_{2:N,T_0+1}\smash{\mathop{\sim}\limits^d} N(\theta_1'e_{2:N,T_0+1},\sigma^2_{v,1})$ is not centered at zero. Conditional coverage under normality is asymptotically determined by $$\Phi\bigg(-\frac{\theta_1'}{\sigma_{v,1}} e_{2:N,T_0+1}+z_{1-\alpha /2}\frac{\sigma_{e,1}}{\sigma_{v,1}}\bigg)- \Phi\bigg(-\frac{\theta_1'}{\sigma_{v,1}} e_{2:N,T_0+1}-z_{1-\alpha /2}\frac{\sigma_{e,1}}{\sigma_{v,1}}\bigg).$$
The conditional bias arising from $\theta_1\ne 0$ and $\sigma_{e,1}\ne \sigma_{v,1}$ will distort inference. Notably, $-\frac{\theta_1'}{\sigma_{v,1}}$ here plays the role of $-\frac{\phi_1}{\sigma_{v,1}}$ in the AR(1) setting. But in contrast to the case of serial correlation, the size distortion does not diminish with $h$. The pup\; prediction is more efficient because $e_{1t}|e_{2:N,t}\smash{\mathop{\sim}\limits^d} (\theta_1'e_{2:N,t},\sigma^2_{v,1})$ which has a smaller variance.
To illustrate, suppose that $N_1=1$ and $e_{1t}=\theta_1 e_{2t}+v_{1t}$, where we assume that $\sigma_{12}$ is the only non-zero cross sectional covariance. In this case, the relevant parameters are $[\Sigma_{00}]_{11}=\sigma^2_{e,2}$, and $[\Sigma_{01}]_{11}=\text{cov}(e_{1t},e_{2t})\equiv \sigma_{12}$. Unlike in the time series case when we only condition on $e_{1,T_0}$, we now condition on $e_{2,T_0+h}$ for each $h$. In the following example, we draw $ e_{2,T_0+h}$ once from the normal distribution using the rndn function in matlab with seed 1234, and $\text{cov}(e_{1t},e_{2t})$ from the uniform distribution using the rand function with seed 57.
The parameterizations result in $\theta_1=\text{sgn} \times 0.72$ where the sign depends on whether $\Sigma_{01}\equiv \sigma_{12}$ is positive or negative. The sign affects the magnitude of the bias which in turn affects the extent of size distortion. However, unlike in the time series case when predictability falls with $h$ by nature of stationarity, the cross-section correlation does not decrease with $h$. As a consequence, the effect on coverage can vary significantly across horizons.
Tables (ref) and (ref) above showed that pup\; corrections reduce bias and mean-squared prediction error significantly. But do these improvements lead to more accurate inference? We first return to Table (ref) when point prediction is based on the simple linear model is $y_t=X_t'\beta+e_t$. The corresponding results for coverage are reported in Table (ref)(b). Note first that there is no gain in unconditional coverage regardless of the error structure because $E[e_{it}]=0$ by assumption. Hence we focus on conditional coverage. In the AR(1) case, coverage is 40% at $h=1$ without pup\;correction, but is at the desired level of 95% with correction. pup\; coverage in the AR(2) case is 91%, which is less accurate, but still better than the OLS coverage of 0.69.
Evaluation of coverage based on $\delta_{it}$ is more involved because the sampling error depends on the estimator. Furthermore, given estimates of $\hat\delta_{i,T_0+h}$ for some $i \in [1,N_1]$, any of the following hypothesis can be considered. \[
\] While $A$ and $B$ are pointwise hypotheses, hypotheses C,D,E concern the average treatment effect of the treated, where the average can be taken over time, over units, or both. The interpretation of each test depends on $T_1$ and $N_1$. Consider $H_0^C$. If $T_1$ is small so that $\Delta_{i,T_0+1:T}=\frac{1}{T_1}\sum_{t=T_0+1}^T \delta_{it}$ is random, we construct a {\em prediction} interval for $\Delta_{i,T_0+1:T}$. When $T_1$ is large, $\Delta_{i,T_0+1:T}\smash{\mathop{\longrightarrow}\limits^p} \Delta_{i,\infty}=E[\delta_{it}]$ is non-random. In this case, we construct a {\em confidence} interval for $\Delta_{1,T_0+1:T}$. Considerations of $N_1$ are likewise needed for testing $H_0^D$.\footnote{ For large $T_1$, $ \frac{1}{T_1}\sum_{h=1}^{T_1} m_{i,T_0+h}-\hat m_{i,T_0+h}+ \frac{1}{T_1}\sum_{h=1}^T e_{i,T_0+h}+ \delta_{i,T_0+h}-E[\delta_{i,T_0+h}] $, which equals $\hat {\Delta}_{i,T_0+1:T}-\Delta_{i,T_0+1:T}+ \delta_{i,T_0+h}-E[\delta_{i,T_0+h}].$ }
To give a flavor of the results, we only consider $H_0^A$ and $H_0^C$ using the same data generating process in Table (ref). We estimate $(F_t,\Lambda_{i})$ using the tall-wide procedure in baing-jasa:21. With this estimator, the asymptotic distribution of $\hat \delta_{i,T_0+h}$ is determined by the distribution of $e_{i,T_0+h}$. Hence the distribution of $\hat \delta_{i,T_0+h}$ is normal only if $e_{it}$ is Gaussian. However, the treatment effect averaged over $T_1$ periods can be asymptotically normal with a convergence rate of $\min(T_0,T_1)$ if the CLT $\frac{1}{\sqrt{T_0}} \sum_{t=1}^{T_0} F_te_{it}\smash{\mathop{\longrightarrow}\limits^d} N(0,\Phi_i)$ holds. In contrast, the average treatment effect over $N_1$ units can be asymptotically normal with a convergence rate of $\min(N_0,N_1)$ if a CLT for $\frac{1}{\sqrt{N_0}}\sum_{i=1}^{N_0}\Lambda_i e_{it}\smash{\mathop{\longrightarrow}\limits^d} N(0,\Gamma_t)$ holds. These distinctions will help understand the coverage results reported in Table (ref)(b).
Turning first to the time series case in the top panel of Table (ref)(b), pup\; coverage for $\delta_{i,T_0+h}$ improves over PCA, but only for $h=1,2$, suggesting that for serially correlated errors, pup\; is most effective when $h$ is small, especially when $T_0$ is small. For cross-sectionally correlated errors reported in the bottom panel, improved coverage of pup\; can occur at any $h$, but is more systematic when $N_0$ is large.
As explained above, the limiting distribution of the baing-jasa:21 estimator of $\hat{\bar\delta_i}$ depends on adequacy of central limit theory, while it is normality of $e_{it}$ that renders $\hat\delta_{it}$ normally distributed. It is thus not surprising that coverage of $\delta_{it}$ does not provide a good guide to the coverage of $\bar \delta_i=\frac{1}{T_1} \sum_{h=1}^{T_1} \delta_{i,T_0+h}$. At $(T_0,N_0)=(50,20)$, pup\; conditional coverage of $\bar\delta_i$ is too low, though no more inaccurate than the standard pca prediction. Whether the errors are serially or cross-sectionally correlated, coverage of $\bar\delta_i$ is more reliable as $T_0$ increases because we average over $T_1$ variables that become increasingly Gaussian. Sampling methods could be useful when $T_0$ is small, see e.g. li-shen-zhou:23.
goldberger:62 shows that if the errors of the prediction model are non-spherical, they can be exploited to improve prediction. Motivated by this result, this paper has suggested pup, a simple way to adjust existing estimates of counterfactual outcomes for dependence in the errors. The adjustment consists of adding a term that exploits the presence of serial or cross sectional correlation in the prediction errors of the model used to obtain the estimated counterfactual outcomes.
We showed that improved mean-squared errors are possible without knowledge of the true error structure, and simple corrections often suffice. We also showed that omitting the pup\; adjustment term when the error is predictable can result in conditional bias, thus leading to prediction intervals that are not conditionally asymptotically valid. In contrast, a prediction interval based on pup\; is conditionally unbiased, resulting in valid inference both conditionally and unconditionally.
While improved predictions are possible, it should be pointed out that when serial correlation is a concern, we can focus on corrections at small $h$ because $\rho^h$ will be small for stationary mixing processes. Furthermore, dependence in the residuals is necessary but not sufficient for improved prediction. This is because the pup\; adjustment depends not just on $\hat \rho_{is}$ or $\hat \rho_{js}$, but also on the values of $\hat e_{i,T_0-s}$ and $\hat e_{j,T_0+j-s}$ which can take on values close to zero. In the German unification example, the cross-section pup\; correction does not make much difference. However, the residuals exhibit serial dependence and the AR(1) pup\;correction changed the imputed growth rate for 1991 from -1.722 to -1.537 and for 1996 from -1.528 to -1.662. If take log GDP as the outcome variable instead of GDP growth, the residuals are still serially correlated with $\hat \rho_1=0.72$. But $(\hat e_{1,T_0},\hat e_{1,T_0-1})=(0.025,-0.0024)$ which are small relative to $y_{1,T_0}=9.26$. Hence in this case, the pup\; corrections did not make appreciable difference. Ultimately, whether the pup\; corrections are large is an empirical matter. Our goal is simply to draw awareness to the possibility of improvements. A practical first step could be to use the in-sample residuals $\hat e_{it}$ to construct an LM test for no dependence using the auxiliary regression \[ \hat e_{it} = X_{it}'\delta_0+\sum_{s=0}^{p_i} \delta_{i,s}\hat e_{i,t-s}+\sum_{j=1}^N \sum_{s=-p_j}^{p_j} \delta_{j,s} \hat e_{j,t-s}.\] We can entertain pup\; corrections if the null hypothesis is rejected, keeping in mind that the significance of the corrections depend not only on the hypothesis to be tested, but also on the estimator used.