%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating predictive posterior distributions
% from the hierarchical normal regression model 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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\mod6_HNRM_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\mod6_HNRM_outage;

r1=size(betafmat,2); %number of draws from gibbs sampler
r2=10; %number of draws of random effects PER original parameter draw
R=r1*r2;

kf=size(betafmat,1);
kr=size(betarmat,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)
% size = 0 (size 2 firm)
% day = 1 (weekday outage)
% time = 1 (daytime outage)
% duration(minutes) = 25

x=[1 0 1 1 25]';
yp1=[]; %will collect predicted draws
yp1alt=[]; %will collect predicted draws, using conditional LH

thoucount=1000;

for i=1:r1;
    E=vec_var(Emat(:,i));
    sig2=sig2mat(i);
    br=betarmat(:,i);
    bf=betafmat(:,i);
    
    %direct approach (using unconditional likelihood)
    muyp=x(1:kf)'*bf + x(kf+1:end)'*br;
    Vyp=x(kf+1:end)'*E*x(kf+1:end)+sig2;
    % draw r2 yi's 
    yp=normrnd(muyp,sqrt(Vyp),r2,1); % r2 by 1
    
   % alternative (and equivalent)approach (using conditional likelihood)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    brp=mvnrnd(br,E,r2); %kr by r2
    muyp=repmat(x(1:kf)'*bf,r2,1) + sum(repmat(x(kf+1:end)',r2,1).*brp,2);
    ypalt=normrnd(muyp,sqrt(sig2),r2,1); % r2 by 1
    
    yp1=[yp1 yp']; %will be 1 by (r1*r2)
    yp1alt=[yp1alt ypalt']; %will be 1 by (r1*r2)
       
        if i== thoucount
          i
          thoucount=thoucount+1000;
        end
    
end


fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics(yp1);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Predicted outage cost, in $1000, 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)
% size = 1 (size 3 firm)
% day = 1 (weekday outage)
% time = 1 (daytime outage)
% duration(minutes) = 25

x=[1 1 1 1 25]';
yp2=[]; %will collect predicted draws
yp2alt=[]; %will collect predicted draws, using conditional LH

thoucount=1000;

for i=1:r1;
    E=vec_var(Emat(:,i));
    sig2=sig2mat(i);
    br=betarmat(:,i);
    bf=betafmat(:,i);
    
    %direct approach (using unconditional likelihood)
    muyp=x(1:kf)'*bf + x(kf+1:end)'*br;
    Vyp=x(kf+1:end)'*E*x(kf+1:end)+sig2;
    % draw r2 yi's 
    yp=normrnd(muyp,sqrt(Vyp),r2,1); % r2 by 1
    
    % alternative (and equivalent)approach (using conditional likelihood)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    brp=mvnrnd(br,E,r2); %kr by r2
    muyp=repmat(x(1:kf)'*bf,r2,1) + sum(repmat(x(kf+1:end)',r2,1).*brp,2);
    ypalt=normrnd(muyp,sqrt(sig2),r2,1); % r2 by 1
    
    yp2=[yp2 yp']; %will be 1 by (r1*r2)
    yp2alt=[yp2alt ypalt']; %will be 1 by (r1*r2)
       
        if i== thoucount
          i
          thoucount=thoucount+1000;
        end
    
end


fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics(yp2);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Predicted outage cost, in $1000, 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)
% size = 1 (size 3 firm)
% day = 1 (weekday outage)
% time = 1 (daytime outage)
% duration(minutes) = 25

x=[1 1 1 1 50]';
yp3=[]; %will collect predicted draws
yp3alt=[]; %will collect predicted draws using cond. LH

thoucount=1000;

for i=1:r1;
    E=vec_var(Emat(:,i));
    sig2=sig2mat(i);
    br=betarmat(:,i);
    bf=betafmat(:,i);
    
    %direct approach (using unconditional likelihood)
    muyp=x(1:kf)'*bf + x(kf+1:end)'*br;
    Vyp=x(kf+1:end)'*E*x(kf+1:end)+sig2;
    % draw r2 yi's 
    yp=normrnd(muyp,sqrt(Vyp),r2,1); % r2 by 1
    
    % alternative (and equivalent)approach (using conditional likelihood)
    brp=mvnrnd(br,E,r2); %kr by r2
    muyp=repmat(x(1:kf)'*bf,r2,1) + sum(repmat(x(kf+1:end)',r2,1).*brp,2);
    ypalt=normrnd(muyp,sqrt(sig2),r2,1); % r2 by 1
    
    yp3=[yp3 yp']; %will be 1 by (r1*r2)
    yp3alt=[yp3alt ypalt']; %will be 1 by (r1*r2)
       
        if i== thoucount
          i
          thoucount=thoucount+1000;
        end
    
end


fprintf(fid,'total number of predictive draws =\t%6.0f \n',R);
fprintf(fid,'\n');
kdiag=klausdiagnostics(yp3);
out=kdiag(:,1:end-3)'; %drop IEF, m* and and CD
fprintf(fid,'Predicted outage cost, in $1000, 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 types via plots, using uncond. LH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[f1,x1]=ksdensity(yp1,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(yp2,'kernel','epanechnikov','npoints',300);
[f3,x3]=ksdensity(yp3,'kernel','epanechnikov','npoints',300);

figure(1);
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('PPD for outage in $1000 using uncond. LH');
legend('ex1','ex2','ex3');
set(gca,'Xlim',[-40 40]);
set(gca,'Ylim',[0 0.12]);

%Compare types via plots, using cond. LH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[f1,x1]=ksdensity(yp1alt,'kernel','epanechnikov','npoints',300);
[f2,x2]=ksdensity(yp2alt,'kernel','epanechnikov','npoints',300);
[f3,x3]=ksdensity(yp3alt,'kernel','epanechnikov','npoints',300);

figure(2);
plot(x1,f1,'-k',x2,f2,'--k',x3,f3,':k');
title('PPD for outage in $1000 using cond. LH');
legend('ex1','ex2','ex3');
set(gca,'Xlim',[-40 40]);
set(gca,'Ylim',[0 0.12]);


save c:\klaus\AAEC6564\mlab\worksp\mod6_HNRM_predict yp1 yp2 yp3;

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;


    