Exact computation of GMM estimators for instrumental variable quantile regression models

TL;DR

This paper introduces an exact computational method for GMM estimators in instrumental variable quantile regression models by formulating the problem as a mixed integer quadratic program, demonstrated through simulations and real data.

Contribution

The paper presents a novel formulation that allows for exact GMM estimation in IVQR models using mixed integer quadratic programming.

Findings

Exact GMM estimators can be computed efficiently.

The method outperforms approximate algorithms in accuracy.

Application to demand for fish illustrates practical usefulness.

Abstract

We show that the generalized method of moments (GMM) estimation problem in instrumental variable quantile regression (IVQR) models can be equivalently formulated as a mixed integer quadratic programming problem. This enables exact computation of the GMM estimators for the IVQR models. We illustrate the usefulness of our algorithm via Monte Carlo experiments and an application to demand for fish.

Tables

Table 1. Table 1: MIQP computation time (CPU seconds)

$\tau$	mean	min	median	max
0.25	94	37	92	197
0.5	348	104	333	989
0.75	86	33	84	186

Table 2. Table 2: Finite-sample performance of the GMM estimator

	mean		median
	bias	RMSE	bias	MAE
${\theta }_0\left(0.25\right)$	0.0109	0.2436	0.0012	0.1643
${\theta }_1\left(0.25\right)$	-0.0327	0.3554	-0.0048	0.2309
${\theta }_2\left(0.25\right)$	0.0003	0.1642	0.0031	0.1008
${\theta }_3\left(0.25\right)$	0.0064	0.2232	-0.0068	0.1522
${\theta }_0\left(0.5\right)$	0.0161	0.2498	-0.0037	0.1724
${\theta }_1\left(0.5\right)$	-0.0412	0.3241	-0.0316	0.2315
${\theta }_2\left(0.5\right)$	-0.0012	0.1561	0.0066	0.1038
${\theta }_3\left(0.5\right)$	0.0012	0.2047	0.0031	0.1396
${\theta }_0\left(0.75\right)$	0.0187	0.3046	0.0055	0.1849
${\theta }_1\left(0.75\right)$	-0.0358	0.3425	-0.0264	0.2235
${\theta }_2\left(0.75\right)$	-0.0022	0.1820	0.0062	0.1181
${\theta }_3\left(0.75\right)$	0.0035	0.2393	-0.0016	0.1538

Table 3. Table 3: Comparison with asymptotic approximation

	standard deviation	asymptotic
	in simulations	standard error
${\theta }_0\left(0.25\right)$	0.2434	0.2297
${\theta }_1\left(0.25\right)$	0.3539	0.3256
${\theta }_2\left(0.25\right)$	0.1642	0.1572
${\theta }_3\left(0.25\right)$	0.2231	0.2059
${\theta }_0\left(0.5\right)$	0.2493	0.2296
${\theta }_1\left(0.5\right)$	0.3215	0.3049
${\theta }_2\left(0.5\right)$	0.1561	0.1474
${\theta }_3\left(0.5\right)$	0.2047	0.1994
${\theta }_0\left(0.75\right)$	0.3040	0.2744
${\theta }_1\left(0.75\right)$	0.3406	0.3400
${\theta }_2\left(0.75\right)$	0.1820	0.1664
${\theta }_3\left(0.75\right)$	0.2393	0.2283

Table 4. Table 4: Coverage probabilities (95% CI)

	$0.8h_{n,HS}$	$h_{n,HS}$	$1.2h_{n,HS}$
${\theta }_0\left(0.25\right)$	0.930	0.940	0.952
${\theta }_1\left(0.25\right)$	0.906	0.914	0.918
${\theta }_2\left(0.25\right)$	0.912	0.924	0.934
${\theta }_3\left(0.25\right)$	0.916	0.926	0.938
${\theta }_0\left(0.5\right)$	0.936	0.944	0.950
${\theta }_1\left(0.5\right)$	0.938	0.944	0.952
${\theta }_2\left(0.5\right)$	0.950	0.958	0.966
${\theta }_3\left(0.5\right)$	0.896	0.916	0.928
${\theta }_0\left(0.75\right)$	0.896	0.916	0.928
${\theta }_1\left(0.75\right)$	0.922	0.938	0.944
${\theta }_2\left(0.75\right)$	0.892	0.896	0.908
${\theta }_3\left(0.75\right)$	0.918	0.928	0.942

