#  Extensions for Inference in Difference-in-Differences with Few Treated
Clusters

Luis Alvarez111Sao Paulo School of Economics - FGV; email:
<EMAIL_ADDRESS>Bruno Ferman222Sao Paulo School of Economics -
FGV; email<EMAIL_ADDRESS>
Abstract

In settings with few treated units, Difference-in-Differences (DID) estimators
are not consistent, and are not generally asymptotically normal. This poses
relevant challenges for inference. While there are inference methods that are
valid in these settings, some of these alternatives are not readily available
when there is variation in treatment timing and heterogeneous treatment
effects; or for deriving uniform confidence bands for event-study plots. We
present alternatives in settings with few treated units that are valid with
variation in treatment timing and/or that allow for uniform confidence bands.

Keywords: inference; difference-in-differences; permutation tests;
heterogeneous treatment effects; uniform confidence bands

JEL Codes: C12; C21; C33

##  1 Introduction

In settings with few treated units, Difference-in-Differences (DID) estimators
are not consistent, and are not generally asymptotically normal, posing
relevant challenges for inference (Donald and Lang,, 2007; Conley and Taber,,
2011). In such settings, relying on standard procedures, such as clustering
the standard errors at the unit level, may lead to severe over-rejection. In
some examples, we may expect rejection rates on the order of more than 60% for
a 5% nominal-level test (Ferman,, 2019). In light of these concerns, some
alternative inference methods have been proposed for settings in which we have
a small number of treated units (Conley and Taber,, 2011; Ferman and Pinto,,
2019; MacKinnon and Webb,, 2019; Hagemann,, 2020).

However, some of these alternatives are not readily available to incorporate
some recent advances in the DID literature. Such advances include (i)
considering settings with heterogeneous treatment effects and variation in the
treatment timing (Goodman-Bacon,, 2021; de Chaisemartin and D’Haultfoeuille,
2020a, ; de Chaisemartin and D’Haultfoeuille, 2020b, ; Callaway and
Sant’Anna,, 2021; Sun and Abraham,, 2020), and (ii) relying on uniform
confidence bands when presenting event-study plots (Freyaldenhoven et al.,,
2021; Callaway and Sant’Anna,, 2021).

In this note, we extend the inference methods proposed by Conley and Taber,
(2011) and by Ferman and Pinto, (2019) so that (i) we can incorporate the
recent recommendations from papers that analyzed settings with heterogeneous
treatment effects and variation in the treatment timing, and/or (ii) we can
allow for uniform confidence bands when considering event-study plots. Conley
and Taber, (2011) discuss in their Section III.A the possibility of extending
their approach to heterogeneous treatment effects. In this note, we formalize
this idea by filling in important implementation details, and we show how
parametric models for heteroskedasticity, as proposed by Ferman and Pinto,
(2019), can be extended to a staggered adoption setting. Conley and Taber,
(2011) also discuss the possibility of computing joint confidence sets based
on the inversion of a test statistic. Differently, we show how one can conduct
uniform inference in this setting by relying on uniform confidence bands,
which have been advocated in the literature due their ease of interpretability
and computation (Montiel Olea and Plagborg-Moller,, 2019). Importantly, due to
the nonstandard nature of the setting, critical values used in the uniform
bands are not based on a known distribution. We show that a specific bootstrap
algorithm is able to recover the required critical values in an asymptotic
framework where the number of treated units is fixed and the number of
controls is large.

##  2 Setting

Let $y_{j,t}(0)$ ($y_{j,t}(1)$) be the potential outcome of unit $j$ at time
$t$ when this unit is untreated (treated) at this period. We consider first
that potential outcomes are given by

$\displaystyle\begin{cases}y_{j,t}(0)=\theta_{j}+\gamma_{t}+\eta_{j,t}\\\
y_{j,t}(1)=\alpha_{j,t}+y_{j,t}(0),\end{cases}$ (1)

where $\theta_{j}$ and $\gamma_{t}$ are, respectively, unit- and time-
invariant unobserved variables, while $\eta_{j,t}$ represents unobserved
variables that may vary at both dimensions; $\alpha_{j,t}$ is the (possibly
heterogeneous) treatment effect on unit $j$ at time $t$. We consider
$\alpha_{j,t}$ as a fixed parameter, which means that we define the target
parameters based on the realized treatment effect of this policy for the
treated units. In this case, uncertainty regarding this parameter comes from
unobserved shocks that may affect the potential outcomes of unit $j$ at time
$t$, such as, for example, weather or economic shocks that are unobserved by
the econometrician ($\eta_{j,t}$). A setting in which treatment assignment and
treatment effects are treated as fixed is common in the literature of DID with
few treated clusters (Conley and Taber,, 2011; Ferman and Pinto,, 2019;
Ferman,, 2020). A similar setting is also considered in other settings in
which the number of treated clusters is fixed, such as in the synthetic
controls literature (Abadie et al.,, 2010; Arkhangelsky et al.,, 2021; Ben-
Michael et al.,, 2021; Chernozhukov et al.,, 2021; Ferman and Pinto,, 2021;
Ferman,, 2021).

Units $j=1,...,N_{1}$ are treated at some point, while units $j=N_{1}+1,...,N$
are never treated, where $N_{0}=N-N_{1}$. We observe information for periods
$t=1,...,T$, and we allow for variation in treatment timing by denoting
$t^{\ast}_{j}\in\\{1,...,T-1\\}$ as the last period before unit $j$ enters
into treatment, for $j\leq N_{1}$. Let $t^{\ast}_{j}=\infty$ for $j>N_{1}$.
Treatment is assumed to be an absorbing state and we assume there is no
anticipation, so observed outcomes are given by
$y_{jt}=\mathbf{1}\\{t>t_{j}^{*}\\}y_{jt}(1)+\mathbf{1}\\{t\leq
t_{j}^{*}\\}y_{jt}(0)$.

We define the target parameter
$\bar{\alpha}=\sum_{i=1}^{N_{1}}\sum_{t=t_{j}^{*}+1}^{T}\omega_{j,t}\alpha_{j,t}$,
which is a linear combination of the treatment effects of different units in
different periods. One example for a target parameter is the average treatment
effects for the treated units in the periods that they were treated.
Alternatively, we can consider, for example, weighted averages depending on
populations sizes.

We might also be interested in a multivariate parameter
$\bar{\boldsymbol{\alpha}}=(\bar{\alpha}_{1},...,\bar{\alpha}_{K})$, in which
case we set
$\bar{\alpha}_{k}=\sum_{i=1}^{N_{1}}\sum_{t=t_{j}^{*}+1}^{T}\omega_{k,j,t}\alpha_{j,t}$.
For example, this may include the treatment effects $\tau$ periods after the
start of the treatment. With some abuse of notation, we can also consider that
some of these $\bar{\alpha}_{k}$ include pre-treatment trends, which are
commonly presented in dynamic DID models as an assessment for the parallel
trends assumption. We provide an example of the latter in the next sections.

