**AUTEUR : Isabelle CADORET-DAVID**


ALL 840
OPEN DATA choixtransport.wks
DATA(FORMAT=wks,ORG=obs) /
*PRINT /

SET modechoix1 = 0.0
SET modechoix2 = 0.0
SET modechoix3 = 0.0
SET modechoix4 = 0.0

***codage de la série mode en 1,2,3,4 selon le choix de l'individu***
COMPUTE nbc = 4      ; *nombre d'alternatives*
COMPUTE nin = 210    ; *le nombre d'individus*
COMPUTE tot = nbc*nin

DO I=1,tot,nbc
COM modechoix1(i) = 1.0
COM modechoix2(i+1) = 1.0
COM modechoix3(i+2) = 1.0
COM modechoix4(i+3) = 1.0
END DO I

*PRINT 1 30  modechoix1 modechoix2 modechoix3 modechoix4

SET modair = modechoix1*mode
SET modbus = modechoix2*mode*2
SET modtrain = modechoix3*mode*3
SET modcar = modechoix4*mode*4


SET mod = modair+modbus+modtrain+modcar


SET selec = %if(mod>0,1,0); *selection des 210 individus, mod indiquant leur choix*
SMPL(series=selec)


DISP '*****************************************************************************'
DISP '**Estimation du premier modèle******'
DISP '*****************************************************************************'

DEC VEC  air(6) bus(6) train(6) car(6)
NONLIN air bus train car

FRML zair = air(1)*ttme+air(2)*invc+air(3)*invt+air(4)*gc+air(5)*hinc+air(6)
FRML zbus = bus(1)*ttme+bus(2)*invc+bus(3)*invt+bus(4)*gc+bus(5)*hinc+bus(6)
FRML ztrain = train(1)*ttme+train(2)*invc+train(3)*invt+train(4)*gc+train(5)*hinc+train(6)
FRML zcar = car(1)*ttme+car(2)*invc+car(3)*invt+car(4)*gc+car(5)*hinc+car(6)

DISP '***estimation en prenant le mode de transport en voiture en référence***'
FRML multilgv = (f1=zair), (f2=zbus), (f3=ztrain), $
                %IF(mod==1,f1,%IF(mod==2,f2,%IF(mod==3,f3,0.0)))-log(1.0+exp(f1)+exp(f2)+exp(f3))
COMPUTE air=%const(0.0), train=%const(0.0),bus=%const(0.0),car=%const(0.0)
MAXIMIZE(METHOD=BFGS,iters=100) multilgv

DISP '***estimation en prenant le mode de transport en avion en référence***'
FRML multilga =  (f2=zbus), (f3=ztrain), (f4=zcar),$
                %IF(mod==4,f4,%IF(mod==3,f3,%IF(mod==2,f2,0.0)))-log(1.0+exp(f2)+exp(f3)+exp(f4))
COMPUTE air=%const(0.0), train=%const(0.0),bus=%const(0.0),car=%const(0.0)
MAXIMIZE(METHOD=BFGS,iters=100) multilga



DISP '*****************************************************************************'
DISP '**Partie 1 : Estimation d"un modèle multinomial logit******'
DISP '*****************************************************************************'


DISP '*****************************************************************************'
DISP '**Estimation avec les caractéristiques des individus******'
DISP '*****************************************************************************'

DEC VEC  air(2) bus(2) train(2) car(2)
NONLIN air bus train car

FRML zair = air(1)+air(2)*hinc
FRML zbus = bus(1)+ bus(2)*hinc
FRML ztrain = train(1)+ train(2)*hinc
FRML zcar = car(1)+car(2)*hinc


DISP '***estimation en prenant le mode de transport en voiture en référence***'
FRML multilgv = (f1=zair), (f2=zbus), (f3=ztrain), $
                %IF(mod==1,f1,%IF(mod==2,f2,%IF(mod==3,f3,0.0)))-log(1.0+exp(f1)+exp(f2)+exp(f3))
COMPUTE air=%const(0.0), train=%const(0.0),bus=%const(0.0),car=%const(0.0)
MAXIMIZE(METHOD=BFGS,smpl=selec,iters=100) multilgv

DISP '***estimation en prenant le mode de transport en avion en référence***'
FRML multilga =  (f2=zbus), (f3=ztrain), (f4=zcar),$
                %IF(mod==4,f4,%IF(mod==3,f3,%IF(mod==2,f2,0.0)))-log(1.0+exp(f2)+exp(f3)+exp(f4))
