Exact Trend Control in Estimating Treatment Effects Using Panel Data with Heterogenous Trends

Exact Trend Control in Estimating Treatment Effects Using Panel Data with Heterogenous Trends

\iffull

\fi

Introduction

In this paper I propose new methods of estimating treatment effects for panel models with heterogenous trends. Two motivational numerical examples are illustrated in Figure (ref) based on simulated data generated by a model considered by Abadie, Diamond and Hainmueller (2010, ADH), with details given in Appendix (ref). Trends are plotted in the figure for the true untreated outcomes, the ADH synthetic control outcomes, and the construction by one of my new methods. In part (a) of Figure (ref), the ADH synthetic control outcomes are far from the truth even for the pre-treatment periods, presumably due to the violation of the convexity or interpolation assumption (ADH, 2010; Gobillon and Magnac, 2016; see also Figure (ref) in the appendix for the generated untreated outcomes). Though not much useful, the ADH results at least do not mislead the researcher as its inappropriateness is unequivocal. In part (b), however, the ADH synthetic control looks flawless for the pre-treatment periods, but the post-treatment synthetic control outcomes are far from the truth. Later developments such as Hsiao, Ching and Wan (2012, HCW hereafter) and Doudchenko and Imbens (2017) suffer from similar biases, while the methods I propose in this paper work well as Figure (ref) shows.

The model considered here is identical to ADH's (2010), and is given by

where $z_i$ and and $y_{it}$ are observed, with $y_{it} = y_{it}^1$ if the $i$th unit is treated in period $t$ and $y_{it} = y_{it}^0$ otherwise. Unit 1 is treated for $t>T_0$, and the rest units ($i=2,\ldots,J+1$) are untreated for all $t$. The unobservable trends $\gamma_t$ and $\delta_t$ are fixed effects that can be dependent on any other random variables. For example, $\gamma_t$ can be small in magnitude in the pre-treatment periods and large in the post-treatment periods. Similarly, $\mu_i$ are arbitrary fixed effects. The unobservable constituents $\mu_i$, $\gamma_t$ and $\delta_t$ should be normalized somehow for identification, but how they are normalized is of no consequence because I take the difference-in-differences (DID) approach. The observed vector $z_i$ contains $K+1$ components including the constant term for common time effects, and the unobservable $h_i$ has $r$ elements, where $K$ and $r$ are typically small. The variables $z_i$ and $h_i$ determine how each unit responds to common shocks $\gamma_t$ and $\delta_t$. The random errors $u_{it}$ are assumed to have zero mean conditional on $z_i$ and $h_i$.

The goal is to find $w_2, \ldots, w_{J+1}$ such that the linear combination $\sum_{j=2}^{J+1} w_j y_{jt}$ forms a sensible counterfactual comparison for the treated unit while controlling for the trends due to the common shocks $\gamma_t$ and $\delta_t$. Unlike Doudchenko and Imbens (2017) I firmly base my analysis on the model given by ((ref)). That is, the goal is to provide weights $w_2, \ldots, w_{J+1}$ such that $y_{1t}^0 - \sum_{j=2}^{J+1} w_j y_{jt}$ is free from confounding trends driven by $\gamma_t'z_i$ and $\delta_t'h_i$ in the model. Identification is sought not by algorithm but by the model and the population distribution of the related random variables.

My approach begins with distinguishing the variables responsible for trend heterogeneity and those which are balanced on in order to enhance comparability. When a set of variables (such as $z_i$ and $h_i$) are responsible for heterogenous trends, they should be exactly balanced on by hard constraints in order to avoid bias due to uncontrolled trends since the components $\gamma_t$, $\delta_t$, $z_i$ and $h_i$ are fixed effects; the balancing covariates such as pre-treatment outcomes, on the other hand, need not be exactly matched on.

The importance of exact matching on $z_i$ and $h_i$ has been overlooked in the literature. As discussed in Section (ref) later, ADH's (2010) algorithm is relevant in a subtle way but their nonnegativity constraint is obstruent. HCW's (2012) regression-based method and Doudchenko and Imbens's (2017) elastic-net proposal do not attend the heterogenous trends $\gamma_t'z_i$ and $\delta_t'h_i$. Consequences of ignoring its significance are visible in Figure (ref) above and in Figure (ref) later in Section (ref).

Exact-balancing on trending covariates does not obliterate the necessity of regularization, especially when the number $J$ of the untreated units is comparable to or larger than the number of balancing covariates as is the case in many applications. Without regularization the weight matrix may not be uniquely identified. Undoubtedly, all the extant methods implement regularization in some ways. ADH (2010) impose the nonnegativity and adding-up constraints as hard restrictions. HCW (2012) select a subset of control groups based on researcher's judgment. Doudchenko and Imbens (2017) implement elastic-net penalties. I also consider regularization, where the penalty term is motivated by ordinary least squares (OLS), rather than given heuristically. My proposal leads to a `constrained ridge' regression and its lasso and elastic-net variants, all of which are now well accepted by the econometric community.

