An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls
Victor Chernozhukov Kaspar Wüthrich Yinchu Zhu
econ.EM stat.ME


1
  
  
 
 

0

  
  
 


An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls
    Victor ChernozhukovMassachusetts Institute of Technology; 50 Memorial Drive, E52-361B, Cambridge, MA 02142, USA; Email: vchern@mit.edu     Kaspar WüthrichDepartment of Economics, University of California San Diego, 9500 Gilman Dr., La Jolla, CA 92093, USA; Email: kwuthrich@ucsd.edu     Yinchu ZhuBrandeis University; 415 South Street, Waltham, MA 02453, USA; Email: yinchuzhu@brandeis.edu
     
=========================================================================================================================================================================================================================================================================================================================================================================================================




We introduce new inference procedures for counterfactual and synthetic control methods for policy evaluation. We recast the causal inference problem as a counterfactual prediction and a structural breaks testing problem. This allows us to exploit insights from conformal prediction and structural breaks testing to develop permutation inference procedures that accommodate modern high-dimensional estimators, are valid under weak and easy-to-verify conditions, and are provably robust against misspecification. Our methods work in conjunction with many different approaches for predicting counterfactual mean outcomes in the absence of the policy intervention. Examples include synthetic controls, difference-in-differences, factor and matrix completion models, and (fused) time series panel data models. Our approach demonstrates an excellent small-sample performance in simulations and is taken to a data application where we re-evaluate the consequences of decriminalizing indoor prostitution. Open-source software for implementing our conformal inference methods is available.



Keywords: permutation inference, model-free validity, difference-in-differences, factor model, matrix completion, constrained Lasso





1.5 
arabic






§ INTRODUCTION

We consider the problem of making inferences on the causal effect of a policy intervention in an aggregate time series setup with a single treated unit. The treated unit is observed for T_0 periods before and T_∗ periods after the intervention occurs. Often, there is additional information in the form of possibly very many untreated units, which can serve as controls. Such settings are ubiquitous in applied research, and there are many different approaches for estimating the causal effect of the policy. Examples include difference-in-differences methods, synthetic control (SC) approaches, factor, matrix completion, and interactive fixed effects (FE) models, and times series models.[We refer to <cit.>, <cit.>, and <cit.> for excellent comparative overviews and reviews.] 
We refer to these methods as counterfactual and synthetic control (CSC) methods. 

This paper provides generic and robust procedures for making inferences on policy effects estimated by CSC methods. We propose a general counterfactual modeling framework that nests and generalizes many traditional and new methods for counterfactual analysis. We focus on methods that are able to generate mean-unbiased proxies, P_t^N, for the counterfactual outcomes of the treated unit in the absence of the policy intervention, Y_1t^N: 

    Y_1t^N=P_t^N+u_t,    E(u_t )=0,    t=1,…,T_0+T_∗.
 
The policy effect in period t is θ_t=Y_1t^I-Y_1t^N, where Y_1t^I is the counterfactual outcome of the treated unit with the policy intervention. We are interested in testing hypotheses about the trajectory of policy effects in the post-treatment period, θ={θ_t}_t=T_0+1^T_0+T_∗. Specifically, we postulate a trajectory θ^0= {θ^0_t}_t=T_0+1^T_0+T_∗ and test the sharp null hypothesis that θ=θ^0. We also consider the problem of testing hypotheses about per-period effects θ_t and propose a simple algorithm for constructing pointwise confidence intervals via test inversion. 

We recast the inference problem as a (counterfactual) prediction and a structural breaks testing problem. This allows us to build on the literature on conformal prediction <cit.> and end-of-sample stability testing <cit.> to construct inference procedures that are provably robust against misspecification and accommodate many classical and modern high-dimensional methods for estimating P_t^N. 

 The basic idea of our testing procedures is as follows. Suppose that there is only one post-treatment period and that P_t^N is known. Under the sharp null that θ_T_0+1=θ_T_0+1^0, we can compute Y_1t^N and u_t=Y_1t^N-P_t^N for all time periods. If the stochastic shock process {u_t} is stationary and weakly dependent, and its distribution is invariant under the intervention, the distribution of the error in the post-treatment period, u_T_0+1, should be the same as the distribution of the errors in the pre-treatment period, {u_t}_t=1^T_0. We operationalize this idea by proposing inference methods in which p-values are obtained by permuting blocks of estimated residuals across the time series dimension.
 

The proposed methods are valid under two different sets of conditions:
 

  (i) Estimator Consistency and Stationary Weakly Dependent Errors

If the data exhibit dynamics, trends, and serial dependence but the stochastic shock sequence {u_t} is stationary and weakly dependent, our inference procedures are approximately valid if the estimator of P_t^N is consistent (pointwise and in prediction norm). Consistency can be verified for many different CSC methods. We provide concrete sufficient conditions for a representative selection of methods, including difference-in-differences, SC, factor models, matrix completion and interactive FE models, linear and nonlinear time series models, and fused time series panel data models.


  (ii) Estimator Stability and Stationary Weakly Dependent Data

In practice, misspecification is an important concern. We show that even if the model for P_t^N is misspecified and the estimator of P_t^N, P̂_t^N, is inconsistent, our procedures are still valid, provided that the data are stationary and weakly dependent and P̂_t^N satisfies a stability condition. This condition requires that P̂_t^N is stable under perturbations in a few observations. It is implied, for instance, if P̂_t^N is consistent for a pseudo-true parameter value but is shown to hold even in high-dimensional settings where consistency results under misspecification are often not available.



The main theoretical results in this paper are finite sample (non-asymptotic) bounds on the size accuracy of our methods; these bounds imply that our methods are exact as T_0→∞. Unlike traditional asymptotic results, which are only informative when the sample size is large enough, our non-asymptotic bounds show how different factors affect the finite sample performance. This feature is relevant in CSC applications where sample sizes are often small.

A key feature of our conformal inference methods is that P_t^N is estimated under the null hypothesis based on data from all T_0+T_∗ periods. Estimation under the null guarantees the exact finite sample validity of our procedures if the data are iid or exchangeable. Even when exchangeability fails, imposing the null for estimation is essential for a good performance in CSC applications where T_0 is often rather small. Figure <ref> plots the empirical rejection probabilities for testing the null that θ_T_0+1=0 when T_0=19, J=50 (as in our empirical application), {u_t} is an AR(1) process, and P_t^N is estimated using SC. The size properties of our method are excellent. By contrast, estimating P_t^N based on the T_0 pre-treatment periods without imposing the null yields substantial size distortions. Figure <ref> suggests that imposing the null continues to improve size accuracy even when exchangeability fails and that these improvements can be substantial in small samples.



[Figure <ref> around here.]




We make two additional contributions that may be of independent interest. First, we introduce the ℓ_1-constrained least squares estimator or constrained Lasso <cit.> as an essentially tuning-free alternative to existing penalized regression estimators and study its theoretical properties. Constrained Lasso nests SC and difference-in-differences, providing a unifying approach for the regression-based estimation of the mean proxies P_t^N. Second, we obtain theoretical consistency results for SC estimators in settings with potentially very many control units.

We develop three extensions of our main results. First, we show that our method can be modified to test hypotheses about average effects over time. Second, we extend our method to settings with multiple treated units. Third, we propose easy-to-implement placebo tests for assessing the credibility of inferences based on our method. 

Monte Carlo simulations suggest that our procedures exhibit excellent size properties and are robust to misspecification. We find that imposing additional constraints (e.g., using SC instead of the more general constrained Lasso) does not improve power when these additional restrictions are correct but can cause power losses when they are not.

Finally, we re-analyze the causal effect of decriminalizing indoor prostitution on sexually transmitted infections. Following <cit.>, we exploit the unanticipated decriminalization of indoor prostitution in Rhode Island in 2003. We find that decriminalizing indoor prostitution significantly decreased the incidence of female gonorrhea.



  
Related Literature. We contribute to the literature on inference procedures for CSC methods with few treated units.  
A popular method is the finite population permutation approach of <cit.>, see also <cit.> and <cit.>. This approach permutes the policy assignment  and relies on permutation distributions for inference. It corresponds to conventional randomization inference <cit.> under random assignment of the policy <cit.>. However, random assignment is not plausible in typical CSC applications, and assignment mechanisms are difficult to model and estimate when there are only few treated units <cit.>.  <cit.> propose randomization tests for settings with staggered treatment adoption, which encompass the approach of <cit.>. The main assumption of <cit.>'s approach is that policy adoption follows a Cox proportional hazards model. We do not model the assignment mechanism. Instead, we exploit stationarity and weak dependence of the errors across time in a repeated sampling framework. One advantage of exploiting the time series dimension is that we only require a suitable model for the potential outcome of the treated unit. By contrast, cross-sectional approaches often require estimating models for all units. That is, we only require a good “local” instead of a good “global” fit, which reduces the risk of model misspecification. On the other hand, our approach requires a large number of pre-treatment periods and relies on invariance of the error distribution under the intervention.

There is also an active literature on asymptotic inference methods for CSC models. Several papers focus on testing hypotheses about average or expected effects over time, requiring T_0 and T_∗ to be large. <cit.>, <cit.>, <cit.>, and <cit.> introduce inference methods based on penalized and constrained regression methods. <cit.> propose inference methods for a version of SC with time and unit weights, which admits a weighted regression formulation. Asymptotic inference methods based on factor and interactive FE models are proposed by <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. 
Here we focus on sharp null hypotheses and permutation distributions and provide non-asymptotic performance guarantees. Our approach is generic and valid with many different methods, including constrained regression, factor models, and interactive FE estimators. <cit.> propose inference methods for difference-in-difference settings with few treated units. They exploit the cross-sectional dimension, relying on weak dependence and stationarity of the error terms across units, which may be hard to justify in typical CSC settings. By contrast, our procedures rely on stationarity and weak dependence of the errors  over time. On the other hand, exploiting the time series dimension, our approach requires T_0 to be large, whereas <cit.> allow T_0 to be fixed.
In related work, <cit.> provide prediction intervals for per-period effects estimated by SC methods. Their key observation is that there is both randomness from estimating the SC weights and from the prediction error. They propose a sampling-based inference method based on non-asymptotic probability bounds that accounts for both types of randomness and is valid with stationary and non-stationary data.


