function clone
%clone  returns the vector MPE that maximizes the log posterior likelihood, using a cloning method.  
%It also returns MP, the value of the log posterior likelihood at MPE. 
%MPE and MP are  of use as a starting point for  estimation using Bayesian techniques in bayes.m. If priors are uniform, then MPE is  the MLE parameter estimate. The model is a latent-variable semi-structural model with transitory and permanent monetary shocks. 
%The  observables  are  Dy, ff-pai, and  Dff, denoting, respectively,  output growth,  the interest-rate-inflation differential, and the change in the  interest rate. 
%The starting point of the estimation is the vector mxpar produced by probe.m.  
%The output (MPE and MP) is saved in the file clone.mat. 
%The estimation procedure adapts the Data Cloning MCMC approach  proposed in
% "Data cloning: easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods," Lele, Subhash R.; Dennis, Brian, & Lutscher, Fritjof (2007), Ecology Letters, vol. 10, pp. 551-563. 
%See also the evaluation of this method presented in  ``Maximum Likelihood Estimation Using Bayesian Monte Carlo Methods,''  by Akbari, Danial Ali,  Master's Thesis, Lund University, Spring 2015. 
%This script is part of the replication package of the paper   
%Replication materials for: ``What Do Long Data Tell Us About the Permanent Component of Inflation?'' by Stephanie Schmitt-Grohe and Martin
%Uribe. 
%© Martín Uribe, May 2017. 

K = 25; %# of clones, which are introduced gradually through the variable k=1, 2, ll, K
ndraw = 1e5; %number of draws of MCMC in the final stage, when the number of clones is K
nburn = 1e5; %burn-in draws
bdraw = 0.5e5; %# of draws per clone block before the number of clons reaches K
bburn = 0.5e4; %burn-in dras within block

create_structure; %creaes the matrices that constitute the structure of the SVAR;

X = 0*Y(:,1)+1; %intercept

ny = size(Y,2); %# of observable variables

load probe.mat mxpar  mxlp %produced by running run_bayes.m (which in turn calls probe.m which saves probe.mat)

param = mxpar; %parameter draw 
posterior = mxlp; %posterior likelihood  associated with param

np = length(param); %# of estimated parameters

%Initialize chains:
param_chain = zeros(np,ndraw+nburn); %Matrix collecting the history of posterior draws

posterior_chain = zeros(ndraw+nburn,1); %Vector collecting the history of evaluations of the log posterior function

accept_chain =  [];

%The following five objects are used to compute the variance-covariance matrix of the proposal distribution in an adaptive fashion
Mparam = 0; %weighted average of parameter draws, with weight lambda, given below
MSparam = 0; %weighted average of square parameters using the weight lambda. 
Macc = 0.23; %weighted  average of accept rate using weight lambda
V = eye(np); %variance-covariance matrix of proposal distribution
Ci = eye(np);  %Cholesky decomposition of V, V=Ci*Ci';

tar = 0.23; %target for acceptance rate
kappa = 0.001; %speed of adjustment of scaler of proposal distribution 
lambda = (1/2)^1e-4; %weight on current draw in computing  Mparam, MSparam, and Macc. This value puts a weight of 1/2 to the last 1e4 draws. 
n = 1; %counter number of iterations
k = 1; %counter of clone blocks

%start drawing from posterior distribution
while n<=ndraw+nburn

%update the variance-covariance matrix of the proposal distribution 
if mod(n,1000)==0
Vrecent = MSparam-Mparam*Mparam';
try
V = (V+Vrecent)/2;
Ci= chol(V,'lower');
catch
end
end

%Draw parameter value
tryparam = param + kappa*Ci * randn(np,1);
trylikelihood = -inf;
exitflag = 0;

tryprior = get_prior(tryparam,PPD);

if tryprior>-inf  %%Make sure the draw is within bounds

evaluate_structure; %assign values to the matrices of the SVAR

if unit_circle==1 %Check that the  SVAR system produces stable dynamics
try
trylikelihood =  log_like_kalman(Y,F1,Q1,Hp1',R,X,Ap');
end
end

tryposterior =trylikelihood + tryprior;

%Determine whether to accept or reject the draw
if  rand <= exp(k*(tryposterior-posterior))

exitflag = 5;  %indicator of draw accepted

param = tryparam;
likelihood = trylikelihood;
prior = tryprior;
posterior = tryposterior;

end %if rand

end %if tryprior>-inf 

%Update the weighted averages of the parameter draws, of their squares. and of the acceptance rate
Mparam = Mparam * lambda + param * (1-lambda);
MSparam = MSparam * lambda  + param*param' * (1-lambda);
Macc = Macc * lambda + (exitflag==5) * (1-lambda);

%Update the MCMC chains
param_chain(:,n) = param;
posterior_chain(n) = posterior;
accept_chain =  [accept_chain  exitflag==5];

%Update the scaler of the proposal distribution
kappa + 0.001*sign(Macc - tar);
kappa = ans *(ans>0) + 0.9*kappa * (ans<0);
tar_satisfied = abs(Macc-tar)<0.02; %indicator of whether the acceptance target is being hit with a tolerance of 2%

%Training
if mod(n,bdraw+bburn)==0 & (k<=K)

if  tar_satisfied %Before adding a clone, 
     [posterior,j] = max(posterior_chain(1:n));
      param = param_chain(:,j);
    %Display summary statistics
    posterior_chain(bburn+1:bburn+bdraw);
mean(ans);
     k = k+1;
end
    n = 0;
end %if mod(n,bdraw+bburn)==0 &...

n = n+1;
end %while n<ndraw+nburn

%remove burn-in draws
param_chain = param_chain(:,nburn+1:ndraw+nburn);
posterior_chain = posterior_chain(nburn+1:ndraw+nburn);
accept_chain = accept_chain(nburn+1:ndraw+nburn);

%Calculate some first and second moments from chains
mean_param = mean(param_chain,2);
varcov_param = cov(param_chain');
var_param = diag(varcov_param);
std_param = sqrt(var_param);
median_param = median(param_chain,2);
[MP,j] = max(posterior_chain);
MPE = param_chain(:,j);
accept_rate = mean(accept_chain);

tryparam = MPE;
evaluate_structure

save clone.mat MPE MP
