@ R2ABQ10.CMD  --  Random A-B probit model with two random
parameters: catch and FCA, and deterministic multiplicative 
differentials for each of the parameters for fish species and
FCA levels. Quadrature is used for approximation.
Sept. 29, 1999  @

@ Standard deviations of catch and FCA parameters vary in 
proportion to their means @

@ TEN evaluation points @

new;
library maxlik;
clear xdif;
output file = R2ABQ10.OUT reset;

"R2ABQ10.CMD. This is the random A-B probit model with two random variables, the
catch & FCA parameters. For different species and FCA levels, a deterministic 
multiplicative component is used, resulting in standard deviations that vary
in proportion to the means for the random parameters.  Quadrature is used
for approximation using TEN evaluation points. Sept. 29, 1999.";

/* Model Selection: choice of explanatory variables */

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

/* 1-4 are the catch rates for four species: perch, salmonids, walleye, bass;
5-12 are the FCA dummies for categories 2-9; 13 is the launch
fee */

/* Choose which parameter(s) to make random. rvec should be the same
length as kvec. All elements with a one will have the same random 
component. Then the random parameters will be scaled up or down 
by a deterministic multiplicative component */

let rvec1 = 1 1 1 1 0 0 0 0 0 0 0 0 0;
let rvec2 = 0 0 0 0 1 1 1 1 1 1 1 1 0;

/* Choose the number of quadrature points (2-10) */

npoints = 10;

/* Choose the starting values */
@ if strtflag eq 0; start with adjusted NRAB parameters @
@ if strtflag eq 1; start with nonrandom probit @

strtflag = 1;

/*  Z is an n x 9 matrix. Its rows are the individuals.
The first eight columns are the results of the paired comparisons (1 or 2).
If the respondent did not choose, then a 0 is entered.
The ninth row contains the version number (1-10).  Z is used to form
Xdif, or the differences in characteristics. */

loadm Z[647, 9] = Z.dat;
n = rows(Z);
Xdif =  zeros(8*n,rows(kvec));
y1 = zeros(10,9);
y2 = zeros(10,9);
w1 = zeros(10,9);
w2 = zeros(10,9);

/*  X is a 160 x 13 matrix of characteristics of the choices. Its rows are
arranged by version (1-10), then choice pair (1-8), then choice (A or B).
Its columns are the four fish catch rates; the FCA dummies (eight
dummies for categories 2-9); and the launch fee.
There are 16 lines for each version. The first two lines of the 16 are
pair 1 alternative A and pair 1 alternative B; the next two lines are 
for pair 2, the next two are for pair 3, and so on... */

loadm X[160, 13] = X.dat;
"
Estimation uses X.dat.";

/*  The following loop creates the matrix Xdif, dimension 8n x 13,
where n is the number of individuals, and there are 13 explanatory
variables in the model. */

ll=0; i=0;
do while i<n;   i=i+1;
  j=0; 
  do while j<8;         j=j+1;  ll=ll+1; 
    c=z[i,j];     @ choice made by i in pair j (1 or 2) @
    v=z[i,9];  @ survey version for individual i @
    cind=16*(v-1)+2*(j-1)+c;
    @ cind is the row of X that has the
      characteristics of the alternative i chose
      in pair j, version v @
    ncind=16*(v-1)+2*(j-1)+3-c;
    @ ncind is the row of X that has the
      characteristics of the alternative i did
      not choose in pair j, version v @
    if c<1; Xdif[ll,.]=zeros(1,13);
    @ for missing choice, Xdif set to 0.
      no effect on likelihood for these pairs @
    else;
    Xdif[ll,.]=X[ncind,kvec] -  X[cind, kvec];
    @ Xdif[ll,.] is the characteristics of the
      alternative not chosen minus the 
      characteristics of the alternative chosen
      for pair ll, ll=1,2,...,nx8. Xdif is arranged
      by i then j @
    endif;
    endo;
  endo;

load bab;
theta0=bab[1]|bab[2:4]/bab[1]|bab[5]|
	      bab[6:12]/bab[5]|bab[13];
@ Starting values are based on NRAB @

names =   "perch"|"salmon fact"|"walleye fact"|"bass fact"|
		 "FCA2 dum"|"FCA3 dum fact"|"FCA4 dum fact"|
		 "FCA5 dum fact"|"FCA6 dum fact"|"FCA7 dum fact"|
		 "FCA8 dum fact"|"FCA9 dum fact"|"launch fee"|
		 "sigmav1"|"sigmav2";
_max_parnames = names[kvec]; @ subset of active parameters @
nparam=rows(_max_parnames);

/* call MAXLIK to use nonrandom probit as starting values */

format 10,4;

rv1=names[1:4];
"First random coefficient set is " $rv1;

rv2=names[5:12];
"Second random coefficient set is " $rv2;

"Number of quadrature points = " npoints;

@ The next section is the nonrandom AB probit model. It only runs if
strtflag eq 1 @