The rest of this paper is organized as follows. Section (ref) presents the new estimators, and Section (ref) compares them with extant estimators. The last section contains concluding remarks. All the proofs are gathered in the appendix, which also contains discussions on establishing asymptotics. Throughout the paper, $Y_t$ and $U_t$ denote the $J\times 1$ vectors of $y_{jt}$ and $u_{jt}$, respectively, for $j\ge 2$, i.e., for the untreated units. The weight vector $(w_2, \ldots, w_{J+1})'$ is denoted by $w$, and $Z$ is the $(K+1)\times J$ matrix $(z_2,\ldots, z_{J+1})$. The exact-balancing restriction for $z_i$ is thus written as $z_1 = Zw$.

Estimation

This section presents the new estimators. Section (ref) considers a model with a common component $\gamma_t'z_i$ but without latent factors in order to motivate the exact-matching constraint $z_1=Zw$ and regularization. Section (ref) considers the same model but introduces balancing covariates. Section (ref) makes an extension to models with unobservable common factors.

Heterogenous trends on observables

To begin with, consider the model in ((ref)) without $h_i$ so the potential untreated outcomes are modeled by $y_{it}^0 = \mu_i + \gamma_t'z_i + u_{it}$, where $u_{it}$ shows no systematic trends if the model is correctly specified. For $s$ and $t$ with $s\le T_0 < t$, where $T_0$ is the last period before treatment, we have

with $I(\cdot)$ denoting the indicator function. In ((ref)), a post-treatment period $t$ is compared with a single pre-treatment period $s$ for the sake of simple exposition. Generalization by changing $y_{is}$ to $T_0^{-1} \sum_{s=1}^{T_0} y_{is}$ or any other weighted average makes no serious differences in the arguments to follow; likewise, $y_{it}$ can be replaced with an average over the post-treatment periods.

An obvious estimator of $\tau_{1t}$ in ((ref)) can be obtained by the OLS regression of $y_{it}-y_{is}$ on $I(i=1)$ and $z_i$ using the $J+1$ cross-sectional observations as the sample. Because the dummy variable $I(i=1)$ has value 1 only for $i=1$, the OLS estimator of $\gamma_t-\gamma_s$ is also obtained by regression $y_{it}-y_{is}$ on $z_i$ using $i\ge 2$, and then $\tau_{1t}$ is estimated as the prediction error for $i=1$. That is, the OLS estimator of $\gamma_t-\gamma_s$ is $(ZZ')^{-1}Z(Y_t-Y_s)$, and \[ \hat\tau_{1t} = (y_{1t}-y_{1s}) - (Y_t-Y_s)'Z'(ZZ')^{-1} z_1. \] With $w_a$ denoting $Z'(ZZ')^{-1}z_1$, this $\hat\tau_{1t}$ is written as $\hat\tau_{1t} = (y_{1t}-y_{1s}) - (Y_t - Y_s)'w_a$, which is the DID estimator using $Y_t'w_a$ as the constructed control group. In the simple case of $z_i=1$, the elements in $w_a$ are uniform, i.e., $w_a = J^{-1} (1,1,\ldots,1)'$, and $Y_t'w_a$ is the unweighted average $y_{jt}$ over the untreated units. In this sense $w_a$ generalizes the unweigted averaging operator. Note that $w_a$ depends on $z_i$ only and choice of $s$ and $t$ is irrelevant.

The weight vector $w_a$ eliminates the confounding trends driven by $z_i$ from $y_{1t} - Y_t'w_a$ because

