Decentralization Estimators for Instrumental Variable Quantile Regression Models

TL;DR

This paper introduces a decomposition approach for IVQR models that simplifies estimation by breaking it into convex quantile regressions, enabling faster and more accessible causal effect estimation.

Contribution

It presents a novel decomposition method that transforms the complex IVQR estimation into convex sub-problems, improving efficiency and implementation.

Findings

Decomposition into convex quantile regressions simplifies IVQR estimation.

New identification results for IVQR models.

Proposed estimators are fast, easy to implement, and tuning-free.

Abstract

The instrumental variable quantile regression (IVQR) model (Chernozhukov and Hansen, 2005) is a popular tool for estimating causal quantile effects with endogenous covariates. However, estimation is complicated by the non-smoothness and non-convexity of the IVQR GMM objective function. This paper shows that the IVQR estimation problem can be decomposed into a set of conventional quantile regression sub-problems which are convex and can be solved efficiently. This reformulation leads to new identification results and to fast, easy to implement, and tuning-free estimators that do not require the availability of high-level "black box" optimization routines.

Figures

Tables

Table 1. Table 1. Algorithms

One endogenous variable
Algorithm	R-Package
Contraction algorithm
Root-finding algorithm	uniroot (R Core Team, 2019)
Profiling algorithm	uniroot (R Core Team, 2019)
IQR
Two endogenous variables
Algorithm	R-Package
Contraction algorithm
Nested root-finding algorithm	uniroot (R Core Team, 2019)
Root-finding algorithm (implemented as optimizer)	optim_sa (Husmann, Lange, and Spiegel, 2017)
Nested profiling algorithm	uniroot (R Core Team, 2019)
IQR (implementation: p.132 in
Chernozhukov, Hansen, and Wüthrich (2017))

Table 2. Table 2. Computation time, 401(k) DGP with one endogenous regressor

$N$	Contr	Brent	Profil	InvQR
1000	0.10	0.03	0.03	0.34
5000	0.44	0.13	0.12	2.44
10000	0.75	0.28	0.30	6.48
Notes: The table reports average computation time in seconds at $\tau =0.5$ over 100 simulation repetitions based on the DGP described in the main text. Contr: contraction algorithm; Brent: root-finding algorithm based on Brent’s method; Profil: profiling estimator based on Brent’s method; InvQR: inverse QR with grid search over 500 grid points. We use 2SLS estimates as starting values.

Table 3. Table 3. Computation time, 401(k) DGP with two endogenous regressor

$N$	Contr	NestBr	SimAnn	Profil	InvQR
1000	0.18	0.41	3.13	0.77	55.40
5000	1.71	2.67	19.05	4.39	282.64
10000	4.50	7.99	54.95	12.50	672.00
Notes: The table reports average computation time in seconds at $\tau =0.5$ over 100 simulation repetitions based on the DGP described in the main text. Contr: contraction algorithm; NestBr: nested algorithm based on Brent’s method; SimAnn: simulated annealing based optimization algorithm; Profil: nested profiling estimator based on Brent’s method; InvQR: inverse QR with grid search over 100$\times 100$ grid points. We use 2SLS estimates as starting values.

Table 4. Table 4. Coverage, 401(k) DGP with one endogenous regressor

	$1-\alpha =0.95$		$1-\alpha =0.9$
$\tau$	Contr	Brent	Contr	Brent
0.15	0.96	0.96	0.93	0.91
0.25	0.95	0.95	0.91	0.92
0.50	0.95	0.95	0.90	0.91
0.75	0.95	0.95	0.90	0.89
0.85	0.94	0.94	0.89	0.89
Notes: Monte Carlo simulation with 1000 repetitions as described in the main text. Contr: contraction algorithm; Brent: root-finding algorithm based on Brent’s method. We use 2SLS estimates as starting values.

Table 5. Table 5. Bias and RMSE, 401(k) DGP with one endogenous regressor

	Bias/${10}^2$				RMSE/${10}^3$
$\tau$	Contr	Brent	Profil	InvQR	Contr	Brent	Profil	InvQR
0.15	-4.43	-6.84	-6.71	-7.26	7.58	7.71	7.64	7.87
0.25	-0.28	-1.86	-1.90	-1.67	4.04	4.10	4.11	4.10
0.50	-1.04	-1.46	-1.49	-1.23	2.06	2.07	2.07	2.08
0.75	-1.60	-1.31	-1.48	-1.14	1.85	1.85	1.86	1.85
0.85	0.61	1.17	0.85	1.24	2.06	2.07	2.07	2.08
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; Brent: root-finding algorithm based on Brent’s method; Profil: profiling estimator based on Brent’s method; InvQR: inverse quantile regression. We use 2SLS estimates as starting values.

Table 6. Table 6. Bias and RMSE, 401(k) DGP with two endogenous regressors

	Bias/${10}^2$					RMSE/${10}^3$
$\tau$	Contr	NestBr	SimAnn	Profil	InvQR	Contr	NestBr	SimAnn	Profil	InvQR
Coefficient on binary endogenous variable
0.15	-4.81	-2.73	3.77	-3.65	-7.77	8.29	7.84	6.74	7.95	8.60
0.25	-3.47	-3.75	-3.17	-3.83	-3.42	4.32	4.31	4.25	4.31	4.32
0.50	0.84	0.56	0.68	0.71	0.74	1.93	1.95	1.95	1.95	1.98
0.75	-0.56	-0.33	-0.37	-0.57	-0.25	1.75	1.74	1.74	1.74	1.78
0.85	-1.03	-0.72	-0.74	-1.28	-0.61	2.18	2.19	2.20	2.19	2.22
Coefficient on continuous endogenous variable
0.15	2.00	0.48	4.72	-0.03	0.23	1.07	1.07	2.20	1.07	1.19
0.25	1.72	-0.09	0.25	-0.48	-0.10	1.00	1.02	1.13	1.03	1.13
0.50	0.75	-0.47	-0.45	-0.54	-0.68	0.89	0.97	0.97	0.97	1.11
0.75	-1.43	-0.49	-0.44	-1.57	-0.15	0.99	1.10	1.11	1.12	1.24
0.85	-2.66	-0.89	-0.91	-2.60	-0.40	1.12	1.27	1.27	1.28	1.32
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; NestBr: nested algorithm based Brent’s method; SimAnn: simulated annealing based optimization algorithm; Profil: nested profiling estimator based on Brent’s method; IQR: inverse quantile regression. We use 2SLS estimates as starting values.

Table 7. Table 7. Bias and RMSE, 401(k) DGP with three endogenous regressors

	Bias/${10}^2$		RMSE/${10}^3$
$\tau$	Contr	Nested	Contr	Nested
Coefficient on $D$
0.15	-3.40	-5.02	7.52	7.62
0.25	-0.98	-1.52	4.11	4.12
0.50	-1.04	-1.42	2.03	2.06
0.75	-1.67	-1.50	1.86	1.86
0.85	0.86	1.03	2.05	2.05
Coefficient on $D_2$
0.15	1.12	-0.13	1.01	1.03
0.25	2.11	0.74	1.01	1.00
0.50	0.80	-0.28	0.94	0.99
0.75	-0.39	0.48	1.00	1.09
0.85	-2.64	-0.83	1.13	1.25
Coefficient on $D_3$
0.15	1.70	-0.25	1.08	1.13
0.25	1.57	-0.21	0.99	1.01
0.50	0.95	-0.06	0.92	0.96
0.75	-1.15	-0.33	1.02	1.11
0.85	-1.01	0.23	1.22	1.37
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; Nested: nested algorithm based on Brent’s method. We use 2SLS estimates as starting values.

Table 8. Table 8. Computation time, 401(k) DGP with three endogenous regressor

$N$	Contr	Nested
1000	0.36	6.29
5000	4.47	42.93
10000	10.11	145.58
Notes: The table reports average computation time in seconds at $\tau =0.5$ over 20 simulation repetitions based on the DGP described in the main text. Contr: contraction algorithm; Nested: nested algorithm based on Brent’s method. We use 2SLS estimates as starting values.

Table 9. Table 9. Bias and RMSE, symmetric design with one endogenous regressor

	$N=500$
	Bias				RMSE
$\tau$	Contr	Brent	Profil	InvQR	Contr	Brent	Profil	InvQR
0.15	0.03	-0.00	-0.02	-0.00	0.10	0.10	0.11	0.10
0.25	0.03	0.00	-0.01	0.00	0.12	0.12	0.12	0.12
0.50	-0.00	-0.00	-0.02	-0.00	0.12	0.14	0.14	0.14
0.75	-0.04	-0.01	-0.03	-0.01	0.13	0.12	0.12	0.12
0.85	-0.04	-0.00	-0.02	-0.00	0.11	0.11	0.11	0.11
	$N=1000$
	Bias				RMSE
$\tau$	Contr	Brent	Profil	InvQR	Contr	Brent	Profil	InvQR
0.15	0.02	0.00	-0.00	0.00	0.07	0.07	0.07	0.07
0.25	0.01	-0.00	-0.01	-0.00	0.08	0.08	0.08	0.08
0.50	-0.01	-0.01	-0.01	-0.01	0.09	0.10	0.10	0.10
0.75	-0.02	-0.00	-0.01	-0.00	0.09	0.08	0.08	0.08
0.85	-0.02	-0.00	-0.01	-0.00	0.08	0.08	0.08	0.08
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; Brent: root-finding algorithm based on Brent’s method; Profil: profiling estimator based on Brent’s method; InvQR: inverse quantile regression. We use 2SLS estimates as starting values.

Table 10. Table 10. Bias and RMSE, asymmetric design with one endogenous regressor

	$N=500$
	Bias				RMSE
$\tau$	Contr	Brent	Profil	InvQR	Contr	Brent	Profil	InvQR
0.15	0.11	0.01	-0.04	-0.00	0.22	0.20	0.21	0.20
0.25	0.07	0.00	-0.02	-0.00	0.17	0.16	0.17	0.16
0.50	0.04	-0.00	-0.02	-0.00	0.13	0.12	0.12	0.12
0.75	0.03	0.00	-0.01	0.00	0.11	0.11	0.11	0.11
0.85	-0.03	-0.01	-0.03	-0.00	0.12	0.11	0.12	0.11
	$N=1000$
	Bias				RMSE
$\tau$	Contr	Brent	Profil	InvQR	Contr	Brent	Profil	InvQR
0.15	0.05	-0.01	-0.03	-0.01	0.16	0.15	0.15	0.15
0.25	0.04	0.00	-0.01	0.00	0.11	0.11	0.11	0.11
0.50	0.03	0.00	-0.01	0.00	0.08	0.08	0.08	0.08
0.75	0.01	-0.01	-0.02	-0.01	0.08	0.08	0.08	0.08
0.85	-0.03	-0.01	-0.02	-0.01	0.09	0.09	0.09	0.09
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; Brent: root-finding algorithm based on Brent’s method; Profil: profiling estimator based on Brent’s method; InvQR: inverse quantile regression. We use 2SLS estimates as starting values.

Table 11. Table 11. Bias and RMSE, symmetric design with two endogenous regressors

	$N=500$
	Bias					RMSE
$\tau$	Contr	NestBr	SimAnn	Profil	InvQR	Contr	NestBr	SimAnn	Profil	InvQR
Coefficient on $D_1$
0.15	0.00	-0.00	0.00	-0.02	-0.01	0.11	0.12	0.14	0.13	0.13
0.25	0.01	-0.00	-0.01	-0.02	-0.01	0.15	0.16	0.17	0.16	0.16
0.50	-0.02	-0.02	-0.02	-0.04	-0.02	0.17	0.19	0.19	0.19	0.20
0.75	-0.04	-0.03	-0.03	-0.05	-0.03	0.21	0.20	0.21	0.20	0.20
0.85	-0.05	-0.03	-0.03	-0.05	-0.03	0.18	0.17	0.18	0.18	0.17
Coefficient on $D_2$
0.15	0.10	-0.01	-0.01	-0.05	-0.02	0.27	0.27	0.29	0.28	0.31
0.25	0.10	-0.00	-0.02	-0.04	-0.02	0.29	0.29	0.30	0.30	0.30
0.50	-0.01	-0.02	-0.02	-0.06	-0.02	0.33	0.38	0.39	0.40	0.39
0.75	-0.15	-0.04	-0.06	-0.08	-0.05	0.40	0.40	0.41	0.41	0.41
0.85	-0.19	-0.05	-0.06	-0.10	-0.07	0.39	0.36	0.40	0.38	0.43
	$N=1000$
	Bias					RMSE
$\tau$	Contr	NestBr	SimAnn	Profil	InvQR	Contr	NestBr	SimAnn	Profil	InvQR
Coefficient on $D_1$
0.15	-0.00	-0.00	-0.01	-0.01	-0.00	0.08	0.09	0.10	0.09	0.10
0.25	-0.00	-0.00	-0.01	-0.01	-0.01	0.10	0.11	0.12	0.11	0.13
0.50	-0.01	-0.01	-0.01	-0.02	-0.01	0.12	0.13	0.13	0.13	0.16
0.75	-0.01	-0.01	-0.01	-0.01	-0.00	0.13	0.13	0.14	0.13	0.14
0.85	-0.02	-0.01	-0.02	-0.02	-0.02	0.12	0.12	0.13	0.12	0.13
Coefficient on $D_2$
0.15	0.05	-0.01	-0.01	-0.02	-0.02	0.19	0.19	0.21	0.19	0.20
0.25	0.05	-0.00	-0.01	-0.02	-0.01	0.22	0.21	0.23	0.22	0.23
0.50	-0.02	-0.02	-0.02	-0.04	-0.03	0.25	0.27	0.27	0.28	0.29
0.75	-0.09	-0.02	-0.02	-0.04	-0.03	0.27	0.25	0.28	0.25	0.26
0.85	-0.09	-0.01	-0.03	-0.04	-0.03	0.26	0.23	0.25	0.24	0.24
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; NestBr: nested algorithm based on Brent’s method; SimAnn: simulated annealing based optimization algorithm; Profil: nested profiling estimator based on Brent’s method; InvQR: inverse quantile regression. We use 2SLS estimates as starting values.

Table 12. Table 12. Bias and RMSE, asymmetric design with two endogenous regressors

	$N=500$
	Bias					RMSE