if strtflag eq 1;
__title = "No random parameters";
{theta,f,g,cov,retcode} = maxprt(MAXLIK(xdif,0,&nrprob,theta0));
endif;
proc nrprob(theta,x); @ start of proc for nonrandom probit @
local beta,p;
beta=(theta[1]*ones(4,1)|theta[5]*ones(8,1)|theta[13]).*
	 (1|theta[2:4]|1|theta[6:12]|1);
p=cdfn(-xdif*beta);
retp(ln(p));
endp;

/* Random probit -- two random parameters
Set quadrature points and weights for npoints 2-10 */

let y1[10,9]=
  0.707106781186548  0.0                               0.524647623275290  
  0.0                              0.436077411927617   0.0                           
  0.381186990207322  0.0                               0.342901327223705
 
 -0.707106781186548  1.224744871391589  -0.524647623275290
  0.958572464613819 -0.436077411927617   0.816287882858965  
 -0.381186990207322  0.723551018752838  -0.342901327223705
 
  0.0                              -1.224744871391589  1.650680123885785
 -0.958572464613819   1.335849074013697 -0.816287882858965  
  1.157193712446780  -0.723551018752838  1.036610829789514

  0.0                               0.0                             -1.650680123885785  
  2.020182870456086  -1.335849074013697  1.673551628767471
 -1.157193712446780   1.468553289216668 -1.036610829789514

  0.0                               0.0                              0.0  
 -2.020182870456086   2.350604973674492 -1.673551628767471  
  1.981656756695843  -1.468553289216668  1.756683649299882

  0.0                               0.0                              0.0  
  0.0                              -2.350604973674492  2.651961356835233
 -1.981656756695843   2.266580584531843 -1.756683649299882

  0.0                               0.0                              0.0  
  0.0                               0.0                             -2.651961356835233
  2.930637420257244  -2.266580584531843  2.532731674232790
  
  0.0                               0.0                              0.0  
  0.0                               0.0                              0.0  
 -2.930637420257244   3.190993201781528  -2.532731674232790

  0.0                               0.0                              0.0
  0.0                               0.0                              0.0
  0.0                              -3.190993201781528  3.436159118837738

  0.0                               0.0                              0.0
  0.0                               0.0                              0.0
  0.0                               0.0                             -3.436159118837738; 

@ e.g., y[1,1] =0.707..., y[1,2]=0, y[1,3]=0.524..., y[2,1]=-0.707...,
the 9 columns are the possible number of quadrature points (2-10),
and the rows are the corresponding values of y. For example, y[3,4]
= -0.9587... is the third y if there are 5 quadrature points, and y[3,1]=
y[4,1]=y[5,1]=0 because the first column corresponds to the case of
only 2 quadrature points. @

let y2[10,9]=
  0.707106781186548  0.0                               0.524647623275290  
  0.0                              0.436077411927617   0.0                           
  0.381186990207322  0.0                               0.342901327223705
 
 -0.707106781186548  1.224744871391589  -0.524647623275290
  0.958572464613819 -0.436077411927617   0.816287882858965  
 -0.381186990207322  0.723551018752838  -0.342901327223705
 
  0.0                              -1.224744871391589  1.650680123885785
 -0.958572464613819   1.335849074013697 -0.816287882858965  
  1.157193712446780  -0.723551018752838  1.036610829789514

  0.0                               0.0                             -1.650680123885785  
  2.020182870456086  -1.335849074013697  1.673551628767471
 -1.157193712446780   1.468553289216668 -1.036610829789514

  0.0                               0.0                              0.0  
 -2.020182870456086   2.350604973674492 -1.673551628767471  
  1.981656756695843  -1.468553289216668  1.756683649299882

  0.0                               0.0                              0.0  
  0.0                              -2.350604973674492  2.651961356835233
 -1.981656756695843   2.266580584531843 -1.756683649299882

  0.0                               0.0                              0.0  
  0.0                               0.0                             -2.651961356835233
  2.930637420257244  -2.266580584531843  2.532731674232790
  
  0.0                               0.0                              0.0  
  0.0                               0.0                              0.0  
 -2.930637420257244   3.190993201781528  -2.532731674232790

  0.0                               0.0                              0.0
  0.0                               0.0                              0.0
  0.0                              -3.190993201781528  3.436159118837738

  0.0                               0.0                              0.0
  0.0                               0.0                              0.0
  0.0                               0.0                             -3.436159118837738; 

let w1[10,9]=

0.8862269254528        1.18635900604               0.8049140900055
0.9453087204829        0.7246295952244           0.8102646175568
0.6611470125582        0.7202352156061           0.6108626337353

0.8862269254528        0.2954089751509           0.8049140900055
0.3936193231522        0.7246295952244           0.4256072526101
0.6611470125582        0.4326515590026           0.6108626337353

0.0                                0.2954089751509           0.08131283544725  
0.3936193231522        0.1570673203229           0.4256072526101
0.2078023258149        0.4326515590026           0.2401386110823

0.0                                0.0                                   0.08131283544725  
0.01995324205905      0.1570673203229           0.05451558281913
0.2078023258149        0.08847452739438         0.2401386110823

