# Imputation of Counterfactual Outcomes when the Errors are Predictable

Sílvia Gonçalves McGill University, Department of Economics, 855 Sherbooke St.
W., Montreal, PQ, H3A,2T7. Email<EMAIL_ADDRESS>Serena Ng
Columbia University and NBER. Department of Economics, 420 W. 118 St. MC 3308,
New York, NY 10027. Email<EMAIL_ADDRESS>
We thank Atsushi Inoue, Simon Lee, Jushan Bai, Kaspar Wüthrich, Marcelo
Medeiros, Bruno Ferman, Yuya Sasaki, Bruce Hansen, Jeff Wooldridge, Pedro
Sant’Anna, Barbara Rossi, and participants at the JBES session at the 2024
ASSA meeting and the Gary Chamberlain Seminar in Econometrics for helpful
comments. This work was supported by an NSERC grant and the National Science
Foundation (Ng SES 018369).

###### Abstract

A crucial input into causal inference is the imputed counterfactual outcome.
Imputation error can arise because of sampling uncertainty from estimating the
prediction model using the untreated observations, or from out-of-sample
information not captured by the model. While the literature has focused on
sampling uncertainty, it vanishes with the sample size. Often overlooked is
the possibility that the out-of-sample error can be informative about the
missing counterfactual outcome if it is mutually or serially correlated.
Motivated by the best linear unbiased predictor (blup) of Goldberger (1962) in
a time series setting, we propose an improved predictor of potential outcome
when the errors are correlated. The proposed pup is practical as it is not
restricted to linear models, can be used with consistent estimators already
developed, and improves mean-squared error for a large class of strong mixing
error processes. Ignoring predictability in the errors can distort conditional
inference. However, the precise impact will depend on the choice of estimator
as well as the realized values of the residuals.

Keywords: treatment effect, synthetic control, missing values, blup.

JEL classification: C22, C23

##  1 Introduction

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. In such problems, a researcher observes
the outcome of unit $i$ after intervention at $t>T_{0}$ (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 1 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))^{\prime}$. The
top-left panel displays actual (log) GDP over the full sample, along with in-
sample fit $\widehat{Y}_{1,1:T_{0}}(0)$, and the counterfactual values
$\widehat{Y}_{1,T_{0}+1:T}(0)$. The effect of reunification on GDP is stark,
but masks the fact that the (in-sample) residuals
$\widehat{e}_{1,1:T_{0}}=\Delta Y_{1,1:T_{0}}(0)-\widehat{\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 $\widehat{e}_{1,1:T_{0}}$ is also cross-correlated with
other errors, though many of the
$\text{corr}(\widehat{e}_{1,1:T_{0}},\widehat{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 2. 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 (1962) 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 first order equivalent to the
infeasible blup and show that 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 then 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 then show that 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.111We 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 conditionally biased which may
distort inference. 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,
$\widehat{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 information in 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 Chernozhukov, Wüthrich, and
Zhu (2021), Fan, Masini, and Medeiros (2022), and Ferman (2023) 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.

##  2 The Econometric Framework

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}=\begin{cases}Y_{it}(0),\quad&i=1,\ldots,N,\;\quad\quad\quad
t=1,\ldots,T_{0}\\\
Y_{it}(0),\quad&i=N_{1}+1,\ldots,N,\;\;\;t=T_{0}+1,\ldots,T\\\
Y_{it}(1)&i=1,\ldots,N_{1},\;\quad\quad\;\;t=T_{0}+1,\ldots,T\end{cases}$

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)}_{\begin{subarray}{c}\text{observed
outcome }\\\ \text{h periods after treatment
}\end{subarray}}-\underbrace{Y_{i,T_{0}+h}(0)}_{\begin{subarray}{c}\text{
outcome without treatment}\\\ \text{in period }T_{0}+h\end{subarray}}.$

The econometrics challenge is that $Y_{it}(0)$ is not observed for $i\leq
N_{1}$ when $t>T_{0}$.

Following the literature, we assume that $Y_{it}(0)$ has a pseudo-true
conditional mean $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$. 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}$,

$\displaystyle Y_{it}(0)$ $\displaystyle=$ $\displaystyle
m_{it}+e_{it},\quad\quad\quad\quad t=1,\ldots,T$ $\displaystyle Y_{it}(1)$
$\displaystyle=$ $\displaystyle m_{it}+\delta_{it}+e_{it},\quad\quad t>T_{0}.$

Let $\widehat{m}_{it}$ be a consistent estimate of $m_{it}$. Then

$Y_{it}(0)=\widehat{m}_{it}+\widehat{e}_{it}.$

Since $\widehat{Y}_{it}(0)=\widehat{m}_{it}$, the treatment effect on unit $i$
at a given $t$ is then estimated by

$\displaystyle\widehat{\delta}_{i,T_{0}+h}$ $\displaystyle=$ $\displaystyle
Y_{i,T_{0}+h}(1)-\widehat{Y}_{i,T_{0}+h}(0)$ $\displaystyle=$ $\displaystyle
Y_{i,T_{0}+h}(1)-Y_{i,T_{0}+h}(0)+(Y_{i,T_{0}+h}(0)-\widehat{m}_{i,T_{0}+h})$
$\displaystyle=$
$\displaystyle\delta_{i,T_{0}+h}+e_{i,T_{0}+h}+(m_{i,T_{0}+h}-\widehat{m}_{i,T_{0}+h}).$

The pointwise imputation/prediction error is

$\displaystyle\widehat{\delta}_{i,T_{0}+h}-\delta_{i,T_{0}+h}$
$\displaystyle=$
$\displaystyle(m_{i,T_{0}+h}-\widehat{m}_{i,T_{0}+h})+e_{i,T_{0}+h}.$

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 $\widehat{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 $e_{it}$ is
unpredictable. 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.

###  2.1 Examples when $e_{it}$ is predictable

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}^{\prime}\beta$. Single equation estimation yields the
imputed value
$\widehat{\delta}_{1,T_{0}+1}=x_{T_{0}+1}^{\prime}\widehat{\beta}$ and
imputation error
$\widehat{\delta}_{1,T_{0}+1}-\delta_{1,T_{0}+1}=-x_{T_{0}+1}^{\prime}(\widehat{\beta}-\beta)+e_{1,T_{0}+1}$
whose variance is

$\text{var}(\widehat{\delta}_{1,T_{0}+1}-\delta_{1,T_{0}+1})=\sigma_{e}^{2}+x_{T_{0}+1}^{\prime}\text{var}(\widehat{\beta})x_{T_{0}+1}.$

Correlation in $x_{t}e_{t}$ may necessitate robust standard errors for
$\widehat{\beta}$, but provided that $E[x_{t}e_{1t}]=0$, $\widehat{\beta}$ is
consistent in the sense that $\text{var}(\widehat{\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.

#### 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}}\neq\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}\neq 0$.

The underlying premise of many imputation procedures is that $m_{it}$ has a
low rank structure. The method of Difference-in-Difference (DID) specifies
$m_{it}=c_{i}+\gamma_{t}$ and is suited for $N_{1}>1$ units. Factor-based
imputation assumes $X_{it}=\lambda_{i}^{\prime}F_{t}+\epsilon_{it}$ where
$F_{t}$ is a vector of $r$ latent common factors with $\lambda_{i}$ as
loadings and $\epsilon_{it}$ is an idiosyncratic error. DID obtains as a
special case with $F_{t}=(1,\gamma_{t})^{\prime}$ and
$\lambda_{i}=(c_{i},1)^{\prime}$.

#### Example misspecification 2:

Suppose that the potential outcome has an interactive fixed effect structure:
$Y_{1t}(0)=\lambda_{1}^{\prime}F_{t}+\epsilon_{1t}$ but a researcher specifies
an additive fixed effect model $Y_{it}(0)=\lambda_{i}+F_{t}+\epsilon_{it}$.
Then $e_{1t}$ will be serially correlated if $F_{t}$ is serially correlated,
even if $\epsilon_{1t}$ is white noise.

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 (2017)
iteratively estimates $F$ and the missing values jointly by principal
components (PCA), Bai and Ng (2021) impute the missing values (or complete the
matrix) using two full-sample applications of PCA. Athey, Bayati, Doudchenko,
Imbens, and Khosravi (2021); Arkhangelsky, Athey, Hirshberg, Imbens, and Wager
(2021) estimate the low rank component using singular value thresholding
(SVT).222Regularization is not necessary to consistently estimate the missing
values, but could give a lower rank common component than the one in Bai and
Ng (2021). Though all consistent estimators of $F$ imply that
$\widehat{F}_{T_{0}+h}$ can be used as though they were observed predictors,
the possibility remains that $\epsilon_{1,T_{0}+h}$ can be predicted by
information available.

#### Example Correlated Idiosyncratic Errors

(from Fan, Masini, and Medeiros (2022)) Suppose that
$Y_{it}(0)=\Lambda_{i}^{\prime}F_{t}+\epsilon_{it}$ and the researcher
correctly assumes $m_{1t}=\Lambda_{1}^{\prime}F_{t}$ but $\epsilon_{it}$ is
correlated with $\epsilon_{jt}$ for $j$ in some index set $C$. Then
$e_{1t}=\epsilon_{1t}+\sum_{j\in C}\theta_{j}\epsilon_{jt}$ is predictable by
those $\epsilon_{jt}$ where $j\in C$.

The method of synthetic control (SC) developed in Abadie and Gardeazabal
(2003) 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\leq
T_{0}$. Abadie, Diamond, and Hainmueller (2010) 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 Arkhangelsky, Athey, Hirshberg, Imbens, and Wager
(2021) also reweights the pre-treatment time periods to balance 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^{*}$.

#### Example Imperfect Fit 1:

(from Ben-Michael, Feller, and Rothstein (2021)) 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.

#### Example Imperfect Fit 2:

(from Ferman and Pinto (2021)) Suppose that
$Y_{it}(0)=c_{i}+\Lambda_{i}^{\prime}F_{t}+\epsilon_{it}$ and one estimates
$\widehat{\beta}=\text{argmin}_{b}\|Y_{1}(0)-X_{0}b)||_{2}^{2}$ where
$X_{0}=(Y_{2},\ldots,Y_{N})$. With $\beta^{*}=\text{plim}\widehat{\beta}$, the
population imputation error is

$e_{1t}=Y_{1t}(0)-X_{0t}^{\prime}\beta^{*}=\bigg{(}c_{1}-\sum_{j=2}^{N}\beta^{*}_{j}c_{j}\bigg{)}+F_{t}^{\prime}\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})\neq(\sum_{j=2}^{N}\beta_{j}^{*}c_{j},\sum_{j=2}^{N}\beta_{j}^{*}\lambda_{j})$.
Ferman and Pinto (2021) 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.

