---
abstract: 'We consider the problem of predicting a response from a set of covariates when the test distribution differs from the training distribution. Here, we consider robustness against distributions that emerge as intervention distributions. Causal models that regress the response variable on all of its causal parents have been suggested for the above task since they remain valid under arbitrary interventions on any subset of covariates. However, in linear models, for a set of interventions with bounded strength, alternative approaches have been shown to be minimax prediction optimal. In this work, we analyze minimax solutions in nonlinear models for both direct and indirect interventions on the covariates. We prove that the causal function is minimax optimal for a large class of interventions. We introduce the notion of distribution generalization, which is motivated by the fact that, in practice, minimax solutions need to be identified from observational data. We prove sufficient conditions for distribution generalization and present corresponding impossibility results. To illustrate the above findings, we propose a practical method, called NILE, that achieves distribution generalization in a nonlinear instrumental variable setting with linear extrapolation. We prove consistency, present empirical results and provide code.'
author:
- |
    Rune Christiansen$^\flat$ Niklas Pfister$^\flat$ Martin Emil Jakobsen$^{\flat}$\
    Nicola Gnecco$^{\sharp}$ Jonas Peters$^\flat$
bibliography:
- 'ref.bib'
title: ' **The Difficult Task of Distribution Generalization in Nonlinear Models** '
---

Introduction {#sec:intro}
============

Large-scale learning systems, particularly those focusing on prediction tasks, have been successfully applied in various domains of application. Since inference is usually done during training time, any difference between training and test distribution poses a challenge for prediction methods [@Quionero2009; @Pan2010; @Csurka2017; @Arjovsky2019]. Dealing with differences in training and test distribution is of great importance in fields such as many environmental sciences, where methods need to extrapolate both in space and time. Tackling this task requires restrictions on how the distributions may differ, since, clearly, generalization becomes impossible if the test distribution may be arbitrary. Given a response $Y$ and some covariates $X$, existing procedures often aim to find a function $f$ which minimizes the worst-case risk $\sup_{P \in \mathcal{N}} \mathbb{E}_P [(Y - f(X))^2]$ across distributions contained in a small neighborhood $\mathcal{N}$ of the training distribution. The neighborhood $\mathcal{N}$ should be representative of the difference between the training and test distributions, and often mathematical tractability is taken into account, too [@abadeh2015distributionally; @sinha2017certifying]. A typical approach is to define a $\rho$-ball of distributions $\cN_\rho(P_0) := \{P: D(P, P_0) \leq \rho\}$ around the training distribution $P_0$, with respect to some divergence measure $D$, such as the Kullback-Leibler divergence or the $\chi^2$ divergence [@hu2013kullback; @ben2013robust; @bertsimas2018data; @lam2019recovering; @duchi2016statistics]. While these divergence functions only consider distributions with the same support as $P_0$, the Wasserstein distance allows to define a neighborhood of distributions around $P_0$ with possibly different supports [@abadeh2015distributionally; @sinha2017certifying; @esfahani2018data; @blanchet2019data]. In our analysis, we do not start from a divergence measure, but we construct a neighborhood of distributional changes by using the concept of interventions [@pearl2009causality; @Peters2017book].

We will see that, depending on the considered setup, one can find models that perform well under interventions which yield distributions that are considered far away from the observational distribution in any commonly used metric. Using causal concepts for the above problem has been motivated by the following observation. A causal prediction model, that uses only the direct causes of the response $Y$ as covariates, is known to be invariant under interventions on variables other than $Y$: the conditional distribution of $Y$ given its causes does not change (this principle is known as invariance, autonomy or modularity) [@Aldrich1989; @Haavelmo1944; @pearl2009causality]. Such a model yields the minimal worst-case prediction error when considering all interventions on variables other than $Y$ [e.g., @Rojas2016 Theorem 1, Appendix]. It has therefore been suggested to use causal models in problems of domain generalization or distributional shifts [@Scholkopf2012; @Rojas2016; @HeinzeDeml17; @Magliacane2018; @Meinshausen2018; @Arjovsky2019; @pfister2019stabilizing]. One may argue, however, that causal methods are too conservative in that the interventions which induce the test distributions may not be arbitrarily strong. As a result, methods which focus on a trade-off between predictability and causality have been proposed for linear models [@rothenhausler2018anchor; @Pfister2019pnas], see also Section \[sec:existingmethods\]. In this work, we consider the problem of characterizing and finding minimax optimal models in a more general, nonlinear framework. 0

Niklas: Maybe use this from below: In cases with more complicated interventions, the causal function is not necessarily a minimax solution anymore. Two interesting such scenarios are the following.

1.  If ${\mathcal{I}}$ consists only of the trivial intervention (does not change anything), then the [[Martin: any]{}]{} conditional mean function, i.e., $f(x)={\mathbb{E}}(Y\vert X=x)$ [[Martin: $x\mapsto{\mathbb{E}}(Y\vert X=x)$]{}]{}, is a minimax solution.

2.  If ${\mathcal{I}}$ consists of interventions that fix the structural form of $X$ apart from adding shifts to $A$, then (under additional assumptions on $A$) the causal function will approximate a minimax solution as if the allowed shift strength increases to infinity, see Section \[sec:learnability\]. [[Martin: If referring to AR, then maybe change to we add shifts to $X$ not $A$, or simply intervene on $A$ such that $Cov(A,{\varepsilon})=0$. ]{}]{}

Contribution {#sec:contr}
------------

We assume that the true data generating process can be described by a model $M$ that belongs to a class of models ${\mathcal{M}}$ and induces an observational distribution ${\mathbb{P}}_{M}$. We then consider the risk of a prediction function ${f_{\diamond}}$ from a function class ${\mathcal{F}}$ under a modified model $M(i)$ that is obtained from $M$ by an intervention $i$, which belongs to a set of interventions ${\mathcal{I}}$. Here, interventions can either act directly on $X$ or indirectly, via an exogenous variable $A$, if the latter exists (precise definitions are provided in Section \[sec:setup\] below). Our work has four main contributions. (1) We analyze the relation between the causal function (defined formally in Section \[sec:setup\]) and the minimizer of $\sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{{M}(i)} [ (Y - {f_{\diamond}}(X))^2 ]$. Our findings go beyond existing results in that the causal function is shown to be minimax optimal already for relatively small intervention classes. We further prove that, in general, the difference between a minimax solution and the causal function can be bounded and that any minimax solution different from the causal function is not robust with respect to misspecification of the intervention class. (2) In practice, we usually have to learn the minimax solution from an observational distribution, in the absence of causal background knowledge. We therefore introduce the concept of distribution generalization, which requires the existence of a prediction model $f^*$ which (approximately) solves the minimax problem $\operatorname*{argmin}_{{f_{\diamond}}\in {\mathcal{F}}} \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}[ (Y - {f_{\diamond}}(X))^2]$ for all $\tilde{M}$ with ${\mathbb{P}}_M = {\mathbb{P}}_{\tilde{M}}$. To the best of our knowledge, the considered setup is novel. (3) We then investigate explicit conditions on $\mathcal{M}$, $\mathcal{I}$ and ${\mathbb{P}}_M$ that allow us to use the observational distribution of $(X, Y, A)$ to identify a function $f^*:{\mathbb{R}}^d\rightarrow{\mathbb{R}}$ that generalizes to ${\mathcal{I}}$, i.e., it (approximately) solves the above minimax problem. We prove several results. E.g., if the interventions are such that the support of $X$ does not increase with respect to the training distribution, then identifiability of the causal function — a well-studied problem in causality — is in general sufficient for generalization. We furthermore give sufficient conditions for generalization to interventions on either $A$ or $X$ that extend the support of $X$. Table \[tab:generalizability\] summarizes some of these results.

We also prove that, without these assumptions, generalization is impossible; (4) In Section \[sec:learning\], we discuss how minimax functions can be learned from finitely many data and explain how existing methodology fits into our framework. We propose a novel estimator, the NILE, that is applicable in a nonlinear instrumental variables (IV) setting and achieves distribution generalization with linear extensions. We prove consistency and provide empirical results. Our code is available as an `R`-package at <https://runesen.github.io/NILE>. Scripts generating all our figures and results can be found at the same url.

Further related work {#sec:related_work}
--------------------

That the causal function is minimax optimal under the set of all interventions on the covariates has been shown by [@Rojas2016], for example, where the additional assumption of no hidden variables is made. In Section \[sec:setup\], we extend this result in various ways. The question of distributional robustness, sometimes also referred to as out-of-distribution generalization, aims to develop procedures that are robust to changes between training and testing distribution. Empirically, this problem is often studied using adversarial attacks, where small digital [@goodfellow2014explaining] or physical [@evtimov2017robust] perturbations of pictures can deteriorate the performance of a model; arguably, these procedures are not yet fully understood theoretically. Unlike the procedures mentioned in Section \[sec:contr\] that aim to minimize the worst-case risk across distributions contained in a neighborhood of the training distribution, e.g., in the Wasserstein metric, [@sinha2017certifying], we assume these neighborhoods to be generated by interventions. To the best of our knowledge, the characterization of distribution generalization that we consider in Section \[sec:generalizability\] is novel.

In settings of covariate shift, one usually assumes that the training and test distribution of the covariates are different, while the conditional distribution of the response given the covariates remains invariant [@daume2006domain; @bickel2009discriminative; @David10; @muandet2013domain]. Sometimes, it is additionally assumed that the support of the training distribution covers the one of the test distribution [@shimodaira2000improving]. In this work, the conditional distribution of the response given the covariates is allowed to change between interventions, due to the existence of a hidden confounder, and we consider settings where the test observations lie outside the training support. Data augmentation methods increase the diversity of the training dataset by changing the geometry and the color of the images (e.g., by rotation, cropping or changing saturation) [@zhang2017mixup; @shorten2019survey]. This allows the user to create models that generalize better to unseen environments [e.g., @volpi2018adversarial]. We view these approaches as a way to enlarge the support of the covariates, which comes with theoretical advantages, see Section \[sec:generalizability\].

Minimizing the worst-case prediction error can also be formulated in terms of minimizing the regret in a multi-armed bandit problem [@lai1985asymptotically; @auer2002finite; @bartlett2008high]. In that setting, the agent can choose the distribution which generates the data. In our setting, though, we do not assume to have control on the interventions and hence on the distribution of the sampled data.

Structure of this work
----------------------

We introduce our framework for generating a collection of intervention distributions in Section \[sec:setup\]. In Section \[sec:robustness\], we formalize the problem considered in this work, namely to find a model that predicts well under a set of intervention distributions. We prove that for a wide range of intervention classes, this is achieved by the causal function. In reality, we are not given the full causal model, but only the observational distribution. This problem is considered in Section \[sec:generalizability\], where we provide sufficient conditions under which distribution generalization is possible and prove corresponding impossibility results. The condition whether the intervened $X$ values are inside the support of the training distribution will play an important role. Section \[sec:learning\] considers the problem of learning models from a finite amount of data. In particular, we propose a method, called NILE, that learns a generalizing model in a nonlinear IV setting. We prove consistency and compare our method to state-of-the art approaches empirically. In Appendix \[sec:causal\_relations\_X\], we comment on the different model classes that are contained in our framework. Appendix \[sec:IVconditions\] summarizes existing results on identifiability in IV models and Appendix \[sec:test\_statistic\] provides details on the test statistic that we use in NILE. Appendix \[sec:additional\_experiments\] contains an additional experiment and all proofs are provided in Appendix \[app:proofs\].

Modeling intervention induced distributions {#sec:setup}
===========================================

We now specify the statistical model used throughout this paper. For a real-valued response variable $Y\in{\mathbb{R}}$ and predictors $X \in {\mathbb{R}}^d$, we consider the problem of estimating a regression function that works well not only on the training data, but also under distributions that we will model by interventions. We require a model that is able to model an observational distribution of $(X,Y)$ (training) and the distribution of $(X,Y)$ under a class of interventions on (parts of) $X$ (testing). We will do so by means of a structural causal model (SCM) [@Bollen1989; @pearl2009causality]. More precisely, denoting by $H \in {\mathbb{R}}^q$ some additional (unobserved) variables, we consider the SCM

$$\begin{aligned}
    H &\coloneqq {{\varepsilon}}_H \; &{ \text{\scriptsize $q$  assignments} }  \\
    X &\coloneqq h_2(H, {\varepsilon}_X)\; &{\text{\scriptsize $d$  assignments} }\\
    Y &\coloneqq f(X) + h_1(H, {\varepsilon}_Y)\; &{\text{\scriptsize $1$  assignment} }
  \end{aligned}$$

Here, $f$, $h_1$ and $h_2$ are measurable functions, the innovation terms ${{\varepsilon}}_X$, ${{\varepsilon}}_Y$ and ${{\varepsilon}}_H$ are independent vectors with possibly dependent coordinates. Two comments are in order. The joint distribution of $(X, Y)$ is constrained only by requiring that $X$ and $h_1({{\varepsilon}}_Y, H)$ enter the equation of $Y$ additively. This constraint affects the allowed conditional distributions of $Y$ given $X$ but does not make any restriction on the marginal distributions of either $X$ or $Y$. Furthermore, we do not assume that the above SCM represents the true causal relationships between the random variables. We do not assume any causal background knowledge of the system. Instead, the SCM is used only to construct the test distributions (by considering interventions on $X$) for which we are analyzing the predictive performance of different methods – similar to how one could have considered a ball around the training distribution. If causal background knowledge exists, however, e.g., in the form of an SCM over variables $X$ and $Y$, it can be made to fit into the above framework. As such, our framework includes a large variety of models, including SCMs in which some of the $X$ are not ancestors but descendants of $Y$ (this requires adapting the set of interventions appropriately), see Appendix \[sec:causal\_relations\_X\] for details. The following remark shows such an example, and may be interesting to readers with a special interest in causality. It can be skipped at first reading.

\[rem:model\] If a priori causal background knowlededge is available, e.g., in form of an SCM, our framework is still applicable after an appropriate transformation. The following example shows a reformulation of an SCM over variables $X_1$, $X_2$ and $Y$.

[ $$\begin{aligned}
        X_1 &\coloneqq {{\varepsilon}}_1 \\
        X_2 &\coloneqq k(Y) + {\varepsilon}_2\\
        Y &\coloneqq f(X_1) + {\varepsilon}_3,
      \end{aligned}$$ with $({\varepsilon}_1,{\varepsilon}_2, {{\varepsilon}}_3)\sim Q$. ]{}

[ $$\xrightarrow{\text{rewrite}}$$]{}

[ $$\begin{aligned}
        H &\coloneqq {{\varepsilon}}_3 \\
        X &\coloneqq h_2(H,({{\varepsilon}}_1,{{\varepsilon}}_2))\\
        Y &\coloneqq f(X_1) + H,
      \end{aligned}$$ with $({\varepsilon}_1,{\varepsilon}_2, {{\varepsilon}}_3)\sim Q$]{}.

Here, $h_2(H,({{\varepsilon}}_1, {{\varepsilon}}_2))\coloneqq({{\varepsilon}}_1, k(f({{\varepsilon}}_1)+H) + {{\varepsilon}}_2)$. Both SCMs induce the same observational distribution over $(X_1,X_2, Y)$ and any intervention on the covariates in the SCM on the left-hand side can be rewritten as an intervention on the covariates in the SCM on the right-hand side. Details and a more general treatment are provided in Appendix \[sec:causal\_relations\_X\].

Sometimes, the vector $X$ contains variables that are independent of $H$ and that enter additively into the assignments of the other covariates. If such covariates exist, it can be useful to explicitly distinguish them from the other covariates. We will denote them by $A$ and call them exogenous variables. Such variables are interesting for two reasons. (i) Under additional assumptions, they can be used as instrumental variables [e.g., @Bowden1985; @greene2003econometric], a well-established tool for ensuring that $f$ can be uniquely recovered from the observational distribution of $(X, Y)$. And (ii), we will see below that in general, interventions on such variables lead to intervention distributions with desirable properties. In the remainder of this article, we will therefore consider a slightly larger class of SCMs that also includes exogenous variables $A$. It contains the SCM presented at the beginning of Section \[sec:setup\] as a special case.[^1]

Model
-----

Formally, we consider a response $Y \in {\mathbb{R}}^1$, covariates $X\in {\mathbb{R}}^d$, exogenous variables $A\in {\mathbb{R}}^r$, and unobserved variables $H \in {\mathbb{R}}^q$. Let further ${\mathcal{F}}\subseteq\{f:{\mathbb{R}}^d\rightarrow{\mathbb{R}}\}$, ${\mathcal{G}}\subseteq\{g:{\mathbb{R}}^r\rightarrow{\mathbb{R}}^d\}$, ${\mathcal{H}}_1\subseteq\{h_1:{\mathbb{R}}^{q+1}\rightarrow{\mathbb{R}}\}$ and ${\mathcal{H}}_2\subseteq\{h_2:{\mathbb{R}}^{q+d}\rightarrow{\mathbb{R}}^d\}$ be fixed sets of measurable functions. Moreover, let $\mathcal{Q}$ be a collection of probability distributions on ${\mathbb{R}}^{d+1+r+q}$, such that for all $Q\in\mathcal{Q}$ it holds that if $({\varepsilon}_X,{\varepsilon}_Y, {{\varepsilon}}_A, {{\varepsilon}}_H)\sim Q$, then ${\varepsilon}_X,{\varepsilon}_Y, {{\varepsilon}}_A$ and ${{\varepsilon}}_H$ are jointly independent, and for all $h_1\in{\mathcal{H}}_1$ and $h_2\in{\mathcal{H}}_2$ it holds that $\xi_Y := h_1({{\varepsilon}}_H, {\varepsilon}_Y)$ and $\xi_X := h_2({{\varepsilon}}_H, {\varepsilon}_X)$ have mean zero.[^2] Let $\mathcal{M}\coloneqq
{\mathcal{F}}\times{\mathcal{G}}\times{\mathcal{H}}_1\times{\mathcal{H}}_2\times\mathcal{Q}$ denote the model class. Every model $M=(f, g, h_1, h_2, Q)\in \cM$ then specifies an SCM by[^3]

$$\begin{aligned}
    A &\coloneqq {{\varepsilon}}_A \; &{\text{\scriptsize $r$  assignments} }   \\
    H &\coloneqq {{\varepsilon}}_H \; &{\text{\scriptsize $q$  assignments} }   \\
    X &\coloneqq g(A) + h_2(H, {\varepsilon}_X)\; &{\text{\scriptsize $d$  assignments} }\\
    Y &\coloneqq f(X) + h_1(H, {\varepsilon}_Y)\; &{\text{\scriptsize $1$  assignment} }
  \end{aligned}$$

with $({\varepsilon}_X,{\varepsilon}_Y, {{\varepsilon}}_A, {{\varepsilon}}_H)\sim Q$. For each model $M = (f,g,h_1,h_2,Q) \in\mathcal{M}$, we refer to $f$ as the *causal function* (for the pair $(X,Y)$) and assume that the entailed distribution has finite second moments. Furthermore, we denote by ${\mathbb{P}}_M$ the joint distribution over the observed variables $(X, Y, A)$ induced by the SCM specified by $M$. If no exogenous variables $A$ exist, one can think of the function $g$ as being a constant function. 0 We now specify the statistical model, we consider in this work. For a real-valued response variable $Y\in{\mathbb{R}}$ and predictors $X \in {\mathbb{R}}^d$, we consider the problem of estimating a regression function that works well not only on the training data, but also under interventions. In order to formally model interventions we will introduce a class of structural causal models (SCMs), over the response $Y \in {\mathbb{R}}^1$, the covariates $X\in {\mathbb{R}}^d$, exogenous variables $A\in {\mathbb{R}}^r$, and hidden variables $H \in {\mathbb{R}}^q$.

Let ${\mathcal{F}}\subseteq\{f:{\mathbb{R}}^d\rightarrow{\mathbb{R}}\}$, ${\mathcal{G}}\subseteq\{g:{\mathbb{R}}^r\rightarrow{\mathbb{R}}^d\}$, ${\mathcal{H}}_1\subseteq\{h_1:{\mathbb{R}}^{q+1}\rightarrow{\mathbb{R}}\}$ and ${\mathcal{H}}_2\subseteq\{h_2:{\mathbb{R}}^{q+d}\rightarrow{\mathbb{R}}^d\}$ be fixed sets of measurable functions. Moreover, let $\mathcal{Q}$ be a collection of probability distributions on ${\mathbb{R}}^{d+1+r+q}$ such that for all $Q\in\mathcal{Q}$ it holds if $({\varepsilon}_X,{\varepsilon}_Y, {{\varepsilon}}_A, {{\varepsilon}}_H)\sim Q$ then ${\varepsilon}_X,{\varepsilon}_Y, {{\varepsilon}}_A$ and ${{\varepsilon}}_H$ are jointly independent and for all $h_1\in{\mathcal{H}}_1$ and $h_2\in{\mathcal{H}}_2$ it holds that $\xi_Y := h_1({{\varepsilon}}_H, {\varepsilon}_Y)$ and $\xi_X := h_2({{\varepsilon}}_H, {\varepsilon}_X)$ have mean zero.[^4] Let $\mathcal{M}\coloneqq
{\mathcal{F}}\times{\mathcal{G}}\times{\mathcal{H}}_1\times{\mathcal{H}}_2\times\mathcal{Q}$ denote the model class, then every model, $M=(f, g, h_1, h_2, Q)\in \cM$, specifies an SCM by[^5]

$$\begin{aligned}
A &\coloneqq {{\varepsilon}}_A \; &{\scriptsize \text{$r$  assignments} }   \\
    H &\coloneqq {{\varepsilon}}_H \; &{\scriptsize \text{$q$  assignments} }   \\
    X &\coloneqq g(A) + h_2(H, {\varepsilon}_X)\; &{\scriptsize \text{$d$  assignments} }\\
    Y &\coloneqq f(X) + h_1(H, {\varepsilon}_Y)\; &{\scriptsize \text{$1$  assignment} }
    \end{aligned}$$

with $({\varepsilon}_X,{\varepsilon}_Y, {{\varepsilon}}_A, {{\varepsilon}}_H)\sim Q$.

We require that these SCMs induce a unique distribution over the variables $(X, Y, A, H)$. For each model $M = (f,g,h_1,h_2,Q) \in\mathcal{M}$, we refer to $f$ as the *causal function* (for the pair $(X,Y)$). Furthermore, we denote by ${\mathbb{P}}_M$ the joint distribution over the observed variables $(X, Y, A)$ induced by the SCM specified by $M$. Each SCM $M\in\mathcal{M}$ can now be modified by the concept of interventions [e.g., @pearl2009causality; @Peters2017book]. An intervention corresponds to replacing one or more of the structural assignments of the SCM. Thus, for each $M\in \cM$ and intervention $i$ we consider the intervened SCM $M(i)$ over the variables $(X^i, Y^i, A^i, H^i)$ given by the SCM of $M$ under intervention $i$. We will henceforth refer to $M(i)$ as an intervened model. We do not require that there exists an $M'\in \cM$ such that the SCM of $M'$ coincides with $M(i)$. We only consider interventions, such that $M(i)$ induces a unique distribution over $(X^i,Y^i,A^i,H^i)$. Let $\cI$ denote a subset of such interventions.

Interventions {#sec:interventions}
-------------

Each SCM $M\in\mathcal{M}$ can now be modified by the concept of interventions [e.g., @pearl2009causality; @Peters2017book]. An intervention corresponds to replacing one or more of the structural assignments of the SCM. For example, we intervene on all covariates $X$ by replacing the $d$ assignments with, e.g., a random variable, which is independent of the other noise variables and has a multivariate Gaussian distribution. Importantly, an intervention on some of the variables does not change the assignment of any other variable. In particular, an intervention on $X$ does not change the conditional distribution of $Y$, given $X$ and $H$ (this is an instance of the invariance property mentioned in Section \[sec:intro\]). More generally, we denote by $M(i)$ the intervened SCM over the variables $(X^i, A^i, Y^i, H^i)$, obtained by performing the intervention $i$ in model $M$. We do not require that the intervened model $M(i)$ belong to the model class $\cM$, but we require that $M(i)$ induces a joint distribution over $(X^i,Y^i,A^i,H^i)$, which has finite second moments. We use ${\mathcal{I}}$ to denote a collection of interventions.

In this work, we only consider interventions on the covariates $X$ and $A$. More specifically, for a given model $M=(f,g,h_1,h_2,Q)\in \cM$ and an intervention $i\in{\mathcal{I}}$, the intervened SCM $M(i)$ takes one of two forms. First, for an intervention on $X$ it is given by $$\begin{aligned}
  A^i := {{\varepsilon}}_A^i, \quad H^i := {{\varepsilon}}_H^i, \quad X^i :=
  \psi^i(g, h_2, A^i,H^i, {{\varepsilon}}^i_X ,I^i), \quad Y^i:= f(X^i) + h_1(H^i,{{\varepsilon}}_Y^i),\end{aligned}$$ and, second, for an intervention on $A$ it is given by $$A^i := \psi^i(I^i, {{\varepsilon}}_A^i), \quad H^i := {{\varepsilon}}_H^i, \quad X^i := g(A^i) + h_2(H^i, {{\varepsilon}}_X^i), \quad Y^i:= f(X^i) + h_1(H^i,{{\varepsilon}}_Y^i).$$ In both cases, $({{\varepsilon}}_X^i, {{\varepsilon}}_Y^i, {{\varepsilon}}_A^i,{{\varepsilon}}_H^i) \sim Q$, the (possibly degenerate) random vector $I^i$ is independent of $({{\varepsilon}}_X^i, {{\varepsilon}}_Y^i, {{\varepsilon}}_A^i,{{\varepsilon}}_H^i)$, and $\psi^i$ is a measurable function, whose arguments are all part of the structural assignment of the intervened variable in model $M$. We will see below that this class of interventions is rather flexible. It does, however, not allow for arbitrary manipulations of $M$. For example, the noise variable ${{\varepsilon}}_Y$ is not allowed to enter the structural assignment of the intervened variable. Interventions on $A$ will generally be easier to analyze than interventions on $X$. We therefore distinguish between the following different types of interventions on $X$. Let $i$ be an intervention on $X$ with intervention map $\psi^i$. The intervention is then called $$\begin{array}{ll}
    &\quad\textit{confounding-preserving}\qquad
    \begin{array}{ll}
      \text{if
      there exists a map ${\varphi}^i$, such that}\\
      \psi^i(g, h_2, A^i ,H^i, {{\varepsilon}}^i_X ,I^i) = {\varphi}^i(A^i,g(A^i),h_2(H^i,
      {{\varepsilon}}^i_X) ,I^i)
    \end{array}\\
    &\text{and it is called}\\
    &\quad\textit{confounding-removing}\phantom{ii}\qquad
    \begin{array}{ll}
      \text{if for all models $M \in {\mathcal{M}}$,}\\
      \psi^i(g, h_2, A^i ,H^i, {{\varepsilon}}^i_X ,I^i) {\perp \!\!\! \perp}H^i\quad\text{under }M(i).
    \end{array}
  \end{array}$$ Furthermore, we call a set of interventions ${\mathcal{I}}$ *well-behaved* either if it consists only of confounding-preserving interventions or if it contains at least one confounding-removing intervention. Confounding-preserving interventions contain, e.g., *shift interventions* on $X$, which linearly shift the original assignment by $I^i$, that is, $\psi^i(g, h_2, A^i,H^i, {{\varepsilon}}^i_X ,I^i) = g(A^i) + h_2(H^i, {{\varepsilon}}_X^i) +
I^i$. The name ‘confounding-preserving’ stems from the fact that the unobserved (confounding) variables $H$ only enter the intervened structural assignment of $X$ via the term $h_2(H^i, {{\varepsilon}}^i_X)$, which is the same as in the original model. Some interventions are confounding-removing and confounding-preserving, but not every confounding-removing intervention is confounding-preserving. For example, the intervention $\psi^i(g, h_2, A^i ,H^i, {{\varepsilon}}^i_X ,I^i)={{\varepsilon}}^i_X$ is confounding-removing but, in general, not confounding-preserving. Similarly, not all confounding-preserving are confounding-removing.

If the context does not allow for any ambiguity, we omit the superscript $i$ and write expressions such as ${\mathbb{E}}_{M(i)}[(Y-f(X))^2]$. The support of random variables under interventions will play an important role for the analysis of distribution generalization. Throughout this paper, ${\mathrm{supp}}^{M}(Z)$ denotes the support of the random variable $Z\in\{A, X, H, Y\}$ under the distribution ${\mathbb{P}}_M$, which is induced by the SCM $M\in\mathcal{M}$. Moreover, ${\mathrm{supp}}_{{\mathcal{I}}}^{M}(Z)$ denotes the union of ${\mathrm{supp}}^{M(i)}(Z)$ over all interventions $i\in{\mathcal{I}}$. We call a collection of interventions on $Z$ *support-reducing* (w.r.t. $M$) if ${\mathrm{supp}}_\cI^M(Z)\subseteq {\mathrm{supp}}^M(Z)$ and *support-extending* (w.r.t. $M$) if ${\mathrm{supp}}_\cI^M(Z)\not \subseteq {\mathrm{supp}}^M(Z)$. Whenever it is clear from the context which model is considered, we may drop the indication of $M$ altogether and simply write ${\mathrm{supp}}(Z)$.

Interventional robustness and the causal function {#sec:robustness}
=================================================

Let ${\mathcal{M}}$ be a fixed model class, let $M = (f,g,h_1, h_2, Q) \in {\mathcal{M}}$ be the true data generating model, and let ${\mathcal{I}}$ be a class of interventions. In this work, we aim to find a function $f^*:{\mathbb{R}}^d\rightarrow{\mathbb{R}}$, such that the predictive model $\hat Y = f^*(X)$ has low worst-case risk over all distributions induced by the interventions in $\cI$. We therefore consider the optimization problem $$\label{eq:minimax_problem}
\operatorname*{argmin}_{{f_{\diamond}}\in {\mathcal{F}}} \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{M(i)} \big[ (Y - {f_{\diamond}}(X))^2 \big],$$ where ${\mathbb{E}}_{M(i)}$ is the expectation in the intervened model $M(i)$. In general, this optimization problem is neither guaranteed to have a solution, nor is the solution, if it exists, ensured to be unique. Whenever a solution $f^*$ exists, we refer to it as a *minimax solution* (for model $M$ w.r.t. ($\mathcal{F},{\mathcal{I}}$)).

If, for example, ${\mathcal{I}}$ consists only of the trivial intervention, that is, ${\mathbb{P}}_M = {\mathbb{P}}_{M(i)}$, we are looking for the best predictor on the observational distribution. In that case, the minimax solution is obtained by any conditional mean function, $f^*: x\mapsto{\mathbb{E}}[Y\vert X=x]$ (provided that $f^*\in\mathcal{F}$). For larger classes of interventions, however, the conditional mean may become sub-optimal in terms of prediction. To see this, it is instructive to decompose the risk under an intervention. Since the structural assignment for $Y$ remains unchanged for all interventions that we consider in this work, it holds for all ${f_{\diamond}}\in {\mathcal{F}}$ and all interventions $i$ on either $A$ or $X$ that $${\mathbb{E}}_{M(i)}[(Y - {f_{\diamond}}(X))^2]
    = {\mathbb{E}}_{M(i)}[(f(X) - {f_{\diamond}}(X))^2]+{\mathbb{E}}_{M}[\xi_Y^2]+2{\mathbb{E}}_{M(i)}[\xi_Y(f(X)-{f_{\diamond}}(X))].$$ Here, the middle term does not depend on $i$ since $\xi_Y = h_1(H, {{\varepsilon}}_Y)$ remains fixed. If $i$ is a confounding-removing intervention, then $\xi_Y{\perp \!\!\! \perp}X$ under ${\mathbb{P}}_{M(i)}$, and, because of ${\mathbb{E}}_{M}[\xi_Y] = 0$, the last term in the above equation vanishes. Therefore, if ${\mathcal{I}}$ consists only of confounding-removing interventions, the causal function is a solution to the minimax- problem . The following proposition shows that an even stronger statement holds: The causal function is already a minimax solution if ${\mathcal{I}}$ contains at least one confounding-removing intervention on $X$.

\[prop:minimax\_equal\_causal\] If ${\mathcal{I}}$ is a set of interventions on $X$ or $A$ and at least one of these is a confounding-removing intervention, then the causal function $f$ is a minimax solution.

One step in the proof of this proposition is to show that the minimal worst-case loss is attained at a confounding-removing intervention. That is, $$\inf_{{f_{\diamond}}\in{\mathcal{F}}}\sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{M(i)} \big[ (Y - {f_{\diamond}}(X))^2 \big]=\inf_{{f_{\diamond}}\in{\mathcal{F}}}\sup_{i \in {\mathcal{I}}_{\text{cr}}} {\mathbb{E}}_{M(i)} \big[ (Y - {f_{\diamond}}(X))^2 \big],$$ where ${\mathcal{I}}_{\text{cr}}{\subseteq}{\mathcal{I}}$ denotes the non-empty subset of confounding-removing interventions. This observation will also be used in the proofs of some of the results that follow below.

We now prove that when restricting ourselves to linear functions only, the causal function is also a minimax solution with respect to the set of all shift interventions on $X$ – interventions that appear in linear IV models and recently gained further attention in the causal community [@rothenhausler2018anchor; @Sani2020]. The proposition below also makes precise in which sense shift interventions are related to linear model classes. Intuitively, when the causal relation between $X$ and $Y$ is linear, shift interventions are sufficient to create unbounded variability in all directions of the covariance matrix of $X$ (more precisely, the unbounded eigenvalue condition below is satisfied if ${\mathcal{I}}$ is the set of all shift interventions on $X$). As the following proposition shows, under this condition, the causal function is a minimax solution.

\[prop:shift\_interventions\] Let ${\mathcal{F}}$ be the class of all linear functions, and let ${\mathcal{I}}$ be a set of interventions on $X$ or $A$ s.t. $\sup_{i\in \cI} \lambda_{\min}\big({\mathbb{E}}_{M(i)}\big[XX^\top\big]\big) =\infty$, where $\lambda_{\min}$ denotes the smallest eigenvalue (assuming that the considered moments exist). Then, the causal function $f$ is the unique minimax solution.

Even if the causal function $f$ does not solve the minimax problem , the difference between the minimax solution and the causal function cannot be arbitrarily large. The following proposition shows that the worst-case $L_2$-distance between $f$ and any function ${f_{\diamond}}$ that performs better than $f$ (in terms of worst-case risk) can be bounded by a term which is related to the strength of the confounding.

\[prop:difference\_to\_causal\_function\] Let ${\mathcal{I}}$ be a set of interventions on $X$ or $A$. Then, for any function ${f_{\diamond}}\in{\mathcal{F}}$ which satisfies that $$\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-{f_{\diamond}}(X))^2]\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-f(X))^2],$$ it holds that $$\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]\leq 4{\operatorname{Var}}_M(\xi_Y).$$

