%probe.m
%Use a Metropolis-Hastings Sampling method to estimate the parameters [PHIM PHIX PHIN PSSI] in the MXN model with tot and interest-rate spread shocks developed in the paper ``How Important Are Terms of Trade Shocks?,'' by Stephanie Schmitt-Grohé and Martín Uribe, submitted to the International Economic Review. The algorithm minimizes  the distance between the empirical and theoretical impulse responses (IRs) targeted  in the estimation.    The IRs  targeted are those associated with a terms-of-trade shock and the country spread shock  in the SVAR model with both shocks  studied in the aforementioned paper and the variables included are the trade balance, output, consumption,   investment and the real exchange rate. The output of probe.m is stored in the file probe.mat. The 38x4 matrix MXPAR contains the estimates of the 4 parameter values [PHIM  PHIX  PHIN  PSSI]  for each of the 38 countries in the panel. The 38x1 vector MXLL contains the minimum distance attained.  This program takes about 2 days to run for NROUND=1e3 in 2016. 
%© Martín Uribe  October 2016.

clear all 
format compact
clc

dir_info = cd %saves the location of this file
datestr(now)
warning off

load z:\uribe\tot\ier\mfiles\sgu_tot_run_baaff.mat  hx  PI country_name  IR IRs
%produced by running 
%sgu_tot_run_baaff.m in
%z:\uribe\tot\ier\mfiles

load z:\uribe\tot\ier\mfiles\bootstrap_ir_baaff.mat IR1 IR2 IRs1 IRs2
%produced by running 
%bootstrap_ir_baaff.m in
%z:\uribe\tot\ier\mfiles

%pick variables in IR's
vir =[3:7];

%Pick period  in IR's
per = 1:5;
IRstd = abs(IR2 - IR1);
IRsstd = abs(IRs2 - IRs1);

%pick out IRs of tb y and i, in this order
LB = 1e-4*ones(4,1); %lower bound
UB = 80*ones(4,1);%upper bound

%standard deviation of stand-in distributioin
[~,~,prior_std] = duniform(LB,LB,UB);

np = length(LB); %# of estimated parameters
ncou = size(hx,3); %# of countries
ndraws = 1e3;%# of draws per country per round
NROUND = 1e3; %# of rounds
nround = 0; %# current of rounds
tar = 0.5;%target acceptance rate 

%Variables where results are stored
if ~exist('probe.mat')
MXLL = -inf(ncou,1); %likelihood
MXPAR = ones(ncou,np)*20; %parameter values
else
load probe.mat MXPAR MXLL nround
end

while nround<NROUND
for k=1:ncou
mxpar =rand(np,1).*UB;
mxll = -inf;
n = 1; %# of iterations
P = cell(0);
P = {hx(1:2,1:2,k), PI(1:2,1:2,k), IR(per,vir,k), IRstd(per,vir,k), IRs(per,vir,k), IRsstd(per,vir,k)};

while n<ndraws 
r1 = 10.^(- randi(3,np,1));
r2 = sign(randn(np,1));
 param = mxpar .*(1+r1.*r2);  %Draw parameter value

flag_like = 1;
if sum(param>=LB & param<=UB)==np

try
[loglikelihood,flag_like] =  get_dist(param,P);
catch
error('An error has occurred, possibly with the qz function. By M. Uribe')
loglikelihood = NaN;
flag_like = 0;
end

if  (rand<exp(loglikelihood-mxll)) & (flag_like==1)
mxll =  loglikelihood;
mxpar = param(:);
end %if loglikelihood>mxll

%Check if objective is improved
if (loglikelihood>MXLL(k)) & (flag_like==1)
MXLL(k) = loglikelihood;
MXPAR(k,:) = param(:)';
end

end %if sum(param>=LB & param<=UB)==np

n = n+flag_like;
end %while  n<ndraws 

end %for k=1:ncou

if rand<0.1
save probe.mat
end
nround = nround + 1;


%display progress
disp(['Number of rounds cmpleted: ' num2str(nround) ' out of ' num2str(NROUND)]);

end %while nround<NROUND
save probe.mat