Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
106,015 characters · 20 sections · 76 citation commands
On Policy Evaluation With Aggregate Time-Series Instruments
\thispagestyle{empty}
\baselineskip=20pt \setcounter{page}{1}
Aggregate events are often used to identify policy-relevant quantities. A classic example is card1990impact, where a single aggregate policy change -- the Mariel Boatlift of 1980 -- allows the author to quantify the effect of immigration on the Miami labor market by comparing it to unaffected locations before and after the shock. Since then, this difference in differences (DiD) strategy has become the standard tool for empirical analysis in similar settings and spurred a vast methodological literature; see arkhangelsky2023causal for a recent survey. However, many problems in empirical economics do not conform with the standard DiD setting because aggregate events occur frequently and affect all units simultaneously. Developing methodology for such applications is the primary goal of this paper.
We focus on environments where the aggregate variation is used as an instrument for unit-specific policy variables. In applications we have in mind the variables of interest are often equilibrium objects, and external instruments are necessary to solve the endogeneity problem. We develop a new estimator for such applications, establish its statistical properties, and show its superiority to the benchmark that dominates current empirical practice -- the two-stage least-squares (TSLS) regression with unit and time fixed effects.
Formally, suppose we observe outcome $Y_{it}$ and treatment $W_{it}$ for $n$ units over $T$ periods. We also have access to aggregate instrument $Z_t$, and our goal is to quantify the contemporaneous effect of $W_{it}$ on $Y_{it}$. For example, in nakamura2014fiscal, $W_{it}$ is the military procurement spending in state $i$ and year $t$, $Y_{it}$ is a measure of state-level growth, $Z_t$ is the national military spending in the US, and the contemporaneous effect describes a local fiscal multiplier. It is natural to view both $W_{it}$ and $Y_{it}$ as jointly determined in equilibrium, forcing the authors to address the endogeneity problem to identify the multiplier. Having access to an aggregate instrument suggests an easy solution: as long as $Z_t$ satisfies conventional assumptions of imbens1994identification, we can establish a causal link between $Y_{it}$ and $W_{it}$ by constructing an instrumental variables (IV) estimator separately for each unit $i$ and reporting a summary of these estimators, e.g., the average.
In practice, however, this naive approach is unlikely to produce a credible answer. Researchers suspect that $Z_t$ correlates with other unobserved aggregate variables that affect the outcomes, and thus, the assumptions of imbens1994identification are not satisfied. For example, in nakamura2014fiscal, the authors are concerned that aggregate military spending is correlated with fiscal and monetary policies, directly affecting local outcomes. Similar problems are evident in simpler environments: a comparison of Miami's labor market before and after the Mariel Boatlift is likely to be corrupted by aggregate confounders.
To address these concerns, researchers are adopting a different strategy. In particular, the following equation is estimated by the TSLS:
using $D_iZ_t$ as an instrument.\footnote{Applications abound from different domains. See, e.g., studies of the effect of Airbnb on the housing market barron2021effect, plantation owners' power on wages and incarceration rates dippel2020outside, commodities prices on violent conflict dube2013commodity, dams on poverty and development duflo2007dams, housing wealth effects guren2018housing, local fiscal multipliers nakamura2014fiscal, foreign aid on conflict nunn2014us and a few others.} Here, $D_i$ is a measure of “exposure” of unit $i$ to the aggregate variable, and $\tau$ is the parameter of interest. The critical part of equation (ref) are time-fixed effects that are meant to capture the unobserved aggregate confounders. To understand the mechanics of this estimator, suppose we have only two units $i =1,2$, and we know that the first unit is affected by $Z_t$, $D_1 =1$, while the second one is unaffected at all, $D_2 =0$. Then, by looking at differences across units, we eliminate the aggregate confounder as long as it affects both units in the same way. This procedure is a natural extension of the conventional DiD strategy, and it is precisely what the TSLS estimation of (ref) amounts to for the case with two units. Arguably, not much else can be done in this setting.
The situation changes once we have access to multiple units. In this case, the TSLS estimator first averages units with high and low values of $D_i$ and then takes the difference between the resulting averages. This particular aggregation scheme is valid and can even be statistically efficient as long as we believe that the unobserved confounder affects all units in the same way. While unavoidable in the case of two units, this assumption becomes increasingly dissatisfying once we have multiple units. Researchers recognize this threat and commonly view it as the main danger to the validity of the TSLS identification strategy (e.g., chodorowstock,guren2018housing). The estimator we propose directly addresses this threat.
To construct our estimator, we reserve a part of the sample to learn weights $\omega^{rob}_i$, which we then use to aggregate the rest of the data and to build a single IV estimator. To produce $\omega^{rob}_i$, we first project out the effect of the aggregate instrument and use residuals to construct a combination of units with higher values of $D_i$, which resembles a combination of units with lower values of $D_i$. The identification assumption that justifies this approach is the existence of a particular aggregation scheme (unit-specific weights), which eliminates the unobserved confounding when applied to the data. This condition is valid whenever $D_i$ cannot be fully explained by systematic unit-level variation in outcomes, i.e., a no-multicollinearity restriction. We show that this weak restriction holds in various models, all of which generalize the one that justifies the TSLS estimation. A similar identification assumption is behind the Synthetic DiD (SDiD) estimator introduced in arkhangelsky2019synthetic (see also the recent analysis in imbens2023identification). The development of the new method that extends the SC ideas to settings with aggregate instruments and thus vastly expands the applicability of the SC approach constitutes the first major contribution of this paper.
The second major contribution lies in the formal framework we use to analyze our method and the statistical properties we establish within this framework. In particular, we combine insights from two strands of the literature: the one on quasi-experimental designs and the one on the properties of the SC-type estimators.
We consider a finite population model where the sole source of uncertainty comes from aggregate variables: the observed instrument and unobserved confounders. This approach takes its roots in the experimental literature neyman1923applications,rubin1977assignment, but it has also been successfully applied to observational data with the recent methodological literature emphasizing the importance of this framework both for estimation borusyak2022quasi,borusyak2023nonrandom and inference abadie2020sampling,abadie2023should,adao2019shift. The key new ingredient we bring to this setup is the presence of unobserved confounders, which we explicitly model. We expect this idea to be useful more broadly because confounding is a natural threat in all observational studies. The finite population framework is especially appealing in our setting for two reasons. First, it allows us to focus explicitly on the key identifying variation, which comes from the aggregate variables. Second, in many relevant empirical applications, it is natural to treat the observed units as the whole population rather than a sample. For example, this is the case when units represent geographic entities, e.g., states or counties.
We use this framework to establish the statistical properties of the new estimator. We analyze it in a high-dimensional regime where $n$ is similar or larger than $T$. In many applications of interest, $n$ and $T$ are comparable, making this regime natural. We prove that our algorithm delivers consistent and $\sqrt{T}$-convergent estimators, even in the presence of aggregate confounding. We also show how to use our method to conduct valid inference as long as there is enough variation in the baseline outcomes. To derive these results, we build on recent insights in SC literature arkhangelsky2019synthetic,ferman2021synthetic, and connect empirical weights $\omega_i^{rob}$ to a specific deterministic analog. However, our analysis is necessarily different for two reasons. First, the policy variable in our environment is endogenous, and the construction of weights relies on aggregate instruments, which are not present in the conventional SC literature. Second, the only source of uncertainty in our model comes from the aggregate variation. To the best of our knowledge, this paper is the first to derive any statistical properties of SC-type estimators in such models. More recently, imbens2023identification follow a related path when analyzing a class of SC estimators.
The TSLS estimator is widely adopted for a reason -- the transparency and flexibility of the method compensate for the underlying strong assumptions. Moreover, practitioners might expect (ref) to provide a good approximation in many applications, rendering confounding concerns less relevant. To investigate this, we demonstrate the benefits of our approach using the data from nakamura2014fiscal. First, we reevaluate their study using our method and find fiscal multipliers larger in magnitude than the original ones. One of the goals of this exercise is to demonstrate to the applied-minded reader that our method is highly transparent and, by using it, one can gain insights about the TSLS estimator itself. We then construct a simulation that mimics the properties of the original dataset of nakamura2014fiscal. We use this simulation to show that our estimator remains competitive in simple designs, can outperform the TSLS even when the latter is consistent, and is a clear winner in more realistic situations with unobserved aggregate confounders. An important practical message that we infer from these results is that there are reasons to prefer our method even in applications where the TSLS is a valid option.
Our approach addresses the major shortcoming of the conventional TSLS estimation of ((ref)): its invalidity in the presence of unobserved aggregate confounders correlated with the instrument. In practice, there are other reasons why equation ((ref)) can be problematic, for example, nonlinearity or dynamic treatment effects. In these cases, the TSLS might be the wrong tool to start with, and by extension, the same holds for our method. As a result, researchers should use our estimator in applications where the TSLS can be a priori reasonable, but they are worried about potential aggregate confounders. These concerns are not unique to our setting and equally apply to methods that aim to relax the DiD assumptions; see arkhangelsky2023causal for a discussion.
Our setup is related to the literature on invalid instruments (e.g., andrews1999consistent,kolesar2015identification,lewbel2012using,windmeijer2019use). In contrast to this literature, we do not need to assume that a sufficient number of instruments or their known combination, such as average, is valid. Instead, we focus on situations with unobserved aggregate confounders, which allows us to impose a factor structure on the omitted-variable bias. This structure connects our setup to the literature on interactive fixed effects (e.g., bai2009panel,moon2015linear), which focuses on exogenous regressors and thus is not directly relevant to applications we have in mind. An important exception is moon2017dynamic, which considers endogenous regressors but does not allow for low-rank instrumental variables we are interested in. More substantially, we depart from this literature by changing both the estimator and the approach to the analysis. First, we do not estimate the underlying factors but look for the appropriate aggregation scheme. Recent theoretical results in the SC literature illustrate the benefits of this approach even in environments where direct factor-based estimators are available imbens2023identification. Second, we establish the properties of our estimator in the finite population framework with aggregate uncertainty.
Our model is also related to the recent econometric literature on shift-share designs (adao2019shift,borusyak2022quasi,goldsmith2020bartik,jaeger2018shift). First, our finite population analysis is motivated by the approach in adao2019shift and the related earlier work by abadie2023should. Second, similar to this literature, we consider situations where an instrument has a particular product structure. However, our goal is quite different: we propose and analyze a new estimator, while this literature has focused on the properties of the standard IV estimator under alternative assumptions. Crucially, we relax the exogeneity assumption made in the shift-share literature and allow for unobserved aggregate variables that affect different units differently.
The paper proceeds as follows: in Section (ref), we discuss the mechanics of TSLS regression ((ref)) in more detail, present our algorithm, apply it to nakamura2014fiscal, and discuss informally when we expect it to be valid. In Section (ref), we introduce the causal model along with statistical restrictions, and we demonstrate the formal properties of our algorithm. Section (ref) discusses possible extensions of our algorithm, heterogeneous treatment effects, and connections to the literature on shift-share designs. Section (ref) demonstrates the properties of our estimator in simulations, and Section (ref) concludes.
We use $\mathbb{E}[\cdot]$ and $\mathbb{V}[\cdot]$ to denote expectation and variance operators, respectively. We use $\|\cdot\|_2$ to denote the $l_2$-norm of a vector, and $\| \cdot\|_{op}$ to denote the operator norm of a matrix. We use $\textbf{tr}\left(A\right)$ to denote the trace of a square matrix $A$. For two sequences $a_k$ and $b_k$, we write $a_k \lesssim b_k$ if $\frac{a_k}{b_k}$ is bounded and $a_k \sim b_k$ if $a_k \lesssim b_k$ and $b_k \lesssim a_k$. We use $O_p(1)$ and $o_p(1)$ for sequences of random variables that are bounded in probability and converge to zero in probability, respectively.
This section introduces our estimator in the context of an empirical example. We describe the data from nakamura2014fiscal and replicate their baseline results. Next, we propose a new estimator and apply it to this dataset. Our estimates are larger in magnitude, though still within the range reported by nakamura2014fiscal in various specifications. We tie our method to a particular econometric model in Section (ref), where we establish its theoretical properties. We also demonstrate the performance of our method in simulations in Section (ref).
In nakamura2014fiscal, the authors investigate the relationship between government spending and state GDP growth. They use state data on total military procurement for 1966 through 2006 and combine it with U.S. Bureau of Economic Analysis state GDP and state employment datasets. The authors complement these data with the oil prices data from the St. Louis Federal Reserve’s FRED database, the state-level inflation series constructed by del1998aggregate, and their own inflation calculations after 1995.
By estimating the growth-spending relationship, nakamura2014fiscal want to capture the open-economy fiscal multiplier. The authors compare different U.S. states and study their reaction to fluctuations in aggregate military spending. They argue that this strategy allows them to control for other aggregate variables (such as monetary policy). It also allows them to account for the potential endogeneity of local procurement spending.
To illustrate their approach, we introduce some notation. For a generic observation -- a state $i$, and a generic period $t$, -- denote per capita output growth in state $i$ from year $t-2$ to $t$ by $Y_{it}$. Similarly, denote two-year growth in per capita military procurement spending in state $i$ and year $t$, normalized by output in year $t-2$ by $W_{it}$. Finally, let $Z_t$ be the change in total national procurement from year $t-2$ to $t$. The final dataset has $n =51$ states and $T =39$ periods.
The main object of interest -- the fiscal multiplier -- is estimated using the TSLS regression ((ref)) with $D_iZ_t$ as the instrument. The authors construct $D_i$ by estimating $n$ state-level first-stage OLS coefficients
and setting $D_i := \hat \pi_{i}$. As expected, for 49 states, $D_i$ is positive, with Mississippi and North Dakota being the exceptions. In the analysis below, we drop these states and the state of Alaska, where the output growth is exceptionally responsive to the changes in national procurement. This choice leaves us with $n=48$ states.\footnote{Results for the whole sample are similar in magnitude but are estimated less precisely. We report the full sample results in Appendix (ref).}
The TSLS estimator for $\tau$ is equal to
and can be interpreted in two different ways. First, it is a combination of the state-level coefficients,
where $\hat \delta_{i}$ is the state-specific reduced-form coefficient, analogous to $\hat \pi_{i}$. Panel A of Figure (ref) plots $\{(\hat \pi_i,\hat\delta_i)\}_{i\le n}$, with the size of each point being proportional to $\left|D_i - \frac{1}{n}\sum_{j\le n}D_j\right|$, and the colors reflecting the sign. Representation ((ref)) shows that $\hat \tau_{TSLS}$ is equal to the slope of the line that connects the centers of mass of points with negative and positive weights (blue triangles). We see that the coefficients vary a lot, but the association is positive, which results in $\hat \tau_{TSLS} = 1.23$.
Alternatively, $\hat \tau_{TSLS}$ is numerically equal to an IV estimator for the aggregate model
where $Y_t := \frac{1}{n}\sum_{i\le n}Y_{it}(D_i - \frac{1}{n}\sum_{j\le n}D_j)$ and $W_t :=\frac{1}{n}\sum_{i\le n}W_{it}(D_i - \frac{1}{n}\sum_{j\le n}D_j)$, and we use $Z_t$ as an instrument. Panel A of Figure (ref) shows the time-series interpretation of $\hat \tau_{TSLS}$, plotting the aggregate data $Y_t$ and $W_t$ vis-a-vis the OLS fit based on $Z_t$. Using the residuals from these regressions, we produce the robust standard error estimate for $\hat \tau_{TSLS}$, resulting in $\hat {\text{s.e.}}(\hat \tau_{TSLS}) = 0.70$. This estimate is based on using Algorithm (ref), which we explain in more detail in Section (ref). The estimates and standard errors are different from the baseline specification in nakamura2014fiscal ($1.43$ and $0.36$, respectively) because we drop the three states and compute a standard error based on aggregate uncertainty instead of clustering at the unit level.\footnote{To construct this error, we use the automatic model selection ARIMA package in R to estimate the model for $Z_t$. See Section (ref) for details.}
Representation ((ref)) shows that $\hat \tau_{TSLS}$ is a weighted combination of the unit-level coefficients $\{(\hat \pi_i,\hat\delta_i)\}_{i \le n}$ with weights proportional to $D_i - \frac{1}{n}\sum_{j\le n}D_j$. These weights sum up to zero, meaning that the TSLS estimator subtracts the weighted average of the units with relatively large exposures from those with relatively small ones. This action is reflected in Panel A of Figure (ref), where we use different colors for states with positive and negative weights.
This particular aggregation scheme is a consequence of the two-way model in (ref). We eliminate the time fixed effects by averaging over the cross-sectional dimension with the weights that sum up to zero. In applications, these effects capture unobserved aggregate confounders potentially correlated with $Z_t$. For example, in nakamura2014fiscal time fixed effects are meant to capture other policy variables that are likely correlated with national procurement.
This strategy is appropriate only if potential confounders affect all cross-sectional units in the same way (or, at least, in a way that is unrelated to $D_i$). Thus, the main threat to the validity of the TSLS estimator is the presence of aggregate variable $H_t$ with heterogeneous coefficients. A natural model that reflects this is given by
where $\epsilon_{it}$ is assumed to be uncorrelated with $D_i$ and $Z_t$. As long as $\theta_i$ is correlated with $D_i$ and $H_t$ is correlated with $Z_t$, the TSLS estimator suffers from the omitted-variable bias (OVB) and is invalid. Researchers recognize this threat (e.g., see discussions in chodorowstock,guren2018housing) and address it by including additional aggregate and unit-specific control variables. Of course, in practice, we cannot guarantee that these controls are sufficient to account for all confounders.
Our estimator complements this strategy by using a more flexible weighting scheme.\footnote{Below, we discuss the basic version of our estimator that does not involve additional controls. We show how to use controls to improve our estimator in Section (ref).} To understand why weighting can help with unobserved confounders, suppose $D_i$ is binary and split all units into two groups accordingly.\footnote{We thank an anonymous referee for suggesting this example.} If the average value of $\theta_i$ varies between these groups, then the TSLS strategy is invalid. However, if there is an overlap in distributions of $\theta_i$ in the two groups, we can correct these differences by reweighting. We cannot follow this strategy directly because $\theta_i$ is unknown. Instead, we can use the observed data to indirectly implement this strategy by searching for weights with specific balancing properties. The algorithm we present next does precisely that.
First, we again compute the state-level first-stage and reduced-form coefficients by estimating the analogs of (ref), but only using data for periods $t =T_0 +1,\dots, T$. Here, $T_0$ is a user-specified parameter, with a default value $T_0 = \ \frac {T}{3}$. We use $\left\{\left(\hat \pi_i^{post},\hat\delta_i^{post}\right)\right\}_{i\le n}$ to denote the corresponding OLS estimates. To estimate the effect, we aggregate these coefficients using weights $\omega^{rob}_i$,
which we define below. Similar to (ref), this estimator is numerically equal to the time-series IV estimator for the equation
where $Y_t^{rob} := \frac{1}{n}\sum_{i\le n}Y_{it}\omega^{rob}_i$ and $W_t^{rob} :=\frac{1}{n}\sum_{i\le n}W_{it}\omega^{rob}_i$, we use $Z_t$ as an instrument, and estimate ((ref)) using data for $t = T_0+1,\dots, T$. Using weights $\omega^{rob}_i$ as opposed to $D_i - \frac{1}{n}\sum_{j\le n}D_j$ is the main conceptual difference between our estimator and the TSLS. Analogously to $\hat \tau_{TSLS}$, one can compute $\hat \tau_{rob}$ by estimating (ref) by the TSLS for periods $t = T_0+1,\dots, T$, using $\omega^{rob}_i Z_t$ as the instrument.
We construct the weights $\omega^{rob}_i$ using the first $T_0$ periods. As discussed before, we want to make units with high values of $D_i$ look similar on average to those with low values of $D_i$. To this end, we residualize the data for each unit with respect to $Z_t$ and look for such a combination. We achieve this goal by solving a quadratic optimization problem:
where $\zeta^2$ is a user-specified regularization parameter, and
As a default value, we use
where $\hat \epsilon^{(y)}$ and $\hat \epsilon^{(w)}$ are $n\times T_0$ matrices of residuals from regressions in ((ref)). The scaling factor $\hat \sigma$ captures the size of the errors and satisfies $\hat \sigma \lesssim \max\{\hat\sigma_{y},\hat\sigma_{w}\} $. The estimation procedure is summarized in Algorithm (ref).
To implement our method, users need to specify two parameters: the size of the learning period $T_0$ and the amount of regularization $\zeta^2$, and we provide the default values for both. Our theoretical results allow for a range of these parameters and only restrict their asymptotic behavior. In practice, we expect the method to be more sensitive to the choice of $T_0$ rather than $\zeta^2$. If the variation in $Z_t$ is particularly large in specific periods, then inclusion (or exclusion) of these periods might meaningfully change the solution of the optimization problem. In such situations, we advise users to report the sensitivity of the final answer to $T_0$. One can use all $T$ periods to learn the weights and estimate the coefficient. There are multiple reasons, though, that make this choice less attractive, and we briefly discuss them in Section (ref). \RestyleAlgo{boxruled} \LinesNumbered
To reevaluate the results from nakamura2014fiscal with our method, we use the original exposures $D_i$, set $T_0 = 10$ (which corresponds to years 1968-1977), and use the default value for $\zeta$ specified in (ref). We then construct $\omega^{rob}_i$ and estimate $\hat \tau_{rob}$ using Algorithm (ref).
Panel B of Figure (ref) plots the cross-sectional representation of our estimator. As before, the points represent the state-level first-stage and reduced-form coefficients but are now estimated using data from 1978 to 2006. The circle size reflects the absolute value of $\omega^{rob}_i$, and the colors reflect the sign. Our estimator $\hat \tau_{rob} = 1.72$ equals the slope of the line connecting two centers of mass (blue triangles) for negative and positive weights. Compared to coefficients from Panel A of Figure (ref), the first-stage coefficients computed for 1981 to 2006 exhibit less variability. By construction, in Panel A, the states with extreme first-stage coefficients have the largest weights (in absolute value), which is no longer the case in Panel B. Aggregating these state-level coefficients, we get a larger multiplier than before, though still within the range nakamura2014fiscal report for alternative specifications.
There are two differences between Panel A and B of Figure (ref). The first is the period we use to construct the state-level coefficients; the second is the weighting scheme. If we only change the period but apply the same weights as before, we get a multiplier of $1.71$, with a standard error of $0.91$. The similarity between point estimates is unsurprising, given visually minor differences between the weights, which we plot in Figure (ref). The differences are mostly in the tails, with the original weights being extreme for several states. As a result, the standard error of the alternative estimator is $26\%$ higher.
We can see this in Panel B of Figure (ref), which demonstrates the time-series representation of our estimator. We plot the aggregate data $Y_{t}^{rob}$ and $W_{t}^{rob}$ vis-a-vis two separate OLS predictions for the years 1968 through 1977 (in blue), and 1978 through 2006 (in red). By changing the aggregation scheme, we reduce the variability: there is an $11\%$ reduction in the standard deviation for $W_t$ and a $24\%$ reduction for $Y_t$. Despite this decrease in variability, the aggregate instrument $Z_t$ remains relevant. Focusing only on 1979--2006 period, the $R^2$ for $W_t$ increases from $54\%$ to $74\%$, and for $Y_t$ it increases from $11\%$ to $20\%$.
The higher -- compared to the original estimate from nakamura2014fiscal -- standard error of our estimator, $\hat {\text{s.e.}}(\hat \tau_{rob}) = 0.74$, is explained by the fact that we compute it using Algorithm (ref) rather than by clustering at the unit level as in the original paper. Note that despite a shorter time span, the standard error is comparable to the one we got for the TSLS estimator.
$\hat \tau_{TSLS}$ and $\hat \tau_{rob}$ rely on aggregating the unit-level first-stage and reduced-form coefficients. As discussed in the introduction, if $Z_t$ is a valid instrument in the sense of imbens1994identification, all these coefficients describe interpretable causal effects. The TSLS and our method build on these coefficients by aggregating them in a particular way. This is why we encourage users to produce analogs of Figure (ref) to investigate the importance of alternative aggregation schemes. In applications where the unit-level coefficients exhibit strong association, the aggregation does not play a major role, resulting in similar estimates. If, however, there is significant variation, then the weighting scheme becomes important.
Our weights are produced by solving the optimization problem ((ref)), and to understand their properties, it is helpful to consider several edge cases. First, if $\zeta$ is equal to infinity, then $\omega^{rob}_i \propto (D_i - \frac{1}{n}\sum_{j\le n}D_j)$, i.e., we get the same aggregate variables as before, in the TSLS. The resulting estimator is similar but not numerically equal to $\hat \tau_{TSLS}$ because we only use periods $t> T_0$ to estimate the coefficients.
To understand what happens when $\zeta \ne \infty$, suppose $D_i$ is binary, $D_i \in \{0,1\}$. As discussed, the original aggregation scheme constructs a difference between average exposed ($D_i = 1$) and not-exposed ($D_i = 0$) units. Aggregation with weights $\omega^{rob}_i$ has a similar flavor but corresponds to taking a weighted average in both groups. This interpretation follows from examining the constraint in ((ref)). Among all possible weighted averages, we select the one that makes aggregate variables $Y_t^{rob}$ and $W_t^{rob}$ as predictable as possible.
Motivation for this choice is evident from looking at the first-stage and reduced-form equations that correspond to ((ref)):
Unobserved confounders make the prediction of $Y_{it}$ and $W_{it}$ by $Z_t$ harder, so the weights that eliminate such factors should also make the prediction easier. Terms $u_{it}^{(y)}$ and $u_{it}^{(w)}$ create a statistical challenge: it is possible that instead of eliminating the confounder, the weights produce a combination of errors that compensates $H_t$. To prevent such overfitting, we include the regularization term in ((ref)), which forces the weights to be as uniform as possible. Using this, we can show that $\omega^{rob}$ are close to deterministic weights $\omega^{\star}$, which minimize
subject to the same constraints. As long as $\frac{1}{n}\sum_{i\le n}\omega^{\star}_i\theta_i^{(w)}$ and $\frac{1}{n}\sum_{i\le n}\omega^{\star}_i\theta_i^{(y)}$ are small, we can expect $\frac{1}{n}\sum_{i\le n}\omega^{rob}_i\theta_i^{(w)}$ and $\frac{1}{n}\sum_{i\le n}\omega^{rob}_i\theta_i^{(y)}$ to be negligible as well.
The expectation in (ref) is not defined until we describe the randomness behind the observed data. One can think of at least three different sources of such randomness (and their combinations): (a) random sampling of units, (b) random variation in the unit-specific errors $(u_{it}^{(y)},u_{it}^{(w)})$, and (c) random variation in aggregate variables $(Z_t, H_t)$. Among these three, we view the last one as the most appropriate for the analysis in nakamura2014fiscal and similar applications. First, the sampling uncertainty is at odds with the fact that we observe all possible states in the US (see the discussion in abadie2020sampling,abadie2023should). Second, probability statements based on $(u_{it}^{(y)},u_{it}^{(w)})$ are disconnected from the original identification idea, which instead relies on the aggregate variation. Guided by this logic, Section (ref) focuses on uncertainty coming from $(Z_t, H_t)$.
In the rest of the paper, we demonstrate the advantages of our algorithm with theory and simulations. However, as with any method, our approach has limitations and should be used carefully. Below, we discuss the main ones, thus defining practical use cases for our estimator.
Our method is designed for applications where the TSLS regression ((ref)) is a priori reasonable, but users worry about potential unobserved confounders. There are multiple reasons why ((ref)) might fail, other than omitted variables. For example, the underlying model can be nonlinear, or the dynamic effects of the past treatments can be sizable. In these cases, the TSLS regression, and by extension, our improvement upon it, might be the wrong tool, and researchers should use other methods. We discuss this issue in more detail in Section (ref).
Our algorithm is data-intensive, particularly in the time dimension. Our formal results require the total number of periods to be large and the number of units $n$ to be at least of a similar order. This assumption can be restrictive in some cases, but we believe it is reasonable for applications currently estimating (ref) using TSLS with aggregate instruments.
Finally, for our weights to improve over the TSLS ones, the environment should be sufficiently stable over time. To construct $\omega^{rob}$, we use the residuals after projecting out $Z_t$, which behave well when the effect of $Z_t$ does not change over time. We return to this point in Section (ref), discussing heterogeneous effects in more detail. For our weights to be useful for the second part of the data, any potential confounder like $H_t$ in ((ref)) should have a similar effect in both data parts. This assumption might be too strong in environments where structural breaks are likely. This limitation can potentially be relaxed by using the whole sample to construct the weights. However, our theoretical results, particularly inference, rely on sample splitting. We also believe that sample splitting is a good general practice that protects from potential abuse.\footnote{See spiess2018optimal for a formalization of this argument.}
Given these limitations, we recommend that researchers use our technique in situations where the number of observed periods is relatively large, structural breaks are unlikely, and the main potential problem is the presence of unobserved aggregate confounders rather than nonlinearity or dynamics. Our formal results and simulations show that our method either dominates the TSLS or performs similarly under this set of assumptions.
This section introduces our formal setup and states three theoretical results.\footnote{All proofs are collected in Appendix (ref).} The first theorem shows that our estimator remains consistent even when the TSLS fails. The second guarantees our estimator is asymptotically unbiased and normal under mild technical assumptions. The third theorem justifies conventional inference in situations with sufficient heterogeneity in the outcomes.
We start this section by discussing a causal model that we use to analyze the properties of our estimator. The model is based on two separate assumptions. The first describes the structural relationship between the observed data and underlying potential outcomes neyman1923applications,rubin1977assignment. The fundamental property of this model is that it allows for both observed and unobserved aggregate inputs. The second assumption describes the joint behavior of the aggregate variables. We close this section by discussing the difference between our setup and the more conventional ones.
We observe $n$ units ($i$ is a generic unit) over $T$ periods ($t$ is a generic period). For each unit, we observe an outcome variable $Y_{it}$, an endogenous policy variable (treatment) $W_{it}$, an aggregate instrument $Z_t$, and a measure of exposure of unit $i$ to this variable $D_i$. We aim to estimate a causal relationship between $Y_{it}$ and $W_{it}$. To formalize causality, we describe a model of the underlying potential outcomes. In addition to $w_{t}$ (potential value of $W_{it}$) and $z_t$ (potential value of $Z_t$), we also introduce $h_t$, an unobserved aggregate confounder that causally affects both the outcome and the treatment variable. We define $w^t := (\dots, w_{1},\dots, w_t)$, $z^t:=(\dots, z_{1},\dots, z_t)$, and $h^t:= (\dots, h_{1},\dots, h_t)$, and we make our first assumption.
The critical part of this assumption is the presence of unobserved aggregate variable $H_t$. The danger such unobservables present for identification is well recognized in applied work chodorowstock,guren2018housing. The standard restriction made in the literature is to assume that $\theta_{i}^{(w)}$ and $\theta_{i}^{(y)}$ do not systematically vary over $i$. We do not make this assumption; instead, we allow for such heterogeneity. Following most empirical applications, we focus on contemporaneous treatment effects and assume that only current quantities affect the outcomes. In practice, this assumption is often justified by looking at particular transformations of the available data nakamura2014fiscal. Finally, to simplify the exposition and focus on a single parameter, we avoid heterogeneity in treatment effects. We relax this in Section (ref), where we discuss when the output of our algorithm can be interpreted as a weighted average of individual treatment effects.
Our next assumption restricts the joint distribution of aggregate variables $\{(Z_t,H_t)\}_{t\le T}$.
The first part of this assumption restricts the means of $Z_t$ and $H_t$, which are assumed to be constant over time. This is without loss of generality for $H_t$ because its mean can be treated as a part of $\alpha_{it}^{(w)}$ and $\alpha_{it}^{(y)}$, but it is restrictive for $Z_t$. This assumption can be relaxed by considering a parametric model for the mean, for example, allowing for seasonality or secular trends. Fundamentally, our approach relies on researchers knowing how to detrend $Z_t$ and exploit random fluctuations $\epsilon_t^{(z)}$, which is straightforward if the mean is constant. This concept is closely linked to borusyak2023nonrandom, who emphasize the importance of demeaning and highlight that the instrument $Z_t$ can be worthless without the ability to demean it.
The second part restricts the distribution of $\epsilon^{(z)}$ and $\epsilon^{(h)}$. These errors are generated by the underlying independent shocks $\nu^{(z)}$ and $\nu^{(h)}$. Since $|\rho_{ag}|$ is less than one, there is variation in $\epsilon^{(h)}_t$ that is not entirely explained by $\epsilon^{(z)}_t$ (and vice versa). Restrictions on the elements of matrices $\Lambda^{(z)}$ and $\Lambda^{(h)}$ exclude persistent cases (e.g., random walks) but allow for time-heterogeneous coefficients, making (ref) particularly appealing. This structure is less flexible than the full linear model:
which would allow $(\epsilon^{(z)},\epsilon^{(h)})$ to be generated by a VAR system with time-heterogeneous coefficients. The additional restrictions we impose with (ref) are needed to have a single parameter $\rho_{ag}$ that governs the correlation, simplifying the exposition and proofs. We also expect our results to hold in the more general model (ref).
The model we describe in this section follows the tradition that goes back at least to neyman1923applications, in which we treat the observed units and their potential outcomes as fixed and focus explicitly on the randomness introduced by external shocks. As discussed in the previous section, this is natural in applications where the observed units exhaust the population of interest, which is often the case with geographic data. The finite population approach is also commonly motivated by the desire to connect the source of randomness to a well-defined shock while keeping the rest of the structure of the model as flexible as possible. This connection is straightforward in experimental applications, with the shock corresponding to the initial treatment assignment, making this type of analysis particularly appealing.
This logic is more complex once we consider observational studies where we cannot explicitly manipulate the treatment. Still, it is natural, especially in environments with aggregate variables. Recall the setup in nakamura2014fiscal discussed in the previous section. In that context, various structural shocks in the US economy (in the sense of ramey2016macroeconomic) affect local-level outcomes through observed and unobserved aggregate variables. Assumption (ref) captures this behavior by modeling $Z_t$ and $H_t$ as functions of the underlying structural shocks $(\nu^{(z)},\nu^{(h)})$. These shocks turn $Z_t$ and $H_t$ into exogenous shifters that, per Assumption (ref), serve as inputs for the economic system that determines the local-level variables. The view that aggregate variables are a natural source of quasi-experimental variation is now widely adopted in the methodological literature on shift-share instruments adao2019shift,borusyak2022quasi,borusyak2023nonrandom. The key assumption distinguishing our model is the possible correlation between the observed and unobserved aggregate variables (e.g., Section 4.2 in adao2019shift).
As another example suppose $(Y_{it},W_{it})$ are determined jointly in the local equilibrium:
This structure arises, for example, in guren2018housing where $Y_{it}$ is the retail employment in location $i$, period $t$, and $W_{it}$ is the house price. The aggregate variables $Z_t$ and $H_t$ correspond to exogenous demand and supply shifters. Substituting $Y_{it}$ in the expression for $W_{it}$, we get the model ((ref)). This example demonstrates the difference between $Z_t$ and $H_t$. The former is a shifter for $W_{it}$ and is excluded from the structural equation for $Y_{it}$. Despite this exclusion restriction, $Z_t$ might be an invalid instrument due to its potential correlation with $H_t$.
The example above also illustrates the nature of the endogeneity problem we aim to solve with our method. The outcomes $(Y_{it}, W_{it})$ are determined in the local equilibrium, and thus any direct regression of $Y_{it}$ on $W_{it}$ or vise versa does not have any economic meaning. This feature describes the identification problem in the model, which is unaffected by any sampling assumptions behind $\alpha_{it}^{(w)},\alpha_{it}^{(y)}$, or any separation of these quantities into systematic fixed effects, and idiosyncratic errors, as in (ref). The key assumption that makes the relationship between $Y_{it}$ and $W_{it}$ meaningful is the exclusion restriction for $Z_t$.
It is also instructive to look at our setup from a more conventional point of view. To this end, consider $Y_{it}(w):= Y_{it}(w,H_t)$ and $W_{it}(z) :=W_{it}(z,H_t)$ -- potential outcomes at realized values of $h_t$. For each $i$, we can view the potential outcomes $\{Y_{it}(w), W_{it}(z)\}_{t\le T}$ as a time-series version of the IV model of imbens1994identification. In particular, Assumption (ref) guarantees that the instrument $Z_t$ satisfies the exclusion restriction for each unit. Assumptions (ref) and (ref), however, do not guarantee that $Z_t$ is independent of $(Y_{it}(w), W_{it}(z))$. The two assumptions together quantify the extent of this dependence at the unit level:
Our setup thus relaxes the independence assumption of imbens1994identification (Condition 1, (i) in the paper) but imposes a product structure on the correlation. Equation ((ref)) emphasizes the difference between $H_{t}$ and $(\alpha_{it}^{(w)},\alpha_{it}^{(y)})$: heterogeneity in the latter does not lead to a violation of the independence assumption and thus is not crucial for identification.
By definition, we do not observe $H_t$ and thus cannot guarantee that it is one-dimensional, potentially making a model with multiple confounders more appropriate. Notationally, this extension is straightforward since we can interpret $(\theta_i^{(y)},\theta_i^{(w)})$ and $H_t$ as $p$-dimensional vectors. When $p$ is large enough, the RHS of ((ref)) can approximate arbitrary covariance between $Y_{it}(w), \, W_{it}(z)$ and $Z_t$. The dimension of $H_t$ does not directly enter Algorithm (ref) but makes its analysis more involved. We focus on the one-dimensional case because it transmits our theoretical insights in the simplest form. In Appendix (ref), we establish our main bound, assuming $H_t$ is a vector, and later specialize it to the scalar case.
This section describes two assumptions that restrict heterogeneity in potential outcomes. As discussed in the previous section, in the absence of $H_t$, one can cast our setup into the conventional IV model from imbens1994identification, which does not require substantial restrictions on heterogeneity in potential outcomes. To understand why restrictions are necessary when $H_t$ is present, suppose $\theta_{i}^{(w)} \propto \pi_i \propto \theta_i^{(y)}$ and researchers observe $\pi_i$ up to proportionality. In each period, we can project the data on $\{\pi_i\}_{i\le n}$ and its orthogonal complement in $\mathbb{R}^n$. The latter projection depends neither on $Z_t$ nor $H_t$, because $\theta_{i}^{(w)},\theta_i^{(y)}$ are proportional to $\pi_i$. As a result, it is useless for identification without additional assumptions on terms $\alpha_{it}^{(w)}, \alpha_{it}^{(y)}$. The projection on $\{\pi_i\}_{i\le n}$ leaves us with a single time series for the outcome and the endogenous variable. By the proportionality assumption, these time series are necessarily corrupted by the aggregate confounder, making them equally useless for identifying the effect.
The exact proportionality of $\pi_i, \theta_{i}^{(w)},\theta_i^{(y)}$ is an extreme restriction, which is unlikely to hold in any practical application. Once we relax this restriction, the identification results are much more positive. As long as the projection of $\pi_i$ on $\theta_{i}^{(w)}$ and $\theta_i^{(y)}$ is not equal to zero, we can use the residuals from this projection to aggregate units cross-sectionally. This aggregation guarantees that the resulting quantity depends on $Z_t$ but not on $H_t$, thus solving the OVB problem. The conditions we impose below essentially amount to this possibility with two important caveats. First, to implement the strategy described above, we need to know $\{\pi_i\}_{i \le n}$ and $\left\{\theta_{i}^{(w)},\theta_i^{(y)}\right\}_{i \le n}$ up to proportionality, and we do not know these quantities. We assume that $\pi_i$ is directly connected to the observed variation $D_i$ to bypass this issue. Second, we need to guarantee that the terms $\alpha_{it}^{(w)}, \alpha_{it}^{(y)}$ do not “mask” important variation, the point we discuss in detail below.
We start by connecting $\pi_i$ to $D_i$.
This assumption guarantees that the observed variation in the unit-level “exposure” $D_i$ is directly connected to the underlying structural parameter $\pi_i$. The linearity is motivated by the empirical practice where researchers often assume that exposures $\pi_i$ are known up to linear transformation (e.g., dube2013commodity,nunn2014us). It can be substantially relaxed: fundamentally, our theoretical results rely on $D_i$ being strongly correlated with $\pi_i$ after adjusting for other unit-specific coefficients. As a result, one can extend Assumption (ref) by explicitly including $\theta_i^{(w)},\theta_i^{(y)}$ or functions of $\{\alpha_{it}^{(w)}, \alpha_{it}^{(y)}\}_{t\le T}$ in the expression for $\pi_i$.
To state our next assumption, we introduce additional notation. For any $T_a >1$, define population analogs of $\hat \sigma_{k,T_0}$ from ((ref)):
For any $T_b > T_a\ge1$, $t\in [T_a,T_b]$, weights $\omega_i$ such that $\sum_{i \le n} \omega_i = 0$, and $k \in \{y,w\}$ define a transformation of the model coefficients:
With this notation, we are ready to introduce our following assumption.
To unpack this condition, first suppose that $D_i = \pi_i$ and $\alpha_{it}^{(w)} = \alpha_{it}^{(y)} = 0$. In this case, Assumption (ref) is satisfied whenever $\pi_i$ is not perfectly predictable by $\theta_i^{(w)},\theta_i^{(y)}$. As we discussed above, this restriction is necessary for identification in our model. When $\alpha_{it}^{(k)} \ne 0$, Assumption (ref) imposes additional restrictions. As we discussed in Section (ref), our statistical results guarantee that the weights $\omega^{rob}$ are close to deterministic oracle weights $\omega^{\star}$ that optimize the expected version of ((ref)):
subject to appropriate constraints. Using Assumptions (ref) -- (ref) we can compute these expectations, and after concentrating $\left\{\eta_{0}^{(k)},\eta_{z}^{(k)}\right\}_{k \in\{y,w\}}$, we get
where $\kappa^2(T_0)$ is strictly positive. Assumption (ref) guarantees that the oracle problem ((ref)) has a well-behaved solution.
Assumption (ref) is a sufficient condition directly motivated by our approach. How restrictive are the additional restrictions we impose on top of those necessary for identification? We do not attempt to answer this question fully and instead ask if our sufficient condition holds in statistical models commonly used in empirical and methodological literature. These models put probability structure on unit-specific quantities that we treat as fixed elsewhere. We use this structure to argue that Assumption (ref) holds for a “generic” dataset, i.e., for most of the datasets generated from the models we consider.
This example covers many familiar cases. First, if $\alpha_i^{(d)} = \alpha^{(d)}$, then $D_i$ is as good as randomly assigned -- a situation that rarely holds in applications but serves as a natural benchmark. In this situation, the TSLS weights satisfy Assumption (ref). Importantly, this is no longer true if we allow heterogeneity in $\alpha_i^{(d)}$. If $L_{it}^{(k)} \equiv 0$, then we recover a conventional two-way model that is commonly used in applications. In practice, we rarely expect the two-way model to hold exactly, and $L_{it}^{(k)}$ can be viewed as an approximation error. If $L_{it}^{(k)}$ has a low-rank structure, then we recover the interactive fixed-effects model (e.g., bai2009panel). However, this structure is not necessary, and $L_{it}^{(k)}$ can vary arbitrarily as long as it remains appropriately bounded.
The result in Proposition (ref) relies on the existence of random $\epsilon_i^{(d)}$, suggesting that this variation plays an important identification role. As we explained at the beginning of this section, identification is impossible in the absence of the residual variation in $\pi_i$ that is orthogonal to $\theta_i^{(w)}, \theta_i^{(y)}$, and Proposition (ref) generates such variation using random errors. Two important details make this variation different from the one in $Z_t$. First, the result in Proposition (ref) holds for a generic parameter realization, allowing us to treat them as fixed elsewhere. Second, the fluctuations $\epsilon_{i}^{(d)}$ are not available to the researcher, because $\alpha_i^{(d)}$ is unknown. Moreover, our analysis does not rely on empirical weights $\omega^{rob}$ converging to $\omega^{\star}_{T_a}$, instead they are converging to a certain projection of these object, see Theorem (ref) for the precise statement.
The setup of Proposition (ref) is motivated by the statistical models commonly used in applications. The class of models it covers is large but by no means exhaustive, and in some applications, one would need to consider a different justification for Assumption (ref). For example, if unit exposures $(\pi_i, \theta_i^{(y)}, \theta_i^{(w)})$ are equilibrium objects in some economic model, then to guarantee Assumption (ref) one needs to restrict the underlying structural parameters of that model. Proposition (ref) might not directly apply to such cases, but we expect the underlying intuition to be helpful.
It is also instructive to compare Assumption (ref) to conventional restrictions from the interactive fixed effects literature, particularly for models with low-rank regressors. For example, Assumption ID (ii) -(iii) in moon2017dynamic requires the terms $\alpha_{it}^{(k)}$ to be uncorrelated with other variables in the model and guarantees the presence of residual variation in the analog of $D_i$ after it is projected on the analogs of $\theta_i^{(w)}, \theta_{i}^{(y)}$ (without explicitly modeling the source of this variation). These restrictions imply the analog of Assumption (ref), because we can use the residual variation to construct $\omega^{\star}_{T_a}$. The model in Proposition (ref) is more general because it allows the terms $\alpha_{it}^{(k)}$ to have a systematic component $L_{it}^{(k)}$ with essentially unrestricted spectral properties (apart from boundedness). See also imbens2023identification for a related discussion.
We now turn to the statistical properties of our estimator. All probability statements in the section except those in Proposition (ref) refer to the joint distribution of $\{(Z_t, H_t)\}_{t\le T}$. We focus on a particular asymptotic regime characterized by the following assumption.
With the first part of this assumption, we restrict the analysis to environments where $n$ is comparable to or larger than $T$, which we expect to hold in many applications. The second part implies that the variability in $D_i$ and $\theta_i^{(k)}$ is present in the limit. For binary $D_i$, this assumption is reasonable if the size of the treated or control group is not too small. The variability in $\theta_i^{(k)}$ implies that $H_t$ is a “strong” factor. While common in the theoretical literature on interactive fixed effects (e.g., bai2009panel,moon2015linear), this assumption can be restrictive in some applications, where researchers expect little variability in $\theta_i^{(k)}$. A version of our results holds in environments where $\sigma^2_{\theta^{(k)}} = 0$, thus allowing for weak factors (see the discussion in Appendix (ref)). The restriction on the correlation is innocuous, as it always holds along a subsequence, and we make it to simplify the exposition. Finally, the restriction on $\alpha^{(k)}_{it}$ guarantees that $Y_{it}$ and $W_{it}$ have finite variances.
Assumption (ref) describes the limit behavior of unit-specific quantities. For the aggregate variables, we define similar objects for fixed $T_b > T_a \ge 1$ and $k\in \{z,h\}$:
As indicated by the indices, these quantities depend on $T_a$ and $T_b$. Assumption (ref) guarantees that $|\rho_{T_a|T_b}| <1$, and $\sigma_{k,T_a|T_b}$ are uniformly bounded from above and below.
Our first result compares probability limits of $\hat \tau_{TSLS}$ and $\hat \tau_{rob}$.
This result demonstrates that $\hat \tau_{rob}$ remains consistent in the regime where $\hat \tau_{TSLS}$ generally fails. It also formalizes a part of the discussion in Section (ref). In our model, the only threat to the validity of the TSLS is the presence of omitted variables. As long as either $\rho_{cs}^{(y)}$ or $\rho_{1|T}$ is equal to zero, the TSLS estimator is consistent. In particular, the latter quantity equals zero as long as $\rho_{ag} = 0$.
Theorem (ref) provides a first justification for using our estimator, but it does not describe its distributional properties in large samples. Under current assumptions, we can provide only relatively weak guarantees on the asymptotic behavior of $\hat \tau_{rob}$. In particular, Assumptions (ref) and (ref) imply that there exist weights such that
however, this property is too weak to achieve asymptotic unbiasedness. To make progress, we impose additional assumptions on $\theta_i^{(y)}$ and $\theta_i^{(w)}$.
This restriction is similar to Assumption (ref) and requires existence of variation in $\theta_i^{(k)}$ that is not captured by $\alpha_{it}^{(w)},\alpha_{it}^{(y)}$, and $\theta_{i}^{(-k)}$. To justify it, we return to the example from Proposition (ref).
To understand why Assumption (ref) can improve the performance of $\hat \tau_{rob}$, it is useful to consider environments where it fails. In particular, if $\theta_i^{(k)}$ is nearly spanned by $\{\alpha_{it}^{(y)},\alpha_{it}^{(w)}\}_{t \le T_0}$ and $\theta_i^{(-k)}$, but the remaining variation is strongly associated with $\pi_i$, then it is very hard to eliminate it by aggregation, which results in a slow rate of convergence. This problem is mitigated when there is enough variability in $\theta_i^{(k)}$, which is not explained by other variables. Our next result demonstrates these gains by characterizing the asymptotic behavior of $\hat \tau_{rob}$. To state it, we define for arbitrary periods $T_b > T_a>0$ a matrix $\Lambda^{(z)}_{T_a|T_b}$ such that
We also use $\alpha^{(k)}_{T_a, T_b}(\omega)$ to denote the vector $ \left(\alpha_{T_b, T_a|T_b}^{(k)}(\omega),\dots, \alpha_{T_a, T_a|T_b}^{(k)}\right)$.
This result implies that our estimator is asymptotically unbiased and normal as long as $\sigma_{n,T}$ remains bounded. This condition can be restrictive if $\alpha_{i,t}^{(y)}$ are close to being degenerate. For example, if we view these quantities as being generated from the normal distribution, i.e., $\alpha_{it}^{(y)} \sim \mathcal{N}\left(\alpha^{(y)}_i + \lambda_t^{(y)},\sigma^2_{\alpha^{(y)}}\right)$, then $\sigma_{n,T} =O_p\left(\frac{1}{\sqrt{n}}\right)$ and higher-order terms in ((ref)) become important. This lack of uniformity is similar to one analyzed in menzel2021bootstrap. In practice, we do not expect the two-way model to hold exactly; rather we view it as an approximation. In this situation, the first term in ((ref)) becomes dominant, and we can use Theorem (ref) for inference.
Theorem (ref) describes the asymptotic behavior of our estimator in the presence of unobserved confounders. If no such variables exist in the structural equation, i.e., $\theta_i^{(y)} \equiv 0$, then our estimator and the standard TSLS estimator are asymptotically normal under mild technical conditions. In this regime, $\sigma_{n, T}$ can be smaller or larger than its TSLS counterpart $\left\| \alpha_{1|T}^{(y)}(\omega^{TSLS}) \Lambda^{(z)}_{1|T}\right\|_2$, depending on the underlying complexity of the potential outcomes and differences in the sample sizes.
We conduct inference in several steps, which we summarized in Algorithm (ref). First, we estimate the variance. We assume that a researcher has access to a consistent estimator for $\Lambda^{(z)}_{T_0+1|T}$, which we denote as $\hat\Lambda^{(z)}_{T_0+1|T}$. In practice, this means that we need to estimate the model for $Z_t$ to conduct inference. In our simulations and empirical analysis, we use off-the-shelf methods to do that, and the resulting estimator performs relatively well.
For $t>T_0$, we construct scaled residuals from the aggregate regression
and estimate the asymptotic standard error of $\hat\tau_{rob}$:
With this quantity, we construct a standard asymptotic confidence interval of level $1-\alpha$:
where $z_{\alpha}$ is $\alpha$-quantile of the standard normal distribution. Our next result characterizes the asymptotic properties of this interval.
\RestyleAlgo{boxruled} \LinesNumbered
As discussed above, we expect this theorem to be useful in practice whenever the two-way model for $\alpha_{it}^{(k)}$ is only approximately correct. This result focuses on the conventional interval ((ref)), which is valid if the first stage is strong, that is, if $\eta_{\pi}$ is large enough. A version of Theorem (ref) also holds for the first-stage and reduced-form coefficients and thus can be used to conduct conventional robust inference (see andrews2019weak for a recent survey on robust inference).
To construct $\hat \sigma_{rob}$, we combine aggregate residuals with the estimator for the parameters of the design model. If $\Lambda^{(z)}$ is diagonal, i.e., the variation in $Z_t$ is independent over time, then $\hat \sigma_{rob}$ corresponds to “clustering at time level.” However, in practice, this assumption can be too restrictive. With general $\Lambda^{(z)}$, we need to consider dependence over time, and $\hat \sigma_{rob}$ does that using $\hat \Lambda^{(z)}_{T_0+1|T}$. In the application, we discussed in Section (ref) as well as in the Monte-Carlo experiments in Section (ref), we use the standard automatic model selection ARIMA package in R to estimate the model for $Z_t$ and then simulate from this model to construct $\hat \sigma_{rob}$. We recommend that researchers use the same approach in practice.
This section discusses three possible extensions of our model and the respective adjustments to the algorithm. We first show how to incorporate covariates in our setting. Secondly, we examine the case of heterogeneous treatment effects as a natural extension. We conclude Section (ref) by connecting our estimator to the literature on shift-share designs.
A typical regression equation estimated in applications will have a more complicated structure than ((ref)):
Here, $X_i$ are observed unit-level attributes, for example, region indicators, and $\tilde H_t$ is a vector of observed aggregate variables that we expect to be correlated with $Z_t$. Equation ((ref)) is estimated by the TSLS using $D_i Z_t$ as an instrument for $W_{it}$ and treating $\alpha_i$ and $\theta_i$ as fixed parameters. Inclusion of $\mu_t(X_i)$ instead of $\mu_t$ and $\tilde \theta_i^\top \tilde H_{t}$ in the equation mitigates the OVB concerns but does not eliminate them.
Our estimator also allows for unit-level covariates and observed aggregate variables.\footnote{To incorporate time-varying covariates $X_{it}$, we can define $X_i := (X_{i1},\dots, X_{iT})$. Alternatively, and more in line with current empirical practice, we can instead residualize $Y_{it}$ and $W_{it}$ with respect to $X_{it}$.} In particular, we suggest estimating equation ((ref)) by TSLS using $\omega^{rob}_iZ_t$ as instrument for $W_{it}$ and data from periods $T_0+1, \dots T$. The weights $\omega^{rob}_i$ then solve an adjusted optimization problem:
The additional constraint guarantees that aggregation eliminates the linear projection of $\theta_i^{(w)}$ and $\theta_i^{(y)}$ on $X_i$. As a typical example, consider a situation where data can be grouped into clusters and where researchers wish to include cluster-specific time fixed effects. This can be achieved using $X_i$ corresponding to cluster indicators, under the natural extension of Assumptions (ref) to (ref), Theorems (ref) and (ref) continue to hold for the weights that solve ((ref)).
In some applications, the unit-level variables $Y_{it}$ and $W_{it}$ have different statistical properties; for example, they are measured using a different number of observations. In such situations, researchers commonly use weighted versions of the TSLS. To achieve the same with our algorithm, researchers can estimate ((ref)) using weighted TSLS with $\omega^{rob}_iZ_t$ as an instrument for $W_{it}$. To construct the weights $\omega^{rob}_i$, we solve the optimization problem ((ref)) but instead of the standard euclidean norm $\|w\|_2^2$, we use a weighted one:
where $A$ is a diagonal matrix, and $(A)_i = a_i^2>0$.
In applications, it is rarely possible to argue that the treatment effects are constant, and thus Assumption (ref) can be too restrictive. To address this, we consider a model with heterogeneous effects:
We also define the reduced form that corresponds to the structural equation above:
where $\tilde\alpha_{it}^{(y)} := \alpha_{it}^{(y)} + \tau_i \alpha_{it}^{(w)}$, and $\tilde \theta_i^{(y)}:= \theta_i^{(y)} + \tau_i \theta_i^{(w)}$. For any estimator $\hat \tau(\omega)$ that averages units with arbitrary weights $\omega$ and constructs the IV ratio from the aggregate regressions, we have
Our goal in this section is to understand when $\tau(\omega):=\frac{\frac{1}{n}\sum_{i\le n} \omega_i\tau_i\pi_i}{\frac{1}{n}\sum_{i\le n} \omega_i\pi_i}$ has a causal interpretation. We thus ignore the errors in ((ref)). Their properties depend on the choice of weights $\omega$ and can be established in the same way as before.
First, we consider a situation where $\pi_i = \eta_{\pi} D_i$, for binary $D_i\in\{0,1\}$. For $\omega^{TSLS}$, we get
which is an average treatment effect for the exposed group. Using $\omega^{rob}$, we get
where $\frac{1}{n}\sum_{i \le n}\omega^{rob}_iD_i =1$. Without additional restrictions, we cannot interpret $\tau_{rob}$ as a convex combination of treatment effects because the weights $\omega^{rob}_i$ can be negative for exposed units. Negative weights lead to extrapolation, which can help with the OVB but at the cost of interpretability.
This problem is easy to address by adding a non-negativity constraint
to the optimization program ((ref)). The resulting $\tau(\omega^{rob})$ is a convex combination of treatment effects by construction. The optimization problem remains convex and can be solved efficiently, even for large datasets. Inequality constraint ((ref)) also acts like a powerful regularizer, improving the statistical properties of the algorithm. To reap these benefits, we need to assume that “good” balancing weights that satisfy ((ref)) exist, for example, by adding this restriction to Assumption (ref). Overall, in applications with binary $D_i$, we recommend imposing ((ref)) unless the user strongly believes that extrapolation is necessary.
Many applications do not have a control group with $D_i$ taking arbitrary values, so non-negativity constraints are harder to motivate. However, one can still interpret $\tau(\omega^{rob})$ with additional assumptions. In particular, suppose that
where $\epsilon_i^{(d)}$ has the same properties as in Proposition (ref). If Assumption (ref) holds, we have
As long as $\omega^{rob}$ converges to $\frac{\epsilon^{(d)}}{\sigma^2_{d}}$, the estimand $\tau(\omega^{rob})$ converges to the average treatment effect. In Appendix (ref), we discuss models where this convergence holds. In this case, our method improves over the TSLS in two ways: it removes the OVB and helps interpretability.
The heterogeneity we consider in this section is restricted in an important way: we do not allow $\tau_i$ and $\pi_i$ to vary over time. Such variation makes it impossible to project $Z_t$ out when constructing the weights $\omega^{rob}$. This problem can be bypassed if the researcher knows that in the initial $T_0$ periods, $\pi_{it} \equiv 0$ for all units. In particular, in applications where $\pi_{it} = (\eta_0+\eta_tD_i)\mathbf{1}_{t>T_0}$, we expect Algorithm (ref) to perform well with both cross-sectional and time-series heterogeneity in treatment effects, as long as $\eta_t >0$.
This section discusses the relationship between our model and models from the shift-share, or Bartik instruments, literature (adao2019shift, borusyak2022quasi, goldsmith2020bartik). We start by considering an extension of our original framework. Assume that instead of a single aggregate variable, we have $|S|$ of them. In a typical application, these will correspond to industry-level shifters. The following equations are now satisfied for all $i$ and $t$:
where $s$ is a generic industry, and we observe $\{\gamma_{is}\}_{i,s}$, $\{W_{it},Y_{it}\}_{it}$, $\{Z_{ts}\}_{t,s}$, and $\sum_{i\le n} \gamma_{is} =1$. It is straightforward to see that our model is a special case of this with $|S| = 1$.
The model typically considered in the shift-share literature is a special case of ((ref)), with $T=1$, and two additional assumptions: (a) $Z_{ts} = \psi_{ts}^{\top}\mu_t + \epsilon_{ts}$, where $\psi_{ts}$ are known, $\mathbb{E}[\epsilon_{ts}] = 0$, and $\epsilon_{ts}$ are uncorrelated over $s$; and (b) for every $t$, $\{H_{ts}\}_{s\in S}$ is uncorrelated with $\{\epsilon_{ts}\}_{s\in S}$. Identification is achieved by exploiting variation over industries (see borusyak2022quasi). In applications, $T$ is usually not equal to $1$, and the model in differences is often considered. At the same time, the identification argument does not exploit the time dimension and focuses on the variation over industries.
Models of the type ((ref)) can be promising because they allow for a combination of two identification arguments: one based on the variation over time and one based on the variation over $s$. In applications, both $|S|$ and $T$ can be modest (especially if we want shifters to be independent over $s$), and thus it is natural to use both sources of variation. The development of strategies that allow for that is an attractive area of future research.
In this section, we illustrate the performance of our estimator in simulations. To make the simulations more realistic, we build them using the dataset from nakamura2014fiscal, which we described in Section (ref). In our experiments, we try to capture the spirit of this empirical exercise and investigate how different features of the data-generating process affect the performance of the algorithms. Formally, our simulations are based on the following model:
Here, parameters $\{\beta_i^{(y)},\beta_i^{(w)},\mu_t^{(y)},\mu_t^{(w)},L_{it}^{(y)},L_{it}^{(w)},\tau, \pi_i, \theta_i^{(w)}, \theta_i^{(y)}\}_{i\le n,t\le T}$ are fixed, while $\epsilon_{it}^{(y)},\epsilon_{it}^{(w)}$, and $\{Z_t,H_t\}_{t\le T}$ are random.
In Appendix (ref), we describe how exactly we use the data to construct $\{L_{it}^{(y)},L_{it}^{(w)}, \pi_i\}_{i\le n,t\le T}$, and the models for $\{Z_t\}_{t\le T}$ and $\{\epsilon_{it}^{(y)},\epsilon_{it}^{(w)}\}_{i\le n,t\le T}$. Heuristically, we extract the components $L_{it}^{(y)}$ and $L_{it}^{(w)}$ using the SVD decompositions of observed data, while for $\pi_i$, we use the estimated $\hat \pi_{i}$ from Section (ref), which we scale to make the instrument relatively strong.\footnote{The median $F$ statistic for $\hat \tau_{TSLS}$ for the fourth design is equal to $78$.} We make these adjustments to focus on the properties of our estimator in the regime covered by Theorems (ref) and (ref). The data are not directly informative about $H_t$ and $\{\theta_{i}^{(w)},\theta_{i}^{(y)}\}_{i\le n}$---we need to make ad hoc choices. We construct $H_t$ as a linear combination of $Z_t$ and an independent random process with the same distribution as $Z_t$. We set $\theta_i^{(w)}$ to be equal to a linear combination of $\hat\pi_i$ and an independent standard normal variable, and we do the same for $\theta_i^{(y)}$.
Note that this simulation design is not perfectly aligned with our theoretical results in Section (ref). We allow for additional randomness in errors $\{\epsilon_{it}^{(y)},\epsilon_{it}^{(w)}\}_{i\le n,t\le T}$. We make this choice to make the simulation less driven by particular features of the dataset, essentially allowing for matrices that are in the “neighborhood” of the original ones. We report the simulation results without these errors in Appendix (ref), and they largely agree with those in Table (ref).
We compare the performance of our estimator (as described by Algorithm (ref)) with the standard TSLS algorithm from Section (ref). In both cases, we use the data to construct $D_i$ by estimating the next equation by OLS, using data for $t\le \frac{T}{3}$:
and set $D_i := \hat \pi_i$. We consider four different designs. In the first design we drop $L_{it}^{(w)},L_{it}^{(y)}$, and $H_t$ from the model ((ref)). In this case, the TSLS algorithm should perform better than ours because it uses the optimal weights. With the second design, we start to increase the complexity, adding $L_{it}^{(w)}$ and $L_{it}^{(y)}$ back to the model. One can think of this design as a DGP for the data from nakamura2014fiscal under which the TSLS approach is justified. Here, we should expect both algorithms to perform well in terms of bias but potentially differ in variance. In the third design, we drop $L_{it}^{(w)},L_{it}^{(y)}$ but add $H_t$. Finally, in the fourth case, we have $L_{it}^{(w)}$ and $L_{it}^{(y)}$, and $H_t$.
In Table (ref), we report results over replication $1000$ for simulations for the case of $\tau = 1.43$, corresponding to the original point estimate obtained by nakamura2014fiscal. The results confirm the above intuition: in the simplest case, our estimator is less precise than $\hat \tau_{TSLS}$, although the difference is small. We see sizable gains in RMSE for the second design. In the third case, our estimator eliminates most of the bias, while the TSLS error is dominated by it. Finally, in the most general design, our estimator is nearly unbiased and dominates the TSLS in terms of RMSE. In Figure (ref), we plot the densities of $\hat \tau - \tau$ over the simulations for the second and fourth designs. These plots demonstrate the gains in variance and bias and show the estimator's overall behavior. Once again, we see that even when TSLS is approximately unbiased, there are gains from using our approach that come from increased precision.
The comparison in Table (ref) emphasizes the benefits of our approach, with an efficiency loss being present only in the most unrealistic first design. Importantly, these results do not reflect the efficiency cost of not using only part of the data, which we discussed in Section (ref). This cost is present for both estimators in our simulation because we use the first third of the data to construct $\{D_i\}_{i\le n}$. Formally, this situation is not covered by our theoretical results, which treat $\{D_i\}_{i \le n}$ as fixed quantities. However, the results in Table (ref) show that the algorithm continues to perform well in this regime, which is relevant for some empirical applications.
We also investigate the performance of our inference approach as described in Algorithm (ref). In Table (ref), we report coverage rates for nominal $95\%$ confidence intervals for $\hat\tau_{rob}$ and $\hat \tau_{TSLS}$. We construct $\hat\Lambda^{(z)}_{T_0+1|T}$ by fitting an ARIMA model to the data $\{Z_t\}_{t\le T}$ using the automatic model selection package in R. We see that the coverage is below nominal for all designs and estimators. This is not surprising, given that the sample size is relatively small, and in the third and fourth designs, both estimators are biased. In relative terms, the coverage for $\hat \tau_{rob}$ is closer to the nominal one.
Aggregate instruments provide a natural source of exogenous variation for unit-level outcomes. As a result, they are frequently used to evaluate the effects of local policies. This exercise has two conceptual steps: aggregating unit-level data into a time series and analyzing the aggregated data. We propose a new algorithm for constructing unit weights to produce aggregate outcomes. We use a flexible statistical model to show that our weights eliminate potential unobserved aggregate confounders, leading to a consistent and asymptotically normal estimator. Using data-driven simulations, we demonstrate the superiority of our proposal over the conventional TSLS estimator in various relevant regimes.