$\tau$	Contr	NestBr	SimAnn	Profil	InvQR	Contr	NestBr	SimAnn	Profil	InvQR
Coefficient on $D_1$
0.15	-0.02	0.02	0.00	-0.03	0.01	0.25	0.26	0.28	0.27	0.26
0.25	-0.05	0.01	0.00	-0.01	-0.00	0.20	0.20	0.21	0.20	0.21
0.50	-0.04	-0.00	0.00	-0.01	0.00	0.16	0.17	0.21	0.17	0.19
0.75	-0.02	-0.02	-0.01	-0.02	-0.02	0.17	0.17	0.18	0.18	0.19
0.85	-0.01	-0.01	-0.02	-0.02	-0.02	0.20	0.19	0.19	0.20	0.19
Coefficient on $D_2$
0.15	0.26	-0.06	-0.11	-0.16	-0.13	0.57	0.52	0.58	0.53	0.59
0.25	0.23	-0.01	-0.02	-0.06	-0.01	0.45	0.41	0.43	0.40	0.44
0.50	0.12	-0.03	-0.04	-0.07	-0.07	0.34	0.32	0.48	0.33	0.73
0.75	0.04	-0.05	-0.06	-0.11	-0.05	0.32	0.31	0.34	0.33	0.34
0.85	-0.13	-0.03	-0.01	-0.08	0.01	0.40	0.34	0.38	0.34	0.36
	$N=1000$
	Bias					RMSE
$\tau$	Contr	NestBr	SimAnn	Profil	InvQR	Contr	NestBr	SimAnn	Profil	InvQR
Coefficient on $D_1$
0.15	-0.03	0.01	-0.00	-0.02	-0.01	0.18	0.19	0.19	0.18	0.19
0.25	-0.04	-0.00	-0.01	-0.02	-0.01	0.15	0.15	0.16	0.15	0.16
0.50	-0.03	-0.01	-0.01	-0.01	-0.01	0.13	0.13	0.14	0.13	0.14
0.75	-0.03	-0.01	-0.01	-0.02	-0.01	0.12	0.12	0.13	0.13	0.14
0.85	0.01	0.00	0.00	-0.01	-0.00	0.14	0.13	0.15	0.14	0.15
Coefficient on $D_2$
0.15	0.15	-0.03	-0.03	-0.07	-0.04	0.37	0.37	0.38	0.37	0.39
0.25	0.10	-0.01	-0.01	-0.05	-0.02	0.28	0.28	0.30	0.28	0.28
0.50	0.05	-0.02	-0.03	-0.04	-0.03	0.22	0.22	0.23	0.22	0.24
0.75	0.06	-0.01	-0.02	-0.03	-0.01	0.24	0.22	0.24	0.23	0.24
0.85	-0.08	-0.03	-0.04	-0.07	-0.03	0.27	0.24	0.26	0.25	0.24
Notes: Monte Carlo simulation with 500 repetitions as described in the main text. Contr: contraction algorithm; NestBr: nested algorithm based on Brent’s method; SimAnn: simulated annealing based optimization algorithm; Profil: nested profiling estimator based on Brent’s method; InvQR: inverse quantile regression. We use 2SLS estimates as starting values.

Table 13. Table 13. Coverage, location-scale DGP with one endogenous regressor

$N=500$
	Symmetric Design				Asymmetric Design
	$1-\alpha =0.95$		$1-\alpha =0.9$		$1-\alpha =0.95$		$1-\alpha =0.9$
$\tau$	Contr	Brent	Contr	Brent	Contr	Brent	Contr	Brent
0.15	0.93	0.98	0.88	0.94	0.89	0.97	0.85	0.95
0.25	0.94	0.96	0.89	0.93	0.93	0.97	0.88	0.95
0.50	0.96	0.96	0.91	0.91	0.95	0.97	0.91	0.93
0.75	0.94	0.97	0.90	0.94	0.95	0.97	0.91	0.93
0.85	0.95	0.98	0.90	0.96	0.96	0.98	0.92	0.96
$N=1000$
	Symmetric Design				Asymmetric Design
	$1-\alpha =0.95$		$1-\alpha =0.9$		$1-\alpha =0.95$		$1-\alpha =0.9$
$\tau$	Contr	Brent	Contr	Brent	Contr	Brent	Contr	Brent
0.15	0.95	0.96	0.90	0.92	0.92	0.97	0.85	0.93
0.25	0.95	0.96	0.91	0.91	0.92	0.95	0.87	0.91
0.50	0.96	0.96	0.90	0.90	0.94	0.96	0.90	0.93
0.75	0.95	0.96	0.90	0.91	0.95	0.95	0.91	0.91
0.85	0.96	0.97	0.93	0.94	0.96	0.95	0.92	0.91
Notes: Monte Carlo simulation with 1000 repetitions as described in the main text. Contr: contraction algorithm; Brent: root-finding algorithm based on Brent’s method. We use 2SLS estimates as starting values.

Equations

$$Y_d=q(d,x,U_d)\text{with}{U}_d\sim U(0,1).$$

$$P\left(Y\le q(D,X,\tau )\mid X,Z\right)=\tau ,\tau \in (0,1).$$

$$q(d,x,\tau )=x^{\prime }{\theta }_X(\tau )+d^{\prime }{\theta }_D(\tau ),$$

$${\Psi }_P\left(\theta (\tau )\right):=E_P\left[\left(1\left\{Y\le X^{\prime }{\theta }_X(\tau )+D^{\prime }{\theta }_D(\tau )\right\}-\tau \right)\left(\begin{array}{c}X\\ Z\end{array}\right)\right].$$

$${\Psi }_P\left({\theta }^{\ast }\right)=0.$$

$${\mathcal{Q}}_N^{GMM}\left(\theta \right)=-\frac{1}{2}{\left(\frac{1}{\sqrt{N}}\sum _{i=1}^N{m}_i\left(\theta \right)\right)}^{\prime }{W}_N\left(\theta \right)\left(\frac{1}{\sqrt{N}}\sum _{i=1}^N{m}_i\left(\theta \right)\right),$$

$$Q_{P,1}\left(\theta \right)$$

$$Q_{P,j}\left(\theta \right)$$

$$\underset{{\tilde{\theta }}_1\in {\mathbb{R}}^{d_X}}{\min }{Q}_{P,1}\left({\tilde{\theta }}_1,{\theta }_{-1}\right),$$

$$\underset{{\tilde{\theta }}_j\in \mathbb{R}}{\min }{Q}_{P,j}\left({\tilde{\theta }}_j,{\theta }_{-j}\right),j=2,\dots ,J.$$

$${\Psi }_{P,1}({\tilde{\theta }}_1^{\ast },{\theta }_{-1})$$

$${\Psi }_{P,j}({\tilde{\theta }}_j^{\ast },{\theta }_{-j})$$

$${\theta }_j^{\ast }\in {\tilde{L}}_j({\theta }_{-j}^{\ast }),j=1,\dots ,J,$$

$$E_P[f_{Y|D,X,Z}\left(D^{\prime }{\theta }_{-1}+X^{\prime }{\theta }_1\right)X{X}^{\prime }]$$

$$E_P[f_{Y|D,X,Z}\left(D^{\prime }{\theta }_{-1}+X^{\prime }{\theta }_1\right)D_{\ell }{Z}_{\ell }],\ell =1,\dots ,d_D,$$

$$R_{-j}:=\{{\theta }_{-j}\in {\Theta }_{-j}:{\Psi }_{P,j}(\theta )=0,\text{ for some }\theta =({\theta }_j,{\theta }_{-j})\in \Theta \}.$$

$${\Psi }_{P,j}\left(L_j({\theta }_{-j}),{\theta }_{-j}\right)=0,\text{for all}{\theta }_{-j}\in R_{-j}.$$

$$K(\theta )=\left(\begin{array}{c}{K}_1(\theta )\\ ⋮\\ K_J(\theta )\end{array}\right)=\left(\begin{array}{c}{L}_1({\theta }_{-1})\\ ⋮\\ L_J({\theta }_{-J})\end{array}\right).$$

$$K(\theta )=\left(\begin{array}{c}{L}_1({\theta }_2)\\ L_2({\theta }_1)\end{array}\right).$$

$$M\left({\theta }_{-1}\right)=\left(\begin{array}{c}{M}_1({\theta }_{-1})\\ M_2({\theta }_{-1})\\ ⋮\\ M_{d_D}({\theta }_{-1})\end{array}\right)=\left(\begin{array}{c}{L}_2\left(L_1({\theta }_{-1}),{\theta }_{-\{1,2\}}\right)\\ L_3\left(L_1({\theta }_{-1}),L_2(L_1({\theta }_{-1}),{\theta }_{-\{1,2\}}),{\theta }_{-\{1,2,3\}}\right)\\ ⋮\\ L_J\left(L_1({\theta }_{-1}),L_2(L_1({\theta }_{-1}),{\theta }_{-\{1,2\}}),\cdots \right)\end{array}\right),$$

$$M({\theta }_2)=L_2\left(L_1\left({\theta }_2\right)\right),$$

$$\displaystyle\tilde{R}_{1}:=\big{\{}\theta_{-1}\in\Theta_{-1}:~{}$$

$$\displaystyle\Psi_{P,2}(\theta_{1},\theta_{2},\pi_{-\{1,2\}}\theta_{-1})=0,~{}\exists(\theta_{1},\theta_{2})\in\Theta_{1}\times\Theta_{2}\big{\}}.$$

$$D_M:=\bigcap _{j=1}^{d_D}{\tilde{R}}_j={\tilde{R}}_{d_D},$$

$${\theta }^{(s+1)}=K\left({\theta }^{(s)}\right),s=0,1,2,\dots ,{\theta }^{(0)}\text{given}.$$

$${\theta }_{-1}^{(s+1)}$$

$$J_K(\theta )=\left(\begin{array}{cc}0& J_{L_1}({\theta }_2)\\ J_{L_2}({\theta }_1)& 0\end{array}\right)\text{and}{J}_M({\theta }_2)=J_{L_2}\left(L_1({\theta }_2)\right)J_{L_1}\left({\theta }_2\right),$$

$$\displaystyle J_{L_{-j}}\left(\theta_{j}\right)=-\left(\frac{\partial\Psi_{P,-j}(\theta_{j},\theta_{-j})}{\partial\theta_{-j}^{\prime}}\bigg{|}_{\theta=(\theta_{j},L_{-j}(\theta_{j}))}\right)^{-1}\frac{\partial\Psi_{P,-j}(\theta_{j},\theta_{-j})}{\partial\theta_{j}^{\prime}}\bigg{|}_{\theta=(\theta_{j},L_{-j}(\theta_{j}))},~{}j=1,2.$$

$$\det \left(J_{{\Psi }_P}(\theta )\right)$$

$$0\ne \det \left(\begin{array}{cc}{I}_{d_1}& -J_{L_1}({\theta }_2^{\ast })\\ -J_{L_2}({\theta }_1^{\ast })& I_{d_2}\end{array}\right)=\det (I_{d_X+d_D}-J_K({\theta }^{\ast })).$$

$$(I_{d_X+d_D}-K)({\theta }^{\ast })=0,$$

$$(I_{d_D}-M)({\theta }_{-1}^{\ast })=0.$$

Full text

**Keywords: instrumental variables, quantile regression, contraction mapping, fixed point estimator, bootstrap.

JEL Codes: C21, C26

**

1. Introduction

Quantile regression (QR), introduced by Koenker and Bassett (1978), 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 Chernozhukov and Hansen (2005, 2006) 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 Chernozhukov and Hansen (2006). 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 Chernozhukov and Hansen (2006). 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 Chernozhukov and Hansen (2004). Based on this application, we perform extensive Monte Carlo simulations, which demonstrate that our estimation and inference procedures have excellent finite sample properties.

1.1. Literature

We contribute to the literature on estimation and inference based on linear IVQR models. Chernozhukov and Hong (2003) have proposed a quasi-Bayesian approach which can accommodate multiple endogenous variables but, as noted by Chernozhukov and Hansen (2013), requires careful tuning in applications. Chernozhukov and Hansen (2006) 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. Chernozhukov and Hansen (2008) and Jun (2008) have studied weak instrument robust inference procedures based on the inversion of Anderson-Rubin-type tests. Chernozhukov, Hansen, and Jansson (2009) have proposed a finite sample inference approach. Andrews and Mikusheva (2016) have developed a general conditional inference approach and derived sufficient conditions for the IVQR model. Kaplan and Sun (2017) and de Castro, Galvao, Kaplan, and Liu (2018) have suggested to use smoothed estimating equations to overcome the non-smoothness of the IVQR estimation problem, although the non-convexity remains. More recently, Chen and Lee (2018) 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, Zhu (2018) 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 Chen and Lee (2018). This procedure typically requires an early termination of the algorithm to ensure computational tractability which is akin to a tuning parameter choice. In addition, Zhu (2018) has proposed a $k$-step approach that allows for estimating models with multiple endogenous regressors based on large datasets, but requires estimating the gradient. Pouliot (2018) 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, Machado and Santos Silva (2018) 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 Chernozhukov, Imbens, and Newey (2007), Horowitz and Lee (2007), Chen and Pouzo (2009), Chen and Pouzo (2012), Gagliardini and Scaillet (2012), and Wüthrich (2019). Chernozhukov and Hansen (2013) and Chernozhukov, Hansen, and Wüthrich (2017) have provided surveys on the IVQR model including references to empirical applications.

Abadie, Angrist, and Imbens (2002) 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 Imbens and Angrist (1994). Their approach has been extended and further developed by Frandsen, Frölich, and Melly (2012), Frölich and Melly (2013), and Belloni, Chernozhukov, Fernandez-Val, and Hansen (2017) among others. We refer to Melly and Wüthrich (2017) for a recent review of this approach and to Wüthrich (2020) for a comparison between this approach and the IVQR model. Identification and estimation in nonseparable models with continuous endogenous regressors have been studied by Chesher (2003), Ma and Koenker (2006), Lee (2007), Jun (2009), Imbens and Newey (2009), D’Haultfoeuille and Février (2015), and Torgovitsky (2015) 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 (Weisberg and Welsh, 1994), general maximum likelihood problems (Smyth, 1996), linear models with high-dimensional fixed effects (e.g., Guimaraes and Portugal, 2010, and the references therein), sample selection models (Marra and Radice, 2013), peer effects models (Arcidiacono, Foster, Goodpaster, and Kinsler, 2012), interactive fixed effects models (e.g., Chen, Fernandez-Val, and Weidner, 2014; Moon and Weidner, 2015), and random coefficient logit demand models (Lee and Seo, 2015). 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.

1.2. Organization of the Paper

