The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
106,458 characters
Distributional robustness of K-class estimators and the PULSE
\title{Distributional robustness of K-class\\estimators and the PULSE\footnote{Published in The Econometrics Journal: \url{https://doi.org/10.1093/ectj/utab031}}}
\author{Martin Emil Jakobsen\\Department of Mathematical Sciences\\University of Copenhagen\\[email removed] \and Jonas Peters\\Department of Mathematical Sciences\\University of Copenhagen\\[email removed]}
\maketitle
\begin{abstract}\normalsize While causal models are robust in that they are prediction optimal under arbitrarily
strong interventions,
they may not be optimal when the interventions are bounded.
We prove that the classical K-class estimator satisfies such optimality by establishing a connection between
K-class estimators and anchor regression.
This connection
further motivates a novel estimator in
instrumental variable settings
that minimizes
the mean squared prediction
error subject to the constraint that the estimator lies in an
asymptotically valid confidence region of the causal coefficient. We call
this estimator PULSE (p-uncorrelated least squares estimator), relate it to work on invariance,
show
that it can be computed efficiently as a data-driven K-class estimator, even though the underlying optimization problem is non-convex, and prove consistency.
We evaluate
the estimators on real data and
perform simulation experiments illustrating
that PULSE suffers from less variability. There are several settings including
weak instrument settings, where it outperforms other estimators.
\end{abstract}
\textbf{Keywords:} Causality, distributional robustness, instrumental variables
\crefname{section}{section}{sections}
\Crefname{section}{Section}{Sections}
\crefname{appendix}{appendix}{appendices}
\Crefname{appendix}{Appendix}{Appendices}
\section{Introduction}
Learning causal parameters from data has been
a key challenge in many scientific fields
and has been a long-studied problem in
econometrics
\citep[e.g.][]{ goldberger1972structural, koopmans1953studies,wold1954causality}.
Many years after the groundbreaking work by
\citet{Fisher35} and \citet{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
\citep[e.g.,][]{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{sec:summary}.
\subsection{Related work}
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 \citep[][]{Pearl2009},
propensity score matching \citep[][]{Rosenbaum83},
inverse probability weighting \citep[][]{Horvitz1952} or
G-computation \citep[][]{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
\citep[e.g.,][]{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 \citet{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 \citep[][]{mariano1972existence}.
\citet{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 \citet{anderson1949estimation}.
\citet{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.\ \citet{MARIANO1975}, \citet{mariano2001simultaneous}. While the LIML does not have moments of any order, it shares the same asymptotic normal distribution with TSLS.
Based on simulation studies,
\citet{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,
\citet{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
\citep[][]{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.
\citet{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$.
\citet{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.
\citet{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.\ \citet{Spirtes2000}.
Score-based methods
assume a certain statistical model and optimize (penalized) likelihood scores; see, e.g.\ \citet{Chickering2002}.
Some methods exploit a simple form of causal assignments, such as
additive noise
(e.g., \citealp{Peters2014jmlr}, and \citealp{Shimizu2006})
and others are based on exploiting invariance statements \citep[e.g.,][]{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.\ \citet{Claassen2013} and \citet{Spirtes1995}.
Recent works in the fields of machine learning and computational statistics \citep[e.g.][]{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 \citep[][]{Lauritzen1996},
for example,
but the underlying fundamental principle
has been discussed already several decades ago;
most prominently
using the terms
`autonomy' or `modularity' (\citealp{Haavelmo1944}, and \citealp{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
\Cref{sec:simsem}.
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,
\begin{align}
\label{eq:optimF}
\operatorname*{arg\,min}_\gamma \sup_{F\in \cF} E_F\left[ (Y-\gamma^\intercal X)^2\right].
\end{align}
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
\citep[][]{MeinshausenMiniMax}.
But estimators of~\Cref{eq:optimF}
may be available even
if we do not observe data
from each distribution in $F$.
The true causal coefficient $\gamma_0$, for example,
minimizes \Cref{eq:optimF}
when $\cF$ is the set of all possible
(hard)
interventions on $X$ \citep[e.g.,][]{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 \citet{AnchorRegression}, see Section~\ref{sec:AR} 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
\Cref{eq:optimF}
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 \citep[][]{Pfister2019stab, rojas2018invariant,Magliacane2018} search over different subsets $S$ and
aim to choose sets that are both invariant and predictive.
\subsection{Summary and contributions} \label{sec:summary}
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.
\subsubsection{Distributional robustness of K-class estimators.}
\label{sec:intro:robustness}
In \Cref{SEC:ROBUSTNESS} 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
\begin{align*}
\begin{bmatrix}
Y & X^\intercal & H^\intercal \end{bmatrix}
:= \begin{bmatrix}
Y &X^\intercal & H^\intercal \end{bmatrix} B + A^\intercal M+ \varepsilon^\intercal,
\end{align*}
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
\begin{align*}
\begin{bmatrix}
Y_v & X_v^\intercal & H_v^\intercal \end{bmatrix} := \begin{bmatrix}
Y_v & X_v^\intercal & H_v^\intercal \end{bmatrix} B + v^\intercal M+ \varepsilon.
\end{align*}
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 \Cref{sec:simsem}
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
\begin{align*}
\fY = \fX \gamma_0 + \fA \beta_0 + \fH \eta_0 + \bm{\varepsilon}_Y = \fX \gamma_0 + \fA \beta_0 + \tilde{\fU}_Y.
\end{align*}
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
\begin{align*}
\hat{\alpha}_{\text{K}}^n(\kappa) = (\fZ_*^\intercal (I-\kappa P_\fA^\perp)\fZ_*)^{-1} \fZ_*^\intercal(I-\kappa P_\fA^\perp)\fY,
\end{align*}
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 \cite{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,
\begin{align*}
\hat{\alpha}_{\text{K}}^n(\kappa) \overset{P}{\underset{n\to\infty }{\longrightarrow}} \operatorname*{arg\,min}_{\alpha}\sup_{v\in C(\kappa)} E^{\mathrm{do}(A:=v)}\left[ (Y - \alpha^\intercal Z_* )^2 \right],
\end{align*}
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 \Cref{fig:DistRobustness} in \Cref{app:DistributionalRobustness}).
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.\ \cite{dhrymes1974}.
\subsubsection{The PULSE estimator} \label{sec:summarypuls}
\Cref{SEC:PULSE} 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
\begin{align*}
\hat{\alpha}^n_{\mathrm{K}}( \lambda^\star_n/(1+\lambda^\star_n) ) = \operatorname*{arg\,min}_{\alpha} l_{\mathrm{OLS}}^n(\alpha) + \lambda^\star_n l_{\mathrm{IV}}^n(\alpha) ,
\end{align*}
where
\begin{align*}
\lambda^\star_n
:= \inf \left\{ \lambda > 0\,:\, \begin{tabular}{c}
\text{testing }\text{Corr}$(A,Y-Z\hat{\alpha}_{\text{K}}^n(\lambda/(1+\lambda)))=0$ \\
\text{yields a }p\text{-value }$ \geq p_{\min}$
\end{tabular}
\right\},
\end{align*}
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 \cite{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
\begin{align*}
\hat{\alpha}_{\text{Pr}}^n(t) := \begin{array}{ll}
\operatorname*{arg\,min}_{\alpha} & l_{\mathrm{OLS}}^n(\alpha) \\
\text{subject to} & l_{\mathrm{IV}}^n(\alpha)\leq t.
\end{array}
\end{align*}
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
\begin{align*}
\begin{array}{ll}
\text{argmin}_\alpha & l_{\mathrm{OLS}}^n(\alpha;\fY,\fZ) \\
\text{subject to} & \alpha \in \cA_n(1-p_{\min}),
\end{array}
\end{align*}
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 \citep[see, e.g.,][]{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 \citep[e.g.][]{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 \Cref{app:underidentifiedexperiment} 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.\ \citet{hahn2004estimation} and \citet{stock2002survey}.
Implementation of PULSE and code for experiments (R) are available on GitHub.\footnote{
\url{https://github.com/MartinEmilJakobsen/PULSE}}
\section{Robustness properties of K-class estimators} \label{SEC:ROBUSTNESS}
\setcounter{equation}{0}
In this section
we consider
K-class
estimators (\citealp{theil1958economic}, and \citealp{nagar1959bias})
and show a connection with anchor regression of \citet{AnchorRegression}.
In \Cref{sec:KclassInIVModelCoincidesWithAr} we establish the connection
in models where we use \textit{a priori} information that there are no included exogenous variables in the target equation of interest.
In \Cref{sec:KclassAsPenalizedRegression} we then show that general
K-class estimators can be written as the solution to a penalized regression problem. In \Cref{sec:intervrobustnessofKclass} we utilize this representation and the ideas of \citet{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 \Cref{sec:RobustnessProofs}.
\subsection{Setup and assumptions} \label{sec:SetupAndAssumptionsRobustness}
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)
\begin{align} \label{ARModel}
\begin{bmatrix}
Y & X^\intercal & H^\intercal \end{bmatrix}
:= \begin{bmatrix}
Y &X^\intercal & H^\intercal \end{bmatrix} B + A^\intercal M+ \varepsilon^\intercal,
\end{align}
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
\begin{align} \label{eq:StructuralEquationOfInterest}
\fY = \fX \gamma_0 + \fA \beta_0 + \fH \eta_0 + \bm{\varepsilon}_Y =
\fZ \alpha_0 + \tilde{\fU}_Y,
\end{align}
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
\Cref{eq:StructuralEquationOfInterest}
for consistent estimation. In the remainder of this work, we disregard any \textit{a priori} (non-sample)
information not concerning the target equation.
The question of identifiability of $\alpha_0$ has been studied extensively \citep[][]{frisch38,Haavelmo1944,measuringtheequationsystems} and more recent overviews can be found in, e.g.,
\citet{Didelez2010}, \citet{fisher1966identification}, and \citet{greene2003econometric}.
We will use the following assumptions concerning the structure of the SEM:
\begin{assumption}[Global assumptions] \label{ass:global}
\begin{enumerate*}[label=(\alph*),ref=\ref{ass:global}.(\alph*)] \justifying
\item
$(Y,X,H,A)$ is generated in accordance with the SEM in \Cref{ARModel};
\label{ass:linearSEM}
\item $\rho(B)<1$ where $\rho(B)$ is the spectral radius of $B$; \label{ass:SpectralRadiusOfBLessThanOne}
\item $\varepsilon$ has jointly independent marginals $\varepsilon_1, \ldots, \varepsilon_{d+1+r}$;\label{ass:epIndependentMarginals}
\item $A$ and $\varepsilon$ are independent; \label{ass:AindepEp}
\item No variable in $Y$, $X$ and $H$ is an ancestor of $A$, that is, $A$ is exogenous;
\label{ass:Aexogenous}
\item $E[\|\varepsilon\|^2_2]$, $E[\|A\|^2_2] < \infty$; \label{ass:SecondMomentEp}
\label{ass:SecondMomentA}
\item $E[\varepsilon] =0$. \label{ass:EpMeanZero}
\item $\text{Var}(A) \succ 0$, i.e., the variance matrix of $A$ is positive definite; \label{ass:VarianceOfAPositiveDefinite}
\item $\fA^\intercal \fA$ is almost surely of full rank; \label{ass:AtAfullRank}
\end{enumerate*}
\end{assumption}
\begin{assumption}[Finite sample assumptions] \label{ass:finiteass}
\begin{enumerate*}[label=(\alph*),ref=\ref{ass:finiteass}.(\alph*)]
\item $\fZ_*^\intercal \fZ_*$ is almost surely of full rank; \label{ass:ZtZfullRank}
\item $\fA^\intercal \fZ_*$ is almost surely of full column rank. \label{ass:AtZfullRank}
\item $\fX^\intercal \fX$ is almost surely of full rank; \label{ass:XtXfullRank}
\end{enumerate*}
\end{assumption}
\begin{assumption}[Population assumptions] \label{ass:popass}
\begin{enumerate*}[label=(\alph*),ref=\ref{ass:popass}.(\alph*)]
\item $\text{Var}(Z_*) \succ 0$, i.e., the variance matrix of $Z_*$ is positive definite;\label{ass:VarianceOfZPositiveDefinite}
\item $E[AZ_*^\intercal]$ is of full column rank. \label{ass:EAZtFullColumnRank}
\end{enumerate*}
\end{assumption}
We will henceforth assume that \Cref{ass:global} 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 \Cref{ass:AtAfullRank} and Assumption \Cref{ass:XtXfullRank} hold if $X$ and $A$ have density with respect to Lebesgue measure, which in turn is guaranteed by Assumption \Cref{ass:AindepEp} if $A$ and $\varepsilon$ have density with respect to Lebesgue measure.
Assumption \Cref{ass:VarianceOfAPositiveDefinite} and \Cref{ass:AtAfullRank} implies that the instrumental variable objective functions introduced below is almost surely well-defined and Assumption \Cref{ass:XtXfullRank} yields that the ordinary least square solution is almost surely well-defined. Assumption \Cref{ass:SecondMomentA,ass:SecondMomentEp} implies that $Y,X$ and $H$ all have finite second moments.
For Assumption \Cref{ass:EAZtFullColumnRank} and \Cref{ass:AtZfullRank} it is necessary that $q\geq \mathrm{dim}(Z_*)$, i.e., that the setup must be just- or over-identified; see \Cref{sec:SetupAndAssumptionsPULSE} below.
\subsection{Distributional robustness of anchor regression} \label{sec:AR}
\citet{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
\begin{align} \label{eq:lossPop}
l_{\mathrm{OLS}}(\gamma ;Y,X) &:= E(Y-\gamma^\intercal X)^2, \\
l_{\mathrm{IV}}(\gamma ;Y,X,A)
&:=E(A(Y-\gamma^\intercal X))^\intercal E(A A^\intercal)^{-1} E(A(Y- \gamma^\intercal X)), \notag \\ \label{eq:lossEmp}
l^n_{\mathrm{OLS}}(\gamma;\fY,\fX) &:= n^{-1} (\fY-\fX \gamma)^\intercal (\fY-\fX \gamma),\\
l^n_{\mathrm{IV}}(\gamma;\fY,\fX, \fA) &:=n^{-1} (\mathbf{Y}-\mathbf{X}\gamma)^\intercal P_\fA (\mathbf{Y}-\mathbf{X}\gamma),
\end{align}
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
\begin{align}
\gamma_{\mathrm{AR}}(\lambda)&:=
\operatorname*{arg\,min}_{\gamma\in \mathbb{R}^d} \{l_{\text{OLS}}(\gamma)+\lambda l_{\text{IV}}(\gamma) \}, \quad
\label{eq:EmpiricalARestimator}
\hat{\gamma}_{\mathrm{AR}}^n(\lambda):=
\operatorname*{arg\,min}_{\gamma\in \mathbb{R}^d} \{l_{\text{OLS}}^n(\gamma)+\lambda l_{\text{IV}}^n (\gamma)\}.
\end{align}
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 \Cref{ass:XtXfullRank}, 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
\begin{align} \label{eq:ARsolution}
\hat{\gamma}_{\text{AR}}^n(\lambda) = [\fX^\intercal (I+\lambda P_\fA)\fX]^{-1}\fX^\intercal(I+\lambda P_\fA)\fY,
\end{align}
which follows from solving the normal equation of \Cref{eq:EmpiricalARestimator}.
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 \citet{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.$
\subsection{Distributional robustness of K-class estimators} \label{sec:RobustnessOfKclass}
We now introduce the limited information estimators known as K-class estimators (\citealp{theil1958economic}, and \citealp{nagar1959bias}) used
for single equation inference.
Suppose that we are given non-sample information about which components of $\gamma_0$ and $\beta_0$, of \Cref{eq:StructuralEquationOfInterest}, 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
\begin{align} \label{eq:KclassEstimatorInSingleLine}
\hat{\alpha}_{\text{K}}^n(\kappa;\fY,\fZ_*,\fA) = (\fZ_*^\intercal (I-\kappa P_\fA^\perp)\fZ_*)^{-1} \fZ_*^\intercal(I-\kappa P_\fA^\perp)\fY,
\end{align}
where
$I-\kappa P_\fA^\perp = I- \kappa(I-P_\fA) = (1-\kappa)I + \kappa P_\fA$.
Comparing \Cref{eq:ARsolution,eq:KclassEstimatorInSingleLine}
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.
\subsubsection{K-class estimators in models with no included exogenous variables}
\label{sec:KclassInIVModelCoincidesWithAr}
Assume that,
in addition to \Cref{ass:global}, 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:
\begin{align} \label{eq:robIV}
\gamma_{\text{K}}(\kappa) &= \gamma_{\mathrm{AR}} \left( \frac{\kappa }{1-\kappa} \right) = \operatorname*{arg\,min}_{\gamma\in \mathbb{R}^d} \sup_{v\in C(\kappa/(1-\kappa))} E^{\mathrm{do}(A:=v)} \left[ \left( Y-\gamma^\intercal X\right) ^2 \right],
\end{align}
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 \citet{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.\ \citet{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~\Cref{eq:robIV}
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.
\subsubsection{The K-class estimators as penalized regression estimators} \label{sec:KclassAsPenalizedRegression}
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 \citet{mcdonald1977k} building upon previous work of \citet{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,
\begin{align} \label{KclassLossFunctionPop}
l_{\mathrm{K}}(\alpha;\kappa,Y,Z_*,A)&= (1-\kappa)l_{\mathrm{OLS}}(\alpha;Y,Z_*) + \kappa l_{\mathrm{IV}}(\alpha;Y,Z_*,A), \\ \label{KclassLossFunctionEmp}
l_{\mathrm{K}}^n(\alpha;\kappa,\fY,\fZ_*,\fA)&= (1-\kappa)l_{\mathrm{OLS}}^n(\alpha;\fY,\fZ_*) + \kappa l_{\mathrm{IV}}^n(\alpha;\fY,\fZ_*,\fA).
\end{align}
\begin{restatable}[]{proposition}{PenalizedKClassSolutionUniqueAndExists}
\label{lm:PenalizedKClassSolutionUniqueAndExists}
Consider one of the following scenarios: 1) $\kappa <1$ and Assumption \Cref{ass:ZtZfullRank} holds, or 2) $\kappa = 1$ and Assumption \Cref{ass:AtZfullRank} holds.
The estimator minimizing the empirical loss function of \Cref{KclassLossFunctionEmp} is almost surely well-defined and coincides with the K-class estimator of \Cref{eq:KclassEstimatorInSingleLine}. That is, it almost surely holds that
\begin{align} \label{eq:kclassminim}
\hat{\alpha}_{\mathrm{K}}^n(\kappa;\fY,\fZ_*,\fA) =\operatorname*{arg\,min}_{\alpha\in \mathbb{R}^{d_1+q_1}} l_{\mathrm{K}}^n(\alpha;\kappa,\fY,\fZ_*,\fA).
\end{align}
\end{restatable}
Assuming $\kappa\not =1$, we can rewrite \Cref{eq:kclassminim} to
\begin{align}
\label{eq:KclassLossFunctionAsPenalizedOLS}
\hat{\alpha}_{\text{K}}^n(\kappa;\fY,\fZ_*,\fA) = \operatorname*{arg\,min}_{\alpha\in \mathbb{R}^{d_1+q_1}} \{ l^n_{\mathrm{OLS}}(\alpha;\fY,\fZ_*) + \frac{\kappa}{1-\kappa} l^n_{\mathrm{IV}}(\alpha;\fY,\fZ_*,\fA) \}.
\end{align}
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 \Cref{ass:finiteass}
with the
corresponding population \Cref{ass:popass}, we get that the minimization estimator of the empirical loss function of \Cref{KclassLossFunctionEmp} 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.
\begin{restatable}[]{proposition}{PopulationPenalizedKClassSolutionUniqueAndExists}
\label{lm:PopulationPenalizedKClassSolutionUniqueAndExists}
Consider one of the following scenarios: 1) $\kappa \in [0,1)$ and Assumption \Cref{ass:VarianceOfZPositiveDefinite} holds, or 2) $\kappa = 1$ and Assumption \Cref{ass:EAZtFullColumnRank} holds.
It holds that $(\hat{\alpha}_{\text{K}}^n(\kappa;\fY,\fZ_*,\fA))_{n\geq 1}$ is an asymptotically well-defined sequence of estimators.
Furthermore, the sequence consistently estimates the well-defined population K-class estimand. That is,
$$
\hat{\alpha}_{\mathrm{K}}^n(\kappa;\fY,\fZ_*,\fA) \overset{P}{\underset{n\to\infty }{\longrightarrow}} \alpha_{\mathrm{K}}(\kappa;Y, Z_*,A) := \operatorname*{arg\,min}_{\alpha\in \mathbb{R}^{d_1+q_1}}l_{\mathrm{K}}(\alpha;\kappa,Y,Z_*,A).
$$
\end{restatable}
\subsubsection{Distributional robustness of general K-class estimators}
\label{sec:intervrobustnessofKclass}
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.
\begin{restatable}[]{theorem}{TheoremIntRobustKclass}
\label{sthm:TheoremIntRobustKclas}
Let \Cref{ass:global} hold.
For any fixed $\kappa \in[0,1)$ and $Z_*=(X_*, A_*)$ with $X_*\subseteq X$ and $A_* \subseteq A$, we have, whenever the population K-class estimand is well-defined, that
\begin{align*}
\alpha_{\mathrm{K}}(\kappa;Y,Z_*,A)&= \operatorname*{arg\,min}_{\alpha\in \mathbb{R}^{d_1+q_1}}\sup_{v\in C(\kappa)} E^{\mathrm{do}(A:=v)}\left[ (Y - \alpha^\intercal Z_* )^2 \right],
\end{align*}
where $
C(\kappa) := \left\lbrace v:\Omega \to \mathbb{R}^{q}: \mathrm{Cov}(v,\varepsilon)=0, \, E [ vv^\intercal ] \preceq \frac{1}{1-\kappa} E[AA^\intercal] \right\rbrace$.
\end{restatable}
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{sec:intro:robustness} and \Cref{sec:simsem})
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 \citet{AnchorRegression}.
Our
proof
follows similar arguments but
additionally allows for $A_* \not = \emptyset$.
The property in Theorem~\ref{sthm:TheoremIntRobustKclas}
has a decision-theoretic interpretation (see \citet{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{sthm:TheoremIntRobustKclas}, 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{sthm:TheoremIntRobustKclas}
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 \Cref{rm:ModelMispecification} in \Cref{app:AddRemarks}.
\section{The p-uncorrelated least square estimator} \label{SEC:PULSE}
\setcounter{equation}{0}
We now introduce
the p-uncorrelated least square estimator (PULSE).
As discussed in \Cref{sec:summary}, 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 \Cref{ass:EAZtfullrank} introduced below, also known as the rank condition for identification; \citet{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
\begin{align} \label{eq:PULSEfirstEQ}
\hat{\alpha}^n_{\mathrm{PULSE}}(p_{\min}) = \begin{array}{ll}
\text{argmin}_\alpha & l_{\mathrm{OLS}}^n(\alpha;\fY,\fZ) \\
\text{subject to} & \text{p-value}(T(\alpha;\fY,\fZ,\fA)) \geq p_{\min},
\end{array}
\end{align}
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{sec:summarypuls} and \Cref{app:underidentifiedexperiment}).
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 \citep[][]{anderson1949estimation}.
While the objective function in \Cref{eq:PULSEfirstEQ}
is quadratic in $\alpha$, the resulting constraint is,
in general, non-convex.
In Section~\ref{sec:DualPULSE}, 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 \Cref{sec:SetupAndAssumptionsPULSE}, we briefly introduce the setup and assumptions.
In \Cref{sec:VanishCorr}, we specify a class of asymptotically consistent tests for $\cH_0(\alpha)$.
In \Cref{sec:PULSEdefi} we formally define the PULSE estimator.
In \Cref{sec:PrimalPULSE}, 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 \Cref{sec:DualPULSE}, 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 \Cref{sec:RemainingProofsOfSecPULSE} unless stated otherwise.
\subsection{Setup and assumptions}\label{sec:SetupAndAssumptionsPULSE}
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 \Cref{ARModel}. 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 \Cref{SEC:ROBUSTNESS},
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 \Cref{ass:global}
from \Cref{sec:SetupAndAssumptionsRobustness} still hold.
Furthermore, we will
make use of the following situational assumptions
\begin{assumption} \label{ass:AIndepUYandMeanZeroA}
\begin{enumerate*}[label=(\alph*),ref=\ref{ass:AIndepUYandMeanZeroA}.(\alph*)]
\item $A \protect\mathpalette{\protect\independenT}{\perp} U_Y$; \label{ass:AIndepUy}
\item $E[A]=0$. \label{ass:MeanZeroA}
\end{enumerate*}
\end{assumption}
\begin{assumption}
$\varepsilon$ has non-degenerate marginals.\label{ass:NonDegenYNoise}
\end{assumption}
\begin{assumption} \label{ass:ZtZfullrankandAtZfullrank}
\begin{enumerate*}[label=(\alph*),ref=\ref{ass:ZtZfullrankandAtZfullrank}.(\alph*)]
\item $\fZ^\intercal \fZ$ is of full rank; \label{ass:ZtZfullrank}
\item $\fA^\intercal \fZ$ is of full rank. \label{ass:AtZfullrank}
\end{enumerate*}
\end{assumption}
\begin{assumption}
$[\fZ \, \, \fY]$ is of full column rank. \label{ass:ZYfullcolrank}
\end{assumption}
\begin{assumption}
$E[AZ^\intercal]$ is of full rank. \label{ass:EAZtfullrank}
\end{assumption}
Assumption \Cref{ass:AIndepUy} 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 \Cref{ass:MeanZeroA} 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 \Cref{rm:AssumptionMeanZero}.
Assumption \Cref{ass:ZtZfullrank} ensures that K-class estimators for $\kappa < 1$ are
well-defined, regardless of the over-identification degree.
In the under-identified setup, Assumption \Cref{ass:AtZfullrank} 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$.
\Cref{ass:ZYfullcolrank} 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. \Cref{ass:NonDegenYNoise,ass:EAZtfullrank} ensure that various limiting arguments are valid. In the just- and over-identified setup \Cref{ass:EAZtfullrank} is known as the rank condition for identification.
\subsection{Testing for vanishing correlation} \label{sec:VanishCorr}
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 \Cref{ass:ZYfullcolrank} holds we can define $T_n^c:\mathbb{R}^{d_1+q_1} \to \mathbb{R}$ by
\begin{align*}
T_n^c(\alpha) := c(n) \frac{l_{\text{IV}}^n(\alpha)}{l_{\text{OLS}}^n(\alpha)} =
c(n) \frac{\|P_\fA (\fY- \fZ \alpha) \|_2^2}{\|\fY- \fZ \alpha\|_2^2},
\end{align*}
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.
\begin{restatable}[Level and power of the test]{lemma}{TheoremTestingVanishingCorr}
\label{prop:TestingVanishingCorr}
Let \Cref{ass:AIndepUYandMeanZeroA,ass:NonDegenYNoise,ass:ZYfullcolrank} hold
and assume that $c(n) \sim n$ as $n\to\infty$. For any $p\in (0,1)$
and any fixed $\alpha$, the statistical test rejecting the null hypothesis $\cH_0(\alpha)$ if $T_n^c(\alpha) > Q_{\chi^2_{q}}(1-p),$ has point-wise
asymptotic level $p$ and point-wise asymptotic power of 1 against all alternatives as $n \rightarrow \infty$.
\end{restatable}
\begin{remark} \label{rm:AssumptionMeanZero}
\textnormal{Assumption \Cref{ass:MeanZeroA}, $E[A]=0$, is important for the test statistic to be asymptotic pivotal under the null hypothesis,
that is,
to ensure that the asymptotic distribution of $T_n^c(\alpha)$
does not depend on the
model parameters except for $q$.
We can drop this assumption
if we change the null hypothesis to $\cH_0(\alpha):E[A(Y-Z^\intercal\alpha)]=0$ and add the assumption that $E[U_Y]=0$. Furthermore, if we are in the just- or over-identified setup and \Cref{ass:EAZtfullrank} holds, both of these hypotheses are under their respective assumptions equivalent to $\tilde{\cH}_0(\alpha): \alpha=\alpha_0$. That is, the test in \Cref{prop:TestingVanishingCorr} becomes an asymptotically consistent test for the causal coefficient.}
\end{remark}
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 (\citealp{anderson1950asymptotic}).
We make this connection precise
in \Cref{rm:ConnectionToAndersonRubinCI}
in \Cref{app:AddRemarks}.
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.\ \citet{staiger1997instrumental} and \citet{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 \cite{dufour1997some}.
\cite{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$.
\subsection{The PULSE estimator} \label{sec:PULSEdefi}
For any level $p_{\min}\in(0,1)$, we formally define the PULSE estimator of \Cref{eq:PULSEfirstEQ} by letting the feasible
set
be given by the acceptance region $\cA_n(1-p_{\min})$ of $\cH_0(\alpha)$ using the test of \Cref{prop:TestingVanishingCorr}. That is, we consider
\begin{align} \label{eq:PULSE}
\hat{\alpha}^n_{\mathrm{PULSE}}(p_{\min}) := \begin{array}{ll}
\operatorname*{arg\,min}_\alpha & l_{\mathrm{OLS}}^n(\alpha) \\
\text{subject to} & T_n(\alpha) \leq Q_{\chi^2_{q}}(1-p_{\min}).
\end{array}
\end{align}
In general,
this
is a non-convex optimization problem (\citealp{boyd2004convex})
as the constraint function is non-convex,
see the blue contours in
\Cref{fig:LevelsetsTestAndOLS}(left).
From \Cref{fig:LevelsetsTestAndOLS}(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 \Cref{sec:PrimalPULSE}, we prove that this is not a coincidence:
\Cref{eq:PULSE}
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 \Cref{sec:DualPULSE}, 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 \cite{gautier2011high}. In \Cref{rm:Pretest} in \Cref{app:AddRemarks} we briefly discuss the connection to pre-test estimators.
Furthermore, any method for inverting the test, see, e.g., \cite{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).
\subsection{Primal representation of PULSE}
\label{sec:PrimalPULSE}
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
\begin{align} \label{PR.t.n}
\begin{array}{ll}
\text{minimize}_\alpha & l_{\mathrm{OLS}}^n(\alpha;\fY,\fZ) \\
\text{subject to} & l_{\mathrm{IV}}^n(\alpha;\fY,\fZ,\fA)\leq t.
\end{array}
\end{align}
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.
\begin{restatable}[Unique solvability of the primal]{lemma}{LemmaPrimalUniqueSolution} \label{lm:PrimalUniqueSolAndSlatersConditions}
Let \Cref{ass:ZtZfullrankandAtZfullrank} hold. It holds that $\alpha \mapsto l_{\mathrm{OLS}}^n(\alpha)$ and $\alpha \mapsto l_{\mathrm{IV}}^n(\alpha)$ are strictly convex and convex, respectively. Furthermore, for any $t > \inf_{\alpha}l_{\mathrm{IV}}^n(\alpha)$ it holds that the constrained minimization problem (Primal$.t.n$) has a unique solution and satisfies Slater's condition. In the under- and just-identified setup the constraint bound requirement
is equivalent to $t>0$ and in the over-identified setup to $t> l_{\mathrm{IV}}^n(\hat{\alpha}^n_{\mathrm{TSLS}})$, where $\hat{\alpha}^n_{\mathrm{TSLS}}= (\fZ^\intercal P_\fA \fZ)^{-1}\fZ^\intercal P_\fA \fY $.
\end{restatable}
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{lm:PrimalUniqueSolAndSlatersConditions} above).
Slater's condition
will play a role
in \Cref{sec:DualPULSE} 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
\begin{align} \label{eq:PrimalProblemSolutionDef}
\hat{\alpha}_{\text{Pr}}^n(t) := \begin{array}{ll}
\operatorname*{arg\,min}_{\alpha} & l_{\mathrm{OLS}}^n(\alpha) \\
\text{subject to} & l_{\mathrm{IV}}^n(\alpha)\leq t.
\end{array}
\end{align}
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
\begin{align} \label{eq.Def.t.star.p}
t_n^\star(p_{\min}) := \sup \{ t \in ( \inf_{\alpha}l_{\mathrm{IV}}^n(\alpha),l_{\mathrm{IV}}^n(\hat{\alpha}_{\mathrm{OLS}}^n) ] : T_n(\hat{\alpha}_{\mathrm{Pr}}^n(t))\leq Q_{\chi^2_{q}}(1-p_{\min})\}.
\end{align}
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.
\begin{restatable}[Primal representation of PULSE]{theorem}{pPULSESolvesPULSE}
\label{thm:pPULSESolvesPULSE}
Let $p_{\min}\in(0,1)$ and \Cref{ass:ZtZfullrankandAtZfullrank,ass:ZYfullcolrank} hold and assume that $t_n^\star(p_{\min}) >-\infty$. If
$T_n(\hat{\alpha}_{\mathrm{Pr}}^n(t_n^\star(p_{\min})))\leq Q_{\chi^2_{q}}(1-p_{\min})$, then the PULSE problem has a unique solution given by the primal PULSE solution. That is,
$
\hat{\alpha}^n_{\mathrm{PULSE}}(p_{\min}) = \hat{\alpha}_{\mathrm{Pr}}^n(t_n^\star(p_{\min})).$
\end{restatable}
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 \Cref{thm:PULSEpPULSEdPULSEEequivalent}. 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. \Cref{fig:LevelsetsTestAndOLSAndIV}
shows
an example of the equivalence in \Cref{thm:pPULSESolvesPULSE}. \Cref{fig:LevelsetsTestAndOLSAndIV}(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})$.
\begin{figure}[t]
\centering
\includegraphics[width=\linewidth-0pt]{Levelsets_Test_OLS_IV_Combined-eps-converted-to.pdf}
\caption{ Illustrations of the level sets of $l_{\mathrm{OLS}}^n$ (red contours), the proposed test-statistic $T_n$ (blue contours) and $l_{\text{IV}}^n$ (green contours) in a just-identified setup. The example is generated with a two dimensional anchor $A = (A_1,A_2)$, one of which is included, and one included endogenous variable $X$, i.e., $Y = \alpha_1 X + \alpha_2 A_1 + H+ \varepsilon_Y$ with $(\alpha_1,\alpha_2)=(1,1)$. Both illustrations
show level sets from the same setup, but they use different scales.
The black text denotes the level of the test-statistic contours. In this setup, the PULSE constraint bound, the rejection threshold of the test with $p_{\min}=0.05$, is $Q_{\chi^2_{2}}(0.95) \approx 5.99$. The blue level sets of $T_n$ are
non-convex. The sublevel set of the test, corresponding to the acceptance region, is illustrated by the blue area. In the right plot, we see that the smallest level set of $l_{\text{OLS}}^n$ that has a non-empty intersection with the $Q_{\chi^2_{q}}(1-p_{\min})$-sublevel set of $T_n$ is a singleton (black dot, $t^*(p)$). This shows that
in this example
the PULSE problem is solvable and has a unique solution. The $l_{\text{IV}}^n$ level set that intersects this singleton is exactly the $t_n^\star(p_{\min})$-level set of $l_{\text{IV}}^n$, illustrating the statement of \Cref{thm:pPULSESolvesPULSE} in that the primal PULSE
with that choice of $t$
solves the PULSE problem. The black line visualizes the solutions
$\{\hat{\alpha}_{\text{Pr}}^n(t) : t \in D_{\text{Pr} }\}$. The black points and corresponding text labels indicates which constraint bound $t$ yields the specific point. In general, the class of primal solutions
does not coincide with the class of convex combinations of the OLS and the TSLS estimators.}
\label{fig:LevelsetsTestAndOLSAndIV} \label{fig:LevelsetsTestAndOLS}
\end{figure}
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))$ \citep[see also][]{AnchorRegression}.
Theorem~\ref{thm:pPULSESolvesPULSE} 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 \Cref{fig:LevelsetsTestAndOLSAndIV}.
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.\ \cite{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.
\subsection{Dual representation of PULSE} \label{sec:DualPULSE}
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
\begin{align} \label{K.lambda.n}
\begin{array}{ll}
\text{minimize} & l_{\mathrm{OLS}}^n(\alpha) + \lambda l_{\mathrm{IV}}^n(\alpha).
\end{array}
\end{align}
Whenever Assumption \Cref{ass:ZtZfullrank} 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 \Cref{lm:PenalizedKClassSolutionUniqueAndExists}. That is,
\begin{align*}
\hat{\alpha}_{\text{K}}^n(\kappa)
= (\fZ^\intercal (\fI+\lambda P_\fA)\fZ)^{-1} \fZ^\intercal(\fI+\lambda P_\fA)\fY
\end{align*}
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 \Cref{ass:AtZfullrank} we have that the minimum of $l_{\text{IV}}^n(\alpha)$ is attainable (see the proof of \Cref{lm:PrimalUniqueSolAndSlatersConditions}). Hence, let the solution space for the minimization problem $\min_\alpha l_{\text{IV}}^n(\alpha)$ be given by
\begin{align}\label{eq:Miv}
\cM_{\mathrm{IV}} := \operatorname*{arg\,min}_\alpha l_{\mathrm{IV}}^n(\alpha) = \{\alpha \in \mathbb{R}^{d_1+q_1} : l_{\mathrm{IV}}^n(\alpha)= \min_{\alpha'}l_{\mathrm{IV}}^n(\alpha')\}.
\end{align}
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.
\begin{restatable}[K-class estimators and TSLS differ]{lemma}{KclassNotEqualToTwoSLS}
\label{lm:KclassNotEqualToTwoSLS}
Assume that we are in the just- or over-identified setup and $n>q$. Furthermore, assume that $\varepsilon$ has density with respect to Lebesgue measure and that the coefficient matrix $B$ of the SEM in \Cref{ARModel} is lower triangular. If the rank conditions of \Cref{ass:ZtZfullrankandAtZfullrank} hold almost surely, then it
almost surely holds, that all K-class estimators with penalty parameter $\lambda\in[0,\infty)$ differ from the
TSLS estimator, i.e., $\hat{\alpha}_{\mathrm{TSLS}}^n \not\in \{\hat{\alpha}_{\mathrm{K}}^n(\lambda): \lambda \geq 0 \}$.
\end{restatable}
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.
\begin{assumption}
No K-class estimator $\hat{\alpha}_{\mathrm{K}}^n(\kappa)$ with $\kappa\in[0,1)$, is a member of $\cM_{\text{IV}}$. \label{ass:KclassNotInIV}
\end{assumption}
Furthermore, when imposing that \Cref{ass:KclassNotInIV} holds we also have that the K-class estimators differ from each other.
\begin{restatable}[K-class estimators differ]{corollary}{KclassSolutionsDistinct}
\label{cor:KclassSolutionsDistinct}
Let
\Cref{ass:ZtZfullrankandAtZfullrank,ass:KclassNotInIV} hold.
If $\lambda_1,\lambda_2\geq 0$ with $\lambda_1\not = \lambda_2$, then $\hat{\alpha}_{\mathrm{K}}^n(\lambda_1) \not = \hat{\alpha}_{\mathrm{K}}^n(\lambda_2)$.
\end{restatable}
The above corollary is proven as \Cref{cor:KclassSolutionsDistinctApp} in \Cref{sec:SomeProofsOfSecPULSE}.
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)$.
\begin{restatable}[Connecting the primal and dual]{lemma}{EquivalenceBetweenKlikeAndPrimal} \label{lm:EquivalenceBetweenKlikeAndPrimal}
If
\Cref{ass:ZtZfullrankandAtZfullrank,ass:ZYfullcolrank,ass:KclassNotInIV} hold, then both of the following statements hold. \textit{(a)} For any $t \in D_{\mathrm{Pr}}$, there exists a unique $\lambda(t) \geq 0$ such that (Primal$.t.n$) and (Dual$.\lambda(t).n$) have the same unique solution. \textit{(b)} For any $\lambda \geq 0$, there exists a unique $t(\lambda) \in D_{\mathrm{Pr}}$ such that (Primal$.t(\lambda).n$) and (Dual$.\lambda.n$) have the same unique solution.
\end{restatable}
\Cref{lm:EquivalenceBetweenKlikeAndPrimal} 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
\begin{align} \label{eq.Def.lambda.star.p}
\lambda_n^\star(p_{\min}) := \inf\{\lambda \geq 0 : T_n(\hat{\alpha}_{\mathrm{K}}^n (\lambda))\leq Q_{\chi^2_{q}}(1-p_{\min}) \}.
\end{align}
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).
$
\begin{restatable}[Dual representation of PULSE]{theorem}{PrimalDualConnectionPvalConstraint} \label{thm:PULSEpPULSEdPULSEEequivalent}
Let $p_{\min}\in(0,1)$ and \Cref{ass:ZtZfullrankandAtZfullrank,ass:ZYfullcolrank,ass:KclassNotInIV} hold.
If $\lambda_n^\star(p_{\min}) <\infty $, then it holds that $t_n^\star(p_{\min})>-\infty$ and $ \hat{\alpha}_{\mathrm{K}}^n(\lambda_n^\star(p_{\min})) = \hat{\alpha}_{\mathrm{Pr}}^n(t_n^\star(p_{\min})) = \hat{\alpha}_{\mathrm{PULSE}}^n(p_{\min})$.
\end{restatable}
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:
\begin{align*}
\hat{\alpha}_{\mathrm{PULSE}}^n(p_{\min}) &= \hat{\alpha}_{\mathrm{K}}^n(\lambda_n^\star(p_{\min})) = (\fZ^\intercal (\fI+\lambda_n^\star(p_{\min}) P_\fA)\fZ)^{-1} \fZ^\intercal(\fI+\lambda_n^\star(p_{\min}) P_\fA)\fY.
\end{align*}
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 \Cref{sec:BinarySearch}, 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.
\begin{restatable}[Infeasibility of the dual representation]{lemma}{LambdaStarFinite}
\label{lm:LamdaStarFiniteIFF}
Let $p_{\min}\in(0,1)$ and \Cref{ass:ZtZfullrankandAtZfullrank,ass:ZYfullcolrank,ass:KclassNotInIV} hold. In the under- and just-identified setup we have that $\lambda_n^\star(p_{\min})<\infty$. In the over-identified setup it holds that
$
\lambda^\star_n(p_{\min}) < \infty \iff T_n(\hat{\alpha}_{\mathrm{TSLS}}^n)< Q_{\chi^2_q}(1-p_{\min}).
$
This is not guaranteed to hold as the event that $\cA_n(1-p_{\min})= \emptyset$ can have positive probability.
\end{restatable}
Thus, under suitable regularity assumptions \Cref{lm:LamdaStarFiniteIFF} 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.
\subsubsection{Binary search for the dual parameter} \label{sec:BinarySearch}
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.
\begin{restatable}[Monotonicity of the losses and the test statistic]{lemma}{MonotonicityOfTestOfLambda}
\label{lm:OLSandIV_Monotonicity_FnctOfPenaltyParameterLambda}
When Assumption \Cref{ass:ZtZfullrank} holds the maps $
[0,\infty)\ni \lambda \mapsto l_{\mathrm{OLS}}^n(\hat{\alpha}_{\mathrm{K}}^n (\lambda) )$ and $ [0,\infty)\ni \lambda \mapsto l_{\mathrm{IV}}^n(\hat{\alpha}_{\mathrm{K}}^n (\lambda) ) $
are monotonically increasing and monotonically decreasing, respectively. Consequently, if \Cref{ass:ZYfullcolrank} holds, we have that the map
$
[0,\infty)\ni \lambda \mapsto T_n (\hat{\alpha}_{\mathrm{K}}^n (\lambda) )
$
is monotonically decreasing. Furthermore, if \Cref{ass:KclassNotInIV} also holds, these monotonicity statements can be strengthened to strictly decreasing and strictly increasing.
\end{restatable}
The above lemma is proven as \Cref{lm:OLSandIV_Monotonicity_FnctOfPenaltyParameterLambdaApp} in \Cref{sec:SomeProofsOfSecPULSE}. 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 \Cref{lm:TestInAlphaLambdaStarEqualsQuantileApp} in \Cref{sec:SomeProofsOfSecPULSE}.
We propose the binary search algorithm
presented in Algorithm~\ref{Binary.Search.Lambda.Star} in \Cref{app:algo},
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 \Cref{eq:LambdaEquality} in the proof of \Cref{lm:LamdaStarFiniteIFF}. This initialization removes the need for the first while loop in (lines 4--6).
We now prove that Algorithm~\ref{Binary.Search.Lambda.Star}
achieves the required precision and
is asymptotically correct.
\begin{restatable}[]{lemma}{BinarySearchLambdaStarConverges}
\label{lm:BinarySearchLambdaStarConverges}
Let $p_{\min}\in(0,1)$ and \Cref{ass:ZtZfullrankandAtZfullrank,ass:ZYfullcolrank} hold. If $\lambda_n^\star(p_{\min})<\infty$, then $\lambda_n^\star(p_{\min})$ can be approximated with arbitrary precision by the binary search \Cref{Binary.Search.Lambda.Star}, that is, $ \mathrm{Binary.Search}(N,p_{\min}) - \lambda_n^\star(p_{\min}) \to 0,$ as $N\to\infty$.
\end{restatable}
\subsection{Algorithm and consistency} \label{sec:AlgoAndConsistency}
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 \citet{Keane2010} and \citet{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
\begin{align*}
\hat{\alpha}_{\text{PULSE}+}^n(p_{\min}) := \left\{\begin{array}{ll}
\hat{\alpha}_{\mathrm{PULSE}}^n(p_{\min}), & \text{if } T_n(\hat{\alpha}_{\text{TSLS}}^n)<Q_{\chi^2_{q}}(1-p_{\min}) \\
\hat{\alpha}^n_{\text{ALT}}, & \text{otherwise}.
\end{array} \right.
\end{align*}
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 \Cref{alg:2} in \Cref{algo:pulseplus} 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$.
\begin{restatable}[Consistency of PULSE$+$]{theorem}{ConsistencyOfPULSE}
\label{thm:ConsistencyOfPULSE}
Consider the just- or over-identified setup and let $p_{\min}\in (0,1)$. If \Cref{ass:AIndepUYandMeanZeroA,ass:ZtZfullrankandAtZfullrank,ass:ZYfullcolrank,ass:EAZtfullrank,ass:KclassNotInIV} hold almost surely for all $n\in \mathbb{N}$ and $\hat{\alpha}_{\mathrm{ALT}}^n$ consistently estimates $\alpha_0$, then $\hat{\alpha}_{\mathrm{PULSE}+}^n(p_{\min})\stackrel{P}{\longrightarrow} \alpha_0$, when $n\to \infty$.
\end{restatable}
We believe that a similar statement also holds in the under-identified setting, see \Cref{app:underidentifiedexperiment}.
\section{Simulation experiments}
In \Cref{sec:Experiments} 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 \Cref{sec:WeakInst}. Here we give a brief overview of the study and the observations.
\subsection{Distributional robustness}
The
theoretical results on distributional robustness
proved in \Cref{SEC:ROBUSTNESS}
translate to finite data. The experiments of \Cref{app:DistributionalRobustness} 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 \cref{app:underidentifiedexperiment}, we furthermore consider an under-identified setting.
\subsection{Estimating causal effects}
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 \Cref{sec:Experiments} for further details and the multivariate simulation experiments mentioned above.
\subsubsection{Univariate model} \label{sec:SimUnivariate}
We first
compared performance measures of the estimators in a univariate instrumental variable model. As seen in \citet{hahn2002new} and \citet{hahn2004estimation}, we consider structural equation models of the form
\begin{align*}
A := A \in \mathbb{R}^q,\quad
X := A^\intercal \bar \xi + U_X \in \mathbb{R}, \quad
Y := X \gamma + U_Y \in \mathbb{R},
\end{align*}
where $A\sim \mathcal{N}(0,I)$ and $A\protect\mathpalette{\protect\independenT}{\perp} (U_X,U_Y) $ with
$
\begin{psmallmatrix}
U_X \\ U_Y
\end{psmallmatrix}\sim \mathcal{N}\left(\begin{psmallmatrix}
0 \\ 0
\end{psmallmatrix}, \begin{psmallmatrix}
1 & \rho \\
\rho & 1
\end{psmallmatrix}\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 \Cref{sec:ExpResultsEstCausPerformanceMeasures}.
\Cref{fig:HahnExpRMSE} 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 \Cref{sec:WeakInst} 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.
\begin{figure}[h]
\centering\includegraphics[width=\linewidth]{Univariate_RMSE_20200716034113-eps-converted-to.pdf}
\caption{ \normalsize Illustrations of the relative change in RMSE.}\label{fig:HahnExpRMSE}
\end{figure}
In \Cref{sec:AppFigs}, 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 \Cref{fig:HahnExpRMSE}, 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.
\section{Empirical applications}
We now consider three classical
instrumental variable applications (see \citet{albouy2012colonial} and \citet{buckles2013season} for discussions on the underlying assumptions).
\begin{itemize}
\item[\textit{(i)}] ``Does compulsory school attendance affect schooling and earnings?'' by \cite{angrist1991does}. This paper investigates the effects of education on wages. The endogenous effect of education on wages are remedied by instrumenting education on quarter of birth indicators.
\item[\textit{(ii)}] ``Using geographic variation in college proximity to estimate the return to schooling'' by \cite{card1993using}. This paper also investigates the effects of education on wages. In this paper education is instrumented by proximity to college indicator.
\item[\textit{(iii)}] ``The colonial origins of comparative development: An empirical investigation'' by \cite{acemoglu2001colonial}. This paper investigates the effects of extractive institutions (proxied by protection against expropriation) on the gross domestic product (GDP) per capita. The endogeneity of the explanatory variables are remedied by instrumenting protection against expropriation on early European settler mortality rates.
\end{itemize}
We have
applied the different estimators OLS, TSLS, PULSE, and Fuller to the classical data sets
\cite{acemoglu2001colonial}, \cite{angrist1991does} and \cite{card1993using}.
All models considered in \cite{angrist1991does} and \cite{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 \Cref{app:EmpricalApp} 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.
\subsection{\cite{acemoglu2001colonial}} \label{sec:mainColonialApplication}
The dataset of \cite{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 \cite{acemoglu2001colonial}. Model M1 is given by the reduced form structural equations
\begin{align*}
\log \mathrm{GDP} = \text{avexpr}\cdot \gamma + \mu_1 + U_1, \quad \text{avexpr} = \log \mathrm{em4}\cdot \delta + \mu_2 + U_2,
\end{align*}
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.
\Cref{tbl:SettlerMortalityCoefficients} 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.
\begin{table}[htp]
\caption{\label{tbl:SettlerMortalityCoefficients}The estimated return of expropriation protection on log GDP per capita. }
\begin{center}
\begin{tabu}to \textwidth {r c c c c c c c}
\toprule \toprule
Model & OLS & TSLS & FUL & PULSE & Message & Test & Threshold \\
\midrule
M1 &0.5221 & 0.9443 & 0.8584 & 0.6583 & -- & 5.9915 & 5.9915\\
M2 &0.4679 & 0.9957 & 0.8457 & 0.5834 & -- & 7.8147 & 7.8147\\
M3 &0.4868 & 1.2812 & 0.9925 & 0.7429 & -- & 5.9914 & 5.9915\\
M4 &0.4709 & 1.2118 & 0.9268 & 0.6292 & -- & 7.8147 & 7.8147\\
M5 &0.4824 & 0.5780 & 0.5573 & 0.4824 & OLS Accepted & 1.1798 & 5.9915\\
M6 &0.4658 & 0.5757 & 0.5476 & 0.4658 & OLS Accepted & 1.1554 & 7.8147\\
M7 &0.4238 & 0.9822 & 0.7409 & 0.4238 & OLS Accepted & 10.7722 & 11.0705\\
M8 &0.4013 & 1.1071 & 0.7059 & 0.4013 & OLS Accepted & 9.7546 & 12.5916\\
\bottomrule\bottomrule
\end{tabu}
\end{center}
\footnotesize
\textbf{Note:}
Point estimates for the return of expropriation protection on log GDP per capita.
The OLS and TSLS values coincide with the ones shown in
\cite{acemoglu2001colonial}.
The right columns show the values of the test statistic (evaluated in the PULSE estimates) and the test rejection thresholds.
The `--' indicates that OLS is not accepted and TSLS is not rejected.
\end{table}
We can also use this example to
illustrate the robustness property of K-class estimators;
see \Cref{sthm:TheoremIntRobustKclas}.
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 \Cref{sec:Colonial} for further details.
\section{Summary and future work}
We have proved that a distributional robustness property similar to the one shown for anchor regression \citep[][]{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 \citep[see also][]{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 \citep[][]{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.
\section*{Acknowledgements}
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.