Debiased/Double Machine Learning for Instrumental Variable Quantile Regressions
Jau-er Chen Chien-Hsun Huang Jia-Jyun Tien






§ INTRODUCTION


Machine learning methods have been actively studied in economic big data settings in recent years, 
cf. <cit.> and <cit.>. 
Most empirical studies in economics aim to understand the program evaluation, or equivalently, the causal effect. Constructing the counterfactual and then estimating causal effects relies on an appropriately chosen identification strategy. In economics, the instrumental variable approach is an extensively used identification strategy for causal inference. Therefore, the machine learning techniques often require the adaptation to exploit the structure of the underlying identification strategy. These adaptations are part of  an emerging research area at the intersection of machine learning and econometrics, which is called the causal machine learning in the economic literature.
Two popular causal machine learning approaches are currently available to estimate treatment effects through adapted machine learning algorithms, and they also provide valid standard errors of an estimated causal parameter of interest, such as the average treatment effect and quantile treatment effect. These two approaches are the double machine learning (DML) cf. <cit.>, and the generalized random forests (GRF) of <cit.>. The GRF estimates heterogeneous treatment effects and explores variable importance accounting for heterogeneity in the treatment effect. The resulting information is crucial for optimal polices mapping from individuals' observed characteristics to treatments. The DML provides a clever and general recipe for use of sample splitting, cross-fitting, and Neyman orthogonalization, to make causal inference possible and allows for almost any machine learner. Furthermore, the DML is feasible for dealing with high-dimensional datasets where researchers observe massive characteristics of the units. For instance, through sample splitting, the DML estimates each of the nuisance function (e.g., the expectations for the target variable and outcome variable given high-dimensional controls) on an auxiliary sample, and then it uses out-of-sample residuals as the basis for the treatment effect estimation. Moreover, the cross-fitting algorithm allows researchers to use all of the data in the final treatment effect estimation instead of throwing an auxiliary sample used in the sample splitting earlier. This procedure in fact follows the Neyman-type orthogonal moment conditions which ensure that the estimation above is insensitive to the first-order perturbations of the nuisance parameter near the true value, and consequently, the regular inference on a low-dimensional causal parameter proceeds. 


