Year: 2021

Abstract:
We consider learning causal relationships under conditional moment restrictions. Unlike causal inference under unconditional moment restrictions, conditional moment restrictions pose serious challenges for causal inference, especially in highdimensional settings. To address this issue, we propose a method that transforms conditional moment restrictions to unconditional moment restrictions through importance weighting, using a conditional density ratio estimator. Using this transformation, we successfully estimate nonparametric functions defined under conditional moment restrictions. Our proposed framework is general and can be applied to a wide range of methods, including neural networks. We analyze the estimation error, providing theoretical support for our proposed method. In experiments, we confirm the soundness of our proposed method.

heading: 1 INTRODUCTION
text: Consider learning the causal relationship between airline ticket prices and demand. As one might expect, prices and demand rise and fall through the seasons, being affected by other events like vacation periods, which are called confounders and may or may not be observable. Due to confounders, naively inferring from this pattern that higher (lower) prices increase (decrease) demand would be incorrect, and potentially detrimental. Thus, controlling for confounding effects is essential. This issue frequently arises in practice, especially when learning causal (structural) relationships is essential to answer counterfactual questions regarding policy intervention and outcome.
One approach to deal with confounding effects (like in the airline example above) is the instrumental variable (IV) approach (Wooldridge, 2002; Greene, 2003). In the IV approach, the conditional moment restriction is defined as such that the causal model satisfies the restriction of zero expected value given IVs, thus conditioning out the confounding effect. The simplest representation of this idea is the two-step least squares (2SLS) for linear models (Wooldridge, 2002; Greene, 2003). However, given the complex nature of the causal effect and confounding effect (and their relation), assuming a linear relation can be too strong. Thus, in this paper, we focus on nonparametric IV (NPIV) regressions, allowing for much more flexible estimation (Newey & Powell, 2003).
NPIV can be viewed as an instance of a more general framework of causal inference under conditional moment restrictions. In this light, the machinery for inference under conditional moment restrictions also applies to NPIV, as well as its shortcomings. One major issue with using the conditional moment restrictions for causal inference is that one must approximate the conditional expectation, which is often difficult to do (see, e.g., Newey (1993); Donald et al. (2003) for parametric and Newey & Powell (2003); Ai & Chen (2003) for nonparametric IVs using sieves). For instance, Lewbel (2007) and Otsu (2007) estimate the conditional expectation by local kernel density estimation. However, local kernel density estimation suffers under high dimensionality. For this problem, recent methods suggest the use of machine learning methods, such as neural networks (Hartford et al., 2017).
In this paper, we propose transforming conditional moment restrictions into unconditional moment restrictions by importance weighting using the conditional density ratio, which is defined as the ratio of the conditional probability density, conditioned on the IVs, to the unconditional probability density. We show that the unconditional expectation of a random variable weighted by the conditional density ratio is equal to the conditional expectation. Further, we show that it is possible to estimate the
conditional density ratio with the least-squares method with a neural network. Once the conditional density ratio is estimated, the usual method of moments, such as GMM, can be used straightforwardly.
The contribution of this paper is as follows: (i) we propose a novel approach to convert conditional moment restrictions to unconditional moment restrictions by importance weighting; (ii) using our proposed transformation, we develop methods for NPIV; (iii) we analyze the estimation error

heading: 2 SETUP AND NOTATION
text: Among various problems of learning causal relationships from conditional moment restrictions, we focus on the NPIV regression for ease of discussion. Note that our proposed method can be applied to more general settings, similar to Ai & Chen (2003).
Suppose that the observations {(Yi, Xi, Zi)}ni=1 are i.i.d., where Yi ∈ Y ⊆ R is an observable scalar random variable, Xi ∈ X ⊆ RdX is a dX dimensional explanatory variable, Zi ∈ Z ⊆ RdZ is a dZ dimensional random variable, called an IV, and X and Z are compact with nonempty interior. We also assume that the probability densities of (Yi, Xi, Zi), (Yi, Xi), Zi exist and denote them by p(y, x, z), p(y, x), and p(z), respectively. Let us define the causal relationships between Yi and Xi as
Yi = f ∗(Xi) + εi,
where f∗ : X → Y is a structural function, εi is the sub-Gaussian error term with mean zero. To learn f∗, suppose the IV Zi satisfies the following conditional moment restrictions:
E [εi|Zi] = 0 ∀i ∈ {1, 2, . . . , n}. (1) Then, we also assume that under the conditional moment restriction, we can uniquely identify f∗. Our goal is to learn f∗ from the conditional moment restrictions in (1). If Zi = Xi, this problem boils down to the estimation of the conditional expectation (regression function) E[Yi|Xi]. However, when E[εi|Xi] 6= 0, E[Yi|Xi] is not equivalent to f∗(Xi); that is, typical regression analysis, such as least squares, may not return the correct estimate of f∗(Xi).

heading: 3 PRELIMINARIES AND LITERATURE REVIEW
text: In this section, we briefly review causal inference under moment restrictions.

heading: 3.1 IV METHOD FOR LINEAR STRUCTURAL FUNCTIONS
text: One of the basic cases of using the IV is when f∗ is a linear model X>i θ ∗ with dX dimensional vector θ∗ and the error term εi is correlated with the explanatory variable Xi. In this case, the parameter θ∗ of the linear model can be estimated if there are IVs of dimension dX or more that satisfy the unconditional moment restrictions, E[Ziεi] = 0dZ , where 0d is a d dimensional zero vector. In the just-identified case (dX = dZ), we can estimate θ∗ as θ̂IV = ( 1 n ∑n i=1 ZiX > i )−1 1 n ∑n i=1 ZiYi. In the over-identified case (dX ≤ dZ), we can estimate it by the 2SLS. In the 2SLS, we first regress Xi by Zi; then using a dX dimensional function ĝ(Zi) ob-
tained from the first stage regression, we estimate θ∗ as θ̂2SLS = ( 1 n ∑n i=1 X̂iX̂ > i )−1 1 n ∑n i=1 X̂iYi,
where X̂d,i = Z>i ( 1 n ∑n j=1 ZjZ > j )−1 1 n ∑n j=1 ZjXd,j and X̂i = (X̂1,i, . . . , X̂dX ,i) >.
More generally, we can formulate the estimation of the linear structural function by unconditional moment restrictions defined using the IV; that is, E[m(θ∗;Yi, Xi, Zi)] = 0dm , where θ∗ ∈ Θ is a parameter representing the causal effect, Θ is the parameter space, and m : Θ×R×X×Z → Rdm is a dm dimensional moment function. Then, the GMM estimator is defined as θ̂GMM = arg minθ∈Θ ( 1 n ∑n i=1m(θ;Yi, Xi, Zi) )> Wdm ( 1 n ∑n i=1m(θ;Yi, Xi, Zi) ) , where Wdm is a dm × dm weight matrix. If Wdm is chosen as W ∗dm ∝ E[m(θ ∗;Yi, Xi, Zi)
>m(θ∗;Yi, Xi, Zi)], the estimator θ̂GMM is efficient under the posited moment conditions. Note that the GMM includes the 2SLS as a special case where εi has the same variance among i ∈ {1, . . . , n} and the zero covariance. Three methods have been proposed to obtain W ∗dm ; two-step GMM, iterative GMM (Hayashi, 2000), and continuous updating GMM (CU estimator; CUE, Hansen et al., 1996). These estimators have the
same asymptotic distribution, but different non-asymptotic properties. In particular, CUE is known as a special case of the generalized empirical likelihood (GEL) estimator (Owen, 1988; Smith, 1997), which plays an important role in estimation with many moment restrictions (Newey & Smith, 2004).

heading: 3.2 NPIV AND LEARNING FROM CONDITIONAL MOMENT RESTRICTIONS
text: In linear structural functions, zero covariance between the IV and error term suffices for identification. However, in nonparametric structural functions, we require a stronger restriction: the error term has conditional expectation zero given the IVs. Then, we can characterize the solution of NPIV as that of an integral equation K(z) = E[Yi|Zi = z] = E[f∗(Xi)|z] = ∫ f∗(x)dF (x|z) with F denoting the conditional c.d.f of x given z. The identification also results in the uniqueness of the solution.
Several estimators have been proposed. Newey & Powell (2003) proposes a nonparametric analogue of the 2SLS for linear structural functions. They first define the linear-in-parameter model as a model to approximate the nonparametric structural function f∗, where they use a linear approximation with basis expansion, such as sieve (series) regression. Let us denote the approximation as f∗(Xi) ≈ ψ(Xi)
>θ∗, where ψ : X → Rdψ , dψ ≤ n is a vector-valued function consisting of outputs of basis functions. Then, they conduct the 2SLS as follows: (i) they define a linear-in-parameter model as an approximation to the nonparametric model of Zi, and estimate E[ψ(Xi)|Zi] using that model; (ii) they regress Yi by E[ψ(Xi)|Zi] to estimate θ∗. In contrast, Ai & Chen (2003) proposes a nonparametric analogue of the GMM. Ai & Chen (2003) also approximate the nonparametric structural function f∗ by a linear-in-parameter model. However, unlike Newey & Powell (2003), Ai & Chen (2003) estimates E[(Yi − f(Xi))|Zi], instead of estimating E[ψ(Xi)|Zi], by another linear-in-parameter model. Using the estimator of E[(Yi−f(Xi))|Zi], Ai & Chen (2003) transforms conditional moment restrictions into unconditional ones and apply the GMM to estimate f∗. In addition to these two typical methods, Darolles et al. (2011), Carrasco et al. (2007) propose their own methods for NPIV.
This is how the NPIV problem is typically cast into the framework of conditional moment restrictions, where a parameter representing the causal relationship satisfies E[m(θ∗;Yi, Xi, Zi)|Zi] = 0dm . As with NPIV, we often define the moment function as a function of the nonparametric function f instead of the parameter θ. We can estimate f∗ defined under conditional moment restrictions using variants of GMM or empirical likelihood, e.g., Ai & Chen (2003); Domínguez & Lobato (2004); Otsu (2007); Lewbel (2007), by transforming the conditional moment restrictions to unconditional ones.

