%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating predictive posterior distributions
% from the tobit model: horse adoption 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_tobit_adoption_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_tobit_adoption;

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

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'Example 1: PPD for expected premium \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

%contents of X
% 1     constant  1
% 2     male      0
% 3     female	  1
% 4     bay       0
% 5     black     0
% 6     brown     1
% 7     buckskin  0
% 8     dun	      0
% 9     gray	  0
% 10    pinto	  0
% 11    redroan	  0
% 12    sorrel    0
% 13    age       2.3 (sample mean)
% 14    training1 0

x1=[1 0 1 0 0 1 0 0 0 0 0 0 2.3 0]';%brown, no training
x2=[1 0 1 0 0 0 0 1 0 0 0 0 2.3 0]';%dun, no training
x3=[1 0 1 0 0 1 0 0 0 0 0 0 2.3 1]';%brown,  training
x4=[1 0 1 0 0 0 0 1 0 0 0 0 2.3 1]';%dun,  training

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


for r=1:R
    br=betamat(:,r);
    sig2=sqrt(sig2mat(r));
    yp1(r)=x1'*br+sig2*(normpdf(x1'*br/sig2)./normcdf(x1'*br/sig2)); 
    yp2(r)=x2'*br+sig2*(normpdf(x2'*br/sig2)./normcdf(x2'*br/sig2)); 
    yp3(r)=x3'*br+sig2*(normpdf(x3'*br/sig2)./normcdf(x3'*br/sig2)); 
    yp4(r)=x4'*br+sig2*(normpdf(x4'*br/sig2)./normcdf(x4'*br/sig2)); 
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,'brown, no training\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,'brown, training\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(yp3);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'dun, no training\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(yp4);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'dun, training\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);
[f3,x3]=ksdensity(yp3,'kernel','epanechnikov','npoints',300);
[f4,x4]=ksdensity(yp4,'kernel','epanechnikov','npoints',300);


figure(1);
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k',x4,f4,'-.k');
title('expected premium for horse types');
legend('brown, no train','dun, no train','brown, train','dun,train');
set(gca,'Xlim',[2 9]);
set(gca,'Ylim',[0 3]);


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

x=mean(X)';
yp1=[];
for r=1:R
    br=betamat(:,r);
    sig2=sqrt(sig2mat(r));
    yp1(r)=br(13)*normcdf(x'*br/sig2); 
    
end

kdiag=klausdiagnostics_greater0(yp1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'mrg. effect of age on expected premium\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 age on premium');
set(gca,'Xlim',[-0.2 0.2]);
set(gca,'Ylim',[0 12]);

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

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

yp1=[];
for r=1:R
    br=betamat(:,r);
    sig2=sqrt(sig2mat(r));
    yp1(r)=(x1'*br+sig2*(normpdf(x1'*br/sig2)./normcdf(x1'*br/sig2)))-...
           (x0'*br+sig2*(normpdf(x0'*br/sig2)./normcdf(x0'*br/sig2))); 
    
end

kdiag=klausdiagnostics_greater0(yp1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'mrg. effect of training on premium\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 training  on premium');
set(gca,'Xlim',[0.5 2.5]);
set(gca,'Ylim',[0 3]);



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;


    