We recast the causal inference problem as a (counterfactual) prediction problem and build on the literature on conformal prediction <cit.> and on the literature on permutation tests <cit.>, which was started by <cit.> in the context of randomization; see <cit.> for a Bayesian justification. On a more general level, our approach is also connected to transformation-based approaches to model-free prediction <cit.>. Let us discuss in more detail the relationship to <cit.>, who extend classical conformal prediction to time series settings. Besides a different focus (prediction intervals for future outcome values vs. inference on policy effects), there are several important differences. First, we rely on permuting residuals, whereas CWZ18 permute blocks of data. Second, we theoretically analyze different types of permutations. In particular, we study the set of all permutations, which yields precise p-values in small samples. This set of permutations cannot be used in the framework of CWZ18 unless the data are iid. Third, we allow for non-stationary data, whereas the prediction methods in CWZ18 are strictly limited to stationary data. Forth, CWZ18 rely on abstract high-level conditions on the test statistics and do not provide any primitive conditions. By contrast, we develop transparent sufficient conditions that can be verified for many traditional and modern CSC methods, and we provide explicit primitive conditions for a large selection of popular approaches. Finally, we establish the validity of our methods with time series data under misspecification and stability, whereas the theoretical results for weakly dependent data in CWZ18 require correct specification and consistency. 
 
 
Finally, we show that the problem of making inferences on policy effects can be recast as a structural breaks testing problem with a known break date. Therefore, we build on and contribute to the literature on testing for structural breaks and, in particular, to the literature on structural breaks testing using permutation approaches <cit.>. Besides a different focus (inference on policy effects vs. testing for structural breaks), our paper differs from the existing literature in that we specifically focus on testing at the end of the sample, allow for a very general class of estimators, including modern high-dimensional methods, and provide non-asymptotic performance guarantees under correct specification and misspecification. Our paper is also related to tests for structural breaks at the end of the sample <cit.>.[<cit.> informally suggest applying a variant of <cit.>'s end-of-sample stability test in the context of SC, and <cit.> use a version of this test in the context of difference-in-differences approaches with few treated groups.] Let us discuss the differences to <cit.>'s end-of-sample instability test based on subsampling in more detail. First, we focus on causal inference, whereas <cit.> is concerned with structural breaks testing. Second, our procedures are exactly valid under exchangeability, and we obtain finite sample bounds under weak conditions on the estimators, while the theoretical properties of <cit.>'s test rely on asymptotic analyses. Third, our methods are valid under misspecification, whereas <cit.> assumes correct specification. Forth, our results under correct specification only require stationarity and weak dependence of {u_t}, while <cit.>'s test assumes stationarity of the data.[ <cit.> briefly comments on page 1681 (comment 4) that his test can be shown to be asymptotically valid under stationary errors but does not provide a formal result.] Finally, our procedures work in conjunction with many modern high-dimensional estimators, whereas <cit.> focuses on low-dimensional GMM models. 
 


  
Notation. For q≥ 1, the ℓ_q-norm of a vector is denoted by ·_q. We use ·_0 to denote the number of nonzero entries of a vector; ·_∞ is used to denote the maximal absolute value of entries of a vector. We use the
notation a≲ b to denote a ≤ cb for some constant c > 0 that does not depend on the sample size. We use the notation a ≍ b to denote a≲ b and b≲ a. For a set A, |A| denotes the cardinality of A.  For any a∈ℝ, we define ⌊ a ⌋ =max{z∈ℤ:z≤ a}  and ⌈ a ⌉ =min{z∈ℤ:z≥ a}, where ℤ is the set of integers. We use ℕ to denote the set of natural numbers.


 



§ A CONFORMAL INFERENCE METHOD




 §.§ The Counterfactual Model


We consider a time series of T outcomes for a treated unit, labeled j=1.  During the first
T_0 periods, the unit is not treated by a policy and, during the remaining T-T_0=T_∗ periods,  it is treated
by the policy. Extensions to more than one treated unit are discussed in the Appendix. Our typical setting is where T_∗ is short compared to T_0. There may be other units that are not exposed to the policy, and they will be introduced below. We denote the observed outcome of the treated unit by Y_1t. We employ the potential (latent) outcomes framework <cit.> and denote potential outcomes with and without the policy as Y_1t^I and Y_1t^N. The effect of the policy intervention in period t is θ_t=Y_1t^I-Y_1t^N. 


Our conformal inference method will rely on the following counterfactual modeling framework, which nests many traditional and new methods for counterfactual policy analysis; see Sections <ref>–<ref> for examples.


[Counterfactual Model] 
Let {P_t^N} be a given sequence of mean-unbiased predictors or proxies for the counterfactual outcomes {Y_1t^N} in the absence of the policy intervention, that is {E( P_t^N)} = {E ( Y_1t^N)}. Let {θ_t} be a fixed policy effect sequence with θ_t = 0 for t ≤ T_0, so that potential outcomes under the intervention are given by {Y_1t^I}  = {Y_1t^N + θ_t}.[In the Appendix, we consider an extension to random policy effects.]
In other words, potential outcomes can be written as

    CMF[       Y_1t^N = P_t^N + u_t; Y_1t^I = P_t^N + θ_t + u_t;                            ] |    E( u_t) = 0,     t=1,…,T ,

where {u_t} is a centered stationary stochastic process. Observed outcomes are related to potential outcomes as Y_1t=Y_1t^N+D_t(Y_1t^I-Y_1t^N), where D_t=1(t>T_0 ). 

Assumption <ref> introduces the potential outcomes, but also postulates an identifying assumption in the form of the existence of mean-unbiased proxies P^N_t such that E ( P_t^N) = E ( Y_1t^N ).  Assumption <ref> allows {P_t^N} to be fixed or random and does not impose any restrictions on the dependence between {P_t^N} and {u_t}. In Sections <ref>–<ref>, we will discuss specific panel data and time series models that postulate (and identify) what P_t^N is under a variety of conditions.  Additional assumptions on the stochastic shock process {u_t} will be introduced later, in essence requiring {u_t} to be either iid or, more generally, a stationary and weakly dependent process. 

Assumption <ref> also postulates that the stochastic shock sequence is invariant under the intervention.  This is the fundamental identifying assumption. It requires that the timing of the policy intervention is independent of factors that change the distribution of {u_t}.[In principle, we can relax this assumption by specifying, for example, the scale and quantile shifts in the stochastic shocks that result from the policy, and then working with the resulting model; we leave this extension to future work.] If the policy changes the distribution of {u_t}, one can either interpret our method as a structural breaks test or view the policy effect as random in which case our method yields valid prediction sets; see the Appendix for details.

Often, there is additional information in the form of untreated units, which can serve as controls. Specifically, suppose that there are J ≥ 1 control units, indexed by  j=2,…,J+1.  We assume that we observe all units for all T periods, although this assumption can be relaxed. Let Y_jt denote the observed outcome for these untreated units.  This observed outcome is equal to  the outcome in the absence of the policy intervention, i.e., Y_jt = Y_jt^N for 2≤ j≤  J+1 and 1≤ t≤ T. For each unit, we may also observe a vector of covariates X_jt. This motivates a variety of strategies for modeling and identifying P_t^N as discussed below.





 §.§ Hypotheses of Interest, Test Statistics, and p-Values


We are interested in testing hypotheses about the trajectory of policy effects in the post-treatment period, θ=(θ_T_0+1,…,θ_T)'. Our main hypothesis of interest is

    H_0: θ=θ^0,

where θ^0=(θ^0_T_0+1,…,θ^0_T)' is a postulated policy effect trajectory. Hypothesis (<ref>) is a sharp null hypothesis. It fully determines the value of the counterfactual outcome in the absence of the intervention in the post-treatment period since Y_1t^N=Y^I_1t-θ_t=Y_1t-θ_t. In the Appendix, we show that our method can also be used to test hypotheses about average effects.

To describe our procedure, we write the data under the null hypothesis as :=(θ^0)=(Z_1,…,Z_T)', where

    Z_t=(Y^N_1t,Y^N_2t,…,Y^N_J+1t,X'_1t,…,X'_J+1t)',     t≤ T_0
    (Y_1t^I-θ_t^0,Y^N_2t,…,Y^N_J+1t,X'_1t,…,X'_J+1t)',     t>T_0.



Using one of the methods described below, we will obtain a counterfactual proxy estimate, P̂^N_t,
based on , and compute the residuals û=(û_1,…,û_T )', where û_t=Y_1t^N-P̂^N_t for 1≤ t≤  T.
Since P_t^N is computed using =(θ^0),  P_t^N is estimated under the null hypothesis, which is essential for a good small sample performance. In “ideal” settings where the data are iid, imposing the null guarantees the model-free exact finite sample validity of our method; see the Appendix for details. By contrast, when P_t^N is estimated based on the pre-treatment data {Z_t}_t=1^T_0 without imposing the null, permuting blocks of residuals does not yield procedures with exact finite sample validity, not even with iid data.



Definition of Test Statistic 𝐒. We consider the following test statistic:

    S  (û) = S_q(û) =  ( 1/√(T_∗)∑_t=T_0+1^T|û_t|^q )^1/q.




Note that S is constructed such that high values indicate rejection. Different choices of q lead to power against different alternatives. For instance, if the intervention has a large but only temporary effect (i.e., if |θ_t| is large for few periods), choosing q=∞ yields high power. On the other hand, if the intervention has a permanent effect (i.e., if θ_t is non-zero for many post-treatment periods), tests using S_1 or S_2 exhibit good power properties. In our application, we will be using S_1, which behaves well under heavy-tailed data. Throughout the paper, when the nature of the statistic is not essential, we write S = S_q. 

 While we focus on S_q, other test statistics can be used as well.  For example, when capturing deviations in the average effect T_*^-1∑_t=T_0+1^T θ_t, it is useful to consider S(û) = T_∗^-1/2 |∑_t=T_0+1^T û_t  |. 

We use (block) permutations to compute p-values. A permutation π is a one-to-one mapping π:{1,…,T}↦{1,…,T}. We denote the set of permutations under study as Π and assume that Π contains the identity map 𝕀. We focus on two different sets of permutations: (i) the set of all permutations, which we call iid permutations,  Π_all, and (ii) the set of all (overlapping) moving block permutations, Π_→.[We can also consider other types of permutations; for example, the “iid block” permutations. Specifically, let {b_1, …, b_K} be a partition of {1,…,T}, then we collect all the permutations π of these blocks, forming the “iid m-block” permutations Π_mb. In our context, choosing m=T_* is natural, though other choices should work as well, similarly to the choice of block size in the time series bootstrap. We refer to CWZ18 for more results on block permutations.] The elements of Π_→ are indexed by j ∈{0,1, …, T-1}, and the permutation π_j is defined as

    π_j(i)=
    i+j    if i+j≤ T
    
    i+j-T   otherwise.

Figure <ref> presents a graphical illustration of Π_all and Π_→.



[Figure  <ref> around here.]



The choice of Π does not matter for the exact finite sample validity of our procedures if the residuals are exchangeable. However, Π_all has more elements than Π_→, allowing us to compute more precise p-values and to test at lower significance levels. For the asymptotic validity under estimator consistency, the choice of Π depends on the assumptions that we are willing to impose on the stochastic shock sequence {u_t} (cf. Section <ref>).



For each π∈Π, let û_π=(û_π(1),…,û_π(T))' denote the vector of permuted residuals.[If the estimator of P_t^N is invariant under permutations of the data {Z_t} across the time series dimension (which is the case for many estimators in Section <ref>), permuting the residuals {û_t} is equivalent to permuting the data {Z_t}.] The permutation p-value is defined as follows.




Definition of p-Value. The p-value is 

    p̂=1- F̂( S(û)),   where F̂( x)=1/|Π|∑_π∈Π1{ S(û_π)<  x}.



We are often interested in testing pointwise hypotheses about θ_t, H_0:θ_t=θ_t^0, and in constructing pointwise confidence intervals for θ_t. Pointwise hypotheses can be tested by defining the data under the null as =(Z_1,…,Z_T_0,Z_t)', provided that P_t^N can be estimated based on . Pointwise (1-α) confidence intervals for θ_t can be constructed via test inversion as described in Algorithm <ref>.

[Pointwise Confidence Intervals] (i) Choose a fine grid of G candidate values Θ̃_t={θ̃_1t^0,…,θ̃_Gt^0}. (ii) For θ̃_t^0∈Θ̃_t, define  for the null hypothesis H_0:θ_t=θ̃_t^0 and compute the corresponding p-value, p̂ (θ̃_t^0), using (<ref>). (iii) Return the (1-α) confidence set 𝒞_1-α(t)={θ̃_t^0∈Θ̃_t:  p̂ (θ̃_t^0)> α}. 






 §.§ Models for Counterfactual Proxies via Synthetic Control and Panel Data


The availability of control units motivates several strategies for modeling the counterfactual mean proxies P_t^N. We estimate P_t^N based on the imputed data under the null hypothesis, (θ^0), and write Y_1t^N instead of Y^I_1t-θ_t^0 to alleviate the exposition. 



  §.§.§ Difference-in-Differences Methods


The difference-in-differences method postulates the following model for the counterfactual mean proxy <cit.>:
P_t^N=μ+1/J∑_j=2^J+1Y^N_jt.
This model automatically embeds the identifying information. The counterfactual mean proxy can be estimated as 
P̂^N_t=1/T∑_s=1^T(Y^N_1s-1/J∑_j=2^J+1Y^N_js)+1/J∑_j=2^J+1Y^N_jt.




  §.§.§ Synthetic Control and Constrained Lasso


The canonical SC method <cit.> postulates the following model:

    P^N_t=∑_j=2^J+1w_jY^N_jt,    where   w≥ 0   and  ∑_j=2^J+1 w_j=1.

We need to impose an identification condition that allows us to identify the weights w, for example:[More generally, other exclusion restrictions and identifying assumptions could be used. See also <cit.>, <cit.> and <cit.>, who study the behavior of SC when the data are generated by a factor model.]
 


  (SC) Assume that the structural shocks u_t for the treated unit are uncorrelated with 
contemporaneous values of the outcomes, namely: E ( u_t Y^N_jt)= 0 for  2≤ j≤ J+1.



 The counterfactual is estimated as P̂^N_t=∑_j=2^J+1ŵ_jY^N_jt.
We focus on the following canonical SC estimator for w:[This formulation of canonical SC
without covariates is due to  <cit.>, who refer to the estimator (<ref>) as “constrained regression”. Note that unlike <cit.>, we estimate w under the null hypothesis based on all the data. We focus on the canonical problem (<ref>) for concreteness. <cit.> consider a more general version that also includes covariates into the estimation of the weights. Our inference method also works in conjunction with more recently proposed modified versions of SC, such as the augmented SC estimator of <cit.>.]

    ŵ=min_w∑_t=1^T(Y^N_1t-∑_j=2^J+1w_jY^N_jt)^2  s.t.   w≥ 0   and  ∑_j=2^J+1w_j=1.


As an alternative, we can consider the more flexible model[The idea to relax the non-negativity constraint on the weights is not new. It first appeared in <cit.>, who compared their factor model approach to SC, and also in <cit.>, who used the cross-validated Lasso to estimate the weights, and in <cit.>,  who used cross-validated Elastic Net for estimation of weights. They do not establish the formal properties of these estimators. Here we emphasize another version of relaxing SC, model (<ref>), which leads to constrained Lasso (<ref>). Constrained Lasso demonstrates an excellent theoretical and practical performance:  it is tuning-free, performs very well empirically and in simulations, and we prove that it is consistent for dependent data without any sparsity conditions on the weights and that it satisfies the estimator stability condition required for validity under misspecification. We emphasize that this estimator generally differs from the cross-validated Lasso estimator.]

    P^N_t=μ + ∑_j=2^J+1w_jY^N_jt,   where   w_1 ≤ 1,

maintaining the same identifying assumption (SC). The counterfactual is estimated as P̂^N_t=μ̂+∑_j=2^J+1ŵ_jY^N_jt 
by the ℓ_1-constrained least squares estimator, or constrained Lasso <cit.>:

    (μ̂,ŵ)=min_(μ,w)∑_t=1^T(Y^N_1t-μ-∑_j=2^J+1w_jY^N_jt)^2   s.t.   w_1≤ 1.


The advantage over other penalized regression methods discussed next is that constrained Lasso is essentially tuning-free, does not rely on any sparsity conditions, and is valid for dependent data under weak assumptions. Moreover, constrained Lasso encompasses both difference-in-differences and canonical SC as special cases (by setting w=(1/J,…,1/J)' and μ=0,w≥ 0, respectively) and, thus, provides a unifying approach for the regression-based estimation of P^N_t. 

Section  <ref> provides primitive conditions that guarantee that the SC and the constrained Lasso estimators are valid in our framework in settings with potentially many control units (large J). Finally, we note that it is straightforward to incorporate (transformations of) covariates X_jt into the estimation problems (<ref>) and (<ref>).



  §.§.§ Penalized Regression Methods


Consider a linear model for P_t^N:
P^N_t=μ+∑_j=2^J+1w_jY^N_jt.
We maintain the identifying assumption (SC). The counterfactual is estimated by P̂^N_t=μ̂+∑_j=2^J+1ŵ_jY^N_jt,
where 

    (μ̂,ŵ)=min_(μ,w)∑_t=1^T(Y^N_1t-μ-∑_j=2^J+1w_jY^N_jt)^2 + 𝒫(w),

and 𝒫(w) is a penalty function that penalizes deviations away from zero.  If it is desired
to penalize deviations away from other focal points w^0, for example, w^0 =(1/J, … , 1/J) used in the difference-in-differences approach, we may always use instead: 𝒫(w) ←𝒫(w- w^0). Note that it is straightforward to incorporate covariates X_jt into the estimation problem (<ref>).  

Different variants of 𝒫(w) can be considered. Examples include: Lasso <cit.>, where 𝒫(w)=λw_1 and λ is a tuning parameter; Elastic Net <cit.>, where 𝒫(w)=λ( (1-α)w_2^2+αw_1 ) and λ and α are tuning parameters; Lava <cit.>, where 𝒫(w)=inf_a + b = wλ((1-α) a_2^2+αb_1 ) and λ and α are tuning parameters.

In the context of CSC methods, Lasso was used by <cit.>, <cit.>, and <cit.>, while <cit.> proposed to use Elastic Net. We will impose only weak requirements on the performance of the estimators (pointwise consistency and consistency in prediction norm), which implies that these estimators are valid in our framework under any set of sufficient conditions that exists in the literature. 



  §.§.§ Interactive Fixed Effects, Factor, and Matrix Completion Models


Consider the following interactive FE model for treated and untreated units:

    [ Y^N_jt=λ_j'F_t+ X_jt'β + u_jt,  for   1≤ j ≤ J+1 and  1≤ t≤ T, ]

where F_t are unobserved factors, λ_j are unit-specific factor loadings, and β is a vector of common coefficients. Model (<ref>) nests the classical factor model when β=0 and also covers the traditional linear FE model, in which λ_j'F_t =  λ_j+F_t. Consider the following assumption.

  

  (FE)  Assume that u_jt is uncorrelated with (X_jt,F_t,λ_j),
as well as other identification conditions in <cit.>. 



The model leads to the following proxy:

    P_t^N=λ_1'F_t + X_1t'β.

Counterfactual proxies are estimated by P̂_t^N=λ̂_1'F̂_t + X_1t'β̂, where λ̂_1 and F̂_t, and β̂ are obtained using the alternating least squares method  applied to the model (<ref>); see, for example, <cit.> and <cit.> for a version with high-dimensional covariates. 


<cit.> appears to the be first work that proposed the use of factor models for predicting the (missing) counterfactual responses specifically in SC settings. <cit.> and <cit.> employ <cit.>'s estimator in this setting, albeit provide no formally justified inference methods. Formal inference results for interactive FE and factor models in SC designs are developed in <cit.> and <cit.> among others.[Factor models are widely used in macroeconomics for causal inference and prediction; see, for example, <cit.> and the references therein. In microeconometrics, factor models are used for estimation of treatment/structural effects; see, for example, <cit.> who use interactive FE models to estimate the effect of gun prevalence on crime.]

Other recent applications to predicting counterfactual responses include <cit.> and <cit.> (using, respectively, singular value thresholding and the nuclear norm penalization).[Note that  <cit.>'s analysis applies to a broader collection of problems with general missing data patterns, nesting SC and difference-in-difference problems as special cases.] Our method delivers a way to perform valid inference for policy effects using any of the factor model estimators used in these proposals applied to the complete data under the null.[Note that in our case the sharp null allows us to impute the missing counterfactual response and apply any of the factor estimators to estimate the factor model
for the entire data, which is then used for conformal inference. Hence our inference approach does 
not provide inference for the counterfactual prediction methods given in those papers. Indeed, there, the missing data entries are being predicted using factor models, whereas in our case the missing data entries are known under the null, and we use any form of low-rank approximation or interactive FE model to estimate the model for the entire data under the null hypothesis.] We shall be focusing on <cit.>'s alternating least squares estimator[We choose to focus on PCA/SVD and the alternating least squares estimator for the following reasons: (1) they are by far the most widely used in practice, (2) the alternating least squares estimator is computationally attractive and easily accommodates unbalanced data.] and on matrix completion via nuclear norm penalization when verifying our conditions.



 §.§ Models for Counterfactual Proxies via Time Series and Fused Models




  §.§.§ Simple Time Series Models
 
If no control units are available, one can use time series models for the single unit exposed to the intervention. For example, consider the following autoregressive model:[We can also add a moving average component for the errors, but we do not do so for simplicity.]

    [         Y_1t^N - μ =  ρ (Y_1(t-1)^N -μ ) + u_t; Y_1t^I  - μ = ρ (Y_1(t-1)^N - μ ) +  θ_t + u_t ] |    E( u_t) = 0,   {u_t} iid,     t=1,…,T.

In model (<ref>), the mean unbiased proxy is given by
P_t^N =  μ + ρ (Y_1(t-1)^N- μ).
Note that the policy effect here is transitory, namely it does not feed-forward itself on the future values of Y_1t^I beyond the current values.[We leave the model with persistent feed-forward effects, 
Y_1t^I = ρ (Y_1(t-1)^I) + θ_t + u_t, to future work. ]  Under the null hypothesis, we can impute the unobserved counterfactual  as Y_1t^N = Y_1t - θ_t  and estimate the model using traditional time series methods, and we can conduct inference by permuting the residuals.



The simplest form of the autoregressive model is the AR(K) process, where the ρ(·) take the form: ρ ( · )  =  ∑_k=0^K ρ_k L^k (· ), where L is the lag operator.  There are many identifying conditions
for these models, see, for example, <cit.> or <cit.>. More generally, we can use a nonlinear function of lag operators, ρ (· ) = m( ·, L^1 (·), …, L^k (· )), as, for example, when applying neural networks to time series data <cit.>, and we refer to the latter for identifying conditions.




  §.§.§ Fused Time-Series/Panel Models


A simple and generic way to combine the insights from the panel data and time series models is as follows. Consider the system
of equations:

    [       Y_1t^N = C_t^N + ε_t; Y_1t^I = C_t^N + θ_t + ε_t;                            ] |  
    [ ε_t = ρ (ε_t-1) + u_t,  {u_t} iid,   E( u_t) = 0,;                  {u_t} is independent of {C_t^N}, ] |    t=1,…,T,

where  C_t^N is a panel model proxy for Y_1t^N, identified by one of the panel data methods.  Note that the model has the autoregressive formulation:
Y_1t^N = C_t^N  + ρ (Y_1(t-1)^N - C_t-1^N) + u_t,
thereby generalizing the previous model. 

Here the mean unbiased proxy for Y_1t^N is given by P_t^N = C_t^N + ρ (ε_t-1).
P_t^N is a better proxy than C_t^N because it provides an additional noise reduction through prediction of the stochastic shock by its lag. The model combines any favorite panel model C_t^N for counterfactuals with a time series model for the stochastic shock model in a nice way:  we can identify C_t^N under the null by ignoring the time series structure, and  then  identify the time series structure of  the residuals Y_1t^N - C_t^N.  Estimation can proceed analogously. This approach will often improve the size accuracy of our inferential procedures.






§ THEORY



When the data are iid (or exchangeable), our procedure is exactly valid in finite samples as shown in the Appendix. In this section, we establish the validity of our inference methods with time series data. Our results are non-asymptotic in nature and, hence, hold in finite samples. Finite sample bounds are provided for the size properties of our procedure; these bounds imply that our approach is exact as T_0→∞. In Section <ref>, we establish the validity of our procedure when the estimator of P_t^N satisfies weak and easy-to-verify small error conditions (pointwise consistency and consistency in the prediction norm). This result accommodates non-stationary data and only requires stationarity and weak dependence of the stochastic shock process {u_t}. In Section <ref>, we consider a setting that accommodates misspecification and inconsistent estimators. We show that if the data are stationary and weakly dependent, our procedure is valid, provided that the estimators are stable.  



 §.§ Approximate Validity under Estimator Consistency



The main condition underlying the results in this section is the following assumption on the stochastic shock process.

[Regularity of the Stochastic Shock Process]
 Assume that the density function of S(u) exists and is bounded, and that the stochastic process {u_t}_t=1^T satisfies one of the following conditions.
   

  * {u_t}_t=1^T are iid,  or  

  * {u_t}_t=1^T are stationary, strongly mixing, with sum of mixing coefficient bounded by M.  



Assumption <ref> allows the data to be non-stationary and exhibit general dependence patterns. Assumption <ref>.<ref> of iid shocks is our first sufficient condition. Under this condition, we will be able to use iid permutations, giving us a precise estimate of the p-value. The iid assumption can be replaced by Assumption <ref>.<ref>, which holds for many commonly encountered stochastic processes such as ARMA and GARCH. It can be easily replaced by an even weaker ergodicity condition, as can be inspected in the proofs. Under this assumption, we will have to rely on the moving block permutations. 

 Assumption <ref> does not rule out conditional heteroscedasticity in the stochastic shock process {u_t}. Unconditional heteroscedasticity is allowed in {Z_t} but not in {u_t}. When we suspect unconditional heteroscedasticity in {u_t}, we can apply another filter or model to obtain “standardized residuals” from {û_t}. This will generally require another layer of modeling assumptions, leading to an overall procedure that reduces the data to “fundamental” shocks that are assumed to be stationary under the null.  

We also impose the following condition on the estimation error under the null hypothesis. Let P^N=(P^N_1,…,P^N_T)' and P̂^N=(P̂^N_1,…,P̂^N_T)'.

[Consistency of the Counterfactual Estimators under the Null] Let there be sequences of constants δ_T and γ_T converging to zero. Assume that with probability 1- γ_T,
   

  * the mean squared estimation error is small, P̂^N - P^N ^2_2/T ≤δ^2_T;

  * for T_0+1≤ t≤ T, the pointwise errors are small, |P̂^N_t- P^N_t| ≤δ_T.



Assumption <ref> imposes weak and easy-to-verify conditions on the performance of the estimators P̂_t^N of the counterfactual
mean proxies P_t^N.  These conditions are readily implied by the existing results for many estimators discussed in Section <ref>. In Section <ref>, we provide explicit primitive conditions and references to primitive conditions implying Assumption <ref>.[While our general results in this section are non-asymptotic, some of the analysis in Section <ref> will not be non-asymptotic in nature.]


  

Assume that T_* is fixed. Suppose that Assumptions <ref> and <ref> hold. Impose  Assumption <ref>.<ref> if Π=Π_all; impose Assumption <ref>.<ref> if Π=Π_→. Assume S(u) has a density function bounded by D under the null. Then, under the null hypothesis, the p-value is approximately unbiased in size:

    |P(p̂≤α)- α| ≤ C ( δ̃_T + δ_T + √(δ_T) + γ_T),

where δ̃_T = (T_*/T_0)^1/4(log T) and the constant C depends on T_*, M and D,  but not on T.


The above bound is non-asymptotic, allowing us to claim uniform validity with respect to a rich variety of data generating processes.  Using simulations and empirical examples, we verify that our tests have good power and generate meaningful empirical results. There are other considerations  that also affect  power. For example, the better the model for P_t^N, the less variance the stochastic shocks will have, subject to assumed invariance to the policy. The smaller the variance of the shocks, the more powerful the testing procedure will be.




 §.§ Approximate Validity under Estimator Stability



Misspecification is an important practical concern, and consistency of the estimators of the counterfactual mean proxies P_t^N may be questionable in certain settings. The classical analysis of misspecification focuses on convergence to pseudo-true values <cit.>. If it is possible to show that the estimator of the counterfactual mean proxy, P̂_t^N, is consistent for some pseudo-true value P_t^N∗ and that {Y_1t^N-P_t^N∗}_t=1^T is stationary and weakly dependent, the theoretical results in Section <ref> imply the validity of our procedure. Pseudo-true consistency can often be verified for low-dimensional models, but consistency results under misspecification remain elusive in high-dimensional settings. Therefore, we consider a notion of approximate exchangeability, which only requires the estimator to be stable instead of consistent for a pseudo-true value. This stability condition does not require P̂_t^N to be consistent for anything, nor does it rely on correct specification of the counterfactual mean proxies. In the Appendix, we illustrate the difference between consistency and stability based on the analytically tractable example of Ridge regression. 

The basic idea underlying the theoretical analysis here is as follows. If the estimators are non-random or independent of the data, then stationarity and weak dependence of the data would mean that p̂ based on moving block permutations approximately has a uniform distribution under the null. This result follows from uniform laws of large numbers for dependent data. However, in practice, the estimators are computed using the data and are thus not independent of the data. Our key insight is that stable estimators are approximately independent of individual observations. 

We now formalize the notion of stability of an estimator. To emphasize the dependence of S(û) on the estimator, with a slight abuse of notation, we write S(,β)=ϕ(Z_T_0+1,…,Z_T_0+T_*;β). 
Let {Z̃_t}_t=1^T be iid from the distribution
of Z_1 and independent of . For any H⊂{1,…,T},
let Z_t,H=Z_t{t∉ H}+Z̃_t{t∈ H},
and _H={Z_t,H}_t=1^T. Hence, _H is a perturbed
version of  under H, i.e.,  with elements in H replaced
by {Z̃_t}_t∈ H. 

By stability, we mean that the estimator computed using  is similar to that computed using  _H for H∈ℍ. Let R∈ℕ and define m=⌊ T_0/R⌋. The class ℍ={_1,…,_R} contains R members with |_j|≤ 3m elements. The plan is to require stability under R≍ T_0/log(T_0) (so |_j|≍log(T_0)). Since log(T_0)≪ T_0, swapping out O(log(T_0))  out of T_0+T_* data points should not cause a large change in the estimator for reasonable estimators. 

We now give precise definitions of sets in ℍ. For j∈{1,…,R}, let H_j={(j-1)m+1,…,jm}. Since the test statistic depends on T_* data points after obtaining the estimator,  defining ℍ to be {H_1,…,H_R } is not enough for technical arguments; we need a “wedge” to ensure that these T_* data points do not cause a problem. To do so, we enlarge H_j as follows.  Let k∈ℕ
satisfy T_*<k<m. We let _j denote the k-enlargement
of H_j, i.e., _j={s: min_t∈ H_j|s-t|≤ k}.
Note that _j={(j-1)m+1-k,…,jm+k} for
2≤ j≤ R-1, _1={1,…,m+k} and _R={(R-1)m+1-k,min{Rm+k,T}}.





[Estimator Stability]
Let Π= Π_→. There exist non-decreasing functions ϱ_T(·)
such that 
P(max_π∈Π|S(^π,β̂())-S(^π,β̂(_H))|≤ϱ_T(|H|))≥1-γ_1,T 
and 

P(max_π∈Π|S(()^π,β̂())-S(()^π,β̂(_H))|≤ϱ_T(|H|))≥1-γ_1,T
for any H∈{_1,…,_R}, where d=
and  is independent of (,{Z̃_t}_t=1^T). 



Assumption <ref> specifies the estimator stability condition. It strengthens the perturb-one sensitivity of <cit.>. When the model is misspecified, Assumption <ref> holds whenever the estimator β̂() is consistent to a pseudo-true parameter value. However, it is more general in that the estimator β̂() need not converge to any non-random quantity as long as it is stable under perturbations in a few observations. This feature is crucial in our setting as it allows us to accommodate high-dimensional CSC methods for many of which consistency results under misspecification are not available. Primitive sufficient conditions for Assumption <ref> are provided in the Appendix. 


Let Ψ(x;β)=P(ϕ(Z_T_0+1,…,Z_T_0+T_*;β)≤ x). Our strategy is to show that, under the null hypothesis,
(ϕ(Z_T_0+1,…,Z_T_0+T_*;())) is approximately uniform on (0,1). We exploit the stability condition in Assumption <ref> and show that (ϕ(Z_T_0+1,…,Z_T_0+T_*;())) can be approximated by Ψ(ϕ(_T_0+1,…,_T_0+T_*;(__R));(__R)), which has the uniform distribution on (0,1). Here (_T_0+1,…,_T_0+T_*) has the same distribution as (Z_T_0+1,…,Z_T_0+T_*) and is independent of __R. This essentially confirms the above intuition that for stable estimators, () is almost independent of the last few observations (Z_T_0+1,…,Z_T_0+T_*).



We impose the following regularity conditions on the data.


[Regularity of the Data]
The data under the null, {Z_t}_t=1^T, are stationary
and β-mixing with coefficient (·) satisfying
(i)≤ D_1exp(-D_2i^D_3) for some constants D_1,D_2,D_3>0.
 For 1≤ j≤ R, there exist  sequences ξ_T>0 and γ_2,T=o(1) such that
P(sup_x∈|∂Ψ(x;(__j))/∂ x|≤ξ_T)≥1-γ_2,T.


Stationarity and β-mixing are commonly
imposed conditions on time series data.  For a large class of Markov chains,  GARCH and various
stochastic volatility models, D_3=1 <cit.>.
Let (Ż_T_0+1,…,Ż_T_0+T_*) be an independent copy of  (Z_T_0+1,…,Z_T_0+T_*) and also independent of (,{Z̃}_t=1^T).    The bounded derivative of Ψ(x;(__j)) condition says that the density of ϕ(Ż_t,…,Ż_t+T_*-1;(__j)) conditional on (__j) is bounded by ξ_T with high probability. 
The bounded density condition states that the distribution of the residual does not collapse into a degenerate one or one with point mass. In many cases, ξ_T=O(1) for continuous distributions. For example, if (Y_t,X_t) is jointly Gaussian and the variance of Y_t given X_t is bounded below by a constant, then for any w satisfying the SC restrictions, the density of Y_t-X_t'w is bounded by a constant that does not depend on w.


The following result states the approximate validity of our testing
procedure. 
 
 Let Π=Π_→. Suppose that Assumptions <ref>
and <ref> hold.  Then, under the null hypothesis, there exists a constant
C_1>0 depending only on D_1, D_2 and D_3 such
that for any R with k<⌊ T_0/R⌋ and R<T_0/2,

    |P(p̂≤α)-α|≤ C_1√(ξ_Tϱ_T(T_0/R+2k))+C_1(T_0^-1R[log(T_0/R)]^1/D_3)^1/4
    
    +C_1exp(-(k-T_*+1)^1/D_3)+C_1√(γ_1,T)+C_1√(γ_2,T).



In the theoretical arguments, we actually show a stronger result. The above bound holds for E|P(p̂≤α|(__R))-α|. Since the stability condition states that (__R)≈(), this means that p̂ conditional on () almost has a uniform distribution on (0,1); with iid or exchangeable data, p̂ conditional on ()  has an exact uniform distribution. Therefore, we can view Theorem <ref> as a result for approximate exchangeability. 

Due to the exponential decay of (·), the bound in Theorem
<ref> tends to zero if  we  choose
k to be a slowly growing sequence and T_0/R to be of the same order. For example, we can choose k and R such that k≍ T_0/R ≍log T_0. Since |_j|=⌊ T_0/R⌋ +2k,
Assumption <ref> only requires that the changes
to S(^π,β̂()) are small if we replace only  log T_0
observations in computing β̂(). Under finite dependence, it suffices to choose k and T_0/R to be large enough constants.  Note that R is only needed in the theoretical arguments;  we do not need to choose R when implementing the proposed procedure.

The theoretical analysis in this section suggests that allowing for both unrestricted patterns of non-stationarity and misspecification is not possible in general. To obtain valid inferences with non-stationary data, one has to either rely on correct specification and consistency or impose assumptions on the particular structure of the non-stationarity, which allow for pre-processing the data to make them stationary.


 




§ SUFFICIENT CONDITIONS FOR CONSISTENT ESTIMATION



In this section, we revisit the representative models of counterfactual proxies introduced in Section <ref>. Primitive conditions are provided to guarantee that the estimation of the counterfactual mean proxies is accurate enough for the asymptotic validity of the proposed procedure. In particular, these conditions can be used to verify Assumption <ref>. The regularity conditions (e.g., bounded moments, weak serial dependence) for different models are stated in the Appendix and are commonly imposed in the literature for these models.  The counterfactual mean proxies P_t^N are estimated based on the imputed data under the null, (θ^0), and we write Y_1t^N instead of Y^I_1t-θ_t^0 to alleviate the exposition. All the results in this section assume that T_0→∞ and J→∞ (if J is present in the model).



 §.§ Difference-in-Differences
 
In Section <ref>, we have seen that the counterfactual mean proxies implied by the canonical difference-in-differences model are:
P_t^N=μ+J^-1∑_j=2^J+1Y^N_jt.
We consider the following estimator:
P̂^N_t=μ̂+1/J∑_j=2^J+1Y_jt, where μ̂=1/T∑_t=1^T(Y^N_1t-1/J∑_j=2^J+1Y^N_jt)=μ+ 1/T∑_t=1^T u_t.
Since P̂_t^N-P_t^N=μ̂-μ, Assumption <ref> holds for the simple difference-in-differences model provided that T^-1∑_t=1^T u_t = o_P(1), which is true under very weak conditions.







 §.§ Synthetic Control and Constrained Lasso
 


Several models in Section <ref> (including SC and constrained Lasso) imply a structure in which  the counterfactual proxy is a linear function of observed outcomes of untreated units. 

To provide a unified framework for these models, we use Y to denote a generic vector of outcomes and X to denote the design matrix throughout this section. For example, in Section <ref>, we set Y=Y^N_1 and  X=(Y^N_2,…,Y^N_J+1), where Y^N_j=(Y^N_j1,…,Y^N_jT)' ∈ℝ^T for 1≤ j≤ J+1. These models can be written as 

    Y=Xw+u,

where u=(u_1,…,u_T)'  ∈ℝ^T. Identification is achieved by requiring that X and u be uncorrelated (cf. Condition (SC)). 

Under the framework in (<ref>), different models correspond to different specifications for the weight vector w. For the SC model in Section <ref>, w is an unknown vector whose elements are nonnegative and sum up to one. More generally, one can simply restrict w to be any vector with bounded ℓ_1-norm. This is the constrained Lasso estimator. 


Since P_t^N is the t-th element of the vector Xw, the natural estimator is P̂_t^N being the t-th element of Xŵ, where ŵ is an estimator for w. The estimation of w depends on the specification. Let 𝒲 be the parameter space for w.  We consider the following version of the original SC estimator

    ŵ =  min_w    Y-Xw_2   s.t.   w  ∈𝒲 = { v ≥0,   v_1 =1}.

The constrained Lasso estimator is 

    ŵ= min_w Y-Xw_2   s.t.   w  ∈𝒲 =  { v:   v_1≤ K },

where K is bounded and K>0. In light of the estimator (<ref>), a natural choice is K=1. 

In general, we choose the parameter space 𝒲 to be an arbitrary subset of an ℓ_1-ball with bounded radius. The following result gives very mild conditions under which the constrained least squares estimators are consistent and satisfy Assumption <ref>.[To simplify the exposition, we do not include an  intercept in Lemma <ref>. Similar arguments could be used to prove an analogous result with an unconstrained intercept.] 


 Consider 
ŵ=min_v Y-Xv_2 s.t. v∈𝒲,
where 𝒲 is a subset of {v:v_1≤ K}
and K is bounded.  Assume w ∈𝒲, the data are β-mixing
with exponential speed, and other assumptions listed
at the beginning of the proof, including the identification condition (SC), then the estimator enjoys
the performance bounds  stated in the proof, in particular:
1/T∑_t=1^T (P̂_t^N-P_t^N )^2=o_P(1) and P̂_t^N-P_t^N=o_P(1), for any T_0+1≤ t≤ T.




Lemma <ref> provides several  features that are important for counterfactual inference in our setup. First, we allow J to be large relative to T. To be precise, we only require log J=o(T^c), where c>0 is a constant depending only on the β-mixing coefficients; see the Appendix for details. This is particularly relevant for settings in which the number of (potential) control units and the number of time periods have a similar order of magnitude as in our empirical application in Section <ref>. Second, Lemma <ref> does not rely on any sparsity assumptions on w, allowing for dense vectors. Third, compared to typical high-dimensional estimators (e.g., Lasso or Dantzig selector), our estimator does rely on tuning parameters that can be difficult to choose in times series settings. Finally, Lemma <ref> provides new theoretical consistency results for the canonical SC estimator in settings with time series data and potentially very many control units.






 §.§ Models with Factor Structures
 
The models for counterfactual proxies introduced in Section <ref> have factor structures. We provide estimation results for pure factor models (without regressors), factor models with regressors (interactive FE models), and matrix completion models. In this subsection, following standard notation, we let  N=J+1.


  §.§.§ Pure Factor Models




Recall from Section <ref> the standard factor model 
Y_jt^N=λ_j'F_t+u_jt,
where F=(F_1,…,F_T)'∈ℝ^T× k and Λ=(λ_1,…,λ_N)'∈ℝ^N× k represent the k-dimensional unobserved factors and their loadings, respectively. The counterfactual proxy for Y_1t^N is P_t^N=λ_1'F_t. We identify P_t^N by imposing the condition that the idiosyncratic terms  and the factor structure are uncorrelated (cf. Condition (FE)). 

We use the standard principal component analysis (PCA) for estimating P_t^N.[Note that PCA amounts to singular value decomposition, which can be computed using polynomial time algorithms, <cit.>. ] Let  Y^N∈ℝ^T× N
be the matrix whose (t,j) entry is Y_jt^N. We compute F̂=(F̂_1,…,F̂_T)'∈ℝ^T× k
to be  the matrix containing the eigenvectors corresponding to the
largest k eigenvalues of Y^N(Y^N)' with F̂'F̂/T=I_k.
Let λ̂_j'
denote the j-th row of  Λ̂=(Y^N)'F̂/T. Let F̂_t' denote the t-th row of F̂.  Our estimate for P_t^N is  P̂_t^N= λ̂_1'F̂_t. The following lemma guarantees the validity of this estimator in our context under mild regularity conditions. 


 
 Assume standard regularity conditions
given in <cit.>, including the identification condition (FE). Consider the factor model and the principal component estimator. Then, for any 1≤ t≤ T, as N →∞ and T →∞, we have P̂^N_t-P^N_t=O_P(1/√(N) + 1/√(T) )
and 1/T∑_t=1^T (P̂^N_t-P^N_t)^2=O_P(1/ N + 1/T).



The only requirement on the sample size is that both N and T need to be large. Similar to Theorem 3 of  <cit.>, we do not restrict the relationship between N and T. This is flexible enough for a wide range of applications in practice as   the number of units is allowed to be much larger than, much smaller than, or similar to the number of time periods. 




  §.§.§ Factor plus Regression Model: Interactive FE Model


Now we study the general form of panel models with interactive FEs. Following Section <ref>, these models take the form 
Y_jt^N=λ_j'F_t+X_jt'β+u_jt,
where X_jt∈ℝ^k_x are observed covariates and F=(F_1,…,F_T)'∈ℝ^T× k and Λ=(λ_1,…,λ_N)'∈ℝ^N× k represent the k-dimensional unobserved factors and their loadings, respectively. The counterfactual proxy for Y_1t^N is P_t^N=λ_1'F_t+X_1t'β. In this model, we identify the counterfactual proxy through the condition that the idiosyncratic terms are  independent of the factor structure and the observed covariates (cf. Condition (FE)). 

 The two most popular estimators are the common correlated effects (CCE) estimator by <cit.> and the iterative least squares estimator by <cit.>. We focus on the iterative least squares approach, but analogous results can be established for CCE estimators. 
 The notations for F_t,  λ_j, F̂_t and λ̂_j are the same as before. We compute 
 
    (F̂,Λ̂,β̂)=F,Λ,βmin∑_t=1^T ∑_j=1^N (Y^N_jt-X_jt'β-F_t'λ_j )^2   s.t.     F'F/T=I_k   Λ'Λ = Diagonal_k.


The estimate for P_t^N is P̂_t^N=λ̂_1'F̂_t+X_1t'β̂. The following result states the validity of applying this estimator in conjunction with our inference method.



 Assume the standard conditions
in <cit.>, including the identification condition (FE). 
Then, for any 1≤ t≤ T, 
P̂^N_t-P^N_t=O_P(1/√(T) + 1/√(N)) and 1/T∑_t=1^T(P̂^N_t-P^N_t)^2=O_P(1/T+ 1/N).



Under the conditions in Theorem 3 of <cit.>,  N is of the same order as T so that rate is really T^-1/2; 
however, the stated bound should hold more generally. 






  §.§.§ Matrix Completion via Nuclear Norm Regularization






Suppose that 

    Y_jt^N=M_jt+u_jt,   for 1≤ j≤ J+1  and 1≤ t≤ T,

where M_jt is the (j,t)-element of an unknown matrix M∈ℝ^(J+1)× T satisfying M_*≤ K, where ·_* denotes the nuclear norm (the sum of singular values). We observe Y_jt^N for (j,t)∈{1,…,T}×{1,…,J+1}\{(1,t):T_0+1≤ t≤ T}. The identifying condition is that E(u| M)=0 and that conditional on M, {u_j}_j=1^J+1 is independent across j, where u_j=(u_j1,…,u_jT)'∈ℝ^T.  The counterfactual proxy is P_t^N=M_1t for 1≤ t≤ T. 

The main challenge is to recover the entire matrix M despite the missing entries {Y_1t^N:T_0+1≤ t≤ T}. The literature on matrix completion considers the model  (<ref>) under the assumption of missingness at random and exploits the assumption that the rank of M is low.[See, for example, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>.] Recently, <cit.> introduce this method to study treatment effects in panel data models and point out  the unobserved counterfactuals correspond to  entries that are missing in a very special pattern, rather than at random. Assuming the usual low rank condition on M, they employ the nuclear norm penalized estimator and provide bounds on the estimation error in the typical setup of causal panel data models. 

We take a different approach here since our main goal is hypothesis testing instead of estimation. The key observation is that under the null hypothesis, there are no missing entries in the data. By imposing the null hypothesis, we  replace the missing entries with the hypothesized values and obtain a dataset that contains {Y_jt^N:1≤ j ≤ J+1, 1≤ t ≤ T }. The estimator for M we examine here is closely related to existing nuclear norm regularized estimators and is defined as

    M̂=    A∈ℝ^N× Tmin∑_t=1^T∑_j=1^N(Y_jt^N-A_jt)^2   s.t.   A_*≤ K,

where K>0 is the bound on the nuclear norm of the true matrix. In principle, it can be a sequence that tends to infinity. When M represents a factor structure with strong factors, K can be shown to grow at the rate √(NT). Clear guidance on how to choose K is still unavailable, but following <cit.> one can use  cross-validation.[The properties of cross-validation remain unknown in these settings.]  Alternatively one can use a pilot thresholded SVD estimator to get a sense of what K is, and use a somewhat larger value of K. The following result guarantees the validity of this estimator in our context under mild regularity conditions. 




Consider  the estimator M̂
defined in (<ref>). Assume that M_* ≤ K. Let the conditions listed at the beginning of the proof hold. Then, for any T_0+1≤ t≤ T,
P̂^N_t-P^N_t=o_P(1) and 1/T∑_t=K+1^T(P̂^N_t-P^N_t)^2=o_P(1).



The result is notable because no sub-Gaussian assumptions are required.  The estimator in (<ref>) does not explicitly require a low-rank condition on M. Instead, we impose a growth restriction on K. When M is generated by a strong factor structure and the null hypothesis contains full information on the missing entries, we can choose K≍√(NT) and our consistency result holds as long as N,T→∞ and   E(|u_jt|^2+c| M) is uniformly bounded for some c>0. In the case of weak factors, we can choose  K ≪√(NT) and obtain consistency. 




 §.§ Time Series and Fused Models

As pointed out in Section <ref>, time series models can be used to model counterfactual proxies with or without control units. We now discuss  low-level conditions under which fitting these models yields estimates good enough for the purpose of our conformal inference approach. 




  §.§.§ Autoregressive  Models

The linear autoregressive model with K lags can be written as
Y_1t^N=ρ_0+∑_j=1^Kρ_jY_1t-j^N+u_t,
where {u_t}_t=1^T is an iid sequence with E(u_t)=0.[Here
the model seems different, but  Section <ref>'s model implies this one with  ρ_0  = μ (1 - ∑_j=1^Kρ_j).] The counterfactual proxy for Y_1t^N is P_t^N= ρ_0+∑_j=1^Kρ_jY_1t-j^N. We write P_t^N as  P_t^N = y_t'ρ, where y_t=(1,Y_1t-1^N,Y_1t-2^N,…,Y_1t-K^N)'∈ℝ^K+1 
and  ρ=(ρ_0,…,ρ_K)'∈ℝ^K+1. The coefficient vector ρ can be estimated using least squares: ρ̂=(∑_t=K+1^Ty_ty_t')^-1(∑_t=K+1^Ty_tY_1t^N). The estimator for P_t^N is P̂_t^N = y_t 'ρ̂. 


Suppose that {u_t}_t=1^T is an iid sequence with E(u_1)=0
and E(u_1^4) uniformly bounded and the roots of   1-∑_j=1^Kρ_jL^j=0
are uniformly bounded away from the unit circle. Then, for any T_0+1≤ t≤ T,
P̂^N_t-P^N_t=o_P(1) and 1/T∑_t=K+1^T(P̂^N_t-P^N_t)^2=o_P(1).


As mentioned in Section <ref>, we can also apply  nonlinear autoregressive models
Y_1t^N=ρ(Y_1t-1^N,Y_1t-2^N,…,Y_1t-K^N)+u_t,
where ρ is a nonlinear function, in which case the counterfactual proxy is P_t^N= ρ(Y_1t-1^N,Y_1t-2^N,…,Y_1t-K^N).

Let ρ̂ be an estimator for ρ and P̂_t^N= ρ̂(Y_1t-1^N,Y_1t-2^N,…,Y_1t-K^N). This estimator can be parametric, semiparametric, or fully nonparametric and is only required to be consistent.



  Suppose  that (1) ρ̂-ρ= O_P(r_T) with r_T =o(1) for some appropriate norm · and  max_K+1 ≤ t ≤ T    |ρ̂(Y_1t-1^N,Y_1t-2^N,…,Y_1t-K^N) -ρ(Y_1t-1^N,Y_1t-2^N,…,Y_1t-K^N)| ≤ℓ_T ρ̂- ρ  for some ℓ_T r_T =o(1). Then, for any T_0+1≤ t≤ T,
P̂^N_t-P^N_t=o_P(1) and 1/T∑_t=K+1^T(P̂^N_t-P^N_t)^2 =o_P(1).



The primitive regularity conditions and the definitions of  the neural network estimators possessing these properties can be found, for example, in <cit.> and <cit.>.





  §.§.§ Fused Panel/Time Series Models with AR Errors




Here we provide generic conditions for the fused panel/time series models described in Section <ref>.  
In particular, AR models can be used to filter the estimated residuals and obtain
near iid errors. In Equation (<ref>) of Section <ref>, we introduce an autoregressive structure in the error terms: 
Y_1t^N=C_t^N+ε_t and ε_t=ρ(ε_t-1)+u_t,
where C_t^N can be specified as a panel data model discussed before. Due to the autoregressive structure in ε_t, the counterfactual proxy is P_t^N=C_t^N+ρ(ε_t-1). 

We estimate P_t^N via a two-stage procedure. In the first stage, we estimate C_t^N using the techniques we considered before and obtain say Ĉ_t^N. In the second stage, we estimate ρ(ε_t-1) by fitting an autoregressive model to the estimated residuals {ε̂_t}_t=1^T, where ε̂_t=Y_1t^N-Ĉ_t^N. For simplicity, we consider a linear model in the second stage estimation. Analogous results can be obtained for more general models. To be specific, assume that  ε_t=x_t'ρ+u_t, where x_t=(ε_t-1,ε_t-2,…,ε_t-K)'∈ℝ^K
and ρ=(ρ_1,ρ_2,…,ρ_K)'∈ℝ^K. 

 Given {ε̂_t }_t=1^T from the first-stage estimation, we define x̂_t=(ε̂_t-1,ε̂_t-2,… ,ε̂_t-K)'∈ℝ^K
and ρ̂=(∑_t=K+1^Tx̂_tx̂_t')^-1(∑_t=K+1^Tx̂_tε̂_t). To compute the p-value, we use  {û_t}_t=K+1^T with  û_t=ε̂_t-x̂_t'ρ̂ in the permutation. By the following result, this procedure is valid under very mild conditions for the first-stage estimation. 



Suppose that {u_t}_t=1^T
is an iid sequence with E(u_t)=0 and E(u_1^4) uniformly
bounded and the roots of 1-∑_j=1^Kρ_jL^j=0 are uniformly
bounded away from the unit circle. We assume that 
(1) ∑_t=1^T(Ĉ^N_t-C^N_t)^2=o_P(T), 
and (2) Ĉ^N_t-C^N_t=o_P(1) for T_0-K+1≤ t≤ T. Then, for any T_0+1≤ t≤ T,
P̂^N_t-P^N_t=o_P(1) and ∑_t=K+1^T(P̂^N_t-P^N_t)^2=o_P(T)



Note that the conditions in Lemma <ref> for the autoregressive part are the same as in Lemma <ref>. Consistency of Ĉ_t^N can be verified using existing results, for example, those in Sections <ref>–<ref>. 





§ EMPIRICAL APPLICATION



We revisit the analysis in <cit.> who study the impact of decriminalizing indoor prostitution. 
They consider the case of Rhode Island, where a judge unanticipatedly decriminalized indoor sex work in July 2003 such that, until the recriminalization in November 2009, Rhode Island had decriminalized indoor and prohibited street prostitution. 

We focus on the effect of legalizing indoor prostitution on female gonorrhea incidence. Our outcome of interest is log female gonorrhea incidence per 100,000. We use the data on gonorrhea cases from the Center for Disease Control (CDC)'s Gonorrhea Surveillance Program previously analyzed by <cit.>; see their Section 3 for a detailed description and descriptive statistics. The female gonorrhea series date back to 1985 such that T_0=19 and T_∗=6. Figure <ref> displays the raw data for Rhode Island and the rest of the U.S. states.



[Figure <ref> around here.]



 We apply three different CSC methods: difference-in-differences, canonical SC, and constrained Lasso with K=1. Recall that constrained Lasso  nests both difference-in-differences and SC. Following <cit.>, the set of potential control units includes all other U.S. states and the District of Columbia (J=50). We choose S_1 as our test statistic and report p-values computed based on moving block and iid permutations.[To keep computation tractable, we randomly sample 10,000 iid permutations with replacement.] All computations were performed in  <cit.>.


Before turning to the main results, we use the placebo tests proposed in the Appendix to assess the plausibility of the underlying assumptions. Specifically, based on the pre-treatment data, we test H_0: θ_2003-τ+1=…=θ_2003=0 for τ∈{1,2,3}. Rejections of this null undermine the credibility of the assumptions underlying our procedure and the inferences on policy effects in the post-treatment period. Table <ref> presents the results. Figure <ref> complements the formal tests with plots of the residuals from fitting the three models to the pre-treatment data. The placebo tests and the residual plots provide evidence in favor of the credibility of our inference method in conjunction with SC and, especially, constrained Lasso, but suggest that the difference-in-differences results need to be interpreted with caution.



[Table <ref> around here.]





[Figure <ref> around here.]




Table <ref> reports p-values from testing the null hypothesis of a zero effect:

    H_0: θ_2004=θ_2005=… =θ_2009=0 .

The null hypothesis (<ref>) is rejected at the 10% level based on both permutation schemes and all three methods. 



[Table <ref> around here.]



Figure <ref> displays pointwise 90% confidence intervals. The results are similar for all three methods. While the effect was not or only marginally significant during the first three years, legalizing indoor prostitution significantly decreased the incidence of female gonorrhea thereafter, corroborating the findings by <cit.>.



[Figure <ref> around here.]



To investigate the robustness of our results, we perform a leave-one-out robustness check <cit.> to assess whether our findings are driven by a single control state. We iteratively exclude from the control group one of the states for which either the SC or constrained Lasso weights estimated based on the pre-treatment data are non-zero and compute the p-values for testing hypothesis (<ref>). Figure <ref> displays the distribution of the resulting p-values. Overall, our results are robust and not driven by a single control state: except for one specification, all results are significant at the 10%-level.



[Figure <ref> around here.]




1


§ ACKNOWLEDGEMENTS

We are grateful to Guido Imbens, Jacopo Diquigiovanni, Bruno Ferman,  the Co-Editor (Matias Cattaneo), anonymous referees, and many seminar and conference participants for valuable comments. We would like to thank Scott Cunningham and Manisha Shah for sharing the data for the empirical application. Wüthrich is also affiliated with CESifo and the Ifo Institute. Victor Chernozhukov gratefully acknowledges funding by the National Science Foundation. All errors are our own.
 
 

0

 


1.0

apalike










§ FIGURES MAIN TEXT


























§ TABLES MAIN TEXT

















  Online Supplemental Appendix to “An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls”




Victor Chernozhukov, Massachusetts Institute of Technology. Email: vchern@mit.edu

Kaspar Wüthrich, University of California San Diego. Email: kwuthrich@ucsd.edu

 Yinchu Zhu, Brandeis University. Email: yinchuzhu@brandeis.edu

 

 
 
 


equationsection
figuresection
tablesection

thmsection
lemsection

[sections]
[sections]l1







 
 
 



§ EXTENSIONS





 §.§ Testing Hypotheses about Average Effects over Time



In addition to testing sharp null hypotheses, researchers are often also interested in testing hypotheses about average effects over time, θ̅=T_∗^-1∑_t=T_0+1^Tθ_t:

    H_0: θ̅=θ̅^0


Hypothesis (<ref>) can be tested by collapsing the data into averages of non-overlapping blocks of T_∗ observations over the time dimension. To simplify the exposition, we assume that T/T_∗ is an integer. Note that Assumption <ref> implies the following model for the average potential outcomes Y̅_1r^N=T_∗^-1∑_t=r^r+T_∗-1Y^N_t and Y̅_1r^I=T_∗^-1∑_t=r^r+T_∗-1Y^I_t: 

    [        Y̅_1r^N = P̅^N_r+u̅_r; Y̅_1r^I = P̅^N_r+ θ̅_r +u̅_r;                              ] |    E( u̅_r) = 0,   r=1,T_∗+1,2T_∗+1,…,T_0-T_∗+1,T_0+1,

where P̅_r^N=T_∗^-1∑_t=r^r+T_∗-1P^N_t and u̅_r=T_∗^-1∑_t=r^r+T_∗-1u_t. Define the aggregated (collapsed) data under the null as =(Z̅_1,…,Z̅_T_0+1)', where

    Z̅_r=(Y̅_1r^N,Y̅_2r^N,…,Y̅^N_J+1r,X̅'_1r,…,X̅'_J+1r)',     r<T_0+1
    (Y̅_1r^I-θ̅^o,Y̅_2r^N,…,Y̅^N_J+1r,X̅'_1r,…,X̅'_J+1r)',     r=T_0+1

and X̅_jr=T_∗^-1∑_t=r^r+T_∗-1X_jt for j=1,…,J+1. Note that testing hypothesis (<ref>) is equivalent to testing a hypothesis concerning a per-period effect based on the aggregated data . Specifically, we estimate the average proxy P̂̅̂_r^N based on the aggregated data  and obtain the residuals û̅̂=(û̅̂_1,û̅̂_T_∗+1,…, û̅̂_T_0+1)', 
where û̅̂_r=Y̅_1r^N- P̂̅̂^N_r. The test statistic is S( û̅̂), and p-values can be obtained based on permutations of û̅̂ as described in Section <ref>. 

The key assumption underlying this procedure is that the average mean proxy P̅_r^N can be identified and consistently estimated based on the aggregated data . This is the case for SC and the other regression-based estimators discussed in Sections <ref>–<ref>, provided that E(u̅_rY̅_jr^N)=0 for 2≤ j≤ J+1 and that the sufficient conditions for consistent estimation in Section <ref> hold for the aggregate data. By contrast, identification and estimation of P̅_r^N may not be possible for nonlinear and dynamic models. Under consistent estimation of P̂̅̂^N_r, the formal properties of the test follow from the results in Section <ref> since stationarity and weak dependence of {u_t} imply stationarity and weak dependence of {u̅_r}. Alternatively, if P̂̅̂^N_r can be shown to be stable and the aggregate data are stationary and weakly dependent, the properties of the test follow from the results in Section <ref>. Finally, we emphasize that the effective sample size is T/T_∗ instead of T, such that T needs to be substantially larger than T_∗.







 §.§ Multiple Treated Units

Our method can be extended to accommodate multiple treated units by collapsing the data into averages across the treated units. Consider a setup with L treated units, indexed by j=1,…,L, and J control units, indexed by j=L+1,…,J+L. Suppose that Assumption <ref> holds for all treated units:

    [        Y_jt^N = P_jt^N + u_jt; Y_jt^I = P_jt^N + θ_jt + u_jt;                               ] |    E( u_jt) = 0,     t=1,…,T ,    j=1,…,L.

Under this assumption, hypotheses about the unit-specific policy effects {θ_jt} can be tested by separately applying the proposed inference procedure to each treated unit. In addition, one is often also interested in conducting inferences about the average treatment effects on the treated units, {θ̅_t}, where θ̅_t=L^-1∑_j=1^Lθ_jt.

Specifically, consider the following null hypothesis:

    H_0: (θ̅_T_0+1,…, θ̅_T)=(θ̅_T_0+1^0,…, θ̅_T^0 ).

To test hypothesis (<ref>), note that if Assumption <ref> holds for all treated units, we have the following model for the average potential outcomes Y̅^N_t=L^-1∑_j=1^LY_jt^N and Y̅^I_t=L^-1∑_j=1^LY_jt^I:

    [        Y̅_t^N = P̅_t^N + u̅_t; Y̅_t^I = P̅_t^N + θ̅_t + u̅_t;                               ] |    E( u̅_t) = 0,     t=1,…,T,

where P̅_t^N=L^-1∑_j=1^LP_jt^N and u̅_t=L^-1∑_j=1^Lu_jt. Define the data under the null as =(Z̅_1,…,Z̅_T)', where

    Z̅_t=(Y̅^N_t,Y^N_L+1t,…,Y^N_J+Lt,X̅'_t,X'_L+1t,…,X'_J+Lt)',     t≤ T_0,
    (Y̅_t^I-θ̅_t^o,Y^N_L+1t,…,Y^N_J+Lt,X̅'_t,X'_L+1t,…,X'_J+Lt)',     t>T_0,

and X̅_t=L^-1∑_j=1^LX_jt. To test hypothesis (<ref>), we compute the estimated average proxy P̂̅̂_t^N based on the aggregated data  and obtain the residuals û̅̂=(û̅̂_1,…, û̅̂_T)', where û̅̂_t=Y̅_t^N- P̂̅̂^N_t for  t=1,…,T. The test statistic is S( û̅̂), and p-values can be obtained based on permutations of û̅̂ as described in Section <ref>.
 The formal properties of this test follow from the results in Section <ref>. 
 


 §.§ Placebo Tests



Here we propose easy-to-implement placebo tests for assessing the credibility of inferences based on our method. We recommend applying these placebo tests when using our inference procedures.

Following <cit.>, the idea is to consider a placebo intervention before the actual intervention took place. For a given τ≥ 1, we use our method to test the null hypothesis 

    H_0: θ_T_0-τ+1=…=θ_T_0=0

based in the pre-treatment data =( Z_1,…,Z_T_0)'. Using an appropriate CSC method, we compute the counterfactual mean proxies P̂_t^N based on  and obtain the residuals 

    û=(û_1,…,û_T_0)',   û_t=Y_1t^N-P̂^N_t,    t=1,…, T_0.

We then apply the proposed inference method, treating {1,…,T_0-τ} as the pre-treatment period and {T_0-τ+1,…,T_0} as the post-treatment period. The theoretical properties of such placebo tests follow directly from the results in Section <ref>. As illustrated in our application, it may be useful to complement the formal testing results with plots of the pre-treatment residuals (û_1,…,û_T_0).

A rejection of the null hypothesis (<ref>) undermines the credibility of the assumptions underlying our procedure and the inferences on the policy effects in the post-treatment period. While non-rejections provide evidence in favor of our method, it is important to emphasize that such non-rejections do not “prove” that our method is valid. Moreover, by construction, the placebo tests cannot be used to assess some of the key assumptions underlying our approach, such as the invariance of the distribution of {u_t} under the intervention.

Given our inference method's genericness, the placebo tests discussed here can be used to assess and compare the credibility of different CSC methods. For example, in our empirical application, the placebo tests provide evidence in favor of SC and constrained Lasso, but suggest that the difference-in-differences results need to be interpreted with caution.



 

§ INTERPRETATION AS A STRUCTURAL BREAKS TEST

 
A key assumption underlying our method is the invariance of the distribution of {u_t} under the intervention (Assumption <ref>). In settings where this assumption fails, following the literature on end-of-sample structural breaks tests <cit.>, our procedure can be used as a test of the null hypothesis that the policy has no impact whatsoever against the alternative hypothesis that θ 0 and or the policy affects the distribution of {u_t}. 

Consider the following testing problem:

    H_0:
      Y_1t^N=P_t^N+u_t,   E(u_t)=0,  t=1,…,T,   θ=0,   and
    {u_t}_t=0^T  is stationary and weakly dependent.

against

    H_1:
      Y_1t^N=P_t^N+u_t,  t=1,…,T,  θ 0  and or the distribution of
    {u_t}_t=T_0+1^T differs from that of {u_t}_t=s^s+T_∗-1 for  s=1,…, T_0-T_∗+1.

Hypothesis (<ref>) can be tested by applying our procedure to the data under the null where Y_1t^N=Y^I_1t=Y_1t for all t=1,…,T. The theoretical size properties of this test follow directly from the results in Section <ref>. This test has power against location shifts induced by θ 0 as well as changes in the distribution of {u_t} that increase the quantiles of S(u) (e.g., scale shifts);  see <cit.> for a related discussion.






§ PREDICTION SETS FOR RANDOM POLICY EFFECTS



In the main text, we assume that the policy effect sequence {θ_t} is fixed (Assumption <ref>). Here we show that our procedure generates valid prediction sets for θ_t when θ_t is assumed to be random, as, for example, in <cit.>.[Following the literature on conformal prediction <cit.>, we use the terminology “prediction set” instead of  “confidence set” because the resulting sets are not conventional confidence sets; see <cit.> for a discussion of the differences between confidence intervals and prediction intervals for SC.] Our analysis here is in the same spirit as the classical conformal prediction literature, which aims at constructing prediction intervals for future values of a (random) target quantity of interest. Specifically, we will show that when θ_t is random, Algorithm <ref> provides unconditionally valid (1-α)-prediction sets with non-asymptotic performance guarantees. 

To state the result, define  ^*=(Z_1^*,…,Z_T^*)' with Z_t^*=(Y^N_1t,Y^N_2t,…,Y^N_J+1t,X'_1t,…,X'_J+1t)' for 1≤ t≤ T_0+T_*. Notice that ^* contains the true counterfactuals.


 
Assume that T_* is fixed. Suppose that Y_1t^N=P_t^N+u_t, 1≤ t≤ T, where {u_t} is a centered and stationary stochastic process, and that the policy effects θ_t:=Y_1t^I-Y_1t^N are random. Suppose that Assumption <ref> holds for P̂_N computed using ^*. Impose  Assumption <ref>.<ref> if Π=Π_all.  Impose Assumption <ref>.<ref>,
if Π=Π_→. Assume the statistic S(u) has a density function bounded by D. Then

    |P(θ_t∈𝒞_1-α(t))- (1-α) | ≤ C ( δ̃_T + δ_T + √(δ_T) + γ_T),

where 𝒞_1-α(t) is defined in Algorithm <ref> and δ̃_T = (T_*/T_0)^1/4(log T). The constant C depends on T_*, M and D,  but not on T.


Let us briefly discuss the interpretation of our model when  θ_t is regarded as random. Suppose that the potential outcomes are Y_1t^I=P^I_t+u_t^I and Y_1t^N=P^N_t+u_t^N. The counterfactual mean proxies and the prediction errors under the policy, {P_t^I} and {u_t^I}, may differ from the counterfactual mean proxies and the errors in the absence of the policy, {P_t^N} and {u_t^N}. Algorithm <ref> provides prediction intervals for 

    θ_t:=Y_1t^I-Y_1t^N=Δ_t^P+Δ_t^u,

where Δ_t^P:=P^I_t-P_t^N and Δ_t^u:=u^I_t-u_t^N. Here  θ_t captures both the effects of the policy on the mean proxy, Δ_t^P, and the effect on distribution of the error, Δ_t^u (e.g., a scale shift).  Therefore, we are only assuming that in the absence of the policy, {u_t^N} is a stationary process. All the changes to the mean, the variance, or other features of the distribution due to the policy are captured by the policy effects {θ_t}.

In sum, with random policy effects, our procedure simply provides a prediction set for the policy effects defined as the difference between the two potential outcomes.




§ MODEL-FREE EXACT VALIDITY UNDER EXCHANGEABILITY


Every permutation procedure that is approximately valid in time series settings should have good properties in “ideal” settings where the data are iid or exchangeable. The following theorem, which is based on standard arguments, shows that under exchangeability of the data, our conformal inference approach achieves exact finite sample size control. This result is model-free in the sense that we do not need to use a correct or consistent estimator for the counterfactual mean proxy. As a result, our procedure controls size under arbitrary forms of misspecification and is fully robust against overfitting.


 Let Π be Π_→ or Π_all. Suppose that the data {Z_t}_t=1^T are iid or exchangeable with respect to Π under the null hypothesis and that û_t =g (Z_t, β̂), where the estimator β̂= β̂({Z_t}_t=1^T) is invariant with respect to any permutation of the data. Then, under the null hypothesis, {û_t}_t=1^T is an exchangeable sequence and the permutation p-value is unbiased in level:

    P(p̂≤α)≤α.

Moreover, if {S(û_π) }_π∈Π has a continuous distribution,

    α-1/|Π|≤ P(p̂≤α).



Theorem <ref> requires that the estimators are invariant under permutations of the data under the null hypothesis. Invariance holds for regression-based estimators such as SC, constrained Lasso, or penalized regression, provided that the null hypothesis is imposed for estimation but may fail for dynamic models such as linear and non-linear autoregressive models.


An inspection of the proof of Theorem <ref> shows that our procedure achieves finite sample size control whenever the residuals {û_t}_t=1^T are exchangeable. We demonstrate that exchangeability of  {û_t}_t=1^T is implied if the data {Z_t}_t=1^T are iid or exchangeable. Exchangeability {û_t}_t=1^T may hold even if the data are not exchangeable.  For example, in the difference-in-difference model, the outcome data can have an arbitrary common trend eliminated by differencing, making it possible for the residuals to be iid or exchangeable with non-iid data.

The finite sample validity under exchangeability is crucial for the robustness of our proposal. While high-dimensional CSC approaches may overfit in small samples, Theorem <ref> shows that the proposed method does not suffer from overfitting or misspecification. The fundamental reason is that we exploit symmetry rather than to completely rely on the consistency of the estimator.  Imposing the null hypothesis when estimating the model leads to invariance of the estimator under permutations, exchangeable residuals, and finite sample validity. As a result, the proposed procedure is more robust than alternative approaches that estimate the proxies based on the pre-treatment data without imposing the null hypothesis.

Even in time series settings where exchangeability fails, imposing the null hypothesis for estimation is crucial for achieving a good performance in typical CSC applications where T_0 is rather small <cit.>. Figure <ref> augments Figure <ref> with results for T=100  (T_0=99,T_∗=1). In the empirically relevant case where T_0=19, estimating P_t^N under the null yields an excellent performance. Irrespective of the degree of persistence, size accuracy is substantially better than when P_t^N is estimated based on pre-treatment data only. Even when T_0=99, which is much larger than the T_0 in many CSC applications, imposing the null yields notable performance improvements. In fact, the size accuracy is better for T_0=19 when P_t^N is estimated under the null than for T_0=99 when P_t^N is estimated based on the pre-treatment data only.







§ SUFFICIENT CONDITIONS FOR ESTIMATOR STABILITY



In this section, we provide sufficient conditions for the estimator stability Assumption <ref>. We first present generic sufficient conditions for low-dimensional models. For high-dimensional models, the theoretical analysis is more difficult and a case-by-case analysis is needed. We are not aware of any theoretical work that establishes Assumption <ref> for any high-dimensional model. Here we verify the stability condition for constrained Lasso; stability of Ridge regression is verified in Appendix <ref>.





 §.§ Generic Sufficient Condition for Low-dimensional Models


Consider β̂()=min_β∈L̂(;β),
where L̂(;β) is a loss function and ⊂ℝ^p
for a fixed p. Let  be a set of subsets of {1,…,T}. Notice that Assumption <ref> only requires  to be a singleton, but in this subsection and the next, we allow  to be a class of subsets. 



	Suppose that the following
	conditions hold:
  

  * sup_β∈|L̂(;β)-L(β)|=o_P(1)
	for some non-random L(·). 

  * max_H∈sup_β∈|L̂(_H;β)-L(β)|=o_P(1).

  * L(·) is continuous at β_*, min_βL(β)
	has a unique minimum at β_* and  is compact.
	
	Then max_H∈()-(_H)_2=o_P(1). 



In the literature of misspecified models, β_* is usually
referred to as the pseudo-true value <cit.>.
In M-estimation with L̂(;β)=T^-1∑_t=1^Tl(Z_t;β),
one can often show sup_β|L̂(;β)-L(β)|=o_P(1)
with L(β)=El(Z_1;β); in GMM models with L̂(;β)=T^-1∑_t=1^Tψ(Z_t;β)_2,
one can often use L(β)=Eψ(Z_1;β)_2. 

The proof of Lemma <ref> shows that
()-β_*_2=o_P(1) and max_H∈(_H)-β_*_2=o_P(1).
In other words, the stability of the estimator arises from the consistency
to the pseudo-true value β_*. Such consistency holds under
very weak conditions. We essentially only require a uniform law of
large numbers. This can be verified for many low-dimensional models
under weakly dependent data. The conclusion of Lemma  <ref> translates to  Assumption <ref> once we derive a bound on  max_π∈Πsup_β_1≠β_2|S(^π;β_1)-S(^π;β_2)|/β_1-β_2_2; this requires knowledge of  the model structure. 

For example, suppose that max_π∈Πsup_β_1≠β_2|S(^π;β_1)-S(^π;β_2)|/β_1-β_2_2 =O_P(1).[We only need S(;β) to be Lipschitz with respect to β. In the simple example of S(;β)=|Y_T_0+1-X_T_0+1'β|, this only requires  X_T_0+1_2 to be bounded.] Then max_π∈Π|S(^π;())-S(^π;(_H))|=o_P(1). Here is how we apply Theorem <ref> to obtain the asymptotic size control. Fix an arbitrary δ>0. We can simply choose the constant function ϱ_T(x)=δ for Assumption <ref>. Since max_π∈Π|S(^π;())-S(^π;(_H))|=o_P(1),  Assumption <ref> holds for some γ_1,T=o(1) (due to the definition of convergence in probability).
Assume that Ψ(x,β) has bounded derivative with respect to x. Then we can choose ξ_T=K_1 for a large constant K_1 and γ_2,T=0. As a result, Theorem <ref> has 

    |P(p̂≤α)-α|
       ≤ C_1√(ξ_Tϱ_T(T_0/R+2k))+C_1(T_0^-1R[log(T_0/R)]^1/D_3)^1/4+C_1exp(-(k-T_*+1)^1/D_3)
           +C_1√(γ_1,T)+C_1√(γ_2,T)
        =C_1√(K_1 δ)+C_1(T_0^-1R[log(T_0/R)]^1/D_3)^1/4+C_1exp(-(k-T_*+1)^1/D_3)+C_1√(o(1)).

Since R≍ T_0/log(T_0), we have T_0^-1R[log(T_0/R)]^1/D_3=o(1). Choosing k≍log(T_0) and assuming that T_* is fixed, we obtain exp(-(k-T_*+1)^1/D_3)=o(1). Therefore, the above display implies 

    |P(p̂≤α)-α|=C_1√(K_1 δ)+o(1).

Since δ>0 is arbitrary, we have |P(p̂≤α)-α|=o_P(1).



 §.§ Constrained Lasso


Here we propose sufficient conditions for estimator stability for constrained Lasso. In contrast to Sections <ref> and <ref>, we do not impose correct specification but study the behavior of the constrained Lasso estimator under potential misspecification. To make this explicit, we use β instead of w to denote the coefficient vector in this subsection.  Here, it is possible that EX_t(Y_t-X_t'β)≠ 0 for any β∈. In practice, this arises when  the relationship between X_t and Y_t is non-linear or when the constraint set  is too small. For example, the true parameter could be non-sparse with exploding ℓ_1-norm (e.g., β=(1,…,1)'/√(J)). 






We first introduce some additional notation. Define Y_t=Y_1t^N and X_t=(Y^N_2t,…,Y^N_J+1t)' and let {(Ỹ_t,X̃_t)}_t=1^T be iid from the
distribution of (Y_1,X_1) and independent of the data {(Y_t,X_t)}_t=1^T. 

The constrained Lasso objective functions based on the data under the original data and after switching out observations with t∈ H are given by

    (β)=1/T∑_t=1^T(Y_t-X_t'β)^2  and  _H(β)=T^-1∑_t=1^T(Y_t,H-X_t,H'β)^2,

where (Y_t,H,X_t,H)=(Y_t,X_t) for t∉ H and (Y_t,H,X_t,H)=(Ỹ_t,X̃_t)
for t∈ H. The corresponding constrained Lasso estimators are

    ()=min_β∈(β)  and  (_H)=min_β∈_H(β),

where ⊆{v∈^J: v_1≤ K} and K>0 is a constant. Furthermore, we define =T^-1∑_t=1^TX_tX_t' and  =T^-1∑_t=1^TX_tY_t. Similarly, for H⊂{1,…,T}, let _H=T^-1∑_t=1^TX_t,HX_t,H' and
_H=T^-1∑_t=1^TX_t,HY_t,H. Finally, let  be a set of subsets of {1,…,T}. 



	Suppose that the following conditions
	hold:
	 
	
	
  * with probability at least 1-γ_1,T, _H-_∞≤ c_T
	and _H-_∞≤ c_T for all H∈.
	
	
  * with probability at least 1-γ_2,T, min_v_0≤ sv' v/v_2^2≥κ_1.
	
	
  * with probability at least 1-γ_3,T, max_H∈(_H)_0≤ s/2
	and ()_0≤ s/2. 
	
  * P(max_1≤ t≤ TX_t_∞≤κ_2)=1.
	
	Let _t=Y_t-X_t'() and _t,H=Y_t-X_t'(_H). Then we have that 
	
    P(max_H∈max_1≤ t≤ T|_t-_t,H|≤2κ_2√(κ_1sc_TK(2K+1)))≥1-γ_1,T-γ_2,T-γ_3,T.



Lemma <ref> provides sufficient conditions for  perturbation stability. Inspecting the proof, we notice that the argument does not require the estimator to converge to anything. To our knowledge, this is the first result of this kind. In the conformal prediction literature, one-observation perturbation stability has been considered in Assumption A3 of <cit.>, who only verify it assuming correct model specification and consistent variable selection. There is also a strand of literature in statistics that considers misspecified models in high dimensions and focuses on the pseudo-true value. For example, for linear models, the pseudo-true value represents the best linear projection and is often assumed to be sparse, making it possible to establish consistency of Lasso to this pseudo-true value <cit.>. We do not make these assumptions. Lemma <ref>  allows the model to be misspecified and the pseudo-true value may or may not be consistently estimated by constrained Lasso.


Lemma <ref> says that when the solution of
constrained Lasso is sparse, the stability of  and 
guarantees the stability of the estimator. When |H|≍log T_0
and the observed variables are bounded, we can choose c_T≍ T_0^-1log(T_0).
The sparse eigenvalue condition can typically be verified whenever
s≤ cT, where c>0 is a constant that depends on the eigenvalues
of E. Thus, Lemma <ref> would guarantee
that when sup_H∈|H|≲log T_0, we have 

    max_H∈max_1≤ t≤ T|_t-_t,H|=O_P(√(sT_0^-1log T_0)).


Therefore, whenever the solutions () and (_H) are
sparse enough with s=o(T_0/log(T_0)), we can expect stability
of the estimated residuals. One implication is that since ()_0
and (_H)_0 are clearly bounded above by J, the stability
should easily hold for J≪ T_0/log(T_0).


We now give an explicit formula for ϱ_T(·) in Assumption <ref>. Suppose that S(,β)=|T_*^-1∑_t=T_0+1^T_0+T_*(Y_t-X_t'β)|. Then the above display implies that max_H |S(;())-S(;(_H))|=O_P(√(sT_0^-1log T_0)). As in Appendix <ref>, we can use Theorem <ref> to obtain the asymptotic size control. 
Assume s=o(T_0/log(T_0)). Let q_T=(sT_0^-1log(T_0))^-1/4. Notice that q_T→∞. Since  q_T→∞ and |H|≲log(T_0), the above result of  max_H |S(;())-S(;(_H))|=O_P(√(sT_0^-1log T_0)) implies that there exists γ_1,T=o(1) such that  Assumption <ref> holds with  ϱ_T(x)=q_T√(sT_0^-1x). Then by the same argument as in Appendix <ref>, Theorem <ref> implies

    |P(p̂≤α)-α|
    ≤ C_1√(K_1 ϱ_T(T_0/R+2k))+o(1),

where K_1 is the same constant as in Appendix <ref>. Since k≤ T_0/R, R≍ T_0/log(T_0) and q_T →∞, we have

    ϱ_T(T_0/R+2k)≤ϱ_T(3T_0/R) = √(3) q_T √(sT_0^-1(T_0/R))≲ q_T √(sT_0^-1log(T_0)) = 1/q_T → 0.




Stability does not imply that the constrained Lasso residuals, _t, are close to ε_*,t=Y_t-X_t'β_*, where β_*=min_β∈E(Y_t-X_t'β)^2 is a pseudo-true value. In <cit.>, we show that -β_* _2=O_P((T_0^-1log J)^1/4). However, this is far from enough to conclude that |X_t'(-β_*)|=o_P(1) due to the high-dimensionality of X_t. The usual Cauchy-Schwarz bound X_t_2 -β_* _2 would not converge to zero; the Hölder bound X_t_∞-β_* _1 does not suffices either since  -β_* _1 does not converge to zero. [When  is a bounded ℓ_1-ball, it is in fact impossible to achieve  -β_* _1=o_P(1) <cit.>.] Even if  β_* is assumed to be sparse, one would still require β_*_0=o(√(T_0)/log J). 
In contrast, our stability condition discussed above only requires the much weaker condition of β_*_0=o(T_0/log T_0). Therefore, stability condition could be satisfied even if the estimated residual does not converge to a pseudo-true target. 








§ CONSISTENCY AND ESTIMATOR STABILITY


In this section, we illustrate the difference between stability and consistency using a simple and analytically tractable example: Ridge regression. We shall show that under correct specification, Ridge may not be consistent, while still satisfying estimator stability. 

To keep theoretical analysis tractable, we work under stylized conditions. Given data Z=(Y,X), we define the ridge estimator

    _λ(Z)=(X'X+λ I_T)^-1X'Y,

where λ is a tuning parameter. Then we can compute the residuals (Z;λ)=Y-X _λ(Z).


	Suppose that Y=Xβ+u, where  X=(X_1,…,X_T)'∈^T× J and E(u| X)=0. Assume that β_2 is bounded away from zero and infinity, E(uu'| X)=σ^2 I_T for  a constant σ>0, J≲  T^κ_0 with κ_0∈(0,2/3) and for some constants κ_1,κ_2>0, 
		
    P( κ_1 T ≤λ_min(X'X)≤λ_max(X'X)≤κ_2 T )≥1-o(1)

	and 
		
    P(  max_1≤ t≤ TX_t≤κ_3 √(J)   and  X'u_2≤κ_4 √(JT))≥1-o(1).

	Let Z_H=(Ỹ,X̃) be the perturbed data with |H|≍log T. For ũ=Ỹ-X̃β, assume that
		
    P( X̃'ũ-X'u_2≤κ_5 √(J|H|))≥1-o(1)

	and 
		
    P(  X̃'X̃-X'X≤κ_6 (|H|+J) )≥1-o(1),

	where κ_5,κ_6>0 are constants. 
	  	If λ≍ T, then E X (_λ(Z)-β)_2^2≳ T and (Z;λ)-(Z_H;λ) _∞=o_P(1)[Notice that consistency is enough to derive asymptotic size control from Theorem <ref>; see the discussion at the end of Appendix <ref> for an explicit formula for quantities in Assumption <ref>. ]. 


Lemma <ref> provides a robustness guarantee for the validity of the procedure. Under the ideal choice of the tuning parameter λ, we would expect consistency of _λ(Z) and hence validity of the procedure.   However, Lemma <ref> states that even when the tuning parameter is badly chosen such that consistency fails, one might still expect the estimator to be stable under perturbations, which is sufficient for the validity of our inference procedure. This is important in practice since computing optimal tuning parameters is often difficult. 

Now we give very simple sufficient (but probably far from necessary) conditions for the assumptions of Lemma <ref>. Suppose that the data is independent across t and rows of X are sub-Gaussian. Then standard results in random matrix theory can be used to verify these assumptions. For example, X'X/T-E(X'X/T) has eigenvalues tending to zero (e.g., Theorem 5.39 and Remark 5.40 in <cit.>). Hence, as long as E(X'X/T) has bounded eigenvalues,  there exist constants κ_1,κ_2 that satisfy the condition in Lemma <ref>. Let us further assume that entries of X are also bounded and entries of u are sub-Gaussian and independent of X. Then κ_3  can be chosen to be a large enough constant. 

The existence of κ_4,κ_5 is a consequence of the Hanson-Wright inequality. We consider X'u_2 conditional on X. By Theorem 6.3.2 in <cit.> (proved using the Hanson-Wright inequality) applied to the conditional probability given X and by the definition of sub-Gaussian variables, we have that with probability one, for any z>0,

    P (|X'u_2-σX_F|>z X| X )≤ C_1 exp(-C_2 z^2 ),

where C_1,C_2>0 are constants depending only on the sub-Gaussian norm of u. Choosing z=√(log(T)), we obtain

    P (X'u_2≤σX_F+ √(log(T))X)≤ o(1).

Notice that EX_F^2=∑_j,tEX_j,t^2=O(JT) and X_F^2-EX_F^2=∑_t(∑_j(X_j,t^2-EX_j,t^2)). By the boundedness of X_j,t, ∑_j(X_j,t^2-EX_j,t^2)=O(J) has mean zero and variance  bounded by O(J^2). It follows by the central limit theorem that  X_F^2-EX_F^2=∑_t(∑_j(X_j,t^2-EX_j,t^2))=√(T)J. Since EX_F^2=O(JT), there exists a constant C_3>0 such that P(X_F^2<C_3 JT)≥ 1-o(1). Again by the random matrix theory (e.g., Theorem 5.39 and Remark 5.40 in <cit.>), X=O_P(√(J+T)). We have that with probability approaching one, X'u_2 ≤σX_F +√(log(T))X≤√(C_3 JT)σ+O(√((J+T)log(T))). Assume J≫log(T). Since (J+T)log(T) ≪ JT, there exists κ_4 satisfying P(X'u_2≤κ_4 √(JT))≥ 1-o(1). The same argument holds for κ_5; we just repeat the same argument with T replaced by |H| once we realize X̃'ũ-X'u=∑_t∈ HX̃_t ũ_t- ∑_t∈ H X_t u_t, where X_t denotes the t-th row of X. (Analogous notations apply to X̃_t,u_t, etc.)

To find κ_6, we observe that X̃'X̃-X'X=∑_t∈ HX̃_t X̃_t' -∑_t∈ H X_t X_t'. By the random matrix theory (e.g., Theorem 5.39 and Remark 5.40 in <cit.>), ∑_t∈ H (X_t X_t' -EX_t X_t')≤ C_4 √(J |H|) with probability approaching one, where C_4>0 is a constant.  Again assuming that EX_t X_t' has eigenvalues bounded by a constant C_5>0, we have that ∑_t∈ H X_t X_t' ≤ C_5 |H|+C_4 √(J |H|) with probability approaching one. Recall the elementary inequality √(J|H|)≤ (J+|H|)/2. Therefore, we can simply set κ_6=4(C_4+C_5) and obtain that with probability approaching one, ∑_t∈ H X_t X_t' ≤ (|H|+J)κ_6/2. Similarly, we can get a bound for ∑_t∈ HX̃_t X̃_t'.

The above analysis is  based on simple sub-Gaussian or boundedness assumptions on (X,u). However, the purpose is to show that even in the simple case, Ridge regression can exhibit stability without being consistent (due to Lemma <ref>). We expect similar results to hold in more general or complicated data-generating processes. We leave this extension for future research.





§ SIMULATION STUDY


This section presents simulation evidence on the finite sample properties of our inference procedures. We consider the three CSC methods used in the empirical application in Section <ref>: difference-in-differences, canonical SC, and constrained Lasso with K=1.

We consider different data generating processes (DGPs) for the treated unit all of which specify the treated outcome as a weighted combination of the control outcomes:

    Y_1t=∑_j=2^J+1w_jY^N_jt+u_t    if t≤ T_0, 
    θ_t+∑_j=2^J+1w_jY^N_jt+u_t    if t >T_0,

where u_t=ρ_uu_t-1+v_t, v_tiid∼ N(0,1-ρ_u^2). Similar to <cit.>, the control outcomes are generated using a factor model:

    Y^N_jt=λ_1j+F_1t+λ_2jF_2t+ϵ_jt,

where λ_1j=(j-1)/J, λ_2j=(j-1)/J, F_1tiid∼ N(0,1), and ϵ_jt=ρ_ϵϵ_jt-1+ξ_jt, ξ_jtiid∼ N(0,1-ρ_ϵ^2). In the simulations, we vary ρ_u, ρ_ϵ, T_0, J, and F_2t. The DGPs differ with respect to the specification of the weights w.


We set T_∗=1 and consider the problem of testing the null hypothesis of a zero effect:

    H_0: θ_T=0.

The p-values are computed using the set of moving block permutations Π_→. The nominal level is α=0.1. 

We first analyze the performance of our procedure with stationary data, letting F_2tiid∼ N(0,1). Table <ref> presents simulation evidence on the size properties of our method when the data are iid (ρ_u=ρ_ϵ=0), which implies exchangeability of the residuals (cf. Theorem <ref>). As expected, our procedure achieves exact size control, irrespective of whether or not the model for P_t^N is correctly specified.  To study the finite sample performance with dependent data, we set ρ_u=ρ_ϵ=0.6. Table <ref> shows that our method exhibits close-to-correct size under correct specification and under misspecification, confirming the theoretical results on the robustness of our procedure under estimator stability and stationarity.


To investigate the performance with non-stationary data, we use trending factors: F_2t∼ N(t,1). Tables <ref> and <ref> show that under correct specification, our method exhibits excellent size properties. However, unlike in stationary settings, misspecification can cause size distortions. This finding is expected given the theoretical results and discussions in Section <ref>.

Figure <ref> displays power curves for a setting where T_0=19 and J=50 as in our empirical application, ρ_u=ρ_ϵ=0.6, and F_2tiid∼ N(0,1). Under correct specification, our method exhibits excellent small sample power properties and comes close to achieving the oracle power bound based on the true marginal distribution of u_t. Moreover, we find that imposing additional constraints when estimating P_t^N (e.g., using SC instead of the more general constrained Lasso) does not improve power when these additional restrictions are correct but can cause power losses when they are not.







§ PROOFS




 §.§ Additional Notation

We introduce some additional notations that will be used in the proofs.
For a,b∈ℝ, a∨ b=max{a,b}.  For two positive sequences a_n, b_n (indexed by n), we use  a_n≪ b_n to denote a_n=o(b_n).  We use Φ(·) to denote the cumulative distribution function of the standard normal distribution. Unless stated otherwise, · denotes the Euclidean norm for vectors or the spectral norm for matrices. We use d= to denote equal in distribution. 




 §.§ Proof of Theorem <ref>

The proof proceeds by verifying the high-level conditions in the following lemma. Let n=|Π| so that n=T! if Π=Π_all and n=T if Π=Π_→.

Let {δ_1n, δ_2n, γ_1n, γ_2n} be sequences of numbers converging to zero. Assume the following conditions.
 

  (E) With
probability 1-γ_1n:  the randomization distribution 

    F̃(x):=1/n∑_π∈Π1{S(u_π)< x},

is approximately ergodic for F(x)=P(S(u)< x), namely 

    sup_x∈ℝ|F̃(x)-F(x)|≤δ_1n,


 

  (A) With
probability 1-γ_2n, estimation errors are small:
 

  * the mean squared error is small, n^-1∑_π∈Π[S(û_π)-S(u_π)]^2≤δ_2n^2; 

  * the pointwise error at π=Identity is small, |S(û)-S(u)|≤δ_2n; 

  * The pdf of S(u) is bounded above by a constant D. 


Suppose in addition that the null hypothesis is true. Then, the approximate p-value obeys for any α∈(0,1) 

    |P(p̂≤α)- α|  ≤ 4δ_1n+ 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) +γ_1n+γ_2n.





With this result at hand, the proof of the theorem is a consequence of following four lemmas, which
verify the approximate ergodicity conditions (E) and 
conditions on the estimation error (A) of Lemma <ref>. 
Putting the bounds together and optimizing the error yields the result of the theorem.


The following lemma verifies approximate ergodicity (E)  (which  allows for large T_*) for the case of moving block permutations. 


Let Π be the moving block
permutations. Suppose that {u_t}_t=1^T is stationary and
strong mixing. Assume the following conditions:
(1) ∑_k=1^∞α_mixing(k) is bounded by a constant
M,  (2) T_0≥ T_*+2, and (3) S(u) has bounded pdf. Then there exists a constant M'>0 depending only on M such that
for any δ_1n>0, 

    P(sup_x∈ℝ|F̃(x)-F(x)|≤δ_1n)≥1-γ_T,

where γ_T=(M'√(T_*/T_0)log T_0+T_*+1/T_0+T_*)/δ_1n. 





The following lemma verifies approximate ergodicity (E)  (which  allows for large T_*) for the case of iid permutations. 



Let Π be
the set of all permutations. Suppose that {u_t}_t=1^T is
iid. Assume that S(u) only depends on the last T_* entries
of u. If T_0≥ T_*+2, then

    P(sup_x∈ℝ|F̃(x)-F(x)|≤δ_1n)≥1-γ_T,

where γ_T=√(π/(2⌊ T/T_*⌋) )/δ_1n. 







The following lemma  verifies the condition on the estimation error (A) for moving block permutations.



  Consider
moving block permutations Π. Let T_* be fixed. Suppose that for
some constant Q>0, |S(u)-S(v)|≤ QD_T_*(u-v)_2
for any u,v∈ℝ^T and D_T_*:= Blockdiag(0_T_*,I_T_*).
Then Condition (A) (1)-(2) is satisfied if there exist sequences γ_T,δ_2n=o(1)
such that with probability at least 1-γ_T, 
    P̂^N-P^N_2/√(T)≤δ_2n and  |P̂^N_t-P_t|≤δ_2n  for   T_0+1≤ t≤ T.
 




The following lemma  verifies the condition on the estimation error (A) for iid permutations.



 Consider
the set of all permutations Π. Let T_* be fixed. Suppose
that for some constant Q>0, |S(u)-S(v)|≤ QD_T_*(u-v)_2
for any u,v∈ℝ^T and D_T_*:= Blockdiag(0,I_T_*). Then Condition (A) (1)-(2)
is satisfied if there exist sequences γ_T,δ_2n=o(1)
such that with probability at least 1-γ_T, 
    P̂^N-P^N_2/√(T)≤δ_2n and  |P̂^N_t-P_t|≤δ_2n  for   T_0+1≤ t≤ T.
 



Now we conclude the proof of Theorem <ref>.

For the moving block permutations, let δ_1n=(T_*/T_0)^1/4. Then we apply Lemma <ref> together with Lemmas  <ref> and <ref>, obtaining

    |P(p̂≤α)- α|    ≤  4δ_1n+ 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) +γ_1n+γ_2n
       ≤ 4δ_1n+ 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) + (M'√(T_*/T_0)log T_0+T_*+1/T_0+T_*)/δ_1n+ γ_2n
       ≤ 4(T_*/T_0)^1/4 + 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) 
            + (M'√(T_*/T_0)log T_0+T_*+1/T_0+T_*)(T_*/T_0)^-1/4 + γ_2n.

The final result for moving block permutations follows by straight-forward computations and the observations that δ_2n =O(√(δ_2n)) (due to δ_2n=o(1)).

For iid permutations, we also use δ_1n=(T_*/T_0)^1/4. Then we apply Lemma <ref>  together with Lemmas <ref> and <ref>, 
 obtaining

    |P(p̂≤α)- α|    ≤   4δ_1n+ 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) +γ_1n+γ_2n
       ≤  4δ_1n+ 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) + √(2π/⌊ T/T_*⌋)/δ_1n+ γ_2n
       ≤  4(T_*/T_0)^1/4  + 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) +√(2π/⌊ T/T_*⌋) (T_*/T_0)^-1/4 + γ_2n
       ≲ (T_*/T_0)^1/4 +δ_2n+ √(δ_2n)+ γ_2n.

This completes the proof for iid permutations. 








  §.§.§ Proof of Lemma <ref>


The proof proceeds in two steps.[The proof follows from the same arguments as the proof of Theorem 2 in <cit.> modified to the problem of permuting residuals and is presented for completeness.] 

Step 1: We bound the difference between the p-value and the oracle p-value, F̂(S(û))-F(S(u)).

Let ℳ be the event that the conditions (A) and (E) hold. By assumption, 

    P(ℳ)≥1-γ_1n-γ_2n.


Notice that on the event ℳ, 

    |F̂(S(û))-F(S(u))|    ≤|F̂(S(û))-F(S(û))|+|F(S(û))-F(S(u))|
        (i)≤sup_x∈ℝ|F̂(x)-F(x)|+D|S(û)-S(u)|
       ≤sup_x∈ℝ|F̂(x)-F̃(x)|+sup_x∈ℝ|F̃(x)-F(x)|+D|S(û)-S(u)|
       ≤sup_x∈ℝ|F̂(x)-F̃(x)|+δ_1n+D|S(û)-S(u)|
       ≤sup_x∈ℝ|F̂(x)-F̃(x)|+δ_1n+Dδ_2n,

where (i) holds by the fact that the bounded pdf of S(u) implies
Lipschitz property for F. 

Let A={π∈Π: |S(û_π)-S(u_π)|≥√(δ_2n)}.
Observe that on the event ℳ, by Chebyshev inequality

    |A|δ_2n≤∑_π∈Π(S(û_π)-S(u_π))^2≤ nδ_2n^2

and thus |A|/n≤δ_2n. Also observe that on the event ℳ,
for any x∈ℝ, 

    |F̂(x)-F̃(x)| 
    
        ≤1/n∑_π∈ A|1{ S(û_π)< x} -1{ S(u_π)< x }|+1/n∑_π∈(Π\ A)|1{ S(û_π)< x} -1{ S(u_π)< x}|
       (i)≤|A|/n+1/n∑_π∈(Π\ A)1{|S(u_π)-x|≤√(δ_2n)}≤|A|/n+1/n∑_π∈Π1{|S(u_π)-x|≤√(δ_2n)}
       ≤|A|/n+P(|S(u)-x|≤√(δ_2n))
            +sup_z∈ℝ|1/n∑_π∈Π1{|S(u_π)-z|≤√(δ_2n)} -P(|S(u)-z|≤√(δ_2n))| 
        =|A|/n+P(|S(u)-x|≤√(δ_2n)) 
           +sup_x∈ℝ|[F̃(z+√(δ_2n))-F̃(z-√(δ_2n))]-[F(z+√(δ_2n))-F(z-√(δ_2n))]| 
       ≤|A|/n+P(|S(u)-x|≤√(δ_2n)) +2sup_z∈ℝ|F̃(z)-F(z)|
       (ii)≤|A|/n+2D√(δ_2n)+2δ_1n(iii)≤δ_1n+ 2 δ_2n+ 2 D√(δ_2n),

where (i) follows by 
the elementary inequality of |1{S(û_π)< x}-1{S(u_π)<x}|≤1{|S(u_π)-x|≤|S(û_π)-S(u_π)|},
(ii) follows by the bounded pdf of S(u) and (iii) follows by |A|/n≤δ_2n.
Since the above display holds for each x∈ℝ, it follows
that on the event ℳ, 

    sup_x∈ℝ|F̂(x)-F̃(x)|≤δ_1n+ 2 δ_2n+ 2D√(δ_2n).


We combine (<ref>) and (<ref>)
and obtain that on the event ℳ, 

    |F̂(S(û))-F(S(u))|≤
    2δ_1n+ 2 δ_2n+ D(δ_2n+2 √(δ_2n)).


Step 2: Here we derive the desired result. Notice that 

    |P(1-F̂(S(û))≤α)-α|
         =|E(1{ 1-F̂(S(û))≤α} -1{ 1-F(S(u))≤α})|
       ≤ E|1{ 1-F̂(S(û))≤α} -1{ 1-F(S(u))≤α}|
       (i)≤P(|F(S(u))-1+α|≤|F̂(S(û))-F(S(u))|)
       ≤ P(|F(S(u))-1+α|≤|F̂(S(û))-F(S(u))|  and ℳ)+P(ℳ^c)
       (ii)≤P(|F(S(u))-1+α|≤
      2δ_1n+ 2 δ_2n+ D(δ_2n+ 2√(δ_2n)))+P(ℳ^c)
       (iii)≤
    4δ_1n+ 4 δ_2n+ 2 D(δ_2n+ 2√(δ_2n)) +γ_1n+γ_2n,

where (i) follows by the elementary inequality |1{1-F̂(S(û))≤α}-1{1-F(S(u))≤α}|≤1{|F(S(u))-1+α|≤|F̂(S(û))-F(S(u))|},
(ii) follows by (<ref>), (iii) follows
by the fact that F(S(u)) has the uniform distribution on (0,1)
and hence has pdf equal to 1, and  by (<ref>).
The proof is complete. 





  §.§.§ Proof of Lemma <ref>

We define 

    s_t=
    (∑_s=t^t+T_*-1|u_s|^q)^1/q    if 1≤ t≤ T_0
    
    (∑_s=t^T|u_s|^q+∑_s=1^t-T_0-1|u_s|^q)^1/q    otherwise.


It is straight-forward to verify that 

    {S(u_π): π∈Π}={ s_t: 1≤ t≤ T} .


Let α̃_mixing be the strong-mixing coefficient for
{s_t}_t=1^T_0. Notice that {s_t}_t=1^T_0
is stationary (although {s_t}_t=1^T is clearly not). Let F̌(x)=T_0^-1∑_t=1^T_01{s_t≤ x}.
The bounded pdf of S(u) implies the continuity of F(·).
It follows, by Proposition 7.1 of <cit.>, that

    E(sup_x∈ℝ|F̌(x)-F(x)|^2)≤1/T_0(1+4∑_k=0^T_0-1α̃_mixing(t))(3+log T_0/2log2)^2.


Notice that α̃_mixing(t)≤2 and that 
α̃_mixing(t)≤α_mixing(max{t-T_*,0})
so that

    ∑_k=0^T_0-1α̃_mixing(t)    =    ∑_k=0^T_*α̃_mixing(t)+∑_k=T_*+1^T_0-1α̃_mixing(t)≤2(T_*+1)+∑_k=1^T_0-T_*-1α_mixing(k)
       ≤    2(T_*+1)+∑_k=1^∞α_mixing(k).

Since ∑_k=1^∞α_mixing(k) is bounded by M, it follows
by (<ref>) that 

    E(sup_x∈ℝ|F̌(x)-F(x)|^2)≤ B_T: =  1+ 4 (2 (T_* +1) + M) /T_0(3+log T_0/2log2)^2.

By Liapunov's inequality, 

    E(sup_x∈ℝ|F̌(x)-F(x)|)≤√(E(sup_x∈ℝ|F̌(x)-F(x)|^2))≤√(B_T).


Since (T_0+T_*)F̃(x)-T_0F̌(x)=∑_t=T_0+1^T_0+T_*1{s_t≤ x},
it follows that 

    sup_x∈ℝ|F̃(x)-F̌(x)|     =sup_x∈ℝ|(T_0/T_0+T_*F̌(x)+1/T_0+T_*∑_t=T_0+1^T_0+T_*1{s_t≤ x})-F̌(x)|
        =sup_x∈ℝ|1/T_0+T_*F̌(x)+1/T_0+T_*∑_t=T_0+1^T_0+T_*1{s_t≤ x}| ≤T_*+1/T_0+T_*,

where the last inequality follows by sup_x∈ℝ|F̌(x)|≤1
and the boundedness of the indicator function. Combining the above
two displays, we obtain that

    E(sup_x∈ℝ|F̃(x)-F(x)|)≤√(B_T) +T_*+1/T_0+T_*.

The desired result follows by Markov's inequality.













  §.§.§ Proof of Lemma <ref>


The proof follows by an argument given by <cit.> for subsampling.
We give a complete argument for our setting here for clarity and completeness.

Recall that Π is the set of all bijections π on {1,...,T}.
Let k_T=⌊ T/T_*⌋. Define the blocks of indices 
    b_i=(T-iT_*+1,T-iT_*+2,...,T-iT_*+T_*)∈ℝ^T_*,         i = 1,...., k_T



Since S(u) only depends on u_b_1,  the last T_* entries of u,
we can define 
    Q(x;u_b_1)=1{S(u)≤ x}-F(x).

Therefore, 

    F̃(x)-F(x)=1/|Π|∑_π∈ΠQ(u_π(b_1);x).

Define  π(b_i):=π_| b_i (b_i) to mean the restriction of the permutation map π:
{1, …, T}→{1 , …, T} to the domain b_i.


Notice that for 1≤ i≤ k_T, the value of ∑_π∈ΠQ(u_π(b_i);x)
does not depend on i. It follows that 

    F̃(x)-F(x)=1/|Π|∑_π∈ΠQ(u_π(b_1);x)     =1/k_T∑_i=1^k_T(1/|Π|∑_π∈ΠQ(u_π(b_i);x))
        =1/|Π|∑_π∈Π[1/k_T∑_i=1^k_TQ(u_π(b_i);x)].


Hence by Jensen's inequality

    E(sup_x∈ℝ|F̃(x)-F(x)|) 
        ≤1/|Π|∑_π∈ΠE(sup_x∈ℝ|1/k_T∑_i=1^k_TQ(u_π(b_i);x)|).


To compute the above expectation,
we observe that for any π∈Π,

    E(sup_x∈ℝ|1/k_T∑_i=1^k_TQ(u_π(b_i);x)|)     =∫_0^1P(sup_x∈ℝ|1/k_T∑_i=1^k_TQ(u_π(b_i);x)|>z)dz
       ≤∫_0^12exp(-2k_Tz^2)dz  < ∫_0^∞2exp(-2k_Tz^2)dz = √(π/(2k_T)),

where the first inequality follows by the Dvoretsky-Kiefer-Wolfwitz inequality (e.g.,
Theorem 11.6 in <cit.>) and the fact that
for any π∈Π, {Q(u_π(b_i);x)}_i=1^k_T is a
sequence of iid random variables (since π is a bijection and
{b_i}_i=1^k_T are disjoint blocks of indices); the last equality follows
from the properties of the normal density. Therefore, the above two
display imply that 

    E(sup_x∈ℝ|F̃(x)-F(x)|)≤√(π/(2k_T)).


The desired result follows by Markov's inequality. 








  §.§.§ Proof of Lemma <ref>

Due to the Lipschitz property of S(·), we have 

    ∑_π∈Π[S(û_π)-S(u_π)]^2   ≤ Q∑_π∈Π‖ D_T_*(û_π-u_π)‖ _2^2  =Q∑_π∈Π∑_t=T_0+1^T_0+T_*(û_π(t)-u_π(t))^2
        =Q∑_t=T_0+1^T_0+T_*∑_π∈Π(û_π(t)-u_π(t))^2=QT_*û-u_2^2 = QT_* P̂^N - P^N^2

where the penultimate equality follows by the observation that for moving block permutation
Π, 
    ∑_π∈Π(û_π(t)-u_π(t))^2=û-u_2^2.

Hence condition (A) (1)  follows with a rescaled value of δ_n. Condition
(A) (2) holds by the Lipschitz property of S(·): 

    |S(û)-S(u)|≤ QD_T_*(û-u)_2≤ Q√(∑_t=T_0+1^T_0+T_*(û_t-u_t)^2)

Hence, Condition (A) (2) follows
since P̂^N_t - P^N_t = |û_t-u_t|≤δ_n for T_0+1≤ t≤ T
with high probability. The proof is complete. 




  §.§.§ Proof of Lemma <ref>

For t,s∈{1,...,T}, we define A_t,s={π∈Π: π(t)=s}.
Recall that Π is the set of all bijections on {1,...,T}.
Thus, |A_t,s|=(T-1)!. It follows that for any t∈{1,...,T},

    ∑_π∈Π(û_π(t)-u_π(t))^2    =∑_s=1^T∑_π∈ A_t,s(û_π(t)-u_π(t))^2
        =∑_s=1^T∑_π∈ A_t,s(û_s-u_s)^2=∑_s=1^T|A_t,s|(û_s-u_s)^2=(T-1)!×û-u_2^2.


Due to the Lipschitz property of S(·), we have that

    1/|Π|∑_π∈Π[S(û_π)-S(u_π)]^2≤Q/|Π|∑_π∈Π‖ D_T_* (û_π-u_π)‖ _2^2
      =Q/|Π|∑_π∈Π∑_t=T_0+1^T_0+T_*(û_π(t)-u_π(t))^2
           = Q/|Π|∑_t=T_0+1^T_0+T_*∑_π∈Π(û_π(t)-u_π(t))^2=Q/|Π|T_*(T-1)!×û-u_2^2=QT^-1T_*û-u_2^2,

where the penultimate equality follows by (<ref>) and 
the last equality  follows by |Π|=T!. Thus, part 1 of Condition (A) follows since
T_* is fixed.

To see part 2 of Condition (A), notice that the Lipschitz property
of S(·) implies 

    |S(û)-S(u)|≤ QD_T_*(û-u)_2≤ Q√(∑_t=T_0+1^T_0+T_*(û_t-u_t)^2).

Hence, part 2 of Condition (A) follows since |û_t-u_t|≤δ_n
for T_0+1≤ t≤ T with high probability. The proof is complete. 





 §.§ Proof of Theorem <ref>


We first state auxiliary results.

	Let {W_t}_t=1^T be a stationary
	and β-mixing sequence with coefficient (·).
	Let G(x)=P(W_t≤ x). 
	Then for any positive integer 1≤ m≤ T/2, we
	have 
	
    E(sup_x∈ℝ|T^-1∑_t=1^T[{W_t≤ x}-G(x)]|)≤  2√(T)(m)+ √(π m/(2T))+(m-1)/T.



The next lemma was derived by <cit.>; see Theorem 16.2.1 in <cit.> and Lemma 7.1 of <cit.> for popular versions. We state the result to make the exposition more self-contained. The proof is omitted. 

	 Let (W,R) be random vectors defined in the same probability space. Let  P_W,R denote the probability distribution of (W,R). Let P_R and P_W denote the probability distributions of R and W, respectively. Let ·_TV denote the total variation metric. Define the β-mixing coefficient β(W,R)=P_W,R-P_W ⊗ P_R_TV/2. Then the probability space can be extended to construct random element  such that (1)  and R are independent, (2)  and W have the same distribution and (3) P(≠ W)≤β(W,R).

















Now we prove Theorem <ref>. In this proof, universal constants refer to constants that depend
only on D_1,D_2,D_3>0. Define (x)=R^-1∑_j=1^R_j(x),
where

    _j(x)=m^-1∑_t∈ H_j{ϕ(Z_t,...,Z_t+T_*-1;β̂())≤ x} .


We note that under  Π=Π_→, F̂(x) can be written as

    (x)=T^-1(∑_t=1^T_0{ϕ(Z_t,...,Z_t+T_*-1;β̂())≤ x} +∑_t=T_0+1^T_0+T_*{ϕ(Z_q(t),...,Z_q(t+T_*-1);β̂())≤ x}),

where q(t)=t{t≤ T}+(t-T){t> T}.

The rest of the proof proceeds in 4 steps. The first three steps bound
sup_x∈|(x)-Ψ(x;(__R))|, where we recall that Ψ(x;β)=P(ϕ(Z_t,...,Z_t+T_*-1;β)≤ x).
The fourth step derives the desired result. 

Step 1: bound sup_x∈ℝ|_j(x)-Ψ(x;β̂(__j))|. 

Let A_j=⋃_t∈ H_j{t,...,t+T_*-1}. Since k>T_*,
we have that A_j⊂_j and min_t∈ A_j, s∈_j^c|t-s|≥ k-T_*+1.
This means that {Z_t}_t∈_j^c and {Z_t}_t∈ A_j
have a gap of at least k-T_*+1 time periods. By Lemma <ref> (applied with W={Z_t}_t∈_j^c and R={Z_t}_t∈ A_j), there exist random elements
{_t}_t∈ A_j (on an enlarged probability space) such
that (1) {_t}_t∈ A_j is independent of {Z_t}_t∈_j^c,
(2) {_t}_t∈ A_jd={Z_t}_t∈ A_j
and (3) P({_t}_t∈ A_j≠{Z_t}_t∈ A_j)≤(k-T_*+1).
Since {Z̃_t}_t∈_j is independent of the data,
we can construct {_t}_t∈ A_j such that it is also
independent of __j. 

Define the event 

    _j={{_t}_t∈ A_j={Z_t}_t∈ A_j}⋂{sup_x∈|∂Ψ(x;(__j))/∂ x|≤ξ_T}
    ⋂{max_π∈Π|S(^π,β̂())-S(^π,β̂(__j))|≤ϱ_T(|_j|)}

 as well as the functions 

    _j(x)=m^-1∑_t∈ H_j{ϕ(_t,...,_t+T_*-1;β̂(__j))≤ x}
    Ḟ_j(x)=m^-1∑_t∈ H_j{ϕ(Z_t,...,Z_t+T_*-1;β̂(__j))≤ x} .


By the construction of {_t}_t∈ A_j and Assumptions
<ref> and <ref>, P(_j^c)≤(k-T_*+1)+γ_1,T+γ_2,T. 

Notice that conditional on (__j), ϕ(_t,...,_t+T_*-1;β̂(__j)) is a stationary β-mixing across t∈ H_j with mixing coefficient β̃_ mixing(i)≤(i-T_*+1)
for i≥ T_*. Moreover, since  (_t,...,_t+T_*-1) is independent of __j, we have

    Ψ(x;β̂(__j))=P(ϕ(_t,...,_t+T_*-1;β̂(__j)) ≤ x |β̂(__j)  ).

  Hence, by Lemma <ref>, we have
that for any m_1≤ m/2,

    E(sup_x∈ℝ|_j(x)-Ψ(x;β̂(__j))|)≤2m^1/2(m_1-T_*+1)+√(π m_1/(2m))+(m_1-1)/m.


We shall choose m_1 later. Observe that on the event _j,
_j(·)=Ḟ_j(·). Therefore, 

    E(a_j)≤2m^1/2(m_1-T_*+1)+√(π m_1/(2m))+(m_1-1)/m+2P(_j^c),

where a_j=sup_x∈ℝ|Ḟ_j(x)-Ψ(x;β̂(__j))|. 

Now we bound sup_x∈|Ḟ_j(x)-_j(x)|. Fix an
arbitrary x∈. Observe that on the event _j, 

    |_j(x)-Ḟ_j(x)|
        =|m^-1∑_t∈ H_j({ϕ(Z_t,...,Z_t+T_*-1;β̂())≤ x} -{ϕ(Z_t,...,Z_t+T_*-1;β̂(__j))≤ x})|
       ≤ m^-1∑_t∈ H_j|{ϕ(Z_t,...,Z_t+T_*-1;β̂())≤ x} -{ϕ(Z_t,...,Z_t+T_*-1;β̂(__j))≤ x}|
       (i)≤m^-1∑_t∈ H_j{|ϕ(Z_t,...,Z_t+T_*-1;β̂(__j))-x|≤ϱ_T(|_j|)}
        = m^-1∑_t∈ H_j{ϕ(Z_t,...,Z_t+T_*-1;β̂(__j))≤ x+ϱ_T(|_j|)}
            -m^-1∑_t∈ H_j{ϕ(Z_t,...,Z_t+T_*-1;β̂(__j))<x -ϱ_T(|_j|)}
        <Ḟ_j(x+ϱ_T(|_j|))-Ḟ_j(x-2ϱ_T(|_j|))
       ≤Ψ(x+ϱ_T(|_j|);β̂(__j))-Ψ(x-2ϱ_T(|_j|);β̂(__j))+2a_j(ii)≤3ξ_Tϱ_T(|_j|)+2a_j,

where (i) follows by the elementary inequality |{x≤ z}-{y≤ z}|≤{|y-z|≤|x-y|}
for any x,y,z∈ and (ii) follows by the definition of _j.
Since the above bound holds for any x∈ and |_j|≤ m+2k,
we have that on the event _j, 

    sup_x∈|_j(x)-Ḟ_j(x)|≤3ξ_Tϱ_T(m+2k)+2a_j.


By the definition of a_j, this means that on the event _j,

    sup_x∈ℝ|_j(x)-Ψ(x;β̂(__j))|≤3ξ_Tϱ_T(m+2k)+3a_j.


By (<ref>) and the fact that _j(·)
and Ψ(·,·) take values in [0,1], we have that for
a universal constant C_1>0, 


    E(sup_x∈ℝ|_j(x)-Ψ(x;β̂(__j))|)
       ≤3ξ_Tϱ_T(m+2k)+3E(a_j)+2P(_j^c)
       ≤3ξ_Tϱ_T(m+2k)+6m^1/2(m_1-T_*+1)+3√(π m_1/(2m))+3(m_1-1)/m+8P(_j^c)
       (i)≤C_1(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T),

where (i) follows by P(_j^c)≤(k-T_*+1)+γ_1,T+γ_2,T
and m_1/m≤√(m_1/m). 

Step 2: bound R^-1∑_j=1^Rsup_x∈|Ψ(x;β̂(__j))-Ψ(x;β̂(__R))|. 

Let ={Ż_t}_t=1^T satisfy that d=
and  is independent of (,{Z̃_t}_t=1^T).
Therefore, for any 1≤ j≤ R, 

    Ψ(x;β̂(__j))     =P(ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))≤ x|(__j))
       (i)=P(ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))≤ x|,{Z̃_t}_t=1^T),

