%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hierarchical normal regression model using outage 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

% load c:\klaus\AAEC6564\mlab\worksp\outage_50firms.txt;
% data=outage_50firms;
% clear outage_50firms;
% save c:\klaus\AAEC6564\mlab\worksp\outage_50firms data;
% 'ok'
% break

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


% Subset of outage data: 50 randomly selected size 2 qnd 3 firms, each with cost reports for 6
% outage scenarios, for a total of 300 obs.

% Contents of data (do "sum" in STATA and copy table - convenient!)

%     Variable |       Obs        Mean    Std. Dev.       Min        Max
% -------------+--------------------------------------------------------
%1      firm_id |       300      110.82    54.66102         24        214
%2     scenario |       300         3.5    1.710679          1          6
%3         cost |       300    1987.033    5493.072          0   50000.01
%4         size |       300         2.6    .4907165          2          3
%5          day |       300    .8333333    .3733007          0          1
%6         time |       300    .8333333    .3733007          0          1
%7      dur_min |       300    190.0056    248.8056   .0333333   720.0001

% There are 28% of zeros, but we'll ignore truncation issues for now.
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)];

% contents of X
%%%%%%%%%%%%%%%%%%%
%1 constant
%2 size dummy
%3 weekday dummy (1=weekday, 0=weekend)
%4 daytime dummy (1=daytime, 0=nighttime)
%5 outage duration in minutes

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 and X
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;

% generic OLS
bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(obs-k);

% 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
% since we're drawing beta first, won't need starting draw

% elements for betar (hierarchical means of random coefficients
mur0=zeros(kr,1); 
Vr0=eye(kr)*100; 
betardraw=bols(kf+1:end);

% elements for E -hierarchical VCOV
v0=kr+1; 
S0=eye(kr);
Edraw=iwishrnd(S0,v0);


% elements for sig2
sig2draw=s2; % use OLS variance as starting draw for Gibbs Sampler
eta0=1/2; %used to be called "v0", but we're already using that for the IW
tau0=1/2; % diffuse prior shape and scale

[betafmat,betarmat,Emat,sig2mat,betarimat]...
    =gs_HNRM_v2(y,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
    mur0,Vr0,betardraw,v0,S0,Edraw,eta0,tau0,sig2draw);


'GS done'


[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod6_HNRM_outage.txt','w');
if fid==-1;
   error('File could not be opened');
   
else;
   disp('File opened successfully');
end;

% put all draws together & run diagnostics
allmat=[betafmat;betarmat;Emat;sig2mat];
kdiag=klausdiagnostics_greater0(allmat);

fprintf(fid,'total number of iterations =\t%6.0f \n',R);
fprintf(fid,'burn-in iterations =\t%6.0f \n',r1);
fprintf(fid,'\n');

% betaf stuff
out=kdiag(1:kf,:)';
fprintf(fid,'Output table for betafs \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

% betar stuff
out=kdiag(kf+1:k,:)';
fprintf(fid,'Output table for betars \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

% E stuff
out=kdiag(k+1:end-1,:)';
fprintf(fid,'Output table for E \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

%sig2 stuff
out=kdiag(end,:)';
fprintf(fid,'Output table for sig2 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


save c:\klaus\AAEC6564\mlab\worksp\mod6_HNRM_outage betafmat betarmat Emat sig2mat betarimat;

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;










