function [ml,pit,score] = covselection(yoos,y,ylag,x,bhatprior,g0,DEMEAN,nkeep,oos)
% Implements the MC(3) algorithm for BMA as in KOOP(2003), page 170.
% INPUTS - y is a Tx1 vector of dependent variable
%          ylag is a Txl matrix of the lagged dependent variable
%          x is a Txk matrix of additional variables
%          bhatprior is a lx1 and includes prior values for the constant and the
%          autoregressive lags. The code adjusts the dimention of it in
%          various regressions to make the size compatible
%          g is the scalar value of the g in the g prior
%          DEMEAN guides whether the x should be demeaned, standardized or not, this is
%          usually done to cancel out the scale effects in various
%          regressors
%          nkeep specified the number of draws that should be accumulated
%          based on the posterior simulation. The burn-in draws are
%          adjusted such that they include 10% of the total number of draws
%          that are being accumulated
%          oos provides the values of the set of regressors the forecasts
%          should be conditioned on
count = 0;
nburn = 0.1*nkeep;
nrep = nburn + nkeep;

l = size(ylag,2);
oosraw = oos(:,1:l);

% The following procedure is in place to get rid of NaNs.
xraw = [];
for i = 1:size(x,2)
    if (sum(isfinite(x(:,i))==1) >= size(x,1) && sum(isfinite(oos(:,i+l))==1) >= size(oos(:,i+l),1))
        xraw = [xraw x(:,i)];
        oosraw = [oosraw oos(:,i+l)];
    end
end
k = size(xraw,2);

if DEMEAN == 0;
    xraw = xraw;
elseif DEMEAN == 1;
    xraw = demean(xraw);
elseif DEMEAN == 2;
    xraw = standard(xraw);
end

% Initialize vectors to keep track of models and their respecive marginal likelihoods, 
% as well as the first and second moment of the normal distribution.
cumdraw = zeros(k,nrep);
cumlprob = zeros(1,nrep);
ml = zeros(nkeep,1);
pit = zeros(nkeep,1);
score = zeros(nkeep,1);
mo1 = zeros(k+l,nkeep);
mo2 = zeros(1,nkeep);

% Initialize the number of regressors in this case by including the first
% regressor in x in the estimation.
ndraw = zeros(k,1);
ndraw(1,1) = 1;
cumdraw(:,1) = ndraw;

selx = [ylag xraw(:,ndraw == 1)];
selk = sum(ndraw)+l;
lprob = modellikel(y,selx,[bhatprior; zeros(selk-l,1)],g0);
cumlprob(:,1) = lprob;

for i=2:nrep
    % Choose at random one of the k potential explanatory variables
    % If it is already in the model, delete it else add it
    % Based on this, make up a candidate model
    covdraw = unidrnd(k,1);
    
    ndrawnew = ndraw;
    
    if (ndraw(covdraw,1) == 1)
        ndrawnew(covdraw,1) = 0;
    else
        ndrawnew(covdraw,1) = 1;
    end
    
    selx = [ylag xraw(:,ndrawnew == 1)];
    selk = sum(ndrawnew) + l;
    lprobnew = modellikel(y,selx,[bhatprior; zeros(selk-l,1)],g0);
    
    % Accept or reject the new model with a probability of alpha detailed
    % in the paper
    if log(rand) < (lprobnew - lprob)
        count = count + 1;
        lprob = lprobnew;
        ndraw = ndrawnew;
        cumdraw(:,i) = ndraw;  
        cumlprob(:,i) = lprob;
    else
        cumdraw(:,i) = ndraw;
        cumlprob(:,i) = lprob;
    end
    
    if i > nburn
        % Calculate the coefficients and the variance to be used to
        % construct the predictive densities
        
        effx = [ylag xraw(:,ndraw == 1)];
        [beta sigma2] = bmareg(y,effx,[bhatprior; zeros(sum(ndraw),1)],g0);
        
        effindex = [ones(size(ylag,2),1); ndraw == 1];
        
        mo1(effindex == 1,i-nburn) = beta;
        mo2(1,i-nburn) = sigma2;
        
        ml(i-nburn,:) = exp(cumlprob(:,i)');
        pit(i-nburn,:) = normcdf(yoos,oosraw*mo1(:,i-nburn),sqrt(mo2(1,i-nburn)));
        score(i-nburn,:) = normpdf(yoos,oosraw*mo1(:,i-nburn),sqrt(mo2(1,i-nburn)));
    end
end