where (i) follows by the fact that  and (,{Z̃_t}_t=1^T)
are independent. (This is the identity that E(f(X;g(Y))| g(Y))=E(f(X;g(Y))| Y)
for any measurable functions f and g if X and Y are independent.
To see this, simply notice that the distribution of X given g(Y)
and the distribution of X given Y are both equal to the unconditional
distribution of X.)

Define the event 

    _j={sup_x∈|∂Ψ(x;(__j))/∂ x|≤ξ_T} .


Clearly, P(_j)≥1-γ_2,T by Assumption <ref>.
Therefore, we have that on the event _j, 

    |Ψ(x;β̂(__j))-Ψ(x;β̂(__R))|
        =|E({ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))≤ x} -{ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))≤ x}|,{Z̃_t}_t=1^T)|
       ≤ E(|{ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))≤ x} -{ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))≤ x}||,{Z̃_t}_t=1^T)
       (i)≤E[{|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))-x|
                   ≤|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))|}|,{Z̃_t}_t=1^T]
        =P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))-x|
                   ≤|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))||,{Z̃_t}_t=1^T]
       ≤ P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))-x|≤2ϱ_T(|_j|)|,{Z̃_t}_t=1^T]
         +P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))|>2ϱ_T(|_j|)|,{Z̃_t}_t=1^T]
        =Ψ(x+2ϱ_T(|_j|);(__j))-Ψ(x-2ϱ_T(|_j|);(__j))
         +P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))|>2ϱ_T(|_j|)|,{Z̃_t}_t=1^T]
       (ii)≤4ξ_Tϱ_T(m+2k)
         +P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))|>2ϱ_T(|_j|)|,{Z̃_t}_t=1^T],

