# Development of Integrated Choice and Latent Variable (ICLV) Models Using Matrix-Based Analytic Approximation and Automatic Differentiation Methods on TensorFlow Platform

Authors: Jie Ma; Xin Ye; Kun Huang
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

This study further explores the multinomial probit-based integrated choice and latent variable (ICLV) models. The LDLT matrix-based analytic approximation methods, including Mendell and Elston (ME) method, bivariate ME (BME) method, and two-variate bivariate screening (TVBS) method, were adapted to calculate the multivariate cumulative normal distribution (MVNCD) function in the ICLV model because of the better performances in accuracy and computational time. Integrated with the composite marginal likelihood (CML) estimation approach, the ICLV model based on high-dimensional integration can be estimated accurately within a reasonable time. In this study, some three-alternative and four-alternative ICLV models are simulated to examine their abilities to recover model parameters. It is found that the parameter estimates and standard error estimates are acceptable for both models and the computational time is expected to decrease using tensor data structures on the TensorFlow platform. For the four-alternative ICLV models, the TVBS method has the highest level of accuracy. The BME method is also a good alternative to TVBS if computational time is of great concern. The application of the automatic differentiation (AD) technique in the model can free researchers from coding analytical gradients of log-likelihood functions and thereby greatly reduce the workload of researchers.


---

## Body

## 1. Introduction

In the discrete choice model (DCM), the choice depends not only on the transport characteristics and the observable characteristics of decision-makers but also on the unobservable characteristics, namely latent variables. Among them, psychological and attitudinal factors as parts of unobservable latent variables have been proved to be significant when being integrated into DCMs [1].One alternative to incorporate latent variables into discrete choice models is the integrated choice and latent variable (ICLV) model [2, 3], which can be viewed as a variation of the traditional structural equation model (SEM) with the advantage to involve the underlying choice process based on psychological concepts. The ICLV model is a popular research topic in the field of transport choice studies since the early 2000s, and the application of this model has been increasing in recent years [3–7]. The earlier ICLV models are generally developed on the basis of the mixed logit model, of which the error terms follow an independent and identically Gumbel distribution. These logit kernel-based ICLV models ignore the correlation among the utility of alternatives and the correlation among the error terms of latent variables. Meanwhile, the models have several deficiencies in the aspects of the estimation method (i.e., maximum simulated likelihood (MSL) approach [8, 9]), the convergence of model estimation with multiple latent variables [7], and computational time of multidimensional integration. Regarding the limitations of the above models, Bhat and Dubey [10] firstly proposed an multinomial probit (MNP) kernel-based ICLV formulation method, which can accommodate a large number of latent variables, and performances in convergence and computational efficiency have been improved. In particular, the maximum approximate composite marginal likelihood (MACML) inference approach is integrated to estimate the model parameters, which further simplifies the estimation process. As the number of latent variables is independent of the integral dimensions of the log-likelihood function and only relates to the number of alternatives of choice models, the computational time of this method is relatively short. Besides, the ordered and continuous response indicators for the latent variables can be involved in the model easily. Up to now, the new MNP-based ICLV model has not been extensively applied [11, 12] and remains a broad prospect for practice.At present, it is critical to compute the multivariate cumulative normal distribution (MVNCD) function of the proposed ICLV model, which can be further improved. In some earlier studies, the simulation techniques are used for evaluating the MVNCD functions, of which the most common method is using a Geweke–Hajivassiliou–Keane (GHK) simulator [13]. However, since there is a large time cost of the simulation method with the same level of accuracy, this method may be replaced by the analytic approximation method. Mendell and Elston (ME) proposed a univariate conditioning approach earlier [14], and then, the researchers made a series of improvements. In the above probit-based ICLV model, the MACML approach of Bhat [15] is adopted, which combines the Switzer–Solow–Joe (SSJ) approximation for the MVNCD function with the CML inference approach for MNP models. Here, we intend to apply a new LDLT matrix-based analytic approximation method proposed by Bhat in 2018 [16] to evaluate the MVNCD function, which is demonstrated to have better performance than the SSJ. The method can effectively reduce difficulties in convergence and covariance matrix computation that occur routinely in the maximum-likelihood estimation of choice models with analytically intractable likelihood functions. Meanwhile, applying the LDLT decomposition method to the algorithms also assists in ensuring no substantial speed reduction in each iteration.In the process of the maximum-likelihood estimation (MLE) for models, programming analytical gradients are generally needed to improve the speed and accuracy at convergence. However, since the derivation and coding of analytical gradients are usually complicated and time-consuming for complex models, the introduction of automatic differentiation (AD) into the model estimation process enables the acquisition of gradients without coding effort, which significantly reduces the workload of researchers. The AD method takes advantage of the fact that every computer program will execute a sequence of elementary arithmetic operations and basic functions, no matter how complicated it is. The MLE procedure usually adopts a reverse mode of AD, which is first published in 1976 by Linnainmaa [17]. Compared with numerical differentiation, it has obvious advantages in terms of calculating higher-order derivatives and partial derivatives with multiple inputs. In the past decades, the theory of AD has been considerably developed. Nevertheless, it has limited applications in the field of econometrics and statistics up to now [18, 19], which provides opportunities for further investigation in this respect.Furthermore, the graph computation technique based on tensors represented as multidimensional arrays in the TensorFlow platform can avoid loop operations based on the two-dimensional matrix data structure used in conventional computational platforms. TensorFlow should have an advantage in computational speed over conventional platforms when the LDLT-based analytic approximation method is applied within a sample. On the one hand, the graph computation technique based on tensors can ensure that all the calculation process in LDLT-based analytic approximation is completely realized at a speed of C language codes, rather than that of scripting language codes to operate matrices repeatedly in multiple loops. On the other hand, in MLE procedures using matrix-based analytic approximations, a tensor with more dimensions is desirable because each individual observation in the sample needs a matrix-based operation for its log-likelihood computation and such computation needs to be repeated as many times as the sample size. Thus, a computational platform with a tensor data structure can avoid repeated matrix operations in a large loop and further contribute to computational speed improvement. In this study, some efforts will be made to explore the feasibility of MNP-based ICLV model estimation using LDLT-based analytic approximation and tensor data structure on the TensorFlow platform.The rest of this study will be organized as follows. The second section elaborates on the MNP-based ICLV model and LDLT-based analytic approximation method and then introduces the AD mechanism. In the third section, the simulation experiments are designed for validating the three-alternative and four-alternative ICLV models. The next section presents the evaluation results of the MVNCD function and the estimation results for the ICLV model with different analytic approximation methods (ME, BME, and TVBS). Conclusions and discussions are summarized in the last section.


## 2. Methodology

### 2.1. MNP-Based ICLV Model

The MNP-based integrated choice and latent variable (ICLV) model proposed by Bhat and Dubey in 2014 [10] contains the following three components. Assume there are L latent variables z∗ (index l=1,2,…,L), G ordinal indicator variables y (index g=1,2,…G), I alternatives’ utilities U (index i=1,2,…I), and Q individuals (index q=1,2,…Q) in this model. To simplify the presentation, the matrix form is taken in the formulas below and the index of the variables is not shown.


#### 2.1.1. Latent Variable Structural Equation Model

(1)z∗=αw+η,η∼N0L,Γ,where w is a D˜×1 vector of observed covariates affecting z∗ (not including a constant) and α is a corresponding L×D˜ matrix of coefficients (i.e., exogenous variable loadings on z∗). η is a L×1 vector of errors assumed to be standard multivariate normally distributed: η∼N0L,Γ, where Γ is a correlation matrix. Each ηl of η is assumed to be the standard normal distribution. For the structural equation model, the same exogenous variables w can be used for different latent variables zl∗.


#### 2.1.2. Latent Variable Measurement Equation Model

(2)y=δ+dz∗+ξ,ψlow<y<ψup,ξ∼N0G,Σ,where y, δ, and ξ are G×1 vectors that represent observed ordered indicator-dependent variables, intercept in measurement equation, and error items in measurement equation, respectively. d is a G×L matrix of coefficients representing the effect of latent variables z∗ on observed indicators y. Also, let Σ be the correlation matrix of ξ=ξ1,ξ2,…,ξG, in which the identification of diagonal is not considered. Here, the normalization on the error terms ξ is needed for identification, as in the usual ordered-response model. The ordered indicator y can be horizontally partitioned into corresponding ordinal categories by setting thresholds. As per Bhat in 2014, ψlow is a G×1 vector that stacks the lower thresholds ψg,ng−1g=1,2,…,G and ψup is another vector that stacks the upper thresholds ψg,ngg=1,2,…,G.


#### 2.1.3. Choice Model (Multinomial Probit Model)

(3)U=xβ+γφz∗+ε=xβ+λz∗+ε,whereλ=γφ,ε∼N0I,Λ,where U is (I×1) vector of alternative utilities, x is (I×D) matrix of exogenous variables in the choice model, and ε is (I×1) utility error vector that followed multivariate normal distribution whose variance-covariance matrix is Λ. In the equation, φ is a ∑i=1INi×L matrix of variables interacting with latent variables and γ is an I×∑i=1INi matrix of coefficients capturing the effects of latent variables and their interactions with exogenous variables. In a case without interactions, λI×L matrix of coefficients can be specified when λ=γφ and the matrix φ is reduced to an identity matrix of size L. Since only the covariance matrix of error differences is estimable, II−1/2−1 elements can be identified in Λ after normalizing the covariance matrix.The MNP-based ICLV model framework can be presented in Figure1.Figure 1
MNP-based ICLV model framework.In the process of the model estimation,z∗ in equations (2) and (3) can be substituted into the right side of equation (1) to further simplify expressions. We assume that the error vectors ξ, ε, and η are independent of each other. The following system can be obtained:(4)y=δ+dαw+η+ξ=δ+dαω+dη+ξ,U=xβ+λz∗+ε=xβ+λαw+η+ε=xβ+λαw+λη+ε.Consider theG+I vector YUG+I=yG′,UI′′. The following is defined:(5)B=δ+dαwxβ+λαw,Ω=dΓd′+ΣdΓλ′λΓd′λΓλ′+Λ.Then,YU∼MVNG+IB,Ω.To estimate the model, we need to define matrixM of size G+I−1×G+I firstly to transform the choice model into the utility difference form. When the alternative i is chosen, the matrix G×G in Mi with the first G rows and first G columns is an identity matrix, and the rest matrix I−1×I in Mi is an identity matrix of I−1 size with an extra column of −1 added as the ith column. For example, when the alternative i=2 is chosen, Mi=1000001000001−10000−11 if G=2 and I=3. Next, let Yu∼MVNG+I−1B˜,Ω˜, where B˜=MB,Ω˜=MΩM′. After that, the matrix dimension is reduced to G+I−1.To integrate the multivariate OP model and MNP model to estimate the parameters, the upper and lower threshold vectors are defined:ψ˜low=ψlow′,−∞I−1′ and ψ˜up=ψup′,0I−1′, where both −∞I−1 and 0I−1 are I−1×1-column vectors that represent negative infinities and zeros, respectively. Thus, the likelihood function of the ICLV model may be written as follows:(6)Lθ=Prψ˜low≤u˜≤ψ˜up=∫Du˜fG+I−1u˜B˜,Ω˜du,where Du¯=u˜:ψ˜low≤u˜≤ψ˜up is the integration domain and fG+I−1. is the multivariate normal density function of dimension G+I−1. The integral dimension of the above likelihood function does not increase with the increase in the latent variables z∗.To further proceed, the composite marginal likelihood (CML) approach can be used to simplify the model estimation process when the full likelihood function may be infeasible to evaluate. The CML estimator is the one that maximizes the compounded probability of all pairwise events, which obtains much easier-to-compute, lower-dimensional marginal likelihood. In the above case, the pairwise marginal likelihood function may be written for the ICLV model as follows:(7)LCMLθ=∏g=1G−1∏g′=g+1GPrjg=ng,jg′=ng′×∏g=1GPrjg=ng,i=m,where jg represents the index for the ordinal outcome category for the gth ordinal variables. Then, consider Sg (similar to M), which is a I×G+I−1 selection matrix to select the gth ordered variables and all alternatives of the choice model. Let ΦE:;Δ∗ be the multivariate standard normal cumulative matrix function of dimension E and correlation matrix Δ∗Δ∗=ωΔ−1ΔωΔ−1, where ωΔ is the diagonal matrix of the standard deviation of the covariance matrix Δ. The LCMLθ function is expanded in detail as follows:(8)LCMLθ=∏g=1G−1∏g′=g+1GΦ2ϑg,up,ϑg′,up,vgg′−Φ2ϑg,up,ϑg′,low,vgg′−Φ2ϑg,low,ϑg′,up,vgg′+Φ2ϑg,low,ϑg′,low,vgg′×∏g=1GΦIωΩ˘g−1ψ˘g,up−B˘g;Ωˇg∗−ΦIωΩ˘g−1ψ˘g,low−B˘g;Ωˇg∗,where ϑg,up=ψ˜upg−B˜g/Ω˜gg,ϑg,low=ψ˜lowg−B˜g/Ω˜gg,vgg′=Ω˜gg′/Ω˜gg∗Ω˜gg′, ψ˘g,low=Sgψlow',0I−1′′,ψ˘g,up=Sgψ˜up,B˘g=SgB˜,Ω˘g=SgΩ˜Sg'.In equations (7) and (8), the CML function can be divided into two components: the first component corresponds to each pair of ordinal indicators (of which the integral dimension of MVNCD is 2), while the second component corresponds to each pair of the choice outcome and an ordinal indicator outcome (of which the integral dimension is equal to I). Accordingly, the maximum integral dimension of the MVNCD function in equation (8) does not exceed I. It is pivotal to solve the multidimensional integral in the process of model estimation. The matrix-based analytic approximation method proposed by Bhat [16] can be applied to evaluate the MVNCD function, as detailed in a later section. At last, accounting for the index q for individuals, the pairwise marginal log-likelihood function is defined as logLCMLθ=∑q=1QlogLCML,qθ.The pairwise estimatorθCML obtained by maximizing the logarithm of the pairwise marginal likelihood function with respect to the vector θ is consistent and asymptotically normally distributed with asymptotic mean θ and covariance matrix given by the inverse of Godambe’s [20] sandwich information matrix Gθ [21]:(9)VCMLθ=Gθ−1=Hθ−1JθHθ−1,where Hθ and Jθ can be estimated straightforwardly at the CML estimate (θ^CML):H^θ^=−∑q=1Q∂2logLCML,qθ/∂θ∂θ′θ^CML and(10)J^θ^=∑q=1Q∂logLCML,qθ∂θ∂logLCML,qθ∂θ′θ^CML.


### 2.2. Matrix-Based Analytic Approximation Method

