% marginal.m
% 
% Marginal Likelihood estimator
% It computes the marginal likelihood of a model given some data
% following the Harmonic Mean procedure proposed by Gelfand and
% Dey (1994). The weighting function follows Geweke (1998). 
%
% It needs the program:
%    fgeweke.m
%
% Notes:
%    1)The normalization constant in section 4 of the program is
%    model specific and need to recomputed by hand in each case.
%    2)fgeweke (the weighting function) needs to be positive only
%    in the support of the prior (up to computationally relevance).
%
% References:
%    Gelfand, A.E. and D.K. Dey (1994), "Bayesian Model Choice:
%    Asymptotics and Exact Calculations"
%    Journal of the Royal Statistical Society B, 56, pp. 501-514.
%    Geweke, J. (1998), "Using Simulation Methods for Bayesian 
%    Econometric Models: Inference, Development, and Communication" 
%    Staff Report 249, Federal Reserve Bank of Minneapolis.
%    
% Pau Rabanal and Juan Rubio-Ramirez
% New York and Minneapolis, 6-21-2001

clear all;

% loading the ouput from EHL
load EHLoutput
EHL_orginal_size=size(Theta_s,1);
% eliminating the parameters that were fixed.
Theta_EHL=[Theta_s(:,1:4),Theta_s(:,7:15)];
loglike_EHL=loglike_s;
logpri_EHL=logpri_s;
loglikepri_EHL=loglike_EHL+logpri_EHL;
% getting the maximun EHL
max_loglikepri_EHL=max(loglikepri_EHL);
% loading the ouput from GGLS
load GGLSoutput
GGLS_orginal_size=size(Theta_s,1);
% eliminating the parameters that were fixed.
Theta_GGLS=[Theta_s(:,1:3),Theta_s(:,5:13)];
loglike_GGLS=loglike_s;
logpri_GGLS=logpri_s;
loglikepri_GGLS=loglike_GGLS+logpri_GGLS;
% getting the maximun GGLS
max_loglikepri_GGLS=max(loglikepri_GGLS);
% getting the global maximun
max_loglikepri=max(max_loglikepri_GGLS,max_loglikepri_EHL);
% we divide all the draw by the maximun draw
loglikepri_EHL=loglikepri_EHL-max_loglikepri;
likepri_EHL=exp(loglikepri_EHL);
loglikepri_GGLS=loglikepri_GGLS-max_loglikepri;
likepri_GGLS=exp(loglikepri_GGLS);

% calculating statistics necesary for geweke function
mu_EHL    = mean(Theta_EHL);
mu_GGLS    = mean(Theta_GGLS);
sigma_EHL = cov(Theta_EHL);
sigma_GGLS = cov(Theta_GGLS);
prec_EHL  = inv(sigma_EHL);
prec_GGLS  = inv(sigma_GGLS);

index=1;
for p=0.1:0.1:0.9;    
    chi_EHL=chi2inv(1-p,size(Theta_EHL,2));
    marginalcon_EHL=((2*pi)^(-.5*size(Theta_EHL,2))*det(sigma_EHL)^(-.5))/p;
    j=0;
    for i=1:size(Theta_EHL,1),
        %calculating geweke function
        geweke_EHL(i,1)=fgeweke(Theta_EHL(i,:),mu_EHL,prec_EHL,chi_EHL,marginalcon_EHL);
        if geweke_EHL(i)>0
            j=j+1;
        end
    end
    bad_EHL(index)=size(Theta_EHL,1)-j;
    geweke=geweke_EHL./likepri_EHL;
    %armonic mean
    marginal1=1/mean(geweke);
    % log marginal likelihood
    logmarginal_EHL(index)=log(marginal1);
    index=index+1;
end

index=1;
for p=0.1:0.1:0.9;    
    chi_GGLS=chi2inv(1-p,size(Theta_GGLS,2));
    marginalcon_GGLS=((2*pi)^(-.5*size(Theta_GGLS,2))*det(sigma_GGLS)^(-.5))/p;
    j=0;
    for i=1:size(Theta_GGLS,1),
        %calculating geweke function
        geweke_GGLS(i,1)=fgeweke(Theta_GGLS(i,:),mu_GGLS,prec_GGLS,chi_GGLS,marginalcon_GGLS);
        if geweke_GGLS(i)>0
            j=j+1;
        end
    end
    bad_GGLS(index)=size(Theta_GGLS,1)-j;
    geweke=(geweke_GGLS./likepri_GGLS);
    % armonic mean
    marginal1=1/mean(geweke);
    % log marginal likelihood
    logmarginal_GGLS(index)=log(marginal1);
    index=index+1;
end

log_marginal_diff=logmarginal_EHL-logmarginal_GGLS
bad_EHL
bad_GGLS