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.
60,925 characters
Inference in High-Dimensional Regression Models without the Exact or $L^p$ sparsity
\title{Inference in High-Dimensional Regression Models \\without the Exact or $L^p$ sparsity}
\thanks{We thank two anonymous referees for very useful comments and suggestions. All the remaining errors are ours. H.D. Chiang is supported by the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation.}
\author[J. Cha]{Jooyoung Cha}
\author[H. D. Chiang]{Harold D. Chiang}
\author[Y. Sasaki]{Yuya Sasaki}
\date{First arXiv version: This version: December 30, 2022}
\address[J. Cha]{
Department of Economics, Vanderbilt University\\
VU Station B \#351819, 2301 Vanderbilt Place,
Nashville, TN 37235-1819, USA.}
\email{[email removed]}
\address[H. D. Chiang]{
Department of Economics, University of Wisconsin-Madison\\
William H. Sewell Social Science Building, 1180 Observatory Drive,
Madison, WI 53706, USA.}
\email{[email removed]}
\address[Y. Sasaki]{
Department of Economics, Vanderbilt University\\
VU Station B \#351819, 2301 Vanderbilt Place,
Nashville, TN 37235-1819, USA.}
\email{[email removed]}
\begin{abstract}
We propose a new inference method in high-dimensional regression models and high-dimensional IV regression models.
The method is shown to be valid without requiring the exact sparsity or $L^p$ sparsity conditions.
Simulation studies demonstrate superior performance of this proposed method over those based on the LASSO or the random forest, especially under less sparse models.
We illustrate an application to production analysis with a panel of Chilean firms.
Our results are closer to the benchmark conventional estimates than the estimates by other machine learning methods.
\bigskip
\noindent {\bf Keywords:}
double/debiased machine learning, high-dimensional Akaike information criterion, orthogonal greedy algorithm, production function.
\medskip
\end{abstract}
\maketitle
\section{Introduction}
The advent of modern machine learning\footnote{\label{foot:ml}We use the phrase ``machine learning'' following the recent related literature, but it is worthy to remark that it is a synonym of semiparametric or high-dimensional `estimation.'} techniques has significantly widened the class of analyzable regression models that include models with high-dimensional controls and/or models with flexible nonlinearity.
Perhaps the most popular and important machine learning approaches to estimation and inference in high-dimensional regression models today are those based on shrinkage and regularization, such as the least absolute shrinkage and selection operation (LASSO).
While they have practically appealing properties, these popular machine learning methods yet rely on a list of assumptions that may not be necessarily mild under certain applications.
In particular, the assumptions of the exact sparsity and the $L^p$ sparsity required for these methods are sometimes controversial and perceived to be strong for some applications.
As such, there still remains room in the literature for further widening the class of analyzable high-dimensional regression models if these assumption can be relaxed by a new machine learning method.
This paper proposes a method of inference in high-dimensional regression models without requiring the exact sparsity or the $L^p$ sparsity.
We set a low-dimensional parameter vector as the object of interest, and treat the remaining high-dimensional parameter vector as a nuisance component.
Under this common setting, our proposed method of inference works as follows.
First, use the orthogonal greedy algorithm \citep[OGA;][]{temlyakov2000weak} to order the high-dimensional regressors in a descending order of explanatory power.
Second, use the high-dimensional Akaike information criterion \citep[HDAIC;][]{ing2020} to select a model among the ordered list of models constructed in the first step.
Third, estimate the models selected in the second step.
Fourth, plug the estimated selected models in a Neyman orthogonal score and estimate the low-dimensional parameter vector of interest.
We take advantage of a number of recent methodological and theoretical developments, namely OGA, HDAIC, and DML, to derive asymptotic statistical properties of this new method of inference.
\citet{ing2020} investigates convergence rate properties of an estimator of high-dimensional regression models based on OGA and HDAIC (hereafter referred to as OGA+HDAIC).
Importantly, the setting imposes assumptions on the high-dimensional parameter vector that are weaker than the exact sparsity and the $L^p$ sparsity.
Furthermore, this approach does not require to impose high-level conditions on the sample Gram matrix, such as a restricted eigenvalue type condition, that are often required in the literature.
Even under these weaker conditions, it is still possible to obtain similar rates of convergence to those based on existing machine learning techniques such as the LASSO.
Given the adequately fast convergence rates of the preliminary nuisance parameter estimators based on OGA+HDAIC, we can apply the post-double-selection approach \citep*{BCH2014RES} or the double/debiased machine learning (DML) framework \citep*{ccddhnr} to in turn obtain a root-$N$ convergence of the low-dimensional parameter vector of interest with a limit normal distribution.
Simulation studies demonstrate that this proposed method performs significantly better than those based on the LASSO or the random forest, especially when the data generating model becomes less sparse.
We apply the proposed method to an analysis of production functions using a panel of Chilean firms.
Our results are closer to the benchmark conventional estimates \citep{levinsohn2003estimating} than the estimates based on other machine learning methods, namely the LASSO and the random forest.
{\bf Relation to the literature:}
This paper is related to several branches of the econometrics and statistics literature.
First, it is closely related to the literature on inference in high-dimensional regression models and high-dimensional IV regression models, e.g., \cite{belloni2012sparse}, \cite{BCH2014RES}, \cite{javanmard2014confidence}, \cite{van2014asymptotically}, \cite{zhang2014confidence}, \cite{caner2018asymptotically}, \cite{caner2018high}, \cite{belloni2018high}, \cite{galbraith2020simple}, \cite{gold2020inference}, \cite{kueck2021estimation} to list a few.
The majority of these papers focus on utilizing LASSO \citep{tibshirani1996regression} and its variants for estimation, and thus rely crucially on the exact sparsity, approximate sparsity\footnote{The approximate sparsity is closely related to the exact sparsity.} or $L^p$ sparsity.
We contribute to this literature, as emphasized above, by relaxing the conventional assumptions of these notions of sparsity.
Second, also related is the literature on model selection methods for high-dimensional models.
For extensive reviews of this vast literature, we refer readers to the monographs of \cite{buhlmann2011statistics}, \cite{giraud2015introduction}, and \cite{hastie2019statistical}.
In particular, we take advantage of the theoretical results of OGA+HDAIC by \cite{ing2020} as one of the main auxiliary steps to our goal as emphasized earlier. Methodologically, our paper benefits from the theoretical studies of various greedy algorithms in \cite{temlyakov2000weak}, \cite{tropp2004greed}, \cite{tropp2007signal}, \cite{ing2011stepwise}, and share ties with other iterated model selection methods such as the least absolute angle regression of \cite{efron2004least}, the $L_2$-boosting of \cite{buhlmann2003boosting}, the test-based forward model selection of \citet{kozbur2017testing,kozbur2020analysis}, and so forth.
Third, this paper is related to the literature on Neyman orthogonal scores or locally robust scores, e.g., \cite{belloni2015uniform}, \cite{chernozhukov2016locally}, \cite{BCCW2018}, and, in particular, the DML \citep*{ccddhnr}.
We contribute to this literature by proposing to add OGA+HDAIC to the library of the list of preliminary estimators.
Fourth, the conditions that we impose in place of the exact sparsity and the $L^p$ sparsity concern the speed at which the absolute size of the regression parameters decays when descendingly ordered. These conditions are analogous to those that are used for model selection problems in autoregressive time series models \citep[e.g.,][]{shibata1980asymptotically,ing2007accumulated} as well as the ordinary- and super-smoothness on probability density functions that are used in the deconvolution literature \citep[e.g.,][]{fan1991optimal,fan1993nonparametric}.
\section{High-Dimensional Linear Regression Models}\label{sec:regression}
\subsection{The Model}
Consider the linear regression model
\begin{align}
Y =& D\theta_0 + X'\Lambda_0 + U,&&E [U|X,D]=0, \label{eq:y}
\end{align}
where $Y$ denotes an outcome variable, $D$ denotes a treatment variable, $X$ denotes a $p$-dimensional vector of controls, and $U$ denotes unobserved factors.
We allow for a high dimensionality in the sense that $p$ can be increasing in $N$ and may be even larger than $N$ -- more details will follow.
In this framework, we are interested in the partial effect $\theta_0$ of $D$ on $Y$.
Also write the linear projection of $D$ on $X$:
\begin{align}
D =& X'\beta_0 + V,&&E [V|X]=0, \label{eq:d}
\end{align}
In Section \ref{sec:approx}, we consider an extended model in which we introduce approximation errors in \eqref{eq:y}--\eqref{eq:d}.
To construct a moment restriction under \eqref{eq:y}--\eqref{eq:d}, consider the orthogonal score function from \cite{Robinson}:
\begin{align}\label{eq:robinson}
\psi(Y,D,X;\theta,\eta) := \left\{Y-X'\gamma-\theta(D-X'\beta) \right\}\left(D-X'\beta \right),
\end{align}
where $X'\gamma_0 = E[Y|X]$ and $\eta=(\gamma,\beta)$.
Note also that $X'\beta_0 = E[D|X]$ follows by construction from \eqref{eq:d}.
\begin{notation}
To proceed, we first fix basic notations.
We use subscripts $i$ and $j$ to denote indices of observations and coordinates, respectively.
Define $X_{I j}=(X_{ij},\,i\in I)$ as a $\left\vert I \right\vert \times 1$ vector, $X_{i J} = (X_{i j},\,j\in J)$ as a $\left\vert J \right\vert\times 1$ vector, and $X_{I J} = (X_{i J},\,i\in I)'$ as a $\left\vert I \right\vert \times \left\vert J \right\vert$ matrix,
where $I$ is a subset of observation indices $\{1,2,...,N\}$, $J$ is a subset of coordinate indices $\mathfrak{P} \equiv \{1,2,\dots,p\}$, and $\left\vert . \right\vert$ denotes the set cardinality.
For any vector, $\left \|.\right \|$ refers to the Euclidean norm.
The $L^q$ norm is defined by $\left \|\xi\right \|_q=(\sum_{j=1}^{p} \xi_j^q)^{1/q}$ for $q < \infty$ and $\left \|\xi\right \|_\infty = \max_{1\le j \le p}\left\vert \xi_j \right\vert$.
\end{notation}
\subsection{The Method}
This section provides an overview of the method.
We propose the following procedure for a root-$N$ consistent estimation and inference about the partial effect $\theta_0$ without assuming sparsity on the high-dimensional parameters, $\beta_0$ or $\gamma_0$.
\begin{algorithm}[OGA+HDAIC with DML for high-dimensional linear models]\label{algorithm:dml_oga_hdaic} \hfill
\begin{enumerate}[label={Step \arabic*.}]
\item Randomly split the sample indices $\{1,...,N\}$ into $K$ folds $(I_k)_{k=1}^K$. For simplicity, let the size of each fold be $n=N/K$ and the size of $I_k^c$ be $n^c$.
\item For each fold $k \in \{1,...,K\}$, perform following procedure using $\{(X_i',D_i)'\}_{i\in I_k^c}$ to get $\widehat\beta_k$.
\begin{enumerate}
\item Compute $\widehat{\mu}_{0,j} = X_{I_k^c j}'D_{I_k^c}/\sqrt{n^c}\lVert X_{I_k^c j} \rVert$. Select the coordinate $\widehat{j}_{1} = \operatorname{argmax}_{1\le j \le p} \left\vert \widehat{\mu}_{0,j} \right\vert$. Define $\widehat{J}_1 = \{\widehat{j}_1\}$.
\item Compute $\widehat{\mu}_{1,j} = X_{I_k^c j}'(I_{n^c} - H_1) D_{I_k^c} / \sqrt{n^c}\lVert X_{I_k^c j} \rVert$, where $H_1 = X_{I_k^c \widehat{j}_1}(X_{I_k^c \widehat{j}_1}'X_{I_k^c, \widehat{j}_1})^{-1}X_{I_k^c, \widehat{j}_1}'$. Select the coordinate $\widehat{j}_2 = \operatorname{argmax}_{1\le j \le p,j\notin \widehat{J}_1} \left\vert \widehat{\mu}_{1,j} \right\vert$. Update $\widehat{J}_2 =\widehat{J}_1\cup \{\widehat{j}_2\}$.
\item Given $m-1$ coordinates $\widehat{J}_{m-1}$ that have been obtained, compute $\widehat{\mu}_{m-1,j} = X_{I_k^c j}'(I_{n^c}- H_{m-1}) D_{I_k^c} / \sqrt{n^c}\lVert X_{I_k^c j} \rVert$, where $H_{m-1} = X_{I_k^c \widehat{J}_{m-1}}(X_{I_k^c \widehat{J}_{m-1}}'X_{I_k^c \widehat{J}_{m-1}})^{-1}X_{I_k^c \widehat{J}_{m-1}}'$. Select the coordinate $\widehat{j}_m = \operatorname{argmax}_{1\le j\le p, j\notin \widehat{J}_{m-1}}\left\vert \widehat{\mu}_{m,j} \right\vert$. Iteractively update $\widehat{J}_m = \widehat{J}_{m-1}\cup \{\widehat{j}_m\}$.
\item Compute $\textup{HDAIC}$ $(\widehat{J}_m) = (1+C^* |\widehat{J}_m|\log p/n^c)\widehat{\sigma}_{m}^2$ for each $m$, where $C^*$ is from (\ref{Cstar}) in Appendix \ref{sec:regression:details} and $\widehat{\sigma}_{m}^2 = 1/n^c D_{I_k^c} '(I-H_m)D_{I_k^c}$. Choose $\widehat{m} = \operatorname{argmin}_{1\le m\le M_n^*}$ $\textup{HDAIC} (\widehat{J}_m)$, where $M_n^*$ is defined in (\ref{Mnstar}) in Appendix \ref{sec:regression:details}.
\item With coordinates $\widehat{J}_{\widehat{m}}$, run OLS of $D_i$ on $X_{i\widehat{J}_{\widehat{m}}}$ to get $\widehat{\beta}_k.$
\end{enumerate}
\item Repeat Step 2 with $\{(X_i',Y_i)'\}_{i\in I_k^c}$ instead of $\{(X_i',D_i)'\}_{i\in I_k^c}$, to get $\widehat{\gamma}_k$ for each fold $k \in \{1,...,K\}$.
\item Obtain $\check{\theta}$ as a solution to $1/K \sum_{k=1}^K 1/n \sum_{i\in I_k}\psi(Y_i,D_i,X_i;\check{\theta},\widehat{\eta}_k)=0$ where $\widehat{\eta}_k=(\widehat{\gamma}_k,\widehat{\beta}_k)$ and $\psi$ is defined in \eqref{eq:robinson}.
\item Compute $\widehat{M} = -1/K\sum_{k=1}^K 1/n\sum_{i\in I_k}(D_i - X_i'\widehat{\beta})^2$. Obtain a variance estimator of $\check{\theta}$ as $\widehat{\Omega} = \widehat{M}^{-1}\frac{1}{K} \sum_{k=1}^K \frac{1}{n}\sum_{i\in I_k} [\psi(Y,D,X;\check{\theta},\widehat{\eta}_k)\psi(Y,D,X;\check{\theta},\widehat{\eta}_k)'](\widehat{M}^{-1})'$.
\end{enumerate}
\end{algorithm}
We highlight three notable elements of this algorithm.
First, the overall procedure (Steps 1--4) uses the cross fitting to remove an over-fitting bias.
Specifically, by using complementary sub-sample $I_k^c$ to estimate the nuisance parameters $\widehat\eta_k=(\widehat\gamma_k,\widehat\beta_k)$ that are in turn evaluated in the $I_k$-mean of the score, we can circumvent a bias that arises from products of dependent factors in the score.
Our combined use of the orthogonal score \eqref{eq:robinson} and this cross-fitting method allows for the high-level theory of the double/debiased machine learning \citep[DML,][]{ccddhnr} to be applicable.
Section \ref{sec:without_cross_fitting} discusses an alternative algorithm that does not rely on the cross fitting at the cost of an additional assumption.
Second, the coordinates $\{\widehat j_1,...,\widehat j_{p}\}$ are ranked in Step 2 (a)--(c) in the order of decreasing importance after successive orthogonalization using OGA as in \cite{ing2020}.
Third, a subset $\widehat J_{\widehat m} = \{\widehat j_1,...,\widehat j_{\widehat m}\}$ of the ordered set $\{\widehat j_1,...,\widehat j_{p}\}$ is selected in Step 2 (d) using HDAIC as in \cite{ing2020}.
Our combined use of these three elements (DML, OGA, and HDAIC) together allows for a novel root $N$ consistent estimation of $\theta_0$ without assuming traditional functional class restrictions (e.g., the sparsity) required by existing popular estimators (e.g., LASSO).
In Section \ref{sec:regression:theory}, we formally present theoretical arguments in support of this claim.
While Algorithm \ref{algorithm:dml_oga_hdaic} provides nearly full details of the proposed method, it omits a couple of details.
Specifically, Step 2 (d) on HDAIC uses two tuning parameters, $C^*$ and $M_n^*$.
We present details about these aspects of the algorithm in Appendix \ref{sec:regression:details}.
\subsection{The Theory}\label{sec:regression:theory}
This section proposes and discusses assumptions under which one can conduct an inference about $\check{\theta}$ based on root-$N$ asymptotic normality using the method described in Algorithm \ref{algorithm:dml_oga_hdaic}.
We use the notations $c,C,\overline{C},\overline{\tau}$ and $q$ for strictly positive constants such that their values can differ depending on the location. Let $q>4$ be a positive integer, $c_q$, $C_q$, $\lambda_1$ be some positive constants and $K_{N,q}$ be a positive sequence of constants such that $K_{N,q}\ge E[\max_{1\le j\le p}\left\vert X_{ij} \right\vert^q]$.
Wherever there is no risk of confusion,
we also use the generic notation $\xi$ to refer to both $\beta$ and $\gamma$ to avoid repetitions. All the random variables and parameter vectors are $N$-dependent unless otherwise specified. We abbreviate the $N$ index for brevity.
\begin{assumption}\label{a:1}
For each $N\in \mathbb N$, it holds that
\begin{enumerate}[(a)]
\item $(Y_i,D_i,X_i')_{i=1}^N$ are i.i.d. copies of $(Y,D,X')$.
\item \eqref{eq:y} and \eqref{eq:d} hold.
\item $E[\left\vert Y \right\vert^q] + E[\left\vert D \right\vert^q] \le C_q$.
\item $E[\left\vert UV \right\vert^2]\ge c_q^2$ and $E[V^2]\ge c_q$.
\item ${\max_{1\le j\le p}E[|X_{ij}|^q]}\le C_q$, $E[\left\vert V \right\vert^q] \le C_q$, and $E[\left\vert U \right\vert^q]\le C_q$. \label{qth_bdd}
\end{enumerate}
Furthermore, it holds asymptotically that
(f) $K_{N,q}^2 \log p/ N^{1-2/q}=o(1).$ \label{maxqth_bdd}
\end{assumption}
Assumption \ref{a:1} (a) requires a random sampling of data.
Assumption \ref{a:1} (b) requires that the correct model is given by \eqref{eq:y} and \eqref{eq:d}.
Assumption \ref{a:1} (c)--(e) requires bounded moments of various variables.
Assumption \ref{a:1} (f) requires constraints on the speed at which the dimensionality as well as the maximal of the covariate vector can grow. Vectors consist of independent subgaussian random variables with bounded variances, for example, are special cases satisfying this restriction, as the expectation of the maximum of $X_{ij}$ is bounded by a factor of $\sqrt{\log p}$.
We emphasize that these assumptions are mild in comparison with the counterpart assumptions made in the high-dimensional regression literature.
\begin{assumption}\label{a:2}
It holds over $N\in \mathbb N$ that
\begin{enumerate}[(a)]
\item $\lambda_{\min} (\Gamma)\ge \lambda_1>0$ and $\lambda_{\max} (\Gamma) \le C_q$, where $\Gamma=E[X X']$. \label{additionalA}
\item Define $\Gamma(J) = E[X_{iJ}X_{iJ}']$ and $d_{\ell}(J) = E[X_{i\ell}X_{iJ}]$ for a set of coordinate indices $J\subseteq \mathfrak{P}$. Then $$\max_{1\le \left\vert J \right\vert\le \overline{C}(N/\log p)^{1/2},\,\ell \notin J} \left\vert \Gamma^{-1}(J)d_\ell(J) \right\vert< C_q.$$ \label{IngA5}
\end{enumerate}
\end{assumption}\vspace{-.6cm}
Assumption \ref{a:2} \ref{additionalA} requires that the minimum eigenvalue $\lambda_{\min}(\Gamma)$ of the Gram matrix $\Gamma$ to be positive, and it is a very common restriction.
Assumption \ref{a:2} \ref{IngA5} is a restriction on the covariance structure of $X_{iJ}$.
Observe that $\Gamma^{-1}(J)d_\ell(J)$ takes the form of regression coefficient of $X_{i\ell}$ on $X_{iJ}$, and so Assumption \ref{a:2} \ref{IngA5} means that $X_{i\ell}$ cannot be strongly correlated with $X_{iJ}$ for $\ell \notin J$.
We remark that these conditions are imposed at the population level; unlike in the LASSO or Dantzig selector \citep{candes2007dantzig}, a restricted eigenvalue type condition for the sample Gram matrix \citep[see e.g.,][]{bickel2009simultaneous} is not required here -- see also the discussion in \citet[Sec. 3.2]{ing2020}.
The following assumption imposes restrictions on the function classes in terms the parameters $\beta_0$ and $\gamma_0$.
We will use the generic notation $\xi_0$ to refer to $\beta_0$ and $\gamma_0$.
Note that $\xi_0 \in \mathbb{R}^p$ in both cases.
Define $\xi(J) = (\xi_j)_{j\in J}$ to be a $\left\vert J \right\vert\times 1$ vector, where recall that $\xi$ is a generic notation to refer to $\beta_0$ and $\gamma_0$.
\begin{assumption}\label{a:3}
It holds over $N\in \mathbb N$ that
for each of $\xi_0 = \beta_0$ and $\gamma_0$, $\xi_0$ follows either (a) or (b) described below.
\begin{enumerate}[(a)]
\item Polynomial decay: $\log p = o(N^{1-2/q})$. Each $\xi_0$ is such that $\left \|\xi_0\right \|_2^2 \le C_0$ for some $C_0>0$ and there exist $\alpha> 1$ such that for any $J \subseteq \mathfrak{P}$,$$\hspace{2cm}\left \|\xi_0(J)\right \|_1 \le C \left( \left \|\xi_0(J)\right \|_2^2\right)^{(\alpha-1)/(2\alpha-1)}.\label{poly} $$
\item Exponential decay: $\log p= o(N^{1/4}).$ Each $\xi_0$ is such that $ \left \|\xi_0\right \|_\infty \le C_0$ for some $C_0>0$ and there exists $C_1>1$ such that for any $J \subseteq \mathfrak{P}$, $$\left \|\xi_0(J)\right \|_1 \le C_1 \left \|\xi_0(J)\right \|_\infty.\label{expo}$$
\end{enumerate}
\end{assumption}
This is a key assumption in this paper, and defines admissible function classes for the high-dimensional linear models.
While the literature on LASSO requires the exact sparsity and the $L^p$ sparsity (including approximate sparsity) conditions, Assumption \ref{a:3} does not impose such conditions.
We remark that, if we rearrange the components of the parameter vector $\xi_0$ by their absolute values in a descending order (denote it again as $\xi_0$ with an abuse of notation), then Condition (a) contains special cases such as the conventional polynomial decay condition
\begin{align*}
Lj^{-\alpha} \le |\xi_{0j}|\le Uj^{-\alpha}, \quad 0<L\le U <\infty,
\end{align*}
as well as the polynomial summability condition
\begin{align*}
\sum_{j=1}^p |\xi_{0j}|^{1/\alpha} < M,\quad M\in(0,\infty),
\end{align*}
following the discussion in \citet[pp. 1962]{ing2020}. On the other hand, Condition (b) implies the conventional exponential decay condition that, for some $\alpha'>0$,
\begin{align*}
L'\exp(-\alpha' j) \le |\xi_{0j}|\le U'\exp(-\alpha' j) ,\quad 0<L'\le U' <\infty,
\end{align*}
so long as the regressors have bounded second moments.
Hence throughout the paper, Conditions (a) and (b) are referred to as the polynomial decay condition and the exponential decay condition, respectively, albeit their extra generality. Clearly, the case of polynomial decay accommodates a larger function class, but we remark that there is a tradeoff in terms of how fast the dimension $p$ can diverge as the sample size $N$ increases.
Following \citet[][pp. 1960]{ing2020}, we now present a concrete example where our Assumption \ref{a:3} holds but the sparsity does not.
Suppose that
\begin{align*}
L j^{-\alpha} \le \left\vert \Lambda_{0(j)}\sigma_{(j)} \right\vert \le U j^{-\alpha},&&j=1,\dots,p,
\end{align*}
for some $\alpha>1$, where $0< L \le U < \infty$, and $|\Lambda_{0(1)}\sigma_{(1)}|\ge |\Lambda_{0(2)}\sigma_{(2)}| \ge \dots \ge |\Lambda_{0(p)}\sigma_{(p)}|$ is a descending reordering of $\{\Lambda_{0,j}\sigma_j\}$ with $\sigma_j^2 = E[X_{ij}^2]$.
In this setting, our Assumption \ref{a:3} holds for the same $\alpha$, but $\sum_{j=1}^p \left\vert \Lambda_{0,j}\sigma_j \right\vert^{1/\gamma}$ is now unbounded as $L(1+\log p)\le \sum_{j=1}^p \left\vert \Lambda_{0,j}\sigma_j \right\vert^{1/\gamma} \le U(1+\log p).$
\begin{theorem}\label{theorem:regression} Let $(\mathcal{P}_N)_{N\in \mathbb N}$ be a sequence of sets of DGPs such that Assumptions \ref{a:1}--\ref{a:3} are satisfied on the model \eqref{eq:y}--\eqref{eq:d}. Then, the estimator $\check{\theta}$ satisfies
\begin{align*}
\sqrt{N} \left(\check{\theta}-\theta_0\right) \xrightarrow{d} N(0,\Omega),
\end{align*}
where $\Omega = (E[V^2])^{-1}E[V^2U^2](E[V^2])^{-1}$.
Define $\widehat{M} :=- 1/K\sum_{k=1}^K 1/n \sum_{i\in I_k} (D_i-X_i'\widehat{\beta})^2$. Then, we can define the variance estimator \begin{align*}
\widehat{\Omega} = \widehat{M}^{-1}\frac{1}{K} \sum_{k=1}^K \frac{1}{n}\sum_{i\in I_k} [\psi(Y,D,X;\check{\theta},\widehat{\eta}_k)\psi(Y,D,X;\check{\theta},\widehat{\eta}_k)'](\widehat{M}^{-1})'
\end{align*}
and the confidence regions with significance level $a\in (0,1 )$ have uniform asymptotic validity:
\begin{align*}
\sup_{P\in \mathcal{P}_N} \left\vert P\left(\theta_0 \in \left[\check{\theta} \pm \Phi^{-1}(1-a/2)\sqrt{\widehat{\Omega}/N}\right]\right)-(1-a) \right\vert =o(1).
\end{align*}
\end{theorem}
A proof is provided in Appendix \ref{sec:theorem:regression}.
This theorem guarantees that the estimator $\check{\theta}$ of $\theta_0$ provided by Algorithm \ref{algorithm:dml_oga_hdaic} converges at the rate of $\sqrt{N}$ and asymptotically follows the normal distribution under Assumption \ref{a:1}--\ref{a:3}.
Furthermore, the sample-counterpart asymptotic variance estimator constructs an asymptotically valid confidence interval.
We emphasize that this result does not rely on the sparsity assumption which is used in the literature on high-dimensional linear models.
\section{Simulation Studies}\label{sec:simulations}
In this section, we investigate the finite sample properties of our proposed estimator $\check{\theta}$ and compare them with those of two existing estimators, namely the LASSO-based DML and random-forest-based DML.\footnote{For these two existing methods, we use the R package ``DoubleML : Double Machine Learning in R.''}
We follow \cite{BCH2014RES} in developing baseline data generating processes (DGPs).
The linear regression model is specified by
\begin{align*}
Y =& D\theta_0 + X'\Lambda_0 + U,
\end{align*}
where $\theta_0=0.5$ and $p = dim(X) = 500$.
Consistently with this specification, data are generated by the system
\begin{align*}
Y =& \theta_0(D-X'\beta_0) + X'\gamma_0 + U, &&U\sim N(0,1),
\\
D =& X'\beta_0 + V, &&V \sim N(0,1),
\end{align*}
where the covariates are in turn generated by $X \sim N(0,\Sigma)$ with $\Sigma_{jk}=(0.5)^{\left\vert k-j \right\vert}$.
For the high-dimensional nuisance parameters, $\eta_0 = (\gamma_0,\beta_0)$, we set $p=500$ throughout and consider a couple of alternative designs.
In the first design, each of $\beta_0$ and $\gamma_0$ has ten coordinates taking the value of 1 and $p-10$ coordinates taking the value of zero, i.e., sparse deign.
In the second design, both $\beta_0$ and $\gamma_0$ decay exponentially.
Specifically, the $j$-th coordinate of each of $\beta_0$ and $\gamma_0$ is set to $e^{-j}$.
The third design has both $\beta_0$ and $\gamma_0$ decaying at polynomial rates.
Specifically, the $j$-th coordinate of each of $\beta_0$ and $\gamma_0$ is set to $j^{-2}$, $j^{-1.75}$, $j^{-1.5}$, $j^{-1.25}$ and $j^{-1}$ for five sets of simulations.
For each of these sets of simulations, we experiment with the two sample sizes $N \in \{500,1000\}$.
\begin{table}
\centering
\scalebox{1}{
\begin{tabular}{ccclcccc}
\hline\hline
$\beta_{0,j},\gamma_{0,j}$ & $N$ & $p$ & Method of Preliminary Estimation & Bias & SD & RMSE & 95\%\\
\hline
Sparse & 500 & 500 & LASSO & 0.020 & 0.044 & 0.049 & 0.937\\
& & & Random Forest & -0.481& 0.000 & 0.498 & 0.000\\
& & & OGA+HDAIC & -0.003 & 0.045 & 0.046 & 0.943\\
\cline{2-8}
&1000 & 500 & LASSO & 0.012 & 0.031 & 0.033 & 0.929\\
& & & Random Forest & -0.347& 0.003 & 0.485 & 0.000\\
& & & OGA+HDAIC & 0.000 & 0.032 & 0.032 & 0.947\\
\hline
$e^{-j}$ & 500 & 500 & LASSO & 0.006 & 0.044 & 0.045 & 0.934\\
& & & Random Forest & 0.008 & 0.044 & 0.045 & 0.940\\
& & & OGA+HDAIC & 0.000 & 0.045 & 0.045 & 0.941\\
\cline{2-8}
&1000 & 500 & LASSO & 0.007 & 0.031 & 0.032 & 0.942\\
& & & Random Forest & 0.006 & 0.031 & 0.031 & 0.950\\
& & & OGA+HDAIC & 0.000 & 0.032 & 0.032 & 0.950\\
\hline
$j^{-2}$ & 500 & 500 & LASSO & 0.010 & 0.044 & 0.046 & 0.934\\
& & & Random Forest & 0.033 & 0.043 & 0.055 & 0.878\\
& & & OGA+HDAIC &-0.002 & 0.045 & 0.046 & 0.938\\
\cline{2-8}
&1000 & 500 & LASSO & 0.009 & 0.031 & 0.032 & 0.939\\
& & & Random Forest & 0.025 & 0.031 & 0.039 & 0.882\\
& & & OGA+HDAIC & 0.001 & 0.032 & 0.032 & 0.945\\
\hline
$j^{-1.75}$& 500 & 500 & LASSO & 0.016 & 0.044 & 0.047 & 0.931\\
& & & Random Forest & 0.046 & 0.043 & 0.063 & 0.818\\
& & & OGA+HDAIC & -0.001 & 0.045 & 0.047 & 0.930\\
\cline{2-8}
&1000 & 500 & LASSO & 0.011 & 0.031 & 0.033 & 0.932\\
& & & Random Forest & 0.035 & 0.031 & 0.046 & 0.805\\
& & & OGA+HDAIC & 0.001 & 0.031 & 0.032 & 0.951\\
\hline
$j^{-1.5}$& 500 & 500 & LASSO & 0.020 & 0.044 & 0.049 & 0.931\\
& & & Random Forest & 0.066 & 0.042 & 0.079 & 0.651\\
& & & OGA+HDAIC & 0.001 & 0.045 & 0.046 & 0.936\\
\cline{2-8}
&1000 & 500 & LASSO & 0.013 & 0.031 & 0.034 & 0.928\\
& & & Random Forest & 0.053 & 0.030 & 0.060 & 0.576\\
& & & OGA+HDAIC & 0.002 & 0.031 & 0.033 & 0.938\\
\hline
$j^{-1.25}$& 500 & 500 & LASSO & 0.028 & 0.044 & 0.053 & 0.904\\
& & & Random Forest & 0.108 & 0.041 & 0.115 & 0.245\\
& & & OGA+HDAIC & 0.006 & 0.044 & 0.047 & 0.933\\
\cline{2-8}
&1000 & 500 & LASSO & 0.018 & 0.031 & 0.036 & 0.909\\
& & & Random Forest & 0.096 & 0.029 & 0.100 & 0.083\\
& & & OGA+HDAIC & 0.004 & 0.031 & 0.034 & 0.923\\
\hline
$j^{-1}$ & 500 & 500 & LASSO & 0.038 & 0.043 & 0.061 & 0.827\\
& & & Random Forest & 0.193 & 0.037 & 0.196 & 0.001\\
& & & OGA+HDAIC & 0.022 & 0.043 & 0.053 & 0.893\\
\cline{2-8}
&1000 & 500 & LASSO & 0.026 & 0.031 & 0.042 & 0.838\\
& & & Random Forest & 0.179 & 0.026 & 0.181 & 0.000\\
& & & OGA+HDAIC & 0.014 & 0.031 & 0.037 & 0.901\\
\hline\hline
\end{tabular}
}
\caption{\setlength{\baselineskip}{5.5mm}Monte Carlo simulation results. Displayed are Monte Carlo simulation statistics including the bias, standard deviation (SD), root mean square error (RMSE), and 95\% coverage frequency.}
\label{tab:simulation}
\end{table}
Table \ref{tab:simulation} summarize simulation results.
Displayed are four Monte Carlo simulation statistics for each set of simulations, including the bias, standard deviation (SD), root mean square error (RMSE), and 95\% coverage frequency.
In the first row group of the table displaying the results the sparse design, both LASSO-based method and our proposed method based on the OGA and HDAIC work well, while that based on Random Forest significantly underperforms.
In the second row group of the table displaying the results under the exponential decay,
all the three machine learning methods yield desired results both in terms of all the displayed statistics.
There are no significant differences across the three methods under this sparse model.
In the subsequent row groups of the table displaying the results for the cases of the polynomial decays, however, observe that the performance varies across the three machine learning methods.
While our proposed method based on the OGA and HDAIC continues to perform well in terms of all the displayed statistics, the LASSO-based method slightly underperforms and the random-forest-based method significantly underperforms.
In particular, these differences in the finite-sample performance widen as the degree of polynomial decay becomes smaller, i.e., as the model becomes less sparse.
These results demonstrate the relative robustness of the method proposed in this paper under less sparse high-dimensional regression models.
We ran many other sets of simulations and present their results in Appendix \ref{sec:additional_simulations}.
In particular, Appendix \ref{sec:simulation:tuning} presents simulation results with various values of the tuning parameters, and demonstrate the robustness of the qualitative patterns observed above.
From these results, we recommend to use the method based on the OGA and HDAIC over the two alternative methods for its robust performance across various designs even including the sparse design.
\section{Extensions} \label{sec:extensions}
\subsection{Models with Approximation Errors}\label{sec:approx}
Extending the baseline model \eqref{eq:y}--\eqref{eq:d}, consider the following partially linear model motivated by \cite{BCH2014RES}:
\begin{align}
Y =& D\theta_0 + f(X) + U,&&E [U|X]=0, \label{eq:y:approx:pl}\\
D =& g(X) + V,&&E [V|X]=0, \label{eq:d:approx:pl}
\end{align}
where $Y$ denotes an outcome variable, $D$ denotes a treatment variable, $X$ denotes a $p$-dimensional vector of controls, and $U$ and $V$ denotes unobserved factors.
We do not directly impose any parametric restriction on $f$ or $g$ unlike the baseline model presented in Section \ref{sec:regression}.
This extension is useful in certain applications, such as the one we present in Section \ref{eq:application}.
In this semi-parametric framework, we are interested in the partial effect $\theta_0$ of $D$ on $Y$.
Now consider the followng reduced form regressions for \eqref{eq:y:approx:pl}--\eqref{eq:d:approx:pl}:
\begin{align}
Y =& \underbrace{X'\gamma_0 + r_Y(X)}_{f(X)+\theta_0 g(X)} + \operatorname{\mathcal{E}},&&E [\operatorname{\mathcal{E}}|X]=0, \label{eq:y:approx}\\
D =& \underbrace{X'\beta_0 + r_D(X)}_{g(X)} + V,&&E [V|X]=0, \label{eq:d:approx}
\end{align}
where $X'\gamma_0$ and $X'\beta_0$ are approximations to $E[Y|X]$ and $E[D|X]$, and $r_{Y}(X)$ and $r_D(X)$ are approximation errors.
The functions $r_{Y}$ and $r_{D}$ are nonparametric as are $f$ and $g$. {\color{black}We will impose conditions on the magnitudes of $r_Y$ and $r_D$ below. Models under these conditions, along with certain sparsity conditions imposed on $\beta_0$ and $\gamma_{0}$, are said to be ``approximate sparse'' in \cite{belloni2012sparse,BCH2014RES}.}
Recall the orthogonal score $\psi(Y,D,X;\theta,\eta)$ defined in \eqref{eq:robinson}.
With this orthogonal score, we propose to obtain $\check\theta$ and $\widehat\Omega$ via Algorithm \ref{algorithm:dml_oga_hdaic} presented in Section \ref{sec:regression} even under the current extended setting with approximation errors.
With the extended model \eqref{eq:y:approx}--\eqref{eq:d:approx}, a different set of assumptions are imposed from those in the baseline model.
First, we slightly modify Assumption \ref{a:1} as follows.
\begin{assumption}\label{a:1:approx}
For each $N\in \mathbb N$, it holds that
\begin{enumerate}[(a)]
\item $(Y_i,D_i,X_i')_{i=1}^N$ are i.i.d. copies of $(Y,D,X')$.
\item \eqref{eq:y:approx} and \eqref{eq:d:approx} hold.
\item $E[\left\vert Y \right\vert^q] + E[\left\vert D \right\vert^q] \le C_q$.
\item $E[\left\vert UV \right\vert^2]\ge c_q^2$ and $E[V^2|(Y,D,X')]\ge c_q$. \label{Vsq_bdd_below}
\item ${\max_{1\le j\le p}E[|X_{ij}|^q]}\le C_q$, $E[\left\vert V \right\vert^q] \le C_q$, and $E[\left\vert \operatorname{\mathcal{E}} \right\vert^q]\le C_q$. \label{qth_bdd}
\end{enumerate}
Furthermore, it holds asymptotically that
(f) $K_{N,q}^2 C\log p/ N^{1-2/q}=o(1).$ \label{maxqth_bdd}
\end{assumption}
\noindent
In part \ref{Vsq_bdd_below}, we require the conditional variance of $V$ given $(Y,D,X')$ to be bounded away from zero whereas the counterpart in the baseline model assumed the unconditional variance to be bounded away from zero.
We continue to use Assumption \ref{a:2} from the baseline model.
However, it should be stressed that we now impose Assumption \ref{a:2} on \eqref{eq:y:approx}--\eqref{eq:d:approx} rather than \eqref{eq:y}--\eqref{eq:d}.
With the approximation errors introduced in the current extended model, we make the following assumption on the approximation error functions $r_Y$ and $r_D$.
\begin{assumption}\label{a:approx}
For $r(X) = r_Y(X)$ and $r_D(X)$, it holds that
\begin{enumerate}[(a)]
\item $E[r^4(X)]\le C$. \label{r4_bdd}
\item $E[r^2(X)]\le C\log p/N$. \label{rsq_bdd}
\item $\max_{1\le j\le p}\left|E[r(X)X_{ij}]\right|\le C_{p,1}\sqrt{\log p}/N^{1/4}$.\label{rX_bdd}
\end{enumerate}
\end{assumption}
\noindent
Assumption \ref{a:approx} \ref{r4_bdd} requires the fourth moment of the approximation error to be bounded, \ref{rsq_bdd} assumes the second moment to be of order $\log p/N$, and \ref{rX_bdd} bounds the maximum cross moment of the approximation and the covariates.
Finally, we focus on the more difficult case, namely the polynomial decay case, for brevity in this section.
\begin{assumption}\label{a:3:approx}
It holds over $N\in \mathbb N$ that
for each of $\xi_0 = \beta_0$ and $\gamma_0$, $\xi_0$ follows
polynomial decay, i.e., $\log p = o(N^{1-2/q})$. Each $\xi_0$ is such that $\left \|\xi_0\right \|_2^2 \le C_0$ for some $C_0>0$ and there exist $\alpha> 1$ and $C_{\alpha} > 0$ such that for any $J \subseteq \mathfrak{P}$,$$\hspace{2cm}\left \|\xi_0(J)\right \|_1 \le C_\alpha \left( \left \|\xi_0(J)\right \|_2^2\right)^{(\alpha-1)/(2\alpha-1)}.$$
\end{assumption}
The following theorem establishes the asymptotic normality of $\check\theta$ along with the asymptotic validity of inference under the extended model with approximation errors.
\begin{theorem}\label{theorem:approx} Let $(\mathcal{P}_N)_{N\in \mathbb{N}}$ be a sequence of sets of DGPs such that Assumptions \ref{a:2} and \ref{a:1:approx}--\ref{a:3:approx} are satisfied on the model \eqref{eq:y:approx:pl}--\eqref{eq:d:approx:pl} entailing the reduced forms \eqref{eq:y:approx}--\eqref{eq:d:approx}. Then, the estimator $\check{\theta}$ defined in Algorithm \ref{algorithm:dml_oga_hdaic} satisfies
\begin{align*}
\sqrt{N} \left(\check{\theta}-\theta_0\right) \xrightarrow{d} N(0,\Omega),
\end{align*}
where $\Omega = (E[V^2])^{-1}E[V^2U^2](E[V^2])^{-1}$.
The confidence regions with significance level $a\in (0,1)$ have uniform asymptotic validity:
\begin{align*}
\sup_{P\in \mathcal{P}_N} \left\vert P\left(\theta_0 \in \left[\check{\theta} \pm \Phi^{-1}(1-a/2)\sqrt{\widehat{\Omega}/N}\right]\right)-(1-a) \right\vert =o(1),
\end{align*}
where $\widehat{\Omega}$ is defined in Section \ref{sec:regression}.
\end{theorem}
\noindent
A proof is presented in Appendix \ref{sec:theorem:approx}.
The same remarks as those presented below the statement of Theorem \ref{theorem:regression} apply here.
{\color{black} We want to stress that Theorem \ref{theorem:approx} is not an immediate consequence given Theorem \ref{theorem:regression}, because the original proofs for convergence rates of OGA+HDAIC in \cite{ing2020} do not permit approximately sparse models. In order to show Theorem \ref{theorem:approx}, we establish convergence rates for the OGA+HDAIC in approximately sparse regression models.}
\subsection{Estimation and Inference without Cross Fitting}\label{sec:without_cross_fitting}
Thus far, our proposed procedures of estimation and inference are based on cross fitting.
A drawback of using the cross fitting is the randomness of estimates given the data.
To overcome this drawback, we provide an alternative procedure of estimation and inference without relying on cross fitting in this section.
However, we stress that this benefit comes with costs in some assumptions as discussed below.
We continue from Section \ref{sec:approx} to consider the partial linear model \eqref{eq:y:approx:pl}--\eqref{eq:d:approx:pl} entailing the reduced forms \eqref{eq:y:approx}--\eqref{eq:d:approx}.
Furthermore, we continue to use the same orthogonal score $\psi(Y,D,X;\theta,\eta)$ defined in \eqref{eq:robinson}.
However, we now replace Algorithm \ref{algorithm:dml_oga_hdaic} by the following algorithm which does not involve the cross-fitting procedure.
Let $[N] = \{1,\dots,N\}$, so that $X_{[N] j} = \{X_{ij}, i\in [N]\}$ and $D_{[N]}=(D_1,\dots,D_N)'$.
\begin{algorithm}[OGA+HDAIC with DML for high-dimensional linear models without cross fitting]\label{algorithm:dml_oga_hdaic:fullsample} \hfill
\begin{enumerate}[label={Step \arabic*.}]
\item Perform following procedure using $\{(X_i',D_i)'\}_{i=1}^N$ to get $\widehat{\beta}$.
\begin{enumerate}
\item Compute $\widehat{\mu}_{0,j} = X_{[N] j}'D_{[N]}/\sqrt{N}\lVert X_{[N] j} \rVert$. Select the coordinate $\widehat{j}_{1} = \operatorname{argmax}_{1\le j \le p} \left\vert \widehat{\mu}_{0,j} \right\vert$. Define $\widehat{J}_1 = \{\widehat{j}_1\}$.
\item Compute $\widehat{\mu}_{1,j} = X_{[N] j}'(I_N- H_1) D_{[N]} / \sqrt{N}\lVert X_{[N] j} \rVert$, where $H_1 = X_{[N] \widehat{j}_1}(X_{[N] \widehat{j}_1}'X_{[N] \widehat{j}_1})^{-1}X_{[N] \widehat{j}_1}'$. Select the coordinate $\widehat{j}_2 = \operatorname{argmax}_{1\le j \le p,j\notin \widehat{J}_1} \left\vert \widehat{\mu}_{1,j} \right\vert$. Update $\widehat{J}_2 =\widehat{J}_1\cup \{\widehat{j}_2\}$.
\item Given $m-1$ coordinates $\widehat{J}_{m-1}$ that have been obtained, compute $\widehat{\mu}_{m-1,j} = X_{[N] j}'(I_N - H_{m-1}) D_{[N]} / \sqrt{N}\lVert X_{[N] j} \rVert$, where $H_{m-1} = X_{[N] \widehat{J}_{m-1}}(X_{[N] \widehat{J}_{m-1}}'X_{[N] \widehat{J}_{m-1}})^{-1}X_{[N] \widehat{J}_{m-1}}'$. Select the coordinate $\widehat{j}_m = \operatorname{argmax}_{1\le j\le p, j\notin \widehat{J}_{m-1}}\left\vert \widehat{\mu}_{m,j} \right\vert$. Iteractively update $\widehat{J}_m = \widehat{J}_{m-1}\cup \{\widehat{j}_m\}$.
\item Compute $\textup{HDAIC}$ $(\widehat{J}_m) = (1+C^* |\widehat{J}_m|\log p/N)\widehat{\sigma}_{m}^2$ for each $m$ and $\widehat{\sigma}_{m}^2 = 1/N D '(I-H_m)D$. Choose $\widehat{m} = \operatorname{argmin}_{1\le m\le M_n^*}$ $\textup{HDAIC} (\widehat{J}_m)$, where $C^*$ and $M_n^*$ are defined in \eqref{Cstar} and \eqref{Mnstar} in Appendix \ref{sec:regression:details}.
\item With coordinates $\widehat{J}_{\widehat{m}}$, run OLS of $D_i$ on $X_{i\widehat{J}_{\widehat{m}}}$ to get $\widehat{\beta}.$
\end{enumerate}
\item Repeat Step 2 with $\{(X_i',Y_i)'\}_{i=1}^N$ instead of $\{(X_i',D_i)'\}_{i=1}^N$, to get $\widehat{\gamma}$.
\item Obtain $\widetilde{\theta}$ as a solution to $1/N \sum_{i=1}^N \psi(Y_i,D_i,X_i;\widetilde{\theta},\widehat{\eta})=0$ where $\widehat{\eta} = (\widehat{\gamma},\widehat{\beta})$ and $\psi$ is defined in \eqref{eq:robinson}.
\end{enumerate}
\end{algorithm}
To establish asymptotic properties for this new estimator $\widetilde\theta$, we continue to impose Assumptions \ref{a:2}, \ref{a:1:approx}, and \ref{a:approx}.
As in Section \ref{sec:approx}, we focus on the more difficult case, namely the polynomial decay case, for brevity.
\begin{assumption}\label{a:3:fullsample}
It holds over $N\in \mathbb N$ that
$\left\vert \theta_0 \right\vert \le C$, and
for each of $\xi_0 = \beta_0$ and $\gamma_0$, $\xi_0$ follows
polynomial decay, i.e.,
each $\xi_0$ is such that $\left \|\xi_0\right \|_2^2 \le C_0$ for some $C_0>0$ and there exist $\alpha> 1$ such that $\log p = o(N^{(\alpha-1)/(3\alpha-1)})$ and for any $J \subseteq \mathfrak{P}$,$$\hspace{2cm}\left \|\xi_0(J)\right \|_1 \le C \left( \left \|\xi_0(J)\right \|_2^2\right)^{(\alpha-1)/(2\alpha-1)}.$$
\end{assumption}
\noindent
Unlike the previous sections, however, we now require $\log p = o(N^{(\alpha-1)/(3\alpha-1)})$ for $\alpha>1$.
The following theorem establishes the asymptotic normality of $\widetilde\theta$ defined without the cross-fitting procedure.
\begin{theorem}\label{theorem:fullsample} Let $(\mathcal{P}_N)_{N\in \mathbb{N}}$ be a sequence of sets of DGPs such that Assumptions \ref{a:2}, \ref{a:1:approx}, \ref{a:approx}, and \ref{a:3:fullsample} are satisfied on the model \eqref{eq:y:approx:pl}--\eqref{eq:d:approx:pl} entailing the reduced forms \eqref{eq:y:approx}--\eqref{eq:d:approx}. Then, the estimator $\widetilde{\theta}$ satisfies
\begin{align*}
\sqrt{N} \left(\widetilde{\theta}-\theta_0\right) \xrightarrow{d} N(0,\Omega),
\end{align*}
where $\Omega = (E[V^2])^{-1}E[V^2U^2](E[V^2])^{-1}$.
\end{theorem}
\noindent
A proof is given in Appendix \ref{sec:theorem:fullsample}.
{\color{black}We stress that, although the proof builds on that of Theorem 1 in \cite{BCH2014RES}, it is far from being trivial as the lack of cross-fitting and $L^p$-sparsity creates extra challenges. Specifically, a key intermediate step is to control the $L^1$ distances between $\beta_0$ and $\widetilde \beta(\widetilde J)$, an oracle regression estimator defined in the proof of Theorem \ref{theorem:fullsample}. Due to the lack of exact or approximate sparsity of $\beta_0$, this is shown via different strategies from those employed in \cite{BCH2014RES}.
}
As emphasized at the beginning of the current subsection, the main advantage of the estimation procedure without cross fitting is that the estimate is now non-random given data.
Besides, this framework without cross fitting offers an additional advantage.
Recall that our main motivation to use the OGA+HDAIC is to weaken the sparsity assumptions required for conventional high-dimensional methods such as the LASSO.
In the current framework without cross fitting, there is another motivation to use the OGA+HDAIC.
Namely, it selects regressors based on their strength in explanatory power by the algorithm.
Hence, we have better interpretations of the model selected by the OGA+HDAIC than the conventional high-dimensional methods under the current framework without cross-fitting.
We highlight this additional advantage of our proposed method.
Appendix \ref{sec:simulation:no_cross_fitting} presents simulation results with this modified method without cross fitting.
The results are similar to those obtained for the baseline model presented in Section \ref{sec:simulations}.
\section{An Empirical Application}\label{eq:application}
In this section, we demonstrate an application of the proposed method to estimation of production functions.
The main challenge in the econometrics of production functions is the simultaneity in the choice of input firms \citep{marschak1944random}.
While early studies of production functions address this simultaneity problem by explicitly modeling rational choice structures of firms, \citet{OlleyPakes1996} more recently propose a novel idea to use the inverse of the reduced-form investment choice function as a control function.
\citet{levinsohn2003estimating} propose to use intermediate input, instead of investment, as a control variable for a number of advantages.
The use of the control function \textit{a la} \citet{levinsohn2003estimating} entails the partial linear estimating equation for the labor elasticity of the form
\begin{equation}\label{eq:production_first_stage}
y_{it} = \ell_{it}\theta + f(k_{it},m_{it}) + u_{it},
\end{equation}
where
$y_{it}$ denotes the logarithm of output,
$\ell_{it}$ denotes the logarithm of labor input,
$k_{it}$ denotes the logarithm of capital input,
$m_{it}$ denotes the logarithm of intermediate input,
$g$ is a nonparametric function that subsumes a part of the production function and the control function, and
$u_{it}$ denotes a mean-orthogonal reduced-form composite error.
See \citet{OlleyPakes1996} and \citet{levinsohn2003estimating} for details.
In light of the partial linear form \eqref{eq:production_first_stage}, \citet{OlleyPakes1996} and \citet{levinsohn2003estimating} propose to use the estimator of \citet{Robinson} which is semiparametric root-$n$ consistent for $\theta$.
Following these seminal papers, numerous researchers have estimated production functions.
That said, many of these subsequent studies follow the Stata command \citep{petrin2004production} which implements estimation of \eqref{eq:production_first_stage} via the parametric third-degree polynomial approximation
\begin{equation}\label{eq:production_first_stage_degree_three}
y_{it} = \ell_{it}\theta + \sum_{\rho_1=0}^3\sum_{\rho_2=0}^{3-\rho_1} \delta_{\rho_1 \rho_2} k_{it}^{\rho_1} m_{it}^{\rho_2} + u_{it}.
\end{equation}
See \citet[][pages 116--118]{petrin2004production}.
To mitigate the approximation bias asymptotically, we consider a higher-dimensional approximation
\begin{equation}\label{eq:production_first_stage_high_dimensional}
y_{it} = \ell_{it}\theta + \underbrace{\sum_{j=1}^p \delta_{j} \phi_j(k_{it},m_{it}) + r_p(k_{it},m_{it})}_{f(k_{it},m_{it})} + u_{it}
\end{equation}
with an error $r_p(k_{it},m_{it})$ in approximation,
where $\phi = (\phi_1,\phi_2,\phi_3,\cdots)$ is a basis and $p$ can be large and increasing with the sample size.
The basis $\phi$ could be defined as the Cartesian product of polynomials, i.e.,
$(\phi_1(k,m),\phi_2(k,m),$ $\phi_3(k,m),\cdots)=(1,k,m,k^2,m^2,km,\cdots)$, as a generalization of the popular estimating equation \eqref{eq:production_first_stage_degree_three} in the Stata command.
More generally, we can define the basis $\phi$ as the tensor product of orthonomal bases.
We employ the tensor product of Hermite bases \citep{gallant1987semi,chen2007large} for our basis $\phi$, and apply our proposed method to \eqref{eq:production_first_stage_high_dimensional} to get an estimate of $\theta$ and its standard error.
Following \citet{levinsohn2003estimating}, we use a plant-level panel of Chilean firms from 1979 to 1986.
See \citet{liu1991entry} for details about the construction of the data.
Among others, we focus on the 3-digit level industry of food products (311) because of its large sample size compared to other industries.
We are interested in the elasticity with respect to unskilled labor input $\ell^u_{it}$ and skilled labor input $\ell^s_{it}$.
The intermediate input variables include electricity $m^e_{it}$, fuels $m^f_{it}$, and materials $m^m_{it}$.
To estimate the elasticity with respect to unskilled labor input $\ell^u_{it}$ using $m^m_{it}$ as a proxy following \citet{levinsohn2003estimating}, we consider the estimating equation of the form
\begin{equation}\label{eq:unskilled}
\underbrace{y_{it}}_{Y} = \underbrace{\ell^u_{it}\theta^u}_{D\theta_0} + \underbrace{\ell^s_{it}\theta^s + m^e_{it}\theta^e + m^f_{it}\theta^f + \sum_{j=1}^p \delta_{j} \phi_j(k_{it},m^m_{it}) + \tau_{t} + r_p(k_{it},m^m_{it})}_{f(X) } + \underbrace{u_{it}}_{U}
\end{equation}
as in \eqref{eq:y:approx:pl}.
To estimate the elasticity with respect to skilled labor input $\ell^s_{it}$, we swap $\ell^u_{it}\theta^u$ and $\ell^s_{it}\theta^s$ in the above estimating equation:
\begin{equation}\label{eq:skilled}
\underbrace{y_{it}}_{Y} = \underbrace{\ell^s_{it}\theta^s}_{D\theta_0} + \underbrace{\ell^u_{it}\theta^u + m^e_{it}\theta^e + m^f_{it}\theta^f + \sum_{j=1}^p \delta_{j} \phi_j(k_{it},m^m_{it}) + \tau_{t} + r_p(k_{it},m^m_{it})}_{f(X) } + \underbrace{u_{it}}_{U}
\end{equation}
as in \eqref{eq:y:approx:pl}.
The term $\tau_{t}$ represents time effects.
Following \citet{levinsohn2003estimating}, we include the indicator for year groups 1979--1981, 1982--1983, and 1984--1986.
For estimation of \eqref{eq:unskilled} using a polynomial basis, we let $X$ consist of
(i) $\ell_{it}^s$
(ii) $m^e_{it}$,
(iii) $m^f_{it}$,
(iv) $k_{it}$, $\ldots$, $k_{it}^{10}$,
(v) $m^m_{it}$, $\ldots$, $(m^m_{it})^{10}$,
(vi) dummy for 1979--1981,
(vii) dummy for 1982--1983, and
(viii) interactions of the terms in (iv) and (v).
We also consider an estimation of \eqref{eq:unskilled} using a Hermite basis $(\psi_0,\ldots,\psi_9)$, we let $X$ consist of
(i) $\ell_{it}^u$
(ii) $m^e_{it}$,
(iii) $m^f_{it}$,
(iv) $\psi_0(k_{it})$, $\ldots$, $\psi_9(k_{it})$,
(v) $\psi_0(m^m_{it})$, $\ldots$, $\psi_9(m^m_{it})$, and
(vi) dummy for 1979--1981,
(vii) dummy for 1982--1983, and
(viii) interactions of the terms in (iv) and (v).
We use a finite-sample adjusted version of the DML estimates following \citet[][Sec. 3.4]{ccddhnr} -- see Appendix \ref{eq:finite} for details.
See Appendix \ref{eq:hermite} for details about the Hermite basis.
We repeat analogous estimation procedures for \eqref{eq:unskilled}.
Table \ref{tab:empirical} summarizes estimation results.
Row (I) copies estimates from \citet{levinsohn2003estimating}.
Rows (II), (III), and (IV) report results based on the DML with LASSO, DML with random forest, and DML with the OGA and HDAIC (the estimator proposed in this paper), respectively.\footnote{For (II) and (III), we use the R package ``DoubleML : Double Machine Learning in R.'' We set the parameters as folds $=10$, num.trees $=100$, min.node.size $= 2$, max.depth $= 5$, and the number of repetitions $= 20$}
The first two columns show results based on the tensor product of polynomial bases, while the last two columns show results based on the tensor product of Hermite bases.
\begin{table}
\centering
\begin{tabular}{clccccc}
\hline\hline
&& \multicolumn{2}{c}{Polynomial Basis} && \multicolumn{2}{c}{Hermite Basis} \\
\cline{3-4}\cline{6-7}
&& Unskilled & Skilled && Unskilled & Skilled \\
&& Labor & Labor && Labor & Labor \\
\hline
(I) &\citet{levinsohn2003estimating}& 0.139 & 0.051 && --- & ---\\
& &(0.010)&(0.009)\\
\hline
(II) & Double Machine Learning with & 0.170 & 0.063 && 0.196 & 0.085\\
& LASSO Preliminary Estimation &(0.011)&(0.008)&&(0.011)&(0.010)\\
\hline
(III) & Double Machine Learning with & 0.185 & 0.061 && 0.189 & 0.062\\
& Random Forest Preliminary Estimation
&(0.013)&(0.010)&&(0.013)&(0.011)\\
\hline
(IV) & Double Machine Learning with & 0.279 & 0.161 && 0.165 & 0.038\\
& OGA+HDAIC Preliminary Estimation &(0.017)&(0.012)&&(0.011)&(0.010)\\
\hline\hline
\end{tabular}
\caption{Estimates of labor elasticities in the 3-digit level industry of food products (311) in Chile.}
\label{tab:empirical}
\end{table}
First, observe that all the three machine learning estimates, (II), (III), and (IV), based on the polynomial basis yield larger point estimates than the low-dimensional estimates (I).
This may indicate a potential bias of the conventional estimator based on a low-dimensional polynomial approximation.
However, it is also worthy of remarking that polynomial bases (including the Legendre bases) are not suitable to approximating functions of variables that have unbounded supports.
Hermite bases, on the other hand, are capable of approximating functions of variables with unbounded supports.\footnote{\label{foot:approximation}Namely, the set of functions spanned by the standard polynomials is dense in the set $C^0(K)$ of continuous functions (or $L^2(K)$ of square integrable functions) defined only on a \textit{compact} support $K \subset \mathbb{R}$, whereas the set of functions spanned by the Hermite polynomials is dense in the set $L^1(\mathbb{R})$ of integrable functions defined on the \textit{entire} real line $\mathbb{R}$. As such, when the regressor(s) are infinitely supported, as is likely the case in the current application, the approximation by the standard polynomial basis is not credible but that by the Hermite polynomial basis is credible.}
We therefore focus on the results based on the Hermite basis.\footnote{\label{foot:not_invariant}In addition to this difference in the approximation theoretic properties between the polynomial and Hermite bases, we also remark that the difference may be due to the fact that the proposed method is not invariant to an invertible linear transformation of the regressors like the LASSO.}
For the unskilled labor coefficient, all the three machine learning methods, (II), (III), and (IV), still yield larger point estimates than that of the low-dimensional method (I), but the estimate based on our proposed method (IV) is relatively smaller and closer to that of (I).
For the skilled labor coefficient, the two machine learning methods, (II) and (III), yield slightly larger estimates than that of (I), while our proposed machine learning method (IV) yields a slightly smaller estimate than that of (I).
In summary, the results of the DML based on the LASSO or the random forest significantly differ from the result of the conventional low-dimensional method, but our proposed method also yields slightly different results from those of the DML based on the LASSO or the random forest, as is also the case in our simulation studies presented in Section \ref{sec:simulations}.
We ran several other estimates for robustness checks.
Appendix \ref{sec:additional_empirical} presents estimation results based on alternative values of the tuning parameters.
Appendix \ref{sec:additional_empirical} also presents results based on the method without cross fitting introduced in Section \ref{sec:without_cross_fitting}.
It turns out that the qualitative pattern of the results summarized above remain robust.
Finally, we conclude this section with a few remarks about the validity of the estimation approach employed in this empirical application following \citet{OlleyPakes1996} and \citet{levinsohn2003estimating}.
It is well known today that the estimating equation \eqref{eq:production_first_stage} fails to identify the parameter $\theta$ in general, as first pointed out by \citet{ackerberg2015identification}.
That said, they also suggest that $\theta$ can be correctly identified by \eqref{eq:production_first_stage} under certain DGPs.
They include DGPs with: (1) i.i.d. optimization error in $\ell_{it}$ and not in $m_{it}$; or (2) i.i.d. shocks to the price of labor or output after $m_{it}$; for instance \citep[][Sec. 3.1]{ackerberg2015identification}.
As such, we stress that the validity of the estimation method present above is contingent on these assumptions about the underlying DGPs.
\section{Summary and Discussions}\label{sec:summary_discussions}
In this paper, we propose a new method of inference in high-dimensional regression models.
The estimation procedure is based on a combined use of the OGA, HDAIC, and DML.
The method of inference about any low-dimensional subvector of high-dimensional parameters is based on a root-$N$ asymptotic normality, which does not require the exact sparsity condition or the $L^p$ sparsity condition.
In stead imposed are conditions on the rate at which the absolute size of parameters decays when descendingly ordered.
We demonstrate through simulation studies superior finite sample performance of this proposed method over those based on two popular alternatives, namely the LASSO and the random forest.
The extent of this outperformance is more prominent under less sparse models characterized by slower polynomial decays.
Finally, we illustrate an application of the method to production analysis using a panel of Chilean firms.
Using the tensor product of Hermite basis as high-dimensional controls, we find that estimates based on our proposed method differ from those based on the LASSO and random forest, similarly to what we observe in the simulation studies.
We close this paper with discussions of limitations, omitted extensions and potential directions for future research.
First, unlike regressions and like the LASSO, the method is not invariant to invertible linear transformation of the regressors $X$.
Practitioners should be aware of this drawback in our proposed method.
Second, as is the case with other DML methods, our proposed method based on DML is subject to random estimates.
To overcome this problem, we present an alternative procedure without cross fitting in Section \ref{sec:without_cross_fitting}, but this comes at the expense of an alternative set of assumptions.
Again, practitioners should be aware of these tradeoffs in choosing an appropriate method.
Third, we focus on high-dimensional linear regression models with exogenous regressors throughout the main text.
We provide an extension to high-dimensional linear IV models in Section \ref{sec:iv}.
Extensions to other important models are left for future research.
\vspace{1cm}