##  3 Prediction with 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.

###  3.1 Goldberger’s BLUP

This subsection summarizes results for best linear unbiased prediction. The
concept seems to be first considered in Henderson (1950) 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.333Robinson (1991) provides a survey of its many derivations, including a
Kalman filter interpretation, see also Spall (1991). Taub (1979) and Baltagi
(2008, 2013) use it in variance components analysis of panel data. Goldberger
(1962) formalizes the idea in a setting where the linear model for predicting
a scalar variable $y_{t}$ is given by

$y_{t}=X_{t}^{\prime}\beta+e_{t}$ (1)

where $X_{t}$ is a $K\times 1$ vector of completely observed predictors
assumed to be fixed in repeated samples, $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 $e$.

Goldberger (1962) 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. 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

$\displaystyle\sigma_{m}^{2}$ $\displaystyle=$ $\displaystyle
E[(y_{m|T_{0}}-y_{m})^{2}]=E[A^{\prime}ee^{\prime}A+e_{m}^{2}-2A^{\prime}e_{m}e]$
$\displaystyle=$ $\displaystyle A^{\prime}\Omega
A+E\left(e_{m}^{2}\right)-2A^{\prime}\omega$

where $\omega\equiv E[e_{m}e]$. 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)

$y_{m|T_{0}}^{\ast}=x_{m}^{\prime}\beta_{GLS}+\omega^{\prime}\Omega^{-1}{e}_{GLS},$

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}\sim(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

$y_{T_{0}+h|T_{0}}^{\ast}=X_{T_{0}+h}^{\prime}\beta_{GLS}+\phi_{1}^{h}{e}_{GLS,T_{0}}$

with prediction error
$e^{*}_{T_{0}+1|T_{0}}=y_{T_{0}+1}-y_{T_{0}+1|T_{0}}^{\ast}=e_{T_{0}+1}-\phi_{1}e_{T_{0}}+o_{p}(1)=v_{T_{0}+1}+o_{p}(1)$.

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.444Cochrane-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}\widehat{\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 and
Orcutt (1949) suggest to improve the standard prediction by incorporating lags
of regressors and the dependent variable. The difference is that
$\widehat{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
$\widehat{\beta}_{GLS}$ instead of $\widehat{\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

$\displaystyle\widehat{y}_{T_{0}+1|T_{0}}^{\ast}$ $\displaystyle=$
$\displaystyle
X_{T_{0}+1}^{\prime}\widehat{\beta}_{GLS}+\widehat{\phi_{1}}\widehat{e}_{GLS,T_{0}}$

has prediction error
$\widehat{e}_{T_{0}+1|T_{0}}^{*}=y_{T_{0}+1}-\widehat{y}_{T_{0}+1|T_{0}}^{*}$,
or

$\displaystyle\widehat{e}_{T_{0}+1|T_{0}}^{\ast}$ $\displaystyle=$
$\displaystyle
v_{T_{0}+1}-\left(X_{T_{0}+1}-\phi_{1}X_{T_{0}}\right)^{\prime}(\widehat{\beta}_{GLS}-\beta)-\left(\widehat{\phi}_{1,GLS}-\phi_{1}\right)\widehat{e}_{GLS,T_{0}}$
(2) $\displaystyle=$ $\displaystyle v_{T_{0}+1}+o_{p}(1)$

where the $o_{p}(1)$ term comes from the fact that the jointly estimated
$\widehat{\beta}_{GLS}$ and $\widehat{\phi}_{1,GLS}$ are $\sqrt{T_{0}}$
consistent for $\beta$ and $\phi_{1}$. Since $\widehat{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

$\displaystyle\widehat{e}_{T_{0}+1|T_{0}}$ $\displaystyle=$ $\displaystyle
e_{T_{0}+1}-X_{T_{0}+1}^{\prime}(\widehat{\beta}_{OLS}-\beta)$ (3)
$\displaystyle=$ $\displaystyle e_{T_{0}+1}+o_{p}(1)$

where the $o_{p}(1)$ term comes from $\sqrt{T_{0}}$ consistency of
$\widehat{\beta}_{OLS}$ for $\beta$. But $\widehat{e}_{T_{0}+1|T_{0}}$ is
asymptotically $e_{T_{0}+1}$ whose variance is
$\sigma^{2}_{e}\geq\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}$.

###  3.2 From BLUP to PLUP

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
$\widehat{\beta}$ denote the least squares estimate of $\beta$. Consider
modifying the (standard) least-squares prediction
$\widehat{y}_{T_{0}+1|T_{0}}=X_{T_{0}+1}^{\prime}\widehat{\beta}$ as follows:

$\widehat{y}^{+}_{T_{0}+1|T_{0}}=\widehat{y}_{T_{0}+1|T_{0}}+\widehat{\rho}_{1}\widehat{e}_{T_{0}}$
(4)

where $\widehat{e}_{T_{0}}=y_{T_{0}}-X_{T_{0}}^{\prime}\widehat{\beta}$ is the
OLS residual, and

$\widehat{\rho_{1}}=\frac{\sum_{t=1}^{T_{0}}\widehat{e}_{t-1}\widehat{e}_{t}}{\sum_{t=1}^{T_{0}}\widehat{e}_{t-1}^{2}}\\\
$ (5)

is the least squares estimate of the first order autocorrelation coefficient
of $\widehat{e}_{t}$. Note that unlike feasible blup which re-estimates
$\beta$ after $\widehat{\rho}_{1}$ is available, we simply adjust the OLS
prediction $\widehat{y}_{T_{0}+1|T_{0}}$ for serial correlation with
$\widehat{\rho}_{1}\widehat{e}_{T_{0}}$. The prediction error
$\widehat{e}_{T_{0}+1|T_{0}}^{+}=y_{T_{0}+1}-\widehat{y}^{+}_{T_{0}+1|T_{0}}=y_{T_{0}+1}-\widehat{y}_{T_{0}+1|T_{0}}-\widehat{\rho}_{1}\widehat{e}_{T_{0}}$
is

$\displaystyle\widehat{e}_{T_{0}+1|T_{0}}^{+}$ $\displaystyle=$
$\displaystyle\widehat{e}_{T_{0}+1|T_{0}}-\widehat{\rho}_{1}\widehat{e}_{T_{0}}$
$\displaystyle=$ $\displaystyle e_{T_{0}+1}-\rho_{1}e_{T_{0}}+o_{p}(1),$

where the last equality follows because
$\widehat{\beta}\smash{\mathop{\longrightarrow}\limits^{p}}\beta$ and
$\widehat{\rho}_{1}\smash{\mathop{\longrightarrow}\limits^{p}}\rho_{1}$. If
$e_{t}$ is indeed an AR(1) model, then $\rho_{1}=\phi_{1}$ and

$\displaystyle\widehat{e}_{T_{0}+1|T_{0}}^{+}$ $\displaystyle=$ $\displaystyle
v_{T_{0}+1}+o_{p}(1),$

which is asymptotically equal to $v_{T_{0}+1}$, the prediction error of the
infeasible blup.

As it turns out, adding the term $\widehat{\rho}_{1}\widehat{e}_{T_{0}}$ to
any consistent prediction $\widehat{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

$\widehat{y}_{T_{0}+1|T_{0}}^{+}=X_{T_{0}+1}^{\prime}\widehat{\beta}+\widehat{\rho}_{1}\widehat{e}_{T_{0}}$

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
$\widehat{e}_{T_{0}+1|T_{0}}^{+}=y_{T_{0}+1}-\widehat{y}_{T_{0}+1|T_{0}}^{+}$,
we assume the following.

Assumption A1

(a)

$E\left|e_{t}\right|^{r}<\infty$ for some $r>2$, for all $t$.

(b)

$\left\\{e_{t}\right\\}$ is a zero mean strictly stationary strong mixing
process with mixing coefficients
$\alpha\left(k\right)=O\left(k^{-\frac{r}{r-2}-\delta}\right)$ for some
$\delta>0$.

We define the strong mixing coefficients as
$\alpha\left(k\right)=\sup_{A,B}\left|P\left(A\cap
B\right)-P\left(A\right)P\left(B\right)\right|$ 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|\psi_{j}\right|<\infty$ and $v_{t}$ is
i.i.d.$\left(0,\sigma_{v}^{2}\right)$ with $E\left|v_{t}\right|^{r}<\infty$,
which includes stationary invertible ARMA(p,q) processes and nonlinear weakly
dependent processes with GARCH and ARCH innovations.

###### Lemma 1

($h=1$): Let $X_{t}$ be predictors and $e_{t}$ be the errors in the model
defined by (1). Suppose that $\\{e_{t}\\}$ satisfies Assumption A1 and that
for $j=0,1$, $E\left(X_{t}e_{t-j}\right)$, $E\left(X_{t-j}e_{t}\right)$ and
$E\left(X_{t}X_{t-j}^{\prime}\right)$ exist such that
$\widehat{\beta}\smash{\mathop{\longrightarrow}\limits^{p}}\beta$. Then as
$T_{0}\rightarrow\infty$,

  * (i)

$\widehat{\rho}_{1}\smash{\mathop{\longrightarrow}\limits^{p}}\rho_{1}\equiv\frac{\gamma_{1}}{\gamma_{0}}$,
where $\gamma_{k}\equiv E\left(e_{t}e_{t-k}\right)$ for all $k$;

  * (ii)

Standard prediction error: $\widehat{e}_{T_{0}+1|T_{0}}=e_{T_{0}+1}+o_{p}(1)$
where $e_{T_{0}+1}\sim(0,\gamma_{0})$;

  * (iii)

plup error:
$\widehat{e}_{T_{0}+1|T_{0}}^{+}=e_{T_{0}+1}-\rho_{1}e_{T_{0}}+o_{p}(1)$ where
$e_{T_{0}+1}-\rho_{1}e_{T_{0}}\sim(0,\gamma_{0}(1-\rho_{1}^{2}))$.

Part (i) shows that $\widehat{\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}$ has mean zero, 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|\rho_{1}\right|\leq 1$,

$\gamma_{0}(1-\rho_{1}^{2})\leq\gamma_{0},$

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
$\widehat{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
$\widehat{\rho}_{1}\smash{\mathop{\longrightarrow}\limits^{p}}\rho_{1}$ as
stated in (i). The plup error is now
$\widehat{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
$\widehat{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|\rho_{1}\right|\leq
1$. Thus the AR(1) correction unambiguously reduces one-step asymptotic mean
squared prediction error.

Next, consider cases when $h>1$. The standard prediction
$\widehat{y}_{T_{0}+h|T_{0}}=X_{T_{0}+h}^{\prime}\widehat{\beta}$ has error
$\widehat{e}_{T_{0}+h|T_{0}}=y_{T_{0}+h}-\widehat{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}$. Iterated plup (or plupi), defined as

$\widehat{y}_{T_{0}+h|T_{0}}^{+I}=X_{T_{0}+h}^{\prime}\widehat{\beta}+\widehat{\rho}_{1}^{h}\widehat{e}_{T_{0}},$

has error
$\widehat{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}}\sim(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 $\widehat{\rho}_{h}=\widehat{\gamma}_{h}/\widehat{\gamma}_{0}$
be the $h^{th}$ order sample autocorrelation coefficient of
$\\{\widehat{e}_{t}:t=1,\ldots,T_{0}\\}$. Direct plup (or plupd) is defined as

$\widehat{y}_{T_{0}+h|T_{0}}^{+d}=X_{T_{0}+h}^{\prime}\widehat{\beta}+\widehat{\rho}_{h}\widehat{e}_{T_{0}}.$

plupd has error
$\widehat{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}}\sim(0,\gamma_{0}\left(1-\rho_{h}^{2}\right)).$

###### Lemma 2

($h\geq 1)$ Under the same assumptions as in Lemma 1, the asymptotic MSE of
plupd is always smaller than or equal to that of plupi and that of the
standard predictor for all $h\geq 1$.

The proof is given in the Appendix. The standard predictor, plupi and plupd
are all asymptotically unbiased provided that $E[e_{t}]=0$. But

$\displaystyle\text{var}(e_{T_{0}+h}-\rho_{h}e_{T_{0}})-\text{var}\left(e_{T_{0}+h}-\rho_{1}^{h}e_{T_{0}}\right)$
$\displaystyle=$
$\displaystyle-\gamma_{0}\left(\rho_{h}-\rho_{1}^{h}\right)^{2}\leq 0$

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}\leq
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\geq 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.

###  3.3 Simulations for Linear Predictions

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}^{\prime}\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}\sim 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 A 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 $\widehat{y}_{T_{0}+h}$ using the least squares estimate
$\widehat{\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 A 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 A 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.

##  4 Imputation of Counterfactual Outcomes

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 and Rubin (2019, Ch. 11) 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 and Lin (1971) 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 and Scanlan (2024)
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
uncorrelated555Brodersen, Gallusser, Koehler, Remy, and Scott (2015) 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, and Medeiros (2018); Masini and Medeiros (2021,
2022) 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:

#### Assumptions A2:

For $h\geq 1$, $\widehat{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(\widehat{m}_{it}-m_{it}\right)^{2}=o_{p}\left(1\right)$.

Assumption A2 is verified in Chernozhukov, Wüthrich, and Zhu (2021) for
estimators including synthetic control, matrix completion, factor-based
methods. Given an asymptotically unbiased $\widehat{m}_{i,T_{0}+h}$ satisfying
Assumptions A1 and A2, the estimated treatment effect

$\displaystyle\widehat{\delta}_{i,T_{0}+h}$ $\displaystyle=$ $\displaystyle
m_{i,T_{0}+h}+\delta_{i,T_{0}+h}+e_{i,T_{0}+h}-\widehat{m}_{i,T_{0}+h}$

has error

$\displaystyle\widehat{\delta}_{i,T_{0}+h}-\delta_{i,T_{0}+h}$
$\displaystyle=$
$\displaystyle(m_{i,T_{0}+h}-\widehat{m}_{i,T_{0}+h})+e_{i,T_{0}+h}.$ (6)

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}=\begin{pmatrix}[l]Y_{1,1:T_{0}}\\\
\vdots\\\ Y_{N,1:T_{0}}\\\ Y_{N_{1}+1,T_{0}+1:T}\\\ \vdots\\\
Y_{N,T_{0}+1:T}\end{pmatrix}\equiv\begin{pmatrix}[l]\mathcal{Y}^{pre}(0)_{(NT_{0}\times
1)}\\\ \\\ \mathcal{Y}^{post}(0)_{(N-N_{1})T_{1}\times 1}\end{pmatrix}$

where $n=(NT-N_{1}T_{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}})^{\prime}$ 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})^{\prime}$.

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,

$\displaystyle\mathcal{Y}(0)$ $\displaystyle=$
$\displaystyle\mathcal{M}+\mathcal{E},\quad\quad\mathcal{E}=\begin{pmatrix}\mathcal{E}^{pre}\\\
\mathcal{E}^{post}\end{pmatrix}\sim(0,\Gamma)$ $\displaystyle\Gamma$
$\displaystyle=$ $\displaystyle
E[\mathcal{E}\mathcal{E}^{\prime}]=\begin{pmatrix}\Gamma_{pre,pre}&\Gamma_{pre,post}\\\
\Gamma_{post,pre}&\Gamma_{post,post}\end{pmatrix}.$

The $n\times n$ matrix $\Gamma$ depends on $\Sigma=E[e_{t}e_{t}^{\prime}]$,
which is the $N\times N$ covariance matrix of
$e_{t}=(e_{1:N_{1},t}^{\prime},e_{N_{1}+1:N,t}^{\prime})^{\prime}$. It also
depends on $\Omega_{i}=E[e_{i}e_{i}^{\prime}]$, 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 when $\mathcal{M}=\mathcal{X}\beta$, where
$\mathcal{X}$ contains observed predictors and consider obtaining the blup
given information on $\mathcal{Y}(0)$ and $\mathcal{X}$. Re-doing Goldberger’s
problem gives the following:

###### Proposition 1

Assume that $\mathcal{Y}(0)=\mathcal{M}+\mathcal{E}$, where
$\mathcal{M}=\mathcal{X}\beta$, and let $\Gamma$ be the $n\times n$ variance-
covariance of $\mathcal{E}$, where $n=NT-N_{1}T_{1}$. The blup of the
counterfactual outcome for unit $i\in[1,N_{1}]$ is

$\mathcal{Y}_{i,T_{0}+h}^{+}(0)=X_{T_{0}+h}^{\prime}\beta_{GLS}+\omega_{ih}^{\prime}\Gamma^{-1}\mathcal{E}_{GLS}$

where $\beta_{GLS}$ is the vector of infeasible GLS estimates and
$\mathcal{E}_{GLS}$ are the corresponding residuals,
$\omega_{ih}=E[\mathcal{E}e_{i,T_{0}+h}]$ is $n\times 1$ vector of covariances
between the unexplained errors in the vector of observed control outcomes and
unit $i$’s counterfactual outcome not explained by the model at $T_{0}+h$.

Proposition 1 provides the individual level best linear unbiased prediction in
a treatment effects setting. Two cases are of special interest.

#### Case 1: serial correlation only:

If $E[e_{\ell t}e_{js}]=0$ for $\ell\neq j$ and for all $t,s$, then for
$i\in[1,N_{1}]$,

$\omega_{ih}^{\prime}\Gamma^{-1}\mathcal{E}_{GLS}=\theta_{i}^{\prime}e_{GLS,i,1:T_{0}}\
$ (7)

where
$\theta_{i}=(E[e_{i,1:T_{0}}e_{i,1:T_{0}}^{\prime}])^{-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 (7) 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}=\begin{pmatrix}E[e_{i,1:T_{0}}e_{i,T_{0}+h}]\\\ 0_{N_{0}T\times
1}\end{pmatrix},$

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.

#### Case 2: cross-section correlation only:

If $E[e_{\ell t}e_{js}]=0$ for $t\neq s$ and for all $\ell,j$, then for
$i\in[1,N_{1}]$,

$\omega_{ih}^{\prime}\Gamma^{-1}\mathcal{E}_{GLS}=\sum_{j=1}^{N-N_{1}}\theta_{i,N_{1}+j}e_{GLS,N_{1}+j,T_{0}+h}.$
(8)

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 (8), 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}=\begin{pmatrix}E(\mathcal{E}^{pre}e_{i,T_{0}+h})\\\
E(\mathcal{E}^{post}e_{i,T_{0}+h})\end{pmatrix}=\begin{pmatrix}0_{NT_{0}\times
1}\\\ E(\mathcal{E}^{post}e_{i,T_{0}+h})\end{pmatrix}.$

