Note

Go to the end to download the full example code.

Serial correlation

Choice model with the latent variable. Mixture of logit, with agent effect to deal with serial correlation. Measurement equation for the indicators. Maximum likelihood (full information) estimation.

author:: Michel Bierlaire, EPFL
date:: Thu Apr 13 18:16:37 2023

import sys

import biogeme.biogeme as bio
import biogeme.biogeme_logging as blog
import biogeme.results as res
from biogeme import models
from biogeme.exceptions import BiogemeError
from biogeme.expressions import (
    Beta,
    bioDraws,
    MonteCarlo,
    Elem,
    bioNormalCdf,
    exp,
    log,
)
from biogeme.data.optima import (
    read_data,
    age_65_more,
    moreThanOneCar,
    moreThanOneBike,
    individualHouse,
    male,
    haveChildren,
    haveGA,
    highEducation,
    WaitingTimePT,
    Envir01,
    Envir02,
    Envir03,
    Mobil11,
    Mobil14,
    Mobil16,
    Mobil17,
    Choice,
    TimePT_scaled,
    TimeCar_scaled,
    MarginalCostPT_scaled,
    CostCarCHF_scaled,
    PurpHWH,
    PurpOther,
    distance_km_scaled,
    ScaledIncome,
)
from read_or_estimate import read_or_estimate

logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example b06serial_correlation.py')

Example b06serial_correlation.py

Read the estimates from the structural equation estimation.

MODELNAME = 'b05latent_choice_full'
try:
    struct_results = res.bioResults(pickle_file=f'saved_results/{MODELNAME}.pickle')
except BiogemeError:
    print(
        f'Run first the script {MODELNAME}.py in order to generate the '
        f'file {MODELNAME}.pickle, and move it to the directory saved_results'
    )
    sys.exit()
betas = struct_results.get_beta_values()

Coefficients.

coef_intercept = Beta('coef_intercept', betas['coef_intercept'], None, None, 0)
coef_age_65_more = Beta('coef_age_65_more', betas['coef_age_65_more'], None, None, 0)
coef_haveGA = Beta('coef_haveGA', betas['coef_haveGA'], None, None, 0)

coef_moreThanOneCar = Beta(
    'coef_moreThanOneCar', betas['coef_moreThanOneCar'], None, None, 0
)
coef_moreThanOneBike = Beta(
    'coef_moreThanOneBike', betas['coef_moreThanOneBike'], None, None, 0
)
coef_individualHouse = Beta(
    'coef_individualHouse', betas['coef_individualHouse'], None, None, 0
)
coef_male = Beta('coef_male', betas['coef_male'], None, None, 0)
coef_haveChildren = Beta('coef_haveChildren', betas['coef_haveChildren'], None, None, 0)
coef_highEducation = Beta(
    'coef_highEducation', betas['coef_highEducation'], None, None, 0
)

Latent variable: structural equation.

Define a random parameter, normally distributed, designed to be used for Monte-Carlo integration.

omega = bioDraws('omega', 'NORMAL')
sigma_s = Beta('sigma_s', betas['sigma_s'], None, None, 0)

Deal with serial correlation by including an error component that is individual specific

error_component = bioDraws('error_component', 'NORMAL')
ec_sigma = Beta('ec_sigma', 10, None, None, 0)

thresholds = [None, 4, 6, 8, 10, None]
betas_thresholds = [
    Beta(
        'beta_ScaledIncome_minus_inf_4',
        betas['beta_ScaledIncome_minus_inf_4'],
        None,
        None,
        0,
    ),
    Beta(
        'beta_ScaledIncome_4_6',
        betas['beta_ScaledIncome_4_6'],
        None,
        None,
        0,
    ),
    Beta(
        'beta_ScaledIncome_6_8',
        betas['beta_ScaledIncome_6_8'],
        None,
        None,
        0,
    ),
    Beta(
        'beta_ScaledIncome_8_10',
        betas['beta_ScaledIncome_8_10'],
        None,
        None,
        0,
    ),
    Beta(
        'beta_ScaledIncome_10_inf',
        betas['beta_ScaledIncome_10_inf'],
        None,
        None,
        0,
    ),
]

formula_income = models.piecewise_formula(
    variable=ScaledIncome,
    thresholds=thresholds,
    betas=betas_thresholds,
)

CARLOVERS = (
    coef_intercept
    + coef_age_65_more * age_65_more
    + formula_income
    + coef_moreThanOneCar * moreThanOneCar
    + coef_moreThanOneBike * moreThanOneBike
    + coef_individualHouse * individualHouse
    + coef_male * male
    + coef_haveChildren * haveChildren
    + coef_haveGA * haveGA
    + coef_highEducation * highEducation
    + sigma_s * omega
    + ec_sigma * error_component
)

