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.
72,794 characters · 17 sections · 76 citation commands
Forward-Selected Panel Data Approach for Program Evaluation
\affil[a]{Department of Economics, the Chinese University of Hong Kong, Sha Tin, New Territories, Hong Kong SAR, China} \affil[b]{Barclays Capital Asia Ltd., Cheung Kong Center, 2 Queen\textquoteright s Road Central, Hong Kong, Hong Kong SAR, China}
Key words: aggressive algorithm, average treatment effect, counterfactual analysis, post-selection inference
JEL code: C13, C21, C23, C38, D73
A controlled experiment compares outcomes of a treatment group with those from a control group. It is the golden standard for scientific research. While the randomized controlled trials are useful in understanding economic mechanisms duflo2007using,banerjee2009experimental, for large-scale questions economists mostly have access to observational datasets only. For example, in economic research we rarely enjoy the luxury to implement a controlled experiment at a national level that would affect millions of people\textemdash such an exercise can be prohibitively expensive or ethically unacceptable. Instead, economists resort to constructing counterfactuals from observational data for policy evaluation. A counterfactual is the potential outcome that can be perceived but has never happened in the real world.
In view of the lack of genuine control groups in many important economic empirical questions, hsiao2012panel (HCW, henceforth) propose the panel data approach (PDA) to exploit the correlation between cross-sectional units in estimating the counterfactual. Consider, for simplicity, that a single treatment (treatment and intervention are used as synonyms) is carried out during the observed time window. PDA is a linear regression on the cross-sectional units in the pre-treatment data, and then these estimated coefficients are used to extrapolate the counterfactual of no policy intervention in the post-treatment period. Its convenience attracts many applications and extensions, for example bai2014property, fujiki2015disentangling, ouyang2015treatment, ke2017china, to name a few. Compared with the popular difference-in-difference method, the combination of control units allows time-varying treatment effect. Alternatively, abadie2003economic and abadie2010synthetic advocate the synthetic control method (SCM). gardeazabal2017empirical and wan2018panel compare PDA and SCM in simulations and empirical applications.
Choices of the control units directly affect PDA's estimation and inference results, and thus a systematic variable selection scheme is of vital practical importance. HCW experiment with the Akaike information criterion (AIC) and the corrected AIC (AICC), and du2015home recommend the latter for consistent variable selection. These conventional variable selection methods compute an information criterion for each candidate model to identify the “best subset”. In PDA the total number of candidate models is $2^{N}$, where $N$ is the number of available potential control units. In spite of the state-of-the-art computing technology, exhaustive search quickly becomes too time-consuming for a moderate $N$. The exhaustive enumeration is inapplicable in the era of big data when data-rich environments offer information at an unprecedented scale.
Furthermore, besides the computational difficulty, a large cross-sectional dimension also challenges PDA's theoretical justification. As PDA is often applied to aggregate data with low-frequency temporal observations, HCW's conventional “fixed $N$, large $T$”\footnote{Here $T$ represents the size of the time dimension in a generic panel data setting. We will elaborate it when introducing the model.} asymptotic framework is unlikely to deliver satisfactory approximation in empirical studies where $N$ is comparable to $T$, or even exceeds $T$. To overcome the high dimensionality, li2016estimation suggest using Lasso tibshirani1996regression but provide no theoretical foundation. carvalho2016arco develop the Lasso theory under a general framework called Artificial Counterfactual (ArCo).
This paper studies estimation and inference of the average treatment effect (ATE) by PDA when a large number of candidate cross-sectional control units are present. Motivated by real data applications, we tackle variable selection in high dimension.
This paper contributes, toward the growing toolkit of program evaluation, a user-friendly procedure with asymptotic guarantee to make it possible to incorporate all potential control units with accessible time series. It alleviates the arbitrariness of control units selection --- an essential step in program evaluation.
In terms of the algorithms, we suggest using forward selection to choose a sequence of control units one by one up to a desirable number. Involving only a series of OLS regressions, forward selection is computationally efficient. For hypothesis testing about the ATE, we advocate calculating a conventional $t$-statistic conditioning on the control variables chosen by forward selection and then comparing it with a quantile of the standard normal distribution. We call this two-step method the forward-selected PDA (fsPDA), as in the title of the paper.
Despite the simplicity of the algorithms, the underlying asymptotic theory for fsPDA nevertheless demands careful development and justification. The environment of independently and identically distributed (i.i.d.) data in which most high-dimensional statistical problems are investigated is too restrictive for economic questions with temporal observations. Accommodating heterogeneous weakly dependent time series, we establish our theory in the asymptotic framework allowing $N/T\to\infty$ when both dimensions diverge.
A unique innovation is that our theory is valid in regression models no matter whether the “true” underlying coefficients are dense or sparse. It differentiates us from the vast literature of high-dimensional statistics that counts on sparse regression models for asymptotic results. Here dense models impose no restrictions on the coefficients \textemdash in principle all coefficients can be simultaneously non-zero.\footnote{Dense model are of rising interest in economics where observed variables are interconnected and little prior knowledge ensures that the driven forces fall into a few of key variables giannone2017economic.} In contrast, sparsity means most of the regression coefficients are exactly zero or too small to matter. As to be discussed in Section (ref), PDA is motivated from a factor model which induces a dense regression in general. We make it clear that the inference of ATE in the post-treatment data for correct test size does not require consistent estimation of the true underlying high-dimensional coefficients in the data generating process (DGP). Instead, it is sufficient if we can recover the linear projection coefficients associated with a small subset of control variables.
We add to the literature of both variable selection and statistical inference in the dense model setting. (i) We show that forward selection is capable of reducing the variance of the regression error as much as that of a computationally infeasible best subset. Although forward selection is by no means a new algorithm, in the past its theory is established either in sparse statistical models Zhu2017forward or dense population models das2011submodular. We are unaware of its theory in dense models with sampling error. (ii) Many asymptotic normality results in high-dimensional models are pointwise asymptotics under a single DGP, for example the so-called oracle property fan2001variable,zou2006adaptive. Our inferential theory is uniform under DGPs that satisfy a set of conditions. In other words, the seemingly naive practice of conventional normal inference is valid and the randomness stemming from the step of variable selection can be safely ignored. This validity comes from the special structure of PDA, in which the pre-treatment and post-treatment periods are naturally separated into two disjoint segments. The control units are selected from the pre-treatment data only, and under weak dependence they become asymptotically independent of the post-treatment data on which the test statistic is based.
This paper fits in the theme of research to better connect settings of economic interest with modern high-dimensional statistics. It transpires that the economic context of PDA underpins our unsophisticated procedure and circumvents statistical challenges encountered by post-selection uniform inference in a single dataset. The PDA environment also helps with other technical building blocks. For instance, the restricted eigenvalue condition bickel2009simultaneous is often viewed as a necessary but somewhat ad hoc assumption in high-dimensional regressions. Rather than imposing it, we argue that a version of the restricted eigenvalue condition is a natural implication of the underlying latent factor model that motivates PDA. The theoretical properties that we established for PDA are supported by extensive Monte Carlo simulations.
Literature Review. Various greedy variable selection algorithms have been studied in operational research, statistics and econometrics. Working with random samples, wang2009forward, Zhu2017forward, and Zhu2019JASA analyze forward selection as a device for model determination in statistical ultrahigh-dimensional sparse regressions. kozbur2017testing, kozbur2018sharp and hansen2018targeted investigate test-based stopping criteria and post-selection inference. Our paper extends the operational research by das2011submodular and das2018approximate who highlight the key role played by submodular ratio in the population model analysis of forward selection. When we carry forward selection over into panel data, we must cope with sampling uncertainty as well as temporal dependence. The greedy nature of forward selection is closely related to the component-wise boosting buhlmann2006,luo2016high which is familiar to econometricians bai2009boosting,shi2015REL,fonseca2018boost,Phillips2019boosting. Alternatively, carvalho2016arco studies asymptotic validity of ATE estimation by Lasso-type methods in sparse models.
Uniform inference after variable selection is a difficult statistical issue, as pointed out by leeb2005model and leeb2006can. Proposed solutions usually resort to non-standard methods, for example berk2013valid, fithian2014optimal and tibshirani2018uniform, when model selection and testing are carried out within a single dataset. Predictive inference, however, provides an amenable environment to work with and leeb2009conditional shows post-selection asymptotic normality. Another related line of uniform inference literature tries to correct the shrinkage bias in high-dimensional regressions belloni2014uniform,belloni2017program,javanmard2018debiasing. In our paper we always use OLS to estimate coefficients so we are free from shrinkage biases caused by penalized estimation.
Notations. Unless explicitly defined otherwise, we use a plain letter, say “$x$”, to denote a scalar, a boldface lowercase letter “$\mathbf{x}$” to denote a column vector, and a boldface uppercase letter “$\mathbf{X}$” to denote a matrix. The square matrix $\mathbf{I}$ is an identity matrix. 1$\left\{ \cdot\right\} $ is the indicator function. For a real number, $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution $N\left(0,1\right)$, $\left\lceil \cdot\right\rceil $ is the ceiling function and $\left\lfloor \cdot\right\rfloor $ is the floor function. For a square matrix, $\left(\cdot\right)^{-}$ is the Moore-Penrose generalized inverse, and $\phi_{\min}\left(\cdot\right)$ and $\phi_{\max}\left(\cdot\right)$ are the minimum eigenvalue and the maximum eigenvalue, respectively. The cardinality of a discrete set $U$ is denoted as $\left|U\right|$. A vector with a discrete set as its subscript, $\mathbf{x}_{U}:=(x_{j})_{j\in U}$, makes a $\left|U\right|$-element subvector of $\mathbf{x}$. $\left\Vert \cdot\right\Vert _{2}$ and $\left\Vert \cdot\right\Vert _{1}$ are the usual $L_{2}$ and $L_{1}$ vector norms, respectively.
Now we introduce PDA's panel data setting. $(N+1)$ cross-sectional units in a panel data are indexed by $\mathcal{N}_{0}:=\left\{ 0,1,\ldots,N\right\} $, in which $j=0$ indexes the sole treated unit whereas $\mathcal{N}:=\left\{ 1,\ldots,N\right\} $ is the index set of the $N$ control units. In the potential outcome framework, let $y_{jt}^{1}$ and $y_{jt}^{0}$ be the outcomes of the unit $j$ at time $t$ with and without a policy intervention, respectively. We cannot witness $y_{jt}^{1}$ and $y_{jt}^{0}$ simultaneously; instead we observe $y_{jt}=y_{jt}^{0}(1-d_{jt})+y_{jt}^{1}d_{jt}$, where $d_{jt}$ is a dummy variable equal to 1 if the $j$-th unit is under intervention at time $t$; otherwise $d_{jt}=0$.
The time dimension of the panel data is shown in the timeline diagram in Figure (ref). The time series extends from $-\infty$ to $\infty$, while we econometricians observe a time window $\left\{ -T_{1},\ldots-1,0,1,\ldots T_{2}\right\} $ for some $T_{1},T_{2}\in\mathbb{N}$. Without loss of generality, a policy intervention occurs at time $t=0$, which partitions the observed interval into two sections: a pre-treatment period $\mathcal{T}_{1}:=\left\{ -T_{1},\ldots,-1\right\} $ and a post-treatment period $\mathcal{T}_{2}:=\left\{ 1,\ldots,T_{2}\right\} $, with lengths $T_{1}=\left|\mathcal{T}_{1}\right|$ and $T_{2}=\left|\mathcal{T}_{2}\right|$. Denote $\mathcal{T}:=\mathcal{T}_{1}\cup\mathcal{T}_{2}$, and $T:=\left|\mathcal{T}\right|$.
The mathematical expectation of a generic random variable $x_{t}$ is denoted as $E\left[x_{t}\right]$. For a heterogeneous time series, we define $\mathcal{E}_{\left(1\right)}\left[x_{t}\right]:=T_{1}^{-1}\sum_{t\in\mathcal{T}_{1}}E\left[x_{t}\right]$ as the average of the expectations $(E[x_{t}])_{t\in\mathcal{T}_{1}}$ in the pre-treatment period, and similarly $\mathcal{E}_{\left(2\right)}\left[x_{t}\right]:=T_{2}^{-1}\sum_{t\in\mathcal{T}_{2}}E\left[x_{t}\right]$ as the average of the expectations after the treatment. We define $\mathbb{E}_{\left(1\right)}\left[x_{t}\right]:=T_{1}^{-1}\sum_{t\in\mathcal{T}_{1}}x_{t}$ as the pre-treatment sample mean, and $\mathbb{E}_{\left(2\right)}\left[x_{t}\right]=T_{2}^{-1}\sum_{t\in\mathcal{T}_{2}}x_{t}$ as the post-treatment sample mean.
Plan. The rest of the paper is organized as follows. Section (ref) introduces PDA and describes forward selection and the post-selection ATE inference. Section (ref) presents asymptotic analysis of fsPDA. Section (ref) reports the simulation results, and Section (ref) carries out an empirical example. All proofs are provided in the appendix. Moreover, an online supplement is prepared with further comparison of methods, one more empirical example, and additional simulation results. Replication data and a documented R package {\tt fsPDA} are hosted at \url{https://github.com/zhentaoshi/fsPDA}.
PDA is motivated from a factor model. For the completeness of the paper, we briefly summarize HCW's proposal. Consider a pure factor model in which all cross-sectional units share at most $K$ common factors:
where $\mathbf{f}_{t}$ is a mean zero $K$-vector of latent factors, $\boldsymbol{\lambda}_{j}$ is a $K$-vector of factor loading, and $e_{jt}$ is a mean zero idiosyncratic error component orthogonal to the factors.\footnote{For simplicity of presentation, we assume $E[y_{jt}^{0}]=0$ for all $j\in\mathcal{N}_{0},\ t\in\mathcal{T}$ and the linear regressions in Section (ref) does not include an intercept. While the presence of the intercept incurs extra notation, incorporating it does not affect the asymptotic theory buhlmann2011statistics. The intercept will be included in the empirical applications.}
HCW assume only one unit ($j=0$) is exposed to a policy intervention, so the intervention does not affect the outcomes of all other units $i\in\mathcal{N}$. After the treatment, the observed outcome for the treated unit is
where $\Delta_{t}$ is the treatment effect at time $t$. PDA is interested in the null hypothesis of, without loss of generality, zero ATE \[ \mathbb{H}_{0}:\ \mathcal{E}_{\left(2\right)}\left[\Delta_{t}\right]=0. \] For only $y_{0t}^{1}$ is observable after the intervention, to evaluate the treatment effect we must estimate the counterfactual $y_{0t}^{0}$ for $t\in\mathcal{T}_{2}$ from the observed data. li2016estimation show that, based on the factor model, there exists an $N$-vector $\boldsymbol{\boldsymbol{\beta}}_{\mathcal{N}}^{0}$ such that the outcome of the treated unit can be written as a linear combination of the outcomes of the control units plus an orthogonal error
where $\mathbf{y}_{\mathcal{N}t}:=(y_{jt})_{j\in\mathcal{N}}$ is an $N$-dimensional vector. Although the regression equation ((ref)) \textemdash PDA's workhorse for estimation and inference \textemdash is derived from the factor model ((ref)), the factor model itself merely serves as a motivation but is irrelevant to PDA's implementation.
With the pre-treatment sub-sample $\mathcal{T}_{1}$, HCW advocate estimating $\widehat{\boldsymbol{\beta}}_{\mathcal{N}}$ by ordinary least squares (OLS) or generalized least squares (GLS), and predicting the counterfactual as $\widehat{y}_{0t}^{0}=\mathbf{y}_{\mathcal{N}t}'\widehat{\boldsymbol{\beta}}_{\mathcal{N}}$ for $t\in\mathcal{T}_{2}$ and the treatment effect $\widehat{\Delta}_{t}=y_{0t}^{1}-\widehat{y}_{0t}^{0}$. The inference is routine. They construct the usual $t$-statistic based on the sample mean of $\{\widehat{\Delta}_{t}\}_{t\in\mathcal{T}_{2}}$ and its long-run variance, and then compare the absolute value of the $t$-statistic with some quantile of $N\left(0,1\right)$, say 1.96 for two-sided test of size 5%, to decide whether the null should be rejected.
The estimate of PDA depends on the choice of the control units. When the number of potential controls is large, HCW's information criterion approach will encounter computational difficulty in exhaustive search. To solve this problem, we propose using the forward selection method hastie2009bible. In the first iteration, we regress $\left(y_{0t}\right)_{t\in\mathcal{T}_{1}}$ on $\left(y_{jt}\right)_{t\in\mathcal{T}_{1}}$ for each $j\in\mathcal{N}$, and choose the one that maximizes the R-squared $\mathscr{R}^{2}\left(\left\{ j\right\} \right)$ of OLS, where $\left\{ j\right\} $ in the parenthesis stresses the set of control units on which the R-squared $\mathscr{R}^{2}$ is based. We denote the index of the maximizer as $\widehat{j}_{1}$ and let $\widehat{U}_{1}=\{\widehat{j}_{1}\}$ be a single-element set. In the $r$-th iteration, where $r=2,..,R$, we run OLS of $\left(y_{0t}\right)_{t\in\mathcal{T}_{1}}$ on $(\mathbf{y}_{\widehat{U}_{r-1}t})_{t\in\mathcal{T}_{1}}$ together with another single $\left(y_{jt}\right)_{t\in\mathcal{T}_{1}}$ for each $j\in\mathcal{N}\backslash\widehat{U}_{r-1}$, choose the one \textemdash denoted as $\widehat{j}_{r}$ \textemdash that maximizes the corresponding R-squared $\mathscr{R}^{2}(\widehat{U}_{r-1}\cup\{j\})$, and incorporate it into the selected set $\widehat{U}_{r}=\widehat{U}_{r-1}\cup\{\widehat{j}_{r}\}$. The total number of iterations, $R$, is a tuning parameter specified by the user. The algorithm is described formally as follows.
The above forward selection procedure is a greedy algorithm that takes the most aggressive direction in each step to increase the R-squared, or equivalently to reduce the sum of squared residuals, conditional on the variables that are already included. Once a variable is selected, there is no mechanism to drop it. Greedy algorithms are popular in modern machine learning. For example, breiman2001random grows regression trees by splitting a single variable each time at the deepest descent, and buhlmann2006's componentwise boosting seeks the most greedy variable without adjusting other coefficients.
The ultimate goal of PDA is statistical inference for the ATE. If we had prior knowledge about an index set $U\subset\mathcal{N}$ of relevant control units, we would naturally carry out the following procedure. We would regress $\left(y_{0t}\right)_{t\in\mathcal{T}_{1}}$ on $(\mathbf{y}_{Ut}=(y_{jt})_{j\in U})_{t\in\mathcal{T}_{1}}$ to obtain the coefficient $\widehat{\boldsymbol{\beta}}_{U}$ and predict the counterfactual $\widehat{y}_{0t}^{0}\left(U\right):=\mathbf{y}_{Ut}'\widehat{\boldsymbol{\beta}}_{U}$. Next, we would estimate the treatment effect based on the set $U$ as \[ \widehat{\Delta}_{Ut}:=y_{0t}^{1}-\widehat{y}_{0t}^{0}\left(U\right)=y_{0t}^{1}-\mathbf{y}_{Ut}'\widehat{\boldsymbol{\beta}}_{U},\ t\in\mathcal{T}_{2}. \] Let $\widehat{\rho}_{\tau U}^{2}:=T_{2}^{-1}\sum_{t,s\in\mathcal{T}_{2}}(\widehat{\Delta}_{Ut}-\bar{\Delta}_{U})(\widehat{\Delta}_{Us}-\bar{\Delta}_{U})\cdot\boldsymbol{1}\left\{ \left|t-s\right|\leq\tau\right\} $ be a heteroskedasticity- and autocorrelation-consistent (HAC) estimator of the long-run variance, where $\tau$ is the number of lags and $\bar{\Delta}_{U}:=\mathbb{E}_{\left(2\right)}[\widehat{\Delta}_{Ut}]$. We calculate the $t$-statistic
which depends on $\tau$ and $T_{2}$ while we suppress them for conciseness. We would reject the null hypothesis at size $a$ when $|\mathcal{Z}_{U}|>\Phi^{-1}\left(1-a/2\right)$, provided that the distribution of $\mathcal{Z}_{U}$ can be approximated by $N\left(0,1\right)$.
In reality we rarely know in advance a set of relevant control units $U$. We suggest using $\widehat{U}_{R}$, the set chosen by forward selection, to substitute the generic $U$ in the $t$-statistic in the above paragraph. That is, we reject the null hypothesis $\mathbb{H}_{0}$ at 5% size if $|\mathcal{Z}_{\widehat{U}_{R}}|>1.96$. Although $\widehat{U}_{R}$ is a random set determined by the pre-treatment data, we will show $N\left(0,1\right)$ is a reasonable approximation to the statistic $\mathcal{Z}_{\widehat{U}_{R}}$ under the null along with mild assumptions of weak temporal dependence.
There are two tuning parameters in fsPDA, $R$ for the total number of selected variables and $\tau$ for the long-run variance estimation. We suggest using wang2009shrinkage's modified Bayesian information criterion (modified BIC) to choose $R$, while the choice of $\tau$ has been well studied in econometrics literature newey1987simple,andrews1991heteroskedasticity.
Before we conclude this section, we emphasize that we do not attempt to directly estimate the factor model due to the following reasons. (i) In the PDA framework the factor model is an abstraction independent of the algorithm based on linear regressions, and this regression approach is also followed by li2016estimation and carvalho2016arco. (ii) To conduct inference literally in the factor model, we will need to estimate the $\left(N+1\right)\times\left(N+1\right)$ covariance matrix for the idiosyncratic noises $\left(e_{jt}\right)_{j\in\mathcal{N}_{0}}$, which involves $\left(N+2\right)\left(N+1\right)/2$ unknown entries so other sparse matrix estimation techniques have to be invoked for dimension reduction.
Section (ref) has introduced forward selection and then the $t$-statistic based on the selected variables. We proceed by establishing asymptotic guarantee for this procedure. After laying out the regularity conditions, we reverse the order by first studying a generic post-selection inference in this context of ATE estimation, and then arguing that the forward selection is a competitive method for variable selection.
We work with a multi-index asymptotic framework. In asymptotic statements, we take the number of cross-sectional units $N\to\infty$, while the number of pre-treatment observations $T_{1}=T_{1}\left(N\right)$ is understood as a deterministic function of $N$ such that $T_{1}\to\infty$ as $N\to\infty$. $N$ is allowed to be larger than $T_{1}$ to accommodate high-dimensional settings, although $\left(\log N\right)/T_{1}\to0$. Similar indexing is applied to the number of the post-treatment observations $T_{2}=T_{2}\left(N\right)\to\infty$ and $\left(\log N\right)/T_{2}\to0$ as $N\to\infty$.
The algorithm of forward selection uses the pre-treatment data only. To study the asymptotic properties of forward selection and post-selection inference, we impose two high-level assumptions. The first one regularizes the minimal eigenvalue of the population Gram matrix. Let $\eta_{r}=\min_{\mathcal{U}_{r}}\phi_{\min}\left(\mathcal{E}_{\left(1\right)}\left[\mathbf{y}_{Ut}\mathbf{y}_{Ut}^{\prime}\right]\right)$ where $\mathcal{U}_{r}:=\left\{ U\subset\mathcal{N}:\left|U\right|\leq\left\lfloor r\right\rfloor \right\} $ for some $r\in\mathbb{\mathbb{R}}^{+}$. A universal constant is a strictly positive finite real number that is independent of sample sizes.
As $R$ diverges to infinity at a rate slower than $\left(T_{1}/\log N\right)^{1/3}$, the sample version of the $u\times u$ Gram submatrix $\mathbb{E}_{\left(1\right)}\left[\mathbf{y}_{Ut}\mathbf{y}_{Ut}^{\prime}\right]$ involving $T_{1}$ time series observations is likely to be of full rank when $u\ll T_{1}$, with the help of the second assumption below about the population second-moment as well as their sample counterpart.
Assumption (ref)(a) postulates a convergence rate of the second moments, and (b) is a common assumption of finite population second moments. With independent observations, belloni2012sparse use the self-normalized Cram\'{e}r-type moderate-deviation theory jing2003self to establish the probabilistic bound in (a). In time series contexts, similar conditions are used in medeiros2016l1, kock2015oracle, and koo2016high under various assumptions of tail bounds and serial dependence.
In the population model, $\boldsymbol{\beta}_{U}^{0}:=(\mathcal{E}_{\left(1\right)}[\mathbf{y}_{Ut}\mathbf{y}_{Ut}^{\prime}])^{-}\mathcal{E}_{\left(1\right)}\left[\mathbf{y}_{Ut}y_{0t}\right]$ is the “true” linear projection coefficient under a given $U$, and the corresponding projection error is $\varepsilon_{Ut}:=y_{0t}-\mathbf{y}_{Ut}^{\prime}\boldsymbol{\beta}_{U}^{0}$. Let $\sigma_{U}^{2}:=\mathcal{E}_{\left(1\right)}[\varepsilon_{Ut}^{2}]$ be the population variance of the projection error under the set $U$, and $\widehat{\sigma}_{U}^{2}$ be the sample variance of $(\widehat{\varepsilon}_{Ut}:=y_{0t}-\mathbf{y}_{Ut}^{\prime}\widehat{\boldsymbol{\beta}}_{U})_{t\in\mathcal{T}_{1}}$. The following lemma shows the sample variance $\widehat{\sigma}_{U}^{2}$ approximates the population counterpart $\sigma_{U}^{2}$, and the OLS estimator $\widehat{\boldsymbol{\beta}}_{U}=(\mathbb{E}_{\left(1\right)}[\mathbf{y}_{Ut}\mathbf{y}_{Ut}^{\prime}])^{-}\mathbb{E}_{\left(1\right)}\left[\mathbf{y}_{Ut}y_{0t}\right]$ approximates $\boldsymbol{\beta}_{U}^{0}$.
Lemma (ref) indicates that over all index sets $U$'s with no more than $\left\lfloor \left(1+\delta_{1}\right)R\right\rfloor $ elements, if $R$ diverges slowly such that $1/R+R/\left(T_{1}/\log N\right)^{1/3}\to0$ as in Assumption (ref), then the difference between $\widehat{\sigma}_{U}^{2}$ and $\sigma_{U}^{2}$ is negligible in probability. Similar approximation holds in the coefficient estimation for $\boldsymbol{\beta}_{U}^{0}$. These are results prepared for the following two subsections.
We first work on the asymptotic property of the post-selection $t$-statistic based on a generic data-driven variable selection method using the pre-treatment data only. It will include the forward selection as a special case.
In inference we must use the post-treatment data, on which we impose a few regularity assumptions.
In the post-treatment subsample, Assumption (ref)(a) is about the convergence rate of the sample mean to the population mean 0, although $y_{0t}^{0}$ is unobservable. (b) is analogous to Assumption (ref)(a) in the pre-treatment period, and the fourth moment in (c) is commonly imposed in studies of inferential procedures for high-dimensional factor models bai2003inferential. The last two items in Assumption (ref) are concerning the long-run variance, where (d) bounds the long-run variance from degeneracy and (e) guarantees the absolute summability of the autocorrelations. (c), (d) and (e) make sure that the self-normalized test statistic behaves well, so that a suitable version of the Berry-Essen bound can be applied to establish the asymptotic normality of the test statistic.
Next, we introduce the time series weak dependence structure. Let $\mathcal{F}_{N}^{t_{1},t_{2}}$ be the smallest $\sigma$-field generated by the Borel sets of the collection $\{\left(\mathbf{f}_{t}',\mathbf{e}_{\mathcal{N}_{0}t}'\right)':t_{1}\leq t\leq t_{2}\}$ from the factor model ((ref)), where it naturally incorporates the fact that no random variables are produced at $t=0$, the calendar date for the treatment. In view of the infinite time series in Figure (ref), for each $k\in\mathbb{N}$ we define
where $\mathbb{Z}$ is the set of all integers. The dependence indicator $\phi_{N}\left(k\right)$ is the uniform strong mixing coefficient davidson1994stochastic. We impose the following Assumption (ref), which is similar to carvalho2016arco's Assumption 3 of geometric strong mixing.
The above assumption is employed for two technical purposes: (i) It allows us to invoke the Berry-Essen bound for heterogeneous time series bentkus1997berry,sunklodas2000approximation. (ii) It implies the asymptotic independence, as $k\to\infty$, between the events in $\mathcal{F}_{N}^{-\infty,-1}$ before the treatment and the events in $\mathcal{F}_{N}^{k,\infty}$ which is $k$ periods after the treatment. The second point is critical for asymptotic normality. If a single dataset is used for model selection and parameter estimation, post-selection inference is in general a very difficult statistical problem that leads to non-standard asymptotic distributions leeb2005model,leeb2006can, and it is a topic of intensive recent research berk2013valid,belloni2014uniform,belloni2017program,hansen2018targeted. However, in conditional (on the selected model from a training sample) predictive inference, post-selection asymptotic normality is achievable leeb2009conditional and the inference can be carried out following standard asymptotically normal procedure.
In our context, the estimated ATE is based on the average involving the predicted outcomes over the post-treatment period $\mathcal{T}_{2}$. Between the two blocks $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$, the observations near the treatment date $t=0$ are essentially dependent. For instance, those with time index $t=1,\ldots,k$ are statistically dependent on the random variables at the end of $\mathcal{T}_{1}$. This dependent episode consists of a smaller and smaller fraction of the post-treatment sample if we devise $k=k\left(N\right)$ such that $k/T_{2}\to\infty$ as $N\to\infty$.
Let $M$ be the DGP that generates $\mathbf{y}_{\mathcal{N}_{0}t}^{0}$ in ((ref)). Let $\check{U}_{R}\in\mathcal{U}_{R}$ be an index set estimated by an arbitrary variable selection method using the pre-treatment dataset only. Let $\mathcal{Z}_{\check{U}_{R}}:=\mathcal{Z}_{U}|_{U=\check{U}_{R}}$ for the $t$-statistic $\mathcal{Z}_{U}$ defined in ((ref)) evaluated at $U=\check{U}_{R}$. Obviously $\mathcal{Z}_{\check{U}_{R}}=\mathcal{Z}_{\check{U}_{R}}\left(M\right)$ depends on the underlying DGP $M$.
Consider a set of DGPs $\mathcal{M}$ such that Assumptions (ref), (ref), (ref) and (ref) hold uniformly. This uniformity requirement strengthens the stochastic orders in Assumptions (ref)(a), (ref)(a) and (b), whereas all other assumptions are already stated with universal constants. The following theorem provides the asymptotic distribution of $\mathcal{Z}_{\check{U}_{R}}$ uniformly over the set of eligible DGPs $M\in\mathcal{M}$.
The restriction $M\in\mathcal{M}$ implicitly imposes Assumptions (ref), (ref), (ref) and (ref) uniformly. Theorem (ref) is established by a Berry-Essen bound for heterogeneous time series sunklodas2000approximation. The key condition that contributes to the uniform asymptotic normality is the weak dependence in Assumption (ref) along with our setting of ATE estimation, in which the policy intervention occurs at time $t=0$ splits the sample into two disjoint subsamples indexed by $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$, respectively. Sampling splitting is a popular approach to achieve uniform inference in statistical machine learning in cross-sectional environments, for example belloni2014uniform and wager2018estimation. The notation of uniform strong mixing is a time series analogy of asymptotic independence.
The uniform asymptotic normality in Theorem (ref) holds regardless of the algorithm that selects a subset of no more than $R$ control variables. Consider an ad hoc non-random way of choosing a sequence of sets. Given an arbitrary ordering of the control units, we may naively choose the first $R$ terms $U_{R}^{\mathrm{naive}}=\left\{ 1,\ldots,R\right\} $ for $R$ satisfying the order regularized by the conditions in Assumption (ref), and we would have $\mathcal{Z}_{U_{R}^{\mathrm{naive}}}\stackrel{d}{\to}N\left(0,1\right)$. It is also applicable to the $t$-statistic based on HCW's best subset method via AIC or AICC. When they developed the asymptotic inference, HCW heuristically took the selected variables, which we denote here as $\check{U}_{R}^{\mathrm{AICC}}$, as if they were fixed. Our result implies $\mathcal{Z}_{\check{U}_{R}^{\mathrm{AICC}}}\stackrel{d}{\to}N\left(0,1\right)$, which helps justify HCW's practice.
Instead of $\check{U}_{R}^{\mathrm{AICC}}$ that is based on exhaustive search over all subsets, we nevertheless advocate the forward selection algorithm for $\widehat{U}_{R}$ in view of its convenience in computation in high-dimensional settings. The asymptotic theory of forward selection is developed in the next subsection.
In order to discuss the efficacy of forward selection, we spell out our target for variable selection. Let $U_{u}^{*}:=\arg\min_{\mathcal{U}_{u}}\sigma_{U}^{2}$ be the best subset of $u$ elements among all $U\subset\mathcal{N}$; the number of elements is no more than some $u\in\mathbb{N}$. Let $\sigma_{u}^{*2}:=\sigma_{U_{u}^{*}}^{2}=\sigma_{U}^{2}|_{U=U_{u}^{*}}$ be the corresponding noise level under this best subset $U_{u}^{*}$. If $U_{u}^{*}$ is not unique, we simply refer to any of them as the best subset and our analysis is not affected no matter $U_{u}^{*}$ is unique or not. It is computationally expensive to locate the best subset $U_{u}^{*}$. Even if $\sigma_{U}^{*2}$ were estimated with no noise, we would exhaustively compare $\sigma_{U}^{*2}$ for ${N \choose u}$ models, which is of exponential order of $N$.
Instead of searching for $U_{u}^{*}$, we seek to identify a subset $U$ on which $\sigma_{U}^{2}$ approximates the optimal variance $\sigma_{u}^{*2}$. Theorem (ref) below states that the greedy forward selection algorithm picks up a set $\widehat{U}_{R}$ with a regression variance asymptotically as small as the desired $u$-element best subset if $R$ dominates $u$ asymptotically. The greedy algorithm only searches among $\sum_{r=1}^{R}\left(N-r+1\right)$ models, which is of linear order of $N$. The latter is computationally much more efficient than exhaustive search.
Since variable selection does not use the post-treatment subsample, only Assumptions (ref) and (ref) are needed for Theorem (ref). The above theorem is a nearly optimal result. It implies with high probability that the computationally feasible sample variance $\widehat{\sigma}_{\widehat{U}_{R}}^{2}$ is asymptotically no worse, up to an arbitrarily small tolerance $\delta_{2}$, than the computationally heavy but theoretically optimal $\sigma_{u}^{*2}$. Such approximation can be achieved by incorporating $R$ units. Though $R$ is of bigger order than $u$ in the asymptotic sense, if we specify $R=\left\lfloor u\log\log N\right\rfloor $, then obviously the number of OLS regressions is fewer than $Nu\log\log N$, and $Nu\log\log N\ll{N \choose u}$ for a non-trivial $u$ and large $N$.
Theorem (ref) holds uniformly if the DGPs under consideration are restricted to $\mathcal{M}$. After forward selection, we use $\mathbf{y}_{\widehat{U}_{R}t}$ to predict the counterfactual $y_{0t}^{0}$ and obtain the time-varying treatment effect $\widehat{\Delta}_{\widehat{U}_{R}t}$ in the post-treatment period. Since $\mathcal{Z}_{\widehat{U}_{R}}\left(M\right)$ is a special case of $\mathcal{Z}_{\check{U}_{R}}\left(M\right)$ in Theorem (ref) when we use forward selection to choose variables, the following corollary is an immediate implication.
We summarize the theoretical results in this paper. Theorem (ref) shows that the $t$-statistic based on a generic variable selection method from the pre-treatment data has correct size. Theorem (ref) highlights that the forward selection algorithm can attain variance for the regression model as small as that of the best subset $U_{u}^{*}$ if $R$ dominates the cardinality $u$ asymptotically. The small $\widehat{\sigma}_{\widehat{U}_{R}}^{2}$ in general improves the statistical efficiency of the hypothesis testing.
We evaluate the finite-sample performance of our proposed procedure by Monte Carlo simulations. We conduct extensive experiments with non-sparse coefficients and sparse ones, and with various degrees of cross-sectional correlation and time dependence.\footnote{Due to the limitations of space, in the main text we present results for a dense underlying linear regression model. In the online supplement, we document the performance of variable selection, parameter estimation and prediction accuracy for a sparse model.} For comparison, we also estimate the model using Lasso. For each DGP, we generate one treated unit $j=0$ along with 100 control units $j=1,\ldots,100$. We run $1000$ replications and check the out-of-sample root mean predicted squared error (RMPSE) as well as the test size or power for $\mathbb{H}_{0}$. For simplicity, we set equal the lengths of the pre-treatment and post-treatment time series, with $T_{1}=T_{2}=50,100$ or 200.
Both forward selection and Lasso need turning parameters: the stopping time $R$ in forward selection and the penalty level $\lambda$ in Lasso. We adopt the modified BIC wang2009shrinkage in choosing the tuning parameters. For forward selection, the stopping time $R$ is determined by \[ \widehat{R}=\arg\min_{r\in\mathbb{N}}\left\{ \log\left(\widehat{\sigma}_{\widehat{U}_{r}}^{2}\right)+\log\log N\cdot r(\log T_{1})/T_{1}\right\} . \] Lasso's tuning parameter is determined by \[ \widehat{\lambda}=\arg\min_{\lambda}\left\{ \log\left(\mathbb{E}_{\left(1\right)}\left[(y_{0t}-\mathbf{y}_{\mathcal{N}t}'\widehat{\boldsymbol{\beta}}_{\lambda}^{las})^{2}\right]\right)+2\log\log N\cdot\Vert\widehat{\boldsymbol{\beta}}_{\lambda}^{las}\Vert_{0}(\log T_{1})/T_{1}\right\} , \] where $\Vert\widehat{\boldsymbol{\beta}}_{\lambda}^{las}\Vert_{0}$ is the number of non-zero coordinates in the vector $\widehat{\boldsymbol{\beta}}_{\lambda}^{las}$. In the second term of the modified BIC, we have the admittedly ad hoc constant 1 for forward selection and 2 for the Lasso, respectively. The difference arises because in our simulations Lasso would select many more variables than forward selection were the same constant shared in the two estimation methods, resulting in unsatisfactory performance. The choice of the constant will be commented in the continued Remark (ref) in this section.
We generate the data via the factor model ((ref)) with 4 common factors. The idiosyncratic shocks $e_{jt}\sim N\left(0,0.5^{2}\right)$ is independent across $j$ and $t$.
The factor loading $\lambda_{jl}$, $l=1,\ldots,4$, is independently drawn from $\mathrm{Uniform}\left(1,2\right)$ if $j=0,\ldots,4$, whereas $\lambda_{jl}\sim\mathrm{Uniform}\left(-0.1,0.1\right)$ if $j=5,\ldots,100$.
For $t\in\mathcal{T}_{2}$, the treated unit $y_{0t}$ is subject to an exogenous shock $\Delta_{t}$. We generate $\Delta_{t}$ by seven DGPs, denoted as $D1$ to $D7$:
The null hypothesis is true under $D1$\textendash $D3$, but false under $D4$\textendash $D7$. The treatment is time-invariant under $D1$, time-varying under $D2$, and serially correlated under $D3$. Mean shifts are introduced to post-treatment outcomes in $D4$ and $D5$, whereas $D6$ and $D7$ add non-zero dynamic treatment effects.
The first two columns of Table (ref) report the number of non-zero coefficients and the empirical RMPSE $\big(\mathbb{E}_{\left(2\right)}[\left(y_{0t}^{0}-\widehat{y}_{0t}\right)^{2}]\big)^{1/2}$, where $\widehat{y}_{0t}$ is the predicted value for $y_{0t}$: forward selection gives $\widehat{y}_{0t}=\widehat{y}_{0t}(\widehat{U}_{R})=\mathbf{y}_{\widehat{U}_{\widehat{R}}t}^{\prime}\widehat{\boldsymbol{\beta}}_{\widehat{U}_{\widehat{R}}}$ and Lasso gives $\widehat{y}_{0t}=\widehat{y}_{0t}(\widehat{\boldsymbol{\beta}}_{\widehat{\lambda}}^{las})=\mathbf{y}_{\mathcal{N}t}'\widehat{\boldsymbol{\beta}}_{\widehat{\lambda}}^{las}$. In both factor structures RMPSE of Lasso are larger than those of forward selection in all cases, and Lasso chooses more variables.
Columns 3\textendash 9 of Table (ref) display the rejection probability of the null hypothesis, that is, the proportion of instances when the test rejects the null. The nominal test size is 5%. As the null hypothesis is true in $D1$\textendash $D3$, the rejection probability is associated with test size; the closer it approaches to $5\%$, the better is the performance. For $D4$\textendash $D7$, on the contrary, the larger is the rejection probability, the more powerful is the test. We observe that as the length of the time series increases, the test size based on forward selection falls down toward $5\%$ under both the static and dynamic factor structures, though there is a slight size inflation in $D3$ when dynamics is present in the factors. This is caused by the relatively imprecise long-run variance estimation. The test is powerful under $D4$\textendash $D7$ when the null is violated. In contrast, the test size of the model selected by Lasso is subject to more severe size inflation and is less powerful. The inferior test performance is caused by Lasso's larger RMPSE, which is further caused by the shrinkage estimation scheme.\footnote{Simulation evidence of Lasso's coefficient estimation bias is shown in Table S4 in the online supplement for a sparse model.}
We plot in Figure (ref) the estimated ATE to facilitate visualization. In each panel, the null hypothesis is true for the first column of subgraphs, whereas the null is violated with $E\left[\Delta_{t}\right]=0.5$ for all $t\in\mathcal{T}_{2}$ in the second column and $E\left[\Delta_{t}\right]=1$ in the last column. We witness in both factor structures that forward selection estimates the counterfactual with little bias and the variance is reduced as the time length grows. Finally, the kernel density of the test statistic $\mathcal{Z}_{\widehat{U}_{\widehat{R}}}$ based on forward selection is shown in Figure (ref). Normality is approximated very well in $D1$ and $D2$, though slightly heavier tails are observed in $D3$. Overall, the $t$-statistic graph is supportive for the theoretical result of asymptotic normality.
In this section, we investigate an empirical example where the number of potential control units overpasses the number of temporal observations. Another application which revisits HCW's original empirical example is included in Section S2 of the online supplement.
China launched an anti-corruption campaign of unprecedented scale in November 2012 shortly after Xi Jinping took office. The campaign aimed at cracking down graft and power abuse in all party apparatus, government bureaucracies and military departments. The influence of the anti-corruption campaign motivates academic research assessing its impact from a multitude of perspectives, for example, stock return lin2016anti,ding2017equilibrium and corporate behavior xu2016does,PAN2017. In this paper, we investigate luxury goods importation.
We use the import data from the United Nations Comtrade Database.\footnote{The United Nations Statistics Division, United Nations Comtrade database. \url{http://comtrade.un.org/}.} The database provides detailed statistics for international commodity trade, and the monthly data for China are available since 2010. We focus on the category named “watches with case of, or clad with, precious metal,” following lan2018swiss who find that Chinese luxury watches import co-moves with leadership transitions and government turnover.
The raw time series of Chinese luxury watch import, plotted as the red curve in the lower subgraph in Figure (ref), dropped sharply around the start of the anti-corruption campaign. However, a seemingly structural break can be the upshot of many factors that influenced the macroeconomic environment, for example, terms of international trade, exchange rate volatility, domestic political attitude. During the period from 2013 to 2015, Chinese economy slowed down and it stirred a turmoil over the global commodity markets. Besides the watches, other commodity importation shrank as well. While the flagging economy would have weakened the imports of a myriad of commodities, we employ PDA to control such overall effect in the hope to better isolate the impact of the anti-corruption campaign.
The dependent variable is set as the monthly growth rate of luxury watch import in US dollars, and the independent variables are chosen by the forward selection out of the import growth rates of 88 commodities.\footnote{To ensure that the control units are insusceptible to the anti-corruption policy, 7 categories commonly consumed as bribe goods or conspicuous consumption are excluded. These 7 categories are (with the UN Comtrade Database code in the parenthesis): Beverages, spirits and vinegar (22), Tobacco and manufactured tobacco substitutes (24), Essential oils, perfumes, cosmetics, toiletries (33), Articles of leather, animal gut, harness, travel goods (42), Fur-skins and artificial fur, manufactures thereof (43), Pearls, precious stones, metals, coins, etc (71), Clocks and watches and parts thereof (91) and Works of art, collectors pieces and antiques (97). As a result, $88$ out of the total 95 categories are left to serve as the pool of control units.} We use the growth rate instead of the level data to avoid time series non-stationarity. January 2013 is regarded as the time of the treatment, which is the month immediately after the Eight-Point Policy announcement. There are 35 pre-treatment observations ranging from February 2010 to December 2012, and 36 post-treatment observations spanning from January 2013 to December 2015. The same automated forward selection algorithm as in the simulation chooses 3 control units.\footnote{The selected categories are “knitted or crocheted fabric,” “cork and articles of cork,” and “salt, sulfur, earth, stone, plaster, lime and cement.”}
With the estimated model, we predict the counterfactuals $\widehat{y}_{0t}^{0}$ and estimate treatment effect. Figure (ref) displays the actual luxury watches import growth (solid line) and its estimated counterparts without anti-corruption campaign (dashed line). The upper subgraph shows the growth rate; the lower one shows the value in US dollars, where the counterfactual in monetary value is constructed according to the predicted growth rate. Before the intervention, the model fits the real data quite well and the R-squared of the selected model is 77.85%. After January 2013, if the anti-corruption policy had not been implemented, the import growth rate would have followed the dashed line, which is visibly higher than the realizations. In particular, in January 2013 the import value slumped by a whopping 42%. In contrast, our counterfactual prediction suggests it would have increased by 1.7%. The ATE over the post-treatment period is \[ \frac{1}{36}\sum_{t\in\mathcal{T}_{2}}\widehat{\Delta}_{\widehat{U}_{\widehat{R}}t}=-3.09\%, \] which means that on average the anti-corruption campaign slowed down the luxury watch import by 3.09% per month. The $t$-statistic is $-2.457$, with a $p$-value $1.40\%$. It rejects the null hypothesis of zero ATE at 5% size. Accumulating such a monthly ATE over 36 months leads to roughly two thirds of reduction in importation ($\left(1-0.0309\right)^{36}=0.323$), which is manifested in the lower subgraph. In December 2015, while the realized import was 29.35 million US dollars, the counterfactual predicts 89.27 million without the campaign. Our empirical evidence suggests that China's anti-corruption has been effective in slashing the luxury watch import.
In this paper, we propose using forward selection to choose control units in PDA and then carrying out standard hypothesis testing. Forward selection method is computationally much more efficient than the exhaustive search for the best subset. We establish asymptotic theory for the nearly optimality of forward selection, and show validity of conducting post-selection inference for the ATE by the $t$-statistic based on the selected set. Our theory is valid no matter the true coefficient in the linear regression model is sparse or dense. These extensions widen the applicability of PDA to real world high-dimensional problems in modern data-rich environments.
Shi acknowledges the financial support from the Hong Kong Research Grants Council No.24614817. We thank Cheng Hsiao, Qi Li, Peter Phillips, Jeffrey Wooldridge and Yinchu Zhu for helpful comments, and Yishu Wang for excellent research assistance. All remaining errors are ours.
\setcounter{footnote}{0} \setcounter{table}{0} \setcounter{figure}{0} \setcounter{equation}{0}