Let $\boldsymbol{\eta}_{j}=(\eta_{j,1},...,\eta_{j_{T}})^{\prime}$.
Uncertainty comes from different realizations of
$\\{\boldsymbol{\eta}_{j}\\}_{j=1}^{N}$, where we treat treatment allocation
as fixed. We also consider the possibility of a vector of observable variables
$Z_{j}$ that may be determinants of the heteroskedasticity, as we discuss
below. These variables are also treated as nonrandom throughout.333We can also
consider the case in which some of these variables enter in the model for
$y_{it}(0)$ in (1).

###### Assumption 1

(Sampling) For each $j=1,...,N$, we observe $(y_{j,1},...,y_{j,T},Z_{j})$,
where $y_{j,t}=y_{j,t}(0)$ if $t<t^{\ast}_{j}$, and $y_{j,t}=y_{j,t}(1)$
otherwise. Potential outcomes are determined by Equation 1. We assume that
$(\boldsymbol{\eta}_{1},...,\boldsymbol{\eta}_{N})$ are mutually independent,
with finite second moments.

We do not need to impose any assumptions on $\theta_{j}$ and $\gamma_{t}$.

##  3 Estimator

An important limitation of the TWFE estimator in this setting is that it may
recover a linear combination of the treatment effects $\alpha_{j,t}$ in which
some of the weights might be negative (de Chaisemartin and D’Haultfoeuille,
2020b, ; Goodman-Bacon,, 2021; Borusyak et al.,, 2021). Most of the solutions
in this case consider alternative estimators that combine simpler $2\times 2$
estimators (Sun and Abraham,, 2020; Callaway and Sant’Anna,, 2021; Borusyak et
al.,, 2021).444Indeed, Borusyak et al., (2021) show that, under Assumptions 1
and 2, all linear-in-outcomes unbiased estimators of
$\boldsymbol{\bar{\alpha}}$ take an “imputation” form, of which aggregation of
simple 2 $\times$ 2 DID estimators constitute a particular case.,555Callaway
and Sant’Anna, (2021) consider a doubly robust estimation method for each of
those $2\times 2$ estimators, instead of considering simple DID estimators.
Even though their estimator is not consistent in our setting, our inference
method remains valid for weighted averages of the treatment effect, where the
weights are given by the estimated propensity score, provided that the outcome
model is correctly specified. Indeed, their estimator is not doubly robust in
our setting with a few number of treated units. We follow a similar approach.
More specifically, for each $\alpha_{j,t}$ in which $\omega_{k,j,t}>0$ for
some $(k,j,t)$, we consider an estimator $\hat{\alpha}_{j,t}$, and then we
aggregate them to get
$\widehat{\bar{\alpha}}_{k}=\sum_{i=1}^{N_{1}}\sum_{t=t_{j}^{*}+1}^{T}\omega_{k,j,t}\hat{\alpha}_{j,t}$.
We note that Conley and Taber, (2011) already recommended estimating each
$\alpha_{j,t}$ separately in settings with heterogeneous treatment effects,
even before the recent papers that pointed out the aggregation problems of the
TWFE estimator.

For a given $(\nu_{1}(j,t),...,\nu_{t^{\ast}_{j}}(j,t))$ that satisfies
$\sum_{\tau=1}^{t^{\ast}_{j}}\nu_{\tau}(j,t)=1$, we consider estimators
$\hat{\alpha}_{j,t}$ of the form

$\displaystyle\hat{\alpha}_{j,t}=\left[y_{j,t}-\sum_{\tau=1}^{t^{\ast}_{j}}\nu_{\tau}(j,t)y_{j,\tau}\right]-\frac{1}{N_{0}}\sum_{j^{\prime}=N_{1}+1}^{N}\left[y_{j^{\prime},t}-\sum_{\tau=1}^{t^{\ast}_{j}}\nu_{\tau}(j,t)y_{j^{\prime},\tau}\right].$
(2)

That is, for the post-pre comparison, we compare period $t$ with a weighted
average of the pre-treatment periods given by the weights $\nu_{\tau}(j,t)$.
Then we compare this post-pre comparison for the treated unit $j$ with the
average of the controls.666We do not consider the possibility of using the
not-yet-treated as controls, because this would be asymptotically irrelevant
when $N_{1}$ is fixed and $N_{0}\rightarrow\infty$, and because this would
make the notation more complicated. Still, it is possible to consider this
alternative.  Two simple examples include

$\displaystyle\hat{\alpha}_{j,t}=\left[y_{j,t}-\frac{1}{t^{\ast}_{j}}\sum_{\tau=1}^{t^{\ast}_{j}}y_{j,\tau}\right]-\frac{1}{N_{0}}\sum_{j^{\prime}=N_{1}+1}^{N}\left[y_{j^{\prime},t}-\frac{1}{t^{\ast}_{j}}\sum_{\tau=1}^{t^{\ast}_{j}}y_{j^{\prime},\tau}\right],$
(3)

where we consider a DID estimator using unit $j$ and the never treated, for
all pre-treatment periods and for period $t$. Alternatively, we can use only
the last period before unit $j$ starts treatment as the pre-period, so that

$\displaystyle\hat{\alpha}_{j,t}^{\prime}=\left[y_{j,t}-y_{i,t^{\ast}_{j}}\right]-\frac{1}{N_{0}}\sum_{j^{\prime}=N_{1}+1}^{N}\left[y_{j^{\prime},t}-y_{j^{\prime},t^{\ast}_{j}}\right].$
(4)

We aggregate these estimators the following way. Let
$\mathbf{Y}_{j}=(y_{j,1},...,y_{j,T})$. For each $j=1,\ldots,N_{1}$, we define
a $(K_{j}\times T)$ matrix $A_{j}$ such that
$A_{j}[\mathbf{Y}_{j}-\frac{1}{N_{0}}\sum_{i=N_{1}+1}^{N}\mathbf{Y}_{i}]$
consists of stacked estimators of $\hat{\alpha}_{jt}$ for a subset of the
periods. We provide examples of choices of $A_{j}$ below. We then aggregate
these building-block estimators onto a $K$-dimensional estimator of
$\bar{\boldsymbol{\alpha}}$ through $(K\times K_{j})$ matrices $B_{j}$. The
resulting estimator is given by
$\widehat{\bar{\boldsymbol{\alpha}}}\equiv\sum_{j=1}^{N_{1}}B_{j}A_{j}\left[\mathbf{Y}_{j}-\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\boldsymbol{Y}_{i}\right]$,
whereas the target parameter may be written as
$\bar{\boldsymbol{\alpha}}=\sum_{j=1}^{N_{1}}B_{j}A_{j}\boldsymbol{\alpha}_{j}$,
where
$\boldsymbol{\alpha}_{j}=(0,\ldots,0,\alpha_{j,t^{*}_{j}+1},\ldots,\alpha_{jT})^{\prime}$
is the vector of identifiable treatment effects for unit $j$.

We consider the following parallel trends assumption.

###### Assumption 2

(Parallel trends) For each $j=1,\ldots,N_{1}$, we assume that:

$B_{j}A_{j}\mathbb{E}[\boldsymbol{\eta}_{j}]=B_{j}A_{j}\mathbb{E}[\boldsymbol{\eta}_{i}],\
\forall i\in\\{N_{1}+1,\ldots,N\\}\,.$

We also assume that each row of $B_{j}A_{j}$ sums up to zero.

