A Generalized Continuous-Multinomial Response Model with a t-distributed Error Kernel

TL;DR

This paper introduces a generalized multinomial response model with a t-distributed error kernel, improving preference estimation especially in class-imbalanced datasets by accounting for decision-uncertainty behavior.

Contribution

It extends the GCM model to include a t-distributed error kernel, providing a new estimator and demonstrating its advantages over the normal distribution in preference modeling.

Findings

GCM-t outperforms GCM-N in preference prediction accuracy.

Accounting for decision-uncertainty affects elasticity and willingness-to-pay estimates.

The model is applied to electric vehicle preferences, influencing policy insights.

Abstract

In multinomial response models, idiosyncratic variations in the indirect utility are generally modeled using Gumbel or normal distributions. This study makes a strong case to substitute these thin-tailed distributions with a t-distribution. First, we demonstrate that a model with a t-distributed error kernel better estimates and predicts preferences, especially in class-imbalanced datasets. Our proposed specification also implicitly accounts for decision-uncertainty behavior, i.e. the degree of certainty that decision-makers hold in their choices relative to the variation in the indirect utility of any alternative. Second, after applying a t-distributed error kernel in a multinomial response model for the first time, we extend this specification to a generalized continuous-multinomial (GCM) model and derive its full-information maximum likelihood estimator. The likelihood involves an open-form expression of the cumulative density function of the multivariate t-distribution, which we propose to compute using a combination of the composite marginal likelihood method and the separation-of-variables approach. Third, we establish finite sample properties of the GCM model with a t-distributed error kernel (GCM-t) and highlight its superiority over the GCM model with a normally-distributed error kernel (GCM-N) in a Monte Carlo study. Finally, we compare GCM-t and GCM-N in an empirical setting related to preferences for electric vehicles (EVs). We observe that accounting for decision-uncertainty behavior in GCM-t results in lower elasticity estimates and a higher willingness to pay for improving the EV attributes than those of the GCM-N model. These differences are relevant in making policies to expedite the adoption of EVs.

Tables

Table 1. Table 1 : Class-imbalance example under Probit and t-distributed error kernel

Demographic profile	t-distributed kernel						Probit kernel
	DOF = 0.1		DOF = 0.5		DOF = 1.0
	P(car)	P(bicycle)	P(car)	P(bicycle)	P(car)	P(bicycle)	P(car)	P(bicycle)
High-income and non-student	0.62	0.38	0.80	0.20	0.88	0.12	0.99	0.01
High-income and student	0.60	0.40	0.76	0.24	0.83	0.17	0.96	0.04
Low-income and non-student	0.58	0.42	0.70	0.30	0.75	0.25	0.84	0.16
Low-income and student	0.54	0.46	0.58	0.42	0.59	0.41	0.62	0.38

Table 2. Table 2 : Simulation Results of GCM-t Model for DOF-I scenario (DOF=1)

Parameter	True Value	MEV	MAB	APB	FSSE	ASE	ASE/ FSSE	CP	Power
Mean effect elements
${\gamma }_{11}$	1	0.98	0.02	1.69	0.06	0.06	1.07	0.93	1
${\gamma }_{12}$	0.5	0.51	0.01	2.06	0.06	0.06	1.12	0.97	1
${\gamma }_{13}$	0.75	0.76	0.01	1.46	0.07	0.06	0.98	0.93	1
${\gamma }_{14}$	-0.5	-0.5	0	0.94	0.06	0.06	1.06	0.97	1
${\beta }_{21}$	-1.5	-1.74	0.24	16.23	0.28	0.16	0.56	0.47	1
${\beta }_{22}$	1	0.91	0.09	8.71	0.11	0.11	0.99	0.89	1
${\beta }_{23}$	0.9	0.81	0.09	10.35	0.13	0.12	0.98	0.9	1
${\beta }_{25}$	1	0.94	0.06	6.46	0.26	0.16	0.63	0.72	0.99
${\beta }_{31}$	-1.3	-1.47	0.17	13.21	0.26	0.2	0.77	0.85	1
${\beta }_{32}$	0.9	0.79	0.11	12.35	0.14	0.12	0.87	0.81	1
${\beta }_{33}$	0.8	0.68	0.12	15.58	0.16	0.13	0.84	0.79	1
${\beta }_{35}$	0.9	0.8	0.1	10.9	0.23	0.17	0.73	0.8	0.99
${\beta }_{41}$	-1.2	-1.36	0.16	13.49	0.27	0.2	0.76	0.79	1
${\beta }_{42}$	0.8	0.67	0.13	16.54	0.15	0.13	0.88	0.76	1
${\beta }_{43}$	0.7	0.57	0.13	18.7	0.19	0.15	0.81	0.8	0.93
${\beta }_{45}$	0.8	0.68	0.12	15.23	0.25	0.18	0.74	0.77	0.92
${\beta }_{51}$	-1	-1.08	0.08	7.72	0.22	0.18	0.83	0.91	1
${\beta }_{52}$	0.7	0.54	0.16	23.33	0.15	0.12	0.85	0.65	0.99
${\beta }_{53}$	0.6	0.43	0.17	27.66	0.17	0.15	0.84	0.76	0.83
${\beta }_{55}$	0.7	0.56	0.14	20.14	0.23	0.18	0.81	0.81	0.8
Covariance matrix elements
${\overline{\Sigma }}_{11}$	1.5	1.62	0.12	8.28	0.1	0.09	0.97	0.75	1
${\overline{\Sigma }}_{21}$	0.3	0.28	0.02	8.22	0.29	0.18	0.61	0.77	0.42
${\overline{\Sigma }}_{31}$	0.4	0.44	0.04	9.78	0.3	0.2	0.66	0.79	0.6
${\overline{\Sigma }}_{33}$	1.1	1.21	0.11	10.41	0.29	0.24	0.83	0.89	1
${\overline{\Sigma }}_{41}$	0.6	0.72	0.12	20.24	0.35	0.24	0.7	0.78	0.76
${\overline{\Sigma }}_{44}$	1.2	1.48	0.28	23.19	0.4	0.33	0.83	0.91	1
${\overline{\Sigma }}_{51}$	0.5	0.6	0.1	19.12	0.34	0.24	0.72	0.82	0.65
${\overline{\Sigma }}_{55}$	1.3	1.48	0.18	14.06	0.41	0.31	0.76	0.9	1
Degree-of-freedom
$\delta$	1	1.08	0.08	8.24	0.03	0.04	1.1	0.43	1
Average			0.14	12.56	0.21	0.16	0.84	0.8	0.93

Table 3. Table 3 : Simulation Results of GCM-t Model for DOF-II scenario (DOF=12)

Parameter	True Value	MEV	MAB	APB	FSSE	ASE	ASE/ FSSE	CP	Power
Mean effect elements
${\gamma }_{11}$	1	1	0	0.38	0.04	0.05	1.09	0.97	1
${\gamma }_{12}$	0.5	0.5	0	0.77	0.05	0.05	0.99	0.97	1
${\gamma }_{13}$	0.75	0.75	0	0.42	0.05	0.05	0.94	0.91	1
${\gamma }_{14}$	-0.5	-0.5	0	0.26	0.05	0.05	0.92	0.91	1
${\beta }_{21}$	-1.5	-1.76	0.26	17.62	0.17	0.15	0.88	0.49	1
${\beta }_{22}$	1	0.89	0.11	11.48	0.11	0.11	0.96	0.8	1
${\beta }_{23}$	0.9	0.79	0.11	12.12	0.12	0.12	0.98	0.87	1
${\beta }_{25}$	1	0.9	0.1	10.28	0.18	0.19	1.08	0.94	1
${\beta }_{31}$	-1.3	-1.42	0.12	9.58	0.27	0.27	1.01	0.95	1
${\beta }_{32}$	0.9	0.76	0.14	16.03	0.13	0.13	1.02	0.77	1
${\beta }_{33}$	0.8	0.67	0.13	16.38	0.15	0.14	0.94	0.77	1
${\beta }_{35}$	0.9	0.74	0.16	17.89	0.18	0.19	1.06	0.83	0.99
${\beta }_{41}$	-1.2	-1.29	0.09	7.49	0.25	0.26	1.04	0.97	1
${\beta }_{42}$	0.8	0.63	0.17	21.82	0.13	0.12	0.92	0.65	1
${\beta }_{43}$	0.7	0.53	0.17	24.39	0.16	0.13	0.83	0.62	0.98
${\beta }_{45}$	0.8	0.64	0.16	20.04	0.19	0.19	1.04	0.87	0.93
${\beta }_{51}$	-1	-1.02	0.02	1.86	0.2	0.22	1.11	0.98	1
${\beta }_{52}$	0.7	0.5	0.2	28.48	0.12	0.11	0.94	0.52	1
${\beta }_{53}$	0.6	0.42	0.18	30.18	0.14	0.12	0.88	0.61	0.92
${\beta }_{55}$	0.7	0.53	0.17	24.78	0.18	0.19	1.02	0.79	0.81
Covariance matrix elements
${\overline{\Sigma }}_{11}$	1.5	1.52	0.02	1.22	0.06	0.06	1.05	0.97	1
${\overline{\Sigma }}_{21}$	0.3	0.17	0.13	43.74	0.21	0.19	0.9	0.85	0.19
${\overline{\Sigma }}_{31}$	0.4	0.34	0.06	15.32	0.23	0.2	0.87	0.9	0.41
${\overline{\Sigma }}_{33}$	1.1	1.1	0	0.04	0.41	0.41	1.01	0.93	1
${\overline{\Sigma }}_{41}$	0.6	0.55	0.05	8.92	0.29	0.23	0.81	0.87	0.67
${\overline{\Sigma }}_{44}$	1.2	1.21	0.01	0.53	0.48	0.48	1.01	0.9	1
${\overline{\Sigma }}_{51}$	0.5	0.43	0.07	13.6	0.26	0.22	0.84	0.89	0.49
${\overline{\Sigma }}_{55}$	1.3	1.2	0.1	7.58	0.4	0.45	1.1	0.91	1
Degree-of-freedom
$\delta$	12	13.48	1.48	12.37	2.89	2.84	0.98	0.97	1
Average			0.18	12.95	0.28	0.27	0.97	0.84	0.91

Table 4. Table 4 : Effect of ignoring non-normality (DOF=2)

Parameter	True Value	MEV	MAB	APB	MEV	MAB	APB
		GCM-t with DOF of 2			GCM-t with DOF of 300
Mean effect elements
${\gamma }_{11}$	1	0.98	0.02	1.53	0.98	0.02	1.6
${\gamma }_{12}$	0.5	0.51	0.01	2.12	0.51	0.01	2.51
${\gamma }_{13}$	0.75	0.76	0.01	0.99	0.76	0.01	1.93
${\gamma }_{14}$	-0.5	-0.5	0	0.67	-0.51	0.01	2.48
${\beta }_{21}$	-1.5	-1.75	0.25	16.91	-1.4	0.1	6.36
${\beta }_{22}$	1	0.9	0.1	10.13	0.66	0.34	33.79
${\beta }_{23}$	0.9	0.81	0.09	9.94	0.6	0.3	32.91
${\beta }_{25}$	1	0.91	0.09	8.95	0.65	0.35	35.18
${\beta }_{31}$	-1.3	-1.49	0.19	14.34	-1.33	0.03	2.44
${\beta }_{32}$	0.9	0.76	0.14	16.1	0.55	0.35	38.33
${\beta }_{33}$	0.8	0.67	0.13	16.66	0.49	0.31	38.6
${\beta }_{35}$	0.9	0.82	0.08	9.33	0.6	0.3	33.77
${\beta }_{41}$	-1.2	-1.36	0.16	13.28	-1.16	0.04	2.98
${\beta }_{42}$	0.8	0.64	0.16	20.45	0.44	0.36	45.38
${\beta }_{43}$	0.7	0.5	0.2	28.62	0.33	0.37	52.55
${\beta }_{45}$	0.8	0.71	0.09	10.75	0.5	0.3	37.95
${\beta }_{51}$	-1	-1.09	0.09	8.82	-2.65	1.65	164.78
${\beta }_{52}$	0.7	0.53	0.17	24.31	0.34	0.36	51.03
${\beta }_{53}$	0.6	0.42	0.18	30.29	0.25	0.35	58.77
${\beta }_{55}$	0.7	0.55	0.15	21.58	0.28	0.42	60.32
Covariance matrix elements
${\overline{\Sigma }}_{11}$	1.5	1.57	0.07	4.42	4.2	2.7	180.26
${\overline{\Sigma }}_{21}$	0.3	0.24	0.06	21.37	-0.37	0.67	222.38
${\overline{\Sigma }}_{31}$	0.4	0.32	0.08	19.77	-0.18	0.58	145.74
${\overline{\Sigma }}_{33}$	1.1	1.09	0.01	0.72	1.54	0.44	39.58
${\overline{\Sigma }}_{41}$	0.6	0.55	0.05	9.13	0.31	0.29	48.77
${\overline{\Sigma }}_{44}$	1.2	1.23	0.03	2.25	1.39	0.19	15.44
${\overline{\Sigma }}_{51}$	0.5	0.51	0.01	1.49	1.07	0.57	113.11
${\overline{\Sigma }}_{55}$	1.3	1.37	0.07	5.48	14.58	13.28	1021.4
Degree-of-freedom
$\delta$	2	2.12	0.12	5.84
Average			0.1	11.59		0.88	88.94
Loglikelihood		-10417.44			-10922.73

Table 5. Table 5 : Effect of ignoring non-normality (DOF=12)

Parameter	True Value	MEV	MAB	APB	MEV	MAB	APB
		GCM-t with DOF of 12			GCM-t with DOF of 300
