MDC Procedure

HEV and Multinomial Probit: Heteroscedastic Utility Function

When the stochastic components of utility are heteroscedastic and independent, you can model the data by using an HEV or a multinomial probit model. The HEV model assumes that the utility of alternative j for each individual i has heteroscedastic random components,

upper U Subscript i j Baseline equals upper V Subscript i j Baseline plus epsilon Subscript i j

where the cumulative distribution function of the Gumbel distributed epsilon Subscript i j is

upper F left-parenthesis epsilon Subscript i j Baseline right-parenthesis equals exp left-parenthesis minus exp left-parenthesis minus epsilon Subscript i j Baseline slash theta Subscript j Baseline right-parenthesis right-parenthesis

Note that the variance of epsilon Subscript i j is one-sixth pi squared theta Subscript j Superscript 2. Therefore, the error variance is proportional to the square of the scale parameter theta Subscript j. For model identification, at least one of the scale parameters must be normalized to 1. The following SAS statements estimate an HEV model under a unit scale restriction for mode "1" (theta 1 equals 1):

/*-- hev with gauss-laguerre method --*/
proc mdc data=newdata;
   model decision = ttime /
            type=hev
            nchoice=3
            hev=(unitscale=1, integrate=laguerre)
            covest=hess;
   id pid;
run;

The results of computation are presented in Figure 14 and Figure 15.

Figure 14: HEV Estimation Summary

The MDC Procedure
 
Heteroscedastic Extreme Value Model Estimates

Model Fit Summary
Dependent Variabledecision
Number of Observations50
Number of Cases150
Log Likelihood-33.41383
Maximum Absolute Gradient0.0000218
Number of Iterations11
Optimization MethodDual Quasi-Newton
AIC72.82765
Schwarz Criterion78.56372


Figure 15: HEV Parameter Estimates

The MDC Procedure
 
Heteroscedastic Extreme Value Model Estimates

Parameter Estimates
ParameterDFEstimateStandard
Error
t ValueApprox
Pr > |t|
ttime1-0.44070.1798-2.450.0143
SCALE210.77650.43481.790.0741
SCALE310.57530.27522.090.0366


The parameters SCALE2 and SCALE3 in the output correspond to the estimates of the scale parameters theta 2 and theta 3, respectively.

Note that the estimate of the HEV model is not always stable because computation of the log-likelihood function requires numerical integration. Bhat (1995) proposed the Gauss-Laguerre method. In general, the log-likelihood function value of HEV should be larger than that of conditional logit because HEV models include the conditional logit as a special case. However, in this example the reverse is true (–33.414 for the HEV model, which is less than –33.321 for the conditional logit model). (See Figure 14 and Figure 3.) This indicates that the Gauss-Laguerre approximation to the true probability is too coarse. You can see how well the Gauss-Laguerre method works by specifying a unit scale restriction for all modes, as in the following statements, since the HEV model with the unit variance for all modes reduces to the conditional logit model:

/*-- hev with gauss-laguerre and unit scale --*/
proc mdc data=newdata;
   model decision = ttime /
            type=hev
            nchoice=3
            hev=(unitscale=1 2 3, integrate=laguerre)
            covest=hess;
   id pid;
run;

Figure 16 shows that the ttime coefficient is not close to that of the conditional logit model.

Figure 16: HEV Estimates with All Unit Scale Parameters

The MDC Procedure
 
Heteroscedastic Extreme Value Model Estimates

Parameter Estimates
ParameterDFEstimateStandard
Error
t ValueApprox
Pr > |t|
ttime1-0.29260.0438-6.68<.0001


There is another option of specifying the integration method. The INTEGRATE=HARDY option uses the adaptive Romberg-type integration method. The adaptive integration produces much more accurate probability and log-likelihood function values, but often it is not practical to use this method of analyzing the HEV model because it requires excessive CPU time. The following SAS statements produce the HEV estimates by using the adaptive Romberg-type integration method:

/*-- hev with adaptive integration --*/
proc mdc data=newdata;
   model decision = ttime /
             type=hev
             nchoice=3
             hev=(unitscale=1, integrate=hardy)
             covest=hess;
   id pid;