COMPUTE air=%const(0.0), train=%const(0.0),bus=%const(0.0),car=%const(0.0)
MAXIMIZE(METHOD=BFGS,smpl=selec,iters=100) multilga

DISP '***Calcul des probabilités estimées***'
SET zairp = air(1)+air(2)*hinc
SET zbusp = bus(1)+ bus(2)*hinc
SET ztrainp = train(1)+ train(2)*hinc
SET zcarp = car(1)+car(2)*hinc


SET pair = 1 / (1+exp(zcarp)+exp(zbusp)+exp(ztrainp))
SET pbus =  exp(zbusp)/(1+exp(zcarp)+exp(zbusp)+exp(ztrainp))
SET ptrain = exp(ztrainp)/(1+exp(zcarp)+exp(zbusp)+exp(ztrainp))
SET pcar = exp(zcarp)/(1+exp(zcarp)+exp(zbusp)+exp(ztrainp))
SET ptot = pair + pbus + ptrain + pcar ;*vérification que la somme est égale à 1**
PRINT 1 30 pair pbus ptrain pcar ptot

DISP '***Calcul des effets marginaux pour un niveau moyen des probabilités***'
STAT(noprint) pcar
COMP pcarm = %mean
STAT(noprint) pbus
COMP pbusm = %mean
STAT(noprint) ptrain
COMP ptrainm = %mean
STAT(noprint) pair
COMP pairm = %mean

COMP pcarhinc = pcarm*(car(2)-(pcarm*car(2)+pbusm*bus(2)+ptrainm*train(2)))
COMP pbushinc = pbusm*(bus(2)-(pcarm*car(2)+pbusm*bus(2)+ptrainm*train(2)))
COMP ptrainhinc = ptrainm*(train(2)-(pcarm*car(2)+pbusm*bus(2)+ptrainm*train(2)))
COMP pairhinc = pairm*(-car(2)-(pairm*(-car(2))+pbusm*(bus(2)-car(2))+ptrainm*(train(2)-car(2))))
COMP effettot =  pcarhinc+pbushinc+ptrainhinc+pairhinc

DISP 'effets marginaux sur pcar, pbus, ptrain'
DISP  pcarhinc   pbushinc    ptrainhinc pairhinc
DISP 'somme des effets ' effettot

DISP '***Calcul des élasticités***'
STAT(noprint) hinc
COMP hincm = %mean
COMP pcarelashinc = pcarhinc*(hincm/pcarm)
COMP pbuselashinc = pbushinc*(hincm/pbusm)
COMP ptrainelashinc = ptrainhinc*(hincm/ptrainm)
DISP 'élasticités sur pcar, pbus, ptrain'
DISP pcarelashinc pbuselashinc ptrainelashinc

DISP '*****************************************************************************'
DISP '**Partie 2 : Estimation d"un modèle logit conditionel******'
DISP '*****************************************************************************'

SMPL
SET gcair = 0.0
SET gcbus = 0.0
SET gctrain = 0.0
SET gccar = 0.0

SET ttair = 0.0
SET ttbus = 0.0
SET tttrain = 0.0
SET ttcar = 0.0

SET invcair = 0.0
SET invcbus = 0.0
SET invctrain = 0.0
SET invccar = 0.0

SET invtair = 0.0
SET invtbus = 0.0
SET invttrain = 0.0
SET invtcar = 0.0

SET moda = 0.0
SET modb = 0.0
SET modt = 0.0
SET modc = 0.0
SET dum  = 0.0

DO I=1,tot,nbc
COM gcair(i) = gc(i)
COM gcbus(i) = gc(i+1)
COM gctrain(i) = gc(i+2)
COM gccar(i) = gc(i+3)

COM ttair(i) = ttme(i)
COM ttbus(i) = ttme(i+1)
COM tttrain(i) = ttme(i+2)
COM ttcar(i) = ttme(i+3)

COM invcair(i) = invc(i)
COM invcbus(i) = invc(i+1)
COM invctrain(i) = invc(i+2)
COM invccar(i) = invc(i+3)

COM invtair(i) = invt(i)
COM invtbus(i) = invt(i+1)
COM invttrain(i) = invt(i+2)
COM invtcar(i) = invt(i+3)

COM moda(i) = modechoix1(i)*mode(i)
COM modb(i) = modechoix2(i+1)*mode(i+1)
COM modt(i) = modechoix3(i+2)*mode(i+2)
COM modc(i) = modechoix4(i+3)*mode(i+3)