Second, the covariance matrix of errors $\Gamma$ has a block diagonal
structure

$\Gamma=\begin{pmatrix}\Sigma\otimes I_{T_{0}}&0\\\ 0&\Sigma_{00}\otimes
I_{T_{1}}\end{pmatrix},\quad\Sigma=\begin{pmatrix}\Sigma_{11}&\Sigma_{10}\\\
\Sigma_{10}^{\prime}&\Sigma_{00}\end{pmatrix},$

where $\Sigma_{11}=E(e_{1:N_{1},t}e_{1:N_{1},t}^{\prime})$,
$\Sigma_{00}=E(e_{N_{1}+1:N,t}e_{N_{1}+1:N,t}^{\prime})$, and
$\Sigma_{10}=E(e_{1:N_{1},t}e_{N_{1}+1:N,t}^{\prime})$. 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.}^{\prime}\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}^{\prime}\Gamma^{-1}\mathcal{E}_{GLS}=([\Sigma_{10}]_{i,.}\Sigma_{00}^{-1}\otimes
J_{h}^{\prime})\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}^{\prime}=(\theta_{i,N_{1}+1},\ldots,\theta_{i,N})=\Sigma_{00}^{-1}[\Sigma_{10}]_{i,.}^{\prime}$,
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})^{\prime}$. Given the
definition of $J_{h}$, we can write the blup correction for unit $i$ in period
$T_{0}+h$ as
$\omega_{ih}^{\prime}\Gamma^{-1}\mathcal{E}_{GLS}=\theta_{i}^{\prime}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}|\leq|\bar{\tau}_{ij}|$ for some $\bar{\tau}_{ij}$ for all $t$,
and $\sum_{j=1}^{N}|\bar{\tau}_{ij}|\leq M\leq\infty$ for all $i$.