Even though the difference can be bounded, it may be non-zero, and one may benefit from choosing a function that differs from the causal function $f$. This choice, however, comes at a cost: it relies on the fact that we know the class of interventions $\mathcal{I}$. In general, being a minimax solution is not entirely robust with respect to misspecification of $\mathcal{I}$. In particular, if the set $\cI_2$ of interventions describing the test distributions is misspecified by a set $\cI_1\neq\cI_2$, then the considered minimax solution with respect to $\cI_1$ may perform worse than the causal function on the test distributions.

\[prop:misspecification\_minimax\] Let ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$ be any two sets of interventions on $X$, and let $f_1^*\in\mathcal{F}$ be a minimax solution w.r.t. ${\mathcal{I}}_1$. Then, if ${\mathcal{I}}_2\subseteq{\mathcal{I}}_1$ it holds that $$\sup_{i\in{\mathcal{I}}_2}{\mathbb{E}}_{M(i)}\big[(Y-f_1^*(X))^2\big] \leq \sup_{i\in{\mathcal{I}}_2}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2\big].$$ If ${\mathcal{I}}_2\not\subseteq{\mathcal{I}}_1$, however, it can happen (even if $\mathcal{F}$ is linear) that $$\sup_{i\in{\mathcal{I}}_2}{\mathbb{E}}_{M(i)}\big[(Y-f_1^*(X))^2\big] > \sup_{i\in{\mathcal{I}}_2}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2\big].$$

The second part of the proposition should be understood as a non-robustness property of non-causal minimax solutions. Improvements on the causal function are possible in situations, where one has reasons to believe that the test distributions do not stem from a set of interventions that is much larger than the specified set.

So far, the optimizer of the minimax problem  depends on the true model $M$. In practice, however, we do not have access to the true model $M$, but only to its observational distribution ${\mathbb{P}}_M$. This motivates the definition of distribution generalization.

Distribution generalization {#sec:generalizability}
===========================

Throughout this section, let ${\mathcal{M}}$ denote a fixed model class, let $M = (f,g,h_1,h_2,Q) \in {\mathcal{M}}$ be the true (but unknown) data generating model, with observational distribution ${\mathbb{P}}_M$, and let ${\mathcal{I}}$ be a set of interventions on $X$ or $A$. Depending on the model class ${\mathcal{M}}$, there may be several models $\tilde{M} \in {\mathcal{M}}$ that induce the observational distribution ${\mathbb{P}}_M$ but do not agree with $M$ on all intervention distributions induced by $\cI$. Each such model induces a potentially different minimax problem, with a potentially different set of solutions. Given knowledge only of ${\mathbb{P}}_M$, it is therefore generally not possible to identify a solution to . In this section, we study conditions on ${\mathcal{M}}$, ${\mathbb{P}}_M$ and ${\mathcal{I}}$, under which this becomes possible. More precisely, we aim to characterize under which conditions $({\mathbb{P}}_M, \mathcal{M})$ generalizes to $\mathcal{I}$.

\[defi:general\] $({\mathbb{P}}_{M},\mathcal{M})$ is said to generalize to ${\mathcal{I}}$ if for every ${\varepsilon}> 0$ there exists a function $f^{*}\in\mathcal{F}$ such that, for all models $\tilde{M}\in\mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}}={\mathbb{P}}_{M}$, it holds that $$\label{eq:def_generalization}
    \left \vert \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
      \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \right \vert \leq {\varepsilon}.$$

Distribution generalization does not require the existence of a minimax solution in ${\mathcal{F}}$ (which would require further assumptions on the function class ${\mathcal{F}}$) and instead focuses on whether an approximate solution can be identified based only on the observational distribution ${\mathbb{P}}_{M}$. If, however, there exists a function $f^* \in {\mathcal{F}}$ which, for every $\tilde{M}\in{\mathcal{M}}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, is a minimax solution for $\tilde{M}$ w.r.t. $({\mathcal{F}}, {\mathcal{I}})$, then, in particular, $({\mathbb{P}}_{M},\mathcal{M})$ generalizes to ${\mathcal{I}}$. As the next proposition shows, generalization is closely linked to the ability of identifying the joint intervention distributions of $(X, Y)$ from the observational distribution.

0 [[Rune: maybe introduce an “${\varepsilon}$-minimax solution” in Section \[sec:robustness\]? That would make the discussion below a bit easier (and the above definition more compact)]{}]{}.

Definition \[defi:general\] says that the observational distribution ${\mathbb{P}}_{M}$ is sufficient for identifying a solution $f^*$ that for all $\tilde{M}\in\mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}}={\mathbb{P}}_{M}$ solves (up to ${\varepsilon}$) the minimax problem corresponding to $\tilde{M}$. The ability to generalize, in particular, ensures that there exists a well-defined mapping [[Jonas: does this not require that every ${\mathbb{P}}_M$ generalizes?]{}]{} $$F:\{{\mathbb{P}}_{\tilde{M}}\,\vert\,\tilde{M}\in\mathcal{M}\}\to {\mathcal{F}}$$ such that $F({\mathbb{P}})$ solves the minimax problem (up to ${\varepsilon}$) for all observational equivalent models (i.e., all models $\tilde{M}\in\mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}}={\mathbb{P}}$). That is, the observational distribution of the variables $(X, Y, A)$ is sufficient to identify a solution that minimizes the worst-case intervention mean squared error under all models with the same observational distribution.

\[prop:suff\_general\] Assume that for all $\tilde{M} \in {\mathcal{M}}$ it holds that[^6] $${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M \quad \Rightarrow \quad {\mathbb{P}}^{(X,Y)}_{\tilde{M}(i)} = {\mathbb{P}}^{(X,Y)}_{M(i)} \quad \forall i \in {\mathcal{I}},$$ where ${\mathbb{P}}^{(X,Y)}_{M(i)}$ is the joint distribution of $(X, Y)$ under $M(i)$. Then, $({\mathbb{P}}_M, {\mathcal{M}})$ generalizes to ${\mathcal{I}}$.

Proposition \[prop:suff\_general\] provides verifiable conditions for distribution generalization, and is a useful result for proving possibility statements. It is, however, not a necessary condition. In Propositions \[prop:genX\_intra\] and  \[prop:genX\_extra\], we give further conditions under which distribution generalization is possible for all well-behaved sets of interventions. In particular, if the set of interventions ${\mathcal{I}}$ contains at least one confounding-removing intervention it can be shown that the causal function always generalizes, even in cases where the interventional marginal of $X$ is not identified. We will see that distribution generalization is closely linked to the relation between the support of ${\mathbb{P}}_M$ and the support of the intervention distributions. Below, we therefore distinguish between support-reducing interventions (Section \[sec:supp\_reducing\_onX\]) and support-extending interventions (Section \[sec:support\_extending\_onX\]) on $X$. In Section \[sec:int\_onA\], we consider interventions on $A$. We will see that parts of the analysis carry over from the interventions on $X$.

Support-reducing interventions on $X$ {#sec:supp_reducing_onX}
-------------------------------------

In order to simplify the following analysis, we will constrain ourselves to cases in which the causal function is identified on the support of $X$. This condition is made precise in the following assumption.

\[ass:identify\_f\] For all $\tilde{M}=(\tilde{f}, \dots) \in \mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_{M}$, it holds that $\tilde{f}(x) = f(x)$ for all $x \in {\mathrm{supp}}(X)$.

Assumption \[ass:identify\_f\] concerns identifiability of the causal function from the observational distribution on the support of $X$. This question has received a lot of attention in literature. In linear instrumental variable settings, for example, one assumes that the functions $f$ and $g$ are linear and the product moment between $A$ and $X$ has rank at least the dimension of $X$ [e.g., @wooldridge2010econometric]. In linear non-Gaussian models, one can identify the function $f$ even if there are no instruments [@Hoyer2008b]. For nonlinear models, restricted structural causal models can be exploited, too. In that case, Assumption \[ass:identify\_f\] holds under regularity conditions if $h_1(H, {\varepsilon}_Y)$ is independent of $X$ [@Zhang2009; @Peters2014jmlr; @Peters2017book] and first attempts have been made to extend such results to non-trivial confounding cases [@Janzing2009uai]. The nonlinear IV setting [e.g., @amemiya1974nonlinear; @newey2013nonparametric; @newey2003instrumental] is discussed in more detail in Appendix \[sec:IVconditions\], where we give a brief overview of identification results for linear, parametric and non-parametric function classes. There is also a technical aspect regarding identifiability: Assumption \[ass:identify\_f\] states that $f$ is identifiable, even on ${\mathbb{P}}_M$-null sets, which is usually achieved by placing further constraints on the function class, such as smoothness. Even though this issue seems technical, it becomes important when considering hard interventions that set $X$ to a fixed value, for example.

Assumption \[ass:identify\_f\] is not necessary for generalization. [@rothenhausler2018anchor] show, for example, that if ${\mathcal{F}}$ and ${\mathcal{G}}$ consist of linear functions it is possible to generalize to a set of bounded interventions on $A$ – even if Assumption \[ass:identify\_f\] does not hold. If, however, Assumption \[ass:identify\_f\] holds, then distribution generalization is possible even in nonlinear settings, under a large class of interventions if these are support-reducing.

\[prop:genX\_intra\] Let ${\mathcal{I}}$ be a well-behaved set of interventions on $X$, and assume that ${\mathrm{supp}}_{{\mathcal{I}}}(X)\subseteq{\mathrm{supp}}(X)$. Then, under Assumption \[ass:identify\_f\], $({\mathbb{P}}_{M}, \cM)$ generalizes to the interventions ${\mathcal{I}}$.

Proposition \[prop:genX\_intra\] states that Assumption \[ass:identify\_f\] is a sufficient condition for generalization when ${\mathcal{I}}$ is a well-behaved set of support-reducing interventions. However, for an arbitrary set of interventions, generalization can become impossible, even if Assumption \[ass:identify\_f\] is satisfied and all interventions are support-reducing.

### Impossibility of generalization under changes in confounding

0

Remove: [[todo: explain example]{}]{} [[todo: turn this into a proposition]{}]{} [[Jonas: Assumption \[ass:identify\_f\] is not satisfied, correct? Isn’t it possible though to adapt the example?]{}]{}

\[ex:impossible\_interpolation\] For $c\in(0,1)$, consider the following models $$M_c:\quad
    \begin{cases}
      &A\coloneqq {\varepsilon}_A\quad H\coloneqq{\varepsilon}_H\\
      &X\coloneqq \sqrt{1-\big(\tfrac{1}{c}\big)^2}{\varepsilon}_X +
      \tfrac{1}{c}H\\
      &Y\coloneqq \beta X + \sqrt{1-c^2}{\varepsilon}_Y + c H,
      \quad\text{with }
{{\color{red}Jonas: {\varepsilon}_A}},      
      {\varepsilon}_X,{\varepsilon}_Y,{\varepsilon}_H \overset{i.i.d.}{\sim}\mathcal{N}(0,1)
    \end{cases}$$ Then, it follows that the joint distribution $\bP_{M_c}$ of $(X,Y,A)$ is Gaussian and independent of $c$ with ${\operatorname{Var}}(X)=1$, ${\operatorname{Var}}(Y)=1+\beta^2+2\beta$ and ${\operatorname{Cov}}(X,Y)=\beta+1$. For every $k\in(0, \infty)$ [[Jonas: can we not simply use $i\in(0, \infty)$? this would save us a subscript]{}]{} define the intervention $i_k$ on $X$ by the mapping $\psi_k(H)=k H$. [[Jonas: $\psi^{i_k}(g, h_2, A^i,H^i, {{\varepsilon}}^i_X ,I^i) = kH$.]{}]{} These interventions are not confounding-preserving [[Jonas: well-behaved?]{}]{}. Using the properties of the ordinary least square estimator we get that $$\operatorname*{argmin}_{\beta}{\mathbb{E}}_{M_c(i_k)}\big[(Y-\beta
    X)^2\big]=\beta^{OLS}_{M_c(i_k)}=\frac{{\operatorname{Cov}}(Y, X)}{{\operatorname{Var}}{X}}
{{\color{red}Jonas: =\frac{{\operatorname{Cov}}_{M_c(i_k)}(Y, X)}{{\operatorname{Var}}_{M_c(i_k)}{X}}}}    
    =\beta+\frac{c}{k}.$$ This implies that for ${\mathcal{I}}\coloneqq\{i_k\}$ the minimax solution can range in the interval $(\beta, \beta+\frac{1}{k})$ depending on the model, i.e., the observational distribution of $(X, Y)$ does not generalize in this setting, i.e., for any [[Jonas: why ‘any’?]{}]{} $\{M_c:c\in(0,1)\}\subseteq\mathcal{M}$ it holds [[Jonas: for any $c\in(0,1)$]{}]{} that $(\bP_{M_{c}},\mathcal{M})$ does not generalize to the intervention $\cI=\{i_k\}$ [[Remove: for any $c\in(0,1)$]{}]{}.

Consider, for example, a one-dimensional linear instrumental variable setting. Let therefore $\mathcal{Q}$ be a class of product distributions on ${\mathbb{R}}^4$, such that for all $Q \in \mathcal{Q}$, the coordinates of $Q$ are non-degenerate, zero-mean with finite second moment. Let ${\mathcal{M}}$ be the class of all models of the form $$\label{eq:modelmc}
A\coloneqq {\varepsilon}_A, \quad 
H\coloneqq\sigma{\varepsilon}_H, \quad 
X\coloneqq \gamma A + {\varepsilon}_X + \tfrac{1}{\sigma}H, \quad
Y\coloneqq \beta X + {\varepsilon}_Y + \tfrac{1}{\sigma}H,$$ with $\gamma, \beta \in {\mathbb{R}}$, $\sigma > 0$ and $({\varepsilon}_A, {\varepsilon}_X,{\varepsilon}_Y,{\varepsilon}_H) \sim Q \in \mathcal{Q}$. Assume that ${\mathbb{P}}_M$ is induced by some (unknown) model $M = M(\gamma, \beta, \sigma, Q)$ from the above model class (here, we slightly adapt the notation from Section \[sec:setup\]). The following proposition shows that if the set of interventions ${\mathcal{I}}$ is not well-behaved, distribution generalization is not always ensured.

\[prop:impossibility\_interpolation\] Assume that ${\mathcal{M}}$ is given as defined above, and let ${\mathcal{I}}\subseteq {\mathbb{R}}_{>0}$ be a compact set of interventions on $X$ defined by $\psi^i(g, h_2, A^i, H^i, {\varepsilon}_X^i, I^i) = iH$, for $i \in {\mathcal{I}}$ (this set of interventions is not well-behaved). Then, $({\mathbb{P}}_M, \cM)$ does not generalize to the interventions in ${\mathcal{I}}$ (even if Assumption \[ass:identify\_f\] is satisfied). In addition, any prediction model other than the causal model may perform arbitrarily bad under the interventions ${\mathcal{I}}$. That is, for any $b \neq \beta$ and any $c > 0$, there exists a model $\tilde{M}\in\cM$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_{M}$, such that $$\abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-bX)^2\big]
      - \inf_{{b_{\diamond}}\in {\mathbb{R}}}\,
      \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{b_{\diamond}}X)^2\big]} \geq c.$$

Support-extending interventions on $X$ {#sec:support_extending_onX}
--------------------------------------

If the interventions in ${\mathcal{I}}$ extend the support of $X$, i.e., ${\mathrm{supp}}_{\mathcal{I}}(X) \not {\subseteq}{\mathrm{supp}}(X)$, Assumption \[ass:identify\_f\] is not sufficient for ensuring distribution generalization. This is because there may exist models $\tilde{M} \in {\mathcal{M}}$ which agree with $M$ on the observational distribution, but whose corresponding causal function $\tilde{f}$ differs from $f$ outside of the support of $X$ In that case, a support-extending intervention on $X$ may result in different dependencies between $X$ and $Y$ in the two models, and therefore induce a different set of minimax solutions. The following assumption on the model class ${\mathcal{F}}$ ensures that any $f\in \cF$ is uniquely determined by its values on ${\mathrm{supp}}(X)$.

\[ass:gen\_f\] For all $\tilde{f}, \bar{f} \in {\mathcal{F}}$ with $\tilde{f}(x) = \bar{f}(x)$ for all $x \in {\mathrm{supp}}(X)$, it holds that $\tilde{f} \equiv \bar{f}$.

We will see that this assumption is sufficient (Proposition \[prop:genX\_extra\]) for generalization with respect to well-behaved interventions on $X$. Furthermore, it is also necessary (Proposition \[prop:impossibility\_extrapolation\]) if $\cF$ is sufficiently flexible. The following proposition can be seen as an extension of Proposition \[prop:genX\_intra\].

\[prop:genX\_extra\] Let ${\mathcal{I}}$ be a well-behaved set of interventions on $X$. Then, under Assumptions \[ass:identify\_f\] and \[ass:gen\_f\], $({\mathbb{P}}_{M}, \cM)$ generalizes to ${\mathcal{I}}$.

Because the interventions may change the marginal distribution of $X$, the preceding proposition includes examples, in which distribution generalization is possible even if some of the considered joint (test) distributions are arbitrarily far from the training distribution, in terms of any reasonable divergence measure over distributions, such as Wasserstein distance or $f$-divergence.

The proposition relies on Assumption \[ass:gen\_f\]. Even though this assumption is restrictive, it is satisifed by several reasonable function classes, which therefore allow for generalization under any set of well-behaved interventions. Below, we give two examples of such a function class.

### Sufficient conditions for generalization {#sec:suff_genX_extrap}

Assumption \[ass:gen\_f\] states that every function in ${\mathcal{F}}$ is globally identified by its values on ${\mathrm{supp}}(X)$. This is, for example, satisfied if $\mathcal{F}$ is a linear space of functions with domain $\mathcal{D} {\subseteq}{\mathbb{R}}^d$ which are linearly independent on ${\mathrm{supp}}(X)$. More precisely, $$\begin{aligned}
\mathcal{F} \text{ is linearly closed}: &\quad 
f_1, f_2 \in \mathcal{F}, c \in {\mathbb{R}}, \implies
f_1 + f_2 \in \mathcal{F}, cf_1 \in \mathcal{F}, \text{ and} \label{eq:linear}\\
\mathcal{F} \text{ is lin.\ ind.\ on }{\mathrm{supp}}(X): 
    &\quad f_1(x) = 0 \quad \forall x \in {\mathrm{supp}}(X) \;\implies \;f_1(x) = 0 \quad \forall x \in \mathcal{D}. \label{eq:linearind}\end{aligned}$$ Examples of such classes include (i) globally linear parametric function classes, i.e., ${\mathcal{F}}$ is of the form $$\mathcal{F}^1 \coloneqq\{{f_{\diamond}}:\mathcal{D}\rightarrow{\mathbb{R}}{\, \vert \,}\text{ there exists } \gamma \in {\mathbb{R}}^k \text{ s.t. } 
  \forall
  x\in\mathcal{D} {\, : \,}{f_{\diamond}}(x)=\gamma^\top \nu (x) \},$$ where $\nu = (\nu_1, \dots, \nu_k)$ consists of real-valued, linearly independent functions satisfying that ${\mathbb{E}}_M[\nu(X) \nu(X)^\top]$ is strictly positive definite, and (ii) the class of differentiable functions that extend linearly outside of ${\mathrm{supp}}(X)$, that is, ${\mathcal{F}}$ is of the form $$\mathcal{F}^2 \coloneqq\{{f_{\diamond}}:\mathcal{D}\rightarrow{\mathbb{R}}\,\vert\,
  {f_{\diamond}}\in C^1 \text{ and } \forall x\in\mathcal{D} \setminus
  {\mathrm{supp}}(X):\, {f_{\diamond}}(x)={f_{\diamond}}(x_b)+\nabla{f_{\diamond}}(x_b)(x-x_b)\},
  $$ where $x_b\coloneqq\operatorname*{argmin}_{z\in{\mathrm{supp}}(X)}\norm{x-z}$ and ${\mathrm{supp}}(X)$ is assumed to be closed with non-empty interior. Clearly, both of the above function classes are linearly closed. To see that ${\mathcal{F}}^1$ satisfies , let $\gamma \in {\mathbb{R}}^k$ be s.t. $\gamma^\top \nu (x) = 0$ for all $x \in {\mathrm{supp}}(X)$. Then, it follows that $0 = {\mathbb{E}}_M[(\gamma^\top \nu(X))^2] = \gamma^\top {\mathbb{E}}_M[\nu(X)
\nu(X)^\top] \gamma$ and hence that $\gamma = 0$. To see that ${\mathcal{F}}^2$ satisfies , let ${f_{\diamond}}\in {\mathcal{F}}^2$ and assume that ${f_{\diamond}}(x) = 0$ for all $x \in {\mathrm{supp}}(X)$. Then, ${f_{\diamond}}(x) = 0$ for all $x \in \mathcal{D}$ and thus ${\mathcal{F}}^2$ uniquely defines the function on the entire domain $\mathcal{D}$. By Proposition \[prop:genX\_extra\], generalization with respect to these model classes is possible for any well-behaved set of interventions. In practice, it may often be more realistic to impose bounds on the higher order derivatives of the functions in ${\mathcal{F}}$. We now prove that this still allows for approximate distribution generalization, see Propositions \[prop:extrapolation\_bounded\_deriv\_cr\] and \[prop:extrapolation\_bounded\_deriv\].

### Sufficient conditions for approximate distribution generalization {#sec:suffapprox}

For differentiable functions, exact generalization cannot always be achieved. Bounding the first derivative, however, allows us to achieve approximate generalization. We therefore consider the following function class $$\label{eq:boundedder}
    \mathcal{F}^3\coloneqq\{{f_{\diamond}}:\mathcal{D}\rightarrow{\mathbb{R}}\,\vert\,
    {f_{\diamond}}\text{ is continuously differentiable with } \norm{\nabla{f_{\diamond}}}_{\infty}\leq K\},$$ for some fixed $K<\infty$, where $\nabla{f_{\diamond}}$ denotes the gradient and $\mathcal{D}\subseteq{\mathbb{R}}^d$. We then have the following result.

\[prop:extrapolation\_bounded\_deriv\_cr\] Let $\mathcal{F}$ be as defined in . Let ${\mathcal{I}}$ be a set of interventions on $X$ containing at least one confounding-removing intervention, and assume that Assumption \[ass:identify\_f\] holds true. (In this case, the causal function $f$ is a minimax solution.) Then, for all $f^*$ with $f^*=f$ on ${\mathrm{supp}}(X)$ and all $\tilde{M}\in\mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}}={\mathbb{P}}_{M}$, it holds that $$\label{eq:gen_cond}
    \abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
      \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]}\leq
    4\delta^2K^2+4\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)},$$ where $\delta :=
  \sup_{x\in{\mathrm{supp}}_{\cI}^{M}(X)}\inf_{z\in{\mathrm{supp}}^{M}(X)}\norm{x-z}$. If ${\mathcal{I}}$ consists only of confounding-removing interventions, the same statement holds when replacing the bound by $4\delta^2K^2$.

Proposition \[prop:extrapolation\_bounded\_deriv\_cr\] states that the deviation of the worst-case generalization error from the best possible value is bounded by a term that grows with the square of $\delta$. Intuitively, this means that under the function class defined in , approximate generalization is reasonable only for interventions that are close to the support of $X$. We now prove a similar result for cases in which the minimax solution is not necessarily the causal function. The following proposition bounds the worst-case generalization error for arbitrary confounding-preserving interventions. Here, the bound additionally accounts for the approximation to the minimax solution.

