%--------------------------------------------------------------------------
% This programs uses a version of the Adaptive Metropolis algorithm of
% Haario et al. (2001) to estiamte the term tructure model in NImark (2010).
%
%--------------------------------------------------------------------------

clear all;
clc;
format short;
warning off all;
global Y 

load('YFB'); %Loads the yield data 
Y=YFB;

load('lastxopt'); %
load('thetamode'); %


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Assign starting values and boundaries to parameter vector THETA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set values here...

% THETA=[0.99,0,0.7,0,0,.6,0.5,0.5,0.2,zeros(1,9),0.1;]';
% LB=[-10*ones(1,6),-15*ones(1,12),1e-8;]';
% UB=[10*ones(1,6),15*ones(1,12),1;]';

lastMCMC=[];

% x=xopt;

% ...or load starting values here
load('lastMCMC');
x=lastMCMC(:,end);
llmax=-NoArb3facLL(thetamode)%Displays likelihood at thetamode

%Define boundary of THETA
THETA=[0.99,0,0.7,0,0,.6,0.5,0.5,0.2,zeros(1,9),0.1;]';
LB=[-10*ones(1,6),-15*ones(1,12),1e-8;]';
UB=[10*ones(1,6),15*ones(1,12),1;]';


%--------------------------------------------------------------------------
% Set control parameters fro Metropolis algorithm
%--------------------------------------------------------------------------
%Number of draws (the ones in the outer loop of the MH)
S=2e4;
%Number of draws (the ones in the inner loop of the MH)
s=1e2; %Every s draw of original chain is thus saved)
%Calibrates vcscale for the first iterinitial iterations of S using an
%identity matrix times epseye.
iterinitial=13;
epseye=1e-8;
% A useful small number
eps=1e-3;
%Values for the constants in the adaptive proposal.
sd=1e-4*(2.38/2)^2; % Scalar that governs the acceptance rate
%Every how many iterations the program display the results
dispiter=100;
%Every how many iterations the program updates the scale matrix
calibiter=100;
%Every how many iterations the MCMC is saved
saveMCMC=1000;

%--------------------------------------------------------------------------
% Initializes the MH algorithm
%--------------------------------------------------------------------------

%Set up so that the first candidate draw is always accepted
lpostdraw = -9e+200;
lpostcan = -9e+200;
bdraw=x;
mean0=x;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set initial covariance of proposal density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if max(size(lastMCMC)) >= (length(x)+1);
    SIGMCMC=cov(lastMCMC(:,1:end)');
    vscale=sd*SIGMCMC+sd*epseye*eye(length(x));
else
    vscale=diag(abs(x))*epseye*100;
end

%Store all draws in the following matrices which are initialized here
bb_=zeros(length(x),s);
OutsideProp=zeros(S,1);
SwitchesProp=zeros(S,1);
%Number of draws outside parameter boundaries
q=0;
%Number of switches (acceptances)
pswitch=0;



%--------------------------------------------------------------------------
% MH algorithm starts here
%--------------------------------------------------------------------------

% Start of the outer MH loop.

for iter=1:S
    tic
% Start of the inner MH loop.
    for iter2=1:s
        bcan = bdraw + norm_rnd(vscale);
        if min(bcan > LB)==1
            if min(bcan < UB)==1
                lpostcan = -NoArb3facLL(bcan);
                laccprob = lpostcan-lpostdraw;
            else
                laccprob=-9e+200;
                q=q+1;
            end
        else
            laccprob=-9e+200;
            q=q+1;
        end
        if lpostcan >= llmax;
            llmax=lpostcan;
            thetamode=bcan;
        end
        %Accept candidate draw with log prob = laccprob, else keep old draw
        if log(rand)<laccprob
            lpostdraw=lpostcan;
            bdraw=bcan;
            pswitch=pswitch+1;
        end

        bb_(:,iter2)=bdraw;

       
    end
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
    toc

    disp('acceptance rate and fraction outside bounds');
    [pswitch/(iter*s) q/(iter*s)]


    lastMCMC=[lastMCMC bb_(:,iter2);];

    if length(lastMCMC(1,:)) >= (length(x)+1);
        SIGMCMC=cov(lastMCMC(:,1:end)');
        vscale=sd*SIGMCMC+sd*epseye*eye(length(x));
        if length(lastMCMC(1,:)) >= (length(x)+101);
        SIGold=cov(lastMCMC(:,1:end-100)');
        oldMCMC=lastMCMC(:,1:end-100);
        else
        SIGold=SIGMCMC; 
        oldMCMC=lastMCMC;
        end
        disp('Current, mode, mean, change in mean last 10 000 draws, s.d., change in s.d. of MCMC last 10 000 draws');
        [bdraw, thetamode, mean(lastMCMC')', (mean(lastMCMC')'-mean(oldMCMC')'), diag(SIGMCMC).^.5, (diag(SIGMCMC).^.5)-(diag(SIGold).^.5);]
         disp('Current likelihood and likelihood at mode');[lpostdraw llmax]
    else
    disp('Current mode');thetamode
    disp('Current likelihood at mode');llmax

    end 
    disp('Total number of draws');length(lastMCMC(1,:))*s
    save('lastMCMC','lastMCMC');
    save('thetamode','thetamode');
 
end 

