Context-dependent Causality (the Non-Nonotonic Case)

Context-dependent Causality (the Non-Monotonic Case)

{Keywords:\,\relax } Counterfactual; Diffusion; Synthetic distribution; Unrestricted support; Contaminated Generative Adversarial Networks; Context; Non-monotonic inverse problem; Unknown kernel function; Neural Networks; Fredholm integral equation

Introduction

Consequences are propagated by actions taken in a specific and often latent contextual platform. As such, context is an intrinsic part of a causal model, playing the role of unobserved heterogeneity. The common practice in the treatment of unobserved heterogeneity is to impose various kinds of structural form restrictions on the relationships among actions and consequences. A widely used structural restriction is the monotonicity assumption, which implies a one-to-one relationship between either the action or the outcome variable and its random disturbance (given the control variables). A key limitation of this assumption however, is that it rules out any source of uncertainty regarding the unobservables when the observables are held fixed. Namely, monotonicity forces a degenerate conditional distribution of the unobservables given the observables. Consequently, the theoretical model is restricted to a subset of structural form functions and thus, may poorly identify and capture the true underlying causal relationship between action and outcome.

We depart from the conventional reliance on monotonicity to achieve identifiability. Yet, we keep the general architecture similar to the conventional modeling. The present model is a nonparametric triangular model with non-separable disturbances, which includes a latent context variable. Context is modeled as a confounder with unobserved heterogeneity given the observables, challenging the conventional causal inference. The main challenge in achieving identifiability is aggravated due to the absence of an adjustment set (observed control variables) to be conditioned on in order to satisfy some conditional independence assumption. This identifiability relies on formulating the structural equations as a system evolving from noise accumulation, known as a diffusion process to account for the influence of the latent context (confounder) on the outcome. This allows us to achieve identifiability without the limitations of monotonicity assumptions. The key result of the offered identifiability approach introduced here importantly emanates from the fact that any possible set of structural functions belonging to an equivalence class of triangular models generating the data yields the same do-interventional counterfactual distribution. This result is achieved by expressing the distributional relationship between a latent context variable and a vector of observables through a system of Fredholm integral equations. The aforementioned set of equations is governed by the set of generator functions and an unknown kernel function, inducing a non-monotonic inverse problem. The role of these generator functions is two-fold: (i) to characterize an admissible kernel function induced by these generator functions and (ii) to ensure that the estimator of the unknown quantity is a continuous function of the data for any given kernel function in the equivalence class. We establish that the interventional distribution is invariant to the choice of the kernel if it belongs to any strongly complete family of densities, e.g., an infinite mixture of Gaussians. This is a significant result, as it shows that it is possible to solve a non-monotonic inverse problem even when the kernel is unknown.

The framework offered is general and may have a significant impact on a wide variety of economic and other phenomena inclusive of those intrinsically exhibiting non-monotonicity in their confounding variables. The relationship between the equivalence class of structural form functions and the causal distribution is formed by the establishment of a newly introduced augmented Fourier series expansion. This novel expansion method is designed to characterize an infinite system of linear regression equations satisfying cross-equation restrictions, induced by an equivalence class of non-separable structural form functions. These cross-equation restrictions are derived from the common triangular model assumptions, without reliance on any variant of monotonicity.

The existing literature using either local average treatment effects or triangular models architecture impose various functional form restrictions, such as conditional quantile restrictions chesher2003identification, strict monotonicity in the case of average local treatment effects angrist1994identification as well as in the case of triangular models imbens2009identification,hoderlein2017corrigendum. Alternative models assume restricted support domain (locality) hoderlein2009identification,angrist1996identification, altonji2005cross, heckman2005structural. These restrictions guarantee the identifiability of the latent confounder up to a strictly monotonic transformation from conditional quantiles of the observables. The main limitation of such treatment however, is to potentially narrow down the support domain of the latent space. This can be detrimental to causal inference as it eliminates a subset of the control (placebo) group, when employing e.g., the common practice of partial means newey1994kernel, interventions pearl2019interpretation and local average treatment effect angrist1996identification. Monotonicity is rather not an innocuous restriction as it is at odds with many phenomena in general and those in economics in particular. The potential bias propagated by relying on the monotonicity assumption has long been estimated, (e.g., klein2010heterogeneous). However, there is a lack of alternative identifiability, inference and an appropriate estimator which are not monotonicity-based. Monotonicity is a strong and perhaps unpalatable for many basic economic phenomena such as demand, supply and consumption behaviour hoderlein2016testing, hoderlein2007identification, multiple equilibria in strategic behavior myerson1999nash, search behavior chetverikov2019testing, agarwal2020searching as well as in prospect theory tversky1979analysis.\footnote{A recent paper dembo2021ever documents strong violations of monotonicity as a feature of expected utility theory.}

In our present framework, we depart from the above mentioned restrictions by allowing for the posterior distribution of the latent context to be a function of the covariates and of the outcome. This generalization incorporates cases in which the context is also a function of the outcome and not only the covariates, unlike the control variable approach \myciteopt{mammen2012nonparametric,mammen2016semiparametric,chernozhukov2020}. The presence of a latent context covariate affects the cause and the outcome simultaneously and thus, is the source of endogeneity in our model. Thus, the endogeneity can be eliminated by allocating each one of the cause and outcome a different and independent realization of the context. Building upon this idea, we develop a novel identification strategy as well as a new estimator for triangular models in the presence of non-separable disturbances. Unlike the common practice, our approach does not rely on the strict monotonicity assumption, but rather on an integral equation characterized by generator functions. These generator functions induce an equivalence class of triangular models, while exploiting the entire support domain of the latent space in a non-monotonic manner. This insight solves the double-hurdle problem of relying on monotonicity as well as on a restricted support domain chernozhukov2020 which may produce inferior results as is shown in our simulations.

The new paradigm presented here yields a synthetic counterfactual distribution, which is entropy-advantageous. We provide a novel feed-forward multi-layer Neural network-based estimator. This is performed by leveraging on a novel extension namely, a Contaminated version we develop and apply to the Generative Adversarial Networks model (CONGAN) goodfellow2014generative. The advantage of the proposed practice is that it enables the Neural network to be trained on the entire support domain of the latent space unlike the common practice, which uncovers counterfactual relationships by averaging available observable data on a restricted support.

The paper proceeds as follows. Section (ref) gives an overview of the common practice of causal inference in triangular models. Section (ref) presents the diffusion problem formulation and its underlying assumptions. Section (ref) shows an alternative formulation of the triangular model as an integral equation to simplify the analysis. Section (ref) provides an illustrative example for identifiability which does not rely on monotonicity assumption. Section (ref) characterizes strong completeness of sufficient statistic as a building block for identifiability. Section (ref) presents the identifiability by relying on strong completeness. The estimator is presented in section (ref). In section (ref) Monte-Carlo simulations are used to verify our theoretical model performance. Section (ref) concludes.

Discussion and literature review