COM dum(i) = 1.0
END DO I

SMPL(series=dum)

SET modca = moda*1+modb*2+modt*3+modc*4

DISP '********************************************************************************************************'
DISP '*******Estimation du modèle en fonction des variables dont la valeur dépend du choix de l"individu******'
DISP '*********************************************************************************************************'


DISP '***estimation en prenant le mode de transport en voiture en référence***'

DEC VEC coef(4)
NONLIN coef
FRML ziair = coef(1)*(gcair-gccar)+coef(2)*(ttair-ttcar)+coef(3)*(invcair-invccar)+coef(4)*(invtair-invtcar)
FRML zibus = coef(1)*(gcbus-gccar)+coef(2)*(ttbus-ttcar)+coef(3)*(invcbus-invccar)+coef(4)*(invtbus-invtcar)
FRML zitrain =coef(1)*(gctrain-gccar)+coef(2)*(tttrain-ttcar)+coef(3)*(invctrain-invccar)+coef(4)*(invttrain-invtcar)

FRML condlogit = (c1=ziair), (c2=zibus), (c3=zitrain),  $
                %IF(modca==1,c1,%IF(modca==2,c2,%IF(modca==3,c3,0)))-log(exp(c1)+exp(c2)+exp(c3)+1)
COMPUTE coef=%const(0.0)
MAXIMIZE(METHOD=BFGS,iters=100) condlogit

DISP '***calcul des probabilités estilées en prenant le mode de transport en voiture en référence***'

SET ziairp = coef(1)*(gcair-gccar)+coef(2)*(ttair-ttcar)+coef(3)*(invcair-invccar)+coef(4)*(invtair-invtcar)
SET zibusp = coef(1)*(gcbus-gccar)+coef(2)*(ttbus-ttcar)+coef(3)*(invcbus-invccar)+coef(4)*(invtbus-invtcar)
SET zitrainp =coef(1)*(gctrain-gccar)+coef(2)*(tttrain-ttcar)+coef(3)*(invctrain-invccar)+coef(4)*(invttrain-invtcar)

SET pcar = 1 / (1+exp(ziairp)+exp(zibusp)+exp(zitrainp))
SET pbus =  exp(zibusp)/(1+exp(ziairp)+exp(zibusp)+exp(zitrainp))
SET ptrain = exp(zitrainp)/(1+exp(ziairp)+exp(zibusp)+exp(zitrainp))
SET pair = exp(ziairp)/(1+exp(ziairp)+exp(zibusp)+exp(zitrainp))
SET ptot = pair + pbus + ptrain + pcar ;*vérification que la somme est égale à 1**
PRINT 1 30 pair pbus ptrain pcar ptot

DISP '***estimation en prenant le mode de transport en avion en référence***'

DEC VEC coef(4)
NONLIN coef
FRML zicar = coef(1)*(gccar-gcair)+coef(2)*(ttcar-ttair)+coef(3)*(invccar-invcair)+coef(4)*(invtcar-invtair)
FRML zibus = coef(1)*(gcbus-gcair)+coef(2)*(ttbus-ttair)+coef(3)*(invcbus-invcair)+coef(4)*(invtbus-invtair)
FRML zitrain =coef(1)*(gctrain-gcair)+coef(2)*(tttrain-ttair)+coef(3)*(invctrain-invcair)+coef(4)*(invttrain-invtair)


FRML condlogit = (c4=zicar), (c2=zibus), (c3=zitrain),  $
                %IF(modca==4,c4,%IF(modca==2,c2,%IF(modca==3,c3,0)))-log(exp(c4)+exp(c2)+exp(c3)+1.0)
COMPUTE coef=%const(0.0)
MAXIMIZE(METHOD=BFGS,iters=100) condlogit

DISP '***Calcul des probabilités estimées en prenant le mode de transport en avion en référence***'

SET zicarp = coef(1)*(gccar-gcair)+coef(2)*(ttcar-ttair)+coef(3)*(invccar-invcair)+coef(4)*(invtcar-invtair)
SET zibusp = coef(1)*(gcbus-gcair)+coef(2)*(ttbus-ttair)+coef(3)*(invcbus-invcair)+coef(4)*(invtbus-invtair)
SET zitrainp =coef(1)*(gctrain-gcair)+coef(2)*(tttrain-ttair)+coef(3)*(invctrain-invcair)+coef(4)*(invttrain-invtair)


