new;

library maxlik;

/*=========================================================================================================*/
/*========================== Loading data in GAUSS: note our ASCII file can be loaded in the same way as this txt file.  =======*/
/*=========================================================================================================*/

/* 10/3 MB These are the paths for nedatap & rmdataps now */
load rdata4[1680,8] = d:\schotter\424cralp\expt3.txt;     
load rdata6[1680,8] = d:\schotter\424cralp\expt1.txt;


/*=========================================================================================================*/
/*==========================Creating some variables--you don't need to do this way ==================================*/
/*=========================================================================================================*/
/* 10/3mb  There are no dummies for elicitation in this likelihood--I don't know where a test was done for the beta0/beta1 rest'n. */
/* Create a data set with both elicitation and no elicitation data & with a dummy for no elicitation.   Use the above models as bases 
for the restricted model, and then write the unrestricted with dummy entries in the b1-4. */

newmdda1 = rdata4 ~ ones(1680,1);
newmdda2 = rdata6[.,1:8] ~ zeros(1680,1);
newmddat = newmdda1 | newmdda2;

/*======================================================================================================*/
/*======= Procedure llhood sets up the likelihood of each observation in a vector.  ====================================*/
/*========= b is the vector of parameters to be estimated, xx is the matrix of data =================================*/
/*========= Note that b should contain all the betas AND the variance coefficient sigma =========================*/
/*=========================================================================================================*/

proc llhood(b,xx);
local beliefs, ii, action_t, action_s, weights, paydiff, pr_red, pr_green, dummy;
local red, green, pr, help, llval;

beliefs = zeros(rows(xx),1);
beliefs[1] = .5;
beliefs[2] = xx[1,6];
ii = 2;
do until ii > rows(xx)-1;
    action_t = xx[ii,6];
    action_s = xx[1:ii-1,6];
    weights = b[5] ^ rev(seqa(1,1,ii-1));
    beliefs[ii+1] = (action_t + sumc(action_s .* weights))/(1 + sumc(weights));
    ii = ii + 1;
endo;

paydiff = xx[.,1] .* (5 * beliefs - 3);
dummy = xx[.,9];

pr_red = exp(b[1] + b[2] * paydiff + b[3] * dummy + b[4] * dummy .* paydiff) ./ (1 + exp(b[1] + b[2] * paydiff + b[3] * dummy + b[4] * dummy .* paydiff));
pr_green = 1 - pr_red;

red = (xx[.,5] .== 1);
green = (xx[.,5] .== 0);

pr = (pr_red .^ red) .* (pr_green .^ green);

/*=======================================================================================================================*/
/*========== retp(pr) outputs the calculated vector of probabilities of the observed outcomes, based on data xx and given parm vector b  =======*/
/*====================================================================================================================*/
retp(pr);
endp;


/*=========================================================================================================*/
/*=========== Procedure llfun takes the sum of the log of the vector of probabilities calculated in llhood. =====================*/
/*======== If you like you can combine the steps creating the probability vector and summing its log into one llfun procedure =======*/
/*=========================================================================================================*/
proc llfun(b,xx);
local pr, pp, r1, r2, data, help, llval;

pr = zeros(rows(xx),1);
pp = 1;
do until pp > 56;
    r1 = 1 + 60 * (pp - 1);
    r2 = 60 + 60 * (pp -1);
    data = xx[r1:r2,.];
    pr[r1:r2] = llhood(b,data);				/* 10/3mb same as loop in above proc but for 54 subject obs's*/
    pp = pp + 1;
endo;

help = pr .* (pr .> .00001) + .00001 * (pr .<= .00001);		
/*=========================================================================================================*/
/*========================== Summing the log probabilities =========================================== =======*/
/*=========================================================================================================*/
llval = sumc(ln(help));

/*=========================================================================================================*/
/*==========================Returning the sum of the log probabilities, which is our log likelihood. =======================*/
/*=========================================================================================================*/
retp(llval);
endp;


/*=========================================================================================================*/
/*========================== This line sets up a vector of initial guesses at the b parameter vector ====================*/
/*=========================================================================================================*/
b5 = 0.5 * ones(5,1);  b3 = 0.5 * ones(3,1);

/*=========================================================================================================*/
/*========================== This line uses the maxlik library to run the likelihood.  ==================================*/
/*========================= It will iterate until it converges to parameter estimates.  ===================================*/
/*=========================================================================================================*/
{bu,fu,gu,hu,ju} = maxprt(maxlik(newmddat,0,&llfun,b5));

end;