Various identification strategies have been designed to uncover counterfactual relationships in triangular models, consisting of non-additively separable disturbances, in the absence of experimental data. The ultimate goal of these approaches is to mimic a “natural experiment” rosenzweig2000natural by detaching a random variable $X$ (an observed action or treatment) from the unobservables affecting some other random variable $Y$ (outcome). The prominent question they answer is what would have been the expected counterfactual outcome of $Y$ if the values of the independent variable $X$ were assigned independently of the disturbances. Generally, in order to achieve this goal, monotonicity is imposed either on the unobservables in the first stage regression \myciteopt{imbens2009identification,blundell2014control, heckman1985alternative,hoderlein2009identification,chernozhukov2020}; or on the unobservables in the outcome's structural equation in the second stage chernozhukov2007instrumental, kim2020partial. In other cases, the monotonicity requirement is imposed on the instrumental variable angrist1994identification,hoderlein2017corrigendum,heckman2018unordered. However, this assumption is at odds with many actual phenomena as discussed earlier. Moreover and perhaps more importantly, it imposes severe limitations on the set of identifiable causal relationships in nonparametric inference, while potentially missing the true causal relationship. The unattended issue in the existing literature is that the direction and strength of causal relationship is context-blind.

An alternative identifiability strategy for causal inference not relying on monotonicity is obtained by assuming ignorability which is known also as an unconfoundedness assumption athey2016recursive or a conditional independence on observables farrell2021deep. This approach limits the scope of causal analysis by requiring a directly observed set of control variables, which are rarely available in reality, in order to introduce conditional independence. The present study departs from this approach by embracing a generative covariate mechanism mammen2012nonparametric,mammen2016semiparametric producing a synthetic sample of the same distribution present in the data. This approach does not necessitate the unconfoundedness assumption, yet allows for the estimation of the counterfactual distribution at any level of accuracy desired (given amount of data available and level of complexity). This methodology also departs from the popular quantile structural estimation, which is intrinsically monotonic chernozhukov2013inference.

We contribute to the existing literature in two specific aspects. The theoretical contribution amounts to developing a new causal identification paradigm, which is non-parametric, monotonicity-free, context-dependent and universal with respect to the identifiable function set. Our new paradigm yields a counterfactual distribution exploiting the unrestricted support domain of the latent space, which is essential for the practice of intervention and counterfactual inference. On the methodological front, we develop a novel estimation strategy referred to here as a Contaminated Generative Adversarial Networks (CONGAN). This approach consists of two generators and one classifier: counterfactual samples generator and a contaminator generator. The classifier determines whether the sample is real or synthetic as in the conventional GAN goodfellow2014generative. The main difference is in the unique role of each one of the generators in the causal inference. The counterfactual samples generator intends to mimic the causal relationship between $X$ and $Y$, which would be the case if $X$ and the unobservables were distributionally independent. However, it cannot be trained directly using the observed data, in which $X$ and the latent context variable are distributionally jointly dependent. Consequently, the role of the contaminator is to minimize the disparity between the observed and the counterfactual samples. As such, the contaminator is a nuisance Neural network which is essential for training the counterfactual generator. Our new paradigm enables to mimic the distribution of $Y$ given $X$ that would have been obtained under random assignment of the latent context variable (independently of $X$). We emphasize at this stage that context is not just a semantic concept but rather an identifiability instrument, in that it permits rendering $Y$ a different contextual regime from the one governing $X$. Once the regime governing $Y$ does not longer comove with $X$, one mimics the effect of an exogenous variation in $X$ on $Y$. The benefit of the contaminated Neural network is in that it precludes the need to invert the generators in order to account for the posterior distribution of the latent context variable given the observables, which is computationally cumbersome. Our CONGAN application utilizes multi-layer feed-forward Neural networks in which, unlike in other series estimators (e.g., sieve), the basis functions themselves are data-driven by combining simple functions. Such a flexible combination is known to be capable of approximating any measurable function to any desired degree of accuracy hornik1989multilayer.\footnote{We note that a specific variant of GANS, referred to as Wasserstein GANS with a penalized gradient, has been recently employed for conducting monte-carlo simulations of potential outcomes under unconfoundedness assumption athey2020. Recent work criticizes this penalized Wasserstein GANS due to disregarding significant parts of the support domain in the data wei2018improving. Note that the utilization of the correct support is at the core of counterfactual analysis.} In the ensuing section we attend to the theoretical triangular model formulation.

Diffusion-based non-monotonic triangular models

Consider the following triangular simultaneous equations model, pioneered by heckman1985alternative. The model includes three observed random variables $X,Y,Z\in\mathbb{R}$, and two unobserved random disturbances $\upeta\in\mathbb{R}$ and $\bm{\epsilon}:=(\epsilon_1,...,\epsilon_d)\in\upchi^d\subset\mathbb{R}^{d}$. The triplet $(Z,\bm{\epsilon},\upeta)$ generates a pair of observed quantities $X,Y\in\mathbb{R}$, as an infinite Gaussian mixtures governed by unknown location as well as variance functional parameters $(\mu_{1,i}, \mu_{2,i})$ and $(\sigma_{1,i}, \sigma_{2,i})$. There are unknown mixing proportion parameters $(\alpha_{1,i}, \alpha_{2,i})$. Let $L$ be the number of components in the Gaussian mixtures of $(X,Y)$, which approaches infinity, depicted as,

and $Y = g(X,\upepsilon)$ satisfying,

The random variable $Y$ is the outcome variable, whereas $X$ represents an action or an endogeneous variable. The random variable $Z$ is an exogenous covariate also known as an instrumental variable. The random variable $Z$ is assumed to be continuous with an unknown density $f_Z(z)$, and is further assumed to satisfy $Z\protect\mathpalette{\protect\independenT}{\perp} (\upeta,\bm{\epsilon})$. The disturbances $(\upeta,\bm{\epsilon})$ are assumed to be continuous and possibly jointly dependent. We denote their unknown joint and marginal distributions by $F_{\upeta,\bm{\epsilon}}$, $F_\upeta$ and $F_{\bm{\epsilon}}$, respectively.

Equivalently, the Gaussian Cholesky decomposition parameterized by the functional parameters $\left\{\upbeta_i\right\}$ with $\upbeta_i(x):=\rho_{i,13}^2(x)$ and a constant $\uptheta_i:=\rho_{12}^2$ can be used to simplify the analytical formula of the data Gaussian mixture data generation process in Eqs. (ref) and (ref),

with independent uniformly distributed random variables $\bm\upnu_i$ and $\bm\upnu_{i,x}$. We define $\alpha_{i,j}^*:=\alpha_{1,i}\alpha_{2,j}$ to express the observed outcome as a function of $(X,Z)$,

where $\Psi_{j}(X,Z)$ for any given noise $\bm\epsilon_0$ is just a one of the infinite potential values of $\upeta$ in the $j$'th component of the mixture for any given $(X,Z)$. Eq. (ref) is a reduced-form representation of $Y$ without any endogenous predictor variables. This formulation shares similarities with diffusion models ho2020denoising by capturing the influence of noise on $Y$.

We attend to the following definition of the support of a continuous random vector $Q$ in $\mathbb{Q}^n$ with a probability density function $f_{Q}$,

We consider a non-parametric setting, whereby the unknown function $h$ is assumed to be continuous, but otherwise of general structural form. The unknown function $g$ can be continuous or piecewise constant, in the latter case leading to a discrete outcome $Y$. As in {imbens2009identification} the functions $h$ and $g$ may be nonseparable in their respective disturbances $\upeta$ and $\bm{\epsilon}$. The model ((ref))-((ref)) implies that $\upeta$ is a confounder or common cause, affecting both $X$ directly via (ref), and $Y$ indirectly through the dependence of $\bm{\epsilon}$ on $\upeta$. Finally, following previous work imbens2009identification, we assume the following standard condition holds:

