Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,950 characters · 11 sections · 70 citation commands
Inference in High-Dimensional Regression Models without the Exact or $L^p$ sparsity
\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]}
The advent of modern machine learning\footnote{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 temlyakov2000weak to order the high-dimensional regressors in a descending order of explanatory power. Second, use the high-dimensional Akaike information criterion 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. 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 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., belloni2012sparse, BCH2014RES, javanmard2014confidence, van2014asymptotically, zhang2014confidence, caner2018asymptotically, caner2018high, belloni2018high, galbraith2020simple, gold2020inference, kueck2021estimation to list a few. The majority of these papers focus on utilizing LASSO 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 buhlmann2011statistics, giraud2015introduction, and hastie2019statistical. In particular, we take advantage of the theoretical results of OGA+HDAIC by 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 temlyakov2000weak, tropp2004greed, tropp2007signal, ing2011stepwise, and share ties with other iterated model selection methods such as the least absolute angle regression of efron2004least, the $L_2$-boosting of buhlmann2003boosting, the test-based forward model selection of kozbur2017testing,kozbur2020analysis, and so forth. Third, this paper is related to the literature on Neyman orthogonal scores or locally robust scores, e.g., belloni2015uniform, chernozhukov2016locally, 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 shibata1980asymptotically,ing2007accumulated as well as the ordinary- and super-smoothness on probability density functions that are used in the deconvolution literature fan1991optimal,fan1993nonparametric.
Consider the linear regression model
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$:
In Section (ref), we consider an extended model in which we introduce approximation errors in (ref)--(ref).
To construct a moment restriction under (ref)--(ref), consider the orthogonal score function from Robinson:
where $X'\gamma_0 = E[Y|X]$ and $\eta=(\gamma,\beta)$. Note also that $X'\beta_0 = E[D|X]$ follows by construction from (ref).
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$.
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 (ref) and this cross-fitting method allows for the high-level theory of the double/debiased machine learning ccddhnr to be applicable. Section (ref) 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 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 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), we formally present theoretical arguments in support of this claim.
While Algorithm (ref) 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).
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). 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.
Assumption (ref) (a) requires a random sampling of data. Assumption (ref) (b) requires that the correct model is given by (ref) and (ref). Assumption (ref) (c)--(e) requires bounded moments of various variables. Assumption (ref) (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.
Assumption (ref) (ref) 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) (ref) 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) (ref) 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 candes2007dantzig, a restricted eigenvalue type condition for the sample Gram matrix bickel2009simultaneous is not required here -- see also the discussion in 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$.
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) 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
as well as the polynomial summability condition
following the discussion in ing2020. On the other hand, Condition (b) implies the conventional exponential decay condition that, for some $\alpha'>0$,
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 ing2020, we now present a concrete example where our Assumption (ref) holds but the sparsity does not. Suppose that
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) 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).$
A proof is provided in Appendix (ref). This theorem guarantees that the estimator $\check{\theta}$ of $\theta_0$ provided by Algorithm (ref) converges at the rate of $\sqrt{N}$ and asymptotically follows the normal distribution under Assumption (ref)--(ref). 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.
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 BCH2014RES in developing baseline data generating processes (DGPs). The linear regression model is specified by
where $\theta_0=0.5$ and $p = dim(X) = 500$. Consistently with this specification, data are generated by the system
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\}$.
Table (ref) 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). In particular, Appendix (ref) 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.
Extending the baseline model (ref)--(ref), consider the following partially linear model motivated by BCH2014RES:
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). This extension is useful in certain applications, such as the one we present in Section (ref). 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 (ref)--(ref):
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 belloni2012sparse,BCH2014RES.}
Recall the orthogonal score $\psi(Y,D,X;\theta,\eta)$ defined in (ref). With this orthogonal score, we propose to obtain $\check\theta$ and $\widehat\Omega$ via Algorithm (ref) presented in Section (ref) even under the current extended setting with approximation errors.
With the extended model (ref)--(ref), a different set of assumptions are imposed from those in the baseline model. First, we slightly modify Assumption (ref) as follows.
In part (ref), 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) from the baseline model. However, it should be stressed that we now impose Assumption (ref) on (ref)--(ref) rather than (ref)--(ref). 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$.
Assumption (ref) (ref) requires the fourth moment of the approximation error to be bounded, (ref) assumes the second moment to be of order $\log p/N$, and (ref) 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.
The following theorem establishes the asymptotic normality of $\check\theta$ along with the asymptotic validity of inference under the extended model with approximation errors.
A proof is presented in Appendix (ref). The same remarks as those presented below the statement of Theorem (ref) apply here.
{\color{black} We want to stress that Theorem (ref) is not an immediate consequence given Theorem (ref), because the original proofs for convergence rates of OGA+HDAIC in ing2020 do not permit approximately sparse models. In order to show Theorem (ref), we establish convergence rates for the OGA+HDAIC in approximately sparse regression models.}
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) to consider the partial linear model (ref)--(ref) entailing the reduced forms (ref)--(ref). Furthermore, we continue to use the same orthogonal score $\psi(Y,D,X;\theta,\eta)$ defined in (ref). However, we now replace Algorithm (ref) 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)'$.
To establish asymptotic properties for this new estimator $\widetilde\theta$, we continue to impose Assumptions (ref), (ref), and (ref). As in Section (ref), we focus on the more difficult case, namely the polynomial decay case, for brevity.
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.
A proof is given in Appendix (ref).
{\color{black}We stress that, although the proof builds on that of Theorem 1 in 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). Due to the lack of exact or approximate sparsity of $\beta_0$, this is shown via different strategies from those employed in 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) 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).
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 marschak1944random. While early studies of production functions address this simultaneity problem by explicitly modeling rational choice structures of firms, OlleyPakes1996 more recently propose a novel idea to use the inverse of the reduced-form investment choice function as a control function. levinsohn2003estimating propose to use intermediate input, instead of investment, as a control variable for a number of advantages.
The use of the control function a la levinsohn2003estimating entails the partial linear estimating equation for the labor elasticity of the form
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 OlleyPakes1996 and levinsohn2003estimating for details.
In light of the partial linear form (ref), OlleyPakes1996 and levinsohn2003estimating propose to use the estimator of 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 petrin2004production which implements estimation of (ref) via the parametric third-degree polynomial approximation
See petrin2004production.
To mitigate the approximation bias asymptotically, we consider a higher-dimensional approximation
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 (ref) 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 gallant1987semi,chen2007large for our basis $\phi$, and apply our proposed method to (ref) to get an estimate of $\theta$ and its standard error.
Following levinsohn2003estimating, we use a plant-level panel of Chilean firms from 1979 to 1986. See 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 levinsohn2003estimating, we consider the estimating equation of the form
as in (ref). 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:
as in (ref).
The term $\tau_{t}$ represents time effects. Following levinsohn2003estimating, we include the indicator for year groups 1979--1981, 1982--1983, and 1984--1986.
For estimation of (ref) 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 (ref) 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 ccddhnr -- see Appendix (ref) for details. See Appendix (ref) for details about the Hermite basis. We repeat analogous estimation procedures for (ref).
Table (ref) summarizes estimation results. Row (I) copies estimates from 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.
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{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 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 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{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).
We ran several other estimates for robustness checks. Appendix (ref) presents estimation results based on alternative values of the tuning parameters. Appendix (ref) also presents results based on the method without cross fitting introduced in Section (ref). 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 OlleyPakes1996 and levinsohn2003estimating. It is well known today that the estimating equation (ref) fails to identify the parameter $\theta$ in general, as first pointed out by ackerberg2015identification. That said, they also suggest that $\theta$ can be correctly identified by (ref) 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 ackerberg2015identification. As such, we stress that the validity of the estimation method present above is contingent on these assumptions about the underlying DGPs.
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), 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). Extensions to other important models are left for future research.