heading: 3.3 APPROACHES IN THE MACHINE LEARNING LITERATURE
text: 2SLS with more flexible models. Hartford et al. (2017) extends the two-stage approach by employing a neural network density estimator, and Singh et al. (2019) does it by conditional mean embedding in RKHS. For example, in Hartford et al. (2017), they first approximate K(z) by estimating F (x|z) with neural networks. Then, for B samples {X̃j}Bj=1 generated from the estimator F̂ (x|z), they train
neural networks by minimizing the empirical risk 1n ∑n i=1 ( Yi − 1B ∑ X̃j∼F̂ (Xj |Zi) f ( X̃j ))2 .
In contrast, Xu et al. (2021a) first approximates f(Xi) by neural networks and regard the last layer as ψ(Xi) in the 2SLS of NPIV. Then, they predict ψ(Xi) by another model ϕ(Zi). Because ψ(Xi) is transitive in the learning process, they alternatively train the models ϕ(Zi) and f(Xi).
Minimax approach. There has also been a recent surge in interest in minimax approaches that reformulate conditional moment restrictions as a minimax optimization problem (Bennett et al., 2019; Bennett & Kallus, 2020; Muandet et al., 2020; Chernozhukov et al., 2020; Liao et al., 2020). For this approach, Dikkala et al. (2020) shows a strong theoretical guarantee on the estimation error bound.

heading: 4 PROPOSED METHOD: NPIV BY IMPORTANCE WEIGHTING
text: Suppose that p(y, x) > 0 for all (y, x) ∈ Y × X . Note that we have already assumed the existence of p(y, x|z) for all (y, x, z) ∈ Y × X × Z by assuming the existence of the conditional moment restrictions. Define the conditional density ratio function r : Y × X × Z → (0, C) as
r∗(y, x|z) = p(y, x|z) p(y, x) = p(y, x, z) p(y, x)p(z) ,
for 0 < C <∞. We assume the existence of the conditional density ratio function. Assumption 1. For all (y, x, z) ∈ Y × X × Z , the conditional density ratio r∗(y, x|z) exists.
Then, we transform conditional moment restrictions into unconditional moment restrictions as E[(Yi − f(Xi))|z] = ∫ ( y − f(x) )p(y, x|z) p(y, x) p(y, x)dydx = E [( Yi − f(Xi) ) r∗(Yi, Xi|z) ] .
Thus, if we know the conditional density ratio function r∗(y, x|z), we can approximate E[(Yi − f(Xi))|z] for z ∈ Z by the following sample average:
1
n n∑ i=1 ( Yi − f(Xi) ) r∗(Yi, Xi|z).
Based on this property, we propose the following two-stage method: (i) estimate the conditional density ratio function r∗; (ii) approximate conditional moment restrictions by the sample average of( Yi − f(Xi) ) using the estimate of the conditional density ratio function r∗. Since we do not know the true value of the conditional density ratio r∗ and cannot calculate the expected value, we consider replacing r∗ with an estimator r̂ and approximating the expected value with the sample mean.

heading: 4.1 CONDITIONAL DENSITY RATIO ESTIMATION
text: First, we consider estimating r∗(y, x|z). While it is possible to estimate the probability density functions of the numerator and the denominator, individually, following Vapnik’s principle (Vapnik, 1998), we should avoid solving more difficult intermediate problems. Sugiyama et al. (2012) summarizes various methods to estimate the density ratio directly. Inspired by least-squares importance fitting (LSIF) of Kanamori et al. (2009), we estimate the conditional density ratio by minimizing
r̃ = arg min r∈R̃
1 2 EY,X
[ EZ [( r∗(Yi, Xi|Zj)− r(Yi, Xi|Zj) )2]] ,
where R̃ denotes a set of measurable functions and for a function g : Y × X × Z → R, EY,X [EZ [g(Yi, Xi, Zj)]] denotes ∫ g(y, x, z)p(y, x)p(z)dydxdz. It is easy to confirm that r̃ = r∗ by taking the first derivative of the risk. Because this risk includes the unknown r∗, it may seem intractable objective function. However, we can obtain the risk that does not include r∗ as
arg min r∈R̃
1
2
( EY,X [ EZ [ r∗2(Yi, Xi|Zj)− 2r∗(Yi, Xi|Zj)r(Yi, Xi|Zj) + r2(Yi, Xi|Zj) ]]) = arg min
r∈R̃
{ −EY,X [EZ [r∗(Yi, Xi|Zj)r(Yi, Xi|Zj)]] + 1
2 EY,X
[ EZ [ r2(Yi, Xi|Zj) ]]} = arg min
r∈R̃
{ −EY,X,Z [r(Yi, Xi|Zi)] + 1
2 EY,X
[ EZ [ r2(Yi, Xi|Zj) ]]} .
Here, we used EY,X [EZ [r∗(Yi, Xi|Zj)r(Yi, Xi|Zj)]] = EY,X [ EZ [ p(Yi,Xi,Zj) p(Yi,Xi)p(Zj) r(Yi, Xi|Zj) ]] = EY,X,Z [r(Yi, Xi|Zi)]. For some hypothesis class R, by approximating the risk with its sample approximation, we estimate the conditional density ratio as
r̂ = arg min r∈R { − 1 n n∑ i=1 r(Yi, Xi|Zi) + 1 2 1 n n∑ j=1 1 n n∑ i=1 r2(Yi, Xi|Zj) } .
For the hypothesis classR, we can use various models, such as linear-in-parameter models and neural networks. Suzuki et al. (2008) also proposes a similar formulation based on maximum likelihood estimation with constraints, which is not easy to solve with neural networks.

heading: 4.2 NPIV REGRESSION BY IMPORTANCE WEIGHTING
text: If we know the conditional density ratio function r∗, we can obtain unconditional moment restrictions EY,X [( ρ(f ;Yi, Xi, Z1, r ∗) ρ(f ;Yi, Xi, Z2, r ∗) · · · ρ(f ;Yi, Xi, Zn, r∗) )>] = 0n,
where ρ(f ;Yi, Xi, z, r∗) = ( Yi − f(Xi) ) r∗(Yi, Xi|z). By replacing r∗ and its expectation with its estimator and the sample average, we obtain the sample vector moment restrictions, 1 n ∑n i=1 ( ρ(f ;Yi, Xi, Z1, r̂) ρ(f ;Yi, Xi, Z2, r̂) · · · ρ(f ;Yi, Xi, Zn, r̂) )> . Once we obtain the sample average, we can apply various methods for learning f∗ from the unconditional moment restriction. For instance, for a hypothesis class F , we estimate f∗ by minimizing
Rn(f, r̂) = 1
n n∑ j=1
( 1
n n∑ i=1 (Yi − f(Xi))r̂(Yi, Xi|Zj)
)2 . (2)
This objective function is closely related with the projected mean squared error (MSE) introduced in Dikkala et al. (2020), defined as
L(f, r∗) := EZ [ (EY,X [(f∗(Xi)− f(Xi))r∗(Yi, Xi|Zj)])2 ] , (3)
The objective function (2) is a special case of the GMM and GEL. In the GMM, we estimate f∗ by minimizing 1n ∑n j=1 ( 1 n ∑n i=1(Yi − f(Xi))r̂(Yi, Xi|Zj) )> wj ( 1 n ∑n i=1(Yi − f(Xi))r̂(Yi, Xi|Zj) ) , where wj > 0 is a weight (Ai & Chen, 2003). Thus, our framework of transforming conditional moment restrictions to unconditional moment restrictions using importance weighting allows us to conduct causal inference using NPIV with conditional moment restrictions using a variety of models, such as linear-in-parameter models with basis functions and neural networks.
Linear-in-parameter models. As an example, we introduce a linear-in-parameter model with some basis functions. Here, let us consider using the Gaussian kernel as the basis function. For x ∈ X , let ϕ(x;σ2) = ( K(x,X1;σ 2), . . . ,K(x,Xn;σ 2) )> , where K(x,Xu;σ2) = exp ( −‖x−Xu‖ 2 2
2σ2
) is
the Gaussian kernel with a hyperparameter σ2 > 0 and ‖ · ‖2 is the L2 norm. Then, we define a linear-in-parameter model as f(x;σ2) = β>ϕ(x;σ2) + β0, where β ∈ Rn, and β0 ∈ R.
Neural networks. We can also use neural networks for approximating f∗. In this case, we need to carefully determine the network structures because we cannot uniquely determine the solution owing to the overparameterization. In existing studies, it is assumed that the network modes are well defined for the complexity of the nonparametric function f∗ and the sample size n.
Advantages of our proposed method. Our method has the following three advantages: (i) our method is applicable to general causal inference problems, and is not limited to the NPIV problem formulated in Section 2, as well as the Ai & Chen (2003), (ii) our method can deal with high dimensional variables, owing to the machine learning technique. This is in contrast to related studies Ai & Chen (2003) and Otsu (2007), which use a sieve and Nadaraya-Watson estimator, respectively, and hence do not work in high dimension settings, and (iii) our method is computationally efficient by the use of the importance weight approach. A similar study Hartford et al. (2017) requires additional sampling of Yi from an estimated conditional density p(y|z, y). Another related study Dikkala et al. (2020) requires a difficult algorithm to solve its minimax optimization problem. Our importance weight approach avoids such computational burden.