Table 5. Table 5: IVQR estimation of demand elasticity

	$\tau =0.25$	$\tau =0.5$	$\tau =0.75$
Specification (5.1)
Estimation method: GMM via the MIQP implementation
${\hat{\theta }}_1\left(\tau \right)$	-1.0880	-0.8876	-0.9755
std. dev.	0.4773	0.5056	0.3027
95% CI	$(-2.0234,-0.1525)$	$(-1.8787,0.1034)$	$(-1.5689,-0.3822)$
Estimation method: Inverse QR
${\hat{\theta }}_1\left(\tau \right)$	-1.3680	-0.8860	-1.2685
std. dev.	0.5704	0.4673	0.3911
95% CI	$(-2.486,-0.250)$	$(-1.802,0.030)$	$(-2.035,-0.502)$
Specification (5.1) augmented with day fixed effects
Estimation method: GMM via the MIQP implementation
${\hat{\theta }}_1\left(\tau \right)$	-0.6915	-0.7152	-1.0904
std. dev.	0.3253	0.4828	0.2465
95% CI	$(-1.3290,-0.0540)$	$(-1.6616,0.2312)$	$(-1.5735,-0.6074)$
Estimation method: Inverse QR
${\hat{\theta }}_1\left(\tau \right)$	-1.3635	-0.5950	-1.1790
std. dev.	0.5304	0.4398	0.3653
95% CI	$(-2.403,-0.324)$	$(-1.457,0.267)$	$(-1.895,-0.463)$

Equations

$$Y=W^{\prime }\theta (U),$$

$$U|X,Z\sim \mathrm{Uniform}\left(0,1\right).$$

$$P\left(Y\le W^{\prime }\theta (\tau )|X,Z\right)=P\left(U\le \tau |X,Z\right)=\tau .$$

$$\hat{\alpha }(\tau )\equiv \arg {\inf }_{\alpha \in \mathcal{A}}{\hat{\gamma }}_{\tau }{\left(\alpha \right)}^{\prime }A{\hat{\gamma }}_{\tau }\left(\alpha \right),$$

$$\left({\hat{\beta }}_{\tau }\left(\alpha \right),{\hat{\gamma }}_{\tau }\left(\alpha \right)\right)\equiv \arg {\inf }_{\left(\beta ,\gamma \right)\in \mathcal{B}\times \mathcal{G}}\frac{1}{n}{\sum }_{i=1}^n{\rho }_{\tau }\left(Y_i-D_i^{\prime }\alpha -X_i^{\prime }\beta -{\Psi }_i^{\prime }\gamma \right),$$

$$\hat{\theta }(\tau )=\left(\hat{\alpha }(\tau ),{\hat{\beta }}_{\tau }\left(\hat{\alpha }(\tau )\right)\right).$$

$$E\left[\left(1\left\{Y\le W^{\prime }\theta \left(\tau \right)\right\}-\tau \right)L\right]=0,$$

$${\hat{\theta }}_{GMM}\left(\tau \right)=\arg {\inf }_{\theta \in \Theta }{s}_{\tau }{\left(\theta \right)}^{\prime }G\hat{Q}{G}^{\prime }{s}_{\tau }\left(\theta \right),$$

$$\underset{e=\left(e_1,...,e_n\right),\theta \in \Theta }{\inf }{\left(e-\tau \right)}^{\prime }G\hat{Q}{G}^{\prime }\left(e-\tau \right)$$

$$e_i\left(-M_i-ϵ\right)<Y_i-W_i^{\prime }\theta \le \left(1-e_i\right)M_i,i\in \{1,...,n\},$$

$$e_i\in \{0,1\},i\in \{1,...,n\},$$