where (i) follows by the elementary inequality |{x≤ z}-{y≤ z}|≤{|y-z|≤|x-y|}
for any x,y,z∈ and (ii) follows by |_j|≤ m+2k and
the definition of _j. Since the above bound does not depend
on x, it holds uniformly in x∈ on the event _j.
Since Ψ(·,·) is also bounded by one, we have that 

    E(sup_x∈|Ψ(x;β̂(__j))-Ψ(x;β̂(__R))|)
       ≤4ξ_Tϱ_T(m+2k)
           +2P(_j^c)+P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))|>2ϱ_T(|_j|)]
       ≤4ξ_Tϱ_T(m+2k)+2P(_j^c)
           +P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__R))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;())|>ϱ_T(|_j|)]
           +P[|ϕ(Ż_T_0+1,...,Ż_T_0+T_*;(__j))-ϕ(Ż_T_0+1,...,Ż_T_0+T_*;())|>ϱ_T(|_j|)]
       (i)≤4ξ_Tϱ_T(m+2k)+2P(_j^c)+2γ_1,T(ii)≤4ξ_Tϱ_T(m+2k)+2γ_1,T+2γ_2,T,

where (i) follows by Assumption <ref> and the fact
that |_j|=|_R| and (ii) follows by P(_j)≥1-γ_2,T.
Since the above bound holds for all 1≤ j≤ R,  we have

    E(R^-1∑_j=1^Rsup_x∈|Ψ(x;β̂(__j))-Ψ(x;β̂(__R))|)≤4ξ_Tϱ_T(m+2k)+2γ_1,T+2γ_2,T.


