Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
83,584 characters · 11 sections · 48 citation commands
Matrix Completion When Missing Is Not at Random and Its Applications in Causal Panel Data Models
\def\spacingset#1{ {#1}} \spacingset{1}
{\it Keywords:} Matrix completion; Missing not at random (MNAR); Weak signal-to-noise ratio; Multiple treatments; Tick size pilot program
\spacingset{1.4}
{8pt} {8pt} \setlength\intextsep{8pt} {4pt}
The problem of noisy matrix completion in which we are interested in reconstructing a low-rank matrix from partial and noisy observations of its entries arises naturally in numerous applications. It has attracted a considerable amount of attention in recent years, and a lot of impressive results have been obtained from both statistical and computational perspectives. See, e.g., candes2010matrix,mazumder:2010, koltchinskii:2011, negahban2012restricted,chen2019inference, chen2020noisy, jin2021factor, xia2021statistical,bhattacharya2022matrix among many others. A common and crucial premise underlying these developments is that observations of the entries are missing at random. Although this is a reasonable assumption for some applications, it could be problematic for many others. In the past several years, there has been growing interest to investigate how to deal with situations where missing is not at random and to what extent the techniques and insights that are initially developed assuming missing at random can be extended to these cases. See, e.g. agarwal2020synthetic, agarwal2021causal, athey2021matrix, bai2021matrix, chernozhukov2021inference,cahan2023factor, xiong2023large among others.
This fruitful line of research is largely inspired by the development of synthetic control methods in causal inference. See, e.g., abadie2003economic, abadie2010synthetic,abadie2021using. The close connection between noisy matrix completion and synthetic control methods for panel data was first made formal by athey2021matrix who showed that powerful matrix completion techniques such as nuclear norm regularization can be very useful for many causal panel data models where missing is not at random. It also helps bring together two complementary perspectives of noisy matrix completion: one focuses on statistical inferences assuming a strong factor structure and the other aims at recovery guarantees with minimum signal strength requirement. The main objective of this work is to further bridge the gap between these two schools of ideas and develop a general and flexible inferential framework for matrix completion when missing is not at random and without the requirement of strong factors.
In particular, we shall follow athey2021matrix and investigate how the technique of nuclear norm regularization can be used to infer individual treatment effects under a variety of missing mechanisms. One of the key observations to our development is the fact that if the number of missing entries is sufficiently small when compared to the panel size, then they can be estimated well even when missing is not at random. For more general missing patterns with an arbitrary proportion of missingness, we can judicially divide the missing entries into smaller groups and leverage this fact by applying the nuclear norm regularization to a submatrix with a small number of missing entries. This is where our approach differs from that of athey2021matrix who suggest applying the nuclear norm regularized estimation to the full matrix. We shall show that subgrouping is essential in producing more accurate estimates and more efficient inferences about individual treatment effects. It is worth noting that it is computationally more efficient to estimate all missing entries together, as suggested by athey2021matrix. But estimating too many missing entries simultaneously can be statistically suboptimal. In a way, our results suggest how to trade-off between the computational cost and statistical efficiency.
Our proposal of subgrouping is similar in spirit to the approach taken by agarwal2021causal who suggested estimating the missing entries one at a time. For estimating a single missing entry, they propose a matching scheme that constructs multiple “synthetic" neighbors and averages the observed outcomes associated with each synthetic neighbor. Separating the observations into different sets of neighbors, however, could lead to a loss in efficiency. For example, when estimating the mean of an $N\times N$ matrix with one missing entry, the estimation error of the approach from agarwal2021causal for the missing entry converges at the rate of $N^{-1/4}$, which is far slower than the rate of $N^{-1/2}$ attained by our method.
Furthermore, we show that, with appropriate debiasing, our proposed estimate is asymptotically normal even with fairly weak signals. More specifically, the asymptotic normality holds if $\psi_{\min}^2\gg \sigma^2N$ where $\psi_{\min}$ is the smallest nonzero singular value of the mean of an $N\times N$ matrix and $\sigma^2$ is the variance of the observed entries. Our development builds upon and complements a series of recent works that show that statistical inference for matrix completion is possible with a low signal-to-noise ratio when the data are missing uniformly at random. See, e.g., chen2019inference, chen2020noisy, xia2021statistical. Our results also draw an immediate comparison with the recent works by bai2021matrix, cahan2023factor who developed an inferential theory for the asymptotic principle component (APC) based approaches when the signal is much stronger, e.g., $\psi_{\min}^2\gtrsim \sigma^2N^2$. It is worth pointing out that the nuclear norm regularization and APC-based approach each has its own merits and requires different treatment. For example, APC-based methods usually assume that the factors are random and impose moment conditions to ensure that the factor structure is strong and identifiable, whereas our development assumes that the factors are deterministic but incoherent and allows for weaker signals.
Our work is motivated by a number of recent studies on the Tick Size Pilot Program, an experiment conducted by the Security and Exchange Commission (SEC) to evaluate the impact of widening the tick size on the market quality of small and illiquid stocks from 2016 to 2018. See, e.g., albuquerque2020price, chung2020tick, werner2022tick. The pilot consisted of three treatment groups with a control group: 1) The first treatment group was quoted in \$0.05 increments but still traded in \$0.01 increments (only Q rule), 2) The second treatment group was quoted and traded in \$0.05 increments (Q+T rule), 3) The third treatment group was quoted and traded in \$0.05 increments, and also subject to the trade-at rule (Q+T+TA rule). The trade-at rule, in general, prevents price matching by exchanges that are not displaying the best price. The control group was quoted and traded in \$0.01 increments. Previous studies chung2020tick on the effects of the quote rule (Q), trade rule (T), and trade rule (TA) on the liquidity measure are based on traditional regression or difference-in-difference methods and assume that the treatment effect is invariant with respect to time and unit. As we shall demonstrate, this assumption is problematic for the Tick Size Pilot Program data and there is significant heterogeneity in the treatment effect across both time and units. Indeed, more insights can be obtained using a potential outcome model with interactive fixed effects to capture such heterogeneity. To do so, we extend our methodology from estimating a single matrix to the simultaneous completion of multiple matrices, accounting for the multiple potential situations.
The remainder of this paper is organized as follows. Section (ref) introduces the method of using the nuclear norm penalized estimation when missing is not at random and provides the convergence rates of the estimator. Section (ref) discusses how to reduce bias and provides inferential theory using the debiased estimator. Section (ref) shows how our proposed methodology can be applied to infer the treatment effect in the Tick Size Pilot Program and presents the empirical findings of our analysis. Section (ref) examines the finite sample performance of our estimators using simulation studies. Finally, we conclude with a few remarks in Section (ref). All proofs are relegated to the Appendix due to the space limit.
In what follows, we use $\|\cdot\|_{\rm F}$, $\|\cdot\|$, and $\|\cdot\|_*$ to denote the matrix Frobenius norm, spectral norm, and nuclear norm, respectively. In addition, $\|\cdot\|_\infty$ denotes the entrywise $\ell_\infty$ norm, and $\|\cdot\|_{2,\infty}$ the largest $\ell_2$ norm of all rows of the matrix, i.e., $\|A\|_{2,\infty}=\max_i (\sum_j a_{ij}^2)^{1/2}$. For any vector $a$, $\|a\|$ denotes its $\ell_2$ norm. For any set $\mathcal{A}$, $|\mathcal{A}|$ is the number of elements in $\mathcal{A}$. We use $\circ$ to denote the Hadamard product or the entry-by-entry product between matrices of conformable dimensions. $a \lesssim b$ means $|a|/|b| \leq C_1$ for some constant $C_1>0$ and $a \gtrsim b$ means $|a|/|b| \geq C_2$ for some constant $C_2>0$. $c \asymp d$ means that both $c/d $ and $d/c$ are bounded. $a \ll b$ indicates $|a| \leq c_1|b|$ for some sufficiently small constant $c_1>0$ and $a \gg b$ indicates $c_2|a| \geq |b|$ for some sufficiently small constant $c_2>0$. In addition, $[K] = \{1, \dots , K\}$.
Consider a panel data setting where $M=(m_{it})_{1\le i\le N, 1\le t\le T}$ is a $N\times T$ matrix of rank $r$ ($\ll \min\{N, T\}$). We use $i$ as the cross-section index and $t$ as the time index. Following the convention of the matrix completion literature, we shall assume that the singular vectors of $M$ are incoherent in that there is a $\mu\geq 1$ such that $\|U_{M}\|_{2,\infty} \leq \sqrt{{\mu r}/{N}}$, $\|V_{M}\|_{2,\infty} \leq \sqrt{{\mu r}/{T}}$ where $U_M$ and $V_M$ denote the left and right singular vectors of $M$, respectively. The incoherence condition requires the singular vectors to be de-localized, in the sense that entries are not dominated by a small number of rows or columns.
Instead of $M$, we observe a subset of the entries of $Y=M+E$ where $E$ is a noise matrix whose entries are independent and identically distributed zero-mean, sub-Gaussian random variable, i.e., $\mathbb{E}[\epsilon_{it}^2] = \sigma^2$, $ \mathbb{E} [\exp(s \epsilon_{it})] \leq \exp(C s^2 \sigma^2)$, $\forall s \in \mathbb{R}$ and some constant $C>0$. Let $\Omega=(\omega_{it})_{1\le i\le N, 1\le t \le T}\in \{0,1\}^{N\times T}$ indicate the observed entries: $\omega_{it}=1$ if and only if $y_{it}$ is observed. The goal of noisy matrix completion is to estimate $M$ from $Y_\Omega:=\{y_{it}: \omega_{it}=1\}$. A popular approach to do so is the nuclear norm penalization: $$ \widetilde{M}=\operatorname*{arg\,min}_{A\in {\mathbb R}^{N\times T}}\left\{\|\Omega\circ (Y-A)\|_{\rm F}^2+\lambda\|A\|_\ast\right\}, $$ where $\lambda\ge 0$ is a tuning parameter. The properties of $\widetilde{M}$ are by now well understood in the case of missing completely at random, especially when the entries of $\Omega$ are independently sampled from a Bernoulli distribution. See, e.g., koltchinskii:2011, chen2020noisy. Instead, we are interested here in the situation where $\Omega$ is not random.
Situations when missing is not at random arise naturally in many causal panel models. Consider, for example, the evaluation of a program that takes effect after time $T_0$ for the last $N-N_0$ units. If $M$ is the potential outcome under the control, then we do not have observations of its entries for $i > N_0$ and $t> T_0$, e.g., $\Omega=1\{t\le T_0 {\rm \ or\ } i\le N_0\}$, yielding a block missing pattern as shown in the left panel of Figure (ref). A more general setting that often arises in causal panel data is the staggered adoption where units may differ in the time they are first exposed to the treatment, yielding a missing pattern as shown in the right panel of Figure (ref). See athey2021matrix,agarwal2021causal for other similar missing patterns that are common in the context of recommendation systems and A / B testing.
Note that if the entries are observed uniformly at random, then $$\|\Omega\circ (Y-A)\|_{\rm F}^2\approx {|\Omega|\over NT}{\mathbb E}\|Y-A\|_{\rm F}^2$$ for sufficiently large $N$ and $T$. The right-hand side is minimized by $M$, which justifies $\widetilde{M}$ as a plausible estimate of $M$. This intuition, however, no longer applies when $\Omega$ is not random and has more structured patterns. Our proposal to overcome this problem is dividing the missing entries into smaller groups and estimating each group via nuclear norm regularization. The main inspiration behind our method is the observation that $\widetilde{M}$ is a good estimate of $M$ when there are only a few missing entries, even if they are missing not at random.
It is instructive to start with a single treated period, e.g., $\Omega=1\{t \leq T-1 {\rm \ or\ } i\le N_0\}$. In this case, the number of missing entries is $|\Omega^c|=N-N_0$. Denote by $\psi_{\max}$ and $\psi_{\min}$ the largest and smallest nonzero singular value of $M$, respectively, and $\kappa = {\psi_{\max}}/{\psi_{\min}}$ its condition number. The following theorem provides bounds for the estimation error of $\widetilde{M}$.
Some immediate remarks are in order. Consider the situation where $\kappa,\mu,r = O(1)$, and $N \asymp T$. Ignoring the logarithmic term, the signal-to-noise ratio requirement given by Assumption (i) reduces to ${\psi_{\min}}\gg \sigma N^{1/2}$ which is significantly weaker than those in the existing literature. More specifically, if there is a single missing entry, e.g., $N_0=N-1$, agarwal2021causal suggest to partition the submatrix $(m_{it})_{1\le i<N, 1\le t<T}$ into $K$ smaller matrices. In particular, their Theorem 2 states that the best estimation error for their estimate is given by $$ |\widehat{m}^{\rm ADSS}_{NT} - m_{NT}| = O_p \left( \frac{1}{N^{1/4}} + \frac{1}{T^{1/4}} \right) $$ by setting $K \asymp N^{1/2}$. In contrast, under the assumptions of agarwal2021causal, $\sigma,\kappa,\mu,r$ are bounded and hence the convergence rate of our estimator is $$ |\widetilde{m}_{NT} - m_{NT}|=O_p \left( \left(\frac{1}{N^{1/2}} + \frac{1}{T^{1/2}}\right)\sqrt{\log(NT)} \right). $$
Theorem (ref) serves as our building block for dealing with more general and common missing patterns, which we shall now discuss in detail.
\paragraph{Single Treated Period.}
Note that Assumption (iii) of Theorem (ref) restricts the number of missing entries not to be large compared to $N$ and $T$. In particular, if $\kappa, \mu, r = O(1)$ and $N\asymp T$, then it requires that $|\Omega^c| = o\left(N\right)$. To deal with a larger number of missing entries, we shall leverage this result by splitting the missing entries into small groups and estimating them separately, as illustrated in Figure (ref).
Specifically, we split the missing entries into small groups, denoted by $\{\mathcal{G}_l \}_{1\leq l \leq L}$, and construct the submatrices $\{Y_l\}_{1\leq l \leq L}$ as illustrated in Figure (ref). For each $1 \leq l \leq L$, we estimate $M_l$, the corresponding submatrix of $M$, using the nuclear norm penalization:
where $N_l = N_0 + |\mathcal{G}_l|$ and $\Omega_l$ is the corresponding submatrix of $\Omega$. We shall then assemble these estimated submatrices into an estimate $\widetilde{M}$ of $M$. Note that each missing entry appears in one and only one of the submatrices and can therefore be estimated accordingly. The entries from $O$ in Figure (ref), e.g., the $N_0 \times (T-1)$ principle submatrix of $M$, on the other hand, are estimated for all groups. We can estimate these entries by averaging all of these estimates. Let the smallest nonzero singular value of $M_O$ be $\psi_{\min,O}$, where $M_O$ is the submatrix of $M$ corresponding to $O$. Denote by $u_i^\top$ and $v_t^\top$ the $i$-th row of $U_M$ and $t$-th row of $V_M$, respectively. We can then derive the following bounds from Theorem (ref).
The main difference from Theorem (ref) lies in Assumptions (iii) and (iv) of Corollary (ref). Assumption (iii) specifies how large a block can be. In principle, we can always take $|\mathcal{G}_l|=1$, that is, recovering one entry at a time so that this condition is trivially satisfied with sufficiently large $N_0$ and $T$. However, there could be enormous computational advantages in creating groups as large as possible because the number of $\widetilde{M}_l$s that need to be computed decreases with increasing group size.
Assumption (iv) can be viewed as an incoherence condition to ensure that the singular vectors of $M$ are not dominated by either the treated or untreated units. It is easy to see that when there are few missing entries, e.g., $N_0\approx N$, the condition is satisfied by virtue of the incoherence of $u_i$s. In general, if $\{u_i\}_{i \in [N]}$ is exchangeable or if the treated units are uniformly selected, then this condition is satisfied with high probability, at least for sufficiently large $N_0$, since $\frac{N}{N_0} \sum_{i \leq N_0} u_{i}u_{i}^\top\approx \sum_{i \leq N} u_{i}u_{i}^\top=I_r$ by means of matrix concentration inequalities tropp2015introduction.
\paragraph{Single Treated Unit.} A similar estimating strategy can also be used to deal with a single treated unit. Without loss of generality, let $\Omega=1\{t \leq T_0 {\rm \ or\ } i\leq N - 1\}$. Then the fully observed submatrix is $O = (y_{it})_{1\leq i \leq N-1, 1\leq t \leq T_0}$. As in the case of a single treated period, we split the missing entries into smaller groups, denoted by $\mathcal{G}_1,\ldots, \mathcal{G}_L$, by periods, and estimate them separately as before. Similar to Theorem (ref), we have the following bounds for the resulting estimate.
\paragraph{General Block Missing Pattern.} We can also apply the grouping and estimating procedure to general block missing structures such as that depicted in the left panel of Figure (ref), e.g., $\Omega=1\{t\le T_0 {\rm \ or\ } i\le N_0\}$, by estimating missing entries one period at a time (or one unit at a time). Denote by $\mathcal{G}_1, \mathcal{G}_2,\ldots, \mathcal{G}_L$ the groups of missing units (or periods). The following result again follows from Theorem (ref):
It is worth noting that both Corollary (ref) and Corollary (ref) can be viewed as special cases of Corollary (ref). It is also of interest to compare the rates of convergence with those of athey2021matrix. athey2021matrix considered a direct application of the nuclear norm penalized estimation to the full matrix. Their Theorem 2 states that $$ \frac{1}{\sqrt{NT}} \left\Vert\widetilde{M} - M\right\Vert_{\rm F} = O_p \left( \sqrt{\frac{T}{N}} + \sqrt{\frac{1}{T}}\right), $$ ignoring the logarithmic factors and $\sigma$, $r$, and $\left\VertM\right\Vert_\infty$. In other words, the estimate could be inconsistent when $N = O(T)$. On the other hand, the convergence rate of our estimator is given by $$\left\Vert\widetilde{M} - M\right\Vert_\infty = O_p \left( \sqrt{\frac{1}{N_0}} + \sqrt{\frac{1}{T_0}} \right), $$ up to a logarithmic factor when we assume $\kappa,\mu = O_p(1)$. Hence, our estimator is consistent as long as $\min\{N_0,T_0\}$ diverges. Furthermore, the simulation results in Section (ref) also show that applying the nuclear norm penalized estimation to the submatrix indeed performs much better than applying it to the full matrix as long as $N_0$ and $T_0$ are not too small.
\paragraph{Staggered Adoption.} More generally, we can take advantage of our estimation strategy for staggered adoption where there are $D$ number of adoption time points, says $T_1< \cdots< T_D$, and $D$ number of corresponding groups of treated units, says $G_1, \dots , G_D$. That is, for each $d \in [D]$, the units in $G_d$ adopt the treatment in the time period $T_d$. We can utilize the strategy for block missing patterns to estimate the missing entries. More specifically, denote by $M_{d,d'}$ the submatrix with missing entries corresponding to units in $G_{d}$ and time periods in $[T_{d'},T_{d'+1})$, with the convention that $T_{D+1}=T+1$, where $d \leq d' \leq D$. To estimate these missing entries, we can assemble a submatrix, denoted by $Y_{d,d'}$, with units untreated prior to $T_{d'+1}$ and time periods in $[1,T_d)\cup [T_{d'},T_{d'+1})$, as well as units in $G_d$ and time periods in $[1,T_d)$. As shown in Figure (ref), $M_{d,d'}$ is now the missing block of $Y_{d,d'}$, and can be estimated as described in the previous case.
Denote by $\mathcal{G}_1, \mathcal{G}_2,\ldots, \mathcal{G}_L$ the groups for missing units in $M_{d,d'}$ such as $\cup_{l\in [L]} \mathcal{G}_l = G_d$, $N_{d'}$ the number of units that are untreated prior to $T_{d'+1}$, and $\psi_{\min,O_{d,d'}}$ the smallest singular value of the submatrix $M_{O_{d,d'}}=(m_{it})_{1\leq i \leq N_{d'}, 1 \leq t \leq T_d}$. The performance of the resulting estimate is given by Corollary (ref).
It is worth comparing the rates of convergence with those of bai2021matrix which apply their TW algorithm to the full matrix. For all missing entries, the convergence rates of the estimators in bai2021matrix are $O_p\left(\frac{1}{\sqrt{N_D}} + \frac{1}{\sqrt{T_1}}\right)$. On the other hand, if we assume $\kappa,\mu = O_p(1)$, the convergence rate of our estimator is $O_p\left(\frac{1}{\sqrt{N_{d'}}} + \frac{1}{\sqrt{T_d}}\right)$ up to a logarithmic factor. Since $N_{d'} > N_D$ and $T_{d} > T_1$ for all $d' < D$ and $d > 1$, our convergence rate is faster than that of bai2021matrix except for the estimation of missing entries in part $M_{1,D}$ for which both estimates have similar rates of convergence. This shows the advantage of exploiting submatrices for the imputation of missing entries.
We now turn our attention to inferences. While the nuclear norm regularized estimator $\widetilde{M}$ enjoys good rates of convergence, it is not directly suitable for statistical inferences due to the bias induced by the penalty. To overcome this challenge, we propose an additional projection step after applying the nuclear norm penalization in recovering missing entries from group $\mathcal{G}_l$:
where $\mathcal{P}_{r}(B) = \operatorname*{arg\,min}_{A:\textrm{rank}(A)\leq r} \left\VertA-B\right\Vert_F$ is the best rank-$r$ approximation of $B$. We now discuss how this enables us to develop an inferential theory for estimating the missing entries. To fix ideas, we shall focus on inferences about the average of a group of entries at a given time period, e.g., $\sum_{i \in \mathcal{G}} m_{it_0}/|\mathcal{G}|$, where $\mathcal{G} \subseteq [N]$.
\paragraph{Block Missing Patterns.} We shall begin with general block missing patterns, e.g., $\omega_{it}=1$ if $t\le T_0$ or $i\le N_0$. Note that both the single treated period and single treated unit examples from the previous section can be viewed as special cases with $T_0=T-1$ and $N_0=N-1$, respectively.
Suppose that we are interested in the inference of the average of a group of entries at the time $t_0$, $\sum_{i\in \mathcal{G}} m_{it_0}/|\mathcal{G}|$, where $\mathcal{G} \subseteq \{ 1 , \cdots, N\}$ and $t_0>T_0$. Similar to before, we split the interesting group, $\mathcal{G}$, into smaller subgroups, denoted by $\{\mathcal{G}_l \}_{0 \leq l \leq L}$ with the convention that $\mathcal{G}_0 = \mathcal{G} \cap \{1, \cdots , N_0\}$, and construct the corresponding submatrices $\{Y_l\}_{1 \leq l \leq L}$ as illustrated in Figure (ref), and construct $Y_0 = [ (y_{it})_{i\leq N_0,t \leq T_0} \ \ (y_{it})_{i\leq N_0,t = t_0}]$ if $\mathcal{G}_0 \neq \emptyset$.
Recall that $\psi_{\min,O}$ is the smallest nonzero singular value of the $N_0 \times T_0$ matrix $M_O = (m_{it})_{1\leq i \leq N_0, 1 \leq t \leq T_0}$. The following theorem establishes the asymptotic normality of the group average estimator, $\sum_{ i \in \mathcal{G}}\widehat{m}_{it_0}/|\mathcal{G}|$.
\paragraph{Staggered Adaption.} More generally, consider the case of staggered adoption when there are $D$ number of adoption time points, $T_1<T_2< \cdots<T_D$, and $D$ number of corresponding groups of treated units, $G_1, \ldots , G_D$. As in the previous situation, suppose that we are interested in inference for the group average at time $t_0$. Denote by $N_{0}$ the number of units that are untreated until $t_0$, and by $T_0$ the number of time periods where $\{1,\dots,N_0\}$ is untreated, respectively.
We proceed by first splitting $\mathcal{G}$ into smaller groups, denoted by $\{\mathcal{G}_l\}_{0\leq l \leq L}$ with the convention that $\mathcal{G}_0 = \mathcal{G} \cap \{1, \cdots , N_0\}$. In doing so, we want to make sure that all units in each subgroup $\{\mathcal{G}_l\}_{1\leq l \leq L}$ have the same adoption time point, e.g., $\mathcal{G}_l \subseteq G_{d_l}$, as illustrated in Figure (ref). Denote by $D_\mathcal{G}=\{d_l: 1\le l\le L\}$ and by $\psi_{\min,O_{d}}$ the smallest singular value of the submatrix $M_{O_d}=(m_{it})_{1\leq i \leq N_0, 1\leq t \leq T_d}$.
\paragraph{Variance Estimation.} In practice, to use the results above for inferences, we also need to estimate the variance. To this end, let $\widetilde{U}_l\widetilde{D}_l\widetilde{V}_l^\top$ be the SVD of $\mathcal{P}_r(\widetilde{M}_l)$. Denote by $\widetilde{X}_l = \widetilde{U}_l \widetilde{D}_l^{1/2}$ and $\widetilde{Z}_l = \widetilde{V}_l \widetilde{D}_l^{1/2}$. They can be viewed as estimates of rescaled left and right singular vectors. However, as such, they are significantly biased and the bias can be reduced by considering instead $$ \widehat{X}_l = \widetilde{X}_l \left(I_r + \lambda_l (\widetilde{X}_l^\top \widetilde{X}_l) \right)^{1/2}, \ \ \widehat{Z}_l = \widetilde{Z}_l \left(I_r + \lambda_l (\widetilde{Z}_l^\top \widetilde{Z}_l) \right)^{1/2}. $$ We can then use $\widehat{X}_l$ and $\widehat{Z}_l$ in place of the left and right singular vector in defining $\mathcal{V}_{\mathcal{G}}$, leading to the following variance estimate
where $\widehat{\bar{X}}_{\mathcal{G}_l} = \frac{1}{|\mathcal{G}_l|} \sum_{j \in \mathcal{G}_l}\widehat{X}_{l,j}$, $\widehat{\sigma}^2 = \frac{1}{N_0T_0} \sum_{i \leq N_0, t \leq T_0} \widehat{\epsilon}_{it}^2$, and $\widehat{\epsilon}_{it} = y_{it} - \widehat{m}_{it}$. The following corollary shows that asymptotic normality established in Theorem (ref) continues to hold if we use this variance estimate.
Since Theorem (ref) is a special case of Theorem (ref), the variance estimator can also be used for Theorem (ref). Specifically, it is enough to change from $T_{d_l}$ in $\widehat{\mathcal{V}}_\mathcal{G}$ to $T_0$ for Theorem (ref).
Our work was motivated by the analysis of the Tick Size Pilot Program, which we shall now discuss in detail to demonstrate how the proposed methodology can be applied in causal panel data models.
\paragraph{Background.} In October 2016, the SEC launched the Tick Size Pilot Program to evaluate the impact of an increase in tick sizes on the market quality of stocks. As noted before, the pilot consisted of a control group and three treatment groups:
This pilot program has attracted considerable attention, and there are a growing number of studies on the impact of these changes on market quality, often represented by a liquidity measure such as the effective spread since its conclusion in 2018. See, e.g., albuquerque2020price, chung2020tick, griffith2019making, rindi2019us, werner2022tick.
\paragraph{Data.} Data for control variables were obtained from the Center for Research in Security Prices (CRSP) and the daily share-weighted dollar effective spread data from the Millisecond Intraday Indicators by Wharton Research Data Services (WRDS). A key control variable introduced by chung2020tick is TBC which measures the extent to which the new tick size (\$0.05) is a binding constraint on the quoted spreads in the pilot periods and is estimated by the percentage of quoted spreads during the day that are equal to or less than 5 cents, which is the new minimum quoted tick size under the Q rule. Specifically, we calculate the percentage of NBBO updates with quoted spread less than or equal to 5 cents for each day. Using the TBC variable, we can check the effect of an increase in the minimum quoted spread (from 1 cent to 5 cents) on the effective spread.
A data-cleaning process similar to chung2020tick yields a total of $N=1,461$ stocks with $N_0=735$ in the control group, $N_1=254$ in the Q group, $N_2=244$ in the Q+T group, and $N_3=228$ in the Q+T+TA group. Following chung2020tick, data from Oct 1, 2015 to Sep 30, 2016 were used as the pre-pilot periods and Nov 1, 2016 to Oct 31, 2017 as the pilot periods, i.e., $T_0 = 253$ and $T_1 = 252$ for daily data. See chung2020tick for further discussion of data collection. As is common in previous studies, we consider the daily effective spread in cents as a measure of liquidity. Denote by $y_{it}^{(d)}$ the potential outcome for stock $i$ at time $t$ under treatment $d$ with the convention that $d=0, 1,2,3$ corresponds to the control, the Q rule, the Q + T rule, and the Q + T + TA rule, respectively. The four matrices $Y^{(d)}=(y_{it}^{(d)})_{1\le i\le N, 1\le t\le T}$ have block missing patterns, as shown in Figure (ref).
\paragraph{Model.} Previous studies of the effects of the quote (Q) rule, the trade (T) rule, and the trade-at (TA) rule on the liquidity measure are usually based on traditional regression or difference-in-difference methods by assuming that the treatment effect is constant across all units and time periods. For instance, chung2020tick postulated $y_{it} = y_{it}^{(d)}$ if unit $i$ receives treatment $d$ at time $t$ where the potential outcomes $$ y_{it}^{(d)}=m_{it}^{(d)}+x_{it}^\top \beta+\epsilon_{it} $$ and
Here, $\mu^{(0)}=0$, $\mu^{(1)}, \mu^{(2)},\mu^{(3)}$, $\alpha_i$s and $\delta_t$s are unknown parameters, and $x_{it}$ is a set of control variables that includes typical stock characteristics like stock prices and trading volumes, and TBC, a variable measuring the extent to which the new tick size (\$0.05) is a binding constraint on the quoted spreads in the pilot period. See Section (ref) in the Appendix for further details. It is worth noting that, in addition to the treatment effects ($\mu^{(1)}$, $\mu^{(2)}$ and $\mu^{(3)}$), their differences $\theta^{(d)}:=\mu^{(d)}-\mu^{(d-1)}$ are also of interest, as they represent the treatment effects of quote rule, trade rule, and trade-at rule, respectively.
However, (ref) fails to account for the significant heterogeneity in the treatment effects across units and time periods. To this end, we shall consider a more flexible model:
where $\zeta_i$ is a $r$-dimensional vector of (latent) unit specific characteristics and $\eta_{t}^{(d)}$ is the corresponding coefficients of $\zeta_i$ at time $t$ in the potential situation $d$. As we shall see later in this section, (ref) allows us to get more insights into the treatment effects of the pilot program.
One of the key assumptions of Model (ref) is that the subspace spanned by the left singular vector of $M^{(d)}=(m_{it}^{(d)})_{1\le i\le N, 1\le t\le T}$ for all $d=1,2,3$ is included in the subspace spanned by the left singular vector of $M^{(0)}$. agarwal2020synthetic propose a subspace inclusion test to check the validity of this assumption. We carried out this test on the pilot data, which confirms this is a reasonable assumption.
We note that similar low-rank models have also been considered by agarwal2020synthetic and chernozhukov2021inference earlier. However, it is unclear how their methodology can be adapted for the analysis of the Tick Size Program. For example, chernozhukov2021inference impose conditions on the missing pattern that are clearly violated by the pilot data; agarwal2020synthetic only study the average treatment effect and so cannot be used to assess the heterogeneity or dynamics of the treatment effects across units and time periods, respectively.
\paragraph{Estimation.} We now discuss how we can apply the methodology in the previous sections to analyze the tick size program, and in particular to estimate and make inferences about (ref). More specifically, we are interested in estimating the group-averaged treatment effects: for an interesting group of treated units $\mathcal{G}$, $$ \mu^{(d)}_t:={1\over |\mathcal{G}|}\sum_{i\in\mathcal{G}} [m^{(d)}_{it}-m^{(0)}_{it}], $$ and their differences: $$ \theta^{(d)}_t:=\mu^{(d)}_t-\mu^{(d-1)}_t, $$ for $t>T_0$. Especially, when $\mathcal{G}$ is a certain unit, it reduces to the individual treatment effect and if $\mathcal{G}$ is the group of all treated units, it becomes the cross-sectional averaged treatment effect. To this end, we shall derive estimates for $m^{(d)}_{it}$ under Model (ref).
First, note that, for this particular application, one of the covariates (TBC) is only present for the pilot periods. Therefore, we cannot hope to estimate the regression coefficient $\beta$ using the pre-pilot data alone, as suggested by bai2021matrix. Nonetheless, under (ref), $y_{it}$s follow an interactive fixed effect model: $$y_{it} = x_{it}^\top \beta + L_{it} + \epsilon_{it}$$ for some low rank components $L_{it}$ and therefore the regression coefficient $\beta$ can be estimated at the rate of $O_p(1/\sqrt{NT})$. See bai2009panel for details. This is much faster than that of the estimates of $m^{(d)}_{it}$. For brevity, we shall, therefore, treat the regression coefficient $\beta$ as known in what follows, without loss of generality.
For $d=0$, we can apply the method proposed in the previous sections to the potential outcome panel $\tilde{Y}^{(0)}_{it}=(y_{it}^{(0)}-x_{it}^\top\beta)_{1\le i\le N, 1\le t\le T}$. As illustrated in Figure (ref), it has a block missing pattern with $\omega_{it}^{(0)}=1$ if and only if $t\le T_0$ or $i\le N_0$. As such, we can derive estimates $\widehat{m}_{it}^{(0)}$ for $t > T_0$.
When $d>0$, we can only observe $y_{it}^{(d)}$ if unit $i$ receives treatment $d$ and $t>T_0$, so our method cannot be applied directly. Instead, we shall combine all observations from prepilot periods and these observations to form a panel $\tilde{Y}^{(d)}$ whose $(i,t)$ entry is $y_{it}^{(d)}-x_{it}^\top\beta$ if $i$ receives treatment $d$ and $t>T_0$, is $y_{it}^{(0)}-x_{it}^\top\beta$ if $t\le T_0$, and is missing otherwise. Let $\tilde{M}^{(d)}$ be a $N\times T$ matrix whose $(i,t)$ entry is $m^{(0)}_{it}$ if $t\le T_0$, and $m^{(d)}_{it}$ otherwise. $\tilde{Y}^{(d)}$ can be viewed as the noisy observation of $\tilde{M}^{(d)}$ with a block missing pattern: $\omega_{it}^{(d)}=1$ if and only if unit $i$ receives treatment $d$ or $t\le T_0$. Under (ref), $\tilde{m}^{(d)}_{it}=\zeta_i^\top\tilde{\eta}^{(d)}_t$ where $\tilde{\eta}^{(d)}_t=\eta_t^{(0)}$ if $t\le T_0$ and $\eta^{(d)}_t$ otherwise. Therefore, we can again apply our method to $\tilde{Y}^{(d)}$ to obtain estimates $\hat{m}_{it}^{(d)}$ for $t>T_0$.
We shall then proceed to estimate the treatment effects by $$ \widehat{\mu}^{(d)}_t:={1\over |\mathcal{G}|}\sum_{i\in \mathcal{G}} [\widehat{m}^{(d)}_{it}-\widehat{m}^{(0)}_{it}] \qquad {\rm and}\qquad \widehat{\theta}^{(d)}_t:={1\over |\mathcal{G}|}\sum_{i\in\mathcal{G}} [\widehat{m}^{(d)}_{it}-\widehat{m}^{(d-1)}_{it}]. $$
\paragraph{Inferences.} We can also use the results from the last section to derive the asymptotic distribution for $\widehat{\mu}^{(d)}_t$ and $\widehat{\theta}^{(d)}_t$. More specifically, let $M$ be a $N\times (T+3T_1)$ matrix that combines all observed outcomes: the first $T$ columns of $M$ consist of the potential outcomes under the control for the whole periods $(m_{it}^{(0)})_{i\leq N, t\leq T}$, the next $T_1$ columns the potential outcomes under the Q rule for the pilot periods $(m_{it}^{(1)})_{i\leq N, t > T_0}$, followed by those under the Q+T rule again for the pilot periods $(m_{it}^{(2)})_{i\leq N, t > T_0}$, and finally those under the Q+T+TA rule $(m_{it}^{(3)})_{i\leq N, t > T_0}$. Note that $M$ is also a rank-$r$ matrix. Let $M=UDV^\top$ be its singular value decomposition. Denote by $u_i^\top$ and $v_t^\top$ the $i$-th row vector of $U$ and $t$-th row vector of $V$, respectively. In addition, denote by $\mathcal{I}_d$ the group of units treated by treatment $d$ with the convention that $\mathcal{I}_0$ is the control group. Then, under suitable conditions, we have
$\mathcal{V}_{\mu} = \mathcal{V}_{\mathcal{G}}(d,0)$ and $\mathcal{V}_{\theta} = \mathcal{V}_{\mathcal{G}}(d,d-1)$ where
Similar to before, the variance can be replaced by its estimate. Due to the space limit, we shall defer the formal statements and proofs, as well as derivations of the variance estimator to the Appendix.
\paragraph{Fixed Effects vs Interactive Effects.} We begin with some exploratory analyses to illustrate the impact of the pilot program. The top left panel of Figure (ref) gives the boxplots of difference in the effective spread, averaged over time, after and before the pilot. There are a few units with differences that are much larger in magnitude than usual. For better visualization, the top right panel zooms in with a difference between -10 cents and 10 cents. Taken together, it is clear that the three treatment groups have a significant impact on the effective spread.
The treatment effect of the pilot, however, differs between units. The bottom panels of Figure (ref) show barplots of the time series of the effective spread of two typical stocks. The impact of the treatment is much clearer for the stock depicted in the bottom right panel.
The difference in treatment effect among the units suggests that the interactive effect model is more suitable than the fixed effect model used in the previous studies. Note that the fixed effect model (ref) can be viewed as a special case of the interactive effect model (ref) with $\zeta_i = [1 \ \ \alpha_i]^\top$, $\eta_t^{(d)} = [\delta_t+ \mu^{(d)} \ \ 1]^\top$. We conducted a Hausman-type model specification test to further show that the fixed effect model is inadequate in capturing the heterogeneity of the treatment effect. More specifically, denote our estimator of $\theta_{it}^{(d)} \coloneqq m_{it}^{(d)} - m_{it}^{(d-1)}$ by $\hat{\theta}^{(d)}_{it}$ and the two-way fixed effect estimator of $\theta^{(d)} \coloneqq \mu^{(d)} - \mu^{(d-1)} (= m_{it}^{(d)} - m_{it}^{(d-1)})$ in Model (ref) by $\tilde{\theta}^{(d)}$. We considered the following test statistic for model specification: $$T-stat_{\rm ms} = \max_{i \in \mathcal{N}_{tr}, T_0 < t \leq T} \max_{1 \leq d \leq 3} |\hat{\tau}^{(d)}_{it}|$$ where $\mathcal{N}_{tr}$ is the group of all treated stocks, $\hat{\tau}^{(d)}_{it} = \hat{\mathcal{V}}_{d,it}^{-1/2} (\hat{\theta}^{(d)}_{it} - \tilde{\theta}^{(d)})$, and $\hat{\mathcal{V}}_{d,it}$ is the estimator of the asymptotic variance of $\hat{\theta}^{(d)}_{it} - \tilde{\theta}^{(d)}$. Moreover, to test whether $\theta_{it}^{(d)}$ is time and unit invariant or not, we also considered the test statistic such that $$ T-stat_{(d)} = \max_{i \in \mathcal{N}_{tr}, T_0 < t \leq T} \left\vert\hat{\mathcal{V}}_{d,it}^{-1/2} (\hat{\theta}_{it}^{(d)} - \bar{\hat{\theta}}^{(d)})\right\vert $$ where $\bar{\hat{\theta}}^{(d)} =\frac{1}{|\mathcal{N}_{tr}|T_{1}}\sum_{i\in \mathcal{N}_{tr}, T_0 < t \leq T}\hat{\theta}_{it}^{(d)}$.
We derived the large sample distributions of the test statistics under the null and corresponding critical values using the Gaussian bootstrap method belloni2018high. And the null hypothesis that Model (ref) is well specified and the null hypotheses that $\{\theta_{it}^{(d)}\}_{1\leq d \leq 3}$ are time and unit invariant are all rejected at 1% significance level, again indicating that Model (ref) is misspecified and $\{\theta_{it}^{(d)}\}_{1\leq d \leq 3}$ are time and unit variant.
To further illustrate the heterogeneity of the treatment effect, we compute the estimated unit-specific treatment effect averaged over time: $\bar{\hat{\theta}}^{(d)}_i:=T_1^{-1}\sum_{t>T_0} \hat{\theta}^{(d)}_{it}$ and Figure (ref) gives the kernel density estimates of these unit-specific treatment effects for the Q rule, T rule and TA rule respectively. It is evident from these density plots that there is considerable amount of variation and skewness among the estimated treatment effects across units.
Note that a key assumption behind the interactive effect model is that the unit specific characteristic $\zeta_i$ remains the same across all treatment groups as well as the control group so that they can be learned from the pre-pilot periods and utilized for the estimation of $m^{(d)}_{it}$ during the pilot period. This amounts to the assumption that the left singular space of $M^{(d)}$ is included in that of $M^{(0)}$. To check the validity of the assumption, we carry out the subspace inclusion test for $d=1,2,3$ introduced in agarwal2020synthetic, and the test statistics are $0.15$, $0.19$ and $0.11$ with corresponding critical values at 95% level $0.43$, $0.48$ and $0.28$. Additionally, we also confirm that the ranks of $(m^{(0)}_{it})_{i \in \mathcal{I}_d , t \leq T_0}$ and $[(m^{(0)}_{it})_{i \in \mathcal{I}_d , t \leq T_0} \ \ (m^{(d)}_{it})_{i \in \mathcal{I}_d , t > T_0}]$ are the same for all $1 \leq d \leq 3$ using the typical rank estimation method ahn2013eigenvalue, which implies the validity of this assumption.
The rank test also indicates that $r=1$ is an appropriate choice for the pilot data. The associated $R^2$ is 0.79. This is to compared with the fixed effect model (ref) whose $R^2$ is 0.67 with the same degrees of freedom. This again suggests that the interactive effect model (ref) is preferable.
\paragraph{Dynamics of Treatment Effects.} Next, we examine the dynamics of the treatment effects of the Q rule, the T rule, and the TA rule.
To better visualize the dynamics, we plot in Figure (ref) the estimated daily treatment effects along with their 95% confidence interval, adjusted with Bonferoni correction. To gain further insights, we also plot in Figure (ref) the weekly average of the estimated daily treatment effects, again with their 95% confidence interval adjusted with Bonferoni correction. Note that to do so, we need to consider the estimator of the form $$\frac{1}{|\mathcal{S}|}\frac{1}{|\mathcal{N}_{tr}|} \sum_{t \in \mathcal{S}} \sum_{i \in \mathcal{N}_{tr}} \hat{\theta}^{(d)}_{it}$$ where $\mathcal{S}$ is a week of interest. We can generalize the inferential theory from the previous section straightforwardly with the new variance: $$ \sum_{\rho \in \{d, d-1 \}} \left[ \frac{\sigma^2}{|\mathcal{S}|} \bar{u}_{\mathcal{N}_{tr}}^\top \left(\sum_{j \in \mathcal{I}_\rho} u_j u_j^\top\right)^{-1} \bar{u}_{\mathcal{N}_{tr}} \right] + \frac{\sigma^2}{|\mathcal{N}_{tr}|} \bar{v}_{\textrm{diff}}^\top \left(\sum_{s \leq T_0} v_{s} v_{s}^\top \right)^{-1} \bar{v}_{\textrm{diff}} , $$ where $$\bar{v}_{\textrm{diff}} = \frac{1}{|\mathcal{S}|} \sum_{t \in \mathcal{S}} v_{(d \cdot T_1+t)}-v_{((d-1) \cdot T_1+t)}.$$
$\theta_{it}^{(2)}$ and $\theta_{it}^{(3)}$ can be interpreted as the treatment effects of the T rule and the TA rule. As expected by theory in the literature, we have the positive treatment effects of T rule most of the time. The T rule has a negative effect on price improvements, as liquidity providers are less likely to offer them when the minimum possible price improvement is larger. For example, if the T rule makes the minimum possible price improvement to be 5 cents, liquidity providers who would have been willing to provide less than 5 cents of price improvements are unlikely to offer any price improvement at all. Since the effective spread is “quoted spread - price improvement”, we can expect that treatment effects of the T rule is positive. Here, we use the following definitions: $\texttt{Quoted Spread}_t = A_t - B_t$, $\texttt{Effective Spread}_{t} = 2(P_t -\frac{A_t + B_t}{2})$, and $\texttt{Price Improvement}_{t} = 2(A_t - P_t)$, where $A_t$ is the national best ask price at time $t$, $B_t$ is the national best bid price at time $t$, and $P_t$ is the transaction price.
Interestingly, one can observe that the periods associated with large effects of the T rule usually correspond to large trading volumes. In particular, there were large trading volumes in November, early and mid-December in 2016, March, mid and late June, early August, early September, and late October in 2017, and, by and large, these periods coincide with periods with larger impact of the T rule. In general, the correlation coefficients between the estimated effect of the T rule and the trading volume is $0.33$. This suggests that the effect of the T rule becomes stronger when transactions are more active. This agrees with the well-known fact that price improvement is more likely to occur when stocks are actively traded, and therefore the effect of the T rule through price improvement will become amplified and strong when trades are active.
Moreover, we find that the treatment effects of the TA rule are negative most of the time. The TA rule increases visible liquidity by exposing hidden liquidity because, under the TA rule, a venue should display the best bid or ask to execute incoming market orders at the NBBO. It implies a decrease in the quoted spread and a smaller room for price improvements. chung2020tick expect that the effect on the quoted spread is likely to be greater than the effect on price improvements, and so the TA rule decreases the effective spread. Our result corroborates with their conjecture. Further discussion about the empirical findings is given in Section (ref) in the Appendix.
To further demonstrate the practical merits and finite sample performance of our methodology, we conducted several sets of simulation experiments.
The first set of simulations was designed to compare the performance of the proposed estimator with that of other existing estimators in a staggered adoption setting. Here, the size of “no adoption” group (G0) was set to 200. There are three adoption groups (G1, G2, G3), and the size of each adoption group was set to 100. The number of time points was 500 with G1 adopting the intervention at the 201st time period, G2 at the 301st time period, and G3 at the 401st time period. The potential outcome under the control follows a low-rank model $y_{it}^{(0)}=\zeta_i^\top\eta_t^{(0)} + \varepsilon_{it}$ where the noise $\varepsilon_{it}$ was sampled independently from the standard normal distribution. The unit specific characteristics $\zeta_i$s were sampled independently from $\mathcal{N}((2.5/\sqrt{2},2.5/\sqrt{2})^\top, I_2)$ for G0, $\mathcal{N}((1/\sqrt{2},1/\sqrt{2})^\top, I_2)$ for G1, $\mathcal{N}((1.5/\sqrt{2},1.5/\sqrt{2})^\top, I_2)$ for G2, and $\mathcal{N}((\sqrt{2},\sqrt{2})^\top, I_2)$ for G3. In addition, the corresponding coefficient $\eta_t^{(0)}$s were sampled independently from $\mathcal{N}((1/\sqrt{2},1/\sqrt{2})^\top, I_2)$.
To fix ideas, we consider estimating the missing potential outcome $m^{(0)}_{it}$ of a randomly chosen unit in G2 during the last time period ($t=500$) using different estimators including ours (\verb+CY+) along with those from bai2021matrix (\verb+BN+), agarwal2021causal (\verb+ADSS+) and athey2021matrix (\verb+ABDIK+). For ADSS, following the recommendation in agarwal2021causal, we set the number of sub-subgroup $K$ to be $K \asymp |AR^{(k)}|_o^{1/3}$. Table (ref) reports the RMSE, summarized from 1,000 simulation runs. The performance of CY, BN, and ADSS are superior to that of ABDIK with CY slightly better than BN and ADSS.
In addition, we recorded the coverage probabilities of the (asymptotic) confidence intervals associated with each method, with the exception of ABDIK for which such inferential tools have not been developed in the literature. From Table (ref), we can see that the coverage probabilities of ADSS are not close to the nominal level, indicating that the asymptotic distributional properties may not provide good approximations in this setting. On the other hand, our method and BN are more accurate, with ours more closely following the target probabilities.
Our next set of simulations mimics the setting of the pilot program studied in the previous section. More specifically, we considered Model (ref) with two treatment groups, $\mathcal{I}_1$ and $\mathcal{I}_2$, and a control group, $\mathcal{I}_0$. Each treatment group receives a different treatment in the pilot periods.
We set $r=2$ and generated the unit specific characteristics from $\zeta_{i} \sim \mathcal{N}((1/\sqrt{2},1/\sqrt{2})^\top, I_2)$, $\varepsilon_{it} \sim \mathcal{N}(0, 1)$, $\eta_{it}^{(0)} \sim \mathcal{N}((1/\sqrt{2},1/\sqrt{2})^\top, I_2)$, $\eta_{it}^{(1)} \sim \mathcal{N}((1.5/\sqrt{2},1.5/\sqrt{2})^\top, I_2)$, and $\eta_{it}^{(2)} \sim \mathcal{N}((\sqrt{2},\sqrt{2})^\top, I_2)$. In addition, two control variables were included: $x_{1,it}$ is generated from $\mathcal{N}(0, 1)$ while $x_{2,it}$ is generated from $\mathcal{N}(0, 1)$ if $t \in \text{Pilot period}$ and 0 otherwise. We set the regression coefficient $\beta = (1,1)^\top$ and estimated it using the interactive fixed effect estimation with data of whole periods. The numbers of $\mathcal{I}_0$, $\mathcal{I}_1$, and $\mathcal{I}_2$ were set to 250 and the numbers of pre-pilot periods and pilot periods were both set to 250.
As before, we estimated $\mu_{it}^{(d)}$ and $\theta_{it}^{(d)}$ for $1 \leq d \leq 2$ of a randomly chosen unit in $\mathcal{I}_2$ at the last period ($t = 500$). Table (ref) reports the coverage probabilities of our methods for $\mu_{it}^{(1)}$, $\mu_{it}^{(2)}$, and $\theta_{it}^{(2)}$, summarized from 1,000 simulation runs. It is evident that our coverage probabilities are quite close to the corresponding target probabilities. This is complemented by Figure (ref) that shows the histograms of the standardized estimates (t-statistics) along with the standard normal distribution, which again confirms the asymptotic normality of our estimates.
Our final experiment is similar to that from agarwal2021causal and athey2021matrix and is based on the tobacco sales data of abadie2010synthetic. In 1988, California introduced the first anti-tobacco legislation in the United States (Proposition 99) and to study the effect of this legislation on tobacco sales, abadie2010synthetic used the per capita cigarette sales data which was collected across 39 U.S. states from 1970 to 2000. We considered the time horizon of $n = 31$ years and restricted our focus to the $m = 38$ untreated states (excluding California) in their dataset. This data was encoded into a $38 \times 31$ matrix, $Y$, where the entry $y_{it}$ represents the potential outcome of per capita cigarette sales (in packs) for state $i$ in year $t$ under control, i.e., without any intervention in place.
To generate MNAR data, we artificially introduced interventions to a subset of states where the probability that a state adopts an intervention (e.g., tobacco control program) depends on their change in cigarette sales pre-1986 and post-1986. More specifically, we considered the following adoption protocol: First, we clustered states into four categories — severe, moderate, mild, and good — based on their percentage change in average cigarette sales during 1986-2000 compared to that during 1970-1985. The severe states are the states where average cigarette sales are hardly reduced ($-0\% \sim -10\%$, MO,WV,SC,AL,AR,TN), and the moderate states are the states whose percentage change is between $-10\%$ and $-15\%$ (KY,DE,GA,IN,OH,MS). The mild states are the states where the percentage change is between $-15\%$ and $-20\%$ (NE,LA,IA,SD,WI,PA). The rest are good states ($-20\% \sim$).
We then designated the timing and probability of intervention for mild, moderate, severe, and good states differently. Half of the severe states adopt an intervention in 1986 and the other half in 1991. Half of the moderate states adopt the intervention at 1991 and the other half in 1996. Half of the mild states adopt the intervention in 1996, and the other half do not adopt the intervention. In addition, the good states do not adopt the intervention at all. This setup reflects the scenario in which a state whose average sales may not be reduced sufficiently without the intervention is more likely to adopt the intervention early.
Table (ref) shows the average RMSE of missing components caused by the intervention in 10 experiments. Here, the missing components mean the potential “control (no adoption)” outcomes in the intervention period. The only randomization lies in the resampling of the observation patterns. We can check that ABDIK performs relatively poorly. In addition, the performance of our estimator is slightly better than that of BN and ADSS.
This article develops an inference framework for the matrix completion when missing is not at random and without the need for strong signals. One of the key observations to our development is that if the number of missing entries is small enough compared to the size of the panel, they can be well estimated even if missing is not at random. We judicially divide the missing entries into smaller groups and use this observation to provide accurate estimates and efficient inferences. Moreover, we showed that our proposed estimate, even with fairly weak signals, is asymptotically normal with suitable debiasing. As an application, we studied the treatment effects in the tick size pilot program, an experiment conducted by the SEC to assess the impact of tick size extension on the market quality of small and illiquid stocks from 2016 to 2018. While previous studies on this program were based on traditional regression or difference-in-difference methods by assuming that the treatment effect is invariant with respect to time and unit, we observed significant heterogeneity in treatment effects and gained further insights about treatment effects in the pilot program using our estimation method. Lastly, we conducted simulation experiments to further demonstrate the practical merits of our methodology.