Assumption 2 guarantees the relevant parallel trends conditions for the
estimator $\widehat{\bar{\alpha}}$ of our choice. If we use all pre-treatment
periods the estimation of each $\alpha_{j,t}$, then Assumption 2 in practice
means that we have parallel trends for all periods.777Assumption 2 would
actually be weaker than that, as it may be satisfied without assuming parallel
trends for all periods. In this case, however, we would need departures from
parallel trends to cancel out, so that this assumption is satisfied.  In
contrast, if we consider the estimator $\hat{\alpha}_{j,t}^{\prime}$, then we
only need parallel trends between period $t$ and the last period before
treatment for unit $j$. The requirement that the rows of $B_{j}A_{j}$ sum up
to zero is made so the estimator removes unit fixed effects $\theta_{j}$, and
it is satisfied for our generic estimator presented in Equation 2, given the
restriction $\sum_{\tau=1}^{t^{\ast}_{j}}\nu_{\tau}(j,t)=1$.

######  Example 1 (label=exa:scalar)

Suppose we are interested in the average post-treatment effect on the treated
units, and that we use all pre-treatment periods in constructing our
estimator. For $j\leq N_{1}$, put $K_{j}=T-t_{j}^{*}$. In this case, the
matrices $A_{j}$ are taken as

$A_{j}=\begin{bmatrix}-\frac{1}{t_{j}^{*}}\mathbf{1}_{K_{j}\times
t_{j}^{*}}&\mathbb{I}_{K_{j}}\end{bmatrix},$ (5)

where $\mathbf{1}_{a\times b}$ is a $a\times b$ matrix of ones and
$\mathbb{I}_{c}$ is a $c\times c$ identity matrix. The matrices $B_{j}$ are
taken as

$B_{j}=\frac{1}{\sum_{i=1}^{N_{1}}K_{i}}\mathbf{1}_{1\times K_{j}}\,,$

and Assumption 2 requires that, for each $j=1,\ldots,N_{1}$,

$\sum_{t=t_{j}^{*}+1}^{T}\frac{\mathbb{E}[\eta_{jt}]}{(T-t_{j}^{*})}-\sum_{t=1}^{t_{j}^{*}}\frac{\mathbb{E}[\eta_{jt}]}{t_{j}^{*}}=\sum_{t=t_{j}^{*}+1}^{T}\frac{\mathbb{E}[\eta_{it}]}{(T-t_{j}^{*})}-\sum_{t=1}^{t_{j}^{*}}\frac{\mathbb{E}[\eta_{it}]}{t_{j}^{*}},\forall
i\in\\{N_{1}+1,\ldots N_{0}\\}\,.$

######  Example 2 (label=exa:multiv)

Suppose that we are interested on average treatment effects by length of
exposure. Suppose that we use all pre-treatment periods in constructing our
estimator. In this case, the $A_{j}$ can be taken as in the previous example.
The $B_{j}$ are constructed as follows. For each $k\in\mathbb{N}$, define
$N_{k}=\sum_{i=1}^{N_{1}}\mathbf{1}_{k\leq T-t_{j}^{*}}$. Let
$K=\max_{i=1,\ldots,N_{1}}K_{i}$. We define $B_{j}$ as:

$B_{j}=\begin{bmatrix}\operatorname{diag}(1/N_{1},1/N_{2},\ldots,1/{N_{K_{j}}})\\\
\boldsymbol{0}_{K-K_{j}\times K_{j}}\,\end{bmatrix}$ (6)

and Assumption 2 requires that, for each $j=1,\ldots,N_{1}$ and
$t^{\prime}>t_{j}^{*}$:

$\mathbb{E}[\eta_{jt^{\prime}}]-\sum_{t=1}^{t_{j}^{*}}\frac{\mathbb{E}[\eta_{jt}]}{t_{j}^{*}}=\mathbb{E}[\eta_{it^{\prime}}]-\sum_{t=1}^{t_{j}^{*}}\frac{\mathbb{E}[\eta_{it}]}{t_{j}^{*}},\forall
i\in\\{N_{1}+1,\ldots N_{0}\\}\,,$

###### Example 3

Suppose we are interested in simunalteously conducting inference on pre- and
post-treatment differential trends between treatment and control groups. We
use the term “post-treatment differential trends” because, when there is
possibly violation of parallel trends, the difference-in-differences estimand
conflates treatment effects with differential trends in the untretated
potential outcome. In this setting, with some abuse of notation, we define the
building-block estimands as, for $j\leq N_{1}$ and $t=1\ldots T$:

$\alpha_{jt}=\mathbb{E}[y_{jt}-y_{jt^{*}_{j}}]-\frac{1}{N_{1}}\sum_{i=N_{1}+1}^{N}\mathbb{E}[y_{it}-y_{it^{*}_{j}}],$

i.e. the differential trend between unit $j$ and the average of controls with
respect to the last pre-treatment period.888In this case, if Model (1) is
valid, and Assumption 2 holds, then this means that $\alpha_{jt}$ is the
treatment effect for $t$ in the post-treatment periods, and it is zero for the
pre-treatment periods. More generally, $\alpha_{jt}$ conflates treatment
effects with violations of the parallel trends assumption. Estimation of these
building-block parameters may be conducted by setting

$A_{j}=\begin{bmatrix}\mathbb{I}_{t_{j}^{*}-1}&-\boldsymbol{1}_{t_{j}^{*}-1\times
1}&\boldsymbol{0}_{t_{j}^{*}-1\times T-t_{j}^{*}}\\\ \boldsymbol{0}_{1\times
t_{j}^{*}-1}&0&\boldsymbol{0}_{1\times T-t_{j}^{*}-1}\\\
\boldsymbol{0}_{T-t_{j}\times t_{j}^{*}-1}&-\boldsymbol{1}_{T-t_{j}^{*}\times
1}&\mathbb{I}_{T-t_{j}^{*}}\end{bmatrix}\,.$

We may then aggregate these parameters by “length of exposure”. Specifically,
we set $L_{j}=(1-t_{j}^{*})$, $U_{j}=(T-t_{j}^{*})$, $L=\min_{j\leq
N_{1}}L_{j}$, $U=\min_{j\leq N_{1}}U_{j}$. For $k\in\mathbb{Z}$, we set
$N_{k}=\sum_{j=1}^{N_{1}}(\mathbf{1}_{T-t_{j}^{*}\leq
k}+\mathbf{1}_{(1-t_{j}^{*})\leq k})$. We then put:

$B_{j}=\begin{bmatrix}\boldsymbol{0}_{(L_{j}-L)\times T}\\\
\operatorname{diag}(1/N_{L_{j}},1/N_{L_{j}+1},\ldots 1/N_{0},\ldots
1/N_{U_{j}})\\\ \boldsymbol{0}_{(U-U_{j})\times T}\,.\end{bmatrix}$

The next proposition shows our estimator is unbiased, albeit inconsistent.

###### Proposition 1

Under Assumptions 1 and 2, the estimator $\widehat{\bar{\alpha}}$ is unbiased
for $\bar{\alpha}$. Moreover,