Step 3: bound sup_x∈|(x)-Ψ(x;(__R))|.

By (<ref>), we notice that 

    sup_x∈|T(x)-m∑_j=1^R_j(x)|=sup_x∈|∑_t=T_0-mR+1^T_0{ϕ(Z_t,...,Z_t+T_*-1;β̂())≤ x}
    
    +∑_t=T_0+1^T_0+T_*{ϕ(Z_q(t),...,Z_q(t+T_*-1);β̂())≤ x}|≤ T_*+(T_0-mR)≤ T_*+R-1.


Moreover, by (<ref>) and (<ref>),
we have that 

    E(sup_x∈|R^-1∑_j=1^R_j(x)-Ψ(x;β̂(__R)))|
    ≤ C_2(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T)

for some universal constant C_2>0. 

The above two displays imply that

    sup_x∈|T/mR(x)-Ψ(x;β̂(__R))|≤T_*+R-1/mR
    
    +C_2(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T).


Since (x)∈[0,1], we have that 

    sup_x∈|(1-T/(mR))(x)|≤T/mR-1≤T-mR/mR≤T_*+R-1/mR.


Since mR≥ T_0/2 (due to R<T_0/2), we have (T_*+R-1)/(mR)≤2T_*T_0^-1+m^-1≲√(m_1/m).
Hence, the above two displays imply that for some universal constant
C_3>0, 

    E(sup_x∈|(x)-Ψ(x;β̂(__R))|)
    ≤ C_3(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T).