For future use, we present two important properties of the above triangular model. First, note that Eq. ((ref)) together with the assumption that $Z$ is independent of $(\upeta,\bm \epsilon)$ imply that $X\protect\mathpalette{\protect\independenT}{\perp} \bm\epsilon|\upeta$. The second property is stated in the following auxiliary lemma (see proof in appendix (ref)).

Given $n$ i.i.d. triplets $(x_i,y_i,z_i)$ from the above triangular model, Eqs. ((ref))-((ref)), the goal is to estimate the counterfactual distribution of the outcome $Y$ if $X$ were set to the value $x$. Depending on the particular application, this quantity is related to policy effect or treatment effect. In the notation of causal inference pearl2019interpretation, the above quantity of interest is known as intervention. It is given by

It is worth emphasizing that under the above triangular model and its assumptions, even exact knowledge of the joint distribution of $(X,Y,Z)$ does not uniquely determine the quantity $T(y|x)$ imbens2007nonadditive. The reason is that even if $g$ were known, under ((ref))-((ref)) the random variables $X$ and $\bm\epsilon$ might be jointly dependent. In contrast, in Eq. ((ref)) the expectation is over the marginal distribution of $\bm\epsilon$, effectively treating $\bm\epsilon$ as independent of $X$.

To guarantee identifiability of $T(y|x)$, several previous works imposed an additional assumption of {\em strict monotonicity} of $h(z,\eta)$ in the second argument $\eta$, required to hold for all $z$ in its support, e.g., blundell2014control,imbens2009identification. Specifically, imbens2009identification defined the following random variable $V:=F_{X|Z}(X)$ and the following functional, both of which can be estimated from the observed data,

The strict monotonicity assumption implies that the random variable $V$ is a one-to-one mapping of the unobserved disturbance $\upeta$. Furthermore, it can be shown that $V$ is a {\em control variable} and that $H_c(y|x) = T(y|x)$. Then, an estimate of $H_c$ directly yields an estimate of $T(y|x)$.

The identifiability strategy described above is based on the ability to uniquely determine $\upeta$, up to a strictly monotonic transformation, for any pair of observed $(X,Z)$. As we now describe, we significantly depart from this approach. Instead, we leverage on the known axiom of continuity and building block of expected utility theory and post the following assumption:\footnote{We are reminded that continuity is an essential part of the axiom of expected utility theory guaranteeing non-intersecting indifference curves. See afriat1967 for the well-known result showing that data cannot be treated as being generated by a utility function, if there are large deviations from rationalizability. \mycite{afriat1967} theorem tells us that continuity is utmost necessary for any finite data set to be rationalizable. }