$\widehat{\bar{\boldsymbol{\alpha}}}\rightarrow_{p}\bar{\boldsymbol{\alpha}}+\sum_{j=1}^{N_{1}}B_{j}A_{j}\boldsymbol{\eta}_{j}$

when $N_{0}\rightarrow\infty$ and $N_{1}$ is fixed.

This proposition is equivalent to Proposition 1 from Conley and Taber, (2011).
While the DID estimator is unbiased, it is not consistent, because the number
of treated units is fixed.

It would also be possible to use the not-yet-treated as part of the control
group. Since the asymptotic theory we consider in this paper considers $N_{1}$
fixed and $N_{0}\rightarrow\infty$, this would not affect any of our
asymptotic results. For the finite-sample result that the estimator
$\widehat{\bar{\alpha}}$ is unbiased, we would have to adjust Assumption 2 to
include parallel trends for the not-yet-treated.

##  4 Inference method

The fact that $\widehat{\bar{\alpha}}$ is inconsistent poses some important
challenges for inference. Conley and Taber, (2011) propose an interesting
alternative, in which the residuals from the control units can be used to
estimate the distribution of the errors of the treated units. In their
standard implementation, they consider a setting in which
$(\boldsymbol{\eta}_{1},...,\boldsymbol{\eta}_{N})$ is iid. Ferman and Pinto,
(2019) analyze the case in which heteroskedasticity has a known structure that
can be estimated from the data. For example, in case the treated units are
state $\times$ time aggregates of individual level observations, then the idea
is to estimate the heteroskedasticity using the residuals from the control
units, and then use this estimated structure to make the residuals of the
controls informative about the distribution of the errors of the
treated.999Conley and Taber, (2011) also consider in their appendix an
alternative that allows for heteroskedasticity (based on a known variable) and
spatial correlation (depending on an observed distance metric). The method
proposed by Ferman and Pinto, (2019) differs from this alternative in that it
does not require parametrization/estimation of the serial correlation
structure, and that it does not assume normality.

Ferman and Pinto, (2019) focus on estimating a unidimensional parameter, and
do not take into account the recent advances in the analysis of staggered
designs with heterogeneous effects. In their Section III.A, Conley and Taber,
(2011) discuss the possibility of extending their approach to settings in
which there are heterogeneous treatment effects. In this note, we formalize a
procedure to draw inference in this setting, taking into account that
estimators for the building blocks $\alpha_{i,t}$ and $\alpha_{i,t^{\prime}}$
are possibly correlated. We also show how parametric models for
heteroskedasticity, as proposed by Ferman and Pinto, (2019), can be extended
to this setting.

Conley and Taber, (2011) also suggest constructing confidence sets on
multidimensional parameters in this setting by inverting a test statistic.
However, confidence sets constructed from quadratic statistics (the typical
choice, e.g. Wald test statistics) will be ellipsoidal, which is difficult to
both compute and visualize even for moderate dimensions of the target
parameter (Montiel Olea and Plagborg-Moller,, 2019). This is particularly
concerning for presenting confidence sets for event-study-like parameters. In
contrast, we follow the recent literature on DID and event-studies in
proposing uniform confidence bands for the target parameter (Freyaldenhoven et
al.,, 2021; Callaway and Sant’Anna,, 2021). We show how such uniform
confidence bands can be computed in a non-standard setting in which the
estimator is not consistent.

The next assumption imposes a parametric model for the heteroskedacity of
$B_{j}A_{j}\boldsymbol{\eta}_{j}$, and can be seen as an extension of the
assumptions considered by Ferman and Pinto, (2019). Observe that Assumption 2
implies that, for each $j\leq N_{1}$, $B_{j}A_{j}\boldsymbol{\eta}_{i}$ has
equal mean for all $i\in\\{j\\}\cup\\{N_{1},\ldots,N\\}$. Without loss, we
assume that this mean is zero and that:

###### Assumption 3

For each $j=1,\ldots N_{1}$, and for every
$i\in\\{j\\}\cup\\{N_{1}+1,\ldots,N\\}$:

$B_{j}A_{j}\boldsymbol{\eta}_{i}=H_{j}(Z_{i};\delta_{j})\xi_{i,j},$

where $\delta_{j}\in\Delta_{j}\subseteq\mathbb{R}^{p_{j}}$ is an unknown
parameter, $H_{j}(\cdot;\cdot)$ is a known function such that
$H_{j}(Z_{i};\delta_{j})$ is a positive definite $K\times K$ matrix for each
$i\in\\{j\\}\cup\\{N_{1}+1,\ldots,N\\}$; and
$\xi_{i,j}\overset{d}{=}\xi_{i^{\prime},j}$ for every
$i,i^{\prime}\in\\{j\\}\cup\\{N_{1}+1,\ldots,N\\}$, where $\xi_{j,j}$ is a
$K\times 1$ random variable with a continuous distribution and
$\mathbb{E}[\xi_{j,j}]=0$.

Note that, if we set $H_{j}(Z_{i};\delta_{j})$ as the identity matrix, then
Assumption 3 would be implied by Assumption 2 from Conley and Taber, (2011),
which states that $\boldsymbol{\eta}_{i}$ is iid. Therefore, all our results
can be considered as extensions to the inference method proposed by Conley and
Taber, (2011) for the particular case in which $H_{j}(Z_{i};\delta_{j})$ is
the identity matrix.

Ferman and Pinto, (2019) consider as a standard example the case in which
outcomes $y_{it}$ come from aggregating data from $Z_{i}$ individuals in unit
$i$ and time $t$. In this case, they show that, under a wide range of
structures on the spatial and serial correlations within unit $i$, the
variance of
$W_{i}=\frac{1}{T-t^{\ast}}\sum_{t=t^{\ast}+1}^{T}\eta_{it}-\frac{1}{t^{\ast}}\sum_{t=1}^{t^{\ast}}\eta_{it}$,
as a function of $Z_{i}$, would have a parametric form given by
$\mathbb{V}[W_{i}]=A+B\frac{1}{Z_{i}}$ for parameters $A,B\geq 0$. In Appendix
A, we show that, for this setting in which $y_{it}$ is the aggregate of
$Z_{i}$ individual-level observations, and under some homogeneity assumptions,
this parametric form generalizes for our Example 1 as $H_{j}(\cdot
Z_{i};\delta_{j})^{2}=\mathbb{V}[B_{j}A_{j}\boldsymbol{\eta}_{i}]=\delta_{0j}+\delta_{1j}\frac{1}{Z_{i}}$,
for $\delta_{0j},\delta_{1j}\geq 0$, and for our Examples 2 and 3 as
$H_{j}(Z_{i};\Lambda_{0j},\Lambda_{1j})^{2}=\Lambda_{0j}+\frac{1}{Z_{i}}\Lambda_{1j}$
for $K\times K$ positive semidefinite matrices $\Lambda_{0j}$ and
$\Lambda_{1j}$.

Assumption 3 suggests the following procedure for approximating the
distribution of $\hat{\bar{\boldsymbol{\alpha}}}-\bar{\boldsymbol{\alpha}}$.

  1. 1.