Step 4: derive the desired result. 

Let A_R be defined as in Step 1 with j=R. Following Step 1,
we can construct random elements {_t}_t∈ A_R (on an
enlarged probability space) such that (1) {_t}_t∈ A_R
is independent of __R, (2) {_t}_t∈ A_Rd={Z_t}_t∈ A_R
and (3) P({_t}_t∈ A_R≠{Z_t}_t∈ A_R)≤(k-T_*+1).

Define (__R)=ϕ(_T_0+1,...,_T_0+T_*;(__R)).
Since {T_0+1,...,T_0+T_*}⊂ A_R, we have that (_T_0+1,...,_T_0+T_*)
is independent of __R, which means that 

    P((__R)≤ x|(__R))=Ψ(x;(__R))    ∀ x∈.


Therefore, 

    conditional on (__R), Ψ((__R);(__R)) has uniform distribution on (0,1).


We also introduce the following notations to simplify the argument:

()=ϕ(_T_0+1,...,_T_0+T_*;())
and G()=ϕ(Z_T_0+1,...,Z_T_0+T_*;()). 

For arbitrary α∈(0,1) and c>0, we observe that 

    |P((G())<α|(__R))-α|
       (i)=|E({(G())<α}|(__R))-E({Ψ((__R);(__R))<α}|(__R))|
       ≤ E(|{(G())<α} -{Ψ((__R);(__R))<α}||(__R))
       (ii)≤E({|Ψ((__R);(__R))-α|≤|Ψ((__R);(__R))-(G())|}|(__R))
        =P(|Ψ((__R);(__R))-α|≤|Ψ((__R);(__R))-(G())||(__R))
       ≤ P(|Ψ((__R);(__R))-α|≤ c|(__R))
           +P(|Ψ((__R);(__R))-(G())|>c|(__R))
       (iii)≤2c+P(|Ψ((__R);(__R))-(G())|>c|(__R))

where (i) follows by (<ref>), (ii) follows
by the elementary inequality |{x<z}-{y<z}|≤{|y-z|≤|x-y|}
for any x,y,z∈ and (iii) follows by (<ref>).
Now we take expectation on both sides, obtaining 

    E|P((G())<α|(__R))-α|
       ≤2c+P(|Ψ((__R);(__R))-(G())|>c)
       ≤2c+P(|Ψ((__R);(__R))-(())|>c)+P({_t}_t∈ A_R≠{Z_t}_t∈ A_R)
       ≤2c+P(|Ψ((__R);(__R))-(())|>c)+(k-T_*+1)
       ≤2c+c^-1E|Ψ((__R);(__R))-(())|+(k-T_*+1)


Define the event 

    ={sup_x∈|∂Ψ(x;(__R))/∂ x|≤ξ_T}⋂{|()-(__R)| ≤ϱ_T(|_R|) } .






By Assumptions <ref> and <ref>,
P(^c)≤γ_1,T+γ_2,T. Therefore, 

    E|Ψ(();(__R))-Ψ((__R);(__R))|
        =E(|Ψ(();(__R))-Ψ((__R);(__R))|×_)
                   +E(|Ψ(();(__R))-Ψ((__R);(__R))|×_^c)
       (i)≤E(|Ψ(();(__R))-Ψ((__R);(__R))|×_)+2P(^c)
       ≤ξ_Tϱ_T(|_R|)+2P(^c)≤ξ_Tϱ_T(m+2k)+2γ_1,T+2γ_2,T,

where (i) follows by the fact that Ψ(·,·)∈[0,1]. Hence,
we have that 

    E|Ψ((__R);(__R))-(())|
       ≤ E|Ψ(();(__R))-(())|+E|Ψ(();(__R))-Ψ((__R);(__R))|
       ≤ Esup_x∈|Ψ(x;(__R))-(x)|+ξ_Tϱ_T(m+2k)+2γ_1,T+2γ_2,T
       (i)≤C_4(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T)

for a universal constant C_4>0, where (i) follows by (<ref>).

Now we combine (<ref>) and the above display.
We also choose 

    c≍√(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T).


Then we can find a universal constant C_5>0 such that 

    E|P((G())<α|(__R))-α|
    ≤ C_5√(ξ_Tϱ_T(m+2k)+m^1/2(m_1-T_*+1)+√(m_1/m)+(k-T_*+1)+γ_1,T+γ_2,T).


Now we choose m_1 satisfying m_1≍(log m)^1/D_3
and m^1/2(m_1-T_*+1)≲ m^-1. Hence, for some
universal constant C_6>0, 

    E|P((G())<α|(__R))-α|
    ≤ C_6√(ξ_Tϱ_T(m+2k))+C_6(m^-1(log m)^1/D_3)^1/4+C_6√((k-T_*+1))+C_6√(γ_1,T)+C_6√(γ_2,T).


Since m≍ T_0/R, the desired result follows once we notice
that p̂≥1-α and (G())<α are
the same event. 





  §.§.§ Proof of Lemma <ref>


	Define K=⌊ T/m⌋ and (x)=m^-1/2∑_r=1^m_r(x),
	where  _r(x)=K^-1/2∑_j=1^K[{W_(j-1)m+r≤ x}-G(x)] for 1≤ r≤ m.
	Let Δ(x)=∑_t=mK+1^T[{W_t≤ x}-G(x)]. Let L_T(x)=T^-1/2∑_t=1^T[{W_t≤ x}-G(x)]. Notice
	that 
	
    √(T)L_T(x)=√(mK)(x)+Δ(x).

	
	Since |{W_t≤ x}-G(x)|≤1, it follows that sup_x∈ℝ|Δ(x)|≤ T-mK≤ m-1
	and thus 
	
    sup_x∈ℝ|√(T)L_T(x)-√(mK)(x)|≤ m-1.

	
	By Berbee's coupling (Lemma <ref>), we
	can enlarge the probability space and define random variables {W̅_t}_t=1^mK
	such that (1) W̅_td=W_t for all 1≤ t≤ mT,
	(2) W̅_(j-1)m+r is independent across 1≤ j≤ K for
	all r and (3) P(⋃_t=1^mK{W̅_t≠ W_t})≤ mK(m)≤ T(m). 
	
	We now define (x)=m^-1/2∑_r=1^m_r(x), where _r(x)=K^-1/2∑_j=1^K[{W̅_(j-1)m+r≤ x}-G(x)]. 
	
	Since {W̅_(j-1)m+r}_j=1^K is independent, it follows
	by Dvoretzky-Kiefer-Wolfowitz inequality that for any z>0, 
	
    P(sup_x∈ℝ|_r(x)|>z)≤2exp(-2z^2).

	    Therefore, we have that 
	
    E(sup_x∈ℝ|_r(x)|)=∫_0^∞P(sup_x∈ℝ|_r(x)|>z) dz ≤ 2 ∫_0^∞exp(-2z^2) dz=√(π/2).

	Hence, we have that 
	
    E(sup_x∈ℝ|(x)|) ≤ m^-1/2∑_r=1^m E(sup_x∈ℝ|_r(x)|) ≤√(π m/2).







	Since (·)=(·) with probability at least 1-T(m), we have that 
	
    E(sup_x∈|(x)-(x)|)≤ 2T(m).

	Therefore, E(sup_x∈ℝ|(x)|)≤  2T(m)+ √(π m/2). 



	
	By (<ref>) and mK/T≤1, we have that 
	
    E(sup_x∈ℝ|L_T(x)|)≤  2T(m)+ √(π m/2)+(m-1)/√(T).

	
	The proof is complete. 




 §.§ Proof of Theorem <ref>



Recall   ^*=(Z_1^*,…,Z_T^*)' with Z_t^*=(Y^N_1t,Y^N_2t,…,Y^N_J+1t,X'_1t,…,X'_J+1t)' for 1≤ t≤ T_0+T_*. Let 

    p̂_^*=1- F̂( S(û(^*))),

where F̂( x;^*)=1/|Π|∑_π∈Π1{ S(û_π(^*) )<  x} and û(^*)=Y^N-P̂^N with P̂^N computed using ^*.
By the proof of Theorem <ref>, we have

    |P(p̂_^*≤α)- α| ≤ C ( δ̃_T + δ_T + √(δ_T) + γ_T),

where δ̃_T = (T_*/T_0)^1/4(log T) and the constant C depends on T_*, M and D,  but not on T.
It follows that 

    |P(p̂_^*> α)- (1-α) | ≤ C ( δ̃_T + δ_T + √(δ_T) + γ_T).

The desired result follows by  observing that θ_t∈𝒞_1-α(t) is the same event as p̂_^*> α; this is simply because (θ^0)=^*, where (θ^0) is defined in Section <ref>.







 §.§ Proof of Theorem <ref>


Let {S^(j)(û)}_j=1^n denoted the non-decreasing rearrangement of {S(û_π): π∈Π}, where n= |Π|, which we refer to as randomization quantiles. The p-value is 

    p̂ = 1/n∑_π∈Π1 ( S(û_π) ≥ S(û)).

Note that 
    1 (p̂≤α) = 1( S(û ) > S^(k)( û) ),

where k = k(α) =  n - ⌊ nα⌋ = ⌈ n(1-α)⌉. 

The proof proceeds in three steps. First, we show that exchangeability of the data implies exchangeability of the residuals. Second, we show that exchangeability of the residuals implies that P(p̂≤α)≤α. Third, we show that if there are no ties, α-1/n≤ P(p̂≤α). The proof follows from standard arguments <cit.>.


Step 1: By the  iid or exchangeability property of data, we have that 

    {û_t}_t=1^T{g(Z_t,β̂({Z_t}_t=1^T))}_t=1^Td={g(Z_π(t),β̂({Z_π(t)}_t=1^T)}_t=1^T.

Since β̂({Z_π(t)}_t=1^T) does not depend on π, we have 

    {g(Z_π(t),β̂({Z_π(t)}_t=1^T)}_t=1^T={û_π(t)}_t=1^T{g(Z_π(t),β̂({Z_t}_t=1^T)}_t=1^T.

Therefore, {û_π(t)}_t=1^Td={û_t}_t=1^T. 


Step 2: Note that Π_all and Π_→ form groups in the sense that Ππ = Π for all π∈Π. Therefore, the randomization quantiles are invariant, 

    S^(k(α))(û_π) = S^(k(α))(û),  for all π∈Π.

Therefore,

    ∑_π∈Π1( S(û_π ) > S^(k(α))(û_π) ) =∑_π∈Π1( S(û_π ) > S^(k(α))(û)) ≤α n.

Since 1( S(û) > S^(k(α))(û) ) is equal in distribution to  1( S(û_π) > S^(k(α))(û_π) )  for any π∈Π by exchangeability (Step 1), we have that

    α≥ E ∑_π∈Π1( S(û_π) > S^(k(α))(û_π) )/n = E  1( S(û) > S^(k(α))(û) ) = E 1 (p̂≤α).



Step 3: By continuity of the distribution of {S(û_π) }_π∈Π, there are no ties with probability one. 
Therefore,  

    ∑_π∈Π1( S(û_π) ≤ S^(k(α))(û))= k(α)≤ n(1-α)+1

Because

    ∑_π∈Π1( S(û_π) ≤  S^(k(α))(û))+∑_π∈Π1( S(û_π) >  S^(k(α))(û))=n,

we have

    ∑_π∈Π1( S(û_π) >  S^(k(α))(û))≥ nα-1.

The result now follows by similar arguments as in Step 2.







 §.§ Proof of Lemma <ref>



Let X_jt denote the (j,t) entry of the matrix X∈ℝ^T× J. We assume the following conditions hold: (1) E(u_tX_jt)=0 for 1≤ j≤ J. (2) there exist constants c_1,c_2>0 such that  E|X_jtu_t|^2≥ c_1 and  E|X_jtu_t|^3≤ c_2 for any 1≤ j≤ J and 1≤ t≤ T;  (3) for each 1≤ j≤ J, the sequence {X_jtu_t}_t=1^T
is β-mixing and the β-mixing coefficient satisfies that
β(t)≤ a_1exp(-a_2t^τ), where a_1,a_2,τ>0
are constants. (4) there exists a constant c_3>0 such that max_1≤ j≤ J∑_t=1^TX_jt^2u_t^2≤ c_3^2T with probability 1-o(1).  (5) log J=o(T^4τ/(3τ+4)) and w ∈𝒲. (6) There exists a sequence ℓ_T>0 such that (X_t'δ)^2≤ℓ_T X δ_2^2/T,  for all  w+ δ∈𝒲 with probability 1-o(1)
for T_0+1≤ t≤ T  and  (7) ℓ_T B_T → 0
for B_T = M[log(T∨ J)]^(1+τ)/(2τ)T^-1/2. 


Then we claim that under conditions (1)-(5)  listed above:


  (1) There exist a constant M>0 depending only on K and the constants
listed above such that with probability 1-o(1)

    X(ŵ-w)^2_2/ T ≤ B_T = M[log(T∨ J)]^(1+τ)/(2τ)T^-1/2


  (2) Moreover, if (6) and (7) also hold, then

    1/T∑_t=1^T(P̂_t^N-P_t^N)^2=o_P(1)  and P̂_t^N-P_t^N=o_P(1),  for any   T_0+1≤ t≤ T.






The following result is useful in deriving the properties of the ℓ_1-constrained
estimator. 

 Suppose that (1)  E(u_tX_jt)=0 for 1≤ j≤ J. (2) max_1≤ j≤ J,1≤ t≤ TE|X_jtu_t|^3≤ K_1
for a constant K_1>0. (3) min_1≤ j≤ J,1≤ t≤ TE|X_jtu_t|^2≥ K_2
for a constant K_2>0.  (4) For each 1≤ j≤ J, {X_jtu_t}_t=1^T is β-mixing
and the β-mixing coefficients satisfy β(s)≤ D_1exp(-D_2s^τ)
for some constants D_1,D_2,τ>0.  Assume log J=o(T^4τ/(3τ+4)). Then there exists a constant
κ>0 depending only on K_1,K_2,D_1,D_2,τ such
that with probability 1-o(1)

    max_1≤ j≤ J|∑_t=1^TX_jtu_t|<κ[log(T∨ J)]^(1+τ)/(2τ)max_1≤ j≤ J√(∑_t=1^TX_jt^2u_t^2)



Define W_j,t=X_jtu_t. Let m=⌊ [4D_2^-1log(JT)]^1/τ⌋
and k=⌊ T/m⌋. For simplicity, we assume
for now that T/m is an integer. Define 

    H_i={ i,m+i,2m+i,...,(k-1)m+i}    ∀1≤ i≤ m.


By Berbee's coupling (Lemma <ref>), there
exist a sequence of random variables {W̃_j,t}_t∈ H_i such that (1) {W̃_j,t}_t∈ H_i
is independent across t, (2) W̃_j,t has the same distribution
as W_j,t for t∈ H_i and (3) P(⋃_t∈ H_i{W̃_j,t≠ W_j,t})≤ kβ(m). 

By assumption, max_j,tE|X_jtu_t|^3 is uniformly bounded
above and min_j,tE|X_jtu_t|^2 is uniformly bounded away
from zero. It follows, by Theorem 7.4 of <cit.>, that
there exist constants C_0,C_1>0 depending on K_1 and K_2
such that for any 0≤ x≤ C_0k^1/6,

    P(|∑_t∈ H_iW̃_j,t/√(∑_t∈ H_iW̃_j,t^2)|>x)≤ C_1(1-Φ(x)),

where Φ(·) is the cdf of N(0,1). Therefore, for any 0≤ x≤ C_0k^1/6, 

    P(|∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2)|>x)≤ P(|∑_t∈ H_iW̃_j,t/√(∑_t∈ H_iW̃_j,t^2)|>x)+P(⋃_t∈ H_i{W̃_j,t≠ W_j,t})
    ≤ C_1(1-Φ(x))+kβ(m).


The Cauchy-Schwarz inequality implies 

    |∑_t=1^TW_j,t|≤∑_i=1^m|∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2)|√(∑_t∈ H_iW_j,t^2)   ≤√(∑_i=1^m(∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2))^2)×√(∑_i=1^m∑_t∈ H_iW_j,t^2)
        =√(∑_i=1^m(∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2))^2)×√(∑_t=1^TW_j,t^2).


