%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Highest Posterior Density Interval for SUR model
% using outage cost data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set random number generator, start stop watch, open output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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\hpdi_outage.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% load SUR results
load c:\klaus\AAEC6564\mlab\worksp\mod4_outage betamat Emat;

% "Test" 1
%%%%%%%%%%%%%%%%%%%%%
% Test if outage effect is same across scenarios 5 and 6 for firm size 2

Rb=betamat(5,:)-betamat(8,:);

%get HPDI using klaus_hpdi
[low,up]=klaus_hpdi(Rb,0.05,1000);

[f1,x1]=ksdensity(Rb,'kernel','epanechnikov','npoints',300);
%add bounds to graph
ml=repmat(low,1,1000);
mly=(0:1:999);
mu=repmat(up,1,1000);
muy=(0:1:999);

figure(1);
plot(x1,f1,'-k',ml,mly,'r',mu,muy,'r');
title('posterior density for $\beta_5=\beta_7$ with HPDI',...
    'interpreter','latex','fontsize',14);
set(gca,'Xlim',[-0.5 0.5]);
set(gca,'Ylim',[0 8]);
text(-0.45,4,['low = ' num2str(low)]);
text(0.1,4,['up = ' num2str(up)]);

% "Test" 2
%%%%%%%%%%%%%%%%%%%%%
% Examine if outage effect is same for size1 and 2 for outage scenario 1

Rb=betamat(2,:)-betamat(1,:);

%get HPDI using klaus_hpdi
[low,up]=klaus_hpdi(Rb,0.05,1000);
%plot posterior with HPDI

[f1,x1]=ksdensity(Rb,'kernel','epanechnikov','npoints',300);
%add bounds to graph
ml=repmat(low,1,1000);
mly=(0:1:999);
mu=repmat(up,1,1000);
muy=(0:1:999);

figure(2);
plot(x1,f1,'-k',ml,mly,'r',mu,muy,'r');
title('posterior density for $\beta_1=\beta_2$ with HPDI',...
    'interpreter','latex','fontsize',14);
set(gca,'Xlim',[-2 2]);
set(gca,'Ylim',[0 2]);
text(-1.5,1.5,['low = ' num2str(low)]);
text(1,1.5,['up = ' num2str(up)]);



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