Measurement equations.

Intercepts.

INTER_Envir01 = Beta('INTER_Envir01', 0, None, None, 1)
INTER_Envir02 = Beta('INTER_Envir02', betas['INTER_Envir02'], None, None, 0)
INTER_Envir03 = Beta('INTER_Envir03', betas['INTER_Envir03'], None, None, 0)
INTER_Mobil11 = Beta('INTER_Mobil11', betas['INTER_Mobil11'], None, None, 0)
INTER_Mobil14 = Beta('INTER_Mobil14', betas['INTER_Mobil14'], None, None, 0)
INTER_Mobil16 = Beta('INTER_Mobil16', betas['INTER_Mobil16'], None, None, 0)
INTER_Mobil17 = Beta('INTER_Mobil17', betas['INTER_Mobil17'], None, None, 0)

Coefficients.

B_Envir01_F1 = Beta('B_Envir01_F1', -1, None, None, 1)
B_Envir02_F1 = Beta('B_Envir02_F1', betas['B_Envir02_F1'], None, None, 0)
B_Envir03_F1 = Beta('B_Envir03_F1', betas['B_Envir03_F1'], None, None, 0)
B_Mobil11_F1 = Beta('B_Mobil11_F1', betas['B_Mobil11_F1'], None, None, 0)
B_Mobil14_F1 = Beta('B_Mobil14_F1', betas['B_Mobil14_F1'], None, None, 0)
B_Mobil16_F1 = Beta('B_Mobil16_F1', betas['B_Mobil16_F1'], None, None, 0)
B_Mobil17_F1 = Beta('B_Mobil17_F1', betas['B_Mobil17_F1'], None, None, 0)

Linear models.

MODEL_Envir01 = INTER_Envir01 + B_Envir01_F1 * CARLOVERS
MODEL_Envir02 = INTER_Envir02 + B_Envir02_F1 * CARLOVERS
MODEL_Envir03 = INTER_Envir03 + B_Envir03_F1 * CARLOVERS
MODEL_Mobil11 = INTER_Mobil11 + B_Mobil11_F1 * CARLOVERS
MODEL_Mobil14 = INTER_Mobil14 + B_Mobil14_F1 * CARLOVERS
MODEL_Mobil16 = INTER_Mobil16 + B_Mobil16_F1 * CARLOVERS
MODEL_Mobil17 = INTER_Mobil17 + B_Mobil17_F1 * CARLOVERS

Scale parameters.

SIGMA_STAR_Envir01 = Beta('SIGMA_STAR_Envir01', 1, None, None, 1)
SIGMA_STAR_Envir02 = Beta(
    'SIGMA_STAR_Envir02', betas['SIGMA_STAR_Envir02'], None, None, 0
)
SIGMA_STAR_Envir03 = Beta(
    'SIGMA_STAR_Envir03', betas['SIGMA_STAR_Envir03'], None, None, 0
)
SIGMA_STAR_Mobil11 = Beta(
    'SIGMA_STAR_Mobil11', betas['SIGMA_STAR_Mobil11'], None, None, 0
)
SIGMA_STAR_Mobil14 = Beta(
    'SIGMA_STAR_Mobil14', betas['SIGMA_STAR_Mobil14'], None, None, 0
)
SIGMA_STAR_Mobil16 = Beta(
    'SIGMA_STAR_Mobil16', betas['SIGMA_STAR_Mobil16'], None, None, 0
)
SIGMA_STAR_Mobil17 = Beta(
    'SIGMA_STAR_Mobil17', betas['SIGMA_STAR_Mobil17'], None, None, 0
)

Symmetric thresholds.

delta_1 = Beta('delta_1', betas['delta_1'], 0, 10, 0)
delta_2 = Beta('delta_2', betas['delta_2'], 0, 10, 0)
tau_1 = -delta_1 - delta_2
tau_2 = -delta_1
tau_3 = delta_1
tau_4 = delta_1 + delta_2

Ordered probit models.