{

The first key contribution of our work is to show that under the triangular model with assumptions (ref)-(ref), the intervention $T(y|x)$ is identifiable. Consequently, the fact that additive and monotonic triangular models are identifiable, are special cases of our more general identifiability result imbens2007nonadditive,imbens2009identification,blundell2014control,heckman1985alternative.

We establish identifiability in a population setting, where assume to have observed an infinite number of triplets $(x_i,y_i,z_i)$ from the triangular model (ref)-(ref). Hence, in what follows we assume that the joint distribution $F_{X,Y,Z}$ as well as various conditional distributions, such as $F_{Y|X=x,Z=z}$ are all perfectly known. Our approach to prove identifiability proceeds as follows. First, in Section (ref), we derive an integral equation relating the intervention of Eq. ((ref)) to known distributions of $(X,Y,Z)$. In general, there may be an infinite number of solutions to this integral equation, which reflects the fact that in the original triangular model, the unknown functions $h$ and $g$ are not identifiable. Yet, in Section (ref) we prove that {\em any} solution of this integral equation gives the same intervention, hence proving identifiability of the intervention.

An alternative representation for the intervention

A key preliminary step is to derive an {\em alternative} representation for the intervention $T(y|x)$. {Consider the following cumulative distribution function (CDF),}

Note that the function inside the integral , with $\upeta$ attaining all possible values in its support, is well defined by assumption (ref). The following lemma is key to our proposed approach, both for proving identifiability and for inference.

By Lemma (ref), instead of estimating $T(y|x)$ which involves a $d$-dimensional integration over $\bm\epsilon$, one may instead estimate the function $H(y|x)$ that depends on a univariate integral w.r.t. $\upeta$. The specific expectation operator in Eq. ((ref)) is known as the partial means estimator newey1994kernel, since the expectation is taken w.r.t. the unconditional distribution $F_{\upeta}$ rather than the conditional one $F_{\upeta|X=x}$. As such, $H(y|x)$ is a {\em counterfactual} distribution.

An Example for a Triangular Data Generation Process

Let's examine the non-monotonic case, $h(z,\eta)\ne F_{X|Z=z}^{-1}(F_{\upeta}(\eta))$.

The coefficient $\rho_{13}(x)$ determines the conditional correlation between $Y$ and $\upeta$ given $X=x, Z=z$.

Recall that $\upeta\protect\mathpalette{\protect\independenT}{\perp} Z$. Yet, $\upeta\not\protect\mathpalette{\protect\independenT}{\perp} Z|X$. Similarly, $Y\protect\mathpalette{\protect\independenT}{\perp} Z|X,\upeta$. In the equation of $X$, each value of $Z$ characterizes a specific conditional Gaussian (given $\upeta$). In the equation of $Y$, each value of $X$ characterizes a specific conditional Gaussian (given $\upeta$). We employ Gaussian Cholesky decomposition to simplify the analytical formula of the data generation process,

where $\upeta(x,z):=\upeta|X=x, Z=z \sim N\left(\mu_{\eta}(x,z), \sigma_{\eta}^2(x,z)\right)$. The unknown parameters characterizing the latent confounder $\upeta$ are: $(\rho_{12}, \rho_{13}(\cdot), \mu_{\eta}, \sigma_{\eta})$. The rest of the parameters can be directly measured from the joint distribution of the observables $(X,Y,Z)$.

Identifying $\rho_{12}$ from the covariance with observable.

This gives,

Given $\rho_{12}$, the coefficient $\rho_{13}(x)$ can be identified by a linear regression of the form,

The counterfactual outcome variable

Define $\upbeta:=\rho_{13}^2(x)$. For any given pair $(x,z)$, $Y$ can be reconstructed from a noise propagated by $\bm\nu$,

The non-contaminated output is defined as follows, $$Y_0:=\Phi^{-1}(\bm\nu).$$ The observed output is a contaminated variant of $Y_0$ is known as a diffusion process ho2020denoising,\footnote{Diffusion models excel at generating high-quality, diverse data for various applications, including creative content generation. They achieve this by adding controlled noise to the data and subsequently learning to reverse the process, effectively capturing underlying structure in datasets.}

It can be seen that the posterior distribution of $\eta$ given $x$ and $z$ does not play any role in $Y_0$. In the present diffusion model we only need to recover the parameters $(\mu_2(x), \sigma_2(x))$ characterizing $Y(x,z)$ for any give pair $(x,z)$. This requires to identify the nuisance parameters $\beta$ and $\rho_{12}$ as depicted above. Eq. (ref) gives us the direct linkage between the identified parameters and the counterfactual output. Namely, by substituting $\eta(x,z)$ in Eq. (ref) with a random draw from the prior distribution of $\upeta$ we get a counterfactual sample from $Y^{\text{CF}}$.

The next building block is formulating a two-regime model describing the relationship among members in a complete family of distributions, facilitating generalizing the identifiability achieved in the simple example above.

Completeness of Sufficient Statistic

In simpler terms, if the expected value of any function of the statistic is zero for all possible values of the parameter $\theta$, then the function itself must be almost everywhere equal to zero.

Now, we attend to a stronger variant of completeness, bridging between regimes nested in the same model.

Informal Definition (Strong completeness). The concept of strong completeness applied to the relationship between regimes (subsets of a parameter space), introduced by alamatsaz1983completeness, asserts the invariance of expected values of functions with respect to a parameter in some regime (for some fixed setting of the remaining parameters), implying invariance with respect to this parameter in any other regime (other fixed settings of the remaining parameters).

In the figure above we present a two-regime model parameterized by $(\theta_1,\theta_2)$ belonging to a strongly complete family of distributions with respect to $\theta_1$. Under strong completeness, family of regime $B$ (each member is a different realization of $\Theta_1$) is informative for policy makers that want to infer about family of regime $A$ because these two families share a common parameter space $\theta_1$. We next attend to a formal definition of strong completeness.

Consequently, in the presence of a model related to a strongly complete family of distributions with respect to common parameters $\theta_1\in\Theta_1$, invariance to these common parameters in one regime $\theta_2=\theta_2^{\prime\prime}$ holds also in the rest of regimes $\theta_2^{\prime}\in\Theta_2$. This enables to inform policy decisions in different real-world settings (regime $A$).

Identifiability

In this section we prove the identifiability of the intervention $T(y|x)$. This is done in two steps. First, we alleviate the issue of the triangular model dependence on multiple unknown quantities. The most prominent one being the dependence of the disturbances on the endogenous covariate. This is implemented by employing an equivalent representation for the structural functions (lemmas (ref) and (ref) to follow) in a canonical form, which is invariant to the specification of the unknown quantities (joint-distribution of the disturbances). Second, we use canonical form representation to characterize an equivalence class of normalized-triangular models and establish that identifiability is universal in that it holds for any structural continuous function in the equivalence class (theorem (ref)).

Normalized triangular models via error-decoupling mechanism

In what follows, ${X}^*\overset{d}{=}X$ indicates that the two random variables $X^*$ and $X$ have the same distribution. Using this notation, the triangular model is reformulated via an error-decoupling mechanism with disturbances that are uniformly distributed. To this end, given random variables $\upomega$ and $\nu_1,\ldots,\nu_d$, define the following transformed random variables,

and recursively for $1\le j\le d$,

The notations in Eqs. (ref)-(ref) are used to perform error-decoupling in the following lemma (see proof in appendix (ref)):

There are two major differences between the normalized triangular model Eqs. (ref)-(ref) and the original Eqs. (ref)-(ref). The first difference is that the disturbances $\upomega$ and $\bm \nu$ are {\em independent} and in fact uniformly distributed. The second difference, or the price to pay for this decoupling, is that now the disturbance $\omega$ appears explicitly both in the equations for $X$ and for $Y$, namely inside the functions $h^*$ and $g^*$. Nevertheless, this price is negligible as it alleviates a major challenge in the original triangular model of Eqs. (ref)-(ref). This challenge stems from the fact that $\upeta$ and $\bm \epsilon$ may in general be dependent, in an unknown fashion.

Identifiability of the interventional distribution

Recall that our primary goal is to estimate the intervention $T(y|x)$ of ((ref)). The key difficulty is that observing $(X,Y,Z)$ does not uniquely determine the functions ${h}^*$ and ${g}^*$. We alleviate this issue by introducing a set of functions, which are all indistinguishable with respect to the observations. To this end, let $\tilde{F}$ be an arbitrary continuous distribution. Define

where $\bm\nu$ and $\widetilde{\upomega}$ are independent. In other words, the set $\widetilde{\cal C}$ is the {\em equivalence class} of the given triangular model. Namely, all pairs of structural functions $\widetilde{h}$ and $\widetilde{g}$ give rise to the same distribution for the observables $(X,Y,Z)$. For any pair $(\tilde{h},\tilde{g})$ there is a corresponding function

Note that any pair of structural functions $(\widetilde{h}, \widetilde{g})\in\mathcal{\widetilde{C}}$ constitutes a distribution,

as well as a counterfactual distribution, \[ F_{X,Y|Z=z}^{\text{CF}}(x,y):=\mathbb{E}\left[\mathds{1}\left\{\widetilde{h}(z, \upomega)<x, \widetilde{g}(x, \upomega,\bm\nu) < y\right\}\right],\] such that the counterfactual density is $f_{X,Y|Z=z}^{\text{CF}}(x,y):=\mathbb{E}\left[\delta(\widetilde{h}(z, \upomega)-x) \delta(\widetilde{g}(x, \upomega,\bm\nu)-y)\right]$. In next steps, this density of interest is shown to be identifiable.

A sufficient condition for identifiability in equivalence class $\mathcal{\widetilde{C}}$ is to uniquely determine the Fourier transform of the counterfactual density from the joint distribution of the triplets $(X,Y,Z)$, namely, having that,

We rely on the following lemma as a building block in establishing that this condition is satisfied under assumptions (ref)-(ref), by formulating the observational equivalence condition in terms of Fourier transform.

The following theorem (ref) states that for any two pairs of functions $(g^*,h^*)$ and $(\widetilde{g}, \widetilde{h})$ satisfying observational equivalence (Eq. (ref)), counterfactual equivalence also holds.

The theorem below relies on theorem (ref) to show that Eq. (ref) holds. Simply put, it states that the interventional distribution is invariant with respect to the choice of the structural form functions in the equivalence class $\widetilde{\mathcal{C}}$ (see proof in appendix (ref)).

The aforementioned may give rise to the following query: can a monotone data generation process produce an equivalent distribution to the one propagated by a non-monotonic data generation process? The answer to such query is definitely negative (as is proven here in lemma (ref), Appendix (ref)). This is so because any monotonic data generation process lacks sufficient variation in the distribution of the outcome variable due to the restriction on any value of context producing the same distribution of $Y$ within a fixed conditional quantile of $X$.

So far we have established identifiability for the interventional distribution. This is a very strong result in that it is universal, neither relying on any structural form restrictions nor on any assumptions regarding the true measure (prior to normalization) of the disturbances. Thus, it lends itself to nonparametric treatment. Building on the above identifiability result of the counterfactual distribution, we next develop the synthetic counterfactual machinery. Theorem (ref) establishes that the kernel can be unknown in Fredholm integral equations, whereas Theorem (ref) establishes the identifiability of the counterfactual distribution based on the insight of theorem (ref).

Synthetic counterfactual machinery

Following theorem (ref), the interventional distribution $H(y|x)$ is the same for any sequential generator functions in the equivalence class $(h^*,g^*)\in\widetilde{\mathcal{C}}$. This result naturally lends itself for an arbitrarily chosen sequential generator functions $(\widetilde{h},\widetilde{g})$ for defining,

which is the observed (O) conditional generated density function of $Y$ and $X$ given $Z=z$. Similarly, the observed (O) conditional generated density function of $X$ given $Z=z$ is defined as,

Note that the key innovation embedded in (ref) is in that the entire support domain $[0,1]\times[0,1]^d$ of the latent space (random disturbances) is used to characterize the conditional joint distribution of $Y$ and $X$ given $Z=z$. This is superior to the concept widely used in the literature for characterizing $F_{Y|X,\upomega}(y|x,\omega)$, the conditional distribution of $Y$ given $X=x$ and the actual latent space (up to normalization), which is restricted by the available data. The present treatment in (ref) however, conveys an entropy superior way to float the interrelationships among the unobservables and the observables, as it leads to the following general characterization of the unobserved counterfactual (CF) conditional joint distribution of $Y$ and $X$ given $Z=z$ (by letting $\widetilde{\upomega}\sim U[0,1]$ satisfying $\widetilde{\upomega}\protect\mathpalette{\protect\independenT}{\perp}\upomega$):

The increase in entropy is due to the presence of two independent contextual covariates $\upomega$ and $\widetilde{\upomega}$ in (ref) introducing independent variations to $X$ and to $Y$. This yields the following counterfactual result:

The result in (ref) is propagated by the unconstrained support domain of (ref) unlike in models employing the additive control function approach heckman1985alternative, the non-additive control function approach newey1994kernel,imbens2009identification,blundell2014control,hoderlein2009identification, the generated covariates mammen2012nonparametric, the non-separable counterfactual distribution based on quantiles chernozhukov2013inference, the treatment of unobservables in schennach2014entropic and the average causal effects by employing do-interventions pearl2019interpretation. Further discussion pertaining to the effect of limited empirical support on counterfactual analysis is detailed in section (ref).

In fact, the context in and by itself makes it feasible to allocate different values for the counterfactual data set of $X$ and $Y$. By doing so, we can generate random draws, independent of each other, from the unconstrained support of the context and the action variables, unlike the case when the observed data set is used. This is demonstrated in Figure (ref):

The observed data set is denoted by the black dots and the counterfactual data set is denoted by the gray dots. By construction, the former is randomly drawn from the joint support of $X$ and $\upeta$, while the latter is randomly drawn from the rectangular support $\text{Supp}(X)\times\text{Supp}(\upeta)$. Consequently, the partial means estimate of $\mathbb{E}[Y|\text{do}(X=x)]$ is calculated by using the observed data $X$ and $Y$ in which $X$ and $\upeta$ are jointly dependent as in Figure (ref).

Our theoretical result (theorem (ref)) implies that both the kernel as well as the integral equation can be completely characterized by generator functions belonging to $\widetilde{\mathcal{C}}$ and being incorporated in our proposed estimator (in section (ref) to follow). To show that, recall the Dirac delta function $\updelta(\cdot)$ and denote, \[ f_{Y,X|Z=z}(y,x) = \int_{0}^{1}\left\{\int_{[0,1]^d}\updelta\left(y-\widetilde{g}\left(\widetilde{h}(z,\omega),\omega,\bm\nu\right)\right) d\bm\nu\right\} \updelta\left(x-\widetilde{h}(z,\omega)\right) d\omega. \]

Next, we propose our estimator and justify the type of Neural network architecture we develop, which involves sequential generator functions. The sequential generator functions estimator is proposed to mimic the underlying data generation process as in a triangular model, where by its very definition, the structural functions are formulated in a sequential manner.

Estimation: contaminated GAN (CONGAN)

The formulation of the system of Fredholm equations in (ref) (in Appendix (ref)) is used only for achieving the identifiability of $H(y|x)$ and not for estimation. This disparity between identification and estimation is related to the fact that such a system of equations might be ill-posed. Namely, it may provide an estimator which is a discontinuous function of the data and thus, inconsistent darolles2011nonparametric, newey2003instrumental. This issue is resolved here by constructing a novel contaminated GAN (CONGAN) estimator, which is a smooth function of the data. The benefit of the contaminated Neural network is in that it precludes the need to invert the generators in order to account for the posterior distribution of the latent context variable given the observables, which is computationally cumbersome.

Consider a probability space $(\Upomega,\mathcal{F},\mathbb{P})$, a measurable output spaces $(\mathcal{X}, \mathcal{Y})$ for the random vector $(X,Y)$, an observed input space $\mathcal{Z}$ for the random variable $Z$ and a latent input space $(\mathcal{U,V})$ for the random vector $(\upomega,\bm\nu)$. The set of paired measurable functions $\mathcal{G}:(\mathcal{X},\mathcal{U},\mathcal{V})\mapsto\mathcal{Y}$ and $\mathcal{H}:(\mathcal{Z},\mathcal{U})\mapsto\mathcal{X}$ map from an input space (observed and latent space) to the output space. These functions are termed generator functions designated to produce random samples. We consider also $\mathcal{D}$ a set of measure functions $d:(\mathcal{X},\mathcal{Y},\mathcal{Z})\mapsto[0,1]$ referred to as discriminators. Their role is classifying the input space as either real or synthetic. Let $\mu$ be a measure on $(\mathcal{X}, \mathcal{Y})$, we make the following assumption singh2018nonparametric:

We define the contaminated GAN (CONGAN) approach as an extension of the original GAN goodfellow2014generative as follows:

The key difference between the original GAN and our proposed CONGAN is entirely related to the sequential construction of the generator functions (the generator function $\widetilde{h}$ belongs to the input of the generator function $\widetilde{g}$) both share the noise $\omega$ as appears under the contamination braces in (ref). Namely, (ref) is reduced to the original GAN when $\widetilde{h}(Z,\omega)$ is substituted with $X$ and the expectation of the CON braces is taken w.r.t $X$, $Z$ and $\bm\nu$ rather than w.r.t $\bm\nu,\omega$ and $Z$.

The optimal discriminator assigns values close to $1$ given data from the joint distribution of $(X,Y,Z)$ and assigns values close to zero given data from the distribution obtained by the generator functions $(\widetilde{h},\widetilde{g})$. Thus, the optimal $d$ given a fixed pair of generator functions $(\widetilde{h},\widetilde{g})$ is:

Then, the optimal pair of generator functions $(\widetilde{h},\widetilde{g})$ is the one solving the $\min\max$ problem:

Define $f_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h}):=\frac{\partial^2}{\partial x\partial y}F_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h})$ and suppose that assumption (ref) holds. Then following Theorem 2.2 in singh2018nonparametric,