###  4.1 From blup to pup

As 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 practical unbiased predictor (henceforth, pup) can be
used with any choice of $\mathcal{M}$ that can be consistently estimated so
that the residuals
$\widehat{\mathcal{E}}=\mathcal{Y}(0)-\widehat{\mathcal{M}}$ can be used to
control predictability. 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 3

Under Assumptions A1 and A2, a standard imputation has error
$\widehat{\delta}_{i,T_{0}+1}-\delta_{i,T_{0}+1}=e_{i,T_{0}+1}+o_{p}(1)$,
whose variance is $\sigma^{2}_{e,i,1}=\gamma_{0,i}$. A pup correction has
error
$\widehat{\delta}_{i,T_{0}+1}^{+}-\delta_{i,T_{0}+1}=e_{i,T_{0}+1}-\rho_{i,1}e_{i,T_{0}}+o_{p}(1)$,
whose variance is $(\sigma^{+}_{e,i,1})^{2}=\gamma_{0,i}(1-\rho^{2}_{i,1})$,
where $\rho_{i,1}$ is the first-order autocorrelation of $e_{it}$. Since
$|\rho_{i,1}|\leq 1$, $(\sigma^{+}_{e,i,1})^{2}\leq\sigma^{2}_{e,i,1}$.

Lemma 3 follows immediately from Lemmas 1 and 2 but the results are presented
in a treatment effect setting. 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,

$\widehat{Y}_{i,T_{0}+h}^{+}(0)=\widehat{Y}_{i,T_{0}+h}(0)+\sum_{s=0}^{p_{i}}\widehat{\rho}_{is}\widehat{e}_{i,T_{0}-s}+\sum_{j\neq
i}^{N}\sum_{s=-p_{j}}^{p_{j}}\widehat{\theta}_{ij,s}\widehat{e}_{j,T_{0}+h-s}.$
(9)

The first correction captures the time series information from
$\widehat{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. Chernozhukov, Wüthrich, and
Zhu (2021) 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. Fan, Masini, and Medeiros (2022) consider a factor model with
observables $W_{it}$ and assume contemporaneously correlated idiosyncratic
errors $\epsilon_{it}=\theta_{i}^{\prime}\epsilon_{-i,t}+v_{it}$, where
$\epsilon_{-i,t}$ is a $N-1$ vector that excludes $\epsilon_{it}$. This
correlation can be controlled by $\widehat{\epsilon}_{jt}$ from estimation of
the factor model. To circumvent overfitting the augmented prediction model
$Y_{it}(0)=\gamma_{i}^{\prime}W_{it}+\lambda_{i}^{\prime}F_{t}+\theta_{1}^{\prime}\widehat{\epsilon}_{-i,t}+v_{it}$
when $N_{0}$ is large, lasso is used to select which of the
$\widehat{\epsilon}_{j,t}$ to keep. Our analysis provides a framework for
thinking about these pup like modifications.

###  4.2 Simulations: $\widehat{\delta}_{it}$ vs $\widehat{\delta}_{it}^{+}$

In Table A, 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:

$\displaystyle Y_{it}(0)$ $\displaystyle=$ $\displaystyle
c_{i}+\Lambda_{i}^{\prime}F_{t}+e_{it},\quad t=1,\ldots,T$ $\displaystyle
Y_{1t}(1)$ $\displaystyle=$ $\displaystyle
Y_{1t}(0)+\delta_{1t},\quad\quad\delta_{1t}=0.1\quad t=T_{0}+1,\ldots,T.$

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}\sim 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}=.5F_{2,t-1}+e^{F}_{2t}$,
$\Lambda_{ik}\sim N(0,1)$, $e^{F}_{1t}\sim N(0,.5)$, $e^{F}_{2t}\sim N(0,.3)$.
The ‘best’ prediction is
$Y_{1,T_{0}+h}=\Lambda_{1}^{\prime}F_{T_{0}+h}+\phi_{1}^{h}e_{1,T_{0}}$ and
the standard prediction is based on the principal components
$\widehat{\Lambda}_{i}^{\prime}\widehat{F}_{T_{0}+h}$ of the demeaned data for
the control group, which will be denoted pca.

