The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
231,019 characters
Context-dependent Causality (the Non-Monotonic Case)
\maketitle
\begin{abstract}
We develop a novel identification strategy as well as a new estimator for context-dependent causal inference in non-parametric triangular models with non-separable disturbances. Departing from the common practice, our analysis does not rely on the strict monotonicity assumption.
Our key contribution lies in leveraging on \textbf{\textit{diffusion models}} to formulate the structural equations as a system evolving from noise accumulation to account for the influence of the latent context (confounder) variable on the outcome.
Our identifiability strategy involves a system of Fredholm integral equations expressing the distributional relationship between a latent context variable and a vector of observables.
These integral equations involve an unknown kernel and are governed by a set of structural form functions, inducing a \textbf{\textit{non-monotonic inverse problem}}.
We prove that if the kernel density can be represented as an infinite mixture of Gaussians, then there exists a unique solution for the unknown function. 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.
On the methodological front we leverage on a novel and enriched \textit{Contaminated} Generative Adversarial (Neural) Networks (\textbf{CONGAN}) which we provide as a solution to the non-monotonic inverse problem.
\end{abstract}
\vspace{.5em}
{\textit{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
\vspace{20em}
\pagebreak
\section{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 \textbf{\textit{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 \textit{unknown kernel function}, inducing a \textbf{\textit{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 \citep{chesher2003identification}, strict monotonicity in the case of average local treatment effects \citep{angrist1994identification} as well as in the case of triangular models \citep{imbens2009identification,hoderlein2017corrigendum}. Alternative models assume restricted support domain (locality) \citep{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 \citep{newey1994kernel}, interventions \citep{pearl2019interpretation} and local average treatment effect \citep{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., \cite{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 \citep{hoderlein2016testing, hoderlein2007identification}, multiple equilibria in strategic behavior \citep{myerson1999nash}, search behavior \citep{chetverikov2019testing, agarwal2020searching} as well as in prospect theory \citep{tversky1979analysis}.\footnote{A recent paper \citep{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 \citep{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 \textit{\textbf{Contaminated}} version we develop and apply to the Generative Adversarial Networks model (\textbf{CONGAN}) \citep{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{Sec:Discussion} gives an overview of the common practice of causal inference in triangular models. Section \ref{Sec:Theoretical:model} presents the diffusion problem formulation and its underlying assumptions.
Section \ref{Sec:alternative_representation} shows an alternative formulation of the triangular model as an integral equation to simplify the analysis. Section \ref{sec:example} provides an illustrative example for identifiability which does not rely on monotonicity assumption. Section \ref{sec:completeness} characterizes strong completeness of sufficient statistic as a building block for identifiability.
Section \ref{Sec:Identification} presents the identifiability by relying on strong completeness. The estimator is presented in section \ref{Sec:Estimation}.
In section \ref{Sec:Simulations} Monte-Carlo simulations are used to verify our theoretical model performance. Section \ref{Sec:Conclusion} concludes.
\section {Discussion and literature review}\label{Sec:Discussion}
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'' \citep{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 \citep{chernozhukov2007instrumental, kim2020partial}. In other cases, the monotonicity requirement is imposed on the instrumental variable \citep{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 \citep{athey2016recursive} or a conditional independence on observables \citep{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 \citep{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 \citep{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 \citep{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 \citep{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 \citep{athey2020}. Recent work criticizes this penalized Wasserstein GANS due to disregarding significant parts of the support domain in the data \citep{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.
\section{Diffusion-based non-monotonic triangular models}
\label{Sec:Theoretical:model}
Consider the following triangular simultaneous equations model, pioneered by \citet{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,
\begin{align}
& X = h(z,\eta) \sim \sum_{i=1}^{L} \alpha_{1,i}\mathcal{N}\left(\mu_{1,i}(z) + \rho_{12}\sigma_{1,i}(z)\left( \frac{\eta-\mu_{\eta}}{\sigma_\eta}\right), \sigma_{1,i}^2(z)(1-\rho_{i,12}^2)\right) \label{X:Structural}
\end{align}
and $Y = g(X,\upepsilon)$ satisfying,
\begin{align}
& g(X,\upepsilon)|\upeta=\eta\sim \sum_{i=1}^{L} \alpha_{2,i}\mathcal{N}\left(\mu_{2,i}(X) + \rho_{i,13}(X)\sigma_{2,i}(X)\left( \frac{\eta-\mu_{\eta}}{\sigma_\eta}\right), \sigma_{2,i}^2(X)(1-\rho_{i,13}^2(X))\right)\label{Y:Structural}.
\end{align}
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. \eqref{X:Structural} and \eqref{Y:Structural},
\begin{align}
& X = \sum_{i=1}^{L}\alpha_{1,i}\left\{\mu_{1,i}(Z) + \sigma_{1,i}(Z)\left(\sqrt{1-\uptheta_i} \Phi^{-1}(\bm\nu_{i,x}) + \sqrt{\uptheta_i}\left( \frac{\upeta-\mu_{\eta}}{\sigma_\eta}\right)\right)\right\},\nonumber\\
& Y = \sum_{i=1}^{L}\alpha_{2,i}\left\{\mu_{2,i}(X) + \sigma_{2,i}(X)\left(\sqrt{1-\upbeta_i(X)} \Phi^{-1}(\bm\nu_i) + \sqrt{\upbeta_i(X)}\left( \frac{\upeta-\mu_{\eta}}{\sigma_\eta}\right)\right)\right\},\nonumber
\end{align}
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)$,
\begin{align}
Y = \sum_{\substack{i,j=1}}^{L}\alpha_{i,j}^*\left\{\mu_{2,i}(X) + \sigma_{2,i}(X)\left(\sqrt{1-\upbeta_i(X)} \Phi^{-1}(\bm\nu_i) + \sqrt{\upbeta_i(X)}\left( \frac{\Psi_{j}(X,Z)+\bm\epsilon_0-\mu_{\eta}}{\sigma_\eta}\right)\right)\right\}, \label{Y:ReducedForm}
\end{align}
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. \eqref{Y:ReducedForm} is a reduced-form representation of $Y$ without any endogenous predictor variables. This formulation shares similarities with diffusion models \citep{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}$,
\begin{definition}[Support]
$\text{Supp}(Q):=\text{The closure of the set of points for which } f_{Q}(q)\ne 0$.
\end{definition}
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 \cite{{imbens2009identification}}
the functions $h$ and $g$ may be nonseparable in their respective disturbances $\upeta$ and $\bm{\epsilon}$.
The model (\ref{X:Structural})-(\ref{Y:Structural}) implies that $\upeta$ is a \textit{confounder} or common cause, affecting both $X$ directly via \eqref{X:Structural},
and $Y$ indirectly through the dependence of
$\bm{\epsilon}$ on $\upeta$.
Finally, following previous work \citep{imbens2009identification}, we assume the following standard condition holds:
\begin{asu}[Common support]\label{asu:common}
The conditional random variable $\upeta|X=x$ has a \textbf{fixed support} independent of $x$.
Namely,
$\text{Supp}(\upeta|X=x) = \text{Supp}(\upeta),$ $\forall x$.
\end{asu}
For future use, we present two important properties of the above triangular model.
First, note that Eq. (\ref{X:Structural}) 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{Appendix:Proof}).
\begin{lemma}\label{lemma:eps:support}
Under assumption \ref{asu:common} it holds that,
\begin{align}
\text{Supp}(\left.\bm\epsilon\right|X=x)=\text{Supp}(\bm\epsilon) \quad \forall x\in\text{Supp}(X).\label{eq:Eps:Support}
\end{align}
\end{lemma}
Given $n$
i.i.d. triplets $(x_i,y_i,z_i)$ from the above triangular model, Eqs. (\ref{X:Structural})-(\ref{Y:Structural}), 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 \citep{pearl2019interpretation}, the above quantity of interest is known as \textit{intervention}. It is given by
\begin{align}
T(y|x) = P\left(Y\le y|\text{do}(X=x)\right):= \mathbb{E}_{\bm\epsilon\sim F_{\bm\epsilon}}\left[\mathds{1}\left\{g(x,\bm\epsilon)\le y\right\}\right].
\label{eq:do_X}
\end{align}
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)$ \citep{imbens2007nonadditive}.
The reason is that even if $g$ were known, under (\ref{X:Structural})-(\ref{Y:Structural}) the random variables $X$ and $\bm\epsilon$ might be jointly dependent.
In contrast, in Eq. (\ref{eq:do_X}) the expectation is over the marginal distribution of $\bm\epsilon$, effectively treating $\bm\epsilon$ as independent of $X$.
\begin{comment}
\textcolor{red}{
1.
\\
the joint distribution of $(X,Y,Z)$ does not uniquely determine the quantity $T(y|x)$.
\\
\\
3. Most works imposed the following two standard conditions, which we shall also assume hold in our work. Describe assumptions 3.1 3.2,
\\
4. These two assumption are still not sufficient for identifiability. EXAMPLE / CITATION ?
\\
5. previous works additional assumptions sufficient for identifiability. $H_c$, control variable...
\\
6.
}
\end{comment}
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., \cite{blundell2014control,imbens2009identification}.
Specifically, \citet{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,
\begin{align}
H_c(y|x) = \int_{0}^1 F_{Y|X=x,V=v}(y)dv
\label{eq:H_c}.
\end{align}
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 \cite{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. }
\begin{asu}[Continuity of $X|Z$ and $\bm{\epsilon}|\upeta$]
\begin{subasu}\label{asu:continuity}
For any $z\in \text{Supp}(Z)$, the random variable $X|Z=z$ is continuous. In particular, the set of values of $\upeta$ where $\partial h(z,\eta)/\partial \eta=0$ is of measure zero.
\end{subasu}
\begin{subasu}\label{asu:cond:continuity}
The random vector $\bm{\epsilon}|\upeta$ is a $d$-dimensional continuous random vector for any $\upeta$.
\end{subasu}
\end{asu}
\begin{asu}[Completeness]
The distribution of the random vector $(X,\upeta)|Z$ belongs to a complete family of distributions with respect to $Z$.
\end{asu}\label{asu:bounded}
\begin{asu}[Square integrability]\label{asu:Integrability}
The density $f(y|X=x,Z=z)$ is square integrable, and is known.
\end{asu}
\begin{comment}
\begin{asu}
For any function $u(\upeta,x)$ with finite expectation with respect to $\upeta$, if
\begin{equation}
\label{eq:E_w_eta_x}
\mathbb{E}_{\upeta|X=x,Z=z}\left[u(\upeta,x)
\right]=0 \quad \forall (x,z)\in\text{Supp}(X,Z),
\end{equation}
then $\mathbb{E}_{\upeta}\left[u(\upeta,x)\right]=0$ for all $x\in \text{Supp}(X)$.
\end{asu}
\end{comment}
{
The first key contribution of our work is to show that under the
triangular model with assumptions \ref{asu:common}-\ref{asu:Integrability},
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 \citep{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 \eqref{X:Structural}-\eqref{Y:Structural}. 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{Sec:integral}, we derive an
integral equation relating the intervention
of Eq. (\ref{eq:do_X}) 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{Sec:Identification} we prove that {\em any}
solution of this integral equation gives the same intervention, hence proving identifiability of the intervention.
\section{An alternative representation for the intervention}\label{Sec:alternative_representation}
A key preliminary step is to derive an {\em alternative} representation
for the intervention $T(y|x)$.
{Consider the following cumulative distribution function (CDF),}
\begin{equation}
H(y|x) = \mathbb{E}_{\upeta}[ F_{Y|X=x,\upeta}(y)] =
\int F_{Y|X=x,\upeta=\eta}(y)f_{\upeta}(\eta)d\eta.
\label{eq:def_H}
\end{equation}
Note that the function inside the integral , with $\upeta$ attaining all possible values in its support, is well defined by assumption \ref{asu:common}.
The following lemma is key to our proposed approach, both for proving identifiability and for inference.
\begin{lemma}[Functional equivalence]
\label{lem:H_T}
Under the triangular model with assumption \ref{asu:common},
the function $H(y|x)$ defined in Eq. \eqref{eq:def_H} satisfies that
$\forall\hspace{0.2em}(x,y)$ in the relevant support,
\begin{eqnarray}
H(y|x) = T(y|x).
\label{eq:H_T}
\end{eqnarray}
\end{lemma}
\begin{proof}
Using the definition of $Y$ in Eq. \eqref{Y:Structural} the function
$F_{Y|X=x,\upeta}(y)$ may be equivalently written as follows
\begin{equation}
F_{Y|X=x,\upeta}(y) = \mathbb{E}_{\bm{\epsilon}|X=x,\upeta}\left[\mathds{1}\left\{Y\le y\right\}|X=x,\upeta\right] = \mathbb{E}_{\bm{\epsilon}|X=x,\upeta}\left[
\mathds{1}\left\{g(x,\bm{\epsilon})<y\right\}\right].
\end{equation}
By conditional independence $X\protect\mathpalette{\protect\independenT}{\perp}\epsilon|\upeta$, the random vector $\bm{\epsilon}|X=x,\upeta$ has the same distribution as
$\bm{\epsilon}|\upeta$. Inserting the resulting expression back into Eq. \eqref{eq:def_H}
gives that
\begin{equation}
H(y|x) = \mathbb{E}_{\upeta\sim F_{\upeta}}\left[\mathbb{E}_{\bm{\epsilon}|\upeta}\left[\mathds{1}\left\{g(x,\bm{\epsilon})\le y\right\}\right]\right].
\end{equation}
By the law of total expectation, the RHS is simply
$\mathbb{E}_{\bm \epsilon}[\mathds{1}\left\{g(x,\bm\epsilon)\le y\right\}]$, which
by Eq. \eqref{eq:do_X} is $T(y|x)$.
\end{proof}
By Lemma \ref{lem:H_T}, 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{eq:def_H}) is known as the \textit{partial means} estimator \citep{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.
\section{An Example for a Triangular Data Generation Process}\label{sec:example}
Let's examine the non-monotonic case, $h(z,\eta)\ne F_{X|Z=z}^{-1}(F_{\upeta}(\eta))$.
\begin{align}
& \bm\Sigma_{y,\eta}(X,Z):=\left[\begin{matrix}
\sigma_y^2(X) & \sigma_2(x)\sigma_{\eta}(X,Z)\rho_{13}(X) \\
\sigma_2(x)\sigma_{\eta}(X,Z)\rho_{13}(X) & \sigma_{\eta}^2(X,Z)
\end{matrix}\right]\nonumber
\end{align}
The coefficient $\rho_{13}(x)$ determines the conditional correlation between $Y$ and $\upeta$ given $X=x, Z=z$.
\begin{align}
& Z \sim F_Z, \quad \upeta\sim N(\mu_{\eta}, \sigma_{\eta}^2) \nonumber \\
& X|Z=z, \upeta=\eta\sim N\left(\mu_1(z) + \rho_{12}\sigma_1(z)\left( \frac{\eta-\mu_{\eta}}{\sigma_\eta}\right), \sigma_x^2(z)(1-\rho_{12}^2)\right)\nonumber\\
& Y|X=x,Z=z\sim N\left(\mu_2(x) + \rho_{13}(x)\sigma_2(x)\left( \frac{\eta(x,z)-\mu_{\eta}}{\sigma_\eta}\right), \sigma_y^2(x)(1-\rho_{13}^2(x))\right)\nonumber
\end{align}
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,
\begin{align}
& X(z,\eta, \bm\varepsilon_x) = \mu_1(z) - \rho_{12}\sigma_1(z)\frac{\mu_{\eta}}{\sigma_\eta} + \sigma_1(z)\rho_{12}\frac{\upeta}{\sigma_\eta} + \sigma_1(z)\sqrt{1-\rho_{12}^2} \Phi^{-1}(\bm\varepsilon_x) \nonumber\\
& Y(x,z,\bm\nu) = \mu_2(x) - \rho_{13}(x)\sigma_2(x)\frac{\mu_{\eta}}{\sigma_{\eta}} + \sigma_2(x)\rho_{13}(x)\frac{\upeta(z,x)}{\sigma_{\eta}} + \sigma_2(x)\sqrt{1-\rho_{13}^2(x)} \Phi^{-1}(\bm\nu)\label{Diffused_outcome}
\end{align}
\begin{align}
&\upeta\frac{\sigma_1(z)\rho_{12}}{\sigma_\eta} = x - \mu_1(z) + \rho_{12}\sigma_1(z)\frac{\mu_{\eta}}{\sigma_\eta} - \sigma_1(z)\sqrt{1-\rho_{12}^2} \Phi^{-1}(\bm\varepsilon_x) \nonumber
\end{align}
\begin{align}
&\upeta =\mu_{\eta} + \frac{\sigma_\eta}{\sigma_1(z)\rho_{12}}(x - \mu_1(z)) - \frac{\sigma_{\eta}\sqrt{1-\rho_{12}^2}}{\rho_{12}} \Phi^{-1}(\bm\varepsilon_x) \nonumber
\end{align}
\begin{align}
& \sigma_{\eta}(x,z):=\frac{\sigma_{\eta}\sqrt{1-\rho_{12}^2}}{\rho_{12}} \nonumber \\
& \mu_{\eta}(x,z):=\mu_{\eta} + \frac{\sigma_{\eta}}{\sigma_1(z) \rho_{12}}\left(x-\mu_1(z)\right) \quad \mu_{\eta}=0 \quad (\text{a normalization}) \nonumber \\
& \upeta(z,x) = \mu_{\eta}(x,z)- \sigma_{\eta}(x,z)\Phi^{-1}(\bm\varepsilon_x). \nonumber
\end{align}
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.
\begin{align}
& \mathbb{E}_{\upeta}\left[\mathbb{COV}\left[X,Z\right]|\upeta\right] = \mathbb{COV}\left[\mu_1(z), Z\right] + \sqrt{1-\rho_{12}^2} \mathbb{COV}\left[\sigma_1(z) \Phi^{-1}(\bm\varepsilon_x), Z\right] \nonumber
\end{align}
This gives,
\begin{align}
& \sqrt{1-\rho_{12}^2}=\frac{\mathbb{COV}\left[X,Z\right]-\mathbb{COV}\left[\mu_1(z), Z\right] }{\mathbb{COV}\left[\sigma_1(z) \Phi^{-1}(\bm\varepsilon_x), Z\right]} \nonumber
\end{align}
Given $\rho_{12}$, the coefficient $\rho_{13}(x)$ can be identified by a linear regression of the form,
\begin{align}
& Y(x,z,\bm\nu) = \overbrace{\mu_2(x)}^{\mathbb{E}[Y|X=x]} - \rho_{13}(x)\left[ \frac{\sigma_2(x)}{\sigma_1(z)\rho_{12}}\overbrace{(x-\mu_1(z))}^{x-\mathbb{E}[X|Z=z]}\right] + \epsilon \nonumber
\end{align}
The counterfactual outcome variable
\begin{align}
& Y^{\text{CF}}(\text{Do}(X=x)) = \mu_2(x) + \rho_{13}(x)\sigma_2(x)\left( \frac{\eta-\mu_{\eta}}{\sigma_\eta}\right) + \sigma_2(x)\sqrt{1-\rho_{13}^2(x)} \Phi^{-1}(\bm\nu) \label{Diffused_CF_outcome}
\end{align}
Define $\upbeta:=\rho_{13}^2(x)$. For any given pair $(x,z)$, $Y$ can be reconstructed from a noise propagated by $\bm\nu$,
\begin{align}
& Y(x,z) = \mu_2(x) + \sqrt{\upbeta}\sigma_2(x)\left( \frac{\eta(x,z)-\mu_{\eta}}{\sigma_\eta}\right) + \sigma_2(x)\sqrt{1-\upbeta} \Phi^{-1}(\bm\nu). \label{Diffused_CF_outcome:beta}
\end{align}
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 \citep{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.}
\begin{align}
\frac{Y(x,z)-\mu_2(x)}{\sigma_2(x)}:=\sqrt{1-\upbeta}Y_{0}+\sqrt{\upbeta}\left(\frac{\eta(x,z)-\mu_{\eta}}{\sigma_\eta}\right).\label{eq:Y_t}
\end{align}
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. \eqref{Diffused_CF_outcome} gives us the direct linkage between the identified parameters and the counterfactual output. Namely, by substituting $\eta(x,z)$ in Eq. \eqref{eq:Y_t} 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.
\section{Completeness of Sufficient Statistic}\label{sec:completeness}
\begin{definition}[Completeness]
Let $X$ be a random variable with density functions parameterized by $\theta$. A statistic $T(X)$ is said to be complete if, for every measurable function $g$, the following holds:
\[ \mathbb{E}[g(T(X))] = 0 \text{ for all } \theta \implies g(T(X)) = 0 \text{ almost everywhere.} \]
\end{definition}
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.
\textbf{Informal Definition (Strong completeness)}.
The concept of strong completeness applied to the relationship between regimes (subsets of a parameter space), introduced by \cite{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).
\begin{figure}[H]{Two Mutually Exclusive Regimes Nested in a Single Model}
\begin{center}
\begin{tikzpicture}[scale=1.5]
\draw[fill=blue, opacity=0.3] (0,0) circle (1);
\draw[fill=green, opacity=0.3] (2,0) circle (1);
\draw[fill=yellow, opacity=0.2] (1,0) circle (2);
\node[align=center, text=black] at (0,0) {Regime A \\ $\theta_2=\theta_2^{\prime}$ \\(\small randomized)};
\node[align=center, text=black] at (2,0) {Regime B \\ $\theta_2=\theta_2^{\prime\prime}$\\(\small observational) };
\node[align=center, text=black] at (1,-1.5) {A Two-regime Model\\ $\theta_1\in\Theta_1$};
\draw[->, thick] (1.5,-0.5) -- (0.5,-0.5);
\end{tikzpicture}
\end{center}
\end{figure}
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.
\begin{definition}[Strong Completeness]
A family of distributions $F:=\left\{ F_{Y|\theta_1, \theta_2} : \theta_1\in \Theta_1\right\}$ parameterized by $(\theta_1, \theta_2)$ of $d$-dimensional random vectors is said to be strongly complete with respect to $\theta_1$, if for every function $m:\mathbb{R}^d \to\mathbb{R}$ satisfying,
\[ \mathbb{E}_{Y\sim F_{Y|\theta_1, \theta_2}}\left[m(\bm Y)\right] = 0 \quad \forall\theta_1 \in \Theta_1, \]
we have that,
\[ \mathbb{E}_{ Y\sim F_{Y|\theta_1, \theta_2}}\left[\mathds{1}\left\{m(\bm Y)=0\right\}\right] = 1 \quad \forall(\theta_1,\theta_2) \in \Theta_1\times \Theta_2. \]
\end{definition}
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$).
\section{Identifiability}\label{Sec:Identification}
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{Lemma:Decoupling} and \ref{lemma:Invariant:Do} 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{Theorem:Identifiability}).
\subsection{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,
\begin{equation}
\Lambda_0 = \Lambda_0(\upomega):=F_{\upeta}^{-1}(\upomega),
\label{eq:Lambda_0}
\end{equation}
and recursively for $1\le j\le d$,
\begin{equation}
\Lambda_j = \Lambda_j(\upomega,\nu_1,...,\nu_j):=
F_{\epsilon_j|(\Lambda_0(\upomega),\Lambda_1(\upomega,\nu_1),\ldots,\Lambda_{j-1}(\upomega,\nu_1,\ldots,\nu_{j-1}))}^{-1}(\nu_j).
\label{eq:Lambda_j}
\end{equation}
The notations in Eqs. \eqref{eq:Lambda_0}-\eqref{eq:Lambda_j} are used to perform error-decoupling in the following lemma (see proof in appendix \ref{Appendix:Proof}):
\begin{lemma}
\label{Lemma:Decoupling}(Error-decoupling)
Let $\upomega\sim U[0,1]$ and $\bm{\nu}:=(\nu_1,\ldots,\nu_d) \sim U[0,1]^d$
be independent of $\upomega$.
Consider the following triangular model,
where $\Lambda_0$ and $\Lambda_j$ are defined in Eqs. (\ref{eq:Lambda_0}) and (\ref{eq:Lambda_j}),
\begin{eqnarray}
X^*&=&h^*(Z,\upomega):=h\left(Z,\Lambda_0(\upomega)\right)
\label{eq:tilde_X},\\
Y^*&=&g^*(X^*,\upomega,\bm{\nu}):=g\left(X^*,\Lambda_1(\upomega,\nu_1),\ldots,\Lambda_d(\upomega,\nu_1,\ldots,\nu_d)\right). \label{eq:tilde_Y}
\end{eqnarray}
Then, $(X^*, Y^*, Z) \overset{d}{=}(X,Y,Z)$.
\end{lemma}
There are two major differences
between the normalized triangular model Eqs. \eqref{eq:tilde_X}-\eqref{eq:tilde_Y}
and the original Eqs. \eqref{X:Structural}-\eqref{Y:Structural}.
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. \eqref{X:Structural}-\eqref{Y:Structural}. This challenge stems from the fact that $\upeta$ and $\bm \epsilon$ may in general be dependent, in an unknown fashion.
\begin{lemma}[Normalized Interventional Distribution]\label{lemma:Invariant:Do}
Suppose that assumption \ref{asu:cond:continuity} holds. Namely, $\bm{\epsilon}|\upeta$ is a $d$-dimensional continuous random vector for any $\upeta$.
In terms of the decoupled system as in Eqs. (\ref{eq:tilde_X})-(\ref{eq:tilde_Y}), the quantity of interest (\ref{eq:do_X})
admits the following form,
\begin{equation}
H(y|x) = \mathbb{E}_{\substack{\upomega\sim U[0,1]\\ \bm\nu\sim U[0,1]^d}}\left[ {\mathds{1}\left \{{g^*(x,\upomega,\bm\nu)\le y} \right\}} \right].
\label{eq:Tyx_nu}
\end{equation}
\end{lemma}
\subsection{Identifiability of the interventional distribution}
Recall that our primary goal is to estimate the intervention $T(y|x)$ of (\ref{eq:do_X}).
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
\begin{align}
&\widetilde{\mathcal{C}}:=\left\{(\widetilde{h},\widetilde{g}): \widetilde{X}=\widetilde{h}(Z,\widetilde{\upomega}), \, \widetilde{Y}=\widetilde{g}(\widetilde{X},\widetilde{\upomega},\bm\nu),\hspace{0.5em} \left(\widetilde{X},\widetilde{Y},Z\right)\overset{d}{=}\left(X^*,Y^*,Z\right), \, \widetilde{\upomega}\sim U[0,1], \bm\nu\sim U[0,1]^d \right\}.
\label{eq:equivalence:class}
\end{align}
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
\begin{equation}
\tilde{H}(y|x) = \mathbb{E}_{\substack{\tilde{\upomega}\sim U[0,1]\\ \bm\nu\sim U[0,1]^d}}\left[ {\mathds{1}\left \{{\tilde{g}(x,\tilde{\upomega},\bm\nu)\le y} \right\}} \right].
\label{eq:Tyx_nu}
\end{equation}
Note that any pair of structural functions $(\widetilde{h}, \widetilde{g})\in\mathcal{\widetilde{C}}$ constitutes a distribution,
\begin{align}
\mathcal{F}_{X,\upomega|Z=z}^{\widetilde{h}}(x,\omega):=\mathbb{E}\left[\mathds{1}\left\{\widetilde{h}(z,\upomega)<x, \upomega<\omega\right\}\right], \label{eq:F:X:omega}
\end{align}
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.
\begin{lemma}[Fourier Transform of Counterfactual Density]\label{lemma:Fourier:CF:representation}
Under assumptions \ref{asu:common}-\ref{asu:Integrability} the Fourier transform of the counter factual distribution admits the representation as a Fredholm integral equation of the first kind for any pair of frequencies $(\xi_1, \xi_2)$ and value of $z$,
\begin{align}
\Gamma_{\xi_1, \xi_2}(z;(\widetilde{h}, \widetilde{g}))=\mathbb{E}_{\substack{\upomega\sim U[0,1]\\X\sim F_{X|Z=z}}}\left[\int \exp(-i\xi_1 \cdot X) \exp(-i\xi_2 \cdot \widetilde{g}(X, \upomega, \bm\nu))d\bm\nu\right].\label{eq:box:counter:X:Y}
\end{align}
\end{lemma}
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,
\begin{align}
\Gamma_{\xi_1, \xi_2}^{\text{CF}}(z; (\widetilde{h}, \widetilde{g})) = \Gamma_{\xi_1, \xi_2}^{\text{CF}}(z; (h^*, g^*)) \quad \forall (z,\xi_1, \xi_2).\label{eq:coeff:obs:pair:CF:Fourier}
\end{align}
We rely on the following lemma as a building block in establishing that this condition is satisfied under assumptions \ref{asu:common}-\ref{asu:Integrability}, by formulating the observational equivalence condition in terms of Fourier transform.
\begin{lemma}[Fourier Transform Coefficients]\label{lemma:Fourier:representation}
Under assumptions \ref{asu:common}-\ref{asu:Integrability}
the following Fredholm integral equation of the first kind uniformly holds for any pair $(\widetilde{g}, \widetilde{h})\in\mathcal{\widetilde{C}}$ constituting a density $\mathcal{F}_{X,\upomega|Z=z}^{h^*}(x,\omega))$ as in Eq. \eqref{eq:F:X:omega},
\begin{align}
\Gamma_{\xi_1, \xi_2}(z)=
\mathbb{E}_{(X,\upomega)\sim \mathcal{F}_{X,\upomega|Z=z}^{\widetilde{h}}}\left[\int \exp(-i\xi_1 \cdot X) \exp(-i\xi_2 \cdot \widetilde{g}(X, \upomega, \bm\nu))d\bm\nu\right] \quad \forall (z, \xi_1, \xi_2) ,\label{eq:coeff:obs:pair:Fourier}
\end{align}
where the LHS is a specific Fourier transform coefficient of the known density $f_{X,Y|Z=z}$ for frequencies $(\xi_1, \xi_2)$.
\end{lemma}
The following theorem \ref{theorem:coefficients} states that for any two pairs of functions $(g^*,h^*)$ and $(\widetilde{g}, \widetilde{h})$ satisfying observational equivalence (Eq. \eqref{eq:coeff:obs:pair:Fourier}), counterfactual equivalence also holds.
\begin{theorem}[Coefficients comparison]\label{theorem:coefficients}
Let $m(x,y):= \exp(-i\xi_1 \cdot x) \exp(-i\xi_2 \cdot y)$. Suppose that the following observational equivalence holds for two mechanisms $(h^*, g^*)$ and $(\widetilde{h}, \widetilde{g})$ generating $(X^*,Y^*)$ and $(\widetilde{X},\widetilde{Y})$, respectively, in regime $\rho\in\left\{0, 1\right\}$,
\begin{align}
& \mathbb{E}\left[ m(\widetilde{X}, \widetilde{Y})|Z=z,\rho=1\right]=\mathbb{E}\left[ m(X^*, Y^*)|Z=z,\rho=1\right] = \varphi(z) \quad \forall z, \label{eq:obs:equiv}
\end{align}
such that, $(X^*,Y^*)|Z, \rho$ and $(\widetilde{X},\widetilde{Y})|Z, \rho$ are strongly complete mixtures of Gaussians with respect to $Z$. Then,
\[\mathbb{E}\left[ m(\widetilde{X}, \widetilde{Y})|Z=z,\rho=0\right]=\mathbb{E}\left[ m(X^*, Y^*)|Z=z,\rho=0\right] \quad\forall z.\]
This implies that in any regime $\rho$ the conditional expectation of $m(\cdot,\cdot)$ is invariant to the choice of the mechanism.
\end{theorem}
The theorem below relies on theorem \ref{theorem:coefficients} to show that Eq. \eqref{eq:coeff:obs:pair:CF:Fourier} 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{Appendix:Proof}).
\begin{theorem}[Identifiability]\label{Theorem:Identifiability}
Given assumptions \ref{asu:common}-\ref{asu:Integrability}, the intervention function is uniquely determined from the joint distribution of the triplets $(X,Y,Z)$. Namely, $\forall (\widetilde{g}, \widetilde{h})\in\widetilde{\mathcal{C}}$ it holds that,
\begin{equation}
H(y|x) = \tilde H(y|x)
\label{eq:H_equals_tilde_H}
\end{equation}
\end{theorem}
\begin{comment}
So far, we have formed the relationship between the observed and the interventional distributions through a system of Fredholm equations with an unknown kernel. In the ensuing section we develop an augmented Fourier series expansion method to show that the moments equality in Eq. \eqref{eq:moments} holds. Then, we show that the interventional distribution is uniquely determined as a linear combination of these moments (in theorem \ref{Theorem:Identifiability} to follow), implying that any pair of structural functions belonging to equivalence class $\widetilde{\mathcal{C}}$ satisfying Eq. \eqref{eq:int:star} yields the same interventional distribution. This will be done without reliance on the monotonicity assumption.
\end{comment}
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{lemma:Auxiliary}, Appendix \ref{Appendix:Proof}). 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{theorem:coefficients} establishes that the kernel can be unknown in Fredholm integral equations, whereas Theorem \ref{Theorem:Identifiability} establishes the identifiability of the counterfactual distribution based on the insight of theorem \ref{Theorem:Identifiability}.
\subsection{Synthetic counterfactual machinery}
Following theorem \ref{Theorem:Identifiability}, 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,
\begin{align}
& F_{Y,X|Z=z}^{\text{O}}(y,x|\widetilde{g},\widetilde{h}):=P\left(Y\le y,X\le x|Y=\widetilde{g}(X,\upomega,\upnu), X=\widetilde{h}(z,\upomega)\right)\label{equ:FYX:givenZ}\\
&=\int_{0}^1\int_{[0,1]^d}\mathds{1}\left\{\widetilde{g}(\widetilde{h}(z,\omega),\omega,\bm\nu)\le y\right\}\mathds{1}\left\{\widetilde{h}(z,\omega)\le x\right\}d\bm\nu d\omega,\nonumber
\end{align}
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,
\begin{align}
& F_{X|Z=z}^{\text{O}}(x|\widetilde{h}):=P\left(X\le x|X=\widetilde{h}(z,\upomega)\right)=\int_{0}^1\mathds{1}\left\{\widetilde{h}(z,\omega)\le x\right\}\hspace{0.25em} d\omega. \label{equ:X:CF:givenZ}
\end{align}
Note that the key innovation embedded in \eqref{equ:FYX:givenZ} 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 \eqref{equ:FYX:CF:givenZ} 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$):
\begin{align}
& F_{Y,X|Z=z}^{\text{CF}}(y,x|\widetilde{g},\widetilde{h}):=P\left(Y\le y,X\le x|Y=\widetilde{g}(X,\widetilde{\upomega},\upnu), X=\widetilde{h}(z,\upomega)\right)\label{equ:FYX:CF:givenZ}\\
&=\int_{0}^1\int_{0}^1\int_{[0,1]^d}\mathds{1}\left\{\widetilde{g}(\widetilde{h}(z,\omega),\widetilde{\omega},\bm\nu)\le y\right\}\mathds{1}\left\{\widetilde{h}(z,\omega)\le x\right\}d\bm\nu d\omega \hspace{0.25em} d\widetilde{\omega}. \nonumber
\end{align}
The increase in entropy is due to the presence of two independent contextual covariates $\upomega$ and $\widetilde{\upomega}$ in \eqref{equ:FYX:CF:givenZ} introducing independent variations to $X$ and to $Y$. This yields the following counterfactual result:
\begin{align}
& \mathbb{E}\left[\mathds{1}\left\{Y\le y\right\}|\text{do}(X=x); \widetilde{g},\widetilde{h}\right] = \mathbb{E}_{z\sim F_Z}\left[\left.\frac{\frac{\partial}{\partial x}F_{Y,X|Z}^{\text{CF}}(y,x;\widetilde{g},\widetilde{h})}{\frac{\partial}{\partial x}F_{X|Z}^{\text{O}}(x;\widetilde{g},\widetilde{h})}\right|Z=z\right].\label{equ:Expected:FY:CF:givenX}
\end{align}
The result in \eqref{equ:Expected:FY:CF:givenX} is propagated by the unconstrained support domain of \eqref{equ:FYX:CF:givenZ} unlike in models employing the additive control function approach \citep{heckman1985alternative}, the non-additive control function approach \citep{newey1994kernel,imbens2009identification,blundell2014control,hoderlein2009identification}, the generated covariates \citep{mammen2012nonparametric}, the non-separable counterfactual distribution based on quantiles \citep{chernozhukov2013inference}, the treatment of unobservables in \cite{schennach2014entropic} and the average causal effects by employing do-interventions \citep{pearl2019interpretation}. Further discussion pertaining to the effect of limited empirical support on counterfactual analysis is detailed in section \ref{Sec:Restricted:Support}.
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{figure:CounterFactualRectSupport}:
\begin{figure}[H]\centering\caption{\label{figure:CounterFactualRectSupport}Counterfactual Machinery vs. Partial Means \\ (Unrestricted vs. Restricted support)}\vspace{2mm}
\includegraphics[width=16cm,height=5cm]{CounterFactualRectSupport}
\end{figure}
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{figure:CounterFactualRectSupport}.
Our theoretical result (theorem \ref{Theorem:Identifiability}) 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{Sec:Estimation} 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.
\section{Estimation: contaminated GAN (CONGAN)}\label{Sec:Estimation}
The formulation of the system of Fredholm equations in \eqref{eq:integral:normalized} (in Appendix \ref{Sec:integral}) 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 \citep{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 \textit{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 \textit{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 \citep{singh2018nonparametric}:
\begin{asu}\label{Generators:assumption}
$\forall(\widetilde{h},\widetilde{g})\in\mathcal{(G,H)}$ the distributions of $\widetilde{g}(X, \omega,\nu)$ and $\widetilde{h}(Z, \omega)$ are absolutely continuous with respect to $\mu$. The joint distribution of $(Y,X,Z)$ is absolutely continuous with respect to $\mu$.\\
\end{asu}
We define the contaminated GAN (CONGAN) approach as an extension of the original GAN \citep{goodfellow2014generative} as follows:
\begin{align}
V(d,(\widetilde{h},\widetilde{g})) = \left\{\mathbb{E}_{\substack{X,Y,Z\sim F_{X,Y,Z}}}\left[\log\left(d(X,Y,Z)\right)\right] + \mathbb{E}_{\substack{\upomega\sim U[0,1]\\\nu\sim U[0,1]\\Z\sim F_Z}}\left[\log\left(1-d\left(\widetilde{h}(Z,\omega), \widetilde{g}\left(\widetilde{h}(z,\omega),\omega,\bm\nu\right)\right),Z\right)\right]\right\}\label{Cont:GAN}
\end{align}
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 \textit{contamination} braces in \eqref{Cont:GAN}. Namely, \eqref{Cont:GAN} 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:
\begin{align}
d=\underset{d\in\mathcal{D}}{\arg\max}\hspace{0.25em} V(\cdot,(\widetilde{h},\widetilde{g}))
\end{align}
Then, the optimal pair of generator functions $(\widetilde{h},\widetilde{g})$ is the one solving the $\min\max$ problem:
\begin{align}
(\widetilde{h},\widetilde{g})=\underset{(\widetilde{h},\widetilde{g})\in\mathcal{(G,H)}}{\arg\min}\hspace{0.25em}\left\{\underset{d\in\mathcal{D}}{\arg\max}\hspace{0.25em} V(d,(\widetilde{h},\widetilde{g}))\right\}
\end{align}
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{Generators:assumption} holds. Then following Theorem 2.2 in \cite{singh2018nonparametric},
\begin{align}
V(d,(\widetilde{h},\widetilde{g})) = \int_{x,y,z}\left\{\log\left(d(x,y,z)\right)f_{Y,X|Z}(y,x|z) + \log\left(1-d(x,y,z)\right)f_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h})\right\}f_z(z)dx \hspace{0.2em}dy\hspace{0.2em} dz
\end{align}
$\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:
\begin{align}
d(x,y,z) = \frac{f_{Y,X|Z}(y,x|z)}{f_{Y,X|Z}(y,x|z) + f_{Y,X|Z;\widetilde{g},\widetilde{h}}^{\text{O}}(y,x|z)}, \hspace{1em} \forall \hspace{0.2em}(x,y,z) \in \text{Supp}(X,Y,Z)
\end{align}
This yields:
\begin{align}
& \underset{d\in\mathcal{D}}{\max}\hspace{0.25em} V(d,(\widetilde{h},\widetilde{g})) = \int_{x,y,z}\left\{f_{Y,X|Z}(y,x|z) \log\left(\frac{f_{Y,X|Z}(y,x|z)}{f_{Y,X|Z}(y,x|z) + f_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h})}\right)\right\}f_z(z)dx \hspace{0.2em}dy\hspace{0.2em} dz\label{Divergence:loss:function}\\
& + \int_{x,y,z}\left\{f_{Y,X|Z}(y,x|z)\log\left(\frac{f_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h})}{f_{Y,X|Z}(y,x|z) + f_{Y,X|Z}^{\text{O}}(y,x|z;\widetilde{g},\widetilde{h})}\right)f_z(z)dx \hspace{0.2em}dy\hspace{0.2em} dz\right\}\nonumber
\end{align}
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 \citep{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:
\begin{figure}[H]\centering\caption{\label{figure:NeuralNetworkSequentialGenerators}Neural Network with sequential generators ($\widetilde{h}$ and $\widetilde{g}$)}\vspace{2mm}
\includegraphics[width=8cm]{NeuralNetworkSequentialGenerators_omega}
\end{figure}
Figure \ref{figure:NeuralNetworkSequentialGenerators} 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 \eqref{equ:FYX:givenZ} as well as the counterfactual distribution in \eqref{equ:FYX:CF:givenZ}. This latent space consists of $\upomega$ and the random disturbance $\nu$.
Next, we characterize explicitly the Neural network architecture to estimate these functions.
\subsection{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:
\begin{align}
\widetilde{g}(X,\omega,\nu,\bm{\uptheta})=\bm{\mathcal{Y}}_{L}, \hspace{1em}\bm{\mathcal{Y}}_{0} = [X,\omega,\nu], \hspace{1em} \bm{\mathcal{Y}}_{\ell} = \upvarphi\left(\bm{W}_{\ell}^Y \bm{\mathcal{Y}}_{\ell-1}\right), \hspace{1em} \ell = 1,...,L
\end{align}
and
\begin{align}
\widetilde{h}(Z,\omega;,\bm{\upbeta})=\bm{\mathcal{X}}_{L}, \hspace{1em}\bm{\mathcal{X}}_{0} = [Z,\omega], \hspace{1em} \bm{\mathcal{X}}_{\ell} = \upvarphi\left(\bm{W}_{\ell}^X \bm{\mathcal{X}}_{\ell-1}\right), \hspace{1em} \ell = 1,...,L
\end{align}
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})$.
\vspace{1em}
It should be emphasized that by adopting a Shannon-divergence minimization criterion function, we depart from the practice of Wasserstein GAN (WGAN) \citep{athey2020}. Although WGANS can be trained more easily, it requires the utilization of a penalty function to alleviate various optimization difficulties \citep{gulrajani2017improved}, which may limit the support domain of the generated data distribution \citep{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 \citep{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:
\begin{align}
\widehat{p}_{\bm{\upsigma}}(x,y,z) = \frac{1}{N}\sum_{i=1}^N \frac{1}{\upsigma_x}K\left(\frac{X_i-x}{\upsigma_x}\right)\frac{1}{\upsigma_y}K\left(\frac{Y_i-y}{\upsigma_y}\right)\frac{1}{\upsigma_z}K\left(\frac{Z_i-z}{\upsigma_z}\right)
\end{align}
\begin{align}
\widehat{q}_{\bm{\upsigma}}(x,y,z) = \frac{1}{N}\sum_{i=1}^N \frac{1}{\upsigma_x}K\left(\frac{\widehat{X}_i-x}{\upsigma_x}\right)\frac{1}{\upsigma_y}K\left(\frac{\widehat{Y}_i-y}{\upsigma_y}\right)\frac{1}{\upsigma_z}K\left(\frac{Z_i-z}{\upsigma_z}\right)
\end{align}
Utilizing the smoothing Jensen-Shannon Divergence emanating from \eqref{Divergence:loss:function} \citep{sinn2018non}:
\begin{align}
\mathcal{K}_n(\bm{\uptheta},\bm{\upbeta},\bm{\upsigma}) = \frac{1}{N}\sum_{i=1}^N\log\frac{\widehat{p}_{\bm{\upsigma}}(X_i,Y_i,Z_i)}{\widehat{p}_{\bm{\upsigma}}(X_i,Y_i,Z_i)+\widehat{q}_{\bm{\upsigma}}(X_i,Y_i,Z_i)} + \frac{1}{N}\sum_{i=1}^N\log\frac{\widehat{q}_{\bm{\upsigma}}(\widehat{X}_i,\widehat{Y}_i,Z_i)}{\widehat{p}_{\bm{\upsigma}}(\widehat{X}_i,\widehat{Y}_i,Z_i)+\widehat{q}_{\bm{\upsigma}}(\widehat{X}_i,\widehat{Y}_i,Z_i)}\label{CONT:GAN:LOSS}
\end{align}
where $\bm{\upsigma}:=(\upsigma_x,\upsigma_y,\upsigma_z)$ is the bandwidth vector. The optimal generator's strategy is to minimize the divergence:
\begin{align}
(\widehat{\bm{\uptheta}}^{(\text{iter+1})},\widehat{\bm{\upbeta}}^{(\text{iter+1})})=\underset{(\bm{\uptheta},\bm{\upbeta})}{\arg\min}\hspace{0.2em} \mathcal{K}_n(\bm{\uptheta},\bm{\upbeta},\bm{\upsigma}^{(\text{iter})})
\end{align}
\begin{comment}
\subsection{Lack of exchangeability between $\nu$ and $\omega$ in the Neural network}
In our model there is a formidable challenge in that the context variable cannot be determined uniquely by the observables, unlike in the common practice in the literature. The latter largely relies on the strict monotonicity of the structural equation of $X$ in the latent confounder. In fact, unlike the case of our latent context variable, monotonicity implies the degeneracy of this latent confounder when all the observables are held fixed. Its degeneracy renders it an identifiable covariate to be plugged-in the regression equation of $Y$ on $X$.
However, this plug-in strategy is no longer valid in our more general case of a non-degenerate latent context variable. In order to alleviate this challenge we develop and utilize a specific Neural network architecture, consisting of two noise variables $\upomega\sim U[0,1]$ and $\nu\sim U[0,1]$ which belong to the input space of the generator functions $\widetilde{g}$ and $\widetilde{h}$. In this architecture $\omega$ is common to both generators and $\nu$ is associated with $\widetilde{g}$ only. This practice induces the question of how one can distinguish between $\omega$ and $\nu$ without being able to plug-in one of them as a control variable in either one of the generators. However, this is possible due to the key role of the excluded variable $Z$ in this network. The rationale is that because $Z$ is excluded from the structural equation of $Y$ it can only affect $Y$ through its affect on $\omega$ when $X$ is held fixed. Consequently, for any arbitrary values $x$ and $y$, any joint variation in $Z$ and $f_{Y|X=x,Z}(y)$ is related to $\omega$.
Another aspect is that if the support of $Y$ for a given $X=x$ and $Z=z$ contains more values than the support of $\omega$ a given $\widehat{X}=x$ and $Z=z$, it follows that the source of variation in $Y$ is at least partially related to $\nu$. The main implications of this analysis is the lack of exchangeability between $\nu$ and $\omega$. This result is essential for uncovering $(\widetilde{g}, \widetilde{h})$.
\end{comment}
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 \eqref{X:Structural} and in \eqref{Y:Structural}, respectively. The latter will be used as a benchmark for verifying the accuracy of the estimate obtained from the synthetic data.
\subsection{Synthetic counterfactual machine vs. partial means estimator \citep{newey1994kernel}}\label{Sec:Restricted:Support}
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 \textit{real} counterfactual (CF) data set from the structural functions $g$ and $h$:
\begin{align}
& X_{\text{real,i}}^{\text{CF}} = h^*\left(Z_i,\omega_{i}^{\text{real}}\right), \hspace{1em} \omega_{i}^{\text{real}}, \widetilde{\omega}_{i}^{\text{real}}, \nu_i^{\text{real}} \sim U[0,1].\\
& Y_{\text{real,i}}^{\text{CF}} = g\left(X_{\text{real,i}},F_{\bm\epsilon|\upeta}^{-1}\left(\nu_i^{\text{real}}|F_{\upeta}^{-1}\left(\widetilde{\omega}_{i}^{\text{real}}\right)\right)\right).
\end{align}
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 \textit{synthetic (synt)} counterfactual (CF) data set from the synthetic (estimated) structural functions $\widetilde{g}$ and $\widetilde{h}$:
\begin{align}
& X_{\text{synt,i}}^{\text{CF}} = \widetilde{h}\left(Z,\omega_{i}^{\text{synt}};\bm{\upbeta}\right), \hspace{1em} \omega_{i}^{\text{synt}}, \widetilde{\omega}_{i}^{\text{synt}}, \nu_i^{\text{synt}} \sim U[0,1]\\
& Y_{\text{synt,i}}^{\text{CF}} = \widetilde{g}\left(X_{\text{synt,i}},\widetilde{\omega}_{i}^{\text{synt}}, \nu_i^{\text{synt}};\bm{\uptheta}\right).
\end{align}
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 \textit{real} observed (O) data set from the structural functions $g$ and $h$:
\begin{align}
& X_{\text{real,i}}^{\text{O}} = h^*\left(Z_i,\omega_{i}^{\text{real}}\right), \hspace{1em} \omega_{i}^{\text{real}}, \nu_i^{\text{real}} \sim U[0,1]\\
& Y_{\text{real,i}}^{\text{O}} = g\left(X_{\text{real,i}},F_{\bm\epsilon|\upeta}^{-1}\left(\nu_i^{\text{real}}|F_{\upeta}^{-1}\left(\omega_{i}^{\text{real}}\right)\right)\right),
\end{align}
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 \textit{synthetic (synt)} observed (O) data set from the synthetic (estimated) structural functions $\widetilde{g}$ and $\widetilde{h}$:
\begin{align}
& X_{\text{synt,i}}^{\text{O}} = \widetilde{h}\left(Z,\omega_{i}^{\text{synt}};\bm{\upbeta}\right), \hspace{1em} \omega_{i}^{\text{synt}}, \nu_i^{\text{synt}} \sim U[0,1]\\
& Y_{\text{synt,i}}^{\text{O}} = \widetilde{g}\left(X_{\text{synt,i}},\omega_{i}^{\text{synt}}, \nu_i^{\text{synt}};\bm{\uptheta}\right).
\end{align}
The estimated conditional expected outcomes given the action in the \textit{real} data is:
\begin{align}
& \widehat{\mathbb{E}}\left[Y_{\text{real}}|X=x\right] = \widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{O}}|X_{\text{real}}^{\text{O}}=x\right] = \frac{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N Y_{\text{real,i}}^{\text{O}}K\left(\frac{X_{\text{real,i}}^{\text{O}}-x}{\upsigma_x}\right)}{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N K\left(\frac{X_{\text{real,i}}^{\text{O}}-x}{\upsigma_x}\right)}.\label{real:O}
\end{align}
The estimated conditional expected outcomes given the action in the \textit{synthetic (synt)} data is:
\begin{align}
& \widehat{\mathbb{E}}\left[Y_{\text{synt}}|X=x\right] = \widehat{\mathbb{E}}\left[Y_{\text{synt}}^{\text{O}}|X_{\text{synt}}^{\text{O}}=x\right] \frac{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N Y_{\text{synt,i}}^{\text{O}}K\left(\frac{X_{\text{synt,i}}^{\text{O}}-x}{\upsigma_x}\right)}{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N K\left(\frac{X_{\text{synt,i}}^{\text{O}}-x}{\upsigma_x}\right)}.\label{synthetic:O}
\end{align}
The estimated counterfactual expected outcomes given the action in the \textit{real} data is:
\begin{align}
& \widehat{\mathbb{E}}\left[Y_{\text{real}}|\text{do}(X_{\text{real}}=x)\right] = \widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{CF}}|X_{\text{real}}^{\text{CF}}=x\right] = \frac{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N Y_{\text{real,i}}^{\text{CF}}K\left(\frac{X_{\text{real,i}}^{\text{CF}}-x}{\upsigma_x}\right)}{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N K\left(\frac{X_{\text{real,i}}^{\text{CF}}-x}{\upsigma_x}\right)}.\label{real:CF}
\end{align}
The estimated counterfactual expected outcomes given the action in the \textit{synthetic (synt)} data is:
\begin{align}
& \widehat{\mathbb{E}}\left[Y_{\text{synt}}|\text{do}(X_{\text{synt}}=x)\right] = \widehat{\mathbb{E}}\left[Y_{\text{synt}}^{\text{CF}}|X_{\text{synt}}^{\text{CF}}=x\right] \frac{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N Y_{\text{synt,i}}^{\text{CF}}K\left(\frac{X_{\text{synt,i}}^{\text{CF}}-x}{\upsigma_x}\right)}{\frac{1}{N}\frac{1}{\upsigma_x}\sum_{i=1}^N K\left(\frac{X_{\text{synt,i}}^{\text{CF}}-x}{\upsigma_x}\right)}.\label{synthetic:CF}
\end{align}
The partial means estimator \citep{newey1994kernel} is presented here as a benchmark, as it is applied to the \textit{observed} data set, unlike the counterfactual expression applied to the \textit{counterfactual} dataset in \eqref{real:CF}:
\begin{align}
& \widehat{\widehat{\mathbb{E}}}\left[Y_{\text{real}}|\text{do}(X_{\text{real}}=x)\right] =\frac{1}{N}\sum_{i=1}^N \widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{O}}|X_{\text{real}}^{\text{O}}=x, \upomega^{\text{real}}=\omega_{i}^{\text{real}}\right]\label{Partial:Means}
\end{align}
with
\begin{align}
& \widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{O}}|X_{\text{real}}^{\text{O}}=x, \upomega^{\text{real}}=\omega\right] = \frac{\frac{1}{N}\frac{1}{\upsigma_x}\frac{1}{\upsigma_{\omega}}\sum_{i=1}^N Y_{\text{real,i}}^{\text{O}}K\left(\frac{X_{\text{real,i}}^{\text{O}}-x}{\upsigma_x}\right)K\left(\frac{\omega_i^{\text{real}}-\omega}{\upsigma_{\omega}}\right)}{\frac{1}{N}\frac{1}{\upsigma_x}\frac{1}{\upsigma_{\omega}}\sum_{i=1}^N K\left(\frac{X_{\text{real,i}}^{\text{O}}-x}{\upsigma_x}\right)K\left(\frac{\omega_i^{\text{real}}-\omega}{\upsigma_{\omega}}\right)}.
\end{align}
\begin{comment}
and similarly,
\begin{align}
& \widehat{\widehat{\mathbb{E}}}\left[Y_{\text{synt}}|\text{do}(X_{\text{synt}}=x)\right] = \frac{1}{N}\sum_{i=1}^N \widehat{\mathbb{E}}\left[Y_{\text{synt}}^{\text{O}}|X_{\text{synt}}^{\text{O}}=x, \upomega^{\text{synt}}=\omega_{i}^{\text{synt}}\right]
\end{align}
with
\begin{align}
& \frac{1}{N}\sum_{i=1}^N \widehat{\mathbb{E}}\left[Y_{\text{synt}}^{\text{O}}|X_{\text{synt}}^{\text{O}}=x, \upomega^{\text{synt}}=\eta\right] = \frac{\frac{1}{N}\frac{1}{\upsigma_x}\frac{1}{\upsigma_{\eta}}\sum_{i=1}^N Y_{\text{synt,i}}^{\text{O}}K\left(\frac{X_{\text{real,i}}^{\text{O}}-x}{\upsigma_x}\right)K\left(\frac{\omega_i^{\text{synt}}-\eta}{\upsigma_{\eta}}\right)}{\frac{1}{N}\frac{1}{\upsigma_x}\frac{1}{\upsigma_{\eta}}\sum_{i=1}^N K\left(\frac{X_{\text{synt,i}}^{\text{O}}-x}{\upsigma_x}\right)K\left(\frac{\omega_i^{\text{synt}}-\eta}{\upsigma_{\eta}}\right)}.
\end{align}
\end{comment}
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 \eqref{Partial:Means} 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)$ \citep{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:
\begin{align}
& \widehat{\widehat{\mathbb{E}}}\left[Y_{\text{real}}|\text{do}(X_{\text{real}}=x)\right] =\frac{1}{N}\sum_{i=1}^N \widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{O}}|X_{\text{real}}^{\text{O}}=x, V_{\text{real}}=V_{\text{real,i}}\right]\label{Partial:Means:Imbens:Newey}
\end{align}
This is the benchmark used for comparison of our results. The reasons for choosing equation \eqref{Partial:Means:Imbens:Newey} 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 \citep{chernozhukov2020}.
\section{Simulations}\label{Sec: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 \citep{Hanoch1965}; the Almost Ideal Demand System (AIDS) \citep{deaton1980almost}; the Constant Elasticity of Substitution (CES) production function \citep{Sato1967} as well as the Translog utility function \citep{christensen1975}, which may also affect the demand and supply curves. A specific model of demand and supply functions have been discussed in \cite{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:
\begin{enumerate}
\item Trans-log utility functions: \label{Translog:model}
\begin{align}
& h_1(Z,\upeta) = \upalpha_0 + \upalpha_1 \log(Z) + \upalpha_2\log(\upeta ) + \upalpha_3\log^2(Z)+ \upalpha_4\log^2(\upeta) + \upalpha_5\log(Z)\log(\upeta)\\
& g_1(X,\bm\epsilon) = \upbeta_0 + \upbeta_1 \log(X) + \upbeta_2\log(\epsilon ) + \upbeta_3\log^2(X)+ \upbeta_4\log^2(\epsilon) + \upbeta_5\log(X)\log(\epsilon)
\end{align}
\item Almost Ideal Demand System (AIDS) functions: \label{AIDS:model}
\begin{align}
& h_2(Z,\upeta) = h_1(Z,\upeta)+ \upalpha_6 (Z\upeta)^{\uprho}\\
& g_2(X,\bm\epsilon) = g_1(X,\bm\epsilon)+ \upbeta_6 (X\epsilon)^{\uprho}
\end{align}
\item Constant Elasticity of Substitution (CES) production functions: \label{CES:model}
\begin{align}
& h_3(Z,\upeta) = \upalpha_1\left[\upalpha_2 Z^{-\uprho}+\upalpha_3 \upeta^{-\uprho}\right]^{-\upalpha_4/\uprho}\\
& g_3(X,\bm\epsilon) = \upbeta_1\left[\upbeta_2 X^{-\uprho}+\upbeta_3 \epsilon^{-\uprho}\right]^{-\upbeta_4/\uprho}
\end{align}
\item Backward-bending supply curves for $\upalpha_1=0, \upalpha_2=1$ and $\upbeta_1=0, \upbeta_2=1$: \label{Hanoch:model}
\begin{align}
& h_4(Z,\upeta) = \exp(Z\upeta - \upalpha_1 Z)-\upalpha_2 (Z\upeta-\upalpha_1 Z)\\
& g_4(X,\bm\epsilon) = \exp(X\epsilon-\upbeta_1 X)-\upbeta_2(X\epsilon-\upbeta_1 X)
\end{align}
\item A generalized cobb-Douglas production function (Cobb-Douglas is obtained for $\upalpha_1=\upalpha_2=0$ and $\upbeta_1=\upbeta_2=0$): \label{CobbDouglas:model}
\begin{align}
& h_5(Z,\upeta) = \upalpha_1 Z + \upalpha_2\upeta +\upalpha_3 Z^{\upalpha_4}\upeta^{\upalpha_5}\\
& g_5(X,\bm\epsilon) = \upalpha_1 X + \upalpha_2\epsilon +\upbeta_3 X^{\upbeta_4}\epsilon^{\upbeta_5}
\end{align}
\item Nonlinearly non-additively-separable error terms \label{tanh:model}
\begin{align}
& h_6(Z,\upeta) = \upalpha_1\tanh(\upalpha_2 Z+\upalpha_3 \upeta)+\upalpha_4\tanh(\upalpha_5 \upeta)\\
& g_6(X,\bm\epsilon) = \upbeta_1\tanh(\upbeta_2 Z+\upbeta_3 \epsilon)+\upbeta_4\tanh(\upbeta_5 \epsilon)
\end{align}
with $\tanh(x)=\frac{\exp(x)-\exp(-x)}{\exp(x)+\exp(-x)}$ commonly known as hyperbolic tangent.
\end{enumerate}
The error terms are generated as follows: $\upeta, \nu\sim N(0,1)$ such that $\upeta\protect\mathpalette{\protect\independenT}{\perp}\nu$,
\begin{align}
\epsilon = \upgamma_1\upeta\nu + \upgamma_2\upeta + \upgamma_3\nu, \hspace{1em} \upgamma_1,\upgamma_2,\upgamma_3\in\left\{-3,-2, -1, 0, 1, 2, 3\right\}
\end{align}
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{Translog:model}; the Almost Ideal Demand System (AIDS) functions in \ref{AIDS:model}; the Constant Elasticity of Substitution (CES) production functions in \ref{CES:model}; the backward-bending supply curves in \ref{Hanoch:model}; a Cobb-Douglas production function in \ref{CobbDouglas:model} and the hyperbolic tangent model in \ref{tanh:model}. 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{tab:CES:Triang:Average}-\ref{tab:BackBending:Average} (appendix \ref{Appendix:Tables}) 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 \citep{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 \citep{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:
\begin{enumerate}
\item \textbf{Observed Similarity} between $\widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{O}}|X_{\text{real}}^{\text{O}}=x\right]$ and $\widehat{\mathbb{E}}\left[Y_{\text{synt}}^{\text{O}}|X_{\text{synt}}^{\text{O}}=x\right]$ in \eqref{real:O} and \eqref{synthetic:O}, respectively.
It measures the resemblance between the estimates of the \textit{observed} expected outcome of $Y$ (at each quantile of $X$) in the synthetic and real data sets.\label{item:Observed:Similarity}
\item \textbf{Counterfactual Similarity} between $\widehat{\mathbb{E}}\left[Y_{\text{real}}^{\text{CF}}|X_{\text{real}}^{\text{CF}}=x\right]$ and $\widehat{\mathbb{E}}\left[Y_{\text{synt}}^{\text{CF}}|X_{\text{synt}}^{\text{CF}}=x\right]$ in \eqref{real:CF} and \eqref{synthetic:CF}, respectively.
It measures the resemblance between the estimates of the \textit{counterfactual} expected outcome of $Y$ (at each quantile of $X$) in the synthetic and real data sets.\label{item:CF:Similarity}
\item \textbf{Strict monotonicity} by comparing the estimated counterfactual results in \eqref{Partial:Means}, which does not rely on monotonicity, to the partial means estimate under monotonicity in \eqref{Partial:Means:Imbens:Newey}. Both of these estimators utilize the restricted empirical support only.
The degree of non-monotonicity is quantified in these simulations by comparing the results that would have been obtained if $\upomega$ were known (the partial means estimator) to the actual results obtained by using $\widehat{\upomega}=$, an estimate of it (the control variable estimator). Consequently, a negligible disparity between the estimates of the control variable and the partial means indicate that the monotonicity assumption could have been imposed.\footnote{Namely, the control variable estimator, $\widehat{V}_i=\widehat{F}_{X|Z}(x_i|z_i)$, is an estimate of the realization of $\upomega$ in the $i$th observation, whereas the partial means estimator, $V_i=F_{\upeta}(\eta_i)$, is the actual realization of $\omega$ in the $i$th observation. The imposition of strict monotonicity is required as in practice $\omega$ is unknown and thus it must be evaluated.} \label{item:MONO:Similarity}
\item \textbf{Restricted support} by comparing our real data counterfactual estimate in \eqref{real:CF} utilizing the unrestricted support to the partial means in \eqref{Partial:Means} utilizing the restricted support. Both of these estimators do not rely on monotonicity. However, the partial means estimator is available only in theory in which latent context is treated as known.
Any disparity between the partial means and the real counterfactual expected outcomes are related to discrepancies between the empirical and theoretical support of the context and action variable.
\label{item:MEANS:Similarity}
\end{enumerate}
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) \citep{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{tab:CES:Triang:Average}-\ref{tab:BackBending:Average} (appendix \ref{Appendix:Tables}). 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{tab:CES:Triang:Average} (appendix \ref{Appendix:Tables}), 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{tab:CES:Triang:Average}, column ($e$)) 13.94 (0.344) for the real counterfactual and (table \ref{tab:CES:Triang:Average}, 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{item:MEANS:Similarity}). Just as an example, in the 90th quantile of $X$, we report an estimate of (table \ref{tab:CES:Triang:Average}, column ($e$)) 15.73 (0.728) for the real counterfactual and (table \ref{tab:CES:Triang:Average}, 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{item:MONO:Similarity}). For instance, in the 5th quantile of $X$, we get the estimates (standard error) of (table \ref{tab:CES:Triang:Average}, column $(g)$) 6.98 (0.338) for the control variable and (table \ref{tab:CES:Triang:Average}, column $(h)$) 6.96 (0.276) for the partial means; in the 25th quantile of $X$, the estimate is (table \ref{tab:CES:Triang:Average}, column $(g)$) 8.52 (0.213) for the control variable and (table \ref{tab:CES:Triang:Average}, column $(h)$) 8.51 (0.188) for the partial means; in the 75th quantile of $X$, it is (table \ref{tab:CES:Triang:Average}, column $(g)$) 15.58 (0.338) for the control variable and (table \ref{tab:CES:Triang:Average}, column $(h)$) 15.59 (0.276) for the partial means.
In table \ref{tab:tranglog:Average} (appendix \ref{Appendix:Tables}) 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{tab:tranglog:Average}, column ($e$)) 16.63 (0.175) for the real counterfactual and (table \ref{tab:tranglog:Average}, 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{tab:tranglog:Average}, column $(g)$) 20.95 (1.087) for the control variable and (table \ref{tab:tranglog:Average}, 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{tab:tranglog:Average}, column ($e$)) 15.88 (0.190) for the real counterfactual and (table \ref{tab:tranglog:Average}, column $(h)$) 13.53 (0.959) for the partial means.
The entries in table \ref{tab:tanhtanh:Average} (appendix \ref{Appendix:Tables}) 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{tab:tanhtanh:Average}, column $(g)$) -0.18 (0.291) for the control variable and (table \ref{tab:tanhtanh:Average}, 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{tab:tanhtanh:Average}, column ($e$)) -0.17 (0.038) for the real counterfactual and (table \ref{tab:tanhtanh:Average}, column $(h)$) -0.70 (0.212) for the partial means.
The entries in table \ref{tab:AIDS:Average} (appendix \ref{Appendix:Tables}) 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{tab:AIDS:Average}, column ($e$)) 16.71 (0.188) for the real counterfactual and (table \ref{tab:AIDS:Average}, 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{tab:AIDS:Average}, column ($e$)) 22.79 (0.669) for the real counterfactual and (table \ref{tab:AIDS:Average}, column $(h)$) 26.58 (0.975) for the partial means; For this quantile, the synthetic counterfactual result (table \ref{tab:AIDS:Average}, column ($f$)) is 21.66 (0.758). For instance, in the 55th quantile of $X$, we report an estimate of (table \ref{tab:AIDS:Average}, column $(g)$) 21.16 (1.101) for the control variable and (table \ref{tab:AIDS:Average}, column $(h)$) 16.67 (0.748) for the partial means; For this quantile, the real counterfactual result is (table \ref{tab:AIDS:Average}, column ($e$)) 16.06 (0.195) and the synthetic counterfactual result is (table \ref{tab:AIDS:Average}, 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{tab:BackBending:Average} (appendix \ref{Appendix:Tables}) 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{tab:BackBending:Average}, column $(g)$) -0.33 (0.135) for the control variable and (table \ref{tab:BackBending:Average}, 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{tab:BackBending:Average}, column $(e)$) is 0.61 (0.040) and the synthetic counterfactual result (table \ref{tab:BackBending:Average}, 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{tab:BackBending:Average}, column ($e$)) 0.47 (0.033) for the real counterfactual and (table \ref{tab:BackBending:Average}, column $(h)$) 0.21 (0.087) for the partial means. The synthetic counterfactual result is (table \ref{tab:BackBending:Average}, 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{item:Observed:Similarity} and \ref{item:CF:Similarity} tables \ref{tab:CES:Triang:Average}-\ref{tab:BackBending:Average} (appendix \ref{Appendix:Tables}) 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{item:MONO:Similarity}, 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{item:MEANS:Similarity}, 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.}
\begin{comment}
\begin{table}[H] \centering\caption{\label{tab:Translog:CobbDouglass}Do-intervention (Trans-log utility function)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & 0.95 & 15.53 & 15.42 & 7.87 & 6.98 & 10.91 \\ \hline
10 & 1.29 & 17.24 & 17.05 & 10.55 & 11.39 & 13.76 \\ \hline
15 & 1.49 & 18.09 & 17.91 & 12.18 & 4.73 & 16.32 \\ \hline
20 & 1.69 & 18.77 & 18.64 & 13.73 & 11.86 & 18.54 \\ \hline
25 & 1.81 & 19.15 & 19.06 & 14.72 & 13.96 & 20.73 \\ \hline
30 & 1.93 & 19.48 & 19.44 & 15.66 & 15.51 & 21.64 \\ \hline
35 & 2.01 & 19.67 & 19.66 & 16.24 & 17.13 & 21.82 \\ \hline
40 & 2.05 & 19.77 & 19.78 & 16.59 & 17.31 & 21.83 \\ \hline
45 & 2.13 & 19.96 & 20.00 & 17.21 & 17.99 & 21.66 \\ \hline
50 & 2.20 & 20.12 & 20.19 & 17.80 & 18.80 & 21.38 \\ \hline
55 & 2.26 & 20.25 & 20.35 & 18.28 & 19.34 & 21.13 \\ \hline
60 & 2.36 & 20.44 & 20.58 & 19.04 & 20.20 & 20.89 \\ \hline
65 & 2.45 & 20.61 & 20.78 & 19.75 & 20.99 & 20.83 \\ \hline
70 & 2.70 & 21.05 & 21.32 & 21.76 & 22.10 & 21.70 \\ \hline
75 & 2.86 & 21.31 & 21.62 & 23.02 & 23.12 & 22.73 \\ \hline
80 & 3.03 & 21.61 & 21.95 & 24.39 & 25.27 & 23.95 \\ \hline
85 & 3.25 & 22.02 & 22.36 & 26.15 & 26.78 & 25.29 \\ \hline
90 & 3.66 & 22.98 & 23.21 & 29.37 & 27.36 & 27.85 \\ \hline
95 & 4.25 & 25.29 & 25.11 & 34.12 & 31.88 & 33.36 \\ \hline
\end{tabular}
\end{table}
In table \ref{tab:Translog:CobbDouglass} we simulate $500$ observations by using model \ref{Translog:model} for $X$ and the generalized Cobb-Douglas model \ref{CobbDouglas:model} for $Y$ settings $(\upalpha_0,\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(0,1,2,0,0.65,-5.25)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(8,6,-1,1,1)$, $(\upgamma_1,\upgamma_2)=(1,1)$ (500 observations):
\begin{align}
& X=h_1\left(\frac{3}{2}\left(\tanh(0.05Z)+1\right), \frac{3}{2}\left(\tanh(0.25\eta)+1\right)\right)\label{X:simulation:Translog:one}\\
& Y=g_5\left(X,\bm\epsilon\right)
\end{align}
In equation \eqref{X:simulation:Translog:one} above the covariates are transformed to positive quantities as they represent units of goods for the utility equation. We note that $\tanh(x)=\frac{\exp(x)-\exp(-x)}{\exp(x)+\exp(-x)}$.
\begin{table}[H] \centering\caption{\label{tab:Translog:CobbDouglass:stochastic}Do-intervention (Trans-log utility function - stochastic gradient descent)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & 1.01 & 15.75 & 17.98 & 8.04 & 6.95 & 8.65 \\ \hline
10 & 1.51 & 18.14 & 19.97 & 12.04 & 12.47 & 17.85 \\ \hline
15 & 1.68 & 18.77 & 20.48 & 13.45 & 13.37 & 20.41 \\ \hline
20 & 1.86 & 19.31 & 20.90 & 14.85 & 14.44 & 23.41 \\ \hline
25 & 1.87 & 19.35 & 20.93 & 14.94 & 14.55 & 23.55 \\ \hline
30 & 1.97 & 19.63 & 21.14 & 15.77 & 15.24 & 24.22 \\ \hline
35 & 1.98 & 19.64 & 21.15 & 15.80 & 15.27 & 24.25 \\ \hline
40 & 1.99 & 19.68 & 21.18 & 15.92 & 15.41 & 24.31 \\ \hline
45 & 2.05 & 19.83 & 21.29 & 16.41 & 15.83 & 24.39 \\ \hline
50 & 2.13 & 20.03 & 21.43 & 17.06 & 16.38 & 24.26 \\ \hline
55 & 2.18 & 20.13 & 21.51 & 17.44 & 16.69 & 24.14 \\ \hline
60 & 2.33 & 20.45 & 21.73 & 18.66 & 17.70 & 23.40 \\ \hline
65 & 2.49 & 20.75 & 21.93 & 19.91 & 19.10 & 23.14 \\ \hline
70 & 2.65 & 21.03 & 22.12 & 21.20 & 21.19 & 23.36 \\ \hline
75 & 2.87 & 21.41 & 22.37 & 22.99 & 22.98 & 23.51 \\ \hline
80 & 3.13 & 21.84 & 22.68 & 25.01 & 24.61 & 23.92 \\ \hline
85 & 3.30 & 22.18 & 22.94 & 26.41 & 26.04 & 25.21 \\ \hline
90 & 3.63 & 22.97 & 23.65 & 29.07 & 28.01 & 27.96 \\ \hline
95 & 4.29 & 25.53 & 26.21 & 34.34 & 30.13 & 32.62\\ \hline
\end{tabular}
\end{table}
\begin{table}[H] \centering\caption{\label{tab:-monotonic:tanh}Do-intervention ($\tanh$-based non-monotonic function)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 2.09 & - 2.13 & - 1.98 & - 1.17 & - 1.22 & - 1.82 \\ \hline
10 & - 1.64 & - 2.02 & - 1.36 & - 0.88 & - 0.85 & - 1.48 \\ \hline
15 & - 1.35 & - 1.91 & - 0.69 & - 0.70 & - 0.79 & - 0.94 \\ \hline
20 & - 1.13 & - 1.73 & - 0.34 & - 0.55 & - 0.67 & - 0.43 \\ \hline
25 & - 0.92 & - 1.47 & - 0.30 & - 0.41 & - 0.43 & - 0.18 \\ \hline
30 & - 0.71 & - 1.14 & - 0.43 & - 0.27 & - 0.12 & - 0.10 \\ \hline
35 & - 0.53 & - 0.82 & - 0.55 & - 0.15 & - 0.15 & - 0.01 \\ \hline
40 & - 0.39 & - 0.56 & - 0.58 & - 0.06 & - 0.13 & 0.04 \\ \hline
45 & - 0.18 & - 0.16 & - 0.55 & 0.08 & - 0.33 & 0.14 \\ \hline
50 & 0.04 & 0.29 & - 0.44 & 0.24 & 0.30 & 0.24 \\ \hline
55 & 0.15 & 0.52 & - 0.34 & 0.31 & 0.29 & 0.29 \\ \hline
60 & 0.32 & 0.88 & - 0.08 & 0.43 & 0.57 & 0.37 \\ \hline
65 & 0.55 & 1.32 & 0.40 & 0.59 & 0.75 & 0.45 \\ \hline
70 & 0.72 & 1.62 & 0.76 & 0.71 & 0.99 & 0.49 \\ \hline
75 & 0.91 & 1.88 & 1.13 & 0.84 & 0.74 & 0.55 \\ \hline
80 & 1.12 & 2.09 & 1.49 & 0.99 & 1.15 & 0.79 \\ \hline
85 & 1.35 & 2.25 & 1.71 & 1.15 & 1.20 & 1.30 \\ \hline
90 & 1.68 & 2.40 & 1.75 & 1.39 & 1.34 & 1.83 \\ \hline
95 & 2.15 & 2.52 & 1.95 & 1.71 & 1.51 & 2.10 \\ \hline
\end{tabular}
\begin{minipage}[t]{.65\linewidth}
\small\vspace{0.5em}The estimator of \cite{imbens2009identification} for the do-interventional distribution capture conditional quantiles of the latent confounder $\eta$.
\end{minipage}
\end{table}
In table \ref{tab:-monotonic:tanh} we simulate $2,000$ observations by using model \ref{tanh:model} for both $X$ as well as $Y$ and setting $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4,\upalpha_5)=(6, 0.25, -0.5, 5, 0.5)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(6, 0.25, -0.5, 10, 0.5)$, $(\upgamma_1,\upgamma_2)=(-3,1)$.
In table \ref{tab:AIDS} the trigonometric and AIDS models \ref{trigonometric:model} and \ref{AIDS:model}, respectively, are used with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(2,3,2,2)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5,\upbeta_6,\uprho)=(50,10,-18.1,-18, 27,-1,-2, 0.25)$, $(\upgamma_1,\upgamma_2)=(1,1)$ (20,000 observations):
\begin{align}
& X=h_7(Z,\eta)\\
& Y=g_1\left(\frac{1}{2}\tanh(0.05x)+1, \frac{1}{2}\tanh(0.25\bm\epsilon)+1\right)\label{Y:simulation:one}
\end{align}
In equation \eqref{Y:simulation:one} above the covariates are transformed to positive quantities as they represent prices for the demand equation. We note that $\tanh(x)=\frac{\exp(x)-\exp(-x)}{\exp(x)+\exp(-x)}$.
\begin{table}[H] \centering\caption{\label{tab:AIDS}Do-intervention (AIDS function \citep{deaton1980almost})}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 4.18 & 12.93 & 13.65 & 15.55 & 15.48 & 11.44 \\ \hline
10 & - 3.44 & 12.58 & 13.13 & 14.20 & 13.94 & 11.26 \\ \hline
15 & - 3.08 & 12.41 & 12.89 & 13.34 & 13.31 & 11.29 \\ \hline
20 & - 2.73 & 12.24 & 12.67 & 12.63 & 12.70 & 11.40 \\ \hline
25 & - 2.29 & 12.03 & 12.41 & 12.09 & 12.02 & 11.53 \\ \hline
30 & - 1.68 & 11.72 & 12.06 & 11.53 & 11.47 & 11.46 \\ \hline
35 & - 1.03 & 11.41 & 11.71 & 10.38 & 10.48 & 11.12 \\ \hline
40 & - 0.45 & 11.14 & 11.42 & 11.20 & 11.14 & 10.94 \\ \hline
45 & 0.22 & 10.87 & 11.12 & 12.00 & 11.91 & 10.84 \\ \hline
50 & 0.86 & 10.64 & 10.86 & 11.28 & 11.36 & 10.75 \\ \hline
55 & 1.42 & 10.47 & 10.65 & 10.95 & 10.98 & 10.60 \\ \hline
60 & 1.87 & 10.35 & 10.49 & 10.73 & 10.75 & 10.35 \\ \hline
65 & 2.33 & 10.23 & 10.34 & 10.51 & 10.52 & 9.89 \\ \hline
70 & 2.69 & 10.15 & 10.22 & 10.37 & 10.35 & 9.35 \\ \hline
75 & 2.93 & 10.09 & 10.15 & 10.26 & 10.25 & 8.96 \\ \hline
80 & 3.15 & 10.05 & 10.08 & 10.15 & 10.10 & 8.61 \\ \hline
85 & 3.34 & 10.00 & 10.02 & 9.98 & 9.89 & 8.32 \\ \hline
90 & 3.56 & 9.95 & 9.95 & 9.58 & 9.52 & 8.06 \\ \hline
95 & 3.94 & 9.87 & 9.84 & 8.71 & 8.75 & 7.83 \\ \hline
\end{tabular}
\end{table}
In table \ref{tab:CES} the trigonometric and CES models \ref{trigonometric:model} and \ref{CES:model}, respectively, are used with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(2,3,2,2)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\uprho)=(6,0.5,0.5, 1.5,0.5)$, $(\upgamma_1,\upgamma_2)=(1,1)$ (20,000 observations):
\begin{align}
& X=h_7(Z,\eta)\\
& Y=g_3\left(\exp(0.5 X), \exp(0.5 \bm\epsilon)\right)\label{Y:simulation:CES}
\end{align}
In equation \eqref{Y:simulation:CES} above the covariates are transformed to positive quantities as they represent inputs for the production equation.
\begin{table}[H] \centering\caption{\label{tab:CES}Do-intervention (CES production function \citep{Sato1967})}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 4.76 & 9.38 & 8.67 & 13.25 & 12.86 & 7.81 \\ \hline
10 & - 3.42 & 10.35 & 9.21 & 6.12 & 6.15 & 9.11 \\ \hline
15 & - 3.11 & 10.65 & 9.49 & 5.84 & 6.00 & 9.58 \\ \hline
20 & - 2.70 & 11.10 & 9.96 & 6.22 & 6.26 & 10.19 \\ \hline
25 & - 2.09 & 11.89 & 10.86 & 7.91 & 7.83 & 10.99 \\ \hline
30 & - 1.17 & 13.42 & 12.65 & 11.13 & 11.09 & 13.02 \\ \hline
35 & - 0.60 & 14.57 & 13.98 & 11.63 & 11.94 & 14.54 \\ \hline
40 & - 0.12 & 15.62 & 15.17 & 14.00 & 13.89 & 15.71 \\ \hline
45 & 0.33 & 16.66 & 16.32 & 14.72 & 14.74 & 16.79 \\ \hline
50 & 0.83 & 17.88 & 17.62 & 14.82 & 15.22 & 17.91 \\ \hline
55 & 1.30 & 19.06 & 18.86 & 17.39 & 17.35 & 18.87 \\ \hline
60 & 1.72 & 20.15 & 19.99 & 18.90 & 19.04 & 19.84 \\ \hline
65 & 2.15 & 21.29 & 21.17 & 21.05 & 21.12 & 21.55 \\ \hline
70 & 2.60 & 22.50 & 22.42 & 23.30 & 23.40 & 23.93 \\ \hline
75 & 2.96 & 23.50 & 23.47 & 24.96 & 25.03 & 25.44 \\ \hline
80 & 3.26 & 24.40 & 24.42 & 26.31 & 26.33 & 26.28 \\ \hline
85 & 3.55 & 25.26 & 25.34 & 27.70 & 27.53 & 27.08 \\ \hline
90 & 4.30 & 27.76 & 27.97 & 34.70 & 34.34 & 40.01 \\ \hline
95 & 5.98 & 36.13 & 35.41 & 50.21 & 50.07 & 53.24\\ \hline
\end{tabular}
\end{table}
\begin{table}[H]\centering\caption{\label{tab:CobbDouglas}Do-intervention (Cobb-Douglas function)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline 5 & 0.33 & 6.54 & 6.57 & 6.63 & 6.64 & 6.86 \\ \hline
10 & 0.42 & 6.64 & 6.66 & 6.80 & 6.79 & 6.96 \\ \hline
15 & 0.53 & 6.78 & 6.84 & 7.00 & 6.99 & 7.10 \\ \hline
20 & 0.65 & 6.97 & 7.11 & 7.22 & 7.20 & 7.26 \\ \hline
25 & 0.77 & 7.16 & 7.40 & 7.43 & 7.37 & 7.44 \\ \hline
30 & 0.87 & 7.33 & 7.64 & 7.60 & 7.54 & 7.58 \\ \hline
35 & 0.97 & 7.51 & 7.86 & 7.76 & 7.72 & 7.72 \\ \hline
40 & 1.08 & 7.71 & 8.10 & 7.92 & 7.87 & 7.88 \\ \hline
45 & 1.20 & 7.94 & 8.36 & 8.09 & 8.03 & 8.04 \\ \hline
50 & 1.33 & 8.19 & 8.62 & 8.27 & 8.19 & 8.20 \\ \hline
55 & 1.46 & 8.44 & 8.83 & 8.43 & 8.36 & 8.36 \\ \hline
60 & 1.60 & 8.70 & 8.98 & 8.59 & 8.55 & 8.52 \\ \hline
65 & 1.76 & 8.99 & 9.08 & 8.76 & 8.76 & 8.70 \\ \hline
70 & 1.94 & 9.30 & 9.15 & 8.92 & 8.95 & 8.88 \\ \hline
75 & 2.13 & 9.62 & 9.21 & 9.07 & 9.12 & 9.05 \\ \hline
80 & 2.36 & 9.98 & 9.34 & 9.22 & 9.33 & 9.24 \\ \hline
85 & 2.62 & 10.35 & 9.52 & 9.37 & 9.60 & 9.50 \\ \hline
90 & 2.88 & 10.66 & 9.84 & 9.48 & 9.80 & 9.73 \\ \hline
95 & 3.32 & 11.04 & 10.99 & 9.62 & 9.99 & 9.99 \\ \hline
\end{tabular}
\end{table}
In table \ref{tab:CobbDouglas} the Cobb-Douglas model \ref{CobbDouglas:model} is used for both $X$ and $Y$ with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4,\upalpha_5)=(0,0,6,1,1)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(12,4,-8,1,1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (20,000 observations):
\begin{align}
& X=h_6\left(\frac{1}{2}\tanh(0.5Z)+1, \frac{1}{2}\tanh(0.5\eta)+1\right)\label{X:simulation:CobbDouglas}\\
& Y=g_6\left(\frac{1}{2}\tanh(0.5X)+1, \frac{1}{2}\tanh(0.5\bm\epsilon)+1\right)\label{Y:simulation:CobbDouglas}
\end{align}
In equations \eqref{X:simulation:CobbDouglas} and \eqref{Y:simulation:CobbDouglas} above the covariates are transformed to positive quantities as they represent inputs for the production equation. We note that $\tanh(x)=\frac{\exp(x)-\exp(-x)}{\exp(x)+\exp(-x)}$.
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\First specification}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 0.16 & 0.38 & 0.91 & - 0.02 & 0.53 & - 0.39 \\ \hline
10 & - 0.12 & 0.35 & 0.88 & 0.00 & 0.41 & - 0.41 \\ \hline
15 & 0.05 & 0.14 & 0.74 & 0.08 & 0.04 & - 0.51 \\ \hline
20 & 0.31 & - 0.30 & 0.61 & 0.21 & - 0.16 & - 0.73 \\ \hline
25 & 0.51 & - 0.66 & 0.66 & 0.31 & - 0.23 & - 1.04 \\ \hline
30 & 0.66 & - 0.97 & 0.49 & 0.38 & 0.61 & - 1.32 \\ \hline
35 & 0.87 & - 1.32 & 0.07 & 0.49 & 0.31 & - 1.65 \\ \hline
40 & 0.99 & - 1.37 & - 0.05 & 0.55 & 0.50 & - 1.73 \\ \hline
45 & 1.01 & - 1.37 & - 0.06 & 0.56 & 0.48 & - 1.73 \\ \hline
50 & 1.15 & - 1.22 & - 0.07 & 0.63 & 0.36 & - 1.66 \\ \hline
55 & 1.39 & - 0.62 & 0.07 & 0.75 & 0.67 & - 1.21 \\ \hline
60 & 1.61 & - 0.09 & 0.29 & 0.86 & 0.89 & - 0.68 \\ \hline
65 & 1.93 & 0.49 & 0.78 & 1.01 & 0.96 & - 0.20 \\ \hline
70 & 2.39 & 1.15 & 1.45 & 1.24 & 1.11 & 0.34 \\ \hline
75 & 2.99 & 1.85 & 1.58 & 1.54 & 1.69 & 0.99 \\ \hline
80 & 3.76 & 2.64 & 1.52 & 1.92 & 2.11 & 1.66 \\ \hline
85 & 4.79 & 3.45 & 1.43 & 2.43 & 2.41 & 2.28 \\ \hline
90 & 6.49 & 4.40 & 1.19 & 3.27 & 3.69 & 3.07 \\ \hline
95 & 10.04 & 5.63 & 1.13 & 5.03 & 5.26 & 4.43 \\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply} the backward-bending supply function \ref{Hanoch:model} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2)=(1,3)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(0.5, 1.2, 1, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (20,000 observations).
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply:two}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\Second specification}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 0.18 & 0.15 & 0.28 & - 0.02 & - 0.48 & - 0.06 \\ \hline
10 & - 0.11 & 0.13 & 0.24 & - 0.01 & - 0.43 & - 0.07 \\ \hline
15 & - 0.00 & 0.08 & 0.15 & - 0.00 & - 0.21 & - 0.08 \\ \hline
20 & 0.16 & - 0.00 & - 0.02 & 0.07 & 0.14 & - 0.12 \\ \hline
25 & 0.34 & - 0.10 & - 0.22 & 0.16 & 0.09 & - 0.19 \\ \hline
30 & 0.55 & - 0.21 & - 0.35 & 0.29 & 0.24 & - 0.29 \\ \hline
35 & 0.72 & - 0.32 & - 0.40 & 0.39 & 0.33 & - 0.42 \\ \hline
40 & 0.85 & - 0.36 & - 0.40 & 0.45 & 0.37 & - 0.50 \\ \hline
45 & 0.95 & - 0.36 & - 0.37 & 0.48 & 0.41 & - 0.52 \\ \hline
50 & 1.06 & - 0.31 & - 0.28 & 0.53 & 0.45 & - 0.50 \\ \hline
55 & 1.22 & - 0.15 & - 0.07 & 0.62 & 0.53 & - 0.38 \\ \hline
60 & 1.45 & 0.17 & 0.27 & 0.73 & 0.66 & - 0.07 \\ \hline
65 & 1.75 & 0.56 & 0.65 & 0.85 & 0.87 & 0.32 \\ \hline
70 & 2.13 & 0.93 & 0.96 & 1.06 & 1.10 & 0.64 \\ \hline
75 & 2.64 & 1.40 & 1.37 & 1.38 & 1.33 & 1.05 \\ \hline
80 & 3.39 & 1.83 & 1.88 & 1.73 & 1.78 & 1.45 \\ \hline
85 & 4.54 & 2.45 & 2.55 & 2.27 & 2.36 & 2.10 \\ \hline
90 & 6.41 & 3.11 & 3.25 & 3.18 & 3.23 & 2.79 \\ \hline
95 & 10.60 & 3.88 & 3.45 & 5.35 & 3.61 & 4.35
\\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply:two} the backward-bending supply function \ref{Hanoch:model} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2)=(1,3)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(0.5, 1.2, 1, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (20,000 observations). Hidden layers are: $[[16,4,3,2],[16,4,3,2]]$.
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply:three}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\(Parameterization II)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 2.35 & - 1.56 & - 1.75 & - 0.53 & - 0.53 & - 1.62 \\ \hline
10 & - 1.82 & - 1.23 & - 1.41 & - 0.42 & - 0.44 & - 1.16 \\ \hline
15 & - 1.50 & - 1.03 & - 1.21 & - 0.35 & - 0.38 & - 0.81 \\ \hline
20 & - 1.21 & - 0.86 & - 1.03 & - 0.29 & - 0.33 & - 0.74 \\ \hline
25 & - 1.04 & - 0.76 & - 0.92 & - 0.25 & - 0.29 & - 0.71 \\ \hline
30 & - 0.88 & - 0.67 & - 0.82 & - 0.22 & - 0.26 & - 0.66 \\ \hline
35 & - 0.68 & - 0.55 & - 0.69 & - 0.18 & - 0.21 & - 0.57 \\ \hline
40 & - 0.45 & - 0.42 & - 0.55 & - 0.12 & - 0.16 & - 0.45 \\ \hline
45 & - 0.20 & - 0.28 & - 0.39 & - 0.07 & - 0.10 & - 0.32 \\ \hline
50 & 0.00 & - 0.17 & - 0.26 & - 0.02 & - 0.06 & - 0.21 \\ \hline
55 & 0.17 & - 0.08 & - 0.16 & 0.01 & - 0.02 & - 0.13 \\ \hline
60 & 0.35 & 0.01 & - 0.05 & 0.05 & 0.03 & - 0.06 \\ \hline
65 & 0.53 & 0.11 & 0.06 & 0.10 & 0.07 & 0.01 \\ \hline
70 & 0.77 & 0.23 & 0.21 & 0.15 & 0.13 & 0.09 \\ \hline
75 & 1.03 & 0.36 & 0.36 & 0.21 & 0.19 & 0.23 \\ \hline
80 & 1.29 & 0.48 & 0.51 & 0.26 & 0.25 & 0.43 \\ \hline
85 & 1.56 & 0.61 & 0.66 & 0.32 & 0.31 & 0.66 \\ \hline
90 & 1.91 & 0.77 & 0.85 & 0.40 & 0.39 & 0.93 \\ \hline
95 & 2.40 & 0.98 & 1.10 & 0.50 & 0.50 & 1.17\\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply:three} the backward-bending supply function (parameterization II) \ref{Hanoch:model:two} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(-2.5084,3.8355,2.4216,0)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(-0.4646, ,-0.1683, -0.6171, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (500 observations). Hidden layers are: $[[80,4,3,2],[80,4,3,2]]$.
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply:five}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\(Parameterization II)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 2.60 & - 1.77 & - 1.82 & - 0.59 & - 0.98 & - 2.29 \\ \hline
10 & - 1.85 & - 1.31 & - 1.30 & - 0.43 & - 0.68 & - 1.88 \\ \hline
15 & - 1.48 & - 1.08 & - 1.06 & - 0.35 & - 0.55 & - 1.54 \\ \hline
20 & - 1.23 & - 0.94 & - 0.90 & - 0.30 & - 0.46 & - 1.23 \\ \hline
25 & - 0.98 & - 0.79 & - 0.75 & - 0.25 & - 0.38 & - 0.92 \\ \hline
30 & - 0.71 & - 0.64 & - 0.59 & - 0.20 & - 0.30 & - 0.66 \\ \hline
35 & - 0.55 & - 0.54 & - 0.49 & - 0.16 & - 0.24 & - 0.48 \\ \hline
40 & - 0.34 & - 0.43 & - 0.37 & - 0.12 & - 0.18 & - 0.33 \\ \hline
45 & - 0.15 & - 0.32 & - 0.27 & - 0.08 & - 0.13 & - 0.21 \\ \hline
50 & 0.02 & - 0.23 & - 0.17 & - 0.05 & - 0.08 & - 0.08 \\ \hline
55 & 0.21 & - 0.13 & - 0.07 & - 0.01 & - 0.02 & 0.05 \\ \hline
60 & 0.39 & - 0.03 & 0.03 & 0.03 & 0.03 & 0.17 \\ \hline
65 & 0.57 & 0.06 & 0.12 & 0.06 & 0.08 & 0.29 \\ \hline
70 & 0.78 & 0.17 & 0.23 & 0.10 & 0.13 & 0.39 \\ \hline
75 & 1.01 & 0.28 & 0.34 & 0.15 & 0.19 & 0.51 \\ \hline
80 & 1.23 & 0.40 & 0.45 & 0.19 & 0.25 & 0.62 \\ \hline
85 & 1.48 & 0.52 & 0.57 & 0.24 & 0.31 & 0.74 \\ \hline
90 & 1.75 & 0.65 & 0.69 & 0.29 & 0.38 & 0.86 \\ \hline
95 & 2.13 & 0.82 & 0.85 & 0.35 & 0.46 & 1.07\\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply:five} the backward-bending supply function (parameterization II) \ref{Hanoch:model:two} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(4.9556, 5.4055, 4.9914,0)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(0.6655,0.2399,0.3445, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (500 observations). Hidden layers are: $[[40,8,3,2],[40,8,3,2]]$.
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply:four}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\(Parameterization II)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 2.54 & - 1.79 & - 1.84 & - 0.46 & - 0.68 & - 1.28 \\ \hline
10 & - 1.83 & - 1.33 & - 1.35 & - 0.32 & - 0.49 & - 1.07 \\ \hline
15 & - 1.50 & - 1.14 & - 1.13 & - 0.25 & - 0.41 & - 0.94 \\ \hline
20 & - 1.25 & - 0.98 & - 0.97 & - 0.20 & - 0.35 & - 0.83 \\ \hline
25 & - 1.00 & - 0.83 & - 0.81 & - 0.15 & - 0.29 & - 0.77 \\ \hline
30 & - 0.78 & - 0.70 & - 0.67 & - 0.11 & - 0.23 & - 0.72 \\ \hline
35 & - 0.57 & - 0.58 & - 0.54 & - 0.07 & - 0.18 & - 0.64 \\ \hline
40 & - 0.38 & - 0.47 & - 0.43 & - 0.03 & - 0.13 & - 0.53 \\ \hline
45 & - 0.17 & - 0.36 & - 0.31 & 0.01 & - 0.08 & - 0.36 \\ \hline
50 & 0.03 & - 0.24 & - 0.20 & 0.05 & - 0.03 & - 0.16 \\ \hline
55 & 0.21 & - 0.15 & - 0.10 & 0.08 & 0.01 & 0.03 \\ \hline
60 & 0.40 & - 0.05 & 0.01 & 0.12 & 0.05 & 0.18 \\ \hline
65 & 0.58 & 0.05 & 0.10 & 0.15 & 0.10 & 0.29 \\ \hline
70 & 0.76 & 0.14 & 0.20 & 0.19 & 0.14 & 0.37 \\ \hline
75 & 0.94 & 0.23 & 0.29 & 0.22 & 0.18 & 0.43 \\ \hline
80 & 1.20 & 0.36 & 0.42 & 0.27 & 0.25 & 0.49 \\ \hline
85 & 1.52 & 0.52 & 0.58 & 0.33 & 0.32 & 0.55 \\ \hline
90 & 1.85 & 0.67 & 0.73 & 0.39 & 0.40 & 0.65 \\ \hline
95 & 2.25 & 0.86 & 0.92 & 0.47 & 0.49 & 0.86\\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply:four} the backward-bending supply function (parameterization II) \ref{Hanoch:model:two} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(-2.5084,3.8355,2.4216,0)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(-0.4646,-0.1683, -0.6171, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (500 observations). Hidden layers are: $[[80,4,3,2],[80,4,3,2]]$.
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply:six}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\(Parameterization II)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 2.38 & - 1.57 & - 1.68 & - 0.47 & - 0.35 & - 1.36 \\ \hline
10 & - 1.80 & - 1.18 & - 1.30 & - 0.35 & - 0.28 & - 1.07 \\ \hline
15 & - 1.43 & - 0.94 & - 1.06 & - 0.28 & - 0.23 & - 0.79 \\ \hline
20 & - 1.16 & - 0.78 & - 0.89 & - 0.23 & - 0.19 & - 0.62 \\ \hline
25 & - 0.92 & - 0.63 & - 0.75 & - 0.18 & - 0.16 & - 0.49 \\ \hline
30 & - 0.71 & - 0.51 & - 0.63 & - 0.14 & - 0.13 & - 0.39 \\ \hline
35 & - 0.51 & - 0.40 & - 0.52 & - 0.10 & - 0.10 & - 0.30 \\ \hline
40 & - 0.34 & - 0.30 & - 0.42 & - 0.07 & - 0.07 & - 0.22 \\ \hline
45 & - 0.16 & - 0.20 & - 0.32 & - 0.04 & - 0.04 & - 0.14 \\ \hline
50 & 0.01 & - 0.11 & - 0.23 & - 0.00 & - 0.01 & - 0.08 \\ \hline
55 & 0.19 & - 0.01 & - 0.14 & 0.03 & 0.02 & - 0.02 \\ \hline
60 & 0.38 & 0.09 & - 0.04 & 0.06 & 0.05 & 0.03 \\ \hline
65 & 0.56 & 0.18 & 0.05 & 0.10 & 0.08 & 0.08 \\ \hline
70 & 0.73 & 0.27 & 0.14 & 0.13 & 0.12 & 0.13 \\ \hline
75 & 0.89 & 0.35 & 0.22 & 0.16 & 0.15 & 0.20 \\ \hline
80 & 1.06 & 0.43 & 0.30 & 0.19 & 0.18 & 0.29 \\ \hline
85 & 1.30 & 0.55 & 0.41 & 0.23 & 0.23 & 0.43 \\ \hline
90 & 1.66 & 0.72 & 0.58 & 0.29 & 0.31 & 0.59 \\ \hline
95 & 2.23 & 0.97 & 0.83 & 0.40 & 0.44 & 0.76\\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply:six} the backward-bending supply function (parameterization II) \ref{Hanoch:model:two} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(9.3728, 4.5853, 1.8393,0)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(-0.4704, -0.1861, -5.089, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (500 observations). Hidden layers are: $[[40,4,3,2],[40,4,3,2]]$.
\begin{table}[H]\centering\caption{\label{tab:BackwardBendSupply:seven}Do-intervention (Backward-bending supply function \citep{Hanoch1965})\\(Parameterization II)}\vspace{10mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Quantile of $X$}} & \multicolumn{1}{c|}{\multirow{2}{*}{$X$}} & \multicolumn{2}{c|}{$\mathbb{E}[Y|X=x]$} & \multicolumn{3}{c|}{$\mathbb{E}[Y|\text{do}(X=x)]$} \\ \cline{3-7}
\multicolumn{1}{|c|}{} & \multicolumn{1}{c|}{} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{real} & \multicolumn{1}{c|}{synthetic (synt)} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}\textcolor{blue}{Imbens} \\ \textcolor{blue}{and} \\ \textcolor{blue}{Newey} \\ \textcolor{blue}{(2009)}\end{tabular}} \\ \hline
5 & - 2.18 & - 1.47 & - 1.27 & - 0.39 & - 0.33 & - 1.61 \\ \hline
10 & - 1.63 & - 1.13 & - 0.93 & - 0.29 & - 0.24 & - 1.17 \\ \hline
15 & - 1.26 & - 0.91 & - 0.72 & - 0.22 & - 0.18 & - 0.70 \\ \hline
20 & - 0.99 & - 0.76 & - 0.57 & - 0.17 & - 0.14 & - 0.40 \\ \hline
25 & - 0.78 & - 0.63 & - 0.44 & - 0.13 & - 0.11 & - 0.23 \\ \hline
30 & - 0.59 & - 0.53 & - 0.34 & - 0.10 & - 0.08 & - 0.11 \\ \hline
35 & - 0.44 & - 0.44 & - 0.26 & - 0.07 & - 0.05 & - 0.01 \\ \hline
40 & - 0.28 & - 0.35 & - 0.17 & - 0.04 & - 0.02 & 0.10 \\ \hline
45 & - 0.10 & - 0.26 & - 0.07 & - 0.01 & 0.01 & 0.20 \\ \hline
50 & 0.08 & - 0.16 & 0.03 & 0.03 & 0.04 & 0.29 \\ \hline
55 & 0.24 & - 0.08 & 0.11 & 0.06 & 0.07 & 0.37 \\ \hline
60 & 0.41 & 0.01 & 0.20 & 0.09 & 0.10 & 0.43 \\ \hline
65 & 0.61 & 0.12 & 0.31 & 0.13 & 0.14 & 0.50 \\ \hline
70 & 0.83 & 0.22 & 0.42 & 0.17 & 0.18 & 0.58 \\ \hline
75 & 1.07 & 0.34 & 0.54 & 0.22 & 0.23 & 0.63 \\ \hline
80 & 1.31 & 0.46 & 0.66 & 0.27 & 0.28 & 0.65 \\ \hline
85 & 1.52 & 0.56 & 0.76 & 0.31 & 0.33 & 0.67 \\ \hline
90 & 1.90 & 0.74 & 0.94 & 0.38 & 0.42 & 0.70 \\ \hline
95 & 2.49 & 1.00 & 1.21 & 0.50 & 0.58 & 0.88\\ \hline
\end{tabular}
\end{table}
In table \ref{tab:BackwardBendSupply:seven} the backward-bending supply function (parameterization II) \ref{Hanoch:model:two} is used for $X$ and the Cobb-Douglass model \ref{CobbDouglas:model} (with interactions) is used for $Y$ with parameter settings $(\upalpha_1,\upalpha_2,\upalpha_3,\upalpha_4)=(0.5229, 1.2882, -3.0170,0)$ and $(\upbeta_1,\upbeta_2,\upbeta_3,\upbeta_4,\upbeta_5)=(2.9507, -3.2428 , 4.7784, 1, 1)$, $(\upgamma_1,\upgamma_2)=(-3,1)$ (500 observations). Hidden layers are: $[[40,4,3,2],[40,4,3,2]]$.
\end{comment}
\begin{comment}
\subsection{Stability of result and divergence measures}
In this section we present our model performance and stability in terms of dissimilarity between the observed distribution of $(X,Y,Z)$ and its synthetic variant. Additionally, we quantify the degree of dissimilarity between the counterfactual distribution of $(X,Y^{\text{CF}},Z)$ and its synthetic variant. We opted for the Jensen-Shannon-Divergence (JSD) criterion function as a dissimilarity measure between distributions due to its direct relation to the objective function in \eqref{CONT:GAN:LOSS}. The $\text{JSD}(\widehat{f}_{X,Y,Z}^\text{O}, \widehat{f}_{X,Y,Z}^\text{Synt})$ equals zero if and only if $\widehat{f}_{X,Y,Z}^\text{O}$ and $\widehat{f}_{X,Y,Z}^\text{Synt}$ are identical.
\begin{table}[H]
\centering\caption{\label{tab:JSD:one}Jensen-Shannon Divergence (JSD) \\ Neural Network architecture\\
$[16, 4, 3, 2]$}\vspace{1mm}
\scaleobj{0.7}{\begin{threeparttable}\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
\hline
\multicolumn{2}{|c|}{\multirow{2}{*}{Model}} & \multirow{2}{*}{N} & \multicolumn{2}{c|}{\begin{tabular}[c]{@{}c@{}}Observed\\ CONGAN\end{tabular}} & \multicolumn{2}{c|}{\begin{tabular}[c]{@{}c@{}}Counterfactual\\ CONGAN\end{tabular}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Imbens and Newey\\ (2009)\end{tabular}} & \multicolumn{2}{c|}{\begin{tabular}[c]{@{}c@{}}No. of\\ iterations\end{tabular}} \\ \cline{4-7} \cline{9-10}
\multicolumn{2}{|c|}{} & & 1 & 5 & 1 & 5 & & 1 & 5 \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}AIDS\\$\scaleobj{0.7}{\begin{matrix}\zeta(t)=&\hspace{-10.5em}\frac{2}{3}\left(1+\tanh(0.15t)\right)\and\\X=&
5 \log(\zeta(Z)) + 10\log(\zeta(\upeta) ) - 26.25\log^2(\zeta(\upeta))\and + 3.25\log(\zeta(Z))\log(\zeta(\upeta)) - 0.15 \sqrt{\zeta(Z)\zeta(\upeta)}\and\\Y=&\hspace{-12.5em}8X+6\epsilon-X\epsilon
\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.022} & \textbf{0.023} & \textbf{0.213} & \textbf{0.091} & \textbf{0.458} & \multirow{2}{*}{\textbf{326.43}} & \multirow{2}{*}{\textbf{653.75}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.006, 0.047{]} & {[}0.006, 0.047{]} & {[}0.169, 0.258{]} & {[}0.043, 0.158{]} & {[}0.368, 0.545{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.016} & \textbf{0.017} & \textbf{0.184} & \textbf{0.082} & \textbf{0.405} & \multirow{2}{*}{\textbf{281.24}} & \multirow{2}{*}{\textbf{578.71}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.004, 0.032{]} & {[}0.004, 0.035{]} & {[}0.149, 0.224{]} & {[}0.038, 0.146{]} & {[}0.331, 0.483{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.006} & \textbf{0.007} & \textbf{0.084} & \textbf{0.037} & \textbf{0.341} & \multirow{2}{*}{\textbf{501.19}} & \multirow{2}{*}{\textbf{486.42}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.002, 0.012{]} & {[}0.002, 0.015{]} & {[}0.067, 0.102{]} & {[}0.017, 0.064{]} & {[}0.277, 0.403{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Backward-bending supply\\$\scaleobj{0.8}{\begin{matrix}X=&
\exp(Z\upeta-Z)-3(Z\upeta-Z)\and\\Y=&0.5X+0.6\epsilon+0.1X\epsilon\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.036} & \textbf{0.040} & \textbf{0.060} & \textbf{0.048} & \textbf{0.427} & \multirow{2}{*}{\textbf{681.90}} & \multirow{2}{*}{\textbf{600.93}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.01, 0.069{]} & {[}0.011, 0.078{]} & {[}0.031, 0.098{]} & {[}0.023, 0.079{]} & {[}0.344, 0.513{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.026} & \textbf{0.028} & \textbf{0.053} & \textbf{0.041} & \textbf{0.380} & \multirow{2}{*}{\textbf{581.20}} & \multirow{2}{*}{\textbf{778.79}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.007, 0.05{]} & {[}0.009, 0.055{]} & {[}0.03, 0.084{]} & {[}0.02, 0.069{]} & {[}0.307, 0.46{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.010} & \textbf{0.011} & \textbf{0.024} & \textbf{0.019} & \textbf{0.321} & \multirow{2}{*}{\textbf{711.83}} & \multirow{2}{*}{\textbf{720.82}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.003, 0.02{]} & {[}0.003, 0.02{]} & {[}0.013, 0.039{]} & {[}0.01, 0.032{]} & {[}0.254, 0.387{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}CES\\$\scaleobj{0.8}{\begin{matrix}\zeta(t)=&\hspace{-4.5em}\frac{2}{3}\left(1+\tanh(0.15t)\right)\and\\X=&
6\left[\frac{1}{2} \zeta(Z)^{-0.5}+\frac{1}{2} \zeta(\upeta)^{-0.5}\right]^{-1}\and\\Y=&\hspace{-4.5em}0.5X+0.6\epsilon+0.1X\epsilon\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.026} & \textbf{0.023} & \textbf{0.208} & \textbf{0.132} & \textbf{0.476} & \multirow{2}{*}{\textbf{835.82}} & \multirow{2}{*}{\textbf{1,064.11}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.008, 0.051{]} & {[}0.006, 0.045{]} & {[}0.16, 0.263{]} & {[}0.092, 0.191{]} & {[}0.381, 0.578{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.019} & \textbf{0.017} & \textbf{0.181} & \textbf{0.117} & \textbf{0.423} & \multirow{2}{*}{\textbf{721.06}} & \multirow{2}{*}{\textbf{1,231.48}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.006, 0.039{]} & {[}0.004, 0.033{]} & {[}0.145, 0.228{]} & {[}0.079, 0.172{]} & {[}0.347, 0.509{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.007} & \textbf{0.006} & \textbf{0.083} & \textbf{0.053} & \textbf{0.360} & \multirow{2}{*}{\textbf{766.22}} & \multirow{2}{*}{\textbf{1,047.15}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.002, 0.015{]} & {[}0.002, 0.013{]} & {[}0.065, 0.105{]} & {[}0.037, 0.079{]} & {[}0.288, 0.432{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Partially linear\\$\scaleobj{0.8}{\begin{matrix}X=&\hspace{-3em}Z\upeta+\upeta\and\\Y=&8X+6\epsilon-X\epsilon\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.033} & \textbf{0.032} & \textbf{0.383} & \textbf{0.201} & \textbf{0.704} & \multirow{2}{*}{\textbf{1,008.29}} & \multirow{2}{*}{\textbf{1,496.78}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.012, 0.061{]} & {[}0.01, 0.066{]} & {[}0.299, 0.46{]} & {[}0.159, 0.258{]} & {[}0.569, 0.829{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.023} & \textbf{0.023} & \textbf{0.328} & \textbf{0.178} & \textbf{0.615} & \multirow{2}{*}{\textbf{1,040.48}} & \multirow{2}{*}{\textbf{1,618.44}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.008, 0.044{]} & {[}0.007, 0.048{]} & {[}0.259, 0.396{]} & {[}0.138, 0.222{]} & {[}0.514, 0.734{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.009} & \textbf{0.009} & \textbf{0.148} & \textbf{0.082} & \textbf{0.517} & \multirow{2}{*}{\textbf{957.03}} & \multirow{2}{*}{\textbf{1,471.66}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.003, 0.016{]} & {[}0.003, 0.021{]} & {[}0.122, 0.18{]} & {[}0.061, 0.101{]} & {[}0.431, 0.62{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Cobb-Douglas\\$\scaleobj{0.8}{\begin{matrix}X=&\hspace{-5em}Z\upeta\and\\Y=&8X+6\epsilon-X\epsilon\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.038} & \textbf{0.028} & \textbf{0.335} & \textbf{0.083} & \textbf{0.643} & \multirow{2}{*}{\textbf{977.12}} & \multirow{2}{*}{\textbf{1,448.76}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.013, 0.086{]} & {[}0.008, 0.062{]} & {[}0.266, 0.413{]} & {[}0.049, 0.136{]} & {[}0.53, 0.769{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.028} & \textbf{0.020} & \textbf{0.298} & \textbf{0.073} & \textbf{0.585} & \multirow{2}{*}{\textbf{880.17}} & \multirow{2}{*}{\textbf{1,446.46}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.009, 0.059{]} & {[}0.005, 0.041{]} & {[}0.239, 0.362{]} & {[}0.042, 0.114{]} & {[}0.479, 0.689{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.011} & \textbf{0.008} & \textbf{0.134} & \textbf{0.033} & \textbf{0.488} & \multirow{2}{*}{\textbf{808.90}} & \multirow{2}{*}{\textbf{1,024.17}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.003, 0.022{]} & {[}0.002, 0.017{]} & {[}0.107, 0.162{]} & {[}0.019, 0.052{]} & {[}0.4, 0.579{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Parabolic-$\tanh$\\$\scaleobj{0.8}{\begin{matrix}X=&
6\tanh(\frac{1}{4}Z-\frac{1}{2}\upeta)+5\tanh(\frac{1}{2}\upeta)\and\\Y=&\hspace{-8.8em}\frac{1}{2}(X-\epsilon)^2
\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.021} & \textbf{0.016} & \textbf{0.105} & \textbf{0.082} & \textbf{0.332} & \multirow{2}{*}{\textbf{383.80}} & \multirow{2}{*}{\textbf{671.47}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.005, 0.044{]} & {[}0.004, 0.032{]} & {[}0.075, 0.137{]} & {[}0.05, 0.11{]} & {[}0.263, 0.397{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.015} & \textbf{0.011} & \textbf{0.092} & \textbf{0.072} & \textbf{0.297} & \multirow{2}{*}{\textbf{631.33}} & \multirow{2}{*}{\textbf{683.52}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.004, 0.032{]} & {[}0.003, 0.023{]} & {[}0.068, 0.121{]} & {[}0.042, 0.096{]} & {[}0.24, 0.356{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.006} & \textbf{0.005} & \textbf{0.042} & \textbf{0.033} & \textbf{0.252} & \multirow{2}{*}{\textbf{403.78}} & \multirow{2}{*}{\textbf{787.55}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.002, 0.013{]} & {[}0.001, 0.009{]} & {[}0.03, 0.055{]} & {[}0.02, 0.044{]} & {[}0.205, 0.302{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Translog\\$\scaleobj{0.7}{\begin{matrix}\zeta(t)=&\hspace{-10.5em}\frac{2}{3}\left(1+\tanh(0.15t)\right)\and\\X=&
5 \log(\zeta(Z)) + 10\log(\zeta(\upeta) ) - 26.25\log^2(\zeta(\upeta))\and + 3.25\log(\zeta(Z))\log(\zeta(\upeta))\and\\Y=&\hspace{-12.5em}8X+6\epsilon-X\epsilon
\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.056} & \textbf{0.038} & \textbf{0.205} & \textbf{0.122} & \textbf{0.568} & \multirow{2}{*}{\textbf{805.88}} & \multirow{2}{*}{\textbf{1,117.95}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.021, 0.1{]} & {[}0.01, 0.081{]} & {[}0.152, 0.271{]} & {[}0.072, 0.187{]} & {[}0.457, 0.673{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.041} & \textbf{0.026} & \textbf{0.182} & \textbf{0.104} & \textbf{0.513} & \multirow{2}{*}{\textbf{624.76}} & \multirow{2}{*}{\textbf{814.52}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.015, 0.077{]} & {[}0.006, 0.053{]} & {[}0.133, 0.235{]} & {[}0.061, 0.163{]} & {[}0.42, 0.601{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.016} & \textbf{0.011} & \textbf{0.082} & \textbf{0.048} & \textbf{0.431} & \multirow{2}{*}{\textbf{688.18}} & \multirow{2}{*}{\textbf{1,168.71}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.006, 0.029{]} & {[}0.003, 0.022{]} & {[}0.062, 0.109{]} & {[}0.028, 0.074{]} & {[}0.349, 0.511{]} & & \\ \hline
\end{tabular}
\end{threeparttable} }
\end{table}
\begin{table}[H]
\centering\caption{\label{tab:JSD:two}Jensen-Shannon Divergence (JSD) \\ Neural Network architecture\\
$[24, 8, 4, 2]$}\vspace{1mm}
\scaleobj{0.7}{\begin{threeparttable}\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
\hline
\multicolumn{2}{|c|}{\multirow{2}{*}{Model}} & \multirow{2}{*}{N} & \multicolumn{2}{c|}{\begin{tabular}[c]{@{}c@{}}Observed\\ CONGAN\end{tabular}} & \multicolumn{2}{c|}{\begin{tabular}[c]{@{}c@{}}Counterfactual\\ CONGAN\end{tabular}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Imbens and Newey\\ (2009)\end{tabular}} & \multicolumn{2}{c|}{\begin{tabular}[c]{@{}c@{}}No. of\\ iterations\end{tabular}} \\ \cline{4-7} \cline{9-10}
\multicolumn{2}{|c|}{} & & 1 & 5 & 1 & 5 & & 1 & 5 \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}AIDS\\$\scaleobj{0.7}{\begin{matrix}\zeta(t)=&\hspace{-10.5em}\frac{2}{3}\left(1+\tanh(0.15t)\right)\and\\X=&5 \log(\zeta(Z)) + 10\log(\zeta(\upeta) ) - 26.25\log^2(\zeta(\upeta))\and + 3.25\log(\zeta(Z))\log(\zeta(\upeta)) - 0.15 \sqrt{\zeta(Z)\zeta(\upeta)}\and\\Y=&\hspace{-12.5em}8X+6\epsilon-X\epsilon\end{matrix}}$\end{tabular}}} & \multirow{2}{*}{1000} & \textbf{0.034} & \textbf{0.030} & \textbf{0.218} & \textbf{0.106} & \textbf{0.451} & \multirow{2}{*}{\textbf{327.04}} & \multirow{2}{*}{\textbf{600.84}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.009, 0.063{]} & {[}0.008, 0.058{]} & {[}0.167, 0.285{]} & {[}0.08, 0.144{]} & {[}0.369, 0.547{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.024} & \textbf{0.008} & \textbf{0.193} & \textbf{0.083} & \textbf{0.408} & \multirow{2}{*}{\textbf{335.69}} & \multirow{2}{*}{\textbf{522.27}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.006, 0.049{]} & {[}0.002, 0.015{]} & {[}0.147, 0.253{]} & {[}0.06, 0.103{]} & {[}0.329, 0.484{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.010} & \textbf{0.008} & \textbf{0.087} & \textbf{0.042} & \textbf{0.339} & \multirow{2}{*}{\textbf{264.72}} & \multirow{2}{*}{\textbf{488.16}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.003, 0.018{]} & {[}0.002, 0.016{]} & {[}0.063, 0.111{]} & {[}0.031, 0.057{]} & {[}0.276, 0.402{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Backward-bending supply\\$\scaleobj{0.8}{\begin{matrix}X=&\exp(Z\upeta-Z)-3(Z\upeta-Z)\and\\Y=&0.5X+0.6\epsilon+0.1X\epsilon\end{matrix}}$\end{tabular}
}} & \multirow{2}{*}{1000} & \textbf{0.036} & \textbf{0.036} & \textbf{0.057} & \textbf{0.036} & \textbf{0.423} & \multirow{2}{*}{\textbf{516.95}} & \multirow{2}{*}{\textbf{580.35}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.011, 0.073{]} & {[}0.015, 0.07{]} & {[}0.028, 0.093{]} & {[}0.015, 0.074{]} & {[}0.343, 0.507{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.025} & \textbf{0.025} & \textbf{0.050} & \textbf{0.032} & \textbf{0.375} & \multirow{2}{*}{\textbf{624.06}} & \multirow{2}{*}{\textbf{656.83}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.007, 0.051{]} & {[}0.01, 0.053{]} & {[}0.025, 0.082{]} & {[}0.013, 0.061{]} & {[}0.309, 0.451{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.010} & \textbf{0.010} & \textbf{0.023} & \textbf{0.014} & \textbf{0.322} & \multirow{2}{*}{\textbf{427.53}} & \multirow{2}{*}{\textbf{774.46}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.003, 0.021{]} & {[}0.004, 0.019{]} & {[}0.012, 0.037{]} & {[}0.006, 0.026{]} & {[}0.258, 0.385{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}CES\\$\scaleobj{0.8}{\begin{matrix}\zeta(t)=&\hspace{-4.5em}\frac{2}{3}\left(1+\tanh(0.15t)\right)\and\\X=&6\left[\frac{1}{2} \zeta(Z)^{-0.5}+\frac{1}{2} \zeta(\upeta)^{-0.5}\right]^{-1}\and\\Y=&\hspace{-4.5em}0.5X+0.6\epsilon+0.1X\epsilon\end{matrix}}$\end{tabular}
}} & \multirow{2}{*}{1000} & \textbf{0.025} & \textbf{0.027} & \textbf{0.058} & \textbf{0.050} & \textbf{0.475} & \multirow{2}{*}{\textbf{997.71}} & \multirow{2}{*}{\textbf{894.44}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.009, 0.051{]} & {[}0.007, 0.054{]} & {[}0.034, 0.087{]} & {[}0.022, 0.096{]} & {[}0.392, 0.567{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.018} & \textbf{0.019} & \textbf{0.051} & \textbf{0.043} & \textbf{0.425} & \multirow{2}{*}{\textbf{727.49}} & \multirow{2}{*}{\textbf{1,878.86}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.007, 0.037{]} & {[}0.005, 0.039{]} & {[}0.031, 0.079{]} & {[}0.019, 0.089{]} & {[}0.345, 0.517{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.007} & \textbf{0.008} & \textbf{0.024} & \textbf{0.020} & \textbf{0.363} & \multirow{2}{*}{\textbf{599.17}} & \multirow{2}{*}{\textbf{1,070.83}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.002, 0.014{]} & {[}0.002, 0.014{]} & {[}0.014, 0.035{]} & {[}0.009, 0.036{]} & {[}0.299, 0.435{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Partially linear\\$\scaleobj{0.8}{\begin{matrix}X=&\hspace{-3em}Z\upeta+\upeta\and\\Y=&8X+6\epsilon-X\epsilon\end{matrix}}$\end{tabular}
}} & \multirow{2}{*}{1000} & \textbf{0.043} & \textbf{0.028} & \textbf{0.447} & \textbf{0.108} & \textbf{0.712} & \multirow{2}{*}{\textbf{1,046.71}} & \multirow{2}{*}{\textbf{1,709.63}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.013, 0.08{]} & {[}0.007, 0.061{]} & {[}0.358, 0.546{]} & {[}0.075, 0.149{]} & {[}0.574, 0.841{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.030} & \textbf{0.020} & \textbf{0.385} & \textbf{0.094} & \textbf{0.625} & \multirow{2}{*}{\textbf{1,233.12}} & \multirow{2}{*}{\textbf{1,550.40}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.01, 0.059{]} & {[}0.005, 0.044{]} & {[}0.307, 0.464{]} & {[}0.063, 0.124{]} & {[}0.515, 0.74{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.012} & \textbf{0.008} & \textbf{0.177} & \textbf{0.043} & \textbf{0.535} & \multirow{2}{*}{\textbf{963.59}} & \multirow{2}{*}{\textbf{1,523.77}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.004, 0.024{]} & {[}0.002, 0.017{]} & {[}0.143, 0.212{]} & {[}0.028, 0.057{]} & {[}0.44, 0.63{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Cobb-Douglas\\$\scaleobj{0.8}{\begin{matrix}X=&\hspace{-5em}Z\upeta\and\\Y=&8X+6\epsilon-X\epsilon\end{matrix}}$\end{tabular}
}} & \multirow{2}{*}{1000} & \textbf{0.045} & \textbf{0.039} & \textbf{0.418} & \textbf{0.053} & \textbf{0.650} & \multirow{2}{*}{\textbf{881.01}} & \multirow{2}{*}{\textbf{1,090.46}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.015, 0.085{]} & {[}0.011, 0.088{]} & {[}0.335, 0.504{]} & {[}0.025, 0.097{]} & {[}0.542, 0.766{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.032} & \textbf{0.026} & \textbf{0.359} & \textbf{0.046} & \textbf{0.570} & \multirow{2}{*}{\textbf{921.45}} & \multirow{2}{*}{\textbf{1,374.14}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.01, 0.06{]} & {[}0.007, 0.054{]} & {[}0.291, 0.443{]} & {[}0.022, 0.084{]} & {[}0.473, 0.676{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.013} & \textbf{0.011} & \textbf{0.166} & \textbf{0.021} & \textbf{0.489} & \multirow{2}{*}{\textbf{757.70}} & \multirow{2}{*}{\textbf{1,107.21}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.004, 0.025{]} & {[}0.003, 0.023{]} & {[}0.135, 0.199{]} & {[}0.01, 0.038{]} & {[}0.397, 0.576{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Parabolic-$\tanh$\\$\scaleobj{0.8}{\begin{matrix}X=&6\tanh(\frac{1}{4}Z-\frac{1}{2}\upeta)+5\tanh(\frac{1}{2}\upeta)\and\\Y=&\hspace{-8.8em}\frac{1}{2}(X-\epsilon)^2\end{matrix}}$\end{tabular}
}} & \multirow{2}{*}{1000} & \textbf{0.026} & \textbf{0.020} & \textbf{0.112} & \textbf{0.072} & \textbf{0.327} & \multirow{2}{*}{\textbf{476.21}} & \multirow{2}{*}{\textbf{732.32}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.009, 0.049{]} & {[}0.006, 0.039{]} & {[}0.081, 0.147{]} & {[}0.043, 0.1{]} & {[}0.266, 0.404{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.018} & \textbf{0.013} & \textbf{0.099} & \textbf{0.061} & \textbf{0.294} & \multirow{2}{*}{\textbf{624.82}} & \multirow{2}{*}{\textbf{846.94}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.006, 0.033{]} & {[}0.004, 0.027{]} & {[}0.075, 0.127{]} & {[}0.04, 0.085{]} & {[}0.237, 0.356{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.007} & \textbf{0.006} & \textbf{0.046} & \textbf{0.028} & \textbf{0.252} & \multirow{2}{*}{\textbf{389.76}} & \multirow{2}{*}{\textbf{670.56}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.002, 0.014{]} & {[}0.002, 0.011{]} & {[}0.034, 0.061{]} & {[}0.018, 0.041{]} & {[}0.198, 0.298{]} & & \\ \hline
\multicolumn{2}{|c|}{\multirow{6}{*}{\begin{tabular}[c]{@{}c@{}}Translog\\$\scaleobj{0.7}{\begin{matrix}\zeta(t)=&\hspace{-10.5em}\frac{2}{3}\left(1+\tanh(0.15t)\right)\and\\X=&5 \log(\zeta(Z)) + 10\log(\zeta(\upeta) ) - 26.25\log^2(\zeta(\upeta))\and + 3.25\log(\zeta(Z))\log(\zeta(\upeta))\and\\Y=&\hspace{-12.5em}8X+6\epsilon-X\epsilon\end{matrix}}$\end{tabular}
}} & \multirow{2}{*}{1000} & \textbf{0.067} & \textbf{0.030} & \textbf{0.199} & \textbf{0.148} & \textbf{0.575} & \multirow{2}{*}{\textbf{843.40}} & \multirow{2}{*}{\textbf{1,032.88}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.021, 0.135{]} & {[}0.008, 0.075{]} & {[}0.144, 0.263{]} & {[}0.095, 0.211{]} & {[}0.469, 0.686{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{5000} & \textbf{0.047} & \textbf{0.021} & \textbf{0.172} & \textbf{0.129} & \textbf{0.508} & \multirow{2}{*}{\textbf{720.56}} & \multirow{2}{*}{\textbf{948.51}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.015, 0.099{]} & {[}0.005, 0.054{]} & {[}0.126, 0.225{]} & {[}0.083, 0.179{]} & {[}0.414, 0.602{]} & & \\ \cline{3-10}
\multicolumn{2}{|c|}{} & \multirow{2}{*}{100000} & \textbf{0.019} & \textbf{0.008} & \textbf{0.078} & \textbf{0.060} & \textbf{0.428} & \multirow{2}{*}{\textbf{712.45}} & \multirow{2}{*}{\textbf{886.32}} \\ \cline{4-8}
\multicolumn{2}{|c|}{} & & {[}0.006, 0.04{]} & {[}0.002, 0.022{]} & {[}0.056, 0.1{]} & {[}0.039, 0.081{]} & {[}0.351, 0.505{]} & & \\ \hline
\end{tabular}
\end{threeparttable} }
\end{table}
\end{comment}
\section{Conclusion}\label{Sec: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 \textit{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.
\pagebreak