Multivariate Ordered Discrete Response Models with Lattice Structures

Multivariate Ordered Discrete Response Models with Lattice Structures

Introduction

Ordered response models are fundamental in empirical economics, used to analyze discrete choices with inherent ordering, such as risk aversion MalmendierNagel2011, political violence BesleyPersson2011, or educational attainment CameronHeckman1998. These models map latent continuous variables to discrete outcomes via thresholds. While univariate models are well-established (CunhaHeckmanNavarro2007, among many others), multivariate extensions can allow researchers to capture joint decisions across multiple dimensions.

Univariate ordered response models, like ordered probit and logit, were formalized by McKelveyZavoina1975 and AndresonPhilips1981. Multivariate extensions, introduced by AshfordSowden1970 for bivariate probit models, account for correlated decisions. Psychometrics and structural-equation literature adopted and extended the latent-variable viewpoint to multiple categorical indicators. In particular, Muthen1984 formalized a structural equation framework that allowed dichotomous and ordered categorical indicators to be treated as manifestations of underlying (multivariate normal) latent variables -- effectively a multivariate ordered model within the SEM tradition. SEM/psychometrics work (e.g., Olsson1979, among others) set out approaches for polychoric/probit models for multiple ordinal indicators. Kim1995 explicitly proposed and implemented a bivariate cumulative probit regression model for ordered categorical margins, with application and numerical estimation details. For detailed coverage of various types of univariate ordered response model we refer the reader to Agresti1990, Boes2006, STEWART2005, greene2010. greene2010 includes a review of recent applications of the bivariate ordered probit model. Applications of trivariate ordered probit models include buliung2007,genius2005,scott2001.