Mean effect elements
${\gamma }_{11}$	1	1	0	0.39	1	0	0.32
${\gamma }_{12}$	0.5	0.49	0.01	1.46	0.49	0.01	1.77
${\gamma }_{13}$	0.75	0.75	0	0.45	0.75	0	0.03
${\gamma }_{14}$	-0.5	-0.5	0	0.26	-0.5	0	0.01
${\beta }_{21}$	-1.5	-1.75	0.25	16.83	-1.7	0.2	13.54
${\beta }_{22}$	1	0.89	0.11	11	0.85	0.15	14.99
${\beta }_{23}$	0.9	0.81	0.09	10.02	0.78	0.12	13.87
${\beta }_{25}$	1	0.89	0.11	10.84	0.85	0.15	14.84
${\beta }_{31}$	-1.3	-1.4	0.1	7.98	-1.36	0.06	4.97
${\beta }_{32}$	0.9	0.77	0.13	13.98	0.74	0.16	18.16
${\beta }_{33}$	0.8	0.72	0.08	10.01	0.69	0.11	13.78
${\beta }_{35}$	0.9	0.71	0.19	20.82	0.68	0.22	24.95
${\beta }_{41}$	-1.2	-1.32	0.12	10.17	-1.31	0.11	9.15
${\beta }_{42}$	0.8	0.63	0.17	21	0.6	0.2	25.02
${\beta }_{43}$	0.7	0.54	0.16	22.19	0.52	0.18	25.77
${\beta }_{45}$	0.8	0.64	0.16	19.57	0.61	0.19	24.17
${\beta }_{51}$	-1	-1.03	0.03	2.95	-1.04	0.04	4.31
${\beta }_{52}$	0.7	0.52	0.18	25.83	0.49	0.21	29.86
${\beta }_{53}$	0.6	0.45	0.15	25.09	0.42	0.18	29.24
${\beta }_{55}$	0.7	0.52	0.18	26.39	0.49	0.21	30.27
Covariance matrix elements
${\overline{\Sigma }}_{11}$	1.5	1.51	0.01	0.33	1.77	0.27	17.8
${\overline{\Sigma }}_{21}$	0.3	0.2	0.1	32.59	0.15	0.15	50.2
${\overline{\Sigma }}_{31}$	0.4	0.4	0	0.16	0.38	0.02	4.28
${\overline{\Sigma }}_{33}$	1.1	1.13	0.03	3.15	1.14	0.04	3.63
${\overline{\Sigma }}_{41}$	0.6	0.57	0.03	5.12	0.59	0.01	1.82
${\overline{\Sigma }}_{44}$	1.2	1.32	0.12	10	1.41	0.21	17.19
${\overline{\Sigma }}_{51}$	0.5	0.47	0.03	5.93	0.48	0.02	4.8
${\overline{\Sigma }}_{55}$	1.3	1.28	0.02	1.77	1.38	0.08	6.22
Degree-of-freedom
$\delta$	12	13.56	1.56	13.04
Average			0.14	11.36		0.12	14.46
Loglikelihood		-8977.45			-8994.84

Table 6. Table 6 : Change in probability due to 50% increase in commute distance (DGP with DOF of 2)

Alternatives (Density per square mile)	GCM-t with the estimable DOF			GCM-t with the fixed DOF of 300			t-stat for the difference in mean elasticity
	Mean	Std. Dev.	t-stat	Mean	Std. Dev.	t-stat
0-99	0.0076	0.0025	3.04	0.0129	0.0030	4.30	9.23
100-499	0.0045	0.0025	1.80	0.0037	0.0024	1.54	4.08
500 – 1499	0.0147	0.0016	9.19	0.0135	0.0030	4.50	27.24
1500-2000	-0.0037	0.0022	-1.68	-0.0135	0.0025	-5.40	12.53
2000+	-0.0232	0.0026	-8.92	-0.0165	0.0015	-11.00	60.78

Table 7. Table 7 : Sample of a choice situation in the discrete choice experiment

	Gasoline version	Electric version
Purchase Price	$19,000	$34,000
Driving cost per 50 miles	$5.5 per 50 miles	$2.5 per 50 miles
Electric driving range		250 miles
EV parking: charging time for 50 miles		90 minutes per 50 miles
EV parking: charging type		On-street
EV parking: time to find space		5 minutes
EC parking: monthly price		$50 per month

Table 8. Table 8 : Descriptive statistics of the sample

Variables	Frequency/Average
Married indicator	36.45%
Indicator for having children	49.55%
Indicator for working full time	59.47%
Indicator for holding master’s or above degree	17.70%
Number of adults in household	2.88
Number of driving-license-holders in household	2.03
Number of vehicles in household	1.70
Indicator for owning hybrid electric vehicle	4.17%
Male indicator	29.51%
Hispanic indicator	9.53%
Walk score	77.33
Population density	21.53
Household annual vehicle miles traveled	14883.60
Race
African-American indicator	24.38%
Asian indicator	3.76%
Caucasian indicator	62.06%
Age category
Baby-boomer indicator	16.41%
GenX indicator	24.58%
Millenial indicator	57.65%

Table 9. Table 9 : Empirical study: comparison of parameter estimates of GCM-t and GCM-N (t-value in parenthesis)

Model fit Statistics	GCM-t Model				GCM-N Model
Sample size	12336				12336
Number of parameters	39				42
Loglikelihood	-22725.81				-22929.19
Bayesian Information Criterion	45611.18				46030.21
Trace of covariance matrix	1.06				3.43
Parameters	Continuous	Unordered			Continuous	Unordered
	Household vehicle miles traveled	Gasoline	Electric	Opt-out	Household vehicle miles traveled	Gasoline	Electric	Opt-out
Intercept	1.983 (81.22)		-1.755 (-5.69)	-4.229 (-9.82)	1.856 (34.64)		-1.599 (-5.66)	-5.027 (-7.83)
Demographic Variables
Married indicator	-0.059 (-4.38)				-0.105 (-6.19)		-0.033 (-1.10)
Indicator for having children			-0.244 (-7.80)		-0.030 (-4.50)		-0.227 (-7.09)
Indicator for working full time	0.174 (13.38)				0.136 (9.16)
Indicator for holding master’s or above degree			0.144 (3.94)		0.025 (1.14)		0.119 (3.44)
Number of adults in household					0.109 (15.05)		0.049 (4.15)
Number of driving license holders in household	0.055 (6.86)		0.076 (3.75)		-0.033 (-2.91)
Number of vehicles in household			-0.076 (-3.78)		0.013 (1.29)		-0.068 (-4.14)
Indicator for owning hybrid electric vehicle			0.580 (8.12)				0.522 (8.02)
Male indicator							0.081 (2.85)
Hispanic indicator			-0.139 (-2.80)		-0.061 (-2.62)		-0.119 (-2.77)
Race (base: Caucasian)
African-American indicator					0.159 (9.62)		-0.071 (-2.25)
Asian indicator					-0.102 (-2.80)		0.274 (4.31)
Age category (base: Millennial)
Baby boomer indicator			-0.473 (-10.30)				-0.436 (-10.10)
GenX indicator			-0.207 (-5.51)				-0.189 (-5.60)
Environmental variables
Walk-score			0.635 (6.46)		0.103 (1.36)		0.538 (6.17)
Population density of neighborhood	0.004 (5.60)				0.003 (2.93)
Alternative specific variables
Price (in $1,000)		-0.037 (-6.41)	-0.019 (-3.48)			-0.037 (-6.77)	-0.021 (-4.13)
Operating cost per 50 miles (in $)		-0.284 (-7.16)	-0.426 (-6.97)			-0.245 (-6.57)	-0.328 (-6.18)
Log of driving range (in 100 miles)			0.604 (8.29)				0.481 (7.50)
Electric vehicle charging time (in hours)			-0.215 (-7.37)				-0.169 (-6.37)
Electric vehicle parking search time (in minutes)			-0.015 (-3.00)				-0.013 (-2.99)
Monthly electric vehicle parking cost (in $100)			-0.224 (-3.05)				-0.282 (-4.39)
Structural effect
Household vehicle miles traveled (in 1000 miles)			-0.019 (-6.77)				-0.007 (-4.85)

Table 10. Table 10 : Degree-of-freedom specification results in GCM-t model (t-value in parenthesis)

Parameter	Estimates (T-value)
Intercept	1.191 (4.67)
Male indicator	0.202 (2.25)
Married indicator	0.128 (1.46)
Indicator for holding master’s or above degree	0.182 (1.67)
Number of adults in household	-0.165 (-5.32)
Number of driving license holders in household	0.132 (2.65)
Walk score	0.372 (1.30)
Race (base: Caucasian)
African-American indicator	-0.396 (-4.64)
Asian indicator	-0.309 (-1.48)

Table 11. Table 11 : Change in Probability of alternatives due to 1% reduction in parking-cost of electric vehicle

DOF lower limit	DOF upper limit	GCM-t model			GCM-N model
		Gasoline	Electric Vehicle	opt-out	Gasoline	Electric Vehicle	opt-out
1.27	2.37	-0.0002	0.0001	0.0000	-0.0003	0.0003	0.0000
2.37	2.83	-0.0002	0.0002	0.0000	-0.0003	0.0002	0.0000
2.83	3.16	-0.0002	0.0002	0.0000	-0.0003	0.0003	0.0000
3.16	3.53	-0.0002	0.0002	0.0000	-0.0003	0.0002	0.0000
3.53	3.86	-0.0002	0.0002	0.0000	-0.0003	0.0003	0.0000
3.86	4.16	-0.0002	0.0002	0.0000	-0.0003	0.0003	0.0000
4.16	4.40	-0.0002	0.0002	0.0000	-0.0004	0.0003	0.0000
4.40	4.83	-0.0002	0.0002	0.0000	-0.0004	0.0003	0.0000
4.83	5.33	-0.0003	0.0002	0.0000	-0.0004	0.0003	0.0000
5.33	9.32	-0.0003	0.0002	0.0000	-0.0004	0.0003	0.0000

Table 12. Table 12 : Change in Probability of alternatives due to 25% reduction in parking-cost of electric vehicle

DOF lower limit	DOF upper limit	GCM-t model			GCM-N model
		Gasoline	Electric Vehicle	opt-out	Gasoline	Electric Vehicle	opt-out
1.27	2.37	-0.0050	0.0038	-0.0005	-0.0087	0.0064	-0.0004
2.37	2.83	-0.0054	0.0041	-0.0006	-0.0088	0.0061	-0.0006
2.83	3.16	-0.0052	0.0044	-0.0005	-0.0084	0.0068	-0.0005
3.16	3.53	-0.0053	0.0043	-0.0004	-0.0084	0.0062	-0.0004
3.53	3.86	-0.0056	0.0050	-0.0006	-0.0089	0.0073	-0.0005
3.86	4.16	-0.0057	0.0051	-0.0005	-0.0088	0.0070	-0.0006
4.16	4.40	-0.0062	0.0053	-0.0004	-0.0094	0.0072	-0.0006
4.40	4.83	-0.0060	0.0052	-0.0007	-0.0091	0.0070	-0.0009
4.83	5.33	-0.0072	0.0053	-0.0004	-0.0104	0.0073	-0.0005
5.33	9.32	-0.0069	0.0046	-0.0005	-0.0100	0.0064	-0.0007

Table 13. Table 13 : Ratio of GCM-t and GCM-N probabilities for chosen alternative for different DOF

DOF lower limit	DOF upper limit	Gasoline	Electric Vehicle	opt-out	Sample size
1.27	2.37	0.96	0.96	1.91	1232
2.37	2.83	0.97	1.00	1.33	1232
2.83	3.16	0.99	0.96	1.30	1232
3.16	3.53	0.98	1.00	1.18	1208
3.53	3.86	1.01	0.97	1.44	1264
3.86	4.16	1.00	1.00	0.98	1216
4.16	4.40	1.01	1.00	0.77	1248
4.40	4.83	1.02	1.00	0.87	1224
4.83	5.33	1.02	1.01	0.96	1240
5.33	9.32	1.03	1.00	0.80	1232
Overall Average		1.00	0.99	1.15

Table 14. Table 14 : Comparison of Brier scores in the five-fold cross-validation analysis

	Model	Gasoline		Electric		Opt-out
		Mean	Std. Dev.	Mean	Std. Dev.	Mean	Std. Dev.
In-Sample	GCM-T	5791.99	21.73	5480.91	22.03	2166.28	31.36
	GCM-N	5803.92	15.83	5503.46	24.74	2194.23	33.31
Out-of-Sample	GCM-T	937.68	6.53	866.83	7.87	176.09	5.74
	GCM-N	934.35	8.09	868.79	9.04	167.36	5.70

Table 15. Table 15 : Covariance matrix (t-value in parenthesis)

GCM-t model				GCM-N model
0.330 (40.42)				0.638 (84.77)
0.109 (7.12)	1.0 (fixed)			0.109 (5.42)	1.0 (fixed)
-0.020 (-1.31)	0.250 (2.35)	0.728 (3.16)		0.029 (1.07)	0.348 (1.58)	2.790 (2.60)

Equations

$$\begin{array}{rl}& \overline{\mathbf{U}}|\mathbf{y}\sim \text{mbox}{MVT}_{I_K-I}\left({\overleftrightarrow{\mathbf{B}}}_{\overline{\mathbf{U}}},{\overleftrightarrow{\mathbf{\Sigma }}}_{\overline{\mathbf{U}}},\overleftrightarrow{\delta }\right),\end{array}$$

$$\begin{array}{rl}& {\overleftrightarrow{\mathbf{B}}}_{\overline{\mathbf{U}}}={\tilde{\mathbf{B}}}_{\overline{\mathbf{U}}}+{\tilde{\mathbf{\Sigma }}}_{\mathbf{y},\overline{\mathbf{U}}}^{\top }{\left({\tilde{\mathbf{\Sigma }}}_{\mathbf{y}}\right)}^{-1}(\mathbf{y}-{\tilde{\mathbf{B}}}_{\mathbf{y}})\\ & {\overleftrightarrow{\mathbf{\Sigma }}}_{\overline{\mathbf{U}}}=\left[\frac{\delta +\alpha }{\delta +H}\right]\left({\tilde{\mathbf{\Sigma }}}_{\overline{\mathbf{U}}}-{\tilde{\mathbf{\Sigma }}}_{\mathbf{y},\overline{\mathbf{U}}}^{\top }{\left({\tilde{\mathbf{\Sigma }}}_{\mathbf{y}}\right)}^{-1}{\tilde{\mathbf{\Sigma }}}_{\mathbf{y},\overline{\mathbf{U}}}\right)\end{array}$$

$$\begin{array}{rl}\overleftrightarrow{\delta }& =\delta +H\\ \alpha & =(\mathbf{y}-{\tilde{\mathbf{B}}}_{\mathbf{y}}{)}^{\top }{\left({\tilde{\mathbf{\Sigma }}}_{\mathbf{y}}\right)}^{-1}(\mathbf{y}-{\tilde{\mathbf{B}}}_{\mathbf{y}})\end{array}$$