$\forall \hspace{0.2em}(\widetilde{h},\widetilde{g})\in\mathcal{(G,H)}$ and $\forall \hspace{0.2em}d\in\mathcal{D}$. Hence, for any fixed pair of generators $(\widetilde{h},\widetilde{g})\in\mathcal{(G,H)}$ and $d\in\mathcal{D}$ maximization of $V(d,(\widetilde{h},\widetilde{g}))$ has the form:

This yields:

The result above represents the Nash equilibrium concept, in the sense that both the discriminator as well as the generator players reach their best responses in maximizing their value functions given the other player action goodfellow2014generative. In practice $f_z(z)$, $f_{Y,X|Z}(y,x|z)$ and $f_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h})$ are estimated by a kernel density estimator as a function of some unknown bandwidth $\bm{\upsigma}=(\upsigma_x,\upsigma_y,\upsigma_z)$. This bandwidth is a parameter of the model satisfying the aforementioned Nash equilibrium. Consequently, it obviates the known problematic cross-validation, which requires formidable successive multiple estimations of the likelihood function.

Next we attend to the characterization of our proposed Neural network to be used in the contaminated GANS methodology. For doing so, we denote a feed-forward Neural network consisting of $L$ layers (such that there are $L-1$ hidden layers). In what follows we explain in detail the formulation of the sequential generator functions in this Neural network, presented in the following figure:

Figure (ref) depicts the sequential generation of $\widetilde{X}$ and $\widetilde{Y}$ induced by a triangular model with error-decoupling. By adopting this sequential data generation process, in which $\widetilde{g}$ takes the output of $\widetilde{h}$ as one of its arguments, one can exploit the unrestricted support domain of the latent space in order to construct the observed distribution in (ref) as well as the counterfactual distribution in (ref). This latent space consists of $\upomega$ and the random disturbance $\nu$. Next, we characterize explicitly the Neural network architecture to estimate these functions.

Multi-layer neural network with sequential generator functions

Denote the input space of its $\ell$'th layer in the generator functions of $Y$ and $X$ by the column vectors $\bm{\mathcal{Y}}_{\ell-1}$ and $\bm{\mathcal{X}}_{\ell-1}$ of sizes $T_{\ell-1}^Y\times 1$ and $T_{\ell-1}^X\times 1$, respectively. The parameters associated with these covariate vectors are referred to as weights (slopes) and biases (intercepts) which are denoted by the matrices $\bm{W}_{\ell}^{Y}$ and $\bm{W}_{\ell}^{X}$ of sizes $(T_{\ell}^Y+1)\times (T_{\ell-1}^Y+1)$ and $(T_{\ell}^X+1)\times (T_{\ell-1}^X+1)$, respectively. Let $\widetilde{g}$ and $\widetilde{h}$ be the generator functions of $Y$ and $X$, respectively, with parameters $\bm{\upbeta}=(\bm{W}_{1}^X,...,\bm{W}_{L}^X)$ and $\bm{\uptheta}=(\bm{W}_{1}^Y,...,\bm{W}_{L}^Y)$, defined as:

and

where $\upvarphi(\cdot):\mathbb{R}\to C\subset\mathbb{R}$ is some differentiable activation function to be applied element-wise. The input space of the first layer is $(T_0^Y, T_0^X)=(3,2)$ and the output of the last layer is $(T_L^Y, T_L^X)=(1,1)$.

Real and synthetic data are denoted by $\left\{(X_i,Y_i,Z_i)\right\}_{i=1}^N$ and $\left\{(\widehat{X}_i,\widehat{Y}_i,Z_i)\right\}_{i=1}^N$, respectively, with $\widehat{X}_i:=\widetilde{h}(Z_i,\omega_{0i},\bm{\upbeta})$ and $\widehat{Y}_i:=\widetilde{g}(\widehat{X}_i,\omega_{0i}, \nu_{i};\bm{\uptheta})$.

It should be emphasized that by adopting a Shannon-divergence minimization criterion function, we depart from the practice of Wasserstein GAN (WGAN) athey2020. Although WGANS can be trained more easily, it requires the utilization of a penalty function to alleviate various optimization difficulties gulrajani2017improved, which may limit the support domain of the generated data distribution wei2018improving.\footnote{This limitation stems from the fact that this gradient penalty term is influenced only by the sampled points without examining a significant parts of the support domain wei2018improving.} As we deal with counterfactual analysis we aspire to utilize the entire support domain of the data distribution function. Thus, we adopt a Shannon-divergence minimization criterion function, in which the density is evaluated at each point in the data to preserve the original distribution support. We define the following density estimators:

Utilizing the smoothing Jensen-Shannon Divergence emanating from (ref) sinn2018non:

where $\bm{\upsigma}:=(\upsigma_x,\upsigma_y,\upsigma_z)$ is the bandwidth vector. The optimal generator's strategy is to minimize the divergence:

We attend to the generation of simulated data for the purpose of causal counterfactual inference by using the estimated structural functions. The synthetic data is constructed from the synthetic (estimated) structural functions $\widetilde{g}$ and $\widetilde{h}$, whereas the real data is constructed from the original structural model $g$ and $h$ in (ref) and in (ref), respectively. The latter will be used as a benchmark for verifying the accuracy of the estimate obtained from the synthetic data.

Synthetic counterfactual machine vs. partial means estimator newey1994kernel

The following equations map an input vector $(Z_i,\upomega_{i},\widetilde{\upomega}_{i},\nu_i)$ to an output vector $(X_{\text{real,i}}^{\text{CF}}, Y_{\text{real,i}}^{\text{CF}})$ in order to generate the real counterfactual (CF) data set from the structural functions $g$ and $h$:

The following equations map an input vector $(Z_i,\upomega_{i},\widetilde{\upomega}_{i},\nu_i)$ to an output vector $(X_{\text{synt,i}}^{\text{CF}}, Y_{\text{synt,i}}^{\text{CF}})$ in order to generate the synthetic (synt) counterfactual (CF) data set from the synthetic (estimated) structural functions $\widetilde{g}$ and $\widetilde{h}$:

Similarly, the following equations map an input vector $(Z_i,\upomega_{i},\nu_i)$ to an output vector $(X_{\text{real,i}}^{\text{O}}, Y_{\text{real,i}}^{\text{O}})$ in order to generate the real observed (O) data set from the structural functions $g$ and $h$:

and the following equations map an input vector $(Z_i,\upomega_{i},\nu_i)$ to an output vector $(X_{\text{synt,i}}^{\text{O}}, Y_{\text{synt,i}}^{\text{O}})$ in order to generate the synthetic (synt) observed (O) data set from the synthetic (estimated) structural functions $\widetilde{g}$ and $\widetilde{h}$:

The estimated conditional expected outcomes given the action in the real data is:

The estimated conditional expected outcomes given the action in the synthetic (synt) data is:

The estimated counterfactual expected outcomes given the action in the real data is:

The estimated counterfactual expected outcomes given the action in the synthetic (synt) data is:

The partial means estimator newey1994kernel is presented here as a benchmark, as it is applied to the observed data set, unlike the counterfactual expression applied to the counterfactual dataset in (ref):

with

We note that $\widehat{\widehat{\mathbb{E}}}\left[Y_{\text{real}}|\text{do}(X_{\text{real}}=x)\right]$ and $\widehat{\mathbb{E}}\left[Y_{\text{real}}|\text{do}(X_{\text{real}}=x)\right]$ are different empirical objects. In the former $Y_{\text{real}}^{\text{O}}$, available only in theory (as $\upeta$ is unknown), is generated from the same realization of $\upomega$ generated $X_{\text{real}}^{\text{O}}$, while in the latter $Y_{\text{real}}^{\text{CF}}$ is generated by $\widetilde{\upomega}$ (another realization of $\upomega$) which is independent of $X_{\text{real}}^{\text{O}}$. Consequently, $\widehat{\mathbb{E}}\left[Y_{\text{real}}|\text{do}(X_{\text{real}}=x)\right]$ is expected to be more accurate by not being limited to the the joint support of $\upomega$ and $X_{\text{real}}^{\text{O}}$.

