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.
91,825 characters · 23 sections · 101 citation commands
Decentralization Estimators for Instrumental Variable Quantile Regression Models
Quantile regression (QR), introduced by KoenkerBassett1978, is a widely-used method for estimating the effect of regressors on the whole outcome distribution. QR is flexible, easy to interpret, and can be computed very efficiently as the solution to a convex problem. However, in many applications, the variables of interest are endogenous, rendering QR inconsistent for estimating causal quantile effects. The instrumental variable quantile regression (IVQR) model of CH2005,CH2006 generalizes QR to accommodate endogenous regressors. Unfortunately, in sharp contrast to QR and other IV estimators such as two-stage least squares (2SLS), estimation of IVQR models is computationally challenging because the resulting estimation problem, formulated as a generalized method of moments (GMM) problem, is non-smooth and non-convex. From an applied perspective, this issue is particularly troublesome as resampling methods are often used to avoid the choice of tuning parameters when estimating the asymptotic variance of IVQR estimators.
In this paper, we develop a new class of estimators for linear IVQR models. The proposed estimators are fast, easy to implement, tuning-free, and do not require the availability of high-level “black box” optimization routines. They are particularly suitable for settings with many exogenous regressors, a moderate number of endogenous regressors, and a large number of observations, which are ubiquitous in applied research. The key insight underlying our estimators is that the complicated and nonlinear IVQR estimation problem can be “decentralized”, i.e., decomposed into a set of more tractable sub-problems, each of which is solved by a “player” who best responds to the other players' actions. Each subproblem is a conventional (weighted) QR problem, which is convex and can be solved very quickly using robust algorithms. The IVQR estimator is then characterized as a fixed point of such sub-problems, which can be viewed as the pure strategy Nash equilibrium of the “game”. Computationally, this reformulation allows us to recast the original non-smooth and non-convex optimization problem as the problem of finding the fixed point of a low dimensional map, which leads to substantial reductions in computation times.
Implementation of our preferred procedures is straightforward and only requires the availability of a routine for estimating quantile regressions and in some cases a univariate root-finder. The resulting estimation algorithms attain significant computational gains. For example, we show that in problems with two endogenous variables, a version of our estimator that uses a contraction algorithm is 149--308 times faster than the most popular existing approach for estimating IVQR models, the inverse quantile regression (IQR) estimator of CH2006. Another version that uses a nested root-finding algorithm, which is guaranteed to converge under a milder condition, is 84--134 times faster than the IQR estimator. Importantly, these computational gains do not come at a cost in terms of the finite sample performance of our procedures, which is very similar to IQR. The computational advantages of our estimators are even more substantial with more than two endogenous variables. The reason is that the dimensionality of the grid search underlying IQR corresponds to the number of endogenous variables, which renders IQR computationally prohibitive in empirically relevant settings whenever the number of endogenous variables exceeds two or three.
The fixed point reformulation also provides new insights into global identification of IVQR models. In particular, it allows us to study identification and stability of the algorithms (at the population level) in the same framework. Exploiting the equivalence of global identification and uniqueness of the fixed point, we give a new identification result and population algorithms based on the contraction mapping theorem. We then compare our identification conditions to those of CH2006. Further, our reformulation is shown to be useful beyond setups where the contraction mapping theorem applies as long as the parameter of interest is globally identified. For such settings, algorithms based on root-finding methods are proposed. Finally, we show that, by recursively nesting fixed point problems, it is possible to recast the IVQR estimation problem as a univariate root-finding problem. While adding nests incurs additional computational costs, our Monte Carlo experiments suggest that estimation procedures based on nesting perform well and are computationally reasonable when the number of endogenous variables is moderate.
We establish consistency and asymptotic normality of the proposed estimators and prove validity of the empirical bootstrap for estimating the limiting laws. The bootstrap is particularly attractive in conjunction with our computationally efficient estimation algorithms as it allows us to avoid the choice of tuning parameters inherent to estimating the asymptotic variance based on analytic formulas. The key technical ingredient for deriving our theoretical results is the Hadamard differentiability of the fixed point map. This result may be of independent interest.
To illustrate the usefulness of our estimation algorithms, we revisit the analysis of the impact of 401(k) plans on savings in CH2004. Based on this application, we perform extensive Monte Carlo simulations, which demonstrate that our estimation and inference procedures have excellent finite sample properties.
We contribute to the literature on estimation and inference based on linear IVQR models. ChernozhukovHong2003 have proposed a quasi-Bayesian approach which can accommodate multiple endogenous variables but, as noted by CH2013, requires careful tuning in applications. CH2006 have developed an inverse QR algorithm that combines grid search with convex QR problems. Because the dimensionality of the grid search equals the number of endogenous variables, this approach is computationally feasible only if the number of endogenous variables is very low. ChernozhukovHansen2008 and Jun2008 have studied weak instrument robust inference procedures based on the inversion of Anderson-Rubin-type tests. Chernozhukovetal2009 have proposed a finite sample inference approach. AndrewsMikusheva2016 have developed a general conditional inference approach and derived sufficient conditions for the IVQR model. KaplanSun2017 and deCastroetal2018 have suggested to use smoothed estimating equations to overcome the non-smoothness of the IVQR estimation problem, although the non-convexity remains. More recently, ChenLee2018 have proposed to reformulate the IVQR problem as a mixed-integer quadratic programming problem that can be solved using well-established algorithms. However, efficiently solving such a problem is still challenging even for low-dimensional settings. By replacing the $\ell_2$ norm by the $\ell_\infty$ norm, Zhu2018 has shown that the problem admits a reformulation as a mixed-integer linear programming problem, which can be computed more efficiently than the quadratic program in ChenLee2018. This procedure typically requires an early termination of the algorithm to ensure computational tractability which is akin to a tuning parameter choice. In addition, Zhu2018 has proposed a $k$-step approach that allows for estimating models with multiple endogenous regressors based on large datasets, but requires estimating the gradient. pouliot2018 proposes a mixed integer linear programming formulation that allows for subvector inference via the inversion of a distribution-free rankscore test and can be modified to accommodate weak instruments. An important drawback of the estimation approaches based on mixed integer reformulations is that they rely on the availability of high-level “black box” optimization routines such as Gurobi and often require careful tuning in applications. Finally, imposing a location-scale model for the potential outcomes, Machado2018 propose moment-based estimators for the structural quantile function.
Compared to the existing literature on the estimation of linear IVQR models, the main advantages of the proposed estimation algorithms are the following. First, by relying on convex QR problems, our estimators are easy to implement, robust, and computationally efficient in settings with many exogenous variables, a moderate number of endogenous variables, and a large number of observations. Second, by exploiting the specific structure of the IVQR estimation problem, our estimators are tuning-free and do not require the availability of high-level “black box” optimization routines. Third, our estimators are based on the original IVQR estimation problem and thus avoid the choice of smoothing bandwidths and do not rely on additional restrictions on the structural quantile function.
Semi- and nonparametric estimation of IVQR models has been studied by ChernozhukovImbensNewey2007, HorowitzLee2007, ChenPouzo2009, ChenPouzo2012, GagliardiniScaillet2012, and Wuthrich2019. CH2013 and Chernozhukov+17handbook have provided surveys on the IVQR model including references to empirical applications.
AAI2002 have proposed an alternative approach to the identification and estimation of quantile effects with binary endogenous regressors, which builds on the local average treatment effects framework of AngristImbens1994. Their approach has been extended and further developed by Frandsenetal2012, FrolichMelly2013jbes, and Bellonietal2017 among others. We refer to MellyWuthrich2017 for a recent review of this approach and to Wuthrich2020 for a comparison between this approach and the IVQR model. Identification and estimation in nonseparable models with continuous endogenous regressors have been studied by Chesher2003, KoenkerMa2006, Lee2007, Jun2009, ImbensNewey2009, DHaultfoeuilleFevrier2015, and Torgovitsky2015 among others.
On a broader level, our paper contributes to the literature which proposes estimation procedures that rely on decomposing computationally burdensome estimation problems into several more tractable subproblems. This type of procedure, which we call decentralization, has been applied in different contexts. Examples include the estimation of single index models with unknown link function WeisbergWelsh94, general maximum likelihood problems Smyth96, linear models with high-dimensional fixed effects GuimaraesPortugal10, sample selection models MarraRadice13, peer effects models Arcidiacono+12, interactive fixed effects models Chen+14panel,MoonWeidner15, and random coefficient logit demand models LeeSeo15ablp. Most of these papers decompose a single estimation problem into two subproblems. The present paper explicitly considers cases in which the number of subproblems may exceed two.
The remainder of the paper is structured as follows. Section (ref) introduces the setup and the IVQR model. Section (ref) shows that the IVQR estimation problem can be decentralized into a series of (weighted) conventional QR problems. In Section (ref), we introduce population algorithms based on the contraction mapping theorem and root-finders. Section (ref) discusses the corresponding sample algorithms. In Section (ref), we establish the asymptotic normality of our estimators and prove the validity of the bootstrap. Section (ref) presents an empirical application. In Section (ref), we provide simulation evidence on the computational performance and the finite sample properties of our methods. Section (ref) concludes. All proofs as well as some additional theoretical and simulation results are collected in the appendix.
Consider a setup with a continuous outcome variable $Y$, a $d_X\times 1$ vector of exogenous covariates $X$, a $d_D\times 1$ vector of endogenous treatment variables $D$, and a $d_Z\times 1$ vector of instruments $Z$. The IVQR model is developed within the standard potential outcomes framework Rubin1974. Let $\{Y_d\}$ denote the (latent) potential outcomes. The object of primary interest is the conditional quantile function of the potential outcomes, which we denote by $q(d,x,\tau)$. Having conditioned on covariates $X=x$, by the Skorokhod representation of random variables, potential outcomes can be represented as
This representation lies at the heart of the IVQR model. With this notation at hand, we state the main conditions of the IVQR model CH2005.
We briefly discuss the most important aspects of Assumption (ref) and refer the interested reader to CH2005,CH2006,CH2013 for more comprehensive treatments. Assumption (ref).(ref) states the Skorohod representation of $Y_d$ and requires strict monotonicity of the potential outcome quantile function, which rules out discrete outcomes. Assumption (ref).(ref) imposes independence between the potential outcomes and the instrument. Assumption (ref).(ref) defines a general selection mechanism. The key restriction of the IVQR model is Assumption (ref).(ref). Rank invariance (a) requires individual ranks $U_d$ to be the same across treatment states. Rank similarity (b) weakens this condition, allowing for random slippages of $U_d$ away from a common level $U$. Finally, Assumption (ref).(ref) summarizes the observables.
The main implication of Assumption (ref) is the following conditional moment restriction CH2005:
In this paper, we focus on the commonly used linear-in-parameters model for $q(\cdot)$ CH2006:
where $\theta(\tau):=(\theta_X(\tau)',\theta_D(\tau)')'\in \mathbb{R}^{d_X+d_D}$ is the finite dimensional parameter vector of interest. The conditional moment restriction (ref) suggests GMM estimators based on the following unconditional population moment conditions:
Our primary goal here is to obtain estimators in a computationally efficient and reliable manner. We therefore focus on just-identified moment restrictions where $d_Z=d_D$, for which the construction of an estimator is straightforward. A potential caveat of this approach is that estimators based on these restrictions do not achieve the pointwise (in $\tau$) semiparametric efficiency bound implied by the conditional moment restrictions (ref). Appendix (ref) provides a discussion of overidentified GMM problems and presents a two-step approach for constructing efficient estimators based on the proposed algorithms.
In what follows, we will often suppress the dependence on $\tau$ to lighten-up the exposition. We then define the true parameter value $\theta^*$ as the solution to the moment conditions, i.e., \[ \Psi_P\left(\theta^*\right)=0. \] The resulting GMM objective function reads
where $m_i\left(\theta\right):=\left( 1\left\{ Y_i\leq X_i^\prime \theta_X + D_i'\theta_D \right\} -\tau\right)(X_i',Z_i')'$ and $W_N\left(\theta\right)$ is a positive definite weighting matrix. Estimation based on (ref) is complicated by the non-smoothness and, most importantly, the non-convexity of $\mathcal{Q}_N^{GMM}$. This paper proposes a new set of algorithms to address these challenges.
Here we describe the basic idea behind our decentralization estimators. To simplify the exposition, we first illustrate our approach with the population problem of finding the true parameter value $\theta^*$ in the IVQR model. Our estimator then adopts the analogy principle, which will be presented in Section (ref). The key insight is that the complicated nonlinear IVQR estimation problem can be “decentralized”, i.e., decomposed into a set of more tractable sub-problems, each of which is solved by a “player” who best responds to other players' actions. Specifically, we first split the parameter vector $\theta$ into $J$ subvectors $\theta_1,\dots,\theta_{J}$, where $J=d_D+1$. We then decompose the grand estimation problem into $J$ subproblems. Each of the subproblems is allocated to a distinct player. For each $j$, player $j$'s choice variable is the $j$-th subvector $\theta_j$. Her problem is to find the value of $\theta_j$ such that a subset of the moment restrictions is satisfied given the other players' actions $\theta_{-j}$, where $\theta_{-j}$ stacks the components of $\theta$ other than $\theta_j$. This reformulation allows us to view the estimation problem as a game of complete information and to characterize $\theta^*$ as the game's pure strategy Nash equilibrium.
We start by describing the parameter subvectors. First, let $\theta_1=\theta_X\in\mathbb R^{d_X}$ denote the vector of coefficients on the exogenous variables. We allocate this subvector to player 1. Similarly, for each $j=2,\dots,J$, let $\theta_j\in\mathbb R$ denote the coefficient on the $(j-1)$-th endogenous variable, which is allocated to player $j$. The coefficient vector for the endogenous variable can therefore be written as $\theta_D=(\theta_2,\dots,\theta_J)'$. For each $\theta\in\mathbb{R}^{d_X+d_D}$, define the following (weighted) QR objective functions:
where $\rho_\tau(u)=u(\tau-1\{u<0\}$) is the “check-function”. We assume that the model is parametrized such that $Z_\ell/D_\ell$ is positive for all $\ell=1,\dots,d_D$. Under our assumptions, we can always reparametrize the model such that this condition is met; see Appendix (ref) for more details.
The players then solve the following optimization problems:
Observe that each player's problem is a weighted QR problem, which is convex in its choice variable. For the sample analogues of these problems fast solution algorithms exist Koenker2017.
For each $j$, let $\tilde L_j(\theta_{-j})$ denote the set of minimizers. Borrowing the terminology from game theory, we refer to these maps as best response (BR) maps. Under an assumption we specify below, the first-order optimality conditions imply that, for each $j$, any element $\tilde\theta^*_j\in \tilde L_j(\theta_{-j})$ of the BR map satisfies
where $D_{-(j-1)}$ stacks as a vector all endogenous variables except $D_{j-1}$. Note that $\Psi_{P}\left(\theta\right)=\left(\Psi_{P,1}(\theta)',\dots, \Psi_{P,J}(\theta)\right)'$ is the set of unconditional IVQR moment conditions. Hence, $\theta^*$ satisfies
which implies that $\theta^*$ is a fixed point of the BR-maps (i.e.\ a Nash equilibrium of the game).
In what follows, we work with conditions that ensure the existence of singleton-valued BR maps $L_j$, $j=1,\dots,J$, such that, for each $j$, $\Psi_{P,j}\left(L_j(\theta_{-j}),\theta_{-j}\right)=0$.\footnote{While it may be interesting to work with set-valued maps, the existence of the BR functions greatly simplifies our analysis of identification and inference.} We say that the IVQR estimation problem admits decentralization if there exist such BR functions defined over domains for which the moment conditions can be evaluated.\footnote{In Appendix (ref), we also provide weaker conditions under which the decentralization holds on a local neighborhood of $\theta^*$. We call such a result local decentralization, which is sufficient for analyzing the (local) asymptotic behavior of the estimator.} To ensure decentralization, we make the following assumption.
For each $j$, let $\Theta_{-j}\subset \mathbb{R}^{d_{-j}}$ denote the parameter space for $\theta_{-j}$. Assumption (ref).(ref) ensures that $\Theta$ is compact. This assumption also ensures that each $\Theta_{-j}$ is a closed rectangle, which we use to show that $L_j$ is well-defined on a suitable domain. Assumptions (ref).(ref) and (ref).(ref) impose standard regularity conditions on the conditional density and the moments of the variables in the model. We assume $D_\ell$ has a compact support, which allows us to always reparameterize the model so that the objective function in (ref) is well-defined and convex (cf.\ Appendix (ref)). The first part of Assumption (ref).(ref) is a standard full rank condition which is a natural extension of the local full rank condition required for local identification and decentralization (cf.\ Assumption (ref) in the appendix). For the second part of Assumption (ref).(ref), it suffices that the model is parametrized such that, for each $\ell\in\{1,\dots,d_D\}$, $D_\ell Z_\ell$ (and $Z_\ell/D_\ell$) is positive with probability one.
For each $j$, define
This is the set of subvectors $\theta_{-j}$ for which one can find $\theta_j\in\Theta_j$ such that $\theta=(\theta_j,\theta_{-j})'$ solve the $j$-th moment restriction. We take this set as the domain of player $j$'s best response function $L_j$.
The following lemma establishes that the IVQR model admits decentralization.
We now introduce maps that represent all players' (joint) best responses. We consider two basic choices of such maps; one represents simultaneous responses and the other represents sequential responses. In what follows, for any subset $a\subset\{1,\dots,J\}$, let $\theta_{-a}$ denote the subvector of $\theta$ that stacks the components of $\theta_j$'s for all $j\notin a$. If $a$ is a singleton (i.e. $a=\{j\}$ for some $j$), we simply write $\theta_{-j}$. For each $j$ and $a\subseteq\{1,\dots,J\}\setminus\{j\}$, let $\pi_{-a}:\Theta_{-j}\to \prod_{k\in \{1,\dots,J\}\setminus(\{j\}\cup a)}\Theta_k$ be the coordinate projection of $\theta_{-j}$ to a further subvector that stacks all components of $\theta_{-j}$ except for those of $\theta_k$ with $k\in a$.
Let $D_K:=\{\theta\in\Theta: \pi_{-j}\theta\in R_{-j},~ j=1,\dots,J\}$. Let $K:D_K\to\mathbb{R}^{d_X+d_D}$ be a map defined by
This can be interpreted as the players' simultaneous best responses to the initial strategy $(\theta_1,\dots,\theta_{J})$. With one endogenous variable, this map simplifies to
Here, $K$ maps $\theta=(\theta_1,\theta_2)$ to a new parameter value through the simultaneous best responses of players 1 and 2.
Similarly, let $D_M\subset \mathbb R^{d_D}$ and let $M:D_M\to\mathbb R^{d_D}$ be a map such that
which can be interpreted as the players' sequential responses (first by player 1, then player 2, etc.) to an initial strategy $\theta_{-1}=(\theta_{2},\dots,\theta_{J})$.\footnote{One may define $M$ by changing the order of responses as well. For theoretical analysis, it suffices to consider only one of them. Once the fixed point $\theta^*_{-1}$ of $M$ is found, one may also obtain $\theta^*_1$ using $\theta^*_1=L_1(\theta^*_{-1})$.} Note that the argument of $M$ is not the entire parameter vector. Rather, it is a subvector of $\theta$ consisting of the coefficients on the endogenous variables. In order to find a fixed point, this feature is particularly attractive when the number of endogenous variables is small. With one endogenous variable (i.e. $\theta_2\in \mathbb R$ is a scalar), the map simplifies to
which is a univariate function whose fixed point is often straightforward to compute.
Define
This is the set on which the map $\theta_{-1}\to L_2\left( L_1\left(\theta_{-1}\right),\pi_{-\{1,2\}}\theta_{-1}\right)$, the first component of $M$, is well-defined. We then recursively define $\tilde R_j$ for $j=2,\dots,d_D$ in a similar manner. A precise definition of these sets is given in Appendix (ref). Now define
where the second equality follows because $\tilde R_{d_D}$ turns out to be a subset of $\tilde R_j$ for all $j\le d_D$. The following corollary ensures that $K$ and $M$ are well-defined on $D_K$ and $D_M$ respectively.
The key insight that we exploit is that, by construction of the BR maps, the problem of finding a solution to $\Psi_P(\theta)=0$ is equivalent to the problem of finding a fixed point of $K$ (or $M$). The following proposition states the formal result.
In view of Proposition (ref), the original IVQR estimation problem can be reformulated as the problem of finding the fixed point of $K$ (or $M$). This reformulation naturally leads to discrete dynamical systems associated with these maps, which in turn provide straightforward iterative algorithms for computing $\theta^*$.
These discrete dynamical systems constitute the basis for our estimation algorithms.\footnote{These discrete dynamical systems can also be viewed as learning dynamics in a game LiBasar87,FudenbergLevine07.}
In this section, we explore the implications of the fixed point reformulation for constructing population-level algorithms for computing fixed points.
We first consider conditions under which $K$ and $M$ are contraction mappings. They ensure that the discrete dynamical systems induced by $K$ and $M$ are convergent to unique fixed points. Moreover, in view of Proposition (ref), (point) identification is equivalent to the uniqueness of the fixed point of $K$ (or $M$). Therefore, the conditions we provide below are also sufficient for the point identification of $\theta^*.$ We will discuss the relationship between our conditions and existing ones in the next section.
For any vector-valued map $E$, let $J_E(x)$ denote its Jacobian matrix evaluated at its argument $x$. For any matrix $A$, let $\|A\|$ denote its operator norm. We provide conditions in terms of the Jacobian matrices of $K$ and $M$, which are well-defined by Corollary (ref).
Under this additional assumption, the iterative algorithms are guaranteed to converge to the fixed point. We summarize this result below.
In the case of a single endogenous variable, the Jacobian matrices of $K$ and $M$ are given by
where
One may therefore check the high-level condition through the Jacobians of the original moment restrictions. In Appendix (ref), we illustrate a simple primitive condition for a local version of Assumption (ref). In practice, we found that violations of Assumption (ref) lead to explosive behavior of the estimation algorithms and, thus, are very easy to detect numerically.
In view of Proposition (ref), identification of $\theta^*$ is equivalent to uniqueness of the fixed points of $K$ and $M$, which is ensured by Proposition (ref). Here, we discuss how the conditions required by Proposition (ref) relate to the ones in the literature.
We start with local identification. The parameter vector $\theta^*$ is said to be locally identified if there is a neighborhood $\mathcal{N}$ of $\theta^*$ such that $\Psi_P(\theta)\ne 0$ for all $\theta\ne \theta^*$ in the neighborhood. Local identification in the IVQR model follows from standard results Rothenberg71,Chen+14. For example, if $\Psi_P(\theta)$ is differentiable, Chen+14 show that full rank of $J_{\Psi_P}(\theta)$ at $\theta^*$ is sufficient for local identification.
It is interesting to compare this full rank condition to Assumption (ref).(ref) in the appendix, which is a local version of Assumption (ref).(ref). Assumption (ref).(ref) requires that $\rho\left(J_K(\theta^*) \right)<1$, where $\rho(A)$ denotes the spectral radius of a square matrix $A$. We highlight the connection in the case with a single endogenous variable. Full rank of $J_{\Psi_P}(\theta^*)$ is equivalent to $\det\left(J_{\Psi_P}(\theta^*)\right)\ne 0$. Observe that, for any $\theta$,
If $\partial\Psi_{P,j}(\theta)/\partial\theta_j'|_{\theta=\theta^*} $ is invertible for $j=1,2$ (which is true under Assumption (ref).(ref)), $J_{\Psi_P}(\theta^*)$ is full rank if and only if
That is, it requires that none of the eigenvalues of $J_K$ has modulus one. Therefore, Assumption (ref).(ref) is sufficient but not necessary for condition (ref) to hold. Specifically, Assumption (ref).(ref) requires all eigenvalues of $J_{K}(\theta^*)$ to lie strictly within the unit circle, while local identification only requires all eigenvalues not to be on the unit circle. In terms of the dynamical system induced by $K$, the former ensures that the dynamical system has a unique asymptotically stable fixed point, while the latter ensures that the system has a unique hyperbolic fixed point, which is a more general class of fixed points galor2007discrete.\footnote{The argument above also applies to settings with multiple endogenous variables. A similar result can also be shown for $M$.} Under the former condition, iteratively applying the contraction map induces convergence, while the latter generally requires a root finding method to obtain the fixed point.
Now we turn to global identification and compare Proposition (ref) to the global identification result in CH2006.
Under Conditions (i)--(iv), which are substantially stronger than the local identification conditions discussed above, the result in Lemma (ref) follows from an application of Hadamard's global univalence theorem (e.g. Theorem 1.8 in ambrosetti1995primer).
Comparing Lemma (ref) to Proposition (ref), we can see that the result in Lemma (ref) establishes identification over the whole parameter space $\Theta$, while Proposition (ref) establishes identification over the sets $\tilde{D}_K$ and $\tilde{D}_M$, which will generally be subsets of $\Theta$. Regarding the underlying assumptions, Conditions (i) and (ii) in Lemma (ref) correspond to our Assumptions (ref).(ref) and (ref).(ref). Moreover, our Assumption (ref).(ref) constitutes an easy-to-interpret sufficient condition for continuity of $J_{\Psi_P}$ as required in Condition (iii). To apply Hadamard's global univalence theorem, CH2006 assume the simple connectedness of the image of $\Psi_P$ (Condition (iv)). By contrast, we use a different univalence theorem by gale1965jacobian (applied to the map $\Xi$ defined in (ref) that arises from each subsystem), which does not require further conditions. However, when establishing global identification based on the contraction mapping theorem, we need to impose an additional condition on the Jacobian (Assumption (ref)). In sum, our conditions are somewhat stronger in terms of restrictions on the Jacobian, but they are relatively easy to check and allow us to dispense with an abstract condition (simple connectedness of the image of a certain map) to apply a global univalence theorem.
Assumption (ref) is a sufficient condition for the uniqueness of the fixed point and the convergence of the contraction-based algorithms. Even in settings where this assumption fails to hold, one may still identify $\theta^*$ and design an algorithm that is able to find it under weaker conditions on the Jacobian. This is the case under the assumptions in the general (global) identification result of CH2006; see Lemma (ref).
Note that, for the simultaneous dynamical system, $\theta^*$ solves
where $I_{{d_X+d_D}}$ is the identity map. Similarly, in the sequential dynamical system, $\theta^*_{-1}$ solves
Therefore, standard root-finding algorithms can be used to compute the fixed point.
For implementing root-finding algorithms, we find that reducing the dimension of the fixed point problem is often helpful. Toward this end, we briefly discuss another class of dynamical systems and associated population algorithms which can be used for the purpose of dimension reduction. Namely, with more than two players, one can construct nested dynamical systems, which induce nested fixed point algorithms. Nesting is useful as it allows for transforming any setup with more than two players into a two-player system.
To fix ideas, consider the case of three players ($J=3$). Fix player 3's action $\theta_3\in\Theta_3\subset \mathbb R$ and consider the associated “sub-game” between players 1 and 2. To describe the sub-game, define $M_{1,2|3}(\cdot\mid \theta_3):\Theta_2\to\Theta_2$ pointwise by
This map gives the sequential best responses of players 1 and 2, while taking player 3's strategy given. Define the fixed point $L_{12}:\Theta_3\to \Theta_1\times\Theta_2$ of the sub-game by
This map then defines a new “best response” map. Here, given $\theta_3$, the players in the sub-game (i.e., players 1 and 2) collectively respond by choosing the Nash equilibrium of the sub-game. The overall dynamical system induced by the nested decentralization is then given by
Hence, we can interpret the nested algorithm as a two-player dynamical system where one player solves an internal fixed point problem. The nesting algorithms require existence and uniqueness of the fixed points in the sub-game between players 1 and 2. In Appendix (ref), we discuss the formal conditions required for the existence and uniqueness of such fixed points.\footnote{For an equilibrium of the sub-game to be well-defined, one may directly assume that CH2006's global identification condition holds within the sub-game, given player 3's action $\theta_3$. Alternatively, if Assumption (ref) holds, we fix $\theta_3$ to a value such that $(\theta_2,\theta_3)\in \tilde D_M$ for some $\theta_2\in\Theta_2$; see Appendix (ref).}
This nesting procedure is generic and can be extended to more than three players by sequentially adding additional layers of nesting.\footnote{In the current example, consider adding player 4 and letting players 1-3 best respond by returning the fixed point of the sub-game through $M_3$ given $\theta_4$. One can repeat this for additional players. This procedure can also be applied to the simultaneous dynamical system induced by $K$.} It follows that any decentralized estimation problem with more than two players can be reformulated as a nested dynamical system with two players: player $J$ and all others $-J$. The resulting dynamical system $M_J(\theta_J)=L_{J}\left(L_{-J}(\theta_J)\right)$ is particularly useful when $M_J$ is not necessarily a contraction map since $\theta_J$ is a scalar such that, as we see below, its fixed point can be computed using univariate root-finding algorithms.
Let $\left\{(Y_i,D_i',X_i',Z_i')\right\}_{i=1}^N$ be a sample generated from the IVQR model. Our estimators are constructed using the analogy principle. For this, define the sample payoff functions for the players as
For each $j=1,\dots, J$, let the sample BR function $\hatL_j(\theta_{-j})$ be a function such that
Assuming that the model is parametrized in such a way that $Z_{\ell,i}/D_{\ell,i}$, $\ell=1,\dots,d_D$, is positive, these are convex (weighted) QR problems for which fast solution algorithms exist. In our empirical applications and simulations, we use the R-package quantreg to estimate the QRs quantreg2018. For example, $\hat L_2(\theta_{-2})$ can be computed by running a QR with weights $Z_{1,i}/D_{1,i}$ in which one regresses $Y_i-X_i'\theta_1-D_{2,i}\theta_3-\dots-D_{d_D,i}\theta_J$ on $D_{1,i}$ without a constant. These sample BR functions also approximately solve the sample analog of the moment restrictions in (ref)--(ref); see Lemma (ref) in the appendix.
We focus on estimators constructed based on the dynamical system $M$. Lemma (ref) in the appendix shows that the estimators based on $M$ and $K$ are asymptotically equivalent. However, in our simulations, we found that, while the convergence properties of contraction algorithms based on $M$ are excellent, the convergence properties of contraction algorithms based on $K$ can be sensitive to the choice of starting value. This may be attributed to the fact that the algorithm based on $K$ requires all $d_X+d_D$ components of the starting value to be in the domain of the contraction map, while the algorithm based on $M$ only requires that the same condition is satisfied by the starting value for $\theta_{-1}$. Our simulations suggest that it is often easier to satisfy this requirement with $M$ since the number of components in $\theta_{-1}$, $d_D$, is typically much smaller than $d_X+d_D$. Also, for root-finding algorithms, we prefer the sequential dynamical system (induced by $M$) because it again leads to a substantial dimension reduction: it reduces the original $(d_X+d_D)$-dimensional GMM estimation problem to a $d_D$-dimensional root-finding problem.
We construct estimation algorithms by mimicking the population algorithms. Let $\hat M$ denote a sample analog of $M$:
where $\theta_1=\hatL_1\left(\theta_{-1}\right)$. This map induces a sample analog of the sequential dynamical system in Section (ref).
The first set of algorithms exploits that, under Assumption (ref), $\hat M$ is a contraction mapping with probability approaching one. In this case, we iterate the dynamical system or (ref) until $\|\theta_{-1}^{(s)}-\hat M(\theta_{-1}^{(s)})\|$ is within a numerical tolerance $e_N$.\footnote{In the next section, we require $e_N=o(N^{-1/2})$, which ensures that the numerical error does not affect the asymptotic distribution.} This iterative algorithm is known to converge at least linearly. The approximate sample fixed point $\hat\theta_N=(\hat\theta_{N,1},\hat\theta_{N,-1})$ that meets the convergence criterion then serves as an estimator for $\theta$.
We construct an estimator $\hat\theta_N$ of $\theta^*$ as an approximate fixed point to the sample problem:
where $\hat\theta_{N,1}=\hatL_1(\hat \theta_{N,-1} )$ and $e_N$ is a numerical tolerance. This problem can be solved efficiently using well-established root-finding algorithms since $\hat M$ is easy to evaluate as the composition of standard QRs. When $d_D=1$, one may use Brent's method Brent1971 whose convergence is superlinear. When $d_D>1$, one could apply the Newton-Raphson method, which achieves quadratic convergence but requires an estimate or a finite difference approximation of the derivative. The corresponding approximation error may affect the performance. Alternatively, one can compute the fixed point by minimizing $\| \hat M(\theta) - \theta \|^2$. The potential issue with this approach is that translating the root-finding problem into a minimization problem can lead to local minima in the objective function. Therefore, it is important to use global optimization strategies.
As described in Section (ref), nesting can be used to reduce the dimensionality even further. In particular, the problem can be reformulated as a one-dimensional fixed point problem, which can be solved using existing methods. We found that Brent's method performs very well in our context. Nesting is suitable when the number of endogenous variables is moderate. While adding nests incurs additional computational costs, our Monte Carlo experiments suggest that they are not excessive when the number of endogenous variables is moderate.\footnote{See Section (ref) and Appendix (ref).}
A key insight underlying our estimation algorithms is that, given $\theta_{-1}$, the estimation problem becomes a standard convex QR problem: \[ \hatL_1\left(\theta_{-1} \right) \in \arg\min_{\tilde\theta_1 \in \mathbb{R}^{d_X}}Q_{N,1}(\tilde\theta_1,\theta_{-1}) \] This insight suggests a profiling estimator.\footnote{We thank an anonymous referee for suggesting this alternative estimator.} In particular, $\hat\theta_{N,-1}$ can be obtained as the approximate root of the function
and $\hat\theta_{N,1}$ can be estimated as $\hat\theta_{N,1}\in \arg\min_{\tilde\theta_1 \in \mathbb{R}^{d_X}}Q_{N,1}(\tilde\theta_1,\hat\theta_{N,-1})$. When there is only one endogenous variable, $f$ is scalar-valued and univariate root-finders can be used. With multivariate endogenous variables, one can either directly apply multivariate root-finders or construct nested algorithms as described in Sections (ref) and (ref).
Relative to the root-finding methods described in Section (ref), the profiling estimator has the advantage that evaluating $f$ only requires estimating one QR. On the other hand, the root-finding algorithms in Section (ref) efficiently exploit the convexity of the subproblems for players $j=2,\dots,J$, and demonstrate a better computational performance than the profiling estimators with multiple endogenous variables (cf.\ Table (ref)).
We define an estimator $\hat\theta_N$ of $\theta^*$ as an approximate fixed point of $\hat M$ in the following sense:
An estimator of $\theta^*_1$ can be constructed by setting
In what follows, we call $\hat\theta_N=(\hat\theta_{N,1},\hat\theta_{N,-1})$ the fixed point estimator of $\theta^*.$ Under the conditions we introduce below, $\hat\theta_N$ is also an approximate fixed point of $\hat K$; see Lemma (ref) in the appendix. This turns out to be useful for stating theoretical results in a concise manner. While we mostly focus on algorithms based on $\hat M$ below, some of our theoretical results will be stated using $K$. $\hat M$ (or $\hat K$) is defined similarly for the nested dynamical system in which one player solves a fixed point problem in a sub-game.
Consistency and parametric convergence rates of $\hat\theta_N$ can be established using existing results. When $\hat M$ is asymptotically a contraction map, one may construct an estimator $\hat\theta_N$ satisfying (ref)--(ref) using the contraction algorithm in Section (ref) with tolerance $e_N=o(N^{-1/2}).$ One may then apply the result of Dominitz:2005aa to obtain the root-$N$ consistency of the estimator.\footnote{Satisfying $e_N=o(N^{-1/2})$ requires the number of iterations to increase as the sample size tends to infinity, which in turn satisfies requirement (ii) in Theorem 2 of Dominitz:2005aa.} For completeness, this result is summarized in Appendix (ref).
More generally, if $\hat M$ is not guaranteed to be a contraction, one may use root-finding algorithms that solve $\theta_{-1}-\hat M(\theta_{-1})=0$ up to approximation errors of $o(N^{-1/2})$. The root-$N$ consistency of $\theta_{N,-1}$ then follows from the standard argument for extremum estimators, in which we take $\mathcal Q_N(\theta_{-1})=\|\theta_{-1}-\hat M(\theta_{-1})\|$ as a criterion function.\footnote{The key conditions for these results, uniform convergence (in probability) of $\hat K$ and its stochastic equicontinuity, are established in Lemma (ref).} Since these results are standard, we omit details and focus below on the asymptotic distribution and bootstrap validity of the fixed point estimators. Our contributions are two-fold. First, we establish the asymptotic distribution of the fixed point estimator without assuming that $\hat M$ is an asymptotic contraction map, which therefore allows the practitioner to conduct inference using the estimator based on the general root-finding algorithm and complements the result of Dominitz:2005aa. Second, to our knowledge, the bootstrap validity of the fixed point estimators is new. These results are established by showing that, under regularity conditions, the population fixed point is Hadamard-differentiable and hence admits the use of the functional $\delta$-method, which may be of independent theoretical interest.
The following theorem gives the limiting distribution of our estimator. For each $w=(y,d',x',z')'$ and $\theta\in\Theta,$ let $f(w;\theta)\in \mathbb R^{d_X+d_D}$ be a vector whose sub-vectors are given by
and let $g(w;\theta)=(g_1(w;\theta)',\dots,g_J(w;\theta))'$ be a vector such that
In the theorem above, the asymptotic variance of $\hat\theta_N$ is characterized in terms of the Jacobian of $K$ and a Gaussian process $\mathbb W$. This can be reformulated to show its asymptotic relationship to a GMM estimator. Let $\tilde\theta_N$ be an estimator that solves the following estimating equations
Let $\Psi(\tau)=(X',Z')'.$ As shown in CH2006 (Theorem 3 and Remark 3), $\sqrt N(\tilde \theta_N-\theta^*)$ converges weakly to a mean zero multivariate normal distribution with variance
where $J_{\Psi_P}(\theta^*)=E\left[f_{\varepsilon(\tau)|X,D,Z}(0)\Psi(\tau) (X',D')\right],$ $\varepsilon(\tau)=Y-X'\theta^*_1-D'\theta^*_{-1}$, and $f_{\varepsilon(\tau)|X,D,Z}$ is $\varepsilon(\tau)$'s conditional density given $(X,D,Z)$. The following corollary shows that the fixed point estimator $\hat\theta_N$ is asymptotically equivalent to $\tilde\theta_N$ in terms of its limiting distribution.
To conduct inference on $\theta^*$, one may employ a natural bootstrap procedure. For this, use in (ref)--(ref) the bootstrap sample instead of the original sample to define the bootstrap analogs $\hat M^*$ and $\hat\theta^*_{N,-1}$ of $\hat M$ and $\hat\theta_{N,-1}$. In practice, the bootstrap can be implemented using the following steps.
The bootstrap is particularly attractive in conjunction with our new and computationally efficient estimation algorithms. By contrast, directly bootstrapping for instance the IQR estimator of CH2006 is computationally very costly. Alternative methods (either an asymptotic approximation or a score-based bootstrap) require estimation of the influence function, which involves nonparametric estimation of a certain conditional density. Directly bootstrapping our fixed point estimators avoids the use of any smoothing and tuning parameters.\footnote{The use of the bootstrap here is for consistently estimating the law of the estimator. Whether one may obtain higher-order refinements through a version of the bootstrap, e.g., the $m$ out of $n$ bootstrap with extrapolation Sakov:2000aa, is an interesting question which we leave for future research.}
The following theorem establishes the consistency of the bootstrap procedure. For this, let $\stackrel{L^*}{\leadsto}$ denote the weak convergence of the bootstrap law in outer probability, conditional on the sample path $\{W_i\}_{i=1}^\infty.$
In this section, we illustrate the proposed estimators by reanalyzing the effect of 401(k) plans on savings behavior as in CH2004. This empirical example constitutes the basis for our Monte Carlo simulations in Section (ref). As explained by CH2004, 401(k) plans are tax-deferred savings options that allow for deducting contributions from taxable income and accruing tax-free interest. These plans are provided by employers and were introduced in the United States in the early 1980s to increase individual savings. To estimate the effect of 401(k) plans ($D$) on accumulated assets ($Y$), one has to deal with the potential endogeneity of the actual participation status. CH2004 propose an instrumental variables approach to overcome this problem. They use 401(k) eligibility as an instrument ($Z$) for the participation in 401(k) plans. The argument behind this strategy, which is due to Poterbaetal1994,Poterbaetal1995,Poterbaetal1998 and Benjamin2003, is that eligibility is exogenous after conditioning on income and other observable factors. We use the same identification strategy here but note that there are also papers which argue that 401(k) eligibility is not conditionally exogenous Engenetal1996.
We use the same dataset as in CH2004. The dataset contains information about 9913 observations from a sample of households from the 1991 Survey of Income and Program Participation.\footnote{The dataset analyzed by CH2004 has 9,915 observations. Here we delete the two observations with negative income.} We refer to CH2004 for more information about the data and to their Tables 1 and 2 for descriptive statistics. Here we focus on net financial assets as our outcome of interest.\footnote{CH2004 also consider total wealth and net non-financial assets.}
We consider the following linear model for the conditional potential outcome quantiles
The vector of covariates $X$ includes seven dummies for income categories, five dummies for age categories, family size, four dummies for education categories, indicators for marital status, two-earner status, defined benefit pension status, individual retirement account participation status and homeownership, and a constant. Because $P(D=0)>0$, we re-parametrize the model by replacing $D$ by $D^{\star}=D+1$ to ensure that $Z/D^{\star}$ is well-defined and positive.
Below, we briefly describe the construction of the sequential response map for this application. First, we allocate $\theta_X$ to player 1 and hence denote this subvector by $\theta_1$. Similarly, we allocate $\theta_D$ to player 2 and denote it by $\theta_2$. For each $\tau$, the sequential response map can be constructed by taking the following steps.
Combining these two steps yields the sequential response map $\hat M(\cdot)=\hat L_2(\hat L_1(\cdot))$. Figure (ref) shows, for each $\tau \in \{0.25,0.50,0.75\}$, the graph of $\theta_2\mapsto\hat{M}(\theta_2)$. For each $\tau$, the intersection between $\hat{M}$ and the 45-degree line (i.e.\ the identity map) is our fixed point estimator $\hat\theta_D(\tau)$. Figure (ref) further provides a straightforward graphical way to check the validity of the sample analog of Assumption (ref). We can see that the sample analog of $J_M$ (i.e. the slope of $M$) is smaller than one at any $\theta_2$. This suggests that the contraction-based sequential algorithm converges at all three quantile levels, which is indeed what we find.
For each $\tau$, the steps for numerically calculating the fixed point estimator are as follows.
We compare our estimators to the IQR estimator of CH2006 based on a grid search over 500 points. IQR provides a very robust benchmark. However, due to the use of grid search, it is computationally expensive. In our empirical Monte Carlo study in Section (ref), we find that, with 10,000 observations and one endogenous regressor, IQR is nine times slower than the contraction algorithm and 23 times slower than the root-finding algorithm. With two endogenous regressors, the computational advantages of our procedures are even more pronounced. IQR based on a grid search over 100$\times$100 points is 149 times slower than the contraction algorithm and 84 times slower than the nested root-finding algorithm.
Figure (ref) displays the estimates of $\theta_D(\tau)$ for $\tau \in \{0.15,0.20,\dots,0.85\}$. We can see that all estimation algorithms yield very similar results. We also note that the contraction algorithm converges for all quantile levels considered.
Figures (ref) depicts pointwise 95% confidence intervals for the proposed estimators obtained using the empirical bootstrap described in Section (ref) with 500 replications. We can see that the resulting confidence intervals are very similar for both algorithms and do not include zero at all quantile levels considered.
In this section, we assess the practical performance of our estimation algorithms in an empirical Monte Carlo study based on the application in Section (ref).
We consider DGPs which are based on the empirical application of Section (ref).\footnote{The construction of our DGPs is inspired by the application-based DGPs in KaplanSun2017.} We focus on a simplified setting with only two exogenous covariates: income and age. The covariates are drawn from their joint empirical distribution. The instrument $Z_i$ is generated as $\text{Bernoulli}\left(\bar{Z} \right)$, where $\bar{Z}$ is the mean of the instrument in the empirical application. We then generate the endogenous variable as $D_i=Z_i\cdot 1\left\{ 0.6\cdot V_i < U_i\right\},$ where $U_i\sim \text{Uniform}(0,1)$ and $V_i\sim \text{Uniform}(0,1)$ are independent disturbances. The DGP for $D_i$ is chosen to roughly match the joint empirical distribution of $(D_i,Z_i)$. The outcome variable $Y_i$ is generated as
The coefficient $\theta_X(\cdot)$ is constant and equal to the IQR median estimate in the empirical application. $\theta_D(U_i)=5000+U_i\cdot 10000$ is chosen to match the increasing shape of the estimated conditional quantile treatment effects in Figure (ref). $G^{-1}(\cdot)$ is the quantile function of a re-centered Gamma distribution, estimated to match the distribution of the IQR residuals at the median. To investigate the performance of our procedure with more than one endogenous variable, we add a second endogenous regressor:
where we set $\theta_{D,2}(U_i)=10000$. The second endogenous variable is generated as \[ D_{2,i}=0.8 \cdot Z_{2,i}+0.2 \cdot \Phi^{-1}(U_i) \] and the second instrument is generated as $Z_{2,i}\sim N(0,1)$. We set $N=9913$ as in the empirical application.
We assess and compare several different algorithms, all of which are based on the dynamical system $\hat M$.\footnote{We do not explore algorithms based on $\hat K$ because, as discussed in Section (ref), the algorithms based on $\hat M$ have advantages in terms of the choice of starting values (for the contraction algorithm) and dimension reduction (for the root-finding algorithm).} For models with one endogenous variable, we consider a contraction algorithm and a root-finding algorithm based on Brent's method. For models with two endogenous variables, we consider a contraction algorithm, a nested root-finding algorithm based on Brent's method, and a root-finding algorithm implemented as a minimization problem based on simulated annealing (SA).\footnote{We have also explored algorithms based on Newton-Raphson-type root-finders. These algorithms are, in theory, up to an order of magnitude faster than the contraction algorithm and the nested algorithms, but, unlike the other algorithms considered here, require an approximation to the Jacobian and are not very robust to the choice of starting values. We therefore do not report the results here.} For all estimators, we use 2SLS estimates as starting values. We compare the results of our algorithms to (nested) profiling estimators based on Brent's method, and to IQR with a grid search over 500 points (one endogenous regressor) and 100$\times$100 points (two endogenous regressors), which serves as a slow but very robust benchmark. Table (ref) presents an overview of the algorithms.
Here we describe the computational performance of the different procedures and the finite sample performance of our bootstrap inference procedure. In Appendix (ref), we further investigate the finite sample bias and root mean squared error (RMSE) of the different methods. We find that the proposed estimation algorithms perform well and exhibit a similar bias and RMSE. Moreover, the finite sample properties of the estimators based on the contraction and root-finding algorithms are comparable to the profiling estimators and IQR. This shows that the computational advantages of our algorithms do not come at a cost in terms of the finite sample performance. Appendix (ref) presents additional simulation evidence, demonstrating that our algorithms perform well and remain computationally tractable with more than two endogenous regressors.\footnote{Specifically, we present simulation results based on the DGP used in this section, augmented with an additional endogenous regressor, generated as $D_{3,i}=0.8 \cdot Z_{3,i}+0.2 \cdot \Phi^{-1}(U_i)$, where $Z_{3,i}\sim N(0,1)$.}
Tables (ref)--(ref) show the average computation time (in seconds) for estimating the model with one and two endogenous variables for different sample sizes. All computations were carried out on a standard desktop computer with a 3.2 GHz Intel Core i5 processor and 8GB RAM.
With one endogenous regressor, the contraction algorithm and the root-finding algorithm based on Brent's method are computationally more efficient than IQR. Specifically, the root-finding algorithm based on Brent's method is 11 to 23 times faster than IQR and the contraction algorithm is 4 to 9 times faster. The root-finding algorithm is as fast as the profiling estimator and about twice as fast as the contraction algorithm.\footnote{Note that the computational speed of the contraction algorithm depends on $|\hat J_M|$ and, thus, will differ across applications.}
The computational advantages of our algorithms become more pronounced with two endogenous variables. Table (ref) shows that IQR's average computation times are around two orders of magnitude slower than those of our preferred procedures. Specifically, the nested root-finding algorithm is 84 to 134 times faster than IQR, while the contraction algorithm is 149 to 308 times faster. This is as expected since, due to the use of grids, IQR's becomes computationally impractical in our setting whenever the number of endogenous variables exceeds two or three.\footnote{Our implementation of IQR with two endogenous variables is inherently slower than the implementation with one endogenous variable, even when the number of grid points is the same. First, there is an additional covariate in the underlying QRs (the second instrument). Second, with one endogenous variable, we choose the grid value that minimizes the absolute value of the coefficient on the instrument. By contrast, with two endogenous regressors, we choose the grid point which minimizes a quadratic form based on the inverse of the estimated QR variance covariance matrix as suggested in Chernozhukov+17handbook, which requires an additional computational step.} The contraction algorithm is almost twice as fast as the nested algorithm, which, in turn, is almost twice as fast as the profiling estimator. Finally, the minimization algorithm based on SA is about one order of magnitude slower than the contraction algorithm, while still being an order of magnitude faster than IQR.
Finally, we analyze the finite sample properties of our bootstrap inference procedure for making inference on $\theta_D$ in the model with a single endogenous variable. Table (ref) shows the empirical coverage probabilities of bootstrap confidence intervals based on the contraction algorithm and the root-finding algorithm based on Brent's method. Both methods demonstrate an excellent performance and exhibit coverage rates that are very close to the respective nominal levels.
The main contribution of this paper is to develop computationally convenient and easy-to-implement estimation algorithms for IVQR models. Our key insight is that the non-smooth and non-convex IVQR estimation problem can be decomposed into a sequence of much more tractable convex QR problems, which can be solved very quickly using well-established methods. The proposed algorithms are particularly well-suited if the number of exogenous variables is large and the number of endogenous variables is moderate as in many empirical applications.
An interesting avenue for further research is to investigate weak identification robust inference within the decentralized model. One may, for example, write the (re-scaled) sample fixed point restriction as $\sqrt N(I_{d_X+d_D}-\hat K)(\theta)=s_N(\theta)+\mathbb W(\theta)+r_N(\theta)$, where $s_N(\theta)=\sqrt N(I_{d_X+d_D}-K)(\theta),$ $\mathbb W$ is a Gaussian process, and $r_N$ is an error that tends to zero uniformly. This paper assumes that $s_N(\theta^*)=0$ uniquely, and outside $N^{-1/2}$-neighborhoods of $\theta^*$, $s_N(\theta)$ diverges and dominates $\mathbb W.$ For a one-dimensional fixed point problem, this requires the BR map to be bounded away from the 45-degree line outside any $N^{-1/2}$-neighborhood of the fixed point. However if $s_N$ fails to dominate $\mathbb W$ over a substantial part of the parameter space, one would end up with weak identification.\footnote{AndrewsMikusheva2016 study weak identification robust inference methods in models characterized by moment restrictions.} How to conduct robust inference in such settings is an interesting question, which we leave for future research.
Finally, we note that, while we study the performance of the proposed algorithms separately, our reformulation and the resulting algorithms are potentially very useful when combined with other existing procedures. For instance, one could choose starting values using an initial grid search over a coarse grid and then apply a fast contraction algorithm.