% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Hierarchical Ordered Probit 
% % re-parameterized a la Nandram & Chen
% % using common pool data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 

tic;  % start stop watch

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

load c:\klaus\aaec6564\mlab\worksp\CommonPoolData;

% variable list
%1  age            years
%2  edu            years of formal education
%3  gender         1 = female
%4  period         round of game (1-20 for each individual)
%5  effort         1-9 for each individual and round
%6  PAC            "1" if region = Pacific
%7  MAG            "1" if region = Magdalena
%8  CAR            "1" if region = Carribean
%9  sub_id         running panel ID for individuals (1-180)
%10 treat1         "1" if treatment = open access
%11 treat2         "1" if treatment = audit with low penalty
%12 treat3         "1" if treatment = audit with high penalty
%13 treat4         "1" if treatment = communication
%14 group          running panel ID for groups (1-36)
%15 runid          running ID (1-3600)
%16 expect         expected sum of others' payoff

%switch CAR and MAG to have a consistent sequencing of "1s"
MAG=data(:,7);
CAR=data(:,8);
data(:,7)=CAR;
data(:,8)=MAG;

N=180; %total number of original individuals
P=20; %total number of periods

f=find(data(:,6)==1);
data=data(f,:); 
y=data(:,5);
% Identify groupings in y by effort level:
f1=find(y==1);
f2=find(y==2);
f3=find(y==3);
f4=find(y==4);
f5=find(y==5);
f6=find(y==6);
f7=find(y==7);
f8=find(y==8);
f9=find(y==9);

l1=length(f1);
l2=length(f2);
l3=length(f3);
l4=length(f4);
l5=length(f5);
l6=length(f6);
l7=length(f7);
l8=length(f8);
l9=length(f9);


X=[data(:,3) data(:,2)/10 data(:,10:13)];
% contents of X:
%1  gender         1 = female
%2  age            years
%3  expect         expected sum of contributions by others
%4  treat1         "1" if treatment = open access
%5  treat2         "1" if treatment = audit with low penalty
%6  treat3
%7  treat4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end data transformations for T-combos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k=size(X,2);
kf=2; %fixed coefficients
kr=4;  %RCs

P=20; %total number of periods per individual
G=12; %number of groups
Gn=5; %individuals per group 
N=G*Gn; %total number of individuals
obs=N*P; %total sample size

% group y and X into individual-specific cells
% this is easier here since we have a perfectly balanced setup
Ximat=mat2cell(X,P*ones(N,1),k); 
yimat=mat2cell(y,P*ones(N,1),1);

% Now group y and X into group-specific cells
Xgmat=mat2cell(X,(P*Gn)*ones(G,1),k);
ygmat=mat2cell(y,(P*Gn)*ones(G,1),1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=3000; % burn-ins
r2=12000; % keepers
R=r1+r2;

nc=7; %number of cutoff points to be estimated (including the
% highest, which will be re-parameterized)

% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(obs-k);

% elements for beta (fixed coefficients)
mub0=zeros(kf,1); %diffuse prior for mean of fixed coefficients
Vb0=eye(kf)*100; % diffuse prior for varcov of fixed coefficients
% since we're drawing beta first, won't need starting draw

% elements for gamma
mug0=zeros(kr,1); 
Vg0=eye(kr)*100; 
gdraw=bols(kf+1:end);

% elements for diagonals of E
% variance; recall: covariances make no sense here...
phi0=1/2;
gam0=1/2;
Edraw=eye(kr);

% elements for del2 (squared inverse of highest cutoff (c7), but also new
% variance of regression error)
eta0=1/2;
kappa0=1/2; % diffuse prior shape and scale
del2draw=s2; % use OLS variance as starting draw for GS (Tobias uses "1" in his code)

% elements for cutoffs (I have 7 unknown cutoffs, but only 6 will be drawn explicitly
% since the highest cutoff is re-parameterized into del2)
c=0.8; %aim for 44% acceptance (single we're drawing from a single dimension in each case
% see Gelman p. 306
stdc1=c*sqrt(1/obs);
stdc2=c*sqrt(1/obs);
stdc3=c*sqrt(1/obs);
stdc4=c*sqrt(1/obs);
stdc5=c*sqrt(1/obs);
stdc6=c*sqrt(1/obs);

Vc=10;

% set starting draws such that cutoff ordering is preserved
% but note: highest re-paramterized cutoff is now 1!
c1draw=1/7;
c2draw=2/7;
c3draw=3/7;
c4draw=4/7;
c5draw=5/7;
c6draw=6/7; % that way we'll stay below 1

[betamat,gmat,Emat,cmat,gimat,acceptmat]...
    =gs_HOP_normal(y,Ximat,yimat,Xgmat,ygmat,k,kf,kr,P,G,Gn,N,obs,nc,r1,r2,mub0,Vb0,...
    mug0,Vg0,gdraw,phi0,gam0,Edraw,eta0,kappa0,del2draw,...
    stdc1,stdc2,stdc3,stdc4,stdc5,stdc6,Vc,f1,f2,f3,f4,f5,f6,f7,f8,f9,...
    l1,l2,l3,l4,l5,l6,l7,l8,l9,c1draw,c2draw,c3draw,c4draw,c5draw,c6draw);

'GS done'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Capture key results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open log file 
[fid]=fopen('c:\klaus\aaec6564\mlab\logs\mod11_HOP_Pac.txt','w');
if fid==-1;
   warning('File could not be opened');
else;
   disp('File opened successfully');
end;

fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'\n');
fprintf(fid,'smallest acceptance rate cor cmat =\t%6.3f \n',min(acceptmat));
fprintf(fid,'highest acceptance rate cor cmat =\t%6.3f \n',max(acceptmat));
fprintf(fid,'mean acceptance rate cor cmat =\t%6.3f \n',mean(acceptmat));

% put all draws together & get diagnostics
allmat=[betamat;gmat;Emat;cmat];%
kdiag=klausdiagnostics_greater0(allmat);

% beta stuff
out=kdiag(1:kf,:)';
fprintf(fid,'Output table for betas \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

%gamma stuff
out=kdiag(kf+1:k,:)';
fprintf(fid,'Output table for gamma \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

%E stuff
out=kdiag(k+1:k+kr,:)';
fprintf(fid,'Output table for RE variances \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

% cstuff
out=kdiag(end-nc+1:end,:)';
fprintf(fid,'Output table for cutoffs \n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);

st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;


save c:\klaus\aaec6564\mlab\worksp\mod11_HOP_Pac betamat gmat ...
    Emat cmat gimat;