Equation (ref) utilizes $\upomega:=F_{\upeta}(\upeta)$ which is unobserved in the data. In case of strict monotonicity the following holds $F_{\upeta}(\upeta)=F_{X|Z}(X|Z)$ imbens2009identification. Thus, one can construct $V_{\text{real,i}}:=\widehat{F}_{X|Z}(X_{\text{real,i}}^{\text{O}}|Z_i)$ to estimate $\upeta$ (up to normalization) and employ the partial means estimator as follows:

This is the benchmark used for comparison of our results. The reasons for choosing equation (ref) for comparison are: (i) to show that our counterfactual goes beyond conditional quantiles of observables captured by the transformation above (violation of monotonicity) and further (ii) to demonstrate our approach performance in the presence of restricted support which is an acute problem in empirical research regardless of the presence or absence of monotonicity chernozhukov2020.

Simulations

In what follows we verify our proposed estimator's performance by employing monte-carlo simulations. The section has two rules: the first is to verify our model performance when the underlying model is characterized by intrinsically non-monotonic structural equations for $X$ and $Y$. The second is to use various examples representing important and widely used economic models of supply and demand which lie in the heart of the economic science, such as: the “backward-bending” supply curve of labor Hanoch1965; the Almost Ideal Demand System (AIDS) deaton1980almost; the Constant Elasticity of Substitution (CES) production function Sato1967 as well as the Translog utility function christensen1975, which may also affect the demand and supply curves. A specific model of demand and supply functions have been discussed in blundell2014control under the assumption of strict monotonicity in the unobserved non-additively-separable error terms. We however, attend to a more general (and non-monotonic) functional form for various economic phenomena.

The various widely used economic models discussed above are estimated from data sets generated by the following data generation processes:

The error terms are generated as follows: $\upeta, \nu\sim N(0,1)$ such that $\upeta\protect\mathpalette{\protect\independenT}{\perp}\nu$,

We consider various specifications for the structural equations of $X$ and $Y$, each one is chosen from the following functions: the Trans-log utility functions in (ref); the Almost Ideal Demand System (AIDS) functions in (ref); the Constant Elasticity of Substitution (CES) production functions in (ref); the backward-bending supply curves in (ref); a Cobb-Douglas production function in (ref) and the hyperbolic tangent model in (ref). For sake of generality, for each of these model specifications, we verify our model performance under different Neural network architectures, such that for each given architecture we generated $100$ data sets. These data sets are used to calculate the means and the standard error of our estimates. Such depth of simulations can be carried out only by parallel computing that can handle very demanding computations at a reasonable time in order to explore many specifications and variations of data generation processes.

Each of tables (ref)-(ref) (appendix (ref)) summarizes the results obtained from Monte-Carlo simulations generating $100$ data sets, each consisting of $10,000$ observations from specified structural model (data generation process) includes the following columns: a specific quantile of the action variable (column ($a$)); the mean value (over data sets) of the action variable in this quantile (column ($b$)); the mean (over data sets) of the estimated expected counterfactual outcome (given the action variable) using the real and synthetic data in (column ($c$)) and (column $(d)$), respectively; the mean (over data sets) of the estimated expected outcome (given the action variable) using the real data in (column $(e)$) and (column $(f)$), respectively; the control variable results imbens2009identification represent the counterfactual expected outcome that would be obtained by evaluating the latent context variable through the imposition of strict monotonicity (column $(g)$); and the partial means results newey1994kernel represent the counterfactual expected outcome that would be obtained if the context were known (column $(h)$);

We quantify the proximity of the synthetic data set generated by the synthetic generator functions $\widetilde{g}$ and $\widetilde{h}$ to the real data set generated by the “true” generator functions, $g$ and $h$. For tractability, our analysis disentangles the restricted support and divergence from strict monotonicity issues by using the following four criteria:

In the theoretical model we have shown that the importance of allowing for measure-preserving transformations when only distributional similarity is required. That is, we are interested in measuring the similarity between the empirical distributions of two sequences of estimates collected from the Monte Carlo simulations. Each sequence belongs to a different estimator. Thus, in the empirical application we opted for the Sliced-Wasserstein Distance (SWD) manole2019minimax, because it is designated for the comparison of two empirical distributions in a non-parametric manner. For robustness, these distributions must be compared over quantiles in the interior of the distribution to avoid the imposition of parametric assumptions regarding the tail behaviour of each distribution. Thus, we specify a scalar $\Lambda\in(0,0.5)$ for which $[\Lambda,1-\Lambda]$ denotes an interval of quantiles in the interior of the distribution. Using the aforementioned (SWD) methodology we append an asterisk and a dagger to indicate sequences of estimates with similar empirical distributions (over the Monte Carlo simulations) to those of the real counterfactual and the real observed results, respectively. Consequently, we perform these similarity tests by comparing distributions using only quantiles in the interval $[0.05,0.95]$ at significance level of $95\%$.

Summary statistics of results are presented in tables (ref)-(ref) (appendix (ref)). Entries represent average (standard error) of estimates obtained from $100$ samples of $10,000$ observations each. The overall picture emerging from the simulations, as appears in the tables, points to the superiority of our proposed casual inference. This is apparent from the similarity (in the SWD sense) between the real (column $(e)$) and synthetic (column ($f$)) counterfactual results, which appears in the various tables representing different known economics models. In particular, when juxtaposed on the much different and inferior results of existing models as shown in columns ($g$) and ($h$). Basically, existing causality models hardly successfully reconstruct the real counterfactual results. Further, even when choosing a monotonic structural form like the CES in table (ref) (appendix (ref)), results when applying conventional casual inference (column ($g$) and ($h$)) are still inaccurate due to an improper restricted support domain. As can be be seen, our model produces real and synthetic counterfactual results (in bold) which obey SWD similarity. For instance, in the 80th quantile of $X$, we report an estimate of (standard error) (table (ref), column ($e$)) 13.94 (0.344) for the real counterfactual and (table (ref), column ($f$)) 14.13 (0.553) for the synthetic counterfactual. The indication that the empirical support domain is restricted is based on the fact that the real counterfactual estimates are SWD different from the partial means estimates (as described in criterion (ref)). Just as an example, in the 90th quantile of $X$, we report an estimate of (table (ref), column ($e$)) 15.73 (0.728) for the real counterfactual and (table (ref), column $(h)$) 20.43 (0.960) for the partial means. The high similarity between the average estimates of the control variable (column $(g)$) and the partial means (column $(h)$) indicates that indeed there is no statistical evidence for the violation of strict monotonicity (as described in criterion (ref)). For instance, in the 5th quantile of $X$, we get the estimates (standard error) of (table (ref), column $(g)$) 6.98 (0.338) for the control variable and (table (ref), column $(h)$) 6.96 (0.276) for the partial means; in the 25th quantile of $X$, the estimate is (table (ref), column $(g)$) 8.52 (0.213) for the control variable and (table (ref), column $(h)$) 8.51 (0.188) for the partial means; in the 75th quantile of $X$, it is (table (ref), column $(g)$) 15.58 (0.338) for the control variable and (table (ref), column $(h)$) 15.59 (0.276) for the partial means.