\[prop:extrapolation\_bounded\_deriv\] Let $\mathcal{F}$ be as defined in . Let ${\mathcal{I}}$ be a set of confounding-preserving interventions on $X$, and assume that Assumption \[ass:identify\_f\] is satisfied. Let ${\varepsilon}> 0$ and let $f^{*}\in\mathcal{F}$ be such that $$\abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
      \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]}\leq {\varepsilon}.$$ Then, for all $\tilde{M}\in\mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}}={\mathbb{P}}_{M}$, it holds that $$\begin{aligned}
    &\abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big]
      - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
      \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]} \\
        &\quad 
         \leq {{\varepsilon}}+ 12 \delta^2 K^2 + 32 \delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} +
            4 \sqrt{2} \delta K \sqrt{{{\varepsilon}}}
    \end{aligned}$$ where $\delta :=
  \sup_{x\in{\mathrm{supp}}_{\cI}^{M}(X)}\inf_{z\in{\mathrm{supp}}^{M}(X)}\norm{x-z}$.

We can take $f^*$ to be the minimax solution if it exists. In that case, the terms involving ${{\varepsilon}}$ disappear from the bound, which then becomes more similar to the one in Proposition \[prop:extrapolation\_bounded\_deriv\_cr\].

### Impossibility of generalization without restrictions on ${\mathcal{F}}$

If we do not constrain the function class ${\mathcal{F}}$, generalization is impossible. Even if we consider the set of all continuous functions ${\mathcal{F}}$, we cannot generalize to interventions outside the support of $X$. This statement holds even if Assumption \[ass:identify\_f\] is satisfied.

\[prop:impossibility\_extrapolation\] Assume that ${\mathcal{F}}= \{{f_{\diamond}}: {\mathbb{R}}^d \to {\mathbb{R}}\mid {f_{\diamond}}\text{ is continuous}\}$. Let ${\mathcal{I}}$ be a well-behaved set of support-extending interventions on $X$, such that ${\mathrm{supp}}_{{\mathcal{I}}}(X) \setminus {\mathrm{supp}}(X)$ has non-empty interior. Then, $({\mathbb{P}}_{M}, \cM)$ does not generalize to the interventions in ${\mathcal{I}}$, even if Assumption \[ass:identify\_f\] is satisfied. In particular, for any function $\bar{f} \in{\mathcal{F}}$ and any $c > 0$, there exists a model $\tilde{M} \in {\mathcal{M}}$, with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_{M}$, such that $$\abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-\bar{f}(X))^2\big]
  - \inf_{{f_{\diamond}}\in {\mathcal{F}}}\,
  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]} \geq c.$$

Interventions on $A$ {#sec:int_onA}
--------------------

We can now derive corresponding results for interventions on $A$, for which, as we will see, parts of the analysis simplify. We will be able to employ several of the above results by realizing that any intervention on $A$ can be written as an intervention on $X$, in which the structural assignment of $X$ is altered in a way that depends on the functional relationship $g$ between $X$ and $A$. The effect of such an intervention on the prediction model is propagated by $g$. More formally, under such an intervention, a model $\tilde{M} = (\tilde{f}, \tilde{g}, \tilde{h}_1, \tilde{h}_2, \tilde{Q})$ with $\tilde{g} \neq g$ may induce a distribution over $(X,Y)$ that differs from the one induced by $M$. Without further restrictions on the function class ${\mathcal{G}}$, this may happen even in cases where $\tilde{M}$ and $M$ agree on the observational distribution. This motivates an assumption on the identifiability of $g$.

\[ass:identify\_g\] For all $\tilde{M} = (\tilde{f}, \tilde{g}, \tilde{h}_1, \tilde{h}_2, \tilde{Q}) \in {\mathcal{M}}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, it holds that $\tilde{g}(a) = g(a)$ for all $a \in {\mathrm{supp}}(A) \cup {\mathrm{supp}}_{\mathcal{I}}(A)$.

Since $g(A)$ is a conditional mean for $X$ given $A$, the values of $g$ are identified from ${\mathbb{P}}_M$ for ${\mathbb{P}}_M$-almost all $a$. If ${\mathrm{supp}}_{\mathcal{I}}(A) {\subseteq}{\mathrm{supp}}(A)$, Assumption \[ass:identify\_g\] therefore holds if, for example, ${\mathcal{G}}$ contains continuous functions only. The point-wise identifiability of $g$ is necessary, for example, if some of the test distributions are induced by hard interventions on $A$, which set $A$ to some fixed value $a \in {\mathbb{R}}^r$. In the case where the interventions ${\mathcal{I}}$ extend the support of $A$, we additionally require the function class ${\mathcal{G}}$ to extrapolate from ${\mathrm{supp}}(A)$ to ${\mathrm{supp}}(A) \cup {\mathrm{supp}}_{\mathcal{I}}(A)$; this is similar to the conditions on ${\mathcal{F}}$ which we made in Section \[sec:support\_extending\_onX\] and requires further restrictions on ${\mathcal{G}}$. Under Assumption \[ass:identify\_g\], we obtain a result corresponding to Propositions \[prop:genX\_intra\] and \[prop:genX\_extra\].

\[prop:genA\] Let ${\mathcal{I}}$ be a set of interventions on $A$ and assume Assumption \[ass:identify\_g\] is satisfied. Then, $({\mathbb{P}}_{M}, \cM)$ generalizes to ${\mathcal{I}}$ if either ${\mathrm{supp}}_{{\mathcal{I}}}(X)\subseteq {\mathrm{supp}}(X)$ and Assumption \[ass:identify\_f\] is satisfied or if both Assumptions \[ass:identify\_f\] and \[ass:gen\_f\] are satisfied.

### Impossibility of generalization without constraints on $\mathcal{G}$

Without restrictions on the model class ${\mathcal{G}}$, generalization to interventions on $A$ is impossible. This holds true even under strong assumptions on the true causal function (such as $f$ is known to be linear). Below, we give a formal impossibility result for hard interventions on $A$, which set $A$ to some fixed value, where ${\mathcal{G}}$ is the set of all continuous functions.

\[prop:impossibility\_intA\] Assume that ${\mathcal{F}}= \{{f_{\diamond}}: {\mathbb{R}}^d \to {\mathbb{R}}{\, \vert \,}{f_{\diamond}}\text{ is linear} \}$ and ${\mathcal{G}}= \{{g_{\diamond}}: {\mathbb{R}}^r \to {\mathbb{R}}^d {\, \vert \,}{g_{\diamond}}\text{ is continuous} \}$. Let $\mathcal{A} {\subseteq}{\mathbb{R}}^r$ be bounded, and let ${\mathcal{I}}$ denote the set of all hard interventions which set $A$ to some fixed value from $\mathcal{A}$. Assume that $\mathcal{A} \setminus {\mathrm{supp}}(A)$ has nonempty interior. Assume further that ${\mathbb{E}}_M[\xi_X \xi_Y] \neq 0$ (this excludes the case of no hidden confounding). Then, ${\mathbb{P}}_M$ does not generalize to the interventions in ${\mathcal{I}}$. In addition, any function other than $f$ may perform arbitrarily bad under the interventions in ${\mathcal{I}}$. That is, for any $\bar{f} \neq f$ and $c > 0$, there exists a model $\tilde{M} \in {\mathcal{M}}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$ such that $$\abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-\bar{f}(X))^2\big] - \inf_{{f_{\diamond}}\in {\mathcal{F}}}\,
  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]} \geq c.$$

0

![[[todo: ]{}]{}[]{data-label="fig:impossibility"}](impossibility_AtoX_nonlinear){width="\linewidth"}

![ The left plot illustrates the straight-forward idea behind the impossibility result in Proposition \[prop:impossibility\_extrapolation\]. The plots in the middle and on the right-hand side illustrate the impossibility result in Proposition \[prop:impossibility\_intA\]. All plots visualize the case of univariate variables. Under well-behaved interventions on $X$ (left; here using confounding-removing interventions) which extend the support of $X$, generalization is impossible without further restrictions on the function class ${\mathcal{F}}$. This holds true even if Assumption \[ass:identify\_f\] is satisfied. Indeed, although the candidate model (blue line) coincides with the causal model (green dashed curve) on the support of $X$, it may perform arbitrarily bad on test data generated under support-extending interventions. Under interventions on $A$ (right and middle), generalization is impossible even under strong assumptions on the function class ${\mathcal{F}}$ (here, ${\mathcal{F}}$ is the class of all linear functions). Any support-extending intervention on $A$ shifts the marginal distribution of $X$ by an amount which depends on the (unknown) function $g$, resulting in a distribution of $(X,Y)$ which cannot be identified from the observational distribution. Without further restrictions on the function class ${\mathcal{G}}$, any candidate model apart from the causal model may result in arbitrarily large worst-case prediction risk. []{data-label="fig:impossibility"}](impossibility_all){width="\linewidth"}

