Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
106,472 characters · 28 sections · 95 citation commands
Distributional robustness of K-class estimators and the PULSE
Keywords: Causality, distributional robustness, instrumental variables
\crefname{section}{section}{sections} \Crefname{section}{Section}{Sections} \crefname{appendix}{appendix}{appendices} \Crefname{appendix}{Appendix}{Appendices}
Learning causal parameters from data has been a key challenge in many scientific fields and has been a long-studied problem in econometrics goldberger1972structural, koopmans1953studies,wold1954causality. Many years after the groundbreaking work by Fisher35 and Peirce83, causality plays again an increasingly important role in machine learning and statistics, two research areas that are most often considered part of mathematics or computer science Imbens2015, Pearl2009, Peters2017,Spirtes2000. Even though current research in mathematics and computer science on the one hand and econometrics on the other, does not develop independently, we believe that there is a lot of potential for more fruitful interaction between these two fields. Differences in the language have emerged, which can make communication difficult, but the target of inference, the underlying principles, and the methodology in both fields are closely related. This paper establishes a link between two developments in these fields: K-class estimation which aims at estimation of causal parameters with good statistical properties and invariance principles that are used to build methods that are robust with respect to distributional shifts. This connection allows us to prove distributional robustness guarantees for K-class estimators and motivates a new estimator, PULSE. We summarize our main results in Section (ref).
Given causal background knowledge, causal parameters can be estimated when taking into account confounding effects between treatment and outcome. Several related techniques have been suggested to tackle that problem, including variable adjustment Pearl2009, propensity score matching Rosenbaum83, inverse probability weighting Horvitz1952 or G-computation Robins1986.
If some of the relevant variables have not been observed, one may instead use exogenous variation in the data to infer causal parameters, e.g., in the setting of instrumental variables angrist1995identification,Newey2013,Wang2016,wright1928tariff. Limited information estimators leverage instrumental variables to conduct single equation inference. An example of such methods is the two-stage least squares estimators (TSLS) developed by theil1953repeated. Instead of minimizing the residual sum of squares as done by the ordinary least square (OLS) estimator, the TSLS minimizes the sample-covariance between the instruments and regression residuals. TSLS estimators are consistent, but are known to have suboptimal finite sample properties, e.g., they only have moments op to the degree of over-identification mariano1972existence. Kadane1971 shows that under suitable conditions, the mean squared error of TSLS might even be larger than the one of OLS if the sample size is small (more precisely and using the notation introduced below, if $0\leq n-q\leq 2(3-(q_2-d_1))$, where $q_2-d_1$ is the degree of overidentification). This result is another indication that under certain conditions, it might be beneficial to use the OLS for regularization. Another method of inferring causal parameters in structural equation models is the limited information maximum likelihood (LIML) estimator due to anderson1949estimation. theil1958economic introduced K-class estimators, which contain OLS, TSLS and the LIML estimator as special cases. This class of estimators is parametrized by a deterministic or stochastic parameter $\kappa\in[0,\infty)$ that depends on the observational data. Under mild regularity conditions a member of this class is consistent and asymptotically normally distributed if $(\kappa-1)$ and $\sqrt{n}(\kappa-1)$ converge, respectively, to zero in probability when $n$ tends to infinity; see, e.g.\ MARIANO1975, mariano2001simultaneous. While the LIML does not have moments of any order, it shares the same asymptotic normal distribution with TSLS. Based on simulation studies, anderson_hildenbrand_1983 argues that, in many practically relevant cases, the normal approximation to a finite-sample estimator is inadequate for TSLS but a useful approximation in the case of LIML. Using Monte Carlo simulations, hahn2004estimation recommend that the no-moment estimator LIML should not be used in weak instrument situations, where Fuller estimators have a substantially smaller MSE. The Fuller estimators fuller1977some form a subclass of the K-class estimators based on a modification to the LIML, which fixes the no-moment problem while maintaining consistency and asymptotic normality. Kiviet2020 proposes a modification to the OLS estimator that makes use of explicit knowledge of the partial correlation between the covariates and the unobserved noise in $Y$. AndrewsArmstrong2017 propose an unbiased estimator that is based on knowledge of the sign of the first stage regression and the variance the reduced form errors and that is less dispersed than TSLS, for example. judge2012minimum consider an affine combination of the OLS and TSLS estimators, which, again, yields a modification in the space of estimators. We prove that our proposed estimator, PULSE, can also be written as a data driven K-class estimator. As such, it minimizes a convex combination of the OLS and TSLS loss functions and can, in general, not be written as a convex combination of the estimators.
All of the above methods exploit background knowledge, e.g., in form of exogeneity of some of the variables. If no such background knowledge is available, it may still be possible, under additional assumptions, to infer the causal structure, e.g., represented by a graph, from observational (or observational and interventional) data. This problem is sometimes referred to as causal discovery. Constraint-based methods assume that the underlying distribution is Markov and faithful with respect to the causal graph and perform conditional independence tests to infer (parts of) the graph; see, e.g.\ Spirtes2000. Score-based methods assume a certain statistical model and optimize (penalized) likelihood scores; see, e.g.\ Chickering2002. Some methods exploit a simple form of causal assignments, such as additive noise (e.g., Peters2014jmlr, and Shimizu2006) and others are based on exploiting invariance statements Peters2016jrssb, Meinshausen2016. Many of such methods assume causal sufficiency, i.e., that all causally relevant variables have been observed, but some versions exist that allow for hidden variables; see, e.g.\ Claassen2013 and Spirtes1995.
Recent works in the fields of machine learning and computational statistics HeinzeDeml17,Pfister2019pnas,Scholkopf2012 investigate whether causal ideas can help to make machine learning methods more robust. The reasoning is that causal models are robust against any intervention in the following sense. Consider a target or response variable $Y$ and covariates $X_1, \ldots, X_p$. If we regress $Y$ on the set $X_S$, $S \subseteq \{1, \ldots, p\}$, of direct causes, then this regression function $x \mapsto E[Y|X_S = x]$ does not change when intervening on any of the covariates (which is sometimes referred to as `invariance'). This statement can be proved using the local Markov property Lauritzen1996, for example, but the underlying fundamental principle has been discussed already several decades ago; most prominently using the terms `autonomy' or `modularity' (Haavelmo1944, and Aldrich1989). As a result, causal models of the form $x \mapsto E[Y|X_S = x]$ may perform well in prediction tasks, where, in the test distribution, the covariates have been intervened on. If, however, training and test distributions coincide, a model focusing only on prediction and the estimand $x \mapsto E[Y|X = x]$ may outperform a causal approach.
The two models described above (OLS and the causal model) formally solve a minimax problem on distributional robustness. Consider therefore an acyclic linear structural equation model (SEM) over $(Y,X)$ with observational distribution $F$. Details on SEMs and interventions can be found in (ref). Assume that the assignment for $Y$ equals $Y = \gamma_0^\intercal X + \varepsilon_Y$ for some $\gamma_0\in \mathbb{R}^{d}$. The variables corresponding to non-zero entries in $\gamma_0^\intercal X$ are called the parents of $Y$, and $\varepsilon_Y$ is assumed to be independent of these parents. Then, the mean squared prediction error when considering the observational distribution is not necessarily minimized by $\gamma_0$, that is, in general, we have $\gamma_0 \neq \gamma_{\text{OLS}} := \operatorname*{arg\,min}_{\gamma} E_F\left[ (Y-\gamma^\intercal X)^2\right]$. Intuitively, we may improve the prediction of $Y$ by including other variables than the parents of $Y$, such as its descendants. When considering distributional robustness, we are interested in finding a $\gamma$ that minimizes the worst case expected squared prediction error over a class of distributions, $\mathcal{F}$, that is,
If we observe data from all different distributions in $\mathcal{F}$ (and know which data point comes from which distribution), we can tackle this optimization directly MeinshausenMiniMax. But estimators of (ref) may be available even if we do not observe data from each distribution in $F$. The true causal coefficient $\gamma_0$, for example, minimizes (ref) when $\cF$ is the set of all possible (hard) interventions on $X$ rojas2018invariant. The OLS solution is optimal when $\mathcal{F}$ only contains the training distribution. In this sense, the OLS solution and the true causal coefficient constitutes the end points of a spectrum of estimators that are prediction optimal under a certain class of distributions.
Intuitively, models trading off causality and predictability may perform well in situations, where the test distribution is only moderately different from the training distribution. Anchor regression by AnchorRegression, see Section (ref) for details, is one approach formalizing this intuition in a linear setup. Similarly to an instrumental variable setting, one assumes the existence of exogenous variables that are called $A$ (for anchor) which may or may not act directly on the target $Y$. The proposed estimator minimizes a convex combination of the residual sum of squares and the TSLS loss function and is shown to be prediction optimal in the sense of (ref) for a class $\mathcal{F}$ containing interventions on the covariates up to a certain strength; this strength depends on a regularization parameter: the weight that is used in the convex combination of anchor regression. Other approaches Pfister2019stab, rojas2018invariant,Magliacane2018 search over different subsets $S$ and aim to choose sets that are both invariant and predictive.
This paper contains two main contributions: A distributional robustness property of K-class estimators with fixed $\kappa$-parameter and a novel estimator for causal coefficients called the p-uncorrelated least squares estimator (PULSE). The following two sections summarize our contributions.
In (ref) we show that anchor regression is closely related to K-class estimators. In particular, we prove that for a restricted subclass of models K-class estimators can be written as anchor regression estimators. For this subclass, this directly implies a distributional robustness property of K-class estimators. We then prove a similar robustness property for general K-class estimators with a fixed penalty parameter, and show that these properties hold even if the model is misspecified.
Consider a possibly cyclic linear SEM over the variables $(Y,X,H,A)$ of the form
subject to regularity conditions that ensure that the distribution of $(Y,X,H,A)$ is well-defined. Here, B and M are constant matrices, the random vectors $A$ and $\varepsilon$ are defined on a common probability space $(\Omega, \cF,P)$, $Y$ is the endogenous target for the single equation inference, $X$ are the observed endogenous variables, $H$ are hidden endogenous variables and $A$ are exogenous variables independent from the unobserved noise innovations $\varepsilon$.
SEMs allow for the notion of interventions, i.e., modeling external manipulations of the system. In this work, we are only concerned with interventions on the exogenous variables $A$ of the form $\text{do}(A:=v)$. Because $A$ is exogeneous, these interventions can be defined as follows: they change the distribution of $A$ to that of a random vector $v$. The interventional distribution of the variables $(Y,X,H,A)$ under the intervention $\text{do}(A:=v)$ is given by the simultaneous distribution of $(X_v, Y_v,H_v,v)$ generated by the SEM
Thus, the intervention does not change any of the original structural assignments of the endogenous variables. Instead, the change in the distribution of the exogeneous variable propagates through the system. We henceforth let $E^{\mathrm{do}(A:=v)}$ denote the expectation with respect to the interventional distribution of the system under the intervention $\mathrm{do}(A:=v)$. More details on interventions can be found (ref)
Let $(\fY,\fX,\fH,\fA)$ consist of $n$ row-wise independent and identically distributed copies of the random vector $(Y,X,H,A)$ and consider the single equation of interest
The K-class estimator with parameter $\kappa$ using non-sample information that only $\fZ_* \subset [\fX \, \, \fA]$ have non-zero coefficients in the target equation of interest is given by
where $P_\fA^\perp$ is the projection onto the orthogonal complement of the column space of $\fA$. For a fixed $\kappa \in [0,1)$ K-class estimators can be represented by a penalized regression problem $\hat{\alpha}_{\text{K}}^n(\kappa) = \operatorname*{arg\,min}_{\alpha} l_{\mathrm{OLS}}^n(\alpha)+ \kappa/(1-\kappa)l_{\mathrm{IV}}^n(\alpha)$, where $l_{\mathrm{OLS}}^n$ and $l_{\mathrm{IV}}^n$ are the empirical OLS and TSLS loss functions, respectively. This representation and the ideas of AnchorRegression allow us to prove that K-class estimator converges to a coefficient that is minimax optimal when considering all distributions induced by a certain set of interventions of $A$. More specifically, we show that for a fixed $\kappa$ and regardless of identifiability,
where $ C(\kappa) := \{ v:\Omega \to \mathbb{R}^{q}: \mathrm{Cov}(v,\varepsilon)=0, \, E [ vv^\intercal ] \preceq \frac{1}{1-\kappa} E[AA^\intercal] \}$. The argmin on the right-hand side minimizes the worst case prediction error when considering interventions up to a certain strength (measured by the set $C(\kappa)$). This objective becomes relevant when we consider a response variable with several covariates and aim to minimize the mean squared prediction error of future realizations of the system of interest that do not follow the training distribution. The above result says that if the new realizations correspond to (unknown) interventions on the exogenous variables that are of bounded strength, K-class estimators with fixed $\kappa\in (0,1)$ minimize the worst case prediction performance and, in particular, outperform the true causal parameter and the least squares solution (see also (ref) in (ref)). For $\kappa$ approaching one, we recover the guarantee of the causal solution and for $\kappa$ approaching zero, the set of distributions contains the training distribution. The above minimax property therefore adds to the discussion whether non-consistent K-class estimators with penalty parameter not converging to one can be useful; see, e.g.\ dhrymes1974.
(ref) contains the second main contribution in this work. We propose a novel data driven K-class estimator for causal coefficients, which we call the p-uncorrelated least square estimator (PULSE). As above, we consider a single endogenous target in an SEM (or simultaneous equation model) and aim to predict it from observed predictors that are with a priori (non-sample) information known to be either endogenous or exogenous. The PULSE estimator can be written in several equivalent forms. It can, first, be seen as a data-driven K-class estimator
where
for some pre-specified level of the hypothesis test $p_{\min}\in(0,1)$. In words, the PULSE estimator outputs the K-class estimator closest to the OLS while maintaining a non-rejected test of uncorrelatedness. In principle, PULSE can be used with any testing procedure. The choice of test, however, may influence the difficulty of the resulting optimization problem. In this paper, we investigate PULSE in connection with a specific class of hypothesis tests that, for example, contain the test of anderson1949estimation. For these hypothesis tests we develop an efficient and provably correct optimization method, that is based on binary line search and quadratic programming.
We show that our estimator can, second, be written as the solution to a constrained optimization problem. To that end, define the primal problems
For the choice $t^\star_n := \sup\{t\,:\, \text{testing } \mathrm{Corr}(A,Y-Z\hat{\alpha}_{\text{Pr}}^n(t))=0$ $\text{yields a}$ $p\text{-value} \geq p_{\min}\}$, we provide a detailed analysis proving that $\hat{\alpha}^n_{\mathrm{K}} (\lambda^\star_n/(1+\lambda^\star_n) ) = \hat{\alpha}_{\text{Pr}}^n(t^\star_n)$.
For the testing procedure proposed in this paper, we show that, third, PULSE can be written as
where $\cA_n(1-p_{\min})$ is the non-convex acceptance region for our test of uncorrelatedness.
This third formulation allows for a simple interpretation of our estimator: among all coefficients (not restricted to K-class estimators) that do not yield a rejection of uncorrelatedness, we choose the one that yields the best prediction. If the acceptance region is empty it outputs a warning indicating a possible model misspecification or an assumption violation to the user (in that case, one can formally output another estimator such as TSLS or Fuller, yielding PULSE well-defined).
In the just-identified setup, the TSLS estimator solves a normal equation which is equivalent to setting a sample covariance between the instruments and the resulting prediction residuals to zero; it then corresponds to $t=0$. For this (and the over-identified) setting, we prove that PULSE is a consistent estimator for the causal coefficient.
The TSLS does not have a finite variance if there is insufficient degree of overidentification, for example. In particular for weak instruments, this usually comes with poor finite sample performance. In such cases, however, the acceptance region of uncorrelatedness is usually large. This yields a weak constraint in the optimization problem and the PULSE will be closer to the OLS, which in certain settings suffers from less variability hahn2004estimation,HAHN05. In simulations we indeed see that, similarly to other data-driven K-class estimators that are pulled towards the OLS, such as Fuller estimators, the PULSE comes with beneficial finite sample properties compared to TSLS and LIML.
Unlike other estimators such as LIML or the classical TSLS, the PULSE is well-defined in under-identified settings, too. Here, its objective is still to find the best predictive solution among all parameters that do not reject uncorrelatedness. Uncorrelatedness to the exogeneous variable is sometimes referred to as invariance. The idea of choosing the best predictive among all invariant models has been investigated in several works Pfister2019stab, rojas2018invariant, Magliacane2018 with the motivation to find models that generalize well (in particular, with respect to interventions on the exogenous variables). Existing methods, however, focus on selecting subsets of variables and then consider least squares regression of the response variable onto the full subset. PULSE can recover such type of solutions if they are indeed optimal. But it also allows to search over coefficients that are different from least squares regression for sets of variables. Consequently, PULSE allows us to find solutions in situations, where the above methods would not find any invariant subsets, which may often be the case if there are hidden variables (see (ref) for an example).
We show in a simulation study that there are several settings in which PULSE outperforms existing estimators both in terms of MSE ordering and several one-dimensional scalarizations of the MSE. More specifically, we show that PULSE can outperform the TSLS and Fuller estimators in weak instrument situations, for example, where Fuller estimators are known to have good MSE properties; see, e.g.\ hahn2004estimation and stock2002survey.
Implementation of PULSE and code for experiments (R) are available on GitHub.\footnote{ \url{https://github.com/MartinEmilJakobsen/PULSE}}
\setcounter{equation}{0} In this section we consider K-class estimators (theil1958economic, and nagar1959bias) and show a connection with anchor regression of AnchorRegression. In (ref) we establish the connection in models where we use a priori information that there are no included exogenous variables in the target equation of interest. In (ref) we then show that general K-class estimators can be written as the solution to a penalized regression problem. In (ref) we utilize this representation and the ideas of AnchorRegression to prove a distributional robustness guarantee of general K-class estimators with fixed $\kappa\in[0,1)$, even under model misspecification and non-identifiability. Proofs of results in this section can be found in (ref).
Denote the random vectors $Y\in \mathbb{R}, X\in \mathbb{R}^{d}, A\in \mathbb{R}^{q},H\in \mathbb{R}^r$ and $\varepsilon\in \mathbb{R}^{d+1+r}$ by the target, endogenous regressor, anchors, hidden and noise variables, respectively. Let further $(Y,X,H)$ be generated by the possibly cyclic structural equation model (SEM)
for some random vectors $\varepsilon\protect\mathpalette{\protect\independenT}{\perp} A$ and constant matrices $B$ and $M$. Let $(\fY,\fX,\fH,\fA)$ consist of $n\geq \min\{d,q\}$ row-wise independent and identically distributed copies of the random vector $(Y,X,H,A)$. Solving for the endogenous variables we get the structural and reduced form equations $[\,\fY\,\, \fX\,\, \fH\,]\, \Gamma = \fA M + \bm{\varepsilon}$ and $[\,\fY\,\, \fX\,\, \fH\,] = \fA \Pi + \bm{\varepsilon} \Gamma^{-1}$, where $\Gamma := I-B$ and $\Pi := M \Gamma^{-1}$. Assume without loss of generality that $\Gamma$ has a unity diagonal, such that the target equation of interest is given by
where $(1, -\gamma_0,- \eta_0)\in \mathbb{R}^{(1+d+r)}$, $ \beta_0\in \mathbb{R}^q$ and $\bm{\varepsilon}_Y$ are the first columns of $\Gamma$, $M$ and $\bm{\varepsilon}\in \mathbb{R}^n$ respectively, $\fZ := [ \fX \, \, \fA]$, $\alpha_0 = (\gamma_0, \beta_0)\in \mathbb{R}^{d+q}$ and $\tilde{\fU}_Y := \fH\eta_0 + \bm{\varepsilon}_Y$.
The possible dependence between the noise $\tilde \fU_Y$ and the endogenous variables, i.e., the influence by hidden variables, generally, renders the standard OLS approach for estimating $\alpha_0$ inconsistent. Instead, one can make use of the components in $A$ that have vanishing coefficient in (ref) for consistent estimation. In the remainder of this work, we disregard any a priori (non-sample) information not concerning the target equation. The question of identifiability of $\alpha_0$ has been studied extensively frisch38,Haavelmo1944,measuringtheequationsystems and more recent overviews can be found in, e.g., Didelez2010, fisher1966identification, and greene2003econometric.
We will use the following assumptions concerning the structure of the SEM:
We will henceforth assume that (ref) always holds. This assumption ensure that the SEM and that the TSLS objectives are well-defined. In the above assumptions, $Z_*$ and $\fZ_*$ are generic placeholders for a subset of endogenous and exogenous variables from $[ X^\intercal \, \, A^\intercal ]^\intercal$ and $[ \fX \, \, \fA ]$, respectively, which should be clear from the context in which they are used. Both Assumption (ref) and Assumption (ref) hold if $X$ and $A$ have density with respect to Lebesgue measure, which in turn is guaranteed by Assumption (ref) if $A$ and $\varepsilon$ have density with respect to Lebesgue measure. Assumption (ref) and (ref) implies that the instrumental variable objective functions introduced below is almost surely well-defined and Assumption (ref) yields that the ordinary least square solution is almost surely well-defined. Assumption (ref) implies that $Y,X$ and $H$ all have finite second moments. For Assumption (ref) and (ref) it is necessary that $q\geq \mathrm{dim}(Z_*)$, i.e., that the setup must be just- or over-identified; see (ref) below.
AnchorRegression proposes a method, called anchor regression, for predicting the endogenous target variable $Y$ from the endogenous variables $X$. The collection of exogenous variables $A$, called anchors, are not included in that prediction model. Anchor regression trades off predictability and invariance by considering a convex combination of the ordinary least square (OLS) loss function and the two-stage least square (IV) loss function using the anchors as instruments. More formally, we define
as the population and finite sample versions of the loss functions. $P_\fA = \fA(\fA^\intercal \fA)^{-1}\fA^\intercal$ is the orthogonal projection onto the column space of $\fA$. To simplify notation, we omit the dependence on $Y$, $X$, $A$, $\fA$, $\fX$ or $\fY$ when they are clear from a given context. For a penalty parameter $\lambda> -1$, the anchor regression coefficients are defined as
The estimator $\hat{\gamma}_{\text{AR}}^n(\lambda)$ consistently estimates the population estimand $\gamma_{\text{AR}}(\lambda)$ and minimizes prediction error while simultaneously penalizing a transformed sample covariance between the anchors and the resulting prediction residuals. Unlike the TSLS estimator, for example, the anchor regression estimator is almost surely well-defined under the rank condition of Assumption (ref), even if the model is under-identified, that is, there are less exogenous than endogenous variables. The solution to the empirical minimization problem of anchor regression is given by
which follows from solving the normal equation of (ref).
The motivation of anchor regression is not to infer a causal parameter. Instead, for a fixed penalty parameter $\lambda$, the estimator is shown to possess a distributional or interventional robustness property: the estimator is optimal when predicting under interventions on the exogenous variables that are below a certain intervention strength. By Theorem 1 of AnchorRegression it holds that $$ \gamma_{\text{AR}}(\lambda) = \operatorname*{arg\,min}_{\gamma\in \mathbb{R}^d} \sup_{v\in C(\lambda)} E^{\text{do}(A:=v)}\left[ \left( Y-\gamma^\intercal X \right) ^2 \right],$$ where $ C(\lambda) := \left\lbrace v:\Omega \to \mathbb{R}^q : \text{Cov}(v,\varepsilon)=0, E(v v^\intercal) \preceq (\lambda+1) E(A A^\intercal ) \right\rbrace.$
We now introduce the limited information estimators known as K-class estimators (theil1958economic, and nagar1959bias) used for single equation inference. Suppose that we are given non-sample information about which components of $\gamma_0$ and $\beta_0$, of (ref), are zero. We can then partition $\fX = [ \fX_* \, \, \fX_{-*} ] \in \mathbb{R}^{n\times(d_1+d_2)}$, $\fA = [\fA_{*} \, \, \fA_{-*}]\in \mathbb{R}^{n\times(q_1+q_2)}$ and $\fZ = [ \fZ_* \, \, \fZ_{-*} ] = [ \fX_{*} \, \, \fA_* \, \, \fX_{-*} \, \, \fA_{-*} ]$ with $ \fZ \in \mathbb{R}^{n\times((d_1+q_1)+(d_2+q_2))}$, where $\fX_{-*}$ and $\fA_{-*}$ corresponds to the variables for which our non-sample information states that the components of $\gamma_0$ and $\beta_0$ are zero, respectively. We call the variables corresponding to $\fA_*$ included exogenous variables. Similarly, we write $\gamma_{0} =(\gamma_{0,*},\gamma_{0,-*}) $, $\beta_0 = (\beta_{0,*},\beta_{0,-*})$ and $\alpha_0 = ( \alpha_{0,*}, \alpha_{0,- *}) =(\gamma_{0,*},\beta_{0,*},\gamma_{0,-*},\beta_{0,-*})$. The structural equation of interest then reduces to $ \fY = \fX_* \gamma_{0,*} + \fX_{-*} \gamma_{0,-*} + \fA_* \beta_{0,*} +\fA_{-*} \beta_{0,-*} + \tilde{\fU}_Y = \fZ_* \alpha_{0,*} + \fU_Y$, where $\fU_Y = \fX_{-*} \gamma_{0,-*} + \fA_{-*} \beta_{0,-*}+ \fH\eta_0 + \bm{\varepsilon}_Y$. In the case that the non-sample information is indeed correct, we have that $\fU_Y = \tilde{\fU}_Y= \fH\eta_0 + \bm{\varepsilon}_Y$. When well-defined, the K-class estimator with parameter $\kappa\in \mathbb{R}$ for a simultaneous estimation of $\alpha_{0,*}$ is given by
where $I-\kappa P_\fA^\perp = I- \kappa(I-P_\fA) = (1-\kappa)I + \kappa P_\fA$.
Comparing (ref) suggests a close connection between anchor regression and K-class estimators for inference of structural equations with no included exogenous variables. In the following subsections, we establish this connection and subsequently extend the distributional robustness property to general K-class estimators.
Assume that, in addition to (ref), we have the non-sample information that $\beta_0=0$, that is, no exogenous variable in $A$ directly affects the target variable $Y$. By direct comparison we see that the K-class estimator for $\kappa<1$ coincides with the anchor regression estimator with penalty parameter $\lambda =\kappa/(1-\kappa)$, i.e., $\hat{\gamma}_{\text{K}}^n(\kappa)= \gamma_{\mathrm{AR}}^n \left( \frac{\kappa }{1-\kappa} \right) $. Equivalently, we have $\gamma_{\mathrm{AR}}^n \left( \lambda \right) = \gamma_{\text{K}}^n \left( \lambda /(1+\lambda) \right) $ for any $\lambda >-1$. As such, the K-class estimator, for a fixed $\kappa$, inherits the following distributional robustness property:
where $ C(\kappa/(1-\kappa)) = \{ v:\Omega \to \mathbb{R}^q : \text{Cov}(v,\varepsilon)=0, E[v v^\intercal] \preceq \frac{1}{1-\kappa} E[A A^\intercal] \}$. This statement holds by Theorem 1 of AnchorRegression.
In an identifiable model with $P \lim_{n\to \infty}\kappa =1$ we have that $\hat{\gamma}^n_{\text{K}}(\kappa)$ consistently estimates the causal parameter; see e.g.\ mariano2001simultaneous. For such a choice of $\kappa$, the robustness above is just a weaker version of what the causal coefficient can guarantee. However, the above result in (ref) establishes a robustness property for fixed $\kappa <1$, even in cases where the model is not identifiable. Furthermore, since we did not use that the non-sample information that $\beta_0=0$ was true, the robustness property is resilient to model misspecification in terms of excluding included exogenous variables from the target equation which generally also breaks identifiability.
We now show that general K-class estimators can be written as solutions to penalized regression problems. The first appearance of such a representation is, to the best of our knowledge, due to mcdonald1977k building upon previous work of basmann1960finite,basmann1960asymptotic. Their representation, however, concerns only the endogenous part $\gamma$. We require a slightly different statement and will show that the entire K-class estimator of $\alpha_{0,*}$, i.e., the simultaneous estimation of $\gamma_{0,*}$ and $\beta_{0,*}$, can be written as a penalized regression problem. Let therefore $l_{\mathrm{IV}}(\alpha;\fY,\fZ_*,\fA)$, $l_{\mathrm{IV}}^n(\alpha;\fY,\fZ_*,\fA)$ and $l_{\mathrm{OLS}}(\alpha;\fY,\fZ_*)$, $l_{\mathrm{OLS}}^n(\alpha;\fY,\fZ_*)$ denote the population and empirical TSLS and OLS loss functions as defined in \Crefrange{eq:lossPop}{eq:lossEmp}. That is, the TSLS loss function for regressing $\fY$ on the included endogenous and exogenous variables $\fZ_*$ using the exogeneity of $\fA$ and $\fA_{-*}$ as instruments and the OLS loss function for regressing $\fY$ on $\fZ_*$. We define the K-class population and finite-sample loss functions as an affine combination of the two loss functions above. That is,
Assuming $\kappa\not =1$, we can rewrite (ref) to
Thus, K-class estimators seek to minimize the ordinary least squares loss for regressing $\fY$ on $\fZ_*$, while simultaneously penalizing the strength of a transform on the sample covariance between the prediction residuals and collection of exogenous variables $\fA$.
In the following section, we consider a population version of the above quantity. If we replace the finite sample (ref) with the corresponding population (ref), we get that the minimization estimator of the empirical loss function of (ref) is asymptotically well-defined. Furthermore, we now prove that whenever the population assumptions are satisfied, then, for any fixed $\kappa \in [0,1]$, $\hat{\alpha}_{\text{K}}^n(\kappa;\fY,\fZ_*,\fA)$ converges in probability towards the population K-class estimand.
We are now able to prove that the general K-class estimator possesses a robustness property similar to the statements above. It is prediction optimal under a set of interventions, now including interventions on all exogenous $A$ up to a certain strength.
Here, $E^{\text{do}(A:=v)}$ denotes the expectation with respect to the distribution entailed under the intervention $\text{do}(A:=v)$ (see Section (ref) and (ref)) and $(\Omega,\cF,P)$ is the common background probability space on which $A$ and $\varepsilon$ are defined.
In words, among all linear prediction methods of $Y$ using $Z_*$ as predictors, the K-class estimator with parameter $\kappa$ has the lowest possible worst case mean squared prediction error when considering all interventions on the exogenous variables $A$ contained in $C(\kappa)$. As $\kappa$ approaches one, the estimator is prediction optimal under a class of arbitrarily strong interventions in the direction of the variance of $A$. (Here, $\kappa$ is arbitrary but fixed; the statement does not cover data-driven choices of $\kappa$, such as LIML or Fuller.) The above result is a consequence of the relation between anchor regression and K-class estimators. The special case $A_* = \emptyset$ is a consequence of Theorem 1 by AnchorRegression. Our proof follows similar arguments but additionally allows for $A_* \not = \emptyset$.
The property in Theorem (ref) has a decision-theoretic interpretation (see Chamberlain2007 for an application of decision theory in IV models based on another loss function). Consider a response $Y$, covariates $Z_*$ and a distribution (specified by $\theta$) over $(Y, Z_*)$, and the squared loss $\ell(Y, Z, \alpha) := (Y-\alpha^\top Z_*)^2$. Then, assuming finite variances, for each distribution the risk $E_{\theta}[(Y-\alpha^\top Z_*)^2]$ is minimized by the (population) OLS solution $\alpha=\alpha_\theta := \mathrm{cov}_\theta(Z_*)^{-1}\mathrm{cov}_\theta(Z_*,Y)$. In the setting of Theorem (ref), we are given a distribution over $(Y,Z_*)$, specified by $\theta$, but we are interested in minimizing the risk $E_{\theta, v}[(Y-\alpha_{\theta}^\top Z_*)^2]$ for another distribution that is induced by an intervention and specified by $(\theta, v)$. The above result states that the K-class estimator minimizes a worst-case risk when considering all $v \in C(\kappa)$.
Theorem (ref) makes use of the language of SEMs in that it yields the notion of interventions.\footnote{In particular, we have not considered the SEM as a model for counterfactual statements.} As such, the result can be rephrased using other causal frameworks. The crucial assumptions are the exogeneity of $A$ and the linearity of the system. Furthermore, the result is robust with respect to several types of model misspecifications that breaks identifiability of $\alpha_0$, such as excluding included endogenous or exogenous predictors or the existence of latent variables; see (ref) in (ref).
\setcounter{equation}{0} We now introduce the p-uncorrelated least square estimator (PULSE). As discussed in (ref), PULSE allows for different representations. In this section we start with the third representation and show the equivalence of the other representations afterwards.
Consider predicting the target $Y$ from endogenous and possibly exogenous regressors $Z$. Let therefore $\cH_0(\alpha)$ denote the hypothesis that the prediction residuals using $\alpha$ as a regression coefficient is simultaneously uncorrelated with every exogenous variable, that is, $\cH_0(\alpha) : \text{Corr}(A,Y-\alpha^\intercal Z) =0$. This hypothesis is in some models under certain conditions equivalent to the hypothesis that $\alpha$ is the true causal coefficient. One of these conditions is the rank condition (ref) introduced below, also known as the rank condition for identification; wooldridge2010econometric.
The two-stage least square (TSLS) estimator exploits the equivalence between the causal coefficient and the zero correlation between the instruments and the regression residuals. Here, one minimizes a sample covariance between the instruments and the regression residuals: we can write $l^n_{\mathrm{IV}}(\alpha;\fY,\fZ, \fA) = \| \widehat{\text{Cov}}_n(A,Y-\alpha^\intercal Z) \|^2_{(n^{-1}\fA^\intercal \fA)^{-1}}$ when $A$ is mean zero.\footnote{$\|\cdot\|_{(n^{-1}\fA^\intercal \fA)^{-1}}$ is the norm induced by the inner product $\langle x,y\rangle = x^\intercal (n^{-1}\fA^\intercal \fA)^{-1} y$.} In the just-identified setup the TSLS estimator yields a sample covariance that is exactly zero and is known to be unstable, in that it has no moments of any order. Intuitively, the constraint of vanishing sample covariance may be too strong.
Let $T(\alpha;\fY,\fZ,\fA)$ be a finite sample test statistic for testing the hypothesis $\cH_0(\alpha)$ and let $\text{p-value}(T(\alpha;\fY,\fZ,\fA))$ denote the p-value associated with the test of $\cH_0(\alpha)$. We then define the p-uncorrelated least square estimator (PULSE) as
where $p_{\min}$ is a pre-specified level of the hypothesis test. In words, we aim to minimize the mean squared prediction error among all coefficients which yield a p-value for testing $\cH_0(\alpha)$ that does not fall below some pre-specified level-threshold $p_{\min} \in (0,1)$, such as $p_{\min}= 0.05$. That is, the minimization is constrained to the acceptance region of the test, i.e., a confidence region for the causal coefficient in the identified setup. Among these coefficient, we choose the solution that is `closest' to the OLS solution.\footnote{ Here, closeness is measured in the OLS distance: We define the OLS norm via $\|\alpha\|_{\text{OLS}}^2 := l_{\mathrm{OLS}}^n(\alpha + \hat{\alpha}^n_{\mathrm{OLS}}) - l_{\mathrm{OLS}}^n(\hat{\alpha}^n_{\mathrm{OLS}}) = \alpha^\top \mathbf{Z}^T \mathbf{Z} \alpha$, where $\hat{\alpha}^n_{\mathrm{OLS}}$ is the OLS estimator. This defines a norm (rather than a semi-norm) if $\mathbf{Z}^T \mathbf{Z}$ is non-degenerate. Minimizing $l_{\mathrm{OLS}}^n(\alpha)=\|\fY-\fZ\alpha\|_2^2 = (\alpha-\hat{\alpha}^n_{\mathrm{OLS}})^\intercal \fZ^\intercal \fZ (\alpha- \hat{\alpha}^n_{\mathrm{OLS}})+\|\fY-\fZ \hat{\alpha}^n_{\mathrm{OLS}}\|_2^2$ is equivalent to minimizing $\|\alpha - \hat{\alpha}^n_{\mathrm{OLS}}\|_{\text{OLS}}^2$. }
Thus, PULSE allows for an intuitive interpretation. We will see in the experimental section that it has good finite sample performance, in particular for weak instruments. Unlike other estimators, such as LIML, the above estimator is well-defined in the under-identified setup, too.\footnote{ The PULSE estimator is defined for finite samples, but the following deliberation may help to build intuition: In an under-identified IV setting, minimizing $l_{\mathrm{OLS}}(\gamma)$ under the constraint that $l_{\mathrm{IV}}(\gamma) = 0$, can be seen as choosing, under all causal models compatible with the distribution, the model with the least amount confounding -- when using $E(Y-\gamma^\top X)^2 - E(Y-\gamma_{\mathrm{OLS}}^\top X)^2$ as a measure for confounding. } In such cases, PULSE extends on existing literature that aims to trade-off predictability and invariance but that so far has been restricted to search over subsets of variables (see Sections (ref) and (ref)). To maintain consistency of the estimator the chosen test must have asymptotic power of one.
In this paper, we propose a class of significance tests, that contains, e.g., the Anderson-Rubin test anderson1949estimation. While the objective function in (ref) is quadratic in $\alpha$, the resulting constraint is, in general, non-convex. In Section (ref), we develop a computationally efficient procedure that provably solves the optimization problem at low computational cost. Other choices of tests are possible, too, but may result in even harder optimization problems.
In (ref), we briefly introduce the setup and assumptions. In (ref), we specify a class of asymptotically consistent tests for $\cH_0(\alpha)$. In (ref) we formally define the PULSE estimator. In (ref), we show that the PULSE estimator is well-defined by proving that it is equivalent to a solvable convex quadratically constrained quadratic program which we denote by the primal PULSE. In (ref), we utilize duality theory and derive an alternative representation which we denote by the dual PULSE. This representation yields a computationally feasible algorithm and shows that the PULSE estimator is a K-class estimator with a data-driven $\kappa$. Proofs of results in this section can be found in (ref) unless stated otherwise.
In the following sections we again let $(\fY,\fX,\fH,\fA)$ consist of $n\geq \min\{d,q\}$ row-wise independent and identically distributed copies of $(Y,X,H,A)$ generated in accordance with the SEM in (ref). The structural equation of interest is $Y = \gamma^\intercal_0 X + \eta^\intercal_0 H + \beta^\intercal_0 A + \varepsilon_{Y}$. Assume that we have some non-sample information about which $d_2=d-d_1$ and $q_2=q-q_1$ coefficients of $\gamma_0$ and $\beta_0$, respectively, are zero. As in (ref), we let the subscript $*$ denote the variables and coefficients that are non-zero according to the non-sample information but to simplify notation, we drop the $*$ subscript from $Z$, $\fZ$ and $\alpha_0$; that is, we write $Z =[X_*^\intercal \; A_*^\intercal]^\intercal \in \mathbb{R}^{d_1+q_1}$, $\fZ=[ \fX_* \, \, \fA_* ]\in \mathbb{R}^{n\times(d_1+q_1)}$ and $\alpha_0 :=(\gamma_{0,*}^\intercal,\beta_{0,*}^\intercal)^\intercal:\in \mathbb{R}^{d_1+q_1}$. That is, $Y = \alpha^\intercal_{0} Z + U_Y$, where $U_Y = \alpha_{0,-*}^\intercal Z_{-*}+\eta^\intercal_0 H + \varepsilon_Y$. If the non-sample information is true, then $U_Y = \eta^\intercal_0 H + \varepsilon_Y$.
We define a setup as being under- just- and over-identified by the degree of over-identification $q_2-d_1$ being negative, equal to zero and positive, respectively. That is, the number of excluded exogenous variables $A_{-*}$ being less, equal or larger than the number of included endogenous variables $X_{*}$ in the target equation.
We assume that the global assumptions of (ref) from (ref) still hold. Furthermore, we will make use of the following situational assumptions
Assumption (ref) holds if our non-sample information is true, and the instrument set $A$ is independent of all unobserved endogenous variables $H_i$ which directly affect the target $Y$. This holds, for example, if the latent variables are source nodes, that is, they have no parents in the causal graph of the corresponding SEM. Assumption (ref) can be achieved by centering the data. Strictly speaking, this introduces a weak dependence structure in the observations, which is commonly ignored. Alternatively, one can perform sample splitting. For more details on this assumption and the possibility of relaxing it, see (ref). Assumption (ref) ensures that K-class estimators for $\kappa < 1$ are well-defined, regardless of the over-identification degree. In the under-identified setup, Assumption (ref) yields that there exists a subspace of solutions minimizing $l_{\text{IV}}^n(\alpha)$. In the just- and over-identified setup this assumption ensures that $l_{\text{IV}}^n(\alpha)$ has a unique minimizer given by the two-stage least squares estimator $\hat{\alpha}_{\text{TSLS}}^n := (\fZ^\intercal P_\fA \fZ)^{-1} \fZ^\intercal P_\fA \fY$. (ref) is used to ensure that the ordinary least square objective function $l_{\text{OLS}}^n(\alpha;\fY,\fZ)$ is strictly positive, such that division by this function is always well-defined. (ref) ensure that various limiting arguments are valid. In the just- and over-identified setup (ref) is known as the rank condition for identification.
We now introduce a class of tests for the null hypothesis $\cH_0(\alpha) : \text{Corr}(A,Y-Z\alpha) =0$ that have point-wise asymptotic level and pointwise asymptotic power. These tests will allow us to define the corresponding PULSE estimator. When (ref) holds we can define $T_n^c:\mathbb{R}^{d_1+q_1} \to \mathbb{R}$ by
where $c(n)$ is a function that will typically scale linearly in $n$. Let us denote the $1-p$ quantile of the central Chi-Squared distribution with $q$ degrees of freedom by $Q_{\chi^2_{q}}(1-p)$. By standard limiting theory we can test $\cH_0(\alpha)$ in the following manner.
Depending on the choice of $c(n)$, this class contains several tests, some of which are well known. With $c(n) = n-q+Q_{\chi^2_q}(1-p_{\min})$, for example, one recovers a test that is equivalent to the asymptotic version of the Anderson-Rubin test (anderson1950asymptotic). We make this connection precise in (ref) in (ref). The Anderson-Rubin test is robust to weak instruments in the sense that the limiting distribution of the test-statistic under the null-hypothesis is not affected by weak instrument asymptotics; see, e.g.\ staiger1997instrumental and stock2002survey.\footnote{Weak instrument asymptotics is a model scheme where the instrument strength tends to zero at a rate of $n^{-1/2}$, i.e., the reduced form structural equation for the endogenous variables is given by $\fX = \fA n^{-1/2} \Pi_X +\bm{\varepsilon} \Gamma^{-1}_X$.} For weak instruments, the confidence region may be unbounded with large probability; see dufour1997some. moreira2009tests show that the test suffers from loss of power in the over-identified setting.
To simplify notation, we will from now on work with the choice $c(n) = n$ and define the acceptance region with level $p_{\min}\in(0,1)$ as $\cA_n(1-p_{\min}) := \{\alpha \in \mathbb{R}^{d_1+q_1}: T_n(\alpha) \leq Q_{\chi^2_{q}}(1-p_{\min})\}$, where $T_n(\alpha)$ corresponds to the choice $c(n) = n$.
For any level $p_{\min}\in(0,1)$, we formally define the PULSE estimator of (ref) by letting the feasible set be given by the acceptance region $\cA_n(1-p_{\min})$ of $\cH_0(\alpha)$ using the test of (ref). That is, we consider
In general, this is a non-convex optimization problem (boyd2004convex) as the constraint function is non-convex, see the blue contours in (ref)(left). From (ref)(right) we see that in the given example the problem nevertheless has a unique and well-defined solution: the smallest level-set of $l_{\text{OLS}}^n$ with a non-empty intersection of the acceptance region $\{\alpha : T_n(\alpha) \leq Q_{\chi^2_{q}}(1-p_{\min})\}$ intersects with the latter region in a unique point. In (ref), we prove that this is not a coincidence: (ref) has a unique solution that coincides with the solution of a strictly convex, quadratically constrained quadratic program (QCQP) with a data-dependent constraint bound. In (ref), we further derive an equivalent Lagrangian dual problem. This has two important implications. (1) It allows us to construct a computationally efficient procedure to compute a solution of the non-convex problem above, and (2), it shows that the PULSE estimator can be written as K-class estimators.
Estimators with similar constraints albeit different optimization objective have been studied by gautier2011high. In (ref) in (ref) we briefly discuss the connection to pre-test estimators. Furthermore, any method for inverting the test, see, e.g., davidson2014confidence, yields a valid confidence set including the proposed point estimator (given that the method outputs the point estimator when the acceptance region is empty).
We now derive a QCQP representation of the PULSE problem, which we call the primal PULSE. For all $t\geq 0$ define the empirical primal minimization problem (Primal.$t.n$) by
We drop the dependence of $\fY$, $\fZ$ and $\fA$ and refer to the objective and constraint functions as $l_{\text{OLS}}^n(\alpha)$ and $l_{\text{IV}}^n(\alpha)$. The following lemma shows that under suitable assumptions these problems are solvable, strictly convex QCQP problems satisfying Slater's condition.
We restrict the constraint bounds to $D_{\text{Pr}}:=(\inf_\alpha l_{\mathrm{IV}}^n(\alpha), l_{\text{IV}}^n(\hat{\alpha}_{\text{OLS}}^n)]$. Considering $t$ that are larger than $\inf_\alpha l_{\mathrm{IV}}^n(\alpha)$ ensures that the problem (Primal$.t.n$) is uniquely solvable and furthermore that Slater's condition is satisfied (see Lemma (ref) above). Slater's condition will play a role in (ref) when establishing a sufficiently strong connection with its corresponding dual problem for which we can derive a (semi-)closed form solution. Constraint bounds greater than or equal to $l_{\text{IV}}^n(\hat{\alpha}_{\text{OLS}}^n)$ yield identical solutions. Whenever well-defined, let $\hat{\alpha}_{\text{Pr}}^n:D_{\text{Pr}}\to \mathbb{R}^{d_1+q_1}$ denote the constrained minimization estimator given by the solution to the (Primal$.t.n$) problem
We now prove that for a specific choice of $t$, the PULSE and the primal PULSE yield the same solutions. Define $t_n^\star(p_{\min})$ as the data-dependent constraint bound given by
If $t^\star_n(p_{\min})>-\infty$ or equivalently $t^\star_n(p_{\min})\in D_{\text{Pr}}$ we define the primal PULSE problem and its solution by (Primal$.t^\star_n(p_{\min}).n$) and $\hat{\alpha}_{\mathrm{Pr}}^n(t_n^\star(p_{\min}))$. The following theorem yields conditions for when the solutions to the primal PULSE and PULSE problems coincide.
We show that $t_n^\star(p_{\min}) >-\infty$ is a sufficient condition for $T_n(\hat{\alpha}_{\mathrm{Pr}}^n(t_n^\star(p_{\min})))\leq Q_{\chi^2_{q}}(1-p_{\min})$ in the proof of (ref). The sufficiency of $t_n^\star(p_{\min}) >-\infty$ is postponed to the latter proof as it easily follows from the dual representation. Hence, we have shown that finding the PULSE estimator, i.e., finding a solution to the non-convex PULSE problem, is equivalent to solving the convex QCQP primal PULSE for a data dependent choice of $t_n^\star(p_{\min})$.\footnote{Given that value, we can use a numerical QCQP solver to calculate the PULSE estimate. } However, $t_n^\star(p_{\min})$ is still unknown. (ref) shows an example of the equivalence in (ref). (ref)(right) shows that the level set of $l_{\text{IV}}(\alpha) = t^{\star}(p_{\min})$ intersects the optimal level curve of $l_{\text{OLS}}^n(\alpha)$ in the same point given by minimizing over the constraint $T_n(\alpha) \leq Q_{\chi^2_q}(1-p_{\min})$.
The set of solutions to the primal problem $\{\hat{\alpha}_{\text{Pr}}^n(t) : t \in D_{\text{Pr} }\}$ can in the just- and over-identified setup be visualized as an (in general) non-linear path in $\mathbb{R}^{d_1+q_1}$ between the TSLS estimator $(t= l_{\mathrm{IV}}^n(\hat{\alpha}_{\text{TSLS}}^n))$ and the OLS estimator $(t= l_{\text{IV}}^n(\hat{\alpha}_{\text{OLS}}^n))$ AnchorRegression. Theorem (ref) yields that the PULSE estimator ($t=t_n^\star(p_{\min})$) then seeks the estimator 'closest' to the OLS estimator along this path that does not yield a rejected test of simultaneous vanishing correlation between the resulting prediction residuals and the exogenous variables $A$, see (ref). The path of possible solutions is not necessarily a straight line (see black line); thus, in general, the PULSE estimator is different from the affine combination of OLS and TSLS estimators studied by e.g.\ judge2012minimum.
In the under-identified setup, the TSLS end point corresponding to $t= \min_\alpha l_{\mathrm{IV}}^n(\alpha)$ is instead given by the point in the IV solution space $\{\alpha \in \mathbb{R}^{d_1+q_1}: l_{\text{IV}}^n(\alpha)=0\}$ with the smallest mean squared prediction residuals.
In this section, we derive a dual representation of the primal PULSE problem which we will denote the dual PULSE problem. This specific dual representation allows for the construction of a binary search algorithm for the PULSE estimator and yields that PULSE is a member of the K-class estimators with stochastic $\kappa$-parameter.
For any penalty parameter $\lambda \geq 0$ we define the dual problem (Dual$.\lambda.n$) by
Whenever Assumption (ref) holds, i.e., $\fZ^\intercal \fZ$ is of full rank, then for any $\lambda\geq 0$ the solution to (Dual$.\lambda.n$) coincides with the K-class estimator with $\kappa = \lambda/(1+\lambda)\in[0,1)$, see (ref). That is,
solves (Dual$.\lambda.n$). Henceforth, let $\hat{\alpha}_{\mathrm{K}}^n (\lambda)$ denote the solution to (Dual$.\lambda.n$), i.e., in a slight abuse of notation we will denote the solution to (Dual$.\lambda.n$) by $\hat{\alpha}_{\text{K}}^n(\lambda)$, such that $\hat{\alpha}_{\text{K}}^n(\lambda)= \hat{\alpha}_{\text{K}}^n(\kappa)$ for $\kappa = \lambda/(1+\lambda)$. We refer to these two representations as the K-class estimator with penalty parameter $\lambda$ and parameter $\kappa$, respectively. The usage of $\kappa$ or $\lambda$ as argument should clarify which notation we refer to.
Under Assumption (ref) we have that the minimum of $l_{\text{IV}}^n(\alpha)$ is attainable (see the proof of (ref)). Hence, let the solution space for the minimization problem $\min_\alpha l_{\text{IV}}^n(\alpha)$ be given by
In the under-identified setup $(q_2<d_1)$, $\cM_{\text{IV}}$ is a $(d_1-q_2)$-dimensional subspace of $\mathbb{R}^{d_1+q_1}$ and in the just- and over-identified setup it holds that $\cM_{\text{IV}}= \{\hat{\alpha}_{\text{TSLS}}^n\}$.
We now prove that, in the generic case, K-class estimators for $\lambda \in [0,\infty)$ are different from the TSLS estimator. This result may not come as a surprise, but we include it as we need the result later and have not found it elsewhere.
We conjecture that the corresponding statement holds in the under-identified setup and without the lower triangular assumption on B, too. That is, $\cM_{\text{IV}}\cap \{\hat{\alpha}_{\mathrm{K}}^n(\lambda): \lambda \geq 0 \}=\emptyset$ holds almost surely. We therefore introduce this as an assumption.
Furthermore, when imposing that (ref) holds we also have that the K-class estimators differ from each other.
The above corollary is proven as (ref) in (ref). We now show that the class of K-class estimators with penalty parameter $\lambda \geq 0$ , i.e., $\kappa\in[0,1)$, coincides with the class of constrained minimization-estimators that minimize the primal problems with constraint bounds $t> \min_{\alpha}l_{\mathrm{IV}}^n(\alpha)$.
(ref) tells us that, under appropriate assumptions, $ \{\hat{\alpha}_{\text{K}}^n(\kappa ): \kappa \in[0,1)\}= \{\hat{\alpha}_{\text{K}}^n(\lambda):\lambda \geq 0\} = \{\hat{\alpha}_{\text{Pr}}^n(t): t\in D_{\text{Pr}} \}. $ In words, we have recast the K-class estimators with $\kappa \in [0,1)$ as the class of solutions to the primal problems previously introduced. That the minimizers of $l_{\text{IV}}^n(\alpha)$ are different from all the K-class estimators with penalty $\lambda \geq 0$ (or $\kappa\in[0,1)$) guarantees that when representing a K-class problem in terms of a constrained optimization problem it satisfies Slater's condition.
We are now able to show the main result of this section. The PULSE estimator $\hat{\alpha}^n_{\text{PULSE}}(p_{\min})$ solves a K-class problem (Dual$.\lambda.n$) and can therefore be seen as a K-class estimator with a data-dependent parameter. To see this, let us define the dual PULSE penalty parameter, i.e., the dual analogue of the primal PULSE constraint $t_n^\star(p_{\min})$ as
If $\lambda^\star_n(p_{\min})<\infty$, we define the dual PULSE problem by (Dual.$\lambda^\star_n(p_{\min}).n)$ with solution $ \hat{\alpha}_{\text{K}}^n (\lambda_n^\star(p)) = \operatorname*{arg\,min}_{\alpha\in \mathbb{R}^{d_1+q_1}} l_{\text{OLS}}^n(\alpha) + \lambda_n^\star(p_{\min}) l_{\text{IV}}^n(\alpha). $
Thus, the PULSE estimator seeks to minimize the K-class penalty $\lambda$, i.e., to pull the estimator along the K-class path $\{\hat{\alpha}_{\text{K}}^n(\lambda):\lambda \geq 0\}$ as close to the ordinary least square estimator as possible. Furthermore, the statement implies that the PULSE estimator is a K-class estimator with data-driven penalty $\lambda_n^\star(p_{\min})$ or, equivalently, parameter $\kappa = \lambda_n^\star(p_{\min})/(1+\lambda_n^\star(p_{\min}))$. Given a finite dual PULSE penalty parameter $\lambda_n^\star(p_{\min})$ we can, by utilizing the closed form solution of the K-class problem, represent the PULSE estimator in the following form:
However, to the best of our knowledge, $\lambda_n^\star(p_{\min})$ has no known closed form, so the above expression cannot be computed in closed-form either. In (ref), we prove that the PULSE penalty parameter $\lambda_n^\star(p_{\min})$ can be approximated with arbitrary precision by a simple binary search procedure.
The following lemma provides a necessary and sufficient (in practice checkable) condition for when the PULSE penalty parameter $\lambda_n^\star(p_{\min})$ is finite.
Thus, under suitable regularity assumptions (ref) yields that our dual representation of the PULSE estimator always holds in the under- and just-identified setup. It furthermore yields a sufficient and necessary condition for the dual representation to be valid in the over-identified setup, namely that the TSLS is in the interior of the acceptance region. Furthermore, this condition is possibly violated in the over-identified setup with non-negligible probability.
The key insight allowing for a binary search procedure for $\lambda_n^\star(p_{\min})$ is that the mapping $\lambda \mapsto T_n(\hat{\alpha}_{\text{K}}^n(\lambda))$ is monotonically decreasing.
The above lemma is proven as (ref) in (ref). If the OLS solution is not strictly feasible in the PULSE problem, then $\lambda_n^\star(p_{\min})$ indeed is the smallest penalty parameter for which the test-statistic reaches a p-value of exactly $p_{\min}$; see (ref) in (ref).
We propose the binary search algorithm presented in Algorithm (ref) in (ref), that can approximate a finite $\lambda_n^\star(p_{\min})$ with arbitrary precision. We terminate the binary search (see line 2) if $\lambda^\star_n(p_{\min})$ is not finite, in which case we have no computable representation of the PULSE estimator. It is possible to improve this algorithm in the under- and just-identified setup, by initializing $\ell_{\max}$ as the quantity given by (ref) in the proof of (ref). This initialization removes the need for the first while loop in (lines 4--6). We now prove that Algorithm (ref) achieves the required precision and is asymptotically correct.
The dual representation of the PULSE estimator is not guaranteed to be well-defined in the over-identified setup. In particular, it is not well-defined if the TSLS is outside the interior of the acceptance region (which corresponds to a p-value of less than or equal to $p_{\min}$). In this case, we propose to output a warning. This can be helpful information for the user since it may indicate a model misspecification. For example, if the true relationship is in fact nonlinear, and one considers an over-identified case (e.g., by constructing different transformations of the instrument), even the TSLS may be rejected when erroneously considering a linear model; see Keane2010 and Mogstad2010. For any $p_{\min}\in(0,1)$ we can still define an always well-defined modified PULSE estimator $\hat{\alpha}_{\text{PULSE}+}^n(p_{\min})$ as $\hat{\alpha}_{\mathrm{PULSE}}^n(p_{\min})$ if the dual representation is feasible and some other consistent estimator $\hat{\alpha}^n_{\text{ALT}}$ (such as the TSLS, LIML or Fuller estimator) otherwise. That is, we define
Similarly to the case of an empty rejection region, we also output a warning for the case when the OLS estimator is accepted. This may, but does not have to, indicate weak instruments. Thus, we have the algorithm presented as (ref) in (ref) for computing the PULSE$+$ estimator.
We now prove that the PULSE$+$ estimator consistently estimates the causal parameter in the just- and over-identified setting. Assume that we choose a consistent estimator $\hat{\alpha}_{\text{ALT}}^n$ (under standard regularity assumptions, this is satisfied for the TSLS).\footnote{This holds as $ \hat{\alpha}_{\text{TSLS}}^n = \alpha_0 + (n^{-1}\fZ^\intercal \fA (n^{-1}\fA^\intercal \fA)^{-1} n^{-1}\fA^\intercal \fZ)^{-1} n^{-1}\fZ^\intercal \fA (n^{-1}\fA^\intercal \fA)^{-1} n^{-1}\fA^\intercal \fU_Y $.} We can then show that, under mild conditions, the PULSE$+$ estimator, too, is a consistent estimator of $\alpha_0$.
We believe that a similar statement also holds in the under-identified setting, see (ref).
In (ref) we conduct an extensive simulation study investigating the finite sample behaviour of the PULSE estimator. The concept of weak instruments is central to our analysis. An introduction to weak instruments can be found in (ref). Here we give a brief overview of the study and the observations.
The theoretical results on distributional robustness proved in (ref) translate to finite data. The experiments of (ref) shows that even for small sample sizes, K-class estimators outperform both OLS and TSLS for a certain range of interventions, matching the theoretical predictions with increasing sample size. In (ref), we furthermore consider an under-identified setting.
When focusing on the estimation of a causal effect in an identified setting, our simulations show that there are several settings where PULSE outperforms the Fuller and TSLS estimators in terms of mean squared error (MSE). In univariate simulation experiments, such settings are characterized by weakness of instruments and weak confounding (endogeneity). The characterization becomes more involved in multivariate settings, but is similar in that PULSE outperforms all other methods for small confounding strengths, an effect amplified by the weakness of instruments. Below we detail the univariate simulation setup and refer the reader to (ref) for further details and the multivariate simulation experiments mentioned above.
We first compared performance measures of the estimators in a univariate instrumental variable model. As seen in hahn2002new and hahn2004estimation, we consider structural equation models of the form
where $A\sim \mathcal{N}(0,I)$ and $A\protect\mathpalette{\protect\independenT}{\perp} (U_X,U_Y) $ with $
\sim \mathcal{N}\left(
,
\right).$ Furthermore, we let $\gamma =1$ and $\bar{\xi}^\intercal=(\xi,....,\xi)\in \mathbb{R}^q$, where $\xi>0$ is chosen according to the theoretical $R^2$-coefficient. We consider the following simulation scheme: for each $q\in \{1,2,3,4,5,10,20,30\}$, $\rho\in\{0.1,0.2,...,0.9\}$, $R^2\in\{0.0001,0.001,0.01,0.1,0.3\}$ and $n\in\{50,100,150\}$, we simulate $n$-samples from the above system and calculate the OLS, TSLS, Fuller(1), Fuller(4) and PULSE ($p_{\min}= 0.05$) estimates; see (ref).
(ref) contains illustrations of the relative change in square-root mean squared error (RMSE) estimated from $15 000$ repetitions. On the horizontal axis we have plotted the average first stage F-test as a measure of weakness of instruments; see (ref) for further details. A test for $H_0:\bar \xi=0$, i.e., for the relevancy of instruments, at a significance level of 5%, has different rejection thresholds in the range $[1.55,4.04]$ depending on $n$ and $q$. The vertical dashed line corresponds to the smallest rejection threshold of 1.55 and the dotted line corresponds to the `rule of thump' threshold of 10. Note that the lowest possible negative relative change is $-1$ and a positive relative change means that PULSE is better.
In (ref), further illustrations of e.g. the relative change in mean bias and variance of the estimators are presented. We also conducted the simulations for setups with combinations of $\gamma\in\{-1,0\}$, components of $\bar \xi$ chosen negatively, with random flipped sign in each coordinate and for negative $\rho$ (not shown but available in the folder 'Plots' in the code repository). The results with respect to MSE are similar to those shown in (ref), while the bias comparison changes depending on the setup.
We observe that there are settings, in which the PULSE is superior to TSLS, Fuller(1) and Fuller(4) in terms of MSE. This is particularly often the case in weak instrument settings ($\hat E_N(G_n)<10$) for low confounding strength $(\rho\leq 0.2)$. Furthermore, as we tend towards the weakest instrument setting considered, we also see a gradual shift in favour of PULSE for higher confounding strengths. In these settings with weak instruments and low confounding we also see that OLS is superior to the PULSE in terms of MSE. However, for large confounding setups PULSE is superior to OLS in terms of both bias and MSE and this superiority increases as the instrument strength increases. The PULSE is generally more biased than the Fuller and TSLS estimators but less biased than OLS. However, in the settings with weak instruments and low confounding the bias of PULSE and OLS is comparable. In summary, the PULSE is in these settings more biased but its variance is so small that it is MSE superior to the Fuller and TSLS estimators.
We now consider three classical instrumental variable applications (see albouy2012colonial and buckles2013season for discussions on the underlying assumptions).
We have applied the different estimators OLS, TSLS, PULSE, and Fuller to the classical data sets acemoglu2001colonial, angrist1991does and card1993using. All models considered in angrist1991does and card1993using, where we estimate the effect on years of education on wages, using quarter of birth and proximity to colleges as instruments, respectively, the OLS estimates are not rejected by our test statistic and PULSE outputs the OLS estimates; see (ref) for futher details. This may be either due to weak endogeneity (weak confounding), or that the test has insufficient power to reject the OLS estimates due to either weak instruments or severe over-identification.
The dataset of acemoglu2001colonial consists of 64 observations, each corresponding to a different country for which mortality rate estimates encountered by the first European settlers are available. The endogenous target of interest is log GDP per capita (in 1995). The main endogenous regressor in the dataset is an index of expropriation protection (averaged over 1985--1995), i.e., protection against expropriation of private investment by the respective governments. The average expropriation protection is instrumented by the settler mortality rates. We consider eight models M1--M8 which correspond to the models presented in column (1)--(8) in Table 4 of acemoglu2001colonial. Model M1 is given by the reduced form structural equations
where avexpr is the average expropriation protection, em4 is the settler mortality rates, $\mu_1$ and $\mu_2$ are intercepts and $U_1$ and $U_2$ are possibly correlated, unobserved noise variables. In model M2 we additionally introduce an included exogenous regressor describing the country latitude. In model M3 and M4 we fit model M1 and M2, respectively, on a dataset where we have removed Neo-European countries, Australia, Canada, New Zealand and the United States. In model M5 and M6 we fit model M1 and M2, respectively, on a dataset where we have removed observations from the continent of Africa. In model M7 and M8 we again fit model M1 and M2, respectively, but now also include three exogenous indicators for the continents Africa, Asia and other.
(ref) shows the OLS and TSLS estimates (which replicate the values from the study), as well as the Fuller(4) and PULSE estimates for the linear effect of the average expropriation protection on log GDP. In model M1, for example, we see that the PULSE estimate suggests that the average expropriation risk linear effect on log GDP is 0.6583 which is 26% larger than the OLS estimate but 34% smaller than TSLS estimate. In models M5--M8, the OLS estimates are not rejected by the Anderson-Rubin test, so the PULSE estimates coincide with the OLS estimates.
We can also use this example to illustrate the robustness property of K-class estimators; see (ref). Even though interventional data are not available, we can consider the mean squared prediction error when holding out the observations with the most extreme values of the instrument. Depending on the degree of generalization, we indeed see that the PULSE and Fuller tend to outperform OLS or TSLS in terms of mean squared prediction error on the held out data; see (ref) for further details.
We have proved that a distributional robustness property similar to the one shown for anchor regression AnchorRegression fully extends to general K-class estimators of possibly non-identifiable structural parameters in a general linear structural equation model that allows for latent endogenous variables. We have further proposed a novel estimator for structural parameters in linear structural equation models. This estimator, called PULSE, is derived as the solution to a minimization problem, where we seek to minimize mean squared prediction error constrained to a confidence region for the causal parameter. Even though this region is non-convex, we have shown that the corresponding optimization problem allows for a computationally efficient algorithm that approximates the above parameter with arbitrary precision using a simple binary search procedure. In the under-identified setting, this estimator extends existing work in the machine learning literature that considers invariant subsets or the best predictive sets among them: PULSE is applicable even in situations when no invariant subsets exist. We have proved that this estimator can also be written as a K-class estimator with data-driven $\kappa$-parameter, which lies between zero and one. Simulation experiments show that in various settings with weak instruments and weak confounding, PULSE outperforms other estimators such as the Fuller(4) estimator. We thus regard PULSE as an interesting alternative for estimating causal effects in instrumental variable settings. It is easy to interpret and automatically provides the user feedback in case that the OLS is accepted (which may be an indication that the instruments are too weak) or that the TSLS is outside the acceptance region (which may indicate a model misspecification). We have applied the different estimators to classical data sets and have seen that, indeed, K-class estimators tend to be more distributionally robust than OLS or TSLS.
There are several further directions that we consider worthwhile investigating. This includes better understanding of finite sample properties and for the identified setups, the study of loss functions other than MSE. It would be helpful, in particular with respect to real world applications, to understand to which extent similar principles can be applied to models allowing for a time structure of the error terms. We believe that the simple primal form of PULSE could make it applicable for model classes that are more complex than linear models ChristiansenTPAMI. Our procedure can be combined with other tests and it could furthermore be interesting to find efficient optimization procedures for tests that are robust with respect to weak instruments, such as Kleibergen's K-statistic Kleibergen2002, for example. In an under-identified setting, the causal parameters are not identified but the solutions obtained by optimizing predictability under invariance might be promising candidates for models that generalize well to distributional shifts.
We are grateful to Trine Boomsma, Peter B\"uhlmann, Rune Christiansen, Steffen Lauritzen, Nicolai Meinshausen, Whitney Newey, Cosma Shalizi, and Nikolaj Thams for helpful discussions. We thank the editor and two anonymous referees for helpful and constructive comments. MEJ and JP were supported by the Carlsberg Foundation; JP was, in addition, supported by a research grant (18968) from VILLUM FONDEN.