$$M_i\equiv {\max }_{\theta \in \Theta }\left|Y_i-W_i^{\prime }\theta \right|,i\in \{1,...,n\}.$$

$$\hat{Q}={\left[\tau \left(1-\tau \right)n^{-1}{\sum }_{i=1}^n{L}_i{L}_i^{\prime }\right]}^{-1}$$

$$\sqrt{n}({\hat{\theta }}_{GMM}\left(\tau \right)-\theta \left(\tau \right))\stackrel{d}{⟶}N(0,\Omega ),$$

$$\Omega =\tau \left(1-\tau \right){\left[{\Sigma }_{WL}{\Sigma }_{LL}^{-1}{\Sigma }_{WL}^{\prime }\right]}^{-1},{\Sigma }_{WL}=E\left[f_{{\epsilon }_{\tau }}\left(0|W,Z\right)W{L}^{\prime }\right],{\Sigma }_{LL}=E\left[L{L}^{\prime }\right].$$

$${\hat{\Sigma }}_{WL}\equiv n^{-1}{\sum }_{i=1}^n\left[K\left({\hat{\epsilon }}_{\tau ,i}/h_n\right)/h_n\right]W_i{L}_i^{\prime },$$

$$Y$$

$$D_1$$

$$V=\left[\begin{array}[]{cccc}1&0.4&0.6&-0.2\\ 0.4&1&0&0\\ 0.6&0&1&0\\ -0.2&0&0&1\end{array}\right].$$

$$Y={\theta }_0(U)+{\theta }_1(U)D_1+{\theta }_2(U)D_2+{\theta }_3(U)D_3,$$

$${\theta }_0\left(\tau \right)=1+0.5F_{\epsilon }^{-1}\left(\tau \right),{\theta }_1\left(\tau \right)=1+F_{\epsilon }^{-1}\left(\tau \right),{\theta }_2\left(\tau \right)=1+0.25F_{\epsilon }^{-1}\left(\tau \right),{\theta }_3\left(\tau \right)=1+0.15F_{\epsilon }^{-1}\left(\tau \right).$$

$$Y={\theta }_0\left(U\right)+{\theta }_1\left(U\right)D$$

Full text

Keywords:* generalized method of moments, instrumental variable, quantile regression, endogeneity, mixed integer optimization*

JEL codes: C21, C26, C61, C63

1 Introduction

The instrumental variable quantile regression (IVQR) and related models have been increasingly popular for studying the impacts of possibly endogenous covariates on the distribution of the outcome of interest. See a recent review by Chernozhukov and Hansen (2013) and references therein for the latest developments in identification, estimation, and inference as well as the list of empirical applications.

The IVQR model admits conditional moment restrictions which can be used to construct the estimating equations for the GMM estimation of the model parameters. However, the sample counterparts of the IVQR estimating equations are discontinuous in the parameters so that the resulting GMM estimation problem becomes a non-convex and computationally non-trivial optimization problem. Honoré and Hu (2004) provided a heuristic for computing the IVQR GMM estimates. Chernozhukov and Hansen (2006) developed the inverse quantile regression (QR) estimator that is not directly a GMM estimator but can be shown to be asymptotically equivalent to the IVQR GMM estimator. Xu and Burer (2017) proposed an alternative algorithm for computing the inverse QR estimator. The Markov chain Monte Carlo (MCMC) based Laplace type estimator of Chernozhukov and Hong (2003) can also be used as an approximation of the IVQR GMM estimator but it requires careful tunning in the MCMC implementation. Kaplan and Sun (2015) proposed a smoothed estimating equation approach which facilitates the GMM computation problem but requires the choice of the smoothing parameter.

In this paper, we are concerned with exact computation of the GMM estimates of the IVQR parameters. As pointed out by Andrews (1997), heuristic algorithms for computation of GMM estimates that do not guarantee to find the exact global optimum or a specific level of approximation to the global optimum may result in extremum estimators which could exhibit statistical behavior that is quite different from that established by theory. This source of computational uncertainty may impact on the empirical results. Hence, as a complement to the previous work on the IVQR computation, our paper provides a method for exact computation of the IVQR estimates within the classical GMM framework.

