% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 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(:,8)==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_Mag.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_Mag betamat gmat ...
Emat cmat gimat;