%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Savage-Dickey Density Ratio for SUR 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;
%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. 2 ($1000)= size1 size2 size3 +eps2
%cost scen. 3 ($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);

% elements for beta
mu0=zeros(k,1); %diffuse prior for mean of betas
V0=eye(k)*100; % diffuse prior for varcov of beta

% elements for E
% E will be 2 by 2
v0=M+1; % number of equations
S0=eye(M); % So these are the diffusest priors possible

load c:\klaus\AAEC6564\mlab\worksp\mod4_outage betamat Emat;

% test restriction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%evaluate prior density under restriction
bigR=[[0 0 0 1 -1 0 -1 1 0];[0 0 0 1 0 -1 -1 0 1]];
q=zeros(2,1);

denom=mvnpdf(q,bigR*mu0,bigR*V0*bigR'); 
%evaluate conditional posterior at all possible values of E
R=size(betamat,2);
nummat=zeros(R,1);
for i=1:R
   Edraw=vec_var(Emat(:,i));%re-construct VCOV E from vector of elements captured in Emat
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=cell2mat(Xcell(j)); % note: round brackets needed for this
       yi=cell2mat(ycell(j));
       int1=int1+Xi'*inv(Edraw)*Xi;
       int2=int2+Xi'*inv(Edraw)*yi;
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   nummat(i)=mvnpdf(q,bigR*mu1,bigR*V1*bigR'); 
   
end
num=mean(nummat);
logSDDR=log(num)-log(denom);


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

fprintf(fid,'logSDDR for equal differences across S5 and S6 for all firm groups=\t%6.4f\n',logSDDR);

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