$$\begin{array}{rl}\mathcal{L}(\mathbf{\theta })=& \text{mbox}Pr(\mathbf{y})\text{mbox}Pr(1_k=1_m,2_k=2_m,\dots ,I_k=I_m|\mathbf{y})\\ =& f_H\left(\mathbf{y}|{\tilde{\mathbf{B}}}_{\mathbf{y}},{\tilde{\mathbf{\Sigma }}}_{\mathbf{y}},\delta \right)\left[\underset{-\infty }{\overset{{\overleftrightarrow{\mathbf{B}}}_{\overline{\mathbf{U}}}}{\int }}{f}_{(I_K-I)}\left(r|{\overleftrightarrow{\mathbf{B}}}_{\overline{\mathbf{U}}},{\overleftrightarrow{\mathbf{\Sigma }}}_{\overline{\mathbf{U}}},\overleftrightarrow{\delta }\right)d\mathbf{r}\right]\end{array}$$

$$T_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )=\underset{\mathbf{a}}{\overset{\mathbf{b}}{\int }}f(\mathbf{r}|\mathbf{\Omega },\delta )d\mathbf{r}=\frac{\Gamma \left(\frac{\delta +p}{2}\right)}{\Gamma \left(\frac{\delta }{2}\right)(\pi \delta {)}^{\frac{p}{2}}|\mathbf{\Omega }{|}^{\frac{1}{2}}}\underset{\mathbf{a}}{\overset{\mathbf{b}}{\int }}{\left(1+\frac{{\mathbf{r}}^{\top }{\mathbf{\Omega }}^{-1}\mathbf{r}}{\delta }\right)}^{-\frac{(\delta +p)}{2}}d\mathbf{r}$$

$$T_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )={\kappa }_{\delta }^p\underset{\mathbf{a}\le {\mathbf{L}}_{\mathbf{\Omega }}\mathbf{w}\le \mathbf{b}}{\int }{\left(1+\frac{{\mathbf{w}}^{\top }\mathbf{w}}{\delta }\right)}^{-\frac{(\delta +p)}{2}}d\mathbf{w}$$

$$\left(1+\frac{{\mathbf{w}}^{\top }\mathbf{w}}{\delta }\right)=\left(1+\frac{w_1^2}{\delta }\right)\left(1+\frac{w_2^2}{\delta +w_1^2}\right)\dots \left(1+\frac{w_p^2}{\delta +\sum _{j=1}^{p-1}{w}_j^2}\right)$$

$$T_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )={\kappa }_{\delta }^p\underset{{\overline{a}}_1}{\overset{{\overline{b}}_1}{\int }}{\left(1+\frac{w_1^2}{\delta }\right)}^{-\frac{(\delta +p)}{2}}\underset{{\overline{a}}_2}{\overset{{\overline{b}}_2}{\int }}{\left(1+\frac{w_2^2}{\delta +w_1^2}\right)}^{-\frac{(\delta +p)}{2}}\cdots \underset{{\overline{a}}_p}{\overset{{\overline{b}}_p}{\int }}{\left(1+\frac{w_p^2}{\delta +\sum _{j=1}^{p-1}{w}_j^2}\right)}^{-\frac{(\delta +p)}{2}}d\mathbf{w}$$

$$\begin{array}{rl}{T}_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )& ={\kappa }_{\delta }^p\sqrt{\left(\frac{\delta }{\delta +1}\frac{\delta }{\delta +2}\dots \frac{\delta }{\delta +p-1}\right)}\dots \\ & \underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}\cdots \underset{{\hat{a}}_p}{\overset{{\hat{b}}_p}{\int }}{\left(1+\frac{u_p^2}{\delta +p-1}\right)}^{-\frac{(\delta +p)}{2}}d\mathbf{u}\end{array}$$

$$\begin{array}{rl}{T}_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )& =\left[{\kappa }_{\delta +1-1}^1\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}d{u}_1\right]\left[{\kappa }_{\delta +2-1}^1\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_1\right]\\ & \dots \left[{\kappa }_{\delta +p-1}^1\underset{{\hat{a}}_p}{\overset{{\hat{b}}_p}{\int }}{\left(1+\frac{u_p^2}{\delta +p-1}\right)}^{-\frac{(\delta +p)}{2}}d{u}_p\right]\end{array}$$

$$T_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )=\underset{d_1}{\overset{e_1}{\int }}\underset{d_2}{\overset{e_2}{\int }}\underset{d_3}{\overset{e_3}{\int }}\cdots \underset{d_p}{\overset{e_p}{\int }}d{z}_{1}d{z}_{2}d{z}_3\dots d{z}_p$$

$$T_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )=(e_1-d_1)\underset{0}{\overset{1}{\int }}(e_2-d_2)\cdots \underset{0}{\overset{1}{\int }}(e_p-d_p)\underset{0}{\overset{1}{\int }}d\mathbf{\varphi }=\underset{(p-1)}{\underbrace{\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1}{\int }}\underset{0}{\overset{1}{\int }}\cdots \underset{0}{\overset{1}{\int }}}}f(\mathbf{\varphi })d\mathbf{\varphi }$$

$${\mathcal{L}}_{CML}(\mathbf{\theta })=f_H\left(\mathbf{y}|{\tilde{\mathbf{B}}}_{\mathbf{y}},{\tilde{\mathbf{\Sigma }}}_{\mathbf{y}},\delta \right)\left(\prod _{i=1}^{I-1}\prod _{j=i+1}^I\text{mbox}Pr(i_k=i_m,j_k=j_m|\mathbf{y})\right)$$

$$\begin{split}\bm{Y}\bm{U}&=\left[\begin{array}[]{c}\bm{y}\\ \hline\cr\bm{U}\end{array}\right]=\left[\begin{array}[]{c}\text{Commute distance}\\ \hline\cr\text{Density 2000+}\\ \text{Density 0-99}\\ \text{Density 100-499}\\ \text{Density 500-1499}\\ \text{Density 1500-2000}\\ \end{array}\right]\\ &=\left[\begin{array}[]{ccccc}1.00&0.50&0.75&-0.50&0.00\\ 0.00&0.00&0.00&0.00&0.00\\ -1.50&1.00&0.90&0.00&1.00\\ -1.30&0.90&0.80&0.00&0.90\\ -1.20&0.80&0.70&0.00&0.80\\ -1.00&0.70&0.60&0.00&0.70\\ \end{array}\right]\left[\begin{array}[]{c}\text{Constant}\\ \text{High-income household}\\ \text{Household with children}\\ \text{Household with a bachelor's degree holder}\\ \text{Commute distance}\\ \end{array}\right]\end{split}$$

$$\overline{\bm{\Sigma}}=\left[\begin{array}[]{c|cccc}\textbf{1.50}&\textbf{0.30}&\textbf{0.40}&\textbf{0.60}&\textbf{0.50}\\ \hline\cr\textbf{0.30}&1.00&0.50&0.50&0.50\\ \textbf{0.40}&0.50&\textbf{1.10}&0.50&0.50\\ \textbf{0.60}&0.50&0.50&\textbf{1.20}&0.50\\ \textbf{0.50}&0.50&0.50&0.50&\textbf{1.30}\\ \end{array}\right]$$

$$\begin{split}\text{Variance-covariance of error difference}\quad\bm{\overline{\Lambda}}&=\left[\begin{array}[]{c c c c}\bm{\overline{\Lambda}}_{1}&\bm{\overline{\Lambda}}_{1,2}&\dots&\bm{\overline{\Lambda}}_{1,I}\\ \bm{\overline{\Lambda}}_{2,1}&\bm{\overline{\Lambda}}_{2}&\dots&\bm{\overline{\Lambda}}_{2,I}\\ \vdots&\vdots&\ddots&\vdots\\ \bm{\overline{\Lambda}}_{I,1}&\bm{\overline{\Lambda}}_{I,2}&\dots&\bm{\overline{\Lambda}}_{I}\end{array}\right]_{(I_{K}-I)\times(I_{K}-I)}\\ \text{where}\quad\bm{\overline{\Lambda}}_{i}&=\left[\begin{array}[]{c c c c}1&\overline{\Lambda}_{1,2}^{i}&\dots&\overline{\Lambda}_{1,i_{K}-1}^{i}\\ \overline{\Lambda}_{2,1}^{i}&\overline{\Lambda}_{2,2}^{i}&\dots&\overline{\Lambda}_{2,i_{K}-1}^{i}\\ \vdots&\vdots&\ddots&\vdots\\ \overline{\Lambda}_{i_{K}-1,1}^{i}&\overline{\Lambda}_{i_{K}-1,2}^{i}&\dots&\overline{\Lambda}_{i_{K}-1,i_{K}-1}^{i}\\ \end{array}\right]_{(i_{K}-1)\times(i_{K}-1)}\\ \text{Variance-covariance of error}\quad\bm{\Lambda}&=\left[\begin{array}[]{c c c c}\bm{\Lambda}_{1}&\bm{\Lambda}_{1,2}&\dots&\bm{\Lambda}_{1,I}\\ \bm{\Lambda}_{2,1}&\bm{\Lambda}_{2}&\dots&\bm{\Lambda}_{2,I}\\ \vdots&\vdots&\ddots&\vdots\\ \bm{\Lambda}_{I,1}&\bm{\Lambda}_{I,2}&\dots&\bm{\Lambda}_{I}\end{array}\right]_{(I_{K})\times(I_{K})}\\ where\quad\bm{\Lambda}_{i}&=\left[\begin{array}[]{c c}0&\bm{0}_{1,i_{K}-1}\\ \bm{0}_{i_{K}-1,1}&\bm{\overline{\Lambda}}_{i}\end{array}\right]_{i_{K}\times i_{K}}\\ and\quad\bm{\Lambda}_{i,j}&=\left[\begin{array}[]{c c}0&\bm{0}_{1,j_{K}-1}\\ \bm{0}_{i_{K}-1,1}&\bm{\overline{\Lambda}}_{i,j}\end{array}\right]_{i_{K}\times j_{K}}\end{split}$$

$$\begin{split}\bm{\overline{\Lambda}}&=\left[\begin{array}[]{c c}\bm{\overline{\Lambda}}_{1}&\bm{\overline{\Lambda}}_{1,2}\\ \bm{\overline{\Lambda}}_{2,1}&\bm{\overline{\Lambda}}_{2}\end{array}\right]=\left[\begin{array}[]{c c|c c c}1&\overline{\Lambda}_{1,2}^{1}&\overline{\Lambda}_{1,1}^{1,2}&\overline{\Lambda}_{1,2}^{1,2}&\overline{\Lambda}_{1,3}^{1,2}\\ \overline{\Lambda}_{2,1}^{1}&\overline{\Lambda}_{2,2}^{1}&\overline{\Lambda}_{2,1}^{1,2}&\overline{\Lambda}_{2,2}^{1,2}&\overline{\Lambda}_{2,3}^{1,2}\\ \hline\cr\overline{\Lambda}_{1,1}^{2,1}&\overline{\Lambda}_{1,2}^{2,1}&1&\overline{\Lambda}_{1,2}^{2}&\overline{\Lambda}_{1,3}^{2}\\ \overline{\Lambda}_{2,1}^{2,1}&\overline{\Lambda}_{2,2}^{2,1}&\overline{\Lambda}_{2,1}^{2}&\overline{\Lambda}_{2,2}^{2}&\overline{\Lambda}_{2,3}^{2}\\ \overline{\Lambda}_{3,1}^{2,1}&\overline{\Lambda}_{3,2}^{2,1}&\overline{\Lambda}_{3,1}^{2}&\overline{\Lambda}_{3,2}^{2}&\overline{\Lambda}_{3,3}^{2}\\ \end{array}\right]\\ \bm{\Lambda}&=\left[\begin{array}[]{c c}\bm{\Lambda}_{1}&\bm{\Lambda}_{1,2}\\ \bm{\Lambda}_{2,1}&\bm{\Lambda}_{2}\end{array}\right]==\left[\begin{array}[]{c c c|c c c c}0&0&0&0&0&0&0\\ 0&1&\overline{\Lambda}_{1,2}^{1}&0&\overline{\Lambda}_{1,1}^{1,2}&0\overline{\Lambda}_{1,2}^{1,2}&\overline{\Lambda}_{1,3}^{1,2}\\ 0&\overline{\Lambda}_{2,1}^{1}&\overline{\Lambda}_{2,2}^{1}&0&\overline{\Lambda}_{2,1}^{1,2}&\overline{\Lambda}_{2,2}^{1,2}&\overline{\Lambda}_{2,3}^{1,2}\\ \hline\cr 0&0&0&0&0&0&0\\ 0&\overline{\Lambda}_{1,1}^{2,1}&\overline{\Lambda}_{1,2}^{2,1}&0&1&\overline{\Lambda}_{1,2}^{2}&\overline{\Lambda}_{1,3}^{2}\\ 0&\overline{\Lambda}_{2,1}^{2,1}&\overline{\Lambda}_{2,2}^{2,1}&0&\overline{\Lambda}_{2,1}^{2}&\overline{\Lambda}_{2,2}^{2}&\overline{\Lambda}_{2,3}^{2}\\ 0&\overline{\Lambda}_{3,1}^{2,1}&\overline{\Lambda}_{3,2}^{2,1}&0&\overline{\Lambda}_{3,1}^{2}&\overline{\Lambda}_{3,2}^{2}&\overline{\Lambda}_{3,3}^{2}\\ \end{array}\right]\\ \bm{D}&=\left[\begin{array}[]{c c c c c}0&0&0&0&0\\ 1&0&0&0&0\\ 0&1&0&0&0\\ 0&0&0&0&0\\ 0&0&1&0&0\\ 0&0&0&1&0\\ 0&0&0&0&1\\ \end{array}\right]\end{split}$$

$$\bm{D}_{m}=\left[\begin{array}[]{c c}\bm{1}_{H\times H}&\bm{0}_{H\times(I_{K}-I)}\\ \bm{0}_{I_{K}\times H}&\bm{D}\end{array}\right]_{(H+I_{K})\times(H+I_{K}-I)}$$

$$\begin{split}\bm{D}_{m}&=\left[\begin{array}[]{c c}\bm{1}_{2\times 2}&\bm{0}_{2\times 5}\\ \bm{0}_{7\times 2}&\bm{D}_{7\times 5}\end{array}\right]\\ &=\left[\begin{array}[]{c c | c c c c c}1&0&0&0&0&0&0\\ 0&1&0&0&0&0&0\\ \hline\cr 0&0&0&0&0&0&0\\ 0&0&1&0&0&0&0\\ 0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0\\ 0&0&0&0&1&0&0\\ 0&0&0&0&0&1&0\\ 0&0&0&0&0&0&1\\ \end{array}\right]\end{split}$$

$$\bm{M}=\left[\begin{array}[]{c c | c c c c c c c}1&0&0&0&0&0&0&0&0\\ 0&1&0&0&0&0&0&0&0\\ \hline\cr 0&0&1&-1&0&0&0&0&0\\ 0&0&0&-1&1&0&0&0&0\\ 0&0&0&0&0&1&0&-1&0\\ 0&0&0&0&0&0&1&-1&0\\ 0&0&0&0&0&0&0&-1&1\\ \end{array}\right]$$

$$\begin{array}{rl}{T}_p(\mathbf{a},\mathbf{b},\mathbf{\Omega },\delta )& ={\kappa }_{\delta }^2\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +2)}{2}}\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_1\left(d{u}_2\sqrt{\frac{\delta +u_1^2}{\delta +1}}\right)\\ & ={\kappa }_{\delta }^2\left[\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}\sqrt{\frac{\delta +u_1^2}{\delta +1}}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +2)}{2}}d{u}_1\right]\left[\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_2\right]\\ & ={\kappa }_{\delta }^2\left[\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}\sqrt{\left(\frac{\delta }{\delta +1}\right)\left(1+\frac{u_1^2}{\delta }\right)}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +2)}{2}}d{u}_1\right]\left[\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_2\right]\\ & ={\kappa }_{\delta }^2\sqrt{\frac{\delta }{\delta +1}}\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_{1}d{u}_2\text{(equation refeq:MVTCD:M3)}\\ & \text{Now deriving equation refeq:MVTCD:M4}\\ & =\frac{\Gamma \left(\frac{\delta +2}{2}\right)}{\Gamma \left(\frac{\delta }{2}\right)(\pi \delta )}\sqrt{\frac{\delta }{\delta +1}}\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_{1}d{u}_2\\ & =\frac{\Gamma \left(\frac{\delta +1}{2}\right)\Gamma \left(\frac{\delta +1+1}{2}\right)}{\Gamma \left(\frac{\delta }{2}\right)(\pi \delta {)}^{\frac{1}{2}}\Gamma \left(\frac{\delta +1}{2}\right)(\pi (\delta +1){)}^{\frac{1}{2}}}\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_{1}d{u}_2\\ & ={\kappa }_{\delta }^1{\kappa }_{\delta +1}^1\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_{1}d{u}_2\\ & =\left[{\kappa }_{\delta +1-1}^1\underset{{\hat{a}}_1}{\overset{{\hat{b}}_1}{\int }}{\left(1+\frac{u_1^2}{\delta }\right)}^{-\frac{(\delta +1)}{2}}d{u}_1\right]\left[{\kappa }_{\delta +2-1}^1\underset{{\hat{a}}_2}{\overset{{\hat{b}}_2}{\int }}{\left(1+\frac{u_2^2}{\delta +2-1}\right)}^{-\frac{(\delta +2)}{2}}d{u}_1\right]\text{(equation refeq:MVTCD:M4)}\end{array}$$