Estimate $\hat{\bar{\boldsymbol{\alpha}}}$. For each $j=1,\ldots,N_{1}$,
compute and store
$\widehat{W}_{i}(j)=B_{j}A_{j}\boldsymbol{Y}_{i}-\frac{1}{N_{0}}\sum_{i^{\prime}=N_{1}+1}^{N}B_{j}A_{j}\boldsymbol{Y}_{i^{\prime}}$
for every $i\in\\{N_{1}+1,\ldots N\\}$.

  2. 2.

Use the $\hat{W_{i}}(j)$ to construct estimators $\hat{\delta}_{j}$ of
$\delta_{j}$, for $j=1,\ldots N_{1}$.

  3. 3.

For each $j=1,\ldots N_{1}$, compute the normalized residuals
$\\{\widetilde{W}_{i}(j)\\}_{i=N_{1}+1}^{N}$, where
$\widetilde{W}_{i}(j)=H_{j}(Z_{i};\hat{\delta}_{j})^{-1}\hat{W}_{i}(j)$.

  4. 4.

Do that $B$ times:

    1. (a)

Draw a sequence of $N_{1}$ number with replacement from $\\{N_{1}+1,...,N\\}$,
$(i_{b}(1),...,i_{b}(N_{1}))$.

    2. (b)

Compute and store $\hat{e}_{b}$, where

$\hat{e}_{b}=\sum_{j=1}^{N_{1}}H_{j}(Z_{j};\hat{\delta}_{j})\widetilde{W}_{i_{b}(j)}(j)$

The next proposition provides conditions under which this approach correctly
estimates the distribution of interest. In what follows, $\lVert\cdot\rVert$
denotes the spectral norm of a matrix, and
$\lambda_{\operatorname{min}}(\cdot)$ is the smallest eigenvalue of a square
matrix. We also denote by $F$ the distribution function of
$\sum_{j=1}^{N_{1}}B_{j}A_{j}\boldsymbol{\eta}_{j}$.

###### Proposition 2

Suppose Assumptions 1, 2 and 3 hold. Suppose that, for each $j\leq N_{1}$, the
estimators $\hat{\delta}_{j}$ of $\delta_{j}$ satisfy, as $N_{0}\to\infty$,
$\max_{i\in\\{j\\}\cup\\{N_{1}+1,\ldots,N_{0}\\}}\lVert
H_{j}(Z_{i};\hat{\delta}_{j})-H_{j}(Z_{i};\delta_{j})\rVert\overset{p}{\to}0$.
Moreover, assume there exist $0<\underline{h}\leq\overline{h}<\infty$ such
that $\underline{h}\leq\min_{j\leq N_{1},i\in\\{N_{1}+1,\ldots
N\\}}\lambda_{\operatorname{min}}(H_{j}(Z_{i};\delta_{j}))$ and $\max_{j\leq
N_{1},i\in\\{j,N_{1}+1,\ldots N\\}}\lVert
H_{j}(Z_{i};\delta_{j}))\rVert\leq\overline{h}$, uniformly as
$N_{0}\to\infty$. We then have that $\hat{F}(c)\overset{p}{\to}F(c)$,
uniformly over $c\in\mathbb{R}^{K}$, where $\hat{F}(c)$ is given by:

$\hat{F}(c)=\frac{1}{N_{0}^{N_{1}}}\sum_{i_{1}={N_{1}+1}}^{N}\sum_{i_{2}={N_{1}+1}}^{N}\ldots\sum_{i_{N_{1}}={N_{1}+1}}^{N}\mathbf{1}\left\\{\sum_{j=1}^{N_{1}}H(Z_{j};\hat{\delta}_{j})\widetilde{W}_{i_{j}}(j)\leq
c\right\\}\,.$

Proof. See Appendix B.1.

###### Remark 1

Assumption 3 requires the matrices $H_{j}(Z_{i};\delta_{j})$ to be positive
definite, whereas in Examples 2 and 3, some of the entries of
$B_{j}A_{j}\boldsymbol{\eta}_{i}$ are degenerate. Our results immediately
follow in these settings by modifying Assumption 3 and the smallest eigenvalue
assumption in Proposition 2 to hold for the submatrix of
$H_{j}(Z_{i};\delta_{j})$ which stores the covariance matrix of the
nondegenerate part of $B_{j}A_{j}\boldsymbol{\eta}_{i}$. We do not state the
assumptions this way for ease of exposition.

###### Remark 2

Propososition 2 requires estimators $\hat{\delta}_{j}$ of the parametric
model. These may be obtained by solving:

$\min_{d_{j}\in\Delta_{j}}\sum_{s=N_{1}+1}^{N}\lVert\hat{W}_{s}(j)\hat{W}_{s}(j)^{\prime}-H_{j}(Z_{i};d_{j})H_{j}(Z_{i};d_{j})^{\prime}\rVert_{F}^{2},$

where $\lVert\cdot\rVert_{F}$ denotes the Froebenius norm. The resulting
estimator satisfies the assumption required in Proposition 2 under standard
conditions. Indeed, in our leading example where the parametric model comes
from aggregation from individual-level data, uniform consistency of the
estimated variances is ensured if
$\lVert\hat{\Lambda}_{0,j}-{\Lambda}_{0,j}\rVert\overset{p}{\to}0$ and
$\lVert\hat{\Lambda}_{1,j}-{\Lambda}_{1,j}\rVert\overset{p}{\to}0$ for all
$j\leq N_{1}$, which can be ensured by relying on standard arguments on least-
squares estimation.

###### Remark 3

We also note that, for our proposed algorithm to work, we require that the
estimated $\hat{H}_{j}(Z_{i},\hat{\delta}_{j})$, $i>N_{1}$, be invertible.
This occurs with probability approaching one in our setting. However, to
ensure invertibility in finite samples, it may be convenient to include a
penalisation on small values of the smallest singular value of
$\hat{H}_{j}(Z_{i},\hat{\delta}_{j})$. Provided the penalty vanishes
asymptotically, the estimator retains the required properties by Proposition
2.

Proposition 2 provides conditions under which the proposed algorithm recovers
the distribution of the test statistic. This can be used to construct
confidence intervals for the target parameters. If the parameter of interest
is unidimensional ($K=1$), an asymptotically valid $(1-\alpha)$ confidence
intervals for $\bar{\alpha}$ can be constructed as:

$[\hat{\bar{\alpha}}-\hat{q}_{1-\alpha},\hat{\bar{\alpha}}+\hat{q}_{1-\alpha}]\,,$

where $\hat{q}_{1-\alpha}$ is the $1-\alpha$ empirical quantile of
$|\hat{e}_{b}|$. Alternatively, if $K>1$,$(1-\alpha)$ uniform confidence bands
may be constructed as:

$\prod_{s=1}^{K}[\boldsymbol{\hat{\bar{\alpha}}}_{s}-\hat{\iota}_{s}\hat{q}_{1-\alpha},\boldsymbol{\hat{\bar{\alpha}}}_{s}+\hat{\iota}_{s}q_{1-\alpha}],$