Table A reports the error in imputing $\widehat{\delta}_{1,T_{0}+h}$ for
$(T,N)=(50,20)$. Because $Y_{1t}(1)$ is observed,
$\widehat{\delta}_{1t}-\delta_{1t}=Y_{1t}(1)-\widehat{Y}_{1t}(0)$ is due
entirely to $\widehat{Y}_{1t}(0)$. The top panel considers
$e_{1t}=0.6e_{1t-1}+\epsilon_{1t}$ but all errors are mutually uncorrelated.
We see that though pca is unconditionally unbiased, its MSEs are larger
compared to the pup ones. In the serial correlation case, 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.5e_{2t}+\epsilon_{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 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 pup in the
next subsection, with the understanding that either direct or iterated
correlation can be used.

##  5 Inference Based on Imputed Counterfactual Outcomes

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

$\displaystyle\widehat{\delta}_{i,T_{0}+1}-\delta_{i,T_{0}+1}$
$\displaystyle\overset{d}{\longrightarrow}$ $\displaystyle e_{i,T_{0}+1}\sim
F$ $\displaystyle\widehat{\delta}_{i,T_{0}+1}^{+}-\delta_{i,T_{0}+1}$
$\displaystyle\overset{d}{\longrightarrow}$ $\displaystyle
e_{i,T_{0}+1}-\rho_{1}e_{i,T_{0}}\sim F^{+}.$

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
$\widehat{\delta}\pm\widehat{\sigma}_{\delta}z_{1-\alpha/2}$ where
$\widehat{\delta}$ is either the standard predictor, or a pup predictor with
variance defined in Lemma 3. These intervals, denoted PIih, and
PI${}^{+}_{ih}$, can be based on asymptotic theory or resampling methods as in
Li (2020) and Cattaneo, Feng, and Titiunik (2021), 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 (1979) 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 3 and assuming that $\left\\{e_{it}\right\\}$ is Gaussian, all three
intervals are asymptotically valid unconditionally. Consider for instance
$\text{{PI}}_{ih}$, an interval for the standard prediction error for unit
$i$. Since
$\widehat{\delta}_{i,T_{0}+h}-\delta_{i,T_{0}+h}=e_{i,T_{0}+h}+o_{p}(1)$ and
$e_{i,T_{0}+h}\sim N\left(0,\gamma_{0,i}\right)$, where
$\gamma_{0,i}\equiv\sigma^{2}_{e,i}$, we have

$\displaystyle P\left(\delta_{i,T_{0}+h}\in\text{{PI}}_{ih}\right)$
$\displaystyle=$ $\displaystyle
P\left(\widehat{\delta}_{i,T_{0}+h}-\widehat{\sigma}_{e,i}z_{1-\alpha/2}\leq\delta_{i,T_{0}+h}\leq\widehat{\delta}_{i,T_{0}+h}+\widehat{\sigma}_{e,i}z_{1-\alpha/2}\right)$
$\displaystyle=$ $\displaystyle
P\left(-z_{1-\alpha/2}\leq\widehat{\sigma}_{e,i}^{-1}(\widehat{\delta}_{i,T_{0}+h}-\delta_{i,T_{0}+h})\leq
z_{1-\alpha/2}\right)$ $\displaystyle=$ $\displaystyle
P\left(-z_{1-\alpha/2}\leq\sigma_{e,i}^{-1}e_{i,T_{0}+h}\leq
z_{1-\alpha/2}\right)+o\left(1\right)$ $\displaystyle=$
$\displaystyle\Phi\left(z_{1-\alpha/2}\right)-\Phi\left(-z_{1-\alpha/2}\right)+o\left(1\right)$
$\displaystyle=$ $\displaystyle 1-\alpha+o\left(1\right).$

In the above, the third equality follows because
$\widehat{\sigma}_{e,i}^{-1}(\widehat{\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.

###  5.1 Conditional Inference

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}\sim 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}\sim 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}}\sim N(0,1)$. But as
$\frac{e_{i,T_{0}+1}}{\sigma_{e,i}}\sim
N(\phi_{i}\frac{e_{i,T_{0}}}{\sigma_{e,i}},\frac{\sigma^{2}_{v,i}}{\sigma^{2}_{e,i}})\neq
N(0,1)$, the standard prediction will not usually have the correct coverage
unless $\phi_{i}$ or $e_{i,T_{0}}$ are zero. Indeed, PIi1 will not have the
correct conditional coverage probability even asymptotically because
$P\left(\delta_{i,T_{0}+1}\in\text{{PI}}_{i1}|e_{i,T_{0}}\right)$ is
asymptotically equal to

$\Phi\left(-\frac{\phi_{i}}{\sigma_{v,i}}{e}_{i,T_{0}}+\frac{\sigma_{e,i}}{\sigma_{v,i}}z_{1-\alpha/2}\right)-\Phi\left(-\frac{\phi_{i}}{\sigma_{v,i}}{e}_{i,T_{0}}-\frac{\sigma_{e,i}}{\sigma_{v,i}}z_{1-\alpha/2}\right).$

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

$e_{i,T_{0}+h}=\phi_{i}^{h}e_{i,T_{0}}+u_{i,T_{0}+h}|e_{i,T_{0}}\sim
N\left(\phi_{i}^{h}e_{i,T_{0}},\omega_{h,i}^{2}\right),$

where $u_{i,T_{0}+h}\sim(0,\omega^{2}_{h,i})$,
$\omega^{2}_{h,i}=\gamma_{0,i}\left(1-\phi_{i}^{2h}\right)$,666This 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\geq 1$,
$P(\delta_{i,T_{0}+h}\in\text{{PI}}_{ih}|e_{i,T_{0}})$ is asymptotically equal
to

$\Phi\left(-\frac{\phi_{i}^{h}}{\omega_{h,i}}e_{i,T_{0}}+z_{1-\alpha/2}\frac{\sigma_{e,i}}{\omega_{h,i}}\right)-\Phi\left(-\frac{\phi_{i}^{h}}{\omega_{h,i}}e_{i,T_{0}}-z_{1-\alpha/2}\frac{\sigma_{e,i}}{\omega_{h,i}}\right)$

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 it 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, bias is
$-\frac{\phi_{i}^{h}}{\omega_{h,i}}e_{i,T_{0}}$. Analytically evaluating
actual coverage for any $i$ for a nominal $95\%$ interval, we have

Standard Prediction with $(e_{1,T_{0}},\sigma_{v,1})=(-2.0,0.5)$