Full text

**A Generalized Continuous-Multinomial Response Model with a t-distributed Error Kernel

**January 18, 2020

Subodh Dubey*

Department of Transport and Planning

Delft University of Technology, Netherlands

S.K.Dubey@tudelft.nl

Prateek Bansal*

**(Corresponding author)

**School of Civil and Environmental Engineering

Cornell University, United States

pb422@cornell.edu

Ricardo A. Daziano

School of Civil and Environmental Engineering

Cornell University, United States

daziano@cornell.edu

Erick Guerra

Department of City and Regional Planning

University of Pennsylvania, United States

erickg@upenn.edu

• These authors contributed equally to this work.

Abstract

In multinomial response models, idiosyncratic variations in the indirect utility are generally modeled using Gumbel or normal distributions. This study makes a strong case to substitute these thin-tailed distributions with a t-distribution. First, we demonstrate that a model with a t-distributed error kernel better estimates and predicts preferences, especially in class-imbalanced datasets. Our proposed specification also implicitly accounts for decision-uncertainty behavior, i.e. the degree of certainty that decision-makers hold in their choices relative to the variation in the indirect utility of any alternative. Second – after applying a t-distributed error kernel in a multinomial response model for the first time – we extend this specification to a generalized continuous-multinomial (GCM) model and derive its full-information maximum likelihood estimation procedure. The likelihood involves an open-form expression of the cumulative density function of the multivariate t-distribution, which we propose to compute using a combination of the composite marginal likelihood method and the separation-of-variables approach. Third, we establish finite sample properties of the GCM model with a t-distributed error kernel (GCM-t) and highlight its superiority over the GCM model with a normally-distributed error kernel (GCM-N) in a Monte Carlo study. Finally, we compare GCM-t and GCM-N in an empirical setting related to preferences for electric vehicles (EVs). We observe that accounting for decision-uncertainty behavior in GCM-t results in lower elasticity estimates and a higher willingness to pay for improving the EV attributes than those of the GCM-N model. These differences are relevant in making policies to expedite the adoption of EVs.

1 Introduction

Discrete-response (e.g., count, ordered, binary, multinomial) models are popular across various disciplines such as applied economics, transportation, marketing, and political science. In these models, error structures are specified at different modeling levels because the researcher does not have full information about the data generating process (DGP). For instance, the total indirect utility in additive random utility maximization (ARUM) (McFadden, 1973) based choice models is specified as the sum of a deterministic index (depending on observables) and a random error term or taste shock. Such error specifications in empirical studies are generally governed by ease of estimation rather than structural appropriateness (Vijverberg and Vijverberg, 2016). In particular, without worrying about the characteristics of the data in hand, Gumbel/extreme-value (logit) or normal (probit) error kernels are commonly used. These traditional logit and probit links have been replaced by a t-distributed error kernel, i.e. a robit link (Liu, 2004), in multi-level modelling applications to handle heavy-tailed error distributions. Moreover, a t-distribution with an estimable degree-of-freedom (DOF) actually generalizes logit and probit.111The t-distribution with about seven and a large (above thirty) DOF approximates logistic and normal distributions, respectively. However, we are not aware of the use of the robit link in modeling multinomial responses, perhaps because the benefits of the generalization are unclear and the estimation of the resulting model is cumbersome.

In this study, we present the first application of a multinomial response model with a t-distributed error kernel, multinomial robit (MNR) model henceforth. Whereas the proposed MNR model retains all merits of the multinomial probit (MNP) model222Similar to MNP, MNR allows for flexible substitution patterns – correlations across indirect utilities of alternatives – without necessarily including variation in parameters across decision-makers., we illustrate three additional advantages of adopting MNR in practice. First, we highlight that the robit link is superior to the probit link in estimating and predicting preferences in unbalanced datasets where one or more alternatives have small shares. Second, we parameterize the DOF of the t-distributed error kernel as a function of demographics and show how this specification can help in capturing decision uncertainty of decision-makers that standard compensatory ARUM models cannot account for. Third, if the error distribution in the true DGP is robit (symmetric heavy-tailed), probit fails to recover the true model parameters. We numerically establish this benefit of MNR over MNP in a Monte Carlo study. Given the growing interest in the joint modeling of mixed datasets across multiple disciplines (see De Leon and Chough, 2013, for applications), we further extend MNR to a generalized continuous-multinomial (GCM) response model with a t-distributed error kernel that facilitates simultaneous consideration of multiple multinomial and multiple continuous dependent variables. We note that the advantages of considering a t-distributed (over a normally-distributed) error kernel are even more evident in the GCM response model because heavy-tailed distributions are more commonly observed in continuous outcomes. Moreover, the joint GCM model offers several advantages, namely: a) statistically efficient estimation, b) easier hypothesis testing and better power of statistical tests, c) avoidance of inconsistencies in a situation when continuous and multinomial outcomes have structural dependence or have common unobserved factors (see Bhat et al., 2015, for a detailed discussion). For example, the GCM response model is appropriate for the joint modeling of commute distance (continuous variable) and choice of residential location (multinominal variable) because both outcomes may affect each other and may share common unobserved factors. Joint modeling is thus desirable, but it is generally intractable due to an absence of convenient distributions to represent the conditional and/or joint relationship between the outcomes. The proposed GCM response model with a t-distributed error kernel (henceforth, GCM-t) could not have been estimated conveniently without an elegant statistical property of the t-distribution: the conditional distribution of a joint/multivariate t-distribution is also a t-distribution (Ding, 2016).

The contribution of this study is thus threefold. First, we illustrate the importance of t-distributed error kernels in multinomial choice modeling. Second, we derive a full-information maximum likelihood procedure to estimate MNR and GCM-t models. The likelihood expressions of both models involve evaluation of high-dimensional multi-variate-t-cumulative density (MVTCD) functions. We adopt the composite marginal likelihood (or paired-likelihood) approach to first decompose the multi-dimensional MVTNCD integral into multiple pairs and thus reduce the integral dimensionality (see Varin et al., 2011; Xu and Reid, 2011). Subsequently, similar to the Geweke-Hajivassiliou-Keane (GHK) method to simulate multivariate normal cumulative density functions (Genz, 1992; Hajivassiliou et al., 1996), we use a separation-of-variables approach to simulate MVTCD functions (Genz and Bretz, 1999). Third, we numerically verify the statistical properties of the maximum likelihood estimator of the GCM-t model in a Monte Carlo study. We also compare the performance of GCM-t and GCM with normally-distributed error kernel (GCM-N) models in a second simulation study. Finally, we validate the simulation results and highlight the advantages of the robit link (over probit) in an empirical study in the context of policies for on-street parking with charging facilities for electric vehicles. The empirical data looks into the adoption of electric vehicles and the outcomes of interest are a household’s vehicle miles traveled (continuous) and vehicle-purchase preferences (multinomial).

The remainder of this paper is organized as follows. We present the contextual literature review in section 2; section 3 details specification of the GCM-t model and derives its estimator; section 4 illustrates advantages of adopting choice models with t-distributed error kernels in practice; section 5 presents a comprehensive Monte Carlo study and highlights the benefits of GCM-t over GCM-N; section 6 validates the simulation-based findings in the empirical study; and, finally, section 7 concludes and discusses avenues of future research.

2 Literature Review

Regarding flexibility in modeling limited dependent variables, there is an extensive literature covering semi-parametric error distributions for binary- and multinomial-response models. However, we restrict our discussion to flexible, parsimonious parametric error specifications, and highlight research gaps that this study addresses.

To overcome constraints from the assumption of symmetric and thin-tailed error kernels in logit and probit models, existing research offers several alternative error specifications with one or two additional shape parameters. In statistical terms, whereas an error kernel with one additional parameter (e.g., t-distribution and skew normal) allows for the trade-off between skewness and kurtosis, two additional parameters (e.g., skew t-distribution) accommodate several choices of both skewness and kurtosis.

The scobit model derived from a Burr-10 distribution (Nagler, 1994), the robit model which considers a t-distribution with estimable DOF (Liu, 2004), and the skew-probit model which assumes a skew-normal distribution with estimable skew parameter (Bazán et al., 2010) are among well-known binary response models that use error kernels with one additional parameter. Error kernels with two additional parameters, namely skewed t-distribution (Kim et al., 2007) and the generalized Tukey lambda (GTL) family of distributions (Vijverberg and Vijverberg, 2016) have also been applied for binary outcomes.

Some of these binary-response models with a flexible error structure have been extended to the multinomial dependent variables. Castillo et al. (2008) and Fosgerau and Bierlaire (2009) derived a multinomial choice model considering a Weibull distribution on the error term in a specification that is also known as the weibit model. Li (2011) proposed a generalized method to construct asymmetric multinomial choice models for a family of error distributions with heteroskedastic variance, which also nests the weibit and logit models. Recently, Nakayama and Chikaraishi (2015) derived a unified multinomial choice model using the q-generalized-extreme-value (q-GEV) distribution with an estimable shape parameter and applied their model to transportation network assignment problems. Brathwaite and Walker (2018) identified that all these flexible multinomial models impose restrictions on the magnitude and/or sign of the index function. To address this concern, Brathwaite and Walker proposed a generalized link function that eliminates the need for such restrictions.

Note that all the multinomial choice models overviewed above have tractable closed-form choice probability expressions but suffer from a major limitation – the joint modeling of multiple types of dependent variables (e.g., continuous, ordinal, and count) and the inclusion of spatial and social dependencies in these models are computationally intractable, if not impossible, due to an increase in the dimensionality of integration (see Guevara et al., 2009, for a discussion on the curse of dimensionality in choice models).

Error kernels with skew-normal, t, or skew-t distributions have appropriate statistical properties for such joint structural modeling, but the use of these flexible distributions in multinomial choice models is cumbersome due to open-form choice probability expressions as well as noted inference issues.333Models with skew-normal and skew-t error kernel can encounter inference problems due to eventual singularity of the Fisher information matrix (when direct parameterization is used) and violation of asymptotic theory for centered parameterization (see Pewsey, 2000; Azzalini and Arellano-Valle, 2013, for a discussion on these issues).

The skew-normal distribution has been used at a few instances in the mixed MNP model but not to specify the error kernel, rather to model unobserved preference heterogeneity (random parameter choice models as in Bhat and Sidharthan, 2012). Since random taste variations and error term heterogeneity are confounded (Brownstone and Train, 1998), the mixed MNP specification can be viewed as an indirect utility with a non-random index function and a skew-normal error term. Bhat et al. (2015) and Bhat et al. (2017) have also used the skew-normal distribution to account for non-normality in latent constructs within an integrated choice and latent variable model (Ben-Akiva et al., 2002) and in the error kernel of an ordered response model, respectively.

However, we are not aware of the application of t- or skew-t-distributed error kernels in multinomial choice models, which are commonly used in binary, linear mixed or multilevel, and censored linear regressions (Pinheiro et al., 2001; Liu, 2004; Koenker and Yoon, 2009; Marchenko and Genton, 2012; Wang et al., 2018) due to their ability to model heavy-tailed distributions (Figure 1). This study particularly contributes to the literature by first illustrating the statistical and behavioral implications of using a heavy-tailed error kernel in multinomial-response models and proposing the first multinomial choice model with the t-distributed error kernel.444The skew-t kernel could additionally account for asymmetry of the error distribution, but we do not consider it because inference in such models would have similar issues as we encounter in models with skew-normal kernel (Azzalini and Arellano-Valle, 2013).

3 Methodology