SET pair = 1 / (1+exp(zicarp)+exp(zibusp)+exp(zitrainp))
SET pbus =  exp(zibusp)/(1+exp(zicarp)+exp(zibusp)+exp(zitrainp))
SET ptrain = exp(zitrainp)/(1+exp(zicarp)+exp(zibusp)+exp(zitrainp))
SET pcar = exp(zicarp)/(1+exp(zicarp)+exp(zibusp)+exp(zitrainp))
SET ptot = pair + pbus + ptrain + pcar ;*vérification que la somme est égale à 1**
PRINT 1 30 pair pbus ptrain pcar ptot

DISP '***Calcul des effets marginaux pour un niveau moyen la probabilité de prendre le bus***'

STAT(noprint) pcar
COMP pcarm = %mean
STAT(noprint) pbus
COMP pbusm = %mean
STAT(noprint) ptrain
COMP ptrainm = %mean

**mesure de l'effet d'une variation des variables gcbus, gccar, gctrain**
COM pbusgcbus = coef(1)*pbusm*(1-pbusm)
COM pbusgctrain = -coef(1)*pbusm*ptrainm
COM pbusgccar = - coef(1)*pbusm*pcarm

**mesure de l'effet d'une variation des variables ttbus, ttcar, tttrain**
COM pbusttbus = coef(2)*pbusm*(1-pbusm)
COM pbustttrain = -coef(2)*pbusm*ptrainm
COM pbusttcar = - coef(2)*pbusm*pcarm

**mesure de l'effet d'une variation des variables invcbus, invccar, invctrain**
COM pbusinvcbus = coef(3)*pbusm*(1-pbusm)
COM pbusinvctrain = -coef(3)*pbusm*ptrainm
COM pbusinvccar = - coef(3)*pbusm*pcarm

**mesure de l'effet d'une variation des variables invtbus, invtcar, invttrain**
COM pbusinvtbus = coef(3)*pbusm*(1-pbusm)
COM pbusinvttrain = -coef(3)*pbusm*ptrainm
COM pbusinvtcar = - coef(3)*pbusm*pcarm

DISP  'effet de gcbus sur pbus' pbusgcbus
DISP  'effet de gctrain sur pbus' pbusgctrain
DISP  'effet de gccar sur pbus' pbusgccar


DISP  'effet de ttbus sur pbus' pbusttbus
DISP  'effet de tttrain sur pbus' pbustttrain
DISP  'effet de ttcar sur pbus' pbusttcar


DISP  'effet de invcbus sur pbus' pbusinvcbus
DISP  'effet de invctrain sur pbus' pbusinvctrain
DISP  'effet de invccar sur pbus' pbusinvccar

DISP  'effet de invtbus sur pbus' pbusinvtbus
DISP  'effet de invttrain sur pbus' pbusinvttrain
DISP  'effet de invtcar sur pbus' pbusinvtcar

DISP '***Calcul des effets marginaux de gcbus sur pbus ptrain pcar pair***'
COM ptraingcbus = -coef(1)*ptrainm*pbusm
COM pcargcbus = -coef(1)*pcarm*pbusm
COM pairbus = -coef(1)*pairm*pbusm
COM effettot = pbusgcbus + ptraingcbus + pcargcbus + pairbus
DISP pbusgcbus ptraingcbus  pcargcbus  pairbus
DISP 'somme des effets' effettot

DISP '*****************************************************************************'
DISP '**Partie 3 : Combinaison des deux modèles******'
DISP '*****************************************************************************'

DISP '********************************************************************************************************'
DISP '**Estimation du modèle en fonction des variables dont la valeur dépend du choix de l"individu******'
DISP '**et des caractéristiques individuelles**'
DISP '*********************************************************************************************************'


DISP '***estimation en prenant le mode de transport en voiture en référence***'
DEC VEC coef(4) air(2) bus(2) train(2) car(2)
NONLIN coef air bus train car

FRML zcair = coef(1)*(gcair-gccar)+coef(2)*(ttair-ttcar)+coef(3)*(invcair-invccar)+coef(4)*(invtair-invtcar) $
              +air(1)+air(2)*hinc
FRML zcbus = coef(1)*(gcbus-gccar)+coef(2)*(ttbus-ttcar)+coef(3)*(invcbus-invccar)+coef(4)*(invtbus-invtcar) $
             +bus(1)+bus(2)*hinc

