The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
46,736 characters
Lasso under multi-way clustering: Estimation and Post-selection Inference
\title[Lasso under multi-way clustering]{Lasso under multi-way clustering:\\
Estimation and Post-selection Inference }
\author[Chiang and Sasaki]{Harold Chiang and Yuya Sasaki}\address[Harold Chiang]{Department of Economics, Vanderbilt University, United States}\address[Yuya Sasaki]{Department of Economics, Vanderbilt University, United States}
\date{First arXiv version: May 6, 2019.
\\
\indent Code files are available upon request from the authors.}
\keywords{cluster robust standard errors, high dimensions, lasso, machine learning, multi-way clustering, post-selection inference.\\\indent \textit{JEL Classification}: C21, C55}
\maketitle
\begin{abstract}
This paper studies high-dimensional regression models with lasso when data is sampled under multi-way clustering.
First, we establish convergence rates for the lasso and post-lasso estimators.
Second, we propose a novel inference method based on a post-double-selection procedure and show its asymptotic validity.
Our procedure can be easily implemented with existing statistical packages.
Simulation results demonstrate that the proposed procedure works well in finite sample.
We illustrate the proposed method with a couple of empirical applications to development and growth economics.
\end{abstract}
\section{Introduction}\label{sec:introduction}
This paper studies a method of estimation and post-selection inference for regression parameters in high-dimensional linear models by lasso under multi-way clustering.
The objective is motivated by recently increasing demands from applied economic research.
On one hand, economists often use multi-way cluster sampled data.
Examples include, but are not limited to, network data, matched employer-employee data, matched student-teacher data, scanner data where observations are double-indexed by stores and products, market share data where observations are double-indexed by market and products, and growth/development data where observations are double-indexed by ethnicity and geographical units -- see Section \ref{sec:empirical illustration} for specific applications of the last example.
On the other hand, researchers also often use machine learning methods of estimation and inference for high-dimensional models in today's big data environments.
There are a number of useful methods in the literature that deal with each of these two issues (multi-way clustering and high dimensionality) separately, but the existing methods do not seem to provide a solution to dealing with both of these practically relevant issues simultaneously.
In this light, we present lasso under multi-way clustering, and propose a post-selection inference method for regression parameters under this sampling assumption.
In the important branch of the literature following the seminal work by \cite{BCCH12}, post-selection inference with lasso has been widely studied under various settings by \cite{BCH14}\footnote{See also \cite{BCH11}.}, \cite{JM14}, \cite{vdG14}, \cite{ZZ14}, \cite{BCK15}, and many others.
For empirical researchers, lasso has become a powerful machine learning tool under data-rich environments.
Most of the papers in this literature assume i.i.d. or independent sampling.
In many empirical applications, it is sometimes more plausible to assume multi-way cluster sampling (e.g., network data, matched employer-employee data, and matched student-teacher data).
Building upon \cite{BCCH12}, \cite{BC13} and \cite{BCH14}, this paper generalizes lasso and post-double-selection procedure by allowing for multi-way cluster sampling.
To our best knowledge, the present paper is the first in the literature of high-dimensional models to consider lasso under multi-way cluster sampling.
The influential work by \cite{CGM11} proposes multi-way cluster-robust inference methods for linear and nonlinear regression models -- also see \citet[][Section V]{CM15} for a survey.
Formal analysis of asymptotic properties and bootstrap validity under multi-way clustering is studied by \cite{menzel17} using the Aldous-Hoover representation -- see \citet[][Chapter 7]{kallenberg05} for example. Under the assumptions of separable exchangeability, the method of \cite{menzel17} covers both degenerate and non-degenerate cases.
Using the same representation, while focusing on the non-degenerate cases, \cite{DDG18} develop empirical process theory under multi-way cluster sampling which applies to a large class of econometric models.
Building upon the asymptotic framework of these two papers, \cite{MNW19} propose several wild bootstrap procedures for linear regression models, and examine their finite-sample performances under several different cluster sampling scenarios.
In this paper, we take advantage of the innovations by these preceding papers to develop a multi-way cluster-robust inference method for high-dimensional models.
To our best knowledge, the present paper is the first in this literature on multi-way clustering to consider high-dimensional models.
The rest of this paper is organized as follows.
Section \ref{sec:model} introduces the model.
Section \ref{sec:overview} presents an overview of the proposed methodology.
Section \ref{sec:asymptotic_theory} discusses a formal asymptotic theory.
Section \ref{sec:heterogeneous_cluster_sizes} presents an extension of the baseline results to cases of heterogeneous cluster sizes.
Section \ref{sec:simulation_studies} presents simulation studies.
Section \ref{sec:empirical illustration} presents an empirical illustration with development and growth economics.
Section \ref{sec:conclusion} concludes.
The appendix contains mathematical proofs and auxiliary lemmas.
\section{The Model}\label{sec:model}
Consider the high-dimensional regression model
\begin{align}
\label{eq:model_y}
Y_{ij}=&D_{ij}\alpha +X'_{ij}\beta + R^Y_{ij} + \varepsilon_{ij},\:
{\mathrm{E}}[\varepsilon_{ij}|D_{ij},X_{ij}]=0,
\end{align}
where
$Y_{ij}$ is an observed outcome variable,
$(D_{ij},X_{ij}')'$ is an observed vector of regressors, and
$R_{ij}^Y$ is an approximation error for the unit of observation with the double index $(i,j)$.
We set $\alpha$ as a scalar parameter of interest.
The dimension $p$ of the nuisance parameter vector $\beta \in \mathbb R^p$ is potentially increasing in the sample size.
Following the literature on high-dimensional post-selection inference \citep[e.g.,][]{BCH14}, we also consider the auxiliary projection
\begin{align}
\label{eq:model_d}
D_{ij}=&X_{ij}\gamma + R^D_{ij} + v_{ij},\:
{\mathrm{E}}[v_{ij}|X_{ij}]=0,
\end{align}
where $R^D$ is an approximation error.
The dimension $p$ of the nuisance parameter vector $\gamma\in \mathbb R^p$ is the same as that of $\beta$, and is potentially increasing in the sample size.
In the absence of two-way clustering, the system (\ref{eq:model_y})--(\ref{eq:model_d}) would be the same as the model considered in \cite{BCH14}.
We first consider two-way clustering where each cell contains one observation.
Section \ref{sec:heterogeneous_cluster_sizes} presents an extension to the case of heterogeneous cluster sizes.
\section{Overview of the Method}\label{sec:overview}
In this section, we present an overview of the proposed method, namely estimation and post-selection inference.
Formal theoretical justifications are discussed in Section \ref{sec:asymptotic_theory}.
A researcher observes a sample $\left\{\left. (Y_{ij},D_{ij},X_{ij}') \right\vert i \in \{1,...,N\}, j \in \{1,...,M\}\right\}$ of size $NM$.
The estimation procedure consists of two steps.
First, define the lasso estimates for (\ref{eq:model_y}) and (\ref{eq:model_d}) by
\begin{align}
(\widehat \alpha,\widehat \beta')'=& \underset{\alpha,\beta }{\rm argmin} \sum_{i=1}^{N}\sum_{j=1}^{M} (Y_{ij} - \alpha D_{ij} - X_{ij}'\beta)^2 +\lambda_1\|(\alpha,\beta')'\|_1,\label{eq:lasso_main}\\
\text{and}\qquad
\widehat \gamma=& \underset{\gamma }{\rm argmin} \sum_{i=1}^{N}\sum_{j=1}^{M} ( D_{ij} - X_{ij}'\gamma)^2 +\lambda_2 \|\gamma\|_1,\label{eq:lasso_second}
\end{align}
respectively, for some regularization parameters $\lambda_1$ and $\lambda_2$, valid choices of which are discussed in the statement of Theorem \ref{theorem:rates} ahead.
Denote the supports of the lasso estimates by $\widehat I_1={\rm support}(\widehat \beta)$ and $\widehat I_2={\rm support}(\widehat \gamma)$, and let $\widehat I= \widehat I_1 \cup \widehat I_2$.
In the second step, define the post-double-selection lasso estimate $\widetilde \alpha$ by
\begin{align}
(\widetilde \alpha,\widetilde \beta')= \underset{{\rm support}(\beta) \subset \widehat I}{\rm argmin} \sum_{i=1}^{N}\sum_{j=1}^{M} (Y_{ij} - \alpha D_{ij} - X_{ij}'\beta)^2.\label{eq:post_lasso}
\end{align}
Let $\underline C =N\wedge M$, $\mu_N=\underline C/N$ and $\mu_M=\underline C/M$.
Under suitable conditions to be formally stated in Section \ref{sec:asymptotic_theory}, we have the asymptotic normality
\begin{align*}
\sigma^{-1}\sqrt{\underline C}(\widetilde \alpha - \alpha) \leadsto N(0,1),
\end{align*}
where the asymptotic variance is given by $\sigma^2= Q^{-1}\Gamma Q^{-1}$ with
\begin{align*}
Q=&{\mathrm{E}}[v_{11}^2],
\\
\Gamma=&\bar\mu_N \Gamma_N + \bar\mu_M \Gamma_M
=\bar\mu_N{\mathrm{E}}[v_{11} \varepsilon_{11} v_{12}\varepsilon_{12} ] + \bar\mu_M {\mathrm{E}}[ v_{11}\varepsilon_{11}v_{21} \varepsilon_{21}],
\end{align*}
and $\bar\mu_N$ and $\bar\mu_M$ denoting the limits of $\mu_N$ and $\mu_M$, respectively.
The asymptotic variance is estimated by the sample counterpart $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$, where
\begin{align*}
\widehat Q=&\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M} \widehat v_{ij}^2,
\\
\widehat \Gamma
=& \frac{\underline C}{(NM)^2}\sum_{i=1}^{N}\sum_{1\le j,j'\le M} \widehat v_{ij} \widehat \varepsilon_{ij} \widehat \varepsilon_{ij'} \widehat v_{ij'}
+ \frac{\underline C}{(NM)^2}\sum_{1\le i,i'\le N}\sum_{j=1}^{M} \widehat v_{ij}\widehat \varepsilon_{ij} \widehat\varepsilon_{i'j}\widehat v_{i'j} ,
\end{align*}
$\widehat v_{ij}= D_{ij}- X_{ij}'\widehat \gamma$, and $\widehat \varepsilon_{ij}=Y_{ij}- \widehat \alpha D_{ij} - X_{ij}'\widehat\beta$.
In summary, we propose to report the post-double-selection lasso estimate $\widetilde\alpha$ as an estimate of $\alpha$ with its standard error given by $\widehat\sigma / \sqrt{\underline C}$.
The $\alpha^\ast$-level confidence interval can be constructed as $\left[\widetilde\alpha + \Phi^{-1}(\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C},\widetilde\alpha + \Phi^{-1}(1-\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C}\right]$, where $\Phi^{-1}$ denotes the quantile function of the standard normal distribution.
\section{Asymptotic Theory}\label{sec:asymptotic_theory}
The two-way sample sizes $(N,M) \in \mathbb{N}^2$ will be index by a single index $n \in \mathbb{N}$ as $(N,M) = (N(n),M(n))$ where $M(n)$ and $N(n)$ are non-decreasing in $n$ and $M(n)N(n)$ is increasing in $n$. For simplicity, each size of intersection $n_{ij}$ is assumed to be uniformly bounded by a positive integer $\bar n$ that is independent of $n$.
With this said, we will suppress the index notation and write $(N,M)$ for simplicity.
We fix a number of notations.
For each $n$, let ${\mathrm{P}}_{n}$ denote the law with respect to sample size $(N,M)$ -- note that we allow the dimension $p$ of $X_{ij}$ to grow with $n$.
Let $a := p\vee (NM)$.
Also recall the notations $\underline C =N\wedge M$, $\mu_N=\underline C/N$, and $\mu_M=\underline C/M$ from Section \ref{sec:overview}.
We use the short-hand notation $[k] = \{1,...,k\}$ and $[k]^c=\mathbb{N} \backslash [k]$ for any $k \in \mathbb{N}$.
For a sequence $(t_{ij})_{i\in[N],j\in [M]}$, denote $\|t_{ij}\|_n = \sqrt{\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} t_{ij}^2 }$.
Thus, $\|X'_{ij}\delta\|_n = \sqrt{\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} \delta' X_{ij}X_{ij}'\delta}$ is the prediction norm of $\delta$.
Let $\|A\|_\infty=\max_{k,l}|A_{k,l}|$ denote the max norm of matrix $A$.
We write $a \lesssim b$ to mean $a \leqslant cb$ for some $c > 0$ that does not depend on $n$.
We also write $a \lesssim_{\mathrm{P}} b$ to mean $a = O_P(b)$.
We write $Z_{ij}=(Y_{ij},D_{ij},X'_{ij})'$ for the $(p+2)$-dimensional random vector in data.
Throughout, we assume that this random vector $Z_{ij}$ is Borel measurable -- see \citet[][pp. 304]{kallenberg05}.
With these notations, we state the following four assumptions.
\begin{assumption}[Sampling]\label{a:sampling}
Suppose that $\underline C \to \infty $, $\mu_N \to \bar \mu_N \geqslant 0$, and $\mu_M \to \bar \mu_M \geqslant 0$.
\begin{enumerate}
\item $(Z_{ij})_{(i,j)\in \mathbbm N^2}$ is an infinite sequence of separately exchangeable $(p+2)$-dimensional random vectors.
That is, for any permutations $\pi_1$ and $\pi_2$ of $\mathbbm N$, we have
\begin{align*}
(Z_{ij})_{(i,j)\in \mathbbm N^2}\overset{d}{=} (Z_{\pi_1(i)\pi_2(j)})_{(i,j)\in \mathbbm N^2}.
\end{align*}
\item $(Z_{ij})_{(i,j)\in \mathbbm N^2}$ is dissociated.
That is, for any $(c_1,c_2)\in \mathbbm N^2$,
$
(Z_{ij})_{i \in [c_1], j \in [c_2]}
$
is independent of
$
(Z_{ij})_{i \in [c_1]^c, j \in [c_2]^c}.
$
\item For each $n$, an econometrician observes $(Z_{ij})_{i\in[N],j\in[M]}$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Moments]\label{a:moments}
There exists a sequence $\{B_n\}_{n=1}^\infty$ of positive constants such that the following conditions hold for all $n \in \mathbb{N}$ for some $q>4$:
\begin{enumerate}
\item
${\mathrm{E}}[|D_{11}|^{2q}]+\max_{k\in [p]}{\mathrm{E}}[|X_{11,k}|^{2q}]+ {\mathrm{E}}[|\varepsilon_{11}|^{2q}|X_{11},v_{11}]+{\mathrm{E}}[|v_{11}|^{2q}|X_{11}]\le K$ a.s.
and
$0 < c\le {\mathrm{E}}[v^{2}_{11}|X_{11}]$ a.s. for positive constants, $c$ and $K$, that are independent of $n$.
\item ${\mathrm{E}}[\|X_{11}\|_\infty^{2q}]\le B_{n}^{2q}$ and $B_n\sqrt{\log a}\lesssim (N\vee M)^{1/2-1/q}$.
\item $\bar \mu_N{\mathrm{E}}[v_{11} \varepsilon_{11} v_{12}\varepsilon_{12} ] + \bar \mu_M {\mathrm{E}}[ v_{11}\varepsilon_{11}v_{21} \varepsilon_{21}] $ and the maximal and minimal eigenvalues of ${\mathrm{E}}[X_{11}X_{11}']$ are bounded and bounded away from zero uniformly in $n$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Sparsity]\label{a:sparsity}
${}$
\begin{enumerate}
\item $\|\beta\|_0+ \|\gamma\|_0 \lesssim s$.
\item $\|R_{ij}^Y\|_n+ \|R_{ij}^D\|_n\le c_s \lesssim_{\mathrm{P}} \sqrt{s /\underline C}$.
\item $\frac{s^2_n(\log (a))^2}{\underline C}=o(1)$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Sparse Eigenvalues]\label{a:sparse_eigenvalues}
There exists a sequence $\{\ell_n\}$ such that $\ell_n \to \infty$ and, with probability at least $1-o(1)$,
\begin{align*}
0<c\le \phi_{{\rm min}}(s\ell_n)\le \phi_{{\rm max}}(s\ell_n)\le c' <\infty
\end{align*}
holds for some constants, $c$ and $c'$, that are independent of $n$, where
\begin{align*}
\phi_{{\rm max}}(m):=
\max_{\substack{1\le \|\delta\|_0 \le m }} \frac{\delta' M\delta}{\|\delta\|^2}
\text{ and } \phi_{{\rm min}}(m):= \min_{\substack{1\le \|\delta\|_0 \le m }} \frac{\delta' M\delta}{\|\delta\|^2},
\end{align*}
with $M=\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} X_{ij}X_{ij}'$, denote the maximal and minimal $m$-sparse eigenvalues.
\end{assumption}
\begin{remark}[Discussion of the Assumptions]\label{remark:discussion_assumptions}
Assumption \ref{a:sampling} is closely related to Assumption 1 of \cite{DDG18}.
The main difference is that we allow $p$ to be changing with $n$.
We remark that the exchangeability assumption is not new in econometrics -- it has been used in \cite{andrews05} and \cite{menzel15} as well as \citet{menzel17}, \citet{DDG18}, and \citet{MNW19}.
Assumption \ref{a:moments} is standard in the literature on post-selection inference with lasso.
Parts (1) and (2) require an existence of higher order moments of key objects.
Note that common assumptions in high-dimensional literature, such as sub-gaussianity or boundedness, are not required.
They can be replaced by some higher level conditions similar to Condition RF of \cite{BCCH12}.\footnote{See their Lemma 3.}
Part (3) of Assumption \ref{a:moments} requires that the asymptotic variance is bounded away from zero.\footnote{Similarly to \cite{DDG18}, we focus on non-degenerate cases in this paper. See \cite{menzel17} for the studies of degenerate cases using a bootstrap-based method.}
Assumption \ref{a:sparsity} is a direct generalization of Condition ASTE (iii) and (iv) of \cite{BCH14}.
Finally, Assumption \ref{a:sparse_eigenvalues} is analogous to Condition SE of \cite{BCH14}, which is standard in the high-dimensional literature.
It only imposes small diagonal submatrices to be well behaved.
$\triangle$
\end{remark}
\subsection{Independentization via H\'ajek Projection}
In this section, we show that an empirical process in multi-way clustered samples can be represented as a sum of independent variables via H\'ajek projection.
Furthermore, its variance can be shown to be approximated by covariances of observed variables.
For any $f:{\rm support}(Z)\to \mathbb R$, we let
\begin{align*}
\mathbb{G}_C f:=\sqrt{\underline C}\Big\{\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} f(Z_{ij}) - {\mathrm{E}}[f(Z_{11})]\Big\}
\end{align*}
denote its empirical process.
\begin{lemma}[Independentization via H\'ajek Projection]\label{lemma:hajek}
If Assumption \ref{a:sampling} holds and $f:{\rm support}(Z)\to \mathbb R$ satisfies ${\mathrm{E}} f^2(Z_{11})<K$ for a finite constant $K$ that is independent of $n$, then there exist i.i.d. uniform random variables $U_{i0}$ and $U_{0j}$ such that the H\'ajek projection $H_n f$ of $\mathbb{G}_C f$ on
$$
\mathcal{G}_n=\Big\{ \sum_{i=1}^{N} g_{i0}(U_{i0}) + \sum_{j=1}^{M} g_{0j}(U_{0j}) : g_{i0}, g_{0j} \in L^2({\mathrm{P}}_{n}) \Big\}
$$
is equal to
\begin{align*}
H_n f=\sum_{i=1}^{N} \frac{\sqrt{\underline C}}{N} {\mathrm{E}}\Big[f(Z_{i1})- {\mathrm{E}} f(Z_{11}) \Big| U_{i0}\Big] + \sum_{j=1}^{M}\frac{\sqrt{\underline C}}{M} {\mathrm{E}}\Big[f(Z_{1j})- {\mathrm{E}} f(Z_{11}) \Big| U_{0j}\Big]
\end{align*}
for each $n$.
Furthermore,
\begin{align*}
V(\mathbb{G}_C f)= V(H_n f)+O(\underline C^{-1})=\bar\mu_N Cov(f(Z_{11}),f(Z_{12})) + \bar\mu_M Cov(f(Z_{11}),f(Z_{21}))+O(\underline C^{-1})
\end{align*}
holds a.s.
\end{lemma}
A proof of this lemma can be found in Appendix \ref{sec:lemma:hajek}.
The first part of the lemma shows that an empirical process $\mathbb{G}_C f$ under multi-way cluster sampling can be represented as a sum of independent unobserved variables via H\'ajek projection $H_n f$.
While $U_{i0}$ and $U_{0j}$ are unobserved, the second part of this lemma in turn shows that the variance of the H\'ajek projection can be approximated by covariances of observed variables.
Note that, since $H_n f$ is a H\'ajek projection, the lemma implies $\frac{\mathbb{G}_C f}{\sqrt{V(\mathbb{G}_C f)}}=\frac{H_n f}{\sqrt{V(H_n f)}}+o_{\mathrm{P}}(1)$ if $\bar\mu_N Cov(f(Z_{11}),f(Z_{12})) + \bar\mu_M Cov(f(Z_{11}),f(Z_{21}))$ is bounded and bounded away from zero uniformly in $n$.
Our Lemma \ref{lemma:hajek} can be seen as an extension to Lemma D.2 in \cite{DDG18}.
Specifically, while \cite{DDG18} consider a fixed data generating process over the sample size $n$, our Lemma \ref{lemma:hajek} allows the data generating process to vary with $n$ in particular for the sake of accommodating the increasing of dimensionality $p$ for high-dimensional models.
The lemma serves as a main building block for all the asymptotic results to be presented ahead.
\subsection{Convergence Rates of Lasso and Post-Lasso under Multi-Way Clustering}
We next show the convergence rates of the lasso estimator $(\widehat\alpha,\widehat\beta',\widehat\gamma')'$ and the post-lasso estimator $(\widetilde\alpha,\widetilde\beta',\widetilde\gamma')'$ under multi-way clustering.
\begin{theorem}[Convergence Rates for Lasso and Post-Lasso under Multi-Way Clustering]\label{theorem:rates}
If Assumptions \ref{a:sampling}, \ref{a:moments} (1)--(2), \ref{a:sparsity} (1)--(2), and \ref{a:sparse_eigenvalues} are satisfied, and $\lambda_1, \lambda_2 =C\sqrt{(NM)^2 \log a/\underline C}$ for some constant $C>1$, then
\begin{align*}
&\|\widehat\eta-\eta\|_1 + \|\widehat \gamma - \gamma\|_1 \lesssim \sqrt{\frac{s^2\log a}{\underline C}},\, \qquad \|W_{ij}'(\widehat\eta-\eta)\|_n + \|X_{ij}(\widehat \gamma - \gamma)\|_n\lesssim \sqrt{\frac{s\log a}{\underline C}},\\
&\|\widetilde\eta-\eta\|_1 + \|\widetilde \gamma - \gamma\|_1 \lesssim \sqrt{\frac{s^2\log a}{\underline C}},\, \qquad \|W_{ij}'(\widetilde\eta-\eta)\|_n + \|X_{ij}(\widetilde \gamma - \gamma)\|_n\lesssim \sqrt{\frac{s\log a}{\underline C}}, \quad\text{and}\\
&\|\widehat\eta-\eta\| + \|\widehat \gamma - \gamma\| + \|\widetilde\eta-\eta\| + \|\widetilde \gamma - \gamma\| \lesssim \sqrt{\frac{s\log a}{\underline C}}
\end{align*}
hold, where $W_{ij}=[D_{ij},X_{ij}']'$ and $\eta=(\alpha,\beta')'$.
\end{theorem}
A proof can be found in Appendix \ref{sec:theorem:rates}, and is based on the previous result (Lemma \ref{lemma:hajek}).
In the multi-way sampling, this lemma can be viewed as a counterpart of Lemma 6 and Lemma 7 in \cite{BCCH12}.
\subsection{Post-Selection-Inference with Post-Lasso under Multi-way Clustering}\label{sec:post_selection_inference}
In this section, we present the main result of this paper.
The limit normal distribution of the post-double-selection lasso estimate $\widetilde\alpha$ is established based on the previous two results (Lemma \ref{lemma:hajek} and Theorem \ref{theorem:rates}).
\begin{theorem}[Asymptotic Normality]\label{theorem:asymptotic_normality}
If Assumptions \ref{a:sampling}, \ref{a:moments}, \ref{a:sparsity} and \ref{a:sparse_eigenvalues} are satisfied, and $\lambda_1$ and $\lambda_2$ are chosen according to the statement of Theorem \ref{theorem:rates}, then
\begin{align*}
\sigma^{-1}\sqrt{\underline C}(\widetilde \alpha - \alpha) \leadsto N(0,1),
\end{align*}
where $\sigma^2= Q^{-1}\Gamma Q^{-1}$, $Q={\mathrm{E}}[v_{11}^2]$ and
\begin{align*}
\Gamma=&\bar\mu_N \Gamma_N + \bar\mu_M \Gamma_M
=\bar\mu_N{\mathrm{E}}[v_{11} \varepsilon_{11} v_{12}\varepsilon_{12} ] + \bar\mu_M {\mathrm{E}}[ v_{11}\varepsilon_{11}v_{21} \varepsilon_{21}].
\end{align*}
\end{theorem}
A proof can be found in Appendix \ref{sec:theorem:asymptotic_normality}.
This result provides a theoretical justification for the asymptotic variance proposed in the overview in Section \ref{sec:overview}.
In practice, we do not know the components, $Q$ and $\Gamma$, of the asymptotic variance.
The following subsection proposes estimators of them.
\subsection{Variance Estimation}
In this section, we propose an analog variance estimator.
The components, $Q$ and $\Gamma$, of the asymptotic variance can be estimated by
\begin{align*}
\widehat Q=&\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M} \widehat v_{ij}^2
\qquad\text{and}
\\
\widehat \Gamma
=& \frac{\underline C}{(NM)^2}\sum_{i=1}^{N}\sum_{1\le j,j'\le M} \widehat v_{ij} \widehat \varepsilon_{ij} \widehat \varepsilon_{ij'} \widehat v_{ij'}
+ \frac{\underline C}{(NM)^2}\sum_{1\le i,i'\le N}\sum_{j=1}^{M} \widehat v_{ij}\widehat \varepsilon_{ij} \widehat\varepsilon_{i'j}\widehat v_{i'j} ,
\end{align*}
respectively, where $\widehat v_{ij}= D_{ij}- X_{ij}'\widehat \gamma$ and $\widehat \varepsilon_{ij}=Y_{ij}- \widehat \alpha D_{ij} - X_{ij}'\widehat\beta$ are the residuals.
With these component estimators, we propose that the asymptotic variance $\sigma^2 = Q^{-1} \Gamma Q^{-1}$ be estimated by $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$.
The following theorem provides a theoretical support for this variance estimator.
\begin{theorem}[Variance Estimation]\label{theorem:variance_est}
If Assumptions \ref{a:sampling}, \ref{a:moments}, \ref{a:sparsity} and \ref{a:sparse_eigenvalues} are satisfied, $\lambda_1$ and $\lambda_2$ are chosen according to the statement of Theorem \ref{theorem:rates},
$\frac{(NM)^{1/q} B_n^2 s^3 (\log a)^2}{\underline C ^2}=o(1)$,
$\frac{(NM)^{1/q} s \log a}{ \underline C}=o(1)$, and
$\|R^D_{ij}R^Y_{ij}\|^2_n=O(1)$,
then the variance estimator $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$ is consistent for $\sigma^2 = Q^{-1} \Gamma Q^{-1}$.
\end{theorem}
A proof is found in Appendix \ref{sec:theorem:variance_est}.
In light of this result, we propose to compute the standard error by $\widehat\sigma / \sqrt{\underline C}$.
Similarly, in light of this result together with Theorem \ref{theorem:asymptotic_normality}, we propose to construct the $\alpha^\ast$-level confidence interval by $\left[\widetilde\alpha + \Phi^{-1}(\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C},\widetilde\alpha + \Phi^{-1}(1-\alpha^\ast/2)\widehat\sigma / \sqrt{\underline C}\right]$, where $\Phi^{-1}$ denotes the quantile function of the standard normal distribution.
\section{Extension: Heterogeneous Cluster Sizes}\label{sec:heterogeneous_cluster_sizes}
Thus far, we focus on the case where each cluster contains one observation.
In this section, we presented an extension of the baseline results to situations where the numbers of observations are heterogeneous across clusters.
Suppose that we have $n_{ij}$ observations for each cell $(i,j) \in [N] \times [M]$, where $n_{ij}$ is a random variable that is allowed to depend on $(X_{ij,\ell})_{\ell\ge 1}$. To deal with the situation of $n_{ij}=0$, for any sequence $(t_\ell)_{\ell \ge 1}$, define $\sum_{\ell=1}^0 t_\ell=0$.
Consider the model
\begin{align*}
Y_{ij,\ell}=&D_{ij,\ell}\alpha +X_{ij,\ell}'\beta + R^Y_{ij,\ell} + \varepsilon_{ij,\ell},\:
{\mathrm{E}}[\varepsilon_{ij,\ell}|D_{ij,\ell},X_{ij,\ell}]=0,
\end{align*}
where
$Y_{ij,\ell}$ is an observed outcome variable,
$(D_{ij,\ell},X_{ij,\ell}')'$ is an observed vector of regressors, and
$R_{ij,\ell}^Y$ is an approximation error for the unit $\ell\in [n_{ij}]$ with the double index $(i,j)$ indicating $i$-th cluster in the first clustering dimension and $j$-th cluster in the second clustering dimension.
Using matrix notations, we can rewrite the model as
\begin{align*}
Y_{ij}=&D_{ij}\alpha +X_{ij}\beta + R^Y_{ij} + \varepsilon_{ij},\:
{\mathrm{E}}[\varepsilon_{ij}|D_{ij},X_{ij}]=0,
\end{align*}
where each of
$Y_{ij}=(Y_{ij,\ell})_{\ell\in [n_{ij}]}$,
$D_{ij}=(D_{ij,\ell})_{\ell\in [n_{ij}]}$,
$R^Y_{ij}=(R^Y_{ij,\ell})_{\ell\in [n_{ij}]}$,
and
$\varepsilon_{ij}=(\varepsilon_{ij,\ell})_{\ell\in [n_{ij}]}$ is of dimension $n_{ij}\times 1$,
and
$X_{ij}=(X_{ij,\ell}')_{\ell\in [n_{ij}]}$ is of dimension $n_{ij}\times p$.
We similarly write the accompanying auxiliary projection as
\begin{align*}
D_{ij}=&X_{ij}\gamma + R^D_{ij} + v_{ij},\:
{\mathrm{E}}[v_{ij}|X_{ij}]=0,
\end{align*}
where $R^D$ is of dimension $n_{ij}\times 1$ representing approximation errors, and $v_{ij}$ is of dimension $n_{ij}\times 1$ representing projection errors.
Under this setting, the first step of estimation procedure consists of
\begin{align}
(\widehat \alpha,\widehat \beta')'=& \underset{\alpha,\beta }{\rm argmin} \sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{\ell=1}^{n_{ij}} (Y_{ij,\ell} - \alpha D_{ij,\ell} - X_{ij,\ell}'\beta)^2 +\lambda_1\|(\alpha,\beta')'\|_1\nonumber\\
\text{and}\qquad
\widehat \gamma=& \underset{\gamma }{\rm argmin} \sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{\ell=1}^{n_{ij}} ( D_{ij,\ell} - X_{ij,\ell}'\gamma)^2 +\lambda_2 \|\gamma\|_1.\nonumber
\end{align}
In turn, the second-step estimates are obtained by
\begin{align}
(\widetilde \alpha,\widetilde \beta')= \underset{{\rm support}(\beta) \subset \widehat I}{\rm argmin} \sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{\ell=1}^{n_{ij}} (Y_{ij,\ell} - \alpha D_{ij,\ell} - X_{ij,\ell}'\beta)^2.\nonumber
\end{align}
The asymptotic variance estimator for $\widetilde\alpha$ is given by $\widehat\sigma^2= \widehat Q^{-1}\widehat\Gamma \widehat Q^{-1}$, where
\begin{align*}
\widehat Q=&\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M} \widehat v_{ij}'\widehat v_{ij},
\\
\widehat \Gamma
=& \frac{\underline C}{(NM)^2}\sum_{i=1}^{N}\sum_{1\le j,j'\le M} \widehat v_{ij}' \widehat \varepsilon_{ij}\widehat \varepsilon_{ij'}' \widehat v_{ij'}
+ \frac{\underline C}{(NM)^2}\sum_{1\le i,i'\le N}\sum_{j=1}^{M} \widehat v_{ij}'\widehat \varepsilon_{ij} \widehat\varepsilon_{i'j}' \widehat v_{i'j},
\end{align*}
$\widehat v_{ij}= D_{ij}- X_{ij}\widehat \gamma$, and $\widehat \varepsilon_{ij}=Y_{ij}- \widehat \alpha D_{ij} - X_{ij}\widehat\beta$.
We now formally state assumptions for the extended theory to support the asymptotic validity of this procedure.
Define
$W_{ij}=(n_{ij},(Z_{ij,\ell})_{\ell \ge 1})$ and $\ddot M=\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M} X_{ij}'X_{ij}$.
\begin{assumption}[Sampling]\label{a:sampling_hetero}
Suppose that $\underline C \to \infty $, $\mu_N \to \bar \mu_N \geqslant 0$, and $\mu_M \to \bar \mu_M \geqslant 0$.
\begin{enumerate}
\item $(W_{ij})_{(i,j)\in \mathbbm N^2}$ is an infinite sequence of separately exchangeable random processes.
\item $(W_{ij})_{(i,j)\in \mathbbm N^2}$ is dissociated.
\item For each $n$, an econometrician observes $((W_{ij,\ell})_{\ell\in [n_{ij}]} )_{i\in[N],j\in[M]}$.
\item ${\mathrm{E}}[n_{ij}]>0$ and $n_{ij}\le \bar n$ for a positive finite constant $\bar n$ independent of $n$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Moments]\label{a:moments_hetero}
There exists a sequence $\{B_n\}_{n=1}^\infty$ of positive constants such that the following conditions hold for all $n \in \mathbb{N}$ for some $q>4$:
\begin{enumerate}
\item
${\mathrm{E}}[\max_{\ell\in [n_{ij}]}|D_{11,\ell}|^{2q}]+\max_{k\in [p]}{\mathrm{E}}[\max_{\ell\in [n_{ij}]}|X_{11,\ell,k}|^{2q}]+ {\mathrm{E}}[\max_{\ell\in [n_{ij}]}|\varepsilon_{11,\ell}|^{2q}|X_{11,\ell},v_{11,\ell}]+{\mathrm{E}}[|v_{11,\ell}|^{2q}|X_{11,\ell}]\le K$ a.s.
and
$0 < c\le {\mathrm{E}}[\max_{\ell\in [n_{ij}]}v^{2}_{11,\ell}|X_{11,\ell}]$ a.s. for positive constants, $c$ and $K$, that are independent of $n$.
\item ${\mathrm{E}}[\max_{\ell\in [n_{ij}]}\|X_{11}\|_\infty^{2q}]\le B_{n}^{2q}$ and $B_n\sqrt{\log a}\lesssim (N\vee M)^{1/2-1/q}$.
\item $\bar \mu_N{\mathrm{E}}[v_{11}' \varepsilon_{11} \varepsilon_{12}'v_{12} ] + \bar \mu_M {\mathrm{E}}[ v_{11}'\varepsilon_{11}\varepsilon_{21}'v_{21} ] $ and the maximal and minimal eigenvalues of ${\mathrm{E}}[X_{11}'X_{11}]$ are bounded and bounded away from zero uniformly in $n$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Sparsity]\label{a:sparsity_hetero}
${}$
\begin{enumerate}
\item $\|\beta\|_0+ \|\gamma\|_0 \lesssim s$.
\item $\sqrt{(NM)^{-1}\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{\ell \in [n_{ij}]}(R_{ij,\ell}^Y)^2}+ \sqrt{(NM)^{-1}\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{\ell \in [n_{ij}]}(R_{ij,\ell}^D)^2}\le c_s \lesssim_{\mathrm{P}} \sqrt{s /\underline C}$.
\item $\frac{s^2_n(\log (a))^2}{\underline C}=o(1)$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Sparse Eigenvalues]\label{a:sparse_eigenvalues_hetero}
There exists a sequence $\{\ell_n\}$ such that $\ell_n \to \infty$ and, with probability at least $1-o(1)$,
\begin{align*}
0<c\le \phi_{{\rm min}}(s\ell_n)\le \phi_{{\rm max}}(s\ell_n)\le c' <\infty
\end{align*}
holds for some constants, $c$ and $c'$, that are independent of $n$, where
\begin{align*}
\phi_{{\rm max}}(m):=
\max_{\substack{1\le \|\delta\|_0 \le m }} \frac{\delta' \ddot M\delta}{\|\delta\|^2}
\text{ and } \phi_{{\rm min}}(m):= \min_{\substack{1\le \|\delta\|_0 \le m }} \frac{\delta'\ddot M\delta}{\|\delta\|^2}.
\end{align*}
\end{assumption}
The following statement provides a theoretical guarantee for the estimation and inference procedure for the extended model outlined above.
\begin{corollary}\label{corollary:hetero}
If Assumptions \ref{a:sampling_hetero}, \ref{a:moments_hetero}, \ref{a:sparsity_hetero} and \ref{a:sparse_eigenvalues_hetero} are satisfied, and $\lambda_1$ and $\lambda_2$ are chosen according to the statement of Theorem \ref{theorem:rates}, then
\begin{align*}
\sigma^{-1}\sqrt{\underline C}(\widetilde \alpha - \alpha) \leadsto N(0,1),
\end{align*}
where $\sigma^2= Q^{-1}\Gamma Q^{-1}$ with
\begin{align*}
Q=&{\mathrm{E}}[v_{11}'v_{11}],
\\
\Gamma
=&\bar\mu_N{\mathrm{E}}[v_{11}' \varepsilon_{11} v_{12}'\varepsilon_{12} ] + \bar\mu_M {\mathrm{E}}[ v_{11}'\varepsilon_{11}v_{21}' \varepsilon_{21}].
\end{align*}
Furthermore, if $\frac{(NM)^{1/q} B_n^2 s^3 (\log a)^2}{\underline C ^2}=o(1)$,
$\frac{(NM)^{1/q} s \log a}{ \underline C}=o(1)$, and
$\|R^D_{ij}R^Y_{ij}\|^2_n=O(1)$,
then the variance estimator $\widehat\sigma^2$ is consistent for $\sigma^2$.
\end{corollary}
A proof of Corollary \ref{corollary:hetero} closely follows that of the results in Section \ref{sec:asymptotic_theory}, and are therefore omitted.
The key difference is that we now apply Aldous-Hoover representation on $W_{ij}$ rather than on $Z_{ij}$.\footnote{For more insights on this extension, see Section 3.1 of \cite{DDG19}.}
\section{Simulation Studies}\label{sec:simulation_studies}
In this section, we present simulation studies of finite-sample performance of the proposed method of estimation and post-selection inference.
We compare the performance of our method against existing alternatives from the lasso literature that do not account for multi-way clustering.
\subsection{Simulation Setup}
We consider the linear model
\begin{align*}
Y_{ij} = D_{ij} \alpha + X_{ij}' \beta + \varepsilon_{ij}.
\end{align*}
The parameter values are fixed at $(\alpha,\beta')' = \left(0.5, 0.5^2, \cdots, 0.5^{\text{dim}(X)+1}\right)'$.
The random vector $(D_{ij}, X_{ij}', \varepsilon_{ij})$ is constructed by
\begin{align*}
\left(D_{ij}, X_{ij}\right)
=&
(1-\omega^{x}_1-\omega^{x}_2) \upsilon^{x}_{ij} + \omega^{x}_1 \upsilon^{x}_{i} + \omega^{x}_2 \upsilon^{x}_{j}
\qquad\text{and}\\
\varepsilon_{ij}
=&
(1-\omega^{\varepsilon}_1-\omega^{\varepsilon}_2) \upsilon^{\varepsilon}_{ij} + \omega^{\varepsilon}_1 \upsilon^{\varepsilon}_{i} + \omega^{\varepsilon}_2 \upsilon^{\varepsilon}_{j}
\end{align*}
with two-way clustering weights $(\omega^{x}_1, \omega^{x}_2)$ and $(\omega^{\varepsilon}_1,\omega^{\varepsilon}_2)$,
where $\upsilon^{x}_{ij}$, $\upsilon^{x}_{i}$, and $\upsilon^{x}_{j}$ are independently generated according to
\begin{align*}
\upsilon^{x}_{ij}, \upsilon^{x}_{i}, \upsilon^{x}_{j} \sim N\left(0,\left(\begin{array}{ccccc}\rho^0 & \rho^1 & \cdots & \rho^{\text{dim}(X)-1} & \rho^{\text{dim}(X)} \\ \rho^1 & \rho^0 & \cdots & \rho^{\text{dim}(X)-2} & \rho^{\text{dim}(X)-1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \rho^{\text{dim}(X)-1} & \rho^{\text{dim}(X)-2} & \cdots & \rho^0 & \rho^1 \\ \rho^{\text{dim}(X)} & \rho^{\text{dim}(X)-1} & \cdots & \rho^1 & \rho^0 \end{array}\right)\right),
\end{align*}
and $\upsilon^{\varepsilon}_{ij}$, $\upsilon^{\varepsilon}_{i}$, and $\upsilon^{\varepsilon}_{j}$ are independently generated according to
\begin{align*}
\upsilon^{\varepsilon}_{ij}, \upsilon^{\varepsilon}_{i}, \upsilon^{\varepsilon}_{j} \sim N(0,1).
\end{align*}
Note that the weights $(\omega^{x}_1, \omega^{x}_2)$ and $(\omega^{\varepsilon}_1,\omega^{\varepsilon}_2)$ specify the extent of dependence in two-way clustering in $(D_{ij},X_{ij}')$ and $\varepsilon_{ij}$, respectively.
Also, the parameter $\rho$ specifies the extent of collinearity among the high-dimensional covariates $(D_{ij},X_{ij}')$.
We set
$(\omega^{x}_1, \omega^{x}_2) = (0.25,0.25)$,
$(\omega^{\varepsilon}_1,\omega^{\varepsilon}_2) = (0.25,0.25)$, and
$\rho = 0.50$.
\subsection{Alternative Variance Estimators}
We compare the performance of our multi-way cluster-robust variance estimator with two existing alternative benchmarks.
One is the heteroskedasticity robust variance estimator (such as the one in \cite{BCH14}) without accounting for cluster sampling, i.e., $\Gamma$ is estimated by
\begin{align*}
\widehat \Gamma_{HC}=\frac{1}{NM} \sum_{i=1}^{N}\sum_{j=1}^{M} \widehat v_{ij}^2 \widehat \varepsilon_{ij}^2.
\end{align*}
We will refer to this variance estimator $\widehat Q^{-1} \widehat \Gamma_{HC} \widehat Q^{-1}$ as the `0-Way' estimator.
The other is the one-way cluster-robust variance estimator (similar to those of \cite{BCHK16} and \cite{kock16}) clustered at one (e.g., second) dimension, i.e., $\Gamma$ is estimated by
\begin{align*}
\widehat \Gamma_{CR}=\frac{1}{NM^2} \sum_{i=1}^{N}\sum_{1\le j,j'\le M} \widehat v_{ij}\widehat \varepsilon_{ij}\widehat v_{ij'}\widehat \varepsilon_{ij'}.
\end{align*}
We will refer to this variance estimator $\widehat Q^{-1} \widehat \Gamma_{CR} \widehat Q^{-1}$ as the `1-Way' estimator.
\subsection{Results}
Table \ref{tab:simulation_results} summarizes simulation results.
The first two columns indicate the two-way sample sizes $(N,M)$.
The third column indicates the dimension (Dim) of $(\alpha,\beta')'$.
The next four columns report simulation statistics for $\widetilde\alpha$.
These statistics include the average (Avg), bias (Bias), standard deviation (SD), and root mean square error (RMSE).
The last three columns report 95\% coverage frequencies of $\alpha$ based on three variance estimators.
The first is the heteroskedasticity robust variance estimator (0-Way).
The second is the one-way cluster-robust variance estimator (1-Way).
The third is our multi-way cluster-robust variance estimator (2-Way).
The results are based on 25,000 Monte Carlo iterations for each row in the table.
In view of the statistics columns, observe that the post-double-selection lasso estimate $\widetilde\alpha$ behaves well in larger sample sizes (e.g., $N, M \geqslant 20$) both in terms of bias and variance.
Next, observe the 95\% coverage frequencies by the three alternative variance estimators.
Both the 0-Way and 1-Way variance estimators significantly underestimate the variances of the post-double-selection lasso estimate $\widetilde\alpha$.
On the other hand, the coverage frequency based on our 2-Way variance estimator approaches the nominal probability (95\%) as the sample size increases.
These results demonstrate that, when the true sampling process entails multi-way clustering, traditional variance estimators may bias the inference and our multi-way cluster-robust variance estimator performs robustly well.
\section{Empirical Illustrations}\label{sec:empirical illustration}
In this section, we illustrate our proposed method with applications to a couple of empirical studies.
There is a sequence of recent growth and development economic studies using empirical data that are clustered at ethnic and geographical levels \citep[e.g.,][]{NunnWantchekon11,Michalopoulos_Papaioannou2013,MP14,MP16,gershman16,anderson18,dickens18}.
The next two subsections present how our method can enrich the model flexibility and robustness of such studies, focusing on the cases of \citet{NunnWantchekon11} and \citet{Michalopoulos_Papaioannou2013}.
\subsection{Slave Trade and Mistrust in Africa}\label{sec:slave trade and mistrust in africa}
\citet{NunnWantchekon11} analyze the effects of slave trade on mistrust in Africa, controlling for various demongraphic and geographical covariates including age, age squared, ethnic fractionalization, gender, urban residence, occupation, religion, and living conditions as well as country fixed effects in their baseline model.
Estimates of these effects are obtained by running regressions with a sample that pools $n_{ij}$ individuals $\ell \in [n_{ij}]$ in ethnic group $i$ and districts $j$ across the cells $(i,j) \in [N] \times [M]$ of $N (=185)$ ethnic groups and $M (=1257)$ districts.
Standard errors are computed by the two-way cluster-robust method of \citet{CGM11} for the ethnic group and district as two ways of clustering.
With our proposed method that is applicable to both high-dimensional models and multi-way clustering, they could consider even more flexible model specifications, for example, allowing for higher orders of age rather than just the quadratic specification and interactions of the age polynomials with various other dummy variables.
We present estimates with standard errors under such extended models with flexible specifications, demonstrate that qualitatively similar results continue to be obtained without substantial loss of statistical significance, and thus confirm further robustness of the main empirical findings by \citet{NunnWantchekon11}.
Consider the model
\begin{align*}
Y_{ij,\ell}=&D_{ij,\ell}\alpha +X_{ij,\ell}'\beta + R^Y_{ij,\ell} + \varepsilon_{ij,\ell},\:
{\mathrm{E}}[\varepsilon_{ij,\ell}|D_{ij,\ell},X_{ij,\ell}]=0,
\end{align*}
where $Y_{ij,\ell}$ denotes a measure of trust, $D_{ij,\ell}$ denotes an intensity measure of slave trade, $X_{ij,\ell}$ contains polynomial basis elements of age up to degree 10, ethnic fractionalization, gender, urban residence, occupation, religion, living conditions, the interactions of the polynomial basis of age with all the dummy variables, and country fixed effects, consisting of 597 dimensions of covariates in total.
Note that the total number of regressors ($p+1=598$) is much larger than the effective sample size ($\underline C = N \wedge M = 185$) of two-way clustering in this extended setting.
Table \ref{tab:results_nunn} summarize the estimates of the effects of slave trade on mistrust as measured by the ``trust of neighbors,'' corresponding to Table 1 of \citet{NunnWantchekon11}.
The last two columns in the table show the original estimates obtained under the prototypical model by \citet[][Table 1]{NunnWantchekon11} and corresponding lasso estimates obtained under more flexible model specification by our method.
Across all the measures of slave exports, the original estimates and our lasso estimates are similar with similar levels of statistical significance.
These results demonstrate that, even for flexible model specifications entailing high-dimensional covariates,
the proposed method allows to produce qualitatively similar results without extensive loss of significance, and we can thus confirm further robustness of the main empirical findings by \citet{NunnWantchekon11}.
\subsection{Pre-Colonial Institutions and Regional Developments in Africa}\label{sec:pre colonial institutions and regional developments in africa}
\citet{Michalopoulos_Papaioannou2013} analyze the effects of pre-colonial institutions on contemporary regional developments in Africa, controlling for various population, locational and geographic covariates including population density, distance to capital, distance to sea coast, distance to border, water area, land area, elevation, land suitable for agriculture, ecological suitability, petrolium, and diamond mine as well as country fixed effects in their baseline model.
Estimates of these effects are obtained by running regressions with a sample that pools $n_{ij}$ populated pixels $\ell \in [n_{ij}]$ in ethnic group $i$ and country $j$ across the cells $(i,j) \in [N] \times [M]$ of $N (=93)$ ethnic groups and $M (=48)$ countries.
Standard errors are computed by the two-way cluster-robust method of \citet{CGM11} for the ethnic group and district as two ways of clustering.
With our proposed method that is applicable to both high-dimensional models and multi-way clustering, they could consider even more flexible model specifications, for example, allowing for interactions of all combinations of geographical covariates and locational covariates.
We present estimates with standard errors under such extended models with flexible specifications, demonstrate that qualitatively similar results continue to be obtained without substantial loss of statistical significance, and thus confirm further robustness of the main empirical findings by \citet{Michalopoulos_Papaioannou2013}.
Consider the model
\begin{align*}
Y_{ij,\ell}=&D_{ij,\ell}\alpha +X_{ij,\ell}'\beta + R^Y_{ij,\ell} + \varepsilon_{ij,\ell},\:
{\mathrm{E}}[\varepsilon_{ij,\ell}|D_{ij,\ell},X_{ij,\ell}]=0,
\end{align*}
where $Y_{ij,\ell}$ denotes a regional development measured by night light intensity, $D_{ij,\ell}$ denotes an intensity measure of pre-colonial ethinic institutions, $X_{ij,\ell}$ contains population density, interactions of all combinations of locational controls (distance to capital, distance to sea coast, and distance to border), interactions of all combinations of geographical controls (water area, land area, elevation, land suitable for agriculture, ecological suitability, petrolium, and diamond mine), and country fixed effects, consisting of 82 or 83 dimensions of covariates in total.
Note that the total number of regressors ($p+1=83$ or 84) is much larger than the effective sample size ($\underline C = N \wedge M = 48$) of two-way clustering in this extended setting.
Table \ref{tab:results_michalopoulos} summarize the estimates of the effects of pre-colonial institutions on regional development as measured by the ``light density,'' corresponding to parts of Table 3 of \citet{Michalopoulos_Papaioannou2013}.
The last two columns in the table show the original estimates obtained under the prototypical model by \citet[][Table 3]{Michalopoulos_Papaioannou2013} and corresponding lasso estimates obtained under more flexible model specification by our method.
Across all the measures of pre-colonial institutions and all specifications, the original estimates and our lasso estimates are similar with similar levels of statistical significance.
These results demonstrate that, even for flexible model specifications entailing high-dimensional covariates,
the proposed method allows to produce qualitatively similar results without extensive loss of significance, and we can thus confirm further robustness of the main empirical findings by \citet{Michalopoulos_Papaioannou2013}.
\section{Conclusion}\label{sec:conclusion}
In this paper, we investigate high-dimensional regression models when data is sampled under multi-way clustering.
We establish the convergence rates for the lasso and post-lasso estimators under multi-way clustering.
We then propose an inference method based on a post-double-selection procedure and show that it is asymptotically valid under multi-way clustering.
Simulation studies demonstrate that the proposed procedure works well in finite sample under multi-way clustering in comparison with existing alternatives.
We demonstrate that our method can enrich the flexibility of regression models and robustness of empirical results through a couple of empirical applications in growth and development economics.
Indeed, both multi-way clustering and high dimensionality are two important issues which concern applied research.
The existing literature provide solutions to each of multi-way clustering and high-dimensionality separately.
To our best knowledge, the literature does not seem to provide a solution to both of these issues simultaneously.
In this paper, we filled this void in the literature.
\newpage