The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
23,976 characters
Bias-robust causal inference for panel data
\maketitle
\vspace{-3em}
\begin{abstract}
\noindent
We develop a bias-robust causal inference method for observational panel data
settings. Such methods typically impute untreated outcomes, so counterfactual
error passes straight into the estimated treatment effect while conventional
standard errors ignore it. We adapt bias-aware minimax methods, developed for
estimating regression coefficients in factor-model panels, to a causal target:
the average effect on the treated, which has to be imputed and may vary across
units and periods. The estimator corrects the imputed counterfactual with
weighted untreated residuals and reports intervals with an explicit allowance
for the error that remains. In simulations the proposed method holds nominal coverage where
alternatives such as the generalized synthetic control have almost none, especially when the factor rank is underfitted, at the cost of wider intervals. By applying the developed methodology to real data the estimated effect
remains significant for counterfactual errors nearly twice the size that the
design's placebos typically exhibit.
\end{abstract}
\section{Introduction}
Most policy settings do not offer random assignment, so researchers recover
missing counterfactuals from observational data. Synthetic control matches a
treated unit's pre-treatment path with a weighted combination of controls
\citep{abadie2010}; interactive-fixed-effects models instead use unit-specific
loadings on common factors \citep{bai2009}, and generalized synthetic control
(GSC) and matrix completion exploit that low-rank structure across many treated
units and adoption dates \citep{gobillon2016,xu2017,samartsidis2020,athey2021}.
This flexibility has a cost. GSC learns factors from control outcomes and each
treated unit's loading from its pre-treatment history, so counterfactual error
enters the treatment effect directly while conventional standard errors ignore
it. Related concerns motivate synthetic difference-in-differences
\citep{arkhangelsky2021}, augmented synthetic control \citep{benmichael2021},
and generalized-Bayes robustness \citep{alexopoulosdemiris2025};
\citet{rambachan2023} and \citet{callaway2021} address parallel trends rather
than counterfactual error.
The contribution of the present paper is to make causal inference bias-robust, using the
bias-aware minimax methods \citet{armstrong2022} developed for regression coefficients in the case of a complete panel dataset.
In the case of a causal target we extend this approach in three ways: the
estimand is defined over untreated potential outcomes that are systematically
missing, it averages effects that may differ across treated units, and the
correction is confined to observed untreated units. More precisely, we keep the GSC
estimand but choose weights on untreated residuals minimizing worst-case bias
plus variance, and report intervals that allow explicitly for the bias that
remains. Inference rests on an assumed bound for the counterfactual error.
Coverage follows from that bound under high-level conditions, but the bound
itself cannot be derived in this design; as in honest
difference-in-differences \citep{rambachan2023} we therefore report it rather
than estimate it, and summarize a conclusion by the largest error it survives.
Section 2 develops the estimator, its bias bound and its interval; Section 3
reports the simulations and Section 4 the real data application. Finally,
Section 5 concludes.
\section{Method}
For $T\times N$ matrices let
$\langle B,G\rangle=\operatorname{tr}(B'G)=\sum_{t,i}B_{ti}G_{ti}$ be the
entrywise inner product, which treats the two matrices as long vectors and takes
their dot product. Let $s_1(B)\geq s_2(B)\geq\cdots$ be the singular values of
$B$, which record how much of it lies along each of its orthogonal directions. Three summaries of
those values appear below: the operator norm $\|B\|_{\mathrm{op}}=s_1(B)$, the
Frobenius norm $\|B\|_F=(\sum_j s_j(B)^2)^{1/2}$ and the nuclear norm
$\|B\|_*=\sum_j s_j(B)$.
With $D_{ti}$ the treatment indicator, $\mathcal O=\{D_{ti}=0\}$ holds all
untreated observations, including the pre-adoption periods of units treated
later, and $\mathcal T=\{D_{ti}=1\}$ the rest.
Outcomes obey
\begin{equation}
Y_{ti}=\Gamma_{ti}+\tau_{ti}D_{ti}+\varepsilon_{ti},\qquad
\Gamma_{ti}=x_{ti}'\beta+\alpha_i+\delta_t+L_{ti},\quad
\operatorname{rank}(L)\leq R .
\label{eq:model}
\end{equation}
Here $\Gamma_{ti}=\operatorname{E}\{Y_{ti}(0)\mid
x_{ti},\alpha_i,\delta_t,L_{ti}\}$ and $\tau_{ti}=Y_{ti}(1)-Y_{ti}(0)$; we
assume consistency, no interference and no anticipation, and
$\operatorname{E}(\varepsilon\mid D,x,\alpha,\delta,L)=0$, so treatment may
depend on the latent untreated-outcome structure, and in particular on the
loadings, but not on the idiosyncratic errors.
Effects may differ across units and periods; with $W=D/|\mathcal T|$ the target
$\tau_W=\langle W,\tau\rangle$ averages them over all treated observations, the first
departure from the regression-coefficient problem. We fit \eqref{eq:model} by
rank-$R$ least squares on $\mathcal O$ alone and complete $\widehat\Gamma$
everywhere. For $A$ supported on $\mathcal O$, that is with $A_{\mathcal T}=0$, let
$\widehat\tau(A)=\langle W,Y-\widehat\Gamma\rangle-\langle
A,Y-\widehat\Gamma\rangle$: the first term is the usual imputation estimate, the
average gap between observed and imputed outcomes on treated observations, and
the second subtracts a weighted average of residuals from untreated
observations, chosen
below to remove as much of the counterfactual error as possible. Write
$E=\Gamma-\widehat\Gamma$ for the counterfactual error matrix and $H=W-A$ for
the matrix that carries it into the estimate, so that
\begin{equation}
\widehat\tau(A)-\tau_W=\langle H,E\rangle+\langle H,\varepsilon\rangle .
\label{eq:decomp}
\end{equation}
\subsection{Worst-case bias and minimax weights}
We choose the weights to minimize the mean squared error
$\operatorname{MSE}(A,E)=\operatorname{E}[\{\widehat\tau(A)-\tau_W\}^2]$. This
involves the counterfactual bias, the first term of \eqref{eq:decomp}, which is
not observed because $E$ is unknown, so the criterion cannot be minimized as it
stands. We therefore maximize it over a class $\mathcal E(C)$, where its radius $C$ bounds the counterfactual errors we are willing to assume.
For errors of that size we take the largest $\operatorname{MSE}(A,E)$ the estimator can have. Importantly, the maximum $\operatorname{MSE}(A,E)$ is free of $E$ and depends on
$A$ alone, so minimizing it is well posed; that minimax problem is what this
subsection develops. More precisely, following \citet{armstrong2022} we take the nuclear class, a ball of radius
$C$ about zero, $\mathcal E(C)=\{E:\|E\|_*\leq C\}$, on which duality gives the exact worst-case
bias
\begin{equation}
\sup_{E\in\mathcal E(C)}|\langle H,E\rangle|=C\,\|H\|_{\mathrm{op}},
\label{eq:duality}
\end{equation}
attained when $E$ is rank one, that is when the whole counterfactual error is a
single time path scaled by a single set of unit loadings. Turning to the second term of \eqref{eq:decomp} we have that $A$ vanishes on $\mathcal T$
and $W$ vanishes on $\mathcal O$, the two have disjoint supports, so
$\|H\|_F^2=\|W\|_F^2+\|A\|_F^2$. Moreover, the errors being uncorrelated with common
variance $\sigma^2$,
\begin{equation}
\operatorname{Var}\langle H,\varepsilon\rangle=\sigma^2\|H\|_F^2
=\sigma^2\bigl(\|W\|_F^2+\|A\|_F^2\bigr),
\label{eq:var}
\end{equation}
of which only $\sigma^2\|A\|_F^2$ depends on $A$. It is clear thus that the correction can reduce the
bias but only by adding variance of its own. Notice also that $A$ is restricted to $\mathcal O$ since only untreated residuals are pure error whereas on the treated observations the residual $Y-\widehat\Gamma$ contains the effect being estimated. If we estimate $\sigma^2$ by the mean
squared residual on $\mathcal O$, $\widehat\sigma^2=|\mathcal O|^{-1}
\sum_{(t,i)\in\mathcal O}(Y_{ti}-\widehat\Gamma_{ti})^2$ and by noting that $\|\cdot\|_{\mathrm{op}}^2$ is convex then,
\begin{equation}
\widehat A=\arg\min_{A:\ A_{\mathcal T}=0}\ \sup_{E\in\mathcal E(C)}
\operatorname{MSE}(A,E)
=\arg\min_{A:\ A_{\mathcal T}=0}
\bigl\{C^2\|W-A\|_{\mathrm{op}}^2+\widehat\sigma^2\|A\|_F^2\bigr\}
\label{eq:minimax}
\end{equation}
is a convex problem. By equation \eqref{eq:duality} the supremum touches only the bias,
giving $C^2\|H\|_{\mathrm{op}}^2$, and by \eqref{eq:var} the variance
contributes $\widehat\sigma^2\|A\|_F^2$ beyond
$\widehat\sigma^2\|W\|_F^2=\widehat\sigma^2/|\mathcal T|$, this is the noise any
estimator of $\tau_W$ carries through the treated observations, fixed by the
design, which shifts the objective without moving its minimizer and is dropped; see in the supplementary material for more details. Moreover, the objective in \eqref{eq:minimax} is homogeneous, so $\widehat A$ depends on $C$ and
$\widehat\sigma$ only through their ratio. At $C=0$ it returns $\widehat A=0$,
the uncorrected imputation; as $C$ grows so does the correction, but only until
weights confined to $\mathcal O$ can flatten the leading direction of $W$ no
further, beyond which the estimate is insensitive to $C$. In a block design
that limit binds at once: if $\mathcal T=P\times Q$ for some set of periods $P$
and units $Q$, then $\widehat A=0$ and $\widehat\tau$ reduces to the imputation
estimator with bias bound $C\|W\|_{\mathrm{op}}$. Since $A$ vanishes on
$\mathcal T$, the restriction of $W-A$ to $P\times Q$ is $W$, and selecting
rows and columns cannot raise an operator norm, so
$\|W-A\|_{\mathrm{op}}\geq\|W\|_{\mathrm{op}}$ throughout, and both terms of
\eqref{eq:minimax} are minimized at $A=0$. The correction draws its power from
staggered adoption, which leaves observations untreated in periods when others
are already treated; one treated unit, or several adopting together, leaves
none.
Let $\widehat H=W-\widehat A$, evaluating \eqref{eq:duality} at $\widehat H$
bounds the counterfactual bias by $b=C\|\widehat H\|_{\mathrm{op}}$, while
\eqref{eq:var} makes the noise standard deviation
$\widehat{\mathrm{se}}=\widehat\sigma\|\widehat H\|_F$. We report
$\widehat\tau(\widehat A)\pm(b+z_{1-\alpha/2}\widehat{\mathrm{se}})$: the
estimate sits off centre by at most $b$ and scatters about that point with
standard deviation $\widehat{\mathrm{se}}$, so allowing for the bias in full
and for the noise at the usual quantile covers whatever the true bias.
A sharper constant is available, because a bias of known size shifts the
sampling distribution rather than widening it \citep{armstrongkolesar2018}; the
saving is negligible here, since $m=b/\widehat{\mathrm{se}}$ is large, so we
report the conservative interval, which is also what \citet{armstrong2022}
implement.
Coverage then rests on two things: that $C$ really does bound the
counterfactual error, and that the noise is approximately normal on the scale
of $\widehat{\mathrm{se}}$. Granting both gives a guarantee at every $C$, which
is what makes the sensitivity analysis of Section 4 meaningful; Section~S4 of
the supplementary material proves the following.
\begin{proposition}[High-level sensitivity coverage]\label{prop:coverage}
Fix $C>0$ and suppose $\Pr(\|E\|_*\leq C)\to1$ and
$\langle\widehat H,\varepsilon\rangle/\widehat{\mathrm{se}}
\overset{d}{\to}N(0,1)$. Then
\[
\liminf\Pr\bigl\{\tau_W\in\widehat\tau(\widehat A)
\pm(b+z_{1-\alpha/2}\widehat{\mathrm{se}})\bigr\}\geq1-\alpha .
\]
\end{proposition}
This adapts the high-level coverage argument of \citet{armstrong2022}. We
maintain rather than verify its Gaussian-approximation and standard-error
conditions for the implemented factor estimator; Section 3 illustrates the
finite-sample coverage that results under the designs stated there, and does
not establish the conditions generally. The homoskedastic form assumed in
\eqref{eq:var} enters twice: in the criterion, where misspecification costs
only efficiency, and in $\widehat{\mathrm{se}}$, where it bears on coverage.
The proposition is stated in terms of $\widehat{\mathrm{se}}$, so a
dependence-robust replacement leaves it intact.
Finally, under their complete-panel assumptions and for their specified
preliminary estimator, \citet{armstrong2022} bound the relevant nuclear-norm
error by $3R\,s_1(\widehat U)(1+\epsilon)$, $\widehat U$ being the residual
matrix; the multiplier is $3R$ for the error $\Gamma-\widehat\Gamma$ and $2R$
for a differently defined target, the version they implement. Their theorem is
proved for an estimator that refits after an initial debiasing pass, and does
not reach a single rank-$R$ fit on $\mathcal O$ in a panel whose treated block
is never observed untreated. We therefore take
$\widehat C=3R\,s_1(\widehat U)(1+\epsilon)$ as a diagnostic reference rather
than a valid radius for this design; Section 3 reports how it compares with the
realized error there, and $C$ is treated as a sensitivity parameter throughout. What Proposition~\ref{prop:coverage}
leaves to be supplied is $C$ itself, so we report it as a sensitivity parameter
and summarize a conclusion by the largest $C$ it survives.
\section{Simulations}
We simulate $Y_{ti}=1.5f_{1t}\lambda_{1i}+\kappa f_{2t}\lambda_{2i}
+\varepsilon_{ti}$ with $N=T=40$ and no treatment effect, $f_1$ a standardized
AR(1) with coefficient $0.6$, $f_2$ a standardized trend, loadings and errors
standard normal. The $15$ units with the largest $\lambda_{2i}+0.25v_i$ adopt in
five three-unit cohorts at dates $20$--$34$, so selection follows exposure to
the trend and $\kappa$ scales the confounding. Over $200$ replications Monte
Carlo standard errors are at most $0.035$ for coverage and $0.04$ for RMSE.
In Panel A of Table~\ref{tab:sim} all methods fit the true rank. GSC coverage
falls from $0.47$ to $0.03$ as the confounding factor weakens, because its
intervals ignore counterfactual error; the proposed intervals hold coverage
throughout, and the minimax correction lowers bias from $0.77$ to $0.69$ at
$\kappa=1$. Panel B fits one factor when two are present. GSC coverage collapses to
$0.00$--$0.04$ while the proposed intervals retain $1.00$, and here the
diagnostic radius does become informative, since the omitted factor leaves
signal above the noise edge. Honest DiD, which fits no factor model, is
unchanged across panels and is narrower than the proposed interval in both.
In the correctly specified design the diagnostic radius exceeds the realized
$\|E\|_*$ in $99.5$--$100\%$ of replications and averages $1.25$--$1.34$ times it, so intervals
run $23$--$31\%$ wider than under the infeasible oracle radius at unchanged
coverage. The constant matters: $2R$ in place of $3R$ bounds $\|E\|_*$ in at
most $12.5\%$ of replications. Since $b/\widehat{\mathrm{se}}$ averages
$42$--$45$, the intervals are almost entirely bias allowance, and the sharp
critical value of Section 2 would shorten them by under one percent.
Figure~\ref{fig:sim} traces both statistics over the full $\kappa$ range. The
parallel-trends comparators are themselves sensitive to aggregation, which the
Supplementary Material documents.
\begin{table}[H]
\centering\footnotesize
\caption{Zero-effect simulations, $N=T=40$, 200 replications. The true ATT is zero, so Bias is the mean estimate. Cov is coverage of the nominal 95\% interval, Len its mean length. Proposed intervals use the oracle radius. Honest DiD returns a robust set, not a point estimate.}
\label{tab:sim}
\begin{tabular}{llcccc}
\toprule
\multicolumn{6}{l}{\textit{Panel A: correct rank ($R=2$)}} \\
Factors & Method & Cov & Len & Bias & RMSE \\
\midrule
\multirow{5}{*}{$\kappa=1$} & Proposed & 1.00 & 6.51 & 0.69 & 0.78 \\
& GSC & 0.47 & 1.49 & 0.77 & 0.86 \\
& C\&S DiD & 0.34 & 2.09 & 1.26 & 1.38 \\
& Honest DiD ($\bar M=0.5$) & 1.00 & 6.23 & -- & -- \\
\cmidrule(l){1-6}
\multirow{5}{*}{$\kappa=0.25$} & Proposed & 1.00 & 6.01 & 0.62 & 0.64 \\
& GSC & 0.03 & 0.65 & 0.72 & 0.74 \\
& C\&S DiD & 0.83 & 2.02 & 0.34 & 0.65 \\
& Honest DiD ($\bar M=0.5$) & 1.00 & 6.32 & -- & -- \\
\midrule
\multicolumn{6}{l}{\textit{Panel B: rank underfitted ($R=1$, truth $R=2$)}} \\
Factors & Method & Cov & Len & Bias & RMSE \\
\midrule
\multirow{5}{*}{$\kappa=1$} & Proposed & 1.00 & 9.82 & 2.38 & 2.42 \\
& GSC & 0.00 & 1.78 & 2.81 & 2.86 \\
& C\&S DiD & 0.34 & 2.09 & 1.26 & 1.38 \\
& Honest DiD ($\bar M=0.5$) & 1.00 & 6.23 & -- & -- \\
\bottomrule
\end{tabular}
\end{table}
\begin{figure}[H]
\centering
\includegraphics[width=0.9\linewidth]{fig_sim_main.pdf}
\caption{Coverage and absolute ATT bias by factor strength; true ATT is zero.
Left column fits the true rank, right fits one factor when two are present.
Dashed line is nominal $95\%$. The two bias panels have separate vertical
scales.}
\label{fig:sim}
\end{figure}
\section{Real Data Analysis}
We revisit the Election Day Registration (EDR) application of \citet{xu2017}:
turnout for $47$ states over $24$ presidential elections, nine adopting in four
cohorts, with mail-in and motor-voter controls, two-way effects and $R=2$. Because the actual effect is unknown, we first measure bias on $50$ confounded donor-state placebos with known zero effect, assigning earlier placebo adoption
to the larger pre-1976 turnout trend. Table~\ref{tab:edr} reports the outcome.
The proposed estimator improves on GSC in bias, MAE and RMSE and covers at
$1.00$; GSC covers at $0.96$ here, so in this design the gain over it lies in
the point estimate rather than in coverage, while Callaway--Sant'Anna is badly
biased and covers $0.28$. The price is intervals three to five times wider.
\begin{table}[H]
\centering\footnotesize
\caption{Fifty confounded donor-state placebos with a known zero effect, so Bias is the mean estimate. Cov is coverage of the nominal 95\% interval, Len its mean length. Monte Carlo standard errors are at most 0.33 for bias, 0.24 for RMSE and 0.06 for coverage.}
\label{tab:edr}
\begin{tabular}{lccccc}
\toprule
Method & Bias & MAE & RMSE & Cov & Len \\
\midrule
Proposed & -0.79 & 1.94 & 2.37 & 1.00 & 37.80 \\
GSC & -0.88 & 2.02 & 2.45 & 0.96 & 10.51 \\
C\&S DiD & 4.62 & 4.62 & 4.93 & 0.28 & 7.36 \\
\bottomrule
\end{tabular}
\end{table}
For the actual policy GSC estimates an ATT of $4.90$ points with interval
$[0.29,9.50]$. The correction leaves the estimate at $4.99$ and, as
Figure~\ref{fig:edr} shows, it barely moves across the sweep, from $4.895$ to
$4.989$ as $C$ runs from $0$ to $400$: the weights saturate, so the radius
governs the interval rather than the estimate.
The half-width grows in $C$
while the estimate does not, the correction saturates, as Section~S3 of the
supplementary material shows, so there is a largest bound whose interval
still excludes zero: the breakdown bound, here found to be $37.9$. Whether this
is large depends on the benchmark, and the comparison is best made on the scale
of the bias itself. A placebo effect is zero, so a placebo estimate is
$\langle H,E\rangle+\langle H,\varepsilon\rangle$; netting out the sampling
part leaves a typical counterfactual bias of $2.33$ points against the $4.15$
the interval absorbs at the breakdown, while the $95$th percentile of the
placebo magnitudes reaches $4.23$. The conclusion therefore tolerates about
$1.8$ times the bias this design usually produces but falls just short of the
error it occasionally delivers; the shaded band of Figure~\ref{fig:edr} spans
the two, and the breakdown sits inside its upper edge.
One placebo in fifty is as large as the estimated $4.99$. At the diagnostic
radius $\widehat C=266.0$, deliberately conservative, the conclusion fails.
Honest
DiD is less encouraging, its breakdown $\bar M$ being $0$, since even under
exact parallel trends the robust set $[-0.98,4.63]$ contains zero.
\begin{figure}[H]
\centering
\includegraphics[width=0.9\linewidth]{fig_edr_sensitivity.pdf}
\caption{The EDR interval against the assumed radius $C$. The shaded band spans
the typical and $95$th-percentile placebo errors, $21.3$ to $38.6$ once divided
by $\|\widehat H\|_{\mathrm{op}}$ to make them commensurate with $C$; the
interval first contains zero at the breakdown bound $37.9$, while
$\widehat C=266.0$ is the diagnostic radius of Section 2.}
\label{fig:edr}
\end{figure}
\section{Conclusion}
We transfer minimax bias-aware inference from panel estimation to causal
inference. The resulting coverage is superior than the one of recently developed methods such as the generalized synthetic control, especially when the factor rank is wrong. The price is a bound
that must be assumed rather than estimated. This is less a concession than a
change of accounting: generalized synthetic control also needs the
counterfactual error to be small and simply never says how small, whereas here
that error is a single number one can vary, benchmark against placebos, and
summarize by the largest value a conclusion survives, the discipline honest
difference-in-differences applies to parallel trends. What remains open is a
bound derived rather than verified. Extending Theorem~2 of
\citet{armstrong2022} to a panel whose treated block is never observed would
close the one step this paper checks by simulation, and a bound proved for this
design should be less conservative than the one we import.
\setlength{\bibsep}{1pt plus 0.3ex}
\begin{thebibliography}{9}
\bibitem[Abadie et al.(2010)]{abadie2010}
Abadie, A., Diamond, A., Hainmueller, J., 2010. Synthetic control methods for
comparative case studies. \emph{J. Am. Stat. Assoc.} 105, 493--505.
\bibitem[Alexopoulos and Demiris(2025)]{alexopoulosdemiris2025}
Alexopoulos, A., Demiris, N., 2025. On robust Bayesian causal inference. arXiv
preprint arXiv:2511.13895.
\bibitem[Arkhangelsky et al.(2021)]{arkhangelsky2021}
Arkhangelsky, D., Athey, S., Hirshberg, D.A., Imbens, G.W., Wager, S., 2021.
Synthetic difference-in-differences. \emph{Am. Econ. Rev.} 111, 4088--4118.
\bibitem[Armstrong and Kolesár(2018)]{armstrongkolesar2018}
Armstrong, T.B., Kolesár, M., 2018. Optimal inference in a class of regression
models. \emph{Econometrica} 86, 655--683.
\bibitem[Armstrong et al.(2026)]{armstrong2022}
Armstrong, T.B., Weidner, M., Zeleneev, A., 2026.
Robust estimation and inference in panels with interactive fixed effects.
\emph{Journal of Political Economy} 134(8), in press.
\bibitem[Athey et al.(2021)]{athey2021}
Athey, S., Bayati, M., Doudchenko, N., Imbens, G.W., Khosravi, K., 2021. Matrix
completion methods for causal panel data models. \emph{J. Am. Stat. Assoc.} 116,
1716--1730.
\bibitem[Bai(2009)]{bai2009}
Bai, J., 2009. Panel data models with interactive fixed effects.
\emph{Econometrica} 77, 1229--1279.
\bibitem[Ben-Michael et al.(2021)]{benmichael2021}
Ben-Michael, E., Feller, A., Rothstein, J., 2021. The augmented synthetic
control method. \emph{J. Am. Stat. Assoc.} 116, 1789--1803.
\bibitem[Callaway and Sant'Anna(2021)]{callaway2021}
Callaway, B., Sant'Anna, P.H.C., 2021. Difference-in-differences with multiple
time periods. \emph{J. Econometrics} 225, 200--230.
\bibitem[Gobillon and Magnac(2016)]{gobillon2016}
Gobillon, L., Magnac, T., 2016. Regional policy evaluation: interactive fixed
effects and synthetic controls. \emph{Rev. Econ. Stat.} 98, 535--551.
\bibitem[Rambachan and Roth(2023)]{rambachan2023}
Rambachan, A., Roth, J., 2023. A more credible approach to parallel trends.
\emph{Rev. Econ. Stud.} 90, 2555--2591.
\bibitem[Samartsidis et al.(2020)]{samartsidis2020}
Samartsidis, P., Seaman, S.R., Montagna, S., Charlett, A., Hickman, M.,
De Angelis, D., 2020. A Bayesian multivariate factor analysis model for
evaluating an intervention by using observational time series data on multiple
outcomes. \emph{J. R. Stat. Soc. Ser. A} 183, 1437--1459.
\bibitem[Xu(2017)]{xu2017}
Xu, Y., 2017. Generalized synthetic control method. \emph{Polit. Anal.} 25,
57--76.
\end{thebibliography}
\clearpage
\setcounter{figure}{0}
\setcounter{table}{0}
\begin{center}
{\large\bfseries Supplementary material}
\end{center}
\medskip