%This is based on the Program for Exercise 3 in Chapter 16 of Bayesian
%Econometric Methods (Koop, Poirier and Tobias)
%Bayesian model averaging program
%Use same setup as Fernandez, Ley and Steel (JAE, 2001)
%this version modifies program to impose some exp vars to always be
%included
%Modified in april 2012 to calculate results for some specific models

load climate_new2.txt;
climate=climate_new2;
[n,kraw]=size(climate);



%Choose the dependent variable
y=climate(:,3);
%choose the explanatory variables that are in all models
xraw1=climate(:,4:13);
kalways=size(xraw1,2);
xraw2 = [climate(:,14:20), climate(:,22:23), climate(:,30:31), climate(:,33:34), climate(:,36:41)]; 
kbma=size(xraw2,2);
xraw=[xraw1 xraw2];
bigk=size(xraw,2);


%subtract mean from all regressors as in FLS
mxraw=mean(xraw);
for i=1:bigk
    xraw(:,i)=xraw(:,i) - mxraw(1,i);
end


%I am creating a vector of length bigk which is 1 if explanatory
%variable is included, else equals zero
%molddraw is initialized to contain all explanatory variables
%with t-stats greater than .5

molddraw=zeros(bigk,1);
molddraw(1:kalways,1)=ones(kalways,1);
xold=[ones(n,1) xraw];
xtxinv=inv(xold'*xold);
bhat=xtxinv*xold'*y;
e=y-xold*bhat;
sse=e'*e;
s2=sse/(n-bigk-1);
bcov=s2*xtxinv;
bt=zeros(bigk+1,1);
for i=2:bigk+1
    bt(i,1)=bhat(i,1)/sqrt(bcov(i,i));
    if abs(bt(i,1))>.5
        molddraw(i-1,1)=1;
    end
end

%the following is an alternative thing I checked for Ross late
molddraw=zeros(bigk,1);
molddraw(1:kalways+2,1)=ones(kalways+2,1);
molddraw(17,1)=1;
molddraw(21,1)=1;
molddraw(22,1)=1;
molddraw(27,1)=1;

%ols results
'ols result and t stat'
[bhat, bt]

%Make up the matrix of explanatory variables for model
xold=ones(n,1);
kold=sum(molddraw)+1;
for i=1:bigk
    if molddraw(i,1)>0
        xold=[xold xraw(:,i)];
    end
end
%specify g0 for the g-prior
if n<=(bigk^2)
    g0=1/(bigk^2);
else
    g0=1/n;
end

yty = (y-mean(y))'*(y-mean(y));
xtxinv=inv(xold'*xold);
ymy = y'*y - y'*xold*xtxinv*xold'*y;
g1=g0/(g0+1);
g2=1/(g0+1);
lprobold = .5*kold*log(g1) -.5*(n-1)*log(g2*ymy + g1*yty);
mstart=molddraw;
lprobstart=lprobold;
lprobmax=-9e200;
lprobmin=9e200;
inccount=zeros(bigk,1);
msize=0;
%I am keeping records for top 10 drawn models
%Here initialize this to initial model
top10mod=[molddraw molddraw molddraw molddraw molddraw ...
       molddraw molddraw molddraw molddraw molddraw];
lprobtop10=lprobold*ones(10,1);
top10count=zeros(10,1);

%calculate first and second moment of all coefficients
%Initialize them here
b1mo=zeros(bigk,1);
b2mo=zeros(bigk,1);
%Number of burnin and kept draws
nburn=1000;
nkeep=10000;

nrep=nburn+nkeep;
for irep=1:nrep
   irep

   %choose at random on of the bigk potential explanatory variables
   %if it is already in the model, delete it else add it
   %Based on this, make up candidate model
indch=kalways+round(kbma*rand);
xnew=xold;
mnewdraw=molddraw;
if indch>kalways
if molddraw(indch,1)==1
     isum=0; 
    for i=1:indch
        isum=isum+molddraw(i,1);
    end
    xnew = [xold(:,1:isum) xold(:,isum+2:kold)];
    mnewdraw(indch,1)=0;
    
else
    isum=0; 
    for i=1:indch
        isum=isum+molddraw(i,1);
    end
    xnew = [xold(:,1:isum+1) xraw(:,indch) xold(:,isum+2:kold)];
    mnewdraw(indch,1)=1;
end
end

knew=sum(mnewdraw)+1;
xtxinv=inv(xnew'*xnew);
ymy = y'*y - y'*xnew*xtxinv*xnew'*y;
lprobnew = .5*knew*log(g1) -.5*(n-1)*log(g2*ymy + g1*yty); 

%Now decide whether to accept candidate draw
if log(rand) < (lprobnew - lprobold)
    xold=xnew;
    lprobold=lprobnew;
    molddraw=mnewdraw;
    kold=knew;
end

if irep>nburn
    
    %If new drawn model better than current top 10, add it to list
    for i=1:10
if lprobold>=lprobtop10(i,1)
  if sum(abs(molddraw - top10mod(:,i)))<.09
      break
  end
    if i<10
       lprobtop10(i+1:10,1)=lprobtop10(i:9,1);
       top10mod(:,i+1:10) = top10mod(:,i:9);
       top10count(i+1:10,1)=top10count(i:9,1);
    end
    lprobtop10(i,1)=lprobold;
    top10mod(:,i)=molddraw;
    top10count(i,1)=0;
    break
  end
 
end

   for i=1:10
       temp1=sum(abs(molddraw-top10mod(:,i)));       
       if temp1<.01
           top10count(i,1)=top10count(i,1)+1;
           break
       end
   end

    inccount = inccount + molddraw;
    msize=msize + kold;
    %calculating posterior properties of coefficients means
    %we have to write out full posterior
    Q1inv = (1+g0)*xold'*xold;
    Q0inv=g0*xold'*xold;
        Q1=inv(Q1inv);
        b1= Q1*xold'*y;
        vs2 = (y-xold*b1)'*(y-xold*b1) + b1'*Q0inv*b1;
        bcov = (vs2/(n-2))*Q1;
        %the next bit of this is awkward, needed to find out if variable is 
        %included in the model and, if so, to find out where it is
           summer1=1;
           bmean1=zeros(bigk,1);
           b2mo1=zeros(bigk,1);
        for i=1:bigk
            bc=zeros(1,kold);        
            if molddraw(i,1)==1
              summer1=summer1+1;
              bc(1,summer1)=1;
              bmean=bc*b1;
            bvar =bc*bcov*bc';
          bmean1(i,1)=bmean;
          b2mo1(i,1)=bvar+bmean^2;
           
            end
          
        end
         b1mo=b1mo + bmean1;
            b2mo=b2mo + b2mo1;
      
        if lprobold>lprobmax
            lprobmax=lprobold;
            mmax=mnewdraw;
            b1max = bmean1;
            bvar1= b2mo1 - b1max.^2;
            bsdmax = sqrt(bvar1);
        end
       if lprobold<lprobmin
            lprobmin=lprobold;
            
        end
end

end

inccount=inccount./nkeep;
'proportion of models visited containing each variable'
varcount=cumsum(ones(bigk,1));
[varcount inccount]
'mean number of regressors in models'
msize/nkeep

'Prob of top 10 models, analytical (after deleting rest of models)'
lt1=lprobtop10 - max(lprobtop10);
exp(lt1)/sum(exp(lt1))

'Prob of top 10 models, numerical (after deleting rest of models)'
top10count./sum(top10count)

'Prob of top 10 models out of total number of models, numerical'
sum(top10count)/nkeep

 b1mo=b1mo./nkeep;
 b2mo=b2mo./nkeep;
 bsd=sqrt(b2mo-b1mo.^2);
 'Posterior mean and stand dev of each coefficient: BMA'
 [varcount b1mo bsd]
 'maximum value for log marginal likelihood'
 lprobmax
 'log marg like for initial model'
 lprobstart
 'minimum value for log marginal likelihood'
 lprobmin
 'Posterior mean and stand dev of each coefficient: Selecting single model'
 [varcount b1max bsdmax]

