%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Post-estimation script for common pool models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;  % start stop watch

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

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

r2=100; % number of secondary draws

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PACIFIC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% prepare data
%%%%%%%%%%%%%%%%%%
load c:\klaus\aaec6564\mlab\worksp\CommonPoolData;
%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,:); 
m1=mean(data(:,3)) %share of female
m2=mean(data(:,2)/10);

f1=find(data(:,4)<6); %period indicator
f2=find((data(:,4)>10) & (data(:,4)<16));
f=sort([f1;f2]);
data(f,:)=[];  %work with period 2 stuff only

clear data;

% load model-specific posterior draws
load c:\klaus\aaec6564\mlab\worksp\mod11_HOP_Mag betamat gmat Emat cmat;

R=size(betamat,2);
kr=size(gmat,1);
kf=size(betamat,1);
%keep every 10th draw
filt=(1:10:R);
betamat=betamat(:,filt);
gmat=gmat(:,filt);
Emat=Emat(:,filt);
cmat=cmat(:,filt);
R=size(betamat,2);

% T1
%%%%%%%%%%%%%%%%%
x1=[1 m2 1 0 0 0]';%female
x2=[0 m2 1 0 0 0]';%male
expect=[];
tier1=[];
tier2=[];
tier3=[];
prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];
prob6=[];
prob7=[];
prob8=[];
prob9=[];