Our computational algorithm is based on the method of mixed integer optimization (MIO). Specifically, we show that the IVQR GMM estimation problem can be equivalently formulated as a mixed integer quadratic programming (MIQP) problem. Thanks to the developments in MIO solution algorithms and fast computing environments, this reformulation allows us to solve for the exact GMM estimates by using the modern efficient MIO solvers. Well-known numerical solvers such as CPLEX and Gurobi can be used to effectively solve large-scale MIQP problems. See Bertsimas, King, and Mazumder (2016, Section 2.1) for discussions on computational advances in solving the MIO problems. See also Florios and Skouras (2008), Kitagawa and Tetenov (2015), Bertsimas, King, and Mazumder (2016), and Chen and Lee (2016) for related but distinct work on solving non-convex optimization problems in statistics and econometrics via the MIO approach.

The rest of this paper is organized as follows. In Section 2, we summarize the setup of the IVQR model and the inverse quantile regression method of Chernozhukov and Hansen (2006). In Section 3, we present the MIQP formulation of the IVQR GMM estimation problem. We conduct a simulation study of the performance of the MIQP based GMM estimates in Section 4 and illustrate the application of our computation approach in a real data exercise concerning the demand estimation in Section 5. We then conclude the paper in Section 6.

2 The instrumental variable quantile regression model

Let $Y$ be an outcome of interest. We consider the quantile regression model under endogeneity, which is characterized by the structural equation

[TABLE]

where $U$ is an unobserved scalar random variable and $W=(D,X)$ is a vector of covariates. The covariates $D$ may not be independent of $U.$ We assume that there is a vector of instrumental variables, denoted as $Z$, which can be excluded from (2.1) but can influence the endogenous variables $D$ such that $\dim(Z)\geq\dim(D)$ and

[TABLE]

Assume that the function $\theta(\cdot)$ in (2.1) is a measurable function such that the mapping $\tau\mapsto W^{\prime}\theta(\tau)$ is strictly increasing in $\tau$ for almost every realization of $W$. Under these assumptions, it follows that

[TABLE]

Given a random sample, $\left(Y_{i},W_{i},Z_{i}\right)_{i=1}^{n}$ of $n$ observations, we are interested in the estimation of $\theta(\tau)$ for some fixed values of $\tau\in(0,1)$.

The model set forth so far is the well known linear IVQR model which has been studied by Chernozhukov and Hansen (2004, 2005, 2006, 2008), Chernozhukov, Hansen, and Jansson (2007, 2009), and Kaplan and Sun (2015) among many others. Note that, when there is no endogenous covariate, this model reduces to the conventional linear quantile regression model of Koenker and Bassett (1978) where $W=X=Z$. In the presence of endogeneity, Chernozhukov and Hansen (2005) provided further modeling assumptions such that the quantile-specific parameter vector $\theta(\tau)$ can be causally interpreted as the structural quantile effect in the setting with counterfactual outcomes.

Chernozhukov and Hansen (2006) developed primitive conditions for the identification of $\theta(\tau)$. They also provided an inverse QR algorithm for the estimation of $\theta(\tau)$. To describe their algorithm, write $\theta=\left(\alpha,\beta\right)$ such that $W^{\prime}\theta(\tau)=D^{\prime}\alpha\left(\tau\right)+X^{\prime}\beta\left(\tau\right)$. Let $\Psi_{i}=\Psi\left(X_{i},Z_{i}\right)$ be a vector of transformations of instruments with $\dim(\Psi_{i})\geq\dim\left(\alpha\right)$. Let $A$ be a given positive definite matrix. The Chernozhukov-Hansen inverse QR procedure proceeds as follows. Let

[TABLE]

where

[TABLE]

$\mathcal{A}$, $\mathcal{B}$ and $\mathcal{G}$ are compact parameter spaces, and the check function $\rho_{\tau}$ is defined by $\rho_{\tau}\left(u\right)=u\left(\tau-1\left\{u<0\right\}\right)$ for $u\in\mathbb{R}$. The inverse QR estimator is then defined by