Learning generalizing models from data {#sec:learning}
======================================

So far, our focus has been on the possibility to generalize, that is, we have investigated under which conditions it is possible to identify generalizing models from the observational distribution. In practice, generalizing models need to be estimated from finitely many data. This task is challenging for several reasons. First, analytical solutions to the minimax problem are only known in few cases. Even if generalization is possible, the inferential target thus often remains a complicated object, given as a well-defined but unknown function of the observational distribution. Second, we have seen that the ability to generalize depends strongly on whether the interventions extend the support of $X$, see Propositions \[prop:genX\_extra\] and \[prop:impossibility\_extrapolation\]. In a setting with a finite amount of data, the empirical support of the data lies within some bounded region, and suitable constraints on the function class ${\mathcal{F}}$ are necessary when aiming to achieve empirical generalization outside this region, even if $X$ comes from a distribution with full support. As we show in our simulations in Section \[sec:experiments\], constraining the function class can also improve the prediction performance at the boundary of the support.

In Section \[sec:existingmethods\], we survey existing methods for learning generalizing models. Often, these methods assume either a globally linear model class ${\mathcal{F}}$ or are completely non-parametric and therefore do not generalize outside the empirical support of the data. Motivated by this observation, we introduce in Section \[sec:nile\] a novel estimator, which exploits an instrumental variable setup and a particular extrapolation assumption to learn a globally generalizing model.

Existing methods {#sec:existingmethods}
----------------

As discussed in Section \[sec:related\_work\], a wide range of methods have been proposed to guard against various types of distributional changes. Here, we review methods that fit into the causal framework in the sense that the distributions that in the minimax formulation the supremum is taken over are induced by interventions.

[&gt;p[1.8cm]{}&gt;p[3.3cm]{}&gt;p[1.9cm]{}&gt;p[1.7cm]{}|&gt;p[3.1cm]{}]{} model class & interventions & ${\mathrm{supp}}_{{\mathcal{I}}}(X)$ & assumptions & algorithm\
$\mathcal{F}$ linear & on $X$ or $A$of which at least one is confounding-removing & – & Ass. \[ass:identify\_f\] & linear IV (e.g., two-stage least squares, K-class or PULSE [@Theil1958; @jakobsen2020distributional])\
${\mathcal{F}}, {\mathcal{G}}$ linear & on A & bounded strength & – & anchor regression [@rothenhausler2018anchor], see also [@Theil1958]\
$\mathcal{F}$ smooth & on $X$ or $A$ of which at least one is confounding-removing & support-reducing & Ass. \[ass:identify\_f\] & nonlinear IV (e.g., NPREGIV [@NPREGIV-CRAN])\
$\mathcal{F}$ smooth and linearly extrapolates & on $X$ or $A$ of which at least one is confounding-removing & – & Ass. \[ass:identify\_f\] & **NILE** Section \[sec:nile\]\

For well-behaved interventions on $X$ which contain at least one confounding-removing intervention, estimating minimax solutions reduces to the well-studied problem of estimating causal relationships. One class of algorithms for this task is given by linear instrumental variable (IV) approaches. They assume that $\mathcal{F}$ is linear and require identifiability of the causal function (Assumption \[ass:identify\_f\]) via a rank condition on the observational distribution, see Appendix \[sec:IVconditions\]. Their target of inference is to estimate the causal function, which by Proposition \[prop:minimax\_equal\_causal\] will coincide with the minimax solution if the set ${\mathcal{I}}$ consists of well-behaved interventions with at least one of them being confounding-removing. A basic estimator for linear IV models is the two-stage least squares (TSLS) estimator, which minimizes the norm of the prediction residuals projected onto the subspace spanned by the observed instruments (TSLS objective). TSLS estimators are consistent but do not come with strong finite sample guarantees; e.g., they do not have finite moments in a just-identified setup [e.g., @mariano2001simultaneous]. K-class estimators [@Theil1958] have been proposed to overcome some of these issues. They minimize a linear combination of the residual sum of squares (OLS objective) and the TSLS objective. K-class estimators can be seen as utilizing a bias-variance trade-off. For fixed and non-trivial relative weights, they have, in a Gaussian setting, finite moments up to a certain order that depends on the sample-size and the number of predictors used. If the weights are such that the OLS objective is ignored asymptotically, they consistently estimate the causal parameter [e.g., @mariano2001simultaneous]. More recently, PULSE has been proposed [@jakobsen2020distributional], a data-driven procedure for choosing the relative weights such that the prediction residuals ‘just’ pass a test for simultaneous uncorrelatedness with the instruments. In cases where the minimax solution does not coincide with the causal function, only few algorithms exist. Anchor regression [@rothenhausler2018anchor] is a procedure that can be used when $\mathcal{F}$ and $\mathcal{G}$ are linear and $h_1$ is additive in the noise component. It finds the minimax solution if the set ${\mathcal{I}}$ consists of all interventions on $A$ up to a fixed intervention strength, and is applicable even if Assumption \[ass:identify\_f\] is not necessarily satisfied.

In a linear setting, where the regression coefficients differ between different environments, it is also possible to minimize the worst-case risk among the observed environments [@meinshausen2015maximin]. In its current formulation, this approach does not quite fit into the above framework, as it does not allow for changing distributions of the covariates. A summary of the mentioned methods and their assumptions is given in Table \[tab:learnability\].

If $\mathcal{F}$ is a nonlinear or non-parametric class of functions, the task of finding minimax solutions becomes more difficult. In cases where the causal function is among such solutions, this problem has been studied in the econometrics community. For example, [@newey2013nonparametric; @newey2003instrumental] treat the identifiability and estimation of causal functions in non-parametric function classes. Several non-parametric IV procedures exists, e.g., NPREGIV [@NPREGIV-CRAN] contains modified implementations of [@darolles2011nonparametric] and [@horowitz2011applied]. Identifiability and estimation of the causal function using nonlinear IV methods in parametric function classes is discussed in Appendix \[sec:IVconditions\]. Unlike in the linear case, most of the methods do not aim to extrapolate and only recover the causal function inside the support of $X$, that is, they cannot be used to predict interventions outside of this domain. In the following section, we propose a procedure that is able to extrapolate when $\mathcal{F}$ consists of functions which extend linearly outside of the support of $X$. In our simulations, we show that such an assumption can improve the prediction performance on the boundary of the support.

NILE {#sec:nile}
----

We have seen in Proposition \[prop:impossibility\_extrapolation\] that in order to generalize to interventions which extend the support of $X$, we require additional assumptions on the function class ${\mathcal{F}}$. In this section, we start from such assumptions and verify both theoretically and practically that they allows us to perform distribution generalization in the considered setup. Along the way, several choices can be made and usually several options are possible. We will see that our choices yield a method with competitive performance, but we do not claim optimality of our procedure. Several of our choices were partially made to keep the theoretical exposition simple and the method computationally efficient. We first consider the univariate case (i.e., $X$ and $A$ are real-valued) and comment later on the possibility to extend the methodology to higher dimensions. Unless specific background knowledge is given, it might be reasonable to assume that the causal function extends linearly outside a fixed interval $[a,b]$. By additionally imposing differentiability on ${\mathcal{F}}$, any function from ${\mathcal{F}}$ is uniquely defined by its values within $[a,b]$, see also Section \[sec:suff\_genX\_extrap\]. Given an estimate $f$ on $[a,b]$, the linear extrapolation property then yields a global estimate on the whole of $\mathbb{R}$. In principle, any class of differentiable functions can be used. Here, we assume that, on the interval $[a,b]$, the causal function $f$ is contained in the linear span of a B-spline basis. More formally, let $B = (B_1, ..., B_k)$ be a fixed B-spline basis on $[a,b]$, and define $\eta := (a,b,B)$. Our procedure assumes that the true causal function $f$ belongs to the function class ${\mathcal{F}}_{\eta} := \{f_\eta(\cdot; \theta) {\, : \,}\theta \in {\mathbb{R}}^k\}$, where for every $x \in {\mathbb{R}}$ and $\theta \in {\mathbb{R}}^k$, $f_\eta(x; \theta)$ is given as $$\label{eq:f_theta}
f_{\eta}(x;\theta) :=
\begin{cases}
B(a)^\top  \theta +  B^\prime (a)^\top  \theta (x - a) & \text{ if } x < a \\
B(x)^\top  \theta & \text{ if } x \in [a, b] \\
B(b)^\top  \theta + B^\prime (b)^\top \theta (x - b) & \text{ if } x >b,\\
\end{cases}$$ where $B^\prime := (B_1^\prime, \dots, B_k^\prime)$ denotes the component-wise derivative of $B$. In our algorithm, $\eta = (a, b, B)$ is a hyper-parameter, which can be set manually, or be chosen from data.

0 Let $({\mathbf{X}}, {\mathbf{Y}}, {\mathbf{A}}) \in {\mathbb{R}}^{n \times 3}$ be the observed data. It its most general form, the Nonlinear Intervention-robust Linear Extrapolator (NILE) solves an optimization problem of the form $$\label{eq:nile_general}
\operatorname*{argmin}_{{f_{\diamond}}\in {\mathcal{F}}} \underbrace{\norm{{\mathbf{Y}} - {f_{\diamond}}({\mathbf{X}})}_2^2}_{\text{prediction}} + 
\lambda \underbrace{\norm{{\mathbf{P}}({\mathbf{Y}} - {f_{\diamond}}({\mathbf{X}}))}_2^2}_{\text{invariance}} + 
\gamma \underbrace{\int ({f_{\diamond}}^{\prime \prime}(x))^2 dx}_{\text{smoothness + lin. extrap.}},$$ where ${\mathcal{F}}$ is some suitable function class, ${\mathbf{P}}$ is the (square) “hat-matrix” for a nonlinear regression of the residuals ${\mathbf{Y}} - {f_{\diamond}}({\mathbf{X}})$ onto ${\mathbf{A}}$, and $\lambda, \gamma  > 0$ are some tuning parameters. This estimator may be seen as a nonlinear version of the PULSE (see Section \[sec:existingmethods\]), with an additional constraint that enforces linear extrapolation. By choosing ${\mathcal{F}}$ as the linear span of a B-spline basis, the smoothness penalty in becomes a quadratic function of the spline coefficients [e.g., @fahrmeir2013regression], and the resulting optimization problem becomes strictly convex.

### Estimation procedure {#sec:estimation}

We now introduce our estimation procedure for fixed choices of all hyper-parameters. Section \[sec:algorithm\] describes how these can be chosen from data in practice. Let $({\mathbf{X}}, {\mathbf{Y}}, {\mathbf{A}}) \in {\mathbb{R}}^{n \times 3}$ be $n$ i.i.d. realizations sampled from a distribution over $(X,Y,A)$, let $\eta = (a,b,B)$ be fixed and assume that ${\mathrm{supp}}(X){\subseteq}[a,b]$. Our algorithm aims to learn the causal function $f_\eta(\cdot ; \theta^0) \in {\mathcal{F}}_\eta$, which is determined by the linear causal parameter $\theta^0$ of a $k$-dimensional vector of covariates $(B_1(X), \dots, B_k(X))$. From standard linear IV theory, it is known that at least $k$ instrumental variables are required to identify the $k$ causal parameters, see Appendix \[sec:IVconditions\]. We therefore artificially generate such instruments by nonlinearly transforming $A$, by using another B-spline basis $C = (C_1, \dots, C_k)$. The parameter $\theta^0$ can then be identified from the observational distribution under appropriate rank conditions, see Section \[sec:consistency\]. In that case, the hypothesis $H_0(\theta) : \theta= \theta^0$ is equivalent to the hypothesis $\tilde H_0(\theta) : {\mathbb{E}}[C(A)(Y - B(X)^\top \theta)] =
  0$. Let ${\mathbf{B}} \in {\mathbb{R}}^{n \times k}$ and ${\mathbf{C}} \in {\mathbb{R}}^{n \times k}$ be the associated design matrices, for each $i \in \{1, \dots, n\}$, $j \in \{1, \dots, k\}$ given as ${\mathbf{B}}_{ij} = B_j(X_i)$ and ${\mathbf{C}}_{ij} = C_{j}(A_i)$. A straightforward choice would be to construct the standard TSLS estimator, i.e., $\hat \theta$ as the minimizer of $\theta \mapsto \norm{{\mathbf{P}} ({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2$, where ${\mathbf{P}}$ is the projection matrix onto the columns of ${\mathbf{C}}$, see also [@hall2005generalized]. Even though this procedure may result in an asymptotically consistent estimator, there are several reasons why it may be suboptimal in a finite sample setting. First, the above estimator can have large finite sample bias, in particular if $k$ is large. Indeed, in the extreme case where $k = n$, and assuming that all columns in ${\mathbf{C}}$ are linearly independent, ${\mathbf{P}}$ is equal to the identity matrix, and $\hat \theta$ coincides with the OLS estimator. Second, since $\theta$ corresponds to the linear parameter of a spline basis, it seems reasonable to impose constraints on $\theta$ which enforce smoothness of the resulting spline function. Both of these points can be addressed by introducing additional penalties into the estimation procedure. Let therefore ${\mathbf{K}} \in {\mathbb{R}}^{k \times k}$ and ${\mathbf{M}} \in {\mathbb{R}}^{k \times k}$ be the matrices that are, for each $i,j \in \{1, \dots, k\}$, defined as ${\mathbf{K}}_{ij} = \int B^{\prime \prime}_i(x) B^{\prime \prime}_j(x) dx$ and ${\mathbf{M}}_{i j} = \int C^{\prime \prime}_{i}(a) C^{\prime \prime}_{j}(a) da$, and let $\gamma, \delta > 0$ be the respective penalties associated with ${\mathbf{K}}$ and ${\mathbf{M}}$. For $\lambda \geq 0$ and with $\mu := (\gamma, \delta, C)$, we then define the estimator $$\label{eq:thetahat}
\hat \theta^n_{\lambda, \eta, \mu}
:= \operatorname*{argmin}_{\theta \in {\mathbb{R}}^{k}} 
\norm{{\mathbf{Y}} - {\mathbf{B}} \theta }_2^2 + \lambda \norm{{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2 + \gamma \theta^\top {\mathbf{K}} \theta, \\$$ where ${\mathbf{P}}_\delta := {\mathbf{C}} ({\mathbf{C}}^\top {\mathbf{C}} + \delta {\mathbf{M}})^{-1} {\mathbf{C}}^\top$ is the ‘hat’-matrix for a penalized regression onto the columns of ${\mathbf{C}}$. By choice of ${\mathbf{K}}$, the term $\theta^\top {\mathbf{K}} \theta$ is equal to the integrated squared curvature of the spline function parametrized by $\theta$. The above may thus be seen as a nonlinear extension of K-class estimators [@Theil1958], with an additional penalty term which enforces linear extrapolation. In principle, the above approach extends to situations where $X$ and $A$ are higher-dimensional, in which case $B$ and $C$ consist of multivariate functions. For example, [@fahrmeir2013regression] propose the use of tensor product splines, and introduce multivariate smoothness penalties based on pairwise first- or second order parameter differences of basis functions which are close-by with respect to some suitably chosen metric. Similarly to , such penalties result in a convex optimization problem. However, due to the large number of involved variables, the optimization procedure becomes computationally burdensome already in small dimensions. Within the function class ${\mathcal{F}}_\eta$, the above defines the global estimate $f_\eta(\cdot; \hat \theta^n_{\lambda, \eta, \mu})$, for every $x \in {\mathbb{R}}$ given by $$\label{eq:fhat_theta}
f_\eta(x; \hat \theta^n_{\lambda, \eta, \mu}) :=
\begin{cases}
 B(a)^\top \hat \theta^n_{\lambda, \eta, \mu} +  B^\prime (a)^\top \theta^n_{\lambda, \eta, \mu} (x - a) & \text{ if } x < a \\
 B(x)^\top \hat \theta^n_{\lambda, \eta, \mu} & \text{ if } x \in [a,b] \\
 B(b)^\top \theta^n_{\lambda, \eta, \mu} + B^\prime (b)^\top \theta^n_{\lambda, \eta, \mu} (x - b) & \text{ if } x > b. \\
\end{cases}$$ We deliberately distinguish between three different groups of hyper-parameters $\eta$, $\mu$ and $\lambda$. The parameter $\eta = (a,b,B)$ defines the function class to which the causal function $f$ is assumed to belong. To prove consistency of our estimator, we require this function class to be correctly specified. In turn, the parameters $\lambda$ and $\mu=(\gamma, \delta, C)$ are algorithmic parameters that do not describe the statistical model. Their values only affects the finite sample behavior of our algorithm, whereas consistency is ensured as long as $C$ satisfies certain rank conditions, see Assumption \[ass:RankCondition\] in Section \[sec:consistency\]. In practice, $\gamma$ and $\delta$ are chosen via a cross-validation procedure, see Section \[sec:algorithm\]. The parameter $\lambda$ determines the relative contribution of the OLS and TSLS losses to the objective function. To choose $\lambda$ from data, we use an idea similar to the PULSE [@jakobsen2020distributional].

### Algorithm {#sec:algorithm}

Let for now $\eta, \mu$ be fixed. In the limit $\lambda \to \infty$, our estimation procedure becomes equivalent to minimizing the TSLS loss $\theta \mapsto \norm{{\mathbf{P}}_\delta ({\mathbf{Y}} - {\mathbf{B}}\theta)}_2^2$, which may be interpreted as searching for the parameter $\theta$ which complies ‘best’ with the hypothesis $\tilde{H}_0(\theta) : {\mathbb{E}}[C(A)(Y - B(X)^\top \theta)] = 0$. For finitely many data, following the idea introduced in [@jakobsen2020distributional], we propose to choose the value for $\lambda$ such that $\tilde{H}_0(\hat \theta^n_{\lambda, \eta, \mu})$ is just accepted (e.g., at a significance level $\alpha = 0.05$). That is, among all $\lambda \geq 0$ which result in an estimator that is not rejected as a candidate for the causal parameter, we chose the one which yields maximal contribution of the OLS loss to the objective function. More formally, let for every $\theta \in {\mathbb{R}}^k$, $T(\theta) = (T_n(\theta))_{n \in {\mathbb{N}}}$ be a statistical test at (asymptotic) level $\alpha$ for $\tilde{H}_0(\theta)$ with rejection threshold $q(\alpha)$. That is, $T_n(\theta)$ does not reject $\tilde{H}_0(\theta)$ if and only if $T_n(\theta) \leq q(\alpha)$. The penalty $\lambda^\star_n$ is then chosen in the following data-driven way $$\begin{aligned}
\lambda^\star_n := \inf \{\lambda \geq 0 : T_n(\hat \theta^n_{\lambda, \eta, \mu})\leq q(\alpha)\}.\end{aligned}$$ In general, $\lambda^\star_n$ is not guaranteed to be finite for an arbitrary test statistic $T_n$. Even for a reasonable test statistic it might happen that $T_n(\hat \theta^n_{\lambda, \eta, \mu} ) > q(\alpha)$ for all $\lambda \geq 0$; see [@jakobsen2020distributional] for further details. We can remedy the problem by reverting to another well-defined and consistent estimator, such as the TSLS (which minimizes the TSLS loss above) if $\lambda^\star_n$ is not finite. Furthermore, if $\lambda \mapsto T_n(\hat \theta^n_{\lambda, \eta, \mu})$ is monotonic, $\lambda^\star_n$ can be computed efficiently by a binary search procedure. In our algorithm, the test statistic $T$ and rejection threshold $q$ can be supplied by the user. Conditions on $T$ that are sufficient to yield a consistent estimator $f_\eta(\cdot , \hat \theta_{\lambda_n^\star, \mu, \eta})$, given that $\mathcal{F}_\eta$ is correctly specified, are presented in Section \[sec:consistency\]. Two choices of test statistics which are implemented in our code package can be found in Appendix \[sec:test\_statistic\]. For every $\gamma \geq 0$, let ${\mathbf{Q}}_\gamma = {\mathbf{B}} ({\mathbf{B}}^\top {\mathbf{B}} + \gamma {\mathbf{K}})^{-1} {\mathbf{B}}^\top$ be the ‘hat’-matrix for the penalized regression onto ${\mathbf{B}}$. Our algorithm then proceeds as follows.\

**input**: data $({\mathbf{X}}, {\mathbf{Y}}, {\mathbf{A}}) \in {\mathbb{R}}^{n \times 3}$ **options**: $k$, $T$, $q$, $\alpha$ **output**: $\hat{f}^n_{\text{NILE}}  := f_{\hat \eta}(\, \cdot \, ; \hat \theta^n_{\lambda_n^\star, \mu^n_{\text{CV}}, \hat \eta})$ defined by

The penalty parameter $\gamma^n_{\text{CV}}$ is chosen to minimize the out-of-sample mean squared error of the prediction model $\hat{{\mathbf{Y}}} = {\mathbf{Q}}_{\gamma} {\mathbf{Y}}$, which corresponds to the solution of for $\lambda = 0$. After choosing $\lambda_n^\star$, the objective function in increases by the term $\lambda_n^\star \norm{{\mathbf{P}}_{\delta_{\text{CV}}^n}({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2$. In order for the penalty term $\gamma \theta^\top {\mathbf{K}} \theta$ to impose the same degree of smoothness in the altered optimization problem, the penalty parameter $\gamma$ needs to be adjusted accordingly. The heuristic update in our algorithm is motivated by the simple observation that for all $\delta, \lambda \geq 0$, $\norm{{\mathbf{Y}} - {\mathbf{B}} \theta}_2^2 + \lambda \norm{{\mathbf{P}}_\delta ({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2 \leq (1 + \lambda) \norm{{\mathbf{Y}} - {\mathbf{B}} \theta}_2^2$.

### Asymptotic generalization (consistency) {#sec:consistency}

We now prove consistency of our estimator in the case where the hyper-parameters $(\eta, \mu)$ are fixed (rather than data-driven), and the function class ${\mathcal{F}}_{\eta}$ is correctly specified. Fix any $a<b$ and a basis $B=(B_1, \dots, B_k)$. Let $\eta_0 = (a,b,B)$ and let the model class be given by ${\mathcal{M}}= {\mathcal{F}}_{\eta_0} \times {\mathcal{G}}\times {\mathcal{H}}_1 \times {\mathcal{H}}_2 \times
\mathcal{Q}$, where ${\mathcal{F}}_{\eta_0}$ is as described in Section \[sec:nile\]. Assume that the data-generating model $M = ( f_{\eta_0}(\, \cdot \,; \theta^0),g,h_1,h_2,Q) \in {\mathcal{M}}$ induces an observational distribution ${\mathbb{P}}_M$ such that ${\mathrm{supp}}^M(X) \subseteq (a,b)$. Let further ${\mathcal{I}}$ be a set of interventions on $X$ or $A$, and let $\alpha \in (0,1)$ be a fixed significance level.

We prove asymptotic generalization (consistency) for an idealized version of the NILE estimator which utilizes $\eta_0$, rather than the data-driven values. Choose any $\delta,\gamma \geq 0$ and basis $C=(C_1,...,C_k)$ and let $\mu=(\delta,\gamma,C)$. We will make use of the following assumptions.

1.  $\forall \tilde{M}\in \cM$ s.t. $\bP_M = \bP_{\tilde M}$ : $\sup_{i\in\cI}{\mathbb{E}}_{\tilde{M}(i)} [X^2], \, \sup_{i\in\cI} \lambda_{\max}({\mathbb{E}}_{\tilde{M}(i)} [B(X)B(X)^\top ])<\infty$. \[ass:MaximumEigenValueBounded\]

2.  ${\mathbb{E}}_M[ B(X)B(X)^\top ]$, ${\mathbb{E}}_M[ C(A)C(A)^\top ]$ and ${\mathbb{E}}_M[ C(A)B(X)^\top  ]$ are of full rank.\[ass:RankCondition\]

<!-- -->

1.  $T(\theta)$ has uniform asymptotic power on any compact set of alternatives. \[ass:ConsistentTestStatistic\]

2.  $\lambda^\star_n := \inf\{\lambda\geq 0 : T_n(\hat \theta^n _{\lambda,\eta_{0},\mu}) \leq q(\alpha)\}$ is almost surely finite. \[ass:LambdaStarAlmostSurelyFinite\]

3.  $\lambda \mapsto T_n(\hat \theta^n _{\lambda,\eta_0,\mu})$ is weakly decreasing and $\theta \mapsto T_n(\theta)$ is continuous. \[ass:MonotonicityAndContinuityOfTest\]

Assumptions \[ass:MaximumEigenValueBounded\]–\[ass:RankCondition\] ensure consistency of the estimator as long as $\lambda^\star_n$ tends to infinity. Intuitively, in this case, we can apply arguments similar to those that prove consistency of the TSLS estimator. Assumptions \[ass:ConsistentTestStatistic\]–\[ass:MonotonicityAndContinuityOfTest\] ensure that consistency is achieved when choosing $\lambda^\star_n$ in the data-driven fashion described in Section \[sec:algorithm\]. In Assumption \[ass:MaximumEigenValueBounded\], $\lambda_{\max}$ denotes the largest eigenvalue. In words, the assumption states that, under each model $\tilde M\in \cM$ with $\bP_M = \bP_{\tilde M}$, there exists a finite upper bound on the variance of any linear combination of the basis functions $B(X)$, uniformly over all distributions induced by ${\mathcal{I}}$. The first two rank conditions of \[ass:RankCondition\] enable certain limiting arguments to be valid and they guarantee that estimators are asymptotically well-defined. The last rank condition of \[ass:RankCondition\] is the so-called rank condition for identification. It guarantees that $\theta^0$ is identified from the observational distribution in the sense that the hypothesis $H_0(\theta):\theta=\theta^0$ becomes equivalent with $\tilde{H}_0(\theta) : {\mathbb{E}}_M[C(A)(Y-B(X)^\top \theta)]=0$. \[ass:ConsistentTestStatistic\] means that for any compact set $K\subseteq {\mathbb{R}}^k$ with $\theta^0 \not \in K$ it holds that $\lim_{n\to \infty} P(\inf_{\theta\in K} T_n(\theta) \leq q(\alpha)) =0$. If the considered test has, in addition, a level guarantee, such as pointwise asympotic level, the interpretation of the finite sample estimator discussed in Section \[sec:algorithm\] remains valid (such level guarantee may potentially yield improved finite sample performance, too). \[ass:LambdaStarAlmostSurelyFinite\] is made to simplify the consistency proof. As previously discussed in Section \[sec:algorithm\], if \[ass:LambdaStarAlmostSurelyFinite\] is not satisfied, we can output another well-defined and consistent estimator on the event $(\lambda^\star_n=\infty)$, ensuring that consistency still holds. Under these conditions, we have the following asymptotic generalization guarantee.

\[thm:consis\] Let $\cI$ be a set of interventions on $X$ or $A$ of which at least one is confounding-removing. If assumptions \[ass:MaximumEigenValueBounded\]–\[ass:RankCondition\] and \[ass:ConsistentTestStatistic\]–\[ass:MonotonicityAndContinuityOfTest\] hold true, then, for any $\tilde{M} \in {\mathcal{M}}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, and any ${\varepsilon}> 0$, it holds that $${\mathbb{P}}_{M} \left( \big\vert \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}_{\eta_0}}\,
      \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \big\vert \leq {\varepsilon}\right) \to 1,$$ as $n\to \infty$. In the above event, only $\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}$ is stochastic.

### Experiments {#sec:experiments}

We now investigate the empirical performance of our proposed estimator, the NILE, with $k = 50$ spline basis functions. To choose $\lambda_n^\star$, we use the test statistic $T_n^2$, which tests the slightly stronger hypothesis $\bar{H}_0$, see Appendix \[sec:test\_statistic\]. In all experiments use the significance level $\alpha = 0.05$. We include two other approaches as baseline: (i) the method NPREGIV (using its default options) introduced in Section \[sec:existingmethods\], and (ii) a linearly extrapolating estimator of the ordinary regression of $Y$ on $X$ (which corresponds to the NILE with $\lambda^\star \equiv 0$). In all experiments, we generate data sets of size $n=200$ as independent replications from $$\label{eq:sim_model}
A:= {\varepsilon}_A, \quad H:= {\varepsilon}_H, \quad X := \alpha_A A + \alpha_H H + \alpha_{\varepsilon}{\varepsilon}_X, \quad Y := f (X) + 0.3  H + 0.2  {\varepsilon}_Y,$$ where $({\varepsilon}_A, {\varepsilon}_H, {\varepsilon}_X, {\varepsilon}_Y)$ are jointly independent with $\mathcal{U}(-1,1)$ marginals. To make results comparable across different parameter settings, we impose the constraint $\alpha_A^2 + \alpha_H^2  +\alpha_{\varepsilon}^2 = 1$, which ensures that in all models, $X$ has variance $1/3$. The function $f$ is drawn from the linear span of a basis of four natural cubic splines with knots placed equidistantly within the $90\%$ inner quantile range of $X$. By well-known properties of natural splines, any such function extends linearly outside the boundary knots. Figure \[fig:overlay\_estimates\_and\_varying\_confounding\] (left) shows an example data set from , where the causal function is indicated in green. We additionally display estimates obtained by each of the considered methods, based on 20 i.i.d. datasets. Due to the confounding variable $H$, the OLS estimator is clearly biased. NPREGIV exploits $A$ as an instrumental variable and obtains good results within the support of the observed data. Due to its non-parametric nature, however, it cannot extrapolate outside this domain. The NILE estimator exploits the linear extrapolation assumption on $f$ to produce global estimates.

![A sample dataset from the model with $\alpha_A = \sqrt{1/3}$, $\alpha_H = \sqrt{2/3}$, $\alpha_{\varepsilon}= 0$. The true causal function is indicated by a green dashed line. For each method, we show 20 estimates of this function, each based on an independent sample from . For values within the support of the training data (vertical dashed lines mark the inner 90% quantile range), NPREGIV correctly estimates the causal function well. As expected, when moving outside the support of $X$, the estimates become unreliable, and we gain an increasing advantage by exploiting the linear extrapolation assumed by the NILE.[]{data-label="fig:overlay_estimates_and_varying_confounding"}](overlay_estimates){width="\linewidth"}

We further investigate the empirical worst-case mean squared error across several different models of the form . That is, for a fixed set of parameters $(\alpha_A, \alpha_H, \alpha_{{\varepsilon}})$, we construct several models $M_1, \dots, M_N$ of the form by randomly sampling causal functions $f_1, \dots, f_N$ (see Appendix \[sec:additional\_experiments\] for further details on the sampling procedure). For every $x \in [0,2]$, let ${\mathcal{I}}_x$ denote the set of hard interventions which set $X$ to some fixed value in $[-x,x]$. We then characterize the performance of each method using the average (across different models) worst-case mean squared error (across the interventions in ${\mathcal{I}}_x$), i.e., for each estimator $\hat{f}$, we consider $$\label{eq:experiments_risk}
\frac{1}{N} \sum_{j=1}^N \sup_{i \in {\mathcal{I}}_x} {\mathbb{E}}_{M_j(i)} \big[ (Y - \hat{f}(X))^2 \big] 
= {\mathbb{E}}[\xi_Y^2] + \frac{1}{N} \sum_{j=1}^N \sup_{\tilde{x} \in [-x,x]}  (f_j(\tilde{x}) - \hat{f}(\tilde{x}))^2,$$ where $\xi_Y :=  0.3  H + 0.2  {\varepsilon}_Y$ is the noise term for $Y$ (which is fixed across all experiments). In practice, we evaluate the functions $\hat{f}$, $f_1, \dots, f_N$ on a fine grid on $[-x,x]$ to approximate the above supremum. Figure \[fig:overlay\_estimates\_and\_varying\_confounding2\] plots the average worst-case mean squared error versus intervention strength for different parameter settings. The optimal worst-case mean squared error ${\mathbb{E}}[\xi_Y^2]$ is indicated by a green dashed line. The results show that the linear extrapolation property of the NILE estimator is beneficial in particular for strong interventions. In the case of no confounding ($\alpha_H = 0$), the minimax solution coincides with the regression of $Y$ on $X$, hence even the OLS estimator yields good predictive performance. In this case, the hypothesis $\bar{H}_0(\hat \theta^n_{\lambda, \delta^n_{\text{CV}}, \gamma^n_{\text{CV}}})$ is accepted already for small values of $\lambda$ (in this experiment, the empirical average of $\lambda^\star_n$ equals 0.015), and the NILE estimator becomes indistinguishable from the OLS. As the confounding strength increases, the OLS becomes increasingly biased, and the NILE objective function differs more notably from the OLS (average $\lambda^\star_n$ of 2.412 and 5.136, respectively). The method NPREGIV slightly outperforms the NILE inside the support of the observed data, but drops in performance for stronger interventions. We believe that the increase in extrapolation performance of the NILE for stronger confounding (increasing $\alpha_H$) might stem from the fact that, as the $\lambda_n^\star$ increases, also the smoothness penalty $\gamma$ increases, see Algorithm \[alg:nile\]. While this results in slightly worse in-sample prediction, it seems beneficial for extrapolation (at least for the particular function class that we consider). We do not claim that our algorithm has theoretical guarantees which explain this increase in performance.

![ Predictive performance under confounding-removing interventions on $X$ for different confounding- and intervention strengths (see alpha values in the grey panel on top). The right panel corresponds to the same parameter setting as in Figure \[fig:overlay\_estimates\_and\_varying\_confounding\]. The plots in each panel are based on data sets of size $n=200$, generated from $N = 100$ different models of the form . For each model, we draw a different function $f$, resulting in a different minimax solution (see Appendix \[sec:additional\_experiments\] for details on the sampling procedure). The performances under individual models are shown by thin lines; the average performance across all models is indicated by thick lines. In all considered models, the optimal prediction error is equal to ${\mathbb{E}}[\xi_Y^2]$ (green dashed line). The grey area indicates the inner 90 % quantile range of $X$ in the training distribution; the white area can be seen as an area of generalization. []{data-label="fig:overlay_estimates_and_varying_confounding2"}](varying_confounding_with_variability){width="\linewidth"}

In the case, where all exogenous noise comes from the unobserved variable ${\varepsilon}_X$ (i.e., $\alpha_A$ = 0), the NILE coincides with the OLS estimator. In such settings, standard IV methods are known to perform poorly, although also the NPREGIV method seems robust to such scenarios. As the instrument strength increases, the NILE clearly outperforms OLS and NPREGIV for interventions on $X$ which include values outside the training data.

Discussion and future work
==========================

In many real world problems, the test distribution may differ from the training distribution. This requires statistical methods that come with a provable guarantee in such a setting. It is possible to characterize robustness by considering predictive performance for distribution that are close to the training distribution in terms of standard divergences or metrics, such as KL divergences or Wasserstein distance. As an alternative view point, we have introduced a novel framework that formalizes the task of distribution generalization when considering distributions that are induced by a set of interventions. Based on the concept of modularity, interventions modify parts of the joint distribution and leave other parts invariant. Thereby, they impose constraints on the changes of the distributions that are qualitatively different from considering balls in a certain metric. As such, we see them as a useful language to describe realistic changes between training and test distributions. Our framework is general in that it allows us to model a wide range of causal models and interventions, which do not need to be known beforehand. We have proved several generalization guarantees, some of which show robustness for distributions that are not close to the training distribution by considering almost any of the standard metrics. We have further proved impossibility results that indicate the limits of what is possible to learn from the training distribution. In particular, in nonlinear models, strong assumptions are required for distribution generalization to a different support of the covariates. As such, methods such as anchor regression cannot be expected to work in nonlinear models, unless strong restrictions are placed on the function class ${\mathcal{G}}$.

![Predictive performance for varying instrument strength. If the instruments have no influence on $X$ ($\alpha_A = 0$), the second term in the objective function is effectively constant in $\theta$, and the NILE therefore coincides with the OLS estimator (which uses $\lambda=0$). This guards the NILE against the large variance which most IV estimators suffer from in a weak instrument setting. For increasing influence of $A$, it clearly outperforms both alternative methods for large intervention strengths. []{data-label="fig:weak_instruments"}](varying_instrument_with_variability){width="\linewidth"}

Our work can be extended into several directions. It may, for example, be worthwhile to investigate the sharpness of the bounds we provide in Section \[sec:suffapprox\] and other extrapolation assumptions on $\mathcal{F}$. While our results can be applied to situations where causal background knowledge is available, via a transformation of SCMs, our analysis is deliberately agnostic about such information. It would be interesting to see whether stronger theoretical results can be obtained by including causal background information. Finally, it could be worthwhile to investigate whether NILE, which outperforms existing approaches with respect to extrapolation, can be combined with non-parametric methods. This could yield an even better performance on estimating the causal function within the support of the covariates.

We view our work as a step towards understanding the problem of distribution generalization. We hope that considering the concepts of interventions may help to shed further light into the question under which assumptions it is possible to generalize knowledge that was acquired during training to a different test distribution.

Acknowledgments {#acknowledgments .unnumbered}
===============

We thank Thomas Kneib for helpful discussions. RC and JP were supported by a research grant (18968) from VILLUM FONDEN.

Transforming causal models {#sec:causal_relations_X}
==========================

As illustrated in Remark \[rem:model\], our framework is able to model cases where causal relations between the observed variables are given explicitly, e.g., by an SCM. The key insight is that most of these causal relations can be absorbed by the hidden confounding $H$ on which we make few restrictions. To show how this can be done in a general setting, let us consider the following SCM $$\begin{aligned}
\label{eq:SCM_full_causal}
  \begin{split}
    A\coloneqq {{\varepsilon}}_A\qquad\qquad
    &X\coloneqq w(X, Y) + g(A) + h_2(H, {{\varepsilon}}_X)\\
    H\coloneqq {{\varepsilon}}_H\qquad\qquad
    &Y\coloneqq f(X) + h_1(H, {{\varepsilon}}_Y).
  \end{split}\end{aligned}$$ Assume that this SCM is uniquely solvable in the sense that there exists a unique function $F$ such that $(A, H, X, Y)=F({{\varepsilon}}_A,{{\varepsilon}}_H,{{\varepsilon}}_X,{{\varepsilon}}_Y)$ almost surely, see [@Bongers2016b] for more details. Denote by $F_X$ the coordinates of $F$ that correspond to the $X$ variable (i.e., the coordinates from $r+q+1$ to $r+q+d$). Assume further that there exist functions $\tilde{g}$ and $\tilde{h}_2$ such that $$\label{eq:linear_decomposition_inverse}
    F_X({{\varepsilon}}_A,{{\varepsilon}}_H,{{\varepsilon}}_X,{{\varepsilon}}_Y)=\tilde{g}({{\varepsilon}}_A) + \tilde{h}_2(({{\varepsilon}}_H, {{\varepsilon}}_Y), {{\varepsilon}}_X).$$ This decomposition is not always possible, but it exists in the following settings, for example: (i) *There are no $A$ variables.* As discussed in Section \[sec:setup\] our framework also works if no $A$ variables exist. In these cases, the additive decomposition becomes trivial. (ii) *There are further constraints on the full SCM.* The additive decomposition holds if, for example, $w$ is a linear function or $A$ only enters the structural assignments of covariates $X$ which have at most $Y$ as a descendant.

Using the decomposition in , we can define the following SCM $$\begin{aligned}
\label{eq:SCM_full_causal_simple}
\begin{split}
  A\coloneqq {{\varepsilon}}_A\qquad\qquad
     &X\coloneqq \tilde{g}(A) + \tilde{h}_2(\tilde{H}, {{\varepsilon}}_X)\\
  \tilde{H}\coloneqq {{\varepsilon}}_{\tilde{H}}\qquad\qquad
     &Y\coloneqq f(X) + h_1(\tilde{H}),
\end{split}\end{aligned}$$ where ${{\varepsilon}}_{\tilde{H}}$ has the same distribution as $({{\varepsilon}}_H, {{\varepsilon}}_Y)$ in the previous model. This model fits the framework described in Section \[sec:setup\], where the noise term in $Y$ is now taken to be constantly zero. Both SCMs  and  induce the same observational distribution and the same function $f$ appears in the assignments of $Y$.

It is further possible to express the set of interventions on the covariates $X$ in the original SCM  as a set of interventions on the covariates in the reduced SCM . The description of a class of interventions in the full SCM may, however, become more complex if we consider them in the reduced SCM . In particular, to apply the developed methodology, one needs to check whether the interventions in the reduced SCM is a well-behaved set of interventions (this is not necessarily the case) and how the support of all $X$ variables behaves under that specific intervention. We now discuss the case that the causal graph induced by the full SCM is a directed acyclic graph (DAG).

*Intervention type.*First, we consider which types of interventions in  translate to well-behaved interventions in . Importantly, interventions on $A$ in the full SCM reduce to regular interventions $A$ also in the reduced SCM. Similarly, performing hard interventions on all components of $X$ in the full SCM leads to the same intervention in the reduced SCM, which is in particular both confounding-removing and confounding-preserving. For interventions on subsets of the $X$, this is not always the case. To see that, consider the following example.

$$\begin{aligned}
\label{eq:ex_causal_1}
  \begin{split}
    A&\coloneqq {{\varepsilon}}_A\\
    X_1 &\coloneqq {{\varepsilon}}_1, \quad
    X_2  \coloneqq Y + {\varepsilon}_2\\
    Y &\coloneqq X_1 + {\varepsilon}_Y
  \end{split}
  \end{aligned}$$

$$\begin{aligned}
\label{eq:ex_causal_simple_1}
  \begin{split}
    A &\coloneqq {{\varepsilon}}_A, \quad
    H \coloneqq {{\varepsilon}}_Y \\
    X &\coloneqq ({{\varepsilon}}_1,  H + {{\varepsilon}}_1 + {{\varepsilon}}_2)\\
    Y &\coloneqq X_1 + H
  \end{split}
  \end{aligned}$$

with ${{\varepsilon}}_A, {\varepsilon}_1, {\varepsilon}_2, {{\varepsilon}}_Y\overset{i.i.d.}{\sim}
\mathcal{N}(0,1)$, where represents the full SCM and corresponds to the reduced SCM using our framework. Consider now, in the full SCM, the intervention $X_1 \coloneqq i$, for some $i\in{\mathbb{R}}$. In the reduced SCM, this intervention corresponds to the intervention $X = (X_1, X_2) \coloneqq (i, H + i + {{\varepsilon}}_2)$, which is neither confounding-preserving nor confounding-removing.[^7] On the other hand, any intervention on $X_2$ or $A$ in the full SCM model corresponds to the same intervention in the reduced SCM. We can generalize these observations to the following statements:

-   *Interventions on A:* If we intervene on $A$ in the full SCM (i.e., by replacing the structural assignment of $A$ with $\psi^i(I^i, {{\varepsilon}}_A^i)$), then this translates to an equivalent intervention in the reduced SCM .

-   *Hard interventions on all X:* If we intervene on all $X$ in the full SCM by replacing the structural assignment of $X$ with an independent random variable $I\in{\mathbb{R}}^d$, then this translates to the same intervention in the reduced SCM which is confounding-removing.

-   *No X is a descendant of $Y$ and there is no unobserved confounding $H$:* If we intervene on $X$ in the full SCM (i.e., by replacing the structural assignment of $X$ with $\psi^i(g, A^i, {{\varepsilon}}^i_X ,I^i)$), then this translates to a potentially different but confounding-removing intervention in the reduced SCM . This is because the reduced SCM does not include unobserved variables $H$ in this case.

-   *Hard interventions on a variable $X_j$ which has at most $Y$ as a descendant:* If we intervene on $X_j$ in the full SCM by replacing the structural assignment of $X_j$ with an independent random variable $I$, then this intervention translates to a potentially different but confounding-preserving intervention.

Other settings may yield well-behaved interventions, too, but may require more assumptions on the full SCM model or further restrictions on the intervention classes. *Intervention support.*A support-reducing intervention in the full SCM can translate to a support-extending intervention in the reduced SCM. Consider the following example.

$$\begin{aligned}
\label{eq:ex_causal_2}
      \begin{split}
        X_1 &\coloneqq {{\varepsilon}}_1 \\
        X_2 &\coloneqq X_1 + \mathbf{1}\{X_1 = 0.5\}\\
        Y &\coloneqq X_2 + {\varepsilon}_Y
      \end{split}
    \end{aligned}$$

$$\begin{aligned}
\label{eq:ex_causal_simple_2}
      \begin{split}
        X &\coloneqq ({{\varepsilon}}_1, {{\varepsilon}}_1 + \mathbf{1}\{{{\varepsilon}}_1 = 0.5\})\\
        Y &\coloneqq X_2 + {{\varepsilon}}_Y,
      \end{split}
    \end{aligned}$$

with ${{\varepsilon}}_1, {{\varepsilon}}_Y\overset{i.i.d.}{\sim} \mathcal{U}(0,1)$. As before, represents the full SCM, whereas corresponds to the reduced SCM converted to fit our framework. Under the observational distribution, the support of $X_1$ and $X_2$ is equal to the open interval $(0, 1)$. Consider now the support-reducing intervention $X_1:= 0.5$ in . Within our framework, such an intervention would correspond to the intervention $X = (X_1, X_2) := (0.5, 1.5)$, which is support-extending. This example is rather special in that the SCM consists of a function that changes on a null set of the observational distribution. With appropriate assumptions to exclude similar degenerate cases, it is possible to show that support-reducing interventions in correspond to support-reducing interventions within our framework .

Sufficient conditions for Assumption 1 in IV settings {#sec:IVconditions}
=====================================================

Assumption \[ass:identify\_f\] states that $f$ is identified on the support of $X$ from the observational distribution of $(Y,X,A)$. Whether this assumption is satisfied depends on the structure of $\cF$ but also on the other function classes $\cG,\cH_1,\cH_2$ and $\mathcal{Q}$ that make up the model class $\cM$ from which we assume that the distribution of $(Y,X,A)$ is generated.

Identifiability of the causal function in the presence of instrumental variables is a well-studied problem in econometrics literature. Most prominent is the literature on identification in linear SCMs [e.g., @fisher1966identification; @greene2003econometric]. However, identification has also been studied for various other parametric function classes. We say that $\mathcal{F}$ is a parametric function class if it can be parametrized by some finite dimensional parameter set $\Theta \subseteq {\mathbb{R}}^p$. We here consider classes of the form $$\begin{aligned}
    \mathcal{F} := \{ f(\cdot ,\theta):{\mathbb{R}}^d \to {\mathbb{R}}\,\vert\, \theta :\Theta \rightarrow {\mathbb{R}}, \theta \mapsto f(x,\theta) \text{ is } C^2 \text{ for all } x\in {\mathbb{R}}^d \}.\end{aligned}$$ Consistent estimation of the parameter $\theta_0$ using instrumental variables in such function classes has been studied extensively in the econometric literature [e.g., @amemiya1974nonlinear; @jorgenson1974efficient; @kelejian1971two]. These works also contain rigorous results on how instrumental variable estimators of $\theta_0$ are constructed and under which conditions consistency (and thus identifiability) holds. Here, we give an argument on why the presence of the exogenous variables $A$ yields identifiability under certain regularity conditions. Assume that ${\mathbb{E}}[h_1(H, {\varepsilon}_Y)|A]=0$, which implies that the true causal function $f(\cdot,\theta_0)$ satisfies the population orthogonality condition $$\begin{aligned}
 \label{Eq:PopOrthCondNonLinearIV}
    {\mathbb{E}}[l(A)^\top (Y-f(X,\theta_0))] = {\mathbb{E}}\big[l(A)^\top {\mathbb{E}}[h_1(H, {\varepsilon}_Y)|A]\big]= 0,\end{aligned}$$ for some measurable mapping $l:{\mathbb{R}}^q\to {\mathbb{R}}^g$, for some $ g \in \mathbb{N}_{>0}$. Clearly, $\theta_0$ is identified from the observational distribution if the map $\theta \mapsto {\mathbb{E}}[l(A)^\top (Y-f(X,\theta))]$ is zero if and only if $\theta=\theta_0$. Furthermore, since $\theta\mapsto f(x,\theta)$ is differentiable for all $x\in {\mathbb{R}}^d$, the mean value theorem yields that, for any $\theta\in \Theta$ and $x\in {\mathbb{R}}^d$, there exists an intermediate point $\tilde{\theta}(x,\theta,\theta_0)$ on the line segment between $\theta$ and $\theta_0$ such that $$f(x,\theta) - f(x,\theta_0) = D_\theta f(x,\tilde{\theta}(x,\theta,\theta_0))(\theta-\theta_0),$$ where, for each $x\in {\mathbb{R}}^d$, $D_\theta f(x,\theta)\in{\mathbb{R}}^{1\times p}$ is the derivative of $\theta\mapsto f(x,\theta)$ evaluated in $\theta$. Composing the above expression with the random vector $X$, multiplying with $l(A)$ and taking expectations yields that $${\mathbb{E}}[l(A)(Y-f(X,\theta_0))] - {\mathbb{E}}[l(A)(Y-f(X,\theta))]=  {\mathbb{E}}[l(A)D_\theta f(X,\tilde{\theta}(X,\theta,\theta_0))](\theta_0-\theta).$$ Hence, if $ {\mathbb{E}}[l(A)D_\theta f(X,\tilde{\theta}(X,\theta,\theta_0))]\in
{\mathbb{R}}^{g\times p}$ is of rank $p$ for all $\theta\in\Theta$ (which implies $g \geq p$), then $\theta_0$ is identifiable as it is the only parameter that satisfies the population orthogonality condition of . As $\theta_0$ uniquely determines the entire function, we get identifiability of $f\equiv f(\cdot,\theta_0)$, not only on the support of $X$ but the entire domain ${\mathbb{R}}^d$, i.e., both Assumptions \[ass:identify\_f\] and \[ass:gen\_f\] are satisfied. In the case that $\theta \mapsto f(x,\theta)$ is linear, i.e.$f(x,\theta) = f(x)^T \theta$ for all $x\in {\mathbb{R}}^d$, the above rank condition reduces to ${\mathbb{E}}[l(A)f(X)^T]\in {\mathbb{R}}^{g\times p}$ having rank $p$ (again, implying that $g \geq p$). Furthermore, when $(x,\theta)\mapsto f(x,\theta)$ is bilinear, a reparametrization of the parameter space ensures that $f(x,\theta)= x^T \theta$ for $\theta\in\Theta \subseteq {\mathbb{R}}^d$. In this case, the rank condition can be reduced to the well-known rank condition for identification in a linear SCM, namely that ${\mathbb{E}}[AX^T] \in {\mathbb{R}}^{q\times p}$ is of rank $p$.

Finally, identifiability and methods of consistent estimation of the causal function have also been studied for non-parametric function classes. The conditions for identification are rather technical, however, and we refer the reader to [@newey2013nonparametric; @newey2003instrumental] for further details.

Choice of test statistic {#sec:test_statistic}
========================

By considering the variables $B(X) = (B_1(X), \dots, B_k(X))$ and $C(A) = (C_1(A), \dots, C_k(A))$ as vectors of covariates and instruments, respectively, our setting in Section \[sec:nile\] reduces to the classical (just-identified) linear IV setting. We could therefore use a test statistics similar to the one propsed by the PULSE [@jakobsen2020distributional]. With a notation that is slightly adapted to our setting, this estimator tests $\tilde{H}_0(\theta)$ using the test statistic $$\begin{aligned}
T^1_n(\theta) = c(n) \frac{\norm{{\mathbf{P}} ({\mathbf{Y}} - {\mathbf{B}}\theta)}_2^2}{\norm{{\mathbf{Y}} - {\mathbf{B}}\theta}_2^2},\end{aligned}$$ where ${\mathbf{P}}$ is the projection onto the columns of ${\mathbf{C}}$, and $c(n)$ is some function with $c(n) \sim n$ as $n\to\infty$. Under the null hypothesis, $T^1_n$ converges in distribution to the $\chi^2_{k}$ distribution, and diverges to infinity in probability under the general alternative. Using this test statistic, $\tilde{H}_0(\theta)$ is rejected if and only if $T^1_n(\theta)> q(\alpha)$, where $q(\alpha)$ is the $(1-\alpha)$-quantile of the $\chi^2_{k}$ distribution. The acceptance region of this test statistic is asymptotically equivalent with the confidence region of the Anderson-Rubin test [@anderson1949estimation] for the causal parameter $\theta^0$. Using the above test results in a consistent estimator for $\theta^0$ [@jakobsen2020distributional Theorem 3.12]; the proof exploits the particular form of $T^1_n$ without explicitly imposing that assumptions \[ass:ConsistentTestStatistic\] and \[ass:LambdaStarAlmostSurelyFinite\] hold.

If the number $k$ of basis functions is large, however, numerical experiments suggest that the above test has low power in finite sample settings. As default, our algorithm therefore uses a different test based on a penalized regression approach. This test has been proposed in [@chen2014note] for inference in nonparametric regression models. We now introduce this procedure with a notation that is adapted to our setting. For every $\theta \in {\mathbb{R}}^k$, let $R_\theta = Y - B(X)^\top \theta$ be the residual associated with $\theta$. We then test the slightly stronger hypothesis $$\bar{H}_0(\theta): \text{ there exists } \sigma_\theta^2>0 \text{ such that }  {\mathbb{E}}[R_\theta {\, \vert \,}A] {\ensuremath{\stackrel{\text{a.s.}}{=}}}0  \text{ and } \text{Var}[R_\theta {\, \vert \,}A] = \sigma_\theta^2$$ against the alternative that ${\mathbb{E}}[R_\theta {\, \vert \,}A] = m(A)$ for some smooth function $m$. To see that the above hypothesis implies $\tilde{H}_0(\theta)$ (and therefore $H_0(\theta)$, see Section \[sec:estimation\]), let $\theta \in {\mathbb{R}}^k$ be such that $\bar{H}_0(\theta)$ holds true. Then, $${\mathbb{E}}[C(A)(Y - B(X)^\top \theta)] = {\mathbb{E}}[C(A) R_\theta] = {\mathbb{E}}[{\mathbb{E}}[C(A) R_\theta {\, \vert \,}A]] = {\mathbb{E}}[C(A) {\mathbb{E}}[R_\theta {\, \vert \,}A]] = 0,$$ showing that also $\tilde{H}_0(\theta)$ holds true. Thus, if $\tilde{H}_0(\theta)$ is false, then also $\bar H_0(\theta)$ is false. As a test statistic $T^2_n(\theta)$ for $\bar{H}_0(\theta)$, we use (up to a normalization) the squared norm of a penalized regression estimate of $m$, evaluated at the data ${\mathbf{A}}$, i.e., the TSLS loss $\norm{{\mathbf{P}}_\delta ({\mathbf{Y}} - {\mathbf{B}}\theta)}_2^2$. In the fixed design case, where ${\mathbf{A}}$ is non-random, it has been shown that, under $\bar{H}_0(\theta)$ and certain additional regularity conditions, it holds that $$\frac{\norm{{\mathbf{P}}_\delta ({\mathbf{Y}} - {\mathbf{B}}\theta)}_2^2 - \sigma_\theta^2 c_n}{\sigma_\theta^2 d_n} \stackrel{\text{d}}{\longrightarrow} \mathcal{N}(0,1),$$ where $c_n$ and $d_n$ are known functions of ${\mathbf{C}}$, ${\mathbf{M}}$ and $\delta$ [@chen2014note Theorem 1]. The authors further state that the above convergence is unaffected by exchanging $\sigma_\theta^2$ with a consistent estimator $\hat{\sigma}_\theta^2$, which motivates our use of the test statistic $$T^2_n(\theta) := \frac{\norm{{\mathbf{P}}_\delta ({\mathbf{Y}} - {\mathbf{B}}\theta)}_2^2 - \hat \sigma_{\theta,n}^2 c_n}{\hat \sigma_{\theta,n}^2 d_n},$$ where $\hat \sigma_{\theta,n}^2 := \frac{1}{n-1} \sum_{i=1}^n \norm{({\mathbf{I}}_n - {\mathbf{P}}_\delta)({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2$. As a rejection threshold $q(\alpha)$ we use the $1-\alpha$ quantile of a standard normal distribution. For results on the asymptotic power of the test defined by $T^2$, we refer to Section 2.3 in [@chen2014note].

In our software package, both of the above tests are available options.

Their result is formalized in a fixed design setup. With a slight change of notation to suit our setting, they consider several instances of a response variable $R_1, \dots, R_n$ and fixed values of a covariate $a_1, \dots, a_n \in [c,d]$. They assume the model $$R_i = m(a_i) + \xi_i, \qquad {\mathbb{E}}[\xi_i] = 0, \qquad \text{Var}(\xi_i) = \sigma^2, \qquad i = 1, \dots, n,$$ for some unspecified smooth function $m$, and consider the null hypothesis $$\tilde{H}_0: m(a) = 0 \qquad \text{for all} \qquad a \in [c,d].$$ To test the above hypothesis, they use a B-spline estimate of $m$ to derive a test statistic which has an asymptotic standard normal distribution under $\tilde{H}_0$. We let $C_1, \dots, C_k$ denote the B-spline basis, and use ${\mathbf{M}}$ to denote the corresponding penalty matrix, defined analogously as in Section \[sec:estimation\].

$$\tilde{H}_0(\theta): 
\begin{cases}
\text{there exists } \xi {\perp \!\!\! \perp}A \text{ with } {\mathbb{E}}[\xi] = 0, \text{Var}(\xi) = \sigma^2 \text{ and } \\ 
\text{a smooth function } m \text{ such that } R_\theta = m(A) + \xi
\end{cases}$$

For every $\theta \in {\mathbb{R}}^k$, let $R_\theta = Y - B(X)^\top \theta$ be the residual associated with $\theta$, and consider the slightly stronger hypothesis

To see that $\tilde{H}_0$ implies $H_0$, let $\theta \in {\mathbb{R}}^k$ be such that $\tilde{H}_0(\theta)$ holds true. Then, $${\mathbb{E}}[C(A)(Y - B(X)^\top \theta)] = {\mathbb{E}}[C(A) R_\theta] = {\mathbb{E}}[{\mathbb{E}}[C(A) R_\theta {\, \vert \,}A]] = {\mathbb{E}}[C(A) {\mathbb{E}}[R_\theta {\, \vert \,}A]] = 0,$$ showing that also $H_0(\theta)$ holds true.

the class of distributions defined by $${\mathcal{P}}:= 
\left \lbrace
\begin{array}{l}
\text{distr. over } (R,A) \in {\mathbb{R}}^2
\text{ s.t. there exists } \xi {\perp \!\!\! \perp}A \text{ with } {\mathbb{E}}[\xi] = 0, \text{Var}(\xi) = \sigma^2 \\ 
\text{and a smooth function } m \text{ such that } R = m(A) + \xi
\end{array}
\right \rbrace.$$ Within this class of distributions, we can define $\tilde{H}_0(\theta): (R_\theta, A) \in {\mathcal{P}}$ with $m = 0$

, and let $m_\theta$ denote a conditional mean function of $R_\theta$ given $A$. The true causal coefficient $\theta^0$ yields the residual $R_{\theta^0} = \xi_Y$, which is independent of $A$. Since $\xi_Y$ has zero mean, it holds that $m_{\theta^0} = 0$ almost surely. For every $\theta$, let $$\tilde{H}_0(\theta) : m_{\theta} = 0 \text{ almost surely.}$$

We use a test for $H_0$ designed for inference in penalized spline regression models. [@chen2014note] propose a test statistic $T_n$ which, up to a normalization depending on $n$, ${\mathbf{C}}$ and ${\mathbf{M}}$, is equal to the TSLS loss $\norm{{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2$, and show that, under $H_0(\theta)$, this test statistic converges to a standard normal distribution [Theorem 1 @chen2014note].

Addition to experiments {#sec:additional_experiments}
=======================

Sampling of the causal function {#sec:exp_sampling}
-------------------------------

To ensure linear extrapolation of the causal function, we have chosen a function class consisting of natural cubic splines, which, by construction, extrapolate linearly outside the boundary knots. We now describe in detail how we sample functions from this class for the experiments in Section \[sec:experiments\]. Let $q_{\min}$ and $q_{\max}$ be the respective $5\%$- and $95\%$ quantiles of $X$, and let $B_1, \dots, B_4$ be a basis of natural cubic splines corresponding to 5 knots placed equidistantly between $q_{\min}$ and $q_{\max}$. We then sample coefficients $\beta_i {\overset{\text{iid}}{\sim}}\mathcal{U}(-1,1)$, $i = 1, \dots, 4$, and construct $f$ as $f = \sum_{i=1}^4 \beta_i B_i$. For illustration, we have included 18 realizations in Figure \[fig:f\_samples\].

![ The plots show independent realizations of the causal function that is used in all our experiments. These are sampled from a linear space of natural cubic splines, as described in Appendix \[sec:exp\_sampling\]. To ensure a fair comparison with the alternative method, NPREGIV, the true causal function is chosen from a model class different from the one assumed by the NILE.[]{data-label="fig:f_samples"}](f_samples){width="\linewidth"}

Violations of the linear extrapolation assumption {#sec:exp_violation_lin_extrap}
-------------------------------------------------

We have assumed that the true causal function extrapolates linearly outside the 90% quantile range of $X$. We now investigate the performance of our method for violations of this assumption. To do so, we again sample from the model , with $\alpha_A = \alpha_H = \alpha_{\varepsilon}= 1/\sqrt{3}$. For each data set, the causal function is sampled as follows. Let $q_{\min}$ and $q_{\max}$ be the $5\%$- and $95\%$ quantiles of $X$. We first generate a function $\tilde{f}$ that linearly extrapolates outside $[q_{\min}, q_{\max}]$ as described in Section \[sec:exp\_sampling\]. For a given threshold $\kappa$, we then draw $k_1, k_2 {\overset{\text{iid}}{\sim}}\mathcal{U}(-\kappa, \kappa)$ and construct $f$ for every $x \in {\mathbb{R}}$ by $$f(x) = \tilde{f}(x) + \frac{1}{2} k_1 ((x-q_{\min})_{-})^2 + \frac{1}{2} k_2 ((x-q_{\max})_{+})^2,$$ such that the curvature of $f$ on $(-\infty, q_{\min}]$ and $[q_{\max}, \infty)$ is $k_1$ and $k_2$, respectively. Figure \[fig:violation\_lin\_extrap\] shows results for $\kappa = 0,1,2,3,4$. As the curvature increases, the ability to generalize decreases.

![Worst-case mean squared error for increasingly strong violations of the linear extrapolation assumption. The grey area marks the inner 90 % quantile range of $X$ in the training distribution. As the curvature of $f$ outside the domain of the observed data increases, it becomes difficult to predict the interventional behavior of $Y$ for strong interventions. However, even in situations where the linear extrapolation assumption is strongly violated, it remains beneficial to extrapolate linearly. []{data-label="fig:violation_lin_extrap"}](varying_extrapolation_curvature_with_variability){width="\linewidth"}

Proofs {#app:proofs}
======

Proof of Proposition \[prop:minimax\_equal\_causal\] {#sec:prop:minimax_equal_causal}
----------------------------------------------------

Assume that $\cI$ is a set of interventions on $X$ with at least one confounding-removing intervention. Let $i \in {\mathcal{I}}$ and ${f_{\diamond}}\in {\mathcal{F}}$, then we have the following expansion $$\label{eq:decomp_equation}
    {\mathbb{E}}_{M(i)}[(Y - {f_{\diamond}}(X))^2]
    = {\mathbb{E}}_{M(i)}[(f(X) - {f_{\diamond}}(X))^2]+{\mathbb{E}}_{M(i)}[\xi_Y^2]+2{\mathbb{E}}_{M(i)}[\xi_Y(f(X)-{f_{\diamond}}(X))],$$ where $\xi_Y=h_1(H, {\varepsilon}_Y)$. For any intervention $i\in \cI$ the causal function always yields an identical loss. In particular, it holds that $$\begin{aligned}
 \label{eq:SupLossCausalFunction}
    \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}[(Y - f(X))^2] =\sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}[\xi_Y^2] = {\mathbb{E}}_{M}[\xi_Y^2],
  \end{aligned}$$ where we used that the distribution of $\xi_Y$ is not affected by an intervention on $X$. The loss of the causal function can never be better than the minimax loss, that is, $$\begin{aligned}
    \label{eq:lowerbdd_prop1}
    \inf_{{f_{\diamond}}\in{\mathcal{F}}}\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y - {f_{\diamond}}(X))^2]
    &\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y - f(X))^2] ={\mathbb{E}}_{M}[\xi_Y^2].
  \end{aligned}$$ In other words, the minimax solution (if it exists) is always better than or equal to the causal function. We will now show that when $\cI$ contains at least one confounding-removing intervention, then the minimax loss is dominated by any such intervention.

Fix $i_0\in \cI$ to be a confounding-removing intervention and let $(X,Y,H,A)$ be generated by the SCM $M(i_0)$. Recall that there exists a map $\psi^{i_0}$ such that $X:= \psi^{i_0}(g, h_2, A, H, {{\varepsilon}}_X ,I^{i_0})$ and that $X{\perp \!\!\! \perp}H$ as $i_0$ is a confounding-removing intervention. Furthermore, since the vectors $A$, $H$, ${{\varepsilon}}_X$, ${{\varepsilon}}_Y$ and $I^{i_0}$ are mutually independent, we have that $(X,H){\perp \!\!\! \perp}{{\varepsilon}}_Y$ which together with $X{\perp \!\!\! \perp}H$ implies $X, H$ and ${{\varepsilon}}_Y$ are mutually independent, and hence $X {\perp \!\!\! \perp}h_1(H,{{\varepsilon}}_Y)$. Using this independence we get that ${\mathbb{E}}_{M(i_0)}[\xi_Y(f(X)-{f_{\diamond}}(X))]={\mathbb{E}}_{M}[\xi_Y]{\mathbb{E}}_{M(i_0)}[(f(X)-{f_{\diamond}}(X))]$. Hence, for the intervention $i_0$ together with the modeling assumption ${\mathbb{E}}_{M}[\xi_Y]=0$ implies that for all ${f_{\diamond}}\in {\mathcal{F}}$, $${\mathbb{E}}_{M(i_0)}[(Y - {f_{\diamond}}(X))^2]
    = {\mathbb{E}}_{M(i_0)}[(f(X) - {f_{\diamond}}(X))^2]+{\mathbb{E}}_{M}[\xi_Y^2]\geq{\mathbb{E}}_{M}[\xi_Y^2].$$ This proves that the smallest loss at a confounding-removing intervention is achieved by the causal function. Denoting the non-empty subset of confounding-removing interventions by $\cI_{\text{cr}}\subseteq \cI$, this implies $$\begin{aligned}
  \label{eq:upperbdd_prop1}
  \inf_{{f_{\diamond}}\in{\mathcal{F}}} \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y - {f_{\diamond}}(X))^2] & \geq 
  \inf_{{f_{\diamond}}\in{\mathcal{F}}} \sup_{i\in \cI_{\text{cr}}} {\mathbb{E}}_{M(i)}[(Y - {f_{\diamond}}(X))^2] \\ \notag
  &\geq 
  \inf_{{f_{\diamond}}\in{\mathcal{F}}} {\mathbb{E}}_{M(i_0)}[(Y - {f_{\diamond}}(X))^2] \\ \notag
  &=  {\mathbb{E}}_{M}[\xi_Y^2].
  \end{aligned}$$ Combining and it immediately follows that $$\begin{aligned}
    \inf_{{f_{\diamond}}\in{\mathcal{F}}} \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y - {f_{\diamond}}(X))^2] =  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}[(Y - f(X))^2],
  \end{aligned}$$ and hence $$f\in\operatorname*{argmin}_{{f_{\diamond}}\in {\mathcal{F}}} \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{M(i)}[(Y-{f_{\diamond}}(X))^2],$$ which completes the proof of Proposition \[prop:minimax\_equal\_causal\].

Proof of Proposition \[prop:shift\_interventions\]
--------------------------------------------------

Let ${\mathcal{F}}$ be the class of all linear functions and let ${\mathcal{I}}$ denote the set of interventions on $X$ that satisfy $$\sup_{i\in \cI} \lambda_{\min}\big({\mathbb{E}}_{M(i)}\big[XX^\top\big]\big) =\infty.$$ We claim that the causal function $f(x)=b^\top x$ is the unique minimax solution of . We prove the result by contradiction. Let $\bar{f}\in\mathcal{F}$ (with $\bar{f}(x)=\bar{b}^{\top}x$) be such that $$\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-\bar{b}^\top X)^2]\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-b^\top X)^2],$$ and assume that $\norm{\bar{b}-b}_2>0$. For a fixed $i\in \cI$, we get the following bound $$\begin{aligned}
    {\mathbb{E}}_{M(i)}[(b^\top X-\bar{b}^\top X)^2]
    =(b-\bar{b})^\top {\mathbb{E}}_{M(i)}[XX^{\top}](b-\bar{b})
    \geq \lambda_{\min}({\mathbb{E}}_{M(i)}[XX^{\top}]) \|b-\bar{b} \|_2^2.
  \end{aligned}$$ Since we assumed that the minimal eigenvalue is unbounded, this means that we can choose $i\in{\mathcal{I}}$ such that ${\mathbb{E}}_{M(i)}[(b^\top X-\bar{b}^\top X)^2]$ can be arbitrarily large. However, applying Proposition \[prop:difference\_to\_causal\_function\], this leads to a contradiction since $\sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{M(i)}[(b^\top X-\bar{b}^\top X)^2]\leq 4\operatorname{Var}_{M}(\xi_Y)$ cannot be satisfied. Therefore, it must holds that $\bar{b}=b$, which moreover implies that $f$ is indeed a solution to the minimax problem $\operatorname*{argmin}_{{f_{\diamond}}\in\mathcal{F}}\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-{f_{\diamond}}(X))^2]$, as it achieves the lowest possible objective value. This completes the proof of Proposition \[prop:shift\_interventions\].

Proof of Proposition \[prop:difference\_to\_causal\_function\] {#sec:prop:difference_to_causal_function}
--------------------------------------------------------------

Let ${\mathcal{I}}$ be a set of interventions on $X$ or $A$ and let ${f_{\diamond}}\in{\mathcal{F}}$ with $$\label{eq:better_than_causal_cond}
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-{f_{\diamond}}(X))^2]\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}[(Y-f(X))^2].$$ For any $i\in{\mathcal{I}}$, the Cauchy-Schwartz inequality implies that $$\begin{aligned}
    &{\mathbb{E}}_{M(i)}[(Y-{f_{\diamond}}(X))^2]
    ={\mathbb{E}}_{M(i)}[(f(X)+\xi_Y-{f_{\diamond}}(X))^2]\\
    &\qquad={\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]+{\mathbb{E}}_{M(i)}[\xi_Y^2]+2{\mathbb{E}}_{M(i)}[\xi_Y(f(X)-{f_{\diamond}}(X))]\\
    &\qquad\geq{\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]+{\mathbb{E}}_{M}[\xi_Y^2]-2\left({\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]{\mathbb{E}}_{M}[\xi_Y^2]\right)^{\frac{1}{2}}.
  \end{aligned}$$ A similar computation shows that the causal function $f$ satisfies $${\mathbb{E}}_{M(i)}[(Y-f(X))^2]={\mathbb{E}}_{M}[\xi_Y^2].$$ So by condition this implies for any $i\in{\mathcal{I}}$ that $${\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]+{\mathbb{E}}_{M}[\xi_Y^2]-2\left({\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]{\mathbb{E}}_{M}[\xi_Y^2]\right)^{\frac{1}{2}}\leq{\mathbb{E}}_{M}[\xi_Y^2],$$ which is equivalent to $$\begin{aligned}
  &{\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]\leq 2\sqrt{{\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]{\mathbb{E}}_{M}[\xi_Y^2]}\\
  &\iff\quad
    {\mathbb{E}}_{M(i)}[(f(X)-{f_{\diamond}}(X))^2]\leq 4{\mathbb{E}}_{M}[\xi_Y^2].
  \end{aligned}$$ As this inequality holds for all $i\in{\mathcal{I}}$, we can take the supremum over all $i\in{\mathcal{I}}$, which completes the proof of Proposition \[prop:difference\_to\_causal\_function\].

Proof of Proposition \[prop:misspecification\_minimax\]
-------------------------------------------------------

As argued before, we have that for all $i \in {\mathcal{I}}_1$, $${\mathbb{E}}_{M(i)}\big[(Y-f(X))^2\big]={\mathbb{E}}_{M(i)}\big[\xi_Y^2\big]={\mathbb{E}}_{M}\big[\xi_Y^2\big].$$ Let now $f_1^*\in\mathcal{F}$ be a minimax solution w.r.t. ${\mathcal{I}}_1$. Then, using that the causal function $f$ lies in ${\mathcal{F}}$, it holds that $$\sup_{i\in{\mathcal{I}}_1}{\mathbb{E}}_{M(i)}\big[(Y-f_1^*(X))^2\big] \leq \sup_{i\in{\mathcal{I}}_1}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2\big] \ = {\mathbb{E}}_{M}\big[\xi_Y^2\big].$$ Moreover, if ${\mathcal{I}}_2\subseteq{\mathcal{I}}_1$, then it must also hold that $$\sup_{i\in{\mathcal{I}}_2}{\mathbb{E}}_{M(i)}\big[(Y-f_1^*(X))^2\big] \leq {\mathbb{E}}_{M}\big[\xi_Y^2\big]=\sup_{i\in{\mathcal{I}}_2}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2\big].$$ To prove the second part, we give a one-dimensional example. Let $\mathcal{F}$ be linear (i.e., $f(x)=b x$) and let ${\mathcal{I}}_1$ consist of shift interventions on $X$ of the form $$X^i\coloneqq g(A^i) + h_2(H^i, {\varepsilon}_X^i)+ c,$$ with $c\in [0, K]$. Then, the minimax solution $f^*_1$ (where $f^*_1(x)=b^*_1 x$) with respect to ${\mathcal{I}}_1$ is not equal to the causal function $f$ as long as ${\operatorname{Cov}}(X, \xi_Y)$ is strictly positive. This can be seen by explicitly computing the OLS estimator for a fixed shift $c$ and observing that the worst-case loss is attained at $c=K$. Now let $\cI_2$ be a set of interventions of the same form as $\cI_1$ but including shifts with $c>K$ such that $\cI_2 \not \subseteq \cI_1$. Since $\cF$ consists of linear functions, we know that the loss ${\mathbb{E}}_{M(i)}\big[(Y-f_1^*(X))^2\big]$ can become arbitrarily large, since $$\begin{aligned}
      &{\mathbb{E}}_{M(i)}\big[(Y-f_1^*(X))^2\big]\\
           &\quad=(b-b^*_1)^2{\mathbb{E}}_{M(i)}[X^2]+{\mathbb{E}}_{M}[\xi_Y^2]+2(b-b^*_1){\mathbb{E}}_{M(i)}[\xi_Y X]\\
      &\quad=(b-b^*_1)^2(c^2+{\mathbb{E}}_{M}[X^2]+2c{\mathbb{E}}_{M}[X])+{\mathbb{E}}_{M}[\xi_Y^2]+2(b-b^*_1)({\mathbb{E}}_{M}[\xi_Y X]+{\mathbb{E}}_{M}[\xi_Y]c),
    \end{aligned}$$ and $(b-b^*)^2>0$. In contrast, the loss for the causal function is always ${\mathbb{E}}_{M}[\xi_Y^2]$, so the worst-case loss of $f^*_1$ becomes arbitrarily worse than that of $f$. This completes the proof of Proposition \[prop:misspecification\_minimax\].

Proof of Proposition \[prop:suff\_general\]
-------------------------------------------

Let ${\varepsilon}> 0$. By definition of the infimum, we can find $f^* \in {\mathcal{F}}$ such that $$\left \vert \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
      \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \right \vert \leq {\varepsilon}.$$ Let now $\tilde{M} \in {\mathcal{M}}$ be s.t. ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$. By assumption, the left-hand side of the above inequality is unaffected by substituting $M$ for $\tilde{M}$, and the result thus follows.

Proof of Proposition \[prop:genX\_intra\]
-----------------------------------------

Let ${\mathcal{I}}$ be a well-behaved set of interventions on $X$. We consider two cases; (A) all interventions in ${\mathcal{I}}$ are confounding-preserving and (B) there is at least one intervention in ${\mathcal{I}}$ that is confounding-removing.

**Case (A):** In this case, we prove the result in two steps: (i) We show that $(A, \xi_X, \xi_Y)$ is identified from the observational distribution ${\mathbb{P}}_M$. (ii) We show that this implies that the intervention distributions $(X^i, Y^i)$, $i \in {\mathcal{I}}$, are also identified from the observational distribution, and conclude by using Proposition \[prop:suff\_general\]. Some of the details will be slightly technical because we allow for a large class of distributions (e.g., there is no assumption on the existence of densities).

We begin with step (i). In this case, ${\mathcal{I}}$ is a set of confounding-preserving interventions on $X$, and we have that ${\mathrm{supp}}_{{\mathcal{I}}}(X)\subseteq{\mathrm{supp}}(X)$. Fix $\tilde{M} =(\tilde{f},\tilde{g},\tilde{h}_1,\tilde{h}_2,\tilde{Q})
  \in {\mathcal{M}}$ such that ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$ and let $(\tilde{X}, \tilde{Y},\tilde{H},\tilde{A})$ be generated by the SCM of $\tilde{M}$. We have that $(X,Y,A) {\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{X}, \tilde{Y},\tilde{A})$ and by Assumption \[ass:identify\_f\], we have that $f \equiv \tilde{f}$ on ${\mathrm{supp}}(X)$, hence $f(X) {\ensuremath{\stackrel{\text{a.s.}}{=}}}\tilde{f}(X)$. Further, fix any $B\in \cB({\mathbb{R}}^p)$ (i.e., in the Borel sigma-algebra on ${\mathbb{R}}^p$) and note that $$\begin{aligned}
    \bE_M[\mathbbm{1}_{B}(A)X|A]
    &=\bE_M[\mathbbm{1}_{B}(A)g(A)+\mathbbm{1}_{B}(A)h_2(H,{{\varepsilon}}_X)|A] \\
    &= \bE_M[\mathbbm{1}_{B}(A)g(A)|A] + \mathbbm{1}_{B}(A)\bE[h_2(H,{{\varepsilon}}_X)] = \mathbbm{1}_{B}(A)g(A),
  \end{aligned}$$ almost surely. Here, we have used our modeling assumption ${\mathbb{E}}[h_2(H, {\varepsilon}_X)] = 0$. Hence, by similar arguments for $\bE_{\tilde{M}}(\mathbbm{1}_{B}(\tilde{A})\tilde{X}|\tilde{A})$ and the fact that $(X,Y,A) {\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{X}, \tilde{Y},\tilde{A})$ we have that $$\begin{aligned}
    \mathbbm{1}_{B}(A)g(A) {\ensuremath{\stackrel{\text{a.s.}}{=}}}\bE_M (\mathbbm{1}_{B}(A)X|A) {\ensuremath{\stackrel{\text{d}}{=}}}\bE_{\tilde{M}}(\mathbbm{1}_{B}(\tilde{A})\tilde{X}|\tilde{A}) {\ensuremath{\stackrel{\text{a.s.}}{=}}}\mathbbm{1}_{B}(\tilde{A})\tilde{g}(\tilde{A}).
  \end{aligned}$$ We conclude that $\mathbbm{1}_{B}(A)g(A){\ensuremath{\stackrel{\text{d}}{=}}}\mathbbm{1}_{B}(\tilde{A})\tilde{g}(\tilde{A})$ for any $B\in \cB({\mathbb{R}}^p)$. Let $\bP$ and $\tilde{\bP}$ denote the respective background probability measures on which the random elements $(X,Y,H,A)$ and $(\tilde{X},\tilde{Y},\tilde{H},\tilde{A})$ are defined. Fix any $F\in \sigma(A)$ (i.e., in the sigma-algebra generated by $A$) and note that there exists a $B\in \cB({\mathbb{R}}^p)$ such that $F=\{A\in B\}$. Since $A {\ensuremath{\stackrel{\text{d}}{=}}}\tilde{A}$, we have that, $$\begin{aligned}
    \int_F g(A) \, \mathrm{d} \bP = \int \mathbbm{1}_{B}(A) g(A) \, \mathrm{d} \bP = \int \mathbbm{1}_{B}(\tilde{A})\tilde{g}(\tilde{A})  \, \mathrm{d} \tilde{\bP} = \int  \mathbbm{1}_{B}(A)\tilde{g}(A)  \, \mathrm{d}\bP = \int_F \tilde{g}(A)  \, \mathrm{d}\bP.
  \end{aligned}$$ Both $g(A)$ and $\tilde{g}(A)$ are $\sigma(A)$-measurable and they agree integral-wise over every set $F\in \sigma(A)$, so we must have that $g(A) {\ensuremath{\stackrel{\text{a.s.}}{=}}}\tilde{g}(A)$. With $\eta(a,b,c)= (a,c-\tilde{f}(b),b-\tilde{g}(a))$ we have that $$\begin{aligned}
    (A,\xi_Y,\xi_X)   {\ensuremath{\stackrel{\text{a.s.}}{=}}}(A,Y-\tilde{f}(X),X-\tilde{g}(A)) 
                    =\eta(A,X,Y) 
                    {\ensuremath{\stackrel{\text{d}}{=}}}\eta(\tilde{A},\tilde{X},\tilde{Y}) 
                    = (\tilde{A},\tilde{\xi}_Y,\tilde{\xi}_X),
  \end{aligned}$$ so $(A,\xi_Y,\xi_X)
  {\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{A},\tilde{\xi}_Y,\tilde{\xi}_X)$. This completes step (i).

Next, we proceed with step (ii). Take an arbitrary intervention $i \in {\mathcal{I}}$ and let ${\varphi}^i, I^i, \tilde{I}^i$ with $I^i {\ensuremath{\stackrel{\text{d}}{=}}}\tilde{I}^i$, $I^i{\perp \!\!\! \perp}({{\varepsilon}}_X^i,{{\varepsilon}}_Y^i,{{\varepsilon}}_H^i,{{\varepsilon}}_A^i) \sim Q$ and $\tilde{I}^i {\perp \!\!\! \perp}(\tilde{{{\varepsilon}}}^i_X,\tilde{{{\varepsilon}}}^i_Y,\tilde{{{\varepsilon}}}^i_H,\tilde{{{\varepsilon}}}^i_A)
  \sim \tilde{Q}$ be such that the structural assignments for $X^i$ and $\tilde{X}^i$ in $M(i)$ and $\tilde{M}(i)$, respectively, are given as $$X^i := {\varphi}^i(A^i,g(A^i), h_2(H^i, {\varepsilon}_X^i), I^i) \quad \text{ and } \quad \tilde{X}^i := {\varphi}^i(\tilde{A}^i,\tilde{g}(\tilde{A}^i), \tilde{h}_2(\tilde{H}^i, \tilde{{\varepsilon}}_X^i), \tilde{I}^i).$$ Define $\xi_X^i := h_2(H^i, {\varepsilon}_X^i)$, $\xi_Y^i := h_1(H^i, {\varepsilon}_Y^i)$, $\tilde{\xi}_X^i := \tilde{h}_2(\tilde{H}^i, \tilde{{\varepsilon}}_X^i)$ and $\tilde{\xi}_Y^i := \tilde{h}_1(\tilde{H}^i,
  \tilde{{\varepsilon}}_Y^i)$. Then, it holds that $$(A^i, \xi_X^i, \xi_Y^i) {\ensuremath{\stackrel{\text{d}}{=}}}(A,\xi_X, \xi_Y)  {\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{A},  \tilde{\xi}_X, \tilde{\xi}_Y)  {\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{A}^i, \tilde{\xi}_X^i, \tilde{\xi}_Y^i),$$ where we used step (i), that $(A^i, \xi_X^i, \xi_Y^i)$ and $(A,\xi_X, \xi_Y)$ are generated by identical functions of the noise innovations and that $ ({{\varepsilon}}_X,{{\varepsilon}}_Y,{{\varepsilon}}_H,{{\varepsilon}}_A) $ and $({{\varepsilon}}_X^i,{{\varepsilon}}_Y^i,{{\varepsilon}}_H^i,{{\varepsilon}}_A^i)$ have identical distributions. Adding a random variable with the same distribution, that is mutually independent with all other variables, on both sides does not change the distribution of the bundle, hence $$\begin{aligned}
    (A^i, \xi_X^i, \xi_Y^i,I^i)  {\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{A}^i, \tilde{\xi}_X^i, \tilde{\xi}_Y^i, \tilde{I}^i).
  \end{aligned}$$ Define $\kappa(a,b,c,d) :=
  ({\varphi}^i(a,\tilde{g}(a),b,d),\tilde{f}({\varphi}^i(a,\tilde{g}(a),b,d))+c)$. As shown in step (i) above, we have that $g(A^i){\ensuremath{\stackrel{\text{a.s.}}{=}}}\tilde{g}(A^i)$. Furthermore, since ${\mathrm{supp}}(X^i) \subseteq{\mathrm{supp}}(X)$ we have that $f(X^i) {\ensuremath{\stackrel{\text{a.s.}}{=}}}\tilde{f}(X^i)$, and hence $$\begin{aligned}
    (X^i, Y^i) &{\ensuremath{\stackrel{\text{a.s.}}{=}}}(X^i,\tilde{f}(X^i)+\xi_Y^i) \\
    &=  ({\varphi}^i(A^i,g(A^i), \xi_X^i , I^i) ,\tilde{f}({\varphi}^i(A^i,g(A^i), \xi_X^i, I^i))+\xi_Y^i) \\
    &{\ensuremath{\stackrel{\text{a.s.}}{=}}}({\varphi}^i(A^i,\tilde{g}(A^i), \xi_X^i , I^i) ,\tilde{f}({\varphi}^i(A^i,\tilde{g}(A^i), \xi_X^i, I^i))+\xi_Y^i) \\
    &=\kappa(A^i, \xi_X^i, \xi_Y^i,I^i) 
                {\ensuremath{\stackrel{\text{d}}{=}}}\kappa(\tilde{A}^i, \tilde{\xi}_X^i, \tilde{\xi}_Y^i, \tilde{I}^i) 
               = (\tilde{X}^i, \tilde{Y}^i).
  \end{aligned}$$ Thus, $\bP_{M(i)}^{(X,Y)} = \bP_{\tilde{M}(i)}^{(X,Y)}$, which completes step (ii). Since $i \in {\mathcal{I}}$ was arbitrary, the result now follows from Proposition \[prop:suff\_general\].

**Case (B):** Assume that the set of interventions ${\mathcal{I}}$ contains at least one confounding-removing intervention. Let $\tilde{M} =(\tilde{f},\tilde{g},\tilde{h}_1,\tilde{h}_2,\tilde{Q})
  \in {\mathcal{M}}$ be such that ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$. Then, by Proposition \[prop:minimax\_equal\_causal\], it follows that the causal function $\tilde{f}$ is a minimax solution w.r.t. $(\tilde{M}, {\mathcal{I}})$. By Assumption \[ass:identify\_f\], we further have that $\tilde{f}$ and $f$ coincide on ${\mathrm{supp}}(X) \supseteq {\mathrm{supp}}_{{\mathcal{I}}}(X)$. Hence, it follows that $$\inf_{{f_{\diamond}}\in {\mathcal{F}}} \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}[(Y - {f_{\diamond}}(X))^2] = \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}[(Y - \tilde{f}(X))^2] = \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}[(Y - f(X))^2],$$ showing that also $f$ is a minimax solution w.r.t. $(\tilde{M}, {\mathcal{I}})$. This completes the proof of Proposition \[prop:genX\_intra\].

Proof of Proposition \[prop:impossibility\_interpolation\]
----------------------------------------------------------

We first show that the causal parameter $\beta$ is not a minimax solution. Let $u := \sup {\mathcal{I}}< \infty$, since ${\mathcal{I}}$ is bounded, and take $b = \beta + 1/(\sigma u)$. By an explicit computation we get that $$\begin{aligned}
 \inf_{{b_{\diamond}}\in {\mathbb{R}}}\,
 \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[(Y-{b_{\diamond}}X)^2\big] & \leq  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[(Y-b X)^2\big]
 = \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[({{\varepsilon}}_Y + \tfrac{1}{\sigma}H - 
 \tfrac{1}{\sigma u}iH)^2\big] \\
 &= \sup_{i\in{\mathcal{I}}} \left[1 + \left(1 - \tfrac{i}{u}\right)^2\right] 
 < 2
 = \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-\beta X)^2\big],
  \end{aligned}$$ where the last inequality holds because $0 < 1 + (1 - i/u)^2 < 2$ for all $i \in {\mathcal{I}}$, and since ${\mathcal{I}}{\subseteq}{\mathbb{R}}_{>0}$ is compact with upper bound $u$. Hence, $$\begin{aligned}
  \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-\beta X)^2\big]
  - \inf_{{b_{\diamond}}\in {\mathbb{R}}}\,
  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[(Y-{b_{\diamond}}X)^2\big] > 0,
  \end{aligned}$$ proving that the causal parameter is not a minimax solution for model $M$ w.r.t. $(\cF, \cI)$. Recall that in order to prove that $(\bP_{M},\cM)$ does not generalize with respect to $\cI$ we have to show that there exists an ${{\varepsilon}}>0$ such that for all $b\in {\mathbb{R}}$ it holds that $$\begin{aligned}
