%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Highest Posterior Density Interval for SUR model
% using simulated 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\mod3_hpdi.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

% load wage regression results
load c:\klaus\AAEC6564\mlab\worksp\mod2application betamat sig2mat;
R=size(betamat,2);
% pick coefficient for "education" 
b4=betamat(4,:); 
%get HPDI using klaus_hpdi
[low,up]=klaus_hpdi(b4,0.05,1000)
%means alpha=0.05, and we want 1000 bins for empirical frequencies
fprintf(fid,'lower bound for HPDI =\t%6.4f \n',low);
fprintf(fid,'upper bound for HPDI =\t%6.4f \n',up);
fprintf(fid,'\n');

%plot posterior with HPDI
[f1,x1]=ksdensity(b4,'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 distributions for $\beta_4$',...
    'interpreter','latex','fontsize',14);
set(gca,'Xlim',[-0.1 0.25]);
set(gca,'Ylim',[0 20]);
text(-0.09,14,['low = ' num2str(low)]);
text(0.15,14,['up = ' num2str(up)]);


% % Informal "Test" for restrictions involving several parameters
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For an example, lets's test if an additional child under 6
% reduces annual earnings 3 times more than an additonal child 6-18.
Rm=[0 0 0 0 1 -3]; %define R-matrix and c
q=0; %threshold value to be examined

%In Bayesian terms, we're interested in the posterior distribution of 
%beta4 - 3*beta5.

%get posterior density for Rb
Rb=zeros(R,1);
for r=1:R
    br=betamat(:,r);
    Rb(r)=Rm*br;
end

%get HPDI using klaus_hpdi
[low,up]=klaus_hpdi(Rb,0.05,1000)
fprintf(fid,'lower bound for HPDI =\t%6.4f \n',low);
fprintf(fid,'upper bound for HPDI =\t%6.4f \n',up);
fprintf(fid,'\n');

%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 distributions for $\beta_4-3*\beta_5$',...
    'interpreter','latex','fontsize',14);
set(gca,'Xlim',[-2 2]);
set(gca,'Ylim',[0 2]);
text(-1.9,1,['low = ' num2str(low)]);
text(1,1,['up = ' num2str(up)]);




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










