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,728 characters
Multiway Cluster Robust Double/Debiased Machine Learning
\title{Multiway Cluster Robust Double/Debiased Machine Learning\thanks{First arXiv date: September 8, 2019\smallskip}}
\author{
Harold D. Chiang\thanks{Harold D. Chiang: [email removed]. Department of Economics, Vanderbilt University, VU Station B \#351819, 2301 Vanderbilt Place, Nashville, TN 37235-1819, USA\smallskip}
\qquad Kengo Kato\thanks{Kengo Kato: [email removed]. Department of Statistics and Data Science, Cornell University, 1194 Comstock Hall, Ithaca, NY 14853, USA\smallskip}
\qquad Yukun Ma\thanks{Yukun Ma: [email removed]. Department of Economics, Vanderbilt University, VU Station B \#351819, 2301 Vanderbilt Place, Nashville, TN 37235-1819, USA\smallskip} \qquad Yuya Sasaki\thanks{Yuya Sasaki: [email removed]. Department of Economics, Vanderbilt University, VU Station B \#351819, 2301 Vanderbilt Place, Nashville, TN 37235-1819, USA\smallskip} \thanks{We benefited from useful comments by seminar participants at Southern Methodist University, Stony Brook University, University of Bristol, and University of Colorado - Boulder, and participants at CeMMAP UCL/Vanderbilt Joint Conference on Advances in Econometrics and CeMMAP Workshop on Causal Learning with Interactions. All remaining errors are ours.\bigskip}
}
\date{}
\maketitle
\begin{abstract}
This paper investigates double/debiased machine learning (DML) under multiway clustered sampling environments.
We propose a novel multiway cross fitting algorithm and a multiway DML estimator based on this algorithm.
We also develop a multiway cluster robust standard error formula.
Simulations indicate that the proposed procedure has favorable finite sample performance.
Applying the proposed method to market share data for demand analysis, we obtain larger two-way cluster robust standard errors for the price coefficient than non-robust ones in the demand model.
\bigskip\\
{\bf Keywords:} double/debiased machine learning, multiway clustering, multiway cross fitting
\bigskip\\
{\bf JEL Codes:} C10, C13, C14
\bigskip\\${}$\bigskip\\${}$
\end{abstract}
\section{Introduction}
We propose a novel multiway cross fitting algorithm and a double/debiased machine learning (DML) estimator based on the proposed algorithm.
This objective is motivated by recently growing interest in use of dependent cross sectional data and recently increasing demand for DML methods in empirical research.
On one hand, researchers frequently use multiway cluster sampled data in empirical studies, such as network data, matched employer-employee data, matched student-teacher data, scanner data where observations are double-indexed by stores and products, and market share data where observations are double-indexed by market and products.
On the other hand, we have witnessed rapidly increasing popularity of machine learning methods in empirical studies, such as random forests, lasso, post-lasso, elastic nets, ridge, deep neural networks, and boosted trees among others.
To date, available DML methods focus on i.i.d. sampled data.
In light of the aforementioned research environments today, a new method of DML that is applicable to multiway cluster sampled data may well be of interest by empirical researchers.
The DML was proposed by the recent influential paper by \citet[CCDDHNR,][]{CCDDHNR18}.
They provide a general DML toolbox for estimation and inference for structural parameters with high-dimensional and/or infinite-dimensional nuisance parameters.
In that paper, the estimation method and properties of the estimator are presented under the typical microeconometric assumption of i.i.d. sampling.
We advance this frontier literature of DML by proposing a modified DML estimation procedure with multiway cross fitting, which accommodates multiway cluster sampled data.
Even for multiway cluster sampled data, we show that the proposed DML procedure works under nearly identical set of assumptions to that of CCDDHNR (\citeyear{CCDDHNR18}).
To our best knowledge, the present paper is the first to consider generic DML methods under multiway cluster sampling.
Another branch of the literature following the seminal work by \cite{CGM11} proposes multiway cluster robust inference methods.
\cite{Menzel17} conducts formal analyses of bootstrap validity under multiway cluster sampling robustly accounting for non-degenerate and degenerate cases.
\cite{DDG18} develop empirical process theory under multiway cluster sampling which applies to a large class of models.
We advance this practically important literature by developing a multiway cluster robust inference method based on DML.
In deriving theoretical properties of the proposed estimator, we take advantage of the Aldous-Hoover representation employed by the preceding papers.
To our knowledge, the present paper is the first in this literature on multiway clustering to develop generic DML methods.
\subsection{Relations to the Literature}
The past few years have seen a fast growing literature in machine learning based econometric methods. For general overviews of the field, see, e.g., \cite{AtheyImbens19} or \cite{MullainathanSpiess17}. For a review of estimation and inference methods for high-dimensional data, see \cite{BCH14review}. For an overview of data sketching methods tackling computationally impractically large number of observations, see \cite{LeeNg19}.
The DML of CCDDHNR (\citeyear{CCDDHNR18}) is built upon \cite{BCK15}, which proposes to use Neyman orthogonal moments for a general class of Z-estimation statistical problems in the presence of high-dimensional nuisance parameters.
This framework is further generalized in different directions by \cite{BCFH17} and \cite{BCCW18}.
CCDDHNR (\citeyear{CCDDHNR18}) combine the use of Neyman orthogonality condition with cross fitting to provide a simple yet widely applicable framework that covers a large class of models under i.i.d. settings.
The DML is also compatible with various types of machine learning based methods for nuisance parameter estimation.
Driven by the need from empiricists, the literature on cluster robust inference has a long history in econometrics.
For recent review of the literature, see, e.g., \cite{CM15} and \cite{MacKinnon2019}.
On the other hand, coping with cross-sectional dependence using a multiway cluster robust variance estimator is a relatively recent phenomenon.
\cite{CGM11} first provide a multiway cluster robust variance estimator for linear regression models without imposing additional parametric assumptions on the intra-cluster correlation structure.
This variance estimator has significantly reshaped the landscape of econometric practices in applied microeconomics in the past decade.\footnote{As of December 31, 2019, \cite{CGM11} has received over 2,500 citations. The majority of such citations came from applied economic papers.}
In contrast to the popularity among empirical researchers, theoretical justification of the validity of this type of procedures was lagging behind.
The first rigorous treatment of asymptotic properties of multiway cluster robust estimators are established by \cite{Menzel17} using the Aldous-Hoover representation under the assumptions of separable exchangeability and dissociation.
The asymptotic theory of \cite{Menzel17} covers both non-degenerate and degenerate cases.
Focusing on non-degenerate situations, \cite{DDG18} further extend this approach to a general empirical process theory.\footnote{See also \cite{DDG19} for further generalization of the empirical process theory for dyadic data under joint exchangeability assumption.}
Using this asymptotic framework, \cite{MacKinnonNielsenWebb2019} study linear regression models under the non-degenerate case and examine the validity of several types of wild bootstrap procedures and the robustness of multiway cluster robust variance estimators under different cluster sampling settings.
Despite of the popularity of both machine learning and cluster robust inference among empirical researchers, relatively limited cluster robust inference results exist for machine learning based methods.
Inference for machine learning based methods with one-way clustering is studied by \cite{BCHK16}, \cite{Kock2016}, \cite{KockTang2018}, \cite{SGCT18} and \cite{HansenLiao19}
for different variations of regularized regression estimators and \cite{AtheyWager19} for random forests.
\cite{ChiangSasaki2019} investigate the performance of lasso and post-lasso in the partially linear model setting of \cite{BCH14} under multiway cluster sampling.
To our best knowledge, there is no general machine learning based procedures with known validity under multiway cluster sampling environments.
\section{Overview}\label{sec:overview}
\subsection{Setup}\label{sec:setup}
Suppose that the researcher observes a sample $\left\{\left. W_{ij} \right\vert i \in \{1,...,N\}, j \in \{1,...,M\}\right\}$ of double-indexed observations of size $NM$.
Let $P$ denote the probability law of $\{W_{ij}\}_{ij}$, and let ${\rm E}_{P}$ denote the expectation with respect to $P$.
Let $\underline C= N \wedge M$ denote the sample size in the smaller dimension.
We consider two-way clustering where each cell contains one observation for simplicity of notations, but results for higher cluster dimensions and random cluster sizes can be obtained at the expense of involved notations -- see Appendix \ref{sec:extension_to_general_multiway_clustering} for a general case.
The structural model is assumed to entail the moment restriction
\begin{align}
{\rm E}_{P}[\psi(W_{11};\theta_0,\eta_0)]=0
\label{eq:existance_condition}
\end{align}
for some score $\psi$ that depends on a low-dimensional parameter vector $\theta \in \Theta \subset \mathbbm R^{d_\theta}$ and a nuisance parameter $\eta \in T$ for a convex subset $T$ of a normed linear space.
The nuisance parameter $\eta$ may be finite-, high-, or infinite-dimensional, and its true value is denoted by $\eta_0 \in T$.
In this setup, the true value of the low-dimensional target parameter, denoted by $\theta_0 \in \Theta$, is the object of interest.
Let $\widetilde T=\{\eta - \eta_0 : \eta \in T\}$, and define the Gateaux derivative map $D_r: \widetilde T \rightarrow \mathbbm R^{d_\theta}$ by
\begin{align*}
D_r[\eta-\eta_0]:=\partial_r \Big\{
{\rm E}_{P}[\psi(W_{11};\theta_0,\eta_0+r(\eta-\eta_0))]\Big\}
\end{align*}
for all $r\in[0,1)$.
Also denote its limit by
\begin{align*}
\partial_\eta{\rm E}_{P}\psi(W_{11};\theta_0,\eta_0)[\eta - \eta_0]:=D_0[\eta-\eta_0].
\end{align*}
We say that the Neyman orthogonality condition holds at $(\theta_0,\eta_0)$ with respect to a nuisance realization set $\mathcal T_n \subset T$ if the score $\psi$ satisfies (\ref{eq:existance_condition}), the pathwise derivative $D_r[\eta-\eta_0]$ exists for all $r\in[0,1)$ and $\eta\in \mathcal T_n$, and the orthogonality equation
\begin{align}
\partial_\eta{\rm E}_{P}\psi(W_{11};\theta_0,\eta_0)[\eta - \eta_0]=0
\label{eq:Neyman_orthogonal_condition}
\end{align}
holds for all $\eta\in \mathcal T_n$.
Furthermore, we also say that the $\lambda_n$ Neyman near-orthogonality condition holds at $(\theta_0,\eta_0)$ with respect to a nuisance realization set $\mathcal T_n\subset T$ if the score $\psi$ satisfies (\ref{eq:existance_condition}), the pathwise derivative $D_r[\eta-\eta_0]$ exists for all $r\in[0,1)$ and $\eta\in \mathcal T_n$, and the orthogonality equation
\begin{align}
\sup_{\eta \in \mathcal T_n}\Big\| \partial_\eta {\rm E}_{P}\psi(W;\theta_0,\eta_0)[\eta-\eta_0] \Big\|\le \lambda_n
\label{eq:Neyman_near_orthogonal_condition}
\end{align}
holds for all $\eta\in \mathcal T_n$
for some positive sequence $\{\lambda_n\}_n$ such that $\lambda_n=o(\underline C^{-1/2})$.
Throughout, we will consider structural models satisfying the moment restriction (\ref{eq:existance_condition}) and either form of the Neyman orthogonality conditions, (\ref{eq:Neyman_orthogonal_condition}) or (\ref{eq:Neyman_near_orthogonal_condition}).
Consider linear Neyman orthogonal scores $\psi$ of the form
\begin{align}
\psi(w;\theta,\eta)=\psi^a(w;\eta)\theta +\psi^b(w;\eta), \text{ for all $w\in \rm{supp}(W)$, $\theta\in\Theta$, $\eta\in T$. } \label{eq:linear_score}
\end{align}
A generalization to nonlinear score follows from linearization with Gateaux differentiability as in Section 3.3 of CCDDHNR (\citeyear{CCDDHNR18}).
We focus on linear scores as they cover a wide range of applications.
\subsection{The Multiway Double/Debiased Machine Learning}\label{sec:multiway_dml}
For the class of models introduced in Section \ref{sec:setup}, we propose a novel $K^2$-fold multiway cross fitting procedure for estimation of $\theta_0$.
For any $r \in \mathbb N$, we use the notation $[r]=\{1,...,r\}$.
With a fixed positive integer $K$, randomly partition $[N]$ into $K$ parts $\{I_1,...,I_K\}$ and $[M]$ into $K$ parts $\{J_1,...,J_K\}$.
For each $(k,\ell) \in [K]^2$, obtain an estimate
$$\widehat \eta_{k\ell}=\widehat \eta\left((W_{ij})_{(i,j)\in ([N]\setminus I_k )\times ([M]\setminus J_\ell)}\right)$$
of the nuisance parameter $\eta$ by some machine learning method (e.g., lasso, post-lasso, elastic nets, ridge, deep neural networks, and boosted trees) using only the subsample of those observations with multiway indices $(i,j)$ in $([N]\setminus I_k ) \times ([M]\setminus J_\ell)$.
In turn, we define $\widetilde \theta$, the multiway double/debiased machine learning (multiway DML) estimator for $\theta_0$, as the solution to
\begin{align}
\frac{1}{K^2}\sum_{(k,\ell)\in [K]^2} \mathbbm E_{n,k\ell}[\psi(W;\widetilde \theta,\widehat \eta_{k\ell})] =0,\label{eq:MDML}
\end{align}
where $\mathbbm E_{n,k\ell} [f(W)] = \frac{1}{|I_k||J_\ell|}\sum_{(i,j)\in I_k\times J_\ell} f(W_{ij})$ denotes the subsample empirical expectation using only the those observations with multiway indices $(i,j)$ in $I_k \times J_\ell$.
We call this procedure the $K^2$-fold multiway cross fitting.
Note that, for each $(k,\ell)\in [K]^2$, the nuisance parameter estimate $\widehat\eta_{k\ell}$ is computed using the subsample of those observations with multiway indices $(i,j) \in ([N]\setminus I_k ) \times ([M]\setminus J_\ell)$, and in turn the score term $\mathbbm E_{n,k\ell}[\psi(W; \cdot,\widehat\eta_{k\ell})]$ is computed using the subsample of those observations with multiway indices $(i,j) \in I_k \times J_\ell$.
This two-step computation is repeated $K^2$ times for every partitioning pair $(k,\ell)\in [K]^2$.
Figure \ref{fig:cross_fitting} illustrates this $K^2$-fold cross fitting for the case of $K=2$ and $N=M=4$, where the cross fitting repeats for $K^2 (= 2^2 = 4)$ times.
\begin{figure}
\caption{An illustration of $2^2$-fold cross fitting.}\label{fig:cross_fitting}
\tikzstyle{my help lines}=[gray,
thick,dashed]
\begin{multicols}{4}
\qquad\\
\begin{tikzpicture}
\draw(4,0) grid (3,1)node[above,black]{Nuisance};
\draw (3,1) grid (2,2);
\draw (3,0) grid (2,1);
\draw (3,2) grid (4,1);
\draw (0,4) grid (1,3)node[above,black]{Score};
\draw (1,4) grid (2,3);
\draw (0,3) grid (1,2);
\draw (1,3) grid (2,2);
\draw[style=my help lines] (2,0) grid (0,2);
\draw[style=my help lines] (4,2) grid (2,4);
\end{tikzpicture}
\qquad\\
\begin{tikzpicture}
\draw(4,2) grid (3,3)node[above,black]{Nuisance};
\draw(3,3) grid (2,4);
\draw(4,3) grid (3,4);
\draw(3,2) grid (2,3);
\draw(2,0) grid (1,1)node[above,black]{Score};
\draw(1,1) grid (0,2);
\draw(2,1) grid(0,2);
\draw(1,0) grid(0,1);
\draw[style=my help lines] (4,0) grid (2,2);
\draw[style=my help lines] (2,2) grid (0,4);
\end{tikzpicture}
\qquad\\
\begin{tikzpicture}
\draw(4,2) [darkgray] grid (3,3)node[above,black]{Score};
\draw(3,3) [darkgray] grid (2,4);
\draw(4,3) [darkgray] grid (3,4);
\draw(3,2) [darkgray] grid (2,3);
\draw(2,0) [darkgray] grid (1,1)node[above,black]{Nuisance};
\draw(1,1) [darkgray] grid (0,2);
\draw(2,1) [darkgray] grid(0,2);
\draw(1,0) [darkgray] grid(0,1);
\draw[style=my help lines] (4,0) grid (2,2);
\draw[style=my help lines] (2,2) grid (0,4);
\end{tikzpicture}
\qquad\\
\begin{tikzpicture}
\draw (4,0) [darkgray] grid (3,1)node[above,black]{Score};
\draw (3,1) [darkgray] grid (2,2);
\draw (3,0) [darkgray] grid (2,1);
\draw (3,2) [darkgray] grid (4,1);
\draw (0,4) [darkgray] grid (1,3)node[above,black]{Nuisance};
\draw (1,4) [darkgray] grid (2,3);
\draw (0,3) [darkgray] grid (1,2);
\draw (1,3) [darkgray] grid (2,2);
\draw[style=my help lines] (2,0) grid (0,2);
\draw[style=my help lines] (4,2) grid (2,4);
\end{tikzpicture}
\end{multicols}
\end{figure}
\begin{remark}
This estimator is a multiway-counterpart of DML2 in CCDDHNR (\citeyear{CCDDHNR18}).
It is also possible to consider the multiway-counterpart of their DML1.
With this said, we focus on this current estimator following their simulation finding that DML2 outperforms their DML1 in most situation settings due to the stability of the score function.
\end{remark}
\begin{remark}[Higher Cluster Dimensions]\label{remark:higher_cluster_dim}
When we have $\alpha$-way clustering for an integer $\alpha>2$, the above algorithm can be easily generalized into a $K^\alpha$-fold multiway DML estimator. See Appendix \ref{sec:extension_to_general_multiway_clustering} for a generalization.
\end{remark}
We propose to estimate the asymptotic variance of $\sqrt{\underline C}(\widetilde\theta-\theta_0)$ by
\begin{align}
\widehat\sigma^2
=&
\widehat J^{-1}
\widehat\Gamma
(\widehat J^{-1})',
\label{eq:variance_estimator}
\end{align}
where $\widehat \Gamma$ and $\widehat J$ are given by
\begin{align}
\widehat \Gamma
=&
\frac{1}{K^2}\sum_{(k,\ell)\in [K]^2} \left\{
\frac{|I|\wedge|J|}{(|I||J|)^2}\sum_{i\in I_k}\sum_{j,j'\in J_\ell} \psi(W_{ij};\widetilde \theta,\widehat\eta_{k\ell}) \psi(W_{ij'};\widetilde \theta,\widehat\eta_{k\ell})'
\right.
\nonumber\\
& \qquad\qquad \ \
+
\left.
\frac{|I|\wedge|J|}{(|I||J|)^2}\sum_{i,i'\in I_k}\sum_{j\in J_\ell} \psi(W_{ij};\widetilde \theta,\widehat\eta_{k\ell}) \psi(W_{i'j};\widetilde \theta,\widehat\eta_{k\ell}) '
\right\}
\qquad\text{and}
\nonumber\\
\widehat J
=&
\frac{1}{K^2}\sum_{(k,\ell)\in [K]^2}\mathbbm E_{n,k\ell}[\psi^a(W;\widehat \eta_{k\ell})],\nonumber
\end{align}
accounting for multiway cluster dependence.
For a $d_\theta$-dimensional vector $r$, the $(1-a)$ confidence interval for the linear functional $r'\theta_0$ can be constructed by
\begin{align*}
\text{CI}_a:=[r'\widetilde \theta\pm \Phi^{-1}(1-a/2)\sqrt{r'\widehat \sigma^2 r/\underline C}].
\end{align*}
\subsection{Example: Partially Linear IV Model with Multiway Cluster Sample}\label{sec:example_partially_linear}
For an illustration, consider as a concrete example the partially linear IV model (cf. Okui, Small, Tan and Robins, \citeyear{OkuiSmallTanRobins2012} ; CCDDHNR, \citeyear{CCDDHNR18}, Section 4.2) adapted to the multiway cluster sample data:
\begin{align}
Y_{ij}=&D_{ij}\theta_0 + g_0(X_{ij})+\epsilon_{ij},\qquad{\rm E}_{P}[\epsilon_{ij}|X_{ij},Z_{ij}]=0,
\label{eq:example:reduced_form}
\\
Z_{ij}=&m_0(X_{ij})+v_{ij},\quad\:\qquad\qquad{\rm E}_{P}[v_{ij}|X_{ij}]=0.
\label{eq:example:projection}
\end{align}
A researcher observes the random variables $Y_{ij}$, $D_{ij}$, $X_{ij}$, and $Z_{ij}$, which are typically interpreted as the outcome, endogenous regressor, exogenous regressors, and instrumental variable, respectively.
The low-dimensional parameter vector $\theta_0$ is an object of interest.
A Neyman orthogonal score $\psi$ for such model is given by
\begin{align}
\psi(w;\theta,\eta)=(y-g_1(x)-\theta(d - g_2(x)))(z-m(x))
\label{eq:IV_Neymand_orthogonal_moment}
\end{align}
as in \cite{OkuiSmallTanRobins2012} and CCDDHNR (\citeyear{CCDDHNR18}),
where $w=(y,d,x,z)$, $\eta=(g_1,g_2,m)$ and $g_1$, $g_2$, $m\in L^2(P)$.
It is straightforward to verify that this score satisfies both the moment restriction (\ref{eq:existance_condition}), ${\rm E}_{P}[\psi(W_{11};\theta_0,\eta_0)]=0$, and the Neyman orthogonality condition (\ref{eq:Neyman_orthogonal_condition}), $\partial_\eta {\rm E}_{P} \psi(W_{11};\theta_0,\eta_0)[\eta - \eta_0]=0$ for all $\eta \in \mathcal{T}_n$ at $\eta_0=(g_{10},g_{20},m_0)$, where $g_{10}(X)={\rm E}_{P}[Y|X]$, $g_{20}(X)={\rm E}_{P}[D|X]$, and $m_0(X)={\rm E}_{P}[Z|X]$.
The following algorithm is our proposed multiway DML procedure introduced in Section \ref{sec:multiway_dml}, specifically applied to this partially linear IV model.
\begin{algorithm}[$K^2$-fold Multiway DML for Partially Linear IV Model with Lasso]\label{algorithm:partial_linear_iv}
${}$
\begin{enumerate}
\item Randomly partition $[N]$ into $K$ parts $\{I_1,...,I_K\}$ and $[M]$ into $K$ parts $\{J_1,...,J_K\}$.
\item For each $(k,\ell)\in[K]^2$:
\begin{enumerate}
\item Run a lasso of $Y$ on $X$ to obtain $\widehat g_{1,k\ell}(x)=x'\widehat\beta_{k\ell}$ using observations from $I_k^c\times J_\ell^c$.
\item Run a lasso of $D$ on $X$ to obtain $\widehat g_{2,k\ell}(x)=x'\widehat\gamma_{k\ell}$ using observations from $I_k^c\times J_\ell^c$.
\item Run a lasso of $Z$ on $X$ to obtain $\widehat m_{k\ell}(x)=x'\widehat \xi_{k\ell}$ using observations from $I_k^c\times J_\ell^c$.
\end{enumerate}
\item Solve the equation
\begin{align*}
\frac{1}{K^2}\sum_{(k,\ell) \in [K]^2}\mathbbm E_{n,k\ell}[(Y_{ij}-X_{ij}'\widehat\beta_{k\ell}-\theta(D_{ij} - X_{ij}'\widehat \gamma_{k\ell}))(Z_{ij}-X_{ij}'\widehat\xi_{k\ell})]=0
\end{align*}
for $\theta$ to obtain the multiway DML estimate $\widetilde\theta$.
\item Let $\widehat \varepsilon_{ij}=Y_{ij}- X_{ij}' \widehat\beta_{k \ell} - \widetilde\theta (D_{ij}-X_{ij}'\widehat\gamma_{k\ell})$, $\widehat u_{ij} = D_{ij} - X_{ij}'\widehat\gamma_{k\ell}$, and $\widehat v_{ij}=Z_{ij}-X_{ij}'\widehat \xi_{k \ell}$ for each $(i,j) \in I_k \times J_\ell$ for each $(k,\ell) \in [K]^2$, and let the multiway DML asymptotic variance estimator be given by
\begin{align*}
\widehat \sigma^2=&\widehat J^{-1} \frac{1}{K^2}\sum_{k=1}^K \sum_{\ell=1}^K\Big\{
\frac{|I|\wedge |J|}{(|I||J|)^2}\sum_{i\in I_k}\sum_{j,j'\in J_\ell}
\widehat \varepsilon_{ij}\widehat v_{ij} \widehat v_{ij'} \widehat \varepsilon_{ij'}
+
\frac{|I|\wedge |J|}{(|I||J|)^2}\sum_{i,i'\in I_k}\sum_{j\in J_\ell} \widehat \varepsilon_{ij}\widehat v_{ij} \widehat v_{i'j} \widehat \varepsilon_{i'j}
\Big\}(\widehat J^{-1})',
\end{align*}
where
\begin{align*}
\widehat J=&-\frac{1}{K^2}\sum_{k=1}^K\sum_{\ell=1}^K\mathbbm E_{n,k\ell}[\widehat u_{ij}\widehat v_{ij} ].
\end{align*}
\item Report the estimate $\widetilde\theta$, its standard error $\sqrt{\widehat\sigma^2/\underline C}$, and/or the $(1-a)$ confidence interval
\begin{align*}
\text{CI}_a:=\left[\widetilde \theta\pm \Phi^{-1}(1-a/2)\sqrt{\widehat \sigma^2 /\underline C}\right].
\end{align*}
\end{enumerate}
\end{algorithm}
For the sake of concreteness, we present this algorithm specifically based on lasso (in the three sub-steps under step 2), but another machine learning method (e.g., post-lasso, elastic nets, ridge, deep neural networks, and boosted trees) may be substituted for lasso.
\begin{example}[Demand Analysis]\label{ex:demand_analysis}
Consider the model of \citet{Berry94}
in which consumer $c$ derives the utility
\begin{align*}
\delta_{ij} + X_{ij}\alpha_c + \varepsilon_{cij}
\end{align*}
from choosing product $i$ in market $j$,
where $\varepsilon_{cij}$ independently follows the Type I Extreme Value distribution, $\alpha_c$ is a random coefficient, and the mean utility $\delta_{ij}$ takes the linear-index form
\begin{align*}
\delta_{ij} = D_{ij}\theta_0 + \epsilon_{ij}.
\end{align*}
In this framework, \citet[][Equation (9)]{LuShiTao19} derive the partial-linear equation
\begin{align*}
Y_{ij} = D_{ij}\theta_0 + g_0(X_{ij}) + \epsilon_{ij}
\end{align*}
for estimation of $\theta_0$, where $Y_{ij} = \log( S_{ij} ) - \log( S_{0j} )$ denotes the observed log share of product $i$ relative to the log of the outside share.
Since $D_{ij}$ usually consists of the endogenous price of product $i$ in market $j$, researchers often use instruments $Z_{ij}$ such that ${\rm E}_{P}[\epsilon_{ij}|X_{ij},Z_{ij}]=0$.
This yields the reduced-form equation (\ref{eq:example:reduced_form}), together with the innocuous nonparametric projection equation (\ref{eq:example:projection}).
Since the random vector $W_{ij} = (Y_{ij},D_{ij},X_{ij},Z_{ij})$ is double-indexed by product $i$ and market $j$, the sample naturally entails two-way dependence.
Specifically, for each product $i$, $\{W_{ij}\}_{j=1}^M$ is likely dependent through a supply shock by the producer of product $i$.
Similarly, for each market $j$, $\{W_{ij}\}_{i=1}^N$ is likely dependent through a demand shock in market $j$.
As such, instead of using standard errors based on i.i.d. sampling, we recommend that a researcher uses the two-way cluster-robust standard error based on Algorithm \ref{algorithm:partial_linear_iv}.
$\triangle$
\end{example}
\section{Theory of the Multiway DML}\label{sec:Theory of the Multiway DML}
In this section, we present formal theories to guarantee that the multiway DML method proposed in Section \ref{sec:overview} works.
We first fix some notations for convenience.
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$.
With this said, we will suppress the index notation and write $(N,M)$ for simplicity.
Let $\{\mathcal P_n\}_n$ be a sequence of sets of probability laws of $\{W_{ij}\}_{ij}$ -- note that we allow for increasing dimensionality of $W_{ij}$ in the sample size $n$.
Let $P=P_{n}\in \mathcal P_n$ denote the law with respect to sample size $(N,M)$.
Throughout, we assume that this random vector $W_{ij}$ is Borel measurable.
Recall the notations $\underline C =N\wedge M$, $\mu_N=\underline C/N$, and $\mu_M=\underline C/M$, and suppose that $\mu_N\to \bar \mu_N$, $\mu_M\to \bar \mu_M$.
We write $a \lesssim b$ to mean $a \leq cb$ for some $c > 0$ that does not depend on $n$.
We also write $a \lesssim_P b$ to mean $a = O_P(b)$.
For any finite dimensional vector $v$, $\|v\|$ denotes the $\ell_2$ or Euclidean norm of $v$.
For any matrix $A$, $\|A\|$ denotes the induced $\ell_2$-norm of the matrix. For any set $B$, $|B|$ denotes the cardinality of the set.
We state the following assumption on multiway clustered sampling.
\begin{assumption}[Sampling]\label{a:sampling}
Suppose $\underline C \to \infty $.
The following conditions hold for each $n$.
\begin{enumerate}[(i)]
\item $(W_{ij})_{(i,j)\in \mathbbm N^2}$ is an infinite sequence of separately exchangeable $p$-dimensional random vectors.
That is, for any permutations $\pi_1$ and $\pi_2$ of $\mathbbm N$, we have
\begin{align*}
(W_{ij})_{(i,j)\in \mathbbm N^2}\overset{d}{=} (W_{\pi_1(i)\pi_2(j)})_{(i,j)\in \mathbbm N^2}.
\end{align*}
\item $(W_{ij})_{(i,j)\in \mathbbm N^2}$ is dissociated.
That is, for any $(c_1,c_2)\in \mathbbm N^2$,
$
(W_{ij})_{i \in [c_1], j \in [c_2]}
$
is independent of
$
(W_{ij})_{i \in [c_1]^c, j \in [c_2]^c}.
$
\item For each $n$, an econometrician observes $(W_{ij})_{i\in[N],j\in[M]}$.
\end{enumerate}
\end{assumption}
Recall that we focus on the linear Neyman orthogonal score of the form
\begin{align*}
\psi(w;\theta,\eta)=\psi^a(w;\eta)\theta +\psi^b(w;\eta), \text{ for all $w\in \rm{supp}(W)$, $\theta\in\Theta$, $\eta\in T$. }
\end{align*}
Let $c_0>0$, $c_1>0$, $s>0$, $q\ge 4$
be some finite constants with $c_0\le c_1$.
Let $\{\delta_n\}_{n\ge 1}$ (estimation errors) and $\{\Delta_n\}_{n\ge 1}$ (probability bounds) be sequences of positive constants that converge to zero such that $\delta_n \ge \underline C^{-1/2}$.
Let $K\ge 2$ be a fixed integer.
Let $W_{00}$ denote a copy of $W_{11}$ that is independent from the data and the random set $\mathcal T_n$ of nuisance realization.
With these notations, we consider the following assumptions.
\begin{assumption}[Linear Neyman Orthogonal Score]\label{a:linear_orthogonal_score}
For $\underline C\ge 3$ and $P\in \mathcal P_n$, the following conditions hold.
\begin{enumerate}[(i)]
\item The true parameter value $\theta_0$ satisfies (\ref{eq:existance_condition}).
\item $\psi$ is linear in the sense that it satisfies (\ref{eq:linear_score}).
\item The map $\eta \mapsto {\rm E}_{P}[\psi(W_{00};\theta,\eta)]$ is twice continuously Gateaux differentiable on $T$.
\item $\psi$ satisfies either the Neyman orthogonality condition (\ref{eq:Neyman_orthogonal_condition}) or more generally
the Neyman $\lambda_n$ near orthogonality condition at $(\theta_0,\eta_0)$ with respect to a nuisance realization set $\mathcal T_n\subset T$ as
\begin{align*}
\lambda_n:=\sup_{\eta \in \mathcal T_n}\Big\| \partial_\eta {\rm E}_{P}\psi(W_{00};\theta_0,\eta_0)[\eta-\eta_0] \Big\|\le \delta_n \underline C^{-1/2}.
\end{align*}
\item The identification condition holds as the singular values of the matrix $J_0:={\rm E}_{P}[\psi^a(W_{11};\eta_0)]$ are between $c_0$ and $c_1$.
\end{enumerate}
\end{assumption}
\begin{assumption}[Score Regularity and Nuisance Parameter Estimators]\label{a:regularity_nuisance_parameters}
For all $\underline C\ge 3$ and $P\in \mathcal P_n$, the following conditions hold.
\begin{enumerate}[(i)]
\item Given random subsets $I\subset [N]$ and $J\subset [M]$ such that $|I|\times |J|=\lfloor NM/K^2\rfloor$, the nuisance parameter estimator $\widehat \eta=\widehat\eta((W_{ij})_{(i,j)\in I^c\times J^c}) $, where the complements are taken with respect to $[N]$ and $[M]$, respectively, belongs to the realization set $\mathcal T_n$ with probability at least $1-\Delta_n$, where $\mathcal T_n$ contains $\eta_0 $.
\item The following moment conditions hold:
\begin{align*}
m_n:=& \sup_{\eta\in \mathcal T_n}({\rm E}_{P}[\|\psi(W_{00};\theta_0,\eta)\|^q])^{1/q} \le c_1,\\
m_n':=& \sup_{\eta\in \mathcal T_n}({\rm E}_{P}[\|\psi^a(W_{00};\eta)\|^q])^{1/q} \le c_1.
\end{align*}
\item The following conditions on the rates $r_n$, $r_n'$ and $\lambda_n'$ hold:
\begin{align*}
r_n:=& \sup_{\eta\in \mathcal T_n}
\|{\rm E}_{P}[\psi^a(W_{00};\eta)]-{\rm E}_{P}[\psi^a(W_{00};\eta_0)]\|\le \delta_n,\\
r_n':=& \sup_{\eta\in \mathcal T_n}
(\|{\rm E}_{P}[\psi(W_{00};\theta_0,\eta)]-{\rm E}_{P}[\psi(W_{00};\theta_0,\eta_0)]\|^2)^{1/2}\le \delta_n,\\
\lambda_n'= & \sup_{r\in (0,1),\eta\in \mathcal T_n}\|\partial^2_r {\rm E}_{P}[\psi (W_{00};\theta_0,\eta_0+r(\eta-\eta_0)) ] \|\le \delta_n/\sqrt{\underline C}.
\end{align*}
\item All eigenvalues of the matrix
\begin{align*}
\Gamma:=\bar\mu_N \Gamma_N + \bar\mu_M \Gamma_M=\bar\mu_N{\rm E}_{P} [\psi(W_{11};\theta_0,\eta_0)\psi(W_{12};\theta_0,\eta_0)']
+ \bar\mu_M{\rm E}_{P}[\psi(W_{11};\theta_0,\eta_0)\psi(W_{21};\theta_0,\eta_0)'].
\end{align*}
are bounded from below by $c_0$.
\end{enumerate}
\end{assumption}
\begin{remark}[Discussion of the Assumptions]
Assumption \ref{a:sampling} is similar to those of the preceding work on multiway cluster robust inference \citep[cf.][]{Menzel17,DDG18,ChiangSasaki2019}.
\citet{Menzel17} does not invoke the dissociation, and follows an alternative approach to inference.
The other papers assume both the separate exchangeability and dissociation, and conduct unconditional inference as in this paper.
See \citet[][Corollary 7.23 and Lemma 7.35]{Kallenberg2006} for representations with and without the dissociation under the separate exchangeability.
Assumption \ref{a:linear_orthogonal_score}
is closely related to Assumptions 3.1 of CCDDHNR (\citeyear{CCDDHNR18}). It requires the score to be
Neyman near orthogonal -- see their Section 2.2.1 for the procedure of orthogonalizing a non-orthogonal score.
It also imposes some mild smoothness and identification conditions.
Assumption \ref{a:regularity_nuisance_parameters} corresponds to Assumption 3.2 of CCDDHNR (\citeyear{CCDDHNR18}). It imposes some high level conditions on the quality of the nuisance parameter estimator as well as the non-degeneracy of the asymptotic variance. This rules out the degenerate cases such as Example 1.6 of \cite{Menzel17}.
\end{remark}
\begin{remark}[Partial Distributions]
Assumptions \ref{a:linear_orthogonal_score} and \ref{a:regularity_nuisance_parameters} state conditions based on $W_{00}$, differently from CCDDHNR (\citeyear{CCDDHNR18}), because of our need to deal with dependent observations in cross fitting in our multiway DML framework.
\end{remark}
The following result presents the main theorem of this paper, establishing the linear representation and asymptotic normality of the multiway DML estimator.
It corresponds to Theorem 3.1 of CCDDHNR (\citeyear{CCDDHNR18}), and is an extension of it to the case of multiway cluster sampling.
\begin{theorem}[Main Result]\label{theorem:main_result_linear}
Suppose that Assumptions \ref{a:sampling}, \ref{a:linear_orthogonal_score} and \ref{a:regularity_nuisance_parameters} are satisfied.
If $\delta_n\ge \underline C^{-1/2}$ for all $\underline C\ge 1$, then
\begin{align*}
\sqrt{\underline C}\sigma^{-1}(\widetilde \theta - \theta_0)=\frac{\sqrt{\underline C}}{NM}\sum_{i=1}^N \sum_{j=1}^M
\bar \psi(W_{ij})+O_P(\rho_n)\leadsto N(0,I_{d_\theta})
\end{align*}
holds uniformly over $P\in\mathcal P_n$, where the size of the remainder terms follows
\begin{align*}
\rho_n :=\underline C^{-1/2} + r_n +r_n' + \underline C^{1/2} \lambda_n + \underline C^{1/2} \lambda_n'\lesssim \delta_n,
\end{align*}
the influence function takes the form $\bar \psi(\cdot):=-\sigma^{-1}J_0^{-1} \psi(\cdot;\theta_0,\eta_0)$,
and the asymptotic variance is given by
\begin{align}
\sigma^2:=J_0^{-1}\Gamma (J_0^{-1})'. \label{eq:population_variance}
\end{align}
\end{theorem}
As is commonly the case in practice, we need to estimate the unknown asymptotic variance.
The following theorem shows the validity of our proposed multiway DML variance estimator.
\begin{theorem}[Variance Estimator]\label{theorem:variance_estimator_linear}
Under the assumptions required by Theorem \ref{theorem:main_result_linear}, we have
\begin{align*}
\widehat \sigma^2=\sigma^2 +O_P(\rho_n).
\end{align*}
Furthermore, the statement of Theorem \ref{theorem:main_result_linear} holds true with $\widehat \sigma^2$ in place of $\sigma^2$.
\end{theorem}
Theorems \ref{theorem:main_result_linear} and \ref{theorem:variance_estimator_linear} can be used for constructing confidence intervals.
\begin{corollary}\label{corollary:inference_t-test}
Suppose that all the Assumptions required by Theorem \ref{theorem:main_result_linear} are satisfied.
Let $r$ be a $d_\theta$-dimensional vector.
The $(1-a)$ confidence interval of $r'\theta_0$ given by
\begin{align*}
\text{CI}_a:=[r'\widetilde \theta\pm \Phi^{-1}(1-a/2)\sqrt{r'\widehat \sigma^2 r/\underline C}]
\end{align*}
satisfies
\begin{align*}
\sup_{P\in\mathcal P_n}|P_P(\theta_0 \in \text{CI}_a)-(1-a)|\to 0.
\end{align*}
\end{corollary}
As in Section 3.4 of CCDDHNR (\citeyear{CCDDHNR18}), we can also repeatedly compute multiway DML estimates and variance estimates $S$-times for some fixed $S\in \mathbbm N$ and consider the average or median of the estimates as the new estimate.
This does not have an asymptotic impact, yet it can reduce the impact of a random sample splitting on the estimate.
\section{Simulation Studies}\label{sec:simulation_studies}
\subsection{Simulation Setup}
Consider the partially linear IV model introduced in Section \ref{sec:example_partially_linear}.
We specifically focus on the following high-dimensional linear representations
\begin{align*}
Y_{ij} =& D_{ij}\theta_0 + X_{ij}'\zeta_0 + \epsilon_{ij}
\\
D_{ij} =& Z_{ij}\pi_{10} + X_{ij}'\pi_{20} + \upsilon_{ij},
\\
Z_{ij} =& X_{ij}'\xi_0 + V_{ij},
\end{align*}
where the parameter values are set to $\theta_0 = \pi_{10} = 1.0$ and $\zeta_0 = \pi_{20} = \xi_0 = (0.5,.0.5^2,\cdots,0.5^{\text{dim}(X)})'$ for some large $\text{dim}(X)$.
The primitive random vector $(X_{ij}',\epsilon_{ij},\upsilon_{ij},V_{ij})'$ is constructed by
\begin{align*}
X_{ij} &= (1-\omega_1^X - \omega_2^X) \alpha_{ij}^X + \omega_1^X \alpha_i^X + \omega_2^X \alpha_j^X,
\\
\epsilon_{ij} &= (1-\omega_1^\epsilon - \omega_2^\epsilon) \alpha_{ij}^\epsilon + \omega_1^\epsilon \alpha_i^\epsilon + \omega_2^\epsilon \alpha_j^\epsilon,
\\
\upsilon_{ij} &= (1-\omega_1^\upsilon - \omega_2^\upsilon) \alpha_{ij}^\upsilon + \omega_1^\upsilon \alpha_i^\upsilon + \omega_2^\upsilon \alpha_j^\upsilon,
\qquad\text{and}\\
V_{ij} &= (1-\omega_1^V - \omega_2^V) \alpha_{ij}^V + \omega_1^V \alpha_i^V + \omega_2^V \alpha_j^V
\end{align*}
with two-way clustering weights $(\omega_1^X,\omega_2^X)$, $(\omega_1^\epsilon,\omega_2^\epsilon)$, $(\omega_1^\upsilon,\omega_2^\upsilon)$, and $(\omega_1^V,\omega_2^V)$, where
$\alpha_{ij}^X$, $\alpha_{i}^X$, and $\alpha_{j}^X$ are independently generated according to
\begin{align*}
\alpha_{ij}^X, \alpha_{i}^X, \alpha_{j}^X
\sim N
\left(0, \left(\begin{array}{ccccc}
s_X^0 & s_X^1 & \cdots & s_X^{\text{dim}(X)-2} & s_X^{\text{dim}(X)-1} \\
s_X^1 & s_X^0 & \cdots & s_X^{\text{dim}(X)-3} & s_X^{\text{dim}(X)-2} \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
s_X^{\text{dim}(X)-2} & s_X^{\text{dim}(X)-3} & \cdots & s_X^0 & s_X^1 \\
s_X^{\text{dim}(X)-1} & s_X^{\text{dim}(X)-2} & \cdots & s_X^1 & s_X^0
\end{array}\right)\right),
\end{align*}
$(\alpha_{ij}^\epsilon,\alpha_{ij}^\upsilon)'$, $(\alpha_{i}^\epsilon,\alpha_{i}^\upsilon)'$, and $(\alpha_{j}^\epsilon,\alpha_{j}^\upsilon)'$ are independently generated according to
\begin{align*}
\left(\begin{array}{c}\alpha_{ij}^\epsilon \\ \alpha_{ij}^\upsilon\end{array}\right),
\left(\begin{array}{c}\alpha_{i}^\epsilon \\ \alpha_{i}^\upsilon\end{array}\right),
\left(\begin{array}{c}\alpha_{j}^\epsilon \\ \alpha_{j}^\upsilon\end{array}\right)
\sim N
\left(0, \left(\begin{array}{cc}
1 & s_{\epsilon\upsilon} \\
s_{\epsilon\upsilon} & 1
\end{array}\right)\right),
\end{align*}
and $\alpha_{ij}^V$, $\alpha_{i}^V$, and $\alpha_{j}^V$ are independently generated according to
\begin{align*}
\alpha_{ij}^V, \alpha_{i}^V, \alpha_{j}^V
\sim
N(0,1).
\end{align*}
The weights $(\omega_1^X,\omega_2^X)$, $(\omega_1^\epsilon,\omega_2^\epsilon)$, $(\omega_1^\upsilon,\omega_2^\upsilon)$, and $(\omega_1^V,\omega_2^V)$ specify the extent of dependence in two-way clustering in $X_{ij}$, $\epsilon_{ij}$, $\upsilon_{ij}$, and $V_{ij}$, respepctively.
The parameter $s_X$ specifies the extent of collinearity among the high-dimensional regressors $X_{ij}$.
The parameter $s_{\epsilon\upsilon}$ specifies the extent of endogeneity.
We set the values of these parameters to $(\omega_1^X,\omega_2^X) = (\omega_1^\epsilon,\omega_2^\epsilon) = (\omega_1^\upsilon,\omega_2^\upsilon) = (\omega_1^V,\omega_2^V) = (0.25, 0.25)$ and $s_X = s_{\epsilon\upsilon} = 0.25$.
\subsection{Results}
Monte Carlo simulations are conducted with 2,500 iterations for each set.
Table \ref{tab:simulation_results} reports simulation results.
The first four columns in the table indicate the data generating process ($N$, $M$, $\underline C$, and dim$(X)$).
The next column indicates the integer $K$ for our $K^2$-fold cross fitting method.
We use $K=2$ and $3$ in the simulations for the displayed results, since $2^2 (\approx 5)$ and $3^2 (\approx 10)$ are close to the common numbers of folds used in cross fitting in practice.
The next column indicates the machine learning method for estimation of $\widehat\eta_{k\ell}$.
We use the ridge, elastic net, and lasso.
The last four columns of the table report Monte Carlo simulation statistics, including the bias (Bias), standard deviation (SD), root mean square error (RMSE), and coverage frequency for the nominal probability of 95\% (Cover).
For each covariate dimension $\text{dim}(X) \in \{100,200\}$, for each choice $K \in \{2,3\}$ for the number $K^2$ of multiway cross fitting, and for each of the three machine learning methods, we observe the following patterns as the effective sample size $\underline C=N \wedge M$ increases: 1) the bias tends to zero; 2) the standard deviation decreases approximately at the $\sqrt{\underline C}$ rate; and 3) the coverage frequency converges to the nominal probability.
These results confirm the theoretical properties of the proposed method.
We ran several other sets of simulations besides those displayed in the table, and this pattern remains the same across different sets.
Comparing the results across the three machine learning methods, we observe that the ridge entails larger bias and smaller variance relative to the elastic net and lasso in finite sample.
This makes the coverage frequency of the ridge less accurate compared with the elastic net and lasso.
This result is perhaps specific to the data generating process used for our simulations.
On one hand, the choice $K=3$ (i.e., $9$-fold) of the multiway cross fitting contributes to mitigating the large bias of the ridge relative to the choice $K=2$, and hence $K=3$ produces more preferred results for the ridge.
On the other hand, the choice $K=2$ tends to yield preferred results in terms of coverage accuracy for the elastic net and lasso.
In light of these results, we recommend the elastic net or lasso along with the use of $2^2$- fold (i.e., $4$-fold) cross fitting.
This number of folds in cross fitting is in fact similar to that recommended by CCDDHNR (\citeyear{CCDDHNR18}) for i.i.d. sampling -- see their Remark 3.1 where they recommend 4- or 5-fold cross fitting.
\section{Empirical Illustration: Demand Analysis with Market Share Data}\label{sec:empirical_illustration}
Let us revisit the demand model of Example \ref{ex:demand_analysis} in Section \ref{sec:example_partially_linear}.
Recall that, for the consumer demand model of \citet{Berry94} introduced in Example \ref{ex:demand_analysis}, \citet[][Equation (9)]{LuShiTao19} derive the partial-linear equation
\begin{align}\label{eq:partial_lienar_demand}
Y_{ij} = D_{ij}\theta_0 + g_0(X_{ij}) + \epsilon_{ij}
\end{align}
for estimation of $\theta_0$, where $Y_{ij} = \log( S_{ij} ) - \log( S_{0j} )$ denotes the observed log share of product $i$ relative to the log of the outside share in market $j$, $D_{ij}$ denotes the log price of product $i$ in market $j$, and $X_{ij}$ denotes a vector of observed attributes of product $i$ in market $j$.
To deal with the likely endogeneity of $D_{ij}$, researchers often use instruments $Z_{ij}$ such that ${\rm E}_{P}[\epsilon_{ij}|X_{ij},Z_{ij}]=0$.
Such instruments often consist of observed attributes of other products in the market.
The implied equation (\ref{eq:partial_lienar_demand}) together with this mean independence assumption yields the reduced-form model (\ref{eq:example:reduced_form}).
Furthermore, we write the innocuous nonparametric projection equation (\ref{eq:example:projection}).
Therefore, we apply Algorithm \ref{algorithm:partial_linear_iv} in Section \ref{sec:example_partially_linear} for the two-way cluster robust DML estimation of $\theta_0$ with a robust standard error.
We present an application of the proposed algorithm to the U.S. automobile data of \citet{BLP95}.
The sample consists of unbalanced two-way clustered observations with $N=557$ models of automobiles and $M=20$ markets.
The observed attributes $X_{ij}$ consist of horsepower per weight, miles per dollar, miles per gallon, and size.
The instrument $Z_{ij}$ is defined as the sum of the values of these attributes of other products.
For the purpose of highlighting the effect of clustering assumptions, we report estimates and standard errors under the zero-way cluster robust DML (based on the i.i.d. assumption) and the one-way cluster robust DML (based on clustering along each of the product and market dimensions), as well as the two-way cluster robust DML (along both of the product and market dimensions).
The number $K=4$ of folds of cross fitting is used for the zero- and one-way cluster robust DML, while the number $K^2=4$ of folds of two-way cross fitting is used for the two-way cluster robust DML following the recommendations from Section \ref{sec:simulation_studies} and those by CCDDHNR (\citeyear{CCDDHNR18}, Remark 3.1).
To mitigate the uncertainty induced by sample splitting, we compute estimates based on the average of ten rerandomized DML following CCDDHNR (\citeyear{CCDDHNR18}, Section 3.4) with variance estimation according to CCDDHNR (\citeyear{CCDDHNR18}, Equation 3.13) adapted to our two-way cluster-robustness.
Table \ref{tab:empirical_results} summarizes the results.
For each of the zero-, one-, and two-way cluster robust DML, both the point estimates and standard errors are similar across all the choices of instrument.
Furthermore, the point estimates are also similar across all of the zero-, one-, and two-way cluster robust DML.
On the other hand, the standard errors tend to increase as the assumed number of ways of clustering increases.
In other words, the zero-way cluster robust DML reports the smallest standard error while the two-way cluster robust DML reports the largest standard error.
To robustly account for possible cross-sectional dependence of observations in such two-way cluster sampled data as this market share data, we recommend that researchers use the two-way cluster robust DML although it may incur larger standard errors as is the case with this application.
\section{Conclusion}\label{sec:conclusion}
In this paper, we propose a multiway DML procedure based on a new multiway cross fitting algorithm. This multiway DML procedure is valid in the presence of multiway cluster sampled data, which is frequently used in empirical research.
We present an asymptotic theory showing that multiway DML is valid under nearly identical reguarity conditions to those of CCDDHNR (\citeyear{CCDDHNR18}).
The proposed method covers a large class of econometric models as is the case with CCDDHNR (\citeyear{CCDDHNR18}), and is compatible with various machine learning based estimation methods.
Simulation studies indicate that the proposed procedure has attractive finite sample performance under various multiway cluster sampling environments for various machine learning methods.
To accompany the theoretical findings, we provide easy-to-implement algorithms for multiway DML.
Such algorithms are readily implementable using existing statistical packages.
There are a couple of possible directions for future research.
First, whereas we focused on linear orthogonal scores that cover a wide range of applications, it may be possible to develop a method and theories for non-linear orthogonal scores as in CCDDHNR (\citeyear{CCDDHNR18}; Section 3.3).
Second, whereas we focused on unconditional moment restrictions, it may be possible and will be important to develop a method and theories for conditional moment restrictions \citep{AiChen2003,AiChen2007,ChenLintonKeilegom2003,ChenPouzo2015}.
We leave these and other extensions for future research.
\newpage