The remainder of the paper is structured as follows. Section 2 introduces the setup and the IVQR model. Section 3 shows that the IVQR estimation problem can be decentralized into a series of (weighted) conventional QR problems. In Section 4, we introduce population algorithms based on the contraction mapping theorem and root-finders. Section 5 discusses the corresponding sample algorithms. In Section 6, we establish the asymptotic normality of our estimators and prove the validity of the bootstrap. Section 7 presents an empirical application. In Section 8, we provide simulation evidence on the computational performance and the finite sample properties of our methods. Section 9 concludes. All proofs as well as some additional theoretical and simulation results are collected in the appendix.

2. Setup and Model

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 (e.g., Rubin, 1974). 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

[TABLE]

This representation lies at the heart of the IVQR model. With this notation at hand, we state the main conditions of the IVQR model (Chernozhukov and Hansen, 2005, Assumptions A1-A5).

Assumption 1**.**

Given a common probability space $(\Omega,F,P)$, the following conditions hold jointly with probability one:

1. (1)

Potential outcomes: Conditional on $X=x$, for each $d$, $Y_{d}=q(d,x,U_{d})$, where $q(d,x,\tau)$ is strictly increasing in $\tau$ and $U_{d}\sim U(0,1)$. 2. (2)

Independence: Conditional on $X=x$, $\{U_{d}\}$ are independent of $Z$. 3. (3)

Selection: $D:=\delta(Z,X,V)$ for some unknown function $\delta(\cdot)$ and random vector $V$. 4. (4)

Rank invariance or Rank similarity: Conditional on $X=x$, $Z=z$,

• (a)

$\{U_{d}\}$* are equal to each other; or, more generally,*

• (b)

$\{U_{d}\}$* are identically distributed, conditional on $V$.* 5. (5)

Observed variables: Observed variables consist of $Y:=q(D,X,U_{D})$, $D$, $X$, and $Z$.

We briefly discuss the most important aspects of Assumption 1 and refer the interested reader to Chernozhukov and Hansen (2005, 2006, 2013) for more comprehensive treatments. Assumption 1.1 states the Skorohod representation of $Y_{d}$ and requires strict monotonicity of the potential outcome quantile function, which rules out discrete outcomes. Assumption 1.2 imposes independence between the potential outcomes and the instrument. Assumption 1.3 defines a general selection mechanism. The key restriction of the IVQR model is Assumption 1.4. 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 1.5 summarizes the observables.

Remark 2.1**.**

Assumption 1 does not impose any restrictions on how the instrument $Z$ affects the endogenous variable $D$. As a consequence, Assumption 1 alone does not guarantee point identification of the structural quantile function. In Sections 3 and 4 we discuss sufficient conditions for global (point) identification. Due to the nonlinearity of the IVQR problem, these conditions are stronger than the usual first stage assumptions in linear instrumental variables models and require the instrument to have a nontrivial impact on the joint distribution of $(Y,D)$. The identification conditions are particularly easy to interpret when $D$ and $Z$ are binary; see Chernozhukov and Hansen (2005, Section 2.4) and Appendix C.2.1 for a further discussion. ∎

The main implication of Assumption 1 is the following conditional moment restriction (Chernozhukov and Hansen, 2005, Theorem 1):

[TABLE]

In this paper, we focus on the commonly used linear-in-parameters model for $q(\cdot)$ (e.g., Chernozhukov and Hansen, 2006):

[TABLE]

where $\theta(\tau):=(\theta_{X}(\tau)^{\prime},\theta_{D}(\tau)^{\prime})^{\prime}\in\mathbb{R}^{d_{X}+d_{D}}$ is the finite dimensional parameter vector of interest. The conditional moment restriction (2.1) suggests GMM estimators based on the following unconditional population moment conditions:

[TABLE]

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 (2.1). Appendix A 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.,

[TABLE]

The resulting GMM objective function reads

[TABLE]

where $m_{i}\left(\theta\right):=\left(1\left\{Y_{i}\leq X_{i}^{\prime}\theta_{X}+D_{i}^{\prime}\theta_{D}\right\}-\tau\right)(X_{i}^{\prime},Z_{i}^{\prime})^{\prime}$ and $W_{N}\left(\theta\right)$ is a positive definite weighting matrix. Estimation based on (2.4) 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.

3. Decentralization

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 5. 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})^{\prime}$. For each $\theta\in\mathbb{R}^{d_{X}+d_{D}}$, define the following (weighted) QR objective functions:

[TABLE]

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 B for more details.

The players then solve the following optimization problems:

[TABLE]

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 (e.g., Koenker, 2017).

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

[TABLE]

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)^{\prime},\dots,\Psi_{P,J}(\theta)\right)^{\prime}$ is the set of unconditional IVQR moment conditions. Hence, $\theta^{*}$ satisfies

[TABLE]

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$.111While 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.222In Appendix C.2, 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.

Assumption 2**.**

The following conditions hold.

1. (1)

$\Theta$* is a closed rectangle in $\mathbb{R}^{d_{X}+d_{D}}$. $\theta^{*}$ is in the interior of $\Theta$.* 2. (2)

$E[|Z_{\ell}|^{2}]<\infty$* for $\ell=1,\dots,d_{D}$. $E[|X_{k}|^{2}]<\infty$ for all $k=1,\dots,d_{X}$. For each $\ell=1,\dots,d_{D}$, $D_{\ell}$ has a compact support;* 3. (3)

The conditional cdf $y\to F_{Y|D,X,Z}(y)$ is continuously differentiable for all $y\in\mathbb{R}$ $a.s.$ The conditional density $f_{Y|D,Z,X}$ is uniformly bounded $a.s.$; 4. (4)

For any $\theta\in\Theta$, the matrices

[TABLE]

and

[TABLE]

are positive definite.

For each $j$, let $\Theta_{-j}\subset\mathbb{R}^{d_{-j}}$ denote the parameter space for $\theta_{-j}$. Assumption 2.1 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 2.2 and 2.3 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 (3.2) is well-defined and convex (cf. Appendix B). The first part of Assumption 2.4 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 4 in the appendix). For the second part of Assumption 2.4, 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

[TABLE]

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})^{\prime}$ 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.

Lemma 1**.**

*** Suppose that Assumptions 1 and 2 hold. Then, there exist functions $L_{j}:R_{-j}\to\mathbb{R}^{d_{j}},j=1,\dots,J$ such that, for $j=1,\dots,J$,*

[TABLE]

Further, $L_{j}$ is continuously differentiable on the interior of $R_{-j}$ for all $j=1,\dots,J$.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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})$.333One 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

[TABLE]

which is a univariate function whose fixed point is often straightforward to compute.

Define

[TABLE]

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 C. Now define

[TABLE]

where the second equality follows because $\tilde{R}_{d_{D}}$ turns out to be a subset of $\tilde{R}_{j}$ for all $j\leq d_{D}$. The following corollary ensures that $K$ and $M$ are well-defined on $D_{K}$ and $D_{M}$ respectively.

Corollary 1**.**

The maps $K:D_{K}\to\mathbb{R}^{d_{X}+d_{D}}$ and $M:D_{M}\to\mathbb{R}^{d_{D}}$ exist and are continuously differentiable on the interior of their domains.

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.

Proposition 1**.**

Suppose that Assumptions 1 and 2 hold. Then,

1. (i)

$\Psi_{P}(\theta^{*})=0$* if and only if $K\left(\theta^{*}\right)=\theta^{*}$* 2. (ii)

$\Psi_{P}(\theta^{*})=0$* if and only if $M\left(\theta^{*}_{-1}\right)=\theta^{*}_{-1}$ and $\theta^{*}_{1}=L_{1}(\theta^{*}_{-1})$.*

In view of Proposition 1, 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^{*}$.

1. (1)

Simultaneous dynamical system:444This algorithm is akin to the Jacobi computational procedure.

[TABLE] 2. (2)

Sequential dynamical system:555Smyth (1996) considers this type of algorithm for $J=2$ and calls it “zigzag” algorithm. It is akin to a Gauss-Seidel procedure.

[TABLE]

where $\theta_{1}^{(s+1)}=L_{1}\left(\theta_{-1}^{(s)}\right)$.

These discrete dynamical systems constitute the basis for our estimation algorithms.666These discrete dynamical systems can also be viewed as learning dynamics in a game (Li and Basar, 1987; Fudenberg and Levine, 2007).

4. Population Algorithms

In this section, we explore the implications of the fixed point reformulation for constructing population-level algorithms for computing fixed points.

4.1. Contraction-based Algorithms

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 1, (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 1.

Assumption 3**.**

There exist open strictly convex sets $\tilde{D}_{K}\subseteq D_{K}$ and $\tilde{D}_{M}\subseteq D_{M}$ such that

1. (1)

$\|J_{K}\left(\theta\right)\|\leq\lambda$* for some $\lambda<1$ for all $\theta\in\tilde{D}_{K}$;* 2. (2)

$\|J_{M}\left(\theta_{-1}\right)\|\leq\lambda$* for some $\lambda<1$ for all $\theta_{-1}\in\tilde{D}_{M}$.*

Under this additional assumption, the iterative algorithms are guaranteed to converge to the fixed point. We summarize this result below.

Proposition 2**.**

Suppose that Assumptions 1 and 2 hold.

1. (i)

Suppose further that Assumption 3.1 holds. Then $K$ is a contraction on the closure of $\tilde{D}_{K}$. The fixed point $\theta^{*}\in cl(\tilde{D}_{K})$ of $K$ is unique. For any $\theta^{(0)}\in\tilde{D}_{K},$ the sequence $\{\theta^{(s)}\}_{s=0}^{\infty}$ defined in (3.15) satisfies $\theta^{(s)}\to\theta^{*}$ as $s\to\infty$. 2. (ii)

Suppose further that Assumption 3.2 holds. Then $M$ is a contraction on the closure of $\tilde{D}_{M}$. The fixed point $\theta^{*}_{-1}\in cl(\tilde{D}_{M})$ of $M$ is unique. For any $\theta^{(0)}_{-1}\in\tilde{D}_{M},$ the sequence $\{\theta^{(s)}_{-1}\}_{s=0}^{\infty}$ defined in (3.16) satisfies $\theta^{(s)}_{-1}\to\theta^{*}_{-1}$ as $s\to\infty$.

In the case of a single endogenous variable, the Jacobian matrices of $K$ and $M$ are given by

[TABLE]

where

[TABLE]

One may therefore check the high-level condition through the Jacobians of the original moment restrictions. In Appendix C.2.1, we illustrate a simple primitive condition for a local version of Assumption 3. In practice, we found that violations of Assumption 3 lead to explosive behavior of the estimation algorithms and, thus, are very easy to detect numerically.

4.2. Connections to the Identification Conditions in the Literature

In view of Proposition 1, identification of $\theta^{*}$ is equivalent to uniqueness of the fixed points of $K$ and $M$, which is ensured by Proposition 2. Here, we discuss how the conditions required by Proposition 2 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)\neq 0$ for all $\theta\neq\theta^{*}$ in the neighborhood. Local identification in the IVQR model follows from standard results (e.g., Rothenberg, 1971; Chen, Chernozhukov, Lee, and Newey, 2014). For example, if $\Psi_{P}(\theta)$ is differentiable, Chen, Chernozhukov, Lee, and Newey (2014, Section 2.1) 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 5.1 in the appendix, which is a local version of Assumption 3.1. Assumption 5.1 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)\neq 0$. Observe that, for any $\theta$,

[TABLE]

If $\partial\Psi_{P,j}(\theta)/\partial\theta_{j}^{\prime}|_{\theta=\theta^{*}}$ is invertible for $j=1,2$ (which is true under Assumption 2.4), $J_{\Psi_{P}}(\theta^{*})$ is full rank if and only if

[TABLE]

That is, it requires that none of the eigenvalues of $J_{K}$ has modulus one. Therefore, Assumption 5.1 is sufficient but not necessary for condition (4.3) to hold. Specifically, Assumption 5.1 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 (e.g. Galor, 2007).777The 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 2 to the global identification result in Chernozhukov and Hansen (2006).

Lemma 2** (Theorem 2 in Chernozhukov and Hansen (2006)).**

Suppose that Assumption 1 holds. Moreover, suppose that (i) $\Theta$ is compact and convex and $\theta^{*}$ is in the interior of $\Theta$; (ii) $f_{Y|D,Z,X}$ is uniformly bounded a.s.; (iii) $J_{\Psi_{P}}(\theta)$ is continuous and has full rank uniformly over $\Theta$; and (iv) the image of $\Theta$ under the mapping $\theta\to\Psi_{P}(\theta)$ is simply connected. Then, $\theta^{*}$ uniquely solves $\Psi_{P}(\theta)=0$ over $\Theta$.

Under Conditions (i)–(iv), which are substantially stronger than the local identification conditions discussed above, the result in Lemma 2 follows from an application of Hadamard’s global univalence theorem (e.g. Theorem 1.8 in Ambrosetti and Prodi (1995)).

Comparing Lemma 2 to Proposition 2, we can see that the result in Lemma 2 establishes identification over the whole parameter space $\Theta$, while Proposition 2 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 2 correspond to our Assumptions 2.1 and 2.3. Moreover, our Assumption 2.3 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, Chernozhukov and Hansen (2006) assume the simple connectedness of the image of $\Psi_{P}$ (Condition (iv)). By contrast, we use a different univalence theorem by Gale and Nikaido (1965) (applied to the map $\Xi$ defined in (E.3) 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 3). 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.

4.3. Root-Finding Algorithms and Nesting

Assumption 3 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 Chernozhukov and Hansen (2006); see Lemma 2.

Note that, for the simultaneous dynamical system, $\theta^{*}$ solves

[TABLE]

where $I_{{d_{X}+d_{D}}}$ is the identity map. Similarly, in the sequential dynamical system, $\theta^{*}_{-1}$ solves

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 C.3, we discuss the formal conditions required for the existence and uniqueness of such fixed points.888For an equilibrium of the sub-game to be well-defined, one may directly assume that Chernozhukov and Hansen (2006)’s global identification condition holds within the sub-game, given player 3’s action $\theta_{3}$. Alternatively, if Assumption 3 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 C.3.

This nesting procedure is generic and can be extended to more than three players by sequentially adding additional layers of nesting.999In 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.

5. Sample Estimation Algorithms

5.1. Sample Estimation Problem

Let $\left\{(Y_{i},D_{i}^{\prime},X_{i}^{\prime},Z_{i}^{\prime})\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

[TABLE]

For each $j=1,\dots,J$, let the sample BR function $\hat{L}_{j}(\theta_{-j})$ be a function such that

[TABLE]

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 (Koenker, 2018). 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}^{\prime}\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 (3.5)–(3.6); see Lemma 10 in the appendix.

