function bayes(param0,posterior0,Ci,kappa_input,nburn,ndraw,nrandpick,progress_report, tar, tartol) 
%bayes(param0,posterior0,Ci,kappa_input,nburn,ndraw,nrandpick,progress_report, tar, tartol)  performs a Bayesian estimation of a latent-variable semi-structural model 
%Inputs:
%param0   initial parameter vector
%posterior0  posterior likelihood associated with param0 (optional, if set to [],  the program computes it).
%Ci Cholesky decomposition of V, V=Ci*Ci' is the variance-covariance matrix of the stand-in distribution. Default Ci=I
%kappa_input = scaler of the stand-in distribution. Default: adjust it to ensure that the acceptance rate equals a target given by the variable tar in the script.
%nburn = number of discarded initial draws. Default 1e5. If the algorithm does not accept any proposed vector  after 1e5 draws, it draws another 1e5 vectors, and so on. 
%ndraw = number of draws after burning. Default 1e6. This chain is not saved and could be somewhat longer or shorter than  ndraw, as the script stops drawing once a chain of length nrandpick has been obtained. 
%nrandpick = length of  the chain of parameeter vectors to be saved. It's a randomly picked subsample of a longer MCMC chain of length ndraw. Default ndraw
%progress_report=1 to display progress report during the program execution, and 0 otherwise (default).
%tar = target accept rate
%tarto = tolerance for deviating from tar
%Output (saved in  bayes.mat): 
%param_chain is a matrix containing nrandpick draws from the posterior distribution of the vector of estimated parameters
%mean_param is the  mean of param_chain
%posterior_chain  chain of posterior likelihood values associated with param_chain 
%varcov_param variance/covariance matrix of the vector of estimated parameters, computed using  param_chain.  
%var_param is a vector equal to the diagonal of  varcov_param. 
%std_param is a vector equal to the square root of var_param
%accept_rate fraction of proposal draws accepted out of ndraw draws. 
%© Martín Uribe and Stephanie Schmitt-Grohé, May 2022. 

format compact

create_structure; %creaes the matrices that constitute the structure of the SVAR;
noposterior0 = 0; %indicator = 1 if posterior0 needs to be calculated
if nargin<1 | isempty(param0)
param0 = PPD(:,1); %prior mean PPD is created by create_structure.m
noposterior0 = 1;
end


if nargin<2 | isempty(posterior0) | noposterior0
posterior0 = -inf;
tryparam = param0;
evaluate_structure
if unit_circle==1
likelihood0 =   log_like_kalman(Y,Fs,Qs,Hps',R,X,Ap',Fs,Qs,Hps',R1,X,Ap',Tbreak);
prior0 = get_prior(param0,PPD)
posterior0 = prior0 + likelihood0;
end  %if unit_circle==1
end %if nargin<2 | isnan(posterior0)

param = param0;
posterior = posterior0;

if nargin<3 | isempty(Ci)
Ci = eye(np); %var-cov matrix of stand-in distribution;  np=# of estimated parameters produced by create_structure.m 
end
Ci_identity = max(max(abs(Ci-eye(size(Ci)))))==0; %indicator = 1 if Ci=I

if nargin<4 | isempty(kappa_input)
kappa = 0.01; %scaler of proposal distribution 
else
kappa = kappa_input;
end

if nargin < 5 | (isempty(nburn))
nburn = 1e5;%burn-in draws
end

if nargin<6 | isempty(ndraw)
ndraw = 1e6; %number of draws of MCMC  chain 
end

if nargin<7 | isempty(nrandpick)
nrandpick = ndraw; %total # of  accepted draw  stored. Default ndraw
end

if nargin<8 | isempty(progress_report)
progress_report = 0;
end

if nargin<9 | isempty(tar)
tar = 0.3; %target accept_rate
end

if nargin<10 | isempty(tartol)
tartol = 0.05; %tolerance target accept_rate
end

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

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

prior_chain = zeros(nrandpick,1); %Vector collecting the history of evaluations of the log posterior function

accept_chain =  zeros(nrandpick,1);
accept_rate = 0;

kappa_chain =  zeros(nrandpick,1);

n = 1; %counter number of iterations
nn = 0; %# of accepted draws randomly picked 

%start drawing from posterior distribution
while   (nn<nrandpick)
%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,Fs,Qs,Hps',R,X,Ap',Fs,Qs,Hps',R1,X,Ap',Tbreak);
end %try
end %if unit_circle

tryposterior =trylikelihood + tryprior;

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

exitflag = 5;  %indicator of draw accepted
param = tryparam;
posterior = tryposterior;
prior = tryprior;
end %if rand

end %if tryprior>-inf 

accept_rate = (n-1)/n * accept_rate + 1/n * (exitflag==5);

%Restart MCMC if accepted 0 draws during burn in

TL(n,1) = trylikelihood;
TP(n,1) = tryprior;
if (n==nburn) & (posterior == posterior0)
n = 1;
disp('n=nburn reached without update')
kappa = 0.9*kappa
end

%Update the MCMC chain
if (n>nburn) & (rand<=nrandpick/ndraw)
nn = nn+1;
param_chain(:,nn) = param;
posterior_chain(nn) = posterior;
prior_chain(nn) = prior;
accept_chain(nn) = (exitflag==5);
kappa_chain(nn) = kappa;

%Progress Report
if rand<1/(nrandpick/100) & progress_report
disp('-------------------------------------------')
disp('Progress report')
disp(['Directory: ' cd])
disp(datestr(now)) 
disp(['completed ' num2str(round(nn/nrandpick*100)) '%']);
disp(['accept_rate = ' num2str(accept_rate)]);
save bayes_temp.mat
end %if rand<1/(nrandpick/100)
end %if  (n>nburn) & (rand<=nrandpick/ndraw)

n = n+1;

%Update kappa if not hitting accept_rate target
if (nn>=nrandpick/10) & (accept_rate> tar+tartol)
kappa = kappa*1.1;
disp(['accept_rate = ' num2str(accept_rate)  ' ==> increased kappa  to '  num2str(kappa)])
n = 1; nn=0; 
elseif (nn>=nrandpick/10) & (accept_rate< tar-tartol)
kappa = kappa/1.1^rand;
disp(['accept_rate = ' num2str(accept_rate)  ' ==> reduced kappa  to '  num2str(kappa)])
n = 1; nn=0; 
end

%Update Ci 
if (nn==nrandpick) & (abs(accept_rate- tar)   <= tartol) & (Ci_identity ==1)
save bayes0.mat
varcov_param = cov(param_chain'); 
Ci = chol(varcov_param,'lower');
disp('Update Ci')
n = 1; nn=0; Ci_identity = 0; kappa = 0.1;
end

end %while n<ntotal

%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);
[mode_posterior,j] = max(posterior_chain);
mode_param = param_chain(:,j); %model of posterior distribution
accept_rate = mean(accept_chain);

param = param_chain(:, nn);
posterior = posterior_chain(nn); 


%evaluate matrices of the SVAR at the mean parameter value
tryparam = mean_param;
evaluate_structure

!del bayes_temp.mat
save bayes.mat
