Trend in Bayesian estimation of Laubach and Williams(2003))

Dear community,

I am trying to estimate the natural rate of interest following Laubach and Williams (2003) using Bayesian methods to examine the degree to which priors shape inference on the natural rate (I follow the Dynare application in https://www.federalreserve.gov/econresdata/feds/2015/files/2015077pap.pdf but with the standard specification in L&W(2003)).

I have checked several posts, and also the guide to observation equations.

The issue is that for different countries, such as Brazil, Colombia, and the US, I find very flat natural rates -under different priors for the parameters- that follow a very flat and small (around zero) estimated unit root component of the natural rate (zz or “other factors” in the original paper). In contrast, the standard L&W(2003) paper finds that the unit root component mostly drives the (decreasing) natural rate in the US.

Below is the model code, and attached is the data and code for the three countries.

I wonder what I am missing? Does it relate with the initial conditions of the (diffuse) filter or the steady state? Would the “filter_initial_state” help? Should I specify another trending variable? Is Dynare implictly imposing something? Does it relate with the observation equation?
Dynare_natural.zip (27.8 KB)

Would be grateful for any suggestion. Regards.

The code:

//////////////////////////////////////////////////////////////////////////
%----------------------------------------------------------------
% 0. Housekeeping
%----------------------------------------------------------------
close all;

%----------------------------------------------------------------
% 1. Variable Definition
%----------------------------------------------------------------

var y ystar ygap gg rr rrstar zz pi;
var dyobs piobs rrobs;

varexo e_pi e_zz e_rr;
varexo e_ystar e_ygap e_gg;
varexo e_yobs;

parameters a1 a2 a3 c1 c2 c3;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

a1 = 1.5 ;
a2 = -0.6 ;
a3 = -0.1 ;
c1 = 0.1 ;
c2 = 0.3 ;
c3 = 0.05 ;

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------
model(linear);

% Measurement, Signal equations
ygap = a1ygap(-1) + a2ygap(-2) + a30.5((rr(-1) - rrstar(-1)) + (rr(-2) - rrstar(-2))) + e_ygap;
pi = c1pi(-1) + c2(pi(-2)+pi(-3)+pi(-4))/3 + (1-c1-c2)(pi(-5)+pi(-6)+pi(-7)+pi(-8))/4 + c3ygap(-1) + e_pi;

%Transition equations
rrstar = gg + zz;
zz= zz(-1) + e_zz;
ystar = ystar(-1) + gg(-1) + e_ystar;
gg = gg(-1) + e_gg;
rr= rr(-1) + e_rr;
//leady=ygap(+1);

% Measurement equations to fit the data with measurement error
y=ystar + ygap + e_yobs;
dyobs = 4*(y-y(-1));
piobs = pi;
rrobs = rr;

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

% Steady state
steady_state_model;
pi=0;
rrstar=0;
rr=rrstar;
gg=0;
ygap=0;
ystar=gg;
y=ystar;
zz=rrstar-gg;
dyobs=0;
piobs=pi;
rrobs=rr;
end;

resid;
steady;
check;
model_diagnostics;

shocks;
var e_yobs; stderr 0;
end;

%----------------------------------------------------------------
% 5. Estimation
%----------------------------------------------------------------

estimated_params;
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF

stderr e_ygap , 0.5 , 1e-4, 10,INV_GAMMA_PDF,0.5 ,2;
stderr e_pi , 1 , 1e-4, 10,INV_GAMMA_PDF,2 ,2;
stderr e_ystar, 2 , 1e-4, 10,INV_GAMMA_PDF,2 ,5;
stderr e_zz , 1 , 1e-4, 10,INV_GAMMA_PDF,1 ,5;
stderr e_rr , 1 , 1e-4, 10,INV_GAMMA_PDF,1 ,5;
stderr e_gg , 0.1 , 1e-4, 10,INV_GAMMA_PDF,0.1 ,5;

a1, 1 , -10, 10, NORMAL_PDF, 1 ,5;
a2, -0.6, -10, 10, NORMAL_PDF, -0.6,5;
a3, -0.3, -10, 10, NORMAL_PDF, -0.3,5;
c1, 0.7 , -10, 10, NORMAL_PDF, 0.7 ,5;
c2, 0.2 , -10, 10, NORMAL_PDF, 0.2 ,5;
c3, 0.1 , -10, 10, NORMAL_PDF, 0.1 ,5;

end;

% Observables
varobs dyobs piobs rrobs;

estimation(graph,order=1,optim=(‘MaxIter’,1000,‘Tolfun’,1.0e-06),datafile=‘dataset_COL’,mode_compute=5,mh_replic=10000,mh_nblocks=2,mh_jscale=0.68,mh_drop=0.2,first_obs=20,lik_init=2,prefilter=1,filtered_vars, mode_check,plot_priors =1) rrstar zz gg;

% Plot natural rate
plot([oo_.SmoothedVariables.rr,oo_.SmoothedVariables.rrstar,oo_.SmoothedVariables.gg,oo_.SmoothedVariables.zz] )