Remark 5.1**.**

The proposed estimators rely on decentralizing the original non-smooth and non-convex IVQR GMM problem into a series of convex QR problems. The quality and the computational performance of our procedures therefore crucially depends on the choice of the underlying QR estimation approach, which deserves some further discussion. The interested reader is referred to Koenker (2017) for an excellent overview on the computational aspects of QR. In this paper, we use the Barrodale and Roberts algorithm which is implemented as the default in the quantreg package and described in detail in Koenker and D’Orey (1987, 1994). This algorithm is computationally tractable for problems up several thousand observations. For larger problems, we recommend using interior point methods, potentially after preprocessing; see Portnoy and Koenker (1997) for a detailed description. These methods are conveniently implemented in the quantreg package. For very large problems, one can resort to first-order gradient descent methods, which are amenable to modern parallelized computation; see Section 5.5 in Koenker (2017) for an introduction and simulation evidence on the performance of such methods. ∎

We focus on estimators constructed based on the dynamical system $M$. Lemma 5 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$:

[TABLE]

where $\theta_{1}=\hat{L}_{1}\left(\theta_{-1}\right)$. This map induces a sample analog of the sequential dynamical system in Section 3.

1. Sample sequential dynamical system:

[TABLE]

where $\theta_{1}^{(s+1)}=\hat{L}_{1}\left(\theta_{-1}^{(s)}\right)$.

5.2. Contraction-based Algorithms

The first set of algorithms exploits that, under Assumption 3, $\hat{M}$ is a contraction mapping with probability approaching one. In this case, we iterate the dynamical system or (5.6) until $\|\theta_{-1}^{(s)}-\hat{M}(\theta_{-1}^{(s)})\|$ is within a numerical tolerance $e_{N}$.101010In 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$.

5.3. Algorithms based on Root-Finders and Optimizers

We construct an estimator $\hat{\theta}_{N}$ of $\theta^{*}$ as an approximate fixed point to the sample problem:

[TABLE]

where $\hat{\theta}_{N,1}=\hat{L}_{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 (Brent, 1971) 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 4.3, 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.111111See Section 8 and Appendix D.

5.4. Profiling Algorithms

A key insight underlying our estimation algorithms is that, given $\theta_{-1}$, the estimation problem becomes a standard convex QR problem:

[TABLE]

This insight suggests a profiling estimator.121212We 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

[TABLE]

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 4.3 and 5.3.

Relative to the root-finding methods described in Section 5.3, the profiling estimator has the advantage that evaluating $f$ only requires estimating one QR. On the other hand, the root-finding algorithms in Section 5.3 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 3).

6. Asymptotic Theory

6.1. Estimators

We define an estimator $\hat{\theta}_{N}$ of $\theta^{*}$ as an approximate fixed point of $\hat{M}$ in the following sense:

[TABLE]

An estimator of $\theta^{*}_{1}$ can be constructed by setting

[TABLE]

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 5 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 (6.1)–(6.2) using the contraction algorithm in Section 5.2 with tolerance $e_{N}=o(N^{-1/2}).$ One may then apply the result of Dominitz and Sherman (2005) to obtain the root-$N$ consistency of the estimator.131313Satisfying $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 and Sherman (2005). For completeness, this result is summarized in Appendix H.

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.141414The key conditions for these results, uniform convergence (in probability) of $\hat{K}$ and its stochastic equicontinuity, are established in Lemma 11. 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 and Sherman (2005). 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.

Remark 6.1**.**

To establish the asymptotic properties, one could try to reformulate our estimator as an estimator that approximately solves a GMM problem. Here, instead of relying on another reformulation, which would require establishing a sample analog version of Proposition 1, we develop and directly apply an asymptotic theory for fixed point estimators. We take this approach as the theory itself contains generic results (Theorem 1 and Lemmas 6–7) surrounding the Hadamard-differentiability of fixed points, which allow for applying the functional $\delta$-method to obtain the asymptotic distribution of $\hat{\theta}_{N}$ and bootstrap validity. These results can potentially be used to analyze decentralized estimators outside the IVQR class. ∎

6.2. Asymptotic Theory and Bootstrap Validity

The following theorem gives the limiting distribution of our estimator. For each $w=(y,d^{\prime},x^{\prime},z^{\prime})^{\prime}$ and $\theta\in\Theta,$ let $f(w;\theta)\in\mathbb{R}^{d_{X}+d_{D}}$ be a vector whose sub-vectors are given by

[TABLE]

and let $g(w;\theta)=(g_{1}(w;\theta)^{\prime},\dots,g_{J}(w;\theta))^{\prime}$ be a vector such that

[TABLE]

Theorem 1**.**

Suppose that Assumptions 1 and 2 hold. Let $\left\{W_{i}\right\}_{i=1}^{N}$ be an i.i.d. sample generated from the IVQR model, where $W_{i}=(Y_{i},D_{i}^{\prime},X_{i}^{\prime},Z_{i}^{\prime})$. Then,

[TABLE]

with

[TABLE]

where $\mathbb{W}$ is a tight Gaussian process in $\ell^{\infty}(\Theta)^{d_{X}+d_{D}}$ with mean zero and the covariance kernel

[TABLE]

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

[TABLE]

Let $\Psi(\tau)=(X^{\prime},Z^{\prime})^{\prime}.$ As shown in Chernozhukov and Hansen (2006) (Theorem 3 and Remark 3), $\sqrt{N}(\tilde{\theta}_{N}-\theta^{*})$ converges weakly to a mean zero multivariate normal distribution with variance

[TABLE]

where $J_{\Psi_{P}}(\theta^{*})=E\left[f_{\varepsilon(\tau)|X,D,Z}(0)\Psi(\tau)(X^{\prime},D^{\prime})\right],$ $\varepsilon(\tau)=Y-X^{\prime}\theta^{*}_{1}-D^{\prime}\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.

Corollary 2**.**

Suppose the conditions of Theorem 1 hold. Then, $V=\tilde{V}$.

To conduct inference on $\theta^{*}$, one may employ a natural bootstrap procedure. For this, use in (6.1)–(6.2) 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.

1. (1)

Compute the fixed point estimator $\hat{\theta}_{N}$ using the original sample. 2. (2)

Draw a bootstrap sample $\{W^{*}_{i}\}_{i=1}^{N}$ randomly with replacement from $P_{N}.$ Use the sequential dynamical system based on $\hat{M}^{*}$ combined with a contraction or root-finding algorithm to compute $\hat{\theta}_{N}^{*}$. 3. (3)

Repeat Step 2 across bootstrap replications $b=1,\dots,B$. Let

[TABLE]

Use $F_{B}$ as an approximation to the sampling distribution of the root $\sqrt{N}(\hat{\theta}_{N}-\theta^{*})$.

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 Chernozhukov and Hansen (2006) 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.151515The 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 and Bickel, 2000), is an interesting question which we leave for future research.

The following theorem establishes the consistency of the bootstrap procedure. For this, let $\stackrel{{\scriptstyle L^{*}}}{{\leadsto}}$ denote the weak convergence of the bootstrap law in outer probability, conditional on the sample path $\{W_{i}\}_{i=1}^{\infty}.$

Theorem 2**.**

Suppose that Assumptions 1 and 2 hold. Let $\left\{W_{i}\right\}_{i=1}^{N}$ be an i.i.d. sample generated from the IVQR model. Then,

[TABLE]

where $V$ is as in (6.5).

7. Empirical Example

In this section, we illustrate the proposed estimators by reanalyzing the effect of 401(k) plans on savings behavior as in Chernozhukov and Hansen (2004). This empirical example constitutes the basis for our Monte Carlo simulations in Section 8. As explained by Chernozhukov and Hansen (2004), 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. Chernozhukov and Hansen (2004) 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 Poterba, Venti, and Wise (1994, 1995, 1998) and Benjamin (2003), 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 (e.g., Engen, Gale, and Scholz, 1996).

We use the same dataset as in Chernozhukov and Hansen (2004). The dataset contains information about 9913 observations from a sample of households from the 1991 Survey of Income and Program Participation.161616The dataset analyzed by Chernozhukov and Hansen (2004) has 9,915 observations. Here we delete the two observations with negative income. We refer to Chernozhukov and Hansen (2004) 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.171717Chernozhukov and Hansen (2004) also consider total wealth and net non-financial assets.

We consider the following linear model for the conditional potential outcome quantiles

[TABLE]

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.

1. (1)

For any given $\theta_{2}$, compute $\theta_{1}=\hat{L}_{1}(\theta_{2})$ by running an ordinary QR (for the $\tau$-th quantile) using $Y_{i}-D_{i}^{\star}\theta_{2}$ as the dependent variable and $X_{i}$ as a regressor vector. 2. (2)

Given $\theta_{1}=\hat{L}_{1}(\theta_{2})$ from the previous step, compute $\hat{L}_{2}(\theta_{1})$ by running an ordinary QR with weights $Z_{i}/D_{i}^{\star}$ using $Y_{i}-X_{i}^{\prime}\theta_{1}$ as the dependent variable, $D_{i}^{\star}$ as the regressor, and omitting the constant.

Combining these two steps yields the sequential response map $\hat{M}(\cdot)=\hat{L}_{2}(\hat{L}_{1}(\cdot))$. Figure 1 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 1 further provides a straightforward graphical way to check the validity of the sample analog of Assumption 3. 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.

1. (3)

Find the fixed point of $\hat{M}$ using one of the following algorithms.

• •

Contraction algorithm: Set the tolerance $e_{N}$.181818In our application and simulations, we set $e_{N}$ equal to the square-root of the machine precision in R: $e_{N}\approx 1.5\cdot 10^{-8}$. We set the initial value $\theta_{2}^{(0)}$ to the 2SLS estimate of the coefficient on $D^{\star}_{i}$. Iterate $\theta_{2}^{(s+1)}=\hat{M}(\theta_{2}^{(s)}),$ for $s=0,1,2,\dots,$ until $|\theta_{2}^{(s+1)}-\theta_{2}^{(s)}|\leq e_{N}$. Report $\theta_{2}^{(s+1)}$ as the estimate $\hat{\theta}_{D}$ of $\theta_{D}$.

• •

Root-finding algorithm: Find the solution to the equation $\theta_{2}-\hat{M}(\theta_{2})=0$ by applying Brent’s method as implemented by the R-package uniroot. Report the solution as the estimate $\hat{\theta}_{D}$ of $\theta_{D}$. 2. (4)

(Optional) Compute $\hat{\theta}_{X}=\hat{L}_{1}(\hat{\theta}_{D})$ (as in Step (2)) if one is also interested in the coefficients on $X$. Recover estimates of the original coefficients following the discussion in Appendix B.

We compare our estimators to the IQR estimator of Chernozhukov and Hansen (2006) 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 8, 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 2 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 3 depicts pointwise 95% confidence intervals for the proposed estimators obtained using the empirical bootstrap described in Section 6.2 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.

8. Empirical Monte Carlo Study

In this section, we assess the practical performance of our estimation algorithms in an empirical Monte Carlo study based on the application in Section 7.

8.1. An Application-based DGP

We consider DGPs which are based on the empirical application of Section 7.191919The construction of our DGPs is inspired by the application-based DGPs in Kaplan and Sun (2017). 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

[TABLE]

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 2. $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:

[TABLE]

where we set $\theta_{D,2}(U_{i})=10000$. The second endogenous variable is generated as

[TABLE]

and the second instrument is generated as $Z_{2,i}\sim N(0,1)$. We set $N=9913$ as in the empirical application.

8.2. Estimation Algorithms

We assess and compare several different algorithms, all of which are based on the dynamical system $\hat{M}$.202020We do not explore algorithms based on $\hat{K}$ because, as discussed in Section 5, 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).212121We 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 1 presents an overview of the algorithms.

8.3. Results

Here we describe the computational performance of the different procedures and the finite sample performance of our bootstrap inference procedure. In Appendix D.1, 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 D.2 presents additional simulation evidence, demonstrating that our algorithms perform well and remain computationally tractable with more than two endogenous regressors.222222Specifically, 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 2–3 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.232323Note 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 3 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.242424Our 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, Hansen, and Wüthrich (2017), 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 4 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.

9. Conclusion

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.252525Andrews and Mikusheva (2016) 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.

**APPENDIX (FOR ONLINE PUBLICATION)

**

Hiroaki Kaido and Kaspar Wüthrich

Table of Contents

\startcontents

[sections] \printcontents[sections]l1

Appendix A Overidentification and Efficiency

In the main text, we focus on just-identified moment restrictions with $d_{Z}=d_{D}$, for which the construction of an estimator is straightforward. If the model is overidentified (i.e. if $d_{Z}>d_{D}$), we can transform the original moment conditions

[TABLE]

into a set of just-identified moment conditions

[TABLE]

where $\tilde{Z}$ is a $d_{D}\times 1$ vector of transformations of $(X,Z)$. A practical choice is to construct $\tilde{Z}$ using a least squares projection of $D$ on $Z$ and $X$.

To achieve pointwise (in $\tau$) efficiency, we can employ the following two-step procedure which is based on Remark 5 in Chernozhukov and Hansen (2006):

Step 1: We first obtain an initial consistent estimate of $\theta^{*}$ using one of our estimators based on a set of just-identified moment conditions such as (A.1). We then use nonparametric estimators to estimate the conditional densities $V(\tau)=f_{\varepsilon(\tau)|X,Z}(0)$ and $v(\tau)=f_{\varepsilon(\tau)|D,X,Z}(0)$, where $\varepsilon(\tau)=Y-X^{\prime}\theta^{*}_{1}(\tau)-D_{1}\theta^{*}_{2}(\tau)-\dots-D_{d_{D}}\theta^{*}_{J}(\tau)$, and the conditional expectation function $E_{P}\left[Dv(\tau)\mid X,Z\right]$.

Step 2: We apply our procedure to obtain a solution to the following moment conditions:

[TABLE]

Consider players $j=1,\dots,J$ solving the following optimization problems:

[TABLE]

where

[TABLE]