i=1
for i=1:R %start main loop
    E=diag(Emat(:,i));
    b=betamat(:,i);
    g=gmat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    c4=cmat(4,i);
    c5=cmat(5,i);
    c6=cmat(6,i);
    c7=cmat(7,i);
    %draw gamgi's
    gimat=repmat(g',r2,1)+ mvnrnd(zeros(kr,1),E,r2); %will be r2 by kr
    %create Xb
    int1=sum(repmat(x1',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    int2=sum(repmat(x2',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    
    tier1=[tier1 (m1*(normcdf(c2-int1))+...
              (1-m1)*(normcdf(c2-int2)))'];
    tier2=[tier2 (m1*(normcdf(c5-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c2-int2)))'];     
    tier3=[tier3 (m1*(1-normcdf(c5-int1))+...
             (1-m1)*(1-normcdf(c5-int2)))'];
    
    prob1=[prob1 (m1*(normcdf(-int1))+...
              (1-m1)*(normcdf(-int2)))']; 
    prob2=[prob2 (m1*(normcdf(c1-int1)-normcdf(-int1))+...
              (1-m1)*(normcdf(c1-int2)-normcdf(-int2)))'];
    prob3=[prob3 (m1*(normcdf(c2-int1)-normcdf(c1-int1))+...
              (1-m1)*(normcdf(c2-int2)-normcdf(c1-int2)))'];
    prob4=[prob4 (m1*(normcdf(c3-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c3-int2)-normcdf(c2-int2)))'];
    prob5=[prob5 (m1*(normcdf(c4-int1)-normcdf(c3-int1))+...
             (1-m1)*(normcdf(c4-int2)-normcdf(c3-int2)))'];
    prob6=[prob6 (m1*(normcdf(c5-int1)-normcdf(c4-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c4-int2)))'];
    prob7=[prob7 (m1*(normcdf(c6-int1)-normcdf(c5-int1))+...
             (1-m1)*(normcdf(c6-int2)-normcdf(c5-int2)))'];
    prob8=[prob8 (m1*(normcdf(c7-int1)-normcdf(c6-int1))+...
             (1-m1)*(normcdf(c7-int2)-normcdf(c6-int2)))'];
    prob9=[prob9 (m1*(1-normcdf(c7-int1))+...
             (1-m1)*(1-normcdf(c7-int2)))'];
         
end
proball= [prob1;prob2;prob3;prob4;prob5;prob6;prob7;prob8;prob9];
clear prob1 prob2 prob3 prob4 prob5 prob6 prob7 prob8 prob9;
expect=sum(proball.*repmat((1:1:9)',1,r2*R));
T1probmat=[expect;tier1;tier2;tier3;proball];

% T2
%%%%%%%%%%%%%%%%%%%%%
x1=[1 m2 0 1 0 0]';%female
x2=[0 m2 0 1 0 0]';%male
expect=[];
tier1=[];
tier2=[];
tier3=[];
prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];
prob6=[];
prob7=[];
prob8=[];
prob9=[];

i=1
for i=1:R %start main loop
    E=diag(Emat(:,i));
    b=betamat(:,i);
    g=gmat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    c4=cmat(4,i);
    c5=cmat(5,i);
    c6=cmat(6,i);
    c7=cmat(7,i);
    %draw gamgi's
    gimat=repmat(g',r2,1)+ mvnrnd(zeros(kr,1),E,r2); %will be r2 by kr
    %create Xb
    int1=sum(repmat(x1',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    int2=sum(repmat(x2',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    
    tier1=[tier1 (m1*(normcdf(c2-int1))+...
              (1-m1)*(normcdf(c2-int2)))'];
    tier2=[tier2 (m1*(normcdf(c5-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c2-int2)))'];     
    tier3=[tier3 (m1*(1-normcdf(c5-int1))+...
             (1-m1)*(1-normcdf(c5-int2)))'];
    
    prob1=[prob1 (m1*(normcdf(-int1))+...
              (1-m1)*(normcdf(-int2)))']; 
    prob2=[prob2 (m1*(normcdf(c1-int1)-normcdf(-int1))+...
              (1-m1)*(normcdf(c1-int2)-normcdf(-int2)))'];
    prob3=[prob3 (m1*(normcdf(c2-int1)-normcdf(c1-int1))+...
              (1-m1)*(normcdf(c2-int2)-normcdf(c1-int2)))'];
    prob4=[prob4 (m1*(normcdf(c3-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c3-int2)-normcdf(c2-int2)))'];
    prob5=[prob5 (m1*(normcdf(c4-int1)-normcdf(c3-int1))+...
             (1-m1)*(normcdf(c4-int2)-normcdf(c3-int2)))'];
    prob6=[prob6 (m1*(normcdf(c5-int1)-normcdf(c4-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c4-int2)))'];
    prob7=[prob7 (m1*(normcdf(c6-int1)-normcdf(c5-int1))+...
             (1-m1)*(normcdf(c6-int2)-normcdf(c5-int2)))'];
    prob8=[prob8 (m1*(normcdf(c7-int1)-normcdf(c6-int1))+...
             (1-m1)*(normcdf(c7-int2)-normcdf(c6-int2)))'];
    prob9=[prob9 (m1*(1-normcdf(c7-int1))+...
             (1-m1)*(1-normcdf(c7-int2)))'];
         
end
proball= [prob1;prob2;prob3;prob4;prob5;prob6;prob7;prob8;prob9];
clear prob1 prob2 prob3 prob4 prob5 prob6 prob7 prob8 prob9;
expect=sum(proball.*repmat((1:1:9)',1,r2*R));
T2probmat=[expect;tier1;tier2;tier3;proball];

% T3
%%%%%%%%%%%%%%%%%%%%%
x1=[1 m2 0 0 1 0]';%female
x2=[0 m2 0 0 1 0]';%male
expect=[];
tier1=[];
tier2=[];
tier3=[];
prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];
prob6=[];
prob7=[];
prob8=[];
prob9=[];

i=1
for i=1:R %start main loop
    E=diag(Emat(:,i));
    b=betamat(:,i);
    g=gmat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    c4=cmat(4,i);
    c5=cmat(5,i);
    c6=cmat(6,i);
    c7=cmat(7,i);
    %draw gamgi's
    gimat=repmat(g',r2,1)+ mvnrnd(zeros(kr,1),E,r2); %will be r2 by kr
    %create Xb
    int1=sum(repmat(x1',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    int2=sum(repmat(x2',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    
    tier1=[tier1 (m1*(normcdf(c2-int1))+...
              (1-m1)*(normcdf(c2-int2)))'];
    tier2=[tier2 (m1*(normcdf(c5-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c2-int2)))'];     
    tier3=[tier3 (m1*(1-normcdf(c5-int1))+...
             (1-m1)*(1-normcdf(c5-int2)))'];
    
    prob1=[prob1 (m1*(normcdf(-int1))+...
              (1-m1)*(normcdf(-int2)))']; 
    prob2=[prob2 (m1*(normcdf(c1-int1)-normcdf(-int1))+...
              (1-m1)*(normcdf(c1-int2)-normcdf(-int2)))'];
    prob3=[prob3 (m1*(normcdf(c2-int1)-normcdf(c1-int1))+...
              (1-m1)*(normcdf(c2-int2)-normcdf(c1-int2)))'];
    prob4=[prob4 (m1*(normcdf(c3-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c3-int2)-normcdf(c2-int2)))'];
    prob5=[prob5 (m1*(normcdf(c4-int1)-normcdf(c3-int1))+...
             (1-m1)*(normcdf(c4-int2)-normcdf(c3-int2)))'];
    prob6=[prob6 (m1*(normcdf(c5-int1)-normcdf(c4-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c4-int2)))'];
    prob7=[prob7 (m1*(normcdf(c6-int1)-normcdf(c5-int1))+...
             (1-m1)*(normcdf(c6-int2)-normcdf(c5-int2)))'];
    prob8=[prob8 (m1*(normcdf(c7-int1)-normcdf(c6-int1))+...
             (1-m1)*(normcdf(c7-int2)-normcdf(c6-int2)))'];
    prob9=[prob9 (m1*(1-normcdf(c7-int1))+...
             (1-m1)*(1-normcdf(c7-int2)))'];
        
end
proball= [prob1;prob2;prob3;prob4;prob5;prob6;prob7;prob8;prob9];
clear prob1 prob2 prob3 prob4 prob5 prob6 prob7 prob8 prob9;
expect=sum(proball.*repmat((1:1:9)',1,r2*R));
T3probmat=[expect;tier1;tier2;tier3;proball];


% T4
%%%%%%%%%%%%%%%%%%%%%
x1=[1 m2 0 0 0 1]';%female
x2=[0 m2 0 0 0 1]';%male
expect=[];
tier1=[];
tier2=[];
tier3=[];
prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];
prob6=[];
prob7=[];
prob8=[];
prob9=[];

i=1
for i=1:R %start main loop
    E=diag(Emat(:,i));
    b=betamat(:,i);
    g=gmat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    c4=cmat(4,i);
    c5=cmat(5,i);
    c6=cmat(6,i);
    c7=cmat(7,i);
    %draw gamgi's
    gimat=repmat(g',r2,1)+ mvnrnd(zeros(kr,1),E,r2); %will be r2 by kr
    %create Xb
    int1=sum(repmat(x1',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    int2=sum(repmat(x2',r2,1).*[repmat(b',r2,1) gimat],2); %will be r2 by 1
    
    tier1=[tier1 (m1*(normcdf(c2-int1))+...
              (1-m1)*(normcdf(c2-int2)))'];
    tier2=[tier2 (m1*(normcdf(c5-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c2-int2)))'];     
    tier3=[tier3 (m1*(1-normcdf(c5-int1))+...
             (1-m1)*(1-normcdf(c5-int2)))'];
    
    prob1=[prob1 (m1*(normcdf(-int1))+...
              (1-m1)*(normcdf(-int2)))']; 
    prob2=[prob2 (m1*(normcdf(c1-int1)-normcdf(-int1))+...
              (1-m1)*(normcdf(c1-int2)-normcdf(-int2)))'];
    prob3=[prob3 (m1*(normcdf(c2-int1)-normcdf(c1-int1))+...
              (1-m1)*(normcdf(c2-int2)-normcdf(c1-int2)))'];
    prob4=[prob4 (m1*(normcdf(c3-int1)-normcdf(c2-int1))+...
             (1-m1)*(normcdf(c3-int2)-normcdf(c2-int2)))'];
    prob5=[prob5 (m1*(normcdf(c4-int1)-normcdf(c3-int1))+...
             (1-m1)*(normcdf(c4-int2)-normcdf(c3-int2)))'];
    prob6=[prob6 (m1*(normcdf(c5-int1)-normcdf(c4-int1))+...
             (1-m1)*(normcdf(c5-int2)-normcdf(c4-int2)))'];
    prob7=[prob7 (m1*(normcdf(c6-int1)-normcdf(c5-int1))+...
             (1-m1)*(normcdf(c6-int2)-normcdf(c5-int2)))'];
    prob8=[prob8 (m1*(normcdf(c7-int1)-normcdf(c6-int1))+...
             (1-m1)*(normcdf(c7-int2)-normcdf(c6-int2)))'];
    prob9=[prob9 (m1*(1-normcdf(c7-int1))+...
             (1-m1)*(1-normcdf(c7-int2)))'];       
end
proball= [prob1;prob2;prob3;prob4;prob5;prob6;prob7;prob8;prob9];
clear prob1 prob2 prob3 prob4 prob5 prob6 prob7 prob8 prob9;
expect=sum(proball.*repmat((1:1:9)',1,r2*R));
T4probmat=[expect;tier1;tier2;tier3;proball];


kdiag=klausdiagnostics(T1probmat);
Tprobmean=kdiag(:,1);
Tprobstd=kdiag(:,2);
Tprobnse1=kdiag(:,3);
Tprobmeanlo=Tprobmean-1.96*Tprobnse1;
Tprobmeanhi=Tprobmean+1.96*Tprobnse1;
Tprobout=[Tprobmean Tprobstd Tprobnse1 Tprobmeanlo Tprobmeanhi]';
fprintf(fid,'T stuff \n\n');
fprintf(fid,'mean\t\tstd\t\tnse\t\tlo mean\thi mean\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',Tprobout);
fprintf(fid,'\n');

kdiag=klausdiagnostics(T2probmat);
Tprobmean=kdiag(:,1);
Tprobstd=kdiag(:,2);
Tprobnse1=kdiag(:,3);
Tprobmeanlo=Tprobmean-1.96*Tprobnse1;
Tprobmeanhi=Tprobmean+1.96*Tprobnse1;
Tprobout=[Tprobmean Tprobstd Tprobnse1 Tprobmeanlo Tprobmeanhi]';
fprintf(fid,'T stuff \n\n');
fprintf(fid,'mean\t\tstd\t\tnse\t\tlo mean\thi mean\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',Tprobout);
fprintf(fid,'\n');

kdiag=klausdiagnostics(T3probmat);
Tprobmean=kdiag(:,1);
Tprobstd=kdiag(:,2);
Tprobnse1=kdiag(:,3);
Tprobmeanlo=Tprobmean-1.96*Tprobnse1;
Tprobmeanhi=Tprobmean+1.96*Tprobnse1;
Tprobout=[Tprobmean Tprobstd Tprobnse1 Tprobmeanlo Tprobmeanhi]';
fprintf(fid,'T stuff \n\n');
fprintf(fid,'mean\t\tstd\t\tnse\t\tlo mean\thi mean\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',Tprobout);
fprintf(fid,'\n');

kdiag=klausdiagnostics(T4probmat);
Tprobmean=kdiag(:,1);
Tprobstd=kdiag(:,2);
Tprobnse1=kdiag(:,3);
Tprobmeanlo=Tprobmean-1.96*Tprobnse1;
Tprobmeanhi=Tprobmean+1.96*Tprobnse1;
Tprobout=[Tprobmean Tprobstd Tprobnse1 Tprobmeanlo Tprobmeanhi]';
fprintf(fid,'T stuff \n\n');
fprintf(fid,'mean\t\tstd\t\tnse\t\tlo mean\thi mean\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',Tprobout);
fprintf(fid,'\n');

% save stuff for graph
T1=T1probmat(1,:);
T2=T2probmat(1,:);
T3=T3probmat(1,:);
T4=T4probmat(1,:);

save c:\klaus\aaec6564\mlab\worksp\mod11_HOP_Mag_forgraph T1 T2 T3 T4;


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