where $\hat{\iota}_{s}$ are normalizing constants, and $\hat{q}_{1-\alpha}$ is
the $1-\alpha$ empirical quantile of
$\max_{s=1,\ldots,K}|\hat{e}_{b,s}/\hat{\iota}_{s}|$. Typical choices include
$\hat{\iota}_{s}=1$ for all $s$, which leads to a constant-width uniform band;
and $\hat{\iota}_{s}=\hat{\sigma}_{s}$, where $\hat{\sigma}_{s}^{2}$ is the
estimator of the variance of $\hat{\bar{\alpha}}_{s}$ using the simulated
draws.

We summarise the discussion in the corollary below.

###### Corollary 1

Suppose Assumptions 1, 2, 3 and the conditions in Proposition 2 hold. Let
$\mathcal{C}$ be a confidence interval constructed as:

$\mathcal{C}=\prod_{s=1}^{K}[\boldsymbol{\hat{\bar{\alpha}}}_{s}-\hat{\iota}_{s}\hat{q}_{1-\alpha},\boldsymbol{\hat{\bar{\alpha}}}_{s}+\hat{\iota}_{s}q_{1-\alpha}],$

where the $\iota_{s}$ are such that
$\hat{\iota}_{s}\overset{p}{\to}\iota_{s}>0$ as $B,N_{0}\to\infty$. We then
have that:

$\lim_{N_{0},B\to\infty}\mathbb{P}[\boldsymbol{{\bar{\alpha}}}\in\mathcal{C}]=1-\alpha.$

Specifically, if $\hat{\iota}_{s}=\hat{\sigma_{s}}$, then
$\hat{\sigma}_{s}^{2}\overset{p}{\to}\sigma_{s}^{2}$ as $B,N_{0}\to\infty$.

Proof. See Appendix B.2.

## References

  * Abadie et al., (2010) Abadie, A., Diamond, A., and Hainmueller, J. (2010).  Synthetic Control Methods for Comparative Case Studies: Estimating the Effect of California’s Tobacco Control Program.  Journal of the American Statiscal Association, 105(490):493–505.
  * Arkhangelsky et al., (2021) Arkhangelsky, D., Athey, S., Hirshberg, D. A., Imbens, G. W., and Wager, S. (2021).  Synthetic difference-in-differences.  American Economic Review, 111(12):4088–4118.
  * Ben-Michael et al., (2021) Ben-Michael, E., Feller, A., and Rothstein, J. (2021).  The augmented synthetic control method.  Journal of the American Statistical Association, 116(536):1789–1803.
  * Borusyak et al., (2021) Borusyak, K., Jaravel, X., and Spiess, J. (2021).  Revisiting event study designs: Robust and efficient estimation.
  * Callaway and Sant’Anna, (2021) Callaway, B. and Sant’Anna, P. H. (2021).  Difference-in-differences with multiple time periods.  Journal of Econometrics, 225(2):200–230.  Themed Issue: Treatment Effect 1.
  * Chernozhukov et al., (2021) Chernozhukov, V., Wüthrich, K., and Zhu, Y. (2021).  An exact and robust conformal inference method for counterfactual and synthetic controls.  Journal of the American Statistical Association, 116(536):1849–1864.
  * Conley and Taber, (2011) Conley, T. G. and Taber, C. R. (2011).  Inference with Difference in Differences with a Small Number of Policy Changes.  The Review of Economics and Statistics, 93(1):113–125.
  * (8) de Chaisemartin, C. and D’Haultfoeuille, X. (2020a).  Difference-in-differences estimators of intertemporal treatment effects.
  * (9) de Chaisemartin, C. and D’Haultfoeuille, X. (2020b).  Two-way fixed effects estimators with heterogeneous treatment effects.  American Economic Review, 110(9):2964–96.
  * Donald and Lang, (2007) Donald, S. G. and Lang, K. (2007).  Inference with Difference-in-Differences and Other Panel Data.  The Review of Economics and Statistics, 89(2):221–233.
  * Ferman, (2019) Ferman, B. (2019).  Assessing Inference Methods.  arXiv e-prints, page arXiv:1912.08772.
  * Ferman, (2020) Ferman, B. (2020).  Inference in differences-in-differences with few treated units and spatial correlation.
  * Ferman, (2021) Ferman, B. (2021).  On the properties of the synthetic control estimator with many periods and many controls.  Journal of the American Statistical Association, 116(536):1764–1772.
  * Ferman and Pinto, (2019) Ferman, B. and Pinto, C. (2019).  Inference in differences-in-differences with few treated groups and heteroskedasticity.  The Review of Economics and Statistics, 0(ja):null.
  * Ferman and Pinto, (2021) Ferman, B. and Pinto, C. (2021).  Synthetic controls with imperfect pretreatment fit.  Quantitative Economics, 0(ja):null.
  * Freyaldenhoven et al., (2021) Freyaldenhoven, S., Hansen, C., Pérez, J. P., and Shapiro, J. M. (2021).  Visualization, identification, and estimation in the linear panel event-study design.  Technical report, National Bureau of Economic Research.
  * Goodman-Bacon, (2021) Goodman-Bacon, A. (2021).  Difference-in-differences with variation in treatment timing.  Journal of Econometrics, 225(2):254–277.  Themed Issue: Treatment Effect 1.
  * Hagemann, (2020) Hagemann, A. (2020).  Inference with a single treated cluster.
  * Lehmann and Romano, (2005) Lehmann, E. L. and Romano, J. P. (2005).  Testing Statistical Hypotheses.  Springer New York.
  * MacKinnon and Webb, (2019) MacKinnon, J. G. and Webb, M. D. (2019).  Randomization Inference for Difference-in-Differences with Few Treated Clusters.  Journal of Econometrics, Forthcoming.
  * Montiel Olea and Plagborg-Moller, (2019) Montiel Olea, J. L. and Plagborg-Moller, M. (2019).  Simultaneous confidence bands: Theory, implementation, and an application to svars.  Journal of Applied Econometrics, 34(1):1–17.
  * Sun and Abraham, (2020) Sun, L. and Abraham, S. (2020).  Estimating dynamic treatment effects in event studies with heterogeneous treatment effects.  Journal of Econometrics.
  * van der Vaart, (1998) van der Vaart, A. W. (1998).  Asymptotic Statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.

##  Appendix A Parametric models for examples

###  A.1 Panel data

Suppose the $y_{it}$ are obtained from aggregating individual-level data from
$Z_{i}$ individuals. In this case, we may write:

$\boldsymbol{\eta_{i}}=\boldsymbol{f}_{i}+\frac{1}{Z_{i}}\sum_{s=1}^{Z_{i}}\boldsymbol{\epsilon}_{s,i},$

where $\boldsymbol{f}_{i}$ are group shocks, $\boldsymbol{\epsilon}_{s,i}$ are
idiosyncratic shocks, and, without loss of generality,
$\operatorname{cov}(\boldsymbol{f}_{i},\boldsymbol{\epsilon}_{s,i})=0$. If we
assume $\boldsymbol{\eta}_{i}$ to be iid across $i$ and
$\boldsymbol{\eta}_{s,i}$ to be iid across $s$ and $i$, we have that:

$\mathbb{V}[\boldsymbol{\eta}_{i}]=V_{0}+\frac{V_{1}}{Z_{i}}$

For positive semidefinite matrices $V_{0}$ and $V_{1}$. For given choices of
$A_{j}$ and $B_{j}$, we thus obtain the model:

$H_{j}(Z_{i};\delta_{j})^{2}=\Lambda_{0j}+\frac{\Lambda_{1j}}{Z_{i}},$

where $\Lambda_{0j}$ and $\Lambda_{1j}$ are semidefinite $K\times K$ matrices.

###  A.2 Repeated cross-sections

The model in the previous section is especially suited for panel data, since
it allows for serial correlation in idiosyncratic shocks. In the case where
the data $y_{jt}$ is constructed from repeated cross-sections with $Z_{jt}$
individuals, a more parsimonious model can be obtained. Indeed, in this case,
we consider the model:

$\boldsymbol{\eta}_{i}=\boldsymbol{f}_{i}+\begin{bmatrix}\frac{1}{Z_{1}}\sum_{s=1}^{Z_{1}}\epsilon_{s,i,1}\\\
\vdots\\\ \ \frac{1}{Z_{T}}\sum_{s=1}^{Z_{T}}\epsilon_{s,i,T}\end{bmatrix}\,.$

where we assume that $\boldsymbol{f}_{i}$ is iid across $i$; for each $t$,
$\epsilon_{s,i,t}$ is iid across $s$ and $i$; and that $\epsilon_{s,i,t}$
$\epsilon_{s^{\prime},i^{\prime},t^{\prime}}$ are independent for $t\neq
t^{\prime}$. This leads to:

$\mathbb{V}[\boldsymbol{\eta}_{i}]=V_{0}+\operatorname{diag}(\sigma^{2}_{1}/Z_{i1},\sigma^{2}_{2}/Z_{i2},\ldots,\sigma^{2}_{T}/Z_{iT}),$

In Example 2, this model leads to the following parametrization:

$H_{j}(Z_{j};\delta_{j})^{2}=\Lambda_{0j}+\frac{\omega_{0j}}{Z_{i,t_{j}^{*}}}\begin{bmatrix}\boldsymbol{1}_{K_{j}\times
K_{j}}&\boldsymbol{0}_{K_{j}\times K-K_{j}}\\\ \boldsymbol{0}_{K-K_{j}\times
K_{j}}&\boldsymbol{0}_{K-K_{j}\times
K-K_{j}}\end{bmatrix}+\operatorname{diag}\left(\frac{\omega_{1j}}{Z_{i,t_{j}^{*}+1}},\frac{\omega_{2j}}{Z_{i,t_{j}^{*}+2}},\ldots,\frac{\omega_{K_{j},j}}{Z_{i,K_{j}}},0,0,\ldots,0\right)$

##  Appendix B Proofs of main results

###  B.1 Proof of Proposition 2

First, we note that, by Weyl’s inequality and the eigenvalue and spectral norm
assumptions in the statement of the theorem,
$\mathbb{P}[\cap_{j=1}^{N_{1}}\cap_{i=N_{1}+1}^{N}\\{H_{j}(Z_{i};\hat{\delta}_{j})^{-1}\text{
exists}\\}]\to 1$. Without loss, we assume these matrices to be always
invertible.101010Otherwise, premultiply random variables that depend on an
inverse by an indicator of the event where all matrices are invertible. The
following argument remains essentially unchanged, except for additional
$o_{p}(1)$ terms appearing in the derivations.  Let:

$\tilde{\boldsymbol{e}}=\sum_{j=1}^{N_{1}}H_{j}(Z_{j};\hat{\delta}_{j})\tilde{W}_{i^{*}_{j}}(j),$

where the $i^{*}_{j}$ are iid draws from
$\operatorname{Uniform}(\\{N_{1}+1,\ldots,N\\})$, independently from
$\\{\boldsymbol{\eta}_{i}\\}_{i}$. Do further define

${\boldsymbol{e}}=\sum_{j=1}^{N_{1}}H_{j}(Z_{j};\delta_{j})\xi_{i^{*}_{j},j}.$

We show $\lVert\hat{\boldsymbol{e}}-\boldsymbol{e}\rVert=o_{p}(1)$. Note that:

$\displaystyle\lVert\hat{\boldsymbol{e}}-\boldsymbol{e}\rVert\leq\sum_{j=1}^{N_{1}}\lVert
H_{j}(Z_{j};\hat{\delta}_{j})-H_{j}(Z_{j};{\delta}_{j})\rVert\lVert\tilde{W}_{i^{*}_{j}}(j)-{\xi}_{i^{*}_{j},j}\rVert+\sum_{j=1}^{N_{1}}\lVert
H_{j}(Z_{j};\hat{\delta}_{j})-H_{j}(Z_{j};{\delta}_{j})\rVert\lVert{\xi}_{i^{*}_{j},j}\rVert+$
(7) $\displaystyle\sum_{j=1}^{N_{1}}\lVert
H_{j}(Z_{j};{\delta}_{j})\rVert\lVert
H_{j}(Z_{j};\hat{\delta}_{j})-H_{j}(Z_{j};{\delta}_{j})\rVert$

Fix $j\leq N_{1}$. First, by the Assumption in the statement of the
Proposition, $\lVert
H_{j}(Z_{j};\hat{\delta}_{j})-H_{j}(Z_{j};{\delta}_{j})\rVert=o_{p}(1)$.
Moreover, by Assumption 3, $\xi_{i_{j}^{*},j}\overset{d}{=}\xi_{j,j}$,
implying $\lVert\xi_{i_{j}^{*},j}\rVert=O_{p}(1)$. Finally, letting
$\mathbb{E}_{*}$ denote expectations with respect to the distribution of
$(i_{1}^{*},\ldots,i_{N_{1}}^{*})$, with $\\{\boldsymbol{\eta}_{i}\\}$ fixed,
we observe that:

$\displaystyle\mathbb{E}_{*}\lVert\tilde{W}_{i^{*}_{j}}(j)-{\xi}_{i^{*}_{j},j}\rVert=\frac{1}{N_{0}}\sum_{i=N_{1}+1}^{N}\lVert\tilde{W}_{i}(j)-{\xi}_{i,j}\rVert\leq$
(8) $\displaystyle\max_{i\in\\{N_{1}+1,\ldots N\\}}\lVert
H_{j}(Z_{i};\delta_{j})^{-1}-H_{j}(Z_{i};\hat{\delta}_{j})^{-1}\rVert\left(\left\lVert\frac{1}{N_{0}}\sum_{i=N_{1}+1}^{N}B_{j}A_{j}\boldsymbol{\eta}_{i}\right\rVert+\frac{1}{N_{0}}\sum_{i=N_{1}+1}^{N}\left\lVert
B_{j}A_{j}\boldsymbol{\eta}_{i}\right\rVert\right)+$
$\displaystyle\max_{i\in\\{N_{1}+1,\ldots N\\}}\lVert
H_{j}(Z_{i};\delta_{j})^{-1}\rVert\left\lVert\frac{1}{N_{0}}\sum_{i=N_{1}+1}^{N}B_{j}A_{j}\boldsymbol{\eta}_{i}\right\rVert=o_{p}(1)$

where we use the law of large numbers, the fact that, for a symmetric
invertible matrix $A$, $\lVert
A^{-1}\rVert=1/\lambda_{\operatorname{min}}(A)^{2}$; and that for invertible
matrices $A$ and $B$,