Hence, 

    |∑_t=1^TW_j,t/√(∑_t=1^TW_j,t^2)|≤√(∑_i=1^m(∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2))^2).


It follows that for any 0≤ x≤ C_0k^1/6√(m),

    P(|∑_t=1^TW_j,t/√(∑_t=1^TW_j,t^2)|>x)  ≤ P(√(∑_i=1^m(∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2))^2)>x)
        =P(∑_i=1^m(∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2))^2>x^2)  ≤∑_i=1^mP(|∑_t∈ H_iW_j,t/√(∑_t∈ H_iW_j,t^2)|>x/√(m))
       (i)≤m[C_1(1-Φ(x/√(m)))+kβ(m)]  (ii)≤C_1m√(m/2π)x^-1exp(-x^2/2m)+D_1kmexp(-D_2m^τ)
        <C_1m^3/2x^-1exp(-x^2/2m)+D_1Texp(-D_2m^τ)

where (i) follows by (<ref>) and (ii)
follows by the inequality 1-Φ(a)≤ a^-1ϕ(a) (with ϕ
being the pdf of N(0,1)) and β(m)≤ D_1exp(-D_2m^τ). 

By the union bound, it follows that for any 0≤ x≤ C_0k^1/6√(m),

    P(max_1≤ j≤ J|∑_t=1^TW_j,t/√(∑_t=1^TW_j,t^2)|>x)≤ C_1Jm^3/2x^-1exp(-x^2/2m)+D_1JTexp(-D_2m^τ).


Now we choose x=2√(mlog(Jm^3/2)). Since log J=o(T^4τ/(3τ+4))
and k≍ T/m, it can be very easily verified that x≪ C_0k^1/6√(m)
and the two terms on the right-hand side of the above display tend
to zero. The desired result follows. 

If T/k is not an integer, then we simply add one observation from
{W_j,t}_t=km+1^T to each of H_i for 1≤ i≤ m.
The bound in (<ref>) holds with C_1
large enough. The proof is complete. 


Now we are ready to prove Lemma <ref>.




Let Δ=ŵ-w. Since w_1≤ K, we have Y-Xŵ_2^2≤Y-Xw_2^2.
Notice that Y-Xw=u and Y-Xŵ=u-XΔ. Therefore, u-XΔ_2^2≤u_2^2,
which means XΔ_2^2≤2u'XΔ. Now we observe that

    XΔ_2^2≤2u'XΔ (i)≤2Xu_∞Δ_1 (ii)≤4KXu_∞,

where (i) follows by Hölder's inequality and (ii) follows by
Δ_1≤2K (since ŵ_1≤ K and w_1≤ K).
By Lemma <ref>, there exists a constant
κ>0 such that 

    P(max_1≤ j≤ J|∑_t=1^TX_jtu_t|>κ[log(T∨ J)]^(1+τ)/(2τ)max_1≤ j≤ J√(∑_t=1^TX_jt^2u_t^2))=o(1).


Since P(max_1≤ j≤ J∑_t=1^TX_jt^2u_t^2≤ c_3^2T)→1,
it follows that 

    P(max_1≤ j≤ J|∑_t=1^TX_jtu_t|>κ c_3[log(T∨ J)]^(1+τ)/(2τ)√(T))=o(1).


Part (1) follows by combining (<ref>) and (<ref>).
Part (2) follows by part (1) and ℓ_TB_T=o(1). 














 §.§ Proof of Lemma <ref>
 We borrow results and notations from <cit.>. Following standard notation, we use i instead of j to denote units.
Here are the regularity conditions  from <cit.>.






Suppose that there exists a constant D_0>0 the following conditions
hold:

 (1) max_1≤ t≤ TEF_t_2^4≤ D_0, max_1≤ j≤ Nλ_j_2^4≤ D_0,
max_jtE|u_jt|^8≤ D_0 and E(u_jt)=0.

  (2) max_sN^-1∑_t=1^T|∑_i=1^NE(u_isu_it)|≤ D_0 and max_i∑_j=1^Nmax_1≤ t≤ T|E(u_itu_jt)|≤ D_0.

(3) (NT)^-1∑_s=1^T∑_t=1^T∑_i=1^N∑_j=1^N|E(u_itu_js)|≤ D_0
and max_s,tE|N^-1/2∑_i=1^N[u_isu_it-E(u_isu_it)]|^4≤ D_0.

(4) N^-1∑_i=1^NE T^-1/2∑_t=1^TF_tu_it_2^2≤ D_0.

(5) max_tE  (NT)^-1/2∑_s=1^T∑_i=1^NF_s[u_isu_it-E(u_isu_it)]_2^2≤ D_0.

(6) E (NT)^-1/2∑_t=1^T∑_i=1^NF_tλ_i'u_it_2^2≤ D_0.