[TABLE]

In the procedure above, the function $\Psi$ and the matrix $A$ can vary across $\tau$ and be replaced by their consistent estimates. Moreover, the QR objective function can be weighted across observations. See Chernozhukov and Hansen (2006) for further details.

For implementation, Chernozhukov and Hansen (2006) proposed to solve the outer optimization problem (2.3) by the grid search method. The inner optimization problem (2.4) is the standard quantile regression problem, which can be solved very efficiently. Thus, when $\dim(\alpha)=1$, the inverse QR method is computationally appealing because its implementation amounts to solving convex optimization sub-problems within a low-dimensional global search procedure. However, this computational merit diminishes rapidly with the increase of the number of endogenous variables. Instead of performing grid search, Xu and Burer (2017) proposed an alternative method to compute the inverse QR estimator. Their approach is based on exact minimization of the quadratic norm as in (2.3) subject to the optimality conditions for the linear programming formulation of the QR problem of (2.4). Xu and Burer (2017) showed that the resulting computational problem reduces to a quadratic programming problem subject to complementarity constraints for which they developed a branch-and-bound algorithm to compute the exact solution.

3 Exact computation of the GMM based IVQR estimator via the mixed

integer optimization approach

The conditional moment restriction (2.2) can be used to form estimating equations for the GMM estimation of $\theta\left(\tau\right)$. That is,

[TABLE]

where $L$ is a vector of instruments consisting of functions of $X$ and $Z$. As noted by Chernozhukov and Hansen (2006), the inverse QR estimator, which is not directly a GMM estimator, can be shown to be asymptotically equivalent to the GMM estimator with the instruments $L_{\text{CH}}\equiv[X^{\prime},\Psi\left(X,Z\right)^{\prime}]^{\prime}.$

In this paper, we provide an algorithm for directly computing the GMM based IVQR estimator using the orthogonality conditions (3.1). Let $s_{\tau}(t)$ denote the vector $(s_{\tau,i}(t))_{i=1}^{n}$, where $s_{\tau,i}(t)\equiv 1\left\{Y_{i}\leq W_{i}^{\prime}t\right\}-\tau$ for $i\in\{1,...,n\}$. Let $G$ be the $n$-by-$\dim(L)$ matrix whose $i$th row vector is $L_{i}^{\prime}$. Let $\widehat{Q}$ be a given positive definite matrix of dimension $\dim(L)$. The GMM based IVQR estimator of $\theta\left(\tau\right)$, denoted by $\widehat{\theta}_{GMM}\left(\tau\right)$, is given by

[TABLE]

where ${\Theta}$ is the compact parameter space of $\theta$.

We now present our computational algorithm, which is based on the method of mixed integer optimization. We note that the optimization problem (3.2) can be equivalently formulated as the following constrained mixed integer quadratic programming (MIQP) problem:

[TABLE]

where $\epsilon$ is a given small and positive real scalar (e.g. $\epsilon=10^{-6}$ as in our simulation study), and

[TABLE]

We now explain the equivalence between (3.2) and (3.3). Note that, for a given value of $\theta\in\Theta$, the sign constraints (3.4) and the dichotomization constraints (3.5) enforce that $e_{i}=1\{Y_{i}\leq W_{i}^{\prime}\theta\}$ for $i\in\{1,...n\}$. Therefore, solving the constrained MIQP problem (3.3) is equivalent to solving the GMM estimation problem (3.2). This equivalence enables us to employ the modern MIQP solvers to exactly compute the GMM estimator $\widehat{\theta}_{GMM}\left(\tau\right)$. For the implementation, note that the values $\left(M_{i}\right)_{i=1}^{n}$ in the inequality constraints (3.4) can be computed by formulating the maximization problem in (3.6) as linear programming problems, which can be efficiently solved by modern optimization solvers. Hence these values can be easily computed and stored as the input to the MIQP formulation (3.3).

Remark 1*.*

