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.
49,546 characters
Synthetic Principal Component Design: Fast Covariate Balancing with Synthetic Controls
\begin{abstract}
The optimal design of experiments typically involves solving an NP-hard combinatorial optimization problem. In this paper, we aim to develop a globally convergent and practically efficient optimization algorithm. Specifically, we consider a setting where the pre-treatment outcome data is available and the synthetic control estimator is invoked. The average treatment effect is estimated via the difference between the weighted average outcomes of the treated and control units, where the weights are learned from the observed data. {Under this setting, we surprisingly observed that the optimal experimental design problem could be reduced to a so-called \textit{phase synchronization} problem.} We solve this problem via a normalized variant of the generalized power method with spectral initialization. On the theoretical side, we establish the first global optimality guarantee for experiment design when pre-treatment data is sampled from certain data-generating processes. Empirically, we conduct extensive experiments to demonstrate the effectiveness of our method on both the US Bureau of Labor Statistics and the Abadie-Diemond-Hainmueller California Smoking Data. In terms of the root mean square error, our algorithm surpasses the random design by a large margin.
\end{abstract}
\section{Introduction}
Estimating the average effects of a binary treatment is one of the main goals of empirical economic and political studies. Randomization in controlled trials is one of the golden rules for estimating average treatment effects (ATE). Suppose the treatment assignment procedure guarantees that the potential outcomes are independent of the treatment status. In that case, a simple difference-in-mean (i.e., average outcomes of the treated and control units) estimator becomes an unbiased estimator of ATE. Nevertheless, a fully randomized experiment may be affected by a significantly high variance in the final estimation. Such variance can be reduced via exploiting the feature information in the observed data. Naturally, we focus on the following question: \emph{can the observed covariates improve the statistical properties of the ATE estimators via experimental designing?} \citep{rubin2008objective,kasy2016experimenters}. This problem is referred to as \emph{covariate balancing}, which restricts the randomization to achieve covariate balance between treatment groups \citep{efron1971forcing,morgan2012rerandomization}.
Covariate balancing has been substantially explored in the literature.
One approach is applying covariate balancing via an adequately designed propensity score~\citep{imai2014covariate,zhao2019covariate}, which requires that we have access to a reasonably large sample of experimental units. However, large samples of experimental units are not always available in practice. Therefore, if only a moderate number of units are accessible to the treatment, \citep{bansal2018gram} solve an NP-hard combinatorial optimization problem to balance the empirical covariates. To address this issue, we aim to design a computationally feasible optimal design of experimental studies.
We specifically consider the optimal design of experiments for the synthetic control estimator as in ~\citep{abadie2003economic,abadie2010synthetic,abadie2015comparative}, which becomes an attractive estimation procedure when only a small number of units can be exposed to the experiment. That is, the experiment designer can observe the pre-treatment panel outcome data for a number of units in a number of time periods. Synthetic control compares treated units with a weighted average of untreated units. The weights are determined via empirical fit on the observed pre-treatment outcome. As an example, consider the study of Proposition 99's effect presented in \citep{abadie2010synthetic}. Proposition 99 is a large-scale anti-smoking legislation program that California implemented in 1988. The policy maker wants to estimate the effect of this piece of legislation. Synthetic control suggests the policy maker to estimate the counterfactual outcomes after 1988 by using the observed outcomes of the states without legislative restrictions. \citep{abadie2010synthetic} produces the synthetic control by a combination of Colorado, Connecticut, Montana, Nevada, and Utah and find out that annual per-capita cigarette sales in California were about 26 packs lower than what they would have been by the year 2000. Beyond this application, the method of synthetic controls has been used in many other empirical policy evaluation problems, including legalized prostitution \citep{cunningham2018decriminalizing}, corporate political connections \citep{acemoglu2016value}, taxation \citep{kleven2013taxation}, to just name a few.
To realize the benefits offered by synthetic control, \citep{doudchenko2021synthetic,abadie2021synthetic} considers how an optimally designed experiment can help the experiment designer to further reduce the variance. \citep{doudchenko2021synthetic,abadie2021synthetic} propose an optimization approach to select the control group based on the observed pre-treatment outcome. Namely, the choice of the treated units aims to balance the weighted average of treated and untreated covariates.
As such, the designer can choose the best non-negative weights. \citep{doudchenko2021synthetic} proves that the underlying optimization problem is still NP-hard and \citep{abadie2021synthetic} relaxes the optimization problem into a canonical Quadratic Constraint Quadratic Program
(QCQP). Nevertheless, the resulting QCQP is rather computationally demanding and applicable algorithms (e.g., SDR(cf. Semidefinite Programming)~\citep{luo2010semidefinite}) are not guaranteed to reach a global optimum.
In this paper, we aim to design the first globally convergent optimization algorithm for the weighted covariate balancing formulation \citep{doudchenko2021synthetic}. Moreover, the algorithm is practically efficient. We achieve this by first removing the so-called units cardinality constraint being treated in \citep{doudchenko2021synthetic}. Surprisingly, we find out that this relaxation can be shown to be equivalent to a phase synchronization problem \citep{singer2011angular}. Although phase synchronization is still an NP-hard problem \citep{zhang2006complex}, many practical algorithms have been recently developed for this nonconvex problem, \citep{bandeira2017tightness,boumal2016nonconvex,liu2017estimation,zhong2018near}. Moreover, as we will argue, this nonconvex problem is polynomial-time solvable under a suitable data generating process (\emph{i.e.} average case complexity). Motivated by this line of research, we propose \emph{Synthetic principal component Design} (SPCD), which optimizes the treatment decision via (a normalized variate of) the generalized power method with spectral initialization \citep{chen2021spectral}. If we assume that the pre-trement data is sampled from the linear fixed-effect model studied in \citep{abadie2010synthetic,xu2017generalized,athey2021matrix,ferman2021properties} and invoke the realizable assumptions considered in \citep{abadie2021synthetic}, we can further establish its \emph{global} optimization guarantee and statistical estimation guarantee for the proposed synthetic control procedure. To the best of our knowledge, this is the first computational paradigm of combinatorial optimization-based experiment design which enjoys a global optimization guarantee.
\paragraph{Paper organization} We organize our paper as follows. In Section \ref{section:setting}, we introduce the setup in \citep{doudchenko2021synthetic,abadie2021synthetic} where the authors consider
covariate balancing under the synthetic control setting and reformulate it as a phase
synchronization problem. In Section \ref{section:GIPW}, we introduce the generalized power
method and provide global a optimization guarantee under the linear factor model
\citep{abadie2003economic,xu2017generalized,athey2021matrix,ferman2021properties} with a realizable
assumption \citep{abadie2021synthetic}. In Section \ref{section:numerical}, we apply our method to
both simulated and real-world data sets. We end with some closing remarks in Section \ref{section:conl}.
\vspace{-0.1in}
\subsection{Related Work}
\paragraph{Synthetic Control} Synthetic control \citep{abadie2003economic,abadie2010synthetic} is one of the leading methods to estimate the causal effect of a binary treatment in the panel data setting. It constructs a weighted combination of untreated groups used as controls, to which the outcome of the treatment group is compared. This construction of the control group significantly improves the performance when there are limited units that can be exposed to the experiment designer.
Statistical consistency and inference properties have been established for the synthetic control under the linear factor models \citep{xu2017generalized,athey2021matrix,ferman2021properties}, which is also the starting point of our theory.
\paragraph{Covariate Balancing} To make the difference-in-mean estimator more precise, experimenters sometimes restrict the treatment assignment to achieve covariate balance between treatment groups, \emph{i.e.,}
\[
\min_{\{D_i\}_{i=1}^n} \left |\left|\frac{1}{\sum
\limits_{i}D_i}\sum_{i:D_i = 1}X_{i} -\frac{1}{\sum
\limits_{i}(1-D_i)}\sum_{i:D_i = 0}X_{i}\right |\right|^2
\]
where $\{X_i\in\mathbb{R}^n\}_{i=1}^n$ are the observed units and $\{D_i\in\{0,1\}\}_{i=1}^n$
denotes the treatment experiments the designer aims to optimize
\citep{efron1971forcing,morgan2012rerandomization,imai2014covariate,harshaw2019balancing}. Indeed,
this is a 0-1 NP-hard combinatorial optimization problem. Various different methods to approximately
solve the problem have been proposed, including rerandomization
\citep{morgan2012rerandomization,kallus2018optimal,kallus2021optimality,li2018asymptotic}, design
particular propensity score \citep{imai2014covariate,zhao2019covariate} and the recently proposed
Gram-Schmidt walk \citep{bansal2019algorithm,harshaw2019balancing}. In this paper, we follow the setting in
\citep{doudchenko2021synthetic,abadie2021synthetic} which extend the covariate balancing problem to the
synthetic control setting. \citep{doudchenko2021synthetic} has proved that this problem is still
NP-hard and \citep{abadie2021synthetic} further relaxes the optimization problem into a canonical
form of QCQP.
\paragraph{Phase Synchronization} Phase synchronization aims to estimate $n$ unknown angles $\theta_i\in [0,2\pi],i\in[N]$ through noisy measurements of their offset $\theta_i-\theta_j$ mod $2\pi$ \citep{singer2011angular}, which has received intense interests in areas such as time synchronization between distributed networks, \citep{giridhar2006distributed}; ranking \citep{cucuringu2016sync}; computer vision \citep{wang2013exact,martinec2007robust}; and optics and inverse problems \citep{rubinstein2001reconstruction,alexeev2014phase,singer2011three}.
To globally address the non-convex optimization problem arising from the inverse problem, \citep{bandeira2017tightness} provides the first global results for the SDR under certain data generating procedures. In this paper, we follow a line of research which takes advantage of generalized power methods \citep{boumal2016nonconvex,liu2017estimation,zhong2018near} to solve the problem in a way that shares a similar global guarantee to that of the SDP relaxation, but is computationally much more amenable. It is worth noting that, although the optimization problem considered in our paper is the same as phase synchronization, the data generating process is different from the considered earlier and this leads to our problem-specific generalized power method.
\subsection{ Main Contributions}
\begin{itemize}
\item We show a surprisingly equivalence between the experiment design with synthetic control \citep{doudchenko2021synthetic,abadie2021synthetic} and phase synchronization problem \citep{singer2011angular,boumal2016nonconvex,liu2017estimation}, where separating the experiment and control group can be transformed to finding the phase of a complex signal. We further reveal the hidden connection between covariate balancing with the smallest eigenvector of the gram matrix
and propose a spectral method for fast experiment design.
\item We propose a novel normalized version of the generalized power method, which enjoys \emph{global} convergence results under certain generative models. The normalization technique can weak the condition for generative model assumptions to guarantee the global optimiality and also consistently improves the empirical results.
\item In terms of the root mean square error, our method empirically surpasses random design by a large margin on both synthetic and real-world datasets. Our performance can even exceed 500000 times of rerandomization over the Abadie-Diamond-Hainmueller smoking legislation data.
\end{itemize}
\section{Mathematical Formulations}
\label{section:setting}
In this section, we introduce the mathematical formulation of synthetic control (SC) and the optimal experiment design problem studied in \citep{doudchenko2021synthetic}.
\paragraph{Problem setup} We aim to estimate the effect of a binary treatment under the panel data setting. Researchers have access to the outcome metric of interest $Y \in\mathbb{R}^{N\times T}$ for $N$ units during $T$ time periods. At time $T$, researchers are required to execute an experiment by assigning a binary treatment described by $D_i\in\{-1,1\}, i\in [N]$ based on the observed pre-treatment data.
If $D_i=1$, then a treatment needs to be applied to unit $i$. After the treatment experiment, furthermore outcomes are observed for additional $S$ time periods $t= T+1,\cdots, T+S$. During this period, every unit $i\in[N]$ in each time period $t$ is associated with the following two random outcomes: $Y_{it}(-1)=\mu_{it}+e_{it}, \text{ and } Y_{it}(1)=Y_{it}(-1)+\tau_i,$ where $\mu_{it}$ is the base outcome, ${\tau}$ is the treatment effect aiming to estimate and $e_{it}$ is the zero mean i.i.d idiosyncratic noise with variance $\text{Var}(\epsilon_{it})=\sigma$. Once treatment $D_i$ is applied, the experimenter is able to realize $Y_{it}=\frac{(D_i+1)}{2}Y_{it}(1)+\frac{(1-D_i)}{2}Y_{it}(-1)$.
Estimating the treatment effect $\tau$ is quite challenging because once we implement a treatment on unit $j$ (\emph{i.e.,} $D_j=1$) and observe the outcome $Y_{j,T+1}(1)$, then counterfactual outcome $Y_{j,T+1}(-1)$ is not observable. With the pre-treatment observation $Y_{iT}$, synthetic control literature \citep{abadie2010synthetic,xu2017generalized} constructs the counterfactual estimate for a treated unit $j$ from a weighted average of other units' outcomes: $\hat Y_{j,T+1}(-1)=\sum_{i:D_i=-1} w_iY_{i,T+1}$. The weights $w_i$ are learned from the pre-treatment observed data via minimizing $\sum_{t=1}^T(Y_{jt}-\sum_{i:D_i=0}w_iY_{it})^2$. Then the treatment effect of unit $j$ we estimate can be written as $\tau_j=Y_{j,T+1}-\hat Y_{j,T+1}(-1)$.
\subsection{Synthetic Design}
In this subsection, we consider the synthetic design objective function proposed for two-way fixed effect in \citep{doudchenko2021synthetic,abadie2021synthetic} and reveal its hidden connection with the phase synchronization problem. \citep{doudchenko2021synthetic,abadie2021synthetic} aim to design treatment assignments $\{D_i=\pm 1\}_{i=1}^N$ and weights $\{w_i\ge0\}_{i=1}^N$ for outcome experiments at time $T+1$. For we aims to estimate a two-way fixed effect where the treatment effects are homogeneous, we can consider the \emph{weighted average treatment effect on the treated (wATET)} $\tau=\sum_{i=1}^N \frac{D_i+1}{2}w_i\tau_i$ instead \citep{bottmer2021design}. wATET can be estimated as a difference in weighted means estimator $\hat\tau=\sum_{i:D_i=1}w_i Y_{i,T+1}-\sum_{i:D_i=-1}w_i Y_{i,T+1}$ with $\sum_{i:D_i=1}w_i=\sum_{i:D_i=-1}w_i=1$. Following \citep{doudchenko2021synthetic}, upon the outcome model, the mean squared error of the difference-in-weighted-means estimator admits the decomposition
$$
\mathbb{E}\left[(\hat \tau-\tau)^2|\{D_i,w_i\}_{i=1}^N\right]=\underbrace{\left(\sum_{i:D_i=1} w_i\mu_{i,T+1}-\sum_{i:D_i=-1}w_i\mu_{i,T+1}\right)^2}_{\text{weighted covariate balancing}}+\sigma\sum_{i=1}^N w_i^2.
$$
The designer aims to design the experiment with a lowest expected mean square error. Thus, \citep{doudchenko2021synthetic} proposed the following mixed-integer programming for experimenting design with Synthetic Control:
\begin{equation}
\begin{aligned}
& \min_{\{D_i,w_i\}_{i=1}^n} \, \frac{1}{T}\sum_{t=1}^T\left(\sum_{i:D_i=1} w_iY_{it}-\sum_{i:D_i=-1}w_iY_{it}\right)^2+\sigma \sum_{i=1}^Nw_i^2\\
& \quad \, \operatorname{s.t.} \quad \,\,\,\, w_i\ge 0, D_i\in\{-1,1\}, \,\, \forall i \in [N], \\
& \quad \quad \quad \quad \sum_{i:D_i=1} w_i=\sum_{i:D_i=-1} w_i=1.
\end{aligned}
\label{eq:originalSD}
\end{equation}
\begin{remark} We remove the constraint $\sum_{i:D_i=1} D_i=K$ for a given integer $K\in\mathbb{N}$ in the mixed-integer programming in \citep{doudchenko2021synthetic} mainly as this constraints is empirically proved not critical in \citep{abadie2021synthetic}. The NP-hard proof in \citep{doudchenko2021synthetic} depends on the constraint $\sum_{i:D_i=1} D_i=K$. In the following discussion, we will show that the problem is also NP-hard even $\sum_{i:D_i=1} D_i=K$ is removed, as the resulting optimization problem can be reformulated as the $\ell_1$-PCA \citep{mccoy2011two,wang2021linear} and phase Synchronization \citep{singer2011angular,boumal2016nonconvex} problems.
\end{remark}
By making a further simplification of the problem \eqref{eq:originalSD}, we introduce a change of
variable $W_i=w_iD_i$. For $w_i\ge 0$, then $D_i=\textnormal{sgn}(W_i)$ and $w_i=|W_i|$. At the same
time, the constraint $\sum_{i:D_i=1} w_i=\sum_{i:D_i=-1} w_i=1$ is equivalent to $\mathbbm{1}^\top
W=0$ and the objective function
\[
\frac{1}{T}\sum_{t=1}^T\left(\sum_{i:D_i=1} w_iY_{it}-\sum_{i:D_i=-1}w_iY_{it}\right)^2+\sigma \sum_{i=1}^Nw_i^2
\]
can be reformulated as $W^\top (YY^\top+\lambda I) W$, where $W=[w_1,\cdots,w_N]^\top$ and
$\mathbbm{1}\in\mathbb{R}^N$ is the all one vector. Thus, \eqref{eq:originalSD} can be
recast into
\begin{equation}
\label{eq:min_l1}
\begin{aligned}
\min_{W\in\mathbb{R}^n, \mathbbm{1}^\top W=0, ||W||_1=1 } W^\top (YY^\top+\sigma I) W.
\end{aligned}
\end{equation}
Although the reformulation \eqref{eq:min_l1} translates the problem into a compact matrix form, it is still a nonconvex problem due to the constraint $||W||_1=1$. To deal with the constraint $\mathbbm{1}^\top
W=0$, we add an extra term $\lambda(\mathbbm{1}^\top W)^2$ to the objective function, where $\lambda$ is a pre-defined hyper-parameter. Although this penalty
method cannot produce the exact global solution, we can still recover the sign of the global
solution (see Theorem \ref{thm:sign}). Once the sign of the global solution is identified, the remaining effort of computing the magnitude reduces to solving a convex problem (\ref{eq:convexdesign}).
\begin{theorem}
For large enough $\lambda$, the global solution $W^\ast$ of (\ref{eq:min_l1}) satisfies
$$
\textnormal{sgn}(W^\ast)=\textnormal{sgn}\left(\mathop{\arg\min}_{W\in\mathbb{R}^n, ||W||_1=1 } W^\top (YY^\top+\sigma I+\lambda \mathbbm{1}\mathbbm{1}^\top) W\right).
$$
\label{thm:sign}
\end{theorem}
The following theorem states that the problem is equivalent to another well-known NP-hard non-convex
problem --- Phase Synchronization \citep{singer2011angular}.
\begin{theorem}[Equivalence between Synthetic Design and Phase Synchronization] If $x^\ast\in\mathbb{R}^n$ is the global solution of $\min_{||x||_1=1} ||Ax||_2^2$ for some matrix $A\in\mathbb{R}^{D\times n}$ ($D>n$) and the matrix $A^\top A\in\mathbb{R}^{n\times n}$ is invertible, then $y^\ast=\textnormal{sgn}(x^\ast)$ is the global solution of $\max_{y\in\{-1,+1\}^n} y^\top((A^\top A)^{-1})^\top y$.
\label{theorem:equaltoPhaseSynch}
\end{theorem}
The proof of Theorem \ref{thm:sign} and Theorem \ref{theorem:equaltoPhaseSynch} is omitted in the main text due to page limit and is shown in Appendix \ref{section:equal}.
\begin{remark} Phase synchronization \citep{singer2011angular,bandeira2017tightness} aims to recover $n$ phases $z_i=e^{i\theta_i}, i\in[n]$ via solving the following optimization problem
\begin{equation}
\begin{aligned}
\max_{|x_1|=\cdots=|x_n|=1}\qquad x^\top C x,\label{eq:phaseSyn}
\end{aligned}
\end{equation}
where $C_{ij}$ is the noisy observation of $z_i\bar{z_j}=e^{i(\theta_i-\theta_j)}$. Our problem is symbolically equivalent to (\ref{eq:phaseSyn}). However, the data generating process is quite different from the Phase synchronization. The matrix we consider is the inverse of the data matrix. In Appendix, we will show that the design we find is actually the first $\ell_1$-principal component \citep{mccoy2011two,wang2021linear}.
\end{remark}
\section{Algorithm description}
\label{section:GIPW}
In this section, we propose the generalized power method with spectral initialization to solve our problem. The method is inspired by the lines of its early efforts for phase synchronization problems \citep{journee2010generalized,boumal2016nonconvex,liu2017estimation,zhong2018near}.
\begin{algorithm}[t]
\caption{\textbf{\underline{Synthetic principal component Design}}}\label{alg:sytheticPCD}
\begin{algorithmic}
\Require Pre-treatment Observations $Y\in\mathbb{R}^{N\times T}$
\State Set initial treatment assignment guess through $y^0=\textnormal{sgn}(v)$, where $v$ is the smallest eigenvector of matrix $(Y Y^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)$, where $\alpha,\lambda$ are two pre-defined hyper-parameter.
\State\Comment{\textbf{Spectral Initialization}}
\While{Converged}
\State Select one of the following two boxes to iterate
\begin{myorangebox}
\State For SPCD, update the design via \Comment{\textbf{Generalized power methods}}
\begin{equation}
\begin{aligned}
{\color{myorange}y^{t+1}=\textnormal{sgn}\left[\left((YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}+\beta I\right) y^{t}\right]},
\end{aligned}
\end{equation}
\State where $\beta$ is a pre-defined hyper-parameter.
\end{myorangebox}
\begin{mybluebox}
\State For NormSPCD, update the design via \Comment{\textbf{Normalize the inverse covariance matrix}}
\begin{equation}
\begin{aligned}
{\color{myblue}y^{t+1}=\textnormal{sgn}\Big[\left[(Y Y^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}+\beta I\right] (y^{t}/d)\Big]},
\end{aligned}
\label{update:normsgd}
\end{equation}
\State where {$d=\sqrt{\textnormal{diag}((Y Y^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1})}$} and $/$ denotes element-wise divide.
\end{mybluebox}
\EndWhile
\State Solve the following \emph{convex} optimization problem
\begin{equation}
\begin{aligned}
\{w_i\}_{i=1}^n=& \mathop{\arg\min}_{\{w_i\}_{i=1}^n} \,\, \frac{1}{T}\sum_{t=1}^T\left(\sum_{i:y(i)=1} w_iY_{it}-\sum_{i:y(i)=-1}w_iY_{it}\right)^2+\sigma \sum_{i=1}^Nw_i^2\\
& \quad \, \operatorname{s.t.} \quad \,\, w_i\ge 0, \,\, \forall i \in [N], \sum_{i:y(i)=1} w_i=\sum_{i:y(i)=-1} w_i=1.
\end{aligned}
\label{eq:convexdesign}
\end{equation}
\State Treat Unit $i$ if $y(i)=-\textnormal{sgn}\left(\sum_{i=1}^N y(i)\right)$ and run the experiment.
\State\Comment{\textbf{To ensure the size of the treated group is smaller than the control group}}
\State Estimate the treatment effect via
$$
\hat\tau=\sum_{t=1}^S\left(\sum_{i:y(i)=-\textnormal{sgn}\left(\sum_{i=1}^N y(i)\right)}w_iY_{i,T+t}-\sum_{i:y(i)=\textnormal{sgn}\left(\sum_{i=1}^N y(i)\right)}w_iY_{i,T+t}\right).
$$
\end{algorithmic}
\end{algorithm}
\vspace{-0.1in}
\subsection{Generalized Power Methods}
Spectral relaxation \citep{singer2011angular} is the first simple and efficient approach to solve the
phase synchronization problem. \citep{singer2011angular} relaxed the $N$ constraints
$|x_i|=1,i\in[N]$ to $||x||_2^2=n$. Then the solution becomes the leading
eigenvector. \citep{liu2017estimation,zhong2018near} showed that the eigenvector estimator is almost
close to the global optima under certain data generating process. Following these works, we take our
initial guess of the optimal experiment to be $\textnormal{sgn}(v)$, where $v$ is the smallest
eigenvector of the matrix $(Y Y^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)$ with
$\alpha,\lambda>0$ as two pre-defined hyper-parameters.
To further improve the experiment assignment, we utilize the generalized power method
\citep{journee2010generalized,luss2013conditional}, which considers the linearization of the
objective function at the current point and moves towards a minimizer of this linear function over
the non-convex set $\mathcal{C}$. It is also referred as to Frank-Wolfe algorithm for non-convex
problems in the literature.
Thus, the generalized power method enjoys monotonic improvement (\citep[Lemma 8]{boumal2016nonconvex}). When $g(x)$ is a quadratic function taking the form as $x^\top A x$, the method will become similar to the power methods for the eigenvector problem. The only difference is the normalization step. The power method normalizes the whole vector but the generalized power method normalizes the individual entries. The generalized power method can be also understood as projected gradient descent \citep{smith1994optimization}. Indeed, the update
\[
y^{t+1}=\textnormal{sgn}[((\frac{1}{\beta} YY^\top+\frac{\sigma}{\beta} I+\frac{\alpha}{\beta } \mathbbm{1}\mathbbm{1}^\top)^{-1}+I) y^{t}]
\]
can be understood as a projection step ($\textnormal{sgn}$) after a gradient descent update with step size
$\frac{1}{\beta}$:
\[
y^{t+1}=\textnormal{sgn}(y^{t}+\frac{1}{\beta}((YY^\top+\sigma I+\alpha \mathbbm{1}\mathbbm{1}^\top)^{-1})y^{t}).
\]
Thus the algorithm shares a sufficient ascent condition for each iteration (shown in Appendix
\ref{subsection:linearrate}). Our algorithm is called Synthetic principal component Design (SPCD) and
it is summarized in Algorithm \ref{alg:sytheticPCD}.
\paragraph{Normalized Variant}
In \citep{boumal2016nonconvex,zhong2018near,liu2017estimation}, the global optimality result is
highly dependent on the assumption that the top eigenvector of the iteration matrix lies in
$\{-1,1\}^N$. However, in our setup, the top eigenvector of the iteration matrix is the smallest
eigenvector of the covariance matrix which may not be a $\{-1,1\}^N$ vector. This is also the case
appearing in the phase retrieval \citep{candes2015phase,chen2015solving} and the degree corrected
stochastic block model \citep{zhao2012consistency,jin2015fast}. Inspired by the SCORE
\citep{jin2015fast,jin2022improvements} method for degree corrected stochastic block models, we
further introduce a normalization step to the generalization power method and call the new algorithm
Normalized SPCD (cf. NormSPCD, see \eqref{update:normsgd} for details). We use the diagonal
component of the inverse covariance as an estimate of the true normalization component. In Appendix
\ref{section:opt}, we show that NormSPCD can be interpreted as a Riemannian gradient descent with a
specific metric. Empirical results show that it is better than the original GPW in Figure
\ref{subfig:differentT}. This normalization technique may be of independent interest in other
applications.
\vspace{-0.1in}
\subsection{Global Guarantee}
In this subsection, we provide the global optimization guarantee for the (normalized) generalized power
method. \citep{bandeira2017tightness,boumal2016nonconvex,liu2017estimation,zhong2018near} have shown
that phase retrieval is globally solvable under certain generative models. We will show that
generalized power method can globally converge under certain data generating processes, which are
quite different from the ones assumed in the previous works. Following \citep{abadie2021synthetic},
we consider a realizable linear factor model (also referred to as "interactive fixed-effects model")
\citep{abadie2010synthetic,xu2017generalized,athey2021matrix,ferman2021properties}, which has already
been commonly employed in the literature as a benchmark model to analyze the properties of synthetic
control estimators \citep{amjad2018robust,li2020statistical}. Recently, \citep{shi2021assumptions}
justify the linear assumption from an independent causal mechanism viewpoint. The linear latent
factor model is stated in the following assumption. .
\begin{assumption}[Linear Latent Factor Model \citep{abadie2010synthetic,xu2017generalized}]
The outcomes are generated via the following linear factor model
\[
Y_{jt}=\delta_t+\frac{D_{jt}+1}{2}\tau+\theta_t^T \mu_j+e_{jt},
\qquad \mathbb{E}[e_{jt}|\delta_t,\mu_j,D_{jt}]=0,
\qquad \textnormal{Var}[e_{jt}|\delta_t,\mu_j,D_{jt}]=\sigma.
\]
Here $\delta_t$ is the time fixed effect; $\mu_j$ is the unobserved common factors; $\theta_t$ is a
vector of unknown factor loading; $e_{jt}$ is the unobserved i.i.d. idiosyncratic noise; $\tau$ is the
treatment effect that we aim to estimate and $D_{jt}$ is the $\{-1,1\}$ variable according to the treatment
assignment to unit $j$ at time $t$. More specifically, in the pre-treatment period, $D_{jt}=-1$ for
all $\forall j\in[N], t\in [T]$.
\label{assumption:factormodel}
\end{assumption}
To obtain the global optimality result, we further make the following realizable assumption that
there is only one realizable experiment (zero error experiment) in population.
\begin{assumption} [Realizable Assumption]
There exists a unique parameter $(w_i,D_i)_{i=1}^n$ that satisfies the following conditions:
\begin{itemize}
\setlength{\itemsep}{0pt}\setlength{\parsep}{0pt}\setlength{\parskip}{0pt}
\item $D_i$ are binary treatments, \emph{i.e.} $D_i\in\{-1,1\}$.
\item $w_i\ge 0$ and $\sum_{i=1}^n D_iw_i=0$.
\item $||w||_2^2=N$ and $\epsilon\le|w_i|\le\frac{1}{\epsilon}$ for all $\forall i \in [N]$.
\item The weights will balance the covariates, \emph{i.e.} $\sum_{i=1}^nw_iD_i\mu_i=0$.
\end{itemize}
\label{assumption:realizable}
\end{assumption}
\begin{remark}
This realizable assumption is similar to \citep{li2020statistical}, \citep[(5)]{shi2021assumptions}
and \citep[Assumption 3]{abadie2021synthetic}. The difference is that \citep[Assumption
3]{abadie2021synthetic} assumes the weight will cancel noisy observation of the untreated
outcome $Y_{jt}(-1)$ which is not realistic when the pre-treatment period $T$ is larger than the
number of units $N$. Our assumption is closer to \citep{li2020statistical} and
\citep[(5)]{shi2021assumptions}, but we further assume the uniqueness of the realizable experiment
that makes the optimization problem easier (in terms of no need to distinguish different
realizable experiments).
\end{remark}
Under Assumptions \ref{assumption:factormodel} and \ref{assumption:realizable}, we can show the following global optimality result and the proof is shown in the Appendix \ref{section:opt}.
\begin{theorem}
(Informal) Suppose that Assumptions \ref{assumption:factormodel} and \ref{assumption:realizable}
hold and that the latent time factor $[\theta_i^\top \delta_i]^\top$ is sampled from a underlying
distribution
with
mean $\Tilde{\theta}$ and covariance $\Tilde{\Sigma}$. Under regularity assumptions (see Appendix
\ref{subsection:generative}), if $\sigma$ is small enough and
$T\ge\textnormal{poly}(N,\frac{1}{\epsilon})$, then
\begin{itemize}
\setlength{\itemsep}{0pt} \setlength{\parsep}{0pt} \setlength{\parskip}{0pt}
\item If $\epsilon>\frac{\sqrt{3}}{2}-1$, then SPCD converges to the global optima.
\item If $\epsilon>0$, then NormSPCD converges to the global optima at a linear rate.
\end{itemize}
\end{theorem}
\section{Numerical Study}
\label{section:numerical}
This section report the numerical tests of our algorithm. In subsection \ref{subsection:simulate},
we validate our algorithm on the latent factor model. In subsection \ref{subsection:realworld}, we
demonstrate our algorithm on two real world datasets. Both experiments have shown the effectiveness
of our proposed experiment design algorithm in terms of the the root-mean-square error (RMSE), where
the squared differences between the true values of the treatment effects and the respective
estimates are computed for each treatment period and averaged.
We first introduce a simplified implementation of (Norm)SPCD, which although not guaranteed optimum but efficient, simple and effective in practice. In the simplified implementation, we don't solve the convex program (\ref{eq:convexdesign}) exactly, but using $w=\frac{2(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}y^\ast}{||(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}y^\ast||_1}$ to approximate instead. From \ref{eq:equalproof}, we know that once the optimal design profile $y^\ast$ is obtained, then $w$ is the optimal design weight. Notice that we don't exactly globally solve the problem (\ref{eq:min_l1}) in the simplified implementation, although we obtained the right experiment profile $y^\ast$ (Theorem \ref{thm:sign}). The weight we obtained here is the solve of the penalized approximation, but empirically it works good. The whole process is described in Algorithm \ref{alg:realalg}. In all the experiment in this paper, we use this simplified implementation.
\begin{figure}[]
\centering \includegraphics[width=3in]{timear.jpg}
\caption{Synthetic principal component design (SPCD) selects the treated units whose features are representative of the whole aggregate market of interest.}
\label{fig:syndesign_timefixedeffect}
\end{figure}
\vspace{-0.1in}
\subsection{Simulated Data}
\label{subsection:simulate}
We generate data from the linear factor model (also referred to as an "interactive fixed-effects
model", \citep{xu2017generalized,athey2021matrix}). The outcome $Y$ comes from
$$Y_{it}= v_t^\top\gamma_i+\tau \frac{D_{it}+1}{2}+e_{it},\,\forall i\in[N],\forall t\in[T+S].$$ where
$\gamma_i$ is a vector of latent unit factor of dimension $L$ generated as a standard Gaussian. We simulate the treatment effect $\tau$ with ground turth $1$. The idiosyncratic
noise is sampled from Normal $(0,\sigma^2)$ with $\sigma=1$. Fixing the test time period $S=10$ and the number of units $N=10$, we simulate different pairs of $(L,T)$ selection. We simulated three different choices of the time factor vector $v_t$
\begin{figure}
\centering
\subfloat[{$L=20,T=9$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallTrandom.jpg}
\label{subfig:arandom}}
\quad
\subfloat[{$L=20,T=20$}]{
\includegraphics[width=0.18\textwidth]{syndatarandom.jpg}\label{subfig:brandom}}
\subfloat[{$L=8,T=9$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallLrandom.jpg}
\label{subfig:crandom}}
\quad
\subfloat[{$L=8,T=20$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallLlargeTrandom.jpg}
\label{subfig:drandom}}
\caption{Treatment estimated via Synthetic Control and Synthetic Principle Design for data generated from pure random latent vector. We run the experiment over 100 runs of different seeds for different selections of $L,T$ on data generated from purely random latent vector.In all cases, Synthetic Principle Design provides more robust estimate of the true treatment effect 1.}
\label{figure:simulatedrandom}
\end{figure}
\begin{figure}
\centering
\subfloat[{\scriptsize$L=20,T=9$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallT.jpg}
\label{subfig:a}}
\quad
\subfloat[{\scriptsize$L=20,T=20$}]{
\includegraphics[width=0.18\textwidth]{syndata.jpg}\label{subfig:b}}
\quad
\subfloat[{\scriptsize$L=8,T=9$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallL.jpg}
\label{subfig:c}}
\quad
\subfloat[{\scriptsize$L=8,T=20$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallLlargeT.jpg}
\label{subfig:d}}
\caption{Treatment estimated via synthetic control (SC) and synthetic principal component design (SPCD) over 100 runs of different seeds for different selections of $L,T$ on data generated from time varying time factor.In all cases, Synthetic Principle Design provides more robust estimate of the true treatment effect 1.}
\label{figure:simulated}
\end{figure}
\begin{figure}
\centering
\subfloat[{$L=20,T=9$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallTar.jpg}
\label{subfig:aar}}
\quad
\subfloat[{$L=20,T=20$}]{
\includegraphics[width=0.18\textwidth]{syndataar.jpg}\label{subfig:bar}}
\subfloat[{$L=8,T=9$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallLar.jpg}
\label{subfig:car}}
\quad
\subfloat[{$L=8,T=20$}]{
\includegraphics[width=0.18\textwidth]{syndatasmallLlargeTar.jpg}
\label{subfig:dar}}
\caption{Treatment estimated via Synthetic Control and Synthetic Principle Design for data generated from an AR(1) process. We run the experiment over 100 runs of different seeds for different selections of $L,T$ on data generated from purely random latent vector.In all cases, Synthetic Principle Design provides more robust estimate of the true treatment effect 1.}
\label{figure:simulatedar}
\end{figure}
\begin{itemize}
\item \textbf{Pure Random Latent Vector} In this experiment, we follow \citep{abadie2021synthetic} and run our algorithm on a synthetic dataset where all the time latent factors is sampled from random Gaussian. We sample the latent unit factor $v\in\mathbb{R}^{N\times T}$ and latent time factor $\gamma\in\mathbb{R}^{N\times T}$ both as random standard Gaussian matrices. We fix the test time period $S=10$ and number of units $N=10$ and simulated different pairs of $L,T$ selection. The final results is shown in Figure \ref{figure:simulatedrandom}.
\item \textbf{Time Varying Factor} In this experiment, we generate the time factor vector $v_t$ as $t-\frac{T+S}{2}+\epsilon_t$, where
$t-\frac{T+S}{2}$ is time trend term and $e_{it}$ are i.i.d. standard Gaussian noise. The final results are
shown in Figure \ref{figure:simulated}.
\item \textbf{AR(1) Process} In this experiment, we follow \citep{abadie2022synthetic} and run our algorithm on a synthetic dataset where the time latent factors is sampled from an AR(1) process. In particular the time factor $\gamma=[\gamma_1,\gamma_2,\cdots,\gamma_T]'\in\mathbb{R}^{N\times T}$ is sampled via
\begin{itemize}
\item $\gamma_1\sim\mathcal{N}(0,I_N),$
\item $\gamma_{t+1}=A\gamma_1+b+\sigma\epsilon,\epsilon\sim\mathcal{N}(0,I_N)$.
\end{itemize}
In our experiment, we take $A=0.7I_N$, $b=\mathbbm{1}$ and $\sigma=1$. The final result is shown in Figure \ref{figure:simulatedar}.
\end{itemize}
Synthetic principal component design (SPCD) selects the treated units
whose features are representative of the whole aggregate market of interest \citep{abadie2021synthetic}. At the same time, SPCD surpasses SC in every setting in Figure \ref{figure:simulatedrandom}, \ref{figure:simulated} and \ref{figure:simulatedar}.
\vspace{-0.1in}
\subsection{Real World Data}
\label{subsection:realworld}
To exam our algorithm on real data, we follow
\citep{arkhangelsky2019synthetic,doudchenko2021synthetic} and utilize the US Bureau of Labor
Statistics and the Anti-Smoking Legislation data to examine the validity of our algorithm. Besides
synthetic control (SC), which randomly selects one unit to implement the treatment, we also
implement an additional random baseline, which randomly select units as control/group with
probability $1/2$. Through estimating the average treatment effect on the treated, we
compare SC and the random assignment baseline with SPCD in terms of the RMSE. The final result is
shown in Table \ref{table:result}. SPCD surpasses SC a large margin on both of the datasets.
\begin{figure}[h]
\centering
\subfloat[{\scriptsize Treatment Selected when $T=25$}]{
\includegraphics[width=0.34\textwidth]{tobacco.png}\label{subfig:usmap}}\quad
\subfloat[{\scriptsize Comparison in terms of RMSE}]{
\includegraphics[width=0.3\textwidth]{norm.jpg}\label{subfig:differentT}}
\subfloat[{\scriptsize Comparison with rerandomization }]{
\includegraphics[width=0.3\textwidth]{rerandom.jpg}\label{subfig:rerand}}
\caption{A typical design selected via synthetic principal component design (SPCD) and its performance.}
\label{figure:treatmentUS}
\end{figure}
\textbf{The Abadie–Diamond–Hainmueller Smoking Data.} \citep{abadie2010synthetic} uses SC to study
the effects of Proposition 99, a large-scale anti-smoking legislation program that California implemented in 1988. \footnote{The Abadie–Diamond–Hainmueller Smoking data is first organized by \citep{abadie2010synthetic}. In this paper, we use the organized version by \citep{arkhangelsky2019synthetic} at \url{https://github.com/synth-inference/synthdid/blob/master/data/california_prop99.csv} which drop the data of minimum wage laws, gun laws to abortion laws in the original data and only considers the smoking outcome data.} To simulate the bias of SC and SPCD on this application, following
\citep{athey2021matrix}, we consider observations for 38 states (excluding California due to Proposition 99) from 1970 through 2000. We regard the first $T$ year as pre-treatment periods to produce the design and use the last $31-T$ years as post-treatment periods to test the performance of the treatment assignment. The final result is shown in Table \ref{table:result} and Figure \ref{subfig:differentT}. Our design surpasses the random design by a large margin on most of the selection of time $T$. We also compare our method with the rerandomization design
\citep{morgan2012rerandomization,kallus2018optimal,li2018asymptotic} in Figure \ref{subfig:rerand}, which shows that our algorithm is still better than 500000 times of rerandomization.
One typical design produced by our algorithm is shown in Figure \ref{figure:treatmentUS}. The
experiment design for different pre-treatment length $T$ is shown in Figure
\ref{fig:differenttusamap}. The plots show that our selection of the control group is robust to
different pre-treatment time period and has the ability to represent all different geographic,
demographic, racial, and social structure of states in the United State.
\textbf{US Bureau of Labor Statistics.} We also apply our algorithm on the unemployment rate of 50 states in 40 months from the US Bureau of Labor Statistics (BLS). \footnote{{The BLS Statistics data is available
from the BLS website. In this paper, we use the organized version by
\citep{arkhangelsky2019synthetic} at
\url{https://github.com/synth-inference/synthdid/blob/master/experiments/bdm/data/urate_cps.csv}. We thank \citep{arkhangelsky2019synthetic}'s authors carefully organize the data and open source it on github.} } We run 50 simulations such that each simulation utilizes a 20-by-$T+S$ matrix sampled from the original 50-by-40 dataset. More specifically, we randomly select 20 units and use the first $T$ time
period to select the synthetic design and synthetic weight. The remaining $S$ time periods are the consecutive months that follow. In our experiment, we fix $S=5$ and run both experiment for $T=5,10$. The final result in terms of the RMSE is shown in Table \ref{table:result}.
\begin{table}[]
\centering
\setlength{\tabcolsep}{6pt}
\caption{Root-mean-square errors of the average treatment effect estimates by both synthetic control (SC) and synthetic principal component design (SPCD) on real data. The random design is simulated 10 times and 95$\%$ confidence interval is demonstrated. The reported RMSE for BLS dataset are multiplied by $10^3$ for readability.}
\begin{tabular}{@{}cccccc@{}}
\hline
\hline
\toprule
\multicolumn{5}{c}{\textbf{US Bureau of Labor Statistics }}
\\
\midrule
\multicolumn{3}{c}{$T=5$} & \multicolumn{3}{c}{$T=10$} \\ \cmidrule(l){1-3} \cmidrule(l){4-6}
\multicolumn{1}{c}{SC} & \multicolumn{1}{c}{Random} & SPCD & \multicolumn{1}{c}{SC} & \multicolumn{1}{c}{Random} & SPCD \\ \cmidrule(l){1-3} \cmidrule(l){4-6}
\multicolumn{1}{c}{14.5} & \multicolumn{1}{c}{7.5} & \textbf{0.9} & \multicolumn{1}{c}{11.6} & \multicolumn{1}{c}{5.6} & \textbf{0.6} \\
\midrule \midrule
\multicolumn{6}{c}{\textbf{Anti-smoking legislation}}\\
\midrule
\multicolumn{3}{c}{\textbf{ $T=15$}} & \multicolumn{3}{c}{\textbf{ $T=25$}} \\ \cmidrule(l){1-3} \cmidrule(l){4-6}
\multicolumn{1}{c}{SC} & \multicolumn{1}{c}{Random} & SPCD & \multicolumn{1}{c}{SC} & \multicolumn{1}{c}{Random} & SPCD \\ \cmidrule(l){1-3} \cmidrule(l){4-6}
\multicolumn{1}{c}{11.65} & \multicolumn{1}{c}{4.32$\pm$0.21} & \textbf{1.14} & \multicolumn{1}{c}{7.89} & \multicolumn{1}{c}{3.13$\pm$0.19} & \textbf{0.98} \\ \bottomrule
\end{tabular}
\label{table:result}
\vspace{-0.15in}
\end{table}
\begin{figure}
\centering
\subfloat[{$T=15$}]{
\includegraphics[width=0.25\textwidth]{T15.jpg}
\label{subfig:a2}}
\quad
\subfloat[{$T=20$}]{
\includegraphics[width=0.25\textwidth]{T20.jpg}
\label{subfig:b2}}
\quad
\subfloat[{$T=30$}]{
\includegraphics[width=0.25\textwidth]{T30.jpg}
\label{subfig:c2}}
\caption{Selection of control and treatment group in the Abadie–Diamond–Hainmueller California Smoking Data when different pre-treatment period length $T$ is available. The experiment design when $T=25$ is shown in Figure \ref{subfig:usmap}.}
\label{fig:differenttusamap}
\vspace{-0.1in}
\end{figure}
\vspace{-0.1in}
\section{Conclusion}
\label{section:conl}
In this paper, we consider the optimal experiment design problem for synthetic control, where an
NP-hard weighted covariate balancing problem is needed to be solve. Surprisingly, we reformulate the
problem as a phase synchronization problem and propose a fast spectral initialized (normalized)
generalized power method to address the resulting optimization problem efficiently. In face of a
realizable linear factor model, we provide the first global optimization results for experiment
design. Empirically, our method surpasses the original Synthetic Control and random design a large
margin in terms of RMSE on various datasets. Additionally, the newly proposed normalization
incorporated in GPM may have separate applications in degree correlated stochastic block models.
\begin{algorithm}
\caption{Empirical Implementation of SPCD}\label{alg:realalg}
\begin{algorithmic}
\Require Pre-treatment Observations $Y\in\mathbb{R}^{T\times N}$
\State Set initial treatment assignment guess through $y^0=\textnormal{sgn}(v)$, where $v$ is the smallest eigenvector of matrix $(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)$, where $\alpha,\lambda$ are two pre-defined hyper-parameter.
\State\Comment{Spectral Initialization}
\While{Converged}
\State Select one of the following two boxes to iterate
\begin{myorangebox}
\State For SPCD, update the design via \Comment{\textbf{Generalized power methods}}
\begin{equation}
\begin{aligned}
{\color{myorange}y^{t+1}=\textnormal{sgn}\left[\left((YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}+\beta I\right) y^{t}\right]},
\end{aligned}
\end{equation}
\State where $\beta$ is a pre-defined hyper-parameter.
\end{myorangebox}
\begin{mybluebox}
\State For NormSPCD, update the design via \Comment{\textbf{Normalize the inverse covariance matrix}}
\begin{equation}
\begin{aligned}
{\color{myblue}y^{t+1}=\textnormal{sgn}\Big[\left[(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}+\beta I\right] (y^{t}/d)\Big]},
\end{aligned}
\end{equation}
\State where {$d=\sqrt{\textnormal{diag}((YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1})}$} and $/$ denotes element-wise divide.
\end{mybluebox}
\EndWhile
\State Once obtained the optimal design $y^\ast$, one can select the design weight $w$ via
\begin{equation}
\begin{aligned}
w=\frac{2(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}y^\ast}{||(YY^\top+\alpha I+\lambda \mathbbm{1}\mathbbm{1}^\top)^{-1}y^\ast||_1}
\end{aligned}
\label{eq:weight}
\end{equation}
\State\Comment{The optimality condition ensures $\textnormal{sgn}(w)=y$.}
\State Treat Unit $i$ if $y(i)=-\textnormal{sgn}\left(\sum_{i=1}^N y(i)\right)$ and run the experiment.
\State Estimate the treatment effect via
$$
\hat\tau=\sum_{t=1}^S\left(\sum_{i=1}^N w(i)Y_{i,T+t}\right)
$$
\end{algorithmic}
\end{algorithm}
\begin{acknowledgement}
Yiping Lu is supported by the Stanford Interdisciplinary Graduate
Fellowship (SIGF). Jose Blanchet is supported in part by the Air Force Office of Scientific Research under award number FA9550-20-1-0397 and NSF grants 1915967, 1820942, 1838576. Lexing Ying is supported by National Science Foundation under award DMS-2011699. Yiping Lu also thanks Yitan Wang for helpful comments and feedback.
\end{acknowledgement}
\printbibliography