%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating predictive posterior distributions
% from the probit model for Fair's extramarital affairs example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

tic;  % start stop watch

[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod5_probit_Fair_predict.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% load posterior results
load c:\klaus\AAEC6564\mlab\worksp\mod5_probit_Fair;

R=size(betamat,2); %number of draws from gibbs sampler
k=size(betamat,1);

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'Example 1: PPD for prob "1" \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

%settings for x1
%1       constant       1
%2       male           1
%3       years married  10
%4       kids           1
%5       religious      0
%6       education      16
%7       happy          0

%x2 is the same except for female
x1=[1 1 10 1 0 16 0]';
x2=[1 0 10 1 0 16 0]';

yp1=[]; %will collect predicted draws
yp2=[];


for r=1:R
    br=betamat(:,r);
    yp1(r)=normcdf(x1'*br); 
    yp2(r)=normcdf(x2'*br);
end

fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics_greater0(yp1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Probability of affair, male\n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

kdiag=klausdiagnostics_greater0(yp2);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Probability of affair, female\n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

%Compare types via plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[f1,x1]=ksdensity(yp1,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(yp2,'kernel','epanechnikov','npoints',300);


figure(1);
plot(x1,f1,'-k',x2,f2,'--k');
title('Prob. of affair');
legend('male','female');
set(gca,'Xlim',[0 1]);
set(gca,'Ylim',[0 9]);


fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'Example 2: Marg effect of Years married \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

x=mean(X)';
yp1=[];
for r=1:R
    br=betamat(:,r);
    yp1(r)=br(3)*normpdf(x'*br); 
    
end

kdiag=klausdiagnostics_greater0(yp1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'mrg. effect of yrs married on Probability of affair\n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

[f1,x1]=ksdensity(yp1,'kernel','epanechnikov','npoints',300);
figure(2);
plot(x1,f1,'-k');
title('Marg. effect of years married on Prob. of affair');
set(gca,'Xlim',[-0.01 0.03]);
set(gca,'Ylim',[0 110]);

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'Example 3: Marg effect of "kids" \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

x0=mean(X)';
x0(4)=0;
x1=mean(X)';
x1(4)=1;

yp1=[];
for r=1:R
    br=betamat(:,r);
    yp1(r)=normcdf(x1'*br)-normcdf(x0'*br); 
    
end

kdiag=klausdiagnostics_greater0(yp1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'mrg. effect of kids on Probability of affair\n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

[f1,x1]=ksdensity(yp1,'kernel','epanechnikov','npoints',300);
figure(3);
plot(x1,f1,'-k');
title('Marg. effect of "kids"  on Prob. of affair');
set(gca,'Xlim',[-0.2 0.3]);
set(gca,'Ylim',[0 10]);


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;


    