%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MC3 application to Fernandez et al's growth regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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\mod8_MC3_growth.txt','w');
if fid==-1;
   warning('File could not be opened');
   return
else;
   disp('File opened successfully');
end;

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

load c:\klaus\AAEC6564\mlab\worksp\growth;
%72 countries (n) with 42 regressors

%Variables in data:
%%%%%%%%%%%%%%%%%%%%%%%%
%1      GDP growth          avg per-capita GDP growth between 1960 and 1992
%2      PrSc Enroll         Primary school enrollment, 1960
%3      Life Exp            Life expectancy
%4      GDPsh560            GDP level in 1960
%5      Mining              Fraction GDP in mining
%6      Eco Org             Economic organization
%7      Yrs Open            Number of years open economy
%8      English pct         Frac. pop. speaking English
%9      Foreign pct         Frac. pop. speaking foreign language
%10     R FEX Dist          Exchange rate distortions
%11     Equip Inv           Equipment investment
%12     NEquip Inv          Non-equipment investment
%13     std(BMP)            std (Black Market Premium)
%14     Outwar Or           Outward orientation
%15     Bl Mkt Pm           Black Market Premium
%16     Area                Area (scale effect)
%17     LatAmerica          Latin American indicator
%18     SubSahara           Sub-Saharan indicator
%19     High Enroll         Higher education enrollment
%20     Publ Edu pct        Public education share
%21     Rev & Coup          Revolutions and coups
%22     War Dummy           War indicator
%23     Pol Rights          Policitcal rights
%24     Civl Lib            Civil liberties
%25     Abs Lat             Absolute latitude
%26     Age                 Age (years of existence as independent nation)
%27     Brit Col            British colony indicator
%28     Buddha              Fraction Buddhist
%29     Catholic            Fraction Catholic
%30     Confucious          Fraction Confucian
%31     EthnoL Frac         Ethnolinguistic fractionalization
%32     French Col          French colony indicator
%33     Hindu               Fraction Hindu
%34     Jewish              Fraction Jewish
%35     Muslim              Fraction Muslim
%36     Pr Exports          Primary exports, 1970
%37     Protestants         Fraction Protestant
%38     Rule of Law         Rule of law
%39     Spanish Col         Spanish colony indicator
%40     Pop g               Population growth
%41     Work/Pop            Ratio workers to population
%42     Lab Force           Size of labor force


y=data(:,1); %avg per-capita GDP growth between 1960 and 1992
n=length(y);
Z=[ones(n,1) data(:,2:42)];

X=Z(:,2:end); %drop constant
%de-mean:
xbar=mean(X,1);% mean over rows, 1 by k
X=X-repmat(xbar,n,1);

k=size(X,2); %size of beta (without constant)
n=size(X,1);

%note: Here, I will call the dimension of beta k by 1 and the dimension of 
%theta = [alpha beta']' k+1 by 1 - it makes the code a
%bit easier to follow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preliminaries, Priors, Starting draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r1=100000; % burn-in
r2=1000000; % keepers
R = r1+r2;

%Prior for Beta
g=k^2; %recommended by Fernandez et al (2001) for cases where n < k^2

del0=0.5; %prior probability that a regressor is in model
% Note: in that case all model priors are
% equal and drop out of the MH formula. So del0 is never actually used in 
% the Gibbs Sampler code

%Starting vector for gamma
%%%%%%%%%%%%%%%%%%%%%%%%%
%A comon strategy is to use all variables that have an absolute t-value exceeding a
%certain value, say 0.5, based on an OLS regression
% generic OLS

bols=inv(X'*X)*X'*y;
res=y-X*bols;
s2=(res'*res)/(n-k);
Vbols=s2*inv(X'*X);
se=sqrt(diag(Vbols));
t=bols./se;

f=find(abs(t)>0.5);
gamstart=zeros(k,1);
gamstart(f)=1;

[amat,betamat,sig2mat,gammat,acceptcount]=...
    gs_MC3(X,y,k,r1,r2,gamstart,g);

'GS done'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Capture key results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Run diagnostics on beta and sigma
allmat=[amat;betamat;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');


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

%sig2 stuff
out=kdiag(k+2,:)';
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');

%gamma stuff
out=[mean(gammat,2)]; %inclusion probabilities for each beta
fprintf(fid,'Inclusion probabilities \n\n');
fprintf(fid,'%6.3f\n',out);
fprintf(fid,'\n');

save c:\klaus\AAEC6564\mlab\worksp\m8_MC3_growth betamat sig2mat gammat X y;

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