Envir01_tau_1 = (tau_1 - MODEL_Envir01) / SIGMA_STAR_Envir01
Envir01_tau_2 = (tau_2 - MODEL_Envir01) / SIGMA_STAR_Envir01
Envir01_tau_3 = (tau_3 - MODEL_Envir01) / SIGMA_STAR_Envir01
Envir01_tau_4 = (tau_4 - MODEL_Envir01) / SIGMA_STAR_Envir01
IndEnvir01 = {
    1: bioNormalCdf(Envir01_tau_1),
    2: bioNormalCdf(Envir01_tau_2) - bioNormalCdf(Envir01_tau_1),
    3: bioNormalCdf(Envir01_tau_3) - bioNormalCdf(Envir01_tau_2),
    4: bioNormalCdf(Envir01_tau_4) - bioNormalCdf(Envir01_tau_3),
    5: 1 - bioNormalCdf(Envir01_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Envir01 = Elem(IndEnvir01, Envir01)


Envir02_tau_1 = (tau_1 - MODEL_Envir02) / SIGMA_STAR_Envir02
Envir02_tau_2 = (tau_2 - MODEL_Envir02) / SIGMA_STAR_Envir02
Envir02_tau_3 = (tau_3 - MODEL_Envir02) / SIGMA_STAR_Envir02
Envir02_tau_4 = (tau_4 - MODEL_Envir02) / SIGMA_STAR_Envir02
IndEnvir02 = {
    1: bioNormalCdf(Envir02_tau_1),
    2: bioNormalCdf(Envir02_tau_2) - bioNormalCdf(Envir02_tau_1),
    3: bioNormalCdf(Envir02_tau_3) - bioNormalCdf(Envir02_tau_2),
    4: bioNormalCdf(Envir02_tau_4) - bioNormalCdf(Envir02_tau_3),
    5: 1 - bioNormalCdf(Envir02_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Envir02 = Elem(IndEnvir02, Envir02)

Envir03_tau_1 = (tau_1 - MODEL_Envir03) / SIGMA_STAR_Envir03
Envir03_tau_2 = (tau_2 - MODEL_Envir03) / SIGMA_STAR_Envir03
Envir03_tau_3 = (tau_3 - MODEL_Envir03) / SIGMA_STAR_Envir03
Envir03_tau_4 = (tau_4 - MODEL_Envir03) / SIGMA_STAR_Envir03
IndEnvir03 = {
    1: bioNormalCdf(Envir03_tau_1),
    2: bioNormalCdf(Envir03_tau_2) - bioNormalCdf(Envir03_tau_1),
    3: bioNormalCdf(Envir03_tau_3) - bioNormalCdf(Envir03_tau_2),
    4: bioNormalCdf(Envir03_tau_4) - bioNormalCdf(Envir03_tau_3),
    5: 1 - bioNormalCdf(Envir03_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Envir03 = Elem(IndEnvir03, Envir03)

Mobil11_tau_1 = (tau_1 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
Mobil11_tau_2 = (tau_2 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
Mobil11_tau_3 = (tau_3 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
Mobil11_tau_4 = (tau_4 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
IndMobil11 = {
    1: bioNormalCdf(Mobil11_tau_1),
    2: bioNormalCdf(Mobil11_tau_2) - bioNormalCdf(Mobil11_tau_1),
    3: bioNormalCdf(Mobil11_tau_3) - bioNormalCdf(Mobil11_tau_2),
    4: bioNormalCdf(Mobil11_tau_4) - bioNormalCdf(Mobil11_tau_3),
    5: 1 - bioNormalCdf(Mobil11_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Mobil11 = Elem(IndMobil11, Mobil11)

Mobil14_tau_1 = (tau_1 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
Mobil14_tau_2 = (tau_2 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
Mobil14_tau_3 = (tau_3 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
Mobil14_tau_4 = (tau_4 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
IndMobil14 = {
    1: bioNormalCdf(Mobil14_tau_1),
    2: bioNormalCdf(Mobil14_tau_2) - bioNormalCdf(Mobil14_tau_1),
    3: bioNormalCdf(Mobil14_tau_3) - bioNormalCdf(Mobil14_tau_2),
    4: bioNormalCdf(Mobil14_tau_4) - bioNormalCdf(Mobil14_tau_3),
    5: 1 - bioNormalCdf(Mobil14_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Mobil14 = Elem(IndMobil14, Mobil14)

Mobil16_tau_1 = (tau_1 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
Mobil16_tau_2 = (tau_2 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
Mobil16_tau_3 = (tau_3 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
Mobil16_tau_4 = (tau_4 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
IndMobil16 = {
    1: bioNormalCdf(Mobil16_tau_1),
    2: bioNormalCdf(Mobil16_tau_2) - bioNormalCdf(Mobil16_tau_1),
    3: bioNormalCdf(Mobil16_tau_3) - bioNormalCdf(Mobil16_tau_2),
    4: bioNormalCdf(Mobil16_tau_4) - bioNormalCdf(Mobil16_tau_3),
    5: 1 - bioNormalCdf(Mobil16_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Mobil16 = Elem(IndMobil16, Mobil16)

Mobil17_tau_1 = (tau_1 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
Mobil17_tau_2 = (tau_2 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
Mobil17_tau_3 = (tau_3 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
Mobil17_tau_4 = (tau_4 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
IndMobil17 = {
    1: bioNormalCdf(Mobil17_tau_1),
    2: bioNormalCdf(Mobil17_tau_2) - bioNormalCdf(Mobil17_tau_1),
    3: bioNormalCdf(Mobil17_tau_3) - bioNormalCdf(Mobil17_tau_2),
    4: bioNormalCdf(Mobil17_tau_4) - bioNormalCdf(Mobil17_tau_3),
    5: 1 - bioNormalCdf(Mobil17_tau_4),
    6: 1.0,
    -1: 1.0,
    -2: 1.0,
}

P_Mobil17 = Elem(IndMobil17, Mobil17)

Choice model.

ASC_CAR = Beta('ASC_CAR', betas['ASC_CAR'], None, None, 0)
ASC_PT = Beta('ASC_PT', 0, None, None, 1)
ASC_SM = Beta('ASC_SM', betas['ASC_SM'], None, None, 0)
BETA_COST_HWH = Beta('BETA_COST_HWH', betas['BETA_COST_HWH'], None, None, 0)
BETA_COST_OTHER = Beta('BETA_COST_OTHER', betas['BETA_COST_OTHER'], None, None, 0)
BETA_DIST = Beta('BETA_DIST', betas['BETA_DIST'], None, None, 0)
BETA_TIME_CAR_REF = Beta('BETA_TIME_CAR_REF', betas['BETA_TIME_CAR_REF'], None, 0, 0)
BETA_TIME_CAR_CL = Beta('BETA_TIME_CAR_CL', betas['BETA_TIME_CAR_CL'], -10, 10, 0)
BETA_TIME_PT_REF = Beta('BETA_TIME_PT_REF', betas['BETA_TIME_PT_REF'], None, 0, 0)
BETA_TIME_PT_CL = Beta('BETA_TIME_PT_CL', betas['BETA_TIME_PT_CL'], -10, 10, 0)
BETA_WAITING_TIME = Beta('BETA_WAITING_TIME', betas['BETA_WAITING_TIME'], None, None, 0)

Definition of utility functions.

BETA_TIME_PT = BETA_TIME_PT_REF * exp(BETA_TIME_PT_CL * CARLOVERS)

V0 = (
    ASC_PT
    + BETA_TIME_PT * TimePT_scaled
    + BETA_WAITING_TIME * WaitingTimePT
    + BETA_COST_HWH * MarginalCostPT_scaled * PurpHWH
    + BETA_COST_OTHER * MarginalCostPT_scaled * PurpOther
    + ec_sigma * error_component
)

BETA_TIME_CAR = BETA_TIME_CAR_REF * exp(BETA_TIME_CAR_CL * CARLOVERS)

V1 = (
    ASC_CAR
    + BETA_TIME_CAR * TimeCar_scaled
    + BETA_COST_HWH * CostCarCHF_scaled * PurpHWH
    + BETA_COST_OTHER * CostCarCHF_scaled * PurpOther
    + ec_sigma * error_component
)

V2 = ASC_SM + BETA_DIST * distance_km_scaled

Associate utility functions with the numbering of alternatives.

V = {0: V0, 1: V1, 2: V2}

Conditional on the random parameters, we have a logit model (called the kernel) for the choice.

condprob = models.logit(V, None, Choice)

Conditional on the random parameters, we have the product of ordered probit for the indicators.

condlike = (
    P_Envir01
    * P_Envir02
    * P_Envir03
    * P_Mobil11
    * P_Mobil14
    * P_Mobil16
    * P_Mobil17
    * condprob
)

We integrate over omega using Monte-Carlo integration

loglike = log(MonteCarlo(condlike))

Read the data

database = read_data()

As the objective is to illustrate the syntax, we calculate the Monte-Carlo approximation with a small number of draws. Therefore, we first define the parameters.

the_biogeme = bio.BIOGEME(database, loglike, number_of_draws=100, seed=1223)
the_biogeme.modelName = 'b06serial_correlation'

Biogeme parameters read from biogeme.toml.

If estimation results are saved on file, we read them to speed up the process. If not, we estimate the parameters.

results = read_or_estimate(the_biogeme=the_biogeme, directory='saved_results')

print(f'Final log likelihood: {results.data.logLike:.3f}')
print(f'Output file: {results.data.htmlFileName}')

Final log likelihood: -18406.350
Output file: b06serial_correlation.html

results.get_estimated_parameters()

	Value	Rob. Std err	Rob. t-test	Rob. p-value
ASC_CAR	0.689343	0.117733	5.855133	4.766288e-09
ASC_SM	0.207845	0.340093	0.611144	5.411045e-01
BETA_COST_HWH	-1.386213	0.332408	-4.170214	3.043140e-05
BETA_COST_OTHER	-0.506482	0.159352	-3.178374	1.481036e-03
BETA_DIST	-1.394661	0.370550	-3.763757	1.673794e-04
BETA_TIME_CAR_CL	-0.931835	0.155398	-5.996450	2.016773e-09
BETA_TIME_CAR_REF	-9.236444	1.880464	-4.911790	9.024891e-07
BETA_TIME_PT_CL	-0.429495	0.113183	-3.794700	1.478222e-04
BETA_TIME_PT_REF	-3.171013	0.807718	-3.925894	8.640831e-05
BETA_WAITING_TIME	-0.019410	0.009485	-2.046372	4.071975e-02
B_Envir02_F1	-0.457658	0.031340	-14.602895	0.000000e+00
B_Envir03_F1	0.485267	0.031978	15.175026	0.000000e+00
B_Mobil11_F1	0.569761	0.041931	13.588143	0.000000e+00
B_Mobil14_F1	0.578074	0.035400	16.329817	0.000000e+00
B_Mobil16_F1	0.530125	0.042462	12.484542	0.000000e+00
B_Mobil17_F1	0.514992	0.041891	12.293621	0.000000e+00
INTER_Envir02	0.456798	0.030753	14.853833	0.000000e+00
INTER_Envir03	-0.365668	0.028703	-12.739868	0.000000e+00
INTER_Mobil11	0.417124	0.036801	11.334545	0.000000e+00
INTER_Mobil14	-0.173064	0.027854	-6.213324	5.187542e-10
INTER_Mobil16	0.144486	0.033413	4.324182	1.530989e-05
INTER_Mobil17	0.138232	0.032669	4.231305	2.323389e-05
SIGMA_STAR_Envir02	0.915628	0.034011	26.921843	0.000000e+00
SIGMA_STAR_Envir03	0.854036	0.034954	24.433443	0.000000e+00
SIGMA_STAR_Mobil11	0.894418	0.040472	22.099739	0.000000e+00
SIGMA_STAR_Mobil14	0.757393	0.032915	23.010639	0.000000e+00
SIGMA_STAR_Mobil16	0.868209	0.039182	22.158210	0.000000e+00
SIGMA_STAR_Mobil17	0.873779	0.038761	22.542752	0.000000e+00
beta_ScaledIncome_10_inf	0.124306	0.037808	3.287789	1.009774e-03
beta_ScaledIncome_4_6	-0.300850	0.118960	-2.529008	1.143855e-02
beta_ScaledIncome_6_8	0.344675	0.140951	2.445359	1.447081e-02
beta_ScaledIncome_8_10	-0.680656	0.160339	-4.245112	2.184846e-05
beta_ScaledIncome_minus_inf_4	0.149449	0.061751	2.420182	1.551273e-02
coef_age_65_more	0.052404	0.080149	0.653838	5.132161e-01
coef_haveChildren	-0.015983	0.059651	-0.267947	7.887404e-01
coef_haveGA	-0.733132	0.097889	-7.489436	6.927792e-14
coef_highEducation	-0.259276	0.068751	-3.771225	1.624483e-04
coef_individualHouse	-0.140217	0.057695	-2.430321	1.508546e-02
coef_intercept	0.388138	0.168012	2.310176	2.087842e-02
coef_male	0.068427	0.055521	1.232465	2.177756e-01
coef_moreThanOneBike	-0.387161	0.071862	-5.387521	7.143610e-08
coef_moreThanOneCar	0.722775	0.067049	10.779820	0.000000e+00
delta_1	0.326445	0.012506	26.103383	0.000000e+00
delta_2	0.984414	0.035152	28.004450	0.000000e+00
ec_sigma	-0.302634	0.222164	-1.362210	1.731316e-01
sigma_s	0.801003	0.096041	8.340198	0.000000e+00

Total running time of the script: (0 minutes 0.796 seconds)

Gallery generated by Sphinx-Gallery