Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
125,736 characters · 14 sections · 29 citation commands
A path-sampling method to partially identify causal effects in instrumental variable models
In recent years, a trend in the literature in econometrics and in particular causal inference has been to obtain bounds on quantities of interest in general instrumental variable models via linear programming approaches, often in connection with capacity- and random set theory (e.g. \citet*{balke1994counterfactual}, \citet*{balke1997bounds}, \citet*{manski2007partial}, \citet*{kitamura2018nonparametric}, \citeauthor*{molchanov2005theory} molchanov2005theory, \citeauthor*{beresteanu2011sharp} beresteanu2011sharp, \citeauthor*{galichon2011set} galichon2011set, \citeauthor*{beresteanu2012partial} beresteanu2012partial, \citeauthor*{molchanov2014applications} molchanov2014applications, \citeauthor*{chesher2017generalized} chesher2017generalized, \citeauthor*{russell2019sharp} russell2019sharp). The arguments put forward in favor of these partial identification approaches are higher flexibility and robustness compared to point-identification approaches (\citet*{chesher2017generalized}, \citet*{manski2003partial}). However, these general partial identification approaches suffer from a severe “curse of cardinality” which limits their broader use in general models, as noted in \citet*{beresteanu2012partial}.
The issue is that the number of restrictions on the identified set grows at least exponentially with the number of points in the support of the endogenous variables. Existing methods are intractable in practical settings with high-cardinal endogenous variables\footnote{We say a random variable is “of high cardinality” if it has significantly more than $2$ points in its support, or is continuous. As we show below, even as few as $5$ support points satisfy this definition in the most general instrumental variable setting already.} and have therefore almost exclusively focused on the case of a binary endogenous variable or have made use of the given particular structure of the problem (e.g. \citet*{aguiar2019prices}, \citet*{chesher2017generalized}, \citet*{chiburis2010semiparametric}, \citet*{cheng2006bounds}, \citet*{demuynck2015bounding}, \citet*{hansen1995econometric}, \citet*{honore2006bounds}, \citet*{honore2006bounds2}, \citet*{manski2007partial}, \citet*{manski2014identification}, \citet*{molinari2008partial}, \citet*{norets2013semiparametric}, \citet*{laffers2015bounding}, \citet*{kamat2017identification}, \citet*{torgovitsky2016nonparametric}, \citet*{kitamura2018nonparametric}, \citet*{mogstad2018using}, \citet*{russell2019sharp}, \citet*{tebaldi2019nonparametric}). A practical and generally applicable method that deals with this curse of cardinality has so far been unavailable.
This paper introduces such a method. It extends the linear programming approach to partial identification of functionals of interest to high-cardinality settings, even allowing for continuous variables. The basic idea is to represent the instrumental variable model as a system of stochastic processes indexed by the unobservable heterogeneity, based on the potential outcome notation for causal models \citep*{rubin1974estimating}. This representation as stochastic processes is useful, as it opens the door for the application of all the tools from time-series analysis and stochastic calculus to this setting. A general, in the case of continuous endogenous variables infinite dimensional, linear program is then constructed on the paths of these processes, generalizing the linear programming approach for binary causal models in \citet*{balke1997bounds}. This approach does not require structural assumptions, but allows to incorporate them seamlessly into the model by ruling out paths of the stochastic processes: the idea is that each path denotes a hypothetical participant in the model, so that the stochastic framework approach is a simple way to explicitly model the unobserved heterogeneity in terms of human behavior.
The main contribution of this article is a computational approach to approximately solve these potentially infinite-dimensional linear programs on path spaces. The idea is to sample paths of the stochastic processes and to solve the linear programs on this sample of paths. This introduced randomness is crucial, because it permits the derivation of probabilistic approximation guarantees of theoretical linear program by its sampled counterpart using standard concentration results (\citeauthor*{vapnik1998statistical} vapnik1998statistical, \citeauthor*{wellner2013weak} wellner2013weak). This approach could be helpful in many other settings with optimization problems on an infinite dimensional state space. The key is that one can straightforwardly define a measure on the paths of these processes and use the theory of stochastic processes to analyze the properties, especially in continuous models. In a nutshell, we use randomization to alleviate the “curse of cardinality” in linear programs on path spaces, which is similar in spirit to using randomization to “break the curse of dimensionality” in classical dynamic optimization problems \citep*{rust1997using}. In this article we argue that our approach is an attractive alternative to other sampling approaches such as \citet*{de2004constraint}, because it works on paths of processes instead of sets of inequalities and is therefore more flexible.
As a demonstration of its capabilities, the method estimates bounds on expenditure differences using the 1995/1996 UK family expenditure survey. This problem is well suited for demonstrating the method's capabilities in practice as it is nonlinear with continuous variables and allows to gauge if the program actually obtains reasonable results.\footnote{In particular, the method should produce results which show that food is a necessity- and leisure is a luxury good, as these are well-established economic facts. This application is actually more challenging than a Monte-Carlo approach, as the method needs to replicate known facts on real data under minimal assumptions (see \citeauthor*{advani2019mostly} advani2019mostly for a recent discussion). A priori, it is not even clear that any approach can deliver informative enough bounds to check its validity. The fact that this method does provide informative bounds is a testament to its potential usefulness.} Surprisingly, the method already seems to provide informative nonparametric bounds on household expenditures under the sole assumption that expenditure on goods is continuous with respect to the budget set, corroborating the nonparametric and semi-nonparametric approaches in \citet*{imbens2009identification}, \citet*{blundell2007semi}, \citet*{de2016nonparametric}, and \citet*{song2018nonseparable}, which assume monotonicity or additive separability in the unobservables in the first- or second stage.
The outline of this article is as follows. Section (ref) introduces the setting and the optimization problem at an intuitive level, highlighting three potential settings for applications: randomized controlled trials with imperfect compliance, time-varying treatment effects in program evaluation, and welfare estimation. Section (ref) is the main section of this article and introduces all theoretical results. It contains introduces the stochastic process representation of the causal model formally (Proposition (ref) in section (ref)), states the main result of this article which quantifies the probabilistic approximation of the optimization problems by path sampling (Theorem (ref) in section (ref)), and derives the large sample asymptotics of the optimization programs (Proposition (ref) in section (ref)). Section (ref) deals with the practical implementation and introduces the practical algorithm. Section (ref) contains empirical results: section (ref) contains a brief simulation exercise and section (ref) contains the application to demand estimation. Section (ref) concludes. The appendix contains all proofs.
The basic model in this article is the most general form of an instrumental variable model. This means it makes no a priori assumptions on the functional relations between the variables and leaves the unobserved heterogeneity unrestricted. It is
where, $Y$ is the outcome variable of interest, $X$ is the treatment variable of interest. $X$ is endogenous in the sense that it depends on the unobservable confounder $W$, which also has an influence on $Y$, in addition to the direct effect of $X$ on $Y$. This additional pathway induces a bias when trying to estimate a causal effect of $X$ on $Y$, the classical problem in causal inference and econometrics. A general solution for this is to use an instrumental variable $Z$, which affects the treatment $X$, but is itself exogenous, i.e. independent of the model. $Z$ satisfies full independence of the unrestricted unobserved heterogeneity $W$, denoted by $Z\protect\mathpalette{\protect\independenT}{\perp} W$, as the production functions $g$ and $h$ are unrestricted.\footnote{$Z\protect\mathpalette{\protect\independenT}{\perp} W$ means that $Z$ is independent of $W$, i.e. $P_{Z,W} = P_ZP_{W}$. This model is often written with two separate unobservable variables $V$ and $U$ in the second- and first stage (e.g. \citeauthor*{imbens2009identification} imbens2009identification). This is an equivalent model to (ref) as one can define the two dependent variables $U$ and $V$ on the same probability space and combine them to one variable $W$.}
In the causal inference literature, it is common to represent model (ref) via a directed acyclic graph (\citeauthor*{balke1994counterfactual} balke1994counterfactual, \citeauthor*{pearl1995causal} pearl1995causal, \citeauthor*{balke1997bounds} balke1997bounds) as in Figure (ref).
A third way to represent model (ref) is via the counterfactual notation \citep*{rubin1974estimating}. In this case, the counterfactual (i.e. unobservable) distribution $P_{Y(x)}$ denotes the distribution of the outcome given that $X$ is exogenous, i.e. that $X\protect\mathpalette{\protect\independenT}{\perp} W$. Due to the endogeneity problem $X\not\protect\mathpalette{\protect\independenT}{\perp} W$ depicted in Figure (ref), the observable distribution $P_{Y|X=x}$ does not coincide with $P_{Y(x)}$, so that one cannot simply use the observable distribution to obtain causal effects of interest. However, using the information of the instrument $Z$, a linear programming approach allows us to partially identify functionals of interest such as average treatment effects (ATE) or probabilities of counterfactual events without making structural assumptions on the relations between the variables in question.
The linear programming approach proposed in this article is the natural generalization of \citet*{balke1994counterfactual} and \citet*{balke1997bounds}, which deal with obtaining sharp bounds on a causal effect (such as the average treatment effect) of $X$ on $Y$ when $Y,X,Z$ are binary, to the settings where all $Y$, $X$, and $Z$ are allowed to be continuous.\footnote{See \citet*{russell2019sharp} for an overview of linear programming approaches to identification in instrumental variable models.} It relies on yet another representation of model (ref), which we now introduce in an intuitive way and will address formally in Section (ref).
This representation is based on Rubin's counterfactual notation in the sense that it considers the counterfactual distributions $P_{Y(x)}$ and $P_{X(z)}$ to be the laws of corresponding stochastic processes $Y_x(w)$ and $X_z(w)$ of the first and second stage. Each element $w\in\mathcal{W}$ indexes one path $Y_x(w)$ and $X_z(w)$ of the processes, respectively. Formally, the process $Y_x$ corresponds to $h(x,W)$, and analogous for $X_z$ and $g(z,W)$. This representation allows one to set up a linear program for obtaining bounds on functionals of interest in model (ref):
The programs (ref) are phrased for the general case of potentially continuous $Y$, $X$, and $Z$, anticipating the formal results in Section (ref), but are otherwise perfectly analogous to the optimization program in \citet*{balke1994counterfactual}. They proceed by finding the optimal distribution $P_W$ which maximizes (for an upper bound) or minimizes (for a lower bound) the functional of interest $E_{P_W}[f(Y_{x},x_0)]$ under the restriction that $P_W$ induces processes $Y_x$ and $X_z$ whose induced joint distribution (defined below) \[F_{[Y,X]^*_z}\coloneqq P_W(Y_{X_z}\leq y,X_z\leq x)\] coincides with the joint observable distribution $F_{Y,X|Z=z}$. The optimization is over some set of probability distributions $\mathscr{P}^*(\mathcal{W})$ over the support $\mathcal{W}$ of $W$. The main contribution of this article is a method to probabilistically approximate the solution to optimization problems of the form (ref). We relegate all formal questions about $\mathscr{P}^*(\mathcal{W})$ and the optimization problem more generally to Section (ref).
The functionals in the objective function of (ref) in general can be any map $\phi:\mathscr{P}^*(\mathcal{W})\to\mathbb{R}$ from probability distributions on $\mathcal{W}$ to the real line. For the sake of conciseness, and because they are ubiquitous in applied research, we focus on linear functionals such as the ATE, probabilities of counterfactual events, and the average derivative, which make (ref) linear programs. In this case the objective functions all take the form $E_{P_W}[f(Y_x,x)]$, where $E_{P_W}[\cdot]$ denotes the expectation with respect to $W$ and where $f(Y_x, x)$ is a function that takes in the position $y$ of the stochastic process $Y_x$ at some value $x$ and provides a real number. Some examples are the following.
Note that in all cases, we can straightforwardly replace the conditioning on values $x_0$ or $x_1$ by more general sets $\mathcal{A}_{x_0}$ and $\mathcal{A}_{x_1}$ without changing anything in the results or proofs. For the sake of notation, we always just consider a point $x_0$ instead of a set.
Let us turn to the constraint of (ref). Since all paths $Y_x(w)$ and $X_z(w)$ are indexed by one realization $w$ of the random variable $W$, $P_W$ must weight these paths in such a way that the joint law $F_{[Y,X]^*(z)}(y,x)$ of the counterfactual stochastic processes $Y_x(w)$ and $X_z(w)$ induced by $P_W(Y_x\leq y, X_z\leq x)$ coincides with the observable distribution $F_{Y,X|Z=z}$.\footnote{Formally, $F_{[Y,X]^*(z)}(y,x)$ is the pushforward measure of $P_W$ via the process $[Y,X]^*_z$, see e.g. \citet*[Definition 3.2]{bauer1996probability} for the definition of a pushforward- or image measure.} The construction of the joint counterfactual process $[Y,X]^*_z$ by the marginal processes $Y_x$ and $X_z$ follows from the structure of model (ref). The idea is that for given paths $Y_x(w)$ and $X_z(w)$ indexed by some $w$, \[[Y,X]^*_z(w) \coloneqq (Y_{X_z(w)}(w), X_z(w)).\] This means that for each value $z$ in the support $\mathcal{Z}$ of $Z$, the path $Y_{X_z(w)}(w)$ is constructed by composing the path $Y_x(w)$ with $X_z(w)$. Specifically, $Y_x(w)$ depends only on the current position in $\mathcal{X}$ of $X_z(w)$ for the given $z$, but no other properties of the respective path $X_z(w)$. This follows directly from the form of model (ref), as $h$ does not depend on $z$, so that intuitively \[Y_{X_z(w)}(w) = h(g(z,w),w).\] The property that $h$ does not depend on $Z$ is called the exclusion restriction. A depiction of the construction of $[Y,X]^*_z$ for two continuous processes $Y_x$ and $X_z$ is provided in Figure (ref).
Figure (ref) also shows how an objective function of the form $E_{P_W}[\mathds{1}_{A_y}(Y_{x_0})]$ is included in this approach. Here, the marginal path $Y_x(w)$ “goes through” the event $A_y$ at the point $x_0$. So for the upper bound (i.e. the maximization of (ref)), the optimal $P_W$ should put a lot of weight on this path. For the lower bound (i.e. the minimization), it should not put a lot of weight on the path.
This representation via $Y_x$ and $X_z$ also allows to conveniently introduce any form of structural assumption into model (ref) by ruling out certain paths. For instance, if the researcher is willing to assume that $h(\cdot,W)$ is continuous in $X$, then she will only require continuous paths $Y_x(w)$. This means that $P_W$ will only put mass on continuous paths and not ones with jumps. Similarly, one can introduce any other form of functional form restriction into the model. Note in general, that all functional form restrictions need to be made with respect to the observable random variables, never with respect to $W$. The unobservable $W$ is simply an index for the respective paths; this is a main difference to other canonical models in econometrics (e.g. \citet*{chernozhukov2005iv}, \citet*{imbens2009identification}), which assume continuity and monotonicity of $h(x,\cdot)$ in $W$.
In this respect, it is also interesting to consider point-identification results from the literature on nonseparable triangular models, as \citet*{imbens2009identification}, \citet*{d2015identification}, or \citet*{torgovitsky2015identification}, which require monotonicity of $h(x,\cdot)$ and $g(z,\cdot)$ in $W$. The idea is that monotonicity makes $h$ and $g$ injective in $v$ and $u$, respectively. In our setting, injectivity of $h$ in $v$ means that the paths of $Y_x(v)$ never intersect, so that for each $(y,x,z)$ there is a unique $w$---this would guarantee that we could point-identify (under some more regularity assumptions) the latent distributions $P_{Y(x)}$ and $P_{X(z)}$. Without monotonicity in the unobservable the joint paths $[Y,X]^*_z$ intersect almost everywhere, so that point-identification is impossible and one has to resort to the partial identification approach using programs (ref). This is another intuitive explanation for why monotonicity is such a staple in the literature on point-identification of nonparametric models.
\paragraph{Intuition in terms of randomized controlled trials with imperfect compliance}\\ The setting of estimating dose-response functions of a treatment through a randomized controlled trial with imperfect compliance is a convenient setting to provide some more intuition of the approach. Let the treatment $X$ be continuous and also assume that the outcome of interest $Y$ is continuous. Then the initially randomly assigned dose $Z$ to each participant is the instrument to $X$. With a continuous treatment $X$, the amount of the treatment actually taken by the participants can deviate from the assigned treatment in an infinite number of ways. The experimenter can only work with the overall distribution $P_{Y,X|Z=z}$ since she knows that the observed distributions $P_{Y|X}$ and $P_{X|Z}$ are for an endogenous $X$ due to the imperfect compliance.
In this setting, the stochastic processes $Y_x$ and $X_z$ have a natural interpretation as the response profiles of hypothetical participants, each being indexed by one element $w\in\mathcal{W}$. A response profile $X_z(w)$ tells the respective hypothetical participant indexed by $w$ which treatment dose $x$ she actually takes for each initially assigned treatment dose $z$. For instance, in analogy to the binary case \citep*{angrist1996identification}, the never taker is the path of the stochastic process which is always zero for each value of $Z$, i.e. this participant's response profile is to take no treatment $x$ for every assigned dose $z$. The path of the complier would be the $45^\circ$ line, since this participant perfectly takes $x=z$ for all $z$.
Analogously a response profile $Y_x(v)$ “tells” a hypothetical participant how much she reacts ($y$) to a certain level of the treatment $x$. In this case the process $Y_x$ models the unobserved heterogeneity in treatment response of each participant. The hypothetical participant who does not respond to the treatment at all has the path $Y_x\equiv 0$ for any value of $x$ for instance. Very intuitively, the idea then is to find the relative optimal composition of hypothetical participants (i.e. the weights $P_W$ on the respective paths of the stochastic processes) that maximizes (for an upper bound) or minimizes (for a lower bound) the relative number of hypothetical participants which on average have the strongest or weakest response to a change of the dose from $x_0$ to $x_1$, subject to the constraint that the composition $P_W$ of the hypothetical participants replicates the joint observable weight $P_{Y,X|Z=z}$. This framework is the natural generalization of the distinction into always taker, never taker, complier, and defier from \citet*{angrist1996identification} to the continuous setting, where there is a continuum of response profiles, not just four as in the binary case.
A way to introduce assumptions in this setting is via $P_W$. In particular, one could assume that all hypothetical participants have continuous responses, i.e. if participant's response profile indexed by $w$ is to take $x_0$ units of the treatment when assigned $z_0$ units, then they will also take $x_0+\delta$ units when assigned $z_0+\delta$ units for some small $\delta>0$. An assumption like this would rule out drop-outs from the trial for instance. Response profiles that allow for drop outs are important in trials where the treatment can have severe side-effects under too high doses, like chemotherapy treatments for instance. One can model this by allowing for processes that jump to zero after a given $z_0$ and stay there for all $z>z_0$. This example should make clear that the stochastic process representation we introduce in this article is convenient for modeling latent behavior in the unobserved heterogeneity.
Let us now turn to examples of applications of model (ref) and the corresponding programs (ref) which will also highlight the curse of cardinality when trying to solve (ref) in practice. We start with the binary analogue to our treatment setting.
\paragraph{Example 1: Curse of cardinality in treatment estimation}\\ Consider the treatment setting as above, but with a binary treatment $X\in\{0,1\}$, binary outcome $Y\in\{0,1\}$, and binary intent-to-treat $Z\in\{0,1\}$. This is the classical setting, see for instance \citet*{angrist1996identification}, \citet*{balke1997bounds}, \citet*{angrist1996identification}, \citet*{manski2003partial}, \citet*{russell2019sharp}.
Suppose we are interested in estimating the ATE of a change from $X=0$ to $X=1$. In this simple case, (ref) reduces to a linear program with $9$ equality constraints and $16$ non-negativity constraints \[P_{Y,X|Z=z}(y,x) = P_W(Y_{X_z}=y, X_z=x)\] over distributions $P_W$ on the support $\mathcal{W}$ of $W$, which in this case consists of $(2^2)^2=16$ elements (\citet*{balke1994counterfactual} and \citet*{balke1997bounds}).
The reason for this is as follows. In both the first- and second stage there are four possible functions mapping $Z$ to $X$ and $X$ to $Y$, respectively. In the first stage, these paths of processes have been given specific names \citep*{angrist1996identification}: the never takers, i.e. the path $X_z=0$ for all $z\in\{0,1\}$; the compliers, i.e. the path $X_0=0$ and $X_1 = 1$; the defiers, i.e. the path $X_0=1$ and $X_1=0$; the always takers, i.e. the path $X_z =1$ for all $z\in\{0,1\}$. The analogous set-up holds for the second stage with the same four types. $W$ therefore indexes every combination between one of the types in the first-and second stage, which is $16$. The number of equality constraints follows from the fact that the observable distribution $P_{Y,X|Z=z}$ can take $2^3=8$ values: each combination of $Y, X$, and $Z$. One more equality constraint is needed to guarantee that all probabilities sum to $1$, and $16$ inequality constraints are required to require that all proabilities are non-negative.
Note that one can arrive at the same program via many different routes: random set theory using Artstein's inequality (\citet*{molchanov2005theory}, \citet*{beresteanu2012partial}, \citet*{chesher2017generalized}, \citet*{russell2019sharp}) or optimal transport theory \citep*{galichon2011set}. \citet*{russell2019sharp} gives a clear account of this setting.
As stated, this problem is a simple linear program and is trivial to compute. However, in settings where $Y$, $X$, and $Z$ are allowed to have more than two points in their support, the complexity of the linear program grows at least exponentially. Consider the case where $Y$ can take $q$ different values, $X$ can take $r$ different values, and $Z$ can take $s$ different values. In this case, the linear program will have $q\cdot r\cdot s +1$ equality constraints. Even more drastically, the number of elements in the support of $\mathcal{W}$ is $r^s\cdot q^r$ with the same number of inequality constraints.
This implies that already in a setting where $Y$, $X$, and $Z$ are allowed to take three variables each, $\mathcal{W}$ would consist of $27^2=729$ points, while the number of equality constraints would be $28$, with another $729$ nonnegativity constraints. Already setting up the linear programming problem, either through Arstein's inequality, optimal transport, or the stochastic process method depicted here, will be complicated to do. In the existing literature, \citet*{cheng2006bounds} circumvent this problem in the setting of three-armed randomized controlled trials with imperfect compliance by making monotonicity assumptions which rule out most of the paths, reducing the linear program to a simpler form.
However, already in the case of where $Y$, $X$, and $Z$ are allowed to take $5$ elements each is it practically impossible to set up the linear program over $(5^5)^2=9,765,625$ elements in $\mathcal{W}$ with the same number of non-negativity constraints and $5^3+1=126$ equality constraints. This means that any existing approach requiring to set up an optimization problem of this form, such as \citet*{balke1994counterfactual}, \citet*{chesher2017generalized}, \citet*{beresteanu2012partial}, \citet*{galichon2011set}, russell2019sharp and others cannot be applied in these practical settings, even though they show how to set-up these linear programs in theory. This constitutes what we call the “curse of cardinality”, and is the main reason for why current practical application of this idea as in \citet*{demuynck2015bounding}, \citet*{laffers2015bounding}, \citet*{mogstad2018using}, or \citet*{chesher2017generalized} focus on the case where $X$ is binary.
The example with randomized controlled trials and a continuous treatment we introduced above is the extreme case of this, where $Y$, $X$, and $Z$ are continuous. In this case, there currently does not exist an approach to solve the programs (ref). In the next section, we introduce a sampling approach that approximates the solutions to the problems (ref) while giving probabilistic approximation guarantees, even if all $Y$, $X$, and $Z$ are continuous. Before we do this, let us give other examples for potential applications.
\paragraph{Example 2: Time-varying treatments in program evaluation}\\ In general, $Z$ in model (ref) is a classical instrument. However, the model and the approach presented in this paper are even more broadly applicable. For instance, in models for the evaluation of a training program \citep*{ashenfelter1978estimating} one can consider $Z$ to be time. In this case, solving the programs (ref) still provides bounds on the causal effect of interest, because the model assumes that the counterfactual distribution $P_{Y(x)}$ does not change over time, i.e. that the actual causal effect of the training $x$ on the outcome of interest is fixed. This follows from the fact that $Y_x\coloneqq h(x,W)$ does not explicitly depend on $t$.
To make this more tangible, consider the hypothetical example of a rehabilitation program (such as a self-help group, a health program, etc.) that runs over several weeks. The treatment $X$ is the time (e.g. days) spent in the program, which makes it a continuous treatment. $X$ being continuous introduces the exact curse of cardinality as described in Example 1. $Y$ could be a binary variable $Y\in \{0,1\}$, indicating if the individual has not relapsed $6$ months after the full program ended. Since participation in the program is voluntary, participants can drop out any time; this introduces an endogeneity bias as participants will choose to opt in or out depending on their expected outcome. One instance of this endogeneity problem is the classical Ashenfelter dip \citep*{ashenfelter1984using}, which shows that the outcome of interest changes to the negative right before they opt to partake in a program.
In this setting, model (ref) implies that the second stage process $Y_x\coloneqq h(x,W)$ is independent of time, while the participation in the treatment $X$ can vary with time. Time $Z=t$ is never influenced by the model; furthermore, the unobservable heterogeneity $W$ is fixed in time (it indexes the response profiles of the hypothetical participants), so that the independence restriction $Z\protect\mathpalette{\protect\independenT}{\perp} W$ is satisfied. Using the linear programming approach (ref), one can make assumptions on the participants in the first stage $X_t=g(t,W)$, by modeling their response profiles including drop-in and drop-out rates. Since $X$ is the time spent in the program, all paths $X_t(w)$ are by definition restricted to the $45$-degree line on a plot between $X$ and $Z$, the only difference being the drop-in, or drop-out rate. To adequately model these drop-ins and drop-outs, the stochastic processes $X_t$ in this case should allow for jumps. Suppose for the sake of exposition that the program allows participants to drop in or drop out on a continuous basis. In this case, we still have a continuum of paths $X_t(w)$, i.e. a continuum of elements $w$. Two examples of paths are depicted in Figure (ref). This continuum of paths is the extreme case of the curse of cardinality, just as in Example 1; in particular, no current approach can deal with a model like this without making stronger assumptions.
This example also shows how one can introduce assumptions into the model via restricting the paths. For instance, it could be the case that the program does not allow people to drop-in, and requires everyone to start at the same time. In this case paths of hypothetical participants that model drop-ins like the dashed path in Figure (ref) should not be part of the model. Another assumption could be that for certain time-periods no drop-ins or drop-outs are allowed, in which case paths with drop-ins or drop-outs during this period would not be part of the model. When solving the programs (ref), such an assumption could be included by requiring the optimal $P_W$ to always put weight of zero on these kinds of paths. It is here where the stochastic process framework is so convenient, as the researcher can model individual behavior by restricting the unobserved heterogeneity in very specific ways. Introducing these assumptions is significantly harder to do in the standard representations of the instrumental variable mentioned above.
In order to model the second stage $Y_x$, one could assume stochastic processes $Y_x$ which are zero everywhere but jump to $1$ after a certain length $[x_0,x_1]$, where $x_0$ is the left-most value of the overall interval. This assumption encapsulates the unobserved heterogeneity as different responses to treatment length: participants might need different time lengths to obtain the required outcome. Figure (ref) depicts two potential paths. Together those marginal paths introduce joint paths as depicted in Figure (ref).
For this model, an objective of interest might be the probability that the participant does not relapse $6$ months after treatment if the participant has continuously stayed in the program for a period $[x_0,x_1]$. The objective function in (ref) then takes the form \[E_{P_W}[f(Y_x,x)]\coloneqq E_{P_W}\left[Y_{x\in [x_0,x_1]}\right],\] where $Y_{x\in [x_0,x_1]}$ means the part of the process $Y_x$ when $x\in[x_0,x_1]$. As we show below, the proposed method below can handle discontinuous objective functions if the discontinuity is of the form $f(P_W;x_0)=E_{P_W}[\mathds{1}_{A_y}(Y_{x_0})]$, where $A_y$ is a Borel subset of the support $\mathcal{Y}$ of $Y$. The method introduced in this article can deal with complex models like this, where the treatment is considered to be more intense over time. This is a case of the classical “dose-response” estimation. Here, the dose is the number of days in the program.
\paragraph{Example 3: Partial identification of average welfare}\\ Unobserved individual heterogeneity is a major factor contributing to the variability in empirical demand estimation, which is why models with general unobserved heterogeneity are particularly important in this setting, as argued in \citep*{hausman2016individual}. In that article, the authors consider the following model for the demand function \[q(x,\eta) = \operatorname*{\arg\!\max}_{q\geq0,a\geq0} U(q,a,\eta)\qquad\text{s.t}\qquad p'q+a\leq y.\] The demand function $q(x,\eta)$ is nonseparable in the unrestricted unobserved heterogeneity $\eta$ and $x\coloneqq(p',y)'$ consisting of the price vector $p$ relative to a numeraire good $a$, and the income level $y$ of the individual relative to the same numeraire. $U(q,a,\eta)$ is a (concave in $\eta$) utility function which the authors assume to be increasing in $q$ and $a$. $p$ and $y$, and therefore $x$, are usually continuous variables, which again introduces the extreme curse of cardinality from Example 1 into the problem.
The demand function corresponds to the second stage of model (ref). Concavity of $U$ guarantees that $q(x,\eta)$ is an actual function in $\eta$ and not a correspondence, so that each path of the process $q(x,\eta)$ is a function. This is the same setting we consider in this article.\footnote{Our approach can straightforwardly be adapted to allow for correspondence-valued processes without any changes in the overall approach.} Allowing for general $\eta$ precludes point-identification without stronger assumptions in general, which is why \citet*{hausman2016individual} provide an informal approach to obtain bounds on the average consumer surplus.
The approach proposed in this article can be applied to this setting and hence provides the formal justification for the approach in \citet*{hausman2016individual}. Furthermore, we also provide the asymptotic properties of the upper- and lower bound, which can be used to perform inference in their setting. \citet*{hausman2016individual} also assume that $x\protect\mathpalette{\protect\independenT}{\perp} \eta$ in their application, a strong assumption in practice. Our approach is designed for handling more general models allowing for $X\not\protect\mathpalette{\protect\independenT}{\perp}\eta$ and instrumenting for $X$.
The stochastic process approach formalized in this article has another important benefit in this setting. One can introduce any form of assumption into the model (the demand function $q(x,\eta)$) by ruling out paths $q(x,\eta)$ which do not conform with this assumption. This is the same idea as in the previous example of response profiles. In particular, each path is a response profile of an actor $\eta$ in the market.
Note that in combination with the idea of \citet*{hausman2016individual}, our approach is a complementary approach to the finitary method proposed in \citet*{kitamura2018nonparametric} for partial identification of average welfare. \citet*{kitamura2018nonparametric} work directly with individual budget sets and define a finite linear programming problem which grows with the number of data points in practice. Our method is an attractive alternative to this in two main ways. First, we can allow for continuous functions and variables. In this setting the programs (ref) become infinite dimensional optimization programs. The practical advantage of this is that we can put a grid on the data and need not evaluate the linear program at each data point, which can be impossible to do when the data-set consists of many oservations. Second, and more importantly, the continuous setting allows for the introduction of assumptions on the demand functions $q(x,\eta)$ like smoothness or differentiability. These assumptions are vacuous in a discrete setting such as \citet*{kitamura2018nonparametric}. Such regularity assumptions follow naturally from the economic models (e.g. the properties of the Slutsky matrix) and should be included as assumptions in the empirical setting, which our continuous setting allows for.
Before providing the solution approach to solving the programs (ref) in practice, we define the stochastic process representation of model (ref) formally. From now on, we let $Y$, $X$, and $Z$ be random variables whose corresponding laws $P_Y$, $P_X$, and $P_Z$ have compact supports $\mathcal{Y}$, $\mathcal{X}$, $\mathcal{Z}\subset\mathbb{R}^d$ for some $d\geq1$. In the following, we set $d=1$ for notational convenience, but all results in higher-dimensional settings.\footnote{One can allow for Polish spaces, i.e. complete, separable and metrizable topological spaces, under some formal complications. The compactness restriction can be dropped at the cost of more mathematical formalism.} Moreover, since we focus on compact supports on $\mathbb{R}$ it is also without loss of generality to assume $\mathcal{Y}=\mathcal{X}=\mathcal{Z}=[0,1]$, which we do for the sake of notation.
As stated previously, the main idea for the approach is to interpret the counterfactual laws $P_{Y(x)}$ and $P_{X(z)}$ of the instrumental variable model as laws of counterfactual stochastic processes $Y_x$ and $X_z$, where the randomness is induced by the unobservable $W$. Each element $w\in\mathcal{W}$ indexes one path $Y_x(w)$ of the process. This representation allows to optimize over the unobserved heterogeneity $W$ directly instead of the laws $P_{Y(x)}$ and $P_{X(z)}$ induced by it, which is the key for making the approach feasible. The first question to answer is in which space these paths $Y_x(w)$ and $X_z(w)$ live. Completely unrestricted functions $h$ and $g$ as allowed for in model (ref) correspond to stochastic processes defined in $\mathbb{R}^{\mathbb{R}}$, i.e. the space of all functions from $\mathbb{R}$ to $\mathbb{R}$. This space is large enough to allow for response profiles such as in Example 2, where candidates are in principle allowed to drop in and drop out uncountably many times.
From a mathematical perspective, $\mathbb{R}^{\mathbb{R}}$ is too big and the Borel $\sigma$-algebra is too coarse for many assumptions of interest, however. For example, many subspaces like $C(\mathbb{R})$, the space of continuous functions on $\mathbb{R}$, are not even measurable in the Borel $\sigma$-algebra on $\mathbb{R}^{\mathbb{R}}$ \citep*[Corollary 38.5]{bauer1996probability}. Also, the response profiles allowed by $\mathbb{R}^{\mathbb{R}}$ are too general for most practical purposes. For instance, allowing for uncountably many drop-ins and drop-outs as one response profile in Example 2 is of no practical interest, as one always observes finitely many drop-ins and drop-outs of a person. This level of generality is hence really not needed for modeling the hypothetical participants. On the other hand, restrictions like continuity which require all response profiles to be continuous, might be too restrictive in many settings. An example is again Example 2, where allowing for jumps is required by the problem.
It is here where the stochastic process representation is useful for the first time, as it allows to introduce very weak but easy to understand assumptions which restrict the unobservable heterogeneity $W$. A rather weak but restrictive enough assumption for all practical settings is to only consider stochastic processes that possess almost surely c\`adl\`ag paths, i.e. paths that are continuous on the right with limits on the left. These paths allow to model all possible practically relevant response profiles, as they allow for finitely many jumps larger than any value $\varepsilon>0$ in the paths. This means that one can realistically model the response profiles of Example 2; moreover, paths that require continuity of the responses are measurable in this space, too.
The standard space for these processes is the Skorokhod space $D(\mathcal{X})$, which we equip with the Skorokhod metric \[d_S(Y_x,Y'_x)\coloneqq \inf_{\lambda\in\Lambda}\left\{\max\{\|\lambda-I\|_\infty, \|Y_x - \lambda Y'_x\|_\infty\}\right\}\] to make it a separable and metrizable space. In particular, a metric similar to the Skorokhod metric induces the same topology on $D(\mathcal{X})$ and makes it a Polish space, so that from now on we consider $D(\mathcal{X})$ with the Skorokhod metric a Polish space.\footnote{The reason for assuming compactness of $\mathcal{Y}$, $\mathcal{X}$, and $\mathcal{Z}$ is to be able to work within the Skorokhod space over a compact set, which is nicer to handle. Furthermore, one can define the Skorokhod space also when $\mathcal{X}$ is higher-dimensional, see \citet*{neuhaus1971weak}, \citet*{straf1972weak}, and \citet*{bloznelis1994central}. The metric which is topologically equivalent to the Skorokhod metric is often called the Billingsley metric \citep*{billingsley1999convergence}.} Here, $\Lambda$ is the set of all strictly increasing and continuous functions $\lambda:\mathcal{X}\to\mathcal{X}$, $I:\mathcal{X}\to\mathcal{X}$ is the identity function on $\mathcal{X}$, and $\|\cdot\|_\infty$ is the supremum norm \citep*[Section 12]{billingsley1999convergence}. The Skorokhod space is perfect as the underlying space for the representation of model (ref) via stochastic processes, as it allows for all practically relevant response profiles while being a Polish space, which provides good measurability properties.
One sufficient and necessary condition for a function $f:\mathcal{X}\to\mathbb{R}$ to lie in the Skorokhod space is that its extended modulus of continuity $\omega'_f(\delta)$ satisfies $\omega'_f(\delta)\to0$ for any $\delta\downarrow 0$ \citep*[p. 123]{billingsley1999convergence}.\footnote{$\omega'_{f}(\delta)\coloneqq \inf_{\{x_i\}}\max_{i\in\mathbb{N}}\sup_{x,x'\in[x_{i-1},x_i)}|f(x)-f(x')|$, where the infimum extends over all $\delta$-partitions $\{x_i\}_{i\in\mathbb{N}}$, i.e. partitions where the values $x_i$ are at a distance of at least $\delta>0$ from each other, see \citet*[p. 122]{billingsley1999convergence}. This definition extends the modulus of continuity for continuous functions.} We therefore require the following Assumption throughout.
Assumption (ref) only allows paths $Y_x(w)$ and $X_z(w)$ that have finitely many discontinuities which exceed any given real number. This assumption is weak enough in all practical settings, as argued above. Under Assumption (ref), one can formally introduce the stochastic process representation of model (ref). Throughout the article, these stochastic processes are defined for “deterministic time”, even though $X$ and $Z$ are random variables, as the supports $\mathcal{X}$ and $\mathcal{Z}$ are fixed. This is without loss of generality in this model; in fact, all of the results in this article go through almost verbatim under the assumption that $X$ and $Z$ have random supports $\mathcal{X}$ and $\mathcal{Z}$. The only difference will be the formal complication of working with stochastic processes in random time.\footnote{In this case, one would need to replace Kolmogorov's extension theorem by the random version provided in Theorem 1 of \citet*{hu1988generalization}. } It is, moreover, not even clear that the random support assumption generates more general models than the deterministic support assumption.
The following proposition provides the formal result for the proposed stochastic process representation of model (ref).
Proposition (ref) states that the counterfactual probability measures $P_{Y(x)}$ and $P_{X(z)}$ are the laws of $Y_x$ and $X_z$, where the randomness stems from the unobservable confounder $W$. Since $W$ is unobserved, one can define it on $([0,1],\mathscr{B}_{[0,1]},P_W)$, i.e. to identify the support $\mathcal{W}$ of $W$ with the unit interval. This follows from standard isomorphism results, see for instance \citet*[Theorem 9.2.2]{bogachev2007measure2}. It is in fact possible to construct these explicit isomorphisms between general Borel measures on Polish spaces and Lebesgue measure on the unit interval based on the approach laid out in \citet*{kuratowski1934generalisation}. One standard application of these isomorphisms is the construction of Wiener measure, see for instance \citet*{hess1982kuratowski}.
This is the main difference to the other existing identification results in the literature: we do not make topological or algebraic assumptions on the support $\mathcal{W}$. $W$ is always infinite dimensional in our setting and completely abstract: only the cardinality of $\mathcal{W}$ is important, which needs to be such that each path $Y_x(w)$ and $X_z(w)$ in the Skorokhod space---a Polish space under the Skorokhod metric---can be indexed by one $w$. By the above isomorphism results, this is already possible if $\mathcal{W}\subset[0,1]$. This is the main reason why we can define $\mathcal{W}$ to be the unit interval while the existing identification approaches in econometrics make structural assumptions on $\mathcal{W}$.
The construction of $W$ on the unit interval allows us to define $\mathscr{P}^*(\mathcal{W})$ as the set of all measures $P_W$ on $([0,1],\mathscr{B}_{[0,1]})$ which induce paths on the Skorokhod space defined in Assumption (ref). The idea is that by a construction like the one in \citet*{kuratowski1934generalisation}, a measure $P_W$ will induce a measure on path space, so that we can optimize over measures on the unit interval in the programs (ref). In this respect, all additional functional form assumptions on the model will be made by shrinking $\mathscr{P}^*(\mathcal{W})$. For instance, making continuity assumptions between $Y$ and $X$ or $X$ and $Z$ in (ref) translates to continuity assumptions on the paths $Y_x$ or $X_z$. In this case, $\mathscr{P}^*(\mathcal{W})$ contains only measures $P_W$ on $\mathcal{W}$ which put positive probability on continuous paths. Since under Proposition (ref) all processes are progressively measurable with respect to their natural filtration, one can even introduce dynamic assumptions this way, like mixing properties. This can be interesting in settings similar to Example 2, where the treatment is time dependent.
This section introduces the probabilistic approach to solve the programs (ref) approximately in practice with probabilistic finite sample guarantees for the validity of the approximation. The method proceeds by sampling bases of the path space and solving semi-infinite dimensional analogues \citep*{anderson1987linear} of (ref) over this sample of basis functions. The randomness introduced by sampling the paths permits the use of concentration results \citep*[Chapter 2.14]{wellner2013weak} which provide probabilistic guarantees for the approximation of (ref) by the random sample.
The idea is to approximate the infinite dimensional problems (ref) by their semi-infinite analogues
for some small $\varepsilon>0$ and all $z\in[0,1]$, and where the $L^2$-norm is taken with respect to $y,x$.\footnote{Rewriting the programs (ref) in their penalized form later on will reveal that $\varepsilon$ fulfills the same purpose as a penalty term for the constraint. In this sense, $\varepsilon$ is a penalty parameter of the program which needs to be chosen appropriately. See the next section for a discussion. Overall, the relaxation of the linear constraints to an $L^2$-constraint is made exclusively based on finite sample considerations. In the population $\varepsilon=0$.} The indicator functions $\mathds{1}_{[0,\cdot]\times[0,\cdot]}\left(\tilde{Y}^\kappa_{\tilde{X}^\kappa_{z}(i)}(i),\tilde{X}^\kappa_{z}(i)\right)$ capture the requirement that $\left(\tilde{Y}^\kappa_{\tilde{X}^\kappa_z(i)}(i),\tilde{X}_z^\kappa(i)\right)\leq(y,x)$ at $z$, where the inequality is considered element-wise. An equivalent way to writing this is $\mathds{1}_{[\tilde{Y}_{\tilde{X}_z^\kappa(i)}^\kappa(i),1]\times [\tilde{X}^\kappa_z(i),1]}(y,x)$.
Underlying the approximations (ref) is a random sample of size $l$ of a number of $\kappa$ basis functions which approximate the path space of the processes $Y_x$ and $X_z$. These approximations are denoted by
for coefficients $\beta,\alpha$ and basis functions $b(x),a(z)$. Of particular convenience are basis functions used in Sieves estimation such as (trigonometric) polynomials, splines, wavelets, etc. \citep*{chen2007large}.
The dependence of $\tilde{Y}_{x}(i)$ and $\tilde{X}_{z}(i)$ on the index $i$ shows that the problems (ref) are indeed semi-infinite dimensional in the sense of \citet*{anderson1987linear}. In particular, the index $i$ now runs over finitely many elements $l$ and replaces the variable $w\in\mathcal{W}$ on the state space of $P_W$.
Let us now go over the most important terms in (ref), $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$, $\tilde{Y}^\kappa_{x}(i)$ and the objective function $f(Y_x,x_0)$, one by one.
\paragraph{1. The representative law}\\ The term $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$ is of fundamental importance. In particular, the empirical sampling law $\hat{P}_0\in\hat{\mathscr{P}}^*(\mathcal{W})$ is one representative law on the paths of stochastic processes, which is used for sampling the basis functions. The optimization is then over all $\hat{P}_W$ which are absolutely continuous with respect to $\hat{P}_0$, so that $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$ is the Radon-Nikodym derivative. This construction arises naturally, as the empirical sampling law $\hat{P}_0$ determines the sample of the $l$ paths, over which an optimal law $\hat{P}_W$ will be chosen to solve the programs (ref).
$\hat{P}_W$ must by construction be absolutely continuous with respect to $\hat{P}_0$ as it can only place positive mass on the $i$ sampled paths, which have been determined via $\hat{P}_0$. In other words, $\hat{\mathscr{P}}^*(\mathcal{W})$ is the set of all probability measures which do not put positive measure on paths other than the $l$ paths sampled via $\hat{P}_0$. This implies a natural assumption on the data-generating process in the population.
Assumption (ref) is the theoretical analogue to the fact that all $\hat{P}_W$ are absolutely continuous with respect to $\hat{P}_0$ by construction. Since by the finite dimensional construction all measures $\hat{P}_W$ are automatically absolutely continuous with respect to $\hat{P}_0$, theoretical measures $P_W$ which are not absolutely continuous to $P_0$ can never be detected. This implies that Assumption (ref) is non-testable.
In practice, it is through $P_0$ that the researcher introduces functional form restrictions into the model. For instance, if one wants to assume continuity in the relation between $Y$ and $X$, one will choose a $P_0$ which only puts positive measure on continuous paths. $\hat{P}_0$ will then only sample continuous paths in practice. This way, one can theoretically introduce any form of functional form restriction into the model. In Example 2 for instance, one would define $P_0$ to be a measure that only puts positive probability on paths $X_t$ on the $45$-degree line with finitely many jumps to zero.
Theorem (ref), the main approximation result, requires some more regularity of the representer. In particular, it requires that $P_W$ does not have atoms. This implies that any constraint in (ref) has zero influence on the overall constraint. Since (ref) allow for a continuum of constraints, they still restrict the problem jointly. Below, right after the statement of Theorem (ref) we argue intuitively what this requirement of absolutely continuous $P_W$ means in other settings.
\paragraph{2. The wavelet basis}\\This articles works with a shape-preserving wavelet basis, which is defined as \[\varphi_{\kappa j}(x) \coloneqq 2^{\frac{\kappa}{2}}\varphi(2^{\kappa}x-j)\qquad \kappa,j\in\mathbb{Z}\] with $\varphi:\mathbb{R}\to[0,1]$ of the form \[\varphi(x) \coloneqq
.\] Based on this the notation for the paths sampled via this wavelet basis for dilation $\kappa$ is
where the sums in the definition are both finite since the unit interval is bounded. This wavelet basis preserves shapes in the sense that an approximation of a monotone (convex) function via this basis will itself be monotone (convex), see \citet*{anastassiou1992monotone} and \citet*{anastassiou1992convex}; this is an important feature when introducing shape assumptions into the model in practice.
Other approaches work equally well. For instance, if one is willing to assume that all paths are smooth, one could simulate paths from a standard diffusion process via the Karhunen-Lo\`eve transform. In this case, one would have to extend Proposition (ref) to the case where $Y$, $X$, and $Z$ are not compact. This can be done based on results in \citet*[section 16]{billingsley1999convergence}.
\paragraph{3. Regularity assumptions on the objective function}\\ One also needs to impose some regularity on the objective function.
The assumption allows for general objective functions, including the ones from Section (ref). Lipschitz continuity is a very weak requirement. Also, allowing for indicator functions is important in order to allow for counterfactual probabilities as argued above.
Before stating the main approximation result, we rewrite (ref) and (ref) as penalized programs. The reason is that in this form, we can use the standard theory of $M$-estimators \citep*[chapter 3.4]{wellner2013weak} to obtain the approximation results.
Lemma (ref) is important, as the alternative programs (ref) and (ref) are written in such a way that we can apply the standard theory of $M$-estimators \citep*[chapter 3.2]{wellner2013weak} to the Radon-Nikodym derivatives $\frac{dP_W}{dP_0}$ in order to obtain the concentration bounds. Note how the $L^2$ norm and the integral with respect to $P_0$ have switched places. We view the programs as $M$-estimators with respect to the random sample $l$ of paths. All classical concentration results are applicable in this setting, the only difference is that we now apply them with respect to the artificially introduced randomness by sampling paths. From now on, we will focus exclusively on the programs (ref) and (ref) when deriving the properties of our approach.
The following theorem is the main result of the article and gives probabilistic guarantees of for the approximation of the infinite programs (ref) by the semi-infinite programs (ref). For this bound, we need to require at Lipschitz-continuity of the objective function, which is not satisfied under part (ii) of Assumption (ref). In this case, we need to approximate the indicator function $\mathds{1}_{[0,y]}(Y_{x_0}(w))$ by a logistic function $\mathcal{S}(Y_{x_0}(w)),y,\eta)$ on $[0,1]$, which takes the form \[\mathcal{S}(Y_{x_0}(w),y,\eta)\coloneqq\frac{1}{1+e^{-\eta(y-Y_x(w))}}.\]
Theorem (ref) jointly provides a lower bound on the number of sampled paths $l$ as well as the number of terms in the wavelet decomposition $\kappa$ for the semi-infinite program to provide a good approximation to the infinite dimensional program for a requested probability $\rho$. In fact, it is a finitary analogue of a classical consistency result, as it also provides the rate of convergence as $l,\kappa\to\infty$. Theorem (ref) holds for a fixed observable distribution $F_{Y,X|Z}$; this observable distribution can be an estimator $\hat{F}_{Y,X|Z;n}$, in which case we consider the number of data points $n$ as fixed. The result follows from concentration results for empirical processes with functions of finite entropy numbers \citep*[chapter 2.14]{wellner2013weak}. The bound on the right hand side is composed of two terms. The first term in brackets captures the nonprobabilistic approximation of the value functions when the paths $Y_x$ and $X_z$ are approximated by the wavelet bases $\tilde{Y}_x^\kappa$ and $\tilde{X}_z^\kappa$, respectively and depends on the extended modulus of continuity $\omega'$ of the paths of either process. The two terms here correspond to the constraint and the objective function, respectively. The second term is the probabilistic concentration bound based on the sample of $l$ paths. The constants here also depend on $\kappa$. In the case where the objective function is an indicator function, we need an additional approximation term which follows from approximating the indicator function by a logistic function.
\paragraph{Intuition for Theorem (ref) and comparison to other approaches}\\ Theorem (ref) is but one example of many different approximation guarantees of similar form, depending on what kind of assumptions one is willing to make on the paths (via $P_0$) and $\frac{dP_W}{dP_0}$. In fact, it is the introduction of the randomness via sampling that permits the use of statistical procedures for the purpose of function approximation. The idea is to perform statistical estimation with respect to the counterfactual elements via sampling. In this respect, Assumption (ref) is important for obtaining the rate of convergence in Theorem (ref): if all Radon-Nikodym derivatives $\frac{dP_W}{dP_0}(w)$ are smooth and have support in all of $[0,1]$, one does not need many draws $l$ to approximate the optimization problem. For more concentrated densities, the approximation can be much worse in the sense that one needs a substantially larger sample $l$ to approximate the optimal density. The approach proposed via Theorem (ref) can hence be seen as a nonparametric estimation of a probability density with infinite dimensional support (i.e. where each data-point is a path of a process) and where one has complete control over the data sample. Different assumptions on the smoothness of the densities will lead to different lower bounds on the number of sampled paths.
Theorem (ref) is a quantitative approximation result which is significantly harder, if not impossible, to achieve with an approach that obtains a complete set of inequalities describing the identified set (e.g. \citeauthor*{chesher2017generalized} chesher2017generalized). In fact, Assumption (ref) has potential analogues in the inequality approach: smoothness of the Radon-Nikodym derivative is similar to the assumption that the identified set can be described by (uncountably) many inequalities which all only contribute minimally to the identified set. If one is not willing to make this assumption, it can happen that the Radon-Nikodym derivatives are much more concentrated, i.e. putting significantly more weight on only a few paths. In this case, the bound in Theorem (ref) implies that one needs to sample more paths depending on the new assumption on $P_0$ and $\frac{dP_W}{dP_0}$. The latter setting corresponds to the case where there are only a few inequalities which “determine the form of the identified set”. In this case, one would also need to sample many more inequalities in order to “stumble upon” the few which determine the shape of the set (see the argument in \citeauthor*{de2004constraint} de2004constraint).
Theorem (ref) hence complements the existing results on the probabilistic approximation of infinite dimensional (linear) programs by semi-infinite dimensional programs proposed in \citet*{girosi1995approximation} and \citet*{de2004constraint}. In particular, the method allows for standard statistical methods to obtain approximation results by introducing the sampling probability $P_0$, which leads to flexible and general quantitative results like Theorem (ref). In contrast, \citet*{de2004constraint} use a learning result (Theorem 8.1.4 in anthony1997computational anthony1997computational) to obtain a probabilistic approximation result, which is not as flexible as using nonparametric estimation results. Intuitively, Theorem (ref) puts a probability measure on the “set of all inequalities” in a classical linear program. Instead of working with inequalities like \citet*{de2004constraint} or \citet*{chesher2017generalized}, we work with stochastic processes and put a probability measure on the paths.
The idea of the sampling approach is to introduce additional randomness into the problem by sampling $l$ paths, i.e. sampling the coefficients of the basis functions. The classical statistical randomness in problems (ref) and (ref) follows from approximating the population distribution $F_{Y,X|Z}$ by a finite-sample estimator $\hat{F}_{Y,X|Z;n}$, potentially smoothed via some bandwidth $h_n$, where $n$ denotes the size of this sample.\footnote{The theoretical results in this section are derived for the standard empirical cumulative distribution function $\hat{F}_{Y,X|Z;n}$. They extend straightforwardly to smoothed estimators $\hat{F}_{Y,X|Z;h_n}$. For this, all one has to do in the proofs is to replace the classical Glivenko-Cantelli and Donsker theorems by analogous versions for smoothed empirical processes. These results (and the corresponding weak assumptions) are contained in \citet*{gine2008uniform} for instance.}
This subsection introduces large sample results which enable the researcher to perform inference on the solution of the programs (ref). These results are only derived for each bound separately, i.e. for the lower bound and the upper bound. In order to derive inference bounds on the whole identified set, it might be possible to use well-established results from the literature, such as \citet*{imbens2004confidence}, \citet*{stoye2009more}, and especially \citet*{kaido2019confidence}.
Even though the theoretical programs (ref) have relaxed constraints, it could potentially still be the case that for very large $\lambda$ there exist data-generating processes $F_{Y,X|Z}$ that no $P_{W}\in\mathscr{P}^*(\mathcal{W})$ can replicate; this is especially true in the case where many additional form restrictions are imposed on the model via $P_0$. In this case, however, the data-generating process directly introduces testable assumptions on the model, as the model under the respective assumptions is not able to replicate the observed data. Since the focus of this article is on estimation of bounds, it is convenient to introduce an assumption on the data-generating process which guarantees that the constraint is non-empty.
In the following, $\mathcal{F}_{Y,X|Z=z}$ denotes a set of all conditional cumulative distribution functions on $[0,1]^2$ satisfying certain assumptions the researcher is comfortable to assume for the given data-generating process. This set is equipped with the $L^\infty([0,1]^2)$-norm.
Assumption (ref) is deliberately high-level, because (i) specific assumptions on the data-generating process $\mathcal{F}_{Y,X|Z=z}$ and the model $\mathscr{P}^*(\mathcal{W})$ usually come from economic theory, (ii) an empty constraint for a given data-generating process corresponds to the existence of testable implications on the model, and (iii) Assumption (ref) is only required to obtain regular and well-behaved asymptotic results, but not for any other results in this article. It is also straightforward to derive a low-level sufficient condition on $\mathcal{F}_{Y,X|Z=z}$ implying Assumption (ref) in the case where only Assumption (ref) but no other shape restrictions are imposed on $\mathscr{P}^*(\mathcal{W})$.
For instance, when the set $\mathcal{F}_{Y,X|Z=z}$ consists only of distribution functions $F_{Y,X|Z=z}$ which are laws to stochastic processes $[Y,X]_z$ whose paths have an extended modulus of continuity $\omega'_{[Y,X]_z}(\delta)$ satisfying $\limsup_{\delta\to0} \omega'_{[Y,X]_z}(\delta) = 0$, then Assumption (ref) is fulfilled. This follows directly from the fact that this condition on the extended modulus of continuity implies that almost all paths of the observable process $[Y,X]_z$ lie in the Skorokhod space defined by Assumption (ref), which is the assumption made on the latent process $[Y,X]^*_z$.
$\hat{F}_{Y,X|Z;n}$ denotes the conditional empirical distribution function, $\hat{\tilde{V}}_{*,l,\kappa}(\hat{F}_{Y,X|Z;n})$ and $\hat{\tilde{V}}^*_{l,\kappa}(\hat{F}_{Y,X|Z;n})$ denote the value functions of (ref) when replacing $F_{Y,X|Z}$ by an estimator $\hat{F}_{Y,X|Z;n}$ and $V_{*}(F_{Y,X|Z})$ and $V^*(F_{Y,X|Z})$ denote their counterparts in the population.
The first result concerns the consistency and is based on the Glivenko-Cantelli theorem \citep*[Theorem 19.1]{van2000asymptotic} which provides the convergence $\hat{F}_{Y,X|Z;n}$ to $F_{Y,X|Z}$ in $L^\infty([0,1]^2)$-norm.
Similar to the consistency result, the derivation of the large sample distribution follows form Donsker's theorem \citep*[Theorem 19.3]{van2000asymptotic} in combination with standard sensitivity arguments in optimization problems \citep*{bonnans2013perturbation} and the functional delta method \citep*[Theorem 2.1]{shapiro1991asymptotic}. This is captured in the following
Here, $dV_{*,F_{Y,X|Z=z}}(F'_{Y,X|Z=z})$ is the directional Hadamard derivative of $V_*$ defined by (ref) at $F_{Y,X|Z=z}$, $\mathbb{G}_{F_{Y,X|Z=z}}$ is a Brownian bridge with covariance function \[\text{Cov}_{\mathbb{G}_{F_{Y,X|Z=z}}} = F_{Y,X|Z=z}(\min\{y,y'\},\min\{x,x'\})-F_{Y,X|Z=z}(y,x)F_{Y,X|Z=z}(y',x')\] for all $(y,x),(y',x')\in[0,1]^2$, $z\in\mathcal{Z}$, and “$\rightsquigarrow$” denotes weak convergence.
The directional Hadamard derivative takes the form \[\delta_{F_{Y,X|Z}}V_{*}(F) = \min_{h\in\mathcal{S}(F_{Y,X|Z})}\int2\lambda\left\langle F, F_{Y,X|Z}-\mathds{1}_{[0,\cdot]\times[0,\cdot]}(Y_{X_z(w)}(w),X_z(w))\frac{dP_W}{dP_0}(w)\right\rangle dP_0(w),\] where $\mathcal{S}(F_{Y,X|Z})$ is the solution set of $V_{*}(F_{Y,X|Z})$, i.e. the set of all $\frac{dP_W}{dP_0}(w)$ with $P_W,P_0\in\mathscr{P}^*(\mathcal{W})$ that solve (ref) for $F_{Y,X|Z}$.\footnote{Proposition (ref) works equally well with a smoothed estimator $\hat{F}_{Y,X|Z=z;h_n}$ of $F_{Y,X|Z=z}$ and bandwidth $h_n$. The only important requirement is that the respective empirical process converges to a Brownian bridge.} $\left\langle f,g\right\rangle$ is the inner product in $L^2([0,1]^2)$, i.e. \[\left\langle f,g\right\rangle\coloneqq\int_{[0,1]^2}f(y,x)g(y,x)dydx.\]
(ref) requires that the sample size of the paths grows faster than the data, and the wavelet approximation $\kappa$ is also sufficiently fast given the continuity properties of the paths. Even though the large sample distribution of the value functions is not a standard Brownian bridge process, it still has a relatively common form from a purely statistical perspective, as it takes the form of the first-order directional Hadamard derivative of the value function taken at $F_{Y,X|Z=z}$ in directions $F'_{Y,X|Z=z}\in\mathcal{F}_{Y,X|Z=z}$.
In addition, there are several results in the literature (\citeauthor*{dumbgen1993nondifferentiable} dumbgen1993nondifferentiable, \citeauthor*{fang2018inference} fang2018inference, \citeauthor*{hong2018numerical} hong2018numerical) which establish bootstrap methods for estimating this type of large sample distribution in practice. In particular, they deal with general directional Hadamard differentiability \citep*{shapiro1991asymptotic}, which conforms with Proposition (ref), so that these subsampling/bootstrap results are directly applicable to the problems (ref). These bootstrap-type arguments are convenient in models with a light computational burden mostly. In more complex models one should use the analytically derived large sample theory.
The programs (ref) are semi-infinite programs \citep*{anderson1987linear}, which naturally reduce to finite dimensional problems in practice by approximating the space $[0,1]^3$ where $Y$, $X$, and $Z$ live. One can do this in two general ways. The first is to simply evaluate $\hat{F}_{Y,X|Z=z;n}$ on the values taken by the sample $(Y_i,X_i,Z_i)_{i=1,\ldots,n}$. The second is to evaluate $\hat{F}_{Y,X|Z=z;n}$ on a finite grid that spans $[0,1]^3$. This article focuses on the latter part as a grid approach gives more flexibility with respect to the computational requirements: one can make the grid coarser or finer, depending on the available memory.\footnote{This is also a difference to the computational approach in \citet*{kitamura2018nonparametric}, who set up their problem based on the observed realizations in the data and not a grid. This is a direct consequence of their finitary approach. In contrast, the infinite approach in this article allows arbitrary discretization of $Y,X$, and $Z$ and is therefore also a potential complementary approach to \citet*{kitamura2018nonparametric} in their setting.} Throughout this section, the index $\iota$ captures the degree of approximation of the grid. For instance, $\iota=11$ means that this approximation decomposes the unit interval into $11$ points $0,0.1,0.2,\ldots,0.9,1$, which will be taken to be equidistant without loss of generality. Also, and without loss of generality, all three intervals for $Y$, $X$ and $Z$ are decomposed in the same way, so that $\iota$ is the only necessary parameter controlling the approximation.
The practical implementation deviates from the theoretical approach in that it uses a smoothed variant $\hat{F}_{Y,X|Z=z;h_n}$ of the empirical conditional cumulative distribution function, where the bandwidth is determined via cross-validation. Heuristically, it seems as though the introduced smoothness gives more robust results compared to the standard empirical cumulative distribution function.\footnote{For the practical estimation of $\hat{F}_{Y,X|Z=z;h_n}$ the method uses the “np”-package in $R$ \citep*{hayfield2008nonparametric} with a standard cross-validated bandwidth.}
Under a given finite approximation, the programs take the form
where $\mu$ is a $l\times 1$ vector which corresponds to the Radon-Nikodym derivative $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$ with row-dimension equal to the number of sampled paths $l$,\footnote{Note that all elements in $\mu$ must lie in $[0,1]$, as $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$ is defined on the finite and discrete space of $l$ paths which were sampled by some $\hat{P}_0$. This means that $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$ can only put non-negative probabilities of at most one on the occurrence of each path. Intuitively, this follows from the fact that $\frac{d\hat{P}_W}{d\hat{P}_0}(i)$ is a probability mass function. These bounds on $\mu$ are included as additional constraints using $\vec{1}$.} $\vec{1}$ denotes the vector of the same dimension as $\mu$ containing all ones, $\Xi'$ is a $1\times l$ vector corresponding to $f(Y_x,x_0)$, and $\|\cdot\|_2$ denotes the Euclidean norm. $A'$ denotes the transpose of the matrix $A$. $\tilde{\Theta}$ is a $\iota^3\times l^2$-matrix which maps the realization of the stochastic processes to the distribution $\hat{F}_{Y,X|Z;h_n}$. The $L^2([0,1]^2)$-norm from (ref) reduces to the Euclidean norm due to the approximation of $[0,1]^3$ by a finite grid.
The choice of the Euclidean norm $\|\cdot\|_2$ for the constraint is convenient, as (ref) can be rewritten as
The programs (ref) are quadratic due to the Euclidean norm used and can easily be solved. This article uses the alternating direction method of multipliers (ADMM) (\citeauthor*{boyd2011distributed} boyd2011distributed and \citeauthor*{parikh2014proximal} parikh2014proximal) for optimization. This algorithm is known to converge rather quickly to reasonable approximations of the optimum, which makes it a perfect tool for this purpose. The algorithm requires two more parameters, the augmented Lagrangian parameter $\rho$ and an over-relaxation parameter $\zeta$, which control the convergence of the ADMM algorithm to the optimum. In practice, an over-relaxation parameter of $\zeta=1.7$ and an augmented Lagrangian parameter of $\rho$ between $100$ and $500$ leads to fast and robust convergence.
The computational bottleneck in a practical implementation is the construction of the matrix $\tilde{\Theta}$, whose dimension grows exponentially with the granulation of the grid $\iota$ and the number of paths sampled $l$. Fortunately, since the penalty terms in (ref) and (ref) include indicator functions, the matrix $\tilde{\Theta}$ takes the form of a binary sparse matrix: for each point $(y_\iota,x_\iota,z_\iota)\in[0,1]^3$ in the grid a given combination of paths $Y_x(i)$ and $X_z(i)$ either gets assigned a $0$ if they jointly “do not go through” the intervals $[0,y_\iota]\times [0,x_\iota]$ for given values $z_i$ or a $1$ if they jointly do. This sparseness is helpful as sparse matrices can be stored efficiently. In addition, the process of setting up $\tilde{\Theta}$ can be parallelized, which abates the computational costs even further if the researcher has access to several cores.
In many cases, however, a researcher might only have access to computational resources with very limited working memory. In such situations, it is still possible to apply the proposed method by a “sampling trick” which trades off memory requirements for time. In particular, the idea is to iteratively (i) sample with replacement a relatively small initial number $l_{0}$ of paths (depending on the available memory), (ii) optimize the programs (ref) on this sample, (iii) obtain the value functions as well as the optimizers $\mu$, (iv) drop all paths which were assigned a probability of (close to) $0$ by the optimizer $\mu$, (v) sample another relatively small number $l_{s}$, add these paths to the already existing paths and go back to (ii). The idea of this “sampling trick” is that paths which were assigned a probability of (close to) $0$ by the optimal $\mu$ do not matter for the optimal value. By dropping these paths before sampling new ones, the memory requirements do not grow or only grow modestly in practice---at the additional cost of having to run this optimization for many iterations.\footnote{Discarding elements ex-post in optimization routines is not new. For instance \citet*{wu2001solving} use discards in solving the general capacity problem on Euclidean state space. Currently, the implemented approach is rather na\"ive as it performs a “random grid search” over the infinite dimensional path space. A more efficient implementation uses ideas from sequential MCMC approaches (see schweizer2012thesis schweizer2012thesis for an overview): in the first iteration, randomly sample paths. Then optimize and discard all paths which were assigned a probability of zero. When sampling new paths, do not just sample randomly, but sample a fraction of paths which are close to the paths that were assigned a positive probability and sample another fraction of paths randomly to find other, different binding constraints. An interesting question in this setting is to optimize this procedure, i.e. to optimize the fraction of randomly sampled paths compared to paths which are close to others. Intuitively, there seems to be an “exploration-exploitation-tradeoff”, which one could optimize.}
We can now present our algorithm for solving the infinite dimensional linear programs approximately in practice.
When applying this sampling trick, the solution will be expressed as a solution path over the sampling iterations.
This solution path in general will be erratic due to the nature of the sampling approach, but has the added benefit over the “static” direct method, where all relevant paths are sampled immediately, that one can gauge if the solutions “converge” to some stable limit after a “burn-in” period. This convergence relies on the choice of the penalty parameter $\lambda$. The larger $\lambda$ is chosen, the more it forces the optimizer $\mu$ to replicate the observable $\hat{F}_{Y,X|Z;h_n}$. In fact, in the limit $\lambda\to+\infty$, the program forces to replicate the constraint perfectly. If $\lambda$ is too low, the program ignores the constraint, which will always result in trivial bounds.
This implies that there exists a range of $\lambda$-values for which constraint and objective function are balanced. Figure (ref) depicts the behavior of the solution paths of this “sampling trick” in a stylized setting of the household expenditure application in the next section, where the coarseness of the approximation of the unit interval makes it possible to sample all possible paths; this is done in order to see whether convergence still occurs if one is able to sample the “whole universe of paths”. The left panel depicts a solution path which converges, while the right panel depicts a case of non-convergence for estimating a lower bound on $F_{Y(X=0.75)}(0.75)$.
In this form, these solution paths are reminiscent of the solution paths of regularized linear programs such as LASSO. The difference, however, is that the paths induced by this program are for a fixed $\lambda$, while the actual LASSO solution paths are traced out while varying $\lambda$. In order to have an analogue of the LASSO solution paths in the current method, one would have to solve the program for many different values of $\lambda$, which would generate a system of solution paths. Then one could choose the largest lambda for which the corresponding solution path converges to a stable value.\footnote{It is an intriguing question how to determine an appropriate $\lambda$ by data-driven methods. Such a data-driven method might open up the way for solving other infinite dimensional programs on path spaces in statistics and mathematics via a “sampling of paths”-approach.}
This section presents practical settings for testing the method. The first is a small Monte-Carlo simulation, where the focus lies on the sensitivity of the method to the choice of the number of basis functions $\kappa$ and the penalty term $\lambda$. The second is an application to real data, where the method manages to obtain informative bounds under minimal assumptions.
The idea for generating the data in this simulation exercise is to (i) generate paths of stochastic processes $Y_x$ and $X_{z}$ for a given set $\{z_i\}_{i=1,\ldots,m}\in[0,1]$, (ii) combine the two processes to a joint process $[Y,X]^*_{z}$, (iii) randomly sample points $z$ and corresponding points $(y,x)$ induced by the paths of $[Y,X]^*_{z}$ to obtain the data $(Y,X)$ for the given set of $\{z_i\}$. The processes $Y_x$ and $X_{z}$ are Gaussian processes with mean $0.5$ and a squared exponential covariance kernel of the form $k_{SE}(z-z')\coloneqq \sigma^2\exp\left(-\frac{(z-z')^2}{\ell^2}\right)$, where the length parameter $\ell$ is $0.5$ for $Y_x$ and $0.2$ for $X_z$, and the variance parameter $\sigma^2$ is $0.2$ for $Y_x$ and $0.15$ for $X_z$. All paths of these processes are restricted to lie within $[0,1]$. $2500$ paths were sampled, generating $5000$ random data points from it on a relatively coarse grid based on a dyadic approximation of order $3$, i.e. grid points at $\frac{k}{8}$ for $k=0,\ldots,8$. The values $z$ were drawn uniformly on the unit interval. This data-set has no additionally introduced randomness, to see if the method can actually obtain correct results in practice.
The goal is to estimate $E[Y(0.5)],$ which is equal to $0.5$ by construction. $\hat{F}_{Y,X|Z=z;h_n}$ is estimated by kernel density methods with a cross-validated bandwidth using the $np$-package in $R$. The only assumption on the paths made for estimation is continuity. The hat-functions defined in (ref) form the basis for different levels of $\kappa$ and $\lambda$. Figure (ref) depicts the convergence of the bounds in the “sampling trick” approach.
The left panel depicts the convergence of the bounds on $E[Y(X=0.5)]$ for $\lambda=5$ and a sieve basis which consists of a sum of the basis functions for $\kappa=1,\ldots,5$. The method seems to converge nicely to bounds which contain the true value and are actually reasonably tight, especially for such a coarse grid and the fact that no functional form restrictions are imposed on the model except continuity. The right panel depicts the same thing, only there the paths $Y_x(w)$ and $X_z(w)$ are constructed for one form of basis functions corresponding to $\kappa=8$. Clearly, the method obtains biased results in this case, as the true value does not lie between the upper- and lower bound. Intuitively, this bias stems from the fact that the paths generated for fixed $\kappa=8$, despite being continuous, are not smooth and very erratic, while the paths generated in the simulation are exceptionally smooth due to the choice of the squared exponential covariance kernel and relatively large length parameters $\ell$. This bias can be understood as a “nonparametric misspecification bias” by sampling paths which are very different from the paths in the true data generating process. Fortunately, the left panel shows that already summing over only a few different basis functions removes this bias in this stylized example, which most likely also holds in more realistic settings \citep*{chen2007large}.
Figure (ref) depicts the sensitivity of the method to choosing different penalty terms $\lambda$ in this setting. In the left panel $\lambda = 0.1$ while in the right $\lambda = 50$. The bounds in the left case are significantly wider than on the right, but both estimators are unbiased. Interestingly, and in contrast to Figure (ref), even larger $\lambda$ does not lead to erratic behavior of the paths. This is most likely a result of the smoothness of the simulated paths in combination with a very high signal-to-noise ratio of the simulated data, which implies that there always exists an optimal $\mu$ in the finite dimensional problems (ref). In contrast, when the signal-to-noise ratio is low, which is often the case in real-world applications, then the behavior of the solution paths will become too erratic as depicted in Figure (ref).
As a demonstration of its capabilities, the method estimates bounds on expenditure differences using the $1995/1996$ UK family expenditure survey. This problem is well-suited as it (i) is nonlinear with continuous variables (\citeauthor*{blundell2007semi} blundell2007semi, \citeauthor*{imbens2009identification} imbens2009identification), (ii) allows to gauge if the program actually obtains reasonable results, and (iii) provides a setting not directly related to causal inference, showing the scope of the proposed method. Therefore, in the following the focus will be on the outcomes food and leisure.
Analogous to \citet*{blundell2007semi} and \citet*{imbens2009identification}, the outcome of interest $Y$ will be the share of expenditure on a commodity and $X$ will be the log of total expenditure, scaled to lie in the unit interval. The instrument used in this setting is gross earnings of the head of the household, which assumes that the way the head of the household earns the money is (sufficiently) independent of the household's expenditure allocation; this instrument is used in both \citet*{blundell2007semi} and \citet*{imbens2009identification}. All three variables are inherently continuous which makes this problem a nice setting for demonstrating the practical implementation of the method.
The sample is restricted to the subset of married and cohabiting couples where the head of the household is aged between 20 and 55, and couples with 3 or more children are excluded. Also excluded are households where the head of the household is unemployed in order to have the instrument available for each observation. The final sample comprises 1650 observations.\footnote{The data used in both \citet*{blundell2007semi} and \citet*{imbens2009identification} is based on the $1994/1995$ UK family expenditure survey. Under the same restriction they end up with $1655$ observations.}
Table (ref) gives a summary of the outcome distributions, showing that the relative expenditure on leisure is much more skewed towards zero but with a higher variance than the distribution for food.
The only shape restriction on the instrumental variable model is continuity, i.e. $h$ and $g$ are continuous functions in $X$ and $Z$, respectively. This is a natural assumptions since Engel curves are usually believed to be continuous. It implies that almost all of the paths of $Y_x$ and $X_z$ lie in $C([0,1])$, the space of continuous functions. No other assumptions are upheld. The most general current approaches either require continuity and strict monotonicity of $g$ \citep*{imbens2009identification} or of $h$ \citep*{blundell2007semi} in the unobservable $W$. In contrast, the proposed method does not require any monotonicity assumptions and hence intuitively gives an indication of how much information is available in the data to solve this problem. Surprisingly, there seems to be a substantial amount of information, as the obtained bounds indicate that food is a necessity- and leisure is a luxury good without any assumptions on the model besides continuity. Furthermore, when introducing monotonicity assumptions in the observable variables, the bounds become significantly tighter, showing the identificatory strength of these assumptions in this setting.
Figure (ref) depicts the solution paths for obtaining bounds on the counterfactual difference $F_{Y(X=0.75)}(0.15) - F_{Y(X=0.25)}(0.15)$ for a reasonably fine approximation of the unit interval into $17$ equidistant points (which corresponds to a dyadic approximation of order $4$).
Here the penalization parameter $\lambda = 1$, which seems to provide reasonable convergence to the solution, in particular for the food data. Remarkably, the estimated bounds in this setting are qualitatively informative for the problem. The left panel depicts the households' expenditures on food and the right depicts their expenditures on leisure. The solid lines are the upper- and the lower bound for a model without further assumptions, while the dashed lines are the upper- and lower bounds for a model with the additional assumption that $Y$ is increasing for leisure and decreasing for food in overall expenditure $X$ and that $X$ is increasing in income of the head of the household $Z$.
Consider the left panel first, which depicts the expenditures on food. Here, the general bounds seem to converge and the average values of the bounds over the last $200$ iterations are $0.84$ and $0.26$ for the general upper- and lower bound, and $0.78$ and $0.50$ for the corresponding upper- and lower bounds for the monotone model. All four bounds are positive, which indicates that families that spend a lot in general ($X=0.75$) and spend up to $15\%$ on food ($Y\in[0,0.15]$) would spend much more on food relatively to overall expenditure if they spent much less overall ($X=0.25$). Put differently, families are much more likely to lie in higher quantiles for expenditure on food if they lie in the lower quartile in overall spending than families that lie in the upper quartile in overall spending, which is the defining characteristic of a necessity good at this given level $y^*=0.15$.
It is rather striking that even the model without monotonicity assumptions produces bounds which reflect this fact via a positive lower bound. In this regard, note that the monotonicity assumptions do not only tighten the bounds, but also shift up the lower bound, indicating that monotonicity has a stronger identificatory content in this setting. In particular, they imply that more families (between $50\%$ and $78\%$) in the upper quartile of overall spending ($X=0.75$) spend only up to $15\%$ on food compared to families in the lower quartile of overall spending ($X=0.25$). As mentioned, without monotonicity assumptions, these differences can be as high as $84\%$ and as low as $26\%$, but all positive.
The results for expenditure on leisure for this given scenario are similar, but more erratic. For a clear indication of a luxury good at the given levels, one would expect both bounds to be negative. In fact, these would imply that families in the upper quartile on overall spending (i.e. $X=0.75$) who spend up to $15\%$ of their overall expenditure on leisure ($Y\in[0,0.15]$) are very likely to spend even less on leisure, relatively, if they had a negative shock to overall spending ($X=0.25$). Put differently, families should be more likely to spend only up to $15\%$ of overall expenditure on leisure ($Y\in[0,0.15]$) if they lie in the lower quartile in overall spending ($X=0.25$) than families that lie in the upper quartile in overall spending ($X=0.75$). The obtained results do reflect this circumstance at this level. In particular, the averages of the last $500$ iterations of the general bounds are $-0.032$ and $-0.31$, implying that typically more families (up to $31\%$) in the lower quartile of overall spending ($X=0.25$) spend only up to $15\%$ of their overall on leisure compared to families in the upper quartile of overall spending.
The overall convergence for the leisure data is more erratic compared to food, which is illuminating for the purpose of this article. In particular, the main culprit for the poor performance is the fact that the data is highly skewed towards zero in the leisure case, while the grid placed uniformly over the unit interval is too coarse to “measure” the behavior of the data around zero. Note that the method did not manage to sample enough monotone paths which correspond to the respective events of interest for this reason, so that no bounds for monotonicity exist in this setting. In order to circumvent this, one can simply put a non-uniform grid on $[0,1]$ which has more points close to zero and fewer points further away.
One can perform this exercise for different levels of $y^*$ in order to gauge the behavior of the bounds at different quantiles. Figure (ref) depicts the same exercise for $y=0.25$, i.e. looking at families that spend up to a quarter of the overall expenditures on food and leisure. Here, the method managed to sample enough relevant monotonic paths in the leisure setting despite the coarser grid around zero. These bounds converge to some value in the more general bounds. One important fact to point out is the behavior of the solution path in the leisure case. At a certain point around the $1000$ iterations mark a path is sampled which forces the solution path to jump up dramatically, i.e. a path which contains a lot of information for the given problem. However, the good news here is that over time the paths approach the previous level as more and more paths are drawn and add to the information from the one path. This is an important property to note, as it shows that the solution paths need not monotonically increase in the number of samples, and that even if paths are sampled which contain a lot of information for the respective program, will the method still converge in the long run. In particular, it is the choice of the penalty term which forces to solution path to approach zero again after the jump.
The results are qualitatively similar to the ones for $y^*=0.15$, which is not surprising as one is now comparing more families, i.e. all families that spend up to a quarter on food/leisure. The more families one compares, the less pronounced the effects become. The following table provides an overview of the results at the levels $y^* = 0.15, 0.25, 0.5, 0.75$, which shows exactly this. At $y=0.75$ one is basically comparing all existing families, so that one does not obtain any effects, another sanity check for the method.
Overall, these estimation results are remarkably informative from a qualitative perspective. Recall that the instrumental variable model allows for general unobserved heterogeneity, in particular measurement error in the treatment variable $X$, which indicates that the ratio of information to noise in the data for answering these questions is rather high. These qualitative results not only corroborate the theoretical predictions for expenditure, but also the previous results obtained in \citet*{blundell2007semi}, \citet*{imbens2009identification}, and \citet*{song2018nonseparable}. During their estimation process \citet*{imbens2009identification} and \citet*{song2018nonseparable} assume a univariate and strictly monotonic production function $g(z,W)$ between $X$ and $W$ for all $z$ and use a control variable approach to estimate the production function $h$; \citet*{blundell2007semi} estimate Engel curves semi-nonparametrically, imposing monotonicity in the second stage, and obtaining similar results. \citet*{de2016nonparametric} work with an additively separable first stage and allow for the outcome $Y$ to be measured with error. This is more general than what the proposed method can handle, which can only encompass measurement error in the dependent variable, but not the outcome. Nonetheless, their results are similar to the ones obtained here. In this sense these qualitative results are a “robustness check” for other non- or semiparametric approaches.
Moreover, this method makes it possible to gauge the identificatory content of monotonicity assumptions in the current model. In all cases is this content rather high. Imposing monotonicity between the observables makes the results much more clear-cut and in turn leads to rather strong implications in the cases considered. In this setting, monotonicity is a plausible assumption based on economic theory, but it is important to be aware of the strength of this assumption in other settings.
This article introduces a novel method for partially identifying and causal effects in instrumental variable models with general heterogeneity. It is the first practically applicable method for the most general models with continuous endogenous variables. The idea is to write the respective instrumental variable model as a system of counterfactual stochastic processes and to solve for an optimal probability measure on the paths of these processes subject to the constraint that the law of the joint processes induced by this probability measure replicates the observable distribution. The resulting optimization problem takes the form of an infinite dimensional linear program on path spaces.
The main contribution of this article is to introduce the “sampling-of-paths” approach to solve these types of infinite dimensional programs. The underlying idea is to reduce the infinite dimensional program to a semi-infinite program \citep*{anderson1987linear} by sampling a subset of the paths over which the program optimizes. Then an approximation of the (finite dimensional) state-space of the random variables leads to a finite dimensional program which can be solved efficiently. The main idea for reducing the infinite dimensional program to a semi-infinite dimensional one is to explicitly introduce randomness by sampling paths. This, in conjunction with large deviation results (\citet*{vapnik1998statistical}, \citet*{wellner2013weak}) allows to obtain probabilistic approximation guarantees of the semi-infinite program to the infinite dimensional program. In particular, these guarantees imply a lower bound on the number of paths required for achieving a good approximation with high probability.
A remaining challenge is to obtain an efficient data-driven method for choosing an appropriate penalty term $\lambda$ for the practical optimization routine. This is a similar challenge to finding good penalty terms in high-dimensional regularized regression estimators, but more general, as the setting here is infinite dimensional in a counterfactual path space. Some heuristic guidelines can be given: one should choose the largest $\lambda$ such that the “solution paths” of the sampling method converge to a fixed value after a “burn-in” period. If the solution path is “too erratic”, then one should lower the value of $\lambda$. Formally establishing what “convergence”, “burn-in period”, and “too erratic” mean would not only solve this issue, but would open up potentially novel approaches for data-driven validation approaches in counterfactual settings. In particular, an analogue to the data-driven method for $\ell1$-regularization in high-dimensional regression models as put forward in \citet*{belloni2011l1} could be valuable.
The current practical implementation of the program works for univariate variables. Moreover, the only currently implemented additional nonparametric restriction which can be placed on the model is monotonicity. The program can straightforwardly be extended to higher dimensional settings, but runs into the curse of dimensionality as the stochastic processes become high-dimensional random fields. One standard way to circumvent the curse of dimensionality is to introduce sparsity- and factor assumptions on the stochastic processes in a higher-dimensional setting. Furthermore, it is also imperative to allow for a wide variety of additional (non-) parametric assumptions in the model, like convexity, bounds, reflection processes, Slutsky-type conditions, first-passage times, martingale properties, etc. Furthermore, the current “sampling-and discarding” approaches are not efficient. From a computational perspective, it would be interesting to obtain more efficient resampling methods. One promising approach is to use ideas from sequential Monte-Carlo approaches (see e.g. \citeauthor*{schweizer2012thesis} schweizer2012thesis), which need to be extended to the infinite dimensional path spaces considered here. Another idea is to directly use diffusion- or Levy-processes in the generation of the paths and use the Karhunen-Lo\`eve transform to generate a basis for the paths.