Moreover, we assume the following conditions:
(7) for each t, N^-1/2∑_i=1^Nλ_iu_it→^dN(0,Γ_t)
as N→∞, where Γ_t=lim_N→∞N^-1∑_i=1^N∑_j=1^Nλ_iλ_j'E(u_itu_jt).
(8) for each i, T^-1/2∑_t=1^TF_tu_it→^dN(0,Φ_i)
as T→∞, where Φ_i=lim_T→∞T^-1∑_s=1^T∑_t=1^TE(F_tF_s'u_isu_it).
(9) N^-1∑_i=1^Nλ_iλ_i'→Σ_Λ
and T^-1∑_t=1^TF_tF_t'=Σ_F+o_P(1) for some
k× k positive definite matrices Σ_Λ and Σ_F
satisfying that Σ_ΛΣ_F has distinct eigenvalues.


What follows below is the proof of the lemma.  We recall some notations used by <cit.>. Let F=(F_1,…,F_T)'∈^T× k and Λ=(λ_1,…,λ_N)'∈^N× k. Define
H=(Λ'Λ/N)(F'F̃/T)V_NT^-1, where V_NT∈ℝ^k× k
is the diagonal matrix with the largest k eigenvalues of Y^N(Y^N)'/(NT)
on the diagonal and F̃ is the normalized F, namely F̃'F̃/T=I_k.

We start with the first equation in the proof of Theorem 3 in <cit.>
(on page 166):

    λ̂_1'F̂_t-λ_1'F_t=(F̂_t-H'F_t)'H^-1λ_1+F̂_t'(λ̂_1-H^-1λ_1).


The rest of the proof proceeds in two steps. We first recall some
results from <cit.> and then derive the desired
result.

Step 1: recall useful results from <cit.>. By Lemma A.1 of <cit.>, 

    ∑_t=1^TF̂_t-H'F_t_2^2=O_P(T/δ_NT^2),

where δ_NT=min{√(N),√(T)}. By definition, F̂'F̂/T=I_k,
which means 

    ∑_t=1^TF̂_t_2^2=∑_t=1^T trace(F̂_tF̂_t')= trace(F̂'F̂)=kT.


By Theorem 2 of <cit.>, 

    λ̂_1=H^-1λ_1+O_P(max{T^-1/2,N^-1} ).


By the proof of part (i) in Theorem 2 of <cit.>,
H converges in probability to a nonsingular matrix; see page 166
of <cit.>. Hence, H^-1=O_P(1). By
assumption, λ_1_2=O(1). Hence, 

    H^-1λ_1_2=O_P(1).


Step 2: prove the desired result.

Therefore, 

    ∑_t=1^T(λ̂_1'F̂_t-λ_1'F_t)^2    (i)≤2∑_t=1^T[(F̂_t-H'F_t)'H^-1λ_1]^2+2∑_t=1^T[F̂_t'(λ̂_1-H^-1λ_1)]^2
       ≤2∑_t=1^TF̂_t-H'F_t_2^2×H^-1λ_1_2^2+2∑_t=1^TF̂_t_2^2×λ̂_1-H^-1λ_1_2^2
        (ii)=O_P(T/δ_NT^2)× O_P(1)+2kT× O_P(max{T^-1,N^-2} )
        =O_P(T/δ_NT^2),

where (i) follows by (<ref>) and the elementary
inequality of (a+b)^2≤2a^2+2b^2 for any a,b∈ℝ
and (ii) follows by (<ref>), (<ref>),
(<ref>) and (<ref>). Since
n=|Π|=T for moving block permutation, we have 

    1/n∑_t=1^T(λ̂_1'F̂_t-λ_1'F_t)^2=O_P(1/min{N,T}).


Finally, notice that Theorem 3 of <cit.> implies λ̂_1'F̂_t-λ_1'F_t=O_P(1/δ_NT). 
The proof is complete.  





 §.§ Proof of Lemma <ref>
 We recite
conditions from <cit.>. Following standard notation, we use i instead of j to denote units.

Suppose that there exists a constant D_1>0 the following conditions hold:

(1) max_i,tEX_it_2^4≤ D_1, max_tEF_t_2^4≤ D_1, max_iEλ_i_2^4≤ D_1 and max_i,tE|u_it|^8≤ D_1. 

(2) N^-1∑_i=1^N∑_j=1^Nmax_t,s|E(u_itu_js)| ≤ D_1and T^-1∑_s=1^T∑_t=1^Tmax_i,j|E(u_itu_js)|≤ D_1.

(3) (NT)^-1∑_i=1^N∑_j=1^N∑_s=1^T∑_t=1^T|E(u_itu_js)|≤ D_1. 

(4) max_t,sE|N^-1/2∑_i=1^N[u_isu_it-E(u_isu_it)]|^4≤ D_1. 

(5) T^-2N^-1∑_t,s,q,v∑_i,j|cov(u_itu_ts,u_jqu_jv)|≤ D_1

 (6) T^-1N^-2∑_t,s∑_i,j,k,q|cov(u_itu_jt,u_ksu_qs)|≤ D_1.

(7) the largest eigenvalue of E(u_iu_i') is bounded by D_1,
where u_i=(u_i1,...,u_iT)'∈ℝ^T. 


Moreover, the following conditions also hold: (8) u=(u_1,…,u_N) is independent of (X,F,Λ). 
(9) F'F/T=Σ_F+o_P(1) and Λ'Λ/N=Σ_Λ+o_P(1)
for some matrices Σ_F and Σ_Λ. 
(10) N/T is bounded away from zero and infinity.
(11) For X_i=(X_i1,...,X_iT)'∈ℝ^T× k_x
and M_F=I_T-F(F'F)^-1F', we have 

    inf_F: F'F/T=I_k1/NT∑_i=1^NX_i'M_FX_i-1/T[1/N^2∑_i=1^N∑_j=1^NX_i'M_FX_jλ_i'(Λ'Λ/N)^-1λ_j]>0.


What follows below is the proof of the lemma.  We introduce some notations used in <cit.>. Let H=(Λ'Λ/N)(F'F̂/T)V_NT^-1,
where V_NT is the diagonal matrix that contains the k largest
eigenvalues of (NT)^-1∑_i=1^N(Y_i^N-X_iβ̂)(Y_i^N-X_iβ̂)'
with  Y_i^N=(Y_i1^N,Y_i2^N,...,Y_iT^N)'∈ℝ^T.
Let δ_NT=min{√(N),√(T)}. The rest of the proof
proceeds in two steps. We first derive bounds for ∑_t=1^T(û_1t-u_1t)^2
and then prove the pointwise result.

Step 1: derive bounds for ∑_t=1^T(û_1t-u_1t)^2. 

Define Δ_β=β̂-β and Δ_F,t=F̂_t-H'F_t.
Denote Δ_F=(Δ_F,1,...,Δ_F,T)'∈ℝ^T× k.
Notice that F̂-FH=Δ_F. As pointed out on page 1237
of <cit.>, 

    λ̂_1=T^-1F̂'(Y_1^N-X_1β̂)=T^-1F̂'(u_1+Fλ_1-X_1Δ_β).

Notice
that

    |û_1t-u_1t|^2  =|F_t'λ_1-F̂_t'λ̂_1-X_1t'Δ_β|^2
        (i)=|F_t'λ_1-T^-1(H'F_t+Δ_F,t)'F̂'(u_1+Fλ_1-X_1Δ_β)-X_1,t'Δ_β|^2
       ≤[|F_t'(I_k-HF̂'F/T)λ_1|+|T^-1Δ_F,t'F̂'Fλ_1|+|T^-1F̂_t'F̂'(u_1-X_1Δ_β)|+|X_1t'Δ_β|]^2
       ≲[F_t'(I_k-HF̂'F/T)λ_1]^2+[T^-1Δ_F,t'F̂'Fλ_1]^2+[T^-1F̂_t'F̂'(u_1-X_1Δ_β)]^2+[X_1t'Δ_β]^2,

where (i) follows by (<ref>) and F̂_t=H'F_t+Δ_F,t.
Therefore, 

    ∑_t=1^T(û_1t-u_1t)^2≲∑_t=1^T[F_t'(I_k-HF̂'F/T)λ_1]^2+∑_t=1^T[T^-1Δ_F,t'F̂'Fλ_1]^2
           +∑_t=1^T[T^-1F̂_t'F̂'(u_1-X_1Δ_β)]^2+∑_t=1^T[X_1t'Δ_β]^2
        (i)=λ_1'(I_k-HF̂'F/T)'(F'F)(I_k-HF̂'F/T)λ_1
           +T^-2(F̂'Fλ_1)'(Δ_F'Δ_F)(F̂'Fλ_1)+T^-1‖F̂'(u_1-X_1Δ_β)‖_2^2+X_1Δ_β_2^2
        (ii)=O_P(TΔ_β_2^2+Tδ_NT^-4)+O_P(TΔ_β_2^2+Tδ_NT^-2)+O_P(1+Tδ_NT^-4+TΔ_β_2^2)+O_P(TΔ_β_2^2)
        =O_P(1+TΔ_β_2^2+Tδ_NT^-2),

where (i) follows by ∑_t=1^TF̂_tF̂_t'=F̂'F̂=TI_k
and (ii) follows by Lemma <ref>, together
with F=O_P(√(T)), λ_1=O(1) and F̂=O_P(√(T)).
Since N≍ T, Theorem 1 of <cit.> implies Δ_β_2=O_P(1/√(NT))=O_P(T^-1).
Therefore, the above display implies 

    ∑_t=1^T(û_1t-u_1t)^2=O_P(1).


Step 2: show the pointwise result.

By (<ref>), we have 

    |û_1t-u_1t|    ≤|F_t'(I_k-HF̂'F/T)λ_1|+|T^-1Δ_F,t'F̂'Fλ_1|+|T^-1F̂_t'F̂'(u_1-X_1Δ_β)|+|X_1t'Δ_β|
        (i)≤F_t_2·λ_1_2· O_P(Δ_β_2+δ_NT^-2)+O_P(TΔ_β_2+Tδ_NT^-2)· T^-1Fλ_1_2
           +T^-1F̂_t_2· O_P(√(T)+Tδ_NT^-2+TΔ_β_2)+X_1t_2·Δ_β_2 (ii)≤O_P(T^-1/2),

where (i) follows by I_k-HF̂'F/T=O_P(Δ_β_2+δ_NT^-2), F̂Δ_F,t=O_P(TΔ_β_2+Tδ_NT^-2), 
and F̂'(u_1-X_1Δ_β)_2 =O_P(√(T)+Tδ_NT^-2+TΔ_β_2)
(due to Lemma <ref>), whereas (ii) follows
by F̂_t_2=O_P(1) (Lemma <ref>),
X_1t_2=O_P(1), F_t_2=O_P(1), λ_1=O(1),
Δ_β_2=O_P(T^-1) and Fλ_1_2=O_P(√(T)). 


Suppose that the assumption of Lemma
<ref> holds. Let δ_NT, H,
Δ_F and u_1 be defined as in the proof of Lemma <ref>.
Then (1) I_k-HF̂'F/T=O_P(Δ_β_2+δ_NT^-2);
(2) Δ_F'Δ_F=O_P(TΔ_β_2^2+Tδ_NT^-2);
(3) ‖F̂'(u_1-X_1Δ_β)‖ =O_P(√(T)+Tδ_NT^-2+TΔ_β_2);
(4) X_1Δ_β_2=O_P(√(T)Δ_β_2);
(5) F̂Δ_F,t_2=O_P(TΔ_β_2+Tδ_NT^-2);
(6) F̂_t_2=O_P(1) for 1≤ t≤ T. 


Proof of part (1). Lemma A.7(i) in <cit.>
implies HH' converges in probability to a nonsingular matrix. Hence,

    H=O_P(1)   and   H^-1=O_P(1).


Notice that 

    I_k-HF̂'F/T (i)=I_k-H(FH+Δ_F)'F/T     =I_k-(HH')(F'F/T)-HΔ_F'F/T
        (ii)=O_P(Δ_β_2)+O_P(δ_NT^-2)-HΔ_F'F/T
        (iii)=O_P(Δ_β_2)+O_P(δ_NT^-2),

where (i) holds by F̂=FH+Δ_F, (ii) holds by I_k-(HH')(F'F/T)=O_P(Δ_β_2)+O_P(δ_NT^-2)
(due to Lemma A.7(i) in <cit.>) and (iii) holds
by (<ref>) and Δ_F'F/T=O_P(Δ_β_2)+O_P(δ_NT^-2)
(due to Lemma A.3(i) in <cit.>). This proves part (1).

Proof of part (2). Part (2) follows by Proposition A.1 of
<cit.>: 

    T^-1Δ_F'Δ_F=O_P(Δ_β_2^2)+O_P(δ_NT^-2).


Proof of part (3). To see part (3), first observe that the
independence between u_1 and F implies that 

    E(F'u_1^2| F)≤∑_t=1^TE(F_t'F_tu_1t^2| F)=∑_t=1^TF_t'F_tE(u_1t^2).

It follows that 

    E(F'u_1^2)≤∑_t=1^TE(F_t'F_t)E(u_1t^2) (i)≲T∑_t=1^TE(u_1t^2)=O(T),

where (i) holds by the uniform boundedness of E(F_t'F_t). This means that 

    F'u_1_2=O_P(√(T)).


Notice that 

    ‖F̂'(u_1-X_1Δ_β)‖_2   ≤‖ H'F'u_1‖_2 +‖(F̂-FH)'u_1‖ +F̂·X_1·Δ_β_2
        (i)=‖ H'F'u_1‖_2 +(O_P(T^1/2Δ_β_2)+O_P(Tδ_NT^-2))+O_P(TΔ_β_2)
        (ii)=O_P(√(T))+(O_P(T^1/2Δ_β_2)+O_P(Tδ_NT^-2))+O_P(TΔ_β_2),

where (i) follows by (F̂-FH)'u_1/T=O_P(T^-1/2Δ_β_2)+O_P(δ_NT^-2)
(due to Lemma A.4 in <cit.>) and the fact that F̂=O(√(T))
and X_1=O_P(√(T)) (see the beginning of Appendix A
in <cit.>), whereas (ii) follows by (<ref>)
and (<ref>). We have proved part (3). 

Proof of part (4). We notice that X_1=O_P(√(T));
see the beginning of Appendix A in <cit.>. Part (4)
follows by X_1Δ_β≤X_1·Δ_β_2.

Proof of part (5). Notice that 

    F̂Δ_F,t_2/T≤F̂Δ_F/T (i)≤O_P(Δ_β_2)+O_P(δ_NT^-2),

where (i) follows by Lemma A.3(ii) in <cit.>. We
have proved part (5). 

Proof of part (6). Notice that 

    T^-1Δ_F,t_2^2≤ T^-1Δ_F'Δ_F  = T^-1F̂'Δ_F-T^-1H'F'Δ_F (i)= O_P(Δ_β_2)+O_P(δ_NT^-2),

where (i) follows by Lemma A.3(i)-(ii) of <cit.>. By Theorem 1 of <cit.> and by the assumption of N ≍ T, we have that Δ_F,t_2=O_P(1). By F̂_t_2≤H'F_t_2+Δ_F,t_2, F_t=O_P(1)
and H=O_P(1), we can see that F̂_t_2=O_P(1). The
proof is complete. 






 §.§ Proof of Lemma <ref>


We start with the assumptions. Recall N=J+1.   Assume that (1) {u_j}_j=1^N
is independent across j conditional on M, (2) max_1≤ j≤ NT^-1∑_t=1^TE(|u_jt|^2κ_1| M)=O_P(1)
for some constant κ_1>1, (3) N^-1∑_j=1^NE(u_ju_j'| M)=O_P(1) and  (4)  there exists a sequence ℓ_T >0 such that

 ℓ_T (NT)^-1K√(N∨(N^1/κ_1Tlog N))=o(1)  and  with probability 1-o(1), T^-1∑_t=1^T(M̂_1t-M_1t)^2≤ℓ_T (NT)^-1∑_t=1^T∑_j=1^N(M̂_jt-M_jt)^2 and (M̂_1t-M_1t)^2≤ℓ_T (NT)^-1∑_t=1^T∑_j=1^N(M̂_jt-M_jt)^2 for T_0+1≤ t≤ T. 



We now prove Lemma <ref>. Define
Δ=M̂-M. Let Y∈ℝ^T× N be the matrix
whose (t,j) entry is Y_jt^N. For  (j,t), define
the matrix Q_jt∈ℝ^N× T by Q_is=1{(i,s)=(j,t)},
i.e, a matrix of zeros except that the (j,t) entry is one. Then
we can write the model as 

    Y_jt^N= trace(Q_jt'M)+u_jt     for   1≤ j≤ N, 1≤ t≤ T.


Notice that the estimator M̂ satisfies 

    M̂_*≤ K

and

    ∑_t=1^T∑_j=1^N(Y_jt^N- trace(Q_jt'M̂))^2≤∑_t=1^T∑_j=1^N(Y_jt^N- trace(Q_jt'M))^2.


Plugging (<ref>) into the above inequality
and rearranging terms, we obtain 

    ∑_t=1^T∑_j=1^N( trace(Q_jt'Δ))^2≤2∑_t=1^T∑_j=1^Nu_jt trace(Q_jt'Δ)     =2 trace([∑_t=1^T∑_j=1^Nu_jtQ_jt]'Δ)
        (i)=2 trace(u'Δ)
        (ii)≤2u·Δ_*
        (iii)≤4Ku,

where (i) follows by ∑_t=1^T∑_j=1^Nu_jtQ_jt=u, (ii) follows
by the trace duality property (see e.g., <cit.>,
<cit.> and <cit.>) and  (iii)
follows by the fact that M̂_*≤ K and M_*≤ K.

To bound u, we apply Lemma <ref>. Note that the conditions of Lemma <ref>
are satisfied by our assumption. Therefore, E(u| M)=O_P(√(N∨(N^1/κ_1Tlog N))).
This and (<ref>) imply that 

    (NT)^-1∑_t=1^T∑_j=1^N( trace(Q_jt'Δ))^2=O_P((NT)^-1K√(N∨(N^1/κ_1Tlog N))).

The desired result follows by Assumption (4) listed at the beginning of the proof. 




Suppose that the following conditions
hold:

(i) {u_j}_j=1^N is independent across j conditional
on M. 

(ii)  max_1≤ j≤ NT^-1∑_t=1^TE(|u_jt|^2κ_1| M)=O_P(1)
for some constant κ_1>1. 

(iii)  N^-1∑_j=1^NE(u_ju_j'| M)=O_P(1). 

Then E(u| M)=O_P(√(N∨ (N^1/κ_1Tlog N))).


Recall the elementary inequality |T^-1∑_t=1^Ta_t|≤[T^-1∑_t=1^T|a_t|^κ]^1/κ
for any κ>1 (due to Liapunov's inequality). It follows that
T^-1∑_t=1^Tu_jt^2≤[T^-1∑_t=1^T|u_jt|^2κ_1]^1/κ_1,
which means 

    (∑_t=1^Tu_jt^2)^κ_1≤ T^κ_1-1∑_t=1^T|u_jt|^2κ_1.


Hence, 

    E([max_1≤ j≤ N∑_t=1^Tu_jt^2]| M)     (i)≤{ E[max_1≤ j≤ N(∑_t=1^Tu_jt^2)^κ_1| M]} ^1/κ_1
       ≤{ E[∑_i=1^N(∑_t=1^Tu_jt^2)^κ_1| M]} ^1/κ_1
        (ii)≤{ E[T^κ_1-1∑_j=1^N∑_t=1^T|u_jt|^2κ_1| M]} ^1/κ_1
       ≤{[NT^κ_1max_1≤ i≤ NT^-1∑_t=1^TE(|u_jt|^2κ_1| M)]} ^1/κ_1
        (iii)={[NT^κ_1O_P(1)]} ^1/κ_1=O_P(N^1/κ_1T),

where (i) follows by Liapunov's inequality, (ii) follows by (<ref>)
and (iii) follows by the assumption that max_1≤ j≤ NT^-1∑_t=1^TE(|u_jt|^2κ_1| M)=O_P(1).
Therefore, it follows, by Theorem 5.48 and Remark 5.49 of <cit.>,
that 

    E(u| M)≤√(E(u^2| M))   ≤E(u'u| M)/N^1/2√(N)
           +O(√(O(N^1/κ_1T)logmin(O(N^1/κ_1T),N)))
        (i)≤O_P(√(N))+O(√(N^1/κ_1Tlog N)),

where (i) holds by the assumption of E(u'u| M)/N=N^-1∑_j=1^NE(u_ju_j'| M)=O_P(1).
The proof is complete. 













 §.§ Proof of Lemma  <ref>



By the analysis on page 215-216 of <cit.> (leading
to Equation (8.2.29) therein), we have that ρ̂-ρ=o_P(1).
Hence, 

    ∑_t=K+1^T(û_t-u_t)^2    =∑_t=K+1^T(y_t'(ρ-ρ̂))^2 =(ρ̂-ρ)'(∑_t=K+1^Ty_ty_t')(ρ̂-ρ) ≤ρ̂-ρ_2^2×‖∑_t=K+1^T y_ty_t'‖.


The analysis on page 215 of <cit.> (leading to Equation
(8.2.26) therein) implies that 
    T^-1∑_t=K+1^T y_ty_t'=E(y_ty_t')+o_P(1),

which means ‖∑_t=K+1^T y_ty_t'‖=O_P(T). Since ρ̂-ρ=o_P(1),
the above display implies that 

    ∑_t=K+1^T(û_t-u_t)^2=o_P(T).


Since û_t-u_t=y_t'(ρ-ρ̂), the pointwise
consistency follows by ρ̂-ρ=o_P(1) and the fact that
y_t=O_P(1) for T_0+1≤ t≤ T (due to the stationarity
property of u_t). 




 §.§ Proof of Lemma <ref>



By assumption, max_K+1 ≤ t ≤ T |P̂_t^N-P_t^N|≤ℓ_T ρ̂- ρ. Therefore, 

    1/T∑_t=K+1^T(P̂^N_t-P^N_t)^2≤ℓ_T^2 ρ̂- ρ

and

    max_T_0+1 ≤ t ≤ T |P̂^N_t-P^N_t|≤ℓ_T ρ̂- ρ .

Since ℓ_T ρ̂- ρ=o_P(1), the desired result follows.





 §.§ Proof of Lemma <ref>



We first derive the following result that is useful in proving Lemma <ref>.






Recall ε_t=x_t'ρ+u_t, where   ρ=(ρ_1,ρ_2,...,ρ_K)'∈ℝ^K
and x_t=(ε_t-1,ε_t-2,...,ε_t-K)'∈ℝ^K.
Suppose that  the following hold: 
(1) {u_t}_t=1^T is an iid sequence with E(u_1^4)
uniformly bounded. 
(2) the roots of 1-∑_j=1^Kρ_jL^j=0 are uniformly
bounded away from the unit circle.

Then we have
(i) (T-K)^-1∑_t=K+1^Tu_t^2=O_P(1);
(ii) (T-K)^-1∑_t=K+1^Tx_tu_t=o_P(1);
(iii) (T-K)^-1∑_t=K+1^Tx_t^2=O_P(1).
(iv) There exists a constant λ_0>0 such that the smallest
eigenvalue of (T-K)^-1∑_t=K+1^Tx_tx_t' is bounded
below by λ_0 with probability approaching one. 


Proof of part (i). Part (i) follows by the law of large numbers;
see e.g., Theorem 3.1 of <cit.>. 

Proof of part (ii). Let ℱ_t be the σ-algebra
generated by { u_s: s≤ t }. First notice that {x_tu_t}_t=K+1^T is a
martingale difference sequence with respect to the filtration {ℱ_t}. Since ε_t is a stationary
process, we have that Ex_tu_t^2=  ∑_j=1^KE(ε_t-j^2u_t^2)= ∑_j=1^KE(ε_t-j^2)E(u_t^2)
is uniformly bounded. Hence, part (ii) follows by Exercise
3.77 of <cit.>. 

Proof of part (iii). To see part (iii), notice that x_t^2 = x_t'x_t=∑_j=1^Kε_t-j^2.
By the analysis on page 215 of <cit.>, for each 1≤ j≤ K,
(T-K)^-1∑_t=K+1^Tε_t-j^2=E(ε_t-j^2)+o_P(1).
Thus, part (iii) follows by

    (T-K)^-1∑_t=K+1^Tx_t^2=(T-K)^-1∑_j=1^K∑_t=K+1^Tε_t-j^2=K(E(ε_t^2)+o_P(1)).

 

Proof of part (iv). Similarly, the analysis on page 215 of
<cit.> implies that 
    (T-K)^-1∑_t=K+1^Tx_tx_t'=o_P(1)+Ex_tx_t'.

By Proposition 5.1.1 of  <cit.>, E(x_tx_t')
has eigenvalues bounded away from zero. Part (iv) follows. 




Now we are ready to prove Lemma <ref>.





 Define δ_t=ε̂_t-ε_t, Δ_t=x̂_t-x_t,
ũ_t=u_t+δ_t-Δ_t'ρ and a_t=ũ_t-u_t.
Notice that 

    ε̂_t=δ_t+ε_t=δ_t+x_t'ρ+u_t=δ_t+(x̂_t-Δ_t)'ρ+u_t=x̂_t'ρ+ũ_t.


Therefore, 

    ρ̂=(∑_t=K+1^Tx̂_tx̂_t')^-1(∑_t=K+1^Tx̂_tε̂_t)     =(∑_t=K+1^Tx̂_tx̂_t')^-1(∑_t=K+1^Tx̂_t(x̂_t'ρ+ũ_t))
        =ρ+(∑_t=K+1^Tx̂_tx̂_t')^-1(∑_t=K+1^Tx̂_tũ_t).


The rest of the proof proceeds in three steps. First two steps show
that (T-K)^-1∑_t=K+1^Tx̂_tx̂_t'
is well-behaved and (T-K)^-1∑_t=K+1^Tx̂_tũ_t=o_P(1).
This would imply ρ̂=ρ+o_P(1). In the third step, we
derive the final result. 

Step 1: show that [(T-K)^-1∑_t=K+1^Tx̂_tx̂_t']^-1=O_P(1). 

It is not hard to see that Δ_t^2=∑_s=t-1^t-Kδ_s^2.
Therefore, 

    ∑_t=K+1^TΔ_t^2=∑_t=K+1^T∑_s=t-1^t-Kδ_s^2≤ K∑_t=1^Tδ_t^2 (i)=o_P(T),

where (i) follows by the assumption of T^-1∑_t=1^Tδ_t^2=o_P(1).
Notice that 

    ‖∑_t=K+1^T(x̂_tx̂_t'-x_tx_t')‖    =‖∑_t=K+1^T(x_tΔ_t'+Δ_tx_t'+Δ_tΔ_t')‖
       ≤2∑_t=K+1^Tx_t·Δ_t+∑_t=K+1^TΔ_t^2
       ≤2√((∑_t=K+1^Tx_t^2)(∑_t=K+1^TΔ_t^2))+∑_t=K+1^TΔ_t^2 (i)=o_P(T),

where (i) follows by (<ref>) and Lemma <ref>.
Thus, 

    ‖1/T-K∑_t=K+1^T(x̂_tx̂_t'-x_tx_t')‖ =o_P(1).


By Lemma <ref>, the smallest eigenvalue of
(T-K)^-1∑_t=K+1^Tx_tx_t' is bounded below by a positive constant with probability approaching one. It follows that 

    [(T-K)^-1∑_t=K+1^Tx̂_tx̂_t']^-1=O_P(1).


Step 2: show that (T-K)^-1∑_t=K+1^Tx̂_tũ_t=o_P(1). 

By Lemma <ref>, we have 

    (T-K)^-1∑_t=K+1^Tx_tu_t=o_P(1).


Notice that 

    ‖1/T-K∑_t=K+1^T(x̂_tũ_t-x_tu_t)‖    =‖1/T-K∑_t=K+1^T(Δ_tu_t+x_ta_t+Δ_ta_t)‖
       ≤1/T-K∑_t=K+1^T(‖Δ_tu_t‖ +‖ x_ta_t‖ +‖Δ_ta_t‖)
       ≤√((1/T-K∑_t=K+1^TΔ_t^2)(1/T-K∑_t=K+1^Tu_t^2))
           +√((1/T-K∑_t=K+1^Tx_t^2)(1/T-K∑_t=K+1^Ta_t^2))
           +√((1/T-K∑_t=K+1^TΔ_t^2)(1/T-K∑_t=K+1^Ta_t^2)).


We observe that

    ∑_t=K+1^Ta_t^2=∑_t=K+1^T(δ_t-Δ_t'ρ)^2   ≤2∑_t=K+1^Tδ_t^2+2∑_t=K+1^T(Δ_t'ρ)^2
       ≤2∑_t=1^Tδ_t^2+2ρ^2∑_t=K+1^TΔ_t^2 (i)=O_P(T),

where (i) follows by (<ref>) and the assumption
of T^-1∑_t=1^Tδ_t^2=o_P(1). Combining (<ref>) with (<ref>) and (<ref>),
we obtain 

    ‖1/T-K∑_t=K+1^T(x̂_tũ_t-x_tu_t)‖
    ≤√(o_P(1)(1/T-K∑_t=K+1^Tu_t^2))+√((1/T-K∑_t=K+1^Tx_t^2)o_P(1))+√(o_P(1)× o_P(1)) (i)=o_P(1),

 where (i) follows by Lemma <ref>. Now we
combine (<ref>) and (<ref>),
obtaining 

    (T-K)^-1∑_t=K+1^Tx̂_tũ_t=o_P(1).


By (<ref>) together with (<ref>)
and (<ref>), it follows that 

    ρ̂-ρ=o_P(1).


Step 3: show the desired result.

Recall that û_t=ε̂_t-x̂_t'ρ̂.
Hence, 

    û_t-u_t=(ε̂_t-x̂_t'ρ̂)-u_t (i)=(x̂_t'(ρ-ρ̂)+ũ_t)-u_t=x̂_t'(ρ-ρ̂)+a_t,

where (i) follows by (<ref>). Therefore, we
have 

    ∑_t=K+1^T(û_t-u_t)^2    =∑_t=K+1^T(x̂_t'(ρ-ρ̂)+a_t)^2
       ≤2∑_t=K+1^T(x̂_t'(ρ̂-ρ))^2+2∑_t=K+1^Ta_t^2
       ≤2ρ̂-ρ^2∑_t=K+1^Tx̂_t^2+2∑_t=K+1^Ta_t^2
        =2ρ̂-ρ^2(∑_t=K+1^T trace(x_tx_t')+∑_t=K+1^T trace(x̂_tx̂_t'-x_tx_t'))+2∑_t=K+1^Ta_t^2
        (i)≤o_P(1)×(O_P(T)+o_P(T))+o_P(T)=o_P(T),

where (i) follows by (<ref>), (<ref>),
(<ref>) and Lemma <ref>. 

To see the pointwise result, we notice that by (<ref>)
and (<ref>), it suffices to verify that a_t=o_P(1)
and x̂_t=O_P(1) for T_0+1≤ t≤ T. 

Since x̂_t-x_t=(δ_t-1,δ_t-2,...,δ_t-K)',
the assumption of pointwise convergence of ε̂_t
(i.e., δ_t=o_P(1) for T_0+1-K≤ t≤ T) implies
that x̂_t-x_t=o_P(1) for T_0+1≤ t≤ T. Since
x_t=O_P(1) due to the stationarity condition, we have x̂_t=O_P(1)
for T_0+1≤ t≤ T. 

Since both δ_t and Δ_t are both o_P(1) for
T_0+1≤ t≤ T, so is a_t=δ_t-Δ_t'ρ.
Hence, we have proved the pointwise result. The proof is complete. 






 §.§ Proof of Lemma  <ref>


	Fix an arbitrary η>0. Define a_η=inf_β-β_*_2≥η(L(β)-L(β_*))/3.
	By the compactness of  and the uniqueness of the minimum of
	L(·), we have a_η>0. 
	
	(Otherwise, one can find a sequence {β_k}_k=1^∞
	with β_k-β_*_2≥η for all k≥1 with
	L(β_k)→ L(β_*). By compactness of 
	implies that some subsequence of β_k converges to a point
	β_**∈. Clearly, β_**-β_*_2≥η.
	The continuity of L(·) implies L(β_**)=L(β_*).
	This contradicts the uniqueness of the minimum of L(·).)
	
	Define the event 
    ={sup_β|L̂(;β)-L(β)|≤ a_η}⋂{max_H∈sup_β|L̂(_H;β)-L(β)|≤ a_η}.

	By the assumption, we know P()=1-o(1). 
	
	Notice that on the event , 
	
    L(())-L(β_*)≤ a_η+L̂(;())-L(β_*)≤2a_η+L̂(;())-L̂(;β_*)≤2a_η.

	
	It follows by the definition of a_η that ()-β_*_2≤η
	on the event . Thus, P(-β_*_2≤η)≥ P()=1-o(1).
	Since η>0 is arbitrary, we have ()-β_*_2=o_P(1). 
	
	By the same analysis, we have that on the event , (_H)-β_*_2≤η
	for all H∈. Thus, on the event , max_H∈(_H)-β_*_2≤η.
	We have that max_H∈(_H)-β_*_2=o_P(1).
	The desired result follows. 



 §.§ Proof of Lemma  <ref>

For notational simplicity, we write =() and  _H=(_H). Define the event =_1⋂_2, where _1={min_v_0≤ mv' v/v_2^2≥κ_0}⋂{_0≤ s/2}
	and 
	
    _2=⋂_H∈({_H-_∞≤ c_T, _H-_∞≤ c_T}⋂{max_H∈_H_0≤ s/2}).

	
	By assumption, P()≥1-γ_1,T-γ_2,T-γ_3,T.
	The rest of the argument are statements on the event . 
	
	Fix H∈, let Δ=_H-. Define ξ=-
	and ξ_H=_H-_H. Since _1≤ K,
	we have 
	
    ξ_H-ξ_∞≤_H-_∞+_H-_∞_1≤ c_T(1+K).

	
	When Δ=0, the result clearly holds. Now we consider the case
	with Δ≠0. 
	
	Step 1: show that on the event , 0≤λ^2Δ'Δ-λΔ'≤ c_TK^2+2c_TK
	for any λ∈[0,1]. 
	
	Recall that (β)=β'β-2'β+T^-1∑_t=1^TY_t^2. Since the term T^-1∑_t=1^TY_t^2 does not affect the minimizer, we modify  by dropping this term. With a slight abuse of notation, we still use the symbol  (β)=β'β-2'β and  _H(β)=β'_Hβ-2_H'β. Therefore,  
	
    (β)-_H(β)=β'(-_H)β-2(-_H)'β.

	
	Since sup_β∈β_1≤ K, we have that on
	the event , 
	
    sup_β∈|(β)-_H(β)|≤ c_TK^2+2c_TK.

	
	Let β̅=+λΔ, where λ∈[0,1].
	Then clearly, β̅=λ_H+(1-λ).
	By definition of , we have that ()≤(β̅),
	which means that 
	
    λ^2Δ'Δ-2λξ'Δ≥0.

	
	Clearly, _H(_H)≤_H(β̅) for any λ∈[0,1]. Hence, _H(+Δ)≤_H(+λΔ) for any λ∈(0,1).	By _H(β)=β'_Hβ-2_H'β, this simplifies to (1+λ)Δ' _HΔ≤ 2ξ_H'Δ. Since  λ can be arbitrarily close to one, this means	Δ' _HΔ≤ξ_H'Δ. It follows that for any λ∈[0,1], we have
	
    λ^2Δ'_HΔ≤2λξ_H'Δ.

	
	Notice that 
	
    0≤2λξ_H'Δ-λ^2Δ'_HΔ=2λξ'Δ-λ^2Δ'Δ+2λ(ξ_H-ξ)'Δ+λ^2Δ'(-_H)Δ
    ≤2λξ'Δ-λ^2Δ'Δ+2λξ_H-ξ_∞Δ_1+λ^2-_H_∞Δ_1^2.

	
	It follows, by (<ref>) and Δ_1≤_H_1+_1≤2K,
	that
	
    λ^2Δ'ΣΔ-2λξ'Δ≤2λξ_H-ξ_∞Δ_1+λ^2-_H_∞Δ_1^2≤4c_T(1+K)K+4c_TK^2.

	
	Since (<ref>) and (<ref>)
	hold for any λ∈[0,1], we have that 
	
    0≤λ^2Δ'Δ-2λξ'Δ≤4c_TK(2K+1)    ∀λ∈[0,1].

	
	Step 2: show the desired result.
	
	Since 0≤λ^2Δ'Δ-2λξ'Δ for
	any λ∈(0,1), we have that ξ'Δ≤λΔ'Δ/2
	for any λ∈(0,1). Thus, 
	
    ξ'Δ≤0.

	
	Hence, by the second inequality in (<ref>),
	for any λ∈[0,1],
	
    λ^2Δ'Δ≤λ^2Δ'Δ-2λξ'Δ≤4c_TK(2K+1).

	
	Now we take λ=1, which implies 
	
    Δ'Δ≤4c_TK(2K+1).

	
	Since Δ_0≤_H_0+_0≤ s
	and Δ_1≤√(Δ_0)Δ_2, it follows
	that 
	
    4c_TK(2K+1)≥Δ'Δ≥κ_1Δ_2^2≥κ_1s^-1Δ_1^2.

	
	Hence, Δ_1≤2√(κ_1sc_TK(2K+1)) and 
	
    |(Y_t-X_t')-(Y_t-X_t'_H)|=|X_t'Δ|≤X_t_∞Δ_1≤2κ_2√(κ_1sc_TK(2K+1)).

	
	On the event , the above bound holds for all H∈.
	The desired result follows by P()≥1-γ_1,T-γ_2,T-γ_3,T. 





 §.§ Proof of Lemma <ref>

	For notational simplicity, we write  instead of _λ
	and c_T=|H|. Let =X'X, =X̃'X̃,
	=X'u and =X̃'ũ. We work on the event
	on which κ_1T≤λ_min()≤λ_max()≤κ_2T,
	-μ̃_2≤κ_3√(Jc_T) and -Ω̃≤κ_4(c_T+J).
	Let q_1≥ q_2≥⋯≥ q_J denote the eigenvalues
	of . 
	
	Step 1: show inconsistency.
	
	Notice that 
	
    =(+λ I)^-1X'(Xβ+u)
        =(+λ I)^-1β+(+λ I)^-1X'u
        =β+[(+λ I)^-1-I]β+(+λ I)^-1X'u
        =β+(+λ I)^-1[-(+λ I)]β+(+λ I)^-1X'u
        =β-λ(+λ I)^-1β+(+λ I)^-1X'u.

	
	Therefore, 
	
	
    E(X(-β)_2^2| X)
        =λ X(+λ I)^-1β_2^2+E(X(+λ I)^-1X'u_2^2| X)
        =λ^2β'(+λ I)^-1(+λ I)^-1β+E(u'X(+λ I)^-1(+λ I)^-1X'u| X)
        =λ^2β'(+λ I)^-1(+λ I)^-1β+E([(+λ I)^-1(+λ I)^-1X'uu'X]| X)
        =λ^2β'(+λ I)^-1(+λ I)^-1β+[(+λ I)^-1(+λ I)^-1].

	
	Notice that 
	
    [(+λ I)^-1(+λ I)^-1]=∑_i=1^Jq_i^2(q_i+λ)^-2≳T^2J/(T+λ)^2

	and 
	
    λ^2β'(+λ I)^-1(+λ I)^-1β≥λ^2β_2^2min_1≤ i≤ Jq_i(q_i+λ)^-2≳λ^2T/(T+λ)^2.

	
	Then we have 
	
    E[T^-1X(-β)_2^2| X]≳TJ+λ^2/(T+λ)^2≳TJ/(T+λ)^2+λ^2/(T+λ)^2.

	
	Since λ≍ T, we have E(X(-β)_2^2)≳ T. 
	
	Step 2: show stability.
	
	Similar to (<ref>), we notice that the perturbed estimator
	would be
	
    =β-λ(+λ I)^-1β+(+λ I)^-1.

	
	Then we bound 
	
    X_t'(-)     =λ X_t'[(+λ I)^-1-(+λ I)^-1]β+X_t'[(+λ I)^-1-(+λ I)^-1]
        =λ X_t'[(+λ I)^-1-(+λ I)^-1]β+X_t'[(+λ I)^-1(-)+((+λ I)^-1-(+λ I)^-1)]
        =X_t'[(+λ I)^-1-(+λ I)^-1](λβ-)+X_t'(+λ I)^-1(-).

	Now we notice that 
	
    (+λ I)^-1-(+λ I)^-1    =‖ (+λ I)^-1[(+λ I)-(+λ I)](+λ I)^-1‖
        =(+λ I)^-1(-)(+λ I)^-1=O_P((J+c_T)/(T+λ)^2)

	and λβ-_2≤λ+_2=O_P(λ+√(JT)).
	This means that 
	
    |X_t'[(+λ I)^-1-(+λ I)^-1](λβ-)|=O_P((J+c_T)/(T+λ)^2(λ+√(JT))√(J)).

	
	Moreover, we also have
	
    |X_t'(+λ I)^-1(-)|≤ O_P(√(J)1/μ+λ√(Jc_T))=O_P(J√(c_T)/T+λ).

	
	It follows that 
	
    max_1≤ t≤ T|X_t'(-)|=O_P((J+c_T)/(T+λ)^2(λ+√(JT))√(J)+J√(c_T)/T+λ).

	
	Since J≪ T and T≳λ≫ J^3/2, we have max_1≤ t≤ T|X_t'(-)|=o_P(1). 





§ TABLES AND FIGURES APPENDIX






