In this section, we first discuss separate models for multiple-continuous and multiple-discrete responses with t-distributed error kernels. We then combine these models to derive a generalized continuous-multinomial response model with a t-distributed error kernel (GCM-t) and outline steps for implementation of its full-information maximum likelihood estimator.555Since multinomial robit (MNR) model is a special case of the GCM-t model, we do not explicitly discuss its estimator. In the proposed model formulation, we do not consider network, social, or spatial effects and therefore the utility of an individual is independent of other individuals in the sample. Thus, without loss of generality, we derive the model and estimation procedure for a single individual.

3.1 Continuous variable model

Consider a standard linear regression setup: $y_{h}=\bm{\gamma_{h}}^{\bm{\top}}\bm{X}_{h}+\xi_{h},$ where $h$ is the index of a continuous outcome $h=\{1,2,\dots,H\}$, $y_{h}$ and $\xi_{h}$ are the corresponding dependent variable and a t-distributed error term with DOF $\delta$, and $\bm{\gamma}_{h}$ and $\bm{X}_{h}$ are respectively $(s\times 1)$ vectors of coefficients and exogenous regressors.

Rewriting the regression equation in matrix form leads to: $\bm{y}=\mbox{diag}(\bm{\gamma}\bm{X}^{\bm{\top}})+\bm{\xi},$ where $(\bm{y})_{H\times 1}=[y_{1},y_{2},\dots,y_{H}]^{\top}$ is a vector of continuous outcomes and $(\bm{\xi})_{H\times 1}=[{\xi}_{1},{\xi}_{2},\dots,{\xi}_{H}]^{\top}$ is a t-distributed error vector with DOF $\delta$, and $(\bm{\gamma})_{H\times s}=\left[\bm{\gamma_{1}},\bm{\gamma_{2}},\dots,\bm{\gamma_{H}}\right]^{\top}$ and $(\bm{X})_{H\times s}=\left[\bm{X_{1}},\bm{X_{2}},\dots,\bm{X_{H}}\right]^{\top}$ are matrices of coefficients and exogenous variables. Note that $\bm{y}\sim\mbox{MVT}_{H}\left[\mbox{diag}(\bm{\gamma}\bm{X}^{\bm{\top}}),\bm{\Xi},\delta\right]$ where $\bm{\Xi}$ is the variance-covariance matrix of $\bm{\xi}$.

We consider the same DOF across all elements of the error vector because the exact distribution of linear or non-linear combinations of two t-distributed random variables with arbitrary DOF values is not known (Ahsanullah et al., 2014). Jones (2002) allows marginal distributions to have an arbitrary DOF in the case of a bivariate t-distributed random variable, but its extension to the multivariate case is not straightforward.

3.2 Choice model

Let $i$ be the index for a nominal outcome $i\in\{1,2,\dots,I\}$, and $k$ be the index of alternatives in each nominal outcome $k\in\{1,2,\dots,i_{K}\}$. Then, we can write the indirect utility of alternative $k$ in the $i^{th}$ nominal variable as $U_{ik}=\bm{\beta}_{ik}^{\bm{\top}}\bm{z}_{ik}+\epsilon_{ik},$ where $\bm{z}_{ik}$ and $\bm{\beta}_{ik}$ are $(g\times 1)$ vectors of exogenous variables and coefficients, and $\epsilon_{ik}$ is a t-distributed error term with DOF $\delta$.

If we define the total number of alternatives $I_{K}=\sum\limits_{i=1}^{I}i_{K}$ , the indirect utility vector $(\bm{U})_{I_{K}\times 1}=\left[\bm{U}_{1},\bm{U}_{2},\dots,\bm{U}_{I}\right]^{\bm{\top}}$ where $\bm{U}_{i}=\left[U_{i1},U_{i2},\dots,U_{ii_{K}}\right]$, the coefficient matrix $(\bm{\beta})_{I_{K}\times g}=\left[\bm{\beta}_{11},\bm{\beta}_{12},\dots,\bm{\beta}_{11_{K}},\dots,\bm{\beta}_{II_{K}}\right]^{\top}$, the exogenuous variable matrix $(\bm{z})_{I_{K}\times g}=\left[\bm{z}_{11},\bm{z}_{12},\dots,\bm{z}_{11_{K}},\dots,\bm{z}_{II_{K}}\right]^{\top}$, and the error vector $(\bm{\epsilon})_{I_{K}\times 1}=\left[\bm{\epsilon}_{1},\bm{\epsilon}_{2},\dots,\bm{\epsilon}_{I}\right]^{\bm{\top}}$ where $\bm{\epsilon}_{i}=\left[\epsilon_{i1},\epsilon_{i2},\dots,\epsilon_{ii_{K}}\right]$, we can write the distribution of the indirect utility as $\bm{U}\sim\mbox{MVT}_{I_{K}}\left[\mbox{diag}(\bm{\beta}\bm{z}^{\bm{\top}}),\bm{\Lambda},\delta\right]$ where $\bm{\Lambda}$ is the variance-covariance matrix of $\bm{\epsilon}$.

Since only differences in utility matter, the difference of error terms $(\overline{\bm{\epsilon}})$ is identifiable after fixing the scale of utility. In fact, we normalize the top diagonal element of the covariance matrix of error differences to 1 to fix scale of utility. We create a transformation matrix $(\bm{D})$ to convert the normalized variance-covariance of error differences $(\overline{\bm{\Lambda}}_{(I_{K}-I)\times(I_{K}-I)})$ into the undifferenced error variance-covariance matrix $(\bm{\Lambda}_{I_{K}\times I_{K}})$ using $\bm{\Lambda}=\bm{D}\overline{\bm{\Lambda}}\bm{D}^{\bm{\top}}$. We provide details of creating the transformation matrix $(\bm{D})$ and an illustration of this operatoor in appendix A.1 (cf. Bhat and Sidharthan, 2012). The indirect utility can thus be written as $\bm{U}\sim\mbox{MVT}_{I_{K}}\left[\mbox{diag}(\bm{\beta}\bm{z}^{\bm{\top}}),\bm{D}\overline{\bm{\Lambda}}\bm{D}^{\bm{\top}},\delta\right]$.666The matrix $\overline{\bm{\Lambda}}$ can be block-diagonal and still the dependencies across alternatives and nominal variables are parsimoniously generated by a single DOF parameter.

3.3 Joint Model Specification

In order to write a joint model of the continuous and nominal variable, we define \bm{YU}=\left[\begin{array}[]{c}\bm{y}\\ \bm{U}\end{array}\right]. Thus, the distribution of $\bm{YU}\sim\mbox{MVT}_{H+I_{K}}\left[\bm{B},\bm{\Sigma},\delta\right]$ where \bm{B}=\left[\begin{array}[]{c}\mbox{diag}(\bm{\gamma}\bm{X}^{\bm{\top}})\\ \mbox{diag}(\bm{\beta}\bm{z}^{\bm{\top}})\end{array}\right] and \bm{\Sigma}=\left[\begin{array}[]{c c}\bm{\Xi}&\text{Cov}(\bm{\xi},\bm{\epsilon})\\ \text{Cov}(\bm{\epsilon},\bm{\xi})&\bm{\Lambda}\end{array}\right]. If $\overline{\bm{\Sigma}}$ is the normalized (up to scale) covariance matrix of the joint differenced error \left[\begin{array}[]{c}\bm{\xi}\\ \overline{\bm{\epsilon}}\end{array}\right], which is identified, the undifferenced full variance-covariance matrix $(\bm{\Sigma})$ can be obtained from $\overline{\bm{\Sigma}}$ using the modified transformation matrix $\bm{D}_{m}$ as follows: $\bm{\Sigma}=\bm{D}_{m}\overline{\bm{\Sigma}}\bm{D}_{m}^{\bm{\top}}$. Appendix A.2 provides the details of creating $\bm{D}_{m}$, together with an example.

3.4 Joint Model Estimation

Similar to MNP estimation, we work with utility differences using the chosen alternative as base. To perform this operation, we construct the utility difference operator $\bm{M}$ of size $(H+I_{K}-I)\times(H+I_{K})$ using the algorithm given in appendix A.3. We transform the original mean and the variance-covariance matrix using $\bm{M}$, and thus derive the distribution of the joint variable $\bm{\widetilde{YU}}$ in utility-differences $(\bm{\overline{U}})$ space. We obtain $\bm{\widetilde{YU}}\sim\mbox{MVT}_{H+I_{K}-I}\left(\bm{\widetilde{B}},\bm{\widetilde{\Sigma}},\delta\right)$, where $\bm{\widetilde{B}}=\bm{M}\bm{B}$ and $\bm{\widetilde{\Sigma}}=\bm{M}\bm{\Sigma}\bm{M}^{\bm{\top}}$.

Consider the partition of $\bm{\widetilde{B}}$ and $\bm{\widetilde{\Sigma}}$ into the continuous and choice (discrete) model (in utility differences) as follows: \bm{\widetilde{B}}=\left[\begin{array}[]{c}\bm{\widetilde{B}}_{\bm{y}}\\ \bm{\widetilde{B}}_{\bm{\overline{U}}}\end{array}\right]_{(H+I_{K}-I)\times 1} and \bm{\widetilde{\Sigma}}=\left[\begin{array}[]{c c}\bm{\widetilde{\Sigma}}_{\bm{y}}&\bm{\widetilde{\Sigma}}_{\bm{y},\bm{\overline{U}}}\\ \bm{\widetilde{\Sigma}}_{\bm{\overline{U}},\bm{y}}&\bm{\widetilde{\Sigma}}_{\bm{\overline{U}}}\end{array}\right]_{(H+I_{K}-I)\times(H+I_{K}-I)}.

The conditional distribution of the utility difference vector is also t-distributed (Ding, 2016):

[TABLE]

where

[TABLE]

[TABLE]

Thus, the joint likelihood can be written as :

[TABLE]

where $i_{m}$ is the chosen alternative corresponding to the $i^{th}$ nominal variable, $f_{H}$ is the probability density function of the $H$-variate t-distribution, and $\bm{\theta}$ is a vector of identified parameters $\{\bm{\gamma},\bm{\beta},\mbox{vec}(\bm{\overline{\Sigma}}),\delta\}$.777Note that $\bm{\overline{\Sigma}}$, the normalized variance-covariance matrix of the joint differenced error, is identified and other matrices ($\bm{\Sigma}$ and $\bm{\widetilde{\Sigma}}$) are derived from it using $\bm{D}_{m}$ and $\bm{M}$ matrices. Moreover, $\mbox{vec}(\bm{\overline{\Sigma}})$ vectorizes the unique element of a matrix $\bm{\overline{\Sigma}}$.

In order to ensure positive definiteness of the normalized covariance matrix of error difference ($\bm{\overline{\Sigma}}$), we work with its Cholesky decomposition in estimation. We effectively consider the normalized Cholesky factorization such that the top diagonal element of the error differenced covariance matrix of every nominal variable ($\overline{\bm{\Lambda}}_{i}$) is fixed to 1. Details of normalization can be found in appendix A.4.

The full-information maximum likelihood estimator requires maximization of the joint likelihood function. All numerical maximization routines require to evaluate the likelihood function at a given parameter vector, which we illustrate in Algorithm 1. This computation involves evaluation of a $H-$dimensional multi-variate-t-probability-density (MVTPD) function and a $(I_{K}-I)-$dimensional multi-variate-t-cumulative-density (MVTCD) function. The MVTCD function does not have a closed-form expression and thus requires the use of simulation-aided inference. We use a separation-of-variables (SoV) approach to compute the MVTCD function (Genz and Bretz, 1999), which is detailed in section 3.4.1. Similar to other simulation-based function evaluation procedures, the SOV approach also suffers from the course of dimensionality coming from the increase in integral dimensionality. As a result, simulation-based evaluation of the function not only loses accuracy, but computation time also becomes unmanageable (Bhat, 2003; Craig, 2008). To reduce the dimension of this integration, we adopt the composite marginal likelihood (CML) method, which we briefly discuss in section 3.4.2.

3.4.1 Separation-of-Variables Approach to evaluate MVTCD Function

The underlined concept of this simulator is similar to the GHK simulator for the MVNCD function (Genz, 1992; Hajivassiliou et al., 1996). In the SOV approach, a $p-$dimensional integral is decomposed into $(p-1)$ unidimensional integrals using the Cholesky decomposition of the covariance matrix. These independent unidimensional integrals are sequentially evaluated based on the realization of all previous integrals. The approach presented below is adopted from Genz and Bretz (1999).

We present this approach for computing the area under the curve888The area under the probability density function between $\bm{a}$ and $\bm{b}$ is the same as the CDF evaluated at point $\bm{b}$ when $\bm{a}$ is $-\bm{\infty}$. of a probability density function of the p-dimensional t-distributed random variable ($\bm{r}\sim\mbox{MVT}_{p}(\bm{0},\bm{\Omega},\delta)$)999We set location parameter to $\bm{0}$ without loss of generality because if $\bm{g}\sim\mbox{MVT}_{p}(\bm{\mu},\bm{\Omega},\delta)$ and $\bm{r}\sim\mbox{MVT}_{p}(\bm{0},\bm{\Omega},\delta)$, then $\bm{g}=\bm{r}+\bm{\mu}$. between $\bm{a}$ and $\bm{b}$:

[TABLE]

Next, $\bm{\Omega}=\bm{L}_{\bm{\Omega}}\bm{L}_{\bm{\Omega}}^{\bm{\top}}$, where $\bm{L}_{\bm{\Omega}}$ is the lower triangular Cholesky factor. Then, by change of variable, $\bm{r}=\bm{L}_{\bm{\Omega}}\bm{w}$ and $\bm{r}^{\bm{\top}}\bm{\Omega}^{-1}\bm{r}=\bm{w}^{\bm{\top}}\bm{w}$, $d\bm{r}=|\bm{L}_{\bm{\Omega}}|d\bm{w}$, $d\bm{r}=|\bm{\Omega}|^{\frac{1}{2}}$. Let $\kappa_{\delta}^{p}=\frac{\Gamma\left(\frac{\delta+p}{2}\right)}{\Gamma\left(\frac{\delta}{2}\right)(\pi\delta)^{\frac{p}{2}}}$ such that equation 3 can be rewritten as:

[TABLE]

where $\bm{a}\leq\bm{L}_{\bm{\Omega}}\bm{w}\leq\bm{b}$ is the same as \left[\begin{array}[]{c}a_{1}\\ a_{2}\\ \vdots\\ a_{p}\end{array}\right]\leq\left[\begin{array}[]{c c c c}L_{11}&0&0&0\\ L_{21}&L_{22}&0&0\\ \vdots&\vdots&\ddots&0\\ L_{p1}&L_{p2}&L_{p3}&L_{pp}\end{array}\right]\left[\begin{array}[]{c}w_{1}\\ w_{2}\\ \vdots\\ w_{p}\end{array}\right]\leq\left[\begin{array}[]{c}b_{1}\\ b_{2}\\ \vdots\\ b_{p}\end{array}\right]