heading: 4.3 LEARNING WITH OVERPARAMETRIZED MODELS
text: In learning from moment restrictions, the structural function f∗ is determined by the equations that the expected value of random variables satisfy. When approximating the structure function f∗ using a model with more parameters than the sample size, the solution cannot be uniquely determined. When minimizing the prediction error directly, such overparameterization may not pose a major problem, and, in fact, may be necessary to improve the generalization error, as a recent finding suggests in the case of linear regression (Bartlett et al., 2020). However, in NPIV, the objective function is the set of conditional moment restrictions. Unlike in direct prediction error minimization, a model trained to minimize empirically approximated conditional moment restrictions does not necessarily minimize the MSE for f∗ (E[(f(Xi) − f∗(Xi))2]). In fact, we empirically confirm that methods using neural network sometimes do not work well partly because of their overfitting to empirical moment minimization, not to the empirical MSE for f∗. Despite these potential problems, there is
a strong motivation to use neural networks, owing to the reported superiority in some applications, such as computer vision (Xu et al., 2021a;b; Yuan et al., 2021), natural language processing (Ash et al., 2019; Chen et al., 2020) tasks. For this reason, we introduce a heuristic to avoid this problem.
As we explained, the parameters are not uniquely determined by conditional moment restrictions alone. Therefore, from the set of parameters satisfying conditional moment restrictions, we need to select a set of parameters that works well in prediction. We consider training a model that has the minimum MSE with Yi while satisfying conditional moment restrictions as follows:
min f∈F
E [ (Yi − f(Xi))2 ] s.t. E [(Yi − f(Xi))|Zj ] = 0 ∀j ∈ {1, 2, . . . , n}.
Since it is difficult to solve constrained optimization with neural networks, we propose to solve the penalized optimization. In the case of linear combination with penalties, we train the model by
min f∈F
1
n n∑ i=1 (Yi − f(Xi))2 + η n∑ j=1
( 1
n n∑ i=1 (Yi − f(Xi))r̂(Yi, Xi|Zj)
)2 , (4)
where η ≥ 0 is a regularization coefficient. The motivation of this heuristic is to select a function f(x) that is the closest to E[Yi|Xi = x], among multiple functions satisfying the conditional moment restriction. This heuristic works well when f∗(x) takes a near value of E[Yi|Xi = x] while f∗(x) 6= E[Yi|Xi = x].

heading: 5 ESTIMATION ERROR ANALYSIS
text: We show the estimation error of the conditional density ratio r∗ and structural function f∗. We denote the sample counterpart of L(f, r∗) as
Ln(f, r∗) = 1
n n∑ j=1
( 1
n n∑ i=1 (f∗(Xi)− f(Xi))r∗(Yi, Xi|Zj)
)2 .
Let us denote the distributions of (Yi, Xi) and Zi by P and Q, respectively, and define the L2 risk of a function g with P and Q as ‖g‖2L2(P×Q) = ∫ ∫ g2(w)dPdQ. Let us denote the distribution of
(Yi, Xi, Zi) by O, and define the L2 risk of a function g with O as ‖g‖2L2(O) = ∫ g2(w)dO. We put the following assumptions on the error term εi. Assumption 2. The error term εi is sub-Gaussian random variables. In addition, the distributions of Xi and Zi have probability densities that are finite and bounded away from zero.
Note that the randomness of Yi depends on εi andXi. Define the hypothesis classes of the conditional density ratio r∗ and structural function f∗ as R and F , respectively. Suppose that the hypothesis classes are Vapnik–Chervonenkis (VC) class (for rigorous definition, see Section 2.6 in van der vaart & Wellner (1996)). This class include the true models, and the hypothesises are bounded. Assumption 3. The hypothesis classR is VC class, includes the true model, r∗ ∈ R, and all r ∈ R are uniformly bounded by B > 0. Assumption 4. The hypothesis class F is VC class, includes the true model, f∗ ∈ F , and all f ∈ R are uniformly bounded by B > 0.
We define a measure of complexity of the hypothesis classes: for F , we define the complexity CB(F) := ∫ B 0 √ logN (δ′,F , ‖ · ‖L∞)dδ′ with a covering number N (δ′,F , ‖ · ‖) := inf{N |{fj}Nj=1 s.t. F ⊂ ∪j=1{f |‖f − fj‖ ≤ δ′}} in terms of a sup-norm ‖f‖L∞ = supx |f(x)|. In this following parts, for the conditional density ratio estimation, we derive the bound of the MSE ‖r̂ − r∗‖L2(P×Q); for the structural function estimation, we derive the bound of the projected MSE L(f̂ , r∗), which corresponds to the upper bound of the MSE of ‖f̂ − f∗‖L2(O) (Section 5.3).

heading: 5.1 MSE OF THE CONDITIONAL DENSITY RATIO
text: First, we consider the MSE of the conditional density ratio. For a multilayer perception with ReLU activation function (Definition 1), we show the following bound. The proof is shown in Appendix A.2.
Lemma 1 (MSE of r∗). Suppose that Assumptions 1–3 hold. Let I(r) be a non-negative function on R and I(r∗) <∞. Define RM = {r ∈ R : I(r) ≤ M} satisfying R = ⋃ M≥1RM . Suppose that there exist c0 > 0 and 0 < γ < 2 such that supg∈RM ‖r − r ∗‖ ≤ c0M and sup r∈RM
‖r−r∗‖L2(P )≤δ ‖r −
r∗‖∞ ≤ c0M for all δ > 0, and that logN (δ,RM , ‖ · ‖L∞) = O (M/δ)γ . Then, ‖r̂ − r∗‖L2(O) = Op ( n−1/(2+γ) ) (5)
5.2 PROJECTED MSE OF STRUCTURAL FUNCTION f∗
Next, we consider bounding the projected MSE of the structural function f∗. To bound the projected MSE, we use a technique associated with U-statistics. We first obtain the following lemma. Let (Yi, Xi) be Wi, and W ′i be the i.i.d. copy of Wi. The proof is shown in Appendix A.3.
Lemma 2. Suppose that Assumptions 1–4 hold. For any f ∈ F and r ∈ R, L(f, r) = L̃(f, r)+op(1) holds as a freely chosen hyper-parameter as n→∞, where
L̃(f, r) = EW,W ′,Z [(f∗(Xi)− f(Xi))r(Wi|Zi)(f∗(X ′i)− f(X ′i))r(W ′i |Zi)] ,
Regarding the form, we define an empirical version of L̃(f, r) in an U-statistic form:
L̃n(f, r) = 1
n n∑ j=1
1
n(n− 1) n∑ i,i′=1,i6=i′ (f∗(Xi)− f(Xi))r(Wi|Zj)(f∗(Xi′)− f(Xi′))r(Wi′ |Zj).
By a property of U-statistics (for example, see Arcones & Giné (1993)), Ln(f, r) = L̃n(f, r) + Op (1/n) clearly holds. Then, we can decompose the projected MSE L(f̂ , r∗) as follows:
L(f̂ , r∗) = L̃(f̂ , r∗)− L̃(f̂ , r̂)︸ ︷︷ ︸ =:∆r + L̃(f̂ , r̂)− L̃n(f̂ , r̂)︸ ︷︷ ︸ =:∆f +L̃n(f̂ , r̂) +Op (1/n) . (6)
We can handle ∆r by the the estimation error of r∗ as shown in Lemma 6 in Appendix A.4. We evaluate the projected MSE by combining (6) with Lemma 1 and Lemma 6: Theorem 1. Assume that the conditions of Lemmas 1–2 hold. Then, for any any δ ∈ (0, 1), there exists a constant c > 0 such that the following inequality holds with n ≥ 1 and with probability at least 1− δ, for any γ ∈ (0, 1):
L(f̂ , r∗) ≤ cCB(F) + CB(R)√ n +O
( max {√ log(1/δ)
n ,
1
n1/(2+γ)
}) .
This result reveals the following two findings: (i) the projected MSE is affected separately by the complexity of F andR, and (ii) the overall convergence is O(1/ √ n) when these complexities are finite. For instance, when using neural networks with Definition 1, the assumptions on the hypothesis classes are satisfied (Lemma 3).
5.3 FROM PROJECTED MSE TO MSE OF f∗
Following Chen & Pouzo (2012) and Dikkala et al. (2020, Appendix C.4), we discuss the derivation of the upper bound of the MSE of f∗ from the projected MSE. Let us define the measure of illposedness (Chen & Pouzo, 2012; Dikkala et al., 2020) with respect to the function class F as τ := supf∈R ‖f−f∗‖L2(O) L(f,r∗) . Then, the MSE of f ∗ is upper bounded as ‖f − f∗‖2L2(O) ≤ τ 2L(f̂ , r∗) (Chen & Pouzo, 2012; Dikkala et al., 2020). Thus, Theorem 1 also implies
‖f̂ − f∗‖2L2(O) ≤ cτ 2 CB(F) + CB(R)√
n +O
( τ2 max {√ log(1/δ)
n ,
1
n1/(2+γ)
}) .
Note that Dikkala et al. (2020) derived 1/n rate on the MSE of f∗. Our method cannot achieve this rate because of the convergence rate of the plug-in estimator. However, this does not mean the inferiority of the proposed method. For instance, from the computational point of view, our method is easily implemented without using minimax optimization, which is often unstable.