and $E_{P}\left[Dv(\tau)\mid X,Z\right]_{j-1}$ is the $j$-th element of $E_{P}\left[Dv(\tau)\mid X,Z\right]$. For each $j,$ the BR function $L_{j}(\theta_{-j}(\tau))$, defined as a member of the set of minimizers of $Q_{P,j}(\cdot,\theta_{-j})$, solves a suitable subset of the moment conditions in (A.2). The optimization problems in (A.3)-(A.4) are convex population QR problems provided that the model is parametrized such that $E_{P}\left[Dv(\tau)\mid X,Z\right]_{j-1}/D_{j-1}$, $j=2,\dots,J$, is positive. Estimation can then proceed by replacing the population QR problems by their sample analogues and applying one of the estimation algorithms discussed in the main text. By Corollary 2, the resulting estimator is asymptotically equivalent to a GMM estimator that uses the optimal instrumental variables and thus achieves pointwise (in $\tau$) efficiency (e.g., Chamberlain, 1987).262626Corollary 2 can be applied to the current setting by replacing the original set of covariates and instruments $\Psi(\tau)=(X^{\prime},Z^{\prime})^{\prime}$ in (6.8) with the optimal instrumental variables $\Psi(\tau)=(V(\tau)^{\prime},E_{P}[D_{1}v(\tau)\mid X,Z]/D_{1},\dots,E_{P}[D_{d_{D}}v(\tau)\mid X,Z]/D_{d_{D}})^{\prime}$.

Appendix B Reparametrization

In the main text, we assume that the model is reparametrized such that $Z_{\ell}/D_{\ell}$ is positive for all $\ell=1,\dots,d_{D}$. This ensures that the weights are well-defined and that the weighted QR problems are convex. However, in empirical applications, the weights may not be well-defined (e.g., if $D_{\ell}$ is an indicator variable with $P(D_{\ell}=0)>0$) or negative in some instances. Assuming that $Z_{\ell}$ is positive, a simple way to reparametrize the model is to add a large enough constant $c$ to $D_{\ell}$.272727Since the IVQR model is characterized by conditional moments (as in (2.1)), one may choose transformations of instruments to generate unconditional moment conditions. In case $Z_{\ell}$ is not positive $a.s.$, one can use a positive transformation (e.g., a logistic function) of $Z_{\ell}$ instead of $Z_{\ell}$ itself. The decentralization and identification results then hold with the transformed instruments as long as they satisfy our assumptions. This transformation is theoretically justified by the compactness of the support of $D_{\ell}$ (Assumption 2.2). To fix ideas, suppose that one is interested in estimating the following linear-in-parameters model with a single endogenous variable:

[TABLE]

where $\theta_{1}=(\theta_{11},\theta_{12}^{\prime})^{\prime}$ and $X=\left(1,\tilde{X}^{\prime}\right)^{\prime}$. Suppose further that the support of $D$ is a compact interval, $\left[d_{\min},d_{\max}\right]\subset\mathbb{R}$, with $d_{\min}<0$. In this case, we can apply the transformation $D^{\star}=D+c$, where $c>|d_{\min}|$. The transformed model reads

[TABLE]

where $\theta^{\star}_{11}=\theta_{11}-c\theta_{2}$. Importantly, one can always back out the original parameters, $\theta=(\theta_{11},\theta_{12}^{\prime},\theta_{2})^{\prime}$, from the parameters in the reparametrized model, $\theta^{\star}=(\theta_{11}^{\star},\theta_{12}^{\prime},\theta_{2})^{\prime}$.

Appendix C Decentralization

C.1. The Domains of $M_{j}$-Maps

Recall that, in (3.13), we defined the set

[TABLE]

Similarly, for $k=2,\dots,d_{D}-1$, define

[TABLE]

For $k=d_{D}$, let

[TABLE]

For each $k$, the map $M_{k}$ is well-defined on $\tilde{R}_{k}$. Note also that $\tilde{R}_{d_{D}}\subset\tilde{R}_{j}$ for all $j\leq d_{D}$.

C.2. Local Decentralization and Local Contractions

We say that an estimation problem admits local decentralization if the BR functions $L_{j}$, $j=1,\dots,J$, and the maps $K$ and $M$ are well-defined over a local neighborhood of $\theta^{*}$. The following weak conditions are sufficient for local decentralization of the IVQR estimation problem.

Assumption 4**.**

The following conditions hold.

1. (1)

The conditional cdf $y\mapsto F_{Y|D,X,Z}(y)$ is continuously differentiable at $y^{*}=d^{\prime}\theta^{*}_{-1}+x^{\prime}\theta^{*}_{1}$ for almost all $(d,x,z)$. The conditional density $f_{Y|D,Z,X}$ is bounded on a neighborhood of $y^{*}$ $a.s.$; 2. (2)

The matrices

[TABLE]

and

[TABLE]

are positive definite.

Assumption 4 is weaker than Assumption 2.3–2.4 as it only requires the conditions, e.g., continuous differentiability of the conditional CDF, at a particular point, e.g. $y^{*}$. Under this condition, we can study the local properties of our population algorithms. For this, the following lemma ensures that the BR maps are well-defined locally.

Lemma 3**.**

Suppose that Assumptions 1, 2.1–2.2, and 4 hold. Then, there exist open neighborhoods $\mathcal{N}_{L_{-j}},j=1,\dots,J$, $\mathcal{N}_{K}$, $\mathcal{N}_{M}$ of $\theta^{*}_{-j}$, $\theta^{*}$, and $\theta^{*}_{-1}$ such that

1. (i)

There exist maps $L_{j}:\mathcal{N}_{-j}\to\mathbb{R}^{d_{j}}$, $j=1,\dots,J$ such that, for $j=1,\dots,J$,

[TABLE]

Further, $L_{j}$ is continuously differentiable for all $j=1,\dots,J$. 2. (ii)

The maps $K:\mathcal{N}_{K}\to\mathbb{R}^{{d_{X}+d_{D}}}$ and $M:\mathcal{N}_{M}\to\mathbb{R}^{d_{D}}$ are continuously differentiable.

Proof.

(i) The proof is similar to that of Lemma 1 (see Appendix E). Therefore, we sketch the argument below for $j=1$. By Assumptions 2.2 and 4.1, $\Psi_{P,1}$ is continuously differentiable on a neighborhood $V$ of $\theta^{*}$. By Assumption 4.2 and the continuity of $\det(\partial\Psi_{P,1}(\theta)/\partial\theta_{1}^{\prime})$, one may choose $V$ so that $\det(\partial\Psi_{P,1}(\theta)/\partial\theta_{1}^{\prime})\neq 0$ for all $\theta=(\theta_{1},\theta_{-1})\in V$. By the implicit function theorem, there is a continuously differentiable function $L_{1}$ and an open set $\mathcal{N}_{-1}$ containing $\theta_{-1}$ such that

[TABLE]

The arguments for $L_{j}$, $j\neq 1$ are similar.

(ii) Let $\mathcal{N}_{K}=\{\theta\in\Theta:\pi_{-j}\theta\in\mathcal{N}_{-j},~{}j=1,\dots,J\}$ and let $\mathcal{N}_{M}$ be defined by mimicking (3.13), while replacing $\Theta_{j}$ with $\mathcal{N}_{j}$ in the definition of $\tilde{R}_{j}$ for $j=1,\dots,J$. The continuous differentiability of $K$ and $M$ follows from that of $L_{j}$, $j=1,\dots,J$. ∎

C.2.1. Local Contractions

Recall that $\rho(A)$ denotes the spectral radius of a square matrix $A$. The following assumption ensures that $K$ and $M$ are local contractions.

Assumption 5**.**

---

1. (1)

$\rho\left(J_{K}\left(\theta^{*}\right)\right)<1$** 2. (2)

$\rho\left(J_{M}\left(\theta^{*}_{2}\right)\right)<1$**

Here, we illustrate a primitive condition for Assumption 5. Consider a simple setup without covariates (i.e. $X=1$), a binary $D$, and a binary $Z$. We only analyze Assumption 5.1. A similar result can be derived for Assumption 5.2. In this setting, the Jacobian of $K$ evaluated at $\theta^{*}$ is given by

[TABLE]

The characteristic polynomial is then given by

[TABLE]

Hence, Assumption 3.1 holds if all eigenvalues (i.e. the roots $\lambda_{K}$ of $p_{K}(\lambda)=0$) have modulus less than one, which holds when

[TABLE]

This condition can be simplified to

[TABLE]

where $f_{Y|d,z}(y):=f_{Y|D=d,Z=z}(y)$ and $p(d|z):=P(D=d\mid Z=z)$. It is instructive to interpret condition (C.1) under the local average treatment effects framework of Imbens and Angrist (1994). Specifically, condition (C.1) holds if (i) their monotonicity assumption is such that there are compliers but no defiers and (ii) the complier potential outcome density functions are strictly positive. Conversely, the condition is violated if there are defiers but no compliers.

Under the local contraction conditions in Assumption 5, we have the following results.

Proposition 3**.**

Suppose that Assumptions 1, 2.1, 2.2, and 4 hold.

1. (i)

Suppose further that Assumption 5.1 holds. Then there exists a closed neighborhood $\bar{\mathcal{N}}_{K}$ of $\theta^{*}$ such that $K(\bar{\mathcal{N}}_{K})\subset\bar{\mathcal{N}}_{K}$ and $K$ is a contraction on $\bar{\mathcal{N}}_{K}$ with respect to an adapted norm. 2. (ii)

Suppose further that Assumption 5.2 holds. Then there exists a closed neighborhood $\bar{\mathcal{N}}_{M}$ of $\theta^{*}_{2}$ such that $M(\bar{\mathcal{N}}_{M})\subset\bar{\mathcal{N}}_{M}$ and $M$ is a contraction on $\bar{\mathcal{N}}_{M}$ with respect to an adapted norm.

Proof.

We only prove the result for $K$, the proof for $M$ is similar. By Lemma 3, $L_{j}$ is continuously differentiable at $\theta^{*}$. Note that $J_{K}$ is given by

[TABLE]

which is continuous at $\theta^{*}$. The desired result now follows, for instance, from Proposition 2.2.19 in Hasselblatt and Katok (2003). ∎

C.3. Nested Algorithms: Existence and Uniqueness of Fixed Points in Sub-games

Here we discuss two different sets of conditions which ensure that the nested fixed point algorithms in Section 4.3 are well-defined. Specifically, we present conditions for the existence and uniqueness of fixed points in the sub-games. Section C.3.1 considers contraction-based identification conditions. In Section C.3.2, we briefly discuss global identification conditions. To illustrate, we consider the case with three players ($J=3$).

C.3.1. Contraction-based Identification.

Suppose that Assumptions 1–3 hold. Consider a sub-game formed by players 1 and 2 given some $\tilde{\theta}_{3}$. We assume that $\tilde{\theta}_{3}$ is chosen so that $\mathsf{D}_{M_{1,2|3}}$ (defined below) is non-empty. Let the moment conditions for the sub-game be defined as

[TABLE]

Note that $\mathsf{\Psi}_{P,1}$ and $\mathsf{\Psi}_{P,2}$ and other objects below depend on $\tilde{\theta}_{3}$. We will often suppress this dependence to alleviate the notation. Moreover, define

[TABLE]

Assumptions 1–2 and Lemma 1 guarantee that BR functions $\mathsf{L}_{1}$ and $\mathsf{L}_{2}$, where

[TABLE]

are well-defined. The $M$ map for the sub-system is

[TABLE]

This map exists and is well-defined on

[TABLE]

Note that, by Assumption 3, for any $(\theta_{2},\theta_{3})\in\tilde{D}_{M}$,

[TABLE]

Given these definitions, we can now investigate the sub-system. Observe that the derivative $J_{M_{1,2|3}}(\theta_{2})$ of $M_{1,2|3}$ with respect to $\theta_{2}$ is a component of $J_{M}.$ In particular, for any $(\theta_{2},\theta_{3})$, $J_{M}$ may be written as

[TABLE]

Let $V_{12}:=\{x\in\mathbb{R}^{2}:x=(x_{1},0)^{\prime},x_{1}\in\mathbb{R}\}$. Then, by the definition of the operator norm (see e.g. Bhatia, 2013),

[TABLE]

Hence, it follows that

[TABLE]

which is an analog of Assumption 3 for the sub-game. Then, Proposition 2 ensures the existence and uniqueness of the fixed point in the sub-game. In this section, we focused on identification conditions based on the dynamical system $M$. Similar arguments can be used to establish identification based on the dynamical system $K$.

C.3.2. Global Identification Conditions.

To ensure the existence and uniqueness of the fixed point in the sub-game, we can alternatively rely on the existing global identification conditions in Chernozhukov and Hansen (2006) (cf. Section 4.2 and Lemma 2) and Proposition 1. For every value $\tilde{\theta}_{3}\in\Theta_{3}$, this requires analogues of the conditions in Lemma 2 to hold for the sub-game between players 1 and 2.

Appendix D Additional Simulation Results

D.1. Bias and RMSE Application-based DGP

Here we report simulation evidence on the finite sample bias and RMSE of the different IVQR estimators based on the application-based DGPs in Section 8. Tables 5–6 present the results. We find that all the proposed algorithms perform well and exhibit comparable bias and RMSE properties. In particular, the finite sample performances of our preferred estimators are comparable to IQR and the profiling estimator, which shows that their computational advantages do not come at a cost in terms of the finite sample performance.

D.2. Three Endogenous Variables

Here we present additional simulation evidence with three endogenous variables. We consider the application-based DGP of Section 8, augmented with an additional endogenous variable:

[TABLE]

where $\theta_{D,3}(U_{i})=10000$, $D_{3,i}=0.8\cdot Z_{3,i}+0.2\cdot\Phi^{-1}(U_{i})$, and $Z_{3,i}\sim N(0,1)$. We only report the results based on the contraction algorithm and the nested fixed point algorithm. We do not report results for IQR, which we found to be computationally prohibitive with three endogenous regressors. Table 7 shows that both methods exhibit similar performances in terms of bias and RMSE, which are comparable to their respective performances with two endogenous regressors. Table 8 displays average computation times. As expected, the computational advantages of the contraction algorithm relative to the nested fixed point algorithm are more pronounced than with two endogenous variables.

D.3. Additional Simulations Simple Location Scale DGP

This section presents some additional simulation evidence based on the following location-scale shift model:

[TABLE]

Here $D_{1,i}$ and $D_{2,i}$ are the endogenous variables of interest and $X_{i}$ is an exogenous covariate. In addition, we have access to two instruments $Z_{1,i}$ and $Z_{2,i}$. For $\gamma_{2}=\gamma_{4}=\gamma_{7}=0$, this model reduces to the model considered in Section 6.1 of Andrews and Mikusheva (2016). We set $\gamma_{1}=\dots=\gamma_{7}=1$. To evaluate the performance of our algorithms with one endogenous variable, we set $\gamma_{4}=\gamma_{7}=0$ and use $Z_{1i}$ as the instrument. Following Andrews and Mikusheva (2016), we consider a symmetric as well as an asymmetric DGP:

[TABLE]

where $\left(\xi_{U,i},\xi_{D_{1},i},\xi_{D_{2},i},\xi_{Z_{1},i},\xi_{Z_{2},i},\xi_{X,i}\right)$ is a Gaussian vector with mean zero, all variances are set equal to one, $Cov(\xi_{U},\xi_{D_{1}})=Cov(\xi_{U},\xi_{D_{2}})=0.5$, $Cov(\xi_{D_{1}},\xi_{Z_{1}})=0.8$, $Cov(\xi_{D_{2}},\xi_{Z_{2}})=0.4$, which allows us to investigate the impact of instrument strength, all other covariances are equal to zero, and $\Phi$ is the cumulative distribution function of the standard normal distribution.

We first investigate the bias and RMSE of the different methods. Tables 9–12 present the results. With one endogenous variable, the performances of the root-finding algorithm using Brent’s method, the profiling estimator, and IQR are similar both in terms of bias and RMSE. The contraction algorithm performs well, but exhibits some bias at the tail quantiles. Turning to the results with two endogenous variables, we can see that the nested algorithm exhibits the best overall performance, both in terms of bias and RMSE. The performances of the SA-based optimization algorithm, IQR, and the profiling estimator are similar and only slightly worse than that of the nested algorithm. The contraction algorithm tends to exhibit some bias at the tail quantiles. However, this bias decreases substantially as the sample size gets larger. Finally, comparing the results for the coefficients on $D_{1}$ and $D_{2}$, we can see that the instrument strength matters for the performance of all estimators (including IQR), suggesting that weak identification can have implications for the estimation of IVQR models.

Table 13 displays the empirical coverage probabilities of the bootstrap confidence intervals. The results show that the our bootstrap procedure exhibits excellent size properties.

Appendix E Proofs of Theoretical Results in Section 3

Proof of Lemma 1.

(i) We first show that $L_{1}$ is well-defined. For a given $\theta_{-1}\in\mathbb{R}^{d_{D}}$, let $\theta_{1}^{*}\in\arg\min_{\tilde{\theta}_{1}\in\mathbb{R}^{d_{X}}}Q_{P,1}(\tilde{\theta}_{1},\theta_{-1})$. Under Assumption 2, the objective function is convex and differentiable with respect to $\tilde{\theta}_{1}$. Therefore, by the necessary and sufficient condition of minimization, $\theta_{1}^{*}$ solves

[TABLE]

In what follows, we show that the map $L_{1}:\theta_{-1}\mapsto\theta^{*}_{1}$ is well-defined on $R_{-1}$ using a global inverse function theorem. Recall that

[TABLE]

This function is continuously differentiable with respect to $\theta$. The Jacobian is given by

[TABLE]

where the second equality follows from Assumption 2 and the dominated convergence theorem. Define a transform $\Xi:\Theta\to\mathbb{R}^{d_{X}+d_{D}}$ by

[TABLE]

We follow Krantz and Parks (2012) (Section 3.3) to obtain an implicit function $L_{1}$ on a suitable domain such that $\theta_{1}=L_{1}(\theta_{2})$ if and only if $\Psi_{P,1}(\theta)=0.$ The key is to apply a global inverse function theorem to $\Xi$. Toward this end, we analyze the Jacobian of $\Xi$, which is given as

[TABLE]

where, for any $d\in\mathbb{N}$, $I_{d}$ denotes the $d\times d$ identity matrix. Let $I\subset\{1,\dots d_{X}+d_{D}\}$.

For any matrix $A$, let $[A]_{I,I}$ denote a principal minor of $A$, which collects the rows and columns of $A$ whose indices belong to the index set $I$. By (E.4), if $I\subset\{1,\dots,d_{1}\}$,

[TABLE]

for a subvector $\tilde{X}$ of $X$, which is positive definite by Assumption 2 and Lemma 4. If $I\subset\{d_{1}+1,\dots,{d_{X}+d_{D}}\}$, $[J_{\Xi}(\theta)]_{I,I}=I_{\ell}$ for some $1\leq\ell\leq d_{D}$ and is hence positive definite. Otherwise, any principal minor is of the following form:

[TABLE]

for some subvector $\tilde{X}$ of $X$ and a $m\times\ell$ matrix $B.$ Note that

[TABLE]

where the last inequality follows again from Assumption 2 and Lemma 4. Hence, $J_{\Xi}(\theta)$ is a $P$-matrix. Note that $\Theta$ is a closed rectangle. By Theorem 4 in Gale and Nikaido (1965), $\Xi$ is univalent, and hence the inverse map $\Xi^{-1}$ is well defined.

Let

[TABLE]

which coincides with the definition in (3.8) with $j=1$. Let $F_{1}=[I_{d_{1}},0_{d_{1}\times d_{-1}}]$. For each $\theta_{-1}\in R_{-1}$, define

[TABLE]

Then, for any $\theta\in\Theta$, $\Psi_{P,1}(\theta)=0$ if and only if $\theta_{-1}\in R_{-1}$ and $\Xi(\theta)=(0,\theta_{-1})$. By the univalence of $\Xi$, this is true if and only if $\theta=\Xi^{-1}(0,\theta_{-1})$, and the first $d_{1}$ components extracted by applying $F_{1}$ is $\theta_{1}$. This ensures $L_{1}$ is well-defined on $R_{-1}.$

Below, for any set $A,$ let $A^{o}$ denote the interior of $A$. Let $R_{-1}^{o}=\{\theta_{-1}\in\mathbb{R}^{d_{-1}}:(0,\theta_{-1})\in\Xi(\Theta^{o})\}.$ Note that $\Psi_{P,1}$ is $\mathcal{C}^{1}$ on $\Theta^{o}$ and, for each $\theta=(\theta_{1},\theta_{-1})\in\Theta$ with $\theta_{-1}\in R_{d_{-1}}^{o}$, $\det(\partial\Psi_{P,1}(\theta)/\partial\theta_{1}^{\prime})\neq 0$. Therefore, by the implicit function theorem, there is a $\mathcal{C}^{1}$-function $\tilde{L}_{1}$ and an open set $V$ containing $\theta_{-1}$ such that

[TABLE]

However, such a local implicit function must coincide with the unique global map $L_{1}$ on $V$. Hence, $L_{1}|_{V}=\tilde{L}_{1}$, and therefore $L_{1}$ is continuously differentiable at $\theta_{-1}$. Since the choice of $\theta_{-1}$ is arbitrary, $L_{1}$ is continuously differentiable for all $\theta_{-1}\in R_{2}^{o}.$

Showing that the conclusion holds for any other $L_{j}$ for $j=2,\dots,J$ is similar, and hence we omit the proof. ∎

Lemma 4**.**

Suppose $E_{P}[f_{Y|D,X,Z}\left(D^{\prime}\theta_{-1}+X^{\prime}\theta_{1}\right)XX^{\prime}]$ is positive definite. Then, for any subvector $\tilde{X}$ of $X$ with dimension $\tilde{d}_{X}\leq d_{X}$, $E_{P}[f_{Y|D,X,Z}\left(D^{\prime}\theta_{-1}+X^{\prime}\theta_{1}\right)\tilde{X}\tilde{X}^{\prime}]$ is positive definite.

Proof.

In what follows, let $W=f_{Y|D,X,Z}\left(D^{\prime}\theta_{-1}+X^{\prime}\theta_{1}\right)$ and let

[TABLE]

Let $\tilde{X}$ be a subvector of $X$ with $\tilde{d}_{X}$ components. Then, there exists a $d_{X}\times d_{X}$ permutation matrix $P_{\pi}$ such that the first $\tilde{d}_{X}$ components of $P_{\pi}X$ is $\tilde{X}$.

Let $B:=E[WP_{\pi}XX^{\prime}P_{\pi}^{\prime}]$ and note that

[TABLE]

by the linearity of the expectation operator and $W$ being a scalar. Let $\lambda$ be an eigenvalue of $B$ such that

[TABLE]

for the corresponding eigenvector $z\in\mathbb{R}^{d_{X}}$. By (E.9)-(E.10),

[TABLE]

Note that $P_{\pi}^{-1}=P_{\pi}^{\prime}$ due to $P_{\pi}$ being a permutation matrix. Letting $y:=P_{\pi}^{\prime}z$ then yields

[TABLE]

which in turn shows that $\lambda$ is an eigenvalue of $A$. For any eigenvalue of $A$, the argument above can be reversed to show that it is also an eigenvalue of $B$. Since the choice of the eigenvalue is arbitrary, $A$ and $B$ share the same eigenvalues.

Now let $C:=E[W\tilde{X}\tilde{X}^{\prime}]$ and note that it is a leading principal submatrix of $B$. Then, by the eigenvalue inclusion principle (Horn and Johnson, 1990, Theorem 4.3.28),

[TABLE]

where the last inequality follows from the positive definiteness of $A$. This completes the claim of the lemma. ∎

Proof of Corollary 1.

The existence of $K$ and its continuous differentiability follows immediately from Lemma 1. For $M$, by the definition of $\tilde{R}_{1}$, for any $\theta_{-1}\in\tilde{R}_{j}$, there exists $(\theta_{1},\theta_{2})\in\Theta_{1}\times\Theta_{2}$ such that

[TABLE]

By (i), one may then write $\theta_{1}=L_{1}(\theta_{-1})$ and $\theta_{2}=L_{2}(L_{1}(\theta_{-1}),\pi_{-\{1,2\}}\theta_{-1})$. Hence, the map $M_{1}:\tilde{R}_{1}\to\Theta_{2}$ below is well-defined:

[TABLE]

Recursively, arguing in the same way, the maps

[TABLE]

are well-defined on $\tilde{R}_{2},\dots,\tilde{R}_{d_{D}}$ respectively. The continuous differentiability of $M$ follows from that of $L_{j}$s and the chain rule. ∎

Proof of Proposition 1.

$\Longrightarrow$: For every solution, $\Psi_{P}(\theta^{*})=0$, $\theta^{*}_{j}=L_{j}\left(\theta^{*}_{-j}\right)$ by construction under Assumptions 1 and 2. It follows that $K\left(\theta^{*}\right)=\theta^{*}$ and $M\left(\theta^{*}_{-1}\right)=\theta^{*}_{-1}$.

$\Longleftarrow$: For the simultaneous response note that $K\left(\bar{\theta}\right)=\bar{\theta}$ implies that $\bar{\theta}_{j}=L_{j}\left(\bar{\theta}_{-j}\right)$ for all $j\in\{1,\dots,J\}$. Thus, $\bar{\theta}$ solves $\Psi_{P}(\bar{\theta})=0$ by Lemma 1. Consider next the sequential response. Let $\tilde{\theta},\bar{\theta}\in\Theta$ be such that $\tilde{\theta}_{j}=L_{j}(\bar{\theta}_{-j})$ for $j=1,\dots,J$. By Lemma 1, they satisfy

[TABLE]

Thus, a fixed point $\tilde{\theta}=\bar{\theta}$ satisfies $\Psi_{P}\left(\bar{\theta}\right)=0$. ∎

Appendix F Proofs of Theoretical Results in Section 4

Proof of Proposition 2.

We prove the result for $K$. By Assumption 3, there exists a strictly convex set $\tilde{D}_{K}$ on which the spectral norm of the Jacobian of $K$ is uniformly bounded by 1. This ensures that $K$ is a contraction map on $cl(\tilde{D}_{K})$, and the claim of the proposition now follows from Theorem 2.2.16 in Hasselblatt and Katok (2003). ∎

Appendix G Proofs of Theoretical Results in Section 6

To state and prove results in a concise manner, we use the population and sample simultaneous response maps $K$ and $\hat{K}$ below to define our estimand $\theta^{*}$ and estimator $\hat{\theta}_{N}$. Namely, $\theta^{*}$ is the fixed point of $K$, and $\hat{\theta}_{N}$ solves

[TABLE]

Note that the fixed point estimator defined in (6.1)-(6.2) is asymptotically equivalent to the estimator above due to Lemma 5.

Proof of Theorem 1.

Let $H:=I_{d_{X}+d_{D}}-K$. A fixed point $\theta^{*}$ of $K$ then satisfies

[TABLE]

Similarly, let $\hat{H}:=I_{d_{X}+d_{D}}-\hat{K}$. The estimator $\hat{\theta}_{N}$ satisfies

[TABLE]

where $r_{N}=o_{p}(N^{-1/2})$. Let $\varphi:\ell^{\infty}(\Theta)^{d_{X}+d_{D}}\times\mathbb{R}\to\mathbb{R}^{d_{X}+d_{D}}$ be a map such that, for each $(H,r)\in\ell^{\infty}(\Theta)^{d_{X}+d_{D}}\times\mathbb{R}$, $\tilde{\theta}=\varphi(H,r)$ is an $r$-approximate solution, which satisfies

[TABLE]

One may then write

[TABLE]

By Lemma 12, $\sqrt{N}(\hat{K}-K)\leadsto\mathbb{W}$ in $\ell^{\infty}(\Theta)^{d_{X}+d_{D}}$, where $\mathbb{W}$ is a Gaussian process defined in Lemma 12. By Lemmas 6-7, Condition $Z$ in Chernozhukov, Fernandez-Val, and Melly (2013)(CFM henceforth) holds, which in turn ensures that one may apply Lemmas E.2 and E.3 in CFM. This ensures

[TABLE]

Hence, we obtain (6.4) with

[TABLE]

Finally, note that $\dot{H}_{\theta^{*}}=I_{d_{X}+d_{D}}-J_{K}(\theta^{*})$ by Lemma 7. This establishes the theorem. ∎

Proof of Theorem 2.

Recall that $\hat{H}=I_{d_{X}+d_{D}}-\hat{K}$. The estimator $\hat{\theta}_{N}$ satisfies

[TABLE]

where $r_{N}=o_{p}(N^{-1/2})$. Similarly, let $\hat{H}^{*}=I_{d_{X}+d_{D}}-\hat{K}^{*}$. Let $P^{*}$ denote the law of $\hat{H}^{*}$ conditional on $\{W_{i}\}_{i=1}^{\infty}.$ The bootstrap estimator $\hat{\theta}_{N}^{*}$ satisfies

[TABLE]

where $r^{*}_{N}=o_{P^{*}}(N^{-1/2})$ conditional on $\{W_{i}\}_{i=1}^{\infty}$.

Using the $r$-approximation, one may therefore write

