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.
69,253 characters
Gaussian Transforms Modeling and the Estimation of Distributional Regression Functions
\title{Gaussian Transforms Modeling and the Estimation of Distributional
Regression Functions}
\author{Richard H. Spady$^\dag$ and Sami Stouli$^\S$}
\begin{abstract}
We propose flexible Gaussian representations for conditional cumulative
distribution functions and give a concave likelihood criterion for
their estimation. Optimal representations satisfy the monotonicity
property of conditional cumulative distribution functions, including
in finite samples and under general misspecification. We use these
representations to provide a unified framework for the flexible Maximum
Likelihood estimation of conditional density, cumulative distribution,
and quantile functions at parametric rate. Our formulation yields
substantial simplifications and finite sample improvements over related
methods. An empirical application to the gender wage gap in the United
States illustrates our framework.
\end{abstract}
\maketitle
\textsc{\small{}Keywords:}{\small{} Conditional density estimation,
conditional distributions, conditional quantiles, maximum likelihood,
misspecification, monotonicity, convexity, gender wage gap.}{\small\par}
\section{Introduction}
The flexible modeling and estimation of conditional distributions
are important for the analysis of various econometric and statistical
problems. Conditional probability density functions (PDF) are used
in the program evaluation literature where the generalized propensity
score takes the form of a conditional density (e.g., \citet{Imbens:2000}).
Conditional cumulative distribution functions (CDF) are core building
blocks in the estimation of nonseparable models with endogeneity where
the control variable takes the form of a conditional CDF (e.g., \citet{Imbens Newey 2009},
\citet{CFNSV}). Both conditional PDFs and CDFs are key components
in counterfactual analysis (e.g., \citet{DiNardoetal}, \citet{Chernozhukov Fern Melly 2013}).
For a continuous outcome variable $Y$ and a vector of explanatory
variables $X$, difficulties arise in the formulation of a flexible
model and in the choice of a loss function for the estimation of the
conditional PDF or CDF of $Y$ given $X$. First, a flexible class
of models may include elements that do not satisfy the defining properties
of PDFs (being positive and integration to one) and CDFs (monotonicity
and range in the unit interval) at each value of $X$. This allows
the chosen loss function to select a model that does not satisfy these
properties in finite samples and/or under misspecification.\footnote{These problems have motivated the development of a variety of post-estimation
methods to correct initial estimators (e.g., \citet{HM:2003}, \citet{GHI:2003},
\citet{Chernozhukov Fern Gali 2010}, \citet{HLL:2020}).} Second, maximum likelihood (ML) formulations are often difficult
to implement because nonconcave and/or unbounded likelihoods naturally
arise in the context of flexible modeling, a fact that has motivated
the development of alternatives to ML.\footnote{\citet[Chapter 1.7.4]{Vapnik:1999} motivates the development of alternatives
to ML by the ``narrowness of the ML method'', which is illustrated
by the unbounded likelihood arising for Gaussian mixture models. The
robust approach of \citet{Huber:1981} provides a convex alternative
to the use of Gaussian likelihoods for the concomitant estimation
of location and scale parameters (cf. \citet{Owen:2007} and \citet{SS:2018b}
for a discussion). Section \ref{subsec:Discussion} further elaborates
on this issue.} Third, alternative nonparametric approaches such as kernel regression
suffer from the curse of dimensionality, a limitation that restricts
their use in practice.
One approach to address these difficulties is to specify a flexible
class of models that satisfy the defining properties of PDFs and CDFs
by construction. Sieve ML approaches specify conditional PDFs as scaled
positive transformations of linear combinations of known functions
of $Y$ and $X$, and use the implied log-likelihood function for
estimation (e.g., \citet{Stone:1994}). Support vector methods specify
conditional PDFs as weighted averages of nonnegative kernel functions
of $Y$ and $X$, with weights restricted to be nonnegative and selected
so that the implied representation is as close as possible to the
empirical data distribution (\citet{Vapnik:1999}). Another approach
is to focus on flexible conditional CDF modeling, while discarding
the monotonicity requirement. Distribution regression (\citet{ForesiPerrachi}, \citet{Chernozhukov Fern Melly 2013})
specifies each level of the CDF of $Y$ given $X$ as a known CDF
transformation of a linear combination of the components of $X$.
The conditional CDF is then estimated at each $Y$ value by a sequence
of binary outcome ML estimators.
In this paper we take a different approach by formulating a flexible
class of Gaussian representations for conditional CDFs, instead of
modeling conditional PDFs or CDFs directly. We first expand the range
of the conditional CDF to the real line by application of a Gaussian
quantile transform, and then specify the resulting object as a linear
combination of known functions of $Y$ and $X$. The coefficients
in this linear combination are characterized by a globally concave
likelihood-based loss function that discards nonmonotone models within
the specified class, and selects the optimal element according to
the Kullback-Leibler Information Criterion (KLIC; \citet{Akaike:1973},
\citet{White:1982}). Implied models for conditional PDFs and CDFs
are flexible and can be estimated efficiently by ML under correct
specification. We provide estimation and inference results allowing
for misspecification, and we derive a dual formulation for our estimator
that we use for implementation.
We make four main contributions to the existing literature. First,
we reformulate the problems of modeling and estimating conditional
PDFs and CDFs in terms of conditional CDF representations with range
the real line. This enables us to specify a flexible class of linear
models for these representations.
Compared to sieve ML formulations,
conditional PDF models implied by our representations preserve flexibility
while avoiding the need for a scaling factor that in general cannot
be calculated in closed-form.
Compared to support vector methods that restrict
model coefficients to be nonnegative, our ML approach avoids discarding
potentially more accurate valid approximations (\citet{WGSVVW:1999}).
Second, we give an information-theoretic criterion for the selection
of a globally monotone model within the class of Gaussian representations
in linear form. Under general misspecification, this formulation characterizes
quasi-Gaussian representations that correspond to well-defined KLIC
optimal conditional PDF and CDF approximations to the true data probability
distribution. These approximations satisfy the defining properties
of PDFs and CDFs at each value of $X$, both in finite samples with
probability approaching one and in the population. Compared to the
distribution regression criterion, our approach allows for the global
approximation of conditional CDFs and affords KLIC optimality for
the selected model.
Compared to sieve ML formulations that specify positive conditional
PDFs, our approach relies on the loss function to rule out conditional
PDFs that are negative or zero with positive probability.
Third, we provide a unified approach to the flexible ML estimation
of conditional PDFs and CDFs at parametric rate. Flexibility is obtained
without the inclusion of infinite-dimensional parameters, in this
way going beyond distribution and quantile regression (\citet{Koenker:Bassett1978}),
as well as parametric location-scale formulations (\citet{He:1997},
\citet[Section 2]{SS:2018a}, \citet{MachadoSantosSilva:2019}). Function-valued
parameters lead to slower than root-$n$ conditional PDF estimation
(\citet{RotheWied:2020}), and location-scale models restrict the
shape of conditional distributions. Our approach alleviates the curse
of dimensionality and hence also provides an alternative to kernel-based
methods for the nonparametric estimation of conditional distributions
(e.g., \citet{LR:2007}). Our approach further applies to conditional
quantile function (CQF) estimation and is thus also related to the
kernel-based estimator of \citet{Matzkin:2003} for nonseparable models.
In contrast, we propose a parametric formulation for their flexible
ML estimation, which also facilitates the imposition of shape constraints
from economic theory.
Fourth, we derive a dual formulation for our ML estimator that demonstrates
considerable computational benefits. Compared to distribution and
quantile regression, the dual formulation provides a convex programming
problem (\citet{BV:2004}) for the one-step estimation of conditional
PDFs and CDFs at each sample points. Compared to dual regression and
its generalization (\citet{SS:2018a}), we find that
the dual formulation has the important advantage of being a mathematical
programming problem with linear constraints.
Taken together, these contributions define a novel, unified framework
for flexible distributional regression analysis. Our framework yields
substantial simplifications over related methods, while providing
all the theoretical guarantees afforded by ML. In numerical simulations,
we find that these features translate into largely improved finite
sample performance compared to kernel methods, the Matzkin estimator,
distribution and quantile regression. A distributional analysis of
the gender wage gap in the United States illustrates the benefits
of our methods for empirical practice.
Our methods extend to a variety of settings, including Logistic transform regression models,
mixed discrete-continuous outcome distributions, and multiple outcomes.
Section \ref{sec:Section2} introduces our modeling framework.
Section \ref{sec:Section3} gives results under misspecification.
Section \ref{sec:Section4} contains estimation and inference results,
and duality theory is derived in Section \ref{sec:Section5}. Section
\ref{sec:Section6} illustrates our methods, Section \ref{sec:Section7}
gives extensions and Section \ref{sec:Section8} concludes. Proofs
of Theorems \ref{thm:Thm1}-\ref{thm:Thm2} and \ref{thm:Thm6} are
given in the Appendix. The Supplemental Material (\citet{SS:2025}) contains proofs
of Theorems \ref{thm:Thm3}-\ref{thm:Thm5}, technical results,
implementation details, and results of numerical simulations.
\section{Gaussian Transforms Modeling\label{sec:Section2}}
Let $Y$ be a continuous outcome variable and $X$ a vector of explanatory
variables. A transformation to Gaussianity of the conditional CDF
$F_{Y|X}(Y|X)$ of $Y$ given $X$ occurs by application of the Gaussian
quantile function $\Phi^{-1}$,
\begin{equation}
e=\Phi^{-1}(F_{Y\mid X}(Y\mid X))\equiv g(Y,X),\label{eq:eyrep}
\end{equation}
where the resulting Gaussian Transform (GT) $e$ is a zero mean and
unit variance Gaussian random variable and is independent from $X$,
by construction. With $y\mapsto F_{Y\mid X}(y|X)$ strictly increasing,
the corresponding map $y\mapsto g(y,X)$ is also strictly increasing,
with well-defined inverse denoted $e\mapsto g^{-1}(e,X)$.
Important statistical objects such as conditional PDFs, CDFs and CQFs
can be expressed as known functionals of $g(Y,X)$. The conditional
PDF of $Y$ given $X$ is
\[
f_{Y\mid X}(Y\mid X)=\phi(g(Y,X))\{\partial_{y}g(Y,X)\},\quad\partial_{y}g(Y,X)\equiv\frac{\partial g(Y,X)}{\partial y},
\]
where $e\mapsto\phi(e)$ is the Gaussian PDF and $\partial_{y}g(y,x)$
is a partial derivative, and the conditional CDF and CQF of $Y$ given
$X$ are
\[
F_{Y\mid X}(Y\mid X)=\Phi(g(Y,X)),\quad Q_{Y\mid X}(u\mid X)=g^{-1}(\Phi^{-1}(u),X),\quad u\in(0,1),
\]
respectively. The GT $g(Y,X)$ thus constitutes a natural modeling
object in the context of distributional regression models for $f_{Y|X}(Y|X)$,
$F_{Y|X}(Y|X)$, and $Q_{Y|X}(u|X)$. We refer to these objects as
the `Distributional Regression Functions' (DRF).
In this paper we consider the class of conditional CDFs with Gaussian
representation $e=g(Y,X)$ in linear form, where $g(Y,X)$ is specified
as a linear combination of known transformations of $Y$ and $X$.
Expanding the range of conditional CDFs from the unit interval to
the real line allows for the formulation of linear yet flexible representations.
The implied models for DRFs are parsimonious
and able to capture complex features of the entire statistical relationship
between $Y$ and $X$. In particular, these models allow for nonlinearity
and nonseparability of this relationship.
\subsection{Gaussian representations in linear form}
Let $W(X)$ be a $K\times1$ vector of known functions of $X$ and
$S(Y)$ a $J\times1$ vector of known functions of $Y$, and write
$\otimes$ for their Kronecker product. Assume that $W(X)$ includes
an intercept, i.e., has first component $1$, and that $S(Y)$ has
first two components $(1,Y)'$ and derivative $dS(Y)/dy=s(Y)$, a
vector of functions continuous on $\mathbb{R}$. Given a random vector
$(Y,X')'$ with support $\mathcal{YX}=\mathcal{Y}\times\mathcal{X}$,
where $\mathcal{Y}=\mathbb{R}$ and $\mathcal{Y}$ and $\mathcal{X}$
are the marginal supports of $Y$ and $X$, respectively, our first
assumption defines a GT regression model.
\begin{assumption}For some $b_{0}\in\mathbb{R}^{JK}$, a GT regression
model takes the form
\begin{equation}
e=b_{0}'T(X,Y),\quad\partial_{y}\{b_{0}'T(X,Y)\}=b_{0}'t(X,Y)>0,\quad e\mid X\sim N(0,1),\label{eq:e}
\end{equation}
where $T(X,Y)\equiv W(X)\otimes S(Y)$ and $t(X,Y)\equiv W(X)\otimes s(Y)$.\label{ass:Ass1}\end{assumption}
The GT $g(Y,X)$ in (\ref{eq:eyrep}) is specified as a linear combination
of the known functions $T(X,Y)$. The linear form of $e$ is preserved
by the derivative function $b_{0}'t(X,Y)$ which is simultaneously
specified as a linear combination of $t(X,Y)$. For bounded $W(X)$
and $s(Y)$ with good approximation properties such as splines or wavelets,
this linear specification is a flexible model for (\ref{eq:eyrep}).\footnote{If some component of $W(X)$ or $s(Y)$ has range the real line over
$\mathcal{X}$ or $\mathcal{Y}$, respectively, then its coefficient
in (2.2) must be zero since there is no $b_{0}\neq0$ such that $b_{0}'t(X,Y)>0$
in that case. For vector $X$, an example of a flexible structure
takes $W(X)=W_{1}(X_{1})\varotimes\cdots\varotimes W_{\dim(X)}(X_{\dim(X)})$,
where $W_{l}(X_{l})$ includes an intercept and a vector of bounded
approximating functions, $l=1,\ldots,\dim(X)$.}$^{,}$\footnote{
When $b_0'T(X,Y)$ is viewed as an approximation, condition $b_0't(X,Y) > 0$
in (\ref{eq:e}) restricts the set of admissible approximating functions,
without compromising approximation quality under regularity conditions and for $W(X)$ rich enough.
For $y \mapsto g(y,x)$ uniformly smooth in $x$ and with $S(Y)$
chosen from a suitable class of approximating functions,
there is $\beta(X)$ with $\beta(X)'s(Y)>0$ and uniform approximation error
of $\beta(X)'S(Y)$ for $g(Y,X)$ of the same order as that of an unconstrained approximation.
This is because the approximation errors for $y\mapsto g(y,X)$ of
each approximation type are of the same order for each $x$
(e.g., \citet{DeVore:1977} for splines and \citet{AY:1992} for wavelets),
and hence uniformly over $x$ if $y\mapsto g(y,x)$ is smooth uniformly in $x$.
Moreover, if $W(X)$ is mean-square spanning
and $S(Y)$ has finite conditional second moment,
then there is $b$ such that $E[\{ b'T(X,Y) - \beta(X)'S(Y) \}^2] \rightarrow 0$ as $K \rightarrow \infty$.
Therefore, for $K$ large enough, $b't(X,Y)>0$, by $\beta(X)'s(Y)>0$, and the orders of constrained and
unconstrained approximations being the same
now implies that taking $b'T(X,Y)$ with $b't(X,Y)>0$ in (\ref{eq:e}) preserves approximation quality.
}
When the nonconstant components of $W(X)$ and $s(Y)$ are specified
as nonnegative spline functions (\citet{CS:1966}, \citet{Ramsay:1988}),
we refer to the implied representations as `Spline-Spline models'.
We note that we do not impose $b_{0}>0$.
Model (\ref{eq:e}) corresponds to a well-defined probability distribution
for $Y$ given $X$ with GT in linear form. The implied forms of the
conditional PDF and CDF are
\begin{equation}
f_{Y\mid X}(Y\mid X)=\phi(b_{0}'T(X,Y))\{b_{0}'t(X,Y)\},\quad F_{Y\mid X}(Y\mid X)=\Phi(b_{0}'T(X,Y)).\label{eq:DRF1}
\end{equation}
The resulting log conditional density forms the basis of our approach:
\begin{equation}
\log f_{Y\mid X}(Y\mid X)=-\frac{1}{2}[\log(2\pi)+\{b_{0}'T(X,Y)\}^{2}]+\log(b_{0}'t(X,Y)).\label{eq:log PDF}
\end{equation}
Viewed as an equation in the GT $b_{0}'T(X,Y)$, (\ref{eq:log PDF})
defines
the problem of characterizing GTs of the form (\ref{eq:e}), ruling
out implied conditional PDFs that are negative or zero with positive
probability, and hence also nonmonotone conditional CDFs.
Formulation (\ref{eq:log PDF})
differs from sieve ML approaches
that specify conditional PDFs as $f_{Y|X}(Y|X)=\varphi(b_{0}'T(X,Y))/\int\varphi(b_{0}'T(X,y))dy$
with $\varphi(\cdot)$ a nonnegative function such as the exponential
or the square functions. The log conditional PDF is
\begin{equation}
\log f_{Y\mid X}(Y\mid X)=\log\varphi(b_{0}'T(X,Y))-\log\int\varphi(b_{0}'T(X,y))dy,\label{eq:sieve ML}
\end{equation}
with scaling factor $\int\varphi(b_{0}'T(X,y))dy$ not available in
closed-form in general, and with implied ML first-order conditions
that are nonlinear in $b_{0}'T(X,Y)$. In contrast, (\ref{eq:log PDF})
avoids the scaling factor and yields first-order conditions linear
in $b_{0}'T(X,Y)$ (cf. (\ref{eq:FOCs}) and (\ref{eq:MMrep}) below).
Formulation (\ref{eq:log PDF}) also
differs from support vector methods that define the problem
of characterizing conditional PDFs as solving the equation
\begin{equation}
F_{YX}(y,x)=\int_{-\infty}^{x}\int_{-\infty}^{y}f(s\mid t)dF_{X}(t)ds\label{eq:svm int eq}
\end{equation}
for $f$, where $f$ is specified as a linear combination of nonnegative
kernel functions of $Y$ and $X$ with nonnegative coefficients, which
is overly restrictive (cf. Section 5 in \citet{WGSVVW:1999}).
Compared to both (\ref{eq:sieve ML}) and (\ref{eq:svm int eq}),
shape constraints from economic theory are also easier to impose using
(\ref{eq:log PDF}), with GT in closed-form.\footnote{Shape constraints often apply to CDFs or CQFs (e.g., \citet{BKM:2014},
\citet{CW:2017}), and they easily translate into restrictions on
the shape of GTs. For example, under Assumption \ref{ass:Ass1}, a
nonseparable demand model $Y=h(X,e)$ is nonincreasing in prices
$X$ if the linear constraints $\partial_{x}t(X,Y)'b_{0}\geq0$ hold
in (\ref{eq:log PDF}). By $\partial_{x}Q_{Y|X}(u|X)=-\left.\partial_{x}F_{Y|X}(y_{0}|X)/\partial_{y}F_{Y|X}(y_{0}|X)\right|_{y_{0}=Q_{Y|X}(u|X)}$
and $\partial_{y}F_{Y|X}(Y|X)>0$, nonincreasing demand is implied
by $\partial_{x}F_{Y|X}(Y|X)=\phi(b_{0}'T(X,Y))\{\partial_{x}t(X,Y)'b_{0}\}\geq0$,
and hence by $\partial_{x}t(X,Y)'b_{0}\geq0$.}
\begin{rem}
For $W(X)=1$, models for marginal PDF and CDF of $Y$ arise as a
particular case of (\ref{eq:e}).
For $J=2$, the Jacobian term $b_{0}'t(X,Y)$ does not depend on $Y$ which restricts $F_{Y\mid X}(Y|X)$
to Gaussianity at all values of $X$.
\begin{rem}
Our modeling framework also applies when $\mathcal{Y}$ is bounded
since $Y$ can always be monotonically transformed to a random variable
with support expanded over the real line (cf. Remark \ref{Rk:Boundedsupp} in Section
\ref{sec:Implementation} of the Supplemental Material).
\end{rem}
\end{rem}
\subsection{Characterization}
For $\Theta=\{b\in\mathbb{R}^{JK}:\Pr[b't(X,Y)>0]=1\}$, we define
the population objective function
\begin{equation}
Q(b)=E\left[-\frac{1}{2}\left(\log(2\pi)+\{b'T(X,Y)\}^{2}\right)+\log\left(b't(X,Y)\right)\right],\quad b\in\Theta.\label{eq:popML}
\end{equation}
This criterion introduces a natural logarithmic barrier function (e.g.,
\citet{BV:2004}) in the form of the log of the Jacobian term $b't(X,Y)$.
This is important because the monotonicity requirement for the conditional
CDF is imposed directly by the objective in the definition of the
effective domain of $Q(b)$, i.e., the region in $\mathbb{R}^{JK}$
where $Q(b)>-\infty$. An equivalent interpretation is that the effective
domain of $Q(b)$ contains the set of parameter values that are admissible
for GT regression models with positive conditional PDF, by virtue
of the presence of both the Gaussian density function and the logarithmic
barrier function in (\ref{eq:popML}).
We characterize the shape and properties of $Q(b)$ under the following
assumption.
\begin{assumption}$E[||T(X,Y)||^{2}]<\infty$, $E[||t(X,Y)||^{2}]<\infty$,
and the smallest eigenvalue of $E[T(X,Y)T(X,Y)']$ is bounded away
from zero.\label{ass: Ass2}\end{assumption}
These conditions restrict the set of dictionaries we allow for, as
well as the probability distribution of $Y$ conditional on $X$.
In particular, because $T(X,Y)$ includes $Y$, Assumption \ref{ass: Ass2}
requires $Y$ to have finite second moment. The moment conditions
in Assumption \ref{ass: Ass2} are also sufficient for the second-derivative
matrix of $Q(b)$,
\begin{equation}
\Gamma(b)\equiv E\left[\gamma(Y,X,b)\right],\quad\gamma(Y,X,b)\equiv-T(X,Y)T(X,Y)'-\frac{t(X,Y)t(X,Y)'}{\{b't(X,Y)\}^{2}},\label{eq:hessian}
\end{equation}
to exist for each $b\in\Theta$. Nonsingularity of $E[T(X,Y)T(X,Y)']$
guarantees that $\Gamma(b)$ is negative definite, and hence that
$Q(b)$ is strictly concave and has a unique maximum.
\begin{thm}
\label{thm:Thm1}If Assumptions \ref{ass:Ass1}-\ref{ass: Ass2} hold
then $Q(b)$ is strictly concave and has a unique maximum in $\Theta$
at $b_{0}$.
\end{thm}
By standard ML theory (e.g., \citet{Newey:McFadden:1994}, p. 2124),
this result implies identification of $b_{0}$, with $b_{0}$ being
the only solution to the first-order conditions
\begin{equation}
E\left[\psi(Y,X,b_{0})\right]=0,\;\psi(Y,X,b)\equiv-T(X,Y)(b'T(X,Y))+\frac{t(X,Y)}{b't(X,Y)},\;b\in\Theta.\label{eq:FOCs}
\end{equation}
When $Y|X\sim N(0,1)$, an interesting connection with the the Stein
equation for standard Gaussian random variables arises (e.g., Lemma
2.1 in \citet{Chen Gold Shao 2010}). In that case, $b_{0}'T(X,Y)=Y$
and $b_{0}'t(X,Y)=1$ satisfy the conditions of model (\ref{eq:e}).
Theorem \ref{thm:Thm1} then implies that $b_{0}=(0,1,0_{JK-2})'$
uniquely solves (\ref{eq:FOCs}):
\begin{align*}
E\left[\psi(Y,X,b_{0})\right]=E\left[-T(X,Y)Y+t(X,Y)\right] & =E[W(X)\otimes\{-S(Y)Y+s(Y)\}]\\
& =E[W(X)]\otimes E[-S(Y)Y+s(Y)]=0,
\end{align*}
since $E[-S(Y)Y+s(Y)]=0$ has the form of the Stein equation, and
hence holds for any vector of continuously differentiable functions
$S(Y)$ with $E[|s_{j}(Y)|]<\infty$, $j\in\{1,\ldots,J\}$. In contrast,
(\ref{eq:FOCs}) holding with $b_{0}\neq(0,1,0_{JK-2})'$ will indicate
deviations of $Y$ from Gaussianity and independence from $X$. Since
$b_{0}$ satisfies (\ref{eq:e}), conditions (\ref{eq:FOCs}) thus
characterize a transformation of $Y$ to Gaussianity at each $X$ value.
Hence, Theorem \ref{thm:Thm1} has the following testable implications for model (\ref{eq:e}).
\begin{cor}
\label{cor:Cor1}If there exists $b_{0}$ such that model (\ref{eq:e})
holds then, for any vectors of functions $\overline{W}(X)$ and of
continuously differentiable functions $\overline{S}(e)$ such that
$\overline{T}(X,e)\equiv\overline{W}(X)\otimes\overline{S}(e)$
and $\overline{t}(X,e)\equiv\partial_{e}\overline{T}(X,e)$ satisfy
Assumption \ref{ass: Ass2} with $T=\overline{T}$, $t=\overline{t}$
and $Y=e$, the following hold: (i) $(0,1,0_{JK-2})'$ uniquely solves
\[
\max_{b\in\overline{\Theta}}E[\log(\phi(b'\overline{T}(X,e))\{b'\overline{t}(X,e)\})],\quad\overline{\Theta}\equiv\{b\in\mathbb{R}^{JK}:\Pr[b'\overline{t}(X,e)>0]=1\},
\]
and (ii) the `Stein score' conditions $E[-\overline{T}(X,e)e+\overline{t}(X,e)]=0$
hold.
\end{cor}
\subsection{Discussion\label{subsec:Discussion}}
The general modeling of $F_{Y|X}(Y|X)$ can be done indirectly by
specifying a representation for $Y$ given $X$,
\begin{equation}
Y=H(X,e),\quad e\mid X\sim F_{e},\label{eq:e form}
\end{equation}
with $H(X,e)$ strictly increasing in $e$, a random variable with
distribution $F_{e}$ and independent of $X$. The specification of
$H$ and $F_{e}$ then determines the form of $F_{Y|X}(Y|X)$:
\begin{equation}
F_{Y\mid X}(y\mid X)=F_{e}(H^{-1}(y,X)),\quad y\in\mathbb{R},\label{eq:eform2}
\end{equation}
where $y\mapsto H^{-1}(y,X)$ denotes the inverse function of $e\mapsto H(X,e)$.
In this approach, for a specified distribution $F_{e}$ the object
of modeling is the function $H(X,e)$.
In Econometrics, (\ref{eq:e form}) is often characterized as `nonlinear
and nonseparable' in order to draw attention to the potentially complex
$Y$\textendash $X$ structure at constant $e$ and the lack of additive
structure in $e$ (e.g., \citet{Chesher:2003}, \citet{Matzkin:2003}).
These are essential features of $H$ that allow for the shape of the
conditional distribution of $Y$ to vary across values of $X$. An
alternative approach to (\ref{eq:e form})-(\ref{eq:eform2}) that
preserves nonlinearity and nonseparability is to model $F_{Y|X}(Y|X)$
directly as
\begin{equation}
F_{Y\mid X}(y\mid X)=F_{e}(g(y,X)),\quad y\in\mathbb{R},\label{eq:Y form}
\end{equation}
for some strictly increasing function $y\mapsto g(y,X)$. In our approach,
the object of modeling is the quantile transform $g(X,Y)=F_{e}^{-1}(F_{Y|X}(Y|X))$
which has distribution $F_{e}$ and is independent of $X$ by construction,
for some specified quantile function $F_{e}^{-1}$.
The difference between modeling the statistical relationship between $Y$ and $X$
according to (\ref{eq:eform2}) or (\ref{eq:Y form})
is not innocuous. With $f_{e}$ denoting the PDF of $e$, the definition
of the conditional PDF of $Y$ given $X$ implied by the indirect
approach (\ref{eq:eform2}),
\begin{equation}
f_{Y\mid X}(y\mid X)=f_{e}(H^{-1}(y,X))\{\partial_{y}H^{-1}(y,X)\},\quad y\in\mathbb{R},\label{eq:dens_eform}
\end{equation}
involves the inverse function of the modeling object $H$. In general
this inverse function does not have a closed-form expression, except
for some simple cases like the location-scale model $H(X,e)\equiv X'\beta_{1}+(X'\beta_{2})e$
with $X'\beta_{2}>0$. Furthermore, expression (\ref{eq:dens_eform})
gives rise to a nonconcave likelihood for even the simplest specifications
of $H$ and $F_{e}$, including the location and location-scale models
with Gaussian $e$ (\citet{Owen:2007}, \citet{SS:2018b}). In contrast,
a major advantage of representation (\ref{eq:Y form}) is that the
corresponding expression for $f_{Y|X}(Y|X)$ circumvents the inversion
step since
\[
f_{Y\mid X}(Y\mid X)=f_{e}(g(Y,X))\{\partial_{y}g(Y,X)\}.
\]
This formulation allows for the direct specification of flexible models
for $g(Y,X)$ that are characterized by a concave likelihood. Hence,
considerable computational advantages accrue in estimation when $e=g(Y,X)$
can be computed in closed-form, as further demonstrated by the duality
analysis in Section \ref{sec:Section5}. This formulation also leads
to well-defined representations for $F_{Y|X}(Y|X)$ under misspecification.
\section{Quasi-Gaussian Representations under Misspecification\label{sec:Section3}}
We study the properties of quasi-Gaussian representations for $F_{Y|X}(Y|X)$
that are generated by maximization of the objective $Q(b)$ under
general misspecification, i.e., when there is no representation of
the form (\ref{eq:e}) that satisfies either the Gaussianity or the
independence properties, or both. We find that the implied approximations
for the true DRFs are well-defined
and KLIC optimal.
\subsection{Existence and uniqueness}
Under Assumption \ref{ass: Ass2} the objective function $Q(b)$ is continuous and
strictly concave over $\Theta$, and hence admits at most one maximizer.
Assumption \ref{ass: Ass2} is also sufficient for the level sets of $Q(b)$ to be compact, and hence for
existence of a maximizer. Compactness of the level sets
is a consequence of the explosive behavior of $Q(b)$ at the boundary of $\Theta$.
By the quadratic term $-\{b'T(X,Y)\}^{2}$
being negative, as $b$ approaches the boundary of $\Theta$ the log
Jacobian term diverges to $-\infty$, and hence so does $-\{b'T(X,Y)\}^{2}/2+\log\{b't(X,Y)\}$
on a set with positive probability.
This is sufficient to conclude that the objective function $Q(b)$
diverges to $-\infty$, and hence that there exists at least one maximizer
to $Q(b)$ in $\Theta$, denoted $b^{*}$.
Under misspecification, to the maximizer $b^{*}$ corresponds the
quasi-Gaussian representation $e^{*}=T(X,Y)'b^{*}\equiv g^{*}(Y,X)$,
where $g^{*}(Y,X)$ is an element of the set
\[
\mathcal{E}\equiv\left\{ m:\Pr[m(Y,X)=b'T(X,Y)]=1\right\}
\]
with $b\in\Theta$. By definition of $\Theta$, $y\mapsto b'T(X,y)$
is strictly increasing for each $b\in\Theta$ with probability one,
and hence each $m\in\mathcal{E}$ has a well-defined inverse function.
We note that nonsingularity of $E[T(X,Y)T(X,Y)']$ implies that $g^{*}(Y,X)$
is unique in $\mathcal{E}$, i.e., there is no $m=g^{*}$ in $\mathcal{E}$
with $m(y,x)=b'T(x,y)$ everywhere and $b\neq b^{*}$.
\begin{sloppy}Define the range of $y\mapsto\Phi(m(y,x))$ as $\mathcal{U}_{x}(m)\equiv\{u\in(0,1):\Phi(m(y,x))=u\textrm{ for some }y\in\mathbb{R}\}$,
for $m\in\mathcal{E}$ and $x\in\mathcal{X}$. To the quasi-Gaussian
representation $g^{*}(Y,X)$ correspond a conditional PDF approximation,
defined as
\begin{equation}
f^{*}(Y,X)\equiv\phi(g^{*}(Y,X))\{\partial_{y}g^{*}(Y,X)\},\label{eq:f*(Y,X)}
\end{equation}
and conditional CDF and CQF approximations, defined as
\[
F^{*}(Y,X)\equiv\Phi(g^{*}(Y,X)),\quad Q^{*}(u,X)\equiv g^{*-1}(\Phi^{-1}(u),X),\quad u\in\mathcal{U}_{X}(g^{*}),
\]
where $e\mapsto g^{*-1}(e,X)$ denotes the inverse of $y\mapsto g^{*}(y,X)$.
These representations are unique in, respectively, the following spaces
\begin{align*}
\mathcal{D} & \equiv\left\{ f:\Pr[f(Y,X)=\phi(m(Y,X))\{\partial_{y}m(Y,X)\}]=1\right\} \\
\mathcal{F} & \equiv\left\{ F:\Pr[F(Y,X)=\Phi(m(Y,X))]=1\right\} \\
\mathcal{Q} & \equiv\left\{ Q:\Pr[Q(u,X)=m^{-1}(\Phi^{-1}(u),X)\;\textrm{for all}\;u\in\mathcal{U}_{X}(m)]=1\right\}
\end{align*}
with $m\in\mathcal{E}$, and where $e\mapsto m^{-1}(e,X)$ denotes the
inverse of $y\mapsto m(y,X)$. Therefore, the DRF approximations are well-defined, with positive conditional
PDF and monotone conditional CDF and CQF approximations.\par\end{sloppy}
\begin{thm}
\label{thm:Thm2}If Assumption \ref{ass: Ass2}
holds then there exists a unique maximum $b^{*}$ to $Q(b)$ in $\Theta$.
Consequently, the quasi-Gaussian representation $g^{*}(Y,X)$ and
the corresponding approximations for the DRFs are unique.
\end{thm}
\subsection{KLIC optimality}
When the elements of $\mathcal{D}$ are proper conditional PDFs, a
further motivation for the use of $Q(b)$ is the KLIC optimality of
the implied DRFs under misspecification (\citet{Akaike:1973}, \citet{White:1982}).
Since each $f\in\mathcal{D}$ satisfies $f>0$ by construction, an
element $f\in\mathcal{D}$ is a proper conditional PDF if it satisfies
$\int_{\mathbb{R}}f(y,X)dy=1$ with probability one. A necessary and
sufficient condition for this to hold is that the boundary conditions
\begin{equation}
\lim_{y\rightarrow-\infty}b'T(X,y)=-\infty,\quad\lim_{y\rightarrow\infty}b'T(X,y)=\infty,\label{eq:Boundary conditions}
\end{equation}
hold with probability one, for all $b\in\Theta$. Given a specified
dictionary such that (\ref{eq:Boundary conditions}) holds, Theorem
\ref{thm:Thm2} implies that the approximation $f^{*}(Y,X)$ in (\ref{eq:f*(Y,X)})
is the unique maximum selected by the population criterion in $\mathcal{D}$,
i.e.,
\[
f^{*}=\arg\max_{f\in\mathcal{D}}E\left[\log f(Y,X)\right],
\]
and hence that $f^{*}(Y,X)$ is the KLIC closest probability distribution
to $f_{Y|X}(Y|X)$. The corresponding $F^{*}$ and $Q^{*}$ are then
the KLIC optimal conditional CDF and CQF approximations for $F_{Y|X}(Y|X)$
and $Q_{Y|X}(u|X)$, respectively.
\begin{thm}
\label{thm:Thm3}If $E[|\log f_{Y|X}(Y|X)|]<\infty$ and (\ref{eq:Boundary conditions})
holds with probability one for all $b\in\Theta$, then $f^{*}$ is
the KLIC closest probability distribution to $f_{Y|X}(Y|X)$ in $\mathcal{D}$,
i.e.,
\[
f^{*}=\arg\min_{f\in\mathcal{D}}E\left[\log\left(\frac{f_{Y\mid X}(Y\mid X)}{f(Y,X)}\right)\right],
\]
where each $f\in\mathcal{D}$ is a proper conditional PDF. Moreover,
$f^{*}$ is related to the KLIC optimal conditional CDF $F^{*}$ in
$\mathcal{F}$ by
\[
F^{*}(y,X)=\int_{-\infty}^{y}f^{*}(t,X)dt,\quad y\in\mathbb{R},
\]
and to the well-defined inverse of $y\mapsto F^{*}(y,X)$, the KLIC
optimal CQF $u\mapsto Q^{*}(X,u)$ in $\mathcal{Q}$ with derivative
\[
\frac{\partial Q^{*}(X,u)}{\partial u}=\frac{1}{f^{*}(Q^{*}(X,u),X)}>0,\quad u\in(0,1),
\]
with probability one.
\end{thm}
Under the boundary conditions (\ref{eq:Boundary conditions}), the
set $\mathcal{F}$ is the space of conditional CDFs with Gaussian
representation in linear form, and the set $\mathcal{Q}$ is the space
of corresponding well-defined CQFs. A necessary and sufficient condition
for (\ref{eq:Boundary conditions}) is obtained, for instance, if
the limits $\lim_{y\rightarrow\pm\infty}|S_{j}(y)|$ are finite, $j\in\{3,\ldots,J\}$.
Under this maintained condition, the varying coefficients representation
of $e$ in (\ref{eq:e}),
\begin{equation}
e=\beta_{1}(X)+\beta_{2}(X)Y+\sum_{j=3}^{J}\beta_{j}(X)S_{j}(Y),\;\beta_{j}(X)=W(X)'b_{0j},\;j\in\{1,\ldots,J\},\label{eq:Varying}
\end{equation}
implies that $\beta_{2}(X)>0$ is necessary for the boundary conditions
(\ref{eq:Boundary conditions}) because otherwise $\lim_{y\rightarrow\infty}\beta(X)'S(y)$
would be finite or $-\infty$, and $\lim_{y\rightarrow-\infty}\beta(X)'S(y)$
would be finite or $\infty$, with $\beta(X)=(\beta_{1}(X),\ldots,\beta_{J}(X))'$.
The linear term $\beta_{2}(X)Y$ in (\ref{eq:Varying}) implies that $\beta_{2}(X)>0$
is also sufficient for (\ref{eq:Boundary conditions}). We note
that $\beta_{2}(X)>0$ is implied by the derivative condition $\beta(X)'s(Y)=\beta_{2}(X)+\sum_{j=3}^{J}\beta_{j}(X)s_{j}(Y)>0$
if the transformations $s_{j}(Y)$, $j\in\{3,\ldots,J\}$, are specified
to be zero outside some compact region of $\mathbb{R}$,\footnote{This and the maintained assumption that $\lim_{y\rightarrow\pm\infty}|S_{j}(y)|<\infty$
are satisfied for instance if, for each $j\in\{3,\ldots,J\}$, the
transformations $S_{j}(Y)$ are defined as $S_{j}(y)\equiv\int_{-\infty}^{y}s_{j}(t)dt$,
for nonnegative spline functions $s_{j}(Y)\neq0$ on a compact subset
of $\mathbb{R}$, as $s_{j}(Y)=0$ outside this region and $S_{j}(Y)$
is then a CDF over the entire real line (\citet{CS:1966}, \citet{Ramsay:1988}).} since the derivative then reduces to $\beta_{2}(X)$ outside this
region. The boundary conditions (\ref{eq:Boundary conditions}) then
effectively hold under a location-scale restriction in the tails of
the distribution of $Y$ given $X$. We also note that (\ref{eq:Boundary conditions})
always holds for $J=2$ since the derivative condition is $\beta_{2}(X)>0$
in that particular case.
\section{Estimation and Inference\label{sec:Section4}}
\subsection{Maximum Likelihood estimation}
We assume that we observe a sample of $n$ independent and identically
distributed realizations $\{(y_{i},x_{i})\}_{i=1}^{n}$ of the random
vector $(Y,X')'$. Using the sample analog of $Q(b)$, we define the
GT regression estimator
\begin{equation}
\widehat{b}\equiv\arg\max_{b\in\Theta}n^{-1}\sum_{i=1}^{n}\left\{ -\frac{1}{2}[\log(2\pi)+\{b'T(x_{i},y_{i})\}^{2}]+\log(b't(x_{i},y_{i}))\right\} .\label{eq:ML}
\end{equation}
We derive the asymptotic properties of $\hat{b}$ under the following
assumptions.
\begin{assumption}(i) $\{(y_{i},x_{i})\}_{i=1}^{n}$ are identically
and independently distributed, and (ii) $E[||T(X,Y)||^{4}]<\infty$.\label{ass:Ass4}\end{assumption}
Assumption \ref{ass:Ass4}(i) can be replaced with the condition that
$\{(y_{i},x_{i})\}_{i=1}^{n}$ is stationary and ergodic (\citet{Newey:McFadden:1994}).
Assumption \ref{ass:Ass4}(ii) is needed for consistent estimation
of the asymptotic variance-covariance matrix of $\widehat{b}$.
\begin{sloppy}Recalling the definitions of $\gamma(Y,X,b)$ and $\Gamma(b)$
in (\ref{eq:hessian}) and $\psi(Y,X,b)$ in (\ref{eq:FOCs}), the
variance-covariance matrix of $\hat{b}$ is $\Gamma^{-1}\Psi\Gamma^{-1}/n$,
where $\Gamma\equiv\Gamma(b^{*})$ and $\Psi\equiv E[\psi(Y,X,b^{*})\psi(Y,X,b^{*})']$.
Estimators of $\Gamma$ and $\Psi$ are defined as $\widehat{\Gamma}=n^{-1}\sum_{i=1}^{n}\gamma(y_{i},x_{i},\hat{b})$
and $\widehat{\Psi}=n^{-1}\sum_{i=1}^{n}\psi(y_{i},x_{i},\hat{b})\psi(y_{i},x_{i},\hat{b})'$,
respectively. An estimator of $\Gamma^{-1}$ is any symmetric generalized
inverse $\widehat{\Gamma}^{-}$ of $\widehat{\Gamma}$. Under Assumptions
\ref{ass: Ass2} and \ref{ass:Ass4}, $\widehat{\Gamma}$ will be
nonsingular with probability approaching one (cf. Lemma \ref{lem:Nonsingular},
Section \ref{sec:Auxiliary-Results} in the Supplemental Material),
and hence $\widehat{\Gamma}^{-}$ will be the standard inverse.\par\end{sloppy}
\begin{thm}
\label{thm:Thm4}If Assumptions \ref{ass: Ass2}-\ref{ass:Ass4} hold,
then (i) there exists $\hat{b}$ in $\Theta$ with probability approaching
one; (ii) $\hat{b}\rightarrow_{p}b^{*}$; and (iii) $n^{\frac{1}{2}}(\hat{b}-b^{*})\rightarrow_{d}N(0,\Gamma^{-1}\Psi\Gamma^{-1})$.
Moreover, $\widehat{\Gamma}^{-}\widehat{\Psi}\widehat{\Gamma}^{-}\rightarrow^{p}\Gamma^{-1}\Psi\Gamma^{-1}$.
\end{thm}
Theorem \ref{thm:Thm4}(i) shows that correct specification is not
required for existence of a globally monotone estimate $\hat{b}'T(X,Y)$
for large enough samples. Theorem \ref{thm:Thm4}(ii) holds without
compactness of $\Theta$, by concavity of the objective (e.g., \citet{Newey:McFadden:1994},
Section 2.6). Under correct specification, $\Gamma=-\Psi$ by the
information matrix equality and the estimator is efficient, with asymptotic
variance-covariance matrix $-\Gamma^{-1}$ (e.g., \citet{Newey:McFadden:1994}).
A valid method for model selection is adaptive Lasso ML (\citet{Zou:2006},
\citet{LGF:2012}, \citet{HN:2020}). The selected model is a sparse
KLIC optimal approximation for $g(Y,X)$.
\begin{rem}
Lasso penalized ML allows for the dimension of $W(X)$ to increase with sample size.
From (\ref{eq:Varying}) it is apparent that the objective function has the multiple-index structure $E[m(W(X)'b_{1},\ldots,W(X)'b_{J},Y)]$,
with strictly convex negative log-likelihood $(z_{1},\ldots,z_{J})\mapsto m(z_{1},\ldots,z_{J},Y)$.
High-dimensional estimation results and penalty selection methods
in \citet{CS:2021} apply to this case.
\end{rem}
\subsection{Estimation of DRFs}
An estimator for $g^{*}(y,x)$ is $\widehat{g}^{*}(y,x)\equiv T(x,y)'\widehat{b}$,
and estimators for DRFs are
\[
\widehat{f}^{*}(y,x)\equiv\phi(\widehat{g}^{*}(y,x))\{\partial_{y}\widehat{g}^{*}(y,x)\},\quad\widehat{F}^{*}(y,x)\equiv\Phi(\widehat{g}^{*}(y,x)),\quad(y,x)\in\mathcal{YX},
\]
and
\[
\widehat{Q}^{*}(x,u)\equiv\{y\in\mathbb{R}:\Phi(\widehat{g}^{*}(y,x))=u,\;\partial_{y}\widehat{g}^{*}(y,x)>0\},\quad x\in\mathcal{X},\quad u\in\mathcal{U}_{x}(g^{*}).
\]
These estimators are known functionals of $\widehat{b}$, and hence
their asymptotic distribution follows by application of the Delta
method.
\begin{thm}
\label{thm:Thm5}Suppose that $\Xi\equiv\Gamma^{-1}\Psi\Gamma^{-1}$
is positive definite. Under Assumptions \ref{ass: Ass2}-\ref{ass:Ass4}
we have: (i) for $(y,x)\in\mathcal{Y}\mathcal{X}$,
\[
n^{\frac{1}{2}}(\widehat{f}^{*}(y,x)-f^{*}(y,x))\rightarrow_{d}N(0,\phi(g^{*}(y,x))^{2}\Delta(x,y)'\Xi\Delta(x,y)),
\]
where $\Delta(x,y)\equiv-g^{*}(y,x)\{\partial_{y}g^{*}(y,x)\}T(x,y)+t(x,y)$,
and
\[
n^{\frac{1}{2}}(\widehat{F}^{*}(y,x)-F^{*}(y,x))\rightarrow_{d}N(0,\phi(g^{*}(y,x))^{2}T(x,y)'\Xi T(x,y));
\]
(ii) for $x\in\mathcal{X}$, $u\in\mathcal{U}_{x}(g^{*})$,
\[
n^{\frac{1}{2}}(\widehat{Q}^{*}(x,u)-Q^{*}(x,u))\rightarrow_{d}N(0,\{\partial_{y}g^{*}(y_{0},x)\}^{-2}T(x,y_{0})'\Xi T(x,y_{0})),
\]
where $y_{0}=Q^{*}(u,x)$.
\end{thm}
Theorem \ref{thm:Thm4} provides an estimator for the asymptotic variance-covariance
matrix $\Xi$.
\begin{rem}
To implement (\ref{eq:ML}) we expand the
parameter space $\Theta$ to the larger space $\Theta_{n}=\{b\in\mathbb{R}^{JK}:b't(x_{i},y_{i})>0,\,i\in\{1,\ldots,n\}\}$,
the effective domain of $Q_{n}(b)$. This implies that there is
$b\in\Theta_{n}$ such that $b't(X,Y)\leq0$ with positive probability.
One can verify that $\widehat{b}\in\Theta$ holds after estimation
by checking the quasi-global monotonicity (QGM) property $\widehat{b}'t(x,y)>0$
on a fine grid of values that covers $\mathcal{Y}\times\mathcal{X}$.\label{Rk:QGM}
\end{rem}
\subsection{Comparison with alternative methods}
Our approach is related to distribution regression estimation of conditional
CDFs, for the (probit) model
\begin{equation}
F_{Y|X}(y\mid X)=\Phi(W(X)'\beta(y)),\quad y\in\mathbb{R},\label{eq:DRmodel}
\end{equation}
where $\beta(y)$ a vector of unknown functions. One sense in which
our approach is flexible is that model (\ref{eq:e}) can approximate
$F_{Y|X}(y|X)$ in (\ref{eq:DRmodel}) arbitrarily well for $S(y)$
rich enough.\footnote{Formally, write $S^{J}(Y)=S(Y)$ and suppose that Assumptions \ref{ass: Ass2}
holds, $E[||\beta(Y)||^{2}]<\infty$ and that for any $K$ vector
of functions $b(Y)$ with $E[||b(Y)||^{2}]<\infty$ there are $J\times1$
vectors $a_{k}^{J}$, $k\in\{1,\ldots,K\}$, such that $E[\sum_{k=1}^{K}\{b_{k}(Y)-S^{J}(Y)'a_{k}^{J}\}^{2}]\rightarrow0$
as $J\rightarrow\infty$. Then, $E[\{F_{Y|X}(Y|X)-\Phi(b'[W(X)\otimes S^{J}(Y)])\}^{2}]\rightarrow0$
as $J\rightarrow\infty$, by an argument similar to Theorem 11 in \citet{NeweyStouli:2025}.} If we define $\beta^{*}(y)=(\beta_{1}^{*}(y),\ldots,\beta_{K}^{*}(y))'$
with $\beta_{k}^{*}(y)\equiv S(y)'b_{k}^{*}$ and $b_{k}^{*}=(b_{k1}^{*},\ldots,b_{kJ}^{*})'$,
$k\in\{1,\ldots,K\}$, then $F^{*}(y,X)=\Phi(W(X)'\beta^{*}(y))$
provides increasingly accurate KLIC optimal approximations to (\ref{eq:DRmodel})
as $J$ increases. When $\beta^{*}(y)=\beta(y)$, our formulation
characterizes $\beta(y)$ for each $y$ in $\mathbb{R}$ simultaneously
whereas distribution regression characterizes $\beta(y)$ pointwise.
While both corresponding estimators are $\sqrt{n}$ consistent, the
ML estimator (\ref{eq:ML}) is efficient. When (\ref{eq:DRmodel})
is misspecified and multiple values of $Y$ are of interest, the distribution
regression criterion does not provide approximation guarantees for
$F_{Y|X}(y|X)$, whereas $F^{*}(y,X)$ is KLIC optimal.
Both conditional CDF models in (\ref{eq:DRF1}) and (\ref{eq:DRmodel})
are particular cases of conditional transformation models of the form
$F(g(Y,X))=F(\sum_{l=1}^{L}g_{l}(Y,X))$ (\citet{HKB:2014}), for
some specified CDF $F$. For estimation, taking this additive structure
as a starting point leads to considering the restricted form
\begin{equation}
F(g(Y,X))=F(\sum_{l=1}^{L}g_{l}(Y,X_{l}))=F(\sum_{l=1}^{L}b_{l}'[W_{l}(X_{l})\otimes S(Y)]),X=(X_{1},\ldots,X_{L})'\label{eq:CTM}
\end{equation}
(\citet[Section 5]{HKB:2014}), with monotonicity constraints on each
partial transformation $g_{l}$ (\citet[Section 4.5]{HMB:2018}, \citet[p. 1362]{CKK:2024}).
Here we find that the additive structure imposed by conditional transformation
models is not required for the formulation of fully flexible models
for conditional CDFs. Model (\ref{eq:DRF1}) gives a flexible generalization
of (\ref{eq:CTM}) and does not impose unnecessary monotonicity
restrictions.
Our approach is also related to quantile regression estimation of
CQFs. For the quantile regression model $Q_{Y|X}(u|X)=W(X)'\beta(u)$,
$u\in(0,1)$, the coefficients $\beta(u)$ are estimated at each $u$
by a sequence of linear programming problems. In general, this model
and (\ref{eq:e}) are not nested, but they coincide for $W(X)$ and
$S(Y)$ rich enough. Compared to our approach, the quantile regression
loss function allows for outcomes with non finite second moment, and
enjoys robustness properties when estimating a specific quantile.
When multiple quantiles are of interest, quantile regression does
not provide approximation guarantees for $Q_{Y|X}(u|X)$, whereas
$Q^{*}(X,u)$ is KLIC optimal.
Quantile and distribution regression can result in both finite sample
estimates and population approximations under misspecification that
do not satisfy the monotonicity properties of CQFs and CDFs, respectively.\footnote{This problem has motivated the development of a variety of methods
to avoid or repair intersecting quantile surfaces (e.g., \citet{He:1997},
\citet{DV:2008}, \citet{Chernozhukov Fern Gali 2010}, \citet{YT:2017},
\citet{SS:2018a}).} When this occurs, rearrangement methods provide improvements
(\citet{Chernozhukov Fern Gali 2010}),
but without optimality guarantees. Here we propose a one-step resolution
to the quantile and probability curves crossing problem through the
information-theoretic selection of globally monotone models within
a specified class. Our criterion performs global model comparisons
within a set of proper conditional PDFs, and hence allows for the
selection of an optimal KLIC approximation. The implied CQF is then
free of crossing and the conditional CDF monotone, both are optimal
in a transparent sense, and they are estimated at parametric rate
together with the conditional PDF.
Compared to nonparametric kernel-based methods, flexible estimation
at parametric rate alleviates the curse of dimensionality when $X$
is a vector. In our simulations we find that it also yields substantial
finite sample improvements for conditional PDF, CDF and nonseparable
model estimation when $X$ is a scalar. The QGM property in Remark \ref{Rk:QGM}
avoids the need for correcting estimates (\citet{GHI:2003}) and multiplicity
of inverses (cf. \citet{Matzkin:2003}, p. 1358). The linear form
of GT estimates makes the QGM property easy to check and has the advantage
that imposing shape constraints from economic theory in estimation
is straightforward. Compared to \citet{Matzkin:2003}, the direct
specification of the distribution of $e$ through the log density
in (\ref{eq:log PDF}) also simplifies the choice and imposition of
a normalization.
Because $Q^{*}(X,u)$ is a proper CQF it can be used for simulation
of data that mimics the true data generating process, using the representation
$\widetilde{Y}=Q^{*}(X,U)$, $U|X\sim U(0,1)$. Thus, our approach
also complements Generative Adversarial Networks (GAN) (\citet{Goodfellow:2014})
used in Econometrics for simulation of complex datasets (\citet{AIMM:2024}).
GANs are flexible predictive methods able to mitigate the curse of
dimensionality. Compared to our approach, GANs do not directly produce
conditional distribution estimates and require solving computationally
challenging nonconvex nonconcave min-max games. Recent proposals substitute
the Wasserstein distance in the objective function and use a gradient-based
penalty for more stable implementation (\citet{AB:2017}). In contrast
our approach uses the KLIC to perform ML estimation of conditional
distributions in closed-form, while global concavity alleviates computational
difficulties.
\section{Duality Theory\label{sec:Section5}}
Considerable computational advantages accrue from our ML criterion
where the GT enters in closed-form. To (\ref{eq:ML}) corresponds
a dual formulation that can be cast into the modern convex programming
framework (\citet{BV:2004}). We derive the dual problem and establish
the properties of the dual solutions.
\begin{thm}
\label{thm:Thm6} If Assumptions \ref{ass: Ass2}-\ref{ass:Ass4}
are satisfied then the following hold.
(i) The dual of (\ref{eq:ML}) is
\begin{align}
\min_{(u,v)\in\mathbb{R}^{n}\times(-\infty,0)^{n}} & -n\left(\frac{1}{2}\log(2\pi)+1\right)+\sum_{i=1}^{n}\left\{ \frac{u_{i}^{2}}{2}-\log\left(-v_{i}\right)\right\} \label{eq:Dual objective}\\
\textrm{subject to}\quad & -\sum_{i=1}^{n}\left\{ T(x_{i},y_{i})u_{i}+t(x_{i},y_{i})v_{i}\right\} =0\label{eq:Dual scores}
\end{align}
the dual GT regression problem, with solution $\widehat{\alpha}=(\widehat{u}',\widehat{v}')'$.
(ii) The dual GT regression program (\ref{eq:Dual objective})-(\ref{eq:Dual scores})
admits the method-of-moments representation
\begin{equation}
\sum_{i=1}^{n}\left\{ -T(x_{i},y_{i})\{b'T(x_{i},y_{i})\}+\frac{t(x_{i},y_{i})}{b't(x_{i},y_{i})}\right\} =0,\label{eq:MMrep}
\end{equation}
the first-order conditions of (\ref{eq:ML}).
(iii) With probability approaching one we have: (a) existence and
uniqueness, i.e., there exists a unique pair $(\widehat{b}',\widehat{\alpha}')'$
that solves (\ref{eq:ML}) and (\ref{eq:Dual objective})-(\ref{eq:Dual scores}),
and
\begin{equation}
\widehat{u}_{i}=\widehat{b}'T(x_{i},y_{i}),\quad\hat{v}_{i}=-\frac{1}{\widehat{b}'t(x_{i},y_{i})},\quad i\in\{1,\ldots,n\};\label{eq:Dual_FOCs}
\end{equation}
(b) strong duality, i.e., the value of (\ref{eq:ML}) equals the value
of (\ref{eq:Dual objective})-(\ref{eq:Dual scores}).
\end{thm}
\begin{sloppy}
The dual formulation in Theorem \ref{thm:Thm6} demonstrates important
computational properties of GT regression. The Hessian matrix of the
dual problem (\ref{eq:Dual objective})-(\ref{eq:Dual scores}) is
\[
\left[\begin{array}{cc}
I_{n} & 0_{n\times n}\\
0_{n\times n} & \textrm{diag}(1/v_{i}^{2})
\end{array}\right],
\]
a positive definite diagonal matrix for all $v\in(-\infty,0)^{n}$,
with $I_{n}$ denoting the $n\times n$ identity matrix and $\textrm{diag}(1/v_{i}^{2})$
the $n\times n$ diagonal matrix with elements $(1/v_{1}^{2},\ldots,1/v_{n}^{2})$.
Thus the dual problem is a strictly convex mathematical program with
sparse Hessian matrix and $JK$ linear constraints. We implement this
computationally convenient formulation using the state-of-the-art
convex programming solvers ECOS (\citet{ECOS:2013}) and SCS (\citet{ODCPB2016}).
\par\end{sloppy}
It is useful to compare problem (\ref{eq:Dual objective})-(\ref{eq:Dual scores})
with the dual formulation of quantile regression. The dual problem
for the linear $\tau$ quantile regression of $y_{i}$ on $W(x_{i})$
is
\[
\max_{u}\left\{ y'u:\sum_{i=1}^{n}W(x_{i})u_{i}=(1-\tau)\sum_{i=1}^{n}W(x_{i}),\quad u\in[0,1]^{n}\right\} ,\quad\tau\in(0,1),
\]
with solutions that take value $0$ or $1$, except for $K$ sample
points that define the fitted quantile regression surface and are
assigned $u$ values that are neither $0$ nor $1$ (cf. Chapter 3.5.4
in \citet{Koenker:2005}). When multiple quantiles are of interest,
violation of monotonicity may arise because each quantile surface
must interpolate $K$ sample points, while each quantile regression
is implemented separately. In contrast, (\ref{eq:Dual objective})-(\ref{eq:Dual scores})
assigns $u$ values in $\mathbb{R}^{n}$, avoiding box constraints
on $u$, while imposing monotonicity since $-1/v_{i}>0$ at a solution.
Implied quantile surfaces are then unrestricted beyond the form of
$u_{i}$ in (\ref{eq:Dual_FOCs}), and strong duality ensures KLIC
optimality.
In addition to KLIC optimality and the logarithmic barrier in the
objective, linearity of the constraints is also an important advantage
of (\ref{eq:Dual objective})-(\ref{eq:Dual scores}) compared to
the alternative generalized dual regression characterization of conditional
CDFs and CQFs (\citet{SS:2018a}) for which the mathematical program
is of the form
\begin{equation}
\max\left\{ y'e:\sum_{i=1}^{n}\mathcal{T}(x_{i},e_{i})=0,\quad e\in\mathbb{R}^{n}\right\} ,\label{eq:DualReg}
\end{equation}
where $\mathcal{T}(x_{i},e_{i})$ is a vector of known functions of
$x_{i}$ and $e_{i}$ including $e_{i}$ and $(e_{i}^{2}-1)/2$, so
that $e$ enters nonlinearly into the constraints. (\ref{eq:DualReg})
has first-order conditions
\begin{equation}
y_{i}=\widehat{b}'\{\partial_{e_{i}}\mathcal{T}(x_{i},e_{i})\},\quad i\in\{1,\ldots,n\},\label{eq:GDR_FOCs}
\end{equation}
where $\widehat{b}$ is the Lagrange multiplier vector for the constraints
in (\ref{eq:DualReg}), but where the solution is now determined by
a system of $n$ nonlinear equations instead of having a closed-form
expression as in (\ref{eq:Dual_FOCs}). This further illustrates the
benefits accruing from closed-form modeling of $e=g(Y,X)$, compared
to direct modeling of $y_{i}$ in (\ref{eq:GDR_FOCs}).
\begin{rem}
Compared to the methods above, (\ref{eq:Dual objective})-(\ref{eq:Dual scores})
provides a one-step estimator for both conditional PDFs and CDFs at
each sample point, formed as $\varPhi(\widehat{u}_{i})$ and $\phi(\widehat{u}_{i})(-1/\hat{v}_{i})$,
respectively. Moreover, uniform convergence of the empirical distribution
of $\{\widehat{u}_{i}\}_{i=1}^{n}$ can be used to establish the validity
of bootstrap methods for uniform inference on DRFs (\citet{Chernozhukov Fern Melly 2013}). Uniform convergence
follows from the method-of-moments representation (\ref{eq:MMrep})
and steps similar to the proof of Theorem 6 in \citet{SS:2018a}.
\end{rem}
\section{Distributional Gender Wage Gap Analysis\label{sec:Section6}}
We illustrate our framework with an application to the distributional
gender wage gap in the United States. We use data on wages, hours
and weeks worked, educational attainment, gender, race, age, industry
and occupation from the 2019 American Community Survey (\citet{Rugglesetal:2025}).\footnote{A large representative
survey that covers 1\% of the U.S. population,
with mandatory participation, and available on the IPUMS-USA website
(https://usa.ipums.org/usa/).} For a sample of white employees working full time (more than $34$
hours a week, at least $50$ weeks a year) in metropolitan areas, we stratified the data
according to industry-occupation pairs, and selected the 41 pairs
for which the support of education is the same across genders, to allow for wage comparisons over
the whole support.
We follow \citet{BCS:2024} and include individuals
aged $25$ to $65$, and we discard individuals below the $7.25$ federal minimum hourly
wage.
We have 199,785 observations in total, with 245 to 32,828 observations per industry-occupation pair.
For each pair, we estimate DRFs implied
by the model
\begin{equation}
e^{*}=g^{*}(Y,X,D)=[\{W(X)\otimes(1,D)'\}\otimes S(Y)]'b^{*},\quad X=(X_{1},X_{2})',\label{eq:e* wages}
\end{equation}
where $Y$ denotes weekly wages, $X_{1}$ years of education, $X_{2}$
experience, and $D$ is a gender dummy variable. This provides
a flexible framework for estimation of distributional and quantile
treatment effects. Differences in wages across genders occur if some
component of $D[W(X)\otimes S(Y)]$ in (\ref{eq:e* wages}) has nonzero
coefficient. The sign, magnitude and location of gender differences
in the wage distribution can be analyzed by estimating the quantile
gender wage gap over the support of $X$ and $u\in(0,1)$,
\begin{equation}
\mathtt{GWG}(X,u)=100\times\left[\log Q^{*}(X,D=\mathtt{Male},u)-\log Q^{*}(X,D=\mathtt{Female},u)\right].\label{eq:qGWG}
\end{equation}
Following \citet{BCS:2024} we treat education and experience as continuous,
and for both $W(X)$ and $S(Y)$ we consider a range of spline transformations.\footnote{For model selection
we implement an adaptive Lasso ML version of our
estimator, and we select the specification that minimizes the Bayes
Information Criterion (BIC). Details of our implementation for the empirical application are given
in Section \ref{sec:Implementation} of the Supplemental material.
All computational procedures can be implemented in the software R
(\citet{R:2024}) using open source packages for convex optimization
such as CVX, and its R implementation CVXR (\citet{CVXR}).} Spline functions satisfy the conditions of our modeling framework
and have been demonstrated to be remarkably effective when applied
to the related problems of log density estimation (\citet{KS:1991})
or monotone regression function estimation (\citet{Ramsay:1988}).
Thus, (\ref{eq:e* wages})-(\ref{eq:qGWG}) is a flexible specification
for the gender wage gap, allowing for both observed and unobserved
heterogeneity.
\begin{figure}[t]
\begin{centering}
\subfloat[Conditional PDF of earnings by $\textrm{Years of Education}\in\{12,16,19,20\}$.\label{fig:PDF}]{\begin{centering}
\includegraphics[width=7.25cm,height=5cm]{main_Fig6_1_A_left}\hfill{}\includegraphics[width=7.25cm,height=5cm]{main_Fig6_1_A_right}
\par\end{centering}
}
\par\end{centering}
\begin{centering}
\subfloat[Conditional CDF of earnings by $\textrm{Years of Education}\in\{12,16,19,20\}$.\label{fig:CDF}]{\begin{centering}
\includegraphics[width=7.25cm,height=5cm]{main_Fig6_1_B_left}\hfill{}\includegraphics[width=7.25cm,height=5cm]{main_Fig6_1_B_right}
\par\end{centering}
}
\par\end{centering}
\caption{Conditional PDF and CDF for $\mathtt{Female}$ (left) and $\mathtt{Male}$
(right), with confidence bands.\label{fig:PDFsCDFs}}
\end{figure}
Distributional gender wage gap analysis is challenging because the
shape of the conditional wage distribution varies across $X$ values,
and because of key features of the dataset such as high skewness of
wages, the presence of large outliers, multiple covariates, and the
effectively discrete measurements of education and experience. Applying
high-dimensional mean regression to similar data, \citet{BCS:2024}
find that accounting for heterogeneity in observables (e.g., education,
industry and occupation) is important in understanding the gender gap.
We complement their analysis by studying distributional
wage differences across genders, and flexibly modeling observed heterogeneity
across education and experience levels. For each industry-occupation
pair, we obtain DRFs over
their entire support, test for their equality across genders,
estimate $\mathtt{GWG}(X,u)$, and give confidence bands for
all objects.
\begin{figure}[t]
\subfloat[CQF with confidence bands, for $u\in\{0.1,0.25,0.5,0.75,0.9\}$.\label{fig:CQF_CI}]{\begin{centering}
\includegraphics[width=7.25cm,height=5cm]{main_Fig6_2_A_left}\hfill{}\includegraphics[width=7.25cm,height=5cm]{main_Fig6_2_A_right}
\par\end{centering}
}
\subfloat[CQF with scatterplots by gender, for $u\in\{0.05,0.10,\ldots,0.95\}$.\label{fig:CQF_Full}]{\includegraphics[width=7.25cm,height=5cm]{main_Fig6_2_B_left}\hfill{}\includegraphics[width=7.25cm,height=5cm]{main_Fig6_2_B_right}
}\caption{CQF for $\mathtt{Female}$ (left) and $\mathtt{Male}$ (right).\label{fig:CQFs}}
\end{figure}
Figures \ref{fig:PDFsCDFs}-\ref{fig:CQFs} show DRF estimates for the legal occupation in the services industry
and for $X_{2}=\widehat{Q}_{X_{2}}(0.5)$, the sample experience median.
This example demonstrates that
parsimonious Spline-Spline models
can capture a wide variety of complex data features. The high
single peaked PDFs at $X_{1}\in\{12,16\}$ differ markedly from the
long-tailed and bimodal PDFs at $X_{1}\in\{19,20\}$.\footnote{$X_{1}\in\{12,16,19,20\}$ indicates completion of high school, bachelor,
professional and doctoral degrees, respectively.} This reveals two different regimes in the conditional wage distribution:
high concentration/relatively low wages at $X_{1}\in\{12,16\}$, and
high dispersion/relatively high wages at $X_{1}\in\{19,20\}$. The
two modes in the second regime are especially apparent for males in
Figure \ref{fig:PDFsCDFs}(A), and are reflected by the two inflection
points in the corresponding CDFs in Figure \ref{fig:PDFsCDFs}(B),
and by the large gap between the low and high CQFs
in Figure
\ref{fig:CQFs}(A). CQF estimates capture both linearity and nonmonotonicity
over the $X_{1}$ support, reflecting substantial heteroskedasticity
and changes in mode locations for the wage distribution.
Lower dispersion of CQFs up to the median reflects distributional asymmetry.
Visual inspection strongly suggests that the DRFs differ across genders. Conditional PDFs for females are
higher for low wages, and those for males are higher
for high wages.
Conditional
CDFs for males stochastically dominate those for females.
Upper quantile CQFs are higher for males. A Wald test of DRFs equality
across genders reinforces this
diagnostic.\footnote{We perform a significance test for the $9$ nonzero coefficients of
$D[W(X)\otimes S(Y)]$ in (\ref{eq:e* wages}), with a test statistic
of $316>16.9$, the critical value at the $5\%$ level.}
Figure \ref{fig:QTE} gives a comprehensive picture of the quantile
gender wage gap. We find substantial heterogeneity, with a statistically
significant gap over the whole support of $X_{1}$ for higher quantiles
(median and above). We also find nonlinearity, with close to a constant
gap up to year 16, and then marked divergence across quantiles.
Overall, we find that parsimonious representations are able to capture
complex distributional data features. In the Supplemental Material we further
illustrate the overall finding that $\mathtt{GWG}(X,u)$ varies across
quantiles and exhibit substantial nonlinearity in $X_{1}$, with heterogeneous
patterns across industry-occupation pairs.
The main features of the
selected Spline-Spline model are well-preserved by models with similar
penalization and BIC.
In our simulations we also find that BIC
performs well in a variety of designs. Thus, although
model selection
in our context is an important topic for future research,
BIC appears to be reliable for practical purposes.
\begin{figure}[t]
\includegraphics[width=9cm,height=5cm]{main_Fig7}
\caption{Quantile gender wage gap.\label{fig:QTE}}
\end{figure}
\section{Extensions\label{sec:Section7}}
\subsection{Logistic transform regression}
This paper has exposited a theory of Gaussian transformation regression
based on the standard normal density. Other probability densities
could be employed instead, with the standard Logistic and Laplace
densities being appealing choices. We provide a brief analysis of
the Logistic case here.
We write as before (cf. (\ref{eq:e})) $e=b_{0}'T(X,Y)$, $\partial_{y}\{b_{0}'T(X,Y)\}=b_{0}'t(X,Y)>0$,
but replace $e\mid X\sim N(0,1)$ with $e\mid X\sim\Lambda$ where
$\Lambda$ is the standard logistic CDF with log-density $e-2\log(1+\exp(e))$.
Objective function (\ref{eq:popML}) is replaced by
\[
Q(b)=E\left[b'T(X,Y)-2\log\left(1+\exp(b'T(X,Y))\right)+\log\left(b't(X,Y)\right)\right],\quad b\in\Theta.
\]
The corresponding first-order conditions that replace (\ref{eq:FOCs})
are
\[
E\left[-T(X,Y)\frac{\exp(b'T(X,Y)-1)}{\exp(b'T(X,Y)+1)}+\frac{t(X,Y)}{b't(X,Y)}\right]=0.
\]
The population objective function remains concave as in the Gaussian
case, since the expected Hessian is now given by:
\begin{equation}
-E\left[2\lambda(b'T(X,Y))\left\{ T(X,Y)T(X,Y)'+\frac{t(X,Y)t(X,Y)'}{\{b't(X,Y)\}^{2}}\right\} \right],\label{eq:hessian-1}
\end{equation}
where $\lambda(\cdot)$ is the standard Logistic density. In the Gaussian
version of (\ref{eq:hessian-1}) the term corresponding to $2\lambda(b'T(X,Y))$
is ``1'' because the log Gaussian density is a simple quadratic
with coefficient $(-1/2)$ so that its second derivative is non-stochastic.
Our limited experience with estimation based on the Logistic density
suggests that it has no particular advantage or disadvantage with
respect to the Gaussian.
\subsection{Mixed discrete-continuous outcomes}
Our formulation extends to the case where $Y$ has both a continuous component with support $\mathcal{Y}_{C}$ and a discrete component with support $\mathcal{Y}_{D}\equiv \{y_{1},\ldots,y_{M}\}$.
For each $y\in\mathcal{Y}_{D}$, the PDF of $Y$ conditional on $X$
can be expressed as $\Pr[Y=y|X]=F_{Y|X}(y|X)-\lim_{z\rightarrow y^{-}}F_{Y|X}(z|X)$ (\citet{Mittelhammer:2013}, p. 67), where
$\lim_{z\rightarrow y^{-}}$ is the limit as $z$ approaches $y$ from below.
For some specified continuous and strictly increasing CDF $F$ with derivative $f$, define $g(Y,X)\equiv F^{-1}(F_{Y|X}(Y|X))$. Then $f_{Y|X}(Y|X)$ can now be expressed as
\[
f_{Y|X}(Y|X) =\{F(g(Y,X))-\lim_{z\rightarrow Y^{-}}F(g(z,X))\}^{1(Y\in\mathcal{Y}_{D})}\,\{f(g(Y,X))\partial_{y}g(Y,X)\}^{1(Y\in\mathcal{Y}_{C})}.
\]
Let $D(Y)$ be a vector of dummy variables $D_{m}$, $m\in \mathcal{M}\equiv \{2,\ldots, M\}$, taking value one if outcome $y_{m}$ occurs and zero otherwise. For $Y\in \mathcal{Y}_{D}$ and $M \geq 2$, we take
$g(Y,X)=d_{0}'T_{D}(X,Y)$, where $T_{D}(X,Y)=W(X)\otimes P(Y)$ with $P(Y)=(1,D(Y)')'$. When $M=1$, we set $P(Y)=1$. For $Y\in \mathcal{Y}_{C}$, we take
$g(Y,X)=c_{0}'T_{C}(X,Y)$, where $T_{C}(X,Y)=W(X)\otimes Q(Y)$ with $Q(Y)=(1,C(Y)')'$, for $C(Y)$ a vector of known functions of $Y$. For $F(g(Y,X))$ to be a valid CDF model, we assume $(d_{0}',c_{0}')'$, and $C(Y)$ are such that $y\mapsto g(y,X)$ is right-continuous,
so that $y \mapsto F(g(y,X))$ also is.
To illustrate, let $\mathcal{Y}_{C} = (y_{M},\infty)$,
$C(Y)$ be such that $C(y_{M})=0$, $d_0=(d_{01}',\ldots, d_{0M}')'$, and $\delta_{m}(X)\equiv d_{0m}'W(X)$, $m\in \{1,\mathcal{M}\}$. For $Y \in \mathcal{Y}_{D}$, $P(Y)=(1,D(Y)')'$, and hence:
\[
g(Y,X)=d_{0}'T_D(X,Y)=d_{01}'W(X) + \sum_{m\in \mathcal{M}}\{d_{0m}'W(X)\}D_{m}=\delta_{1}(X) +
\sum_{m\in \mathcal{M}}\delta_{m}(X)D_{m}.
\]
Similarly, for $J\equiv \dim(Q)$, write $c_{0}=(c_{01}',\ldots, c_{0J}')'$ and let $\xi_{j}(X)\equiv c_{0j}'W(X)$, $j\in \{1,\mathcal{J}\}$, $\mathcal{J}\equiv \{2,\ldots,J\}$.
For $Y \in \mathcal{Y}_{C}$, $Q(Y)=(1,C(Y)')'$, and hence: $g(Y,X)=\xi_{1}(X) + \sum_{j\in \mathcal{J}}\xi_{j}(X)C_{j-1}(Y)$.
By $\sum_{j\in \mathcal{J}}\xi_{j}(X)C_{j-1}(y_{M})=0$ and $\sum_{m\in \mathcal{M}}\delta_{m}(X)D_{m}=\delta_{M}(X)$ when $Y=y_{M}$, setting $\delta_1(X) + \delta_{M}(X) = \xi_{1}(X)$ implies right-continuity of $y\mapsto g(y,X)$ at $y_{M}$.
If $E[W(X)W(X)']$ is nonsingular, this holds for $d_{01}+d_{0M}=c_{01}$.
\begin{rem}
For $M=1$ and $Q(Y)=(1,Y)'$, an important particular case is the Tobit model (\citet{Tobin:1958}). With censoring at 0, the standard model has $Y=W(X)'\beta+\sigma e$ if $W(X)'\beta+\sigma e>0$, $e|X\sim\Phi$, and $Y=0$ otherwise. For $Y>0$ this is the same as $e=-W(X)'(\beta/\sigma) + (1/\sigma)Y$, which is of the form $c_{0}'T_{C}(X,Y)$ above for $F=\Phi$ and $c=(c_{01}', c_{02})'\equiv (-\beta'/\sigma, 1/\sigma)'$. This is \citet{Olsen:1978}'s concave proposal commonly used in practice, and our formulation provides a nonlinear and nonseparable generalization.
\begin{rem}
With $Y$ discrete, we have $\lim_{z\rightarrow y_{1}^{-}}F_{Y|X}(z|X)=0$, $\lim_{z\rightarrow y_{m}^{-}}F_{Y|X}(z|X)=F_{Y|X}(y_{m-1}|X)$, $m\in \mathcal{M}$, and $F_{Y|X}(y_M|X)=1$.
Moreover, specifications with $S(Y)=(1,D(Y)')'$ are not restrictive in the $Y$ dimension, and coincide with the semiparametric distribution regression model $F_{Y|X}(Y|X)=F(W(X)'\beta(Y))$.
In contrast with distribution regression, monotonicity is obtained directly from
the implied objective function, with the terms $\log(F(b'T(X,y_{m}))-F(b'T(X,y_{m-1})))$
ruling out nonincreasing CDFs. When also $X$ is discrete, with support $\{x_{1},\ldots,x_{L}\}$, a nonparametric formulation is achieved with $W(X)=(1,Z(X)')'$, for $Z(X)$ a vector of dummy variables $Z_{l}$, $l=2,\ldots,L$, taking value one when $X=x_{l}$ occurs and zero otherwise.
\end{rem}
\end{rem}
\subsection{Multiple outcomes}
With multiple outcomes $(Y_{1},\ldots,Y_{M})'\equiv Y$, $M\geq2$,
writing $\mathbb{Y}_{m}\equiv(Y_{1},\ldots,Y_{m})'$, a compact generalization
of (\ref{eq:e}) is the recursive formulation
\begin{align*}
e_{m} & =T_{m}(X,\mathbb{Y}_{m})'b_{0,m},\quad e_{m}\mid X,\mathbb{Y}_{m-1}\sim N(0,1),\quad m\in\{2,\ldots,M\},\\
e_{1} & =T_{1}(X,Y_{1})'b_{0,1},\quad e_{1}\mid X\sim N(0,1),
\end{align*}
where $T_{m}(X,\mathbb{Y}_{m})\equiv T_{m-1}(X,\mathbb{Y}_{m-1})\otimes S_{m}(Y_{m})$
and $T_{1}(X,Y_{1})\equiv W(X)\otimes S_{1}(Y_{1})$, with
\begin{align*}
\partial_{y_{m}}\{T_{m}(X,\mathbb{Y}_{m})'b_{0,m}\} & =t_{m}(X,\mathbb{Y}_{m})'b_{0,m}>0,\quad m\in\{1,\ldots,M\},
\end{align*}
where $t_{m}(X,\mathbb{Y}_{m})\equiv t_{m-1}(X,\mathbb{Y}_{m-1})\otimes s_{m}(Y_{m})$,
$m\in\{2,\ldots,M\}$, and $t_{1}(X,Y_{1})\equiv W(X)\otimes s_{1}(Y_{1})$.
By construction, the $e_{m}$'s are jointly Gaussian and mutually
independent, with variance-covariance the identity matrix. This is a Gaussian version of \citet{Rosenblatt:1952}'s multivariate
probability transformation. The conditional
CDF is
\[
F_{Y\mid X}(y_{1},\ldots,y_{M}\mid X)=\int_{-\infty}^{y_{1}}\ldots\int_{-\infty}^{y_{M}}f_{Y\mid X}(t_{1},\ldots,t_{M}\mid X)dt_{1}\ldots dt_{M},
\]
where the conditional PDF takes the form
\[
f_{Y\mid X}(\overline{y}_{M}\mid X)=\prod_{m=1}^{M}\phi(T_{m}(X,\overline{y}_{m})'b_{0,m})\{t_{m}(X,\overline{y}_{m})'b_{0,m}\},\quad\overline{y}_{m}\equiv(y_{1},\ldots,y_{m})\in\mathbb{R}^{m}.
\]
\section{Conclusion\label{sec:Section8}}
The formulation of flexible models for the GT $e=g(Y,X)$ leads to
a unifying information-theoretic framework for the global estimation
of DRFs. The implied convex programming
formulation is easy to implement and the linear form of the proposed
GT regression models also constitutes a good starting point for nonparametric
estimation. In this paper we have considered a few extensions to our
original formulation such as misspecification, quantile treatment
effects, Logistic transform regression,
mixed discrete-continuous distributions and multiple outcomes.
Two important further extensions for future work are
sample selection and endogenous regressors.