heading: 6 EXPERIMENTS
text: We implement the following three methods based on our proposed method: first, we use neural networks to predict f∗ and train the model by penalized least-squares in (4) (IW-LS); second, we use neural networks to predict f∗ and train the model by minimizing the sum of approximated moment restrictions in (2) (IW-MM), which is the same as IW-LS except for the penalized term in the IW-LS; third, we use a linear-in-parameter model with the Gaussian kernel to predict f∗ and train the model by GMM (IW-Krnl). For all cases, we use neural networks for estimating r∗. We compare our proposed methods with four methods: DeepGMM (Bennett et al. (2019)), DFIV (Xu et al. (2021a)), DeepIV (Hartford et al. (2017)), and KIV (Singh et al. (2019)). We use the datasets proposed in Newey & Powell (2003), Ai & Chen (2003), and Hartford et al. (2017). In addition, we train neural networks by simple least squares (LS), ignoring the dependency between Xi and εi, as a comparison. To fairly evaluate the performances, for DeepGMM, DFIV, DeepIV, and KIV, we use the code and hyperparameters used in Xu et al. (2021a)1. For IW-LS, IW-MM, IW-Krnl, and LS, we also follow the model and hyperparameters of the code as possible. More details are shown in Appendix B.

heading: 6.1 EXPERIMENTS WITH DATASETS OF NEWEY & POWELL (2003) AND AI & CHEN (2003)
text: First, we investigate the performances of the proposed methods using econometric settings of Newey & Powell (2003) and Ai & Chen (2003). These settings have simpler structures than recently proposed settings, such as in Hartford et al. (2017). When using complex and high-dimensional datasets, there is an inherent difficulty in learning due to its complexity, separate from the problem setting. Therefore, we use simple datasets to check whether the proposed method can actually learn f∗.
Newey & Powell (2003) generates {(Yi, Xi, Zi)}ni=1 as follows: first, they generate {(εi, Ui, Zi)}ni=1
from the multivariate normal distribution N (( 0 0 0 ) , (( 1 0.5 0 0.5 1 0 0 0 1 ))) ; then, they generate Xi = Zi + Ui and Yi = f∗(Xi) + εi, where f∗(Xi) = ln(|Xi − 1|+ 1) sgn(Xi − 1). 1https://github.com/liyuan9988/DeepFeatureIV
Ai & Chen (2003) generates {(Yi, Xi, Zi)}ni=1 as follows: first, they generate {(εi, X1i, Vi, Ui)}ni=1 as εi ∼ N ( 0, X21i + V 2 i ) , X1i i.i.d.∼ Unif[0, 1], Vi i.i.d.∼ Unif[0, 1], and Ui i.i.d.∼ N ( 0, X21i + V 2 i ) ; second, they generate X2i = X1i + Vi + R × εi + Ui and Y1 = X1iγ0 + h0 (X2i) + εi, where h0 (X2i) = exp (X2i) / (1 + exp (X2i)) and R is chosen as 0.9; then, obtain Xi = (X1i X2i)> and Zi = (X1i Vi). Here, εi and Ui are unobservable, and f∗(Xi) = X1iγ0 + h0 (X2i), where the function h0 and γ0 are unknown.
We run each algorithm 20 times on each dataset with n = 1, 000 and calculate the mean squared error (MSE). In the left graph of Figure 1 and Figure 2, we report the MSEs. In the dataset of Newey & Powell (2003), DFIV and our proposed IW-MM produce smaller MSE. The dataset only has a one-dimensionalXi, which mitigates the identification problem of neural networks. In contrast, in the dataset of Ai & Chen (2003), IW-LS leads to smaller MSE. In Appendix B.3, using the experimental setting of Newey & Powell (2003), we also show additional experimental results on the empirical convergence of the MSE of the IWMM and a comparison between the IWMM and classical 2SLS.
Because the dimensions of the IVs are low in the original settings, we add more IVs to the original settings and investigate the performances of the algorithms. The detailed settings are described in Appendix B.1. The results are shown in the right graph of Figures 1 and Figure 3.

heading: 6.2 SIMULATION STUDIES USING DEMAND DESIGN DATASETS
text: To investigate the performance with more complicated and high-dimensional datasets, we use the demand design dataset for synthetic airline ticket sales proposed by Hartford et al. (2017). In this dataset, we observe (Yi, Pi, Ti, Si, Ci), where Yi is sales, Pi is price, Ti is time, Si is consumer’s emotion, and Ci is cost to use as IV. Here, Ti and Si are covariates, that is, Xi = (Ti, Si), and the IV is Zi = (Ci, Ti, Si). The sales Yi is generated as Yi = 100 + (10 + Pi)Sih(Ti) − 2Pi + εi, where h(t) = 2 ( (t−5)4 600 + exp ( −4(t− 5)2 ) + t10 − 2 ) . Since Pi is an endogenous variable and
correlated with the IV, we generate Pi to contain Ci as Pi = 25 + (Ci + 3)h(Ti) + Vi. In the simulation, we assume εi ∼ N ( ρVi, 1− ρ2 ) , Vi ∼ N (0, 1), Ci ∼ N (0, 1), Ti is sampled from the uniform distribution with the continuous support [0, 10], and Si is sampled from the uniform distribution with the discrete support {1, . . . , 7}. Here note that E[εi|Ci, Ti, Si] = E[εi|Zi] = 0. The extent of correlation between Pi and Vi is controlled by ρ ∈ {0.25, 0.75} and the larger the ρ is, the more severe the correlation problem becomes.
Figure 4 shows the results with 1, 000 samples. The results under other settings are reporetd in Appendix B.2. In this setting, IW-LS and LS outperform the other methods. We consider this is because f∗(Xi) takes large values compared to the error term. Under this situation, a model trained to predict Yi may perform well because the influence of E[εi|Xi] 6= 0 is limited. However, the purpose of using NPIV in the first place is because the influence is large, or else effects of E[εi|Xi] 6= 0 can be ignored. Therefore, this dataset is often used in the existing studies, it may not be suitable.

heading: 7 CONCLUSION
text: This paper proposed a method for learning causal relationships from conditional moment restrictions. Our method is based on importance weighting using the conditional density ratio. The proposed method showed superior performance in experiments. We point out potential problems in recently proposed methods concerning identification and empirical performances.

heading: A PROOFS IN SECTION 5

heading: A.1 MATHEMATICAL TOOLS
text: Concentration inequalities for empirical processes Given a probability distribution P and a random variable g(X), we denote the expectation of g(X) under P by ∫ gdP . Given samples X1, X2, . . . , Xn from P , the empirical distribution is denoted by Pn, and the empirical mean is denoted by ∫ gdPn; that is, ∫ gdPn = 1 n ∑n i=1 g(Xi). We also denote ∫ gdP − 1n ∑n i=1 g(Xi) by∫
gd ( P − Pn ) .
We also define
ρ2K(g) = 2K 2
∫ ( exp ( |g|/K − 1− |g|/K )) dP, K > 0. (7)
Let G satisfy sup g∈G ρK(g) ≤ R (8)
Then, we summarize some tools used in the theoretical analysis. Proposition 1 (Theorem 5.11 in van de Geer (2000)). C is a (sufficiently large) universal constant, whereas a, C0, and C1 may be chosen, but do have to satisfy the following conditions:
• a ≤ C1 √ nR2/K; • a ≤ 8 √ nR;
• a ≥ C0 ( max {∫ R
a/(26 √ n)
√ logN (u,G, ‖ · ‖)du,R }) • C20 ≥ C2(C1 + 1).