Our MIQP based computational approach can be used to find the exact global solution in the IVQR GMM estimation problem. Modern MIQP solvers employ branch-and-bound type algorithms which maintain along the solution process both the feasible solutions and lower bounds on the optimal objective function value. Therefore, for computationally demanding applications, this feature enables us to solve for an approximate IVQR GMM estimator with a guaranteed approximation error bound, thus facilitating the design of an early stopping rule as described in Chen and Lee (2016, Section 4.3).

We can perform inference on $\theta\left(\tau\right)$ using the GMM estimator $\widehat{\theta}_{GMM}\left(\tau\right)$. As noted by Chernozhukov, Hansen, and Jansson (2009), we can take

[TABLE]

as a convenient and natural choice of the GMM weight matrix. By (2.2), this weight matrix equals the inverse of the variance of $n^{-1/2}\sum\nolimits_{i=1}^{n}s_{\tau,i}(\theta\left(\tau\right))L_{i}$ conditional on $\left(L_{i}\right)_{i=1}^{n}$. Let $\varepsilon_{\tau}\equiv Y-W^{\prime}\theta\left(\tau\right)$. In the GMM estimation (3.2) with $\widehat{Q}$ given by (3.7), it is straightforward to establish via empirical process theory (see e.g., Pakes and Pollard, 1989) that

[TABLE]

where the asymptotic variance matrix $\Omega$ is given by

[TABLE]

We can estimate $\Sigma_{LL}$ by the sample analog $\widehat{\Sigma}_{LL}\equiv n^{-1}\sum\nolimits_{i=1}^{n}L_{i}L_{i}^{\prime}$. Let $\widehat{\varepsilon}_{\tau,i}\equiv Y_{i}-W_{i}^{\prime}\widehat{\theta}_{GMM}\left(\tau\right)$. Following Powell (1986), $\Sigma_{WL}$ can be consistently estimated by

[TABLE]

where $K\left(\cdot\right)$ is a kernel function and $h_{n}$ is a bandwidth sequence satisfying that $h_{n}\longrightarrow 0$ and $\sqrt{n}h_{n}\longrightarrow\infty.$ Specific rule-of-thumb choices of $h_{n}$ can be based on Koenker (1994). See also Chernozhukov and Hansen (2006, Section 3.4) and Chernozhukov and Hansen (2008, Section 4.4) for the estimation of the IVQR variance components. Based on these results, it is therefore straightforward to construct the confidence interval estimates of $\theta\left(\tau\right)$ within the GMM framework.

4 Simulation study

In this section, we study the performance of the GMM estimator $\widehat{\theta}_{GMM}\left(\tau\right)$ in finite-sample simulations. We used the MATLAB implementation of the Gurobi Optimizer (version 7.0) to solve the MIQP problems for all numerical results of this paper.111The MATLAB codes for the computation of $\widehat{\theta}_{GMM}\left(\tau\right)$ are available from the authors. This implementation requires the Gurobi Optimizer, which is freely available for academic purposes. All computations were done on a desktop PC (Windows 7) equipped with 32 GB RAM and a CPU processor (Intel i7-5930K) of 3.5 GHz.

We generated $n=100$ observations from the following simple location scale model:

[TABLE]

where $\Phi$ denotes the cdf of the standard normal random variable, $Z_{1}$, $Z_{2}$ and $Z_{3}$ are independent standard normal random variables, and $\left(\varepsilon,v_{1},v_{2},v_{3}\right)$ is generated independently of $\left(Z_{1},Z_{2},Z_{3}\right)$ from multivariate normal distribution with mean zero and variance $0.25V$ where

[TABLE]

By Skorohod representation, we can rewrite the model (4.1) as

[TABLE]

where $U=F_{\varepsilon}\left(\varepsilon\right)$ with $F_{\varepsilon}$ being the cdf of the unobservable $\varepsilon$, and

[TABLE]

