%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Chib 1995 marg LH routine for the Hierarchical Normal Regression Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
load c:\klaus\AAEC6564\mlab\worksp\outage_50firms.txt;
data=outage_50firms;
clear outage_50firms;
n=50;
J=6;
obs=n*J;
y=data(:,3)/1000; %outage cost in $1000
% create a 0/1 dummy for size==3
s3=data(:,4);
f=find(s3==2);
s3(f)=0;
f=find(s3==3);
s3(f)=1;
X=[ones(obs,1) s3 data(:,5:7)];
kf=2; %let's specify fixed coefficients for the constant and the size dummy
kr=3; %random coeffs for all others
k=kf+kr;
% break y, Xf, and Xr into individual-specific cells
yimat=mat2cell(y,J*ones(n,1),1);
Ximat=mat2cell(X,J*ones(n,1),k);
% Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=10000; % burn-in
r2=5000; % keepers
R=r1+r2;
% elements for betaf (fixed coefficients)
muf0=zeros(kf,1); %diffuse prior for mean of fixed coefficients
Vf0=eye(kf)*100; % diffuse prior for varcov of fixed coefficients
% elements for betar (hierarchical means of random coefficients
mur0=zeros(kr,1);
Vr0=eye(kr)*100;
% elements for E -hierarchical VCOV
v0=kr+1;
S0=eye(kr);
Edraw=iwishrnd(S0,v0);
% elements for sig2
eta0=1/2;
tau0=1/2; % diffuse prior shape and scale
load c:\klaus\AAEC6564\mlab\worksp\mod6_HNRM_outage;
[logpmarg,lp,llf,lpost]...
=gs_HNRM_chib(y,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
mur0,Vr0,v0,S0,eta0,tau0,betafmat,betarmat,betarimat,Emat,sig2mat);
'GS done'
[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod6_HNRM_outage_chib.txt','w');
if fid==-1;
error('File could not be opened');
else;
disp('File opened successfully');
end;
fprintf(fid,'logpmarg=\t%6.4f \n',logpmarg);
fprintf(fid,'lp=\t%6.4f \n',lp);
fprintf(fid,'llf=\t%6.4f \n',llf);
fprintf(fid,'lpost=\t%6.4f \n',lpost);
save c:\klaus\AAEC6564\mlab\worksp\mod6_HNRM_outage_chib logpmarg lp llf lpost;
st=fclose(fid);
if st==0;
disp('File closed successfully');
else;
warning('Problem with closing file');
end;