0.0                                0.0                                   0.0  
0.01995324205905      0.004530009905509       0.05451558281913
0.01707798300741      0.08847452739438         0.03387439445548

0.0                                0.0                                   0.0
0.0                                0.004530009905509       0.0009717812450995
0.01707798300741      0.004943624275537       0.03387439445548

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0009717812450995
0.0001996040722114  0.004943624275537       0.001343645746781

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0
0.0001996040722114  0.00003960697726326   0.001343645746781

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0
0.0                                0.00003960697726326   0.000007640432855233

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0
0.0                                0.0                                   0.000007640432855233;

@ Note that the sum of the weights is always the square root of pi.
e.g., .88622+.88622 = (pi)^.5 @

let w2[10,9]=

0.8862269254528        1.18635900604               0.8049140900055
0.9453087204829        0.7246295952244           0.8102646175568
0.6611470125582        0.7202352156061           0.6108626337353

0.8862269254528        0.2954089751509           0.8049140900055
0.3936193231522        0.7246295952244           0.4256072526101
0.6611470125582        0.4326515590026           0.6108626337353

0.0                                0.2954089751509           0.08131283544725  
0.3936193231522        0.1570673203229           0.4256072526101
0.2078023258149        0.4326515590026           0.2401386110823

0.0                                0.0                                   0.08131283544725  
0.01995324205905      0.1570673203229           0.05451558281913
0.2078023258149        0.08847452739438         0.2401386110823

0.0                                0.0                                   0.0  
0.01995324205905      0.004530009905509       0.05451558281913
0.01707798300741      0.08847452739438         0.03387439445548

0.0                                0.0                                   0.0
0.0                                0.004530009905509       0.0009717812450995
0.01707798300741      0.004943624275537       0.03387439445548

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0009717812450995
0.0001996040722114  0.004943624275537       0.001343645746781

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0
0.0001996040722114  0.00003960697726326   0.001343645746781

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0
0.0                                0.00003960697726326   0.000007640432855233

0.0                                0.0                                   0.0
0.0                                0.0                                   0.0
0.0                                0.0                                   0.000007640432855233;

_max_parnames=names;
nparam=rows(_max_parnames);

if strtflag eq 0; theta0=bab[1]|bab[2:4]/bab[1]|
			bab[5]|bab[6:12]/bab[5]|
			bab[13]|-0.8941|-1.2031; endif; 
@ Start with NRAB if strtflag = 0 @
if strtflag eq 1; theta0=theta|-0.8941|-1.2031; endif; 
@ Start with nonrandom probit, and values for the standard
deviations based on an earlier unweighted model estimated
with an incomplete sample if strtflag = 1 @

"Starting values are: " theta0';

__title="R2ABQ10.cmd: A-B probit 
with two random parameters (catch, FCA), 
uses TEN evaluation points";
{theta,ff,g,cov,retcode} = maxprt(MAXLIK(Xdif,0,&r2prob,theta0));
"Sigmav1 = " exp(theta[nparam-1]);
"Sigmav2 = " exp(theta[nparam]);

proc r2prob(theta,X);  @ Start of proc for two random parameters @
local beta, bscale, p, sigv1, sigv2, nquad1, nquad2, i, j, 
sum1, sum2, f;
beta=(theta[1]*ones(4,1)|theta[5]*ones(8,1)|theta[13]);
bscale=(1|theta[2:4]|1|theta[6:12]|1);
sigv1=exp(theta[nparam-1]);
sigv2=exp(theta[nparam]);

@ Quadrature loop: npoints is the number of evaluation points. @
i=0;
p=0;
do while i<n;   @ looping over individuals @ 
	i=i+1;
	@ outer loop @
	nquad1=0;
	sum1=0;
	do while nquad1 < npoints;
	@ Hermite quadrature excluding pi or f
	  for i where p is the joint probability of 
	  AB choices over the pairs @
		nquad1=nquad1+1;
		@ inner loop @
		nquad2=0;
		sum2=0;
		do while nquad2<npoints;
			nquad2=nquad2+1;
			f=1;
			j=0;
			do while j<8;
			@ Looping over the pairs for i @
				j=j+1;
				f = f*cdfn(-xdif[8*(i-1)+j,.]*(beta.*bscale)
				- 1.414214*xdif[8*(i-1)+j,.]*(rvec1.*bscale)
				*sigv1*y1[nquad1,npoints-1]
				- 1.414214*xdif[8*(i-1)+j,.]*(rvec2.*bscale)
				*sigv2*y2[nquad2,npoints-1]);
				@ product over J @
			endo;
			sum2=sum2 + w2[nquad2,npoints-1]*f;
			@ sum over pts (2) @
		endo;
		sum1 = sum1 + w1[nquad1,npoints-1]*sum2;
		@ sum over pts (1) @
	endo;
	p = p + ln(sum1/pi);   @ pi, not sqrt(pi) @
	endo;
retp(p);
endp;
end;