Let $\overline{a}_{i}=\frac{\left(a_{i}-\sum\limits_{j=1}^{i-1}L_{ij}w_{i}\right)}{L_{ii}}$ and $\overline{b}_{i}=\frac{\left(b_{i}-\sum\limits_{j=1}^{i-1}L_{ij}w_{i}\right)}{L_{ii}}$. Also note that

[TABLE]

Thus, equation 4 can be rewritten as :

[TABLE]

Consider $w_{i}=u_{i}\sqrt{\frac{\delta+\sum\limits_{j=1}^{i-1}w_{j}^{2}}{\delta+i-1}}$. We can rewrite equation 6 by substituting $w_{i}$ as follows:

[TABLE]

where $\hat{a}_{i}=\overline{a}_{i}\sqrt{\frac{\delta+i-1}{\delta+\sum\limits_{j=1}^{i-1}w_{j}^{2}}}$ and $\hat{b}_{i}=\overline{b}_{i}\sqrt{\frac{\delta+i-1}{\delta+\sum\limits_{j=1}^{i-1}w_{j}^{2}}}$. Further, equation 7 can be rewritten as

[TABLE]

The derivation of equation 8 from equation 6 is illustrated for a bivariate t-cumulative distribution function $(p=2)$ in appendix B. Next, we substitute $u_{i}=t_{\delta+i-1}^{-1}(z_{i})$ where $t_{\delta+i-1}(z_{i})=\kappa_{\delta+i-1}^{1}\int\limits_{-\infty}^{z_{i}}\left(1+\frac{s^{2}}{\delta+i-1}\right)^{-\frac{\delta+i}{2}}ds$ is the cumulative density function (CDF) of the univeriate t-distribution with DOF $\delta+i-1$, and thus $dz_{i}=k_{\delta+i-1}^{1}\left[1+\frac{u_{i}^{2}}{\delta+i-1}\right]^{-\frac{\delta+i}{2}}du_{i}$.

Finally, equation 8 can thus be rewritten as:

[TABLE]

where $d_{i}$ and $e_{i}$ are CDF of t-distributed random variable with DOF $\delta+i-1$ at points $\hat{a}_{i}$ and $\hat{b}_{i}$, respectively. After the change of variables $z_{i}=d_{i}+\phi_{i}(e_{i}-d_{i})$, equation 9 becomes:

[TABLE]

which is an integral of $f(\bm{\phi})=\prod\limits_{i=1}^{p}(e_{i}-d_{i})$ over the $(p-1)$-dimensional unit hyper-cube. This integral can be evaluated using different Quasi-Monte-Carlo or randomized lattice rule methods.

3.4.2 Composite Marginal Likelihood Approach

Computation of the joint likelihood function in equation 2 requires the evaluation of a high-dimensional integral. The dimension of this integral grows with the number of alternatives per nominal variable and also with the number of nominal variables. For example, if there are ten nominal variables, each with six alternatives, the integral would be fifty-dimensional.

Rather than directly evaluating such high-dimensional integrals, we use the CML approach (also known as paired-likelihood approach) for simplification. The CML method breaks down the joint likelihood function into multiple pairs, decreasing the dimension of the integral. More specifically, if choices made by an individual across all nominal variables is an event, this event is represented as pairwise observations in CML:

[TABLE]

where $i_{m}$ represents the chosen alternative for the $i^{th}$ nominal variable. In the CML expression above, the first term is the same as the MVTPD function, but the second term corresponds to the pairing between nominal variables with the highest dimension of integration being equal to $2[max(i_{K}\forall i)]$. Thus, by employing the CML approach, the dimension of integration in the above example can be reduced to ten from fifty. Of course, the CML approach is only applicable if there are three or more nominal variables.

A comprehensive discussion on the CML approach is outside the scope of this paper. Readers can refer to Varin and Vidoni (2005) and Varin et al. (2011) for a detailed discussion on the CML approach and Bhat (2014) for its derivation and application in the context of discrete choice models. Apart from well-established asymptotic properties (Bhat, 2014), Bhat and co-authors have tested finite sample properties of CML by applying it to complex econometric models and have observed satisfactory results (Bhat and Sidharthan, 2012; Bhat et al., 2015, 2017).

4 Implications of using GCM-t in practice

4.1 Class imbalance

Class imbalance, a very high market share of few alternatives relative to others, is often encountered in choice modeling applications such as residential location choice, travel mode choice, and credit card ownership. We take an example of binary-response data to illustrate the importance of using the robit link in datasets with class imbalance. Consider a scenario where a commuter chooses an alternative between car ($C$) and bicycle ($B$). We use a data generating process with a much higher share of car than bicycle, which is reflected in the index functions: $V_{C}=0.5HI-0.2S$ and $V_{B}=-1-0.9HI+0.5S$, where $HI$ and $S$ are indicators for high income and student commuters, respectively.

In these class imbalance situations, accurate prediction of choices is challenging. A good model should ideally predict a higher probability of choosing bicycle (the alternative with a very low market share) when a commuter actually chooses bicycle. We compare predicted choice probabilities for all four demographic groups under normally-distributed and t-distributed error kernels with varying DOFs (Table 1). The predicted probabilities of choosing bicycle by a high-income non-student commuter using probit and robit (with 0.1 DOF) links are 0.01 and 0.38, respectively. These values are 0.38 and 0.46 for a low-income student commuter. Clearly, improvement in the predicted probability of bicycle using robit is higher when the difference between the index function values of the alternatives is higher.101010The difference between the index function values of alternatives is the highest for a high-income non-student commuter and is the lowest for a low-income student commuter.

In sum, the fat-tailed nature of the t-distributed kernel is able to better predict the probability of choosing bicycle for all demographic groups, but the extent of the improvement over probit is higher if the difference in index values lies more toward the tails of the error distribution. Thus, t-distributed error kernels increase the likelihood of predicting correct travel mode assignment in such imbalanced datasets.

[Insert Table 1 here]

4.2 Behavioral implications

Intuitively, decision rules can be viewed as a mapping from attributes of alternatives to observed choices, which is governed by a latent construct. The latent construct includes an index function (or deterministic utility) and an idiosyncratic error term. Adoption of different decision rules manifests into different latent constructs and, thus, different mappings. Most of the previous studies have modeled different decision rules by modifying the index function. In fact, it turns out that changing the linear index function to non-linear functions in the latent construct of the fully-compensatory RUM-based model can translate it into a model with non-compensatory decision rules (Swait, 2001; Elrod et al., 2004; Martínez et al., 2009).111111Whereas the decision-maker is assumed to trade all attributes of alternatives and to choose an alternative with the maximum indirect utility in the compensatory RUM framework, non-compensatory rules allow the decision-maker to choose or reject an alternative based on the value of even a single attribute (Schoemaker, 2013).

We argue that a variation in DOF of the t-distributed error kernel across decision-makers provides a flexible way to model “decision uncertainty" – a degree of (un)certainty that decision-makers hold in their choices relative to the variation in the indirect utility of any alternative – without modifying the index function in the latent construct. We illustrate this choice behavior using a plot of the choice probability vs. the utility at different DOF in a binary choice scenario (see Figure 1(b)). Since the relationship between the choice probability and utility is steeper for the respondents with a higher DOF, choice of these decision-makers is more sensitive to a variation in utility, but in a narrower range. More specifically, decision uncertainty of the respondents with the DOF of 10 is much larger as compared to those with the DOF of 0.1.

In discrete choice experiments, decision uncertainty can depend on the familiarity and experience of the decision-maker about the situation encountered during the experiment, among many other factors. Decision-makers belonging to a specific demographic group can be just more certain about their choices than others. Such choice behaviors cannot be modeled using standard logit and probit links, even when accounting for preference heterogeneity. Also, neglecting decision uncertainty results into underestimation of welfare measures (Dekker et al., 2016). A few studies have quantified the decision uncertainty by asking follow-up questions after each choice task and incorporating these self-reported responses as explanatory variables or in other structural forms (Lundhede et al., 2009; Olsen et al., 2011; Beck et al., 2013; Börger, 2016).

What we propose is to capture individual-specific decision uncertainty by parameterizing the DOF of the error kernel as a function of characeteristics of the decision maker. Since the DOF of each individual is obtained as a byproduct of estimation, GCM-t implicitly captures decision-uncertainty behavior without imposing additional econometric structure and also reduces the cognitive burden of respondents by avoiding the need of asking additional questions.

[Insert Figure 1 here]

5 Monte Carlo study and results

Since this is the first study to use a t-distributed error kernel in multinomial response models, we conduct a simulation study to numerically assess the statistical properties (e.g., recovery of true parameters) of the maximum likelihood estimator of the GCM-t model. In another simulation study, we compare the performances of GCM-t and GCM-N models under thick- and thin-tailed error distributions and highlight the consequences of misspecified error kernels.

The context of our simulation study is joint modeling of the commute distance and residential location choice, i.e. integrated land-use transportation modeling. For residential location, we consider a nominal variable with five population-density-based alternatives: 0 – 99, 100-499, 500 – 1499, 1500 – 1999 and 2000 or more households/square mile. The considered joint model is shown in equation 12.

To generate all three exogenous indicators, we take a draw from a standard uniform distribution. If the sampled value is higher than 0.5, the indicator takes a value 1; otherwise, it takes a value 0. That is, for a certain household if the sampled value is 0.64, 0.32, and 0.75, then the generated household is a high-income household with no children, but with a bachelor’s degree holder. All the assumed parameter values and their directions are intuitive and consistent with the literature. For example, a high-income household with children prefers to live in a low-density area (Paleti et al., 2013). Similarly, assuming that the commute distance precedes the choice of residential location (Clark et al., 2003; Rashidi et al., 2012), the likelihood of living in high-density areas decreases with the increase in the commute distance.

Equation 13 shows the considered covariance matrix ($\overline{\bm{\Sigma}}$) of the joint-differenced error, which is normalized up to scale. After taking a draw from this covariance matrix, response variables are generated. Since the commute distance has a structural relationship with the choice of residential location, it is first obtained using the continuous variable model and is then used as an explanatory variable in the choice model to determine the residential location of a household.

[TABLE]

[TABLE]

We initially estimated all elements of the covariance matrix, but such flexible specification resulted in convergence issues – a few elements of the matrix converged to values near zero, leading to a non-positive-definite covariance matrix. This empirical identification concern is commonly encountered, forcing researchers to adopt a diagonal error-covariance matrix. Such concerns are generally not reported, but we are aware of only handful of studies with a non-diagonal error-covariance structure in multinomial choice models with open-form expressions of choice probabilities (e.g., Paleti et al., 2013; Bhat, 2015).

In the simulation study, we estimate diagonal elements of the error-covariance matrix of the choice model and the elements representing the correlation between continuous and multinomial parts of the joint model. The estimated elements of the covariance matrix are in bold in equation 13.121212The fixed elements of the covariance matrix do not retain their true values during optimization due to two reasons. First, we work in the Cholesky space to ensure the positive definitiveness of the variance-covariance matrix. The non-active elements of the variance-covariance matrix are the function of active Cholesky elements which change in optimization iterations. Second, the Cholesky factor is adjusted to normalize the top diagonal element of the variance-covariance matrix of the nominal variable. This is standard practice in probit models. We note that the diagonal error-covariance in a choice model with a t-distributed error kernel does not translate into independence across alternatives as long as the common DOF is finite.

5.1 Statistical properties of GCM-t estimator

To test the statistical properties of the GCM-t model, we consider the above DGP with DOF($\delta$) 1 (DGP-I) and 12 (DGP-II). Whereas a DOF value of 1 represents a distribution far from normal with flat probability curves and long thick tails, a value of 12 mimics the normal error kernel with steep probability curves and thin tails. In both DGPs, we take a sufficiently large sample of 3000 individuals to circumvent the effect of sample size on parameter recovery. Further, we generate 150 resamples for each DGP to ensure that statistical properties are not affected by the choice of the number of repetitions. To evaluate the MVTCD function, we use 200 Halton draws in the earlier discussed separation-of-variables approach. We have also tested sensitivity of parameter estimates and standard errors by increasing the number of draws from 200 to 500 for a few resamples, but have not observed any improvement in results.131313200 Halton draws have proven to be sufficient to estimate up to eight-dimensional integrals in the GHK simulator for the MNP estimation (Patil et al., 2017). We compute the following performance measures across resamples:

Mean estimated value (MEV): Average value of the estimated parameter across all resamples.

Mean absolute bias (MAB): Average bias (|true value-MEV|) of the parameter estimates.

Absolute percentage bias (APB): |MAB/true value|$\times$100.

Finite sample standard error (FSSE): Standard deviation of the parameter estimates across all resamples.

Asymptotic standard error (ASE): Average of standard error values across all resamples, which are obtained using the sandwich estimator. For the sufficient estimator, FSSE and ASE values are close to each other.

Coverage probability (CP): Proportion of times the 95% confidence interval contains the true value.

Power of the test: Proportion of times the null hypothesis is rejected (i.e., |t-statistic| is greater than 1.96).

A lower APB, a ratio of ASE and FSSE closer to 1, higher CP, and higher power are desirable for a better statistical performance of the estimator (Koehler et al., 2009).

The resulting performance measures for DGP-I and DGP-II are summarized in Tables 2 and 3, respectively. We first highlight important insights related to parameter recovery (i.e., bias). The mean APB values of the parameter vector ($\bm{\gamma}$) in the continuous variable model are 1.54% and 0.46% for DGP-I and DGP-II, respectively. These results not only indicate excellent recovery of $\bm{\gamma}$ vector, but also suggest improvement in its recovery with the increase in the DOF. However, parameters $(\bm{\beta})$ associated with the nominal response model appear to be rather difficult to recover accurately irrespective of the DOF. The mean APB values of the $\bm{\beta}$ vector are 14.79% and 16.90% for DGP-I and DGP-II, respectively. These values for the covariance matrix ($\overline{\bm{\Sigma}}$) are 14.16% and 11.37%, respectively, suggesting that the recovery of the error-covariance matrix is better for the DGP with a higher DOF value. The bias in parameter estimates is substantially higher than the ones obtained when the specification is kept the same but the normally-distributed error replaces the t-distributed error kernel in the DGP and in the estimation. More specifically, the mean APB values of $\bm{\beta}$ and $\overline{\bm{\Sigma}}$ in the corresponding GCM-N model are 4.88% and 7.30%, respectively.141414We have also conducted a Monte Carlo study for the corresponding GCM-N model. We used the GHK simulator with 200 Halton to evaluate the MVNCD function. The detailed results are available upon request. This comparison suggests that allowing for flexibility in the parametric distribution of the error kernel beyond mean and variance parameters might lead to a deterioration in the recovery of model parameters. Bhat and Sidharthan (2012) also observed a higher bias in the parameter estimates of the mixed multinomial choice model with a skew-normal distribution than those of the corresponding mixed MNP model. The recovery of the DOF parameter ($\delta$) in both cases is good, with APB values of 8.24% and 12.37%, respectively.