heading: Then,
text: P (
sup g∈G ∣∣∣∣√n ∫ gd(Pn − P )∣∣∣∣ ≥ a) ≤ C exp(− a2C2(C1 + 1)R2 )
(9)
Proposition 2 (Lemma 5.14 in van de Geer (2000)). Let G ⊂ L2(P ) be a function class and the map I(g) be a complexity measure of g ∈ G, where I is a non-negative function onR and I(g0) <∞ for a fixed g0 ∈ G. We now define GM = {g ∈ G : I(g) ≤M} satisfying G = ⋃ M≥1 GM . Suppose that there exist c0 > 0 and 0 < γ < 2 such that sup g∈GM ‖g − g0‖ ≤ c0M, sup r∈GM
‖g−g0‖L2(P )≤δ
‖g − g0‖∞ ≤ c0M, for all δ > 0,
and that logN (δ,GM , P ) = O (M/δ)γ . Then, we have
sup g∈G ∣∣∫ (g − g0)d(P − Pn)∣∣ D(g) = Op (1) , (n→∞),
where D(g) is defined by
D(g) = max ‖g − g0‖ 1−γ/2 L2(P ) I(g) γ/2 √ n , I(g) n2/(2+γ)  . Lemma 3 (Lemma 9 of Kato & Teshima (2021)). Let ` : (br, Br)→ R be a ν-Lipschitz continuous function. Let logN ( δ,F , ‖ · ‖L2(P ) ) denote the bracketing entropy ofF with respect to a distribution P . Then, for any distribution P , any γ > 0, any M ≥ 1, and any δ > 0, we have
logN ( δ, ` ◦ HM , ‖ · ‖L2(P ) ) ≤ (s+ 1)(2ν) γ
γ
( M
δ
)γ .
Moreover, there exists c0 > 0 such that for any M ≥ 1 and any distribution P , sup
`◦r∈`◦HM ‖` ◦ r − ` ◦ r∗‖L2(P ) ≤ c0νM,
sup `◦r∈`◦HM
‖`◦r−`◦r∗‖L2(P )≤δ
‖` ◦ r − ` ◦ r∗‖∞ ≤ c0νM, for all δ > 0.
Concentration inequalities for U-statistics. Given a probability distributions P and Q and a random variable g(X,Z), we denote the expectation of g(X,Z) under P and Q by ∫ ∫ gdPdQ; that
is, EX [ EZ [ g(X,Z) ]] = ∫ ∫
gdPdQ. Given samples X1, X2, . . . , Xn from P and Z1, Z2, . . . , Zn from Q, the empirical distributions are denoted by Pn and Qn, and the empirical mean is denoted by∫ ∫ gdPndQn; that is, ∫ ∫ gdPndQn = 1 n ∑n j=1 1 n ∑n i=1 g(Xi, Zj). We also denote ∫ ∫ gdPdQ− 1 n ∑n j=1 1 n ∑n i=1 g(Xi, Zj) by ∫ ∫ gd ( P − Pn ) d ( Q−Qn ) .
We also define
ρ2K(g) = 2K 2
∫ ( exp ( |g|/K − 1− |g|/K )) dP, K > 0. (10)
Let G satisfy sup g∈G ρK(g) ≤ R (11)
Then, we summarize some tools used in the theoretical analysis. Proposition 3 (Concentration inequality on empirical processes with U-statistics: Theorem 5 in Arcones (1995), adjusted to our setting). Suppose that S1, ..., Sn are S-valued i.i.d. random variables and consider a class of symmetric functionsH ⊂ L2(S ×S). Also, suppose that any function inH is uniformly bounded by b > 0 and define σ2 = suph∈HVarS(ES′ [h(S, S′)]) where S′ is an i.i.d. copy of S. If the N (ε,H, ‖ · ‖2) ≤ (A/ε)ν for some A, ν > 0, for any t ≥ c with some c > 0, we obtain
P n1/2 sup h∈H  1n(n− 1) n∑ i,i′=1,i6=i′ h(Si, Si′)− ES,S′ [h(S, S)]  ≥ t 
≤ 8 exp(−t2/217(σ2 + tbn−1/2)) + 8A2ν(σ2 + 2tbn−1/2)−ν exp(−n(σ2 + tbn−1/2/2)/2b2) + 2 exp(−t2/(211bc(σ2 + tbn−1/2))).
Proposition 4 (Bernsteins’ inequality for U-statistics (Hoeffding, 1963; Arcones & Giné, 1993)). Let ‖g‖∞ < c, ∫ ∫ gdPdQ = 0, and σ2 = ∫ ∫ g2dPdQ. Then, for any a′ > 0:
P (∫ ∫ gdPndQn > a ′ ) ≤ exp ( na ′2/2
2σ2 + (2/3)ca′
) .
This implies that for any a > 0, P (∣∣∣∣∫ ∫ g − g0d(P − Pn)d(Q−Qn)∣∣∣∣ > a√n ) ≤ exp ( a2/2 6σ2 + 2ca/ √ n ) .
As well as Suzuki et al. (2008), by applying this result, we can obtain the following result. Proposition 5 (From Proof of Theorem 1 in Suzuki et al. (2008)). Let G ⊂ L2(P ×Q) be a function class and the map I(g) be a complexity measure of g ∈ G, where I is a non-negative function on R and I(g0) < ∞ for a fixed g0 ∈ G. We now define GM = {g ∈ G : I(g) ≤ M} satisfying G = ⋃ M≥1 GM . Suppose that there exist c0 > 0 and 0 < γ < 2 such that
sup g∈GM ‖g − g0‖ ≤ c0M, sup r∈GM
‖g−g0‖L2(P×Q)≤δ
‖g − g0‖∞ ≤ c0M, for all δ > 0,
and that logN (δ,GM , P ×Q) = O (M/δ)γ . Then, we have
sup g∈G ∣∣∫ ∫ (g − g0)d(P − Pn)d(Q−Qn)∣∣ D(g) = Op (1) , (n→∞),
where D(g) is defined by
D(g) = max ‖g − g0‖ 1−γ/2 L2(P×Q)I(g) γ/2 √ n , I(g) n2/(2+γ)  . sup g∈R ∣∣∣∣∫ ∫ (g − g0)d(P − Pn)d(Q−Qn)∣∣∣∣ = Op (1/√n) , (n→∞). (12)

heading: Complexities of neural networks.
text: Definition 1 (ReLU neural networks; Schmidt-Hieber, 2020). For L ∈ N and p = (p0, . . . , pL+1) ∈ NL+2,
F(L, p) :={f : x 7→WLσvLWL−1σvL−1 · · ·W1σv1W0x : Wi ∈ Rpi+1×pi , vi ∈ Rpi(i = 0, . . . , L)},
where σv(y) := σ(y − v), and σ(·) = max{·, 0} is applied in an element-wise manner. Then, for s ∈ N, F ≥ 0, L ∈ N, and p ∈ NL+2, define
H(L, p, s, F ) := {f ∈ F(L, p) : L∑ j=0 ‖Wj‖0 + ‖vj‖0 ≤ s, ‖f‖∞ ≤ F},
where ‖ · ‖0 denotes the number of non-zero entries of the matrix or the vector, and ‖ · ‖∞ denotes the supremum norm. Now, fixing L̄, p̄, s ∈ N as well as F > 0, we define
IndL̄,p̄ := {(L, p) : L ∈ N, L ≤ L̄, p ∈ [p̄]L+2}, and we consider the hypothesis class
H̄ := ⋃
(L,p)∈IndL̄,p̄
H(L, p, s, F )
H := {r ∈ H̄ : Im(r) ⊂ (br, Br)}.
Moreover, we define I1 : IndL̄,p̄ → R and I : H → [0,∞) by
I1(L, p) := 2|IndL̄,p̄| 1 s+1 (L+ 1)V 2,
I(r) := max ‖r‖∞, min(L,p)∈IndL̄,p̄ r∈H(L,p,s,F ) I1(L, p)  , where V := ∏L+1 l=0 (pl + 1), and we define
HM := {r ∈ H : I(r) ≤M}. Lemma 4 (Lemma 5 in Schmidt-Hieber (2020)). For L ∈ N and p ∈ NL+2, let V := ∏L+1 l=0 (pl + 1).