We focus on a particular class of multivariate ordered response models with a lattice structure, where decision makers narrowly bracket their choices, treating dimensions in isolation, in line with the behavioral economics framework of narrow bracketing ReadLoewensteinRabin1999. The lattice structure is characterized by functionally independent decision thresholds across dimensions, producing a grid-like latent space -- hence our terminology `lattice models.”\footnote{This terminology is our own and is not standard in the literature.} In practice, lattice models (often coupled with some parametric assumptions on the distribution of unobservables) have often served as the default and most straightforward extension of univariate ordered response models in applied work. kmnonlattice explicitly adopt the term “lattice models” to distinguish these restricted structures from more general multivariate formulations.

In this paper, we develop a formal and rigorous semiparametric framework for lattice ordered response models, where latent processes are specified as linear combinations of covariates and unobserved errors. We derive identification conditions for regression parameters, thresholds, and the joint distribution of the unobservables. The literature on univariate ordered models models provides several foundational insights that aid some identification results in multivariate settings as narrow bracketing allows us to isolate decision making across different dimensions. Namley, under full independence between unobservables and covariates, identification of index parameters and thresholds can rely on single-index methodologies just like in the univariate case. More general semiparametric approaches have allowed weaker conditions: lee1992 studied median independence following manski1975,manski1985,manski1988 work on maximum score, while lewbel2000 and chenkhan2003 allowed heteroskedastic unobservables with the latter focusing on multiplicative heteroskedasticity. In our analysis, we maintain full stochastic independence between unobservables and covariates to primarily focus on the identification and estimation of the joint cumulative distribution function (c.d.f.) which is a topic largely unexplored in the literature, even for lattice models.

Identification of the joint c.d.f. of unobservables in semiparametric models is a core theoretical contribution of this paper. Understanding this joint distribution is crucial for policy analysis. In lattice models, complementarity and substitutability in decision structures are not directly modeled. Thus, all dependence in observed decisions (conditional on covariates) is captured by the dependence among unobservables. This dependence structure is central to policy design involving joint outcomes such as household decisions on healthcare and education investments where the correlation between latent factors determines whether bundled interventions reinforce or crowd out each other. Semiparametric identification avoids restrictive parametric assumptions (e.g., joint normality) that can distort estimated policy effects if misspecified MalmendierNagel2011.

From an estimation perspective, we outline how existing semiparametric estimation methods can recover index parameters and thresholds in a sample, and we discuss how one could approach the estimation of the joint c.d.f. after those parameters are estimated at the $\sqrt{n}$-rate. We also describe how the approach of Coppejans2007 can be extended to jointly estimate all unknown components in one step.

For the parametric case, we focus on the multivariate normal specification, which conveniently captures varying degrees of dependence.\footnote{Alternative parametric specifications for the joint c.d.f. in bivariate ordered response models include Forcina2008 and Ferdous2010.} Because the lattice structure allows identification results for thresholds and indices to extend from univariate models, our attention centers on identifying the correlation parameters. We provide several sufficient conditions for identification in the bivariate case, including (i) configurations where one latent index is pinned at zero, (ii) variation in sign of index–threshold differences across subgroups, and (iii) the presence of exclusive covariates that shift one margin but not the other.

In short, this paper provides a rigorous foundation for lattice ordered response models, establishing semiparametric identification and outlining estimation strategies that make these models suitable for empirical applications where narrow bracketing is plausible such as consumer preference formation Train2009 and policy evaluation HeckmanVytlacil2007.

The remainder of the paper is structured as follows. Section (ref) introduces the general multivariate lattice model. Section (ref) develops the semiparametric specification, identification results and also discusses various approaches to estimation includiing those that utilize existing estimation techniques for univariate models, Section (ref) details the parametric model focusing on multivariate normal errors and identification of correlation coefficients. Section (ref) presents simulation evidence, and Section (ref) provides an empirical application estimating a joint ordered response model for health and happiness rankings. Section (ref) concludes. The Appendix collects proofs of the main theoretical results.

Model formulation

We model a single agent’s decisions across $D \geq 2$ dimensions, mapping a $D$-variate latent continuous metric $\left(Y^{* c_1}, \ldots, Y^{* c_D}\right)$ to a discrete metric $\left(Y^{c_1}, \ldots, Y^{c_D}\right)$. Discrete responses in dimension $d$ are $y_j^{(d)}, j=1, \ldots, M_d$, with ordering $y_1^{(d)} < \ldots < y_{M_d}^{(d)}$.

Thresholds $\alpha_{j_d}^{(d)}$ depend only on $j_d$, ensuring functionally independent decision rules across dimensions. The intersections of these threshodls across different dimensions form a lattice in $\mathbb{R}^D$. This reflects narrow bracketing ReadLoewensteinRabin1999, with intervals $\mathcal{I}_{j_d}^{(d)}$ partitioning $\mathbb{R}$ and rectangles $\times_{d=1}^D \mathcal{I}_{j_d}^{(d)}$ partitioning the latent space.

Semiparametric specification

The $d^{\text{th}}$ latent process is \[ Y^{* c_d} = x_d \beta_d + \varepsilon_d, \quad d=1, \ldots, D, \] where $x_d$ is a row vector of covariates, $\beta_d$ a column vector of parameters, and $\varepsilon_d$ an error term. Errors in $\left(\varepsilon_1, \ldots, \varepsilon_D\right)$ may be correlated, allowing latent processes $Y^{*cd}$ to be correlated conditional on observables.

Let $x = (x_1, \dots, x_D)$ and $\varepsilon = (\varepsilon_1, \dots, \varepsilon_D)'$ combine full vectors of covariates and unobservables, respectively. Denote the joint c.d.f. of $\varepsilon$ as $F$ and the marginal c.d.f. of $\varepsilon_d$ as $F_d$, $d=1,\dots,D$. The length of vector $x_d$ is $k_d$, $d=1,\dots,D$. Let $\mathcal{X}_d$ denote the support of $x_d$ and for each $d$, define \[ S^{(d;j)} = \{x_d \in \mathcal{X}_d \mid P(Y^{(d)} \leq y^{(d)}_j | x_d) \in (0,1)\}, \quad j=1, \ldots, M_d, \] and $S^{(d)} = \cup_{j=1}^{M_d} S^{(d;j)}$. Let $x_{d,m}$ denote the $m$th component of $x_d$ and $x_{d,-m}$ denote the subvector of $x_d$ excluding the $m$th component, with similar notations for $\beta$. $S^{(d)}_m$ denotes the projection of $S^{(d)}$ on $x_{d,m}$ with $S^{(d)}_{-m}$ being the projection on $x_{d,-m}$.

Identification

We derive identification conditions for $\beta_d$, thresholds $\alpha_{j_d}^{(d)}$ and the joint c.d.f. of unobservables under certainm assumptions. We start with Assumption (ref).

In univariate ordered response models, the assumption of independence between the unobservable and covariates is common, being used in klein2002, Coppejans2007, among many others.\footnote{Some papers (see e.g. chenkhan2003) on univariate ordered response allow for heteroskedasticity. In our framework, this would correspond to $\sigma_d(x_d,\theta_0) \varepsilon_d$ with independent $\varepsilon_d$. Some other papers further deviate from the setting of independence. lee1992 considers ordered response under the median independence assumption from manski1975, manski1985. In a recent paper, wangchen take a partial identification approach and consider a generalized maximum score estimator when regressors are interval measured. All of these settings are beyond the scope of this paper and provide interesting avenues for extensions of our work.} We formulate an analogue of a rank condition in the form of Assumption (ref).

Identification of threshold differences or gaps requires additional conditions to those assumed in Theorem (ref). This is given in Theorem (ref).

The new condition of Theorem (ref) would be guaranteed if for sets $S^{(d;j)}$ and $S^{(d;j+1)}$ the intersection of the sets of probabilities $\left\{P\left(Y^{c_d} \leq y^{(d)}_{j} \, | \, x_d \right): x_d \in S^{(d;j)} \right\}$ and $\left\{P\left(Y^{c_d} \leq y^{(d)}_{j+1} \, | \, x_d \right): x_d \in S^{(d;j+1)} \right\}$ contains an interval $(\underline{p}_j,\overline{p}_{j})$. Large support conditions would, e.g, ensure that this interval is $(0,1)$.

Figure (ref), which shows a bivariate lattice model, presents an intuitive summary of the identification strategy in the models with lattice structures. We consider each dimension individually and, within that dimension, express probabilities of discrete values up to certain points in terms of the marginal c.d.f. of the unobservable in that dimension and the index in that dimension. Theorem (ref) is based on consideting just one shaded area for many different $x_d$ -- either the one the left panel or the one on the right panel in Figure (ref). Theorem (ref) requires the computation of both shaded regions for many different $x_d$.

The result of Theorem (ref) immediately implies conditions for identification of marginal distributions of $\varepsilon_d$, $d=1, \ldots, D$.

Condition ((ref)) ensures that any point in the support of $\varepsilon_d$ corresponds to the underlying $\alpha^{(d)}_j-x_d\beta_d$ for some $j$ and $x_d$. Condition (i) explicitly normalizes one threshold (the identification of values of the other thresholds then immediately follows from Theorem (ref)), whereas condition (ii) enforces a normalization of one threshold in an indirect way.

The result of Theorem (ref) does not guarantee identification of the joint distribution of unobservables, even if the conditions of this corollary hold for every $d=1,\ldots,D$. The reason is two-fold. First, Assumption (ref) does not give any information about how the vector $\varepsilon$ relates to $x_h$, $h \neq d$. Under a full stochastic independence of the vector $\varepsilon$ from the whole vector $x$, the identification process easier as $P(\varepsilon_1 \leq e_1, \ldots, \varepsilon_D \leq e_D|x)$ does not depend on $x$ and we only need to identify one $D$-variate c.d.f. $F(e_1,\ldots, e_D)=P(\varepsilon_1 \leq e_1, \ldots, \varepsilon_D \leq e_D)$. The main channel in which we can proceed with identification of $F$ is considering observed probabilities $$P\left(Y^{(1)} \leq y^{(1)}_{j_1}, \ldots, Y^{(D)} \leq y^{(D)}_{j_D}|x\right) = F(\alpha^{(1)}_{j_1} - x_1\beta_1, \ldots, \alpha^{(D)}_{j_D} - x_D\beta_D)$$ but then the question becomes of whether the data provides enough joint variation in indices $(x_1\beta_1, \ldots, x_D\beta_D)$ to identify $F$ on the whole support $\mathcal{E}$ of $\varepsilon$. The issue is that some (potentially each) $x_{d}$ could share all its covariates with another process. In this case $(\alpha^{(1)}_{j_1} - x_1\beta_1, \ldots, \alpha^{(D)}_{j_D} - x_D\beta_D)'$ could take values only in a proper subset of $\mathcal{E}$ and could vary only in certain directions as we vary the values of covariates. Since at this identification stage $(\alpha^{(1)}_{j_1} - x_1\beta_1, \ldots, \alpha^{(D)}_{j_D} - x_D\beta_D)'$ is observed, one could try and assess whether this vector covers the whole support $\mathcal{E}$. What we do is present conditions under which this is guaranteed. The illustration of our idea is given in Figure (ref) for $D=2$. In one dimension (e.g. for $\varepsilon_1$) we ensure that $\alpha^{(1)}_{j_1} - x_1\beta_1$ can cover the whole marginal support of $\varepsilon_1$ (can be checked using conditions of Theorem (ref)). In the other dimension (e.g. for $\varepsilon_2$) we can require an exclusive covariate with non-zero coefficient -- without a loss of generality $x_{2,1}$ -- that can provide enough own variation in $\alpha^{(2)}_{j_2} - x_2\beta_2$ while keeping $x_1\beta_1$ fixed. In Figure (ref) this variation is shown using vertical arrows. Once $x_1\beta_1$ is fixed, this variation can be checked to cover both lower and upper boundaries of $\mathcal{E}$ (either finite or infinite) by checking whether $\sup_{j_2} \sup_{x_{2,1}} F(\alpha^{(1)}_{j_1} - x_1\beta_1, \alpha^{(2)}_{j_2} - x_2\beta_2)$ coincides with $F(\alpha^{(1)}_{j_1} - x_1\beta_1)$ (upper boundary) and whether $\inf_{j_2} \inf_{x_{2,1}} F(\alpha^{(1)}_{j_1} - x_1\beta_1, \alpha^{(2)}_{j_2} - x_2\beta_2)$ is 0 (lower boundary). For general $D$, this identification strategy can be translated into the requirements on exclusive covariates in $D-1$ processes.

Conditions ((ref)) and ((ref)) guarantee that $(\alpha^{(1)}_{j_1} - x_1\beta_1, \ldots, \alpha^{(D)}_{j_D} - x_D\beta_D)'$ for some $j_1,\ldots,j_D$ when taken in any direction $\lambda$ in $\mathbb{R}^D$ can reach the boundary of $\mathcal{E}$ in both positive and negative directions of $\lambda$.

To illustrate the progressive restrictiveness of the identification conditions outlined in Theorems (ref) through (ref), we construct four nested data-generating processes (DGPs) for a bivariate ($D=2$) lattice model, each building sequentially on its predecessor. Each latent process contains a two-dimensional covariate vector associated with $\beta_1 = \beta_2 = (1, 0.5)'$. In each dimension, there are three ordered responses and the threshold differences are 2, Suppose the vector of unobservables is independent of covarioates and has a joint normal distribution.

In DGP 1, covariates are defined as $x_1 = x_2=(x_{common1}, x_{common2})$, where $x_{common1}, x_{common2} \sim \text{Uniform}[-0.5, 0.5]$ and $x_{common1}, x_{common2}$ are not perfectly linearly related. This DGP provides limited support for $x_1\beta_1$ and $x_2\beta_2$ (it is within $[-0.75,0.75]$). Theorem 1 is satisfied, which ensures identification of $\beta_{1}$, $\beta_{2}$ up to scale, but fails to meet the conditions of Theorem (ref) as it lacks overlaps in choice probabilities for threshold differences, Indeed, $P(Y^{c_d} \leq y^{(d)}_1|x_d) \in [\alpha^{(d)}_1-0.75,\alpha^{(d)}_1+0.75]$ whereas $P(Y^{c_d} \leq y^{(d)}_2|x_d) \in [2+\alpha^{(d)}_1-0.75,2+\alpha^{(d)}_1+0.75]=[\alpha^{(d)}_1 +1.25,\alpha^{(d)}_1+2.75]$ with $[\alpha^{(d)}_1-0.75,\alpha^{(d)}_1+0.75]$ and $[\alpha^{(d)}_1 +1.25,\alpha^{(d)}_1+2.75]$ obviously not overlapping. The narrow range of the indices precludes the probability matching required by Theorem (ref).

DGP 2 extends the first by widening the support of covariates: $x_{common1} \sim \text{Uniform}[-2, 2]$, $x_{common2} \sim \text{Uniform}[-0.5, 0.5]$ enabling overlaps in conditional probabilities (e.g., $P(Y^{c_d} \leq y^{(d)}_1|x_d) \in [\alpha^{(d)}_1-2.25,\alpha^{(d)}_1+2.25]$ and $P(Y^{c_d} \leq y^{(d)}_2|x_d) \in [\alpha^{(d)}_1-0.25,\alpha^{(d)}_1+4.25]$). This satisfies the conditions up to Theorem 2 but falls short of Theorem (ref), as the support, while sufficient for probability matching, does not fully cover the interval (0, 1). The added restrictiveness stems from the need for broader support to align probabilities, yet the coverage remains incomplete.

DGP 3 further extends the second by setting $x_{common1} \sim \text{Laplace}$, $x_{common2} \sim \text{Uniform}[-0.5, 0.5]$ to ensure full probability coverage over (0, 1), and incorporates a normalization $F_d(0) = 0.5$. This setup satisfies the conditions up to Theorem (ref) but fails Theorem (ref). as the absence of exclusive covariates prevents independent shifting of dimensions to capture joint dependence.

DGP 4 builds on the third by defining $x_{excl1} \sim Laplace$, $x_{excl2} \sim Laplace$, $x_{common2} \sim \text{Uniform}[-0.5, 0.5]$ (the suppose of the distribution of $(x_{excl1},x_{excl2},x_{common2})$ has an interior in $\mathbf{R}^3$, This allows independent shifting of dimensions 1 abd 2, satisfying the requirement of Theorem (ref) (note this theorem only requires independent shifting of one dimension but for simplicity we allow that in both dimensions). All the parameters including the joint c.d.f. can then be identified fully.

Estimation

In what follows, we briefly outline some possibilities for estimating parameters in semiparametric models. A theme of this section is to outline existing univariate ordered response estimation methods that generalize to lattice models.

\paragraph*{Two-step approach}

The idea of this method is to (i) first use existing estimation approaches for semiparametric univariate ordered response models to estimate index and threshold parameters at a suitable rate (albeit suboptimally as the dependence of the latent processes is ignored), and (ii) second, construct estimates of the joint c.d.f. using some well known statistical methods.

We start by discussing which estimation approaches in the literature can be utilized in the first step.

lewbel2000 develops a semiparametric estimator for qualitative response models (binary, ordered, multinomial) allowing for unknown heteroscedasticity in the latent errors with respect to regressors, or instrumental variables for endogeneity. The method relies on a “special regressor” $v$ that is conditionally independent of the error $\varepsilon$ given other regressors $x$ (i.e., $F_{\varepsilon|v,x}(\varepsilon \mid v, x) = F_{\varepsilon|x}(\varepsilon \mid x)$), with large support. The estimator resembles OLS or 2SLS on a transformed response $y^* = [y - I(v < 0)] / f(v \mid x)$, where $f$ is the conditional density of $v$ given $x$, yielding for ordered response models $\sqrt{n}$-consistent and asymptotically normal estimates for coefficients $\beta$ and thresholds (for ordered models).

To generalize lewbel2000 to lattice models, we need to have $x_d=(v_d,w_d)$ with a continuous special regression $v_d$ with large support per dimension $d$ -- this would effectively extend our Assumption ... (and accommodating heteroscedasticty by allowing $\text{Var}(\varepsilon_d \mid x_d)$ to be arbitrary). The estimator would proceed marginally per dimension using lewbel2000 ordered method to recover $\beta_d$ and thresholds $\alpha^{(d)}_j$. lewbel2000 consider a univariate ordered response, hence the question of joint c.d.f. does not arise (note, however, that for multinomial choice the estimation of joint c.d.f. of unobservables is relevant but lewbel2000 does not address it).

klein2002 approach analyzes the univariate model, estimates the index parameter in the first stage using kernel density estimates of the conditional probability of choosing below a certain level. In the second stage, the approach estimates threshold parameters using shift restrictions. We can extend this approach to multivariate lattice models because the functional independence of thresholds across dimensions allows us to apply stages 1 and 2 marginally for each $d=1,\dots,D$, using univariate techniques and our Assumption (ref) which mirrors a cre assumption of independence in klein2002 and leads to $P(Y_c^d \leq y^{(d)}_j | x_d) = F_d(\alpha^{(d)}_j - x_d \beta_d)$. The estimators or index and threshold parameters obtained from this stage are $\sqrt{n}$-consistent and asymptotically normal.

chenkhan2003 derives rates of convergence for estimating index parameters in heteroskedastic discrete response models, assuming multiplicative heteroskedasticity $\varepsilon_i = \sigma(x_i) \cdot u_i$, where $u_i$ is homoskedastic and independent of $x_i$. For ordered response models with at least three categories, $\sqrt{n}$-consistent estimators are possible. To generalize chenkhan2003 to lattice models, we can consider each dimesn ion $d$ separately and consider at least three responses in that dimension. at the same time, we can generalize it to multiplicative heteroskedasticity per dimension: $\varepsilon_d = \sigma_d(x_d) \cdot u_d$, where $u_d$ is homoskedastic and independent of $x_d$. The chenkhan2003 estimator for index parameters and thresholds proceeds marginally per dimension. Marginal stages inherit rates from chenkhan2003: $\sqrt{n}$-consistent $\hat{\beta}_d$, $\hat{\alpha}^{(d)}_j$ for $M_d \geq 3$.

Liu2024 proposes two simple semiparametric estimators for univariate ordered response models with an unknown error distribution $F_0$, achieving $\sqrt{n}$-consistent and asymptotically normal estimators of the index parameters and thresholds. The first method (binary choice-based) constructs nonparametric maximum likelihood estimates (NPMLE) of $F_0$ from recast binary data, then uses moment conditions index and threshold parameters. The second method (full ordered data) extends this by incorporating all outcomes via a weighted NPMLE. Both enforce monotonicity of $F_0$ and use bootstrap for inference. In lattice models, one can apply Liu2024 methods marginally per dimension to estimate $\beta_d$ and thresholds $\alpha^{(d)}_j$ (up to scale/location). All these estimators will be $\sqrt{n}$-consistent and asymptotically normal.

Thus, all these approaches are suitable when one's goal is to estimate index and thresholds parameters. Given these estimates, one can now proceed with the estimation of the joint c.d.f. $F$ in the second stage (this, of course, is not addressed in the papers mentioned above due to the univariate nature of the problem there). Let us now discuss some specific approaches that can be used to obtain $\widehat{F}$.

One possible approach is the grid inversion method that discretizes the error space and solves a constrained optimization problem. It is a direct, computationally intensive non-parametric method. Let us outline it for $D=2$. Its idea is based on the fact that given $x_i=(x_{i1},x_{i2})$ and $(Y^{(c_1)}_{j_1}=y^{(1)}_{j_1}, Y^{(c_2)}_{j_2}=y^{(2)}_{j_2})$, the latent pair $(\varepsilon_{1i}, \varepsilon_{2i})$ lies in the rectangle $R_i= \times_{d=1}^2 \left( \alpha^{(d)}_{j_d-1} - x_{id}\beta_d, \, \alpha^{(d)}_{j_d} - x_{id}\beta_d \right]$ and, hence, due to independence of errors from covariates,

In the sample each observation $i$ implies a rectangular interval $ \widehat{R}_i=\widehat{\mathbf{\varepsilon}}_i \in \times_{d=1}^2 \left( \hat{\alpha}^{(d)}_{j_{d}(i)-1} - x_{di}\widehat{\beta}_d, \hat{\alpha}^{(d)}_{j_{d} (i)} - x_{di}\widehat{\beta}_d \right] $ for the residual $\widehat{\mathbf{\varepsilon}}_i$, where $j_{d}(i)$ is the observed category in dimension $d$ for $i$ (with $-\infty, +\infty$ boundaries). Let $\mathcal{G} = \{ (e_{1,k}, e_{2,\ell}) : k=1,\dots,K_1,\ \ell=1,\dots,K_2 \}$ be the set of unique lower/upper bounds from all such implied sample rectangles. Let $ \phi = \big( F(e_{1,k}, e_{2,\ell}) \big)_{k,\ell} \in \mathbb{R}^{K_1 K_2} $ collect the unknown c.d.f. values on this grid. We want to find the probability mass assigned to each grid point such that the implied probabilities for each cell match the empirical probabilities in the data as closely as possible. To do this, for each distinct covariate pattern $x_g$ (group), define the empirical cell probabilities \[ \widehat{\pi}_{j_1 j_2}(x_g) \equiv \widehat{P}(Y^{(c_1)}_{i}=y^{(1)}_{j_1(i)}, Y^{(c_2)}_{i}=y^{(2)}_{j_2(i)} \mid X=x_i) = \frac{ \sum_{i : x_i = x_g} \mathbf{1}\{ Y^{(1)}_{i} = y^{(1)}_{j_1}, Y^{(2)}_{i} = y^{(2)}_{j_2} \}}{ \sum_{i }: 1(x_i = x_g) }. \]

Then for each $(j_1,j_2,g)$, \[ \widehat{\pi}_{j_1 j_2}(x_g) = \sum_{k,\ell} A_{j_1 j_2,g}(k,\ell)\, \phi_{k\ell} + u_{j_1 j_2,g}, \] where $A_{j_1 j_2,g}(k,\ell) \in \{-1,0,1\}$ encodes which c.d.f. corner terms enter each rectangle probability using ((ref)). Stacking over all $(j_1,j_2,g)$ yields $A \phi = \widehat{\pi} + u, $ where $\widehat{\pi}$ collects all empirical cell probabilities. We can estimate $\phi$ by solving $\widehat{\phi} = \arg\min_{\phi \in \mathbf{\Phi}} \|A \phi - \widehat{\pi}\|^2$, where the feasible set $\mathbf{\Phi}$ enforces the defining properties of a c.d.f.: \[ \mathbf{\Phi} = \left\{ \phi : 0 \le \phi_{k\ell} \le 1, \, \phi_{k\ell} \text{ nondecreasing in } k \text{ and in } \ell \right\}. \] Optionally we can include a smoothness penalty and optimize $\min_{\phi \in \mathbf{\Phi}} \|A \phi - \widehat{\pi}\|^2 + \lambda \|D \phi\|^2$, where $D$ is a finite-difference matrix. The estimator provides \[ \widehat{F}(e_{1,k}, e_{2,\ell}) = \widehat{\phi}_{k\ell}, \] which can be extended to a continuous surface by bilinear interpolation.

There are some variations of this method. E.g., instead of the grid determined by the implied rectangular regions, one can consider a completely exogenous sample-free grid.

Another possible approach is the kernel smoothing approach. Just like the inversion grid method it uses the fact that $\mathbf{\varepsilon}_{i}\in R_i$ given $x_i=(x_{i1},x_{i2})$ and $(Y^{(c_1)}_{j_1}=y^{(1)}_{j_1}, Y^{(c_2)}_{j_2}=y^{(2)}_{j_2})$, and with with $R_i$ defined in the same way as in the grid inversion method. In the sample each observation $i$ implies $ \widehat{\mathbf{\varepsilon}}_i \in \widehat{R}_i $. We can implement a simulated kernel density estimator, where for each observation $i$ we draw $S$ random samples $(\widetilde{\varepsilon}^{(s)}_1,\widetilde{\varepsilon}^{(s)}_2)$ uniformly from its rectangle $\widehat{R}_i$. We then pool all these $N \times S$ simulated points together. We then perform a standard bivariate kernel density estimation on this large pooled sample. The resulting density is an estimate of $f(\varepsilon_1,\varepsilon_2)$. We can then integrate this estimated density numerically to get the estimated c.d.f. .

Other possible approaches include nonparametric sieve estimator subject to suitable choice of base (for monotonicity-preserving properties) and nonparametric maximum likelihood estimator. We have implemented the grid inversion and the simulated kernel density estimator in simulations biut not the other approaches.

\paragraph*{One-step approach}

If one is interested in estimating the joint c.d.f of unobservable $\varepsilon$ (for purposes of analysing policy intervention or other counterfactuals), then one could extend Coppejans2007 originally developed for univariate ordered response models under independence of the error and covariates. In what follows, we extend it to multivariate ordered response models, describing the bivariate case for illustrational simplicity. Suppose we have a random sample $\left\{(y^{(1)(i)}, y^{(2)(i)},x_1^{(i)},x_2^{(i)}) \right\}_{i=1}^N$. The idea is to maximize the log-likelihood function

for joint c.d.f. of unobservables $F$. Coppejans2007 uses a quadratic B-spline to estimate the c.d.f of unobservables. The multivariate analogy is tensor-product B-splines. For instance, in the bivariate case the tensor-product basis consists of $S_1\cdot S_2$ products of polynomials $\mathcal{R}$ in the form

here calculated for specific values of $e_1$ and $e_2$, with $q_d$ denoting the degree of B-spline in dimension $d=1,2$. A general tensor-product B-spline, which approximates $F(e_1,e_2)$, is a linear combination of these base tensor-product polynomials with coefficients $\{h_{s_{1}s_{2}}\}$, $s_{d}=1,\ldots ,S_{d}$, $d=1,2$:

The linear constraints

guarantee monotonicity of the tensor-product B-spline in each dimension. Additionally, the linear constraints $$ 0 \leq h_{s_{1},s_{2}} \leq 1, \quad \forall \, s_1, s_2 $$ guarantee natural c.d.f. bounds of 0 and 1.\footnote{For more details on shape constraints in tensor-product B-splines, see bk2022.} Linear equality constraints on $h_{s_1s_2}$ can impose normalization restrictions on $F_d$: ...

Parametric specification

In practice a researcher may choose a parametric family to model the distribution of unobservables conditional on covariates. On one hand, choosing a parametric family may allow researcher to explicitly model the distribution of observables as that depending on $x$ and be able to identify all the primitives given the assumed (potentially complicated) dependence structure. The exact identification strategy and assumption behind it will depend on the assumed structure. On the other hand, a researcher may still opt for independent errors and covariates and rely on less stringent data requirements for identification that those given in Section (ref) and a simpler estimation approach.

We illustrate the latter case focusing on the lattice ordered probit (Gaussian errors) case.

Identification

As expected, due to our ability to view decisions rules across different dimensions index and threshold parameters can be identified using the same rank condition commonly employed in univariate ordered probit models. This is formally presented in Theorem (ref) below. Its proof is well known and we replicate it in the Appendix purely for completeness.

Identification of correlation coefficients $\rho_{d_1,d_2}$ in the multivariate lattice setting does not follow from any readily available results in the literature. We can carry out this identification in the pairwise manner under supplementary variation/exclusion conditions. They are collected in Theorem (ref) below

Condition (a) requires a covariate configuration where the latent index in one dimension ($d_1$) is exactly at some threshold. It creates a “pivot” where the error $\varepsilon_{d_1}$ is symmetrically distributed around zero, making joint probabilities with dimension 2 purely a function of $\rho_{d_1,d_2}$'s influence on $\varepsilon_{d_2}$. Condition (b) requires sign-flipping covariates. Namely, it assumes covariate variation creating “same-sign” indices in one dimension (both above or both below the median threshold) but “sign-flipping” in the other. There are other ways to formulate related sufficient conditions in this spirit but we have opted to present this one. Condition (c) is an IV-style exclusion: a covariate (or subvector) affects dimension $d_1$'s outcome (via assocated nonzero $\beta_{d_1}$'s) but not dimension $d_2$'s directly (exclusion from $x_{d_2}$). By having a variable that affects only one outcome, we can trace out how joint probabilities shift when one margin's latent index moves while the other stays fixed. This variation rotates the joint probability surface (enough to do it once), letting us solve for $\rho_{d_1,d_2}$.

Estimation in the parametric model

Estimation in the parametric model is standard via maximum likelihood. The log-likelihood function is equal to that in equations ((ref)) and ((ref)) with a specified cumulative distribution function $F$. We use bivariate normal as the natural example of $F$, as in Assumption (ref), so that $\boldsymbol{\varepsilon}=(\varepsilon_1,\varepsilon_2)'$ is jointly normal with mean \((0,0)'\), unit variances, and correlation \(\rho\). Given a random sample $\left\{(y^{(1)(i)}, y^{(2)(i)},x_1^{(i)},x_2^{(i)}) \right\}_{i=1}^N$ and collecting $\beta_1,\beta_2,\rho$ and all the thresholds in $\alpha$ in one parameter vector $\theta$, we can construct the log-likelihood function

$$\text{ with } \quad \ell^{(i)}_{j_{1},j_{2}} = \sum_{t_1=0}^{1} \sum_{t_2=0}^{1} (-1)^{t_1 + t_2} \Phi_{2}\left( \alpha_{j_1 - t_1, j_2}^{(1)} - x_1^{(i)}\beta_1, \alpha_{j_1, j_2 - t_2}^{(2)} - x_2^{(i)}\beta_2; \rho \right)$$ where $\Phi(\cdot,\cdot;\rho)$ denotes the standard bivariate normal c.d.f. with correlation parameter $\rho$.

The maximum likelihood estimator (MLE) $\hat{\theta}$ solves the optimization problem $\max_{\theta} \mathcal{L}(\theta)$.\footnote{One can impose inequality constraints on $\alpha$ and $\rho$ and maximize a constrained likelihood, or, more straightforwardly, re-parameterize the likelihood to estimate ($\alpha^{(j)}_{1},\sqrt{\alpha^{(j)}_{2}-\alpha^{(j)}_{1}},\dots,\sqrt{\alpha_{M_{d}}^{(j)}-\alpha_{M_{d}-1}^{(d)}})$ and $\tanh{^-1}(\rho)$ so that constraints are enforced automatically.} Under the typical MLE regularity conditions (e.g. \citep*{newey1994}), we have $\sqrt{N}(\hat{\theta}-\theta_{0}) \overset{d}{\longrightarrow} \mathcal{N}(0,V)$, $V = \mathbb{E}\left[\frac{\partial \log(\ell^{(i)}(\theta_0))}{\partial \theta}\frac{\partial \log(\ell^{(i)}(\theta_0))}{\partial \theta'}\right]$. The natural plug-in sample-analogue estimator of $V$ provides a consistent estimator for the variance-covariance matrix.

Simulations

We consider a bivariate ordered response model with

Semiparametric model

Here we focus on two-step approaches. We take the index and threshold parameters to be known (in reality, they would have been estimated consistently at $\sqrt{n}$ rate) and just focus on the estimation of the joint c.d.f. given these parameters. This allows us to compare the performance of different second-stage approaches in their pure form without first-stage inference.\footnote{Moreover, first-stage estimation approaches come from already existing literature.} We take both $x_{1i}$ and $x_{2i}$ to be univariate with respective $\beta_1 = 0.8$, $\beta_2 = -0.5$. We adopt $3 \times 3$ categorical outcomes with the thresholds determining the decision structure given by $\alpha^{(1)}_0=-\infty$, $\alpha^{(1)}_1=-1$, $\alpha^{(1)}_2=1$, $\alpha^{(1)}_3=+\infty$ for dimension 1 and $\alpha^{(2)}_0=-\infty$, $\alpha^{(2)}_1=-0.8$, $\alpha^{(2)}_2=0.8$, $\alpha^{(2)}_3=+\infty$. for dimension 2. We take ${x}_1 \sim N(0,1)$, ${x}_2 \sim 0.5N(0,1) + 0.3{x}_1$, and $\boldsymbol{\varepsilon}$ is bivariate normal with mean zero, unit variances and the correlation coefficient 0.6. We draw $S=10$ points from the rectanbgle associated with observation $i$.

We compare the performance of our estimators on an $80 \times 80$ evaluation grid over $[-2.5, 2.5]^2$: We use $G=6,400$ to denote the number of points in the evaluation grid and $g$ to denote a particular point on this grid. As criteria we use Root Mean Square Error $\sqrt{\frac{1}{G}\sum_{g=1}^n (\hat{F}(g) - F(g))^2}$ (RMSE), Kolmogorov-Smirnov (KS) distance $\max_g |\hat{F}(g) - F(g)|$, Cramér-von Mises (CvM) distance $\frac{1}{G}\sum_{g=1}^G (\hat{F}(g) - F(g))^2$ and correlation $\text{corr}(\hat{{F}}, {F})$.

Table (ref) presents simulation results comparing both approaches on the average of the four metrics in 400 simulations.

On the basis of these results, the kernel smoothing method dominates the grid inversion method on the four metrics. We do not pursue further with regard to how well these methods do with regard to various regions (central vs tail ones) or which method performs better with regard to some specific distributional characteristics such as entropy, or tail mass but one could of course pursue this type of simulation analysis as well. It may very well be the case that the grid inversion method may perform better in other criteria as it explicitly incorporates the ordinal response structure through the design matrix $A$ and its associated constrained least squares formulation provides a globally optimal solution for the discrete approximation. Moreover. this

Parametric model

We now examine Monte Carlo simulations for the parametric case with normal errors. For the purposes of the simulations, we rewrite Equations ((ref)) and ((ref)) as $$Y_{i}^{*c_{1}} = x_{i}\beta_{1} + w_{1i}\gamma_{1} + \varepsilon_{1i}, \qquad Y_{i}^{*c_{2}} = x_{i}\beta_{2} + w_{2i}\gamma_{2} + \varepsilon_{2i}$$ to distinguish exclusive ($w$) and non-exclusive ($x$) covariates. We explore a first scenario with no exclusive covariates ($\gamma_{1}=\gamma_{2}=0$), a second scenario with an exclusive covariate in one latent process, and a third scenario with exclusive covariates in both latent processes. Each simulation design uses 400 independent random samples of size 1,000. A summary of the following results is that in all models, which vary in their number of discrete values $M$, type of regressors (discrete or continuous) and exclusivity of regressors, all parameters are estimated with essentially no bias; threshold and index parameters are estimated more precisely than the correlation parameter.

Parametric Design 1: 2$\times$2 structure, no excluded regressors

We investigate parametric estimation without exclusive covariates by setting $\gamma_1=\gamma_2=0$, removing $w_1$ and $w_2$. We set $\beta_1=3$, $\beta_2=2.5$, $\rho=0.33$, and use a $2\times2$ non-lattice structure with thresholds $\alpha^{(1)}_{1}=1$ and $\alpha_{1}^{(2)}=1.25$ (see Figure (ref)). The common regressor $x$ follows a uniform $[-4,4]$ distribution.

Table (ref) Panel 1 reports mean and standard deviation of parameter estimates. The method estimates all parameters with minimal bias. Estimates of $\rho$ are less precise due to the absence of excluded regressors.

Parametric Design 2: 4$\times$3 with one excluded covariate

In the second simulation design, we extend the number of discrete values $M_{d}$ in both dimensions. The discrete dependent variable $Y^{c_1}$ can take four values and $Y^{c_2}$ can take three values. This generates a 4$\times$3 structure, illustrated in Figure (ref). The common covariate $x$ follows a uniform $[-2,2]$ distribution (alternatively we could have taken it to be discrete). The covariate $w_{1}$ is a discrete random variable taking values -2.5, -1.5, -0.5 and 0.5 with equal probability. We set $\gamma_2=0$ thus effectively removing $w_{2}$ in the second equation. The parameter values are $\beta_{1} = 2, \gamma_{1} = -3, \beta_{2} = 3$ and $\rho = 0.25$.

Table (ref) Panel 2 lists the across-simulation means and standard deviations of the index parameters, thresholds, and the correlation coefficient. The bivariate ordered probit method estimates all parameters with no bias. The correlation parameter remains the least precise estimate across parameters.

Parametric Design 3: 6$\times$2

We consider a design that creates a 6$\times$2 structure on the latent variable space. Figure (ref) illustrates the threshold structure and the values of thresholds. In this design, the common regressor $x$ is drawn from uniform $[-2,2]$ and both latent equations have excluded regressors $w_{1},w_{2} \overset{iid}{\sim} t_{7}$. We also include an additional regressor $z_{2}$ in equation 2, drawn from a logistic (3,2) distribution. The parameter corresponding to $z_{2}$ is denoted $\delta_{2}$, so that the latent equations read

The parameter values are $\beta_1 = 1.75$, $\beta_2 = 2.5$, $\gamma_1 = -2.75$, $\gamma_2=-4$, and $\delta_{2} = 2$. Table (ref) presents the results for the index parameters, thresholds, and the correlation coefficient. Generally, all parameters are estimated with low bias, though slightly more bias in the index parameters in dimension two than one.

Application

Now we study the main factors driving self-reported health and happiness as well co-movement in unobservables driving them. To do this, we pool data on the United States of America and Canada from six waves of the World Values Survey inglehart2014. The results are also fully robust to the use of a different dataset: the National Health and Nutrition Examination Survey (NHANES). A full description of these two datasets and variable construction is provided in the Appendix.

The specification considered follows the standard setup described in the paper. Namely, for latent physical health ($p$) and sadness ($m$) variables $Y^{*}_{p}$ and $Y^{*}_{m}$ respectively, we have

with common row of covariates $x$ and exclusive covariates $w_{p}$ and $w_{m}$ in those two processes. The discrete dependent variable for health we use takes three values: 0, 1, and 2. The value 0 represents a self-reported “State of health” as “fair”, “poor”, or “very poor”. The value 1 represents a report of “good”, and the value 2 a report of “very good”. The dependent variable for happiness again takes values 0, 1, and 2. In this case, a value of 0 represents a self-reported “Feeling of happiness” as “very happy”. A value of 1 represents a reporting of “quite happy”, and a value of 2 a reporting of either “not very happy” or “not at all happy”. It is coded so that higher values reflect lower self-reported happiness.

The common set of regressors $x$ includes the variabless: male, white, a college education dummy, age, regional dummies, and 5 dummies for income brackets. The 5 income brackets are: (1) less than \$20,00; (2) between \$20,000 and \$35,000; (3) between \$35,000 and \$50,000; (4) between \$50,000 and \$ 75,000; and (5) greater than \$100,000. We include no excluded health regressors, so that $w_{p} = 0$, but include dummies for employment status and living with a partner as excluded happiness regressors, $w_{m}$.

The coefficients from the lattice bivariate ordered probit regression are provided in Table (ref). The signs of the regression coefficients are as expected. Positive partial effects on the probability of a better health (above a certain level) are given by varaiables that include dummies for white ethnicity and college education, and higher income brackets. The variables that have a positive impact on the probability of a higher happiness (or, equivalently, lower sadness) include living with a partner and higher income brackets. The effect of employment status on the probability of a higher happiness appears negative but is not statistically significant.

The estimated value of the correlation between the two model errors is negative, which is consistent with our expectations that shocks increasing health would tend to decrease sadness and shocks decreasing health would tend to increase sadness. The thresholds produced are shown in Figures (ref).

Detailed estimation results for index parameters are given in Table

Conclusion

We formulate lattice ordered response models for narrow bracketing, identifying parameters, thresholds, and the joint c.d.f. in a semiparametric framework. In the bivariate probit case, we separately identify $\beta_d$ and $\alpha_{j_d}^{(d)}$ using marginal probabilities, and $\rho$ using joint probabilities with an exclusive covariate. The lattice structure simplifies estimation, suitable for empirical applications. Future work could develop estimation methods.