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.
59,747 characters · 15 sections · 20 citation commands
Dual regression
keywords: Conditional distribution; Duality; Monotonicity; Quantile regression; Method of moments; Mathematical programming; Convex approximation.
Let $Y$ be a continuously distributed random variable and $X$ a random vector. Then the conditional distribution function of $Y$ given $X$, written $U=F_{Y\mid X}(Y\mid X)$, has three properties: $U$ is standard uniform, $U$ is independent of $X$, and $F_{Y\mid X}(y\mid x)$ is strictly increasing in $y$ for any value $x$ of $X$. We will refer to these three properties as uniformity, independence and monotonicity. For some specified mean zero and unit variance distribution function $F$ with support the real line and inverse function $F^{-1}$, define $\varepsilon=F^{-1}\{F_{Y\mid X}(Y\mid X)\}$. Then $\varepsilon$ satisfies independence and monotonicity, has distribution $F$, and is transformed to uniformity by taking $U=F(\varepsilon)$.
The use of dual is thus motivated by the general observation that the estimation problem for a conditional distribution function $F_{Y\mid X}$ indexed by a parameter $\theta$ is usually formulated in terms of a procedure that obtains $\theta$ directly and $F_{Y\mid X}$ as a byproduct that follows from a calculation from the representation evaluated at a specific value of $\theta$. A classical example is the linear location shift model $F_{Y\mid X}(y_{i}\mid x_{i})=F\{(y_{i}-\beta^\textrm{T} x_{i})/\sigma\}$, for which the parameter vector $\theta=(\beta,\sigma)^{\textrm{T}}$ needs to be estimated in order to obtain the $n$ values $\varepsilon_{i}=(y_{i}-\beta^\textrm{T} x_{i})/\sigma$. Here we turn that process around, obtaining $\varepsilon_{i}$ first (from a mathematical programming problem) and backing out $\theta$ afterwards, if at all.
In its simplest form, dual regression augments the median regression dual programming problem KB:1978 with global second moment orthogonality constraints, while expanding the support of parameter values from the unit interval to the real line. Adding further global orthogonality constraints gives rise to a sequence of augmented, generalized dual regression programs. Although each of these programs seeks only to find the $n$ values $\varepsilon_{i}=F^{-1}(u_{i})$, their first-order conditions show that the assignment of these $n$ values corresponds to a sequence of augmented location-scale representations, the simplest element of which is a linear heteroscedastic model. Moreover, their second-order conditions are equivalent to monotonicity, so optimal dual regression solutions are free of quantile-crossing problems.
To each element of the sequence of dual programs corresponds a convex primal problem, both nontrivial to determine and difficult to implement, the convexity of which guarantees uniqueness of optimal dual regression solutions. For a given specification of $F_{Y\mid X}(Y\mid X)$, the first-order conditions of the corresponding primal problem also describe necessary and sufficient conditions for independence of the associated dual solutions. Thus our dual formulation reveals a sequence of convex optimization problems, gives a feasible and direct implementation of each of them, and uniquely characterizes the family of associated globally monotone representations, which can then be used as complete estimates of a flexible class of conditional distribution functions.
We introduce the basic principles underlying our general method by first providing a new characterization of the conditional distribution function $F_{Y\mid X}(Y\mid X)$ associated with the linear location-scale model
where $X$ is a $K\times1$ vector of explanatory variables including an intercept, and $F$ a mean zero and unit variance cumulative distribution function over the real line.
Suppose that we observe a sample of $n$ identically and independently distributed realizations $\{(y_{i},x_{i})\}_{i=1}^{n}$ generated according to model ((ref)). The primary population target of our analysis is \[ \varepsilon_{i}=\frac{y_{i}-\beta_{1}^\textrm{T} x_{i}}{\beta_{2}^\textrm{T} x_{i}}=F^{-1}\{F_{Y\mid X}(y_{i}\mid x_{i})\} \quad(i=1,\ldots,n), \] knowledge of which is equivalent to knowledge of the $n$ values $F_{Y\mid X}(y_{i}\mid x_{i})$ up to the monotone transformation $F$.
Let $\lambda=(\lambda_{1},\lambda_{2})^{\textrm{T}}\in\mathbb{R}^{2\times K}$ and $e_{o}\in\mathbb{R}^{n}$ satisfy the system of $n$ equations and $n$ inequality constraints
where $e_{o}$ further satisfies the $2\times K$ orthogonality conditions $\sum_{i=1}^{n}x_{i}e_{oi}=0$ and $\sum_{i=1}^{n}x_{i}(e_{oi}^{2}-1)=0$. Since $x_{i}$ includes an intercept, the sample moments of $e_{o}$ and $e_{o}^{2}$ are $0$ and $1$, and $e_{o}$ and $e_{o}^{2}$ are orthogonal to each column of the $n\times K$ matrix $(x_{1},\ldots,x_{n})^{\textrm{T}}$ of explanatory variables. We propose a characterization of the sequence of vectors $e_{o}$ that satisfy representation ((ref)) and the associated orthogonality constraints for each $n$. The corresponding sequence of empirical distribution functions then provides an asymptotically valid characterization of the conditional distribution function $F_{Y\mid X}(Y\mid X)$ corresponding to the data-generating process ((ref)). As a by-product of this approach, we simultaneously obtain a characterization of the parameter vector $\lambda$ in ((ref)), which then provides a consistent estimator of the population parameter $\beta$ in ((ref)).
For each $x_{i}$, with the scale function $\lambda_{2}^\textrm{T} x_{i} > 0$, $y_{i}$ is an increasing function of $e_{oi}$, and to representation ((ref)) corresponds a convex function \[ C(x_{i},e_{oi},\lambda)=\int_{0}^{e_{oi}}\{\lambda_{1}^\textrm{T} x_{i}+(\lambda_{2}^\textrm{T} x_{i})s\}ds=(\lambda_{1}^\textrm{T} x_{i})e_{oi}+\frac{1}{2}(\lambda_{2}^\textrm{T} x_{i})e_{oi}^{2} \quad(e_{oi}\in\mathbb{R}), \] and whose quadratic form corresponds to a location-scale representation for $F_{Y\mid X}(Y\mid X)$. Letting $y$ be the $n\times1$ vector of dependent variable values, and assuming knowledge of $\lambda$ and $e_{o}$, we consider assigning a value $e_{i}$ to each observation in the sample by maximizing the correlation between $y$ and $e=(e_{1},\ldots,e_{n})^{\textrm{T}}$ subject to a constraint that embodies the properties of $e_{o}$:
Problem ((ref)) describes the assignment of $e$ values to $y$ values in a sample generated according to a location-scale model, and it admits $e=e_{o}$ as its only solution. Since $e_{o}$ and $\lambda$ are unknown, the assignment problem ((ref)) is infeasible: we thus introduce the equivalent, feasible formulation \[ \textrm{(D)}\quad\max_{e\in\mathbb{R}^{n}}\left\{ y^{\textrm{T}}e : \sum_{i=1}^{n}x_{i}e_{i} =0, \; \frac{1}{2}\sum_{i=1}^{n}x_{i}(e_{i}^{2}-1) =0\right\}, \] the dual regression program.
The solution to (D) is easily found from the Lagrangian \[ \mathscr{L}=\sum_{i=1}^{n}y_{i}e_{i}-\lambda_{1}\sum_{i=1}^{n}x_{i}e_{i}-\frac{1}{2}\lambda_{2}\sum_{i=1}^{n}x_{i}(e_{i}^{2}-1). \] Differentiating with respect to $e_{i},$ we obtain $n$ first-order conditions:
Upon rearranging we obtain the closed-form solution
which is of the location-scale form $e_{i}=\{y_{i}-\mu(x_{i})\}/\sigma(x_{i})$, with $\mu(x_{i})$ and $\sigma(x_{i})$ linear in $x_{i}$.
Another view is obtained by writing the first-order conditions as
a linear location-scale representation, with corresponding quantile regression representation
Program (D) thus provides a complete characterization of linear representations of the form ((ref)) and ((ref)), as they arise from its first-order conditions. Moreover, the parameters of these representations are the Lagrange multipliers $\lambda_{1}$ and $\lambda_{2}$ of an optimization problem with solution $e=e_{o}$.
The quantile regression representation of the first-order conditions of (D) sheds additional light on the monotonicity property of dual regression solutions, when there are no repeated $X$ values. For $u,u'\in(0,1)$, $u'>u$, the no-crossing property of conditional quantiles requires \[ \beta(u')^\textrm{T} x_{i}-\beta(u)^\textrm{T} x_{i}>0,\quad(i=1,\ldots,n), \] which is satisfied if $\lambda_{2}^\textrm{T} x_{i}$ is strictly positive for each $i$, and coincides with the $n$ second-order conditions of program (D): \[ \frac{\partial^{2}\mathscr{L}}{\partial e_{i}\partial e_{i}}=-\lambda_{2}^\textrm{T} x_{i}<0,\quad(i=1,\ldots,n). \] Therefore, an optimal $e$ solution that violates the monotonicity property is ruled out by the requirement that for an observation with $X$ value $x_{i}$, the ordering of the $Y$ values $\beta(u')^\textrm{T} x_{i}$ and $\beta(u)^\textrm{T} x_{i}$ must correspond to the ordering of the $u$ values. Hence the correlation criterion of system (D) suffices to impose monotonicity, with optimality of a solution then being equivalent to monotonicity at the $n$ sample points. Dual regression is thus able to incorporate this property in the estimation procedure, which facilitates extrapolation beyond the empirical support of $X$, and yields significant finite-sample improvements in the estimation of conditional quantile functions as illustrated by our simulations in $\mathsection$(ref).
By Lagrangian duality arguments (BV:2004, Chapter 5), the objective function of the dual of problem (D) is \[ Q_{n}(\lambda)=\sup_{e\in\mathbb{R}^{n}}y^{\textrm{T}}e-\sum_{i=1}^{n}\left\{ C\left(x_{i},e_{i},\lambda\right)-C\left(x_{i},e_{oi},\lambda\right)\right\} , \] defined for all $\lambda\in\Lambda_{0}$, where $\Lambda_{0}=\Lambda_{1}\times\Lambda_{2}$, with $\Lambda_{1}=\mathbb{R}^{K}$ and $\Lambda_{2}=\{\lambda_{2}\in\mathbb{R}^{K}:\;\inf_{i\leq n}\lambda_{2}^\textrm{T} x_{i}>0\}$. Under the conditions of Theorem (ref) below, $Q_{n}(\lambda)$ has a closed-form expression, is strictly convex over $\Lambda_{0}$, and minimizing $Q_{n}(\lambda)$ over $\Lambda_{0}$ is equivalent to solving (D). Given a vector $\omega\in\mathbb{R}^{n}$, we let $\textrm{diag}(\omega_{i})$ denote the $n\times n$ diagonal matrix with diagonal elements $\omega_{1},\ldots,\omega_{n}$.
Theorem (ref) summarizes our finite-sample analysis of dual regression. The proofs of all formal results in the paper are given in the Supplementary Material.
Theorem (ref) establishes formal duality of our initial assignment problem under first and second moment orthogonality constraints and the global $M$-estimation problem (P). Convexity of (P) guarantees that to a unique assignment of $e$ values corresponds a unique linear representation of the form ((ref)). Uniqueness further implies that if $e_{o}$ satisfies independence, then the orthogonality conditions in ((ref)) are both necessary and sufficient for the dual solution $e^{*}$ to satisfy independence.
The primal problem (P) is a locally heteroscedastic generalization of a simultaneous location-scale estimator proposed by Huber:1981 and further analyzed in Owen:2001. The linear heteroscedastic model of equation ((ref)) has been previously encountered in the quantile regression literature: see Koenker:Zhao:1994 and He:1997. The former consider the efficient estimation of ((ref)) via $L$-estimation while the latter develops a restricted quantile regression method that prevents quantile crossing. Compared to these quantile-based methods, dual regression trades local estimation and the convenient linear programming formulation of quantile regression for simultaneous estimation of location and scale parameters.
The dual problem of the linear $0\cdot5$ quantile regression of $Y$ on $X$ is (cf., Koenker, 2005, p. 87, equation 3.12):
The solution to problem ((ref)) produces values of $u$ that are largely 0 and 1, with $K$ sample points being assigned $u$ values that are neither 0 nor 1. The points that are assigned 1 fall above the median quantile regression; the points receiving 0's fall below; and the remaining points fall on the median quantile regression plane. One direction of extension of ((ref)) is to replace the 1/2 with values $\alpha$ that fall between 0 and 1 to obtain the $\alpha$ quantile regression.
Another extension is to augment problem ((ref)) by adding $K$ more constraints:
It is apparent that the solution to ((ref)) does not satisfy ((ref)): the variance of $u$ around 0 in the solution to ((ref)) is approximately $1/2$, not $1/3$. To satisfy program ((ref)), the $u$'s have to be moved off $0$ and $1$. Since $x_{i}$ contains an intercept, the sample moments of $u$ and $u^{2}$ will be $1/2$ and $1/3$; $u$ and $u^{2}$ will be orthogonal to the columns of the matrix $(x_{1},\ldots,x_{n})^{\textrm{T}}$, relations that are necessary but not sufficient for uniformity and independence.
Both systems ((ref)) and ((ref)) demand monotonicity by maximally correlating $y$ and $u$. A violation of monotonicity requires there to be two observations that share the same $X$ values but have different $y$ values, with the lower of the two $y$ values having the weakly higher value of $u$. However, a solution characterized by such a violation could be improved upon by exchanging the $u$ assignments. Thus violation of monotonicity in program ((ref)) arises because the set of admissible exchanges in $u$ assignments is overly restricted: ((ref)) is dual to a linear program well-known to have solutions at which $K$ observations are interpolated when $K$ parameters are being estimated, i.e., the hyperplanes obtained by regression quantiles must interpolate $K$ observations.
By reformulating program ((ref)) into a constrained optimization problem over $\mathbb{R}^{n}$, program (D) further expands the set of admissible exchanges in $u$ assignments, since $u$ is restricted to $[0,1]^{n}$. Doing this, the problem corresponding to ((ref)) becomes the dual regression program (D), where $e$ can take on any real value. It is then natural to take $u_{i}^{*}=F_{n}(e_{i}^{*}),$ the empirical cumulative distribution function of the dual regression solution $e^{*}$, thereby imposing uniformity to high precision even at small $n$.
The dual regression characterization of location-scale conditional distribution functions via the monotonicity element, the objective, and the independence element, the constraints, can be exploited to characterize more flexible representations. Similarly to the approach introduced in $\mathsection$(ref), we first analyze the infeasible assignment problem for a general representation of the stochastic structure of $Y$ conditional on $X$:
where $F$ is a specified cumulative distribution function with support the real line, and for each value $x$ of $X$, the derivative $H_{x}'(\varepsilon)$ of $H_{x}(\varepsilon)$ is strictly positive. Representation ((ref)) always exists with $H_{x}$ defined as the composition of the conditional quantile function of $Y$ given $X=x$ and the distribution function $F$.
To each monotone function $H_{x}$ also corresponds a convex function $\widetilde{H}_{x}$ defined as \[ \widetilde{H}_{x}(e)\equiv\int_{0}^{e}H_{x}(s)ds \quad(e\in\mathbb{R}). \] The monotonicity of $H_{x}(\varepsilon)$ guarantees the convexity of $\widetilde{H}_{x}(\varepsilon)$. The slope of this function gives the value of $Y$ corresponding to a value $e$ of $\varepsilon$ at $X=x$. Thus $F_{Y\mid X}(Y\mid X)$ corresponds to a collection of convex functions, with one element of this collection for each value of $X,$ together with a single random variable whose distribution is common to all the convex functions.
Equipped with $\widetilde{H}_{X}$, suppose we are tasked with assigning a value $e_{i}$ to each of the $n$ realizations $\{(y_{i},x_{i})\}_{i=1}^{n}$. Then, for $S_{n}=\sum_{i=1}^{n}\widetilde{H}_{x_{i}}(\varepsilon_{i})$, solving the infeasible problem \[ \textrm{(IGD)} \quad \max_{e\in\mathbb{R}^{n}} \left\{ y^{\textrm{T}}e : \sum_{i=1}^{n}\widetilde{H}_{x_{i}}(e_{i})=S_{n}\right\} , \] generates the correct $y-e$ assignment: writing the Lagrangian \[ \mathscr{L}=y^{\textrm{T}}e-\Lambda\left\{ \sum_{i=1}^{n}\widetilde{H}_{x_{i}}(e_{i})-S_{n}\right\} , \] the $n$ associated first-order conditions are
Problem (IGD) is infeasible because neither $\widetilde{H}_{x_{i}}$ nor $S_{n}$ is known. However, Theorem (ref) motivates a feasible approach once $H_{X}$ and $F$ are specified. Denote the components of $X$ without the intercept by $\widetilde{X}$, so that $X=(1,\widetilde{X})^{\textrm{T}}$. Without loss of generality, let $\widetilde{X}$ be centered, denoted $\widetilde{X}^{c}$, and let $X^{c}=(1,\widetilde{X}^{c})^{\textrm{T}}$. With $h_{1}(\varepsilon)=1$ and $h_{2}(\varepsilon)=\varepsilon$, we specify $H_{X}$ by a linear combination of $J$ basis functions $h(\varepsilon)=\{h_{1}(\varepsilon),\ldots,h_{J}(\varepsilon)\}^{\textrm{T}}$, the coefficients of which depend on $X$:
and we assume that $H_{X}$ is linear in $X$ and set:
Finally, we specify a zero mean and unit variance distribution for $\varepsilon$ by imposing $E(\varepsilon)=0$ and $E(\varepsilon^{2}-1)/2=0$, and setting $\alpha_{j}=0$ for $j=3,\ldots,J$, in ((ref)).
With $\alpha_{2}+\beta_{2}^\textrm{T}\widetilde{X}^{c}>0$, our normalization and ((ref))--((ref)) together yield the augmented, generalized dual regression model
Equation ((ref)) admits of the following interpretation. When $\widetilde{X}^{c}=0,$ $Y=\alpha_{1}+\alpha_{2}\varepsilon$ and $\varepsilon=(Y-\alpha_{1})/\alpha_{2}$, so that $\varepsilon$ is just a re-scaled version of the distribution of $Y$ at $\widetilde{X}^{c}=0$. Since $\varepsilon$ is independent of $X$, transformations of this shape of $\varepsilon$ must suffice to produce $Y$ at other values of $X$. The first two transformations, $\beta_{1}^\textrm{T}\widetilde{X}^{c}$ and $(\beta_{2}^\textrm{T}\widetilde{X}^{c})\varepsilon$, are translations of location and scale which do not essentially affect the shape of $Y$'s response to changes in $\varepsilon$ at all. The additional terms $(\beta_{j}^\textrm{T}\widetilde{X}^{c})h_{j}(\varepsilon)$ achieve that end.
Suppose that we observe a sample of $n$ identically and independently distributed realizations $\{(y_{i},x_{i})\}_{i=1}^{n}$ generated according to model ((ref)). Define $x_{ij}^{c}=x_{i}^{c}$ for $j=1,2$, and $x_{ij}^{c}=\widetilde{x}_{i}^{c}$ for $j=3,\ldots,J$, and let $(\gamma,\lambda)\in\mathbb{R}^{2+J(K-1)}$ and $e_{o}\in\mathbb{R}^{n}$ satisfy the system of $n$ equations and $2n$ inequality constraints
where $e_{o}$ further satisfies $\sum_{i=1}^{n}x_{ij}^{c}\widetilde{h}_{j}(e_{oi})=0$ $(j=1,\ldots,J)$, with $\widetilde{h}_{1}(e_{oi})=e_{oi}$, $\widetilde{h}_{2}(e_{oi})=(e_{oi}^{2}-1)/2$, and $\widetilde{h}_{j}(e_{oi})=\int_{0}^{e_{oi}}h_{j}(s)ds$ $(j=3,\ldots,J)$. These relations reduce to the linear heteroscedastic representation of $\mathsection$(ref) for $J=2$, and impose that $e_0$ be a zero mean and unit variance vector satisfying the augmented set of orthogonality conditions $\sum_{i=1}^{n}\widetilde{x}_{i}^{c}\widetilde{h}_{j}(e_{oi})=0$ $(j=1,\ldots,J)$. The sequence of vectors $e_{o}$ that satisfies the generalized dual regression representation ((ref)) as well as the associated orthogonality constraints for each $n$ then provides an asymptotically valid characterization of the data-generating process ((ref)).
Each element of this sequence is characterized by the assignment problem
where $\widetilde{H}_{x_{i}}(e_{i};\theta)=\int_{0}^{e_{i}}H_{x_{i}}(s;\theta)ds$, and $\theta=(\theta_{1},\ldots,\theta_{J})^{\textrm{T}}$, with $\theta_{j}=(\gamma_{j},\lambda_{j})^{\textrm{T}}\in\mathbb{R}^{K}$ for $j=1,2$, and $\theta_{j}=\lambda_{j}\in\mathbb{R}^{K-1}$ for $j=3,\ldots,J$. Since $e_{0}$ and $\theta$ are unknown, problem ((ref)) is infeasible; we thus formulate an equivalent, feasible implementation of problem (IGD): \[ \textrm{(GD)}\quad\max_{e\in\mathbb{R}^{n}} \left\{ y^{\textrm{T}}e : \sum_{i=1}^{n}x_{ij}^{c}\widetilde{h}_{j}(e_{i})=0\quad(j=1,\ldots,J)\right\} , \] the generalized dual regression program. (GD) then uniquely characterizes representation ((ref)).
In order to state the properties of (GD) formally, we define the parameter space $\Theta_{n}$, which specifies parameter values compatible with monotone representations: \[ \Theta_{n}=\left\{ \theta\in\Theta_{0,n}:\textrm{there exists }e\in\mathbb{R}^{n}:y_{i}=H_{x_{i}}(e_{i};\theta)\textrm{ and }\inf_{e\in\mathbb{R}}H'_{x_{i}}(e;\theta)>0\quad(i=1,\ldots,n)\right\} , \] with $\Theta_{0,n}=\{\theta\in\mathbb{R}^{2+J(K-1)}:\;\inf_{i\leq n}\theta_{2}^\textrm{T} x_{i}^{c}>0\}$. For $\theta\in\Theta_{n}$, let $e(y_{i},x_{i},\theta)$ denote the inverse function of $H_{x_{i}}(e_{i};\theta)$, which is well-defined for each $x_{i}$. We assume that the basis functions $h$ and the pair $(\theta,e_{0})$ satisfy the following conditions.
Let $\phi(\theta)=[H_{x_{1}}'\{e(y_{1},x_{1},\theta);\theta\},\dots,H_{x_{n}}'\{e(y_{n},x_{n},\theta);\theta\}]^{\textrm{T}}$. Theorem (ref) summarizes our finite-sample analysis of generalized dual regression.
Problem (GD) augments the set of orthogonality constraints in (D) and generates increasingly flexible representations of the form ((ref)). For each element of this sequence, (GD) then provides a feasible formulation of the generalized $y-e$ assignment problem (IGD) with optimality condition $-H'_{x_{i}}(e_{i}^{*};\theta^{*})<0$ equivalent to monotonicity. Theorem (ref) also states the form of the corresponding primal problem, whose convexity guarantees that (GD) and (GP) uniquely and equivalently characterize representation ((ref)). Uniqueness further implies that if $e_{o}$ satisfies independence, then the orthogonality conditions in ((ref)) are both necessary and sufficient for the dual solution $e^{*}$ to satisfy independence as well. Theorem (ref) thus characterizes and establishes the duality between specification of orthogonality constraints on $e$ and specification of a globally monotone representation for $Y$ conditional on $X$.
Formally, (GP) is the restriction of the dual of (GD) to $\Theta_{n}$. The existence Condition (ref) and the form of (GD) optimality conditions together ensure that (GD) does not admit a global maximum with associated multipliers outside $\Theta_{n}$. Implementing (GP) thus requires the imposition of inequality constraints with $e_{i}$ only implicitly defined in the specification of (GP) for $J>2$, and problem (GD) therefore provides a greatly simplified dual implementation. The special case of dual regression corresponds to $J=2$, and imposing $\sum_{i=1}^{n}\widetilde{h}_{j}(e_{i})=0$, for $j=1,2$, is a normalization. The simple basis $\{e_{i},(e_{i}^{2}-1)/2\}$ is obviously impoverished for the space of all convex functions, although quite practical for many applications once the flexibility in the distribution of $e$ is taken into account.
An alternative approach is to specify $F$ to a known distribution, and alter representation ((ref)) and the corresponding problem accordingly. If $F$ is specified to be the standard uniform distribution, then ((ref)) in $\mathsection$(ref) can be generalized as
For $u_{i}\in[0,1]$, let $m^{J}(u_{i})=\{m_{J1}(u_{i}),\ldots,m_{JJ}(u_{i})\}^{\textrm{T}}$, with $m_{Jj}(u_{i})=j^{-1}\{u_{i}^{j}-(j+1)^{-1}\}$. With $\otimes$ denoting the Kronecker product, the large-sample form of program ((ref)) is
Letting $J$ increase, both the distributional and the orthogonality constraints get strengthened. Because $X$ includes an intercept, the distribution of $U_{J}$ approaches uniformity, while simultaneously satisfying an increasing sequence of orthogonality constraints. In the limit, a uniformly distributed random variable $U$ satisfying the full set of orthogonality constraints is thus specified. Since $E\{X\otimes m^{J}(U)\}=0$ for all $J$ is equivalent to the mean-independence property $E(X\mid U)=E(X)$ and the uniformity constraint $U\sim U(0,1)$, in the large $J$ limit program ((ref)) coincides with the scalar quantile regression problem proposed in independent work by CCG:2016 (cf., equation 19, p. 1180)
which provides an optimal transport formulation of quantile regression (we are grateful to an anonymous referee for highlighting this connection). Program ((ref)) is directly amenable to a linear programming implementation which maintains and exploits the full specification of the marginal distribution of $U$ to a known distribution, whereas ((ref)) provides a sequential nonlinear programming characterization of $U$ which relaxes uniformity for finite $n$ and $J$.
For $e_{i}\in\mathbb{R}$, let $\widetilde{h}^{J}(e_{i})=\{\widetilde{h}_{1}(e_{i}),\ldots,\widetilde{h}_{J}(e_{i})\}^{\textrm{T}}$. The large-sample form of program (GD) is
Program ((ref)) relaxes the support constraint in ((ref)) and only specifies first and second moments of $e_{J}$, while the centering of $X$ ensures that this is sufficient for $e_{J}$ to be uniquely determined. The empirical distribution of solutions of the finite-sample analog (GD) of ((ref)) then provides an asymptotically valid characterization of the distribution of $e_{J}$.
Letting $J$ increase, orthogonality constraints in ((ref)) are strengthened, and $e_{J}$ gets closer and closer to satisfying the mean-independence property $E(\widetilde{X}^{c}\mid e_{J})=0$. It follows that for $J$ large enough, ((ref)) is equivalent to
the limiting generalized dual regression problem. Theorem (ref) summarizes this discussion.
We apply our framework to the estimation of a $J$--term generalized dual regression model of the form ((ref)). Denote the support of $X$ by $\mathcal{X}$, and, for some finite constant $C_{\theta}$, define $\Theta_{0}=\{\theta\in\mathbb{R}^{2+J(K-1)}:||\theta||\leq C_{\theta}\,\textrm{and}\:\inf_{x\in\mathcal{X}}\theta_{2}^\textrm{T} x^{c}>0\}$. Letting $\mathcal{C}^{1}(\mathbb{R})$ denote the space of continuously differentiable functions on $\mathbb{R}$, define the space of strictly increasing functions indexed by $X$ values, $\mathcal{M}(X)=\{e_{X}\in\mathcal{C}^{1}(\mathbb{R}):\inf_{y\in\mathbb{R}}e'_{x}(y)>0\;\textrm{for all}\,x\in\mathcal{X}\}$. The large-sample analog of $\Theta_{n}$ is then the space of vectors in $\Theta_{0}$ such that there exists a corresponding optimal generalized dual regression representation: \[ \Theta=\left\{ \theta\in\Theta_{0}:\textrm{there exists }e_{X}\in\mathcal{M}(X)\,\textrm{with}\,\Pr[Y=H_{X}\{e_{X}(Y);\theta\}]=1\right\} . \] For any $\theta\in\Theta$, denote $e_{X}$ in $\mathcal{M}(X)$ such that $\Pr[Y=H_{X}\{e_{X}(Y);\theta\}]=1$ by $e(Y,X,\theta)$.
These conditions are used to establish existence and consistency of dual regression solutions, and Condition (ref)(ii) is needed for asymptotic normality of estimates of $\theta_{0}$. In view of uniqueness stated in part (iii) of Theorem (ref), these properties are shared by $\theta_{n}$ and $\theta^{*}$, which we denote by $\widehat{\theta}$ for notational simplicity. We also denote both $e_{i}^{*}$ and indirect estimates $e(y_{i},x_{i},\theta_{n})$, constructed after solving (GP), by $\hat{e}_{i}$, with empirical distribution function $F_{n}(e)=n^{-1}\sum_{i=1}^{n}1(\hat{e}_{i}\leq e)$, $e\in\mathbb{R}$. Furthermore, part (ii) of Theorem (ref) shows that while the solution $e^{*}$ is obtained directly by solving the mathematical program (GD), knowledge that the solution obeys representation ((ref)) can be exploited to write estimating equations for $\hat{\theta}$ in the form of system ((ref)). The computation of the asymptotic distribution of $\hat{\theta}$ follows from this characterization.
Knowledge of the statistical properties of $\hat{\theta}$ can be used to establish the limiting behaviour of the empirical distribution of $\hat{e}$. Define the empirical dual regression process \[ \mathbb{U}_{n}(e)=n^{1/2}\{F_{n}(e)-F(e)\}\quad(e\in\mathbb{R}). \] Theorem (ref) establishes weak convergence of the empirical distribution of $\hat{e}$ and the limiting behaviour of $\mathbb{U}_{n}$, accounting for its dependence on the distribution of $n^{1/2}(\hat{\theta}-\theta_{0})$.
Theorems (ref) and (ref) together establish that the pair $(\hat{\theta},\hat{e})$ provides an asymptotically valid characterization of the generalized dual regression representation specified in Condition (ref). When $\varepsilon$ is independent of $X$, Theorem (ref) further implies that the empirical distribution of $\hat{e}$ provides an asymptotically valid estimator of the conditional distribution of $Y$ given $X$. For $u\in(0,1)$, estimates of the $X$ coefficients in quantile regression form can then be constructed as $\sum_{j=1}^{J}\hat{\lambda}_{j}h_{j}\{F_{n}^{-1}(u)\}$, exploiting the structure of the conditional quantile function of $Y$ given $X$ implied by representation ((ref)). Theorem (ref) also establishes asymptotic normality of the empirical dual regression process. The form of the covariance function of $\mathbb{U}$ reflects the influence of imposing sample orthogonality constraints in (GD) on the empirical distribution of $e^{*}$, or equivalently, of sample variability of parameter estimates $\theta_{n}$ on the empirical distribution of $e(y_{i},x_{i},\theta_{n})$, as expected from the classical result of Durbin:1973.
Theorem (ref) can be applied to perform pointwise inference on the conditional distribution function of $Y$ conditional on $X$. However, simultaneous inference over regions of the joint support of $Y$ and $X$ is typically of interest in practice. Several approaches for uniform inference in the presence of non-pivotal limit processes have been considered in the literature (e.g., Koenker:Xiao:2002, and Parker:2013), including simulation methods CFG:2013. Extension of existing results to dual regression is beyond the scope of this paper but they provide a natural direction for future study of uniform inference on the empirical dual regression process.
The classical dataset collected by Engel consists of food expenditure and income measurements for 235 households, and has been studied by means of quantile regression methods Koenker:2005. We illustrate dual regression methods by estimating the statistical relationship between food expenditure and income, with household income as a single regressor and food expenditure as the outcome of interest.
We specify the vector of basis functions by means of trigonometric series. Alternative choices such as splines and shape-preserving wavelets (e.g., DeVore:1977, and CSS:2007). In order to choose $J$, we first implement program (GD) for $J=2$, which we then augment sequentially adding one pair of cosine and sine basis at a time, up to a representation of order $J=8$. At each step, we compute a Schwarz Information Criterion Schwarz:1978 applied to the primal generalized dual regression problem, exploiting the strong duality result of Theorem (ref) in order to compute its value as $y^{\textrm{T}}e^{*}+\{2+J(K-1)\}\log n$. Our procedure selects the location-scale representation $J=2$. In the Supplementary Material, we describe the procedure and report results from the augmented specifications, which show that our results are robust to the number of terms included. In order to test for the validity of the selected model, a complementary procedure that should be explored in future research is to test for independence of dual regression solutions and explanatory variables. The test for multivariate independence proposed by GQR:2007 constitutes an interesting starting point for such a development.
All computational procedures can be implemented in the software R R:2017 using open source software packages for nonlinear optimization such as Ipopt or Nlopt, and their R interface Ipoptr and Nloptr developed by Jelmer Ypma. Quantile regression procedures in the package quantreg have been used to carry our comparisons.
Figure (ref) illustrates our results and plots the estimated distribution of food expenditure conditional on household income. Estimates $\{u_{i}^{*}\}_{i=1}^{n}$, where $u_{i}^{*}=F_{n}(e_{i}^{*})$, are used in order to plot each observation in the $xyu$-space with predicted coordinates $(x_{i},y_{i},u_{i}^{*})$, and the solid lines give the $u$-level sets for a grid of values $\{0\cdot1,\ldots,0\cdot9\}$. Although nonstandard, this representation relates to standard quantile regression plots since the levels of the distribution function give the conditional quantiles of food expenditure for each value of income. These are the plotted shadow solid lines corresponding for each $u$ to dual regression estimates of conditional quantile functions of food expenditure given household income.
Figure (ref) shows that the predicted conditional distribution function obtained by dual regression is indeed endowed with all desired properties. Of particular interest is the fact that the estimated function is monotone in food expenditure. Also, our estimates satisfy some basic smoothness requirements across probability levels, in the food expenditure values. This feature does not typically characterize estimates of the conditional quantile process by quantile regression methods, as conditional quantile functions are then estimated sequentially and independently of each other. The decreasing slope of the distribution function across values of income provides evidence that the data indeed follow a heteroscedastic generating process. This is the distributional counterpart of quantile functions having increasing slope across probability levels, a feature characterizing the conditional quantile functions on the $xy$ plane and signalling increasing dispersion in food expenditure across household income values.
Figure (ref) gives the more familiar quantile regression plots. The plots presented show scatterplots of Engel's data as well as conditional quantile functions obtained by dual and quantile regression methods. The rescaled plots in the right panels of Fig. (ref) highlight some features of the two procedures. The fitted lines obtained from dual regression are not subject to crossing in this example, whereas several of the fitted quantile regression lines actually cross for small values of household income. Last, the more evenly spread dual regression conditional quantile functions illustrate the effect of specifying a functional form for the quantile regression coefficients, while preserving asymmetry in the conditional distribution of food expenditure.
Figure (ref) compares our estimates of intercept and income coefficients in quantile regression form, with estimates obtained by quantile regression. For interpretational purposes, we follow Koenker:2005 and estimate the functional coefficients after having recentered household income. This avoids having to interpret the intercept as food expenditure for households with zero income. After centering, the intercept coefficient can be interpreted as the $u$-th quantile of food expenditure for households with mean income. Fig. (ref) shows the estimated quantile regression coefficients as a function of $u$. It illustrates the fact that the flexible structure imposed by dual regression yields estimates that are indeed smoother than their quantile regression counterpart, the latter having a somewhat erratic behaviour around our estimates.
We give a brief summary of the results of a Monte Carlo simulation in order to assess the finite-sample properties dual regression. The data-generating process is
with parameter values calibrated to the empirical application, from which $4999$ samples are simulated. As a benchmark, we compare generalized dual regression estimates of the values $F_{Y\mid X}(y_{i}\mid x_{i})$ $(i=1,\ldots,n)$, to those obtained by applying the inversion procedure of CFG:2010 to the quantile regression process. For each simulation, the estimation and selection procedures are identical to those implemented in the empirical application.
Table (ref) reports a first set of results regarding the accuracy of conditional distribution function estimates. We report average estimation errors across simulations of dual regression and quantile regression estimators, respectively, and their ratio in percentage terms. Estimation errors are measured in $L^{p}$ norms $\left\Vert \cdot\right\Vert _{p}$, for $p=1,2$, and $\infty$, where for $f:\mathbb{R}\mapsto[0,1]$, $\left\Vert f\right\Vert _{p}=\left\{ \int_{\mathbb{R}}\left|f(s)\right|^{p}ds\right\} ^{1/p}$, and are computed with $e^{*}$ the solutions to the selected generalized dual regression program. Correct model selection ranges from $75\%$ of the simulations for $n=100$ to $90\%$ for $n=1000$, providing encouraging evidence about the validity of the proposed criterion. The results show that for this setup our estimates systematically outperform quantile regression-based estimates, with the spread in performance increasing with sample size. Whereas the reduction in average estimation error is between $8\%$ and $17\%$, depending on the norm, for $n=235$, estimation error is reduced up to $30\%$ when $n=1000$. The larger reduction in average errors in $L^{\infty}$ norm reflects the higher accuracy in estimation of extreme parts of the distribution.
In the Supplementary Material, we describe the experiment in detail, and report results on estimation of quantile regression coefficients and the distribution of selected models across simulations. We also include additional simulations that illustrate the empirical performance of dual regression with multiple covariates and show that it performs well relative to the noncrossing quantile regression method proposed by BRW:2010.
If we designate problems such as (D) and (GD) as already dual, then their solutions reveal a corresponding primal. Typically, the Lagrange multipliers of the dual appear as parameters in the primal, and the primal has an interpretation as a data-generating process. So perhaps not surprisingly the constraints on the construction of the stochastic elements have shadow values that are parameters of a data-generating representation. In this way the relation between identification and estimation is made perspicuous: a parameter of the data-generating process is the Lagrange multiplier of a specific constraint on the construction of the stochastic element, so to specify that some parameters are non-zero and others are zero is to say that some constraints are in the large-sample limit binding and others are not.
Another way of expressing this is to say that when a primal corresponds to the data-generating process, additional moment conditions are superfluous: they will in the limit attract Lagrange multiplier values of zero and consequently not affect the value of the program nor the solution. In a sense, this is obvious: the parameters of the primal can typically be identified and estimated through an $M$--estimation problem that will generate $K$ equations to be solved for the $K$ unknown parameters. Nonetheless, the recognition that the only moment conditions that contribute to enforcing the independence requirement are those whose imposition simultaneously reduces the objective function while providing multipliers that are coefficients in the stochastic representation of $Y$ suggests the futility of portmanteau approaches (e.g., those based on characteristic functions) to imposing independence. The dual formulation reveals that to specify the binding moment conditions is to specify an approximating data-generating process representation, which then can be extrapolated to provide estimates of objects of interest beyond the $n$ explicitly estimated values of $\varepsilon_{i}$ that characterize the sample and the definition of the mathematical program.
As is well understood in mathematical programming, dual solutions provide lower bounds on the values obtained by primal problems. In the generic form of the problems we have considered here there is no gap between the primal and dual values; hence in econometrics these problems are said to display point identification. We conjecture that the problems without point identification do have gaps between their dual and primal values, and that this characterization will enhance our understanding.