heading: Then, for any δ > 0,
text: logN (δ,H(L, p, s,∞), ‖ · ‖∞) ≤ (s+ 1) log(2δ−1(L+ 1)V 2). Definition 2 (Derived function class and bracketing entropy). Given a real-valued function class F , define ` ◦ F := {` ◦ f : f ∈ F}. By extension, we define I : ` ◦ H → [1,∞) by I(` ◦ r) = I(r) and ` ◦ HM := {` ◦ r : r ∈ HM}. Note that, as a result, ` ◦ HM coincides with {` ◦ r ∈ ` ◦ H : I(` ◦ r) ≤M}.
Notations. Let us denote a pair of random variables (Yi, Xi) by Wi. We denote the distribution of (Wi, Zi) by R and its empirical distribution as On. Besides, we denote the distributions of Wi and Zi by P and Q, and their empirical distributions by Pn and Qn. For a function g(X,Z), we define the L2 risk over the distribution P and Q as∥∥g∥∥
L2(P×Q) =
√ EX [ EZ [ g2(X,Z) ]] = √∫ ∫ g2dPdQ.

heading: A.2 PROOF OF LEMMA 1: ESTIMATION ERROR IN CONDITIONAL DENSITY RATIO ESTIMATION
text: Lemma 5 (Decomposition of MSE).
‖r̂ − r∗‖2L2(O) ≤ ∣∣∣∣∣EW,Z [r∗(W |Z)− r̂(W |Z)]− 1n n∑ i=1 ( r∗(Wi|Zi)− r̂(Wi|Zi) )∣∣∣∣∣ + 1
2 ∣∣∣∣∣EZ [EW [r∗2(W |Z)− r̂2(W |Z)]]− 1n n∑ j=1 1 n n∑ i=1 ( r∗2(Wi|Zj)− r̂2(Wi|Zj) )∣∣∣∣∣.
Proof of Lemma 5. Since The estimator r̂ is the minimizer of the empirical risk, it satisfies the inequality
− 1 n n∑ i=1 r̂(Wi|Zi) + 1 2 1 n n∑ j=1 1 n n∑ i=1 r̂2(Wi|Zj) ≤ − 1 n n∑ i=1 r∗(Wi|Zi) + 1 2 1 n n∑ j=1 1 n n∑ i=1 r∗2(Wi|Zj).
(13)
Additionally, we consider an expectation of squared residuals as 1 2 EZ [ EW [( r̂(W |Z)− r∗(W |Z) )2]] = 1
2 EZ [ EW [ r̂2(W |Z)− 2r̂(W |Z)r∗(W |Z) + r∗2(W |Z) ]] = 1
2 EZ [ EW [ r̂2(W |Z) ]] − EW,Z [r̂(W |Z)] + 1 2 EW,Z [r∗(W |Z)]
+ EW,Z [r∗(W |Z)]− EW,Z [r∗(W |Z)]
= 1 2 EZ [ EW [ r̂2(W |Z) ]] − EW,Z [r̂(W |Z)]− 1
2 EW,Z [r∗(W |Z)]︸ ︷︷ ︸ =EZ [EW [r∗2(W |Z)]] +EW,Z [r∗(W |Z)]
= −1 2 EZ [ EW [ r∗2(W |Z)− r̂2(W |Z) ]] + EW,Z [r∗(W |Z)− r̂(W |Z)] . (14)
Taking sum of the both hand sides of (13) and (14), then we obtain the following by subtracting − 1n ∑n i=1 r̂(Wi|Zi) + 1 2 1 n ∑n j=1 1 n ∑n i=1 r̂
2(Wi|Zj) from the both side as 1 2 EW,Z [( r̂(W |Z)− r∗(W |Z) )2] ≤ −1
2 EZ [ EW [ r∗2(W |Z)− r̂2(W |Z) ]] + EW,Z [r∗(W |Z)− r̂(W |Z)]
+ 1
2
1
n n∑ j=1 1 n n∑ i=1 ( r∗2(Wi|Zj)− r̂2(Wi|Zj) ) − 1 n n∑ i=1 ( r∗(Wi|Zi)− r̂(Wi|Zi) ) .
Then, we obtain the statement.
Our remaining task is to bound the following target values:
sup r∈R ∣∣∣∣∣∣EZ [EW [r2(W |Z)]]− 1n n∑ j=1 1 n n∑ i=1 r2(Wi|Zj) ∣∣∣∣∣∣ , and (15) sup r∈R ∣∣∣∣∣EW,Z [r(W |Z)]− 1n n∑ i=1 r(Wi|Zi)
∣∣∣∣∣ . (16) Proof of Lemma 1. Since 0 < γ < 2, we can apply Propositions 2 and 5 in combination with Lemma 3 to obtain
sup r∈H
|EZ [ EW [ r2(W |Z) ]] − 1n ∑n j=1 1 n ∑n i=1 r
2(Wi|Zj)| D1(r) = Op (1) ,
sup r∈H ∣∣EW,Z [r(W |Z)]− 1n∑ni=1 r(Wi|Zi)∣∣ D2(r) = Op (1) ,
where
D1(r) = max ‖r 2 − r∗2‖1−γ/2L2(P×Q)I(r 2)γ/2 √ n , I(r2) n2/(2+γ)  , and D2(r) = max ‖r − r ∗‖1−γ/2L2(O) I(r) γ/2 √ n , I(r) n2/(2+γ)  .
Noting that supr∈H I(r) <∞ and supr∈H I(r2) <∞. Then, for any 0 < γ < 2, we have∥∥∥∥∥∥EZ [EW [r̂2(W |Z)]]− 1n n∑ j=1 1 n n∑ i=1 r̂2(Wi|Zj) ∥∥∥∥∥∥ 2
L2(P×Q)
. Op max ‖r̂ 2 − r∗2‖1−γ/2L2(P×Q)√ n , 1 n2/(2+γ)   , and
∥∥∥∥∥EW,Z [r̂(W |Z)]− 1n n∑ i=1 r̂(Wi|Zi) ∥∥∥∥∥ 2
L2(P×Q)
. Op max ‖r̂ − r ∗‖1−γ/2L2(O)√ n , 1 n2/(2+γ)   .
These inequalities imply ‖r̂ − r∗‖L2(P×Q) ≤ Op ( n− 1 2+γ ) .

heading: A.3 PROOF OF LEMMA 2
text: We start with the expected loss function for a learner f : X → R and a conditional density function r : X × Y × Z → R as
L(f, r) := EZ [ (EW [(f∗(X)− f(X))r(W |Z)])2 ] , (17)
where W = (Y,X). For T ∈ T, we consider T i.i.d. random variables Wi = (Xi, Yi), i = 1, ..., n generated from a conditional distribution P and approximate the inner expectation term as
EW [(f∗(X)− f(X))r(W |Z)] = n−1 n∑ i=1 (f∗(Xi)− f(Xi))r(Wi|Z) + κn(f, r),
where κn(f, r) := EW [(f∗(X)− f(X))r(W |Z)]− n−1 n∑ i=1 (f∗(Xi)− f(Xi))r(Wi|Z)
is a residual. We substitute this form into (3) and obtain
L(f, r) = EZ (n−1 n∑ i=1 (f∗(Xi)− f(Xi))r(Wi|Z) + κn(f, r) )2 = EZ
n−2 T∑ k,k′=1 (f∗(Xi)− f(Xi))r(Wi|Z)(f∗(Xi′)− f(Xi′))r(Wi′ |Z)  + 2EZ [ κn(f, r)n
−1 n∑ i=1 (f∗(Xi)− f(Xi))r(Wi|Z) ] ︸ ︷︷ ︸
=:T1
+EZ [κn(f, r)2]︸ ︷︷ ︸ =:T2
= EZ n−2 n∑ i,i′=1,i′ 6=i (f∗(Xi)− f(Xi))r(Wi|Z)(f∗(Xi′)− f(Xi′))r(Wi′ |Z)  ︸ ︷︷ ︸
=:T ∗
+ EZ [ n−2
n∑ i=1
(f∗(Xi)− f(Xi))2r(Wi|Z)2 ]
︸ ︷︷ ︸ =:T0
+T1 + T2.
For T0, the uniformly bounded property of f∗, f and r implies that it is T0 = O(1/n). For T1 and T2, the VC-class property and the Glivenko-Cantelli theorem (Section 2.8.1 in van der vaart & Wellner (1996)), we obtain κn(f, r) → 0 in probability as T → ∞ uniformly on F × R. Associate with the boundedness of f∗, f and r, we obtain that T1 = op(1) and T2 = op(1). About T ∗, which is an U-statistic, a convergence theorem (Theorem 3.1 in Arcones & Giné (1993)) as n→∞ implies that the U-statistic converges to EW,W ′|Z [(f∗(X)− f(X))r(W |Z)(f∗(X ′)− f(X ′))r(W ′|Z)], where W ′ = (Y ′, X ′) is an i.i.d. copied random element of W . Hence, we obtain
L(f, r) = EZ [ EW,W ′|Z [(f∗(X)− f(X))r(W |Z)(f∗(X ′)− f(X ′))r(W ′|Z)] ] + op(1)
= EW,W ′,Z [(f∗(X)− f(X))r(W |Z)(f∗(X ′)− f(X ′))r(W ′|Z)] + op(1),
for any f ∈ F and r ∈ R as T →∞.
A.4 ESTIMATION ERROR BOUND OF ∆r
We can bound ∆r as follows:
Lemma 6. Suppose that Assumptions 1 holds. Then,
∆r = O(‖r̂ − r∗‖L2(P×Q)) (18)
Proof. Let us define ξf̂ := f̂ − f ∗. We bound ∆r by the Lipschitz continuity of L(f̂ , r) in r. For any r, r′ ∈ R, we obtain
L(f̂ , r)− L(f̂ , r′) = E [ ξ2 f̂ (X)ξ2 f̂ (X ′)(r(W |Z)r(W ′|Z)− r′(W |Z)r′(W ′|Z)) ] = E [ ξ2 f̂ (X)ξ2 f̂ (X ′)r(W ′|Z)(r(W ′|Z)− r′(W ′|Z))
] + E [ ξ2 f̂ (X)ξ2 f̂ (X ′)r′(W |Z)(r(W ′|Z)− r′(W ′|Z))
] . 32B ′5‖r − r′‖L2(O),
by the Cauchy-Schwartz inequality and the uniformly bounded property over F andR. As a result, if the conditional density ratio p(w,z)p(w)p(z) is bounded, we obtain
∆r = O(‖r̂ − r∗‖2L2(O)) = O(‖r̂ − r ∗‖2L2(P×Q)).