[TABLE]

Let $E_{P^{*}}$ denote the conditional expectation with respect to $P^{*}$. Let $BL_{1}$ denote the space of bounded Lipschitz functions on $\mathbb{R}^{d_{X}+d_{D}}$ with Lipschitz constant 1. Then, for any $\epsilon>0$,

[TABLE]

By Lemma 12, $\sqrt{N}(\hat{H}^{*}-\hat{H})=-\sqrt{N}(\hat{K}^{*}-\hat{K})\stackrel{{\scriptstyle L^{*}}}{{\leadsto}}-\mathbb{W}\stackrel{{\scriptstyle d}}{{=}}\mathbb{W}$. Noting that $h\circ\varphi^{\prime}_{H,0}\in BL_{1}(\ell^{\infty}(\Theta)\times\mathbb{R})$ and $r_{N}=o_{p}(N^{-1/2})$, it follows that

[TABLE]

with probability approaching 1 due to $r_{N}=o_{P}(N^{-1/2})$. Hence, for the conclusion of the theorem, it suffices to show that the right hand side of (G.10) tends to 0 in probability.

For this, as shown in the proof of Theorem 1, $\varphi$ is Hadamard differentiable at $(H,0)$. Hence, by Theorem 3.9.4 in Van der Vaart and Wellner (1996),

[TABLE]

Take the difference of the left and right hand sides respectively and note that $\varphi^{\prime}_{H,0}$ is linear. This implies the right hand side of (G.10) tends to 0 in probability. This ensures

[TABLE]

∎

Proof of Corollary 2.

Note that $V$ and $g$ may be written as

[TABLE]

where $R(\theta^{*})$ is a $d_{X}+d_{D}$-by-$d_{X}+d_{D}$ matrix given by

[TABLE]

Further, by Lemma 1 and the form of $J_{L_{-j}}(\theta^{*}_{j})$ given in (4.2),

[TABLE]

The form of $R(\theta^{*})$ and $J_{K}(\theta^{*})$ imply

[TABLE]

where $J_{\Psi_{P}}$ is the Jacobian of the estimating equations. Eq. (6.5)-(G.14) and (G.19) ensure that one may also write

[TABLE]

As shown in Chernozhukov and Hansen (2006), the Jacobian of $\Psi_{P}$ is given by

[TABLE]

where $\Psi(\tau)=(X^{\prime},Z^{\prime})^{\prime}.$ Furthermore,

[TABLE]

Hence, (G.20)-(G.22) show that $V$ coincides with the asymptotic variance of the estimator that solves the estimating equations in (6.7). ∎

Lemma 5**.**

Suppose Assumptions 1-2 hold. (i) Let $\hat{\theta}_{N}$ be an estimator of $\theta^{*}$ that satisfies (G.23). Then, it also satisfies (6.1)-(6.2); (ii) Let $\hat{\theta}_{N}$ be an estimator of $\theta^{*}$ that satisfies (6.1)-(6.2). Then, it also satisfies

[TABLE]

Proof.

(i) Consider the case $j=2$. Note that, by (G.23),

[TABLE]

where $r_{N,1}=o_{p}(N^{-1/2})$, and the second equality follows from the definition of $\hat{\theta}_{N,2}$. (G.25) can be written as

[TABLE]

where the last equality follows from the stochastic equicontinuity of $\mathcal{L}_{N}$ shown in the proof of Lemma 11 and $L_{2}$ being Lipschitz since $L_{2}$ is continuously differentiable with a derivative that is uniformly bounded on the compact set $\Theta$. By (G.24)-(G.26), it holds that $\hat{\theta}_{N,j}=M_{j}(\hat{\theta}_{N,-1})+o_{p}(N^{-1/2})$ for $j=2$. Repeat the same argument sequentially for $j=3,\dots,J$. The first conclusion of the lemma then follows.

(ii) Suppose now that $r_{N,1}:=\hat{\theta}_{N,1}-\hat{L}_{1}(\hat{\theta}_{N,-1})\neq o_{P}(N^{-1/2}).$ Then, there is a subsequence $k_{N}$ along which, for any $\eta>0$, $\sqrt{k}_{N}r_{k_{N},1}>\eta$ for all $k_{N}$ with positive probability. Then, the $O_{P}(r_{k_{N},1})$-term in (G.26) is not $o_{p}(k_{N}^{-1/2})$, which therefore implies $\hat{\theta}_{N,j}\neq M_{j}(\hat{\theta}_{N,-1})+o_{p}(N^{-1/2})$ for $j=2$. The second conclusion of the lemma then follows. ∎

Lemma 6**.**

Let $\Lambda\subset\mathbb{R}^{p}$ be a compact set, and let $K:\Lambda\to\mathbb{R}^{p}$ be a map that has a unique fixed point $\lambda_{0}\in\Lambda$. let $H:\Lambda\to\mathbb{R}^{p}$ be defined by $H(\lambda):=\lambda-K(\lambda)$. Then $H^{-1}(x)=\{\lambda\in\Lambda:H(\lambda)=x\}$ is continuous at $x=0$ in Hausdorff distance.

Proof.

For any $x$, write

[TABLE]

Let $x_{n}\to 0$. Since $\lambda_{0}$ is the unique fixed point of $K$, $H^{-1}(0)=\{\lambda_{0}\}$. Therefore,

[TABLE]

Hence, it suffices to show that $\sup_{\lambda\in H^{-1}(x_{n})}\|\lambda-\lambda_{0}\|=o(1)$. We show this by contradiction. Suppose that there is a sequence $\{\lambda_{n}\}\subset\Lambda$ and $\delta>0$ such that $\lambda_{n}\in H^{-1}(x_{n})$ for all $n$ and $\{\lambda_{n}\}$ has a subsequence $\{\lambda_{k_{n}}\}$ such that $\|\lambda_{k_{n}}-\lambda_{0}\|>\delta$ for all $n$. $\lambda_{k_{n}}\in\Lambda$ is a sequence in a compact space, and hence there is a further subsequence $\lambda_{h_{n}}$ such that $\lambda_{h_{n}}\to\lambda^{*}$ for some $\lambda^{*}\in\Lambda$ with $\lambda^{*}\neq\lambda_{0}$. By the continuity of $K$, one then has

[TABLE]

By $\lambda_{h_{n}}-K(\lambda_{h_{n}})=x_{n}$ and $x_{n}\to 0$, it must hold that

[TABLE]

However this contradicts the fact that $\lambda_{0}$ is the unique fixed point, and hence the conclusion follows. ∎

Lemma 7**.**

Suppose $H=I-K$ and $K:\mathbb{R}^{p}\to\mathbb{R}^{p}$ is continuously differentiable at $\lambda_{0}$. Suppose further that $\det(I-J_{K}(\lambda_{0}))\neq 0$. Let $\dot{H}_{\lambda_{0}}:=I-J_{K}(\lambda_{0})$. Then,

[TABLE]

and

[TABLE]

Proof.

Let $\{h_{n}\}\subset\mathbb{S}^{p}$ be a sequence on the unit sphere. Then,

[TABLE]

where $\bar{\lambda}_{n}$ is a mean value between $\lambda_{0}+th_{n}$ and $\lambda_{0}$. Therefore, by the Cauchy-Schwarz inequality,

[TABLE]

where we used $\|h_{n}\|=1$, $\bar{\lambda}_{n}\to\lambda_{0}$, and the continuity of the Jacobian.

For the second claim, note that

[TABLE]

and $h\mapsto\|(I-J_{K}(\lambda_{0}))h\|$ is continuous. Since the domain of $h$ is compact, there is $h^{*}\in\mathbb{S}^{p}$ such that $\inf_{\|h\|=1}\|\dot{H}_{\lambda_{0}}h\|=\|(I-J_{K}(\lambda_{0}))h^{*}\|.$ Let $q=(I-J_{K}(\lambda_{0}))h^{*}$ and note that $I-J_{K}(\lambda_{0})$ is linearly independent (due to $\det(I-J_{K}(\lambda_{0}))\neq 0$), and hence $q\neq 0$. Hence $\inf_{\|h\|=1}\|\dot{H}_{\lambda_{0}}h\|=\|q\|>0$. Hence, the second conclusion follows. ∎

The following result is a slight extension of Lemma E.1 in CFM.

Lemma 8**.**

Suppose that $\Lambda\subset\mathbb{R}^{p}$ and $\mathcal{U}$ is a compact and convex set in $\mathbb{R}^{q}$. Let $\mathcal{I}$ be an open set containing $\mathcal{U}$. Suppose that $\Psi:\Lambda\times\mathcal{I}\to\mathbb{R}^{p}$ is continuous and $\lambda\mapsto\Psi(\lambda,u)$ is the gradient of a convex function in $\lambda$ for each $u\in\mathcal{U}$; (b) for each $u\in\mathcal{U}$, $\Psi(\lambda_{0}(u),u)=0$; (c) $\frac{\partial}{\partial(\lambda^{\prime},u^{\prime})}\Psi(\lambda,u)$ exists at $(\lambda_{0}(u),u)$ and is continuous at $(\lambda_{0}(u),u)$ for each $u\in\mathcal{U}$ and $\dot{\Psi}_{\lambda_{0}(u),u}:=\frac{\partial}{\partial\lambda^{\prime}}\Psi(\lambda,u)|_{\lambda_{0}(u)}$ obeys $\inf_{u\in\mathcal{U}}\inf_{\|h\|=1}\|\dot{\Psi}_{\lambda_{0}(u),u}h\|>c_{0}>0$. Then, Condition $Z$ in CFM holds and $u\mapsto\lambda_{0}(u)$ is continuously differentiable with derivative $J_{\lambda_{0}}(u)=-\dot{\Psi}^{-1}_{\lambda_{0}(u)u}\frac{\partial}{\partial u^{\prime}}\Psi(\lambda_{0}(u),u)$.

Proof.

The proof is the same as that of Lemma E.1 in CFM, in which $\mathcal{U}$ is a compact interval in $\mathbb{R}$. A slight modification is needed when one computes the derivative of $\lambda_{0}(u)$ with respect to $u$. Since $u$ is allowed to be multidimensional, the implicit function theorem gives

[TABLE]

which is uniformly bounded and continuous in $u$ by condition (c), which ensures continuous differentiability of $u\mapsto\lambda_{0}(u)$. Note that for any $\delta>0$ and $\lambda\in B_{\delta}(\lambda_{0}(u))$, there is $\eta>0$ and $u^{\prime}$ such that $\|u^{\prime}-u\|\leq\eta$ so that

[TABLE]

Since $\mathcal{U}$ is compact (and hence totally bounded), there is a finite set $\{u_{j}\}_{j=1}^{J}\subset\mathcal{U}$ such that $\mathcal{U}\subset\bigcup_{j}B_{\eta}(u_{j})$. The argument above then shows that $\mathcal{N}=\bigcup_{u\in\mathcal{U}}B_{\delta}(\lambda_{0}(u))\subset\bigcup_{j}B_{2\delta}(\lambda_{0}(u_{j}))$, which ensures that $\mathcal{N}$ is totally bounded. Since $\mathcal{N}$ is a subset of a Euclidean space (equipped with a complete metric), it follows that $\mathcal{N}$ is compact. This ensures condition $Z$ (i) in CFM. The rest of the proof is essentially the same as the case, in which $\mathcal{U}$ being a compact interval. ∎

Lemma 9**.**

Suppose Assumption 2 holds. Let $w=(y,d^{\prime},x^{\prime},z^{\prime})$ and let $\tau\in(0,1)$. Define

[TABLE]

Then, $\mathcal{M}$ is a Donsker-class.

Proof.

The proof is standard, and hence we give a brief sketch for the first component of $f$, $f_{1}(w;\theta)=(1\{y\leq d^{\prime}\theta_{-1}+x^{\prime}\theta_{1}\}-\tau)x$. Note that $w\mapsto 1\{y\leq d^{\prime}\theta_{-1}+x^{\prime}\theta_{1}\}-\tau$ belongs to Type I-class in Andrews (1994), and the map $w\mapsto x$ does not depend on the parameter. By Theorems 2 and 3 in Andrews (1994), this function then satisfies the uniform entropy condition with the envelope function $\bar{M}(w)=x$, which is square integrable by assumption. Similar arguments apply to the other components of $f$. By Theorem 1 in Andrews (1994), the empirical process: $\mathbb{G}_{n}f$ is stochastically equicontinuous, and $\mathbb{G}_{n}f(\cdot,\theta)$ obeys the classical central limit theorem for each $\theta\in\Theta$. Hence, we conclude that $\mathcal{M}$ is Donsker. ∎

Below, let $g(w;\theta)=(g_{1}(w;\theta)^{\prime},\dots,g_{J}(w;\theta))^{\prime}$ be a vector such that

[TABLE]

Let \rho(\theta,\tilde{\theta}):=\big{\|}\text{diag}\big{(}E_{P}\big{[}(g(W;\theta)-E_{P}[g(W;\theta)])(g(w;\tilde{\theta})-E_{P}[g(w;\tilde{\theta})])^{\prime}\big{]}\big{)}\big{\|} be the variance semimetric. Let $W_{i}=(Y_{i},D_{i}^{\prime},X_{i}^{\prime},Z_{i}^{\prime}),i=1,\dots,N$ be an i.i.d. sample generated from the IVQR model. Define

[TABLE]

Similarly, let $W_{i}^{*}=(Y_{i}^{*},D_{i}^{*\prime},X^{*\prime}_{i},Z^{*\prime}_{i})^{\prime},i=1,\dots,N$ be an bootstrap sample from the empirical distribution $\hat{P}_{N}$ of $\{W_{i}\}$. Define

[TABLE]

where $\hat{L}^{*}_{j}$ is the sample best response map of player $j$, which is defined as in (5.3)-(5.4) while replacing $W_{i}$ with the bootstrap sample $W_{i}^{*}$ in (5.1)-(5.2).

Lemma 10 below shows that the sample BR functions approximately solve sample estimating equations and Lemma 11 characterizes the limiting distributions of $\mathcal{L}_{N}$ and $\mathcal{L}_{N}^{*}$.

Lemma 10**.**

Let the sample BR functions be $\hat{L}_{j}(\theta_{-j})\in\textrm{argmin}_{\tilde{\theta}_{j}}Q_{N,j}(\tilde{\theta}_{j},\theta_{-j}),j=1,\dots,J$. Let $\hat{L}_{j}(\theta_{-j})^{*}$ be an analog of $\hat{L}_{j}(\theta_{-j})$ for the bootstrap sample. Then, (i) the sample BR functions satisfy