In table (ref) (appendix (ref)) we apply the simulations to the known Translog function which is non-monotonic. The synthetic counterfactual results are not statistically different than the real counterfactual results. For instance, in the 25th quantile of $X$, we report an estimate of (table (ref), column ($e$)) 16.63 (0.175) for the real counterfactual and (table (ref), column ($f$)) 16.61 (0.407) for the synthetic counterfactual. However, the common practice results are extremely inaccurate in magnitudes (columns $(g)$ and $(h)$). In fact, the non-monotonicity is expressed by the disparity between the control function and the partial means estimates. Just as an example, in the 20th quantile of $X$, the estimate is (table (ref), column $(g)$) 20.95 (1.087) for the control variable and (table (ref), column $(h)$) 14.71 (0.857) for the partial means. In this specific quantile there is also a disparity between the empirical and restricted supports, captured by the gap between real counterfactual and the partial means estimates. For instance, in the 15th quantile of $X$, it is (table (ref), column ($e$)) 15.88 (0.190) for the real counterfactual and (table (ref), column $(h)$) 13.53 (0.959) for the partial means.

The entries in table (ref) (appendix (ref)) the estimates of the hyperbolic tangent form where the strict monotonicity is violated. Applying the our offered model produces very accurate results and exhibits insignificant statistical difference between the real and synthetic interventional distributions (columns $(e)$ and $(f)$). However, the common practice results are extremely inaccurate in both magnitudes and signs (columns $(g)$ and $(h)$). The violation of monotonicity is expressed by a strong disparity between the control function and the partial means estimates. For instance, the 20th quantile of $X$, we report an estimate of (table (ref), column $(g)$) -0.18 (0.291) for the control variable and (table (ref), column $(h)$) -0.56 (0.189) for the partial means. There is also a disparity between the empirical and restricted supports, captured by the gap between real counterfactual and the partial means estimates. E.g., in the 15th quantile of $X$, we report an estimate of (table (ref), column ($e$)) -0.17 (0.038) for the real counterfactual and (table (ref), column $(h)$) -0.70 (0.212) for the partial means.

The entries in table (ref) (appendix (ref)) exhibit estimates of the well-known demand function (AIDS). This very well-known functional form is intrinsically non-monotonic. Applying the our offered model produces very accurate results and exhibits insignificant statistical difference between the real and synthetic interventional distributions (columns $(e)$ and $(f)$). For instance, in the 65th quantile of $X$, we report an estimate of (table (ref), column ($e$)) 16.71 (0.188) for the real counterfactual and (table (ref), column ($f$)) 16.65 (0.348) for the synthetic counterfactual. However, as can be seen results largely demonstrate the failure to mimic the real counterfactual distribution when applying conventional casual inference (column ($g$) and ($h$)). For instance, in the 90th quantile of $X$, we report an estimate of (table (ref), column ($e$)) 22.79 (0.669) for the real counterfactual and (table (ref), column $(h)$) 26.58 (0.975) for the partial means; For this quantile, the synthetic counterfactual result (table (ref), column ($f$)) is 21.66 (0.758). For instance, in the 55th quantile of $X$, we report an estimate of (table (ref), column $(g)$) 21.16 (1.101) for the control variable and (table (ref), column $(h)$) 16.67 (0.748) for the partial means; For this quantile, the real counterfactual result is (table (ref), column ($e$)) 16.06 (0.195) and the synthetic counterfactual result is (table (ref), column ($f$)) 16.09 (0.344). There is also a disparity between the empirical and restricted supports, captured by the gap between real counterfactual and the partial means estimates.

The entries in table (ref) (appendix (ref)) represent those emerging from the estimation of backward-bending supply curve. In this instance, monotonicity is intrinsically violated. This is attested to the disparity between the control function and the partial means estimates in both magnitude as well as in sign, which in reality may lead to wrong public policy. For instance, in the 35th quantile of $X$, we report an estimate of (table (ref), column $(g)$) -0.33 (0.135) for the control variable and (table (ref), column $(h)$) 0.33 (0.091) for the partial means. Applying our offered model produces very accurate results and exhibits insignificant statistical difference between the real and synthetic interventional distributions (columns $(e)$ and $(f)$). Just as an example, in the 60th quantile of $X$, the real counterfactual result (table (ref), column $(e)$) is 0.61 (0.040) and the synthetic counterfactual result (table (ref), column $(f)$) is 0.59 (0.093). Further, below the 40’th quantile it can be seen that there is a disparity between the empirical and restricted supports, captured by the gap between real counterfactual and the partial means estimates. In another instance, in the 25th quantile of $X$, we report an estimate of (table (ref), column ($e$)) 0.47 (0.033) for the real counterfactual and (table (ref), column $(h)$) 0.21 (0.087) for the partial means. The synthetic counterfactual result is (table (ref), column ($f$)) 0.44 (0.109) indicating that it is not statistically different than the real result.

To sum up, the observed and counterfactual similarity criteria (ref) and (ref) tables (ref)-(ref) (appendix (ref)) developed in this paper demonstrate the high accuracy of the synthetic counterfactual machinery relative to the real data counterfactual result. In terms of sensitivity to the violation of the strict monotonicity assumption, as captured by criterion (ref), our proposed estimator produces accurate results regardless of the absence or presence of monotonicity and outperforms control variable existing approach. In terms of sensitivity to the violation of unrestricted support domain, as captured by criterion (ref), our proposed estimator outperforms both the partial means as well as the control variable estimators, by mimicking the results obtained for the real data counterfactual outcome. It emphasized the results generated by the aforementioned existing models produce results that are inaccurate both in magnitude and sign of the estimates.\footnote{In the supplementary material we provide different stability measures.}

Conclusion

Different latent contexts may affect the relations between cause and outcome and hence can have significant ramification for causal relations and inference. We develop a novel identification strategy as well as a new estimator for triangular models in the presence of non-separable disturbances, which unlike the common practice, does not rely on the strict monotonicity assumption. The key result of this identifiability approach is the explicit characterization of the distributional relationship between the latent context variable and the vector of observables through Fredholm integral equations governed by generator functions and an unknown kernel function, inducing a non-monotonic inverse problem. The role of these generator functions is two-fold: (i) to characterize the unknown kernel function induced by these generator functions and (ii) to ensure that the estimator of the unknown quantity is a continuous function of the data given this unknown kernel function. This very formulation facilitates the establishment of uniqueness of the interventional distribution given the observables. Furthermore, we develop a novel CONGAN estimator based on a feed-forward Neural network architecture generating a synthetic counterfactual distribution. This synthetic distribution represents various combinations of actions, outcomes and contexts, rarely available in finite samples. In the simulations, the proposed estimator's performance in finite samples has been validated in several aspects by testing various data generation processes each associated with the commonly employed structural functions: Cobb-Douglas, AIDS, CES, Translog and backward-bending supply models. It can be seen that our model generates synthetic data mimicking the behavior in the real data (in terms of SWD similarity). This comparison has been done for different quantiles of the action variable. Our results are also compared to both the monotonic control variable as well as partial means estimators. In some of these models, the counterfactual results obtained by the conventionally used partial means estimator, largely deviate from the benchmark real data results in quantity and sign. This may in reality have important ramifications as to the proper public policy.

In future work, our framework can be incorporated in dynamic models as well as panel models to capture the dependence of the action and outcome on past events through a different Neural network architecture, such as recurrent Neural networks.