@ LCJ_ 4.CMD 05/28/2010 @
@ This is the Gauss .cmd file to estimate the lcj model with 8 GB SP choice pairs and the likert-scale responses to the 15 preference statements. Assuming 4 classes, and 4 angler types
to be used as covariates. @
@The code is easily modified to decrease, or increase, the assumed number of classes @
@two data sets are read in, zs.fmt (a gauss matrix) and x.dat, a text-delimited .tex file @
@ zs is the individual-specifc data, and .fmt file, 640 rows x 121 columns - there are 640 individuals in the data set.@
@This matrix is already sorted: 117 women, 191 non-retired men income <$50K, 286 non-retired men income > $50K, 46 retired men @
@Columns 30-121 indentify each individual's answer to each preference statement. The questions in order are Q10a,b,c,Q13a,b,Q135a,b,c,d,e,f,g,h,i,and Q1
All of these questions have five levels, plus a no response option, except for Q1, which has seven levels plus no response @
@ So, for example, column 32 is a 1 if the individual choose level 2 to question Q10a @
@ column 27 is the survey version the individual saw @
@ z[i,j], j=1-8 is a 1 if individual i choose alternative j in SP pair j, and a 2 if he choose alternative j @
@ describe what is in the other columns @
@ x.dat is a text file (space delimited), not a guass data set. It dimensions are 160 x 13 @
@ It specifies the attribute levels for each GB alternative in each pair in each version (1-10 versions) of the survey @
@ first 8 rows correspond are version 1, alternative A, next 8 lines version 1 alternative B, ect. etc. @
@ the first 8 columns are the eight dummies for FCA levels (one 1 and seven zeros), the next four are the catch times by species, and the last column is the price variable @
@-------------------------------------------------@
@ First we create a Xdif vector, dimension (640 x 8) x 13, where 8 is the # of choice pairs and 13 is the number of attributes @
@ each of the 13 columns corresponds to the attribute level of the alternative the individual did not choose minus its level in the chosen alternative @
@ For example, row 3 corresponds to individual 1, pair 3, and row 13 corresponds to individual 2, pair 5, where row 3 column 13 is, for indiv 1, the fee for the alternative
not chosen minus the fee for the alternative chosen @
new;
loadm Zs;
n = rows(Zs); @ n=640, the number of individuals in the sample @
loadm X[160, 13] = X.dat; @ loading attribute data @
@ create xdiff and initially fill it with zeros @
@xc and xnc are matrices used to creat xdiff. "c" refers to "chosen" and "nc" refers to "not chosen" @
clear xdif, xc, xnc;
Xdif=zeros(8*n,13);
Xc=zeros(8*n,13);
Xnc=zeros(8*n,13);
@ the matix xc identifies the attributes of the alternative the individual chose @
@ the maxrix xnc identifies the attributes of the alternative the individual did not chose @
ll=0; i=0; @ ll is a scalar indicator for the individual @
do while i> @
j=0;
do while j<8; j=j+1; ll=ll+1; @ << iterate over each choice pair >> @
c=zs[i,j]; @ choice made by i in pair j (1 or 2); thatis z[i,j] has a 2 in it if individual i chose alternative 2 i pair j @
v=zs[i,27]; @ survey version for individual i@
cind=16*(v-1)+2*(j-1)+c; @cind=choice index, it identifies the row number in x.dat that has the attribute levels for the alternative indiv i choose in pair j @
@ for example if it is pair 5 in version 3, and the individual chose alternative 2, it is 16(3-1)+2(5-1)+2= 42@
ncind=16*(v-1)+2*(j-1)+3-c;
@ncind=choice index, it identifies the row number in x.dat that has the attribute levels for the alternative indiv i did not choose in pair j @
@for example if it is pair 5 in version 3, and the individual chose alternative 2, it is 16(3-1)+2(5-1)+1=41 @
if c<1; Xc[ll,.]=zeros(1,13); @ if didn't make a choice, set all attribute levels = 0 @
else; Xc[ll,.]=x[cind,1 2 3 4 5 6 7 8 9 10 11 12 13]; endif; @ if the individual answered pair j, it plugs into row ll the 13 attribute levels of the alternative he choose @
if c<1; Xc[ll,13]=0; @erm: resets the price level to zero if he did not choose an alternative in this pair @
else; Xc[ll,13]=Xc[ll,13]/10; endif; @ rescaling the price by dividing it by 10 @
@ so at this point have filled in one row of Xc, so, for one individual for one pair, the next set of lines does the same of Xnc @
if c<1; Xnc[ll,.]=zeros(1,13);
else; Xnc[ll,.]=x[ncind,1 2 3 4 5 6 7 8 9 10 11 12 13]; endif;
if c<1; Xnc[ll,13]=0;
else; Xnc[ll,13]=Xnc[ll,13]/10; endif;
Xdif[ll,.]=Xnc[ll,.]-Xc[ll,.]; @ this files in one row of Xdiff, for individual i, pair j @
endo; @ end of j loop @
endo; @ end of the i loop @
@ so at this point the Xdif matrix has been created. It will be called everytime we calculate the prob that the individual chose the SP alternatives that he chose @
@------------------------------------------------------------------------------------------ @
library maxlik; @loading the packages for maxlik, so maxlik can be use below @
maxset;
@ everything below this point is in an E-M iteration loop@
@ each iteration starts with working estimates for the conditional class membership probabilities,
and uses them, with the data to come up with better estimates of the utility parameters .......@
iter=1; do while iter<199; @ 200 iterations sufficient for convergence @
clearg m, j, beta1, beta2, beta3, beta4, p, p1a, p2a, p3a, p4a, p1, p2, p3, p4, prob, lik,
att, pi1, pi2, pi3, pi4, pis, probatt1, probatt2, probatt3, probatt4;
output file = LCJ_4.out reset;
@ STEP 1 in E-M iteration: load the conditional class-membership probabilities: from last iteration, (guesses if first iteration) @
@ E.g. cond14 is dimension 640 x 1 and a .fmt file. The 4 refers to the 4-class model, the 1 refers to class 1@
@ condit14 is a vector of the conditional probabilities of being a member of class 1 @
@ we are reading into the iteration current best estimates for the 4 unconditional class-membership probabilities @
@To start, assume that conditional (and therefore unconditional) probabilities are the same - 0.25 - for everyone @
cond14=0.25*ones(640,1);
cond24=cond14;
cond34=cond14;
cond44=1-cond14-cond24-cond34;
uncond14=cond14;
uncond24=cond14;
uncond34=cond14;
uncond44=1-uncond14-uncond24-uncond34;
@ Step 2 in the E-M iteration @
@ Based on these working estimates of the conditional and unconditional class-membership probabilities, we calculate best estimates of the (pi)qs:c,
where (pi)qs:c is the probability that that an individual in class c will choose level s in preference statement q.@
@These are calculated using Eq 6 @
att=Zs[1:640,30:121];
@ This is 640 x 92matrix. "att" refers to attitudes @
@ the "att" are the xiqs? ones and zeros wrt whether individual i (1-640) choose level s to question q? For each q there is one 1 and the rest are zeros.See the explan of z above @
pi1=sumc(cond14.*att)/(meanc(uncond14)*640); @ This is 92 x 1: for class 1 there are 92 (Pi)qs for class1 @
pi2=sumc(cond24.*att)/(meanc(uncond24)*640); @ note that the denominator is a scalar, the estimated number of people in class 2 @
pi3=sumc(cond34.*att)/(meanc(uncond34)*640);
pi4=sumc(cond44.*att)/(meanc(uncond44)*640);
@ these are the current best estimates of the (pi)qs:c @
@the (pi)c are so simple because there are no covariates, they are simply the proportion of times those in class c choose level s to question q @
@------------------------------------------------------------------------- @
@ Step 3 in the E-M algorithm @
@create here the joint probability of obser the individual's responses to the attitudinal questions, conditional on the (pi)c just calculated @
@ this will be used in proc that defines the conditional lik function @
Probatt1=prodc(((pi1').^att)'); @Eq. 3, 640 x1 vector, conditional on being in class 1 @
Probatt2=prodc(((pi2').^att)'); @ conditional on being in class 2 @
Probatt3=prodc(((pi3').^att)');
Probatt4=prodc(((pi4').^att)');
@----------------------------------------------------------------------------@
@ Step 4 in the E-M iteration: now we have to estimate the utility parameters, the beta, conditional on all of the above estimates @
@ Step 4 does not estimate the utility parameters but rather creates the lnlik proc that will be called by maxlik @
@ Note that the probability of choosing an alternative, conditional on class, is specified a a probit (one could change it to a logit probability) @
@ create a proc that contains the ln lik function, Eq. 1, conditional on the above estimates @
@ It caclulates, for each individual, the ln of the joint probability of their choices and answers.
That is, it calculates thei contribution to the ln lik function, conditional on the current estimates of the class-membership probabilities @
proc abprob(b,X);
Prob=zeros(n,1); @640 rows and 1 column, row i will become, below the joint probability of observing indiv i's SP choices and attitudinal answers, conditional on the above estimates @
p1=zeros(n,1); @640 rows and 1 column, row i will become, below the joint probability of observing only indiv i's SP choices, conditonal on them being in class 1 @
p2=zeros(n,1);
p3=zeros(n,1);
p4=zeros(n,1);
@ The objective is to create the prob vector, which is a function the data,the beta, and the most current estimates of the conditional class-membership probabilities and (pi)qs:c @
@ the vector is filled in row by row using the the following loop over individuals and subloop over pairs @
@ start angler loop @
@ << These betas represent the betas on the attributes, by class. >> @
m=1; do while m<=n;
p1a=1; @ p1a is a scalar, initially specified as a 1, it will become the joint probability of the SP choices for individual m, conditional on being in class 1
(as a function of the beta for class 1 @
@ after it is calculated as a function of the betas it will become a row in p1 @
p2a=1;
p3a=1;
p4a=1;
j=0; do while j<8; j=j+1; @ for the given m we are now looping over the j pairs. @
@before this program is first written, one must create b.fmt, a 52 x 1 vector, of starting values for the beta @
beta1=b[1:13]; @ beta1 is class 1's betas, read from the first 13 row of the vector b: the betas assuming class 1 @
beta2=b[14:26];
beta3=b[27:39];
beta4=b[40:52];
@ these next lines are a bit tricky. Note that the first is p1a multiplied by p1a multipled by a term that is a function of j. It is building up a product across the 8 j @
p1a = p1a*cdfn(-Xdif[(m-1)*8+j,.]*beta1); @ when it finishes the j loop, this will be P(y|c=1): the prob that individual m made the SP choices he made conditional
on the betas for class 1 @
p2a = p2a*cdfn(-Xdif[(m-1)*8+j,.]*beta2); @ << !! Choice model assumes probit. Needs to be specified in paper !! >> @
p3a = p3a*cdfn(-Xdif[(m-1)*8+j,.]*beta3);
p4a = p4a*cdfn(-Xdif[(m-1)*8+j,.]*beta4);
endo; @ the end of the j loop @
@note again how Pca is the the joint probability of only the SP choices for indiv m, conditional on class but this only after the j loop is completed @
@ the next line fills in the indiv m row in each of the pc vectors @
p1[m]=p1a; p2[m]=p2a; p3[m]=p3a; p4[m]=p4a; m=m+1; @for the individual put is in the correct p1 row. @
endo; @ this is the end of the m indiv loop @
@ we now have the pc vectors as a function of the betas, conditional on class @
@ the next line is individual i joint contribution to the like function, equation 2 @
Prob=uncond14.*p1.*probatt1 + uncond24.*p2.*probatt2 + uncond34.*p3.*probatt3 + uncond44.*p4.*probatt4;
@ row i is the prob of indiv i's choices and answers as a function of the beta parameters and conditonal on the current uncond class-membership probabilities, 640x1 vector, and
the current best estim of the joint probability of the individual's answer to the attitudinal questions @
@ more specifically: it consists of 4 terms, one for each class, each term is the unconditional class membershp prob, multipled by the joint probability
of their SP choices (conditional on being in that class multiplied by the joint prob of observ the attitudina responses (conditional on being in that class) @
retp(ln(Prob)); @ this lns each Prob and returns their sum, which is the ln likeilhood function, conditonal on .... @
endp; @ << Ending procedure >> @
@-------------------------------------------------------------------------------------@
@Step 5 of an E-M iteration @
@ now maklik calls the proc to estimate the the beta, the utility parameters @
_max_parnames="perch G1"|"salmonG1"|"walleyG1"|"bass G1"|
"FCA2dum1"|"FCA3dum1"|"FCA4dum1"|"FCA5dum1"|"FCA6dum1"|"FCA7dum1"|
"FCA8dum1"|"FCA9dum1"|"fee G1"|
"perch G2"|"salmonG2"|"walleyG2"|"bass G2"|
"FCA2dum2"|"FCA3dum2"|"FCA4dum2"|"FCA5dum2"|"FCA6dum2"|"FCA7dum2"|
"FCA8dum2"|"FCA9dum2"|"fee G2"|
"perch G3"|"salmonG3"|"walleyG3"|"bass G3"|
"FCA2dum3"|"FCA3dum3"|"FCA4dum3"|"FCA5dum3"|"FCA6dum3"|"FCA7dum3"|
"FCA8dum3"|"FCA9dum3"|"fee G3"|
"perch G4"|"salmonG4"|"walleyG4"|"bass G4"|
"FCA2dum4"|"FCA3dum4"|"FCA4dum4"|"FCA5dum4"|"FCA6dum4"|"FCA7dum4"|
"FCA8dum4"|"FCA9dum4"|"fee G4";
nparam=rows(_max_parnames);
theta0=0.0*ones(52,1); @ starting values of the preference parameters @
__output=2;
__title="LCJ 4-Class with covariates";
{b,f,g,c,retcode} = maxprt(MAXLIK(Xdif,0,&abprob,theta0)); @ << Estimate beta based on the conditional lik function, conditional on the unconditional class-membership prob @
save beta4=b; @ max lik outputs new estimates of the preference parameters @
@---------------------------------------------------------------------------------------- @
@ STEP 6 of E-M iteration: given we now have new beta estimates these can used to update/compute new estimates of the conditional class-membership probabilities @
@ Specifically, use Bayes theorem -Equation 7 - to update the estimates of the conditional class-membership probabilities
which depend on the current best estimates of the (pi),the the old unconditional class-membership probabilities, the joint prob of the attitud responses, and
the new joint prob of observing the sp choices (this last one the only thing that depends on the proc and maxlik. @
@ Note that condc4 is simply rearranging the stuff in Prob @
cond14=uncond14.*Probatt1.*p1./prob; @ if I understand Gauss? p1 and prob should be the last ones produced in the proc by running max like @
cond24=uncond24.*Probatt2.*p2./prob;
cond34=uncond34.*Probatt3.*p3./prob;
cond44=1-cond14-cond24-cond34;
@ << Save conditional and unconditional class membership probabilities for next iteration >> @
save cond14;
save cond24;
save cond34;
uncond14=meanc(cond14[1:117,1])*ones(117,1)|meanc(cond14[118:308])*ones(191,1)|meanc(cond14[309:594])*ones(286,1)|meanc(cond14[595:640])*ones(46,1);
uncond24=meanc(cond24[1:117,1])*ones(117,1)|meanc(cond24[118:308])*ones(191,1)|meanc(cond24[309:594])*ones(286,1)|meanc(cond24[595:640])*ones(46,1);
uncond34=meanc(cond34[1:117,1])*ones(117,1)|meanc(cond34[118:308])*ones(191,1)|meanc(cond34[309:594])*ones(286,1)|meanc(cond34[595:640])*ones(46,1);
uncond44=1-uncond14-uncond24-uncond34;
save uncond14; save uncond24; save uncond34; save uncond44;
@-----------------------------------------------------------@
iter=iter+1; endo;
output off; end;