and thus the DID estimator $\hat\tau_{1t}$ satisfies \[ \hat\tau_{1t} = \tau_{1t} + [ (u_{1t}-u_{1s}) - (U_t - U_s)'w_a]. \] We clearly have $\operatorname{E}(\hat\tau_{1t}) = \tau_{1t}$ because $w_a$ is a function of $z_1,\ldots, z_{J+1}$, provided that the random disturbances $u_{jt}$ have zero mean for all $t$ conditional on the trending covariates $z_1, \ldots, z_{J+1}$.

The above $w_a$ is not the only weight vector that gives an unbiased estimator of $\tau_{1t}$ by DID. Any $w$ satisfying $z_1=Zw$ and $\operatorname{E}(U_t'w)=0$ works because then $y_{1t}^0 - Y_t'w = (\mu_1-\boldsymbol{\mu}'w) + (u_{1t} - U_t'w)$. Given the arbitrariness of $\gamma_t$, unbiased estimation of $\tau_{1t}$ requires $z_1=Zw$ as a minimal condition, which is the exact balancing constraint emphasized in the introduction, and which $w_a$ turns out to satisfy.

It is noteworthy that $w_a=Z'(ZZ')^{-1} z_1$ is the solution to the constrained $\ell_2$ minimization

(See the appendix for a proof that $w_a$ solves ((ref)).) That is, $w_a$ is the smallest (in terms of Euclidean norm) of those satisfying $z_1=Zw$. Under the $iid$ assumption for $u_{jt}$, $w_a$ also minimizes the sampling variability in the constructed counterfactual outcomes conditional on $z_1,\ldots, z_{J+1}$, since $\operatorname{var}(U_t'w) = \sigma_u^2 w'w$ for nonrandom $w$. In plain words, $Y_t'w_a$ would exhibit least fluctuations over time while satisfying $z_1=Zw_a$.

It is subtle to discuss how a weight $w$ is defined for the model $y_{it}^0 = \mu_i + \gamma_t'z_i + u_{it}$. For a given $w$, let $\hat\tau_{1t}(w) = (y_{1t} - Y_t'w) - (y_{1s}-Y_s'w)$ for $s\le T_0<t$, which is the DID estimator using $Y_t'w$ as the constructed comparison group. The restriction that $\hat\tau_{1t}(w)$ should be unbiased for $\tau_{1t}$ alone does not identify a $w$ in the population since $z_1=Zw$ and $\operatorname{E}(u_{jt}w)=0$ are satisfied by infinitely many $w$'s, if $J>K+1$. For example, when $z_i=1$, any $J\times 1$ vector of fixed numbers that sum up to 1, such as the uniform weights, uneven weights like $w=(0.2,0.8,0,\ldots,0)'$, non-convex weights like $w=(-0.3,1.3,0,\ldots,0)'$, and infinitely many others, allows $\hat\tau_{1t}(w)$ to be unbiased for $\tau_{1t}$ if the model is correctly specified so that $\operatorname{E}(u_{it}|z_1,\ldots,z_{J+1})=0$ for all $t$. The weight $w_a = Z'(ZZ')^{-1} z_1$ is just one particular choice that generalizes the uniform weights. The identification of $w_a$ requires further the minimization of $w'w$ in ((ref)) on top of the unbiasedness requirement ($z_1=Zw$).

A natural alternative to $w'w$ in ((ref)) is the $\ell_1$ norm $\Vert w\Vert_1 = \sum_{j=2}^{J+1} |w_j|$, which leads to

a constrained $\ell_1$ minimization problem, also known as the basis pursuit minimization (see Mallat, 2009, Chapter 12). Algorithms using Alternating Direction Method of Multipliers (ADMM) are available for this problem (the R package ADMM). The minimization problem ((ref)) can also be written as the standard quadratic programming

because $w=w^+-w^-$ and $\Vert w\Vert_1 = w^+ + w^-$ for $w_j^+ = \max(w_j,0)$ and $w_j^- = -\min(w_j,0)$. Note that the $\ell_1$ minimization problem does not necessarily have a unique solution (e.g., when $z_i=1$), in which case we can minimize $\varepsilon w'w + \Vert w\Vert_1$ instead of $\Vert w\Vert_1$ for some small positive constant $\varepsilon$ such as $10^{-4}$ to achieve uniqueness (see Gains et al., 2018, p. 863). The elastic-net style loss function $\frac{1-\alpha}{2} w'w + \alpha \Vert w\Vert_1$ using other $\alpha$ parameter values can also be used. The elastic-net minimization algorithm can be implemented as a constrained lasso using $\alpha \Vert w\Vert_1$ as penalty, the zero vector as the response vector, and $[(1-\alpha)/2]^{1/2} I_J$ as the feature matrix. See James et al. (2019) for a fast algorithm for constrained lasso and its implementation by the R package PACLasso.

Given a weight vector $w$, the presence of systematic trends in the prediction error $y_{1t}-Y_t'w$ can be tested for the pre-treatment periods by regressing it on $t$, unless $T_0$ is too small. There is no `generated regressors' problem if $w$ is a function of $z_1, \ldots, z_{J+1}$. In addition, the mutual compatibility of two estimated weight vectors, $w_{(1)}$ and $w_{(2)}$, say, can be tested by regressing $Y_t'w_{(1)} - Y_t'w_{(2)}$ on $t-T_0$, $\mathit{after}_t$ and $\mathit{after}_t (t-T_0)$ using all the observations, where $\mathit{after}_t$ is the dummy variable for $t>T_0$. Overall significance can be interpreted as an evidence of model misspecification, although overall insignificance does not necessarily imply correct model specification because $U_t'[ w_{(1)} - w_{(2)}]$ can show no systematic trends while some $u_{it}$'s still do. If $z_i$ contains pre-treatment outcomes (e.g., ADH, 2010), the estimated $w$ is not necessarily exogenous, and the generated regressors problem applies. In that case, the testing results should be taken only as a diagnostic summary statistic. In all cases decision by human intuition using visual examination rather than formal testing is a promising alternative.

With regard to how to present the estimated counterfactual outcomes, if $z_i$ contains no pre-treatment dependent variables, then $Y_t'w$ and $y_{1t}$ may have systematically different levels just like in the standard DID framework. The counterfactual outcomes are, thus, better presented by $c+Y_t'w$ such that the intercept $c$ deals with the pre-treatment level difference. For example, $c$ can be the average of $y_{1s} - Y_s'w_a$ over the pre-treatment periods. This modification does not change anything about the estimation of treatment effects but only helps presentation.

Balancing covariates

We have thus far considered controlling for heterogenous trends driven by $\gamma_t'z_i$ by imposing the exact-matching constraints that $z_1=Zw$. In most applications the number $K$ of the nonconstant variables in $z_i$ is much smaller than the number $J$ of untreated units, and the restrictions $z_1=Zw$ do not identify a unique $w$. As a supplementary means to identify a single vector, we have considered minimizating the $\ell_2$, the $\ell_1$, or an elastic-net norm of $w$.

Now, beside the trending covariates $z_i$, the researcher may also want some other variables to be balanced on in pursuit of robustness against outliers or local misspecification. Typical balancing covariates include pre-treatment outcomes or their deviations from the pre-treatment average, while other exogenous features such as post-treatment controls can also be taken into consideration. Unlike the trend predictors $z_i$, these balancing covariates need not be matched on exactly.

Let $q_i$ denote the $m\times 1$ vector of such balancing covariates, e.g., $q_i = (y_{i1}, \ldots, y_{iT_0})'$, where $m$ can be larger or smaller than $J$. Let $Q$ be the $m\times J$ matrix of $q_i$ for the untreated units, i.e., $Q=(q_2, \ldots, q_{J+1})$. Matching seeks to make $(q_1-Qw)'(q_1-Qw)$ as small as possible, which leads to a natual extension of ((ref)) to

for a user-specified tuning parameter $\lambda \ge 0$ (and $\lambda>0$ if $Q'Q$ is singular). This is a constrained ridge (CRIDGE) regression of $q_1$ on $Q$ with penalty $\lambda w'w$ and constraints $z_1=Zw$. The shrinkage parameter $\lambda$ inversely relates to the desired matching quality relative to the magnitude $w'w$. If $\lambda=0$ (allowed if $Q'Q$ is nonsingular), we pursue best matching without shrinkage. If $\lambda=\infty$, we give up on balancing and pursue maximal shrinkage, leading to $w_a$ in the previous section. A finite positive $\lambda$ is a compromise. In all cases, we explicitly impose the restrictions that $z_1=Zw$, and thus heterogenous trends due to different $z_i$ are perfectly controlled for.

Given $\lambda$, the solution to ((ref)) is

where $\tilde{w}_{ridge} = G_{\lambda}^{-1} Q'q_1$ is the unconstrained ridge estimator (see the appendix for a proof). Note that $G_{\lambda}$ is invertible if $\lambda>0$ whether or not $Q'Q$ is, and thus $\hat{w}$ is well defined if $Z$ is of full row-rank and $\lambda>0$. The resulting treatment effect estimators are obtained by DID using $Y_t'\hat{w}$ as the constructed control group.

There is a more revealing expression for $\hat{w}$ than ((ref)). To derive it, let us first partial out $z_i$ from $Q$ and from $q_1$. Precisely, let $B = QZ' (ZZ')^{-1}$, the matrix of the OLS estimators from the regression of the rows of $Q$ on $Z'$, and let $\tilde{Q} = Q - BZ$ and $\tilde{q}_1 = q_1 - Bz_1$, the prediction errors. Then $\hat{w}$ is decomposed as follows:

which is the sum of the maximum shrinkage estimator $w_a$ subject to $z_1=Zw$ and the unconstrained ridge estimator $\hat{w}_b$ for balancing on the covariates orthogonal to $Z$ (proved in the appendix). Note that ((ref)) does not hold if the variables are automatically normalized in the ridge regression procedure, but whether to normalize $q_i$ or not is not critical under $z_1 = Z\hat{w}$, according to experiments. See Doudchenko and Imbens (2017) for more on normalization without the constraints.

By substituting the $\ell_1$ norm for the squared $\ell_2$ norm $w'w$ in ((ref)), we have the constrained lasso (CLASSO) version

where $\lambda$ is again a user-specified parameter. A fast optimization algorithm is available (James et al., 2019; see also Gaines et al., 2018). CRIDGE and CLASSO both shrink the parameters but only CLASSO achieves variable selection. Though a simple decomposition like ((ref)) is not available for CLASSO, the modified constrained lasso \[ \min_w\, \tfrac12 (\tilde{q}_1-\tilde{Q}w)' (\tilde{q}_1-\tilde{Q}w) + \lambda\Vert w\Vert_1\text{~~subject to~~} z_1 = Zw \] after partialing out $z_i$ from $q_1$ and $Q$ is identical to the original problem ((ref)). Again, if the balancing covariates are to be scaled within the optimization algorithm, the original variables and the variables after partialing-out give different results, naturally.

When usual constrained-lasso algorithms fail, one can again modify $\ell_1$ to a nominal elastic-net norm as Gaines et al. (2018) remark. The elastic-net objective function is $\frac12 (q_1-Qw)'(q_1-Qw) + \lambda (\frac{1-\alpha}{2} w'w + \alpha \Vert w \Vert_1)$, which equals the lasso objective function \[ \tfrac12 (q_1^{aug} - Q^{aug} w)' (q_1^{aug}-Q^{aug} w) + \lambda \alpha \Vert w\Vert_1, \] where $q^{aug} = (q_1',0)'$ and $Q^{aug} = [Q', \sqrt{\lambda (1-\alpha)} I_J]'$; see Gaines et al. (2018). Doudchenko and Imbens (2017) propose a cross-validation method of selecting $\lambda$ (and $\alpha$). I propose comparison by visualization after trying several different $\lambda$ values.

For the model given by ((ref)), ADH (2010) treat the constant term in $z_i$ and the nonconstant terms differently, where the constant term is exactly matched on by the adding-up constraint and the nonconstant terms appear in minimization. My approach, on the other hand, treats all the terms in $z_i$ identically by exact matching. Balancing $q_1$ and $Qw$ is a different issue; they are matched by the minimization of $(q_1-Qw)'(q_1-Qw)$ without requiring exact balancing. The roles of trending covariates and balancing covariates are different, which is natural considering that $z_i$ appears in the model as the drivers of nuisance trends and $q_i$ is introduced to enhance comparability.

A practical remark on the selection of $\lambda$ is worth making. If the model is correctly specified so that $u_{it}$ shows no systematic trends, i.e., if $u_{it}$ has zero mean conditional on $z_1,\ldots, z_{J+1}$ for all $t$, then any $w$ satisfying $z_1=Zw$ will eliminate confounding systematic trends in $y_{1t} - Y_t'w$. When it happens, the choice of $\lambda$ would not make much difference in principle. On the other hand, systematicity in $y_{1t} - Y_t'w$ in the pre-treatment periods would be an evidence of possible misspecification of the model for some $i$ or all, in which case matching on variables such as pre-treatment outcomes will hopefully mitigate the problem. Since a larger $\lambda$ deteriorates the matching quality and increases the variability in $Y_t'w$, it would be an acceptable practice to enlarge $\lambda$ while keeping the discrepancy between $q_1$ and $Qw$ within a tolerable range. Though fuzzy theoretically, the acceptability is usually clear to human eyes as the time-series of $y_{1s}$ and $Y_s'w$ in the pre-treatment periods can be visually compared without difficulty. Also, the constrained shrinkage estimators are continuous in $\lambda$ (except at $\lambda=0$ for which $Q'Q$ may be singular) for given data, and small changes in $\lambda$ will lead to only small changes in the trend of $Y_t'w$.

Unobservable factors

We have thus far considered the case $h_i$ is empty in ((ref)). In many application, a few variables in $z_i$ would be sufficient as the driving force of trend heterogeneity. Besides, soft matching on the lagged dependent variables often obliterates the necessity of unobservable common factors. In some cases, however, researchers may want to allow for unobservable $h_i$, especially if no observable trending covariates are available. In this section, we discuss how to handle $h_i$.

Because $h_i$ makes heterogenous trends, it is again essential to have $h_1$ and $Hw$ exactly balanced, where $H=(h_2,\ldots,h_{J+1})$. But this is infeasible since $h_i$ are not observed. ADH (2010) replace $h_1 = Hw$ with the sufficient condition that $y_{1s} = Y_s'w$ for all $s\le T_0$, which is not attainable unless $T_0$ is smaller than $J$. But even when $J$ is large enough for $y_{1s}=Y_s'w$ for all $s\le T_0$, the nonnegativity of $w_j$ imposed by ADH (2010) does not necessarily guarantee $z_1=Zw$ and $y_{1s} = Y_s'w$ at the same time. Adverse examples have been illustrated in Figure (ref).

When $h_i$ are unobserved, an obvious strategy is to estimate them rather than attempting to find a detour. If $\check{h}_i$ denotes the initial estimator of $h_i$ and $\check{H} = (\check{h}_2, \ldots, \check{h}_{J+1})$, the corresponding $\ell_2$ optimization problem is \[ \min_w \; (q_1-Qw)'(q_1-Qw) + \lambda w'w \text{~~subject to~~} z_1=Zw \text{ and } \check{h}_1 = \check{H}w. \] There are total $1+K+r$ constraints, which are generally satisfied by nonempty parameters if $J> K+1+r$, which holds in usual applications. If $J$ is too small, the researcher would try to reduce $K$ or $r$ or both; it is not very sensible to have more common factors than the number of untreated units in applications.

A convenient way of estimating $h_i$ is to use least squares using the pre-treatment data: \[ \min_{\substack{\mu_1, \ldots, \mu_{J+1}\\ \gamma_1, \ldots, \gamma_{T_0}\\ \delta_1,\ldots,\delta_{T_0}\\ h_1, \ldots, h_{J+1}}} \; \sum_{i=1}^{J+1} \sum_{t=1}^{T_0} (y_{it} - \mu_i - \gamma_t'z_i - \delta_t'h_i)^2, \] or in matrix notations \[ \min_{\mu^*, \Gamma, F, H^*} \operatorname{tr} \Big\{ (Y^* - 1\mu^{*\prime }- \Gamma Z^* - \delta H^*)' (Y^* - 1\mu^{*\prime }- \Gamma Z^* - \delta H^*) \Big\}, \] where $Y^*$ is the $T_0\times (J+1)$ matrix of $y_{it}$ for $i=1,\ldots, J+1$ (columns) and $t=1,\ldots,T_0$ (rows), $\mu^*$ is the $(J+1)\times 1$ vector of $\mu_i$, $i=1,\ldots,J+1$, $\Gamma = (\gamma_1, \ldots, \gamma_{T_0})'$, $Z^* = (z_1,Z)$, $\delta=(\delta_1,\ldots, \delta_{T_0})'$, and $H^* = (h_1, H)$. The concentrated loss function is

where $M_1 = I_{T_0} - T_0^{-1}11'$ and $M_{Z^{*\prime }} = I - Z^{*\prime }(Z^* Z^{*\prime })^{-1} Z^*$. Let $A=M_1 Y^* M_{Z^{*\prime }}$. The common factors in $A$ are estimated as $\sqrt{T_0}$ times the orthonomal eigenvectors of $AA'$ corresponding to the $r$ largest eigenvalues, and the associated factor loading estimators are $(\tilde{h}_1, \ldots, \tilde{h}_{J+1}) = T_0^{-1} \tilde{\delta}'A$, where $\tilde\delta$ is the matrix of estimated common factors. Note that the estimated common factors correspond to $M_1\delta$ rather than $\delta$ itself, and the estimated factor loadings to $H_*^{\dag} = H^* M_{Z^{*\prime }} = H^* - H^* Z^{*\prime }(Z^*Z^{*\prime })^{-1} Z^* = [h_1, H] - H^* Z^{*\prime }(Z^* Z^{*\prime })^{-1} [ z_1, Z]$ rather than $H^*$ itself. But, given that $z_1=Zw$, we have $h_1 = H w$ if and only if $h_1^{\dag} = H^{\dag} w$, where $H_*^{\dag} = [h_1^{\dag}, H^{\dag}]$. We can therefore use the estimated factor loadings $\tilde{h}_i$ in the constrained ridge, lasso and elastic-net optimization. Although $h_1^{\dag}$ and $H^{\dag} \hat{w}$ are not exactly balanced due to the discrepancy of $\tilde{h}_i$ and $h_i^{\dag}$ (after rotation),

It is nuisance that the constrained estimator vector $\hat{w}$ satisfies $z_1 = Z\hat{w}$ and $\tilde{h}_1 = \tilde{H} \hat{w}$, but not $h_1^{\dag} = H^{\dag} \hat{w}$ or $h_1 = H\hat{w}$. Thus, $y_1 - Y_t\hat{w}$ still contains a remaining trend term as shown in \[ y_1 - Y_t\hat{w} = (\mu_1 - \boldsymbol{\mu}'\hat{w}) + \delta_t' (h_1 - H\hat{w}) + (u_{1t} - U_t' \hat{w}). \] But, given that $z_1=Z\hat{w}$ and $\tilde{h}_1 = \tilde{H} \hat{w}$, we have \[ \delta_t' (h_1 - H\hat{w}) = \delta_t'B^{-1} [ (Bh_1^{\dag} - \tilde{h}_1) - (BH^{\dag} - \tilde{H}) \hat{w} ] \]

If some common factors are observed (e.g., incidental linear or quadratic trends), then they can be partialed out by replacing the $M_1$ matrix in ((ref)) with an appropriate projection matrix. For example, if $y_{it}^0 = \gamma_t' z_i + g_t' \mu_i + \delta_t'h_i + u_{it}$, where $g_t$ is observable and the fixed effects are subsumed in $g_t'\mu_i$, then $M_1$ is to be replaced with $M_{[1,g]}$, say, where $g = (g_1,\ldots, g_{T_0})'$. Finally, the number $r$ of common factors may be chosen exogenously by the researcher or by using an automatic selection procedure. I recommend the former method. Specifically, increasing $r$ starting from zero and plotting the estimated counterfactual outcomes will give the researcher clear ideas how the results change as more hidden factors are allowed for in the model.

Comparison with extant estimators

This section compares the new methods with ADH (2010), HCW (2012), and Doudchenko and Imbens (2017).

Comparison to ADH (2010)

ADH's (2010) synthetic control algorithm consists of two layers of optimization, which I call the `inner' and `outer' optimization loops. The inner loop finds an optimal $\hat{w}(V)$ for a given $V$ by minimizing $(z_1-Zw)' V (z_1-Zw)$ subject to the adding-up and nonnegativity constraints (called the `ADH constraints' in short in this subsection), and the outer loop finds an optimal diagonal positive semidefinite $V$ by minimizing $\sum_{s=1}^{T_0} [y_{1s} - Y_s'\hat{w}(V) ]^2$. The final weight estimator is $\hat{w} = \hat{w}(\hat{V})$. ADH (2010) also discuss using a user-specified $V$.

For a given $V$, if there exists a $w$ satisfying the ADH constraints and the exact-balancing condition $z_1=Zw$ simultaneously, the inner-loop loss function $(z_1-Zw)' V(z_1-Zw)$ attains zero at such a $w$. Even in that case, however, a unique $w$ is not identified in general because the constraints are linear in $w$. For example, if $z_1=0$, a scalar, and $(z_2, z_3, z_4, z_5) = (-2,-1,1,2)$, any symmetric kernels such as $w=(\frac14, \frac14, \frac14, \frac14)'$, $w=(0,\frac12, \frac12,0)'$, etc., minimize the loss function for the inner optimization loop. In such a case a particular weight will be chosen arbitrarily by the numerical procedure used for the optimization. In contrast, if no $w$ satisfies both the ADH constraints and the exact-balancing condition simultaneously, then ADH's algorithm sacrifices exact balancing to abide by the ADH constraints. The consequences of abandoning exact balancing to save the ADH constraints are illustrated in Figure (ref) as discussed repeatedly.

The $V$-weight is determined by the outer-loop minimization for balancing the pre-treatment outcomes. (If a fixed $V$ is used, balancing the pre-treatment outcomes is irrelevant.) For the finally chosen $\hat{V}$, no matter whether it is the outcome of the outer-loop optimization or given exogenously, the solution $\hat{w}=\hat{w}(\hat{V})$ need not be unique nor satisfy $z_1 = Z\hat{w}$. Notably, the selection of $V$ is blind to whether $z_1=Z\hat{w}(V)$ because $V$ is chosen by the outer loop involving only the pre-treatment outcomes. For example, if some $V$ allows for $z_1=Z\hat{w}(V)$ and others do not, the ADH algorithm does not necessarily choose the one that allows for $z_1=Z\hat{w}(V)$ since $V$ is determined by minimizing $\sum_{s=1}^{T_0} [y_{1s} - Y_s'\hat{w}(V)]^2$, which does not necessarily minimize $[z_1-Z\hat{w}(V)]' [z_1-Z\hat{w}(V)]$.

The nonnegativity and adding-up constraints provide attractive interpretations to practitioners, but the benefits come with nontrivial costs. First, ADH's (2010) two-layer optimization procedure may fail to converge or give a suboptimal choice of synthetic control. For example, Abadie and Gardeazabal (2003) find the `Synthetic Basque' of $0.851\times \text{Cataluna} + 0.149 \times \text{Madrid}$' in their study on the political turmoil in Spain. But a thorough investigation reveals that a lower root mean squared prediction error can be achieved by an alternative synthetic Basque of $0.633\times \text{Cataluna} + 0.148\times \text{Madrid} + 0.219 \times \text{Baleares}$. (Finding this weight vector requires more direct use of the Karush-Kuhn-Tucker theorem. Neither the Stata `synth' package nor the R `Synth' package identifies this synthetic control.) This suggests that researchers should not be overly confident about the meaningfulness of the estimated $w$ weights.

The second issue involves the nonnegativity, and is more subtle. The nonnegativity constraint may violate $z_1=Zw$, i.e., $\mathbb{R}^J_+ \cap \{ w\colon z_1=Zw \} = \emptyset$, in which case trends in $y_{1t} - Y_t'w$ due to $z_1-Zw$ may confound the treatment effects if $w$ is forced to be in $\mathbb{R}_+^J$. The importance of nonnegativity can be controversial, but it is noteworthy that a discrepancy between $z_1$ and $Zw$ can lead to a nonnegligible confounding trend in $y_{1t}^0-Y_t'w$ while a negative $w_j$ only affects interpretation. If one wishes, the nonnegativity restriction can be made soft by, for example, the constrained lasso \[ \min_{w^+, w^-} \tfrac12 \Vert q_1-Qw^++Qw^-\Vert_2^2 + \lambda \sum_{j=2}^{J+1} (w^+_j + \kappa w^-_j), \] for some large positive $\kappa$, subject to the constraints that $z_1 = Zw^{+} - Zw^{-}$, $w^+_j\ge 0$ and $w^-_j\ge 0$ for all $j$, which modifies a generalized version of ((ref)). The above soft nonnegativity will allow $w_j<0$ for some $j$ if hard nonnegativity is incompatible with $z_1=Zw$, but will try to keep $w_j$ as close to the nonnegative domain as possible. However, the benefit looks only minor because the appealing interpretation attached to nonnegativity is lost anyway if some $w_j$ are negative.

Comparison to HCW (2012)

HCW (2012) take an alternative approach of regressing $y_{1t}$ on $Y_t$ for a selected subset of the untreated units using the pre-treatment observations to estimate the intercept $c$ and the slope vector $w$. Then the counterfactual outcomes are formed as $\hat{c} + Y_t'\hat{w}$ for $t>T_0$, where $\hat{c}$ and $\hat{w}$ are the OLS estimators. As Li and Bell (2017) derive, this estimator is justified under mean stationarity. If the unobserved trends show mean nonstationarity, HCW's (2012) method needs modification.

To see the source of bias and its remedy, let us take a simple example with $z_i=1$. Given the OLS estimators $\hat{c}$ and $\hat{w}$, the estimated treatment effects are

where $\ddot\gamma_t - \gamma_t - \bar{\gamma}_{\ifthenelse{\equal{}{}}{}{,}pre}$, $\ddot{\delta}_t = \delta_t - \bar{\delta}_{\ifthenelse{\equal{}{}}{}{,}pre}$, $\ddot{u}_{1t} = u_{1t} - \bar{u}_{\ifthenelse{\equal{1}{}}{}{1,}pre}$, and $\ddot{U}_t = U_t - \bar{U}_{\ifthenelse{\equal{}{}}{}{,}pre}$, with $\bar{\xi}_{\ifthenelse{\equal{}{}}{}{,}pre}$ denoting $T_0^{-1} \sum_{t=1}^{T_0} \xi_t$ for variable $\xi_t$.

The OLS regression of $y_{1t}$ on $Y_t$ for $t\le T_0$ may give systematic biases in $\hat\tau_{1t}$ for this model due to the $\ddot\gamma_t (1-1'\hat{w})$ term among others, because the stated OLS regression does not guarantee $1'\hat{w} \xrightarrow{p} 1$. The origin of this failure is in fact endogeneity. Example (ref) below demonstrates that $1'\hat{w}<1$ asymptotically (as $T_0\to\infty$) if $y_{1t}$ is regressed on $Y_t$ for $t\le T_0$ for a model with $z_i=1$ and empty $h_i$, so that systematic changes in trend ($\gamma_t$) may confound the treatment effects.

An obvious solution to the problem is to impose the restrictions that $1'w=1$ in case $z_i=1$ as in Example (ref) and that $z_1=Zw$ for general $z_i$, which is exactly our exact-balancing constraint. If $h_i$ is nonempty in ((ref)), then $h_i$ can be estimated and the constraints that $\tilde{h}_1 = \tilde{H} w$ can be added as explained in Section (ref). Because the number of common factors are typically small, there exist almost certainly some $w$ vectors that satisfy the restrictions. This modified HCW method is a special case of the constrained ridge regressions proposed in this paper corresponding to $\lambda=0$.

The above constrained OLS is easy to implement, but it requires $T_0>J-K-1$. If there are many untreated units ($J$ large), HCW (2012) select a sufficiently small subset a priori by the researcher's judgment, which is sometimes arbitrary but often acceptable as long as rationales are provided. The constraints that $z_1=Zw$ and that $h_1=Hw$ are always crucial.

Doudchenko and Imbens's (2017) elastic net

Doudchenko and Imbens (2017) propose minimizing the elastic-net loss function $\sum_{s=1}^{T_0} (y_{1t} - c - Y_t'w)^2 + \lambda (\frac{1-\alpha}{2} \Vert w\Vert_2^2 + \alpha \Vert w\Vert_1)$ without constraints. Their proposal (elastic net, and no constraints) can be understood as a modification of ADH (2010) and also a modification of HCW (2012) to an elastic-net framework. When signal is strong in the pre-treatment period such that matching on the observed pre-treatment outcomes deals with trends adequately, this elastic-net solution may work well (though bias may still exist due to the endogeneity reason explained in Section (ref)), but otherwise there is no device to control for heterogenous trends in the outcomes in the post-treatment periods.

Let us take numerical examples. Figure (ref) is obtained by applying Doudchenko and Imbens's (2017) proposal to the two simulated data sets considered for Figure (ref). The elastic-net mixing parameter is set to $\alpha = 0.9$ (close to lasso), and the tuning parameter is $\lambda=0.01$, a value that gives a visually appealing pre-treatment matching; larger $\lambda$ values such as 0.1 and 1 are poor in reproducing the trend in the pre-treatment outcomes. The results are compromised for both data sets in the post-treatment periods, which seems to be due to the endogeneity bias discussed in Section (ref). Imposing $1'w=1$ as a hard restriction controls for common time effects, and $z_1=Zw$ for more general models, gives the elastic-net version of what the present paper proposes.

It is noteworthy that Doudchenko and Imbens (2017) do not refer to an explicit model; see their introduction. In other words, their aim is not at controlling for heterogenous trends for models like ((ref)) but at estimating counterfactual trends based on regularized matching on pre-treatment outcomes (identification by regularization).

Conclusion

For model ((ref)) considered by ADH (2010), I propose new estimators of treatment effects by treating the trending variables ($z_i$ and $h_i$ in the model) and other balancing covariates (denoted $q_i$ in this paper) differently. Without further assumptions on the time-varying coefficients ($\gamma_t$ and $\delta_t$ in the model), exact-balancing of the trend predictors as hard restrictions is crucial for properly dealing with heterogenous trends driven by the trending covariates. The adverse consequences of making the exact matching soft are illustrated in Figures (ref) and (ref), where all the extant estimators exhibit compromised behaviors for data generated by ((ref)) without hidden factors. The new estimators proposed in this paper work well.

\iffull

References

\begingroup \setlength\labelsep{0pt}

\endgroup

\fi