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.
110,355 characters · 24 sections · 62 citation commands
Learning non-smooth models: instrumental variable quantile regressions and related problems
\global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long \global\long
\setcounter{page}{1}
The linear instrumental variables quantile model (IVQR) formulated by chernozhukov2005iv,chernozhukov2006instrumental,Chernozhukov2008 has found wide applications in economics. The basic moment condition can be written as follows
where $\tau\in(0,1)$ is the quantile of interest, $Y_{i}\in\mathbb{R}$, $X_{i}\in\mathbb{R}^{p}$ and $Z_{i}\in\mathbb{R}^{L}$ are i.i.d observed variables and $\beta_{*}\in\mathbb{R}^{p}$ is the unknown model parameter. Assume that $p$ and $L$ are fixed with $L\geq p$. The typical setup is that only one (or few) component of $X_{i}$ is endogenous and other components of $X_{i}$ are contained in $Z_{i}$. In the policy evaluation setting, the variable denoting the status of treatment is usually considered endogenous. If this variable only enters the regression equation as one endogenous regressor, then we can apply the existing methods (e.g., the popular method by chernozhukov2006instrumental) for estimating the treatment effect. When multiple endogenous regressors enter the regression equation, it imposes enormous (or even prohibitive) computational challenges to common estimation strategies, which typically involve solving nonconvex and non-smooth optimization problems.
In empirical studies that investigate the causal effects of an endogenous treatment, multiple endogenous regressors arise naturally. For example, the interactions between the treatment variable and other variables are often included in the regression to study the heterogeneity of the treatment effects. In the IVQR setting, this leads to multiple endogenous variables in the regression equation. Consider the randomized training experiment conducted under the Job Training Partnership Act (JTPA). JTPA training services are randomly offered to people, who can then choose whether to participate in the program. One key policy question is whether this program has an effect on earnings. Of course, the baseline question is whether the program has a positive effect overall. In addition, one might ask questions such as whether the effect of the program differs by participants' race, age, etc. These questions can be answered in the regression setting by learning the coefficients for interactions between the treatment status and other variables denoting race, age, etc.
We show that the efficient computations for estimating IVQR from a purely algorithmic perspective might not exist. In particular, we prove that optimization of the usual GMM or similar formulation of the IVQR problem is NP (non-deterministic polynomial-time) hard. To the best of our knowledge, this is the first formal statement on the NP-hardness of the IVQR problem. Due to this result, it is unlikely to find a computationally efficient algorithm that minimizes the GMM criterion function. Notice that this is different from some existing methods that cast the GMM problem as NP-hard problems, e.g., mixed integer program or quadratic program with complementarity constraints burer2012milp,chen2017exact. These papers only imply that these NP-hard problems are more difficult than IVQR GMM-estimation; in other words, algorithms that solve certain NP-hard problems can be used for IVQR estimation. However, the question of whether there exist simple algorithms for IVQR GMM-estimation is still open. In this paper, we provide a negative answer by showing that IVQR estimation is at least as difficult as the set partition problem, which is known to be NP-hard since karp1972reducibility. Moreover, we show that it is also very difficult to find an approximate solution in that obtaining a solution within a constant factor from the global solution is also NP-hard. Therefore, there is inherent complexity that is unlikely to be fully handled by convex relaxation, smoothing or other pure computational remedies for handling non-convex and non-smooth optimization problems.
The NP-hardness suggests that serious efforts are needed to address the computational challenge, but it does not necessarily mean that the model has to be shunned by practitioners. In fact, some of the most popular methods, such as decision trees and $k$-means clustering laurent1976constructing,MAHAJAN201213, are NP-hard. One key observation is that NP-hardness is a statement on the difficulty of the worst-case situation, but the data that arise in practice might not be these “bad” situations. Although optimization algorithms derived from a pure computational point of view are not guaranteed to exploit “probabilistic” properties of the data, there are many inspiring precedents in machine learning that exploit statistical properties to overcome computational difficulties. For example, finding a hyperplane to classify binary targets is a common approach. Whereas using the misclassification rate as the loss function leads to an NP-hard optimization problem, one can use a convex loss function (e.g., hinge loss or logistic loss) and still obtain nice statistical properties, see bishop2006pattern,shalev2014understanding.
Here, we propose to apply the same principle to general non-smooth models. Instead of focusing only on clever computational tricks for the usual GMM formulation, we exploit the statistical aspect of the model and provide an alternative estimation and inference strategy that admits fast computation algorithms. Although estimates from any feasible optimization programs might be of questionable quality, our plan is to use a computationally cheap method and turn a dubious estimate into one with nice statistical properties. Of course, this final estimate might not be the global solution of any optimization problem but good statistical properties would often suffice for estimation and inference purposes.
Our proposal is a $k$-step framework that exploits the identification strength of a general (possibly nonsmooth) GMM model.\footnote{In a sense, many methods for high-dimensional sparse estimation can be also viewed as exploiting identification in order to avoid computational difficulties. Due to the NP-hardness of computing “$\ell_{0}$-regularized” estimators chen2014complexity, many $\ell_{1}$-regularized estimators have been considered and, under strong identification (in the form of sparse eigenvalue conditions or similar conditions), have been shown to achieve optimal statistical properties candes2007dantzig,RaskuttiWainwrightYu2011,cai2017confidence.} This framework starts with an initial estimate and updates it iteratively, where each iteration involves only matrix multiplications and no optimization at all. We show that when the identification is strong, only few iterations are needed to obtain an estimator that is asymptotically equivalent to the GMM estimator.
One key finding is that the initial estimator does not need to be consistent. This is a new result as far as we know. We show that there is a neighborhood of the true parameter value whose area is determined by the identification strength and any point in this neighborhood can be used as an initial estimate; when the identification is strong in the usual sense for GMM models (i.e., Jacobian of the population moment has well-behaved singular values), this neighborhood is fixed and its area does not shrink to zero with sample size $n\rightarrow\infty$. Allowing for inconsistent initial estimates provides valuable convenience in practice. When the initial estimates are computed via complicated optimization problems, computationally feasible solutions might not be consistent since they are often not global solutions even under strong identification. Moreover, even if initial estimates are easily obtained via efficient algorithms, these algorithms might not be feasible when the sample size becomes extremely large. For large datasets that arise in modern research, even common convex problems often encounter scalability issues, see e.g., yang2013quantile,lin2017distributed. Since the initial estimate does not need to be even consistent, one can simply compute it using a very small fraction of the entire sample. Efficient algorithms are available for $k$-step iterations as they only require matrix multiplications. This leads to a fast procedure for handling massive datasets.
We also provide methods and theory for two important aspects of implementing the $k$-step framework. First, we provide a tuning-free estimate for the Jacobian matrix of the moment condition. This estimate is needed for the $k$-step correction and can be tricky to obtain in practice. The formula defining the Jacobian typically contains the (conditional) density of a variable and one typical estimation method is to use kernel with a bandwidth choice. However, the bandwidth might be quite difficult to choose in practice. We propose to generalize a fast bootstrap scheme for density estimation. The basic idea is as follows. Consider the estimation of the density of a variable at a point. Since the asymptotic variance of the bootstrapped sample quantile of this variable is directly related to the density, bootstrapping sample quantile can be used to obtain a consistent density estimate. We generalize this approach to estimation of Jacobian of general nonsmooth GMM models. Notice that the Jacobian matrix can be estimated entry by entry. In the case of IVQR, this has a closed-end solution with computation burden similar to bootstrapping sample quantiles; in our experiments, with a sample size $n=5000$, computing 2000 bootstrap samples takes 0.5 second. Since the entire procedure can be used for general GMM settings and does not require choices such as bandwidth, we believe it is of independent interest. For example, one can use it as a tuning-free density estimator, which can be used to compute the optimal bandwidth for a second-stage estimation. We establish the consistency for this estimate.
Second, we present a formulation that yields initial estimators for IV quantile models and related problems via mixed integer linear programs (MILP). Although MILP is also an NP-hard problem, it is one of the most well-studied and well-understood hard problems and so much progress has been made in recent years that many believe large-scale problems are now feasible, see e.g., bertsimas2005optimization,bixby2007progress,junger200950,linderoth2010milp. As pointed out by bertsimas2016best, the speed of finding global solutions for mixed integer optimization improved approximately 450 billion times between 1994 and 2015. In our experiments, we deliver good estimates for coefficients of 20 endogenous variables within 5 seconds. The high-dimensional version can handle regression equations with 500 endogenous variables within minutes. Since we do not need the initial estimate to be consistent, one can terminate MILP optimization before a global solution is found. We also provide an early-termination rule and establish its theoretical validity. The MILP constructions are not unique to low-dimensional IV quantile regressions. We outline how MILP can be used for related problems, including high-dimensional IV quantile regressions, censored regressions and censored IV quantile regressions.
This paper is related to the literature of analyzing computational complexity of statistical learning methods. One strand of this literature is based on the NP-completeness theory dating back to at least cook1971complexity,karp1972reducibility,levin1973universal. The classical textbook by johnson1979computers documented hundreds of NP-complete problems. The efficient solution of any of these problems still remains elusive today and is now widely considered impossible. Since NP-hard problems are at least as hard as these difficult problems, finding out whether a statistical procedure is computationally NP-hard has been an important task for analyzing the implementability, see e.g., laurent1976constructing,berthet2013optimal,chen2014complexity,chen2017strong. Another recent active area of research is to find efficient approximation algorithms for NP-hard problems and characterize theoretical possibilities of such approximations, see arora2009computational,williamson2011design for excellent textbook treatments and surveys.
An early version of this paper was inspired by the fascinating literature of applying mixed integer programming to statistical learning. Recent progress has drastically improved the speed of mixed integer optimizations, which are now considered a feasible tool for some high-dimensional problems. Most of the advancement concerns high-dimensional linear models; see bertsimas2014least,liu2016global,bertsimas2016best,mazumder2017discrete. The main argument for considering these nonconvex algorithms is that they, compared to convex regularized methods, enjoy more desirable statistical properties. zubizarreta2012using proposed using mixed integer programming for matching estimators in causal inference.
Our work contributes to the fast growing literature of IV quantile regression. The IV quantile regression extends the advantage of quantile regression (Koenker1978) to the settings with endogenous regressors. The conceptual framework and identification of the IV quantile models has been studied by abadie2002instrumental, chernozhukov2005iv and imbens2009identification; see wuthrich2014comparison,wuthrich2019closed, melly2016local and chernozhukov2017instrumental for more discussions. Semiparametric and nonparametric specifications have been studied in horowitz2007nonparametric,chernozhukov2007instrumental,chen2009efficient,chen2012estimation,gagliardini2012nonparametric. The GMM estimation approach applies the classical GMM method for the moment condition in ((ref)). The computational burden of minimizing a nonconvex and non-smooth objective function is known to be challenging for larger dimensional models. The quasi-Bayesian approach of chernozhukov2003mcmc has been suggested, but could be difficult to tune it to sufficiently explore the entire parameter space. In an interesting paper, chen2017exact proposed formulating the original GMM problem as a mixed integer quadratic program (MIQP). Smoothing the GMM objective function has also been considered by kaplan2017smoothed and deCastroGalvaoKaplanLiu2018. The so-called inverse quantile regression by chernozhukov2006instrumental,Chernozhukov2008 takes a different route and reduces the dimension of the space over which the optimization is needed. lee2007endogeneity considers a control function approach but deviates from the model ((ref)). 1kaido2018 propose a decomposition of the IVQR problem into a set of convex sub-problems. pouliot2019 develops mixed integer programs that avoid nonparametric density estimation and allow for weak identification.
Our work is also related to the $k$-step estimator in the econometrics and statistics literature. The classical references include robinson1988stochastic and andrews2002equivalence. The main difference in assumption is that our results do not assume that the sample version of the moment condition is differentiable. We provide a general theory in this setting, which might be of independent interest. Moreover, we show that a consistent starting point is not necessary.
We will use $E_{n}$ to denote the sample average $n^{-1}\sum_{i=1}^{n}$. The $\ell_{q}$-norm of a vector will be denoted by $\|\cdot\|_{q}$ for $q\geq1$; $\|\cdot\|_{\infty}$ denotes the maximum absolute value of a vector, i.e., the $\ell_{\infty}$-norm. Hence, $\|\cdot\|_{2}$ denotes the Euclidean norm. We use $\|\cdot\|$ to denote the spectral norm of a matrix. The indicator function is denoted by $\mathbf{1}\{\}$. For any positive integer $r$, we use $\mathbf{1}_{r}$ to denote the $r$-dimensional vector of ones and $I_{r}$ to denote the $r\times r$ identity matrix. We use $\lambda_{\max}(\cdot)$ and $\lambda_{\min}(\cdot)$ to denote the maximal and the minimal eigenvalues of symmetric matrices. We use $\log$ to denote the natural logarithm. The rest of the paper is organized as follows. Section (ref) provides a general theory of $k$-step correction for non-smooth problems and outlines the details of implementation for IVQR. Section (ref) presents the MILP formulation of IVQR and related problems; we also derive theroetical results for early termination of the algorithm without needing to find a global solution. Section (ref) presnets a new tuning-free methodology for Jacobian estimation, which is of independent interest. Monte Carlo simulations are presented in Section (ref). Section (ref) considers the JTPA example. The proofs of theoretical results are in the appendix.
We start by considering the computational complexity of the GMM formulation of IVQR. The moment condition for ((ref)) is \[ EZ_{i}\left(\mathbf{1}\{Y_{i}-X_{i}'\beta\}-\tau\right)=0. \]
Of course, we can replace $Z_{i}$ with transformations of $Z_{i}$; here, we assume that $Z_{i}$ is already the desired transformation of the instruments. Then the GMM estimator is \[ \hat{\beta}_{GMM}=\underset{\beta\in\mathbb{R}^{p}}{\arg\min}\left\Vert \hat{\Omega}^{1/2}\sum_{i=1}^{n}Z_{i}\left(\mathbf{1}\{Y_{i}-X_{i}'\beta\}-\tau\right)\right\Vert _{2}, \] where $\hat{\Omega}$ is a weighting matrix, either estimated from the data or pre-determined. For simplicity, we use the identity $\hat{\Omega}=I_{r}$. Then the computation can be summarized as follows.
This seemingly routine GMM estimation turns out to be NP-hard. The characterization of NP-hardness is an important assessment on the computational complexity of a problem. Loosely speaking, the class of NP is the set of problems for which a potential answer can be verified efficiently (in polynomial time). Consider the set partition problem: given $\{a_{i}\}_{i=1}^{m}$ integers, decide whether there exists a subset $B\subset\{1,....m\}$ such that $\sum_{i\in B}a_{i}=\sum_{i\notin B}a_{i}$. For any given set $B$, we can easily verify whether $\sum_{i\in B}a_{i}=\sum_{i\notin B}a_{i}$ holds; as a result, the set partition problem is in NP. However, the question of whether such a set $B$ exists is a notoriously difficult problem. In fact, the set partition problem belongs to the set of the hardest NP problems, the so-called NP-completeness class. All the problems in the NP-completeness class are equivalent to each other and solving any NP-complete problem would also solve every NP problem. The NP-completeness class contains many problems that are considered difficult, such as Boolean satisfiability problem, traveling salesman problem, set partition problem, vertex cover problem, etc. No algorithms that can solve an NP-complete problem in polynomial time have been found and whether such an algorithm exists remains an open question, one of the most fundamental questions in computer science and modern mathematics. A detailed treatment on NP-hardness can be found in the classical textbook johnson1979computers or chapter 7 of sipser2012introduction.
The proof of Theorem (ref) establishes that any algorithm solving Problem (ref) can also be used to solve the set partition problem, which is an NP-complete problem. From the proof, one can see that the difficulty arises from the dimension $p$ since the problem is already NP-hard when both $n$ and $L$ are equal to $p$. Since there are still no fast algorithms for any NP-complete problem, we do not expect to magically find a fast algorithm for Problem (ref) and thereby solve all the NP-complete problems.
Since finding the exact solution for NP-hard problems is difficult (if not impossible), one important question is whether it is possible to find a good approximation. To be specific, we introduce the following definition adapted from chapter 3 of ausiello2012complexity.
In terms of the above criteria, there has been tremendous success in finding efficient approximation algorithms for many problems, but satisfactory enough approximation might not be possible. For example, the $k$-center problem admits a simple $2$-approximation and the traveling sales problem has a fast $3/2$-approximation; however, it is NP-hard to find a $\rho$-approximation of the former for any $\rho<2$ and of the latter for any $\rho<220/219$, see williamson2011design. Unfortunately, there are also problems, such as the set covering problem, for which an approximation solution is as hard as the exact solution since a $\rho$-approximation is NP-hard for any $\rho>1$, see lund1994hardness. We now show that the GMM formulation of IVQR belongs to this class of very difficult optimization problems.
The proof of Theorem (ref) is built on the idea that any algorithm solving the GMM formulation in Problem (ref) can be used to solve the so-called minimum unsatisfiability problem of linear systems. The latter problem is known to have no polynomial-time $\rho$-approximation for any $\rho>1$. In the machine learning literature, it is sometimes referred to as the half-space learning problem or linear perceptron learning problem, which finds the best hyperplane to classify binary outcomes. A direct consequence of the non-approximability of this problem is that minimizing the misclassification rate is not a computationally efficient way of training classifiers, see e.g., ben2001efficient; for this reason, convex functions, such as the hinge loss or logistic loss, often serve as the objective function in these learning problems.
In light of Theorems (ref) and (ref), it appears unrealistic to guarantee adequate econometric properties by relying exclusively on algorithms that attempt to accurately solve the GMM formulation. Instead, we focus on designing a procedure that is based on the statistical properties and is as computationally cheap as possible.
We now present our proposal and its theoretical justification. Let $\{W_{i}\}_{i=1}^{n}$ be i.i.d observations. Let $G(\beta)=Eg(W_{i};\beta)$ be an GMM model, where $g$ is an $\mathbb{R}^{L}$-valued function that is possibly non-smooth in $\beta\in\mathcal{B}$. The true parameter value $\beta_{*}$ is assumed to be uniquely defined by $G(\beta_{*})=0$. Let $\Gamma_{*}=(\partial G(\beta)/\partial\beta)(\beta_{*})\in\mathbb{R}^{L\times p}$, $G_{n}(\beta)=n^{-1}\sum_{i=1}^{n}g(W_{i};\beta)$ and $H_{n}(\beta)=\sqrt{n}(G_{n}(\beta)-G(\beta))$.
We view our proposal as performing two tasks: $\sqrt{n}$-estimation and inference by asymptotic normality.
Suppose that we have an initial estimator $\bar{\beta}$ for $\beta_{*}$ and an estimator $\hat{\Gamma}$ for $\Gamma_{*}$. Computing the initial inputs $\bar{\beta}$ and $\hat{\Gamma}$ will be addressed in Sections (ref) and (ref), respectively. Neither $\bar{\beta}$ nor $\hat{\Gamma}$ is assumed to be consistent. Now consider the following one-step correction estimator. We define the one-step correction operator by
We also define the sequence $\{\mathcal{A}_{k}(v,Q)\}_{k=1}^{\infty}$ recursively by \[ \mathcal{A}_{k+1}(v,Q)=\mathcal{A}(\mathcal{A}_{k}(v,Q),Q), \] where $\mathcal{A}_{1}(v,Q)=\mathcal{A}(v,Q)$. The following result allows us to compare the estimation errors of $\bar{\beta}$ and its one-step correction.
Lemma (ref) depicts the basic intuition that underlies the $k$-step estimator. Suppose that \[ \rho:=c\|(\hat{\Gamma}'\hat{\Gamma})^{-1}\hat{\Gamma}'\|\sup_{\beta\in\mathcal{B}_{0}}\|\beta-\beta_{*}\|_{2}+\|(\hat{\Gamma}'\hat{\Gamma})^{-1}\hat{\Gamma}\|\cdot\|\hat{\Gamma}-\Gamma_{*}\|<1. \]
Then Lemma (ref) implies that \[ \|\mathcal{A}(\beta,\hat{\Gamma})-\beta_{*}\|_{2}\leq\rho\|\beta-\beta_{*}\|_{2}+n^{-1/2}\|(\hat{\Gamma}'\hat{\Gamma})^{-1}\hat{\Gamma}'\|\sup_{v\in\mathcal{B}_{0}}\|H_{n}(v)\|_{2}. \]
Define $T_{*}=n^{-1/2}\|(\hat{\Gamma}'\hat{\Gamma})^{-1}\hat{\Gamma}'\|\sup_{v\in\mathcal{B}_{0}}\|H_{n}(v)\|_{2}/[2(1-\rho)]$. If $\|\beta-\beta_{*}\|_{2}\geq T_{*}$, then we have $\|\mathcal{A}(\beta,\hat{\Gamma})-\beta_{*}\|_{2}\leq\bar{\rho}\|\beta-\beta_{*}\|_{2}$, where $\bar{\rho}=(\rho+1)/2<1$. If $\|\mathcal{A}(\beta,\hat{\Gamma})-\beta_{*}\|_{2}\geq T_{*}$, then $\|\mathcal{A}_{2}(\beta,\hat{\Gamma})-\beta_{*}\|_{2}\leq\bar{\rho}^{2}\|\beta-\beta_{*}\|_{2}$. As we iterate, the estimation error shrinks exponentially (due to $\bar{\rho}<1$) until it becomes smaller than $\text{threshold}_{*}$. Typically, we can establish $T_{*}=O_{P}(n^{-1/2})$. This means that $k$-step iteration would yield an $\sqrt{n}$-consistent estimator very quickly.
By induction, we can invoke Lemma (ref) and obtain the following result on the rates of convergence for $\mathcal{A}_{K}(\beta,\hat{\Gamma})$. For now, we do not consider the randomness yet so the following result serves as a formalization of the intuition and a finite-sample result for establishing the final statistical properties.
Theorem (ref) has two important implications. First, the starting point $(\bar{\beta},\hat{\Gamma})$ does not need to be a consistent estimator for $(\beta_{*},\Gamma_{*})$. By Theorem (ref) , whenever we start from a small enough neighborhood (i.e., small enough $c_{1},c_{2}$), $\|\mathcal{A}_{K}(\bar{\beta},\hat{\Gamma})-\beta_{*}\|_{2}$ decays exponentially with $K$ until it reaches the parametric rate $n^{-1/2}$. The only requirement is that $c_{5}\leq c_{1}c_{3}^{1/2}(1-\rho_{*})\sqrt{n}$. A sufficient condition is $c_{1}=c_{3}^{1/2}/(2c_{4})$, $c_{2}=c_{3}^{1/2}/2$ and $n\geq16c_{3}^{-2}c_{4}^{2}c_{5}^{2}$. Notice that this requirement on $c_{1}$ and $c_{2}$ depends on $c_{3}$, which measures the identification strength. Since $c_{3},c_{4},c_{5}$ are bounded away from zero and infinity under strong identification, we allow $c_{1}$ and $c_{2}$ to be bounded away from zero and only require $n$ to be large enough (instead of tending to infinity).
Second, the parametric rate $n^{-1/2}$ is guaranteed after $O(\log n)$ iterations. Notice that $\rho_{*}<1$. This means that $K\geq(\log n)/(2\log\rho_{*}^{-1})$, we have that \[ \|\mathcal{A}_{K}(\bar{\beta},\hat{\Gamma})-\beta_{*}\|_{2}\leq n^{-1/2}\left(c_{1}+\frac{c_{3}^{-1/2}c_{5}}{1-\rho_{*}}\right). \] This is computationally quite attractive. Even if the starting point is not consistent or its rate of convergence can be arbitrarily slow, we only need a few iterations to obtain a $\sqrt{n}$-consistent estimator. Since there is no optimization in each iteration, this can be done extremely fast. In fact, this is the key property we shall exploit when dealing with massive samples. Now we state the result under commonly imposed regularity conditions.
We have the following result.
By Corollary (ref), if we have strong identification, bounded Hessian for $G(\cdot)$ and the empirical process $H_{n}(\cdot)$ is a Donsker class, then after $2\log n$ iterations, we will obtain a $\sqrt{n}$-consistent estimator as long as $\hat{\Gamma}$ and $\bar{\beta}$ lie in a small enough but fixed neighborhood of the true parameters with high probability. Once we obtain a $\sqrt{n}$-consistent estimator for $\beta_{*}$, we can use it to construct a consistent estimator for $\Gamma_{*}$; it turns out that the consistency of $\Gamma_{*}$ is needed to obtain asymptotic normality.
We now derive the asymptotic normality for the $k$-step estimator. We also address an important robustness issue. Obviously, the output of the $k$-step iteration depends on the number of iterations $K$ and the initial estimators $\bar{\beta}$ and $\hat{\Gamma}$. We now explicitly express such dependence and address the issue of sensitivity with respect to $(K,\bar{\beta},\hat{\Gamma})$.
Theorem (ref) provides the main tool for inference. It says that as long as $\hat{\beta}$ and $\hat{\Gamma}$ are consistent (easily achieved by first iterating $2\log n$ times from the initial estimator as shown in Section (ref)), we have $\|\mathcal{A}_{K}(\hat{\beta},\hat{\Gamma})-\beta_{*}+n^{-1/2}(\Gamma_{*}'\Gamma_{*})^{-1}\Gamma_{*}'H_{n}(\beta_{*})\|_{2}=o_{P}(n^{-1/2})$ for any $K\geq1+2\log n$. Commonly imposed regularity conditions would require that $H_{n}(\beta_{*})\rightarrow^{d}N(0,\Omega_{*})$ for some matrix $\Omega_{*}\in\mathbb{R}^{L\times L}$. Hence, we obtain \[ \sqrt{n}(\mathcal{A}_{K}(\hat{\beta},\hat{\Gamma})-\beta_{*})\rightarrow^{d}N(0,(\Gamma_{*}'\Gamma_{*})^{-1}\Gamma_{*}'\Omega_{*}\Gamma_{*}(\Gamma_{*}'\Gamma_{*})^{-1}). \]
Moreover, Theorem (ref) also provides a robustness guarantee on the asymptotic approximation. Since we are taking a supreme in ((ref)), the approximation of $\mathcal{A}_{K}(\beta,\Gamma)-\beta_{*}$ by $-n^{-1/2}(\Gamma_{*}'\Gamma_{*})^{-1}\Gamma_{*}'H_{n}(\beta_{*})$ holds uniformly in $(K,\beta,\Gamma)$. This means this approximation is robust to choices of $(K,\beta,\Gamma)$. For example, one can run the $k$-step, update the initial estimates, use the output to run the $k$-step iterations again, update the initial estimates, and repeat these steps arbitrarily many times. Theorem (ref) says that by doing so, one should not expect to change the inference results. We now summarize the entire procedure for estimation and inference in Algorithm (ref).
In many empirical applications, the sample size $n$ can be enormous. Notice that in Algorithm (ref), all the steps require only matrix multiplication, except Step 1 and potentially Step 3. In obtaining the initial estimation in Step 1, the MILP formulation presented in Section (ref) would be computationally very costly when the sample size $n$ exceeds 1000. However, since we do not require the consistency of $\bar{\beta}$ in Step 1, one can simply run the MILP on a randomly selected subsample of size $m$, where $m\ll n$. In simulations, we find that for $n=5\times10^{6}$, using $m=500$ yields decent performance. In this case, we only use $m/n=0.01\%$ of the data for initial estimation and Algorithm (ref) takes less than 15 seconds! This is a massive reduction in computing time because even linear programs can be slow in such massive scale. Hence, Algorithm (ref) can be used for large-scale quantile regressions.
The $k$-step correction in Section (ref) requires an estimate for the Jacobian of the population moment condition. We now provide a tuning-free option. To fix ideas, we recall $G(\beta)=Eg(W_{i};\beta)$. Now we are interested in estimating the Jacobian \[ \Gamma_{0}:=\Gamma(b_{0}), \] where $\Gamma(\beta)=\frac{\partial G(\beta)}{\partial\beta'}\in\mathbb{R}^{L\times p}$ and $b_{0}$ is an observed quantity, either random (e.g., GMM estimator or any initial estimate) or a deterministic quantity.
In these examples, estimation of $\Gamma(\beta)$ is typically done by a kernel method that requires a choice of bandwidth. Unfortunately in many applications, the choice of bandwidth can be quite difficult to determine and such a choice can often significantly affect the final estimation and inference. Here, we propose a tuning-free estimation scheme. Since the need of Jacobian estimation also arises in other situations (e.g., estimating asymptotic variance), we think this proposal is of independent interest. Moreover, for both examples, our proposal can be easily impolemented by essentially closed-end formulas.
At a high level, the mechanism that delivers the tuning-free property is closely related to the bootstrapped standard errors for quantile regressions. There are two popular methods for computing the standard errors for quantile regressions. One is to use the explicit formula and replace unknown density in this formula with a kernel estimate, which requires a bandwidth choice. The other is to bootstrap the quantile estimate. Notice that the latter avoids explicitly choosing a tuning parameter because the bootstrapped estimates naturally generate perturbations in a local neighborhood that allow us to learn the slope; see Section (ref) for more discussions. Here, we develop a computationally simple scheme implementing a similar mechanism for estimating the Jacobian for general GMM models. Our proposed estimator for $\Gamma(\beta)$ is based on the following two observations:
Let $\{\xi_{i}\}_{i=1}^{n}$ be random variables simulated independent of the data with $E(\xi_{i})=1$, say i.i.d $N(1,1)$ or binary variables in $\{0,2\}$ with equal probability. Recall $G_{n}(\beta)=n^{-1}\sum_{i=1}^{n}g(W_{i};\beta)$ and $H_{n}(\beta)=\sqrt{n}(G_{n}(\beta)-G(\beta))$. We define $G_{n}^{*}(\beta)=n^{-1}\sum_{i=1}^{n}\xi_{i}g(W_{i};\beta)$ and $H_{n}^{*}(\beta)=n^{-1/2}\sum_{i=1}^{n}(\xi_{i}-1)g(W_{i};\beta)$. We have some flexibility in generating the multipliers; for example, one can also simulate $(\xi_{1},...,\xi_{n})$ from a multinomial distribution with parameter $n$ and probabilities $(n^{-1},...,n^{-1})$, leading to the empirical bootstrap.
To estimate $\Gamma(b_{0})$, we solve the one-dimensional problem
where $\bar{C}>0$ is a large enough constant. If there are multiple solutions, we choose $b_{*}$ to be the solution closest to $b_{0}$. The following result provides the basic intuition.
Consider the case in which $G(\cdot)$ has a bounded second-order derivative, $H_{n}(\cdot)$ is Donsker and $b_{*}-b_{0}=o_{P}(1)$. Lemma (ref) essentially says \[ -n^{-1/2}H_{n}^{*}(b_{*})\approx\Gamma_{0}\times(b_{*}-b_{0}). \]
Since both $-n^{-1/2}H_{n}^{*}(b_{*})$ and $b_{*}-b_{0}$ are directly observed, we can simulate enough of them and regress $-n^{-1/2}H_{n}^{*}(b_{*})$ on $b_{*}-b_{0}$ via OLS without intercept. Formally, we use simulation to approximate
where $\mathcal{F}$ is the $\sigma$-algebra generated by the data. We now give the formal result.
In Theorem (ref), we allow $b_{0}$ to be random, which means $\Gamma_{0}=\Gamma(b_{0})$ can be random; when $b_{0}$ is non-random, $P(\rho_{2}<|\Gamma_{0}|<\rho_{3})$ is either one or zero. The assumptions on $H_{n}(\cdot)$ are satisfied when $H_{n}(\cdot)$ is Donsker. The conditions on $H_{n}^{*}(\cdot)$ can be easily verified when $\{\xi_{i}\}_{i=1}^{n}$ are sub-Gaussian multipliers. In this case, $H_{n}^{*}(\cdot)$ conditional on $\mathcal{F}$ is a sub-Gaussian process and one can invoke the usual maximal inequalities and tail bounds in van1996weak. Under these weak regularity conditions, Theorem (ref) establishes the consistency of $\hat{\Gamma}_{0}$. We now present simple computational methods for implementing this estimation strategy.
In Example (ref), the problem in ((ref)) becomes \[ \text{find}\ b_{*}\qquad s.t.\qquad\sum_{i=1}^{n}\mathbf{1}\{Y_{i}\leq X_{i}b_{*}\}Z_{i}\xi_{i}=d, \] where $d$ is a known number $d=\sum_{i=1}^{n}Z_{i}\left(\mathbf{1}\{Y_{i}\leq X_{i}b_{0}\}+(\xi_{i}-1)\tau\right)$. Since $X_{i}$ and $Z_{i}$ are scalars, we can easily obtain a closed-end solution. We summarize the idea below.
By Lemma (ref), we just need to use the cumulative sum of $\tilde{w}_{i}$. We summarize the details in Algorithm (ref). To solve ((ref)), we can simply apply Algorithm (ref) with $a=\sum_{i=1}^{n}Z_{i}\left(\mathbf{1}\{Y_{i}\leq X_{i}\beta_{0}\}+(\xi_{i}-1)\tau\right)$ and obtain a solution that satisfies ((ref)) with error $O_{P}(n^{-1}\max_{1\leq i\leq n}|Z_{i}\xi_{i}|)$.
One can also view Algorithm (ref) as a smart grid since we are essentially saying that the solution has to be one of $\{\tilde{y}_{i}\}_{i=1}^{n}$. Therefore, the number of grid points and the location of the grid points are completely determined by the data. Since the algorithm essentially only consists of sorting variables, it can be implemented efficiently. In our experiment, for a sample size $n=5102$, bootstrapping 2000 samples takes only 0.5 second.
We can use the method discussed above (with $X_{i}=Z_{i}=1$) for the computation of Example (ref). As mentioned before, $\{\xi_{i}\}_{i=1}^{n}$ can be generated from a multinomial distribution and in this case the computation procedure of $b_{*}$ described in Algorithm (ref) becomes computing the empirical quantile in bootstrap samples. Therefore, our proposed estimate $\hat{\Gamma}_{0}$ has a natural connection with bootstrap standard errors for quantile regression.
We would like to point out that bootstrapping quantile regressions in this case can be viewed as an alternative application of Lemma (ref). In Example (ref), we have$\Gamma_{0}=f_{Y}(b_{0})$, where $f_{Y}(\cdot)$ is the probability density function of the scalar variable $Y_{i}$ and $b_{0}$ is the $\tau$-quantile of the sample. Another way of writing ((ref)) is $H_{n}^{*}(b_{*})\approx-\Gamma_{0}\sqrt{n}(b_{*}-b_{0})$. The proposed $\hat{\Gamma}_{0}$ regresses $H_{n}^{*}(b_{*})$ onto $\sqrt{n}(b_{*}-b_{0})$ without intercept. Alternatively, since $\Gamma_{0}$ is known to be positive, we can simply estimate $\Gamma_{0}$ using $\sqrt{E[(H_{n}^{*}(b_{*}))^{2}\mid\mathcal{F}]/E[n(b_{*}-b_{0})^{2}\mid\mathcal{F}]}.$ Since $H_{n}^{*}(\cdot)$ is approximating a Brownian bridge, we know that $E[(H_{n}^{*}(b_{*}))^{2}\mid\mathcal{F}]\approx\tau(1-\tau)$ and simplify this formula as \[ \sqrt{\frac{\tau(1-\tau)}{E[n(b_{*}-b_{0})^{2}\mid\mathcal{F}]}}. \]
Notice that this is exactly the density estimate that we implicitly use when we use the bootstrap standard error. The asymptotic variance formula for the sample quantile is $\tau(1-\tau)/\hat{f}_{Y}^{2}(b_{0})$ for some density estimate $\hat{f}_{Y}(b_{0})$. If we equate this with the bootstrap variance $E[n(b_{*}-b_{0})^{2}\mid\mathcal{F}]$, we obtain that the implicitly used $\hat{f}_{Y}(b_{0})$ is given by the above formula. Therefore, bootstrapping sample quantiles can serve as a tuning-free density estimator. Our estimator $\hat{\Gamma}_{0}$, which is an OLS estimator, extends this strategy to Jacobian estimation for general GMM models.
The methodology we have presented so far does not assume a particular choice of the initial estimate. Therefore, one can choose from methods in the existing literature. In this section, we provide an MILP approach, which can be used for IVQR and related problems. We argue that this is an attractive alternative since it capitalizes on recent advances in mixed integer programming, which has been intensively studied over the past decades and is starting to be applied in large-scale problems. Here, we do not need to wait for MILP to find a global solution; theoretical results on the statistical property of early termination are provided.
In this section, we consider the IV quantile model in ((ref)). Our proposal is a method of moment approach:
where $\mathcal{B}\subseteq\mathbb{R}^{p}$ is a convex set. In practice, we can choose $\mathcal{B}=\mathbb{R}^{p}$ or a bounded rectangular subset of $\mathbb{R}^{p}$. The above estimator is based on the fact that $EZ_{i}(\mathbf{1}\{y_{i}-X_{i}'\beta\leq0\}-\tau)=0$ for $\beta=\beta_{*}$. Of course we can replace $Z_{i}$ with transformations of $Z_{i}$. The idea of the estimator is to find a value $\beta$ to minimize the “magnitude” of the empirical version $E_{n}Z_{i}(\mathbf{1}\{y_{i}-X_{i}'\beta\leq0\}-\tau)$.
The estimator ((ref)) differs from GMM in that we use the $\ell_{\infty}$-norm, instead of the $\ell_{2}$-norm. The choice of $\ell_{\infty}$-norm over $\ell_{2}$-norm is due to computational reasons. As we shall see, the formulation with $\ell_{\infty}$-norm in ((ref)) can be cast as an MILP. If we use $\ell_{2}$-norm instead, then the optimization problem would become a mixed integer quadratic program (MIQP), which is the formulation in chen2017exact.\footnote{In their Appendix C3, an MILP formulation is provided, but it requires much more binary variables. Their formulation needs $n+n(n\text{\textminus}1)/2$ binary variables, while our formulation requires $n$ binary variables.} However, as pointed out in hemmecke2010nonlinear,burer2012milp,mazumder2017discrete, it is quite well known in the integer programming community that current algorithms for MILP problems are a much more mature technology than MIQP. For this reason, we use the formulations in ((ref)).
We now show that the estimator ((ref)) can be cast as an MILP. The key is to introduce $n$ binary variables and use constraints to force them to represent $\mathbf{1}\{Y_{i}-X_{i}'\beta\leq0\}$.
Let $\xi_{i}\in\{0,1\}$. Suppose that $M>0$ is an arbitrary number such that $\max_{1\leq i\leq n}|Y_{i}-X_{i}'\hat{\beta}|\leq M$. Notice that this is not a statistical tuning parameter since we can choose any large enough $M>0$. The key insight is to realize that imposing the constraint $-M\xi_{i}\leq Y_{i}-X_{i}'\beta\leq M(1-\xi_{i})$ will force $\xi_{i}$ to behave like $\mathbf{1}\{Y_{i}-X_{i}'\beta\leq0\}$. To see this, consider the following two cases (ignoring the case of $Y_{i}-X_{i}'\beta=0$): (1) $Y_{i}-X_{i}'\beta<0$ and (2) $Y_{i}-X_{i}'\beta>0$. In Case (1), $\xi_{i}=1$ is the only possibility to make $-M\xi_{i}\leq Y_{i}-X_{i}'\beta\leq M(1-\xi_{i})$ hold. Similarly, in Case (2), $\xi_{i}=0$ is the only choice of $\xi_{i}$ in $\{0,1\}$ to satisfy the constraint. Hence, we need to consider variables $\xi_{i}\in\{0,1\}$ and $\beta\in\mathbb{R}^{p}$ such that $-M\xi_{i}\leq Y_{i}-X_{i}'\beta\leq M(1-\xi_{i})$.
In order to minimize $\|E_{n}Z_{i}(\xi_{i}-\tau)\|_{\infty}$, we introduce an auxiliary variable $t\geq0$ with the constraint $-t\leq E_{n}Z_{i,j}(\xi_{i}-\tau)\leq t$ for $j\in\{1,...,L\}$, where $Z_{i,j}$ is the $j$th component of $Z_{i}$. By minimizing $t$, we equivalently achieve minimizing $\|E_{n}Z_{i}(\xi_{i}-\tau)\|_{\infty}$. To summarize, the final MILP formulation reads
In the case of $Y_{i}-X_{i}'\beta=0$, we have an indeterminancy since both $\xi_{i}=0$ and $\xi_{i}=1$ would satisfy $-M\xi_{i}\leq Y_{i}-X_{i}'\beta\leq M(1-\xi_{i})$. However, for most of the design matrices, $\{i:\ Y_{i}-X_{i}'\beta=0\}$ is empty. If we encounter a lot of zeros for $Y_{i}-X_{i}'\beta$ in the solution, we can simply incorporate a small wedge to solve the indeterminacy: $-M\xi_{i}+D\leq Y_{i}-X_{i}'\beta\leq M(1-\xi_{i})$, where $D>0$ is a very small number, such as machine precision tolerance. In our experience, this is not necessary and does not make a difference in the solution.
We now derive the rate of convergence of $\hat{\beta}$. We also discuss how the rate is affected if we terminate MILP before a global solution is reached. A practical guide for early termination is provided and its theoretical validity is also established.
We start with the following simple high-level condition for identification. Let us introduce the following notations. Recall the notations $G(\beta)=EZ_{i}(\mathbf{1}\{Y_{i}-X_{i}'\beta\leq0\}-\tau)$, $G_{n}(\beta)=n^{-1}\sum_{i=1}^{n}Z_{i}(\mathbf{1}\{Y_{i}-X_{i}'\beta\leq0\}-\tau)$ and $H_{n}(\beta)=\sqrt{n}(G_{n}(\beta)-G(\beta))$.
Assumption (ref) guarantees the identification of $\beta_{*}$ and can be verified using primitive conditions similar to Assumption 2 in chernozhukov2006instrumental. In this paper, we do not consider the case with weak identification.\footnote{Inference under potentially weak instruments is quite challenging even for linear IV models. For joint inference on the entire vector $\beta$ or all the coefficients of the endogenous variables, we can rely on the method proposed in Chernozhukov2008. However, for subvector inference (inference only on part of endogenous variables), it is quite challenging even in the linear IV models, for which weak identification has been thoroughly understood only for homoscedastic errors; see e.g., guggenberger2012asymptotic. } We also assume that the empirical process for $Z_{i}(\mathbf{1}\{Y_{i}-X_{i}'\beta\leq0\}-\tau)$ is globally Glivenko-Cantelli and locally Donsker.
Assumption (ref) is not difficult to verify. For example, straight-forward arguments using Lemmas 2.6.15 and 2.6.18 in van1996weak imply that under enough moments of $\|Z_{i}\|_{2}$, the entropy condition in Theorem 2.14.1 therein holds, which means that $E\sup_{\|\beta-\beta_{*}\|_{2}\leq c}\|H_{n}(\beta)\|_{2}=O(1)$. Since we typically terminate the MILP algorithm before a global solution is found, we would like to consider the properties of estimations from early termination.
Theorem (ref) says that when $\|G_{n}(\hat{\beta})\|_{\infty}$ is small, the rate for $\|\hat{\beta}-\beta_{*}\|_{2}$ is $\|G_{n}(\hat{\beta})\|_{\infty}+n^{-1/2}$. Notice that we observe $\|G_{n}(\hat{\beta})\|_{\infty}$ in the MILP algorithm. Hence, we can terminate it once it reaches certain threshold. A natural threshold is $\|G_{n}(\beta_{*})\|_{\infty}$. Although we cannot really compute $\|G_{n}(\beta_{*})\|_{\infty}$ in practice, we can provide a finite-sample bound for it using the moderate deviation result for self-normalized sums. Let $Z_{1,j}$ denote the $j$-th component of $Z_{i}\in\mathbb{R}^{L}$.
In practice, we can simply take $\alpha=1/n$ and thus Lemma (ref) tells us that for $n$ not too small, we have \[ P\left(\|G_{n}(\beta_{*})\|_{\infty}>Q_{*}\right)\leq4Ln^{-2}, \] where $Q_{*}=\Phi^{-1}(1-n^{-2})n^{-1}\sqrt{\max_{1\leq j\leq L}\sum_{i=1}^{n}Z_{i,j}^{2}}$. Notice that $Q_{*}$ can be explicitly computed from the data. Moreover, we know that $Q_{*}=O_{P}(\sqrt{n^{-1}\log n})$. Therefore, if we stop the MILP algorithm once $\|G_{n}(\hat{\beta})\|_{\infty}\leq Q_{*}$, Lemma (ref) and Theorem (ref) imply that $\|\hat{\beta}-\beta_{*}\|_{2}=O_{P}(\sqrt{n^{-1}\log n})$. As we have seen in Section (ref), this is more than enough for the $k$-step correction to yield an estimator that is asymptotically equivalent to GMM.
Now we provide simulation results to illustrate this point. We find that the MILP algorithm reaches $Q_{*}$ within seconds. Let $p=20$. We generate $Y_{i}=X_{i}'\theta+(X_{i}'\gamma)U_{i}$, where $X_{i}$ and $U_{i}$ are generated from the uniform distribution on $(0,1)$. Entries of $\theta$ and $\gamma$ are randomly generated from the uniform distribution on $(0,1)$. We set $\tau=0.7$. The starting point of the MILP algorithm is generated from $N(0,I_{p})$. In Table (ref), we report the frequency of $\|G_{n}(\hat{\beta})\|_{\infty}\leq Q_{*}$ based on 1000 simulations.
As we can see from Table (ref), we only need to run the algorithm for 10 seconds to ensure that $\|G_{n}(\hat{\beta})\|_{\infty}\leq Q_{*}$, which implies $\|\hat{\beta}-\beta_{*}\|_{2}=O_{P}(\sqrt{n^{-1}\log n})$.
When $p\gg n$ and $\beta$ is a sparse vector, the model ((ref)) becomes a high-dimensional IV quantile model. Although our $k$-step correction framework does not cover the case of $p$ growing with $n$, we still present the formulation of estimating high-dimensional IV quantile regression to demonstrate the generality of MILP. In high dimensions, successful estimation relies on proper regularization on $\beta$. Similar to the regularization in Dantzig selector for linear models (candes2007dantzig), we propose
where $\lambda\asymp\sqrt{n^{-1}\log p}$ is tuning parameter.
Similar to the formulation in Section (ref), we can cast the above problem as an MILP. To account for the $\ell_{1}$-norm in the objective function, we decompose each entry of $\beta$ into the positive and negative part: we write $\beta_{j}=\beta_{j}^{+}-\beta_{j}^{-}$ with $\beta_{j}^{+},\beta_{j}^{-}\geq0$. Then the above problem can be rewritten as
The censored regression proposed by powell1986censored reads
where $\rho_{\tau}(x)=x(\tau-\mathbf{1}\{x\leq0\})$ is the “check” function for a given $\tau\in(0,1)$ and $\{(Y_{i},X_{i})\}_{i=1}^{n}$ is the observed data. Notice that this is a nonconvex and non-smooth optimization problem. Computationally it might not be very attractive, especially when the dimensionality is large. The literature has seen alternative estimators that explicitly model the probability of being censored; see e.g., buchinsky1998alternative,chernozhukov2002three. Recently, there is work in high-dimensional statistics (e.g., muller2016censored) studying the statistical properties of
where $\lambda\asymp\sqrt{n^{-1}\log p}$ is a tuning parameter. However, discussions regarding the computational burden for the above estimator are not common. Here, we case the problem ((ref)) as a MILP. Since problem ((ref)) is a special case of problem ((ref)) with $\lambda=0$, our framework can be used for the computation of both ((ref)) and ((ref)).
We introduce variables $\zeta_{i}^{+},\zeta_{i}^{-}\geq0$ to denote the positive and negative parts of $Y_{i}-\max\{X_{i}'\theta,0\}$: $Y_{i}-\max\{X_{i}'\theta,0\}=\zeta_{i}^{+}-\zeta_{i}^{-}$. Similarly, we introduce $r_{i}^{+},r_{i}^{-}\geq0$ such that $X_{i}'\theta=r_{i}^{+}-r_{i}^{-}$; also, let $\theta_{j}^{+},\theta_{j}^{-}\geq0$ satisfy $\theta_{j}=\theta_{j}^{+}-\theta_{j}^{-}$. As in Section (ref), we use $\xi_{i}\in\{0,1\}$ to represent $\mathbf{1}\{X_{i}'\theta<0\}$ by imposing $-\xi_{i}M\leq X_{i}'\theta\leq(1-\xi_{i})M$, where $M>0$ is any number satisfying $\|X\hat{\theta}\|_{\infty}\leq M$.
Notice that $\max\{X_{i}'\theta,0\}=r_{i}^{+}$ if we can force one of $r_{i}^{+}$ and $r_{i}^{-}$ to be exactly zero. The key idea to achieve this is to impose $0\leq r_{i}^{+}\leq M(1-\xi_{i})$ and $0\leq r_{i}^{-}\leq M\xi_{i}$. If $X_{i}'\theta>0$, then $\xi_{i}=0$, which forces $r_{i}^{-}=0$; if $X_{i}'\theta<0$, then $\xi_{i}=1$, which forces $r_{i}^{+}=0$. Now we write down the MILP formulation for ((ref)):
Consider the following moment condition: \[ P\left(Y_{i}\leq\max\{X_{i}'\beta,C_{i}\}\mid Z_{i}\right)=\tau, \] where we observe i.i.d $\{(Y_{i},X_{i},Z_{i},C_{i})\}_{i=1}^{n}$. chernozhukov2015quantile proposed an estimator strategy that uses a control variable. Here, we consider a direct approach based on the above moment condition:
Now we rewrite ((ref)) as an MILP. Similar to Section (ref), we shall introduce binary variables for the max function. Then we use additional binary variables for the indicator function.
We start by introducing $r_{i}^{+},r_{i}^{-}\geq0$ and $\xi_{i}\in\{0,1\}$ such that $X_{i}'\beta-C_{i}=r_{i}^{+}-r_{i}^{-}$, $-\xi_{i}M\leq X_{i}'\beta-C_{i}\leq(1-\xi_{i})M$, $0\leq r_{i}^{+}\leq M(1-\xi_{i})$ and $0\leq r_{i}^{-}\leq M\xi_{i}$, where $M>0$ is a large enough number. As explained in Section (ref), these constraints will force $\xi_{i}$ to behave like $\mathbf{1}\{X_{i}'\beta<C_{i}\}$ and ensure that one of $r_{i}^{+}$ and $r_{i}^{-}$ is exactly zero, thus $r_{i}^{+}=\max\{X_{i}'\beta-C_{i},0\}$. Hence, $Y_{i}-\max\{X_{i}'\beta,C_{i}\}\leq0$ becomes $Y_{i}-C_{i}-r_{i}^{+}\leq0$.
Now we introduce $q_{i}\in\{0,1\}$ such that $-Mq_{i}\leq Y_{i}-C_{i}-r_{i}^{+}\leq M(1-q_{i})$. Again, this constraint would would make $q_{i}$ behave like $\mathbf{1}\{Y_{i}-C_{i}-r_{i}^{+}\leq0\}$. Therefore, we only need to introduce an extra variable $t\geq0$ to serve as $\|E_{n}Z_{i}(\xi_{i}-\tau)\|_{\infty}$. The final formulation reads
We now conduct simulations for Algorithm (ref) and set $m=500$, . As discussed above, when we run the MILP on the subsample of size $m$, we can expect the rate of convergence to be $\sqrt{m^{-1}\log m}$. We report the coverage probabilities of 95% confidence intervals for $\beta_{j}(\tau)$ for $1\leq j\leq p$.
The results provide quite favorable evidence for the proposed estimator. The empirical coverage probability is close to the nominal level of confidence intervals. This is quite impressive for large $n$. When $n=5\times10^{6}$ and $m=500$, we only use $0.01\%$ of the data for MILP. This still yields good performance in terms of coverage probability of confidence intervals.
We generate \[ Y_{i}=\alpha_{0}+\alpha D_{i}+\sum_{j=1}^{q}W_{i,j}\gamma_{j}+\sum_{j=1}^{q}D_{i}W_{i,j}\theta_{j}+\left(2\sqrt{3}q+\sum_{j=1}^{q}W_{i,j}\lambda_{j}+\sum_{j=1}^{q}D_{i}W_{i,j}\pi_{j}\right)V_{i}, \] where $W_{i,j}\sim\text{uniform}(-\sqrt{3},\sqrt{3})$, $V_{i}\sim N(0,1)$ and $(D_{i},S_{i})$ are mutually independent. We generate $(S_{i},D_{i})\in\{(1,1),(1,0),(0,0)\}$ with probability 0.42, 0.25 and 0.33; these frequencies are estimated from the JTPA data. We set $\alpha_{0}=\alpha=\gamma_{j}=\theta_{j}=\lambda_{j}=\pi_{j}=1$ and $q=10$. We set $X_{i}=(1,D_{i},W_{i,1},...,W_{i,q},D_{i}W_{i,1},...,D_{i}W_{i,q})\in\mathbb{R}^{p}$ and $Z_{i}=(1,S_{i},W_{i,1},...,W_{i,q},S_{i}W_{i,1},...,S_{i}W_{i,q})'\in\mathbb{R}^{p}$ with $p=2q+2=22$. One can easily verify that for any $\tau\in(0,1)$, \[ P\left(Y_{i}\leq X_{i}'\beta(\tau)\mid Z_{i}\right)=\tau, \] with $\beta(\tau)=(1+2\sqrt{3}+c_{\tau},1+c_{\tau},1+c_{\tau},...,1+c_{\tau})'\in\mathbb{R}^{p}$ and $c_{\tau}$ is the $\tau$-th quantile of $N(0,1)$.
We use Gurobi 8.1 for mixed integer programming and implement it in Matlab R2019a. For the starting values, we first randomly generate $\beta_{{\rm start}}$ from $N(0,I_{p})$ and compute the starting points for $\xi_{i}$ and $t$ using $\mathbf{1}\{Y_{i}-X_{i}'\beta_{{\rm start}}\leq0\}$ and $\|E_{n}Z_{i}(\mathbf{1}\{Y_{i}-X_{i}'\beta_{{\rm start}}\leq0\}-\tau)\|_{\infty}$, respectively. We run the MILP program on a randomly selected subsample of size $m=500$ and terminate the optimization algorithm after 5 seconds although a strict guarantee for global solutions would typically take a few hours. The output of the MILP program is used as the initial estimate for the $k$-step iteration. We report the performance of the $k$-step corrected estimator in Table (ref), which is based on 800 random samples. In the experiments, we choose $n\in\{5\times10^{3},5\times10^{4},10^{6}\}$, which means $m/n\in\{10\%,1\%,0.05\%\}$. Five quantiles $\tau\in\{0.15,0.25,0.5,0.75,0.85\}$ are considered. Since the usual $\chi^{2}$-test involves inverting the $p\times p$ asymptotic variance matrix, the finite-sample performance of the test might not be ideal as $p$ increases. One main reason is that the asymptotic variance matrix is estimated and could be badly conditioned even if $p$ is only 20. For this reason, we consider the following non-pivotal test statistic that avoids inverting a large matrix.
In Table (ref), we consider inference of $\beta(\tau)=(\beta_{1}(\tau),...,\beta_{p}(\tau))'$. Notice that $\beta_{2}(\tau)$ corresponds to the coefficient for $D_{i}$, $(\beta_{3}(\tau),...,\beta_{2+q}(\tau))$ for $(W_{i,1},...,W_{i,q})$ and $(\beta_{3+q}(\tau),...,\beta_{p}(\tau))$ for $(D_{i}W_{i,1},...,D_{i}W_{i,q})$. We consider the coverage probability of confidence intervals (sets) for various components (and subvectors) of $\beta(\tau)$. We see that when the sample size is 5000, ellipsoid confidence sets, which requires inverting a $22\times22$ estimated asymptotic variance, have some undercoverage while rectangular confidence sets still have accurate coverage probabilities. For larger sample size, all the confidence intervals (sets) have quite accurate performance. We include the results for $n=10^{6}$ to emphasize the point that for extremely large samples, our method is still quite fast as it takes less than 15 seconds to compute one sample. In this case, the main factor that limits the speed is the memory.
Consider $Y_{i}=X_{i}+Z_{i}\varepsilon_{i}$ and $X_{i}=Z_{i}V_{i}$, where $\varepsilon_{i}\sim\text{Exp}(\lambda)$,\footnote{$\text{Exp}(\lambda)$ denotes the exponential distribution with parameter $\lambda$. The mean of this distribution is $\lambda^{-1}$.} $V_{i}\sim{\rm uniform}(0,1)$ and $Z_{i}\sim{\rm uniform}(0,2)$ are mutually independent. Let $g(W_{i};\beta)=Z_{i}\mathbf{1}\{Y_{i}\leq X_{i}\beta\}$. We can explicitly compute the population derivative: for $\beta>1$, \[ \Gamma(\beta)=\frac{dEg(W_{i};\beta)}{d\beta}=\frac{1-[\lambda(\beta-1)+1]\exp\left(\lambda(1-\beta)\right)}{\lambda(\beta-1)^{2}}. \]
In Table (ref), we evaluate the performance of two estimators in terms of root-mean-squared error (RMSE). The first estimator $\hat{\Gamma}_{\text{tuning-free}}$ is the one proposed in Section (ref) computed using $\sqrt{n}$ bootstrap samples. The second estimator $\hat{\Gamma}_{\text{kernel}}$ is the kernel estimator with Gaussian kernel and bandwidth given by Silverman's rule of thumb. We see that the tuning-free estimator is a good alternative to well-tuned kernel estimators. Since theoretically $\hat{\Gamma}_{\text{kernel}}$ has a faster rate of convergence, we see that the constants in the bandwidth choice can be quite important for performance in practice. The tuning-free estimator is attractive in that one does not have to make these tricky choices.
In Section (ref), we mentioned the problem of investigating the effect of JTPA. We now provide more details. The participation status will be denoted by $D_{i}\in\{0,1\}$, where $D_{i}=1$ means that individual $i$ participates in the program. The random offers, denoted by $S_{i}\in\{0,1\}$, will be used as instruments, where $S_{i}=1$ means that individual $i$ has an offer to participate. Following Chernozhukov2008, we also consider other 13 exogenous variables denoted by $\{W_{i,j}\}_{j=1}^{13}$.\footnote{The data is downloaded from Christian Hansen's \href{http://faculty.chicagobooth.edu/christian.hansen/research/sampdata.zip}{website}.} The outcome variable $Y_{i}$ is earnings. We consider the following model:
where $\tau\in(0,1)$. Under our notation ((ref)), we have \[
\] where $p=L=28$. We rescale $Z_{i}$ such that $n^{-1}\sum_{i=1}^{n}Z_{i,2:28}Z_{i,2:28}'=I_{27}$ for $j\in\{1,...,L\}$, where $Z_{i,2:28}$ denotes the vector $Z_{i}$ with the first component removed. We are interested in $\alpha(\tau)$, which denotes the overall effect of JTPA, as well as $\theta_{j}(\tau)$, which measures the heterogeneity of the effect. Following Chernozhukov2008, we consider $\tau\in\{0.15,0.25,0.5,0.75,0.85\}$. We use the proposed methods to test the model parameters. In Table (ref), we report the test statistics for three hypotheses at these 5 quantiles.
We find strong evidence for heterogeneity in treatment effects. The coefficient $\alpha(\tau)$ for the treatment status is not significant\footnote{Since we include interaction terms in ((ref)), the estimates for $\alpha(\tau)$ would be not represent the “average” effect if there is heterogeneity in the treatment effect; hence, these results should not be directly comparable to results reported in abadie2002instrumental and Chernozhukov2008.}; moreover, the controls by themselves are insignificant or barely significant at 5% or 10% level. However, we reject the insignificance of $\theta_{j}(\tau)$ coefficients, which means that the covariates are very informative on the magnitude of treatmente effects.
We also present two prominent patterns in the heterogeneity. In Figure (ref), we see that the treatment effects of JTPA depend on the marriage status and prior employment records. On the very left tail, the the treatment effect does not depend on the marriage status, whereas among higher-income individuals, married ones benefit more from JTPA than the unmarried. Moreover, although prior employment status is not related to the effect of JTPA for lower-income individuals, we find that on the upper level of the income spectrum, those with insufficient prior employment records tend to benefit less than those with at least 13 weeks of employment in the past year.
In this paper, we provide an alternative estimation and inference method for IVQR and related problems. This alternative approach is needed as we show that the GMM formulation of IVQR (or even a reasonable approximation) is computationally NP-hard. Our proposal focuses on obtaining good statistical properties instead of trying to improve optimization algorithms for GMM. The proposed estimator can transform, in a computationally efficient manner, an inconsistent initial estimator into one that is asymptotically equivalent to GMM. The initial estimator is obtained via MILP. We also propose a tuning-free method for estimating the Jacobian of the moment condition in a non-smooth GMM model. We illustrate our proposal in simulated and empirical data.