$\lVert A^{-1}-B^{-1}\rVert\leq\left(\frac{1}{\lVert
B\rVert}\right)\frac{\lVert A-B\rVert}{\lVert B\rVert-\lVert A-B\rVert}\,.$

By the conditional Markov inequality,
$\mathbb{P}_{*}[\lVert\tilde{W}_{i^{*}_{j}}(j)-{\xi}_{i^{*}_{j},j}\rVert>\epsilon]=o_{p}(1)$
for every $\epsilon>0$. Iterated expectations and the bounded convergence
theorem imply
$\lVert\tilde{W}_{i^{*}_{j}}(j)-{\xi}_{i^{*}_{j},j}\rVert=o_{p}(1)$. We have
thus shown that, for fixed $j\leq N_{1}$, the sum of terms on the left-hand
side of (7) associated with such $j$ is $o_{p}(1)$. Since the chosen $j$ was
arbitrary, we conclude the whole term is $o_{p}(1)$, which proves that
$\lVert\tilde{\boldsymbol{e}}-\boldsymbol{e}\rVert=o_{p}(1)$.

Next, define

$\tilde{\boldsymbol{g}}=\sum_{j=1}^{N_{1}}H_{j}(Z_{j};\hat{\delta}_{j})\tilde{W}_{s^{*}_{j}}(j),$

where $s^{*}_{j}$ are iid draws from
$\operatorname{Uniform}(N_{1}+1,\ldots,N)$, independently from both the
$i^{*}_{j}$ and the $\boldsymbol{\eta}_{i}$. Similarly, let:

$\boldsymbol{g}=\sum_{j=1}^{N_{1}}\xi_{s^{*}_{j},j}.$

Note that $\lVert\tilde{\boldsymbol{g}}-\boldsymbol{g}\rVert=o_{p}(1)$. Next,
we claim that
$(\boldsymbol{e}^{\prime},\boldsymbol{g}^{\prime})^{\prime}\overset{d}{\to}F\otimes
F$. To see this, observe that the event $E=\\{\exists j\neq
j^{\prime}:i_{j}^{*}=i_{j^{\prime}}^{*}\text{ or
}s_{j}^{*}=s_{j^{\prime}}^{*}\\}\cup\\{\exists
j,j^{\prime}:i_{j}^{*}=s_{j^{\prime}}^{*}\\}$ is such that $\mathbb{P}[E]\to
0$. We thus have that, for any $c_{1},c_{2}\in\mathbb{R}^{K}$:

$\displaystyle\mathbb{P}[\boldsymbol{e}\leq c_{1},\boldsymbol{g}\leq
c_{2}]=F(c)^{2}\mathbb{P}[E^{c}]+\mathbb{P}[\\{\boldsymbol{e}\leq
c_{1}\\}\cap\\{\boldsymbol{g}\leq c_{2}\\}\cap E]=F(c)^{2}+o(1)\,.$

Application of Slutsky lemma then yields
$(\tilde{\boldsymbol{e}}^{\prime},\tilde{\boldsymbol{g}}^{\prime})^{\prime}\overset{d}{\to}F\otimes
F$. It then follows from Theorem 15.2.3 of Lehmann and Romano, (2005) that
$\hat{F}(c)\overset{p}{\to}F(c)$ for every $c\in\mathbb{R}^{K}$. Uniform
convergence follows from the fact that $F$ is continuous (van der Vaart,,
1998, p. 339).

###  B.2 Proof of Corollary 1

Let $\hat{F}_{B}$ denote the empirical distribution obtained from running the
proposed algorithm $B$ times. Denoting by $\mathbb{E}_{*}$ the expectation
with respect to the randomness of the algorithm and with the data fixed, we
have that, for each $c\in\mathbb{R}^{K}$:

$\mathbb{E}_{*}[\hat{F}_{B}(c)]=\hat{F}(c)\,,$

$\mathbb{V}_{*}[\hat{F}_{B}(c)]=\frac{\hat{F}(c)(1-\hat{F}(c))}{B}\,.$

It then follows from Proposition 2 that, for each $\epsilon>0$:

$\lim_{N_{0},B\to\infty}\mathbb{P}_{*}[|\hat{F}_{B}(c)-F(c)|>\epsilon]=0.$

The bounded convergence theorem and the fact that $F$ is continuous thus imply
that:

$\sup_{c\in\mathbb{R}^{K}}|\hat{F}_{B}(c)-F(c)|=o_{p}(1)$

Let $\hat{\iota}=\operatorname{diag}(\hat{\iota}_{1},\ldots\hat{\iota}_{s})$.
Since $F$ is continuous and $\hat{\iota}\overset{p}{\to}\iota$, we have that:

$\sup_{c\in\mathbb{R}^{K}}|\hat{F}_{B}(\hat{\iota}c)-F(\iota c)|=o_{p}(1),$

and it follows from the continuous mapping theorem that:

$\sup_{c\in\mathbb{R}}|\hat{G}_{B}(c)-G(c)|=o_{p}(1)\,,$

where $\hat{G}_{B}$ is the empirical distribution of
$\max_{s=1,\ldots,K}|\hat{e}_{s,b}|/|\hat{\iota}_{s}|$, and
$G(c)=\int_{\boldsymbol{[-c,c]^{K}}}F(\iota dx)$. By Lemma 21.2 of van der
Vaart, (1998), it follows that $\hat{q}_{1-\alpha}$ converges in proability to
the $(1-\alpha)$ quantile of $G$. Sjnce by Slutsky theorem
${\hat{\iota}}^{-1}(\boldsymbol{\hat{\bar{\alpha}}-\boldsymbol{\bar{\alpha}}})$
converges weakly to $F(\iota c)$, we have, by the continuous mapping theorem:

$\displaystyle\mathbb{P}[\boldsymbol{{\bar{\alpha}}}\in\mathcal{C}]=\mathbb{P}[\mathbf{1}\hat{q}_{1-\alpha}\leq{\hat{\iota}}^{-1}(\boldsymbol{\hat{\bar{\alpha}}}-\boldsymbol{\bar{\alpha}})\leq\mathbf{1}\hat{q}_{1-\alpha}]=$
(9)
$\displaystyle\mathbb{P}[{\hat{\iota}}^{-1}(\boldsymbol{\hat{\bar{\alpha}}}-\boldsymbol{\bar{\alpha}})\leq\mathbf{1}\hat{q}_{1-\alpha}]-\sum_{\boldsymbol{\kappa}\in\\{-1,1\\}^{K}:\exists\boldsymbol{\kappa}_{s}=1}\mathbb{P}[\operatorname{diag}(\boldsymbol{\kappa}){\hat{\iota}}^{-1}(\boldsymbol{\hat{\bar{\alpha}}}-\boldsymbol{\bar{\alpha}})\leq-\operatorname{diag}(\boldsymbol{\kappa})\hat{q}_{1-\alpha}]\to$
$\displaystyle\int_{[-q_{1-\alpha}\boldsymbol{1},q_{1-\alpha}\boldsymbol{1}]}F(\iota
dx)=1-\alpha\,.$

The last assertion of the corollary is immediate and therefore not proved.