run;

The results are displayed in Figure 17 and Figure 18.

Figure 17: HEV Estimation Summary Using Alternative Integration Method

The MDC Procedure
 
Heteroscedastic Extreme Value Model Estimates

Model Fit Summary
Dependent Variabledecision
Number of Observations50
Number of Cases150
Log Likelihood-33.02598
Maximum Absolute Gradient0.0001202
Number of Iterations8
Optimization MethodDual Quasi-Newton
AIC72.05197
Schwarz Criterion77.78803


Figure 18: HEV Estimates Using Alternative Integration Method

The MDC Procedure
 
Heteroscedastic Extreme Value Model Estimates

Parameter Estimates
ParameterDFEstimateStandard
Error
t ValueApprox
Pr > |t|
ttime1-0.45800.1861-2.460.0139
SCALE210.77570.42831.810.0701
SCALE310.69080.33842.040.0412


With the INTEGRATE=HARDY option, the log-likelihood function value of the HEV model, –33.026, is greater than that of the conditional logit model, –33.321. (See Figure 17 and Figure 3.)

When you impose unit scale restrictions on all choices, as in the following statements, the HEV model gives the same estimates as the conditional logit model. (See Figure 19 and Figure 6.)

/*-- hev with adaptive integration and unit scale --*/
proc mdc data=newdata;
   model decision = ttime /
            type=hev
            nchoice=3
            hev=(unitscale=1 2 3, integrate=hardy)
            covest=hess;
   id pid;
run;

Figure 19: Alternative HEV Estimates with Unit Scale Restrictions

The MDC Procedure
 
Heteroscedastic Extreme Value Model Estimates

Parameter Estimates
ParameterDFEstimateStandard
Error
t ValueApprox
Pr > |t|
ttime1-0.35720.0776-4.60<.0001


For comparison, the following statements estimate a heteroscedastic multinomial probit model by imposing a zero restriction on the correlation parameter, rho 31 equals 0. The MDC procedure requires normalization of at least two of the error variances in the multinomial probit model. Also, for identification, the correlation parameters associated with a unit normalized variance are restricted to be zero. When the UNITVARIANCE= option is specified, the zero restriction on correlation coefficients applies to the last choice of the list. In the following statements, the variances of the first and second choices are normalized. The UNITVARIANCE=(1 2) option imposes additional restrictions that rho 32 equals rho 21 equals 0. The default for the UNITVARIANCE= option is the last two choices (which would have been equivalent to UNITVARIANCE=(2 3) for this example). The result is presented in Figure 20.

The utility function can be defined as

upper U Subscript i j Baseline equals upper V Subscript i j Baseline plus epsilon Subscript i j

where

bold-italic epsilon Subscript i Baseline tilde upper N left-parenthesis bold 0 comma Start 3 By 3 Matrix 1st Row 1st Column 1 2nd Column 0 3rd Column 0 2nd Row 1st Column 0 2nd Column 1 3rd Column 0 3rd Row 1st Column 0 2nd Column 0 3rd Column sigma 3 squared EndMatrix right-parenthesis
/*-- mprobit estimation --*/
proc mdc data=newdata;
   model decision = ttime /
            type=mprobit
            nchoice=3
            unitvariance=(1 2)
            covest=hess;
   id pid;
   restrict RHO_31 = 0;
run;

Figure 20: Heteroscedastic Multinomial Probit Estimates

The MDC Procedure
 
Multinomial Probit Estimates

Parameter Estimates
ParameterDFEstimateStandard
Error
t ValueApprox
Pr > |t|
Parameter Label
ttime1-0.32060.0920-3.490.0005 
STD_311.69130.69062.450.0143 
RHO_31000   
Restrict111.18541.54900.770.4499*Linear EC [ 1 ]

* Probability computed using beta distribution.



Note that in the output the estimates of standard errors and correlations are denoted by STD_i and RHO_ij, respectively. In this particular case the first two variances (STD_1 and STD_2) are normalized to one, and corresponding correlations (RHO_21 and RHO_32) are set to zero, so they are not listed among parameter estimates.

Last updated: June 19, 2025