heading: A.5 PROOF OF THEOREM 1: ESTIMATION ERROR BOUND
text: We use the following lemma to show Theorem 1.
Lemma 7. L̃n(f̂ , r̂) ≤ Γ(f̂ , r̂) +Op (1/n) .
Proof. We start with the following basis inequality following the definition of f̂ ; since it is an minimizer of the empirical risk with r̂, we obtain
1
n n∑ j=1
( 1
n n∑ i=1 (Yi − f̂(Xi))r̂(Wi|Zj) )2 ≤ 1 n n∑ j=1 ( 1 n n∑ i=1 (Yi − f∗(Xi))r̂(Wi|Zj) )2 .
Since Yi = f∗(Xi) + ε, the above inequality updated as
1
n n∑ j=1
( 1
n n∑ i=1 (εi − ξf̂ (Xi))r̂(Wi|Zj) )2 ≤ 1 n n∑ j=1 ( 1 n n∑ i=1 εir̂(Wi|Zj) )2 , (19)
where we define ξf̂ := f̂ − f ∗. Since the left-hand side is expanded as
1
n n∑ j=1
( 1
n n∑ i=1 (εi − ξf̂ (Xi))r̂(Wi|Zj)
)2
= 1
n n∑ j=1
( 1
n n∑ i=1 εir̂(Wi|Zj)
)2 + 1
n n∑ j=1
( 1
n n∑ i=1 ξf̂ (Xi)r̂(Wi|Zj) )2 ︸ ︷︷ ︸
=Ln(f̂ ,r̂)+Op(1/n)
− 2 1 n n∑ j=1 1 n2 n∑ i,i′=1,i6=i′
εiεi′ξf̂ (Xi)ξf̂ (Xi′)r̂(Wi|Zj)r̂(Wi′ |Zj)︸ ︷︷ ︸ =:Γ(f̂ ,r̂)
− 2 1 n3 n∑ i,j=1
ε2i ξf̂ (Xi) 2r̂(Wi|Zj)2︸ ︷︷ ︸
=Op(1/n)
.
Substituting this expansion into (19), we obtain
Ln(f̂ , r̂) ≤ Γ(f̂ , r̂) +Op (1/n) ≤ sup f∈F,r∈R Γ(f, r) +Op (1/n) .
Therefore,
L̃n(f̂ , r̂) = Ln(f̂ , r̂) +Op (1/n) = Γ(f̂ , r̂) +Op (1/n) ≤ sup f∈F,r∈R Γ(f, r) +Op (1/n) .
Then we show Theorem 1.
Step (i): Bound of L̃n(f̂ , r̂). A goal of this step is to show that the definition of f̂ implies that for some ≥ 0, with high probability,
L̃n(f̂ , r̂) ≤ .
From Lemma 7,
L̃n(f̂ , r̂) ≤ Γ(f̂ , r̂) +Op (1/n) ≤ sup f∈F,r∈R Γ(f, r) +Op (1/n) .
A rest of this step is to bound supf∈F,r∈R Γ(f, r). We bound the tail-probability by the Talagrand’s inequality (Theorem 3.3.9 in Giné & Nickl (2021)) for U-statistics (Theorem 5 in Arcones (1995), or Theorem 1 in Li et al. (2014)) as the following inequality; for any δ > 0, with probability at least 1− δ, we obtain
sup f∈F,r∈R
Γ(f, r) . ED [ sup
f∈F,r∈R Γ(f, r)
] +O (√ log(1/δ)
n
) ,
provided the uniformly bounded and finite variance properties over F ×R.
To the end, we study an expectation of Γ(f̂ , r̂) in terms of the dataset D = {(Yi, Xi, Zi)}ni=1 and obtain
ED [ sup
f∈F,r∈R Γ(f, r)
] ≤ ED [Γ(f∗, r∗)] + E [ sup
f,f ′∈F,r,f ′∈R |Γ(f, r)− Γ(f ′, r′)|
]
. 0 + n−1/2 ∫ B
0
√ logN (δ,F ×R, ‖ · ‖)dδ.
Here, for the joint set F ×R, we define a distance ‖(f, r)‖ := ‖f‖+‖r‖. The last inequality follows Γ(f∗, r∗) = 0 and the sub-Gaussianity inequality (Corollary 2.2.8 in van der vaart & Wellner (1996)) associated with the sub-Gaussianity of εi and the Lipschitz continuity of Γ(f, r) in (f, r).
By these results, we obtain
Ln(f̂ , r̂) ≤ n−1/2 ∫ B
0
√ logN (δ,F ×R, ‖ · ‖)dδ +O
(√ log(1/δ)
n
)
. n−1/2 ∫ B
0
√ logN (δ′,F , ‖ · ‖)dδ′ + n−1/2 ∫ B 0 √ logN (δ′,R, ‖ · ‖)dδ′ +O (√ log(1/δ) n ) ,
which follows logN (δ,F ×R, ‖ · ‖) . logN (δ,F , ‖ · ‖) + logN (δ,R, ‖ · ‖).
Step (ii): Bound of ∆f . To bound ∆f = L̃(f̂ , r∗) − L̃n(f̂ , r∗), we apply the concentration inequality on empirical processes (displayed in Proposition 3, which is originally developed in Theorem 5 in Arcones (1995)). We regard h in Proposition 3 as (y, x, y′, x′) 7→ (y − f∗(x))(y′ − f∗(x′))∆f (x)∆f (x ′)r(y, x|Zj)r(y′, x′|Zj) and obtain
∆f = EW,W ′,Z [(f∗(X)− f(X))r(W |Z)(f∗(X ′)− f(X ′))r(W ′|Z)]
− 1 n n∑ j=1
1
n(n− 1) n∑ i,i′=1,i6=i′ (f∗(Xi)− f(Xi))r(Wi|Zj)(f∗(Xi′)− f(Xi′))r(Wi′ |Zj)
.
√ log(1/δ)
n ,
with probability at least δ for any δ ∈ (0, 1).
Step (iii): Conclusion. Finally, by combining the above results, with probability at least 1− δ with δ > 0,
L(f̂ , r∗) = L̃(f̂ , r∗)− L̃n(f̂ , r∗)︸ ︷︷ ︸ =:∆f + L̃n(f̂ , r∗)− L̃n(f̂ , r̂)︸ ︷︷ ︸ =:∆r +L̃n(f̂ , r̂)
≤ ∆f + ∆r + n−1/2 ∫ B
0
√ logN (δ′,F ×R, ‖ · ‖)dδ′ +O
(√ log(1/δ)
n
)
. O(‖r̂ − r∗‖) + n−1/2 ∫ B
0
√ logN (δ,F , ‖ · ‖)dδ
+ n−1/2 ∫ B
0
√ logN (δ,R, ‖ · ‖)dδ +O
(√ log(1/δ)
n
)
= n−1/2 ∫ B
0
√ logN (δ,F , ‖ · ‖)dδ + n−1/2 ∫ B 0 √ logN (δ,R, ‖ · ‖)dδ
+O ( max {√ log(1/δ)
n ,
1
n1/(2+γ)
}) (from Lemma 1).

heading: B EXPERIMENTAL SETTINGS AND ADDITIONAL RESULTS

