Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,253 characters · 19 sections · 91 citation commands
Gaussian Transforms Modeling and the Estimation of Distributional Regression Functions
Keywords:{ Conditional density estimation, conditional distributions, conditional quantiles, maximum likelihood, misspecification, monotonicity, convexity, gender wage gap.}
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., 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., Imbens Newey 2009, CFNSV). Both conditional PDFs and CDFs are key components in counterfactual analysis (e.g., DiNardoetal, 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., HM:2003, GHI:2003, Chernozhukov Fern Gali 2010, 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{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 Huber:1981 provides a convex alternative to the use of Gaussian likelihoods for the concomitant estimation of location and scale parameters (cf. Owen:2007 and SS:2018b for a discussion). Section (ref) 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., 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 (Vapnik:1999). Another approach is to focus on flexible conditional CDF modeling, while discarding the monotonicity requirement. Distribution regression (ForesiPerrachi, 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; Akaike:1973, 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 (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 (Koenker:Bassett1978), as well as parametric location-scale formulations (He:1997, SS:2018a, MachadoSantosSilva:2019). Function-valued parameters lead to slower than root-$n$ conditional PDF estimation (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., LR:2007). Our approach further applies to conditional quantile function (CQF) estimation and is thus also related to the kernel-based estimator of 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 (BV:2004) for the one-step estimation of conditional PDFs and CDFs at each sample points. Compared to dual regression and its generalization (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) introduces our modeling framework. Section (ref) gives results under misspecification. Section (ref) contains estimation and inference results, and duality theory is derived in Section (ref). Section (ref) illustrates our methods, Section (ref) gives extensions and Section (ref) concludes. Proofs of Theorems (ref)-(ref) and (ref) are given in the Appendix. The Supplemental Material (SS:2025) contains proofs of Theorems (ref)-(ref), technical results, implementation details, and results of numerical simulations.
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}$,
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.
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.
The GT $g(Y,X)$ in ((ref)) 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)).\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)) 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., DeVore:1977 for splines and 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)) preserves approximation quality. } When the nonconstant components of $W(X)$ and $s(Y)$ are specified as nonnegative spline functions (CS:1966, Ramsay:1988), we refer to the implied representations as `Spline-Spline models'. We note that we do not impose $b_{0}>0$.
Model ((ref)) 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
The resulting log conditional density forms the basis of our approach:
Viewed as an equation in the GT $b_{0}'T(X,Y)$, ((ref)) defines the problem of characterizing GTs of the form ((ref)), ruling out implied conditional PDFs that are negative or zero with positive probability, and hence also nonmonotone conditional CDFs.
Formulation ((ref)) 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
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)) avoids the scaling factor and yields first-order conditions linear in $b_{0}'T(X,Y)$ (cf. ((ref)) and ((ref)) below). Formulation ((ref)) also differs from support vector methods that define the problem of characterizing conditional PDFs as solving the 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 WGSVVW:1999). Compared to both ((ref)) and ((ref)), shape constraints from economic theory are also easier to impose using ((ref)), with GT in closed-form.\footnote{Shape constraints often apply to CDFs or CQFs (e.g., BKM:2014, CW:2017), and they easily translate into restrictions on the shape of GTs. For example, under Assumption (ref), 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)). 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$.}
For $\Theta=\{b\in\mathbb{R}^{JK}:\Pr[b't(X,Y)>0]=1\}$, we define the population objective function
This criterion introduces a natural logarithmic barrier function (e.g., 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)).
We characterize the shape and properties of $Q(b)$ under the following 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) requires $Y$ to have finite second moment. The moment conditions in Assumption (ref) are also sufficient for the second-derivative matrix of $Q(b)$,
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.
By standard ML theory (e.g., Newey:McFadden:1994, p. 2124), this result implies identification of $b_{0}$, with $b_{0}$ being the only solution to the first-order conditions
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 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)). Theorem (ref) then implies that $b_{0}=(0,1,0_{JK-2})'$ uniquely solves ((ref)):
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)) 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)), conditions ((ref)) thus characterize a transformation of $Y$ to Gaussianity at each $X$ value. Hence, Theorem (ref) has the following testable implications for model ((ref)).
The general modeling of $F_{Y|X}(Y|X)$ can be done indirectly by specifying a representation for $Y$ given $X$,
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)$:
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)) 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., Chesher:2003, 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))-((ref)) that preserves nonlinearity and nonseparability is to model $F_{Y|X}(Y|X)$ directly as
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)) or ((ref)) 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)),
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)) 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$ (Owen:2007, SS:2018b). In contrast, a major advantage of representation ((ref)) 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). This formulation also leads to well-defined representations for $F_{Y|X}(Y|X)$ under misspecification.
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)) 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.
Under Assumption (ref) the objective function $Q(b)$ is continuous and strictly concave over $\Theta$, and hence admits at most one maximizer. Assumption (ref) 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^{*}$.
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 (Akaike:1973, 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
hold with probability one, for all $b\in\Theta$. Given a specified dictionary such that ((ref)) holds, Theorem (ref) implies that the approximation $f^{*}(Y,X)$ in ((ref)) 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.
Under the boundary conditions ((ref)), 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)) 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)),
implies that $\beta_{2}(X)>0$ is necessary for the boundary conditions ((ref)) 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)) implies that $\beta_{2}(X)>0$ is also sufficient for ((ref)). 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 (CS:1966, Ramsay:1988).} since the derivative then reduces to $\beta_{2}(X)$ outside this region. The boundary conditions ((ref)) then effectively hold under a location-scale restriction in the tails of the distribution of $Y$ given $X$. We also note that ((ref)) always holds for $J=2$ since the derivative condition is $\beta_{2}(X)>0$ in that particular case.
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
We derive the asymptotic properties of $\hat{b}$ under the following assumptions.
Assumption (ref)(i) can be replaced with the condition that $\{(y_{i},x_{i})\}_{i=1}^{n}$ is stationary and ergodic (Newey:McFadden:1994). Assumption (ref)(ii) is needed for consistent estimation of the asymptotic variance-covariance matrix of $\widehat{b}$.
Theorem (ref)(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)(ii) holds without compactness of $\Theta$, by concavity of the objective (e.g., 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., Newey:McFadden:1994). A valid method for model selection is adaptive Lasso ML (Zou:2006, LGF:2012, HN:2020). The selected model is a sparse KLIC optimal approximation for $g(Y,X)$.
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.
Theorem (ref) provides an estimator for the asymptotic variance-covariance matrix $\Xi$.
Our approach is related to distribution regression estimation of conditional CDFs, for the (probit) model
where $\beta(y)$ a vector of unknown functions. One sense in which our approach is flexible is that model ((ref)) can approximate $F_{Y|X}(y|X)$ in ((ref)) arbitrarily well for $S(y)$ rich enough.\footnote{Formally, write $S^{J}(Y)=S(Y)$ and suppose that Assumptions (ref) 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 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)) 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)) is efficient. When ((ref)) 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)) and ((ref)) are particular cases of conditional transformation models of the form $F(g(Y,X))=F(\sum_{l=1}^{L}g_{l}(Y,X))$ (HKB:2014), for some specified CDF $F$. For estimation, taking this additive structure as a starting point leads to considering the restricted form
(HKB:2014), with monotonicity constraints on each partial transformation $g_{l}$ (HMB:2018, 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)) gives a flexible generalization of ((ref)) 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)) 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., He:1997, DV:2008, Chernozhukov Fern Gali 2010, YT:2017, SS:2018a).} When this occurs, rearrangement methods provide improvements (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) avoids the need for correcting estimates (GHI:2003) and multiplicity of inverses (cf. 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 Matzkin:2003, the direct specification of the distribution of $e$ through the log density in ((ref)) 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) (Goodfellow:2014) used in Econometrics for simulation of complex datasets (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 (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.
Considerable computational advantages accrue from our ML criterion where the GT enters in closed-form. To ((ref)) corresponds a dual formulation that can be cast into the modern convex programming framework (BV:2004). We derive the dual problem and establish the properties of the dual solutions.
It is useful to compare problem ((ref))-((ref)) 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 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))-((ref)) 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)), 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))-((ref)) compared to the alternative generalized dual regression characterization of conditional CDFs and CQFs (SS:2018a) for which the mathematical program is of the form
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)) has first-order conditions
where $\widehat{b}$ is the Lagrange multiplier vector for the constraints in ((ref)), but where the solution is now determined by a system of $n$ nonlinear equations instead of having a closed-form expression as in ((ref)). This further illustrates the benefits accruing from closed-form modeling of $e=g(Y,X)$, compared to direct modeling of $y_{i}$ in ((ref)).
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 (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 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
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)) 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)$,
Following 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) of the Supplemental material. All computational procedures can be implemented in the software R (R:2024) using open source packages for convex optimization such as CVX, and its R implementation CVXR (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 (KS:1991) or monotone regression function estimation (Ramsay:1988). Thus, ((ref))-((ref)) is a flexible specification for the gender wage gap, allowing for both observed and unobserved heterogeneity.
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, 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.
Figures (ref)-(ref) 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)(A), and are reflected by the two inflection points in the corresponding CDFs in Figure (ref)(B), and by the large gap between the low and high CQFs in Figure (ref)(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)), with a test statistic of $316>16.9$, the critical value at the $5\%$ level.}
Figure (ref) 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.
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))) $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)) 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)) 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:
where $\lambda(\cdot)$ is the standard Logistic density. In the Gaussian version of ((ref)) 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.
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)$ (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}$.
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)) is the recursive formulation
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
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 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}. \]
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.