new;

load x[1522,56]=c:\tokyo\ensyu2.csv;

rows(x);

@e1=x[.,3].==200;@
e1=x[.,1].==256;
e2=x[.,3].==800;
e3=x[.,5].==-999;
e4=sumc((x[.,9:11].==-999)');
e5=sumc((x[.,14:16].==-999)');
e6=x[.,17].==0;
e7=x[.,19].>=500;
e8=(x[.,22].==-999).or(x[.,22].<=10);
e9=x[.,30].==-999;
e10=(x[.,25].==0).*(x[.,19].==200)+(x[.,41].==0).*(x[.,19].==240)+(x[.,23].==0).*(x[.,19].==100)+(x[.,35].==0).*(x[.,19].==410)+(x[.,38].==0).*(x[.,19].==420);
@e11=x[.,35].<=3;@
@e12=x[.,38].<=3;@
@e13=x[.,23].<=5;@
@e14=x[.,25].<=5;@
@e15=x[.,41].<=5;@
@e16=((x[.,34].>5000).*(x[.,19].==410));@
@e17=((x[.,37].>5000).*(x[.,19].==420));@
@e18=x[.,40].<=10;@



e=e1+e2+e3+e4+e5+e6+e7+e8+e9+e10@+e11+e12+e13+e14+e15+e16+e17+e18@;
x=delif(x,e.>0);

@x=selif(x,x[.,37].<=5000);@

rows(x);

library maxlik;
#include maxlik.ext;
maxset;

output file=c:\tokyo\report1.txt on;


__output=2;

clearg prob,selut,ut,rail,bus,car,bicycle,walk,para;

proc li(b,x);


@*****MNL*****@

@
rail     =(x[.,32]./=0).*exp(                     b[1]+b[5]*(x[.,26]/1000).*(x[.,55].>=0.3) +b[6]*(x[.,26]/1000).*(1-(x[.,55].>=0.3))           +b[7]*((x[.,25]+x[.,29]+x[.,31])/60));
bus    =(x[.,47]./=0).*(x[.,42]./=0).*exp( b[2]+b[5]*(x[.,42]/1000).*(x[.,55].>=0.3) +b[6]*(x[.,42]/1000).*(1-(x[.,55].>=0.3))           +b[7]*((x[.,41]+x[.,44]+x[.,46])/60));
car    =(x[.,23]./=0).*exp(                     b[3]+b[14]*(x[.,22]/10000*0.15)                                                                                          +b[8]*(x[.,23]/60)+b[11]*(x[.,18].==1) +b[14]*(x[.,3] .==500));
bicycle=(x[.,35]./=0).*exp(                     b[4]                                                                                                                                   +b[9]*(x[.,35]/60)+b[12]*(x[.,17].<=30));
walk   =(x[.,38]./=0).*exp(                                                                                                                                                               b[10]*(x[.,38]/60)+b[13]*(x[.,17].>=50)+b[14]*(x[.,51] .>=1));
@

rail     =(x[.,32]./=0).*exp(                     b[1]+b[5]*(x[.,26]/1000).*(x[.,55].>=0.3) +b[6]*(x[.,26]/1000).*(1-(x[.,55].>=0.3)) +b[7]*x[.,25]/60+b[15].*((x[.,29]+x[.,31])/60).*(x[.,50].>=2)+b[16].*((x[.,29]+x[.,31])/60).*(1-(x[.,50].>=2)));
bus    =(x[.,47]./=0).*(x[.,42]./=0).*exp( b[2]+b[5]*(x[.,42]/1000).*(x[.,55].>=0.3) +b[6]*(x[.,42]/1000).*(1-(x[.,55].>=0.3)) +b[7]*x[.,41]/60+b[15].*((x[.,44]+x[.,46])/60).*(x[.,50].>=2)+b[16].*((x[.,44]+x[.,46])/60).*(1-(x[.,50].>=2)));
car    =(x[.,23]./=0).*exp(                     b[3]+b[14]*(x[.,22]/10000*0.127)                                                                      +b[8]*(x[.,23]/60)+b[11]*(x[.,18].==1));
bicycle=(x[.,35]./=0).*exp(                     b[4]                                                                                                                 +b[9]*(x[.,35]/60)+b[12]*(x[.,17].<=30)+b[17]*(x[.,51].<=0.5).*(x[.,55].>=0.5));
walk   =(x[.,38]./=0).*exp(                                                                                                                                             b[10]*(x[.,38]/60)+b[13]*(x[.,17].>=50)+b[17]*(x[.,51].<=0.5).*(x[.,55].>=0.5));






@*****Probability*****@
selut=(x[.,19].==200).*rail+(x[.,19].==240).*bus+(x[.,19].==100).*car
+(x[.,19].==410).*bicycle+(x[.,19].==420).*walk;

ut=(x[.,24].==1).*rail+(x[.,39].==1).*bus+(x[.,21].==1).*car+(x[.,33].==1).*bicycle+(x[.,36].==1).*walk;
prob=selut./ut;

para=b;

retp(ln(prob));


endp;

let _max_parnames=
"rconst" 
"bconst"
"cconst"
"jconst" 
"cost_r" 
"cost_nr"
"ptime" 
"Ctime"
"Btime"
"Wtime"
"male"
"age30"
"age50"
@"ab5mm"@
"costcar"
"acig_2m"
"acig_0m"
"ame"
;



start=zeros(17,1);

{b,ff,gg,cov,retcode}=maxlik(x,0,&li,start);
call maxprt(b,ff,gg,cov,retcode);



sumc(((x[.,24].==1).*rail./ut)~((x[.,39].==1).*bus./ut)~((x[.,21].==1).*car./ut)~((x[.,33].==1).*bicycle./ut)~((x[.,36].==1).*walk./ut))/rows(x);

Lb=sumc(li(para,x));
L0=sumc(li(start,x));


rail     =(x[.,32]./=0).*exp(                     b[1]+b[5]*(x[.,26]/1000).*(x[.,55].>=0.3)*0.1 +b[6]*(x[.,26]/1000).*(1-(x[.,55].>=0.3)) +b[7]*x[.,25]/60+b[15].*((x[.,29]+x[.,31])/60).*(x[.,50].>=2)+b[16].*((x[.,29]+x[.,31])/60).*(1-(x[.,50].>=2)));
bus    =(x[.,47]./=0).*(x[.,42]./=0).*exp( b[2]+b[5]*(x[.,42]/1000).*(x[.,55].>=0.3)*0.1 +b[6]*(x[.,42]/1000).*(1-(x[.,55].>=0.3)) +b[7]*x[.,41]/60+b[15].*((x[.,44]+x[.,46])/60).*(x[.,50].>=2)+b[16].*((x[.,44]+x[.,46])/60).*(1-(x[.,50].>=2)));
car    =(x[.,23]./=0).*exp(                     b[3]+b[14]*(x[.,22]/10000*0.127)                                                                      +b[8]*(x[.,23]/60)+b[11]*(x[.,18].==1));
bicycle=(x[.,35]./=0).*exp(                     b[4]                                                                                                                 +b[9]*(x[.,35]/60)+b[12]*(x[.,17].<=30)+b[17]*(x[.,51].<=0.5).*(x[.,55].>=0.5));
walk   =(x[.,38]./=0).*exp(                                                                                                                                             b[10]*(x[.,38]/60)+b[13]*(x[.,17].>=50)+b[17]*(x[.,51].<=0.5).*(x[.,55].>=0.5));


selut=(x[.,19].==200).*rail+(x[.,19].==240).*bus+(x[.,19].==100).*car
+(x[.,19].==410).*bicycle+(x[.,19].==420).*walk;


ut=(x[.,24].==1).*rail+(x[.,39].==1).*bus+(x[.,21].==1).*car+(x[.,33].==1).*bicycle+(x[.,36].==1).*walk;
prob=selut./ut;

sumc(((x[.,24].==1).*rail./ut)~((x[.,39].==1).*bus./ut)~((x[.,21].==1).*car./ut)~((x[.,33].==1).*bicycle./ut)~((x[.,36].==1).*walk./ut))/rows(x);


print "p^2";(L0-Lb)/L0;
print "p^2_bar";(L0-(Lb-rows(start)))/L0;



output off;
end;