We now discuss performance measures related to model inference – ratio of asymptotic standard error and finite sample standard error (ASE/FSSE), coverage probability (CP), and power. ASE/FSEE values of all parameters $\{\bm{\gamma},\bm{\beta},\overline{\bm{\Sigma}},\delta\}$ are close to 1 for both DGPs, suggesting that the proposed estimator is sufficient. However, consistent with the findings of the parameter recovery, mean values of these ratios are much closer to 1 for DGP-II $\{0.99,0.98,0.95,0.98\}$ as compared to those of the DGP-I $\{1.06,0.81,0.76,1.10\}$. Whereas the mean CP values of the $\bm{\gamma}$ vector are close to 0.95 for both DGPs, these values are much lower for other parameters. Particularly, the mean CP value of $\bm{\beta}$ vector is 0.78 for DGPs with DOF 1 and 12, and these numbers are 0.83 and 0.90 for the error-covariance matrix $\overline{\bm{\Sigma}}$.151515CP of the DOF parameter ($\delta$) in DGP-I is 0.43 (see table 2). Even if the estimated DOF is close to the true DOF (low APB), such low CP occurs due to low standard errors. Similarly, as expected, the power value of all elements of the $\bm{\gamma}$ vector is 1 for both DGPs. The estimator also provides very low variance in power values across the vector $\bm{\beta}$ with the mean values of 0.97 and 0.98 for both DGPs. However, the power variance across error-covariance matrix ($\overline{\bm{\Sigma}}$) slightly increases with the increase in DOF from 1 to 12, but the mean power value decreases from 0.80 to 0.72. Importantly, the power value of all diagonal elements of $\overline{\bm{\Sigma}}$ is 1 across both DGPs, but it is relatively lower for the off-diagonal elements.

[Insert Tables 2 and 3 here]

5.2 Effect of modeling heavy-tailed data with normal distribution

We conduct another Monte Carlo study to understand the consequences of assuming a normally-distributed error kernel when the actual DGP has a heavy-tailed error distribution. To accomplish this, we adopt the same model specification as described in the first Monte Carlo study, but consider a DGP with the DOF of 2 to represent a heavy-tailed distribution. We generate 50 resamples using this DGP, and estimate the GCM-t model with the estimable DOF (correct model specification) and the GCM-t model with the fixed DOF of 300 (which is equivalent to the GCM-N model) for each resample. To benchmark this analysis, we further generate 50 resamples with a thin-tailed error distribution (DOF value of 12) and estimate GCM-t with the estimable DOF and the fixed DOF of 300. The estimation results for DGPs with DOF values of 2 and 12 are provided in Tables 4 and 5, respectively.

We obtain several insights from this second simulation study. First, modeling thick/heavy-tailed data (DOF <= 7) with a normally-distributed error kernel can introduce a large bias in the parameter estimates. The mean APB value (across all parameters) of the GCM-t model with the fixed DOF of 300 (equivalent to GCM-N model) increases from 14.46% to 88.94% with the decrease in DOF of the DGP from 12 to 2. These values are 11.36% and 11.59% for the GCM-t model with the correct error specification. Second, we also observe a high degree of deterioration in the goodness of fit due to error misspecification. Under the heavy-tailed DGP (DOF=2), the loglikelihood value of the GCM-t model with the DOF 2 is 505 points higher than that of the GCM model with the fixed DOF of 300 (Table 4). This difference in the loglikelihood values at convergence is just around 17 points under the DGP with the DOF of 12.

We also replicate this analysis for the true DGP with DOF=1. When we estimate GCM-t with DOF=300 (i.e., equivalent GCM-N model) on resamples of this DGP, some elements of the error-covariance matrix explode and also result into the wrong direction of the parameter estimates, even when only diagonal elements are estimated. This behavior might be a manifestation of the GCM-N model’s attempt to fit the index function and then broaden the support of the normal distribution around the fitted index function to put the probability mass on tails to mimic the underlying DGP. This phenomenon can potentially explain the common convergence concerns faced by researchers in estimating the error-covariance matrix of probit-based choice models. Thus, problems in recovering the error-covariance matrix can indicate the possibility of a misspecified error kernel (i.e. the underlying assumption of normally distributed error kernel might be incorrect, and a t-distributed error kernel may perform better).

We finally compute changes in choice probabilities to reside in neighborhoods with different density levels due to a 50% increase in commute distance for the DGP with the DOF of 2. Table 6 shows the mean and standard deviation of the elasticity estimates across 50 resamples for true GCM-t specification and GCM-t with the fixed DOF of 300 (i.e., GCM-N). We use these estimates to test the null hypothesis that differences in the mean elasticity estimates of both models are zero. The p-values of the two-tailed t-test indicate that the null hypothesis is rejected at the 5% significance level for all alternatives. Thus, discrepancies in parameter estimates of GCM-t and GCM-N also translate into differences in the elasticity estimates. Since GCM-t recovers true parameters much better than GCM-N, its elasticity estimates are more reliable than those of GCM-N.

[Insert Tables 4, 5, and 6 here]

6 Empirical study

We now present statistical and behavioral insights from comparison of the GCM-t and GCM-N models in an empirical setting that jointly models individual’s vehicle-miles-traveled (VMT) and stated vehicle purchase preferences. We also illustrate how decision-uncertainty behavior of different demographic groups can be captured by the GCM-t model.

6.1 Data description

We conducted a stated preference survey of 1542 individuals in 2018 to examine the policies related to on-street parking with charging facilities for battery electric vehicles (EV) in the city of Philadelphia, Pennsylvania. Philadelphia has substantial variation in residents’ socioeconomic attributes, driving patterns, and neighborhood characteristics. The city also has the kinds of short trips and stop-and-go traffic that are well-suited to EVs. Charging locations and neighborhood parking, however, remain a substantial barrier to the adoption of EVs. In 2017, the City Council put a moratorium on a ten-year old policy to permit on-street EV charging stations and parking spaces after installing fewer than one hundred of them throughout the city.

The main goal of the survey was to better understand preferences for adoption of electric vehicles by Philadelphians, within the context of availability of dedicated public parking with charging stations. We thus asked survey respondents to imagine themselves in a situation where they had to buy a new vehicle, had settled on a make and model, and must make a choice about whether to buy an electric version, gasoline version, or no car replacement (opt-out option) given a set of vehicle attributes. In a discrete choice experiment (DCE), each survey participant responded to purchase preferences in eight choice situations based on varying purchase prices, operating costs, electric vehicle performance in terms of driving range, and EV parking characteristics (monthly cost for access, time to recharge the battery as a proxy for type of the charging station, average time to find a parking spot as a proxy for availability, and on/off-street location). Table 7 presents a sample of a choice situation.

[Insert Table 7 here]

The experimental shares of gasoline vehicles, electric vehicles, and the opt-out alternative in the sample are 64%, 32%, and 4%, respectively. The sample summary statistics of key socio-demographic and built-environment attributes are reported in Table 8. Comparing our survey sample to Census micro data of persons over 17 in households with one or more cars, our respondents are substantially more likely to be female (70% vs 52%), white non-Hispanic (64% vs. 36%), younger (39 vs. 44 years old on average), and well-educated (61% vs. 34% with a BA or higher). The income and housing type of respondents are generally representative of Philadelphia’s adult population with cars and those who commute to work by car. For example, 63% live in row homes and 20% live in multi-unit buildings. Although, we had respondents from all Philadelphia Zip Codes, respondents are disproportionately from predominantly white neighborhoods in Northeast Philadelphia. In addition to the vehicle choice experiment, respondents were also asked to report their annual household vehicle-miles-traveled (VMT). In the empirical analysis, we jointly model the household’s annual VMT and choice of vehicle as a function of socio-demographic, built-environment, and alternative-specific attributes.

[Insert Table 8 here]

6.2 Results and discussion

Table 9 summarizes the results of GCM-t and GCM-N models, which we discuss in the next subsections.

[Insert Table 9 here]

6.2.1 Socio-demographic and built-environment factors

In this section, we discuss the relation of individual’s or household’s characteristics with household’s annual VMT and vehicle purchase preferences. Several demographic characteristics have a significant association with the household’s VMT in both GCM-t and GCM-N models. The results indicate that married households tend to have lower VMT as compared to unmarried households, ceteris paribus. This association can be a manifestation of the way different types of households spend time on their activities. For example, whereas married (or larger) households are likely to spend more time together at home in joint activities (Fang, 2008; Spissu et al., 2009), unmarried individuals may drive more frequently for various leisure activities (e.g., going to a club, theater, or eating out) besides commuting trips.

As expected, a household with a full-time working individual is likely to drive more than those with part-time or non-workers (Lee and McDonald, 2003; McQuaid and Chen, 2012). Surprisingly, households in higher density Zip Codes tend to drive more than those in lower density ones. This stands in sharp contrast with findings from recent meta-analyses (Ewing and Cervero, 2010; Stevens, 2017), as well as findings from the Philadelphia region more specifically (Klein et al., 2018). The counter-intuitive finding may relate to our sample excluding non-drivers and drawing disproportionately from higher-educated, white, women drivers in relatively dense neighborhoods of Northeast Philadelphia that are somewhat far from Philadelphia’s major employment centers.

We now discuss the factors determining the vehicle purchase preferences of Philadelphia residents. We discuss the effect of socio-demographic attributes, followed by built-environment variables. First, married households with children are less inclined toward the purchase of EVs as compared to single individuals, possibly due to the compact size of most current electric cars. Second, households with highly educated individuals (Masters or post-graduate degree) are more likely to purchase electric vehicles as compared to less educated households, perhaps as a sign of highly responsible behavior toward the planet. The literature also confirms that highly educated individuals exhibit a higher degree of environmental consciousness in terms of recycling, purchase of organic food, and awareness about factors contributing to global warming (Fisher et al., 2012). Third, the number of driving license holders has a positive association with the likelihood of purchasing an EV. This is because households with more drivers can better perceive the long-term benefits of buying fuel-efficient EVs as they are likely to have higher VMT. Fourth, as expected, households with a higher number of vehicles are less likely to purchase EVs. This is perhaps because car-lovers generally have higher income and may view compact design and lower range of EVs as major barriers to adoption (Jakobsson et al., 2016). Fifth, households who already own a hybrid EV also exhibit a higher inclination toward purchasing an EV. Sixth, male respondents are more likely to purchase EVs than females. Whereas most studies have found that females are more environmental-friendly than males, some studies have also reported males to be greener than females (Fisher et al., 2012). Seventh, the results indicate that Asians are more likely to purchase EVs followed by Caucasians. However, low-income groups with generally less access to opportunities (such as African-American and Hispanics) are more likely to stick to traditional gasoline vehicles. Finally, as expected, millennials are more likely to purchase EVs than baby boomers and generation-X individuals. Millennials grew up in an era of prevalent information on topics such as global warming and environmental sustainability, which perhaps made them more inclined to adopt technology aimed at benefiting the environment.

Both built-environment (walk-score and population density of the neighborhood) factors have expected and intuitive directions of effects on vehicle purchase preferences. In terms of structural effects, households with a higher vehicle mileage are less likely to purchase EVs, possibly due to limitations in driving range or lower concern for the environment.

6.2.2 Willingness to pay estimates

We now analyze the results related to alternative-specific explanatory variables in the choice model. The directions of effects are the same and intuitive in both GCM-t and GCM-N models. Since the magnitude of marginal-utility parameters is not directly comparable, we derive willingness to pay (WTP) for improving various EV attributes.

Figure 2 presents the maximum WTP as premium for the purchase of an electric car with a marginal improvement in driving range.161616Since driving range enters the utility specification logarithmically, the WTP estimates are non-linear and decreases with the driving range. The results of GCM-t and GCM-N models indicate that Philadelphians are willing to pay additional $212 (GCM-t) and$153 (GCM-N) in purchase price, respectively, to increase driving range by one mile for a car offering 150 miles with a full battery. These WTP estimates decrease to $64 (GCM-t) and$46 (GCM-N) when the driving range of the car reaches 500 miles, which is close to parity with gasoline cars. Clearly, GCM-t estimates are higher than those of GCM-N, but the difference decreases with the increase in driving range.

[Insert Figure 2 here]

Since monthly cost of having access to public parking with EV charging was an experimental attributes, we also estimate WTP measures for characteristics of the parking spot. For instance, we derive the WTP to reduce average search time for available EV parking. According to GCM-t and GCM-N estimates, Philadelphians are willing to add $6.7 and$4.6 in their monthly parking cost to reduce parking search time by a minute, respectively. This conforms to findings that on-street residential parking is frequently underpriced (Shoup, 2005). In a neighborhood where drivers spend just five minutes searching for parking on average, the average stated willingness to pay to avoid that search time is eight to eleven times higher than the price of a parking permit for a single vehicle.

WTP to reduce the EV charging time by an hour ranges from $96 (GCM-t) to$60 (GCM-N). The higher WTP estimates of GCM-t are aligned with the finding of the study by Dekker et al. (2016), who also observed an increase in WTP for flood risk reductions after controlling for the behavioral responses to decision uncertainty in an integrated choice and latent variable model.

Finally, because vehicle purchase and parking costs happen at different times, the annual subjective discount rate of parking cost is estimated at 10.7% (GCM-t) and 9.3% (GCM-N), which is slightly above market interest rates in the automotive industry.

6.2.3 Behavioral Insights

