@   ABGBNP.CMD Aug 4, 2004 -- A-B choice and frequency data, observed GB
days, and nonparticipation @ 

new;

library maxlik;
maxset;

clear xdif;
clearg beta,V0,V0r,p,e1,e2,rho,p0,lnp,lnL,lnr,pnr,sigmae,soe,sog;

output file = abgbnp.out reset;

"ABGBNP.CMD";

let kvec = 1 2 3 4 5 6 7 8 9 10 11 12 13;

loadm Z2[647, 28] = Z2.dat;
n = rows(Z2);

"Number of observations: " n;

@ Number of occasions@
oc=50;

loadm X[160, 13] = X.dat;
X=X[.,1:4]/10~X[.,5:13];

Xdif=zeros(8*n,rows(kvec));
Xc=zeros(8*n,rows(kvec));
Xnc=zeros(8*n,rows(kvec));
tps=zeros(8*n,1);
gend=zeros(8*n,1);
age=zeros(8*n,1);
ret=zeros(8*n,1);
boat=zeros(8*n,1);
reg1=zeros(8*n,1);
reg2=zeros(8*n,1);
tcc=zeros(8*n,1);
tcnc=zeros(8*n,1);
tot1=oc*ones(8*n,1); @ number of total occasions for SP portion @

ng = z2[.,21]; @ actual GB trips @
dist = z2[.,18]; @ closest distance to GB @
income = z2[.,27];
county = z2[.,28];
tot2=oc*ones(n,1); @ number of total occasions for RP portion @
inc1 = zeros(n,1); 
reg1_2=zeros(n,1);
reg2_2=zeros(n,1);

q=0;
do while q<n; q=q+1; 
        if income[q] eq 1; inc1[q] = 5000; 
        elseif income[q] eq 2; inc1[q] = 15000;
        elseif income[q] eq 3; inc1[q] = 25000;
        elseif income[q] eq 4; inc1[q] = 35000;
        elseif income[q] eq 5; inc1[q] = 45000;
        elseif income[q] eq 6; inc1[q] = 55000;
        elseif income[q] eq 7; inc1[q] = 70000;
        elseif income[q] eq 8; inc1[q] = 90000;
        elseif income[q] eq 9; inc1[q] = 125000;
        elseif income[q] eq 10; inc1[q] = 150000; endif;
endo;

r=0;
do while r<n; r=r+1;
        if county[r] eq 2; reg1_2[r] = 1; endif;
        if county[r] eq 3; reg2_2[r] = 1;
        elseif county[r] eq 4; reg2_2[r] = 1; 
        elseif county[r] eq 5; reg2_2[r] = 1; endif;
endo;

ll=0; i=0;
do while i<n; i=i+1;
  j=0;
  do while j<8;     j=j+1; ll=ll+1;
    c=z2[i,j];    @ choice made by i in pair j (1 or 2) @
    t=z2[i,j+8];  @ GB trips that would be associated with that choice @
    d=z2[i,18];  @ nearest distance from home or cabin to GB launch @
    v=z2[i,17];    @ survey version for individual i @
    inc = z2[i,27]; @ income @
        if inc eq 1; inc = 5000; 
        elseif inc eq 2; inc = 15000;
        elseif inc eq 3; inc = 25000;
        elseif inc eq 4; inc = 35000;
        elseif inc eq 5; inc = 45000;
        elseif inc eq 6; inc = 55000;
        elseif inc eq 7; inc = 70000;
        elseif inc eq 8; inc = 90000;
        elseif inc eq 9; inc = 125000;
        elseif inc eq 10; inc = 150000; endif;
    county = z2[i,28]; @ county @
        if county eq 2; reg1[ll] = 1; endif;
        if county eq 3; reg2[ll] = 1; 
        elseif county eq 4; reg2[ll] = 1; 
        elseif county eq 5; reg2[ll] = 1; endif;
    gend[ll]=z2[i,23];
    age[ll]=z2[i,24];
    ret[ll]=z2[i,25];
    boat[ll]=z2[i,26];
 
    cind=16*(v-1)+2*(j-1)+c;
    ncind=16*(v-1)+2*(j-1)+3-c; 
    
if c<1; Xc[ll,.]=zeros(1,13);
    else; Xc[ll,.]=x[cind,kvec]; endif;