[TABLE]

and

[TABLE]

where $\sup_{\theta_{-j}\in\Theta_{-j}}|r_{N,j}(\theta_{-j})|=o_{P}(N^{-1/2})$ for all $j$; (ii) the sample BR functions $\hat{L}_{j}^{*}(\theta_{-j}),j=1,\dots,J$ in the bootstrap sample satisfy (G.33)-(G.34) while replacing $(Y_{i},D_{i},X_{i},Z_{i})$ with a bootstrap sample $(Y_{i}^{*},D_{i}^{*},X_{i}^{*},Z_{i}^{*})$, each $\hat{L}_{j}$ with $\hat{L}_{j}^{*}$, and each $r_{N,j}$ with $r^{*}_{N,j}$ such that $\sup_{\theta_{-j}\in\Theta_{-j}}|r^{*}_{N,j}(\theta_{-j})|=o_{P^{*}}(N^{-1/2})$.

Proof.

(i) For $j\geq 2,$ the subgradient of $Q_{N,j}$ is

[TABLE]

and hence by the property of the subgradient, for any $v\in\mathbb{R}$, one has

[TABLE]

where $\nabla_{\theta_{j}}Q_{N,j}(\hat{L}_{j}(\theta_{-j}),\theta_{-j},v)$ is the directional derivative of $Q_{N,j}(\theta_{j},\theta_{-j})$ with respect to $\theta_{j}$ toward direction $v\in\mathbb{R}$ evaluated at $(\hat{L}_{j}(\theta_{-j}),\theta_{-j}).$ Note that the directional derivative is given by

[TABLE]

where

[TABLE]

Observe that $-\nabla_{\theta_{j}}Q_{N,j}(\hat{L}_{j}(\theta_{-j}),\theta_{-j},-v)\leq\xi v\leq\nabla_{\theta_{j}}Q_{N,j}(\hat{L}_{j}(\theta_{-j}),\theta_{-j},v)$. This implies

[TABLE]

Noting that $\sum_{i=1}^{N}1\{Y_{i}=(X^{\prime}_{i},D^{\prime}_{i,-(j-1)})^{\prime}\theta_{-j}+D_{i,j}^{\prime}\hat{L}_{j}(\theta_{-j})\}=\text{dim}(\theta_{j})=1$ and taking $v=1$, we obtain,

[TABLE]

uniformly in $\theta_{-j}$, where the last equality is due to $E[|Z_{i,j-1}|^{2}]<\infty$ by Assumption 2.2. Therefore, for some $r_{N,j}$ satisfying the assumption of the lemma, we may write

[TABLE]

The proof for $j=1$ is similar. Also, (ii) can be shown by mimicking the argument above. ∎

Lemma 11**.**

Suppose that Assumptions 1 and 2 hold. Then, (i) $\mathcal{L}_{N}:=(\mathcal{L}_{N,1},\dots,\mathcal{L}_{N,J})$ defined in (G.31) satisfies

[TABLE]

where $\mathbb{W}$ is a tight Gaussian process in $\ell^{\infty}(\Theta)^{d_{X}+d_{D}}$ with the covariance kernel

[TABLE]

$\mathcal{L}_{N}$* is stochastically equicontinuous with respect to the variance semimetric $\rho$; (ii) $\mathcal{L}^{*}_{N}:=(\mathcal{L}^{*}_{N,1},\dots,\mathcal{L}^{*}_{N,J})$ satisfies*

[TABLE]

(iii) $\rho$ satisfies $\lim_{\delta\downarrow 0}\sup_{\|\theta-\tilde{\theta}\|<\delta}\rho(\theta,\tilde{\theta})\to 0.$

Proof.

(i) We first work with $\mathcal{L}_{N,1}$. For this, we establish that $L_{1}$ is Hadamard differentiable. Note that $\theta_{1}=L_{1}(\theta_{-1})$ solves

[TABLE]

Take $\mathcal{U}=\Theta_{-1}$, $\Xi=\Theta_{1}$, $\psi(\lambda,u)=E_{P}[(1\{Y\leq Du+X^{\prime}\lambda\}-\tau)X]$. Define $\phi:\ell^{\infty}(\Xi\times\mathcal{U})^{k_{b}}\times\ell^{\infty}(\mathcal{U})\to\ell^{\infty}(\mathcal{U})$, which maps $(\psi,r)$ to a solution $\phi(\psi,r)=\lambda(\cdot)$ such that

[TABLE]

Then, one may write $L_{1}(\cdot)=\phi(\psi,0).$ We then show that $\psi$ satisfies the conditions of Lemma 8. Note first that $\mathcal{U}$ and $\Xi$ are compact. $\psi$ is continuous and $\lambda\mapsto\psi(\lambda,u)$ is the gradient of the convex function $\lambda\mapsto E_{P}\left[\rho_{\tau}(Y-Du-X^{\prime}\lambda)\right].$ The function $L_{1}(u)=\lambda_{0}(u)$ is defined as the exact solution of $\psi(\lambda,u)=0$. Note also that, by Assumption 2,

[TABLE]

where the second equality follows from the dominated convergence theorem, and the last display is well-defined by the square integrability of $X$. Similarly,

[TABLE]

Hence, the derivative

[TABLE]

exists and is continuous by Assumption 2. By Assumption 2.4, $\dot{\Psi}_{\lambda_{0}(u),u}=\frac{\partial^{2}}{\partial\theta_{1}\partial\theta_{1}^{\prime}}Q_{P,1}(L_{1}(\theta_{-1}),\theta_{-1})$ obeys

[TABLE]

Then, by Lemma 8 and Lemma E.2 in CFM, $\phi$ is Hadamard differentiable tangentially to $\mathcal{C}(\mathcal{N}\times\mathcal{U})^{K}\times\{0\}$ with the Hadamard derivative (of $L_{1}$)

[TABLE]

where $(z,0)\mapsto\phi^{\prime}_{\Psi,0}(z,0)$ is continuous over $z\in\ell^{\infty}(\Theta)^{K}$.

For $j\geq 2$, the argument is similar. For example, for $j=2$, one may take $\mathcal{U}=\Theta_{-2}$, $\Xi=\Theta_{2}$ and $\psi(\lambda,u)=E_{P}[(1\{Y\leq D_{2}\theta_{2}+(D_{1},X)^{\prime}u\}-\tau)Z_{2}]$ and write $L_{2}(\cdot)=\phi(\psi,0)$. The rest of the argument is the same.

Continuing with $j=1$, by Lemma 10, one may write $\hat{L}_{j}(\cdot)=\phi(\psi_{N},r_{N,1})$ with $\psi_{N}(\lambda,u)=\frac{1}{N}\sum_{i=1}^{N}1\{Y_{i}\leq D^{\prime}_{i}u+X^{\prime}_{i}\lambda\}X_{i}$ and $\sup_{\theta_{-1}\in\Theta_{-1}}|r_{N,1}(\theta_{-1})|=o_{p}(N^{-1/2}).$ By Lemma 9 and applying the $\delta$-method (as in Lemma E.3 in CFM), we obtain

[TABLE]

where $\mathbb{W}=(\mathbb{W}_{1}^{\prime},\dots,\mathbb{W}_{J}^{\prime})^{\prime}$ is a tight Gaussian process in $\ell^{\infty}(\Theta)^{d_{X}+d_{D}}$, where for each $j$, $\mathbb{W}_{j}\in\ell^{\infty}(\Theta_{-j})^{d_{j}}$ is given pointwise by

[TABLE]

Hence, its covariance kernel is as given in (G.43). By Lemma 1.3.8. in Van der Vaart and Wellner (1996), $\{\mathcal{L}_{N}\}$ is asymptotically tight, which in turn means that $\{\mathcal{L}_{N}\}$ is stochastically equicontinuous with respect to $\rho$ by Theorem 1.5.7 in Van der Vaart and Wellner (1996).

(ii) For each $j$, let $\mathcal{L}^{*}_{N,j}\in\ell^{\infty}(\Theta_{-j})^{d_{j}}$ be defined pointwise by

[TABLE]

Below, again we work with the case $j=1$. Using $\phi$ (the solution to (G.46)) and applying Lemma 10, we may write

[TABLE]

where $\hat{\psi}_{N}(\lambda,u)=N^{-1}\sum_{i=1}^{N}(1\{Y_{i}\leq D_{i}u+X_{i}^{\prime}\lambda\}-\tau)X_{i}$, and $\hat{\psi}_{N}^{*}$ is defined similarly for the bootstrap sample. Let $E_{P^{*}}$ denote the conditional expectation with respect to $P^{*}$, the law of $\{W^{*}_{i}\}_{i=1}^{N}$ conditional on the sample path. Let $BL_{1}$ denote the space of bounded Lipschitz functions on $\mathbb{R}^{d_{1}}$ with Lipschitz constant 1. Then, for any $\epsilon>0$,

[TABLE]

By Lemma 9 and Theorem 3.6.2 in Van der Vaart and Wellner (1996), $\sqrt{N}(\hat{\psi}^{*}_{N}-\hat{\psi}_{N})\stackrel{{\scriptstyle L^{*}}}{{\leadsto}}\mathbb{G}f_{1}$. Noting that $h\circ\phi^{\prime}_{\Psi,0}\in BL_{1}(\ell^{\infty}(\Theta_{-1})^{d_{1}}\times\mathbb{R})$ and $r_{N}=o_{p}(N^{-1/2})$, it follows that

[TABLE]

with probability approaching 1 due to $r_{N}=o_{P}(N^{-1/2})$. Hence, for the conclusion of the theorem, it suffices to show that the second term on the right hand side of (G.55) tends to 0.

As shown in the proof of (i), $\phi$ is Hadamard differentiable at $(\psi,0)$. Hence, by Theorem 3.9.4 in Van der Vaart and Wellner (1996),

[TABLE]

Take the difference of the left and right hand sides respectively and note that $\phi^{\prime}_{\Psi,0}$ is linear. This implies the right hand side of (G.55) tends to 0 in probability. This, together with (G.55)-(G.56), ensures

[TABLE]

where $\mathbb{W}_{1}(\theta_{-1})=-\frac{\partial^{2}}{\partial\theta_{1}\partial\theta_{1}^{\prime}}Q_{P,1}(L_{1}(\theta_{-1}),\theta_{-1})^{-1}\mathbb{G}f_{j}(\cdot;L_{1}(\theta_{-1}),\theta_{-1}).$ The analysis for any $j\neq 1$ is similar, and one may apply the arguments above jointly across $j=1,\dots,J$, which yields the second claim of the lemma.

(iii) Consider the first submatrix of $E_{P}[(g(W;\theta)-E_{P}[g(W;\theta)])(g(w;\tilde{\theta})-E_{P}[g(w;\tilde{\theta})])^{\prime}]$. It is given by

[TABLE]

Note that $\Theta$ is compact and \theta_{-1}\mapsto\big{(}\frac{\partial^{2}}{\partial\theta_{1}\partial\theta_{1}^{\prime}}Q_{P,1}(L_{1}(\theta_{-1}),\theta_{-1})\big{)}^{-1} is continuous by Lemma 1, which implies that this map is uniformly continuous. Therefore, it remains to show the uniform continuity of $\theta\mapsto\text{Var}(f_{1}(w;\theta))$. Note that

[TABLE]

The right hand side of the display above is continuous on the compact domain $\Theta$, and hence it is uniformly continuous. One can argue the same way for the other subcomponents of \text{diag}\big{(}E_{P}[(g(W;\theta)-E_{P}[g(W;\theta)])(g(w;\tilde{\theta})-E_{P}[g(w;\tilde{\theta})])^{\prime}]\big{)}. This completes the proof. ∎

Lemma 12**.**

Suppose that Assumptions 1 and 2 hold. (i) Let $W_{i}=(Y_{i},D_{i}^{\prime},X_{i}^{\prime},Z_{i}^{\prime})^{\prime},i=1,\dots,N$ be an i.i.d. sample generated from the IVQR model. Then,

[TABLE]

(ii) Let $W_{i}^{*}=(Y_{i}^{*},D_{i}^{*\prime},X^{*\prime}_{i},Z^{*\prime}_{i})^{\prime},i=1,\dots,N$ be an bootstrap sample from the empirical distribution $\hat{P}_{N}$ of $\{W_{i}\}_{i=1}^{N}$. Then,

[TABLE]

Proof.

(i) By Lemma 11, it follows that

[TABLE]

Note that, by the definition of $\hat{L}$ and $L$, one has

[TABLE]

The conclusion of the lemma then follows. The proof of (ii) is similar, and is therefore omitted. ∎

Appendix H Consistency of the Contraction Estimator

Below, we adopt the framework of Dominitz and Sherman (2005). Let $(\mathcal{X},d)$ be a metric space. For a contraction map $F:\mathcal{X}\to\mathcal{X}$, let $c_{F}$ be the modulus of contraction such that

[TABLE]

for any $x,x^{\prime}\in\mathcal{X}.$ As discussed in Section 5 the fixed point estimator $\hat{\theta}_{N}$ can be computed using the sample sequential dynamical system (in (5.6)) or the following sample simultaneous dynamical system.

[TABLE]

Lemma 13**.**

Suppose Assumptions 1, 2, and 3 hold. Let $\hat{\theta}_{N}$ be an estimator constructed by iterating the dynamical system in (H.1) or (in (5.6)) $s_{N}$ times, where $s_{N}\geq-\frac{1}{2}\ln N/\ln c_{K}$. Then,

[TABLE]

Proof.

We show the result by applying Theorem 1 in Dominitz and Sherman (2005) to the estimator obtained from the simultaneous dynamical system. The argument for the sequential system is similar.

By Assumption 3, $K$ is a contraction map on $D_{K}$. Let $\theta^{(s)}$ be obtained from iterating $s$-times the population dynamical system in (3.15). The iteration on the dynamical system is covergent at least linearly (Bertsekas and Tsitsiklis, 1989, Proposition 1.1). Under the condition on $s_{N}$, arguing as in (Dominitz and Sherman, 2005, p.842), it follows that $N^{1/2}\|\theta^{(s_{N})}-\theta^{*}\|\leq\|\theta^{(0)}-\theta^{*}\|$. Finally, by Lemma 12 and tightness of $\mathbb{W}$, $N^{1/2}\sup_{\theta\in D_{K}}\|\hat{K}(\theta)-K(\theta)\|=O_{p}(1)$. These imply the conditions of Theorem 1 in Dominitz and Sherman (2005) with $\delta=1/2$. The claim of the lemma then follows. ∎