FRML zctrain =coef(1)*(gctrain-gccar)+coef(2)*(tttrain-ttcar)+coef(3)*(invctrain-invccar)+coef(4)*(invttrain-invtcar)  $
		+train(1)+train(2)*hinc


FRML clogit = (c1=zcair), (c2=zcbus), (c3=zctrain),  $
                %IF(modca==1,c1,%IF(modca==2,c2,%IF(modca==3,c3,0)))-log(exp(c1)+exp(c2)+exp(c3)+1)
COMPUTE coef=%const(0.0), air=%const(0.0), train=%const(0.0), bus=%const(0.0), car=%const(0.0)
MAXIMIZE(METHOD=BFGS,iters=100) clogit


DISP '***estimation en prenant le mode de transport en avion en référence***'
DEC VEC coef(4) air(2) bus(2) train(2) car(2)
NONLIN coef air bus train car

FRML zccar = coef(1)*(gccar-gcair)+coef(2)*(ttcar-ttair)+coef(3)*(invccar-invcair)+coef(4)*(invtcar-invtair) $
              +car(1)+car(2)*hinc
FRML zcbus = coef(1)*(gcbus-gcair)+coef(2)*(ttbus-ttair)+coef(3)*(invcbus-invcair)+coef(4)*(invtbus-invtair) $
             +bus(1)+bus(2)*hinc

FRML zctrain =coef(1)*(gctrain-gcair)+coef(2)*(tttrain-ttair)+coef(3)*(invctrain-invcair)+coef(4)*(invttrain-invtair)  $
		+train(1)+train(2)*hinc


FRML clogit = (c4=zccar), (c2=zcbus), (c3=zctrain),  $
                %IF(modca==4,c4,%IF(modca==2,c2,%IF(modca==3,c3,0)))-log(exp(c4)+exp(c2)+exp(c3)+1)
COMPUTE coef=%const(0.0), air=%const(0.0), train=%const(0.0), bus=%const(0.0), car=%const(0.0)
MAXIMIZE(METHOD=BFGS,iters=100) clogit

DISP '***Calcul des probabilités estimées***'

SET zccarp = coef(1)*(gccar-gcair)+coef(2)*(ttcar-ttair)+coef(3)*(invccar-invcair)+coef(4)*(invtcar-invtair) $
              +car(1)+car(2)*hinc
SET zcbusp = coef(1)*(gcbus-gcair)+coef(2)*(ttbus-ttair)+coef(3)*(invcbus-invcair)+coef(4)*(invtbus-invtair) $
             +bus(1)+bus(2)*hinc

SET zctrainp =coef(1)*(gctrain-gcair)+coef(2)*(tttrain-ttair)+coef(3)*(invctrain-invcair)+coef(4)*(invttrain-invtair)  $
		+train(1)+train(2)*hinc



SET pair = 1 / (1+exp(zccarp)+exp(zcbusp)+exp(zctrainp))
SET pbus =  exp(zcbusp)/(1+exp(zccarp)+exp(zcbusp)+exp(zctrainp))
SET ptrain = exp(zctrainp)/(1+exp(zccarp)+exp(zcbusp)+exp(zctrainp))
SET pcar = exp(zccarp)/(1+exp(zccarp)+exp(zcbusp)+exp(zctrainp))
SET ptot = pair + pbus + ptrain + pcar ;*vérification que la somme est égale à 1**
PRINT 1 30 pair pbus ptrain pcar ptot

DISP '***Calcul des effets marginaux pour un niveau moyen la probabilité de prendre le bus***'

STAT(noprint) pcar
COMP pcarm = %mean
STAT(noprint) pbus
COMP pbusm = %mean
STAT(noprint) ptrain
COMP ptrainm = %mean

**mesure de l'effet d'une variation des variables gcbus, gccar, gctrain et de hinc**
COM pbusgcbus = coef(1)*pbusm*(1-pbusm)
COM pbusgctrain = -coef(1)*pbusm*ptrainm
COM pbusgccar = - coef(1)*pbusm*pcarm
COM pbushinc = pbusm*(bus(2)-(pcarm*car(2)+pbusm*bus(2)+ptrainm*train(2)))

DISP  'effet de gcbus sur pbus' pbusgcbus
DISP  'effet de gctrain sur pbus' pbusgctrain
DISP  'effet de gccar sur pbus' pbusgccar
DISP  'effet de hinc sur pbus' pbushinc