We used 500 simulation repetitions for all simulation experiments. In the GMM estimation, we took $W=\left(1,D_{1},D_{2},D_{3}\right)$ and $L=\left(1,Z_{1},Z_{2},Z_{3}\right)$. The GMM weight matrix $\widehat{Q}$ was constructed based on (3.7). We set the parameter space $\Theta$ in the MIQP problem (3.3) to be the product of the intervals $[\widehat{\theta}_{j,2SLS}-10\widehat{\sigma}_{j,2SLS},\widehat{\theta}_{j,2SLS}+10\widehat{\sigma}_{j,2SLS}]$, where for $j\in\{0,1,2,3\}$, $\widehat{\theta}_{j,2SLS}$ and $\widehat{\sigma}_{j,2SLS}$, respectively denote the parameter estimate and its estimated heteroskedasticity-robust standard error from the two-stage least square regression of $Y$ on the covariates $W$ using $L$ as the instruments. The value of $\epsilon$ in (3.4) was set to be $10^{-6}$.

We now present the simulation results. First, we report the computational performance of our MIQP algorithm for computing the IVQR GMM estimator. Table 1 gives the summary statistics of the MIQP computation time in CPU seconds across simulation repetitions. From this table, we can see that the MIQP problems (3.3) were solved very efficiently in these simulations which incorporated three endogenous covariates. For the two cases with $\tau\in\{0.25,0.75\}$, the computation time was comparable. Both cases could be easily solved with the mean and median computation time not exceeding 100 seconds and the maximum time below 200 seconds. The case of $\tau=0.5$ appeared to be the most computationally demanding but its maximum time remained capped within 17 minutes.

We now study the statistical performance of the IVQR GMM estimator. In Table 2, we report the mean and median biases, root mean squared error (RMSE) and median absolute error (MAE) of the GMM estimators $\widehat{\theta}_{GMM}\left(\tau\right)$ for $\tau\in\{0.25,0.5,0.75\}$. From these results, we find that the GMM estimators performed quite well in terms of estimation bias. Across the three quantile cases, the estimators for $\theta_{1}\left(\tau\right)$ appeared to have larger dispersion in terms of both RMSE and MAE.

It is also interesting to assess how well the finite-sample behavior of the IVQR GMM estimator can be approximated by asymptotic theory. For this purpose, our exact GMM estimator can be used to eliminate the unquantified uncertainty on the solution inaccuracy that might emerge in a heuristic optimization procedure. In Table 3, we calculated the asymptotic standard error based on the formula (3.9) evaluated at true parameter values of the simulation design. This quantity was then compared to standard deviation of $\widehat{\theta}_{GMM}\left(\tau\right)$ in simulations. The results of Table 3 indicate that the finite-sample standard error of the GMM estimator in this simulation setup, though being slightly larger, can be well approximated by the asymptotic standard error.

In practice, for carrying out inference, the asymptotic variance of the GMM estimator has to be estimated. We used the Gaussian kernel in the estimation of $\Sigma_{WL}$. The bandwidth sequence $h_{n}$ in (3.10) was based on the Hall-Sheather bandwidth choice, which was suggested by Koenker (1994) and also used by Chernozhukov, Hansen, and Jansson (2009). We also checked the sensitivity of the inference results with respect to this bandwidth choice. Specifically, we reported in Table 4 the finite-sample cover probabilities of the 95% confidence interval (CI) estimates for $\theta\left(\tau\right)$, which were constructed based on the normal approximation theory described in Section 3 with three different bandwidth choices: $h_{n}\in\{0.8h_{n,HS},h_{n,HS},1.2h_{n,HS}\}$, where $h_{n,HS}$ denotes the Hall-Sheather bandwidth sequence. From Table 4, we find that the coverage probabilities results were not very sensitive across bandwidth values although the CI estimates were slightly under-sized. We also notice that the CI estimates based on taking $h_{n}=h_{n,HS}$ or $h_{n}=1.2h_{n,HS}$ performed quite well in terms of overall performance.

5 An illustrative empirical example: estimating the demand for fish

We illustrate usefulness of our method for exact computation of the IVQR GMM estimator in an empirical study of the demand for fish. We used the dataset constructed by Graddy (1995) on the transactions of whiting in the Fulton fish market in New York. The data were also previously studied in Chernozhukov and Hansen (2008) and Chernozhukov, Hansen, and Jansson (2009) to illustrate the econometric methods developed for quantile regression models with endogeneity. In what follows, we mainly focused on analyzing the results estimated by the MIQP approach and comparing them to the inverse QR estimation results.