We parameterize the DOF of the t-distributed error kernel as a function of demographics in GCM-t to capture decision-uncertainty behavior of decision-makers. Table 10 presents the relation between the DOF and demographics in GCM-t, which is used to obtain the DOF of each respondent. These results also provide insights about decision uncertainty of different socio-demographic groups. For example, the positive relationship between married males and the DOF indicates that married males are likely to have a higher DOF, and thus are less certain about their choices as compared to unmarried females. We are not aware of any study which has measured decision-uncertainty in the context of the vehicle purchase. Therefore, beyond what observed in the current empirical analysis, we do not provide comparisons with the existing literature.

[Insert Table 10 here]

We also compute elasticity estimates with respect to 1% and 25% reduction in the EV parking price for both models. Tables 11 and 12 present elasticity estimates for 1% and 25% reduction, respectively, for ten different ranges of the DOF such that each bin contains 10% respondents of the sample.171717Note that there is no concept of DOF in the GCM-N model. To facilitate the comparison between GCM-t and GCM-N models, DOF ranges are obtained from the GCM-t model and then the same set of individuals are used for respective calculations in the GCM-N models. Whereas the difference between elasticity estimates of GCM-t and GCM-N models are not apparent for the 1% reduction, they are quite stark for the 25% reduction in parking price (see Table 12). In general, the magnitude of elasticity estimates for gasoline and electric alternatives in the GCM-N model is higher (1.5 times, on average) than that of the corresponding GCM-t model. Specifically, this ratio is higher for respondents with lower DOF values. This supports the our earlier observation that individuals with lower DOF are more certain about their choices and are less sensitive to changes in the utility of alternatives.

[Insert Tables 11 and 12 here]

Finally, we plot the change in the probability of choosing electric and gasoline vehicles in GCM-t and GCM-N model as a function of the change in the utility of electric vehicle in Figures 3(a) and 3(b), respectively. For illustration, we only make plots for three DOF ranges. These plots offer many interesting insights. First, plots for three DOF ranges in the GCM-N model are virtually the same. Second, these change-in-probability plots become steeper with the increase in the DOF in the GCM-t model, capturing the varied decision-uncertainty behavior in the sample. Finally, curves for the lower DOFs are not asymptotic to the x-axis, demonstrating a phenomenon that it is implausible to achieve 100% market share for any alternative in spite of being far superior to other available options.

[Insert Figure 3 here]

6.2.4 Prediction performance

The dataset used in this case study is a good example of class-imbalance as the opt-out alternative has a very small sample share. This is an appropriate scenario to compare the predicted probability of the chosen alternative in GCM-t and GCM-N model. Table 13 presents the ratio of GCM-t and GCM-N choice probabilities for a group of respondents with a varying range of DOF.

The probability ratios for two alternatives (Gasoline and Electric vehicle which have 96% sample share in total) are close to 1 across all DOF ranges. However, the average ratio of 1.15 for the opt-out option (4% sample share) confirms our earlier observation that GCM-t model is indeed better than a GCM-N model in modeling such imbalanced datasets. These ratios for the opt-out option are significantly higher than 1 in lower ranges of the DOF ($\leq$ 3.79). This finding is aligned with our discussion in subsection 4.1 on the ability of GCM-t to better predict the preference for the low-share alternatives under heavy-tailed error distributions.

Surprisingly, GCM-N outperforms GCM-t for DOF ranges 4.16-4.40 and 5.33-9.32 with probability ratios 0.77 and 0.80, respectively. Since GCM-t can approximate GCM-N with the higher DOF, we would not expect this pattern. We speculate that such anomalies might be a manifestation of the strict constraint to use the same DOF expression across all alternatives of the choice model and the continuous response model. Nonetheless, the superior performance of GCM-N over GCM-t for certain DOF values can be leveraged by developing a latent class model with two classes where classes differ based on the distributional assumptions on the error kernel (i.e., t and normal). We discuss another flexible way to relax the “same DOF” assumption in section 7.

We also conduct a five-fold cross-validation analysis to compare the predictive performance of GCM-t and GCM-N. We use a variant of the Brier score as a performance measure, which is a function of the predicted choice probability distribution. A lower Brier score implies better predictive performance. Table 14 shows the mean and the standard deviation of Brier scores for in-sample (70% of total sample) and out-of-sample (30% of total sample) predictions across five samples. In terms of the mean estimates of the Brier score, GCM-t performs better than GCM-N in in-sample predictions, but GCM-N outperforms GCM-t in out-of-sample predictions of Gasoline and opt-out alternatives. However, we cannot conclude much from this empirical study because these differences are not statistically significant due to high standard deviations. Therefore, a detailed Monte Carlo study and more empirical applications are necessary to evaluate the relative predictive performances of GCM-t and GCM-N in practice.

[Insert tables 13 and 14 here]

6.2.5 Model selection

We first compare the goodness-of-fit measures of GCM-t and GCM-N models. The log-likelihood value at convergence of GCM-t is better than the corresponding GCM-N model by 200 points (see Table 9). Further, the lower Bayesian information criterion (BIC) value of GCM-t also suggests that the DGP of this empirical study is better modeled using the GCM-t model (with 39 parameters) than the GCM-N model (with 42 parameters).

Apart from these measures, we also compare the trace of the error-covariance matrix of GCM-t and GCM-N to assess the amount of (un)explained variance in these models (see Table 11). The trace values of the error-covariance matrix in GCM-t and GCM-N models are 1.06 and 3.44, respectively. This result is as expected because parameterization of the DOF in GCM-t model enables it to explain more variation in the DGP than the GCM-N model.

[Insert Table 15 here]

Piatek and Gensowski (2017) did a similar effort to reduce the unexplained variance in the MNP model with latent factors. The authors split the error-covariance matrix into two parts – the allocation (or factor loading) matrix which can be specified to model a covariance structure across alternatives, and the idiosyncratic error matrix which captures the correlation in residual. The key idea behind the allocation matrix is to allow the researcher to specify driving factors behind the choice variation beyond the observables included in the deterministic part of the utility. However, the allocation matrix is not straightforward to define in empirical studies and, beyond a mathematical inconvenience, is difficult to parametrize as a function of other variables. Thus, in contrast to splitting the error-covariance matrix, parametrization of the DOF in the GCM-t model not only offers insights about choice behavior but also provides a new parsimonious specification to reduce the unexplained variance.

7 Conclusions and future work

In this study, we have proposed and applied for the first time a t-distributed error kernel within random-utility-maximization (RUM) multinomial choice models. We have not only discussed statistical advantages of using this t-distributed kernel over a normally-distributed one in class-imbalance situations, but also illustrated how the t-distributed error kernel can capture decision-uncertainty behavior (i.e., how certain the decision-makers are about the choices that they make). Furthermore, we have extended this model to a generalized continuous-multinomial response model with a t-distributed error kernel (GCM-t). Using the composite marginal likelihood method, separation-of-variables approach, and properties of the t-distribution, we have derived the full-information maximum likelihood estimator of the GCM-t model and tested its statistical properties in a Monte Carlo study. Finally, we have compared the statistical performance of GCM-t and GCM with normally-distributed error kernel (GCM-N) in a simulation study, and validated behavioral hypotheses about the advantages of GCM-t over GCM-N in an empirical study related to the adoption of electric vehicles in the city of Philadelphia.

The results of the simulation study indicate that GCM-t and GCM-N perform equally well in terms of recovering model parameters, goodness-of-fit, and model inference when the true error distribution is thin-tailed. However, GCM-t either outperforms GCM-N by a significant margin or GCM-N may not even converge in datasets with heavy-tailed error distributions. In the case study about adoption of electric vehicles, GCM-t is more accurate than GCM-N in estimating the choice of the alternative with a small sample share. Moreover, accounting for decision-uncertainty behavior in GCM-t manifests into lower elasticity estimates (relative to those of GCM-N). These lower elasticity estimates of GCM-t, as well as the evidence we found for underestimation in the GCM-N model of individual’s willingness to pay to increase the driving range of electric vehicles and for reducing the parking search time for a parking spot with charging are all relevant in the context of planning policies for broader adoption of electric vehicles.

We now discuss key limitations of this study and possible research avenues to overcome them in the future. First, unlike GCM-N, incorporating random parameters with parametric heterogeneity distributions in GCM-t is possible, but not straightforward because the resulting distribution – the sum of t-distributed random variables or sum of t-distributed random variable and other parametric distributions – is not of a known form. Therefore, even inclusion of parametric unobserved preference heterogeneity in the GCM-t model requires an additional layer of simulation to compute the likelihood function of the model. This necessity of adding a simulation layer would motivate researchers to rather incorporate flexible semi-parametric (instead of parametric) preference heterogeneity distributions in the GCM-t model (see Train, 2016; Vij and Krueger, 2017; Bansal, Daziano and Achtnicht, 2018, for recent developments in semi-parametrics). We leave this extension for future work. Second, we considered the same degrees of freedom across all alternatives of the choice model and the continuous variable model because the exact distribution of linear or non-linear combinations of two t-distributed random variables with a arbitrary degrees-of-freedom values is not known. Extending GCM-t to a more flexible model with varying degrees of freedom across alternatives using a copula-based approach (Bhat and Eluru, 2009) is a potential direction for future developments. This copula-based extension would allow researchers to capture decision-uncertainty behavior across alternatives, not just across decision-makers. Third, the proposed GCM-t model cannot handle asymmetric error distributions. From an estimation standpoint, extending GCM-t model to a GCM model with a skew-t-distributed error kernel is not challenging because all relevant statistical properties are well-established in the literature (Azzalini and Genton, 2008). However, practical concerns related to model inference and asymptotic properties of this new estimator require further investigation (Azzalini and Arellano-Valle, 2013). Fourth, the proposed model can be extended to accommodate both situations – structural dependence and common unobserved errors in multinomial and continuous outcomes. For instance, the specification can be modified to simultaneously incorporate multiple causality directions by embedding a latent class framework where each class represents a causality direction. Fifth, the current specification implicitly assigns the appropriate DOF to decision-makers, but it only captures deterministic heterogeneity in DOF. The DOF can also have variation within each demographic group. The model can be extended to capture this unobserved heterogeneity in the DOF by introducing a stochastic term in its parameterization. This generalization is straightforward but can lead to a loss in computational efficiency due to an additional layer of simulation in the estimation.

Finally, the full-information maximum likelihood estimation of GCM-t is slow because it relies on a numerical gradient in the absence of tractable expressions of the analytical gradient and Hessian matrix. The computational efficiency of this estimator can be improved by adopting tools like PyTorch or TensorFlow, which use automatic differentiation to provide analytic gradients and Hessians without customized derivations. The possibilities of developing a Bayesian estimator (Kim et al., 2007) or an expectation-maximization procedure (Liu, 2004; Bansal, Daziano and Guerra, 2018) can also be explored in the future to speed up the estimation of the GCM-t model.

Acknowledgements

PB and RAD are thankful to the National Science Foundation CAREER Award CBET-1253475 for financially supporting methodological contributions of this research. PB is also thankful to Prof. Joan Walker and Prof. Kenneth Train for their guidance during his visit to UC Berkeley under the Exchange Scholar Program. RAD and EG also acknowledge support from a faculty research grant from the Kleinman Center for Energy Policy of the University of Pennsylvania, which supported collection of the data used as empirical case study. The authors are thankful to three anonymous reviewers for their detailed comments which significantly improved the manuscript.

Appendix A Matrix transformations

A.1 Transformation matrix $(\bm{D})$ to compute $\bm{\Lambda}$ from $\overline{\bm{\Lambda}}$

We show below the relation between the variance-covariance matrix of the error $\left(\bm{\Lambda}_{I_{K}\times I_{K}}\right)$ and the normalized variance-covariance matrix of the error difference $\left(\bm{\overline{\Lambda}}_{(I_{K}-I)\times(I_{K}-I)}\right)$:

[TABLE]

Note that $\bm{\Lambda}=\bm{D}\bm{\overline{\Lambda}}\bm{D}^{\bm{\top}}$, where the transformation matrix $\bm{D}$ is constructed based on algorithm 2.

We provide an example to illustrate the transformation from $\bm{\overline{\Lambda}}$ to $\bm{\Lambda}$ using $\bm{D}$. We consider a case of $I=2$, with $1_{K}=3$ and $2_{K}=4$. The $\bm{\overline{\Lambda}}_{5\times 5}$, $\bm{\Lambda}_{7\times 7}$, and the transformation matrix $\bm{D}_{7\times 5}$ would be:

[TABLE]

A.2 Modified transformation matrix $(\bm{D}_{m})$ to compute $\bm{\Sigma}$ from $\overline{\bm{\Sigma}}$

The modified transformation matrix $(\bm{D}_{m})$ can be easily computed from the transformation matrix $(\bm{D})$ by appending an identity matrix as follows:

[TABLE]

where $H$ is the number of continuous outcomes, $\bm{1}_{H\times H}$ is an identity matrix of size $H\times H$, and $\bm{D}$ is obtained from algorithm 2. We expand on the example of appendix A.1 which considers $I=2$ with $1_{K}=3$ and $2_{K}=4$. If $H=2$, then the modified transformation matrix $\bm{D}_{m}$ would be:

[TABLE]

A.3 Utility difference generator $(M)$ to compute $\bm{\widetilde{\Sigma}}$ from $\bm{\Sigma}$

We consider the same example of appendix A.2 which considers $H=2$ and $I=2$ with $1_{K}=3$, $1_{m}=2$, $2_{K}=4$, and $2_{m}=3$. We use algorithm 3 to construct $\bm{M}$ matrix for this example:

[TABLE]

A.4 Reparametrization of the Cholesky decomposition of $\bm{\overline{\Sigma}}$

Consider $\bm{L}\bm{L}^{\bm{\top}}=\bm{\overline{\Sigma}}$, where $\bm{L}$ is the lower triangular Cholesky matrix of size $(H+I_{K}-I)\times(H+I_{K}-I)$. We reparametrize all rows of the $\bm{L}$ matrix for which the diagonal element of $\bm{\overline{\Sigma}}$ is normalized to 1. We compute $a_{i}=\sqrt{1+\sum\limits_{j=1}^{i-1}L_{i,j}^{2}}$ for the $i^{th}$ row and modify non-diagonal elements $L_{i,r}=\frac{L_{i,r}}{a_{i}}\forall r\in\{1,2,\dots,i-1\}$ and the diagonal element $L_{i,i}=\frac{1}{a_{i}}$.

Appendix B MVTNCD illustration

The derivation of equation 8 from equation 6 can be understood easily based on a transformation applied on a bivariate t-cumulative distribution function. For $p=2$ , equation 6 after substitution can be written as (note that $w_{1}=u_{1}$ ):

[TABLE]

Figures

Tables