%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bayesian 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_parallel.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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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+2; % number of equations
S0=eye(M); % So these are the diffusest priors possible
Edraw=iwishrnd(S0,v0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parpool('local',2); %start parallel
% my PC has 12 cores, so I'm using them all.
% You might need to change this to 4,6,8 etc
% depending on your computer
% you can find the number of cores by typing "feature('numCores')" in the
% command window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[betamat,Emat]=gs_sur_parallel(X,Xcell,y,ycell,k1,k2,k3,k,M,n,r1,r2,mu0,...
V0,betadraw,v0,S0,Edraw);
'GS done'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete(gcp); %terminates parallel pool
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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,'number of parallel cores =\t%6.0f \n',feature('numcores'));
fprintf(fid,'\n');
fprintf(fid,'true betas \t%6.3f\n',[b1true;b2true;b3true]);
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'true E \t%6.3f\n',var_vec(Etrue));
fprintf(fid,'\n');
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');
save c:\klaus\AAEC6564\mlab\worksp\mod4_sur_parallel 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;