\sup_{\tilde M: \bP_{\tilde{M}}= \bP_M}\big|   \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-b X)^2\big]
- \inf_{{b_{\diamond}}\in {\mathbb{R}}}\,
\sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{b_{\diamond}}X)^2\big]\big| \geq {{\varepsilon}}.\end{aligned}$$ Thus, it remains to show that for all $b \not = \beta$ there exists a model $\tilde M \in \cM$ with $\bP_M = \bP_{\tilde M}$ such that the generalization loss is bounded below uniformly by a positive constant. We will show the stronger statement that for any $b \neq \beta$, there exists a model $\tilde{M}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, such that under $\tilde{M}$, $b$ results in arbitrarily large generalization error. Let $c > 0$ and $i_0 \in {\mathcal{I}}$. Define $$\begin{aligned}
  \tilde\sigma \coloneqq \frac{\operatorname{sign}{((\beta - b)i_0)}\sqrt{1 + c} - 1}{
  (\beta - b)i_0} > 0,
  \end{aligned}$$ and let $\tilde M \coloneqq M(\gamma, \beta, \tilde{\sigma}, Q)$. By construction of the model class ${\mathcal{M}}$, it holds that ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_{M}$. Furthermore, by an explicit computation we get that $$\begin{aligned}
\label{eq:proof_imp_a1}
  \begin{split}
  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-bX)^2\big] 
  &\geq {\mathbb{E}}_{\tilde{M}(i_0)}\big[(Y-bX)^2\big]
  = {\mathbb{E}}_{\tilde{M}(i_0)}\big[((\beta-b)i_0H+{{\varepsilon}}_Y+
  \tfrac{1}{\tilde\sigma}H)^2\big]\\
  &={\mathbb{E}}_{\tilde{M}(i_0)}\big[([(\beta-b)i_0\tilde{\sigma}+1] {{\varepsilon}}_H+{{\varepsilon}}_Y )^2\big] 
  = [(\beta-b)i_0\tilde{\sigma}+1]^2 + 1 \\
  &= ((\beta - b)i_0\tilde\sigma)^2 + 2(\beta - b)i_0\tilde\sigma +
  2\\
  &= (\operatorname{sign}{((\beta-b)i_0)}\sqrt{1+c}-1)^2 + 2\operatorname{sign}{((\beta-b)i_0)}\sqrt{1+c}\\
  &= c + 2.
  \end{split}
  \end{aligned}$$ Finally, by definition of the infimum, it holds that $$\begin{aligned}