The computation of the MVNCD function can be applied for the multidimensional integral of the ICLV model above. In the previous ICLV model of Bhat, the MACML approach in 2011 [15] is adopted. In this section, a new matrix-based analytic approximation method proposed by Bhat in 2018 [16] can replace the MACML approach to evaluate the MVNCD function for acquiring better estimation results [22, 23].There are five methods for approximating the MVNCD function mentioned and summarized in the literature of Bhat in 2018. As per Bhat, the TVBS method is recommended as the evaluation approach for the MVNCD function with the highest accuracy, and the BME method has a significant speed advantage over these evaluation approaches. Besides, as the foundation of the series of analytic approximation methods, the ME method can also be included in the evaluation system in this study. Thus, these three typical methods (ME, BME, and TVBS) will be selected for elaboration below. Conceptually, the ME method is based on univariate conditioning, while BME and TVBS are based on a bivariate conditioning mechanism.Lete=e1,e2,…,eK be a multivariate normally distributed random vector with zero means, variances of 1, and a correlation matrix Σ=1ρ12⋯ρ1Kρ121⋯ρ2K⋮⋮⋱⋮ρ1Kρ2K⋯1 (K > 2), where K is the dimensionality of the MVNCD function. The K-dimensional MVNCD function can be expressed as follows: ΦKe1,e2,…,eK:Σ=Pe<V=Pe1<V1,e2<V2,…,eK<VK. To facilitate the understanding, assume that t represents e<V and tn represents e1<V1, and then, ΦKe:Σ=Pe<V=Pt=Pt1,t2,…,tK.The joint probability of the MVNCD function of the three methods (ME, BME, and TVBS) may be expressed as follows. All the joint probability may be written as the product of a marginal probability and conditional probabilities, and the normal distribution can be used to approximate the skew normal distribution for conditional probabilities. In the equations below, the superscriptn of tn is the number of truncations; nmax is the maximum number of truncations, and a lower nmax means a higher accuracy of the analytic approximation method. Let N be an intermediate parameter used to measure K: K is equal to N for the ME method, and K is discussed as being odd K=2N−1 or even K=2N for the BME (or TVBS) method. The expected values μk and variances σk2 of single-sided truncations in the following equations can be obtained by the properties of truncated multivariate normal distribution, which is not detailed here (refer to Section 2.1 of Bhat in 2018).


#### 2.2.1. The ME MethodN≥2

(11)Pt=Pt1,t2,…,tN,=Pt1×∏n=1N−1Ptn+1|t1t2,…,tn,≈Pt1×∏n=1N−1Ptn+1n,≈ΦV1×∏n=1N−1ΦVn+1−μn+1σn+1.The method relies on a single-sided truncation of a multivariate normal distribution in which some variables are truncated, while others are not. After a single variable is truncated, the residual variables follow a skew normal distribution in the ME method.Ptn+1n is defined as the univariate conditional probabilities of tn+1 after n times truncations, which follows a univariate normal distribution (being used for approximating the skew normal distribution). Thus, the above Pt can be transformed into the product of recursive univariate normal cumulative distribution function (CDF) (represented by Φ.). In equation (11), K=N and nmax=N−1.For instance, whenK=6, the probability of multivariate normal distribution can be written as follows: Pt1,t2,t3,t4,t5,t6≈Pt1×Pt21×Pt32×Pt43×Pt54×Pt65=Pe1<V1×Pe2∗<V2|e1<V1×Pe3∗<V3|e2∗<V2×Pe4∗<V4|e3∗<V3×Pe5∗<V5|e4∗<V4×Pe6∗<V6|e5∗<V5, where ek∗ follows the skew normal distribution after k−1 times truncations. Thereinto, t2 was truncated once and Pt21=Pe2∗<V2|e1<V1≈ΦV2−μ2/σ2. z1=e1|e1<V1 is defined, and it can be deduced that the expected values μ2=ρ12Ez1 and variances σ22=ρ122Varz1+1−ρ122, where Ez1=μ1+σ1λ1=λ1, Varz1=σ121+λ1V1−λ12=1+λ1V1−λ12, and λ1=−ϕV1/ΦV1. Similar to the calculation principle of ΦV2−μ2/σ2, ΦVn+1−μn+1/σn+1 in equation (11) can be calculated by iteratively updating the expected value and covariance matrix. In the equation above, the maximum number of truncations nmax is 5.


#### 2.2.2. The Bivariate ME (BME) MethodN≥2

WhenK=2N is even,(12)Pt=Pt1,t2,…,t2N,=Pt1,t2×∏n=1N−1Pt2n+1,t2n+2|t1,t2,…,t2n,≈Pt1,t2×∏n=1N−1Pt2n+12n,t2n+22n,≈Φ2V1,V2,ρ12×∏n=1N−1Φ2V2n+1−μ2n+1σ2n+1,V2n+2−μ2n+2σ2n+2,ρ2n+1,2n+2∗.WhenK=2N−1 is odd,(13)Pt=Pt1,t2,…,t2N−1,=Pt1,t2×∏n=1N−2Pt2n+1,t2n+2|t1,t2,…,t2n×Pt2N−1|t1,t2,…,t2N−2,≈Pt1,t2×∏n=1N−2Pt2n+12n,t2n+22n×Pt2N−12N−2,≈Φ2V1,V2,ρ12×∏n=1N−2Φ2V2n+1−μ2n+1σ2n+1,V2n+2−μ2n+2σ2n+2,ρ2n+1,2n+2∗×ΦV2N−1−μ2N−1σ2N−1.The BME method is in accordance with a bivariate truncation scheme, assuming that the marginal distribution of the third and fourth untruncated variables, given the first two are truncated, remains bivariate normal, and so on. In the equation above, similar to the ME method,Pt2n+12n,t2n+22n uses the bivariate normal distribution to approximate the bivariate skew normal distribution. The conditional probability Pt of this method for MVNCD function is a product of bivariate normal CDF Φ2. when the dimension is even and a product of Φ2. and Φ. when the dimension is odd. Here, ρ2n+1,2n+2∗ is the updated correlation coefficient between error terms after bivariate truncation, and the calculation principle of the expected values μk and variances σk2 is similar to the ME method. As the extension of the ME method, the BME method truncates two dimensions each time and nmax times in total (also does the TVBS method), where nmax=N−1 in equations (12) and (13).Similarly, whenK=6, the Pt using the BME method can be expressed as follows: Pt1,t2,t3,t4,t5,t6≈Pt1,t2×Pt32,t42×Pt54,t64. Thereinto, t3 and t4are truncated once n=1, and t5 and t6 are truncated twice n=2. The maximum number of truncations nmax is 2, which is smaller than the nmax of the ME method. This indicates that the accuracy of the bivariate truncation scheme is better than the univariate truncation scheme.


#### 2.2.3. The Two-Variate Bivariate Screening (TVBS) Method

WhenK=2NN≥2 is even,(14)Pt=Pt1,t2,…,t2N,=Pt1,t2,t3,t4×∏n=1N−2Pt2n+1,t2n+2,t2n+3,t2n+4|t1,t2,…,t2nPt2n+1,t2n+2|t1,t2,…,t2n,≈Pt1,t2,t3,t4×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32n,t2n+42nPt2n+12n,t2n+22n,=Pt1,t2,t3×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32nPt2n+12n×Pt2N−12N−2,t2N2N−2Pt2N−12N−2,≈Φ3V1,V2,V3:0,Σ1:3,1:3×∏n=1N−2Φ3V2n+1,V2n+2,V2n+3:μ˜,Ω˜ΦV2n+1−μ2n+1/σ2n+1×Φ2V2N−1−μ2N−1/σ2N−1,V2N−μ2N/σ2N,ρ2N−1,2N∗ΦV2N−1−μ2N−1/σ2N−1.WhenK=2N−1N≥3 is odd,(15)Pt=Pt1,t2,…,t2N−1,=Pt1,t2,t3,t4×∏n=1N−3Pt2n+1,t2n+2,t2n+3,t2n+4|t1,t2,…,t2nPt2n+1,t2n+2,t1,t2,…,t2n×Pt2N−3,t2N−2,t2N−1|t1,t2,…,t2N−4Pt2N−3,t2N−2|t1,t2,…,t2N−4,≈Pt1,t2,t3,t4×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32n,t2n+42nPt2n+12n,t2n+22n×Pt2N−32N−4,t2N−22N−4,t2N−12N−4Pt2N−32N−4,t2N−22N−4,=Pt1,t2,t3×∏n=1N−3Pt2n+12n,t2n+22n,t2n+32nPt2n+12n,≈Φ3V1,V2,V3:0,Σ1:3,1:3×∏n=1N−2Φ3V2n+1,V2n+2,V2n+3:μ˜,Ω˜ΦV2n+1−μ2n+1/σ2n+1.Similarly, the TVBS method combines the bivariate truncation scheme, in which the bivariate distribution is assumed to be rectangle-screened normal (RSN) distribution [24]. The corresponding conditional CDF at each step can be expressed as a ratio of a four-variate normal CDF Φ4. and a two-variate normal CDF Φ2.. Since the calculation of four-variate normal CDF is time-consuming, the TVBS method adopts an approximation method by taking the trivariate CDF of the first three variates P123. Given the first three variables, the CDF of the skew normal distribution of the fourth variable is denoted as P4|123. Accordingly, a four-variate normal CDF is obtained: P1234=P123∗P4|123=P123∗P34|123/P3|12. Further, the MVNCD function is composed of the product of Φ3. (denotes trivariate normal CDF) and Φ3./Φ. when K is an odd number, and the product of Φ3., Φ3./Φ., and Φ2./Φ. when K is odd. In equations (14) and (15), μ˜ and Ω˜ represent the expected value and covariance matrix updated after bivariate truncation, respectively. Additionally, nmax=N−1 (K is even) and nmax=N−2 (K is odd) for the TVBS method.WhenK=6, the Pt using the TVBS method can be expressed as follows: Pt1,t2,t3,t4,t5,t6≈Pt1,t2,t3×Pt32,t42,t52/Pt32×Pt54,t64/Pt54, where N=3 and nmax=2. Here, t3, t4, and t5 conducted one truncation n=1 in Pt32,t42,t52/Pt32, while t5 and t6 conducted two truncations n=2 in Pt54,t64/Pt54.Furthermore, the matrix-based analytic approximation method applies an LDLT decomposition method for the correlation matrix followed by rank 1 (ME) or rank 2 (BME and method) updates of the factorization. The implementation is much easier to code and more computationally efficient.As we know from the above discussion, the multidimensional integral for the analytical approximation method of MNVCD can be completely converted into integrals with no more than three dimensions. At present, the CDF of bivariate normal (BVN) distribution and trivariate normal (TVN) distribution is generally realized by numerical methods based on the Legendre–Gauss quadrature integral approximation. This numerical integration method adopted in most existing packages is relatively old with low levels of accuracy and speed. In this study, referring to the study of Genz in 2004 [25], a new numerical integration method is used to replace the original method. For BVN probabilities, we adopted the algorithm proposed by Drezner and Wesolowsky [26] with a higher-order Ox4 Taylor expansion; for TVN probabilities, the second Plackett formula method studied by Gassmann [27] is applied (not be detailed here). Accordingly, after the above algorithms are applied along with matrix operations, the performance of both functions is improved in terms of accuracy and speed.


### 2.3. Automatic Differentiation (AD)

