function usg_noise_bayes
%Bayesian estimation of the open economy model with stationary, nonstationary, and noise shocks and imperfect information presented in the book Open Economy Macroeconomics, by Martín Uribe and Stephanie Schmitt-Grohé (Princeton University Press) in  the chapter `Business Cycles in Emerging Countries: Productivity Shocks Versus Financial Frictions.'' 
%
%output (save in the file usg_noise_bayes.mat):
%param_chain is a matrix containing 1 million draws from the posterior distribution of the vector of estimated parameters
%mean_param mean of param_chain
%logposterior_chain vector containing 1 million draws from the posterior log lilelihood function
%var_param variance/covariance matrix of the vector of estimated parameters, computed using the 1 million MCMC chain. 
%std_param  lower choleski decomposition of std_param
%std_param0 lower choleski decomposition of the  variance matrix of the stand-in distribution. 
%accept_rate fraction of draws accepted
%exitflag_chain is a 1 million vector of 0s, 1s, 2s, 3s, 4s, or 5s:
%0 the parameter draw delivered no equilibrium in the vicinity of the steady state
%1 the draw was rejected even though the equilibrium was unique (the key feature of the mcmc algorithm)
% 2 the draw was rejected because it implied multiple equilibria
%3 the draw was rejected because it resulted in the message ``invertibility violated.''
%4 the draw was rejected because it fell outside of the support of the prior distribution
%5 the draw was accepted.
%Running time: On a top-of-the-line desktop computer in vintage 2011, this program ran 55 minutes to produce an MCMC of 1 million draws. 
%© M. Uribe and S. Schmitt-Grohé, January 2014. 

%Number of draws of the MCMC chain
ndraw = 1e+6;

%scaling factor of var-cov matrix of normal distribution from which parameter values are drawn
k = xlsread('k.xls'); %this xls file is created by running usg_noise_run.m

%load data
Y1 = xlsread('mexico_data.xls');
%columns: date, growth rate of output, growth rate of consumption, growth rate of investment, trade-balance-to-gdp ratio.
%mexico_data.xls was created by running data_level2growth.m in
%z:\uribe\book\noise_shocks\mfiles
T = size(Y1,1);  %number of observations

%Demean the data and remove date column
Y = Y1(:,2:end) - ones(T,1)*mean(Y1(:,2:end));

%The vector of parameters to be estimated, param,  has the following order: 
%SIGMAX = param_estim(1);
%RHOX = param_estim(2);
%SIGMAZ = param_estim(3);
%RHOZ = param_estim(4);
%SIGMAN = param_estim(5);
%RHON = param_estim(6);
%PHI = param_estim(7);
%STDmey = param_estim(8);
%STDmec = param_estim(9);
%STDmeiv = param_estim(10);
%STDmetby = param_estim(11);


%Upper and lower bounds of supports of parameter distributions
get_d
%The last 4 elements of the second column of D  were obtained as std(Y)'*0.25;  This implies that the upper bound of the prior distribution of the standard  deviation of measurement errors is 25 percent of the standard deviation of the observables.  Equivalently, measurement errors are allowed to explain at most 6.25 percent of the variance of the observable. 

LB = D(:,1); %lower bound
UB = D(:,2);%upper bound

%Initial guess for estimated parameter
param0 = xlsread('param0.xls'); %this xls file is created by  usg_noise_run.m

%Number of estimated parameters
np = length(param0);

%Initial var/cov matrix  of stand-in normal density
std_param0 = xlsread('std_param0.xls'); %this xls file is created by usg_noise_run.m

pick_variables = [1 2 3 4]; %picks the elements of the policy function (rows of gx) corresponding to the observables.

%Initial Prior likelihood
logprior0 = log(prod( duniform(param0,LB,UB)));
 
%Calculate the initial log likelihood of the model for param0
loglikelihood0 =  log_likelihood(param0,pick_variables,Y);
 
%calculate initial posterior density value
logposterior0 =  loglikelihood0 + logprior0;

%Initialize the vectors used in the chain: 

%Matrix collecting the history of posterior draws
param_chain = zeros(np,ndraw);

%Vector collecting history of log posterior distribution
logposterior_chain = zeros(ndraw,1);

exitflag_chain = zeros(ndraw,1); %Vector collecting record of reason for rejection of draw: 0=no equilibrium, 1=draw rejected; 2=multiple equilibriu; 3=invertibility violated; 4=bound violated; 5=draw accepted

%The following two variables initialize the computatiion of the posterior mean and variance of the vector of estimated parameters
mean_param = 0;
var_param = 0;
accept_rate = 0;

%start drawing from posterior distribution

param_current = param0; %latest param eter draw of the chain 

logposterior_current = logposterior0; %latest draw of log posterior likelihood in the chain

for j=1:ndraw

exitflag = 4;

%Draw parameter value
param = param_current + k * std_param0 * randn(np,1);

%Make sure the draw is within bounds

if (param>=LB) & (param <= UB) 


 [LL,exitflag] =  log_likelihood(param,pick_variables,Y);


%Make sure eq'm is locally unique
if exitflag==1

logposterior =LL + logprior0;

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


exitflag = 5; 

param_current = param;

logposterior_current = logposterior;

end %if rand

end %if exitflag

end %^if param>=LB ...

param_chain(:,j) = param_current;

logposterior_chain(j) = logposterior_current;

exitflag_chain(j) = exitflag;

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

mean_param = mean_param * (j-1)/j + param_current/j;

var_param = var_param * (j-1)/j + param_current*param_current'/j;

end %for j=1:ndraw

var_param = var_param - mean_param*mean_param';

std_param = chol(var_param,'lower');

median_param = median(param_chain,2);

eval(['save usg_noise_bayes.mat param_chain logposterior_chain exitflag_chain param0 std_param0  mean_param var_param median_param std_param accept_rate LB UB k Y1 Y pick_variables'])