new;

load x[1522,55]=c:\tokyo1\ensyu2.csv;

rows(x);

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:\tokyo1\report.txt on;


__output=2;

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

proc li(b,x);


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


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


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

para=b;

retp(ln(prob));


endp;

let _max_parnames=
"rconst" 
"bconst"
"cconst"
"bconst" 
"cost" 
@"c_cost"@
"time" 
"c_time"
"Btime"
"Wtime"
"h_max"
"otoko"
"30"
"k_50"
"50-"
"h_5"
"k_30a"
@"ame"@
"c_cost"
@"r_hare" @

@
"holiday"
"shop"
@
@
"jdansei"
"cabe50 "
"job"    
"shop"  
"age"
@
;



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

selut=(x[.,19].==200).*rail+(x[.,19].==240).*bus+(x[.,19].==100).*car
+(x[.,19].==410).*bicycle+(x[.,19].==420).*walk;
ut=rail+bus+car+bicycle+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;




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


sumc(((x[.,20].==1).*rail./ut)~((x[.,21].==1).*bus./ut)~((x[.,22].==1).*car./ut))/rows(x);


output off;


end;@