h | coverage | bias | coverage | bias
---|---|---|---|---
| AR(1): $\phi_{1}=0.8$ | MA(1): $\theta_{1}=0.8$
1 | 0.84 | -2.26 | 0.58 | -1.77
2 | 0.87 | -1.41 | 0.95 | -0.00
3 | 0.90 | -1.01 | 0.95 | -0.00
4 | 0.92 | -0.76 | 0.95 | -0.00
5 | 0.93 | -0.59 | 0.95 | -0.00

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}\neq 0$ in the AR(1)
case, bias decreases with $\sigma^{2}_{v,1}$ and $h$. When $\theta_{1}\neq 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
$\widehat{\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 $\widehat{\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}^{+}=\widehat{\delta}_{i,T_{0}+h}^{+}\pm\widehat{\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{{PI}}_{i1}^{+}|e_{i,T_{0}}\right)$ is
asymptotically

$\Phi\left(\frac{\widehat{\sigma}_{i,1}^{+}}{\sigma_{v,i}}z_{1-\alpha/2}\right)-\Phi\left(-\frac{\widehat{\sigma}_{i,1}^{+}}{\sigma_{v,i}}z_{1-\alpha/2}\right)=1-\alpha+o_{p}\left(1\right),$

where the last equality uses the fact that
$\widehat{\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}]\neq 0$ for
at least one $j\neq 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}^{\prime}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
$\widehat{\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
$\widehat{\delta}^{+}_{1,T_{0}+1}-\delta_{1,T_{0}+1}=e_{1,T_{0}+1}-\theta_{1}^{\prime}e_{2:N,T_{0}+1}+o_{p}(1)$,
whose asymptotic variance is
$\sigma^{2}_{v,1}=\sigma^{2}_{e,1}-\Sigma_{01}^{\prime}\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}|e_{2:N,T_{0}+1}\sim
N(\theta_{1}^{\prime}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}^{\prime}}{\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}^{\prime}}{\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}\neq 0$ and
$\sigma_{e,1}\neq\sigma_{v,1}$ will distort inference. Notably,
$-\frac{\theta_{1}^{\prime}}{\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}\sim(\theta_{1}^{\prime}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.

Standard Prediction when
$e_{1,T_{0}+h}=\theta_{1}e_{2,T_{0}+h}+v_{1,T_{0}+h}$, $e_{2,T_{0}+h}$ fixed

$\Sigma=\text{diag}(\sigma^{2}_{e,1},\Sigma_{00})=\text{diag}(0.5,0.841)$.

|  | $\Sigma_{01}=0.613$ | $\Sigma_{01}=-0.613$
---|---|---|---
h | $\bar{e}_{2,T_{0}+h}$ | coverage | bias | coverage | bias
1 | 0.68 | 0.63 | -1.64 | 0.89 | -0.70
2 | -0.83 | 0.81 | 1.07 | 0.86 | 0.85
3 | -0.92 | 0.93 | -0.38 | 0.84 | 0.95
4 | 0.09 | 0.95 | 0.17 | 0.95 | -0.09
5 | 0.86 | 0.94 | 0.32 | 0.86 | -0.89

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.

###  5.2 Simulations: PI vs PI+

Tables A and A 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 A when point prediction is based
on the simple linear model is $y_{t}=X_{t}^{\prime}\beta+e_{t}$. The
corresponding results for coverage are reported in Table A(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
$\widehat{\delta}_{i,T_{0}+h}$ for some $i\in[1,N_{1}]$, any of the following
hypothesis can be considered.

$\begin{array}[]{lllll}H_{0}^{A}&:\delta_{i,T_{0}+h}&=\delta_{i,T_{0}+h}^{0}\\\
H_{0}^{B}&:\delta_{i,T_{0}+h}&=0&\forall h\\\
H_{0}^{C}&:\Delta_{i,T_{0}+1:T}&=\Delta_{i,T_{0}+1:T}^{0}&\Delta_{i,T_{0}+1:T}&=\frac{1}{T_{1}}\sum_{h=1}^{T_{1}}\delta_{i,T_{0}+h}\\\
H_{0}^{D}&:\Delta_{1:N_{1},T_{0}+h}&=\Delta_{1:N_{1},T_{0}+h}^{0},&\Delta_{1:N_{1},T_{0}+h}&=\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}\delta_{i,T_{0}+h}\\\
H_{0}^{E}&:\Delta_{1:N_{1},T_{0}+1:T}&=\Delta_{1:N_{1},T_{0}+1:T_{1}}^{0},&\Delta_{1:N_{1},T_{0}+1:T}&=\frac{1}{T_{1}}\frac{1}{N_{1}}\sum_{h=1}^{T_{1}}\sum_{i=1}^{N_{1}}\delta_{i,T_{0}+h}\end{array}$

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 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 confidence interval for
$\Delta_{1,T_{0}+1:T}$. Considerations of $N_{1}$ are likewise needed for
testing $H_{0}^{D}$.777 For large $T_{1}$,
$\frac{1}{T_{1}}\sum_{h=1}^{T_{1}}m_{i,T_{0}+h}-\widehat{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
$\widehat{\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 A. We estimate
$(F_{t},\Lambda_{i})$ using the tall-wide procedure in Bai and Ng (2021). With
this estimator, the asymptotic distribution of $\widehat{\delta}_{i,T_{0}+h}$
is determined by the distribution of $e_{i,T_{0}+h}$. Hence the distribution
of $\widehat{\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_{t}e_{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 A(b).

Turning first to the time series case in the top panel of Table A(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 Bai and Ng (2021)
estimator of $\widehat{\bar{\delta}_{i}}$ depends on adequacy of central limit
theory, while it is normality of $e_{it}$ that renders $\widehat{\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, and Zhou (2023).

##  6 Conclusion

Goldberger (1962) 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
$\widehat{\rho}_{is}$ or $\widehat{\rho}_{js}$, but also on the values of
$\widehat{e}_{i,T_{0}-s}$ and $\widehat{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 $\widehat{\rho}_{1}=0.72$. But
$(\widehat{e}_{1,T_{0}},\widehat{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 $\widehat{e}_{it}$ to construct an LM test for no dependence
using the auxiliary regression

$\widehat{e}_{it}=X_{it}^{\prime}\delta_{0}+\sum_{s=0}^{p_{i}}\delta_{i,s}\widehat{e}_{i,t-s}+\sum_{j=1}^{N}\sum_{s=-p_{j}}^{p_{j}}\delta_{j,s}\widehat{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.

##  Appendix A Appendix: Proofs

Proof of Lemma 1. The proof of part (i) follows from standard arguments, see
e.g., Hayashi (2000), p. 145. We provide a brief sketch here. Write
$\widehat{\rho_{1}}=\frac{\widehat{\gamma}_{1}}{\widehat{\gamma}_{0}}$, with
$\widehat{\gamma}_{j}=\frac{1}{T_{0}}\sum_{t=1}^{T_{0}}\widehat{e}_{t}\widehat{e}_{t-j}$,
for $j=0,1$. Let
$\tilde{\gamma}_{j}=\frac{1}{T_{0}}\sum_{t=1}^{T_{0}}e_{t}e_{t-j}$ and note
that $\tilde{\gamma}_{j}\rightarrow_{p}\gamma_{j}$ under Assumption A1. Since
$\widehat{e}_{t}=y_{t}-X_{t}^{\prime}\widehat{\beta}=e_{t}+X_{t}^{\prime}(\beta-\widehat{\beta})$,
we can write

$\widehat{\gamma}_{j}=\tilde{\gamma}_{j}-T_{0}^{-1}\sum_{t=1}^{T_{0}}(X_{t-j}e_{t}+X_{t}e_{t-j})^{\prime}(\widehat{\beta}-\beta)+(\widehat{\beta}-\beta)^{\prime}(T_{0}^{-1}\sum_{t=1}^{T_{0}}X_{t}X_{t-j}^{\prime})^{-1}(\widehat{\beta}-\beta).$

The last two terms are $o_{p}\left(1\right)$ under the assumption that
$\widehat{\beta}-\beta=o_{p}\left(1\right)$ and assuming that
$E(X_{t-j}e_{t})$, $E(X_{t}e_{t-j}^{\prime})$ and
$E\left(X_{t}X_{t-j}^{\prime}\right)$ are all finite for $j=0,1$. It follows
that
$\widehat{\gamma}_{j}=\tilde{\gamma}_{j}+o_{p}(1)=\gamma_{j}+o_{p}\left(1\right)$,
which implies the result. To prove (ii), note that
$\widehat{e}_{T_{0}+1|T_{0}}=y_{T_{0}+1}-X_{T_{0}+1}^{\prime}\widehat{\beta}=e_{T_{0}+1}-X_{T_{0}+1}^{\prime}(\widehat{\beta}-\beta)=e_{T_{0}+1}+o_{p}\left(1\right)$
under the assumption that $\widehat{\beta}-\beta=o_{p}(1)$. It follows
immediately that $E(e_{T_{0}+1})=0$ and
$Var(e_{T_{0}+1})\equiv\gamma_{0}=\sigma_{e}^{2}$. For (iii), note that we can
write
$\widehat{e}_{T_{0}+1|T_{0}}^{+}=e_{T_{0}+1}-\rho_{1}e_{T_{0}}+o_{p}(1)$,
given (i) and the assumption that $\widehat{\beta}-\beta=o_{p}(1)$. Since
$E\left(e_{t}\right)=0$ for all $t$, we have $E(e_{T_{0}+1}-\rho_{1}e_{T})=0$.
In addition,
$Var(e_{T_{0}+1}-\rho_{1}e_{T})=\gamma_{0}+\rho_{1}^{2}\gamma_{0}-2\rho_{1}\gamma_{1}=\gamma_{0}-\gamma_{1}^{2}/\gamma_{0}\leq\gamma_{0}$,
where the second equality follows by replacing
$\rho_{1}=\gamma_{1}/\gamma_{0}$.

Proof of Lemma 2. The standard OLS prediction error is

$\widehat{e}_{T_{0}+h|T_{0}}=y_{T_{0}+h}-\widehat{y}_{T_{0}+h|T_{0}}=e_{T_{0}+h}-X_{T_{0}+h}^{\prime}(\widehat{\beta}-\beta)=e_{T_{0}+h}+o_{p}\left(1\right)\text{.}$

The mean of $e_{T_{0}+h}$ is zero and the variance is $\gamma_{0}$. The
$h$-step ahead direct plup based on an AR(1) approximation is

$\widehat{e}_{T_{0}+h|T_{0}}^{+d}=y_{T_{0}+h}-\widehat{y}_{T_{0}+h|T_{0}}^{+d}=e_{T_{0}+h}-\widehat{\rho}_{h}\widehat{e}_{T_{0}}+o_{p}\left(1\right)=e_{T_{0}+h}-\rho_{h}e_{T_{0}}+o_{p}\left(1\right)\text{,}$

where $\widehat{e}_{T_{0}}$ is the OLS residual for observation $T_{0}$, and
$\widehat{\rho}_{h}$ is the $h$-th order sample autocorrelation of
$\widehat{e}_{t}$. It follows immediately that the asymptotic MSE of the
direct plup is smaller or equal than that of OLS since the variance of
$e_{T_{0}+h}-\rho_{h}e_{T_{0}}$ is $\gamma_{0}(1-\rho_{h}^{2})\leq\gamma_{0}$,
since $\rho_{h}^{2}\leq 1$. To compare direct plup with iterated plup , note
that the $h$-step ahead iterated plup based on an AR(1) approximation is

$\widehat{y}_{T_{0}+h|T_{0}}^{+I}=\widehat{y}_{T_{0}+h|T_{0}}+\widehat{\rho}_{1}^{h}\widehat{e}_{T_{0}},$

where $\widehat{\rho}_{1}$ is defined in (5) and $\widehat{e}_{T_{0}}$ is the
OLS residual for observation $T_{0}$. The iterated plup error is

$\widehat{e}_{T_{0}+h|T_{0}}^{+I}=y_{T_{0}+h}-\widehat{y}_{T_{0}+h|T_{0}}^{+I}=e_{T_{0}+h}-\widehat{\rho}_{1}^{h}\widehat{e}_{T_{0}}+o_{p}\left(1\right)=e_{T_{0}+h}-\rho_{1}^{h}e_{T_{0}}+o_{p}\left(1\right)\text{.}$

Both versions of plup have mean zero but the variance of iterated plup is
larger or equal than that of direct plup since we can show that the difference
between these two variances is $\gamma_{0}(\rho_{h}-\rho_{1}^{h})^{2}\geq 0$.
A comparison between iterated plup and the standard prediction reveals that
both prediction errors have mean zero, but the variance of the standard OLS
prediction is $Var\left(e_{T_{0}+h}\right)=\gamma_{0}$ whereas that of the
iterated plup error is equal to

$\displaystyle Var\left(e_{T_{0}+h}-\rho_{1}^{h}e_{T_{0}}\right)$
$\displaystyle=$
$\displaystyle\gamma_{0}+\rho_{1}^{2h}\gamma_{0}-2\rho_{1}^{h}\gamma_{h}$
$\displaystyle=$
$\displaystyle\gamma_{0}[1+\rho_{1}^{2h}-2\rho_{1}^{h}(\frac{\gamma_{h}}{\gamma_{0}})]$
$\displaystyle=$
$\displaystyle\gamma_{0}[1+\rho_{1}^{2h}-2\rho_{1}^{h}\rho_{h})],$

since $\rho_{h}=\frac{\gamma_{h}}{\gamma_{1}}$, where
$\gamma_{h}=E\left(e_{t+h}e_{t}\right)$. When
$1+\rho_{1}^{2h}-2\rho_{1}^{h}\rho_{h}\leq 1$ we can conclude that the mean
square prediction error of plup is smaller than the mean square prediction
error of OLS.

Figure 1: Effect of German Unification, GDP growth

The data are downloaded from https://doi.org/10.7910/DVN/24714

Figure 2: Effect of German Unification on level of GDP

Table 1: Prediction Errors

$\displaystyle y_{t}$ $\displaystyle=$ $\displaystyle
X_{1t}\beta_{1}+X_{2t}\beta_{2}+e_{t},\quad\beta_{1}=\beta_{2}=1,\quad T=50$
$\displaystyle X_{1t}$ $\displaystyle=$ $\displaystyle
0.6X_{1,t-1}+u_{1t},\quad u_{1t}\sim N(0,1),\quad\quad X_{2t}=u_{2t},\quad
u_{2t}\sim N(0,1)$ $\displaystyle e_{t}$ $\displaystyle=$
$\displaystyle\phi_{1}e_{t-1}+\phi_{2}e_{t-2}+v_{t},\quad v_{t}\sim N(0,0.05)$

Unconditional

h | best | noadj | ols | plupI | plupD | best | noadj | ols | plupI | plupD
---|---|---|---|---|---|---|---|---|---|---
| bias | mse
|  ($\phi_{1},\phi_{2}$)=(0.80,0.00)
1 | -0.00 | -0.01 | -0.00 | -0.00 | -0.00 | 0.05 | 0.14 | 0.14 | 0.05 | 0.05
2 | 0.00 | -0.00 | 0.00 | 0.00 | 0.00 | 0.05 | 0.14 | 0.14 | 0.08 | 0.08
5 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.05 | 0.13 | 0.14 | 0.12 | 0.12
10 | 0.00 | -0.00 | -0.00 | -0.00 | 0.00 | 0.05 | 0.14 | 0.14 | 0.14 | 0.14
avg | -0.00 | -0.00 | -0.00 | -0.00 | -0.00 | 0.00 | 0.07 | 0.08 | 0.06 | 0.06
|  ($\phi_{1},\phi_{2}$)=(1.30,-0.40)
1 | -0.00 | -0.01 | -0.00 | -0.00 | -0.00 | 0.05 | 0.43 | 0.43 | 0.06 | 0.06
2 | 0.00 | -0.00 | -0.00 | -0.00 | -0.00 | 0.05 | 0.43 | 0.43 | 0.16 | 0.16
5 | 0.00 | -0.00 | 0.00 | -0.00 | 0.00 | 0.05 | 0.42 | 0.43 | 0.37 | 0.35
10 | 0.00 | -0.01 | -0.00 | -0.00 | -0.00 | 0.05 | 0.43 | 0.44 | 0.48 | 0.44
avg | -0.00 | -0.00 | -0.00 | -0.00 | -0.00 | 0.00 | 0.28 | 0.29 | 0.21 | 0.20

Conditional on $e_{T_{0}}=1.0$ and $e_{T_{0}-1}=0.5$

bias mse ($\phi_{1},\phi_{2}$)=(0.80,0.00) 1 -0.00 0.80 0.79 0.01 0.01 0.05
0.68 0.68 0.05 0.05 2 0.00 0.64 0.64 0.01 0.02 0.05 0.49 0.49 0.09 0.09 5 0.00
0.33 0.33 0.01 0.02 0.05 0.23 0.23 0.13 0.13 10 0.00 0.11 0.10 0.00 0.02 0.05
0.15 0.15 0.14 0.15 avg -0.00 0.36 0.35 0.01 0.02 0.00 0.18 0.18 0.06 0.06
($\phi_{1},\phi_{2}$)=(1.30,-0.40) 1 -0.00 1.10 1.09 0.18 0.18 0.05 1.25 1.26
0.09 0.09 2 0.00 1.03 1.03 0.18 0.24 0.05 1.19 1.20 0.17 0.20 5 0.00 0.62 0.62
-0.04 0.20 0.05 0.72 0.73 0.34 0.38 10 0.00 0.21 0.21 -0.23 0.08 0.05 0.47
0.48 0.49 0.45 avg -0.00 0.61 0.61 -0.04 0.17 0.00 0.56 0.56 0.19 0.22

5000 replications using Matlab R2019a, seed = rng(1234,’twister’)

]

& T=200

Unconditional Conditional ($\phi_{1},\phi_{2}$)=(0.80,0.00) 1 0.96 0.95 0.95
0.95 0.95 0.96 0.39 0.40 0.95 0.95 2 0.95 0.94 0.94 0.94 0.94 0.96 0.63 0.63
0.94 0.93 5 0.95 0.95 0.94 0.94 0.94 0.96 0.87 0.86 0.93 0.92 10 0.95 0.94
0.93 0.93 0.92 0.96 0.93 0.92 0.94 0.91 avg 0.94 0.99 0.98 0.92 0.92 0.95 0.93
0.91 0.91 0.89 ($\phi_{1},\phi_{2}$)=(1.30,-0.40) 1 0.96 0.95 0.94 0.95 0.95
0.96 0.73 0.69 0.91 0.91 2 0.95 0.94 0.94 0.94 0.94 0.96 0.72 0.70 0.93 0.90 5
0.95 0.94 0.93 0.94 0.93 0.96 0.85 0.84 0.94 0.91 10 0.95 0.93 0.92 0.93 0.91
0.96 0.93 0.91 0.93 0.91 avg 0.94 0.98 0.97 0.94 0.91 0.95 0.92 0.90 0.94 0.89

Table 2: Errors in Estimated Treatment Effect : $T=50,N=20$

& $e_{1t}=0.6e_{1t-1}+\epsilon_{1t}$, $e_{T_{0}}$=1.00 Conditional bias of
$\widehat{\delta}_{i,T_{0}+h}$ mse 1 0.00 0.60 0.56 0.15 0.15 0.25 0.61 0.66
0.36 0.36 2 -0.01 0.35 0.32 0.12 0.14 0.33 0.45 0.53 0.45 0.47 5 -0.00 0.21
0.19 0.09 0.13 0.37 0.42 0.51 0.48 0.52 10 -0.01 0.12 0.09 0.04 0.09 0.39 0.40
0.51 0.50 0.54 avg -0.00 0.15 0.12 0.04 0.09 0.12 0.14 0.19 0.17 0.19

$e_{1t}=0.5e_{2t}+\epsilon_{1t}$ Unconditonal bias of
$\widehat{\delta}_{i,T_{0}+h}$ mse 1 0.00 -0.00 0.00 0.00 0.00 0.15 0.20 0.33
0.28 0.28 2 -0.01 -0.01 -0.00 -0.01 -0.01 0.15 0.21 0.35 0.29 0.29 5 0.00
-0.00 0.01 0.01 0.01 0.15 0.21 0.35 0.30 0.30 10 -0.01 -0.01 -0.01 0.00 0.00
0.15 0.21 0.34 0.29 0.29 avg -0.00 -0.00 0.00 0.01 0.01 0.01 0.02 0.05 0.04
0.04 $e_{1t}=0.5e_{2t}+\epsilon_{1t}$, $e_{2T_{0}+1}=0.27539$ Conditional bias
of $\widehat{\delta}_{i,T_{0}+h}$ mse 1 0.00 0.14 0.14 0.04 0.04 0.15 0.17
0.29 0.27 0.27 2 -0.01 0.49 0.49 0.12 0.12 0.15 0.38 0.52 0.35 0.35 5 0.00
0.47 0.48 0.13 0.13 0.15 0.37 0.52 0.35 0.35 10 -0.01 -0.45 -0.44 -0.14 -0.14
0.15 0.35 0.48 0.34 0.34 avg -0.00 -0.12 -0.12 -0.06 -0.06 0.01 0.03 0.06 0.05
0.05

5000 replications using Matlab 2019a with seed rng(456,’twister’). For cross-
section correlation, $e_{2,T_{0}:1:T}$=(0.6361, -0.6386, 0.7118, -1.7044,
-1.2992, 1.6402, 0.1395, 0.9348, 0.5051, 0.9692).

& ($T_{0},N_{0}$)=(200,50): Unconditional Conditional 1 0.95 0.95 0.94 0.95
0.95 0.96 0.90 0.89 0.95 0.95 2 0.91 0.95 0.94 0.94 0.94 0.91 0.92 0.92 0.94
0.94 5 0.89 0.95 0.94 0.93 0.93 0.89 0.95 0.94 0.94 0.93 10 0.89 0.95 0.94
0.94 0.93 0.90 0.95 0.94 0.94 0.93 avg 0.99 0.93 0.92 0.92 0.92 0.99 0.92 0.90
0.92 0.91

$e_{1t}=\theta_{1}e_{2t}+v_{1t}$

($T_{0},N_{0}$)=(50,20): Unconditional Conditional 1 0.95 0.95 0.94 0.93 0.93
0.95 0.97 0.95 0.93 0.93 2 0.95 0.95 0.93 0.93 0.93 0.95 0.86 0.87 0.90 0.90 5
0.95 0.94 0.93 0.91 0.91 0.95 0.91 0.91 0.90 0.90 10 0.95 0.95 0.93 0.88 0.88
0.95 0.98 0.96 0.89 0.89 avg 0.88 0.88 0.79 0.78 0.78 0.88 0.81 0.75 0.77 0.77

($T_{0},N_{0}$)=(200,50): Unconditional Conditional 1 0.95 0.96 0.95 0.95 0.95
0.95 0.94 0.94 0.95 0.95 2 0.95 0.95 0.94 0.94 0.94 0.95 0.94 0.93 0.94 0.94 5
0.95 0.95 0.95 0.95 0.95 0.95 0.75 0.78 0.93 0.93 10 0.95 0.95 0.95 0.94 0.94
0.95 0.86 0.87 0.93 0.93 avg 0.94 0.94 0.92 0.92 0.92 0.94 0.92 0.91 0.92 0.92

## References

  * (1)
  * Abadie, Diamond, and Hainmueller (2010) Abadie, A., A. Diamond, and J. Hainmueller (2010): “Synthetic Control Methods for Comparative Case Studies: Estimating the Effect of California’s Tobacco Control Program,” _Journal of the American Statistical Association_ , 105:490, 493–505, DOI: 10.1198/jasa.2009.ap08746.
  * Abadie and Gardeazabal (2003) Abadie, A., and J. Gardeazabal (2003): “The Economic Costs of Conflict: A Case Study of the Basque Country,” _American Economic Review_ , 93:1.
  * Arkhangelsky, Athey, Hirshberg, Imbens, and Wager (2021) Arkhangelsky, D., S. Athey, D. A. Hirshberg, G. W. Imbens, and S. Wager (2021): “Syntheticx Difference-in-Differences,” _American Economic Review_ , 111-(12), 4088–4118, DOI: 10.1257/aer.20190159.
  * Athey, Bayati, Doudchenko, Imbens, and Khosravi (2021) Athey, S., M. Bayati, N. Doudchenko, G. Imbens, and K. Khosravi (2021): “Matrix Completion Methods for Causal Panel Data Models,” _Journal of the American Statistical Association_ , 116:536, 1716–1730, DOI: 10.1080/01621459.2021.1891924.
  * Bai and Ng (2021) Bai, J., and S. Ng (2021): “Matrix Completion, Counterfactuals, and Factor Analysis of Missing Data,” _Journal of the American Statistical Association_ , 116:536, 1746–1763.
  * Baltagi (2008) Baltagi, B. (2008): “Forecasting with Panel Data,” _Journal of Forecasting_ , 27, 153–173.
  * Baltagi (2013) (2013): “Panel Data Forecasting,” _Handbook of Forecasting_ , 2B, 995–2021.
  * Ben-Michael, Feller, and Rothstein (2021) Ben-Michael, E., A. Feller, and J. Rothstein (2021): “The Augmented Synthetic Control Method,” _Journal of the American Statistical Associaton_ , 116:536, 1789–1803.
  * Brodersen, Gallusser, Koehler, Remy, and Scott (2015) Brodersen, K., F. Gallusser, J. Koehler, N. Remy, and S. Scott (2015): “Inferring Causal Impact Using Bayesian Structural Time-Series Models,” _Annals of Applied Statistics_ , 9:1, 247–274.
  * Carvalho, Masini, and Medeiros (2018) Carvalho, C., R. Masini, and M. Medeiros (2018): “ArCo: An Artificial Counterfactual Approach for High-Diemsnional Panel Time Series Data,” _Journal of Econometrics_ , 207, 352–380.
  * Cattaneo, Feng, and Titiunik (2021) Cattaneo, M. D., Y. Feng, and R. Titiunik (2021): “Prediction Intervals for Synthetic Control Methods,” _Journal of the American Statistical Association_ , 116:536, 1865–1880, DOI: 10.1080/01621459.2021.1979561.
  * Chernozhukov, Wüthrich, and Zhu (2021) Chernozhukov, V., K. Wüthrich, and Y. Zhu (2021): “An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls,” _Journal of the American Statistical Association_ , 116:536, 1849–1864.
  * Chow and Lin (1971) Chow, G. C., and A. Lin (1971): “Best Linear Unbiased Interpolation, Distribution, and Extrapolation of Time Series by Related Series,” _The Review of Economics and Statistics_ , 53:4, 372–375.
  * Cochrane and Orcutt (1949) Cochrane, D., and G. H. Orcutt (1949): “Application of Least Squares Regression to Relationships Containing Autocorrelated Error Terms,” _Journal of the American Statistical Assocation_ , 44, 32–64.
  * Fan, Masini, and Medeiros (2022) Fan, J., R. Masini, and M. Medeiros (2022): “Do We Exploit All Information for Counterfactual Analysis? Benefits of Factor Models and Idiosyncractic Correction,” _Journal of the American Statistical Association_ , 117:538, 574–590.
  * Ferman (2023) Ferman, B. (2023): “Inference in Difference-in-Differences: How Much Should Should We Trust in Independent Clusters,” _Journal of Applied Econometrics_ , 38:3, 270–449.
  * Ferman and Pinto (2021) Ferman, B., and C. Pinto (2021): “Synthetic Controls with Imperfect Pretreatment Fit,” _Quantitative Economics_ , 12, 1197–1221.
  * Goldberger (1962) Goldberger, A. (1962): “Best Linear Unbiased Prediction in the Generalized Linear Regression Model,” _Journal of the American Statistical Association_ , 57, 369–375.
  * Henderson (1950) Henderson, C. R. (1950): “Estimation of Genetic Parameters,” _Annals of Mathematical Statistics_ , 21, 309–310.
  * Hsiao, Ching, and Wan (2012) Hsiao, C., H. Ching, and S. Wan (2012): “A Panel Data Approach for Program Evaluation: Measuring the Benefits of Political and Eocnomic Integration of Hong Kong with Mainland China,” _Journal of Applied Econometrics_ , 27, 705–740.
  * Lei and Candés (2021) Lei, L., and E. Candés (2021): “Conformal Inference of Counterfactuals and Individual Treatment Effects,” _Journal of the Royal Statistical Association Series B_ , 83, 911–938.
  * Li (2020) Li, K. (2020): “Statistical Inference for Average Treatment Effects Estimated by Synthetic Control Methods,” _Journal of the American Statistical Association_ , 115:532, 2068–2083.
  * Li, Shen, and Zhou (2023) Li, X., Y. Shen, and Q. Zhou (2023): “Confidence Intervals of Treatment Effects in Panel Data Models with Interactive Fixed Effects,” _Journal of Econometrics_ , 240:1, Article 105684.
  * Little and Rubin (2019) Little, R., and D. B. Rubin (2019): _Statistical Analysis with Missing Data_. Wiley, New York, third edn.
  * Masini and Medeiros (2022) Masini, R., and M. Medeiros (2022): “Counterfactual Analysis and Inference with Nonstationary Data,” _Journal of Business and Economic Statistics_ , 40:1, 227–239.
  * Masini and Medeiros (2021) Masini, R., and M. C. Medeiros (2021): “Counterfactual Analysis With Artificial Controls: Inference, High Dimensions, and Nonstationarity,” _Journal of the American Statistical Association_ , 116:536, 1773–1788, DOI: 10.1080/01621459.2021.1964978.
  * Ng and Scanlan (2024) Ng, S., and S. Scanlan (2024): “Constructing High Frequency Economic Indicators by Imputation,” _Econometrics Journal_ , 27:1, 1–30.
  * Phillips (1979) Phillips, P. (1979): “The Sampling Distribution of Forecasts From a First-Order Autoregression,” _Journal of Econometrics_ , 9, 241–261.
  * Robinson (1991) Robinson, G. K. (1991): “That BLUP is a Good Thing: The Estimation of Random Effects,” _Statistical Science_ , 6(1), 15 – 32, https://doi.org/10.1214/ss/1177011926.
  * Spall (1991) Spall, J. (1991): “Comment: The Kalman Filter and BLUP,” _Statistical Science_ , 6:1, 39–41.
  * Taub (1979) Taub, A. (1979): “Prediction in Context of Variance-Components Model,” _Journal of Econometrics_ , 10, 103–107.
  * Xu (2017) Xu, Y. (2017): “Generalized Synthetic Control Methods: Causal Inference with Interactive Fixed Effects Models,” _Political Analysis_ , pp. 56–75.