if c<1; Xc[11,13]=0;
    else; Xc[ll,13]=((inc*0.5/2080*1/3*(2*d)/45)+0.26*(2*d*0.5)+Xc[ll,13])/100; endif;

if c<1; Xnc[ll,.]=zeros(1,13);
    else; Xnc[ll,.]=x[ncind,kvec]; endif;
if c<1; Xnc[11,13]=0;
    else; Xnc[ll,13]=((inc*0.5/2080*1/3*(2*d)/45)+0.26*(2*d*0.5)+Xnc[ll,13])/100; endif;

Xdif[ll,.]=Xnc[ll,.]-Xc[ll,.]; 
    if t le -98; tps[ll,.]=0;
    else; tps[ll,.]=t; endif; 

    if tps[ll] ge 50; tps[ll]=50;
    endif;
    endo;
endo; 
 
Xinitial=zeros(n,rows(kvec));

perch0 = .75;
trout0 = 19.4;
walleye0 = 7.4;
bass0 = 15.0;
FCA20 = 0;
FCA30 = 0;
FCA40 = 1;
FCA50 = 0;
FCA60 = 0;
FCA70 = 0;
FCA80 = 0;
FCA90 = 0;
fee0 = 3;
tc0=(((inc1*0.5/2080.*(1/3).*(2*dist)/45)+0.26*(2*dist*0.5)+fee0)/100);

Xinitial=(perch0~trout0~walleye0~bass0~FCA20~FCA30~FCA40
      ~FCA50~FCA60~FCA70~FCA80~FCA90).*ones(647,1)
      ~tc0;

_max_parnames="perch   "|"salmon  "|"walleye "|"bass    "|
        "FCA2 dum"|"FCA3 dum"|"FCA4 dum"|"FCA5 dum"|"FCA6 dum"|"FCA7 dum"|
        "FCA8 dum"|"FCA9 dum"|"tripcost"|"constant"|"gend    "|"age     "|
        "ret     "|"boat    "|"s0e     "|"s0g     "|"reg1     "|"reg2    ";

theta0=zeros(22,1);
load babogb; theta0[1:4]=babogb[1:4]*10; theta0[5:13]=babogb[5:13]; 
theta0[14:22]=zeros(9,1); 
load babgbnp; theta0=babgbnp; theta0[20]=5;
_max_active=ones(13,1)|ones(9,1);_max_active[14]=1; _max_active[19]=1; _max_active[20]=1;
_max_GradTol=1e-5;
_max_covpar=0;
"Tolerance = " _max_Gradtol;
__output=2;

__title="A-B-GB-NP-sigmae^2 = 1/2, V(sum)'s estimated, trip costs ";

{b,f,g,c,retcode} = maxprt(MAXLIK(Xdif,0,&nrprob,theta0));

proc nrprob(theta,X);

beta=theta[1:13];
soe=exp(theta[19]);
sog= exp(theta[20]);
sigmae=.5 ;

V0=(theta[14]+theta[15]*gend+theta[16]*age/100
   +theta[17]*ret+theta[18]*boat+theta[21]*reg1+theta[22]*reg2);
V0r=(theta[14]+theta[15]*z2[.,23]+theta[16]*z2[.,24]/100
   +theta[17]*z2[.,25]+theta[18]*z2[.,26]+theta[21]*reg1_2+theta[22]*reg2_2);

p = cdfn(-Xdif*beta/sqrt(2*sigmae));

e1 = (Xc*beta-V0)/sqrt(soe);

e2 = (-Xdif*beta/sqrt(2*sigmae));

rho = sigmae/sqrt(2*sigmae*(soe));

p0 = cdfbvn(e1, e2, rho)./p;

lnp = sumc( tps.*ln(p0)+(tot1 - tps).*ln(1-p0)+ln(p) );

pnr = cdfn((Xinitial*beta - V0r)/sqrt(sog));
lnr = sumc(ng.*ln(pnr) + (tot2-ng).*ln(1 - pnr));
lnL = lnp + lnr; 
retp(lnL);

endp;

save babgbnp=b; @ estimated parameters @
save cabgbnp=c; @ estimated covariance matrix @

"% correct predictions (p > .5): " sumc(p .gt .5)/(5176-138);

output off; end;