\label{eq:proof_imp_b1}
    \begin{split}
      \inf_{{b_{\diamond}}\in {\mathbb{R}}}\,\sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{b_{\diamond}}X)^2\big] \leq \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-\beta
      X)^2\big] = 2.
    \end{split}
  \end{aligned}$$ Combining  and  yields that the generalization error is bounded below by $c$. That is, $$\begin{aligned}
  \big|   \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-b X)^2\big]
  - \inf_{{b_{\diamond}}\in {\mathbb{R}}}\,
  \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{b_{\diamond}}X)^2\big]\big| \geq c.
  \end{aligned}$$ The above results make no assumptions on $\gamma$, and hold true, in particular, if $\gamma \neq 0$ (in which case Assumption \[ass:identify\_f\] is satisfied, see Appendix \[sec:IVconditions\]). This completes the proof of Proposition \[prop:impossibility\_interpolation\].

Proof of Proposition \[prop:genX\_extra\]
-----------------------------------------

Let $\tilde{M} \in {\mathcal{M}}$ be such that ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$. By Assumptions \[ass:identify\_f\] and \[ass:gen\_f\], it holds that $f \equiv \tilde{f}$. The proof now proceeds analogously to that of Proposition \[prop:genX\_intra\].

Proof of Proposition \[prop:extrapolation\_bounded\_deriv\_cr\] {#sec:prop:extrapolation_bounded_deriv_cr}
---------------------------------------------------------------

By Assumption \[ass:identify\_f\], $f$ is identified on ${\mathrm{supp}}^{M}(X)$ by the observational distribution $\bP_{M}$. Let $\cI$ be a set of interventions containing at least one confounding-removing intervention. For any $\tilde{M}=(\tilde{f},\tilde{g},\tilde{h}_1,\tilde{h}_2,\tilde{Q})\in
  \cM$, Proposition \[prop:minimax\_equal\_causal\] yields that the causal function is a minimax solution. That is, $$\begin{aligned}
 \notag
    \inf_{{f_{\diamond}}\in\mathcal{F}}\sup_{i\in\cI}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]
    &= \sup_{i\in\cI}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-\tilde{f}(X))^2\big] = \sup_{i\in \cI }{\mathbb{E}}_{\tilde{M}(i)}[\xi_Y^2] \\
    &={\mathbb{E}}_{\tilde{M}}[\xi_Y^2], \label{eq:propboundedderiv_causalfunctionsolvesminimax}
  \end{aligned}$$ where we used that any intervention $i\in \cI$ does not affect the distribution of $\xi_Y=\tilde{h}_2(H,{{\varepsilon}}_Y)$. Now, assume that $\tilde{M}=(\tilde{f},\tilde{g},\tilde{h}_1,\tilde{h}_2,\tilde{Q})\in
  \cM$ satisfies $\bP_{\tilde{M}} = \bP_M$. Since $({\mathbb{P}}_M,\cM)$ satisfies Assumption \[ass:identify\_f\], we have that $f \equiv \tilde{f}$ on ${\mathrm{supp}}^M(X)={\mathrm{supp}}^{\tilde{M}}(X)$. Let $f^*$ be any function in ${\mathcal{F}}$ such that $f^*=f$ on ${\mathrm{supp}}^M(X)$. We first show that $\norm{\tilde{f} - f^*}_{\cI,\infty} \leq 2\delta K$, where $\|f\|_{\cI,\infty} := \sup_{x\in{\mathrm{supp}}_{\cI}^M(X)}\|f(x)\|$. By the mean value theorem, for all ${f_{\diamond}}\in{\mathcal{F}}$ it holds that $\abs{{f_{\diamond}}(x) - {f_{\diamond}}(y)} \leq K\norm{x - y}$, for all $x, y\in \mathcal D$. For any $x \in {\mathrm{supp}}^M_\cI(X)$ and $y \in {\mathrm{supp}}^M(X)$ we have $$\begin{aligned}
    \abs[\big]{\tilde{f}(x) - f^*(x)} 
    &=  \abs[\big]{\tilde{f}(x) - \tilde{f}(y) + f^*(y) - f^*(x)}\\
    &\leq  \abs[\big]{\tilde{f}(x) - \tilde{f}(y)} + \abs[\big]{f^*(y) - f^*(x)}\\
    &\leq 2 K \norm{x - y},
  \end{aligned}$$ where we used the fact that $\tilde{f}(y)= f(y) =f^*(y)$, for all $y\in {\mathrm{supp}}^M(X)$. In particular, it holds that $$\begin{aligned}
\label{eq:unif_norm}
    \begin{split}
      \norm{\tilde{f} - f^*}_{\cI,\infty} 
      &=\sup_{x\in {\mathrm{supp}}^M_\cI(X)} \abs[\big]{\tilde{f}(x) - f^*(x)}\\
      &\leq  2K\sup_{x\in {\mathrm{supp}}^M_\cI(X)} \inf_{y\in{\mathrm{supp}}^M(X)}\norm{x - y}\\
      &= 2\delta K.
    \end{split}
  \end{aligned}$$ For any $i\in{\mathcal{I}}$ we have that $$\begin{aligned}
    {\mathbb{E}}_{\tilde{M}(i)}\big[(Y- f^*(X))^2\big]
    &={\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)+\xi_Y-f^*(X))^2\big]\nonumber\\
    &={\mathbb{E}}_{\tilde{M}}\big[\xi_Y^2\big] +
      {\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f^*(X))^2\big]\nonumber\\
    &\qquad\qquad + 2{\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y(\tilde{f}(X)-f^*(X))\big].\label{eq:starting_pt_split}
  \end{aligned}$$ Next, we can use Cauchy-Schwarz, and in to get that $$\begin{aligned}
    &\left|\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big]-\inf_{{f_{\diamond}}\in{\mathcal{F}}}\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]\nonumber \right|\\
    &\quad= \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big]-{\mathbb{E}}_{\tilde{M}}[\xi_Y^2]\nonumber \\
    &\quad = \sup_{i\in{\mathcal{I}}} \left(
      {\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f^*(X))^2\big]+2{\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y(\tilde{f}(X)-f^*(X))\big] \right) \nonumber\\
    &\quad\leq 4\delta^2K^2+4\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)},\label{eq:ineq_prt}
  \end{aligned}$$ proving the first statement. Finally, if ${\mathcal{I}}$ consists only of confounding-removing interventions, then the bound in can be improved by using that ${\mathbb{E}}[\xi_Y]=0$ together with $H{\perp \!\!\! \perp}X$. In that case, we get that ${\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y(\tilde{f}(X)-f(X))\big]=0$ and hence the bound becomes $4\delta^2 K^2$. This completes the proof of Proposition \[prop:extrapolation\_bounded\_deriv\_cr\].

Proof of Proposition \[prop:extrapolation\_bounded\_deriv\] {#sec:prop:extrapolation_bounded_deriv}
-----------------------------------------------------------

By Assumption \[ass:identify\_f\], $f$ is identified on ${\mathrm{supp}}^{M}(X)$ by the observational distribution $\bP_{M}$. Let $\cI$ be a set of confounding-preserving interventions. For a fixed ${\varepsilon}>0$, let $f^*\in\mathcal{F}$ be a function satisfying $$\label{eq:fstarDiffEpsilon}
    \abs{\sup_{i\in\cI}{\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2)\big]-\inf_{{f_{\diamond}}\in\mathcal{F}}\sup_{i\in\cI}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2)\big]}\leq{\varepsilon}.$$ Fix any secondary model $\tilde{M}=(\tilde{f},\tilde{g},\tilde{h}_1,\tilde{h}_2,\tilde{Q})\in
        \cM$ with $\bP_{\tilde{M}} = \bP_M$. The general idea is to derive an upper bound for $\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}[(Y-f^*(X))^2]$ and a lower bound for $\inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
        \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}[(Y-{f_{\diamond}}(X))^2]$ which will allow us to bound the absolute difference of interest.

Since $({\mathbb{P}}_M,\cM)$ satisfies Assumption \[ass:identify\_f\], we have that $f \equiv \tilde{f}$ on ${\mathrm{supp}}^M(X)={\mathrm{supp}}^{\tilde{M}}(X)$. We first show that $\norm{\tilde{f} - f}_{\cI,\infty} \leq 2\delta K$, where $\|f\|_{\cI,\infty} := \sup_{x\in{\mathrm{supp}}_{\cI}^M(X)}\|f(x)\|$. By the mean value theorem, for all ${f_{\diamond}}\in{\mathcal{F}}$ it holds that $\abs{{f_{\diamond}}(x) - {f_{\diamond}}(y)} \leq K\norm{x - y}$, for all $x, y\in \mathcal D$. For any $x \in {\mathrm{supp}}^M_\cI(X)$ and $y \in {\mathrm{supp}}^M(X)$ we have $$\begin{aligned}
    \abs[\big]{\tilde{f}(x) - f(x)} 
    &=  \abs[\big]{\tilde{f}(x) - \tilde{f}(y) + f(y) - f(x)}\\
    &\leq  \abs[\big]{\tilde{f}(x) - \tilde{f}(y)} + \abs[\big]{f(y) - f(x)}\\
    &\leq 2 K \norm{x - y},
    \end{aligned}$$ where we used the fact that $\tilde{f}(y)=f(y)$, for all $y\in {\mathrm{supp}}_M(X)$. In particular, it holds that $$\begin{aligned}
