%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_op_paddlers_predict.txt','w');
if fid==-1;
   warning('File could not be opened');
else;
   disp('File opened successfully');
end;

load c:\klaus\aaec6564\mlab\worksp\paddlers;
y=data(:,4); %let's take a closer look at "beaches"
n=length(y);

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

R=size(betamat,2);
k=size(betamat,1);

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'PPD for example 1 \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');
% define set of characteristics (don't forget the constant term)
% 1     constant    1
% 2     yrs_padd	10
% 3     gender	    0 (male)
% 4     age	        43			
% 5     resident	0 (visitor)
% 6     inc000	    105 (HH income in 1000's - sample mean)
% 7     kids	    1 (kids present in household)
% 8     college	    1 (college education or higher)

x=[1 10 0 43 0 105 1 1]';

prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];

i=1
for i=1:R %start main loop
    b=betamat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    %create Xb
    int=x'*b; %will be r2 by 1
    prob1=[prob1 (normcdf(-int))']; 
    prob2=[prob2 (normcdf(c1-int)-normcdf(-int))'];
    prob3=[prob3 (normcdf(c2-int)-normcdf(c1-int))'];
    prob4=[prob4 (normcdf(c3-int)-normcdf(c2-int))'];
    prob5=[prob5 (1-normcdf(c3-int))'];
           
end

probmat1=[prob1;prob2;prob3;prob4;prob5];

fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics(probmat1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Score probabilities, example 1\n\n');
fprintf(fid,'mean\t\tstd\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'PPD for example 2 \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');
% define set of characteristics (don't forget the constant term)
% 1     constant    1
% 2     yrs_padd	10
% 3     gender	    1 (female)
% 4     age	        43			
% 5     resident	0 (visitor)
% 6     inc000	    105 (HH income in 1000's - sample mean)
% 7     kids	    1 (kids present in household)
% 8     college	    1 (college education or higher)

x=[1 10 1 43 0 105 1 1]';

prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];

i=1
for i=1:R %start main loop
    b=betamat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    %create Xb
    int=x'*b; %will be r2 by 1
    prob1=[prob1 (normcdf(-int))']; 
    prob2=[prob2 (normcdf(c1-int)-normcdf(-int))'];
    prob3=[prob3 (normcdf(c2-int)-normcdf(c1-int))'];
    prob4=[prob4 (normcdf(c3-int)-normcdf(c2-int))'];
    prob5=[prob5 (1-normcdf(c3-int))'];
           
end

probmat2=[prob1;prob2;prob3;prob4;prob5];

fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics(probmat2);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Score probabilities, example 2\n\n');
fprintf(fid,'mean\t\tstd\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'PPD for example 3 \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');
% define set of characteristics (don't forget the constant term)
% 1     constant    1
% 2     yrs_padd	10
% 3     gender	    0 (male)
% 4     age	        43			
% 5     resident	1 
% 6     inc000	    105 (HH income in 1000's - sample mean)
% 7     kids	    0 (kids present in household)
% 8     college	    1 (college education or higher)

x=[1 10 0 43 1 105 0 1]';

prob1=[];
prob2=[];
prob3=[];
prob4=[];
prob5=[];

i=1;
for i=1:R %start main loop
    b=betamat(:,i);
    c1=cmat(1,i);
    c2=cmat(2,i);
    c3=cmat(3,i);
    %create Xb
    int=x'*b; %will be r2 by 1
    prob1=[prob1 (normcdf(-int))']; 
    prob2=[prob2 (normcdf(c1-int)-normcdf(-int))'];
    prob3=[prob3 (normcdf(c2-int)-normcdf(c1-int))'];
    prob4=[prob4 (normcdf(c3-int)-normcdf(c2-int))'];
    prob5=[prob5 (1-normcdf(c3-int))'];
           
end

probmat3=[prob1;prob2;prob3;prob4;prob5];

fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics(probmat3);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Score probabilities, example 3\n\n');
fprintf(fid,'mean\t\tstd\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


%Compare pr(highest score) across types via plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yp1=probmat1(end,:);
yp2=probmat2(end,:);
yp3=probmat3(end,:);

[x1,f1]=epanech_klaus(yp1,min(yp1),max(yp1),300);
[x2,f2]=epanech_klaus(yp2,min(yp2),max(yp2),300);
[x3,f3]=epanech_klaus(yp3,min(yp3),max(yp3),300);
figure(2);
xl=0;
xu=1;
yl=0;
yu=7;
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('PPD for Pr(score=5), beaches, 3 types');
legend('vis, male, kids','vis, fem, kids','res, male, nokids');
axis([xl xu yl yu]);




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









