%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bayesian SUR model with 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


% %%
% %Load & prepare data
% %%%%%%%%%%%%%%%%%%%%
% load c:\klaus\AAEC6564\mlab\worksp\outage.txt;
% data=outage;
% clear outage;
% save c:\klaus\AAEC6564\mlab\worksp\outage;
% 'ok'
% break
%
load c:\klaus\AAEC6564\mlab\worksp\outage;
%contents of data

%1  size                firm size (categories 1 to 4)
%2  outage scenario     scen1= 1 hour, daytime, weekday
%                       scen2= moment, daytime, weekday
%                       scen3= 4hrs,   daytime , weekday
%                       scen4= 12hrs, daytime, weekday
%                       scen5= 1 hr, nightitme, weekday
%                       scen6= 1 hr, daytime, weekend
%3  outage cost (1998 dollars)
%4  firm id

% We have 1046 firms, each with 6 outage scenarios for a total of 6276
% observations; data are already sorted by firm and scenario

% We'll focus on outage scenarios 1, 5, and 6 (all 1 hr-scenarios), and
% firm groups 1-3

f=find(data(:,2)==2 | data(:,2)==3 | data(:,2)==4);
data(f,:)=[];
f=find(data(:,1)==4);
data(f,:)=[];

% this leaves us with 964 firms, for a total of 2892 obs

M=3; 
n=964;

y=data(:,3)/1000; %costs in $1000

%generate dummies for firm size
size1=zeros(M*n,1);
f=find(data(:,1)==1);
size1(f)=data(f,1);

size2=zeros(M*n,1);
f=find(data(:,1)==2);
size2(f)=data(f,1);

size3=zeros(M*n,1);
f=find(data(:,1)==3);
size3(f)=data(f,1);

f1=find(data(:,2)==1);
f2=find(data(:,2)==5);
f3=find(data(:,2)==6);

X1=[size1(f1) size2(f1) size3(f1)]; 
X2=[size1(f2) size2(f2) size3(f2)];
X3=[size1(f3) size2(f3) size3(f3)];

y1=y(f1);
y2=y(f2);
y3=y(f3);

k1=size(X1,2);
k2=size(X2,2);
k3=size(X3,2);

k=k1+k2+k3;
M=3; % number of equations

% here is our full SUR model, in summary:
% equ.1:

%cost scen. 1 ($1000)= size1 size2 size3 +eps1
%cost scen. 5 ($1000)= size1 size2 size3 +eps2
%cost scen. 6 ($1000)= size1 size2 size3 +eps3

% create individual-specific panels - useful for Gibbs Sampler
Xcell=cell(n,1);
ycell=cell(n,1);
for i=1:n
    Xcell{i}=[X1(i,:) zeros(1,k2+k3); zeros(1,k1) X2(i,:) zeros(1,k3);...
        zeros(1,k1+k2) X3(i,:)];
    ycell{i}=[y1(i);y2(i);y3(i)];
end
X=cell2mat(Xcell); %convert cell array to full matrix
y=cell2mat(ycell);

% Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% starting values, priors, and tuners
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% general elements
r1=10000; % burn-in
r2=5000; % keepers
R=r1+r2;

% generic OLS
bols1=inv(X1'*X1)*X1'*y1;
bols2=inv(X2'*X2)*X2'*y2;
bols3=inv(X3'*X3)*X3'*y3;

% elements for beta
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta
betadraw=[bols1;bols2;bols3];

% elements for E
% E will be 2 by 2
v0=M+1; 
S0=eye(M);
Edraw=iwishrnd(S0,v0);


[betamat,Emat]=gs_sur(X,Xcell,y,ycell,k1,k2,k3,k,M,n,r1,r2,mu0,...
    V0,betadraw,v0,S0,Edraw);
'GS done'


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

% put all draws together & run diagnostics
allmat=[betamat;Emat];
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');

% beta stuff
out=kdiag(1:k,:)';
fprintf(fid,'Output table for betas \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');

out=kdiag(k+1:end,:)';
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');

%output in correlation form
rho12=Emat(2,:)./sqrt(Emat(1,:).*Emat(3,:));
rho13=Emat(4,:)./sqrt(Emat(1,:).*Emat(end,:));
rho23=Emat(5,:)./sqrt(Emat(3,:).*Emat(end,:));

allmat=[rho12;rho13;rho23];
kdiag=klausdiagnostics_greater0(allmat);
out=kdiag';
fprintf(fid,'Output table for E, correlations \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\mod4_outage betamat Emat;

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;










