%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating predictive posterior distributions
% from the normal regression model with independent priors
% Wetlands application (refined)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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\mod3_wetlands_predictions.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\mod2_wetlands2 betamat sig2mat;
R=size(betamat,2); %number of draws from the GS

%For this application we can interpret predictions of wetland values
%as "Benefit Transfer (BT)" for a specific policy context (wetland type);

fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'BT for 4 surrounding Counties \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');
% Recall contents of X:
% 
% 1     constant                
% 2     log(income)
% 3     users (% of population)
% 4     wetland acres, in units of $1000

x=[1 10.64189 5 4.44]';
yr1=zeros(1,R);

for r=1:R
    br=betamat(:,r);
    sig2r=sig2mat(r);
    yr1(r)=normrnd(x'*br,sqrt(sig2r)); %careful - normrnd requires 
    %the standard deviation as input, not the variance
end
expyr1=exp(yr1);

kdiag=klausdiagnostics(expyr1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Output table for predictions in $$s, 4 Counties\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,'BT for Nevada \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

x=[1 10.81525 1 4.44]';
yr2=zeros(1,R);

for r=1:R
    br=betamat(:,r);
    sig2r=sig2mat(r);
    yr2(r)=normrnd(x'*br,sqrt(sig2r)); %careful - normrnd requires 
    %the standard deviation as input, not the variance
end
expyr2=exp(yr2);

kdiag=klausdiagnostics(expyr2);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Output table for predictions in $$s, Nevada\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,'BT for Nevada & Utah \n'); 
fprintf(fid,'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

x=[1 10.8307 1 4.44]';
yr3=zeros(1,R);

for r=1:R
    br=betamat(:,r);
    sig2r=sig2mat(r);
    yr3(r)=normrnd(x'*br,sqrt(sig2r)); %careful - normrnd requires 
    %the standard deviation as input, not the variance
end
expyr3=exp(yr3);

kdiag=klausdiagnostics(expyr3);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Output table for predictions in $$s, Nevada & Utah\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 BT distributions via plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[f1,x1]=ksdensity(yr1,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(yr2,'kernel','epanechnikov','npoints',300);
[f3,x3]=ksdensity(yr3,'kernel','epanechnikov','npoints',300);


figure(1);
subplot(1,2,1)
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('PPD for wetland values in log form');
legend('4 Counties','NV','NV & UT');
set(gca,'Xlim',[-15 15]);
set(gca,'Ylim',[0 0.35]);


%%%%%%%%%%%%%%%%%%%
% predictions in dollars
%%%%%%%%%%%%%%%%%%%%

[f1,x1]=ksdensity(expyr1,'kernel','epanechnikov','npoints',500);
[f2,x2]=ksdensity(expyr2,'kernel','epanechnikov','npoints',500);
[f3,x3]=ksdensity(expyr3,'kernel','epanechnikov','npoints',200);

subplot(1,2,2)
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('PPD for wetland values in $$s');
legend('4 Counties','NV','NV & UT');
set(gca,'Xlim',[0 50]);
set(gca,'Ylim',[0 0.1]);





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;


    