\label{eq:unif_norm2}
    \begin{split}
    \norm{\tilde{f} - f}_{\cI,\infty} 
    &= \sup_{x\in {\mathrm{supp}}^M_\cI(X)} \abs[\big]{\tilde{f}(x) - f(x)}\\
    &\leq  2K\sup_{x\in {\mathrm{supp}}^M_\cI(X)} \inf_{y\in{\mathrm{supp}}^M(X)}\norm{x - y}\\
    &= 2\delta K.
    \end{split}
    \end{aligned}$$ Let now $i\in{\mathcal{I}}$ be fixed. The term $\xi_Y = h_1(H, {\varepsilon}_Y)$ is not affected by the intervention $i$. Furthermore, ${\mathbb{P}}^{(X,\xi_Y)}_{M(i)}={\mathbb{P}}^{(X,\xi_Y)}_{\tilde{M}(i)}$ since $i$ is confounding-preserving (this can be seen by a slight modification to the arguments from case (A) in the proof of Proposition \[prop:genX\_intra\]). Thus, for any ${f_{\diamond}}\in \cF$ we have that $$\begin{aligned}
    &{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \nonumber\\
    &\qquad={\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)+\xi_Y-{f_{\diamond}}(X)+f(X)-f(X))^2\big]\nonumber\\
    &\qquad={\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y^2\big] +
    {\mathbb{E}}_{\tilde{M}(i)}\big[(f(X)-{f_{\diamond}}(X))^2\big]+
    {\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f(X))^2\big]\nonumber\\
    &\qquad\qquad\qquad + 2{\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y(f(X)-{f_{\diamond}}(X))\big]\nonumber\\
    &\qquad\qquad\qquad +
    2{\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f(X))(f(X)-{f_{\diamond}}(X))\big]\nonumber\\
    &\qquad\qquad\qquad +
    2{\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y(\tilde{f}(X)-f(X))\big]\nonumber\\
            &\qquad={\mathbb{E}}_{M(i)}\big[\xi_Y^2\big] +
    {\mathbb{E}}_{M(i)}\big[(f(X)-{f_{\diamond}}(X))^2\big]+
    {\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))^2\big]\nonumber\\
    &\qquad\qquad\qquad +2{\mathbb{E}}_{M(i)}\big[\xi_Y(f(X)-{f_{\diamond}}(X))\big]\nonumber\\
    &\qquad\qquad\qquad +
    2{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))(f(X)-{f_{\diamond}}(X))\big]\nonumber\\
    &\qquad\qquad\qquad + 2{\mathbb{E}}_{M(i)}\big[\xi_Y(\tilde{f}(X)-f(X))\big]\nonumber\\
    &\qquad={\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big] + L_1^i(\tilde{f}) +
    L_2^i(\tilde{f},{f_{\diamond}}) + L_3^i(\tilde{f}),\label{eq:decomposition_intoLs}
    \end{aligned}$$ where, we have made the following definitions $$\begin{aligned}
    L_1^i(\tilde{f})&\coloneqq {\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))^2\big],\\
    L_2^i(\tilde{f},{f_{\diamond}})&\coloneqq
    2{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))(f(X)-{f_{\diamond}}(X))\big],\\
    L_3^i(\tilde{f})&\coloneqq 2{\mathbb{E}}_{M(i)}\big[\xi_Y(\tilde{f}(X)-f(X))\big].
    \end{aligned}$$ Using it follows that $$\label{eq:estimate_L1}
    0\leq L_1^i(\tilde{f}) \leq 4\delta^2 K^2,$$ and by the Cauchy-Schwarz inequality it follows that $$\begin{aligned}
    \abs[\big]{L_3^i(\tilde{f})} &\leq 2 \sqrt{{\operatorname{Var}}_M(\xi_Y)4\delta^2K^2} =4\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)}.  \label{eq:estimate_L3}
    \end{aligned}$$ Let now ${f_{\diamond}}\in{\mathcal{F}}$ be any function such that $$\label{eq:ConditionBetterThanTildeCausal}
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde M(i)}\big[(Y-{f_{\diamond}}(X))^2)\big]\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde M(i)}\big[(Y-\tilde{f}(X))^2)\big],$$ then by , the Cauchy-Schwarz inequality and Proposition \[prop:difference\_to\_causal\_function\], it holds for all $i\in \cI$ that $$\begin{aligned}
 \nonumber
    L_2^i(\tilde{f},{f_{\diamond}}) &= 2{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))(f(X)-{f_{\diamond}}(X))\big]\\ 
    &=2{\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f(X))(f(X)-{f_{\diamond}}(X))\big]\nonumber\\
    &= -2{\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f(X))^2\big]+2{\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f(X))(\tilde{f}(X)-{f_{\diamond}}(X))\big]\nonumber\\
    &\geq -8\delta^2K^2 -
    2\sqrt{4\delta^2K^2}\sqrt{4{\operatorname{Var}}_M(\xi_Y)}\nonumber\\
    &=-8\delta^2K^2- 8\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)}, \label{eq:estimate_L2_BetterThanTildeCausal}
    \end{aligned}$$ where, in the third equality, we have added and subtracted the term $2{\mathbb{E}}_{\tilde{M}(i)}\big[(\tilde{f}(X)-f(X))\tilde{f}(X)\big]$. Now let $\cS := \{{f_{\diamond}}\in \cF :
        \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \leq
        \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-\tilde{f}(X))^2\big]
        \}$ be the set of all functions satisfying . Due to , , and we have the following lower bound of interest $$\begin{aligned}
 \notag
    &\inf_{{f_{\diamond}}\in{\mathcal{F}}}\, \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \\\notag
    =& \inf_{{f_{\diamond}}\in\cS}\,   \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \\\notag
    =& \inf_{{f_{\diamond}}\in \cS}\,
    \sup_{i\in{\mathcal{I}}}\big\{ {\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big] + L_1^i(\tilde{f}) +
    L_2^i(\tilde{f},{f_{\diamond}}) + L_3^i(\tilde{f}) \big\} \\\notag
    \geq& \inf_{{f_{\diamond}}\in \cS}\,
    \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]  -8\delta^2K^2- 8\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} - 4\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} \\ \label{eq:lowerboundoninf}
    \geq & \inf_{{f_{\diamond}}\in \cF}\,
    \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]  - 8 \delta^2 K^2 - 12 \delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)}.
    \end{aligned}$$ Next, we construct the aforementioned upper bound of interest. To that end, note that $$\begin{aligned}
 \notag
    &\sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big] \\ \label{eq:QuantityForUpperBound}
    &\quad = \sup_{i\in{\mathcal{I}}} \left\{ {\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big] + L_1^i(\tilde{f}) +
    L_2^i(\tilde{f},f^*) + L_3^i(\tilde{f}) \right\},
    \end{aligned}$$ by . We have already established upper bounds for $L_1^i(\tilde f)$ and $L_3^i(\tilde{f})$ in and , respectively. In order to control $L_2^i(\tilde{f},f^*)$ we introduce an auxiliary function. Let $\bar{f}^*\in \cF$ satisfy $$\label{eq:BarfStarBetterThanCausal}
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2)\big]\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2)\big],$$ and $$\label{eq:epsilonineqforBarfStar}
    \abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
        \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]}\leq {\varepsilon}.$$ Choosing such a $\bar f^*\in \cF$ is always possible. If $f$ is an ${{\varepsilon}}$-minimax solution, i.e., it satisfies , then choose $\bar
        f^*=f$. Otherwise, if $f$ is not a ${{\varepsilon}}$-minimax solution, then choose any $\bar f^*\in \cF$ that is an ${{\varepsilon}}$-minimax solution (which is always possible). In this case we have that $$\begin{aligned}
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \leq {{\varepsilon}},
    \end{aligned}$$ and $$\begin{aligned}
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \geq {{\varepsilon}},
    \end{aligned}$$ which implies that is satisfied. We can now construct an upper bound on $L_2^i(\tilde{f},f^*)$ in terms of $L_2^i(\tilde{f},\bar f^*)$ by noting that for all $i\in\cI$ $$\begin{aligned}
 \notag
    \abs[\big]{L_2^i(\tilde{f},f^*)}
    &= 2\big|{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))(f(X)-f^*(X))\big] \big|  \\
    \notag
    &\leq 2\big|{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))(f(X)-\bar{f}^*(X))\big]\big| \\ \notag
    &\qquad \qquad +2{\mathbb{E}}_{M(i)}\big|(\tilde{f}(X)-f(X))(\bar{f}^*(X)-f^*(X))\big| \\ \notag
    &= \abs[\big]{L_2^i(\tilde{f},\bar f^*)} +2{\mathbb{E}}_{M(i)}\big|(\tilde{f}(X)-f(X))(\bar{f}^*(X)-f^*(X))\big| \\  \notag
    &\leq \abs[\big]{L_2^i(\tilde{f},\bar f^*)} +2\sqrt{{\mathbb{E}}_{M(i)}\left[(\tilde{f}(X)-f(X))^2\right]{\mathbb{E}}_{M(i)}\left[(\bar{f}^*(X)-f^*(X))^2 \right]} \\ \label{eq:UpperboundOfL2TildefStarf}
    &\leq \abs[\big]{L_2^i(\tilde{f},\bar f^*)} +4\delta K\sqrt{{\mathbb{E}}_{M(i)}\left[(\bar{f}^*(X)-f^*(X))^2 \right]},
    \end{aligned}$$ where we used the triangle inequality, Cauchy-Schwarz inequality and . Furthermore, and together with Proposition \[prop:difference\_to\_causal\_function\] yield the following bound $$\begin{aligned}
    |L_2^i(\tilde{f},\bar f^*)|
    &=2\big|{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))(f(X)-\bar f^*(X))\big] \big|\nonumber\\
    &= 2 \sqrt{{\mathbb{E}}_{M(i)}\big[(\tilde{f}(X)-f(X))^2\big] {\mathbb{E}}_{M(i)}\big[(f(X)-\bar f^*(X))^2\big]}\nonumber\\
    &\leq  2\sqrt{4\delta^2K^2} \sqrt{4{\operatorname{Var}}_M(\xi_Y)}\nonumber \\
    &=  8\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)}, \label{eq:estimate_L2}
    \end{aligned}$$ for any $i\in\cI$. Thus, it suffices to construct an upper bound on the second term in the final expression in . Direct computation leads to $$\begin{aligned}
    {\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big] &={\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2\big] \\
    &\quad \quad+{\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big] \\
    &\quad \quad +2{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))(\bar{f}^*(X)-f^*(X))\big].
    \end{aligned}$$ Rearranging the terms and applying the triangle inequality and Cauchy-Schwarz results in $$\begin{aligned}
    &{\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big] \\
    &\quad ={\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big]-{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2\big]\\
    &\quad \quad \quad -2{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))(\bar{f}^*(X)-f^*(X))\big]\\
    &\quad\leq \big|{\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big]- \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]\big| \\
    &\quad \quad \quad+ \big| \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]-{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2\big] \big| \\
    & \quad \quad\quad + 2{\mathbb{E}}_{M(i)}\big|(Y-\bar{f}^*(X))(\bar{f}^*(X)-f^*(X))\big|\\
    &\quad\leq 2{\varepsilon}+2\sqrt{{\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2\big]}\sqrt{{\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big]}\\
    &\quad\leq 2{\varepsilon}+2\sqrt{{\operatorname{Var}}_{M}(\xi_Y)}\sqrt{{\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big]},
    \end{aligned}$$ for any $i\in \cI$. Here, we used that both $f^*$ and $\bar f^*$ are ${{\varepsilon}}$-minimax solutions with respect to $M$ and that $\bar f^*$ satisfies which implies that $$\begin{aligned}
    {\mathbb{E}}_{M(i)}\big[(Y-\bar{f}^*(X))^2)\big]\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f(X))^2)\big] = \sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{M(i)}\big[\xi_Y^2\big] = {\operatorname{Var}}_{M}(\xi_Y),
    \end{aligned}$$ for any $i\in \cI$, as $\xi_Y$ is unaffected by an intervention on $X$. Thus, ${\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big]$ must satisfy $\ell({\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big])\leq 0$, where $\ell:[0,\infty)\to {\mathbb{R}}$ is given by $\ell(z)=z-2{{\varepsilon}}-2\sqrt{{\operatorname{Var}}_{M}(\xi_Y)} \sqrt{z}$. The linear term of $\ell$ grows faster than the square root term, so the largest allowed value of ${\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big]$ coincides with the largest root of $\ell(z)$. The largest root is given by $$\begin{aligned}
    C^2:=2{{\varepsilon}}+2{\operatorname{Var}}_{M}(\xi_Y) + 2\sqrt{{\operatorname{Var}}_{M}(\xi_Y)^2+2{{\varepsilon}}{\operatorname{Var}}_{M}(\xi_Y)},
    \end{aligned}$$ where $(\cdot)^2$ refers to the square of $C$. Hence, for any $i\in \cI$ it holds that $$\label{eq:LastFactorUpperBound}
    {\mathbb{E}}_{M(i)}\big[(\bar{f}^*(X)-f^*(X))^2\big]\leq C^2.$$ Hence by , and we have that the following upper bound is valid for any $i\in \cI$. $$\begin{aligned}
    \abs[\big]{L_2^i(\tilde{f},f^*)} &\leq  8\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} + 4\delta K C \label{eq:estimate_L2_fstar}.
    \end{aligned}$$ Thus, using with , and , we get the following upper bound $$\begin{aligned}
 \notag
    &\sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big] \\ 
            &\quad  \leq \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big] + 4\delta^2K^2+4\delta K C+ 12\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)}.\label{eq:upperboundsup}
    \end{aligned}$$ Finally, by combining the bounds and together with we get that $$\begin{aligned}
    &\abs[\Big]{\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f^*(X))^2\big]
        -
        \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
        \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]}\nonumber\\ \nonumber 
    &\quad\leq\sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-f^*(X))^2\big]
    -
    \inf_{{f_{\diamond}}\in{\mathcal{F}}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{M(i)}\big[(Y-{f_{\diamond}}(X))^2\big]\nonumber\\ \nonumber
    &\qquad\qquad+  4\delta^2K^2+4\delta K C+ 12\delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)}\\ \nonumber
    &\qquad \qquad  + 8 \delta^2 K^2 + 12 \delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} \\
    & \quad \leq {{\varepsilon}}+ 12 \delta^2 K^2 + 24 \delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} + 4 \delta K C. \label{eq:CorrectUpperBound}
    \end{aligned}$$ Using that all terms are positive, we get that $$C = \sqrt{{\operatorname{Var}}_{M}(\xi_Y)} + \sqrt{{\operatorname{Var}}_{M}(\xi_Y) + 2 {{\varepsilon}}}
    \leq 2\sqrt{{\operatorname{Var}}_{M}(\xi_Y)} + \sqrt{2 {{\varepsilon}}}$$ Hence, is bounded above by $$\begin{aligned}
          {{\varepsilon}}+ 12 \delta^2 K^2 + 32 \delta K \sqrt{{\operatorname{Var}}_M(\xi_Y)} +
            4 \sqrt{2} \delta K \sqrt{{{\varepsilon}}}.
    \end{aligned}$$ This completes the proof of Proposition \[prop:extrapolation\_bounded\_deriv\].

Proof of Proposition \[prop:impossibility\_extrapolation\]
----------------------------------------------------------

Let $\bar{f}\in {\mathcal{F}}$ and $c > 0$. By assumption, ${\mathcal{I}}$ is a well-behaved set of support-extending interventions on $X$. Since ${\mathrm{supp}}_{{\mathcal{I}}}^{M}(X) \setminus {\mathrm{supp}}^{M}(X)$ has non-empty interior, there exists an intervention $i_0 \in {\mathcal{I}}$ and ${{\varepsilon}}>0$ such that ${\mathbb{P}}_{M(i_0)}(X\in B) \geq {\varepsilon}$, for some open subset $B \subsetneq \bar{B}$, such that $\mathrm{dist}(B,{\mathbb{R}}^d \setminus \bar{B})>0$, where $\bar{B}:={\mathrm{supp}}_{{\mathcal{I}}}^{M}(X) \setminus {\mathrm{supp}}^{M}(X)$. Let $\tilde{f}$ be any continuous function satisfying that, for all $x \in B\cup({\mathbb{R}}^d\setminus\bar B)$, $$\begin{aligned}
    \tilde{f}(x) =
    \begin{cases}
      \bar{f}(x) + \gamma, \quad & x \in B\\
      f(x), \quad & x \in {\mathbb{R}}^d \setminus \bar{B},
    \end{cases}
  \end{aligned}$$ where $\gamma := {\varepsilon}^{-1/2} \left\{(2 {\mathbb{E}}_{\tilde{M}}[\xi_{Y}^2] +
    c)^{1/2} + ({\mathbb{E}}_{\tilde{M}}[\xi_{Y}^2])^{1/2}\right\}$. Consider a secondary model $\tilde{M} = (\tilde{f}, g, h_1, h_2, Q)\in\cM$. Then, by Assumption \[ass:identify\_f\], it holds that ${\mathbb{P}}_{M} = {\mathbb{P}}_{\tilde{M}}$. Since ${\mathcal{I}}$ only consists of interventions on $X$, it holds that ${\mathbb{P}}_{M(i_0)}(X\in B) = {\mathbb{P}}_{\tilde{M}(i_0)}(X\in B)$ (this holds since all components of $\tilde{M}$ and $M$ are equal, except for the function $f$, which is not allowed to enter in the intervention on $X$). Therefore, $$\begin{aligned}
\label{eq:proof_imp_a}
  {\mathbb{E}}_{\tilde{M}(i_0)}\big[(Y-\bar{f}(X))^2\big]
  &\geq   {\mathbb{E}}_{\tilde{M}(i_0)}\big[(Y-\bar{f}(X))^2 \mathbbm{1}_{B}(X)\big] \nonumber\\
  &= {\mathbb{E}}_{\tilde{M}(i_0)}\big[(\gamma + \xi_Y)^2 \mathbbm{1}_{B}(X)\big] \nonumber\\
  &\geq   \gamma^2{\varepsilon}+ 2 \gamma {\mathbb{E}}_{\tilde{M}(i_0)}\big[\xi_Y
  \mathbbm{1}_{B}(X) \big]\nonumber\\
  &\geq  \gamma^2{\varepsilon}- 2 \gamma \left({\mathbb{E}}_{\tilde{M}}\big[\xi_Y^2\big]
  {\varepsilon}\right)^{1/2} \nonumber\\
  &=  c + {\mathbb{E}}_{\tilde{M}}[\xi_{Y}^2],
  \end{aligned}$$ where the third inequality follows from Cauchy–Schwarz. Further, by the definition of the infimum it holds that $$\begin{aligned}
\label{eq:proof_imp_b}
  \begin{split}
  \inf_{{f_{\diamond}}\in {\mathcal{F}}}\,\sup_{i\in{\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]
  \leq \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-\tilde{f}(X))^2\big] 
  = {\mathbb{E}}_{\tilde{M}}[\xi_{Y}^2].
  \end{split}
  \end{aligned}$$ Therefore, combining  and , the claim follows.

Proof of Proposition \[prop:genA\]
----------------------------------

We prove the result by showing that under Assumption \[ass:identify\_g\] it is possible to express interventions on $A$ as confounding-preserving interventions on $X$ and applying Propositions \[prop:genX\_intra\] and \[prop:genX\_extra\]. To avoid confusion, we will throughout this proof denote the true model by $M^0 = (f^0, g^0, h_1^0, h_2^0, Q^0)$. Fix an intervention $i\in{\mathcal{I}}$. Since it is an intervention on $A$, there exist $\psi^i$ and $I^i$ such that for any $M = (f,g,h_1, h_2, Q) \in {\mathcal{M}}$, the intervened SCM $M(i)$ is of the form $$A^i := \psi^i(I^i, {{\varepsilon}}_A^i), \quad H^i := {{\varepsilon}}_H^i, \quad X^i := g(A^i) + h_2(H^i, {{\varepsilon}}_X^i), \quad Y^i:= f(X^i) + h_1(H^i,{{\varepsilon}}_Y^i),$$ where $({{\varepsilon}}^i_X, {{\varepsilon}}^i_Y, {{\varepsilon}}^i_A, {{\varepsilon}}^i_H) \sim Q$. We now define a confounding-preserving intervention $j$ on $X$, such that, for all models $\tilde{M}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, the distribution of $(X,Y)$ under $\tilde{M}(j)$ coincides with that under $\tilde{M}(i)$. To that end, define the intervention function $$\bar{\psi}^j(h_2, A^j, H^j, {{\varepsilon}}^j_X ,I^j) \coloneqq g^0(\psi^i(I^j, A^j)) +
    h_2(H^j, {{\varepsilon}}_X^j),$$ where $g^0$ is the fixed function corresponding to model $M$, and therefore not an argument of $\bar{\psi}^j$. Let now $j$ be the intervention on $X$ satisfying that, for all $M = (f,g,h_1, h_2, Q) \in {\mathcal{M}}$, the intervened model $M(j)$ is given as $$A^{j} := {{\varepsilon}}_A^{j}, \quad H^{j} := {{\varepsilon}}_H^{j}, \quad X^{j} := \bar{\psi}^j(h_2, A^{j}, H^{j}, {{\varepsilon}}^{j}_X ,I^{j}), \quad Y^{j}:= f(X^{j}) + h_1(H^{j},{{\varepsilon}}_Y^{j}),$$ where $({{\varepsilon}}^j_X, {{\varepsilon}}^j_Y, {{\varepsilon}}^j_A, {{\varepsilon}}^j_H) \sim Q$ and where $I^j$ is chosen such that $I^j{\ensuremath{\stackrel{\text{d}}{=}}}I^i$. By definition, $j$ is a confounding-preserving intervention. Let now $\tilde{M} = (\tilde{f}, \tilde{g}, \tilde{h}_1, \tilde{h}_2,
  \tilde{Q})$ be such that ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, and let $(\tilde{X}^i, \tilde{Y}^i)$ and $(\tilde{X}^j, \tilde{Y}^j)$ be generated under $\tilde{M}(i)$ and $\tilde{M}(j)$, respectively. By Assumption \[ass:identify\_g\], it holds for all $a \in {\mathrm{supp}}(A) \cup {\mathrm{supp}}_{{\mathcal{I}}}(A)$ that $\tilde{g}(a) = g^0(a)$. Hence, we get that $$\begin{aligned}
    (\tilde{X}^i, \tilde{Y}^i) 
    &{\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{g}(\psi^i(I^i, \tilde{{{\varepsilon}}}_A^i)) + \tilde{h}_2(\tilde{{{\varepsilon}}}^i_H, \tilde{{{\varepsilon}}}_X^i), 
      \tilde{f}(\tilde{g}(\psi^i(I^i, \tilde{{{\varepsilon}}}_A^i)) + \tilde{h}_2(\tilde{{{\varepsilon}}}^i_H, \tilde{{{\varepsilon}}}_X^i)) + \tilde{h}_1(\tilde{{{\varepsilon}}}_H^i, \tilde{{{\varepsilon}}}^i_Y)) \\
    &=
      (g^0(\psi^i(I^i, \tilde{{{\varepsilon}}}_A^i)) + \tilde{h}_2(\tilde{{{\varepsilon}}}^i_H, \tilde{{{\varepsilon}}}_X^i), 
      \tilde{f}(g^0(\psi^i(I^i, \tilde{{{\varepsilon}}}_A^i)) + \tilde{h}_2(\tilde{{{\varepsilon}}}^i_H, \tilde{{{\varepsilon}}}_X^i)) + \tilde{h}_1(\tilde{{{\varepsilon}}}_H^i, \tilde{{{\varepsilon}}}^i_Y)) \\
    &{\ensuremath{\stackrel{\text{d}}{=}}}(g^0(\psi^i(I^j, \tilde{{{\varepsilon}}}_A^j)) + \tilde{h}_2(\tilde{{{\varepsilon}}}^j_H, \tilde{{{\varepsilon}}}_X^j), 
      \tilde{f}(g^0(\psi^i(I^j, \tilde{{{\varepsilon}}}_A^j)) + \tilde{h}_2(\tilde{{{\varepsilon}}}^j_H, \tilde{{{\varepsilon}}}_X^j)) + \tilde{h}_1(\tilde{{{\varepsilon}}}_H^j, \tilde{{{\varepsilon}}}^j_Y)) \\
    &{\ensuremath{\stackrel{\text{d}}{=}}}(\bar{\psi}^j(\tilde{h}_2, \tilde{{{\varepsilon}}}_A^j, \tilde{{{\varepsilon}}}_H^j, \tilde{{{\varepsilon}}}^j_X ,I^j), 
      \tilde{f}(\bar{\psi}^j(\tilde{h}_2, \tilde{{{\varepsilon}}}_A^j, \tilde{{{\varepsilon}}}_H^j, \tilde{{{\varepsilon}}}^j_X ,I^j)) + \tilde{h}_1(\tilde{{{\varepsilon}}}_H^j, \tilde{{{\varepsilon}}}^j_Y)) \\
    &{\ensuremath{\stackrel{\text{d}}{=}}}(\tilde{X}^j, \tilde{Y}^j),
  \end{aligned}$$ as desired. Since $i \in {\mathcal{I}}$ was arbitrary, we have now shown that there exists a mapping $\pi$ from ${\mathcal{I}}$ into a set $\mathcal{J}$ of confounding-preserving (and hence a well-behaved set) of interventions on $X$, such that for all $\tilde{M}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$, ${\mathbb{P}}^{(X,Y)}_{\tilde{M}(i)} = {\mathbb{P}}_{\tilde{M}(\pi(i))}^{(X,Y)}$. Hence, we can rewrite Equation  in Definition \[defi:general\] in terms of the set $\mathcal{J}$. The result now follows from Propositions \[prop:genX\_intra\] and \[prop:genX\_extra\].

Proof of Proposition \[prop:impossibility\_intA\]
-------------------------------------------------

Let $b \in {\mathbb{R}}^d$ be such that $f(x) = b^\top x$ for all $x \in {\mathbb{R}}^d$. We start by characterizing the error ${\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big]$. Let us consider models of the form $\tilde{M} = (f, \tilde{g}, h_1, h_2, Q) \in {\mathcal{M}}$ for some function $\tilde{g} \in {\mathcal{G}}$ with $\tilde{g}(a) = g(a)$ for all $a \in {\mathrm{supp}}_M(A)$. Clearly, any such model satisfies that ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$. For every $a \in \mathcal{A}$, let $i_a \in {\mathcal{I}}$ denote the corresponding hard intervention on $A$. For every $a \in \mathcal{A}$ and ${b_{\diamond}}\in {\mathbb{R}}^d$, we then have $$\label{eq:linear_mse}
\begin{aligned}
&\ {\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y - {b_{\diamond}}^\top X)^2\big]  \\
= &\ {\mathbb{E}}_{\tilde{M}(i_a)}\big[(b^\top X + \xi_Y - {b_{\diamond}}^\top X)^2\big] \\
                                                                    = &\ (b - {b_{\diamond}})^\top {\mathbb{E}}_{\tilde{M}(i_a)}[X X^\top] (b - {b_{\diamond}}) + 2(b - {b_{\diamond}})^\top {\mathbb{E}}_{\tilde{M}(i_a)}[X \xi_Y] + {\mathbb{E}}_{\tilde{M}(i_a)}\big[\xi_Y^2] \\
                                                                = &\ (b - {b_{\diamond}})^\top \underbrace{(\tilde{g}(a) \tilde{g}(a)^\top + {\mathbb{E}}_{M}[\xi_X \xi_X^\top])}_{=: K_{\tilde{M}}(a)} (b - {b_{\diamond}}) + 2(b - {b_{\diamond}})^\top {\mathbb{E}}_{M}[\xi_X \xi_Y] + {\mathbb{E}}_{M}\big[\xi_Y^2],
\end{aligned}$$ where we have used that, under $i_a$, the distribution of $(\xi_X, \xi_Y)$ is unaffected. We now show that, for any $\tilde{M}$ with the above form, the causal function $f$ does not minimize the worst-case mean squared error across interventions in ${\mathcal{I}}$. The idea is to show that the worst-case mean squared error  strictly decreases at ${b_{\diamond}}= b$ in the direction $u := {\mathbb{E}}_{M}[\xi_X \xi_Y] / \norm{{\mathbb{E}}_{M}[\xi_X \xi_Y]}_2$. For every $a \in \mathcal{A}$ and $s \in {\mathbb{R}}$, define $$\ell_{\tilde{M},a}(s) := {\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y-(b + s u)^\top X)^2\big] = u^\top K_{\tilde{M}}(a) u \cdot s^2  - 2 u^\top {\mathbb{E}}_{M}[\xi_X \xi_Y] \cdot s + {\mathbb{E}}_{M}\big[\xi_Y^2].$$ For every $a$, $\ell_{\tilde{M},a}^\prime (0) = -2 \norm{{\mathbb{E}}_{M}[\xi_X \xi_Y]}_2 < 0$, showing that $\ell_{\tilde{M},a}$ is strictly decreasing at $s=0$ (with a derivative that is bounded away from 0 across all $a \in \mathcal{A}$). By boundedness of $\mathcal{A}$ and by the continuity of $a \mapsto \ell_{\tilde{M},a}^{\prime \prime } (0)  = 2 u^\top K_{\tilde{M}}(a) u$, it further follows that $\sup_{a \in \mathcal{A}} \card{\ell^{\prime \prime}_{\tilde{M},a} (0)} < \infty$. Hence, we can find $s_0 > 0$ such that for all $a \in \mathcal{A}$, $\ell_{\tilde{M}, a}(0) > \ell_{\tilde{M}, a}(s_0)$. It now follows by continuity of $(a, s) \mapsto \ell_{\tilde{M},a}(s)$ that $$\sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-b^\top X)^2\big] = \sup_{a \in \mathcal{A}} \ell_{\tilde{M}, a}(0) > \sup_{a \in \mathcal{A}} \ell_{\tilde{M}, a}(s_0) = \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-(b+s_0 u)^\top X)^2\big],$$ showing that $b+s_0 u$ attains a lower worst-case mean squared error than $b$.

We now show that all functions other than $f$ may result in an arbitrarily large error. Let $\bar{b} \in {\mathbb{R}}^d \setminus \{ b \}$ be given, and let $j \in \{1, \dots, d\}$ be such that $b_j \neq \bar{b}_j$. The idea is to construct a function $\tilde{g} \in {\mathcal{G}}$ such that, under the corresponding model $\tilde{M} = (f, \tilde{g}, h_1, h_2, Q) \in {\mathcal{M}}$, some hard interventions on $A$ result in strong shifts of the $j$th coordinate of $X$. Let $a \in \mathcal{A}$. Let $e_j \in {\mathbb{R}}^d$ denote the $j$th unit vector, and assume that $\tilde{g}(a) = n e_j$ for some $n \in {\mathbb{N}}$. Using , it follows that $$\begin{aligned}
&{\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y-\bar{b}^\top X)^2\big] \\
&\qquad = n^2 (\bar{b}_j - b_j)^2 + (\bar{b} - b)^\top {\mathbb{E}}_{M}[\xi_X \xi_X^\top] (\bar{b} - b) + 
2
(\bar{b} - b)^\top {\mathbb{E}}_{M}[\xi_X \xi_Y] + {\mathbb{E}}_{M}\big[\xi_Y^2].\end{aligned}$$ By letting $n \to \infty$, we see that the above error may become arbitrarily large. Given any $c > 0$, we can therefore construct $\tilde{g}$ such that ${\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y-\bar{b}^\top X)^2\big] \geq c
+{\mathbb{E}}_{M}\big[\xi_Y^2]$. By carefully choosing $a \in \text{int}(\mathcal{A} \setminus {\mathrm{supp}}_M(A))$, this can be done such that $\tilde{g}$ is continuous and $\tilde{g}(a) = g(a)$ for all $a \in {\mathrm{supp}}_M(A)$, ensuring that $\bP_{\tilde M}= \bP_M$. It follows that $$\begin{aligned}
c   &\leq {\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y-{\bar{b}}^\top X)^2\big]  - {\mathbb{E}}_{M}\big[\xi_Y^2] \\
    &= {\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y-{\bar{b}}^\top X)^2\big]  -  \sup_{i \in {\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)} \big[(Y- b^\top X)^2] \\
    &\leq {\mathbb{E}}_{\tilde{M}(i_a)}\big[(Y-{\bar{b}}^\top X)^2\big]  - \inf_{{b_{\diamond}}\in {\mathbb{R}}^d} \sup_{i \in {\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)} \big[(Y-{b_{\diamond}}^\top X)^2]  \\
    &\leq \sup_{i \in {\mathcal{I}}} {\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{\bar{b}}^\top X)^2\big]  - \inf_{{b_{\diamond}}\in {\mathbb{R}}^d} \sup_{i \in {\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)} \big[(Y-{b_{\diamond}}^\top X)^2],\end{aligned}$$ which completes the proof of Proposition \[prop:impossibility\_intA\].

Proof of Proposition \[thm:consis\]
-----------------------------------

By assumption, $ \cI$ is a set of interventions on $X$ or $A$ of which at least one is confounding-removing. Now fix any $$\tilde M=(f_{\eta_0}(x;\tilde{\theta}),\tilde g,\tilde h_1,\tilde h_2,\tilde Q)\in \cM,$$ with $\bP_{M}=\bP_{\tilde M}$. By Proposition \[prop:minimax\_equal\_causal\], we have that a minimax solution is given by the causal function. That is, $$\begin{aligned}
\inf_{{f_{\diamond}}\in{\mathcal{F}}_{\eta_0}}\,
\sup_{i\in {\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] = \sup_{i\in {\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-
f_{\eta_0}(X;\tilde{\theta})
)^2\big] = {\mathbb{E}}_{M}[\xi^2_Y],\end{aligned}$$ where we used that $\xi_Y$ is unaffected by an intervention on $X$. By the support restriction ${\mathrm{supp}}^M(X) \subseteq (a,b)$ we know that $$\begin{aligned}
f_{\eta_0}(x;\theta^0)=B(x)^\top \theta^0, \quad 
f_{\eta_0}(x;\tilde \theta)=B(x)^\top \tilde \theta,  \quad 
f_{\eta_0}(x;\hat \theta_{\lambda^\star_n,\eta_0,\mu}^n)=B(x)^\top  \hat \theta_{\lambda^\star_n,\eta_0,\mu}^n,\end{aligned}$$ for all $x\in {\mathrm{supp}}^M(X)$. Furthermore, as $Y=B(X)^\top \theta^0+\xi_Y$ $\bP_{M}$-almost surely, we have that $$\begin{aligned}
 \label{eq:ExpansionOfECY_0}
{\mathbb{E}}_M\left[ C(A) Y \right] = {\mathbb{E}}_M\left[ C(A) B(X)^\top \theta^0\right] + {\mathbb{E}}_M\left[C(A) \xi_Y \right] = {\mathbb{E}}_M\left[ C(A) B(X)^\top \right]\theta^0,\end{aligned}$$ where we used the assumptions that ${\mathbb{E}}\left[ \xi_Y \right] =0$ and $A {\perp \!\!\! \perp}\xi_Y$ by the exogeneity of $A$. Similarly, $$\begin{aligned}
 {\mathbb{E}}_{\tilde M}\left[ C(A) Y \right] = {\mathbb{E}}_{\tilde M}\left[ C(A) B(X)^\top \right]\tilde\theta.\end{aligned}$$ As $\bP_M = \bP_{\tilde M}$, it holds $
{\mathbb{E}}_M[ C(A) Y ] = {\mathbb{E}}_{\tilde M}[ C(A) Y ]$ and ${\mathbb{E}}_M[ C(A) B(X)^\top ]= {\mathbb{E}}_{\tilde{M}}[ C(A) B(X)^\top ]$, hence $$\begin{aligned}
{\mathbb{E}}_{ M}\left[ C(A) B(X)^\top \right]\tilde\theta = {\mathbb{E}}_{M}\left[ C(A) B(X)^\top \right]\theta^0 \iff \tilde \theta = \theta^0,\end{aligned}$$ by assumption \[ass:RankCondition\], which states that $ {\mathbb{E}}[ C(A) B(X)^\top ]$ is of full rank (bijective). In other words, the causal function parameterized by $\theta^0$ is identified from the observational distribution. Assumptions \[ass:identify\_f\] and \[ass:gen\_f\] are therefore satisfied. Furthermore, we also have that $$\begin{aligned}
    &\sup_{i\in  {\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] \\
    &\quad =  \sup_{i\in\cI} \big\{ {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] +  {\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y^2\big] \\
    & \qquad \qquad  + 2 {\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y (f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu})) \big]  \big\} \\
    &\quad \leq  \sup_{i\in\cI} \big\{ {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] +  {\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y^2\big] \\
    & \qquad \qquad  + 2 \sqrt{ {\mathbb{E}}_{\tilde{M}(i)}\big[\xi_Y^2 \big] {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2 \big] } \big\} \\
    &\quad \leq  \sup_{i\in\cI}  {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] +  {\mathbb{E}}_{M}\big[\xi_Y^2\big] \\
    & \qquad \qquad  + 2 \sqrt{ {\mathbb{E}}_{M}\big[\xi_Y^2 \big] \sup_{i\in\cI}  {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2 \big] },
    \end{aligned}$$ by Cauchy-Schwarz inequality, where we additionally used that $ {\mathbb{E}}_{\tilde M(i)} [ \xi_Y^2 ] = {\mathbb{E}}_{M} [ \xi_Y^2 ]$ as $\xi_Y$ is unaffected by interventions on $X$. Thus, $$\begin{aligned}
    &\big\vert \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}_{\eta_0}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \big\vert \\
