%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Chib 1995 marg LH routine for the SUR 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

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

% Load & prepare data
%%%%%%%%%%%%%%%%%%%%%%
load c:\klaus\AAEC6564\mlab\worksp\sur_data;
n=length(y1); %since all equations have same # obs's

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

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

% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=10000; % won't be needed, no burn-in necessary for chib method
r2=5000; % needs to be same setting as for the original GS

% 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_sur betamat Emat;

[logpmarg,lp,llf,lpost]...
    =gs_sur_chib(X,Xcell,y,ycell,k1,k2,k3,k,M,n,r1,r2,mu0,...
    V0,v0,S0,betamat,Emat);

'GS done'

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\sur_chib logpmarg lp llf lpost;

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









