%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;

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
phi0=1/2; 
gam0=1/2;

% elements for sig2
eta0=1/2;
tau0=1/2; % diffuse prior shape and scale

load c:\klaus\AAEC6564\mlab\worksp\mod6_HNRM_outage_nocovs;

[logpmarg,lp,llf,lpost]...
    =gs_HNRM_nocovs_chib(y,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
    mur0,Vr0,phi0,gam0,eta0,tau0,betafmat,betarmat,betarimat,Emat,sig2mat);

'GS done'


[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod6_HNRM_outage_nocovs_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_nocovs_chib logpmarg lp llf lpost;

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