With the identification strategy of selection on observables (aka. unconfoundedness), 
empirical applications have been investigated by using the aforementioned two approaches, 
including the works by <cit.> and <cit.>. 
When it comes to the identification strategy of selection on unobservables, 
few empirical papers that use  causal machine learning can be found in the existing literature. 
Those empirical applications very often lack important observed control variables or involve reverse causality,
and thus,  researchers resort to the instrumental variable approach.
In this study, we investigate the estimation and inference of a low-dimensional causal parameter in the presence of high-dimensional controls in an instrumental variable quantile regression. In particular, we build on a previous study <cit.> and then further concretize the econometric procedure. To the best of our knowledge, this study is the first to investigate Monte Carlo performance and empirical studies based on the DML procedure within the framework of instrumental variable quantile regressions. We also make our R codes available on the GitHub repository[The R scripts conducting the estimation and inference of the Double Machine Learning for Instrumental Variable Quantile Regressions  can be found at https://github.com/FieldTien/DML-IVQR/tree/master/example] so that other researchers can benefit from the proposed estimation method.




<cit.> investigated the instrumental variable quantile regression in the context of GRF. 
Their econometric procedure yielded a measure of variable importance 
in terms of characterizing heterogeneity in the treatment effect. 
They proceeded by empirically investigating the distributional effect of 401(k) participation on net financial assets. They demonstrated that income, age, education, 
and family size are the first four important variables in explaining treatment effect heterogeneity. 
In contrast to our study, their GRF-based estimator is not designed for high-dimensional settings. 
With the same dataset, 
we also apply the proposed procedure to empirically investigate the
distributional effects of the 401(k) participation on net financial assets. 
Empirical results signify that 
the 401(k) participants with low savings propensity are more associated with the nonlinear income effect, 
which complements the findings in studies conducted by <cit.> and <cit.>. 
In addition, nonlinear transformations of the four aforementioned variables are also identified as important variables 
in the current context of DML-based instrumental variable quantile regression 
with high-dimensional observed characteristics. 


The rest of the paper is organized as follows. 
The model specification and practical algorithm are introduced in Section 2, which includes detailed descriptions of a general recipe for the DML. 
Section 3 presents finite-sample performances of the estimator through Monte Carlo experiments. 
Section 4 reinvestigates an empirical study on quantile treatment effects: the effect of 401(k) participation on wealth.
Section 5 concludes the paper. 


 






§ THE MODEL AND ALGORITHM



In this study, we use the instrumental variable quantile regression (IVQR) of <cit.> and <cit.> to identify the quantile treatment effect.
In Section 2.1, we briefly review the DML procedure developed in <cit.>.
In Section 2.2, we briefly review the conventional IVQR based on the exposition in <cit.>. 
In Section 2.3, we present DML-IVQR within the framework of high-dimensional controls.




 §.§ The Double Machine Learning



In this section, we briefly review the plain DML procedure.
Let us consider the following canonical example of estimating treatment effect α_0 in a partial linear regression under the identification strategy of selection on observables.

    Y = Dα_0 + h_0(X) + U,     𝔼[U|X, D]=0

where Y is the outcome variable, D is the target variable, and X is a high-dimensional vector of controls.
X are control variables in the sense that

    D= m_0(X) + V,

where m_0(·)≠ 0 and 𝔼[V|X]=0.
Note that h_0(X) and m_0(X) are nuisance functions because they are not the primary objects of interest.
<cit.> develop the DML procedure for estimating α_0, which is outlined in the following three steps.


  I. [Sample splitting]  Split the data into K random and roughly equally sized folds. For k=1, …, K,  a machine learner is used to fit the high-dimensional nuisance functions, 𝔼̂_(-k)[Y|X] and 𝔼̂_(-k)[D|X], using all data except for the kth fold.

  II. [Cross-fitting and residualizing]  Calculate out-of-sample residuals for these fitted nuisance functions on the kth fold; that is, Ŷ_(k) = Y_(k) - 𝔼̂_(-k)[Y|X] and D̂_(k) = D_(k) - 𝔼̂_(-k)[D|X].

  III. [Treatment effect estimation and inference]  Collect all of the out-of-sample residuals from the cross-fitting stage, and use the ordinary least squares to regress Ŷ on D̂ to obtain α̌, the estimator of α_0.
The resulting α̌ estimate can be paired with heteroskedastic consistent standard errors to obtain a confidence interval for the treatment effect.


Because estimating the nuisance functions through machine learners induces regularization biases,
the cross-fitting step was used to refrain from its biasing the treatment effect estimate.
The procedure is random due to the sample splitting. 
Different researchers with the same data set but making different random splits will obtain distinct estimators. This randomness can be reduced by using a larger value of K, but this increases computation cost. K≥ 10 is recommended.
In fact, the DML procedure follows a unified approach in terms of moment conditions and the Neyman orthogonality condition, cf. <cit.>.
In a nutshell, we seek to find moment conditions

    𝔼[g(Y, D, X, α_0, η_0 ) ] = 0

such that the following Neyman orthogonality condition holds

    .∂_η𝔼[g(Y, D, X, α_0, η_0 ) ]|_η=η_0 = 0,

where η_0 are nuisance functions with the true values.
Equation (4) is insensitive to the first-order perturbations of the nuisance function η near the true value. This property allows the estimation of  η_0 using regularized estimators (machine learners) η̂. Without this property, regularization may have too much effect on the estimator of α_0 for regular inference to proceed. 
The estimator α̌ of α_0 solves the empirical analog of the equation (3):

    1/n∑_i=1^n g(y_i, d_i, x_i, α̌, η̂)  = 0,

where we have plugged in the estimator η̂ for the nuisance function. 
Owing to the Neyman orthogonality property, the estimator is first-order equivalent to the infeasible estimator 
α̃ solving

    1/n∑_i=1^n g(y_i, d_i, x_i, α̃, η_0 )  = 0,

where we use the true value of η. 

Therefore, we recast the canonical example set by equations (1) and (2) into the moment conditions that guide the DML procedure outlined above.

    𝔼[g(Y, D, X, α_0, η_0 ) ] 
       =𝔼[[ (Y-E[Y|X] ) - (D-E[D|X]  )α_0 ] ×(D-E[D|X] )]
       =𝔼[[(Dα_0 + h_0(X) + U - (m_0(X)α_0 + h_0(X))  ) -Vα_0 ] × V]
       =𝔼[ (Dα_0 - m_0(X)α_0 + U -Vα_0 )× V ]
       =𝔼[ m_0(X)α_0 V + V^2α_0 - m_0(X)α_0 + UV - V^2α_0 ]
       =𝔼[  UV ] = 0,

where η_0 = [ 𝔼[Y|X]   𝔼[D|X] ].
It is easy to see that the corresponding Neyman orthogonality condition holds 

    .∂_η𝔼[g(Y, D, X, α_0, η_0 ) ]|_η=[ 𝔼[Y|X]   𝔼[D|X] ]
           = .∂_η𝔼[[ (Y-E[Y|X] ) - (D-E[D|X]  )α_0 ] ×(D-E[D|X] )]|_η=[ 𝔼[Y|X]   𝔼[D|X] ] =0.








 §.§ The Instrumental Variable Quantile Regression



Based on the exposition in <cit.>, the following conditional moment restriction yields an IVQR estimator:

    ℙ[Y ≤ q(τ,D,X)|X,Z ]= τ,

where q(·) is the structural quantile function, τ is the quantile index,
Y is the outcome variable, D is the target (endogenous) variable,
and X and Z are control variables and instruments, respectively.
Equation (1) and linear structural quantile specification lead to the following unconditional moment restriction

    𝔼[(τ - 1(Y - D'α - X'β≤ 0)Ψ]=0

where 
    Ψ:= Ψ(X, Z)
 is a vector of the  function of the instruments and control variables,
and (α', β')' are the unknown parameters. In particular, α is a causal parameter of interest.
The parameters depend on the quantile of interest, 
but we suppress  τ associated with α and β for simplicity of presentation.
Equation (2) leads to a particular moment condition for residualization. That is


    g_τ(α;β,δ) = (τ - 1(Y ≤ D'α + X'β))Ψ(α,δ(α))

with the instrument

    Ψ(α,δ(α)) := (Z - δ(α)X)


    δ(α) = M(α)J^-1(α),

where δ is a matrix parameter for weighting the least square Z on the X coefficient,

    M(α) = 𝔼[ZX'f_ε(0|X,Z)],   J(α) = 𝔼[XX'f_ε(0|X,Z)]

and f_ε(0|X,Z) is the conditional density of ϵ = Y - D'α - X'β(α)
with β(α)  defined by

    𝔼[(τ - 1(Y ≤ D'α + X'β(α))X]=0.


First, we construct the grid search interval for α  and then profile out the coefficient for each α in the interval on the exogenous variable using equation (5). Specifically,

    β̂(a) = min_b ∈ℬ1/N∑_i=1^N ρ_τ(Y_i - D_i'a - X_i'b).

By substituting these  estimates into the sample counterpart of the moment restriction (2), we obtain

    ĝ_N(a) = 1/N∑_i=1^N g(a,β̂(a),δ̂(a)),

where

    δ̂(a) = M(a)J^-1(a)

with

    M(a) = 1/Nh_N∑_i=1^N Z_i X_i'K_h_N(Y_i - D_i'a-X_i'β̂(a))


    J(a) = 1/Nh_N∑_i=1^N X_i X_i'K_h_N(Y_i - D_i'a-X_i'β̂(a))

where K_h_N is a kernel function with bandwidth h_N.
In the Monte Carlo simulations, we assume that we know the density function according to our data generation process. Thus, we can solve for the parameters by optimizing the criterion function of generalized method of moments (GMM) as follows: 

    α̂(τ) = min_a ∈𝒜 Nĝ_N(a)'Σ(a,a)^-1ĝ_N(a),

where

    Σ(a_1,a_2) = 1/N∑_i=1^N g(a_1,β̂(a_1))g(a_2,β̂(a_2))'

is a weighting matrix used in the GMM estimation.
Note that the estimator α̂ based on the inverse quantile regression (i.e., IVQR) of <cit.> is the first-order equivalent to the estimator defined by the GMM above.






 §.§ Estimation with High-dimensional Controls



We modify the procedure presented in Subsection 2.2  to deal with a dataset of high-dimensional control variables. 
To this end, we construct the grid search interval for α and profile the coefficients on exogenous variables using the L_1-norm penalized quantile regression estimator of <cit.>:

    β̂(a) = min_b ∈ℬ1/N∑_i=1^N ρ_τ(Y_i - D_i'a - X_i'b) + λ∑_j=1^dim(b)σ̂_j
    |b_j|,

where ρ(·) is the check function and σ̂^2_j = (1/n)∑_i=1^n x^2_ij.
The penalty level λ is chosen as follows.

    λ = 2 ·Λ(1-α | X),

where Λ(1-α | X):= (1-α)-quantile of Λ conditional on X.
The random variable 

    Λ = n sup_u∈𝒰max_1≤ j≤ dim(b)| 1/n∑_i=1^n [ x_ij(u - 𝕀{u_i ≤ u} )/σ̂_j √(u(1-u))] |,

where u_1, …, u_n are i.i.d. uniform (0, 1) random variables that are independently distributed from the controls x_1, …, x_n. The random variable Λ has a pivotal distribution conditional on 
X = [x_1, …, x_n]'. Therefore, we compute Λ(1-α | X) using simulation of Λ.
<cit.> show that the aforementioned choice for the penalty level λ leads to the optimal rates of convergence for the L_1-norm penalized quantile regression estimator.
Namely, 
the choice of the penalization parameter λ based on (13) is theoretically grounded and feasible.
In high-dimensions setting, K-fold cross-validation is very popular in practice.
However, computational cost is roughly proportional to K.
The recently derived non-asymptotic error bounds in <cit.> imply that the K-fold cross-validated  Lasso estimator has nearly optimal convergence rates. 
Although their theoretical guarantees do not directly apply to the L_1-norm penalized quantile regression estimator,  it still sheds some light on the use of cross-validation as an alternative to determine the penalty level λ in our analysis.[We have conducted Monte Carlo experiments indicating that the choice of λ based on (13) or 5-fold cross-validation leads to similar finite sample performances of our proposed procedure in terms of root-mean-square error, mean absolute error, and bias.  Simulation findings are tabulated in Section 3. When there are many binary control variables, the L_1-norm penalized quantile regression may suffer singularity issues in estimation. If this is the case, researchers can utilize the algorithm developed by <cit.> using the Huber loss function to approximate the quantile loss function.]

In addition, we estimate

    M(a) = 1/Nh_N∑_i=1^N Z_i X_i'K_h_N(Y_i - D_i'a-X_i'β̂(a))

and

    J(a) = 1/Nh_N∑_i=1^N X_i X_i'K_h_N(Y_i - D_i'a-X_i'β̂(a)).


We also perform dimension reduction on J because of the large dimension of X.
In particular, we implement the following regularization.

    δ̂_j(a) = min_δ1/2δ'Ĵ(a)δ - M̂_j(a)δ + ϑ||δ||_1.

The regularization above does a weighting Lasso for each instrument variable on control variables.
Consequently, the L_1 norm optimization obeys the Karush-Kuhn-Tucker condition

    ||δ̂_j(a)'Ĵ(a) - M̂_j(a)||_∞≤ϑ,   ∀ j.

More importantly, the aforementioned procedure is the double machine learning algorithm for the IVQR, which satisfies the Neyman orthogonality condition as follows.
Let us present the IVQR as a first-order-equivalent GMM estimator.
To this end, we define

    g(α, η) = (τ - 1(Y ≤ D'α + X'β)) ( Z - δ(α) X )

where η = [β(α)'  δ(α)' ]' are high-dimensional nuisance parameters in the DML setting discussed in Section 2.1 with true values 
η_0 = [β(α_0)'  δ(α_0)' ]'.
Therefore,

    𝔼[g(α_0, η_0)] 
       =𝔼[(τ - 1(Y ≤ D'α_0 + X'β_0)) ( Z - δ(α_0) X )]
       =𝔼[ 𝔼[ . τ - 1(Y ≤ D'α_0 + X'β_0) | X, Z ] ( Z - δ(α_0) X )]=0.

We then calculate

    .∂_η𝔼[g(α_0, η)] |_η=η_0
    = .[ .∂_β𝔼[g(α_0, η)] |_η=η_0; .∂_δ𝔼[g(α_0, η)] |_η=η_0 ]. .

Specifically,

    .∂_β𝔼[g(α_0, η)] |_η=η_0   = ∂_β𝔼[ 𝔼[ . τ - 1(Y ≤ D'α_0 + X'β_0) | X, Z ] ( Z - δ(α_0) X )]
       = ∂_β𝔼[ ( . τ - F(Y ≤ D'α_0 + X'β_0  | X, Z) ) ( Z - δ(α_0) X )]
       = 𝔼[ZX' f_ϵ(0|X, Z) ] - δ(α_0) 𝔼[XX' f_ϵ(0|X, Z) ]
       = M(α_0) -δ(α_0) J(α_0)
       = M(α_0) -M(α_0)J^-1(α_0) J(α_0)=0.


    .∂_δ𝔼[g(α_0, η)] |_η=η_0   = ∂_δ𝔼[(τ - 1(Y ≤ D'α_0 + X'β_0)) ( Z - δ(α_0) X )]
       = - 𝔼[(τ - 1(Y ≤ D'α_0 + X'β_0)) X ] =0.

We thus verify that .∂_η𝔼[g(α_0, η)] |_η=η_0=0,
which indicates the Neyman orthogonality condition holds.




After implementing the DML outlined above,
we solve for the low-dimensional causal parameter α by optimizing the GMM defined as follows.
The sample counterpart of the moment condition

    ĝ_N(a) = 1/N∑_i=1^N (τ - 1(Y_i - D_i'a -X_i'β̂(a) ≤ 0))Ψ(a,δ̂(a)).

Accordingly,

    α̂ = min_a ∈𝒜 Nĝ_N(a)'Σ(a,a)^-1ĝ_N(a).



<cit.> show that the key condition enabling us to perform valid inference on α_0 is the adaptivity condition:
√(N)(ĝ(α_0, η̂) -  ĝ(α_0, η_0) )
ℙ_N⟶ 0.
In particular, each element ĝ_j of ĝ=(ĝ_j )_j=1^k
can be expanded as
√(N)(ĝ_j(α_0, η̂) -  ĝ_j(α_0, η_0) )
=T_1, j + T_2, j +T_3, j, which are formally defined on page 663 in their paper.
The term T_1, j vanishes precisely because of orthogonality, that is, T_1, j=0.
However, the terms T_2, j and T_3, j do not vanish.
The T_2, j and T_3, j vanish when cross-fitting and sample splitting are implemented. 
These two terms are also asymptotically negligible when we impose a further structure on the problem:
such as using a sparsity-based machine learner (e.g., L_1-norm penalized quantile regression)
under approximate sparsity conditions. 
In our procedure, equations (12) and (17) are sparsity-based machine learners.
Therefore, we use no cross-fitting in the DML-IVQR algorithm.

Theoretically speaking, based on equation (19), the approach can be applied to machine learners other than the Lasso. The chief difficulty in implementing an estimation based on Equation (19) is that 
the function being minimized is both non-smooth and non-convex,
and any machine learners are used to dealing with a functional response variable in this context,
cf. <cit.>.
In addition, the corresponding DML with non-linear equations is difficult. 
Therefore, our practical strategy is to implement the DML-IVQR procedure 
described in equations (12) – (17) and (21) – (22), which is equivalent to the Neyman orthogonality condition defined in (19) and (20).









 §.§ Weak-Identification Robust Inference



Under the regularity conditions listed in <cit.>, 
asymptotic normality for the GMM estimator with a non-smooth objective function is guaranteed.
We have

    √(n)ĝ_N(a) d⟶  N(0,Σ(a,a)).

Consequently, it leads to

    Nĝ_N(a)'Σ(a,a)^-1ĝ_N(a) d⟶χ^2_dim(Z).

We define

    W_N≡ Nĝ_N(a)'Σ(a,a)^-1ĝ_N(a).

It  follows that a valid
(1-p) percent confidence region for the true parameter, α_0, can be constructed as the set

    CR := {α∈𝒜:W_N (α) ≤ c_1-p},

where c_1-p is the critical point such that

    P[χ^2_dim(Z) > c_1-p]=p,

and 𝒜 can be numerically approximated by the grid {α_j,j=1,...,J}.






§ MONTE CARLO EXPERIMENTS



We evaluate the finite-sample performance, in terms of 
mean bias (BIAS), mean absolute error (MAE) and root-mean-square error (RMSE)
of the DML-IVQR through 1000 simulations.
The following data generating process is modified from that considered in <cit.>.



    [ u_i; ϵ_i;     ]∼ N(0,[   1 0.3; 0.3   1;     ])


    [ x_ji; z_1i; z_2i; v_1i; v_2i;      ]∼ N (0,I)


    Z_1i = z_1i + x_2i+ x_3i+  x_4i+ v_1i


    Z_2i = z_2i + x_7i+ x_8i + x_9i+ x_10i+ v_2i


    D_i = Φ(z_1i+z_2i +ϵ_i)


    X_ji = Φ(x_ji)


    Y_i =1+ D_i + 5X_1i + 5X_2i + 5X_3i + 5X_4i + 5X_5i + 5X_6i + 5X_7i + D_i× u_i,

where Φ(·) is the cumulative distribution function of a Normal random variable;
i=1, 2, …, n; j= 1, 2, …, p; p is the dimension of controls X, and p=100.
There are ten relevant controls:  X_1i, …, X_10i. 
The instrumental variable is Z.
The target variable is D.
Consequently,

    α(τ)=1+F_ϵ^-1(τ),

where τ is the quantile index and F_ϵ (·) is the cumulative distribution function of the random variable ϵ. 
Therefore, the median treatment effect α(0.5)=1.




 §.§ Residualizing Z on X


We focus on comparing the BIAS, MAE and RMSE resulting from different procedures under the exact specification (10 control variables).  res-GMM represents residualizing Z on X.
GMM stands for doing no residualizing Z on X.
Table 1 shows that residualizing Z on X leads to an efficiency gain across quantiles especially when the sample size is moderate. 










 §.§ IVQR with High-dimensional Controls


We now evaluate the finite-sample performance of the IVQR with high-dimensional controls.
The data generating process involves 100 control variables with an approximate sparsity structure.
In particular, the exact model (true model) depends only on 10 relevant control variables out of the 100 controls. 

Let's fix the name of different estimators first.
The full-GMM uses 100 control variables without regularization. 
The oracle-GMM knows the identity of the true controls and then uses the ten relevant variables.
The DML-IVQR is our proposed estimator.
Table 2 shows that the RMSE  stemmed from the DML-IVQR are close to those from the oracle estimator. The numbers in parentheses are ratios of RMSE or MAE of any estimator to those of the oracle-GMM. 
The BIAS and MAE indeed signify that the DML-IVQR achieves a lower bias in the simulation study.  
In addition,  Figure 1 plots the distributions of the IVQR estimator with and without double machine learning.
The DML-IVQR stands for the double machine learning for the IVQR with high-dimensional controls.
Histograms signify that the DML-IVQR estimator is more efficient and less biased than IVQR using many control variables. Because a weak-identification robust inference results naturally from the IVQR, we construct the robust confidence regions for the full-GMM, oracle-GMM and the DML-IVQR estimators. 
In Figures 2-4, the vertical axis displays the value of the test statistic W_N(α) which is defined in Section 2.4. The horizontal line in gray is the 95% critical value from χ^2_dim(Z).
<cit.> robust confidence region is all values of α such that
the W_N(α) lies below the horizontal line. 
The robust inferential procedure is still valid when identification is weak or fails partially or completely.
Thus Figures 2-4 show that, across quantiles, the robust confidence region based on the DML-IVQR is relatively sharp compared to those of the full-GMM.
In addition, the confidence regions based on the DML-IVQR are remarkably close to those obtained by the oracle estimator.

As to the choice of penalty parameter, 
researchers can chose λ based on Equation (13) proposed by <cit.>  
or based on the K-fold cross-validation.
Both methods of choosing λ  lead to similar finite sample performances of DML-IVQR in terms of the RMSE, MAE and BIAS.  Simulation findings are summarized in Table 3.
  






 


 
 
 
 















 















§ AN EMPIRICAL STUDY: QUANTILE TREATMENT EFFECTS OF 401(K) PARTICIPATION ON ACCUMULATED WEALTH


In this section, 
we reinvestigate an empirical study on quantile treatment effects: 
the effect of 401(k) participation on wealth, cf. <cit.>.
Not only does this conduct data-driven robustness checks on the econometric results, 
but the DML-IVQR sheds light on the treatment effect heterogeneity among the control variables. 
This complements the existing empirical findings. 
In addition, we compare our empirical results with those from <cit.> that conduct
the IVQR estimation by using generalized random forest approach,
which is an alternative in causal machine learning literature. 


Examining the effects of 401(k) plans on accumulated wealth is an issue of long-standing empirical interest. 
For example, based on the identification of selection on observables, 
<cit.> and <cit.> suggest that 
the income nonlinear effect exists in the 401(k) study. 
Nonlinear effects from other control variables are identified as well. 


Based on DML-IVQR, we reinvestigate the impact of the 401(k) participation on accumulated wealth.
Total wealth (TW) or net financial assets (NFTA) is the outcome variable Y.
The treatment variable D is a binary variable that stands for participation in the 401(k) plan.
Instrument Z is an indicator of eligibility to enroll in the 401(k) plan.
The vector of covariates X consists of income, age, family size, marriage, an IRA individual retirement account, a defined benefit status indicator, a home ownership indicator and the different education-year indicator variables. 
The data consists of 9915 observations.


Following the regression specification set up in <cit.>,
Table 3 presents quantile treatment effects obtained from different estimation procedures which 
have been defined in the previous sections including
IVQR, res-GMM and GMM. 
The resulting estimates are similar.
As to the high-dimensional analysis, we create 119 technical control variables 
including those constructed by polynomial bases, interaction terms, and cubic splines (thresholds). 
To ensure each basis has equal length, 
we utilize the minimax normalization for all technical control variables.
Consequently, we use the plug-in method to determine the penalty value when performing the Lasso under the moment condition, and tune the penalty in the quantile L1-norm objective function based on
the Huber approximation by 5-fold cross-validation.
The DML-IVQR also implements feature normalization of the outcome variable for  computational efficiency. 
To make the estimated treatment effects across different estimation procedures roughly comparable, 
Table 5 shows that the effect obtained through the DML-IVQR is multiplied by the standard deviation of the outcome variable.
Weak identification/instrument robust inference on quantile treatment effects are depicted in 
Figures 5 and 6.
However, the robust confidence interval widens as the sample size decreases at the upper quantiles.
Estimated quantile treatment effects are significantly different from zero.
We can use the result from the DML-IVQR as a data-driven robustness check on 
those summarized in the Table 4.

Tables 6 and 7 present the selected important variables across different quantiles.
The approximate sparsity is asymmetric across the conditional distribution in the sense that
the number of selected variables decreases as the quantile index τ increases,
although it hinges on a relatively small number of observations at the upper quantiles.
In this particular example, 
τ captures the rank variable that governs the unobservable heterogeneity: savings propensity.
Small values of τ  represent participants with low savings propensity.
Our empirical results thus signify that the 401(k) participants with low savings propensity are more associated with the nonlinear income effect than those with high savings propensity, which complements the results concluded in previous studies
<cit.> and <cit.>.
The nonlinear income effects, across quantiles ranging from (0, 0.5], are picked up by the selected variables, such as
max(0, inc-0.2), max(0, inc^2-0.2),max(0, inc^3-0.2) and etc.
Technical variables in terms of age, education, family size, and income are more frequently selected in Tables 6 and 7.
In addition, these four variables are also identified as important variables in the 
context of the generalized random forests, cf.  <cit.>.










































§ CONCLUSION




In this study, we investigate the performance of a debiased/double machine learning algorithm within the framework of high-dimensional IVQR. The simulation results indicate that our procedure performs more efficiently than those based on conventional estimators with many controls. Furthermore, we evaluate the corresponding weak identification robust confidence interval of the low-dimensional causal parameter. 
Given many technical controls, we reinvestigate  thequantile treatment effects of the 401(k) participation on accumulated wealth and then highlight the non-linear income effects across the savings propensity.









 










Three authors contributed equally to the paper.


This research was funded by the JSPS KAKENHI (Grant no. JP20K01593), 
the personal research fund from Tokyo International University, and it was
financially supported by the Center for Research in Econometric Theory and Applications (Grant no. 107L900203 ) from The Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education (MOE) in Taiwan.



We are grateful to Tsung-Chih Lai and Hsin-Yi Lin for discussions and comments. This paper has benefited from presentations at Ryukoku University and the 2nd International Conference on Econometrics and Statistics (EcoSta 2018). The usual disclaimer applies.


The authors declare no conflict of interest. 



The following abbreviations are used in this manuscript:



DML     Double machine learning

GMM     Generalized method of moments

GRF     Generalized random forests

IVQR     Instrumental variable quantile regression

Lasso     Least absolute shrinkage and selection operator




no 







References













999



[Athey(2017)]ref-Athey2017
Athey, S. 2017. 
Beyond prediction: Using big data for policy problem.
Science, 355, 483–485.



[Athey and Imbens(2019)]ref-AtheyImbens2019
Athey, S. and G. Imbens. 2019.
Machine learning method that economists should know about.
Annual Review of Economics, 11, 685–725.



[Athey, Tibshirani and Wager(2019)]ref-Athey2019
Athey, S., Tibshirani, J., and S. Wager. 2019.
Generalized random forests.
Annals of Statistics, 47(2), 1148–1178.



[Belloni and Chernozhukov(2011)]ref-BelloniChernozhukov2011
Belloni, A. and V. Chernozhukov. 2011.
l_1-penalized quantile regression in high-dimensional sparse models.
Annals of Statistics, 39(1), 82–130.


[Belloni, Chernozhukov,  Fernández-Val and Hansen(2017)]ref-BCFH2017
Belloni, A., Chernozhukov, C., Fernández-Val, I. and C. Hansen. 2017.
Program evaluation and causal inference with high-dimensional data.
Econometrica, 85(1), 233–298.



[Chen and Hsiang(2019)]ref-ChenHsiang2019
Chen, J.-E. and C.-W. Hsiang. 2019. 
Causal random forests model using instrumental variable quantile regression.
Econometrics, 7, 1-22.


[Chen and Lee(2018)]ref-ChenLee2018
Chen, L.-Y. and S. Lee.  2018.
Exact computation of GMM estimators for instrumental variable quantile regression models.
Journal of Applied Econometrics, 33(4), 553–567.



[Chernozhukov et al.(2018)]ref-Chernozhukov2018
Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and J. Robins.  2018.
Double/debiased machine learning for treatment and structural parameters. 
The Econometrics Journal, 21(1), C1–C68.


[Chernozhukov and Hansen(2004)]ref-ChernozhukovHansen2004
Chernozhukov, V. and C. Hansen. 2004.
The effects of 401(k) participation on the wealth distribution: An instrumental quantile regression analysis. 
Review of Economics and Statistics, 86(3), 735–751.



[Chernozhukov and Hansen(2005)]ref-ChernozhukovHansen2005
Chernozhukov, V. and C. Hansen. 2005.
An IV model of quantile treatment effects. 
Econometrica, 73(1), 245–261.


[Chernozhukov and Hansen(2008)]ref-ChernozhukovHansen2008
Chernozhukov, V. and C. Hansen. 2008.
Instrumental variable quantile regression: A robust inference approach. 
Journal of Econometrics, 142, 379–398.


[Chernozhukov and Hansen(2013)]ref-ChernozhukovHansen2013
Chernozhukov, V. and C. Hansen. 2013.
NBER 2013 Summer Institute: Econometric methods for high-dimensional data. 
http://www.nber.org/econometrics_minicourse_2013/.


[Chernozhukov, Hansen, and Spindler(2015)]ref-ChernozhukovHansenSpindler2015
Chernozhukov, V., Hansen, C. and M. Spindler. 2015.
Valid post-selection and post-regularization inference: An elementary, general approach. 
Annual Review of Economics, 7, 649–688.



[Chernozhukov, Hansen, and Wüthrich(2018)]ref-Chernozhukovetal2018
Chernozhukov, V., Hansen C., and K. Wüthrich. 2018.   
Instrumental variable quantile regression. 
Handbook of Quantile Regression, Chapman & Hall/CRC.


[Chetverikov, Liao and Chernozhukov(2021)]ref-Chetverikovetal2021
Chetverikov, D., Liao, Z. and V. Chernozhukov. 2021.
On cross-validated lasso in high dimensions.
Annals of Statistics, forthcoming




[Chiou, Chen, and Chen(2018)]ref-ChiouChenChen2018
Chiou, Y.-Y., Chen, M.-Y., and J.-E. Chen. 2018.
Nonparametric regression with multiple thresholds: estimation and inference.
Journal of Econometrics, 206(2), 472–514.






[Davis and Heller(2017)]ref-Davis2017
Davis, J. and S.B. Heller. 2017.
Using causal forests to predict treatment heterogeneity: An application to summer jobs. 
American Economic Review, 107(5), 546–550.



[Gilchrist and Sands(2016)]ref-Gilchrist2016
Gilchrist, D.S. and E.G. Sands. 2016.
Something to talk about: Social spillovers in movie consumption. 
Journal of Political Economy, 124(5), 1339–1382.




[Yi and Huang(2017)]ref-YiHuang2017
Yi, C. and J. Huang. 2017.
Semismooth Newton coordinate descent algorithm for elastic-net penalized Huber loss regression and quantile regression.
Journal of Computational and Graphical Statistics, 26(3) 547–557.