&\quad \leq  \sup_{i\in\cI}  {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big]  \\
& \qquad \qquad  + 2 \sqrt{ {\mathbb{E}}_{M}\big[\xi_Y^2 \big] \sup_{i\in\cI}  {\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2 \big] } .
    \end{aligned}$$ For the next few derivations let $\hat \theta =\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}$ for notational simplicity. Note that, for all $x \in {\mathbb{R}}$, $$\begin{aligned}
(f_{\eta_0}(x;\theta^0)-f_{\eta_0}(x;\hat \theta ))^2 &\leq (\theta^0 -\hat \theta)^\top B(x)B(x)^\top (\theta^0 -\hat \theta) \\
& \quad + (B(a)^\top (\theta^0-\hat \theta ) + B'(a)^\top (\theta^0-\hat \theta )(x-a))^2 \\
& \quad + (B(b)^\top (\theta^0-\hat \theta ) + B'(b)^\top (\theta^0-\hat \theta )(x-b))^2.\end{aligned}$$ The second term has the following upper bound $$\begin{aligned}
& (B(a)^\top (\theta^0-\hat \theta ) + B'(a)^\top (\theta^0-\hat \theta )(x-a))^2 \\
 &\quad  = (\theta^0 -\hat \theta)^\top B(a)B(a)^\top (\theta^0 -\hat \theta) \\
 & \qquad +  (x-a)^2(\theta^0 -\hat \theta)^\top B'(a)B'(a)^\top (\theta^0 -\hat \theta)  \\
 & \qquad + 2(x-a) (\theta^0 -\hat \theta)^\top B'(a)B(a)^\top (\theta^0 -\hat \theta) \\
 &\quad \leq \lambda_{\max}(B(a)B(a)^\top) \|\theta^0 - \hat \theta\|_2^2 \\
 & \qquad + (x-a)^2 \lambda_{\max}(B'(a)B'(a)^\top ) \|\theta^0 - \hat \theta\|_2^2 \\
 & \qquad + 2(x-a) \lambda_{\max}((B'(a)B(a)^\top+B(a)B'(a)^\top )/2 ) \|\theta^0 - \hat \theta\|_2^2,\end{aligned}$$ where $\lambda_{\max}$ denotes the maximum eigenvalue. An analogous upper bound can be constructed for the third term. Thus, by combining these two upper bounds with a similar upper bound for the first term, we arrive at $$\begin{aligned}
&{\mathbb{E}}_{\tilde{M}(i)} \big[(f_{\eta_0}(X;\theta^0)-f_{\eta_0}(X;\hat{\theta}))^2\big] \\
&\quad  \leq \lambda_{\max}({\mathbb{E}}_{\tilde{M}(i)} [B(X)B(X)^\top])\|\theta^0 -\hat \theta\|_2^2 \\
&\qquad + \lambda_{\max}(B(a)B(a)^\top) \|\theta^0 - \hat \theta\|_2^2 \\
& \qquad + {\mathbb{E}}_{\tilde{M}(i)} [(X-a)^2] \lambda_{\max}(B'(a)B'(a)^\top ) \|\theta^0 - \hat \theta\|_2^2 \\
& \qquad + 2{\mathbb{E}}_{\tilde{M}(i)} [X-a] \lambda_{\max}((B'(a)B(a)^\top+B(a)B'(a)^\top )/2 ) \|\theta^0 - \hat \theta\|_2^2\\
&\qquad + \lambda_{\max}(B(b)B(b)^\top) \|\theta^0 - \hat \theta\|_2^2 \\
& \qquad + {\mathbb{E}}_{\tilde{M}(i)} [(X-b)^2] \lambda_{\max}(B'(b)B'(b)^\top ) \|\theta^0 - \hat \theta\|_2^2 \\
& \qquad + 2{\mathbb{E}}_{\tilde{M}(i)} [X-b] \lambda_{\max}((B'(b)B(b)^\top+B(b)B'(b)^\top )/2  ) \|\theta^0 - \hat \theta\|_2^2.\end{aligned}$$ Assumption \[ass:MaximumEigenValueBounded\] imposes that $\sup_{i\in\cI}{\mathbb{E}}_{\tilde{M}(i)} [X^2]$ and $\sup_{i\in\cI}\lambda_{\max}({\mathbb{E}}_{\tilde{M}(i)} [B(X)B(X)^\top])$ are finite. Hence, the supremum of each of the above terms is finite. That is, there exists a constant $c>0$ such that $$\begin{aligned}
    &\big\vert \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-f_{\eta_0}(X;\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}))^2\big] - \inf_{{f_{\diamond}}\in{\mathcal{F}}_{\eta_0}}\,
    \sup_{i\in{\mathcal{I}}}{\mathbb{E}}_{\tilde{M}(i)}\big[(Y-{f_{\diamond}}(X))^2\big] \big\vert \\
&\quad \leq  c \|\theta^0 - \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} \|_2^2 + 2 \sqrt{ {\mathbb{E}}_{M}\big[\xi_Y^2 \big] c }\|\theta^0 - \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} \|_2 .
    \end{aligned}$$

It therefore suffices to show that $$\begin{aligned}
\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} \underset{n\to\infty}{\stackrel{P}{\longrightarrow}} \theta^0,\end{aligned}$$ with respect to the distribution induced by $M$. To simplify notation, we henceforth drop the $M$ subscript in the expectations and probabilities. Note that by the rank conditions in \[ass:RankCondition\], and the law of large numbers, we may assume that the corresponding sample product moments satisfy the same conditions. That is, for the purpose of the following arguments, it suffices that the sample product moment only satisfies these rank conditions asymptotically with probability one.

Let $B:= B(X)$, $C:= C(A)$, let $\fB$ and $\fC$ be row-wise stacked i.i.d. copies of $B(X)^\top$ and $C(A)^\top $, and recall the definition $\fP_\delta := \fC \left( \fC^\top \fC + \delta \fM \right)^{-1}
\fC^\top$. By convexity of the objective function we can find a closed form expression for our estimator of $\theta^0$ by solving the corresponding normal equations. The closed form expression is given by $$\begin{aligned}
\hat \theta^n_{\lambda, \eta, \mu}
: &= \operatorname*{argmin}_{\theta \in {\mathbb{R}}^{k}} 
\norm{{\mathbf{Y}} - {\mathbf{B}} \theta }_2^2 + \lambda \norm{{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2 + \gamma \theta^\top {\mathbf{K}} \theta, \\
&=  \left(  \frac{\fB^\top  \fB}{n}   + \lambda^\star_n  \frac{\fB^\top \fP_\delta  \fP_\delta   \fB}{n} + \frac{\gamma \fK}{n} \right)^{-1}  \left( \frac{\fB^\top \fY }{n} + \lambda^\star_n \frac{\fB^\top \fP_\delta  \fP_\delta   \fY}{n} \right),\end{aligned}$$ where we used that $\lambda^\star_n \in [0,\infty)$ almost surely by \[ass:LambdaStarAlmostSurelyFinite\]. Consequently (using standard convergence arguments and that $n^{-1} \gamma \fK $ and $n^{-1} \delta \fM$ converges to zero in probability), if $\lambda^\star_n$ diverges to infinity in probability as $n$ tends to infinity, then $$\begin{aligned}
\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} &\stackrel{P}{\to} \left( {\mathbb{E}}\left[ BC^\top   \right]{\mathbb{E}}\left[ CC^\top   \right]^{-1}{\mathbb{E}}\left[ CB^\top   \right] \right)^{-1} {\mathbb{E}}\left[ BC^\top \right] {\mathbb{E}}\left[ CC^\top   \right]^{-1} {\mathbb{E}}\left[ C Y   \right]  \\
&= \theta^0.\end{aligned}$$ Here, we also used that the terms multiplied by $\lambda^\star_n$ are the only asymptotically relevant terms. These are the standard arguments that the K-class estimator (with minor penalized regression modifications) is consistent as long as the parameter $\lambda^\star_n$ converges to infinity, or, equivalently, $\kappa_n^\star= \lambda^\star_n/(1+\lambda^\star_n)$ converges to one in probability.

We now consider two cases: *(i)* ${\mathbb{E}}[B\xi_Y]\not =0$ and *(ii)* ${\mathbb{E}}[B\xi_Y]=0$, corresponding to the case with unmeasured confounding and without, respectively. For *(i)* we show that $\lambda^\star_n$ converges to infinity in probability and for *(ii)* we show consistency by other means (as $\lambda^\star_n$ might not converge to infinity in this case).

**Case (i):** The confounded case ${\mathbb{E}}[B\xi_Y]\not =0$. It suffices to show that $$\lambda^\star_n := \inf\{\lambda\geq 0 : T_n(\hat \theta^n
_{\lambda,\eta_{0},\mu})\leq q(\alpha)\}\underset{n\to\infty}{\stackrel{P}{\longrightarrow}} \infty.$$ To that end, note that for fixed $\lambda \geq 0$ we have that $$\begin{aligned}
 \label{eq:ThetaLambdaConsistentEstimator}
\hat{\theta}^n_{\lambda,\eta_0,\mu} & \underset{n\to\infty}{\stackrel{P}{\longrightarrow}} \theta_\lambda,\end{aligned}$$ where $$\begin{aligned}
 \label{eq:ThetaLambdaFullRepresentation}
\theta_\lambda \, &:= \left( {\mathbb{E}}\left[ BB^\top \right] + \lambda {\mathbb{E}}\left[ BC^\top   \right]{\mathbb{E}}\left[ CC^\top   \right]^{-1}{\mathbb{E}}\left[ CB^\top   \right] \right)^{-1} \\ \notag
& \qquad \qquad \times \left( {\mathbb{E}}\left[ B Y  \right]  + \lambda {\mathbb{E}}\left[ BC^\top \right] {\mathbb{E}}\left[ CC^\top   \right]^{-1} {\mathbb{E}}\left[ C Y   \right] \right). \end{aligned}$$ Recall that states that $
{\mathbb{E}}\left[ C Y \right] = {\mathbb{E}}\left[ C B^\top \right]\theta^0$. Using and that $Y=B^\top \theta^0 + \xi_Y$ $\bP_{M}$-almost surely, we have that the latter factor of is given by $$\begin{aligned}
&{\mathbb{E}}\left[ B Y  \right]  + \lambda {\mathbb{E}}\left[ BC^\top \right] {\mathbb{E}}\left[ CC^\top   \right]^{-1} {\mathbb{E}}\left[ C Y   \right] \\
& \quad = {\mathbb{E}}\left[ B B^\top   \right]\theta^0 + {\mathbb{E}}\left[ B \xi_Y   \right] + \lambda {\mathbb{E}}\left[ BC^\top \right] {\mathbb{E}}\left[ CC^\top   \right]^{-1} {\mathbb{E}}\left[ C B^\top    \right] \theta^0 \\
&\quad = \left( {\mathbb{E}}\left[ B B^\top   \right] + \lambda {\mathbb{E}}\left[ BC^\top \right] {\mathbb{E}}\left[ CC^\top   \right]^{-1} {\mathbb{E}}\left[ C B^\top    \right] \right)\theta^0 +  {\mathbb{E}}\left[ B \xi_Y   \right]\end{aligned}$$ Inserting this into we arrive at the following representation of $\theta_\lambda$ $$\begin{aligned}
 \label{eq:ThetaLambdaInTermsOfTrueTheta}
\theta_\lambda  &= \theta^0 +  \left( {\mathbb{E}}\left[ BB^\top \right] + \lambda {\mathbb{E}}\left[ BC^\top   \right]{\mathbb{E}}\left[ CC^\top   \right]^{-1}{\mathbb{E}}\left[ CB^\top   \right] \right)^{-1}  {\mathbb{E}}\left[ B \xi_Y  \right].\end{aligned}$$ Since ${\mathbb{E}}\left[ B \xi_Y  \right]\not =0$ by assumption, the above yields that $$\begin{aligned}
 \label{eq:ThetaTrueNotEqualThetaLambda}
\forall \lambda \geq 0 : \quad \quad \theta^0 \not = \theta_\lambda.\end{aligned}$$ Now we prove that $\lambda^\star_n$ diverges to infinity in probability as $n$ tends to infinity. That is, for any $\lambda \geq 0$ we will prove that $$\begin{aligned}
\lim_{n \to \infty }\bP (\lambda^\star_n \leq \lambda ) =0. \end{aligned}$$ We fix an arbitrary $\lambda\geq 0$. By we have that $\theta^0 \not = \theta_{\lambda}$. This implies that there exists an ${{\varepsilon}}>0$ such that $\theta^0\not \in \overline{B(\theta_{\lambda},{{\varepsilon}})}$, where $\overline{B(\theta_{\lambda},{{\varepsilon}})}$ is the closed ball in ${\mathbb{R}}^k$ with center $\theta_{\lambda}$ and radius ${{\varepsilon}}$. By the consistency result , we know that the sequence of events $(A_n)_{n\in {\mathbb{N}}}$, for every $n \in {\mathbb{N}}$, given by $$A_n:= (|\hat \theta_{\lambda,\eta_0,\mu}^n -\theta_{\lambda}|\leq {{\varepsilon}}) = (\hat \theta_{\lambda,\eta_0,\mu}^n \in\overline{B(\theta_{\lambda},{{\varepsilon}})}),$$ satisfies $\bP(A_n)\to 1$ as $n\to \infty$. By assumption \[ass:MonotonicityAndContinuityOfTest\] we have that $$\begin{aligned}
\tilde \lambda\mapsto T_n(\theta^n_{\tilde \lambda, \eta_0,\mu} ), \qquad \text{and} \qquad \theta \mapsto T_n(\theta ),\end{aligned}$$ are weakly decreasing and continuous, respectively. Together with the continuity of $\tilde \lambda \mapsto \hat \theta_{\tilde \lambda,\eta_0,\mu }^n$, this implies that also the mapping $\tilde \lambda \mapsto T_n(\hat \theta_{\tilde \lambda,\eta_0,\mu }^n)$ is continuous. It now follows from Assumption \[ass:LambdaStarAlmostSurelyFinite\] (stating that $\lambda^\star_n$ is almost surely finite) that for all $n \in {\mathbb{N}}$, $\bP(T_{n}( \hat \theta^{n}_{\lambda^\star_{n},\eta_0,\mu} ) \leq q(\alpha))=1$. Furthermore, since $\tilde \lambda\mapsto T_n(\theta^n_{\tilde \lambda, \eta_0,\mu} )$ is weakly decreasing, it follows that $$\begin{aligned}
\bP (\lambda^\star_{n} \leq \lambda ) &= \bP( \{\lambda^\star_{n} \leq \lambda \}\cap \{T_{n}( \hat \theta^{n}_{\lambda^\star_{n}, \eta_0,\mu} ) \leq q(\alpha)\} )  \\
&\leq  \bP( \{\lambda^\star_{n} \leq \lambda \} \cap \{T_{n}( \hat \theta^{n}_{\lambda, \eta_0,\mu} ) \leq q(\alpha) \}) \\
& =\bP( \{\lambda^\star_{n} \leq \lambda \}\cap \{T_{n}( \hat \theta^{n}_{\lambda, \eta_0,\mu} ) \leq q(\alpha)\} \cap A_{n}) \\
& \qquad \qquad + \bP( \{\lambda^\star_{n} \leq \lambda \} \cap \{T_{n}( \hat \theta^{n}_{\lambda, \eta_0,\mu} ) \leq q(\alpha)\} \cap A_{n}^c) \\
&\leq  \bP( \{\lambda^\star_{n} \leq \lambda \} \cap  \{ T_{n}( \hat \theta^{n}_{\lambda, \eta_0,\mu} ) \leq q(\alpha)\} \cap \{ |\hat \theta_{\lambda,\eta_0,\mu}^n -\theta_{\lambda}|\leq {{\varepsilon}}\} ) + \bP(A_n^c).\end{aligned}$$ It now suffices to show that the first term converges to zero, since $\bP(A_{n}^c)\to 0$ as $n\to \infty$. We have $$\begin{aligned}
&\bP( \{\lambda^\star_{n} \leq \lambda \} \cap \{ T_{n}( \hat \theta^{n}_{\lambda, \eta_0,\mu} ) \leq q(\alpha)\} \cap \{ |\hat \theta_{\lambda,\eta_0,\mu}^n -\theta_{\lambda}|\leq{{\varepsilon}}\} ) \\
& \quad \leq \bP \Big( \{ \lambda^\star_{n} \leq \lambda \} \cap \Big\{ \inf_{\theta\in \overline{B(\theta_{\lambda},{{\varepsilon}})}}T_{n}( \theta ) \leq q(\alpha) \Big\} \cap \{ |\hat \theta_{\lambda,\eta_0,\mu}^n -\theta_{\lambda}|\leq{{\varepsilon}}\} \Big)  \\
& \quad \leq \bP \Big(  \inf_{\theta\in \overline{B(\theta_{\lambda},{{\varepsilon}})}}T_{n}( \theta ) \leq q(\alpha) \Big)\\ 
& \quad \stackrel{P}{\to}0, \end{aligned}$$ as $n\to \infty$, since $\overline{B(\theta_{\lambda},{{\varepsilon}})}$ is a compact set not containing $\theta^0$. Here, we used that the test statistic $(T_n)$ is assumed to have compact uniform power \[ass:ConsistentTestStatistic\]. Hence, $\lim_{n\to \infty} \bP (\lambda^\star_{n} \leq \lambda ) = 0 $ for any $\lambda \geq 0$, proving that $\lambda^\star_n$ diverges to infinity in probability, which ensures consistency.

**Case (ii):** the unconfounded case ${\mathbb{E}}[B(X)\xi_Y]=0$. Recall that $$\begin{aligned}
 \notag
\hat{\theta}^n_{\lambda,\eta_0,\mu} \,:&= \operatorname*{argmin}_{\theta \in {\mathbb{R}}^{k}}  
\norm{{\mathbf{Y}} - {\mathbf{B}} \theta }_2^2 + \lambda \norm{{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2 + \gamma \theta^\top {\mathbf{K}} \theta  \\ \label{eq:ThetaLambdaMinimizesObjectiveFunction}
&=\operatorname*{argmin}_{\theta \in {\mathbb{R}}^{k}}   l_{\text{OLS}}^n(\theta) + \lambda l_{\text{TSLS}}^n(\theta) + \gamma l_{\text{PEN}}(\theta) ,\end{aligned}$$ where we defined $l_{\text{OLS}}^n(\theta):= n^{-1}\norm{{\mathbf{Y}} - {\mathbf{B}} \theta }_2^2$, $l_{\text{TSLS}}^n(\theta) := n^{-1}\norm{{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \theta)}_2^2$, and $l_{\text{PEN}}(\theta) := n^{-1} \theta^\top {\mathbf{K}} \theta$. For any $0\leq \lambda_1 < \lambda_2$ we have $$\begin{aligned}
&l_{\text{OLS}}^n(\hat{\theta}^n_{\lambda_1,\eta_0,\mu}) + \lambda_1 l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_1,\eta_0,\mu}) + \gamma l_{\text{PEN}}(\hat{\theta}^n_{\lambda_1,\eta_0,\mu})   \\
&\quad \leq l_{\text{OLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) + \lambda_1 l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) + \gamma l_{\text{PEN}}(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) \\
&\quad = l_{\text{OLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) + \lambda_2 l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) + \gamma l_{\text{PEN}}(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) + (\lambda_1-\lambda_2) l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}) \\
&\quad \leq l_{\text{OLS}}^n(\hat{\theta}^n_{\lambda_1,\eta_0,\mu}) + \lambda_2 l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_1,\eta_0,\mu}) + \gamma l_{\text{PEN}}(\hat{\theta}^n_{\lambda_1,\eta_0,\mu}) + (\lambda_1-\lambda_2) l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}),\end{aligned}$$ where we used . Rearranging this inequality and dividing by $(\lambda_1 - \lambda_2)$ yields $$\begin{aligned}
 l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_1,\eta_0,\mu}) \geq   l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda_2,\eta_0,\mu}),\end{aligned}$$ proving that $\lambda \mapsto l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda,\eta_0,\mu})$ is weakly decreasing. Thus, since $\lambda^\star_n \geq 0$ almost surely, we have that $$\begin{aligned}
 \label{eq:unconfounded_IVinOLSboundedByConvergenceToZero}
l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}) \leq l_{\text{TSLS}}^n(\hat{\theta}^n_{0,\eta_0,\mu}) =n^{-1} ({\mathbf{Y}} - {\mathbf{B}} \hat{\theta}^n_{0,\eta_0,\mu})^{\top } {\mathbf{P}}_\delta{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \hat{\theta}^n_{0,\eta_0,\mu}).\end{aligned}$$ Furthermore, recall from that $$\begin{aligned}
 \label{eq:unconfoundedOLSconvergence}
\hat{\theta}^n_{0,\eta_0,\mu}  \underset{n\to\infty}{\stackrel{P}{\longrightarrow}} \theta_0 = \theta^0,\end{aligned}$$ where the last equality follows from using that we are in the unconfounded case ${\mathbb{E}}[B(X)\xi_Y]=0$. By expanding and deriving convergence statements for each term, we get $$\begin{aligned}
 \notag
&({\mathbf{Y}} - {\mathbf{B}} \hat{\theta}^n_{0,\eta_0,\mu})^{\top } {\mathbf{P}}_\delta{\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \hat{\theta}^n_{0,\eta_0,\mu}) \\ \notag
&\quad \underset{n\to\infty}{\stackrel{P}{\longrightarrow}}  ({\mathbb{E}}[ YC^\top ] - \theta_0{\mathbb{E}}[B C^\top ]) {\mathbb{E}}[C^\top C]^{-1} (
{\mathbb{E}}[CY]
- {\mathbb{E}}[CB^\top] \theta_0 ) \\
& \quad = 0, \label{eq:unconfoundedTSLSinOLSConvpZero}\end{aligned}$$ where we used Slutsky’s theorem, the weak law of large numbers, and . Thus, by and it holds that $$\begin{aligned}
l_{\text{TSLS}}^n(\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}) = n^{-1}\| {\mathbf{P}}_\delta({\mathbf{Y}} - {\mathbf{B}} \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu})  \|_2^2 \underset{n\to\infty}{\stackrel{P}{\longrightarrow}} 0.\end{aligned}$$ For any $z\in {\mathbb{R}}^n$ we have that $$\begin{aligned}
\| {\mathbf{P}}_\delta z  \|_2^2 &  = 
z^\top  \fC ( \fC^\top \fC + \delta \fM )^{-1}
\fC^\top\fC ( \fC^\top \fC + \delta \fM )^{-1}
\fC^\top z \\
& =
z^\top  \fC ( \fC^\top \fC + \delta \fM )^{-1}
(\fC^\top\fC)^{1/2}(\fC^\top\fC)^{1/2} ( \fC^\top \fC + \delta \fM )^{-1}
\fC^\top z \\ 
&= \| (\fC^\top\fC)^{1/2} ( \fC^\top \fC + \delta \fM )^{-1}
\fC^\top z  \|_2^2,\end{aligned}$$ hence $$\begin{aligned}
 \label{eq:consitUconfoundedDifferenceConvPtoZero}
\norm{H_n-G_n \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} }_2^2 &= \| n^{-1/2}(\fC^\top\fC)^{1/2} ( \fC^\top \fC + \delta \fM )^{-1}
\fC^\top ({\mathbf{Y}} - {\mathbf{B}} \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu})\|_2^2 \\ \notag
&\stackrel{P}{\to } 0,\end{aligned}$$ where for each $n \in {\mathbb{N}}$, $G_n\in {\mathbb{R}}^{k\times k}$ and $H_n \in {\mathbb{R}}^{k\times 1}$ are defined as $$\begin{aligned}
G_n &:= n^{-1/2}(\fC^\top\fC)^{1/2} ( \fC^\top \fC + \delta \fM )^{-1} \fC^\top \fB, \text{ and } \\
H_n &:= n^{-1/2}(\fC^\top\fC)^{1/2} ( \fC^\top \fC + \delta \fM )^{-1} \fC^\top \fY.\end{aligned}$$ Using the weak law of large numbers, the continuous mapping theorem and Slutsky’s theorem, it follows that, as $n \to \infty$, $$\begin{aligned}
G_n \stackrel{P}{\to}  G    &:= E[CC^\top]^{1/2} E[CC^\top]^{-1} E[CB^\top], \text{ and }\\
H_n \stackrel{P}{\to}  H    &:= E[CC^\top]^{1/2} E[CC^\top]^{-1} E[CY] \\
                                        & = E[CC^\top]^{1/2} E[CC^\top]^{-1} E[CB^\top ]\theta^0 \\
                                        &= G\theta^0,\end{aligned}$$ where the second to last equality follows from . Together with , we now have that $$\begin{aligned}
\norm{G_n \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}  - G \theta^0}_2^2 
&\leq \norm{G_n \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu}  - H_n}_2^2 + 
\norm{H_n  - G \theta^0}_2^2 \underset{n\to\infty}{\stackrel{P}{\longrightarrow}} 0.\end{aligned}$$ Furthermore, by the rank assumptions in \[ass:RankCondition\] we have that $G_n\in {\mathbb{R}}^{k\times k}$ is of full rank (with probability tending to one), hence $$\begin{aligned}
\|\hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} -\theta^0\|_2^2 &= \|G_n^{-1}G_n( \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} -\theta^0)\|_2^2 \\
& \leq \|G_n^{-1}\|_{\text{op}}^2 \|G_n( \hat{\theta}^n_{\lambda^\star_n,\eta_0,\mu} -\theta^0) \|_2^2 \\
&\stackrel{P}{\to}\|G^{-1}\|_{\text{op}}^2 \cdot 0 \\
&=0,\end{aligned}$$ as $n\to \infty$, proving the proposition.

[^1]: This follows from choosing $A$ as an independent noise variable and a constant $g$.

[^2]: This can be assumed without loss of generality if ${\mathcal{F}}$ and ${\mathcal{G}}$ are closed under addition and scalar multiplication, and contain the constant function.

[^3]: For appropriate choices of $h_2$, the model includes settings in which (some of) the $A$ directly influence $Y$.

[^4]: This is wlog if ${\mathcal{F}}$ and ${\mathcal{G}}$ are closed under addition and scalar multiplication, and contain the constant function.

[^5]: For an appropriate choice of $h_2$, the model includes settings in which (some of) the $A$ directly influence $Y$.

[^6]: It is in fact sufficient if the marginal distribution of $X$, ${\mathbb{E}}_{\tilde{M}(i)}[Y\,\vert\, X]$ and ${\mathbb{E}}_{\tilde{M}(i)}[Y^2\,\vert\, X]$ remain fixed for all $\tilde{M}\in\mathcal{M}$ with ${\mathbb{P}}_{\tilde{M}} = {\mathbb{P}}_M$.

[^7]: This may not come as a surprise since without the help of an instrument, it is impossible to distinguish whether a covariate is an ancestor or a descendant of $Y$.
