%==========================================================================
%                       DSGE MODEL ESTIMATION 
%                   Metropolis-Hastings Algorithm
% 
%
% Author: Luigi Bocola         lbocola@sas.upenn.edu
% Date  : 06/16/2013
%==========================================================================


%=========================================================================
%                              HOUSEKEEPING
%=========================================================================
clc
clear all
close all
delete *.asv

tic

l = path;

path('Matlabfiles',path);
path('Minimization Routines',path);
path('Rational Expectations Solution',path);
path('Datafiles',path);


disp('                                                                  ');
disp('    BAYESIAN ESTIMATION OF DSGE MODEL: METROPOLIS-HASTINGS        ');
disp('                                                                  ');

%=========================================================================
%                  METROPOLIS-HASTINGS ALGORITHM 
% (Report the Acceptance Rate and Recursive Averages Every 1000 draws) 
%=========================================================================



Nsim          = input('How Many Posterior Draws?:  ');
disp('                                                                  ');
disp('                                                                  ');

mode      = zeros(Nsim,18);
mu      = zeros(Nsim,18);
mu(1,:) = [0.9 0.1 0.4 4 2 1 10 0.7 0.7 0.7 0 0.3 0.8 0.5 0.1 0 0 0];

Sigma(:,:,1) = eye(18);
% Sigma = zeros(18,18,Nsim);
% Sigma(1,1,1) = 0.25;
% Sigma(2,2,1) = 0.0225;
% Sigma(3,3,1) = 0.04;
% Sigma(4,4,1) = 4;
% Sigma(5,5,1) = 1;
% Sigma(6,6,1) = 0.25;
% Sigma(7,7,1) = 25;
% Sigma(8,8,1) = 0.01;
% Sigma(9,9,1) = 0.01;
% Sigma(10,10,1) = 0.01;
% Sigma(11,11,1) = 0.16;
% Sigma(12,12,1) = 0.0256;
% Sigma(13,13,1) = 0.04;
% Sigma(14,14,1) = 0.2704;
% Sigma(15,15,1) = 0.01;
% Sigma(16,16,1) = 1;
% Sigma(17,17,1) = 1;
% Sigma(18,18,1) = 1;

c             = 0.2;
Nburn         = int32(0.10*Nsim);

mode(1,:) = mu(1,:);
accept        = 0;
obj           = dsgeliki(mode(1,:)) + prior(mode(1,:));
counter       = 0;
logposterior  = obj*ones(Nsim,1);

for i=1:Nsim
    
    Thetac = mvnrnd(mu(i,:),c*Sigma(:,:,i));

prioc=prior(Thetac);    

if prioc==-Inf
    mode(i+1,:)   = mode(i,:);
    logposterior(i+1) = obj;
else

likic = dsgeliki(Thetac);
objc  = prioc+likic;
alpha = exp(objc-obj);


if alpha>1
    mode(i+1,:)   = Thetac;
    accept            = accept+1;
    obj               = objc;
    logposterior(i+1) = objc;
else
    
    mode(i+1,:)   = mode(i,:);
    logposterior(i+1) = obj;
    
end

mu(i+1,:) = mu(i,:) + ((1/i)*(Thetac - mu(i,:)));
Sigma(:,:,i+1) = Sigma(:,:,i) + mu(i,:)*mu(i,:)' - mu(i+1,:)*mu(i+1,:)' + ((1/i)*(Thetac*Thetac' - Sigma(:,:,i) - mu(i,:)*mu(i,:)'));

acceptancerate     = accept/i;
counter            = counter + 1;

if counter==500
disp('                                                                  ');
disp(['                               DRAW NUMBER:', num2str(i)]         );
disp('                                                                  ');
disp('                                                                  ');    
disp(['                           ACCEPTANCE RATE:', num2str(acceptancerate)]);
disp('                                                                  ');
disp('                                                                  ');    
disp('                            RECURSIVE AVERAGES                    ');
disp('                                                                  ');
% disp('   psi1       psi2     rhoR    pistar   rstar    kappa    tau      rhoG     rhoZ     rhoGZ      sigmaR    sigmaG      sigmaZ  ');
% disp(num2str(mean(Thetasim(1:i,:))));
% disp('                                                                  ');

counter = 0;
end
end
end

MODE    = mode(Nsim+1,:)
LOGOFPOSTERIOR= logposterior(Nsim+1);