Automatic differentiation (AD) is a chain rule-based technique for evaluating the derivatives with respect to the input variables of functions defined by a computer program [28]. There are two different AD strategies: the forward mode and the reverse mode. The advantage of the reverse mode over the forward mode will be highlighted when the number of parameters is large, which is more suitable for the MLE procedure.For a simple example, consider a binary probit model:U=a+bx+ε,y=IU>0, and ε∼N0,1. The logarithmic likelihood function of the model can be written as LLθ=LLa,b=lnΦ2y−1a+bx. In the reverse mode of AD, the derivative is computed with respect to each sub-expression recursively in the chain rule:(16)∂LL∂θ=∂LL∂v1∂v1∂θ=∂y∂v2∂v2∂v1∂v1∂θ=∂y∂v3∂v3∂v2∂v2∂v1∂v1∂θ=….The derivative ofLL with respect to vi is defined as v¯i=∂LL/∂vi, and the reverse accumulation with computation graph is illustrated in Figure 2. In this figure, the partial derivatives are estimated as follows: ∂LL/∂a=a¯=2y−1∗φ2y−1∗a+bx/Φ2y−1∗a+bx and ∂LL/∂b=b¯=x∗2y−1∗φ2y−1∗a+bx/Φ2y−1∗a+bx. Thereinto, φ. has a closed form, and Φ. can be approximated by the numerical integration method.Figure 2
Example of reverse-mode AD in the computation graph.TensorFlow (https://www.tensorflow.org/), an open-source framework for large-scale distributed numerical computation, is widely applied in the areas of machine learning and deep neural network during early explorations. The main feature of the TensorFlow framework is based on a data flow graph and uses tensor as the core component to denote data. The framework builds a gradient graph corresponding to the original computation graph by backward propagation and chain rule of derivative values, which is more mature and flexible than other programming modes in the application of automatic gradient. Due to the advantages of TensorFlow, we applied it for developing econometric models in this study.AD has obvious advantages when it is applied to the MLE of parameters in econometric models. Using this technique, modelers do not need to derive and code the complex analytic gradient expression of the log-likelihood function, which greatly reduces the workload. When the CML approach is used to estimate model parameters, the Hessian matrix and Jacobian matrix of estimated parameters need to be calculated. The errors and complexity of the numerical differentiation method increase when second-order derivatives are computed, and the speed is low in calculating partial derivatives of a function with respect to many inputs. Therefore, the use of automatic differentiation instead of numerical differentiation can aid in attaining more accurate estimation results with less computational time.


## 2.1. MNP-Based ICLV Model

The MNP-based integrated choice and latent variable (ICLV) model proposed by Bhat and Dubey in 2014 [10] contains the following three components. Assume there are L latent variables z∗ (index l=1,2,…,L), G ordinal indicator variables y (index g=1,2,…G), I alternatives’ utilities U (index i=1,2,…I), and Q individuals (index q=1,2,…Q) in this model. To simplify the presentation, the matrix form is taken in the formulas below and the index of the variables is not shown.


### 2.1.1. Latent Variable Structural Equation Model

(1)z∗=αw+η,η∼N0L,Γ,where w is a D˜×1 vector of observed covariates affecting z∗ (not including a constant) and α is a corresponding L×D˜ matrix of coefficients (i.e., exogenous variable loadings on z∗). η is a L×1 vector of errors assumed to be standard multivariate normally distributed: η∼N0L,Γ, where Γ is a correlation matrix. Each ηl of η is assumed to be the standard normal distribution. For the structural equation model, the same exogenous variables w can be used for different latent variables zl∗.


### 2.1.2. Latent Variable Measurement Equation Model

(2)y=δ+dz∗+ξ,ψlow<y<ψup,ξ∼N0G,Σ,where y, δ, and ξ are G×1 vectors that represent observed ordered indicator-dependent variables, intercept in measurement equation, and error items in measurement equation, respectively. d is a G×L matrix of coefficients representing the effect of latent variables z∗ on observed indicators y. Also, let Σ be the correlation matrix of ξ=ξ1,ξ2,…,ξG, in which the identification of diagonal is not considered. Here, the normalization on the error terms ξ is needed for identification, as in the usual ordered-response model. The ordered indicator y can be horizontally partitioned into corresponding ordinal categories by setting thresholds. As per Bhat in 2014, ψlow is a G×1 vector that stacks the lower thresholds ψg,ng−1g=1,2,…,G and ψup is another vector that stacks the upper thresholds ψg,ngg=1,2,…,G.


### 2.1.3. Choice Model (Multinomial Probit Model)

(3)U=xβ+γφz∗+ε=xβ+λz∗+ε,whereλ=γφ,ε∼N0I,Λ,where U is (I×1) vector of alternative utilities, x is (I×D) matrix of exogenous variables in the choice model, and ε is (I×1) utility error vector that followed multivariate normal distribution whose variance-covariance matrix is Λ. In the equation, φ is a ∑i=1INi×L matrix of variables interacting with latent variables and γ is an I×∑i=1INi matrix of coefficients capturing the effects of latent variables and their interactions with exogenous variables. In a case without interactions, λI×L matrix of coefficients can be specified when λ=γφ and the matrix φ is reduced to an identity matrix of size L. Since only the covariance matrix of error differences is estimable, II−1/2−1 elements can be identified in Λ after normalizing the covariance matrix.The MNP-based ICLV model framework can be presented in Figure1.Figure 1
MNP-based ICLV model framework.In the process of the model estimation,z∗ in equations (2) and (3) can be substituted into the right side of equation (1) to further simplify expressions. We assume that the error vectors ξ, ε, and η are independent of each other. The following system can be obtained:(4)y=δ+dαw+η+ξ=δ+dαω+dη+ξ,U=xβ+λz∗+ε=xβ+λαw+η+ε=xβ+λαw+λη+ε.Consider theG+I vector YUG+I=yG′,UI′′. The following is defined:(5)B=δ+dαwxβ+λαw,Ω=dΓd′+ΣdΓλ′λΓd′λΓλ′+Λ.Then,YU∼MVNG+IB,Ω.To estimate the model, we need to define matrixM of size G+I−1×G+I firstly to transform the choice model into the utility difference form. When the alternative i is chosen, the matrix G×G in Mi with the first G rows and first G columns is an identity matrix, and the rest matrix I−1×I in Mi is an identity matrix of I−1 size with an extra column of −1 added as the ith column. For example, when the alternative i=2 is chosen, Mi=1000001000001−10000−11 if G=2 and I=3. Next, let Yu∼MVNG+I−1B˜,Ω˜, where B˜=MB,Ω˜=MΩM′. After that, the matrix dimension is reduced to G+I−1.To integrate the multivariate OP model and MNP model to estimate the parameters, the upper and lower threshold vectors are defined:ψ˜low=ψlow′,−∞I−1′ and ψ˜up=ψup′,0I−1′, where both −∞I−1 and 0I−1 are I−1×1-column vectors that represent negative infinities and zeros, respectively. Thus, the likelihood function of the ICLV model may be written as follows:(6)Lθ=Prψ˜low≤u˜≤ψ˜up=∫Du˜fG+I−1u˜B˜,Ω˜du,where Du¯=u˜:ψ˜low≤u˜≤ψ˜up is the integration domain and fG+I−1. is the multivariate normal density function of dimension G+I−1. The integral dimension of the above likelihood function does not increase with the increase in the latent variables z∗.To further proceed, the composite marginal likelihood (CML) approach can be used to simplify the model estimation process when the full likelihood function may be infeasible to evaluate. The CML estimator is the one that maximizes the compounded probability of all pairwise events, which obtains much easier-to-compute, lower-dimensional marginal likelihood. In the above case, the pairwise marginal likelihood function may be written for the ICLV model as follows:(7)LCMLθ=∏g=1G−1∏g′=g+1GPrjg=ng,jg′=ng′×∏g=1GPrjg=ng,i=m,where jg represents the index for the ordinal outcome category for the gth ordinal variables. Then, consider Sg (similar to M), which is a I×G+I−1 selection matrix to select the gth ordered variables and all alternatives of the choice model. Let ΦE:;Δ∗ be the multivariate standard normal cumulative matrix function of dimension E and correlation matrix Δ∗Δ∗=ωΔ−1ΔωΔ−1, where ωΔ is the diagonal matrix of the standard deviation of the covariance matrix Δ. The LCMLθ function is expanded in detail as follows:(8)LCMLθ=∏g=1G−1∏g′=g+1GΦ2ϑg,up,ϑg′,up,vgg′−Φ2ϑg,up,ϑg′,low,vgg′−Φ2ϑg,low,ϑg′,up,vgg′+Φ2ϑg,low,ϑg′,low,vgg′×∏g=1GΦIωΩ˘g−1ψ˘g,up−B˘g;Ωˇg∗−ΦIωΩ˘g−1ψ˘g,low−B˘g;Ωˇg∗,where ϑg,up=ψ˜upg−B˜g/Ω˜gg,ϑg,low=ψ˜lowg−B˜g/Ω˜gg,vgg′=Ω˜gg′/Ω˜gg∗Ω˜gg′, ψ˘g,low=Sgψlow',0I−1′′,ψ˘g,up=Sgψ˜up,B˘g=SgB˜,Ω˘g=SgΩ˜Sg'.In equations (7) and (8), the CML function can be divided into two components: the first component corresponds to each pair of ordinal indicators (of which the integral dimension of MVNCD is 2), while the second component corresponds to each pair of the choice outcome and an ordinal indicator outcome (of which the integral dimension is equal to I). Accordingly, the maximum integral dimension of the MVNCD function in equation (8) does not exceed I. It is pivotal to solve the multidimensional integral in the process of model estimation. The matrix-based analytic approximation method proposed by Bhat [16] can be applied to evaluate the MVNCD function, as detailed in a later section. At last, accounting for the index q for individuals, the pairwise marginal log-likelihood function is defined as logLCMLθ=∑q=1QlogLCML,qθ.The pairwise estimatorθCML obtained by maximizing the logarithm of the pairwise marginal likelihood function with respect to the vector θ is consistent and asymptotically normally distributed with asymptotic mean θ and covariance matrix given by the inverse of Godambe’s [20] sandwich information matrix Gθ [21]:(9)VCMLθ=Gθ−1=Hθ−1JθHθ−1,where Hθ and Jθ can be estimated straightforwardly at the CML estimate (θ^CML):H^θ^=−∑q=1Q∂2logLCML,qθ/∂θ∂θ′θ^CML and(10)J^θ^=∑q=1Q∂logLCML,qθ∂θ∂logLCML,qθ∂θ′θ^CML.


## 2.1.1. Latent Variable Structural Equation Model

(1)z∗=αw+η,η∼N0L,Γ,where w is a D˜×1 vector of observed covariates affecting z∗ (not including a constant) and α is a corresponding L×D˜ matrix of coefficients (i.e., exogenous variable loadings on z∗). η is a L×1 vector of errors assumed to be standard multivariate normally distributed: η∼N0L,Γ, where Γ is a correlation matrix. Each ηl of η is assumed to be the standard normal distribution. For the structural equation model, the same exogenous variables w can be used for different latent variables zl∗.


## 2.1.2. Latent Variable Measurement Equation Model

(2)y=δ+dz∗+ξ,ψlow<y<ψup,ξ∼N0G,Σ,where y, δ, and ξ are G×1 vectors that represent observed ordered indicator-dependent variables, intercept in measurement equation, and error items in measurement equation, respectively. d is a G×L matrix of coefficients representing the effect of latent variables z∗ on observed indicators y. Also, let Σ be the correlation matrix of ξ=ξ1,ξ2,…,ξG, in which the identification of diagonal is not considered. Here, the normalization on the error terms ξ is needed for identification, as in the usual ordered-response model. The ordered indicator y can be horizontally partitioned into corresponding ordinal categories by setting thresholds. As per Bhat in 2014, ψlow is a G×1 vector that stacks the lower thresholds ψg,ng−1g=1,2,…,G and ψup is another vector that stacks the upper thresholds ψg,ngg=1,2,…,G.


## 2.1.3. Choice Model (Multinomial Probit Model)

(3)U=xβ+γφz∗+ε=xβ+λz∗+ε,whereλ=γφ,ε∼N0I,Λ,where U is (I×1) vector of alternative utilities, x is (I×D) matrix of exogenous variables in the choice model, and ε is (I×1) utility error vector that followed multivariate normal distribution whose variance-covariance matrix is Λ. In the equation, φ is a ∑i=1INi×L matrix of variables interacting with latent variables and γ is an I×∑i=1INi matrix of coefficients capturing the effects of latent variables and their interactions with exogenous variables. In a case without interactions, λI×L matrix of coefficients can be specified when λ=γφ and the matrix φ is reduced to an identity matrix of size L. Since only the covariance matrix of error differences is estimable, II−1/2−1 elements can be identified in Λ after normalizing the covariance matrix.The MNP-based ICLV model framework can be presented in Figure1.Figure 1
MNP-based ICLV model framework.In the process of the model estimation,z∗ in equations (2) and (3) can be substituted into the right side of equation (1) to further simplify expressions. We assume that the error vectors ξ, ε, and η are independent of each other. The following system can be obtained:(4)y=δ+dαw+η+ξ=δ+dαω+dη+ξ,U=xβ+λz∗+ε=xβ+λαw+η+ε=xβ+λαw+λη+ε.Consider theG+I vector YUG+I=yG′,UI′′. The following is defined:(5)B=δ+dαwxβ+λαw,Ω=dΓd′+ΣdΓλ′λΓd′λΓλ′+Λ.Then,YU∼MVNG+IB,Ω.To estimate the model, we need to define matrixM of size G+I−1×G+I firstly to transform the choice model into the utility difference form. When the alternative i is chosen, the matrix G×G in Mi with the first G rows and first G columns is an identity matrix, and the rest matrix I−1×I in Mi is an identity matrix of I−1 size with an extra column of −1 added as the ith column. For example, when the alternative i=2 is chosen, Mi=1000001000001−10000−11 if G=2 and I=3. Next, let Yu∼MVNG+I−1B˜,Ω˜, where B˜=MB,Ω˜=MΩM′. After that, the matrix dimension is reduced to G+I−1.To integrate the multivariate OP model and MNP model to estimate the parameters, the upper and lower threshold vectors are defined:ψ˜low=ψlow′,−∞I−1′ and ψ˜up=ψup′,0I−1′, where both −∞I−1 and 0I−1 are I−1×1-column vectors that represent negative infinities and zeros, respectively. Thus, the likelihood function of the ICLV model may be written as follows:(6)Lθ=Prψ˜low≤u˜≤ψ˜up=∫Du˜fG+I−1u˜B˜,Ω˜du,where Du¯=u˜:ψ˜low≤u˜≤ψ˜up is the integration domain and fG+I−1. is the multivariate normal density function of dimension G+I−1. The integral dimension of the above likelihood function does not increase with the increase in the latent variables z∗.To further proceed, the composite marginal likelihood (CML) approach can be used to simplify the model estimation process when the full likelihood function may be infeasible to evaluate. The CML estimator is the one that maximizes the compounded probability of all pairwise events, which obtains much easier-to-compute, lower-dimensional marginal likelihood. In the above case, the pairwise marginal likelihood function may be written for the ICLV model as follows:(7)LCMLθ=∏g=1G−1∏g′=g+1GPrjg=ng,jg′=ng′×∏g=1GPrjg=ng,i=m,where jg represents the index for the ordinal outcome category for the gth ordinal variables. Then, consider Sg (similar to M), which is a I×G+I−1 selection matrix to select the gth ordered variables and all alternatives of the choice model. Let ΦE:;Δ∗ be the multivariate standard normal cumulative matrix function of dimension E and correlation matrix Δ∗Δ∗=ωΔ−1ΔωΔ−1, where ωΔ is the diagonal matrix of the standard deviation of the covariance matrix Δ. The LCMLθ function is expanded in detail as follows:(8)LCMLθ=∏g=1G−1∏g′=g+1GΦ2ϑg,up,ϑg′,up,vgg′−Φ2ϑg,up,ϑg′,low,vgg′−Φ2ϑg,low,ϑg′,up,vgg′+Φ2ϑg,low,ϑg′,low,vgg′×∏g=1GΦIωΩ˘g−1ψ˘g,up−B˘g;Ωˇg∗−ΦIωΩ˘g−1ψ˘g,low−B˘g;Ωˇg∗,where ϑg,up=ψ˜upg−B˜g/Ω˜gg,ϑg,low=ψ˜lowg−B˜g/Ω˜gg,vgg′=Ω˜gg′/Ω˜gg∗Ω˜gg′, ψ˘g,low=Sgψlow',0I−1′′,ψ˘g,up=Sgψ˜up,B˘g=SgB˜,Ω˘g=SgΩ˜Sg'.In equations (7) and (8), the CML function can be divided into two components: the first component corresponds to each pair of ordinal indicators (of which the integral dimension of MVNCD is 2), while the second component corresponds to each pair of the choice outcome and an ordinal indicator outcome (of which the integral dimension is equal to I). Accordingly, the maximum integral dimension of the MVNCD function in equation (8) does not exceed I. It is pivotal to solve the multidimensional integral in the process of model estimation. The matrix-based analytic approximation method proposed by Bhat [16] can be applied to evaluate the MVNCD function, as detailed in a later section. At last, accounting for the index q for individuals, the pairwise marginal log-likelihood function is defined as logLCMLθ=∑q=1QlogLCML,qθ.The pairwise estimatorθCML obtained by maximizing the logarithm of the pairwise marginal likelihood function with respect to the vector θ is consistent and asymptotically normally distributed with asymptotic mean θ and covariance matrix given by the inverse of Godambe’s [20] sandwich information matrix Gθ [21]:(9)VCMLθ=Gθ−1=Hθ−1JθHθ−1,where Hθ and Jθ can be estimated straightforwardly at the CML estimate (θ^CML):H^θ^=−∑q=1Q∂2logLCML,qθ/∂θ∂θ′θ^CML and(10)J^θ^=∑q=1Q∂logLCML,qθ∂θ∂logLCML,qθ∂θ′θ^CML.


## 2.2. Matrix-Based Analytic Approximation Method

The computation of the MVNCD function can be applied for the multidimensional integral of the ICLV model above. In the previous ICLV model of Bhat, the MACML approach in 2011 [15] is adopted. In this section, a new matrix-based analytic approximation method proposed by Bhat in 2018 [16] can replace the MACML approach to evaluate the MVNCD function for acquiring better estimation results [22, 23].There are five methods for approximating the MVNCD function mentioned and summarized in the literature of Bhat in 2018. As per Bhat, the TVBS method is recommended as the evaluation approach for the MVNCD function with the highest accuracy, and the BME method has a significant speed advantage over these evaluation approaches. Besides, as the foundation of the series of analytic approximation methods, the ME method can also be included in the evaluation system in this study. Thus, these three typical methods (ME, BME, and TVBS) will be selected for elaboration below. Conceptually, the ME method is based on univariate conditioning, while BME and TVBS are based on a bivariate conditioning mechanism.Lete=e1,e2,…,eK be a multivariate normally distributed random vector with zero means, variances of 1, and a correlation matrix Σ=1ρ12⋯ρ1Kρ121⋯ρ2K⋮⋮⋱⋮ρ1Kρ2K⋯1 (K > 2), where K is the dimensionality of the MVNCD function. The K-dimensional MVNCD function can be expressed as follows: ΦKe1,e2,…,eK:Σ=Pe<V=Pe1<V1,e2<V2,…,eK<VK. To facilitate the understanding, assume that t represents e<V and tn represents e1<V1, and then, ΦKe:Σ=Pe<V=Pt=Pt1,t2,…,tK.The joint probability of the MVNCD function of the three methods (ME, BME, and TVBS) may be expressed as follows. All the joint probability may be written as the product of a marginal probability and conditional probabilities, and the normal distribution can be used to approximate the skew normal distribution for conditional probabilities. In the equations below, the superscriptn of tn is the number of truncations; nmax is the maximum number of truncations, and a lower nmax means a higher accuracy of the analytic approximation method. Let N be an intermediate parameter used to measure K: K is equal to N for the ME method, and K is discussed as being odd K=2N−1 or even K=2N for the BME (or TVBS) method. The expected values μk and variances σk2 of single-sided truncations in the following equations can be obtained by the properties of truncated multivariate normal distribution, which is not detailed here (refer to Section 2.1 of Bhat in 2018).


### 2.2.1. The ME MethodN≥2

(11)Pt=Pt1,t2,…,tN,=Pt1×∏n=1N−1Ptn+1|t1t2,…,tn,≈Pt1×∏n=1N−1Ptn+1n,≈ΦV1×∏n=1N−1ΦVn+1−μn+1σn+1.The method relies on a single-sided truncation of a multivariate normal distribution in which some variables are truncated, while others are not. After a single variable is truncated, the residual variables follow a skew normal distribution in the ME method.Ptn+1n is defined as the univariate conditional probabilities of tn+1 after n times truncations, which follows a univariate normal distribution (being used for approximating the skew normal distribution). Thus, the above Pt can be transformed into the product of recursive univariate normal cumulative distribution function (CDF) (represented by Φ.). In equation (11), K=N and nmax=N−1.For instance, whenK=6, the probability of multivariate normal distribution can be written as follows: Pt1,t2,t3,t4,t5,t6≈Pt1×Pt21×Pt32×Pt43×Pt54×Pt65=Pe1<V1×Pe2∗<V2|e1<V1×Pe3∗<V3|e2∗<V2×Pe4∗<V4|e3∗<V3×Pe5∗<V5|e4∗<V4×Pe6∗<V6|e5∗<V5, where ek∗ follows the skew normal distribution after k−1 times truncations. Thereinto, t2 was truncated once and Pt21=Pe2∗<V2|e1<V1≈ΦV2−μ2/σ2. z1=e1|e1<V1 is defined, and it can be deduced that the expected values μ2=ρ12Ez1 and variances σ22=ρ122Varz1+1−ρ122, where Ez1=μ1+σ1λ1=λ1, Varz1=σ121+λ1V1−λ12=1+λ1V1−λ12, and λ1=−ϕV1/ΦV1. Similar to the calculation principle of ΦV2−μ2/σ2, ΦVn+1−μn+1/σn+1 in equation (11) can be calculated by iteratively updating the expected value and covariance matrix. In the equation above, the maximum number of truncations nmax is 5.


### 2.2.2. The Bivariate ME (BME) MethodN≥2

WhenK=2N is even,(12)Pt=Pt1,t2,…,t2N,=Pt1,t2×∏n=1N−1Pt2n+1,t2n+2|t1,t2,…,t2n,≈Pt1,t2×∏n=1N−1Pt2n+12n,t2n+22n,≈Φ2V1,V2,ρ12×∏n=1N−1Φ2V2n+1−μ2n+1σ2n+1,V2n+2−μ2n+2σ2n+2,ρ2n+1,2n+2∗.WhenK=2N−1 is odd,(13)Pt=Pt1,t2,…,t2N−1,=Pt1,t2×∏n=1N−2Pt2n+1,t2n+2|t1,t2,…,t2n×Pt2N−1|t1,t2,…,t2N−2,≈Pt1,t2×∏n=1N−2Pt2n+12n,t2n+22n×Pt2N−12N−2,≈Φ2V1,V2,ρ12×∏n=1N−2Φ2V2n+1−μ2n+1σ2n+1,V2n+2−μ2n+2σ2n+2,ρ2n+1,2n+2∗×ΦV2N−1−μ2N−1σ2N−1.The BME method is in accordance with a bivariate truncation scheme, assuming that the marginal distribution of the third and fourth untruncated variables, given the first two are truncated, remains bivariate normal, and so on. In the equation above, similar to the ME method,Pt2n+12n,t2n+22n uses the bivariate normal distribution to approximate the bivariate skew normal distribution. The conditional probability Pt of this method for MVNCD function is a product of bivariate normal CDF Φ2. when the dimension is even and a product of Φ2. and Φ. when the dimension is odd. Here, ρ2n+1,2n+2∗ is the updated correlation coefficient between error terms after bivariate truncation, and the calculation principle of the expected values μk and variances σk2 is similar to the ME method. As the extension of the ME method, the BME method truncates two dimensions each time and nmax times in total (also does the TVBS method), where nmax=N−1 in equations (12) and (13).Similarly, whenK=6, the Pt using the BME method can be expressed as follows: Pt1,t2,t3,t4,t5,t6≈Pt1,t2×Pt32,t42×Pt54,t64. Thereinto, t3 and t4are truncated once n=1, and t5 and t6 are truncated twice n=2. The maximum number of truncations nmax is 2, which is smaller than the nmax of the ME method. This indicates that the accuracy of the bivariate truncation scheme is better than the univariate truncation scheme.


### 2.2.3. The Two-Variate Bivariate Screening (TVBS) Method

WhenK=2NN≥2 is even,(14)Pt=Pt1,t2,…,t2N,=Pt1,t2,t3,t4×∏n=1N−2Pt2n+1,t2n+2,t2n+3,t2n+4|t1,t2,…,t2nPt2n+1,t2n+2|t1,t2,…,t2n,≈Pt1,t2,t3,t4×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32n,t2n+42nPt2n+12n,t2n+22n,=Pt1,t2,t3×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32nPt2n+12n×Pt2N−12N−2,t2N2N−2Pt2N−12N−2,≈Φ3V1,V2,V3:0,Σ1:3,1:3×∏n=1N−2Φ3V2n+1,V2n+2,V2n+3:μ˜,Ω˜ΦV2n+1−μ2n+1/σ2n+1×Φ2V2N−1−μ2N−1/σ2N−1,V2N−μ2N/σ2N,ρ2N−1,2N∗ΦV2N−1−μ2N−1/σ2N−1.WhenK=2N−1N≥3 is odd,(15)Pt=Pt1,t2,…,t2N−1,=Pt1,t2,t3,t4×∏n=1N−3Pt2n+1,t2n+2,t2n+3,t2n+4|t1,t2,…,t2nPt2n+1,t2n+2,t1,t2,…,t2n×Pt2N−3,t2N−2,t2N−1|t1,t2,…,t2N−4Pt2N−3,t2N−2|t1,t2,…,t2N−4,≈Pt1,t2,t3,t4×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32n,t2n+42nPt2n+12n,t2n+22n×Pt2N−32N−4,t2N−22N−4,t2N−12N−4Pt2N−32N−4,t2N−22N−4,=Pt1,t2,t3×∏n=1N−3Pt2n+12n,t2n+22n,t2n+32nPt2n+12n,≈Φ3V1,V2,V3:0,Σ1:3,1:3×∏n=1N−2Φ3V2n+1,V2n+2,V2n+3:μ˜,Ω˜ΦV2n+1−μ2n+1/σ2n+1.Similarly, the TVBS method combines the bivariate truncation scheme, in which the bivariate distribution is assumed to be rectangle-screened normal (RSN) distribution [24]. The corresponding conditional CDF at each step can be expressed as a ratio of a four-variate normal CDF Φ4. and a two-variate normal CDF Φ2.. Since the calculation of four-variate normal CDF is time-consuming, the TVBS method adopts an approximation method by taking the trivariate CDF of the first three variates P123. Given the first three variables, the CDF of the skew normal distribution of the fourth variable is denoted as P4|123. Accordingly, a four-variate normal CDF is obtained: P1234=P123∗P4|123=P123∗P34|123/P3|12. Further, the MVNCD function is composed of the product of Φ3. (denotes trivariate normal CDF) and Φ3./Φ. when K is an odd number, and the product of Φ3., Φ3./Φ., and Φ2./Φ. when K is odd. In equations (14) and (15), μ˜ and Ω˜ represent the expected value and covariance matrix updated after bivariate truncation, respectively. Additionally, nmax=N−1 (K is even) and nmax=N−2 (K is odd) for the TVBS method.WhenK=6, the Pt using the TVBS method can be expressed as follows: Pt1,t2,t3,t4,t5,t6≈Pt1,t2,t3×Pt32,t42,t52/Pt32×Pt54,t64/Pt54, where N=3 and nmax=2. Here, t3, t4, and t5 conducted one truncation n=1 in Pt32,t42,t52/Pt32, while t5 and t6 conducted two truncations n=2 in Pt54,t64/Pt54.Furthermore, the matrix-based analytic approximation method applies an LDLT decomposition method for the correlation matrix followed by rank 1 (ME) or rank 2 (BME and method) updates of the factorization. The implementation is much easier to code and more computationally efficient.As we know from the above discussion, the multidimensional integral for the analytical approximation method of MNVCD can be completely converted into integrals with no more than three dimensions. At present, the CDF of bivariate normal (BVN) distribution and trivariate normal (TVN) distribution is generally realized by numerical methods based on the Legendre–Gauss quadrature integral approximation. This numerical integration method adopted in most existing packages is relatively old with low levels of accuracy and speed. In this study, referring to the study of Genz in 2004 [25], a new numerical integration method is used to replace the original method. For BVN probabilities, we adopted the algorithm proposed by Drezner and Wesolowsky [26] with a higher-order Ox4 Taylor expansion; for TVN probabilities, the second Plackett formula method studied by Gassmann [27] is applied (not be detailed here). Accordingly, after the above algorithms are applied along with matrix operations, the performance of both functions is improved in terms of accuracy and speed.


## 2.2.1. The ME MethodN≥2

(11)Pt=Pt1,t2,…,tN,=Pt1×∏n=1N−1Ptn+1|t1t2,…,tn,≈Pt1×∏n=1N−1Ptn+1n,≈ΦV1×∏n=1N−1ΦVn+1−μn+1σn+1.The method relies on a single-sided truncation of a multivariate normal distribution in which some variables are truncated, while others are not. After a single variable is truncated, the residual variables follow a skew normal distribution in the ME method.Ptn+1n is defined as the univariate conditional probabilities of tn+1 after n times truncations, which follows a univariate normal distribution (being used for approximating the skew normal distribution). Thus, the above Pt can be transformed into the product of recursive univariate normal cumulative distribution function (CDF) (represented by Φ.). In equation (11), K=N and nmax=N−1.For instance, whenK=6, the probability of multivariate normal distribution can be written as follows: Pt1,t2,t3,t4,t5,t6≈Pt1×Pt21×Pt32×Pt43×Pt54×Pt65=Pe1<V1×Pe2∗<V2|e1<V1×Pe3∗<V3|e2∗<V2×Pe4∗<V4|e3∗<V3×Pe5∗<V5|e4∗<V4×Pe6∗<V6|e5∗<V5, where ek∗ follows the skew normal distribution after k−1 times truncations. Thereinto, t2 was truncated once and Pt21=Pe2∗<V2|e1<V1≈ΦV2−μ2/σ2. z1=e1|e1<V1 is defined, and it can be deduced that the expected values μ2=ρ12Ez1 and variances σ22=ρ122Varz1+1−ρ122, where Ez1=μ1+σ1λ1=λ1, Varz1=σ121+λ1V1−λ12=1+λ1V1−λ12, and λ1=−ϕV1/ΦV1. Similar to the calculation principle of ΦV2−μ2/σ2, ΦVn+1−μn+1/σn+1 in equation (11) can be calculated by iteratively updating the expected value and covariance matrix. In the equation above, the maximum number of truncations nmax is 5.


## 2.2.2. The Bivariate ME (BME) MethodN≥2

WhenK=2N is even,(12)Pt=Pt1,t2,…,t2N,=Pt1,t2×∏n=1N−1Pt2n+1,t2n+2|t1,t2,…,t2n,≈Pt1,t2×∏n=1N−1Pt2n+12n,t2n+22n,≈Φ2V1,V2,ρ12×∏n=1N−1Φ2V2n+1−μ2n+1σ2n+1,V2n+2−μ2n+2σ2n+2,ρ2n+1,2n+2∗.WhenK=2N−1 is odd,(13)Pt=Pt1,t2,…,t2N−1,=Pt1,t2×∏n=1N−2Pt2n+1,t2n+2|t1,t2,…,t2n×Pt2N−1|t1,t2,…,t2N−2,≈Pt1,t2×∏n=1N−2Pt2n+12n,t2n+22n×Pt2N−12N−2,≈Φ2V1,V2,ρ12×∏n=1N−2Φ2V2n+1−μ2n+1σ2n+1,V2n+2−μ2n+2σ2n+2,ρ2n+1,2n+2∗×ΦV2N−1−μ2N−1σ2N−1.The BME method is in accordance with a bivariate truncation scheme, assuming that the marginal distribution of the third and fourth untruncated variables, given the first two are truncated, remains bivariate normal, and so on. In the equation above, similar to the ME method,Pt2n+12n,t2n+22n uses the bivariate normal distribution to approximate the bivariate skew normal distribution. The conditional probability Pt of this method for MVNCD function is a product of bivariate normal CDF Φ2. when the dimension is even and a product of Φ2. and Φ. when the dimension is odd. Here, ρ2n+1,2n+2∗ is the updated correlation coefficient between error terms after bivariate truncation, and the calculation principle of the expected values μk and variances σk2 is similar to the ME method. As the extension of the ME method, the BME method truncates two dimensions each time and nmax times in total (also does the TVBS method), where nmax=N−1 in equations (12) and (13).Similarly, whenK=6, the Pt using the BME method can be expressed as follows: Pt1,t2,t3,t4,t5,t6≈Pt1,t2×Pt32,t42×Pt54,t64. Thereinto, t3 and t4are truncated once n=1, and t5 and t6 are truncated twice n=2. The maximum number of truncations nmax is 2, which is smaller than the nmax of the ME method. This indicates that the accuracy of the bivariate truncation scheme is better than the univariate truncation scheme.


## 2.2.3. The Two-Variate Bivariate Screening (TVBS) Method

WhenK=2NN≥2 is even,(14)Pt=Pt1,t2,…,t2N,=Pt1,t2,t3,t4×∏n=1N−2Pt2n+1,t2n+2,t2n+3,t2n+4|t1,t2,…,t2nPt2n+1,t2n+2|t1,t2,…,t2n,≈Pt1,t2,t3,t4×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32n,t2n+42nPt2n+12n,t2n+22n,=Pt1,t2,t3×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32nPt2n+12n×Pt2N−12N−2,t2N2N−2Pt2N−12N−2,≈Φ3V1,V2,V3:0,Σ1:3,1:3×∏n=1N−2Φ3V2n+1,V2n+2,V2n+3:μ˜,Ω˜ΦV2n+1−μ2n+1/σ2n+1×Φ2V2N−1−μ2N−1/σ2N−1,V2N−μ2N/σ2N,ρ2N−1,2N∗ΦV2N−1−μ2N−1/σ2N−1.WhenK=2N−1N≥3 is odd,(15)Pt=Pt1,t2,…,t2N−1,=Pt1,t2,t3,t4×∏n=1N−3Pt2n+1,t2n+2,t2n+3,t2n+4|t1,t2,…,t2nPt2n+1,t2n+2,t1,t2,…,t2n×Pt2N−3,t2N−2,t2N−1|t1,t2,…,t2N−4Pt2N−3,t2N−2|t1,t2,…,t2N−4,≈Pt1,t2,t3,t4×∏n=1N−2Pt2n+12n,t2n+22n,t2n+32n,t2n+42nPt2n+12n,t2n+22n×Pt2N−32N−4,t2N−22N−4,t2N−12N−4Pt2N−32N−4,t2N−22N−4,=Pt1,t2,t3×∏n=1N−3Pt2n+12n,t2n+22n,t2n+32nPt2n+12n,≈Φ3V1,V2,V3:0,Σ1:3,1:3×∏n=1N−2Φ3V2n+1,V2n+2,V2n+3:μ˜,Ω˜ΦV2n+1−μ2n+1/σ2n+1.Similarly, the TVBS method combines the bivariate truncation scheme, in which the bivariate distribution is assumed to be rectangle-screened normal (RSN) distribution [24]. The corresponding conditional CDF at each step can be expressed as a ratio of a four-variate normal CDF Φ4. and a two-variate normal CDF Φ2.. Since the calculation of four-variate normal CDF is time-consuming, the TVBS method adopts an approximation method by taking the trivariate CDF of the first three variates P123. Given the first three variables, the CDF of the skew normal distribution of the fourth variable is denoted as P4|123. Accordingly, a four-variate normal CDF is obtained: P1234=P123∗P4|123=P123∗P34|123/P3|12. Further, the MVNCD function is composed of the product of Φ3. (denotes trivariate normal CDF) and Φ3./Φ. when K is an odd number, and the product of Φ3., Φ3./Φ., and Φ2./Φ. when K is odd. In equations (14) and (15), μ˜ and Ω˜ represent the expected value and covariance matrix updated after bivariate truncation, respectively. Additionally, nmax=N−1 (K is even) and nmax=N−2 (K is odd) for the TVBS method.WhenK=6, the Pt using the TVBS method can be expressed as follows: Pt1,t2,t3,t4,t5,t6≈Pt1,t2,t3×Pt32,t42,t52/Pt32×Pt54,t64/Pt54, where N=3 and nmax=2. Here, t3, t4, and t5 conducted one truncation n=1 in Pt32,t42,t52/Pt32, while t5 and t6 conducted two truncations n=2 in Pt54,t64/Pt54.Furthermore, the matrix-based analytic approximation method applies an LDLT decomposition method for the correlation matrix followed by rank 1 (ME) or rank 2 (BME and method) updates of the factorization. The implementation is much easier to code and more computationally efficient.As we know from the above discussion, the multidimensional integral for the analytical approximation method of MNVCD can be completely converted into integrals with no more than three dimensions. At present, the CDF of bivariate normal (BVN) distribution and trivariate normal (TVN) distribution is generally realized by numerical methods based on the Legendre–Gauss quadrature integral approximation. This numerical integration method adopted in most existing packages is relatively old with low levels of accuracy and speed. In this study, referring to the study of Genz in 2004 [25], a new numerical integration method is used to replace the original method. For BVN probabilities, we adopted the algorithm proposed by Drezner and Wesolowsky [26] with a higher-order Ox4 Taylor expansion; for TVN probabilities, the second Plackett formula method studied by Gassmann [27] is applied (not be detailed here). Accordingly, after the above algorithms are applied along with matrix operations, the performance of both functions is improved in terms of accuracy and speed.


## 2.3. Automatic Differentiation (AD)

Automatic differentiation (AD) is a chain rule-based technique for evaluating the derivatives with respect to the input variables of functions defined by a computer program [28]. There are two different AD strategies: the forward mode and the reverse mode. The advantage of the reverse mode over the forward mode will be highlighted when the number of parameters is large, which is more suitable for the MLE procedure.For a simple example, consider a binary probit model:U=a+bx+ε,y=IU>0, and ε∼N0,1. The logarithmic likelihood function of the model can be written as LLθ=LLa,b=lnΦ2y−1a+bx. In the reverse mode of AD, the derivative is computed with respect to each sub-expression recursively in the chain rule:(16)∂LL∂θ=∂LL∂v1∂v1∂θ=∂y∂v2∂v2∂v1∂v1∂θ=∂y∂v3∂v3∂v2∂v2∂v1∂v1∂θ=….The derivative ofLL with respect to vi is defined as v¯i=∂LL/∂vi, and the reverse accumulation with computation graph is illustrated in Figure 2. In this figure, the partial derivatives are estimated as follows: ∂LL/∂a=a¯=2y−1∗φ2y−1∗a+bx/Φ2y−1∗a+bx and ∂LL/∂b=b¯=x∗2y−1∗φ2y−1∗a+bx/Φ2y−1∗a+bx. Thereinto, φ. has a closed form, and Φ. can be approximated by the numerical integration method.Figure 2
Example of reverse-mode AD in the computation graph.TensorFlow (https://www.tensorflow.org/), an open-source framework for large-scale distributed numerical computation, is widely applied in the areas of machine learning and deep neural network during early explorations. The main feature of the TensorFlow framework is based on a data flow graph and uses tensor as the core component to denote data. The framework builds a gradient graph corresponding to the original computation graph by backward propagation and chain rule of derivative values, which is more mature and flexible than other programming modes in the application of automatic gradient. Due to the advantages of TensorFlow, we applied it for developing econometric models in this study.AD has obvious advantages when it is applied to the MLE of parameters in econometric models. Using this technique, modelers do not need to derive and code the complex analytic gradient expression of the log-likelihood function, which greatly reduces the workload. When the CML approach is used to estimate model parameters, the Hessian matrix and Jacobian matrix of estimated parameters need to be calculated. The errors and complexity of the numerical differentiation method increase when second-order derivatives are computed, and the speed is low in calculating partial derivatives of a function with respect to many inputs. Therefore, the use of automatic differentiation instead of numerical differentiation can aid in attaining more accurate estimation results with less computational time.


## 3. Simulation Experiments

To evaluate the ability of the proposed ICLV model to recover parameters from finite samples, it is necessary to conduct Monte Carlo simulation experiments. In the study of Bhat [10], four ordered indicator variables and three alternatives are concerned. In this section, referring to Bhat’s experiment and further increasing the number of alternatives, we applied the LDLT matrix-based analytic approximation method for a higher integral dimension in simulation experiments.Here, we design simulation experiments and discuss the identifiability of the model system based on two assumptions (excluding those mentioned in the last section): (1) the number of indicatorsL will not be less than the number of latent variables G (i.e., each latent variable at least has one indicator, which only loads on that variable but does not load on any other latent variable); (2) there is no common variable between each row of x matrix and the w vector. To explain the model more clearly and concisely, the following ICLV model is designed as an example. Consider that the ICLV model is composed of four latent variables zl∗, four ordered indicator variables yg, and three alternatives’ utilities Ui, where the structural equation model of the latent variables in equation (1) can be represented as follows:(17)z1∗z2∗z3∗z4∗=α10.00.0α3α20.00.00.0α40.0α50.00.0α60.00.00.00.00.0α7w1w2w3w4w5+η1η2η3η4.Assume that each ordered variable has only one thresholdψg (i.e., ψg=ψg,2), and the four-variate ordered probit model for the latent variable measurement model of equation (2) is as follows:(18)y1y2y3y4=δ1δ2δ3δ4+d10.00.0d20.00.00.00.00.00.00.00.0d30.00.0d4z1∗z2∗z3∗z4∗+ξ1ξ2ξ3ξ4Then, the choice model with three alternatives of equation (3) can be specified as follows:(19)U1U2U3=001001tt1tc1tt2tc2tt3tc3β2β3βttβtc+0.00.0λ10.00.0λ20.00.0λ30.00.0λ4z1∗z2∗z3∗z4∗+ε1ε2ε3,where tti and tci are the exogenous explanatory variables in the choice model.In the above model, the correlation matrixΣ of the multivariate ordered probit model is fixed as an identity matrix, and the correlation matrix Γ of latent variables and the covariance matrix Λ of the choice model can be defined as follows:(20)Σ=1000010000100001,Γ=1r12r13r14r121r23r24r13r231r34r14r24r341,Λ=σ11σ12σ13σ12σ22σ23σ13σ23σ33.There areL∗L−1/2 parameters in Γ and I∗I−1/2−1 parameters in Λ that need to be estimated. The covariance matrix Ω of the whole model after considering z∗ can be expressed as follows:(21)Ω=dΓd′+ΣdΓλ′λΓd′λΓλ′+Λ=d12+1d1d2r12d1d3r130d1λ1+d1λ3r13d1λ2r12+d1λ4r14λ12+λ32+2λ1λ3r13+σ22d2d4r24d3d4r140d2λ1r12+d2λ3r23d2λ2+d2λ4r24λ1λ2r12+λ1λ4r14+λ2λ3r23+λ3λ4r34+σ23d32+1d3d4r240d3λ1r13+d3λ3d3λ2r23+d3λ4r34λ22+λ42+2λ2λ4r24+σ33d2d3r23d3d4r340d4λ1r14+d4λ3r34d4λ2r24+d4λ4d22+1⋅d42+1σ11σ12σ13.Applying the CML approach of in equations (7) and (8), the number of identifiable parameters in the Ω is 16 (6 + 8+2) in the above sample. Thereinto, G∗G−1/2 parameters of the matrix dΓd′+Σ can be identified (the covariance matrix needs to normalize the diagonal elements to 1 in advance), G∗I−1 parameters of dΓλ′ can be identified where I−1 parameters can be identified in each row at most, and I∗I−1/2−1 parameters can be identified in the matrix λΓλ′+Λ as mentioned in the section of methodology.As shown in Figure3, we conduct the simulation experiments of two models: one contains three alternatives of the choice model and the other contains four alternatives, with the same settings in the latent variable structural equation model and measurement equation model in which r12, r14, r23, and r34 are fixed at the value of zero. For model 1 of three alternatives, the simulation experiment is undertaken 50 times for a sample size of 10,000 with different realizations of YU vector to generate 50 different datasets. The estimator is then applied to each dataset to estimate data specific values for the 31 parameters, and the numerical integration method for calculating the CDF of TVN distribution is used for the MVNCD function; for model 2 of four alternatives, we perform the above data generation progress once with a large sample size of 100,000 and then carry out the Monte Carlo simulation experiments on the same dataset using ME, BME, and TVBS methods, respectively. The analytic approximation method of four-variate normal CDF Pt1,t2,t3,t4 can be written as follows: Pt1×Pt21×Pt32×Pt43 for the ME method, Pt1,t2×Pt32,t42 for the BME method, and Pt1,t2,t3×Pt32,t42/Pt32 for the TVBS method. The 37 parameters can be estimated in the ICLV model.Figure 3
Parameter settings of the ICLV model (for three and four alternatives).


## 4. Model Estimation Results

### 4.1. The MVNCD Function Evaluation Results

To evaluate the accuracy of computing the MVNCD function for different dimensions, we generate 1000 groups of random numbers, which consist of four sets of 250 MVNCD evaluations each, corresponding to (a) low correlations, high MVNCD values, (b) low correlations, low MVNCD values, (c) high correlations, high MVNCD values, and (d) high correlations, low MVNCD values. On the one hand, to generate the datasets of low and high MVNCD values, the upper integration limits can be drawn from0,+K for the cases of high MVNCD values, while it can be drawn from −K/2,+K for the cases of low MVNCD values. On the other hand, to generate the datasets of low correlations and high correlations, the positive-definite covariance matrices C can be calculated by C=RR′+δ×diagru, where R represents a matrix of K×K random univariate standard normal variates and diagru is a diagonal matrix with the K×1 vector ru of standard uniform random variates on the diagonal. δ is a scalar to determine the relative magnitude of the diagonal elements relative to the non-diagonal elements. We can use the two different values of δ to distinguish the matrix with low correlations δ=10 and high correlations δ=0.The MVNCD evaluation results from the different methods are presented in Tables1 and 2. Here, we use the Gaussian quadrature integral based on 400 nodes as the exact value when K ≤ 3 and use the value calculated by the CDFMVNe function of the Gauss (up to an accuracy of 1e − 6) as the exact value when K > 3. Meanwhile, the mean absolute error (MAE) and the mean absolute percentage error (MAPE) are chosen to evaluate the accuracy level of the MVNCD evaluation results. The “%MAE > 0.005” represents the percentage in which the MAE was over 0.005, and “%MAPE > 2” represents the percentage in which the MAPE was over 2.Table 1
BVN and TVN evaluation results for the different methods.
DimMeasureNumerical integration methodsAnalytic approximation methods (TensorFlow)GaussTFMEBMETVBSK = 2 (BVN)MAE4.8∗e − 101.9∗e − 150.0037——%MAE > 0.0050.00.018.8——MAPE8.6∗e − 81.3∗e − 121.28——%MAPE > 20.00.015.2——Time (s)0.0010.0020.003——K = 3 (TVN)MAE8.8∗e − 106.1∗e − 110.00270.0018—%MAE > 0.0050.00.018.112.4—MAPE1.9∗e − 72.9∗e − 81.341.00—%MAPE > 20.00.016.711.0—Time (s)0.0050.0050.0060.003—Note. “—“ indicates that the corresponding analytic approximation methods cannot calculate the values of BVN or TVN.Table 2
MVNCD evaluation results for the different analytic approximation methods.
DimMeasureMEBMETVBSGaussTFGaussTFGaussTFK = 4MAE0.002260.002260.001520.001520.001170.00117MAPE0.590.590.400.400.260.26%MAE > 0.00515.215.29.89.85.85.8%MAPE > 218.018.012.812.88.48.4Time (s)0.0520.0100.0330.0040.0390.007K = 5MAE0.001710.001710.001420.001420.001080.00108MAPE0.540.540.440.440.330.33%MAE > 0.0059.29.27.07.04.84.8%MAPE > 218.818.815.415.413.013.0Time (s)0.0770.0160.0600.0070.0650.011K = 7MAE0.001200.001200.001040.001040.000930.00093MAPE0.450.450.380.380.330.33%MAE > 0.0055.95.94.74.73.63.6%MAPE > 221.521.520.220.218.318.3Time (s)0.1410.0310.0970.0140.1050.021K = 10MAE0.000640.000640.000570.000570.000540.00054MAPE0.320.320.290.290.290.29%MAE > 0.0051.61.61.21.21.11.1%MAPE > 225.025.023.123.122.622.6Time (s)0.2580.0610.1560.0230.1820.033K = 15MAE0.000410.000410.000400.000400.000370.00037MAPE0.290.290.280.280.250.25%MAE > 0.0050.20.20.20.20.10.1%MAPE > 234.134.132.632.632.032.0Time (s)0.5580.1460.3170.0580.3520.076From Table1, the following findings can be obtained. Compared with the calculated results by the cdfBvn and cdfTvn functions in the GAUSS (https://www.aptech.com/wp-content/uploads/2016/05/GAUSS_16_PDF_UG.pdf), the accuracy levels (shown as the MAE) of BVN and TVN probability (based on the improvement in Section 2.3) calculated with the assistance of TensorFlow have been improved from e−10 to e−15 and from e−10 to e−11, respectively, with a slight loss in time. When the integration dimension K is 2 or 3, as shown in Table 1, the accuracy level of the MAE for the ME and BME methods is e−3, which is far lower than that of bivariate and trivariate normal distribution function calculated by numerical integration. Similar conclusion can be obtained by choosing the MAPE as an index. Therefore, numerical integration should be adopted to calculate BVN or TVN values.WhenK ≥ 4, we employed the LDLT matrix-based analytic approximation method for the evaluation on Gauss and TensorFlow platforms. From Table 2, when the number of dimensions “K” increases, the performance of the three methods can be indicated by “MAE,” “MAPE,” “%MAE > 0.005,” and “%MAPE > 2.” The TVBS method is taken as an example, and there is a variation as the dimension “K” goes from 4 to 15, with “MAE” reducing from 0.00117 to 0.00037, “%MAE > 0.005” reducing from 5.8% to 0.1%, and “%MAPE > 2” increasing from 8.4% to 32.0%. As shown in Figure 4(a), “%MAE > 0.005” gradually decreases and “%MAPE > 2” gradually increases with the increase in dimension of all three methods. To explain the intrinsic reason, we checked the data and found that the overall average value of 1000 MVNCD values calculated by the method of generating random number would decrease with the increase in dimension. Due to the lower MVNCD values in the higher dimension, the relative error (MAPE) may become larger with the higher dimension. Meanwhile, as expected, with the increase in K, the time required by the analytic approximation gradually increased for all the methods. In addition, when we compare the three analytical approximation methods with each other at the same dimension, it can be found that the time consumed by the BME method is the shortest. In general, the time consumed by BME is less than half of that by the ME method and that by the TVBS method falls in the middle. Compared with the ME method, the error of the BME and TVBS methods is relatively low, as both of them are based on the bivariate truncation scheme. Overall, the TVBS is the best in terms of accuracy in evaluating individual MVNCD functions. These findings are consistent with Bhat’s study (2018) and mutually corroborate.Figure 4
Comparison of accuracy and time with different dimensions. (a) “%MAE > 0.005” and “%MAPE > 2” (b) time.
(a)(b)Furthermore, there are some findings by comparing the calculation error and calculation time of the different analytic approximation methods (i.e., ME, BME, and TVBS) in Gauss and TensorFlow. In terms of calculation accuracy, it can be seen in Table2 that “MAE,” “MAPE,” “%MAE > 0.005,” and “%MAPE > 2” of the three analytic approximation methods are basically the same. Therefore, although the new method improves accuracy levels of BVN and TVN probability in TensorFlow, its contribution to accuracy improvement is limited for high-dimensional analytic approximation; in terms of calculation time, it can be seen in Table 2 and Figure 4(b) that under the same dimension and algorithm, the calculation speed of TensorFlow platform is much higher than that of Gauss platform, and the average running time of Gauss is about 5-6 times that of TensorFlow. As shown in Table 1, there is no significant difference in the calculation time of BVN and TVN values on the two platforms. Thus, the declines in calculation time are mainly due to the obvious time advantage of TensorFlow with tensor data structure and graph computation technique.


### 4.2. Simulation Results for ICLV Model

Following Section3, the two ICLV models are designed to evaluate the ability for recovering the model parameters. For the three-alternative ICLV model with the maximum integral dimension being 3, the performance of the improved trivariate normal integral approach in estimating the parameters and the standard errors of the ICLV model is evaluated in Table 3. The evaluation system can be adopted as follows:(1)
For the parameter estimates, the mean estimate for each model parameter across the 50 datasets can be obtained. Then, the absolute bias of the estimator is computed asAbs.Bias=estimate value−true value and the absolute percentage bias as APB=estimate value−true value/true value×100.(2)
For the standard error estimates, the standard deviation of each estimated parameter across the 50 datasets can be computed, and this is labeled as the finite sample standard deviation or FSSD (essentially, this is the empirical standard error); also, we can compute the mean standard deviation for each model parameter across the 50 datasets and denote this as the asymptotic standard error or ASE (the standard error here is computed using the Godambe (sandwich) matrix at convergence for a theoretical approximation to the FSSD). Next, to evaluate the accuracy of the ASE formula for the finite sample size used, the relative efficiency of the estimator is computed asrelative efficiency=ASE/FSSD. As per Bhat [10], the closeness between ASE and FSEE may be examined as captured by the 0.75–1.25 range for the relative efficiency value.Table 3
Simulation results for the three-alternative ICLV model.
Parm.True valueParameter estimatesStandard error estimatesMean est.Abs. biasAPB (%)ASEFSSDR. effα10.50.53000.03005.99240.09610.08071.1908α20.60.64110.04116.85060.11770.10031.1730α30.50.50220.00220.44320.20370.14931.3648α40.60.60170.00170.27750.24510.18451.3288α50.30.31250.01254.16300.04080.04660.8740α60.30.30270.00270.91190.13360.11181.1950α7−0.4−0.40610.00611.51860.17730.14551.2182s11——————r120——————r130.60.68550.085514.24850.25670.25161.0206r140——————s21——————r230——————r240.60.62300.02303.84110.48000.39511.2149s31——————r340——————s41——————δ1−1−1.00130.00130.13120.03660.03341.0974δ2−1−1.02440.02442.44180.06800.06381.0664δ3−1−1.00310.00310.31170.04410.04530.9720δ4−1−1.06650.06656.65130.16290.16091.0124d10.30.28800.01203.99110.05080.04411.1524d20.40.44590.045911.48470.15500.13551.1441d30.50.49240.00761.52960.07540.08320.9064d40.60.71570.115719.27580.30460.26191.1633ψ11.51.49710.00300.19670.02750.02651.0375ψ21.51.53540.03542.35780.08730.08231.0601ψ31.51.49490.00510.34310.05040.05490.9191ψ41.51.61410.11417.60540.23460.21141.1101β20.50.48630.01372.74090.09540.08881.0741β3−1−0.99720.00280.27840.13600.15300.8893βtt−1−0.99800.00200.19710.04190.04340.9667βtc−0.8−0.79780.00220.27720.03380.03430.9852λ10.50.48070.01933.86940.08850.08111.0907λ20.20.21290.01296.44770.06610.05171.2771λ30.50.49290.00711.41460.07180.07230.9929λ40.20.21120.01125.61360.06220.05431.1441σ111——————σ120——————σ130——————σ221——————σ230.60.60650.00651.07940.13740.11961.1483σ331.361.37530.01531.12860.27420.25711.0668Overall mean value0.02363.79400.13080.11691.1190Note: “--” indicates that the corresponding parameter is fixed at the true value.The following aspects can be observed in Table3 (the experiment is repeated 50 times with a sample size of 10,000). Firstly, the overall mean value across parameters of APB is 3.7940%, with the individual parameter APB ranging from 0.1% to 19%, which are close to those in Bhat. Besides, the simulation errors are also similar and therefore considered reasonable and acceptable; secondly, the mean APB values for the d vector (9.0%), λ vector (4.3%), and r in Γ matrix (9.0%) are relatively higher than the overall mean APB. This indicates that the coefficients on the latent variable vector z∗ in the measurement equation and the choice model are somewhat more difficult to recover than other parameters. It is expected, since these elements enter into the covariance matrix Ω in a nonlinear fashion as shown in equation (6), and Ω itself enters into the composite likelihood function (equation (8)) in a complex manner. Lastly, as for simulation error estimates, the overall mean value of relative efficiency (RE) across the parameters is 1.1190 (and the RE of almost all parameters is between 0.75 and 1.25), which suggests that the FSSE and ASE values are close to each other and the asymptotic formula is performing well in estimating the finite sample standard errors.The simulation results (with a sample size of 100,000) of the four-alternative ICLV model estimated by the three analytical approximation methods, namely the ME, BME, and TVBS, are presented in Tables4 –6, respectively. In these tables, Abs.Bias=estimate value−true value and APB=estimate value−true value/truevalue×100. As shown in the experimental results, it is ensured that all the identifiable parameters can be significantly estimated and seemingly consistent with their true values. We can see that the values of absolute percentage bias (APB) fall in the ranges of 0–15% in the three tables, of which the APB of four parameters exceeds 10% (i.e., α3, α4, d2, and σ34). It is a similar finding showing that the parameter estimators related to z∗ are less consistent. Additionally, in these three tables (also shown in Figure 5), it can be observed that among the three methods (i.e., ME, BME, and TVBS), the overall mean value of absolute bias is 0.0266, 0.0262, and 0.0259, and the average value of APB is 4.269%, 4.240%, and 4.204%, respectively. It is noticed that the values of these two indicators slightly decrease. On the other hand, the standard errors of all methods are about 0.037. It shows that the standard errors of the four-alternative ICLV model seem reasonable and acceptable. Overall, we have validated the ability for recovering parameters of these three LDLT-based analytic approximation methods for the ICLV model, and we can conclude that TVBS is superior to BME and BME is superior to ME. These findings are consistent with the evaluation results with respect to individual MVNCD functions. Although model parameters can be estimated significantly with a sample of large size, there are issues of empirical identification, namely whether the data can support the model specification. If the empirical data cannot allow for significantly estimating all the theoretically identifiable parameters, some parameters may be fixed and only parameters of particular interest need to be estimated.Table 4
Simulation results for the four-alternative ICLV model with the ME method.
Parm.True valueEstimateAbs. biasAPB (%)Std. err.t-testP valueα10.50.48210.01793.58060.017527.4871≤0.001α20.60.57340.02664.43450.021027.2569≤0.001α30.50.56230.062312.45820.054510.3235≤0.001α40.60.66590.065910.98580.064510.3299≤0.001α50.30.30390.00391.28330.010329.4528≤0.001α60.30.30590.00591.96930.030610.0031≤0.001α7−0.4−0.40990.00992.48170.0408−10.0422≤0.001s11——————r120——————r130.60.59820.00180.30820.046013.0071≤0.001r140——————s21——————r230——————r240.60.64750.04757.92420.10076.4324≤0.001s31——————r340——————s41——————δ1−1−1.00680.00680.68180.0115−87.6878≤0.001δ2−1−0.96710.03293.29330.0161−59.9248≤0.001δ3−1−1.00810.00810.80700.0131−76.7456≤0.001δ4−1−0.98750.01251.25310.0373−26.4410≤0.001d10.30.31470.01474.91600.012525.2686≤0.001d20.40.35080.049212.30200.03819.2070≤0.001d30.50.49960.00040.08280.021023.7820≤0.001d40.60.59100.00901.50230.07937.4520≤0.001ψ11.51.51200.01200.80040.0081186.2783≤0.001ψ21.51.47920.02081.38350.018878.6552≤0.001ψ31.51.49990.00010.00340.0145103.2049≤0.001ψ41.51.51580.01581.05530.053828.1516≤0.001β20.50.45470.04539.05600.030414.9540≤0.001β3−1−1.04160.04164.16160.0476−21.8600≤0.001β4−0.8−0.87280.07289.10370.0445−19.5951≤0.001βtt−1−1.01150.01151.14900.0158−63.8326≤0.001βtc−0.8−0.81070.01071.34240.0128−63.5246≤0.001λ10.50.53800.03807.60640.023422.9537≤0.001λ20.20.18580.01427.07700.018610.0106≤0.001λ30.20.18770.01236.14750.018710.0234≤0.001λ40.50.48910.01092.17560.023221.0530≤0.001λ50.50.51030.01032.06320.021224.0984≤0.001λ60.20.19970.00030.15650.02049.7718≤0.001σ111——————σ120——————σ130——————σ140——————σ221——————σ230.60.59560.00440.73650.052611.3276≤0.001σ240.60.64070.04076.77720.046013.9313≤0.001σ331.361.40080.04082.99990.107413.0425≤0.001σ340.60.69390.093915.64220.07798.9014≤0.001σ441.361.47230.11238.25420.100714.6142≤0.001Overall mean value0.02664.26910.0371Note: “—” indicates that the corresponding parameter is fixed at the true value.Table 5
Simulation results for the four-alternative ICLV model with the BME method.
Parm.True valueEstimateAbs. biasAPB (%)Std. err.t-testP valueα10.50.48210.01793.57600.017527.4890≤0.001α20.60.57340.02664.43070.021027.2577≤0.001α30.50.56220.062212.43020.054010.4099≤0.001α40.60.66580.065810.95830.063910.4166≤0.001α50.30.30440.00441.45370.010329.4308≤0.001α60.30.30600.00602.01600.03069.9954≤0.001α7−0.4−0.41010.01012.52800.0409−10.0344≤0.001s11——————r120——————r130.60.59900.00100.16620.046112.9991≤0.001r140——————s21——————r230——————r240.60.64770.04777.95000.10046.4540≤0.001s31——————r340——————s41——————δ1−1−1.00680.00680.68160.0115−87.6890≤0.001δ2−1−0.96710.03293.29030.0161−60.1612≤0.001δ3−1−1.00760.00760.76470.0131−76.8664≤0.001δ4−1−0.98730.01271.26820.0373−26.4510≤0.001d10.30.31470.01474.91100.012525.2700≤0.001d20.40.35090.049112.27750.03789.2838≤0.001d30.50.49850.00150.29240.021023.7831≤0.001d40.60.59060.00941.56250.07937.4485≤0.001ψ11.51.51200.01200.80000.0081186.2896≤0.001ψ21.51.47930.02071.38040.018779.1925≤0.001ψ31.51.49930.00070.04560.0145103.4676≤0.001ψ41.51.51560.01561.03930.053828.1629≤0.001β20.50.45460.04549.07900.030414.9559≤0.001β3−1−1.04170.04174.17480.0476−21.8780≤0.001β4−0.8−0.87190.07198.99350.0445−19.5883≤0.001βtt−1−1.01120.01121.12150.0158−63.9360≤0.001βtc−0.8−0.81060.01051.31870.0127−63.6256≤0.001λ10.50.53790.03797.57200.023422.9583≤0.001λ20.20.18590.01417.04550.018410.0864≤0.001λ30.20.18770.01236.16000.018610.1024≤0.001λ40.50.48820.01182.36540.023221.0692≤0.001λ50.50.50940.00941.87680.021124.1120≤0.001λ60.20.19950.00050.26650.02049.7651≤0.001σ111——————σ120——————σ130——————σ140——————σ221——————σ230.60.59520.00480.79270.052511.3353≤0.001σ240.60.63770.03776.28870.045713.9456≤0.001σ331.361.40020.04022.95340.107313.0519≤0.001σ340.60.69300.093015.49500.07768.9351≤0.001σ441.361.46280.10287.55510.100114.6143≤0.001Overall mean value0.02624.24000.0370Note: “—” indicates that the corresponding parameter is fixed at the true value.Table 6
Simulation results for the four-alternative ICLV model with the TVBS method.
Parm.True valueEstimateAbs. biasAPB (%)Std. err.t-testP valueα10.50.48200.01803.60060.017527.4956≤0.001α20.60.57330.02674.45500.021027.2643≤0.001α30.50.56210.062112.42480.054110.3956≤0.001α40.60.66570.065710.95320.064010.4023≤0.001α50.30.30420.00421.39130.010429.3686≤0.001α60.30.30580.00581.93330.030510.0157≤0.001α7−0.4−0.40980.00982.44530.0408−10.0548≤0.001s11——————r120——————r130.60.59860.00140.22630.046013.0038≤0.001r140——————s21——————r230——————r240.60.64720.04727.86280.10026.4570≤0.001s31——————r340——————s41——————δ1−1−1.00680.00680.68440.0115−87.6833≤0.001δ2−1−0.96710.03293.28980.0161−60.1169≤0.001δ3−1−1.00780.00780.78070.0131−76.7399≤0.001δ4−1−0.98760.01241.24080.0373−26.4477≤0.001d10.30.31480.01484.94070.012525.2747≤0.001d20.40.35090.049112.27280.03799.2709≤0.001d30.50.49890.00110.21540.021023.7248≤0.001d40.60.59130.00871.45550.07937.4595≤0.001ψ11.51.51200.01200.80260.0081186.2650≤0.001ψ21.51.47930.02071.37990.018779.0895≤0.001ψ31.51.49960.00040.02990.0145103.1776≤0.001ψ41.51.51600.01601.06770.053828.1597≤0.001β20.50.45440.04569.11440.030414.9438≤0.001β3−1−1.04250.04254.24790.0476−21.8930≤0.001β4−0.8−0.87280.07289.10250.0445−19.6065≤0.001βtt−1−1.01160.01161.16070.0158−63.9212≤0.001βtc−0.8−0.81090.01091.35670.0127−63.6051≤0.001λ10.50.53820.03827.64040.023422.9612≤0.001λ20.20.18600.01406.98350.018510.0749≤0.001λ30.20.18780.01226.11800.018610.0900≤0.001λ40.50.48870.01132.25980.023221.0338≤0.001λ50.50.50990.00991.98140.021224.0589≤0.001λ60.20.19970.00030.15100.02049.7806≤0.001σ111——————σ120——————σ130——————σ140——————σ221——————σ230.60.59230.00771.28320.052411.3058≤0.001σ240.60.63820.03826.37020.045813.9236≤0.001σ331.361.38950.02952.16590.107412.9430≤0.001σ340.60.68750.087514.58380.07758.8672≤0.001σ441.361.46290.10297.56430.100114.6142≤0.001Overall mean value0.02594.20370.0370Note. “—” indicates that the corresponding parameter is fixed at the true value.Figure 5
Comparison of three methods for the four-alternative ICLV model: (a) overall mean value of Abs. bias and (b) overall mean value of APB (%).
(a)(b)Furthermore, more attempts were made to discuss the operational efficiency of models using the automatic differential technique. The convergence time of the three-alternative ICLV model in Table3 has a median value of about 8 minutes for the case of 10,000 observations, which is much shorter than that of the literature of Bhat in 2014 (about one hour for the case of 1000 observations). It can also be inferred that the computational time of the improved four-alternative ICLV model based on TensorFlow is reduced. This can be attributed to two factors: one is that the LDLT-based analytic approximation method (instead of the MACML approach) can shorten the calculation time of the ICLV model, as mentioned in Bhat. The improvement on bivariate normal CDF and trivariate normal CDF also facilitates this process; the other is based on the time advantage of TensorFlow in calculating tensor and multidimensional arrays. In addition, in the absence of coding the analytic gradient artificially, the convergence speed of the model can be improved using automatic differentiation to obtain the gradient of the likelihood function. For the four-alternative ICLV model with a sample size of 100,000, the convergence times of three analytical approximation methods are 3 h7 min (for ME), 1 h14 min (for BME), and 2 h6 min (for TVBS), respectively. Among them, the BME algorithm has obvious advantages in terms of time consumption. Hence, if the accuracy level can be lowered under some situations (e.g., in a variable screening process), the BME algorithm can also be chosen to save computational time. In this section, it is not necessary to evaluate the ICLV model of Bhat for comparison, because Bhat has explicitly shown that the LDLT matrix-based analytic approximation method runs much faster than other methods (such as the SSJ method) in evaluating MVNCD functions and estimating MNP model. It is worth noting that all simulation experiments are conducted on a desktop computer with an Intel Core™ i5-4460@3.20 GHz Processor and 16 GB of RAM.


## 4.1. The MVNCD Function Evaluation Results

To evaluate the accuracy of computing the MVNCD function for different dimensions, we generate 1000 groups of random numbers, which consist of four sets of 250 MVNCD evaluations each, corresponding to (a) low correlations, high MVNCD values, (b) low correlations, low MVNCD values, (c) high correlations, high MVNCD values, and (d) high correlations, low MVNCD values. On the one hand, to generate the datasets of low and high MVNCD values, the upper integration limits can be drawn from0,+K for the cases of high MVNCD values, while it can be drawn from −K/2,+K for the cases of low MVNCD values. On the other hand, to generate the datasets of low correlations and high correlations, the positive-definite covariance matrices C can be calculated by C=RR′+δ×diagru, where R represents a matrix of K×K random univariate standard normal variates and diagru is a diagonal matrix with the K×1 vector ru of standard uniform random variates on the diagonal. δ is a scalar to determine the relative magnitude of the diagonal elements relative to the non-diagonal elements. We can use the two different values of δ to distinguish the matrix with low correlations δ=10 and high correlations δ=0.The MVNCD evaluation results from the different methods are presented in Tables1 and 2. Here, we use the Gaussian quadrature integral based on 400 nodes as the exact value when K ≤ 3 and use the value calculated by the CDFMVNe function of the Gauss (up to an accuracy of 1e − 6) as the exact value when K > 3. Meanwhile, the mean absolute error (MAE) and the mean absolute percentage error (MAPE) are chosen to evaluate the accuracy level of the MVNCD evaluation results. The “%MAE > 0.005” represents the percentage in which the MAE was over 0.005, and “%MAPE > 2” represents the percentage in which the MAPE was over 2.Table 1
BVN and TVN evaluation results for the different methods.
DimMeasureNumerical integration methodsAnalytic approximation methods (TensorFlow)GaussTFMEBMETVBSK = 2 (BVN)MAE4.8∗e − 101.9∗e − 150.0037——%MAE > 0.0050.00.018.8——MAPE8.6∗e − 81.3∗e − 121.28——%MAPE > 20.00.015.2——Time (s)0.0010.0020.003——K = 3 (TVN)MAE8.8∗e − 106.1∗e − 110.00270.0018—%MAE > 0.0050.00.018.112.4—MAPE1.9∗e − 72.9∗e − 81.341.00—%MAPE > 20.00.016.711.0—Time (s)0.0050.0050.0060.003—Note. “—“ indicates that the corresponding analytic approximation methods cannot calculate the values of BVN or TVN.Table 2
MVNCD evaluation results for the different analytic approximation methods.
DimMeasureMEBMETVBSGaussTFGaussTFGaussTFK = 4MAE0.002260.002260.001520.001520.001170.00117MAPE0.590.590.400.400.260.26%MAE > 0.00515.215.29.89.85.85.8%MAPE > 218.018.012.812.88.48.4Time (s)0.0520.0100.0330.0040.0390.007K = 5MAE0.001710.001710.001420.001420.001080.00108MAPE0.540.540.440.440.330.33%MAE > 0.0059.29.27.07.04.84.8%MAPE > 218.818.815.415.413.013.0Time (s)0.0770.0160.0600.0070.0650.011K = 7MAE0.001200.001200.001040.001040.000930.00093MAPE0.450.450.380.380.330.33%MAE > 0.0055.95.94.74.73.63.6%MAPE > 221.521.520.220.218.318.3Time (s)0.1410.0310.0970.0140.1050.021K = 10MAE0.000640.000640.000570.000570.000540.00054MAPE0.320.320.290.290.290.29%MAE > 0.0051.61.61.21.21.11.1%MAPE > 225.025.023.123.122.622.6Time (s)0.2580.0610.1560.0230.1820.033K = 15MAE0.000410.000410.000400.000400.000370.00037MAPE0.290.290.280.280.250.25%MAE > 0.0050.20.20.20.20.10.1%MAPE > 234.134.132.632.632.032.0Time (s)0.5580.1460.3170.0580.3520.076From Table1, the following findings can be obtained. Compared with the calculated results by the cdfBvn and cdfTvn functions in the GAUSS (https://www.aptech.com/wp-content/uploads/2016/05/GAUSS_16_PDF_UG.pdf), the accuracy levels (shown as the MAE) of BVN and TVN probability (based on the improvement in Section 2.3) calculated with the assistance of TensorFlow have been improved from e−10 to e−15 and from e−10 to e−11, respectively, with a slight loss in time. When the integration dimension K is 2 or 3, as shown in Table 1, the accuracy level of the MAE for the ME and BME methods is e−3, which is far lower than that of bivariate and trivariate normal distribution function calculated by numerical integration. Similar conclusion can be obtained by choosing the MAPE as an index. Therefore, numerical integration should be adopted to calculate BVN or TVN values.WhenK ≥ 4, we employed the LDLT matrix-based analytic approximation method for the evaluation on Gauss and TensorFlow platforms. From Table 2, when the number of dimensions “K” increases, the performance of the three methods can be indicated by “MAE,” “MAPE,” “%MAE > 0.005,” and “%MAPE > 2.” The TVBS method is taken as an example, and there is a variation as the dimension “K” goes from 4 to 15, with “MAE” reducing from 0.00117 to 0.00037, “%MAE > 0.005” reducing from 5.8% to 0.1%, and “%MAPE > 2” increasing from 8.4% to 32.0%. As shown in Figure 4(a), “%MAE > 0.005” gradually decreases and “%MAPE > 2” gradually increases with the increase in dimension of all three methods. To explain the intrinsic reason, we checked the data and found that the overall average value of 1000 MVNCD values calculated by the method of generating random number would decrease with the increase in dimension. Due to the lower MVNCD values in the higher dimension, the relative error (MAPE) may become larger with the higher dimension. Meanwhile, as expected, with the increase in K, the time required by the analytic approximation gradually increased for all the methods. In addition, when we compare the three analytical approximation methods with each other at the same dimension, it can be found that the time consumed by the BME method is the shortest. In general, the time consumed by BME is less than half of that by the ME method and that by the TVBS method falls in the middle. Compared with the ME method, the error of the BME and TVBS methods is relatively low, as both of them are based on the bivariate truncation scheme. Overall, the TVBS is the best in terms of accuracy in evaluating individual MVNCD functions. These findings are consistent with Bhat’s study (2018) and mutually corroborate.Figure 4
Comparison of accuracy and time with different dimensions. (a) “%MAE > 0.005” and “%MAPE > 2” (b) time.
(a)(b)Furthermore, there are some findings by comparing the calculation error and calculation time of the different analytic approximation methods (i.e., ME, BME, and TVBS) in Gauss and TensorFlow. In terms of calculation accuracy, it can be seen in Table2 that “MAE,” “MAPE,” “%MAE > 0.005,” and “%MAPE > 2” of the three analytic approximation methods are basically the same. Therefore, although the new method improves accuracy levels of BVN and TVN probability in TensorFlow, its contribution to accuracy improvement is limited for high-dimensional analytic approximation; in terms of calculation time, it can be seen in Table 2 and Figure 4(b) that under the same dimension and algorithm, the calculation speed of TensorFlow platform is much higher than that of Gauss platform, and the average running time of Gauss is about 5-6 times that of TensorFlow. As shown in Table 1, there is no significant difference in the calculation time of BVN and TVN values on the two platforms. Thus, the declines in calculation time are mainly due to the obvious time advantage of TensorFlow with tensor data structure and graph computation technique.


## 4.2. Simulation Results for ICLV Model

Following Section3, the two ICLV models are designed to evaluate the ability for recovering the model parameters. For the three-alternative ICLV model with the maximum integral dimension being 3, the performance of the improved trivariate normal integral approach in estimating the parameters and the standard errors of the ICLV model is evaluated in Table 3. The evaluation system can be adopted as follows:(1)
For the parameter estimates, the mean estimate for each model parameter across the 50 datasets can be obtained. Then, the absolute bias of the estimator is computed asAbs.Bias=estimate value−true value and the absolute percentage bias as APB=estimate value−true value/true value×100.(2)
For the standard error estimates, the standard deviation of each estimated parameter across the 50 datasets can be computed, and this is labeled as the finite sample standard deviation or FSSD (essentially, this is the empirical standard error); also, we can compute the mean standard deviation for each model parameter across the 50 datasets and denote this as the asymptotic standard error or ASE (the standard error here is computed using the Godambe (sandwich) matrix at convergence for a theoretical approximation to the FSSD). Next, to evaluate the accuracy of the ASE formula for the finite sample size used, the relative efficiency of the estimator is computed asrelative efficiency=ASE/FSSD. As per Bhat [10], the closeness between ASE and FSEE may be examined as captured by the 0.75–1.25 range for the relative efficiency value.Table 3
Simulation results for the three-alternative ICLV model.
Parm.True valueParameter estimatesStandard error estimatesMean est.Abs. biasAPB (%)ASEFSSDR. effα10.50.53000.03005.99240.09610.08071.1908α20.60.64110.04116.85060.11770.10031.1730α30.50.50220.00220.44320.20370.14931.3648α40.60.60170.00170.27750.24510.18451.3288α50.30.31250.01254.16300.04080.04660.8740α60.30.30270.00270.91190.13360.11181.1950α7−0.4−0.40610.00611.51860.17730.14551.2182s11——————r120——————r130.60.68550.085514.24850.25670.25161.0206r140——————s21——————r230——————r240.60.62300.02303.84110.48000.39511.2149s31——————r340——————s41——————δ1−1−1.00130.00130.13120.03660.03341.0974δ2−1−1.02440.02442.44180.06800.06381.0664δ3−1−1.00310.00310.31170.04410.04530.9720δ4−1−1.06650.06656.65130.16290.16091.0124d10.30.28800.01203.99110.05080.04411.1524d20.40.44590.045911.48470.15500.13551.1441d30.50.49240.00761.52960.07540.08320.9064d40.60.71570.115719.27580.30460.26191.1633ψ11.51.49710.00300.19670.02750.02651.0375ψ21.51.53540.03542.35780.08730.08231.0601ψ31.51.49490.00510.34310.05040.05490.9191ψ41.51.61410.11417.60540.23460.21141.1101β20.50.48630.01372.74090.09540.08881.0741β3−1−0.99720.00280.27840.13600.15300.8893βtt−1−0.99800.00200.19710.04190.04340.9667βtc−0.8−0.79780.00220.27720.03380.03430.9852λ10.50.48070.01933.86940.08850.08111.0907λ20.20.21290.01296.44770.06610.05171.2771λ30.50.49290.00711.41460.07180.07230.9929λ40.20.21120.01125.61360.06220.05431.1441σ111——————σ120——————σ130——————σ221——————σ230.60.60650.00651.07940.13740.11961.1483σ331.361.37530.01531.12860.27420.25711.0668Overall mean value0.02363.79400.13080.11691.1190Note: “--” indicates that the corresponding parameter is fixed at the true value.The following aspects can be observed in Table3 (the experiment is repeated 50 times with a sample size of 10,000). Firstly, the overall mean value across parameters of APB is 3.7940%, with the individual parameter APB ranging from 0.1% to 19%, which are close to those in Bhat. Besides, the simulation errors are also similar and therefore considered reasonable and acceptable; secondly, the mean APB values for the d vector (9.0%), λ vector (4.3%), and r in Γ matrix (9.0%) are relatively higher than the overall mean APB. This indicates that the coefficients on the latent variable vector z∗ in the measurement equation and the choice model are somewhat more difficult to recover than other parameters. It is expected, since these elements enter into the covariance matrix Ω in a nonlinear fashion as shown in equation (6), and Ω itself enters into the composite likelihood function (equation (8)) in a complex manner. Lastly, as for simulation error estimates, the overall mean value of relative efficiency (RE) across the parameters is 1.1190 (and the RE of almost all parameters is between 0.75 and 1.25), which suggests that the FSSE and ASE values are close to each other and the asymptotic formula is performing well in estimating the finite sample standard errors.The simulation results (with a sample size of 100,000) of the four-alternative ICLV model estimated by the three analytical approximation methods, namely the ME, BME, and TVBS, are presented in Tables4 –6, respectively. In these tables, Abs.Bias=estimate value−true value and APB=estimate value−true value/truevalue×100. As shown in the experimental results, it is ensured that all the identifiable parameters can be significantly estimated and seemingly consistent with their true values. We can see that the values of absolute percentage bias (APB) fall in the ranges of 0–15% in the three tables, of which the APB of four parameters exceeds 10% (i.e., α3, α4, d2, and σ34). It is a similar finding showing that the parameter estimators related to z∗ are less consistent. Additionally, in these three tables (also shown in Figure 5), it can be observed that among the three methods (i.e., ME, BME, and TVBS), the overall mean value of absolute bias is 0.0266, 0.0262, and 0.0259, and the average value of APB is 4.269%, 4.240%, and 4.204%, respectively. It is noticed that the values of these two indicators slightly decrease. On the other hand, the standard errors of all methods are about 0.037. It shows that the standard errors of the four-alternative ICLV model seem reasonable and acceptable. Overall, we have validated the ability for recovering parameters of these three LDLT-based analytic approximation methods for the ICLV model, and we can conclude that TVBS is superior to BME and BME is superior to ME. These findings are consistent with the evaluation results with respect to individual MVNCD functions. Although model parameters can be estimated significantly with a sample of large size, there are issues of empirical identification, namely whether the data can support the model specification. If the empirical data cannot allow for significantly estimating all the theoretically identifiable parameters, some parameters may be fixed and only parameters of particular interest need to be estimated.Table 4
Simulation results for the four-alternative ICLV model with the ME method.
Parm.True valueEstimateAbs. biasAPB (%)Std. err.t-testP valueα10.50.48210.01793.58060.017527.4871≤0.001α20.60.57340.02664.43450.021027.2569≤0.001α30.50.56230.062312.45820.054510.3235≤0.001α40.60.66590.065910.98580.064510.3299≤0.001α50.30.30390.00391.28330.010329.4528≤0.001α60.30.30590.00591.96930.030610.0031≤0.001α7−0.4−0.40990.00992.48170.0408−10.0422≤0.001s11——————r120——————r130.60.59820.00180.30820.046013.0071≤0.001r140——————s21——————r230——————r240.60.64750.04757.92420.10076.4324≤0.001s31——————r340——————s41——————δ1−1−1.00680.00680.68180.0115−87.6878≤0.001δ2−1−0.96710.03293.29330.0161−59.9248≤0.001δ3−1−1.00810.00810.80700.0131−76.7456≤0.001δ4−1−0.98750.01251.25310.0373−26.4410≤0.001d10.30.31470.01474.91600.012525.2686≤0.001d20.40.35080.049212.30200.03819.2070≤0.001d30.50.49960.00040.08280.021023.7820≤0.001d40.60.59100.00901.50230.07937.4520≤0.001ψ11.51.51200.01200.80040.0081186.2783≤0.001ψ21.51.47920.02081.38350.018878.6552≤0.001ψ31.51.49990.00010.00340.0145103.2049≤0.001ψ41.51.51580.01581.05530.053828.1516≤0.001β20.50.45470.04539.05600.030414.9540≤0.001β3−1−1.04160.04164.16160.0476−21.8600≤0.001β4−0.8−0.87280.07289.10370.0445−19.5951≤0.001βtt−1−1.01150.01151.14900.0158−63.8326≤0.001βtc−0.8−0.81070.01071.34240.0128−63.5246≤0.001λ10.50.53800.03807.60640.023422.9537≤0.001λ20.20.18580.01427.07700.018610.0106≤0.001λ30.20.18770.01236.14750.018710.0234≤0.001λ40.50.48910.01092.17560.023221.0530≤0.001λ50.50.51030.01032.06320.021224.0984≤0.001λ60.20.19970.00030.15650.02049.7718≤0.001σ111——————σ120——————σ130——————σ140——————σ221——————σ230.60.59560.00440.73650.052611.3276≤0.001σ240.60.64070.04076.77720.046013.9313≤0.001σ331.361.40080.04082.99990.107413.0425≤0.001σ340.60.69390.093915.64220.07798.9014≤0.001σ441.361.47230.11238.25420.100714.6142≤0.001Overall mean value0.02664.26910.0371Note: “—” indicates that the corresponding parameter is fixed at the true value.Table 5
Simulation results for the four-alternative ICLV model with the BME method.
Parm.True valueEstimateAbs. biasAPB (%)Std. err.t-testP valueα10.50.48210.01793.57600.017527.4890≤0.001α20.60.57340.02664.43070.021027.2577≤0.001α30.50.56220.062212.43020.054010.4099≤0.001α40.60.66580.065810.95830.063910.4166≤0.001α50.30.30440.00441.45370.010329.4308≤0.001α60.30.30600.00602.01600.03069.9954≤0.001α7−0.4−0.41010.01012.52800.0409−10.0344≤0.001s11——————r120——————r130.60.59900.00100.16620.046112.9991≤0.001r140——————s21——————r230——————r240.60.64770.04777.95000.10046.4540≤0.001s31——————r340——————s41——————δ1−1−1.00680.00680.68160.0115−87.6890≤0.001δ2−1−0.96710.03293.29030.0161−60.1612≤0.001δ3−1−1.00760.00760.76470.0131−76.8664≤0.001δ4−1−0.98730.01271.26820.0373−26.4510≤0.001d10.30.31470.01474.91100.012525.2700≤0.001d20.40.35090.049112.27750.03789.2838≤0.001d30.50.49850.00150.29240.021023.7831≤0.001d40.60.59060.00941.56250.07937.4485≤0.001ψ11.51.51200.01200.80000.0081186.2896≤0.001ψ21.51.47930.02071.38040.018779.1925≤0.001ψ31.51.49930.00070.04560.0145103.4676≤0.001ψ41.51.51560.01561.03930.053828.1629≤0.001β20.50.45460.04549.07900.030414.9559≤0.001β3−1−1.04170.04174.17480.0476−21.8780≤0.001β4−0.8−0.87190.07198.99350.0445−19.5883≤0.001βtt−1−1.01120.01121.12150.0158−63.9360≤0.001βtc−0.8−0.81060.01051.31870.0127−63.6256≤0.001λ10.50.53790.03797.57200.023422.9583≤0.001λ20.20.18590.01417.04550.018410.0864≤0.001λ30.20.18770.01236.16000.018610.1024≤0.001λ40.50.48820.01182.36540.023221.0692≤0.001λ50.50.50940.00941.87680.021124.1120≤0.001λ60.20.19950.00050.26650.02049.7651≤0.001σ111——————σ120——————σ130——————σ140——————σ221——————σ230.60.59520.00480.79270.052511.3353≤0.001σ240.60.63770.03776.28870.045713.9456≤0.001σ331.361.40020.04022.95340.107313.0519≤0.001σ340.60.69300.093015.49500.07768.9351≤0.001σ441.361.46280.10287.55510.100114.6143≤0.001Overall mean value0.02624.24000.0370Note: “—” indicates that the corresponding parameter is fixed at the true value.Table 6
Simulation results for the four-alternative ICLV model with the TVBS method.
Parm.True valueEstimateAbs. biasAPB (%)Std. err.t-testP valueα10.50.48200.01803.60060.017527.4956≤0.001α20.60.57330.02674.45500.021027.2643≤0.001α30.50.56210.062112.42480.054110.3956≤0.001α40.60.66570.065710.95320.064010.4023≤0.001α50.30.30420.00421.39130.010429.3686≤0.001α60.30.30580.00581.93330.030510.0157≤0.001α7−0.4−0.40980.00982.44530.0408−10.0548≤0.001s11——————r120——————r130.60.59860.00140.22630.046013.0038≤0.001r140——————s21——————r230——————r240.60.64720.04727.86280.10026.4570≤0.001s31——————r340——————s41——————δ1−1−1.00680.00680.68440.0115−87.6833≤0.001δ2−1−0.96710.03293.28980.0161−60.1169≤0.001δ3−1−1.00780.00780.78070.0131−76.7399≤0.001δ4−1−0.98760.01241.24080.0373−26.4477≤0.001d10.30.31480.01484.94070.012525.2747≤0.001d20.40.35090.049112.27280.03799.2709≤0.001d30.50.49890.00110.21540.021023.7248≤0.001d40.60.59130.00871.45550.07937.4595≤0.001ψ11.51.51200.01200.80260.0081186.2650≤0.001ψ21.51.47930.02071.37990.018779.0895≤0.001ψ31.51.49960.00040.02990.0145103.1776≤0.001ψ41.51.51600.01601.06770.053828.1597≤0.001β20.50.45440.04569.11440.030414.9438≤0.001β3−1−1.04250.04254.24790.0476−21.8930≤0.001β4−0.8−0.87280.07289.10250.0445−19.6065≤0.001βtt−1−1.01160.01161.16070.0158−63.9212≤0.001βtc−0.8−0.81090.01091.35670.0127−63.6051≤0.001λ10.50.53820.03827.64040.023422.9612≤0.001λ20.20.18600.01406.98350.018510.0749≤0.001λ30.20.18780.01226.11800.018610.0900≤0.001λ40.50.48870.01132.25980.023221.0338≤0.001λ50.50.50990.00991.98140.021224.0589≤0.001λ60.20.19970.00030.15100.02049.7806≤0.001σ111——————σ120——————σ130——————σ140——————σ221——————σ230.60.59230.00771.28320.052411.3058≤0.001σ240.60.63820.03826.37020.045813.9236≤0.001σ331.361.38950.02952.16590.107412.9430≤0.001σ340.60.68750.087514.58380.07758.8672≤0.001σ441.361.46290.10297.56430.100114.6142≤0.001Overall mean value0.02594.20370.0370Note. “—” indicates that the corresponding parameter is fixed at the true value.Figure 5
Comparison of three methods for the four-alternative ICLV model: (a) overall mean value of Abs. bias and (b) overall mean value of APB (%).
(a)(b)Furthermore, more attempts were made to discuss the operational efficiency of models using the automatic differential technique. The convergence time of the three-alternative ICLV model in Table3 has a median value of about 8 minutes for the case of 10,000 observations, which is much shorter than that of the literature of Bhat in 2014 (about one hour for the case of 1000 observations). It can also be inferred that the computational time of the improved four-alternative ICLV model based on TensorFlow is reduced. This can be attributed to two factors: one is that the LDLT-based analytic approximation method (instead of the MACML approach) can shorten the calculation time of the ICLV model, as mentioned in Bhat. The improvement on bivariate normal CDF and trivariate normal CDF also facilitates this process; the other is based on the time advantage of TensorFlow in calculating tensor and multidimensional arrays. In addition, in the absence of coding the analytic gradient artificially, the convergence speed of the model can be improved using automatic differentiation to obtain the gradient of the likelihood function. For the four-alternative ICLV model with a sample size of 100,000, the convergence times of three analytical approximation methods are 3 h7 min (for ME), 1 h14 min (for BME), and 2 h6 min (for TVBS), respectively. Among them, the BME algorithm has obvious advantages in terms of time consumption. Hence, if the accuracy level can be lowered under some situations (e.g., in a variable screening process), the BME algorithm can also be chosen to save computational time. In this section, it is not necessary to evaluate the ICLV model of Bhat for comparison, because Bhat has explicitly shown that the LDLT matrix-based analytic approximation method runs much faster than other methods (such as the SSJ method) in evaluating MVNCD functions and estimating MNP model. It is worth noting that all simulation experiments are conducted on a desktop computer with an Intel Core™ i5-4460@3.20 GHz Processor and 16 GB of RAM.


## 5. Conclusions and Discussion

On the basis of the research advances of Bhat in 2018 and the automatic differentiation (AD) technique on the TensorFlow platform, this study integrates the LDLT matrix-based analytic approximation method into the probit kernel-based ICLV formulation, which can further improve the performance of the ICLV model in terms of accuracy level and computational speed. The main achievements of this study can be summarized as follows.Firstly, the individual MVNCD function was evaluated with different dimensions. By improving the numerical integration algorithm of CDF of BVN and TVN, the accuracy level of BVN and TVN functions was enhanced in this study (compared with algorithms in GAUSS). Based on the TensorFlow platform, we also tested the errors and operation time of three types of matrix-based analytic approximation methods in various dimensions. It is found that the operation time of BME is the shortest and the accuracy level of TVBS is the highest in the same dimension. By comparing the calculation of MVNCD functions on Gauss and TensorFlow platforms, it shows that the calculation speed of TensorFlow is much higher than that of Gauss.Subsequently, we designed the simulation experiments of the three-alternative ICLV model and four-alternative ICLV model. The estimates and standard errors of both models are within a reasonable range, of which the parameters related to latent variables may be more difficult to recover. For the three-alternative ICLV model, the numerical integration method of calculating TVN value was employed and the simulation experiments were repeated 50 times with a sample size of 10,000. The simulation results show that the overall mean value of APB is 3.79% and the relative efficiency is 1.1190, exhibiting a desirable capability of parameter recovery. After that, this study applies three analytic approximation methods (i.e., ME, BME, and TVBS) in the MVNCD function calculation of the ICLV model and conducts simulation experiments with a sample size of 100,000 regarding the four-alternative ICLV model, and significant results of consistent estimation value can be obtained. The average values of absolute bias and APB decreased in the order of ME, BME, and TVBS, indicating the increased accuracy. To sum up, when we estimate the ICLV model with high-dimensional integration, the TVBS method can be chosen with the highest accuracy level. Also, the BME method is an alternative to save computational time.Meanwhile, with the assistance of the AD technique, the researchers are freed from coding analytical gradients. Compared with the probit-based ICLV model proposed by Bhat, the computational time of the proposed model is reduced remarkably, which is mainly attributed to two reasons: one is adopting the LDLT matrix-based analytic approximation methods to evaluate MVNCD function; the other is based on the advantages of TensorFlow in calculating tensor (represented as multidimensional array). The operation of a three-alternative ICLV model with a sample size of 10,000 only takes 8 min with the time advantage of LDLT-based analytic approximation methods, which can facilitate the applications of the proposed method in practice.In the future work, we will further explore the practical application of the improved ICLV model with multiple alternatives and incorporate more psychological attitude variables into the model so as to explore the potential influencing mechanism of travel behaviors by using real data. If possible, further comparisons in terms of parameter estimation results between the ICLV model and previous models will be discussed.



---

*Source: 6556282-2022-08-02.xml*