The data consist of 111 observations on the price and quantity of whiting transactions aggregated by day. The outcome variable $Y$ is the logarithm of total amount of whitings sold on each day and the endogenous explanatory variable $D$ is the logarithm of the average daily price. The exogenous explanatory variables include the indicators ($Monday$, $Tuesday$, $Wednesday$ and $Thursday$) for days of the week. The instrumental variables are indicators ($Stormy$ and $Mixed$) for weather conditions at sea. These instruments capture the wave height and wind speed, which should affect the supplied quantity of fish and hence the price in the market but should not influence the demand for fish. See Graddy (1995, 2006) for further details on the operation of the Fulton fish market, and the data and variables used for this study.

Following Chernozhukov, Hansen, and Jansson (2009), we considered the simple demand equation

[TABLE]

for the estimation of $\theta_{1}$, the price elasticity of the demand, which may vary across the demand level $U$. We also augmented the specification (5.1) by incorporating the day effect variables as additional controls, and then performed the estimation. Table 5 presents the estimation results for $\theta_{1}\left(\tau\right)$ under these two different specifications. For GMM estimation results, we took $L=\left(1,Stormy,Mixed\right)$ as instruments and configurated the MIQP setting in the same fashion as in Section 4. We used the Gaussian kernel and the Hall-Sheather bandwidth choice for estimating the standard deviation of the GMM estimator and constructing the 95% CI for $\theta_{1}\left(\tau\right)$. We also performed some sensitivity check and found that the results were not very sensitive to the bandwidth choice. Moreover, we also extracted the inverse QR and the corresponding 95% asymptotic CI estimation results provided by Chernozhukov, Hansen, and Jansson (2009, Table

1. on the same estimating model specifications and listed them in Table 5 for comparison.

We now summarize the results in Table 5. First, we find that, for both model specifications, the point estimates of the demand elasticity were all negative but the magnitudes varied across quantile indices. Moreover, both the GMM and inverse QR estimates of $\theta_{1}(\tau)$ were of similar values under the basic specification (5.1). When the day effect variables were included as additional controls, the values of $\widehat{\theta}_{1}(\tau)$ across these two estimation methods differed to a larger extent in the case of $\tau=0.25$. Furthermore, we note that the CI results based on both the GMM and inverse QR methods indicate that the negativity of $\theta_{1}\left(\tau\right)$ was significant for $\tau\in\{0.25,0.75\}$ but we could not reject the case of $\theta_{1}\left(\tau\right)$ being zero at $\tau=0.5$.

6 Conclusions

In this paper, we have proposed a mixed integer quadratic programming approach for estimating the IVQR model within the GMM framework. One possible application of our approach is panel data quantile regression for group-level treatments (Chetverikov, Larsen, and Palmer, 2016). To deal with group-level unobservables, the estimation procedure in Chetverikov, Larsen, and Palmer (2016) consists of group-by-group quantile regression followed by two-stage least squares. They mention (in their footnote 10) that the latter step could be replaced by an IV median regression, if one is willing to replace the usual assumption that the group-level errors are uncorrelated with instruments with median uncorrelation (Komarova, Severini, and Tamer, 2012). This alternative step can be computed using our computation algorithm. It is an interesting topic for future research to fully develop this alternative to IV quantile regression for group-level treatments.

Our approach is limited to GMM estimators for parametric IVQR models. One may consider semiparametric models with endogeneity. For example, Chen, Linton, and Van Keilegom (2003) considered partially linear median regression with some endogenous regressors as one of their examples. Their proposed estimator consists of a two-step procedure: in the first step, nonparametric median regression is carried out given the parameter of interest and in the second step, GMM estimation is implemented with the first step estimates as inputs. Our proposed algorithm is not directly applicable because of the first nonparametric step. It is another interesting topic for future research to develop an algorithm to compute this kind of two-step semiparametric quantile IV estimators.