heading: B.1 SIMULATION STUDIES USING ECONOMICS DATASETS
text: Here, we report the details of the model and hyperparameters used in our experiments. Since the network structures used in Ai & Chen (2003) and Newey & Powell (2003) have many common features, we describe the network structure in Newey & Powell (2003).
For DeepGMM, DFIV, DeepIV, and KIV, we use the same model and hyperparameters published by Xu et al. (2021a). For IW-LS, we use the following two networks for estimating density ratio and predicting f∗ respectively : FC(4,128)-FC(128,128)-FC(128,1) and FC(1,128)-FC(128,128)FC(128,1). Each fully-connected layer (FC) is followed by leaky ReLU activations with leakiness
α = 0.2. A regularization coefficient η is set to 0.001 as a result of cross-validation. For LS, we use the same network structure in IW-LS for predicting f∗. For IW-MM, we use the same network structure as IW-LS. For IW-Krnl, the same network structure as IW-LS to estimate the density ratio and ζ and σ2 are selected via cross-validation. In Ai & Chen (2003) experiment, we change only the first layer in the network structure to match the dimension of Xi, and the rest of the network structure is the same.
In addition to the original settings of Newey & Powell (2003) and Ai & Chen (2003), we investigate cases where we add more IVs. In Newey & Powell (2003), we generate {(Yi, Xi, Zi,Wi,1,Wi,2,Wi,3)}ni=1, where Wi,1,Wi,2,Wi,3 are additional IVs, as follows: first, they generate {(εi, Ui, Zi,Wi,1,Wi,2,Wi,3)}ni=1 from the multivariate normal distribution
N   0 0 0 0 0 0  ,  1 0.5 0 0 0 0 0.5 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1  ; then, they generateXi = Zi+ ∑3 j=1Wi,j+Ui and Yi = f∗(Xi)+εi, where f∗(Xi) = ln(|Xi−1|+1) sgn(Xi−1). Here, εi and Ui are unobservable. In Ai & Chen (2003), they generate {(Yi, Xi, Zi,Wi,1,Wi,2,Wi,3)}ni=1, where Wi,1,Wi,2,Wi,3 are additional IVs, as follows: first, they generate {(εi, X1i, Vi, Ui)}ni=1 as εi ∼ N ( 0, X21i + V 2 i ) , X1i i.i.d.∼
Unif[0, 1], Vi i.i.d.∼ Unif[0, 1], ( Wi,1 Wi,2 Wi,3 ) i.i.d.∼ N (( 0 0 0 ) , (( 1 0.3 0.3 0.3 1 0.3 0.3 0.3 1 ))) , Wi = ∑3 j=1Wi,j , and Ui i.i.d.∼ N ( 0, X21i + V 2 i + |Wi| ) ; second, they generate X2i = X1i + Vi +Wi +R× εi + Ui and Y1 = X1iγ0 + h0 (X2i) + εi, where h0 (X2i) = exp (X2i) / (1 + exp (X2i)) and R is chosen as 0.9; then, obtain Xi = (X1i X2i)> and Zi = (X1i Vi Wi,1 Wi,2 Wi,3). Here, εi and Ui are unobservable, and f∗(Xi) = X1iγ0 + h0 (X2i), where the function h0 and γ0 are unknown.

heading: B.2 SIMULATION STUDIES USING DEMAND DESIGN DATASETS
text: For DeepGMM, DFIV, DeepIV, and KIV, we use the exact same model and hyperparameters in Xu et al. (2021a). For LS and IW-LS, we use the network structure based on the economics datasets network. η is 0.001 as a result of cross-validation. For IW-Krnl, we use the same network structure as IW-LS to estimate the density ratio and select ζ and σ2 via cross-validation.
Figure 5 shows the results with 5000 samples. First, our proposed IW-LS outperforms the existing methods and minimizes the MSE. Second, LS, which is a naive nonlinear regression without IV, outperforms the other existing methods using IV. In contrast, IW-LS outperforms LS.
We consider the case where endogeneity causes a larger change in outcome. Since the original demand design dataset has a small effect of the bias, we slightly change Vi in the price equation: Pi = 25 + (Ci + 3)h(Ti) + 10Vi. In the original dataset, Yi is generated as Yi = 100 + (10 + Pi)Sih(Ti) − 2Pi + εi. However, because the impact of εi on the price is very limited because the variance of εi is relatively small compared with that of Yi. This is the reason why the LS also performs well in the previous result, that is, we can obtain good performance even when ignoring the endogeneity (because f∗(Xi)) is close to E[Yi|Xi]). For this reason, this dataset is not appropriate for investing the performances of the methods, although it is used in existing studies. Here, we consider the different model Yi = 100+(10 +Pi)Sih(Ti)−2Pi+100εi, in which we multiply the error term by 100. In this case, the bias has a more serious impact on the values of price and outcome. Figure 6 shows the results with 1, 000 samples, and Figure 7 shows the results with 5, 000 samples. Regardless of the sample size and ρ, our proposed method outperforms existing methods. Even in this experiment, the LS still performs well. We consider that this is because the correlation between Xi and εi is not strong. In the dataset, Vi, which is a cause of the bias, follows a standard normal distribution and has a limited impact on price due to the constant term 25 (Pi = 25 + (Ci + 3)h(Ti) + Vi), that is, in the variation of Pi, Vi has a small effect compared with the other variables.
In this dataset, f∗(Xi) takes large values compared to the error term. Under this situation, a model trained to predict Yi may perform well because the influence of E[εi|Xi] 6= 0 is limited. However, the purpose of using NPIV in the first place is because the latter influence is large, or else effects of E[εi|Xi] 6= 0 can be ignored. As expected, in our experiments using Hartford et al. (2017), the least-squares method also performs well, even though the training process ignores the problem of NPIV. To investigate the performance for causal inference, we recommend using simpler datasets, before using more complicated datasets.
We also show how the MSE of the IWMM decreases as the sample size grows in Appendix B.3

heading: B.3 SIMULATION STUDIES ON THE RATE OF THE MSE
text: We show experimental results on the MSEs of the NPIV and the classical 2SLS under various sample sizes. Here, we have two goals: (i) we validate the theoretical MSE derived in Section 5.3, and (ii) we show the relative performance of the NPIV method against the classical 2SLS in the learning problem with conditional moment restrictions.
In Figure 8, we show the MSEs of the IWMM in the setting of Newey & Powell (2003) with various sample sizes. We follow the same setting used in Section 6.1, except for the choices of sample sizes. We investigate whether the empirical MSE of IWMM converges to 0 with Op(1/ √ n) for sample size n, which is theoretically provided in Section 5.3 under some assumptions. We compute the empirical MSEs with sample sizes {500, 1000, 1500, 2000, 3500, 4000, 4500, 5000} with 10 trials, then obtain its mean and standard deviation. In Figure 8, we compare the empirical MSEs with the line that represents 1/ √ n. As in our theoretical results, the MSEs decays in Op(1/ √ n).
With the same setting, we also compare the MSEs of IWMM with the classical 2SLS in Figure 9, to understand the effectiveness of the NPIV method. We consider three models for the 2SLS:
Yi = β0 + β1Xi + εi,
Yi = β0 + β1Xi + β2X 2 i + εi,
Yi = β0 + β1Xi + β2X 2 i + β3X 3 i + εi.
The difference between these three models lies in the choice of the polynomial basis. For instance, the last model approximates the structure function by a cubic function. If a model introduce an infinite number of polynomial bases, they can approximate various smooth functions by the Series expansion. In other words, as the number of the polynomial basis is increased, the classical 2SLS approaches the NPIV. This method of introducing a basis is called sieve regression in econometrics, and Newey & Powell (2003) proposed using it to solve NPIV. The experimental result in Figure 9 also shows that the second and third model, with the polynomial bases, X2i and X 3 i , perform closer to the NPIV method than the first model using only Xi.
It is important to note that the 2SLS with a polynomial basis function is difficult to implement in high-dimensional situations. In the setting of Newey & Powell (2003), both X and Z have one dimension, so polynomial approximation is effective. However, if the dimension increases, the Series expansion of multivariate functions requires a huge number of basis functions, hence approximation becomes very difficult. In machine learning, for instance, Singh et al. (2019) proposes to introduce RKHS to avoid this difficulty.

Review summaries:
The present paper proposes an importance weighting approach to address the issue of regression under conditional moment restrictions in the context of non-parametric instrumental variables. The authors provide error bounds on the learned structural function. They show that their approach has a convergence rate of $\mathcal{O}(1/\sqrt{n})$. The main contributions of this work are presented in Sections 4 and 5, where they first introduce their method based on re-weighting by an estimated conditional density ratio function and then provide an estimation error analysis. The theoretical results are complemented with two synthetic experiments which show that the proposed method can compete with and in certain cases improves upon state of the art.

The paper presents a method for estimating causal effects (or more generally non/parametric functions) under conditional moment restrictions, focusing in this case on nonparametric IVs. The main idea is casting these conditional restrictions to an unconditional version through importance resampling using a conditional density estimator (e.g. based on a least-squares method with a NN). The paper also provides a theoretical analysis of the estimation error, providing am error bound (Thm1). The method is evaluated based on three econometrics datasets from literature, showing a smaller MSE than the compared methods.