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.
69,652 characters · 16 sections · 11 citation commands
Algorithmic Subsampling under Multiway Clustering
In the era of big data, it is not uncommon that data sets are so large that researchers may not need to use the whole sample for statistical inference to draw informative conclusions. Furthermore, computational bottlenecks in time and/or memory may even prohibit econometric analyses with such large data sets. The recent econometrics literature \citep*[e.g.,][]{LeeNg2020ARE,LeeNg2020sketching} suggests methods to deal with these circumstances that started to arise in today's data rich environments. The algorithmic subsampling or data sketching explored by these authors paves the way for econometric and statistical analyses based on random subsampling of big data.
The existing study of the algorithmic subsampling focuses on i.i.d. cases, and it has been “silent about how to deal with data that are dependent over time or across space” LeeNg2020ARE. On the other hand, some of big data may exhibit cross-sectional statistical dependence, such as multiway clustering. For instance, common scanner data (leading examples of big data) are clustered in two ways by markets and products. Common demand shocks within a market may induce statistical dependence among different products within that market. Similarly, common supply shocks by a producer may induce statistical dependence among different markets within the product produced by that producer. “A natural stochastic framework for the regression model with multiway clustered data is that of separately exchangeable random variables.” \citep*{mackinnon2019wild}
In this paper, we propose a novel method of algorithmic subsampling for separately exchangeable random variables, which we will refer to as the multiway algorithmic subsampling, and develop asymptotic statistical properties of this method. We first establish basic theories for the multiway algorithmic subsample means, namely their uniform weak law of large numbers and central limit theorem, which differ from the standard ones in a meaningful way. In particular, the form of the central limit theorem that is unique to the multiway algorithmic subsampling entails a practically useful property of robustness against potential degeneracy. Researchers do not know ex ante how the data in use are affected by the cluster sampling. In case that the cluster-specific shocks have no mean effect on the data, the standard multiway-cluster-robust asymptotic distribution without the algorithmic subsampling would suffer from degeneracy, which can lead to either a Gaussian limiting distribution with a faster convergence rate or a non-Gaussian limiting distribution menzel2017bootstrap, and invalidates the statistical inference based on standard multiway cluster-robust standard errors. On the other hand, we show that the multiway algorithmic subsampling allows for a non-degenerate asymptotic distribution regardless of whether the data are dependent or not. In other words, the algorithmic subsampling ensures a robustness against potential degeneracy, and even non-Gaussian degeneracy, of the asymptotic distribution, thereby allowing researchers to robustly enjoy valid statistical inference without knowing whether data are dependent or not. This finding about the additional practical advantage of the algorithmic subsampling is novel in the literature to our knowledge. We emphasize that these advantages of robustness come at the cost of efficiency and power loss due to the algorithmic subsampling.
Once these basic asymptotic statistical theories are established, we apply them to common econometric frameworks. Specifically, we propose a multiway algorithmic subsampling generalized method of moments (GMM) estimator, and derive asymptotic theories for it, including the consistency, asymptotic normality, and consistent variance estimation. This multiway algorithmic subsampling GMM estimator also enjoys the aforementioned property of robustness against potential degeneracy. Likewise, we also propose a multiway algorithmic subsampling M-estimator, and derive similar asymptotic theories for it, including the consistency, asymptotic normality, and consistent variance estimation.
{\bf Relation to the Literature:} This paper intersects with two branches of the literature, namely the algorithmic subsampling and multiway clustering. In econometrics, the algorithmic subsampling and its properties are first studied by \citet*{LeeNg2020ARE,LeeNg2020sketching}. This literature has focused on random (i.i.d.) sampling as emphasized earlier. Robust variances under multiway clustering have been proposed by \citet*{cameron2012robust}, \citet*{thompson2011simple}, and \citet*{cameron2014robust}. Asymptotic statistical properties under multiway clustering have been rigorously investigated by \citet*{davezies2019empirical}, \citet*{chiang2020inference}, \citet*{mackinnon2019wild}, \citet*{menzel2017bootstrap}, and \citet*{chiang2019multiway} under various contexts. This literature has not considered the algorithmic subsampling. To our knowledge, this paper is the first to study the properties of algorithmic subsampling under multiway clustering, and therefore, is the first to propose the aforementioned advantage of the algorithmic subsampling for robustness against potential degeneracy, including non-Gaussian degeneracy, in the asymptotic distribution under multiway clustering. We take advantage of the asymptotic distributional theory for incomplete one-sample U-statistics \citep*{Janson1984} to develop parts of our basic theoretical results. The method of algorithmic subsampling is also closely related to the general scheme of resampling methods for clustered data, which has been studied for one or multiway clustering by \citet*{mackinnon2017wild}, \citet*{mackinnon2018wild}, \citet*{djogbenou2019asymptotic}, \citet*{davezies2019empirical}, \citet*{mackinnon2019wild}, \citet*{chiang2020inference}, and \citet*{menzel2017bootstrap}, to name but a few.
{\bf Organization:} The rest of this paper is organized as follows. Section (ref) introduces the multiway algorithmic subsampling and presents its asymptotic statistical theories. Sections (ref) and (ref) demonstrate applications to the GMM and M-estimation frameworks, respectively. Section (ref) presents Monte Carlo simulation studies. Section (ref) presents an empirical application to scanner data. Section (ref) concludes. All mathematical proofs and details are collected in the appendix.
Suppose that a researcher observes data $\{W_{ij} : 1 \le i \le N, 1 \le j \le M\}$, where $N$ and $M$ are the sample sizes in the first and second cluster dimensions. For instance, $N$ and $M$ are the number of markets and the number of products, respectively, in scanner data. The data may be two-way clustered, in the sense that we allow for arbitrary statistical dependence of $W_{ij}$ across $j \in \{1,...,M\}$ within each market $i$ (due to a common demand shock) and also allow for arbitrary statistical dependence of $W_{ij}$ across $i \in \{1,...,N\}$ within each product $j$ (due to a common supply shock). A formal assumption of this sampling process will be stated as Assumption (ref) ahead. Using such two-way cluster sampled data, we are interested in (uniformly) consistent estimation of and statistical inference about $E[f(W_{ij})]$ based on standard econometric techniques. Since scanner data are very big, however, computational bottlenecks in time and/or memory may limit or even prohibit implementation of standard econometric analysis. \citet*{LeeNg2020ARE,LeeNg2020sketching} therefore suggest the algorithmic subsampling of big data to alleviate computational burdens.
Adapting the ideas of \citet*{LeeNg2020ARE,LeeNg2020sketching} to our framework of two-way clustered data, we propose the following multiway algorithmic subsampling procedure. (We remark that the algorithmic subsampling is different from the subsampling as a resampling method.) Let $p_{NM}$ denote the probability of subsample selection that may depend on the current sample size $(N,M)$. Generate i.i.d Bernoulli$\left(p_{NM}\right)$ random variables $\{Z_{ij}: 1 \le i \le N, 1 \le j \le M\}$ independently from data. Let $\widehat L = \sum_{i=1}^N \sum_{j=1}^M Z_{ij}$ denote the number of non-zero elements. Note that $\widehat L$ follows Binomial$\left(NM,p_{NM}\right)$, and thus $L \equiv E[\widehat L]=NMp_{NM}$ in particular. In fact, this formulation of the algorithmic subsampling is called the Bernoulli subsampling, and is one of the alternative approaches to subsampling proposed by \citet*{LeeNg2020ARE}. We focus on this Bernoulli subsampling in this paper for simplicity as well as its desired property of the aforementioned robustness against degeneracy. That said, we remark that it is also feasible to use alternative subsampling methods (namely, the uniform subsampling with and without replacement) proposed by \citet*{LeeNg2020ARE}.
We use $ \widehat L^{-1}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}f\left(W_{ij}\right) $ to estimate and make inference about $E[f(W_{ij})]$. To this end, we first develop the uniform weak law of large numbers and the central limit theorem under this setting of the multiway algorithmic subsampling in Sections (ref) and (ref), respectively. We then apply these basic theories in turn to establish the consistency and the asymptotic normality for the generalized method of moments and M-estimation in Sections (ref) and (ref), respectively. Hereafter for conciseness of notations, $p_{NM}$ will be abbreviated as $p$. We let $[k]$ denote the set $\{1,...,k\}$ for any $k \in \mathbb{N}$, and let $[k]^c=\mathbb{N}\backslash[k]$ for any $k\in \mathbb{N}.$ We use the short-hand notation $\underline{C}=\min\left\{N,M\right\}$. Throughout the paper, the asymptotics is understood as $\underline{C}\to \infty$. For a vector $v\in \mathbbm R^k$, let $\|v\|$ denote the Euclidean norm of $v$.
We first formally state the assumption of two-way cluster sampling.
Part (i) requires a form of the identical distribution condition in separate permutations of the $i$ index and the $j$ index. Although we relax the independent sampling, we maintain a form of the identical distribution as such. Part (ii) requires that sets of observations are independent if they do not share the same $i$ index or the same $j$ index, i.e., $\left(W_{ij}: i \in \{1,...,c_1\},j \in \{1,...,c_2\}\right)$ and $\left(W_{ij}: i \in \{c_1+1,...\},j \in \{c_2+1,...\}\right)$ are assumed to be independent for any $(c_1,c_2)\in \mathbb{N}^2$. However, any observations are sharing either the same $i$ index or the same $j$ index, then they are allowed to be arbitrarily dependent. For example, in the scanner data, two observations in the same market $i$ may be dependent due to a common demand shock, and likewise two observations in the same product $j$ may also be dependent due to a common supply shock.
Assumption (ref) consists of a sufficient condition for what we actually need. These conditions can be relaxed to the assumption that the data $(W_{ij})_{(i,j)\in\mathbb{N}^2}$ are generated via the process $W_{ij} = f(\alpha_i,\beta_j,\varepsilon_{ij})$ for some Borel-measurable function $f$, where $(\alpha_i)_{i \in \mathbb{N}}$, $(\beta_j)_{j \in \mathbb{N}}$, and $(\varepsilon_{ij})_{(i,j) \in \mathbb{N}^2}$ are mutually independent, and each of $(\alpha_i)_{i \in \mathbb{N}}$, $(\beta_j)_{j \in \mathbb{N}}$, and $(\varepsilon_{ij})_{(i,j) \in \mathbb{N}^2}$ is i.i.d. This data generating process, or so-called the Aldous-Hoover-Kallenberg representation, is implied by Assumption (ref). We can interpret $\alpha_i$ and $\beta_j$ as $i$- and $j$-specific effects, respectively, while $\varepsilon_{ij}$ is an idiosyncratic effect. This representation is also consistent with the data generating processes considered in the simulation studies.
For convenience of stating the next assumption, we introduce additional notations and definitions. Let $(T,d)$ be pseudometric space\footnote{That is, $d(x,y)=0$ does not imply $x=y$}. For $\varepsilon>0$, an $\varepsilon$-net of $T$ is a subset $T_\varepsilon$ of $T$ such that for every $t\in T$ there exists $t_\varepsilon\in T_\varepsilon$ with $d(t,t_\varepsilon)\le \varepsilon$. We define the $\varepsilon$-covering number $N(T,d,\varepsilon)$ of $T$ by
For any probability measure $Q$ on a measurable space $(S,\mathcal S)$ and any $q\ge 1$, define $\|f\|_{Q,q}=\left\{\int |f|^q dQ\right\}^{1/q}$ and let $L^q(S) = \{f:S\to \mathbbm R : \|f\|_{Q,q} < \infty\}$. A function $G:S\to\mathbbm R$ is an envelope of a class of functions $\mathcal G\ni g$, $g:S\to\mathbbm R$, if $\sup_{g\in\mathcal G}|g(s)|\le G(s)$ for all $s\in S$. With these notations and definitions, we state the following assumption regarding the function class where $f$ resides.
Part (i) is a location normalization (centering) and is therefore without loss of generality. Although this part will not be needed in the short run (Lemma (ref)), we state it here as this Assumption (ref) collects requirements about the function space where $f$ resides. Part (ii) is a regularity condition imposed to establish a uniform weak law of large numbers. Part (iii) is a technical requirement that is used to avoid measurability issues. At this moment, we are stating these high-level assumptions for the sake of generality, but we will provide the standard lower-level primitive conditions in the contexts of the application to the generalized method of moments presented in Section (ref) and the application to the M-estimation presented in Section (ref).
Under these assumptions, we can establish the uniform weak law of large numbers for multiway algorithmic subsample means as formally stated in the lemma below.
This result is not very surprising, but we state above as Lemma (ref) and prove it (in Appendix (ref)) for the following two reasons. First, this is nonetheless the first time it is stated and proved in the literature to the best of our knowledge. Second, more importantly, this lemma serves as a useful auxiliary device for other results to be presented ahead that are practically more relevant.
To establish the central limit theorem under the multiway algorithmic subsampling, we augment our assumptions with the following additional condition. Recall the short-hand notation $\underline{C}=\min\left\{N,M\right\}$.
It entails that there exist constants $\lambda_1 \geq 0$ and $\lambda_2\geq 0$ such that $\underline{C}/N\rightarrow \lambda_1$, $\underline{C}/M\rightarrow \lambda_2.$
To facilitate the subsequent discussions, we introduce a notion of degenerate asymptotic distribution. For any scalar-valued sequence of random variables $\left(X_{ij}\right)_{(i,j)\in \mathbb{N}^2}$ satisfying Assumption (ref), we say the asymptotic distribution is degenerate if $\operatorname{\text{Var}}\left( (\sqrt{\underline C} /NM)\sum_{i=1}^N \sum_{j=1}^M X_{ij}\right)\to 0$ as $\underline C\to\infty$. The following theorem establishes the central limit theorem under the multiway algorithmic subsampling.
A proof is provided in Appendix (ref). This central limit theorem entails a novel and useful feature of the algorithmic subsampling for multiway clustered data in practice. Notably, the asymptotic variance $\Gamma$ consists of a sum of two components, $\Gamma_A$ and $\Lambda \Gamma_B$. This is in contrast with the algorithmic subsampling for independent data, where the asymptotic variance consists of only one term. The first part, $\Gamma_A$, in fact coincides with the asymptotic variance that we would get without algorithmic subsampling. Specifically, $\Gamma = \Gamma_A$ if $p = 1$. More generally, if $p$ is a constant, as the sample sizes $(N,M)$ increase, then $\Lambda = 0$ so that $\Gamma = \Gamma_A$. On the other hand, if $p$ is chosen so that $\Lambda = \lim_{N,M \rightarrow \infty} (\underline{C}/NM)((1-p)/p) > 0$, then the second part, $\Lambda \Gamma_B$, is also present. Furthermore, $\Gamma_B$ is nonzero whenever the distribution of $f(W_{ij})$ is non-degenerate, and this feature provides a practically useful property of the robustness in inference against possible events of no cross sectional dependence.
In practice, a researcher may not ex ante know whether data exhibit cross sectional dependence ($E[f(W_{11})f^T(W_{12})] \neq 0$ or $E[f(W_{11})f^T(W_{21})] \neq 0$) or not. If a researcher knew the true dependence structure, he or she could set the correct cluster dimension to conduct valid inference. However, this premise is implausible. In case where there is no cross sectional dependence, then $\Gamma_A=0$ and the statistical inference based on the asymptotic normality without the algorithmic subsampling would suffer from the degeneracy problem. Because of the algorithmic subsampling, however, we can robustly safeguard against such degenerate asymptotic distributions without requiring a prior knowledge of the researcher about the presence/absence of cross sectional dependence in data. This result is novel in the literature, and also uncovers an additional useful property of the algorithmic subsampling in practice.\footnote{Indeed, the method of inference by \citet*{mackinnon2019wild} as well as \citet*{cameron2012robust} adapts to specific classes of degenerate asymptotic distributions. However, these restrict to the cases of Gaussian degeneracy, where the convergence rate is $\sqrt{NM}$ yet the asymptotic distribution is still Gaussian. On the other hand, these existing methods of inference by \citet*{cameron2012robust} and \citet*{mackinnon2019wild} do not adapt to the class of non-Gaussian degenerate asymptotic distributions.} Simulation studies presented in Section (ref) support this practically relevant property of the multiway algorithmic subsampling.
Intuitively, the algorithmic subsampling with smaller $p$ makes it less likely that multiple observations from the same row $i$ or same column $j$ are selected. Thus, it results in placing relatively more weights on the variance $E[f(W_{ij}) f(W_{ij})]$ than on the covariances, $E[f(W_{ij}) f(W_{ij'})]$ and $E[f(W_{ij}) f(W_{i'j})]$, in the asymptotic distribution. Hence, the part $E[f(W_{ij}) f(W_{ij})]$ of the asymptotic variance becomes dominant in the case of degenerate covariances, and this feature of the algorithmic subsampling prevents the degeneracy problem.
We can apply these theoretical results to a number of common frameworks of econometric analysis. Two of the most frequently used classes of econometric methods are the generalized method of moments (GMM) and the M-estimation. Therefore, we will demonstrate applications of these basic theories of the uniform weak law of large numbers (Lemma (ref)) and the central limit theorem (Theorem (ref)) to establish the consistency and the asymptotic normality of the GMM and M-estimators under the multiway algorithmic subsampling in Sections (ref) and (ref).
We conclude this section with a remark on alternative subsampling methods. As mentioned earlier, our method is based on the Bernoulli subsampling, and is one of the alternative approaches to subsampling proposed by \citet*{LeeNg2020ARE}. Besides the Bernoulli subsampling on which we focus in this paper, they propose the uniform subsampling with replacement, the uniform subsampling without replacement, and the leverage score subsampling. Among these alternative methods, it is also feasible to use the uniform subsampling with replacement and the uniform subsampling without replacement. Similar asymptotic properties will follow through similar lines of the argument following Janson1984 to those in the proof of Theorem (ref). See Appendix (ref) for details.
This section demonstrates an application of the basic theories to the ordinary least squares (OLS) estimator. Consider the linear regression model
where $Y_{ij}$ is a response variable, $X_{ij}$ is a vector of $d$ covariates and $u_{ij}$ is an error satisfying $E\left[u_{ij}|X_{ij}\right]=0.$ Let $W_{ij}=\left(Y_{ij},X_{ij}^T\right)^T$, and we apply the proposed multiway algorithmic subsampling to $W_{ij}$. The parameter of interest is the vector of linear projection coefficients
and the multiway algorithmic subsampling OLS estimator is
Applications of Lemma (ref) and Theorem (ref) yield the following limit distribution property about $\widehat\beta$.
See Appendix (ref) for a proof of this corollary.
In this section, we apply the basic methods and theories presented in Section (ref) to the multiway algorithmic subsampling generalized method of moments (GMM). Suppose that an economic model implies moment restrictions $E\left[g\left(W_{ij},\theta^0\right)\right]=0$ for a true parameter vector $\theta^0=(\theta_1^0,...,\theta_k^0)^T \in \Theta$, $\Theta \subset \mathbb{R}^k,$ where $g =\left(g_1 ,...,g_m \right)^T$, and $m\geq k$. With a Bernoulli sample $\{Z_{ij} : 1 \le i \le N, 1 \le j \le M\}$, the algorithmic subsample moment evaluated at $\theta=(\theta_1,...,\theta_k)^T \in \Theta$ is given by $\widehat g_{NM}\left(\theta\right)=\widehat L^{-1}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}g\left(W_{ij},\theta\right)$. Let $\widehat V$ be a positive semi-definite random matrix, which may depend on $\theta$. We define the multiway algorithmic subsampling GMM estimator $\widehat\theta$ as the solution to $$ \max_{\theta \in \Theta}\widehat Q_{NM}\left(\theta\right), $$ where $\widehat Q_{NM}\left(\theta\right)=-\widehat g_{NM}\left( \theta\right)^T\widehat V\widehat g_{NM}\left(\theta\right).$ The true parameter vector $\theta^0\in \Theta$ is assumed to uniquely solve the population problem $\max_{\theta \in \Theta}-E[g(W_{ij},\theta)]^TVE[g(W_{ij},\theta)],$ where $V$ is positive semi-definite and $\widehat V \overset{P}{\rightarrow}V$.
To establish the consistency and asymptotic normality for the multiway algorithmic subsampling GMM estimator $\widehat\theta$, we make the following assumption. For concisely stating the following assumption, we introduce one additional definition regarding Lipschitz continuity. A function $g: \mathbb{R}^k \to\mathbbm R,$ is Lipschitz with a universal Lipschitz constant, if there exists a positive constant $M$ such that $\left|g\left(w,\theta\right)-g\left(w,\theta'\right)\right|\leq M \left\|\theta-\theta'\right\|$ for all $w\in\rm{supp}(W_{ij})$.
Assumption (ref) is analogous to the conditions required for Theorem 2.6 and Theorem 3.4 in \citet*{Newey1994}, which state the consistency and asymptotic normality, respectively, of the GMM estimator under the conventional random sampling.
We first state the consistency of the multiway algorithmic subsampling GMM estimator.
A proof is provided in Appendix (ref). It follows from combining the arguments in the proofs of Newey1994 with our uniform weak law of large numbers for the multiway algorithmic subsampling (Lemma (ref)) presented in Section (ref).
We next state the asymptotic normality of the multiway algorithmic subsampling GMM estimator.
A proof is provided in Appendix (ref). It follows from combining the arguments in the proofs of Newey1994 with our central limit theorem for the multiway algorithmic subsampling (Theorem (ref)) presented in Section (ref).
The components, $G$ and $\Omega$ in Theorem (ref), of the asymptotic variance of the multiway algorithmic subsampling GMM estimator can be estimated by $$ \widetilde G=\frac{1}{\widehat L}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}\nabla_{\theta}g\left(W_{ij},\widehat \theta\right) $$ and $$ \widetilde \Omega=\widetilde \Gamma_1+\Lambda \widetilde \Gamma_2, $$ respectively, where $$ \widetilde \Gamma_1=\frac{\underline{C}}{\widehat L^2}\sum_{i=1}^{N}\sum_{1\leq j,j'\leq M}Z_{ij}Z_{ij'}g\left(W_{ij},\widehat \theta\right)g^T\left(W_{ij'},\widehat \theta\right)+ \frac{\underline{C}}{\widehat L^2}\sum_{1\leq i, i'\leq N}\sum_{j=1}^MZ_{ij}Z_{i'j}g\left(W_{ij},\widehat \theta\right)g^T\left(W_{i'j},\widehat \theta\right) $$ and $$ \widetilde \Gamma_2= \frac{1}{\widehat L}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}g\left(W_{ij},\widehat \theta\right)g^T\left(W_{ij},\widehat \theta\right). $$
We propose to estimate the asymptotic variance $\left(G^TVG\right)^{-1}G^TV\Omega VG\left(G^TVG\right)^{-1}$ by the sample counterpart $\left(\widetilde G^T\widehat V\widetilde G\right)^{-1}\widetilde G^T\widehat V\widetilde \Omega \widehat V\widetilde G\left(\widetilde G^T\widehat V\widetilde G\right)^{-1}.$ To guarantee that this algorithmic subsampling variance estimator works asymptotically, we make the following assumption in addition.
With this additional assumption, $\left(\widetilde G^T\widehat V\widetilde G\right)^{-1}\widetilde G^T\widehat V\widetilde \Omega \widehat V\widetilde G\left(\widetilde G^T\widehat V\widetilde G\right)^{-1}$ is consistent for the asymptotic variance $\left(G^TVG\right)^{-1}G^TV\Omega VG\left(G^TVG\right)^{-1}$, as formally stated in the following theorem.
A proof is provided in Appendix (ref). It follows by combining Lemma (ref) and similar lines of arguments to those in the proofs of Lemma (ref) and Theorem (ref).
In this section, we apply the basic methods and theories presented in Section (ref) to the multiway algorithmic subsampling M-estimation. Let $\Theta \subset \mathbb{R}^k$ be a parameter space and define the class $\mathcal{Q}=\left\{q\left(\cdot, \theta\right): \theta \in \Theta\right\}$ of functions $q(\cdot,\theta)$ indexed by $\theta$. With a Bernoulli sample $\{Z_{ij}, 1 \leq i \leq N, 1 \leq j \leq M\}$, we define the multiway algorithmic subsampling M-estimator $\widehat \theta$ as the solution to $$ \max_{\theta \in \Theta}-\frac{1}{\widehat L}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}q\left(W_{ij},\theta\right). $$ The true parameter vector $\theta^0=(\theta_1^0,...,\theta_k^0)^T \in \Theta$ is assumed to uniquely solve the population maximization problem $\max_{\theta \in \Theta}-E\left[q\left(W_{ij},\theta\right)\right]$, in the sense that $E\left[q\left(W_{ij},\theta^0\right)\right]<E\left[q\left(W_{ij},\theta\right)\right]$ holds for all $\theta=(\theta_1,...,\theta_k)^T \in \Theta$ and $\theta \neq \theta^0$. For each $\theta \in \Theta$, let $-\widehat L^{-1}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}q\left(W_{ij}, \theta\right)$ and $-E\left[q\left(W_{ij},\theta\right)\right]$ be denoted by $\widehat Q_{NM}\left(\theta\right)$ and $Q_0\left(\theta\right)$, respectively, for conciseness.
To establish the consistency and asymptotic normality for the multiway algorithmic subsampling M-estimator $\widehat\theta$, we make the following assumption.
Assumption (ref) is analogous to the conditions required for Theorem 2.1 and Theorem 3.1 in \citet*{Newey1994}, which state the consistency and asymptotic normality, respectively, of the M-estimator under the conventional random sampling.
We first state the consistency of the multiway algorithmic subsampling M-estimator.
A proof is provided in Appendix (ref). It follows from combining the arguments in the proof of Newey1994 with our uniform weak law of large numbers for the multiway algorithmic subsampling (Lemma (ref)) presented in Section (ref).
We next state the asymptotic normality of the multiway algorithmic subsampling M-estimator.
A proof is provided in Appendix (ref). It follows from combining the arguments in the proof of Newey1994 with our central limit theorem for the multiway algorithmic subsampling (Theorem (ref)) presented in Section (ref).
The components, $H$ and $\Sigma$ in Theorem (ref), of the asymptotic variance of the multiway algorithmic subsampling M-estimator can be estimated by $$ \widetilde H=-\frac{1}{\widehat L}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}\nabla_{\theta\theta^T}q\left(W_{ij},\widehat \theta\right) $$ and $ \widetilde \Sigma=\widetilde \Sigma_1+\Lambda \widetilde \Sigma_2, $ respectively, where
and $$ \widetilde \Sigma_2= \frac{1}{\widehat L}\sum_{i=1}^{N}\sum_{j=1}^{M}Z_{ij}\nabla_{\theta}q\left(W_{ij},\widehat \theta\right)\nabla_{\theta}q\left(W_{ij},\widehat \theta\right)^T. $$ Thus, we propose to estimate $H^{-1}\Sigma H^{-1}$ by the sample counterpart $\widetilde H^{-1}\widetilde \Sigma \widetilde H^{-1}.$ To guarantee that this asymptotic variance estimator works, we use the following assumption in addition.
With this additional assumption, $\widetilde H^{-1}\widetilde \Sigma \widetilde H^{-1}$ is consistent for the asymptotic variance $H^{-1}\Sigma H^{-1}$, as formally stated in the following theorem.
A proof is provided in Appendix (ref). It follows by combining Lemma (ref) and similar lines of arguments to those in the proofs of Lemma (ref) and Theorem (ref).
As emphasized in Section (ref), we discovered a new advantage of the algorithmic subsampling that it allows for robustness in inference against potential degeneracy of the asymptotic distribution under multiway clustering. In this section, we use Monte Carlo simulations to demonstrate this robustness property. Following menzel2017bootstrap, we consider two broad categories of designs, namely additively separable designs (Section (ref)) and nonseparable designs (Section (ref)). For each of these two broad categories, we experiment with a design that leads to a non-degenerate asymptotic distribution and another design that leads to a degenerate asymptotic distribution if the algorithmic subsampling were not to be employed. In total, we consider four designs. The multiway algorithmic subsampling will be shown to yield more accurate finite sample coverage results than conventional methods robustly across all the four cases, thereby supporting the aforementioned theoretical discovery by this paper.
First, we generate the two-way clustered array $\{Y_{ij}\}_{i \in [N], j \in [M]}$ according to the additively separable model
where $\beta_j$ and $\varepsilon_{ij}$ are i.i.d. standard normal, and $\alpha_i = (\zeta_i - \mu_\zeta)/\sigma_\zeta$ for $\log(\zeta_i) \stackrel{\text{i.i.d.}}{\sim} N(0,1)$, $\mu_\zeta = E[\zeta_i]$, and $\sigma_\zeta^2 = \text{Var}(\zeta_i)$. With this basic setup, we consider two designs:
Note that Design 2, without $i$-specific randomness or $j$-specific randomness, would lead to a degenerate asymptotic distribution if the algorithmic subsampling were not employed.
Table (ref) reports simulation results for $N=M=40$, $80$, $160$, $320$, and $640$. The top panel reports results for Design 1 (non-degenerate case), and the bottom panel reports results for Design 2 (degenerate case). Each panel contains results based on no algorithmic subsampling (i.e., $p=1$)\footnote{The 95% coverage is computed based on our asymptotic variance formula as the special case with $p=1$.} and results based on the algorithmic subsampling (with the subsampling probabilities of $p=1\underline{C}/(NM)$ and $p=2\underline{C}/(NM)$) for estimation of the mean. The asymptotic variance is estimated using a random subsample of ten percent of the sample. The displayed statistics are the bias (Bias), the standard deviation (SD), the root mean square error (RMSE), and the 95% coverage (95% Cover).
Observe that the 95% coverage frequencies are closer to the nominal probability of 95% with a use of the algorithmic subsampling than without a use of it. This observation is robustly true in both Design 1 (non-degenerate case) and Design 2 (degenerate case). For Design 2 or the degenerate case, in particular, the coverage frequency moves away from the nominal probability as the sample size increases if the algorithmic subsampling were not used. On the other hand, the coverage frequency approaches the nominal probability as the sample size increase if the algorithmic subsampling is used. These results demonstrate the aforementioned robustness property of the multiway algorithmic subsampling against potential degeneracy of the asymptotic distribution. We also experimented with additional simulation settings with much larger $N$ and $M$ and other subsampling probabilities for the algorithmic subsampling variance estimation, but we observe the same qualitative patterns in the results under these alternative settings.
On the one hand, $p=1$ leads to more precision, as quantified by smaller RMSE. On the other hand, $p=1$ leads to larger coverage as observed above. These two phenomena may appear contradictory at first glance. The relevant issues are with the variance estimation, and not with the point estimates. These results precisely highlight the cases of degeneracy. The asymptotic normality with the $\sqrt{\underline{C}}$-rate fails under the degeneracy if we do not use the algorithmic subsampling, i.e., if $p=1$. Therefore, the standard errors are misleadingly larger compared to the actual RMSE of the estimator and the simulated coverage rates exceed the nominal coverage probability in the degenerate case with $p=1$. This is the main reason why we propose to use the algorithmic subsampling (i.e., $p < 1$) to have inference with estimated variance robust against the degeneracy.
Second, we generate the two-way clustered array $\{Y_{ij}\}_{i \in [N], j \in [M]}$ according to the non-additive model
where $\alpha_i$, $\beta_j$ and $\varepsilon_{ij}$ are i.i.d. standard normal. With this basic setup, we consider two designs:
Note that Design 4 would lead to a degenerate asymptotic distribution that is a Gaussian chaos, which is non-Gaussian menzel2017bootstrap, if the algorithmic subsampling were not employed.
Table (ref) reports simulation results for $N=M=40$, $80$, $160$, $320$, and $640$. The top panel reports results for Design 3 (non-degenerate case), and the bottom panel reports results for Design 4 (degenerate case). Each panel contains results based on no algorithmic subsampling (i.e., $p=1$) and results based on the algorithmic subsampling (with the subsampling probabilities of $p=1\underline{C}/(NM)$ and $p=2\underline{C}/(NM)$) for estimation of the mean. The asymptotic variance is estimated using a random subsample of ten percent of the sample. The displayed statistics are the bias (Bias), the standard deviation (SD), the root mean square error (RMSE), and the 95% coverage (95% Cover).
Similarly to the case with the additively separable design, observe that the 95% coverage frequencies are closer to the nominal probability of 95% with a use of the algorithmic subsampling than without a use of it. This observation is robustly true in both Design 3 (non-degenerate case) and Design 4 (degenerate case). For Design 4 or the degenerate case, in particular, the coverage frequency moves away from the nominal probability as the sample size increases if the algorithmic subsampling were not used. On the other hand, the coverage frequency approaches the nominal probability as the sample size increase if the algorithmic subsampling is used. As before, these results demonstrate the aforementioned robustness property of the multiway algorithmic subsampling against potential degeneracy of the asymptotic distribution. We also experimented with additional simulation settings with much larger $N$ and $M$ and other subsampling probabilities for the algorithmic subsampling variance estimation, but we observe the same qualitative patterns in the results under these alternative settings.
In this section, we demonstrate an application of our proposed method to an analysis of demand for differentiated products using scanner data from the Dominick's Finer Foods (DFF) retail chain.\footnote{We thank James M. Kilts Center, University of Chicago Booth School of Business for allowing us to use this data set. It is available at https://www.chicagobooth.edu/research/kilts/datasets/dominicks.} Scanner data may be subject to two-way cluster dependence, as mentioned in Section (ref). Specifically, common demand shocks within a market may induce statistical dependence among different products within that market. Similarly, common supply shocks by a producer may induce statistical dependence among different markets within the product produced by that producer. In this light, a researcher would like to use a two-way cluster robust variance estimate for inference about the model parameters. However, the scanner data from the Dominick's Finer Foods (DFF) retail chain are too large, and today's computational resources will not permit the two-way cluster robust variance estimation in reasonable lengths of time. A simple way to overcome this problem is to use the full sample for parameter estimation and to use a subsample for variance estimation, but this approach fails to deliver robustly valid inference. Hence, we use our proposed multiway algorithmic subsampling method for estimation and two-way cluster robust inference about the key demand model parameter.
Following the literature nevo2000practitioner on analysis of demand for differentiated products with an additive Type-I-Extreme-Value error, we use the GMM approach with the moment restriction
where $i$ indexes products (universal product code, hereafter referred to as UPC), $j$ indexes markets (store $\times$ week), $S_{ij}$ denotes the share of product $i$ in market $j$, $P_{ij}$ denotes the price, $X_{ij}$ denotes a vector of controls (the UPC fixed effects and a time trend), $\zeta_{ij}$ denotes instruments, and $W_{ij} = (S_{ij},P_{ij},X_{ij}^T,\zeta_{ij}^T)^T$.\footnote{In case where the model involves product fixed effects, the algorithmic subsampling can be applied to within-transformation. This operation incurs additional computational costs, although this is a common issue in fixed-effect methods in general. In case a model involves two-way fixed effects, two-way differencing may induce a more complicated dependence structure especially under unbalanced panels. An alternative approach may be to use instrumental variables. We leave rigorous treatments of such a variety of extensions to fixed-effect models for future research.} In addition to the elements in $X_{ij}$, the instrument vector includes $\zeta_{ij}$ as an excluded variable the wholesale costs, which are calculated by inverting the gross margin. We drop those observations for which $\ln(S_{ij}) - \ln(S_{0j})$ is not finite,\footnote{In other words, we drop observations with the zero market share. Dropping these observations may generally incur a trimming bias. We adopt this trimming as it is a standard practice in the literature of demand analysis for differentiated products markets, and we consider the possibly biased estimand as our pseudo-true value.} as well as those observations with missing values. The parameter vector in the model consists of $\theta = (\theta_1,\theta_{-1}^T)^T$, and we are in particular interested in the price coefficient $\theta_1$.
We consider four product categories: beer, oats, snacks, and canned tuna. Table (ref) summarizes the sizes of the original data in terms of various dimensions. It first shows the number of UPCs, the number of weeks, and the number of stores for each product category. As we define a product as that identified by the UPC, the number of products $N$ coincides with the number of UPCs. We define a market as the unique combination of the week and the store. Therefore, the number of markets $M$ is close to, but is generally smaller than, the product of the number of weeks and the number of stores. It is smaller than the na\"ive product because of the unbalancedness in data. Finally, the bottom row shows the total number of observations, which is again smaller than the na\"ive product $NM$ because of the unbalancedness in data.
We now apply our multiway algorithmic subsampling GMM with the moment function defined in (ref) for each of the four product categories. Table (ref) summarizes the estimation results. The table displays the probability $p$ of algorithmic subsampling, the corresponding estimates and their standard errors for the price coefficient, and computational time in seconds for each of parameter estimation and asymptotic variance estimation.
First, observe that the estimates of the price coefficient are negative, as expected, and are statistically significant at the level of 95% for each column except for tuna despite efficiency loss due to the algorithmic subsampling and despite the two-way cluster robustness in the asymptotic variance. As emphasized in Sections (ref) and (ref), the algorithmic subsampling with $p \propto \underline{C}/(NM)$ allows these standard errors to have asymptotically accurate coverage robustly against potential degeneracy, unlike the conventional two-way cluster robust standard errors without the algorithmic subsampling.
Second, the computational time for parameter estimation is within about a dozen of seconds for each column, given that the algorithmic subsampling extracts only the proportions, $p \approx 0.003--0.009$, of the original sample sizes. However, it is the asymptotic variance estimation that costs more computational time under multiway cluster dependence. Focusing on the beer product category, for instance, even the algorithmic subsampling that extracts only the $p \approx 0.004$ portion of the original sample size requires 1223 seconds of computation for variance estimation. When the proportion doubles to $p \approx 0.009$, then the computational time nearly quadruples to 4458 seconds. A na\"ive calculation implies that the use of the full sample without the algorithmic subsampling would require about three years.
In this paper, we propose a novel method of algorithmic subsampling for multiway cluster dependent data. We develop asymptotic statistical properties of this proposed method. Specifically, we develop a new uniform weak law of large numbers and a new central limit theorem for the multiway algorithmic subsample means. As a consequence of the new central limit theorem, we show that the algorithmic subsampling allows for robustness against potential degeneracy of the asymptotic distribution under multiway clustering at the cost of efficiency and power loss due to the algorithmic subsampling. Applying these basic asymptotic statistical theories, we derive the consistency and the asymptotic normality for the multiway algorithmic subsampling generalized method of moments estimator and for the multiway algorithmic subsampling M-estimator.
Our main finding that the algorithmic subsampling allows for the robustness against degeneracy in the asymptotic distribution is novel in the literature on multiway clustering. Indeed, the method of inference by \citet*{mackinnon2019wild} as well as \citet*{cameron2012robust} adapts to the Gaussian degeneracy. However, these existing methods do not adapt to the class of non-Gaussian degenerate asymptotic distributions. In contrast, the asymptotic distribution under the algorithmic subsampling adapts even to the non-Gaussian degeneracy as well. The bootstrap method of \citet*{menzel2017bootstrap} is robust against the non-Gaussian degeneracy. Our proposed method via the algorithmic subsampling leads to the exact limit distribution, and thus non-conservative inference, unlike the method of \citet*{menzel2017bootstrap}. With these said, we once again emphasize that these merits come at the cost of efficiency and power loss by disposing parts of big data.
Finally, we shed some light on possible future directions. In this paper we consider non-nested multiway clustering \citep*[as in][]{cameron2012robust}. In practice, the researcher may be interested in applications with nested clustering in one or more cluster dimensions. Under the current framework, one could take the coarsest levels of clustering. Handling it in a more efficient way is a useful topic but is out of the scope of this paper. In addition, in \citet*{mackinnon2020testing}, formal theory is developed for testing the correct level of (one-way) clustering. One could consider to generalize such test for multiway nested clustering, which is also left for future research.