Commit bbd95b3a authored by MichelJuillard's avatar MichelJuillard
Browse files

adding Ed Herbst fast implementation of the Kalman filter and test

cases with timing. Still needs preprocessor interface (option) and documentation.
parent 78f8f045
function [LIK, LIKK, a, P] = kalman_filter_fast(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P)
% [LIK, LIKK, a, P] =
% kalman_filter_fast(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P)
% computes the likelihood of a stationnary state space model using Ed
% Herbst fast implementation of the Kalman filter.
%@info:
%! @deftypefn {Function File} {[@var{LIK},@var{likk},@var{a},@var{P} ] =} kalman_filter_fast(@var{Y}, @var{start}, @var{last}, @var{a}, @var{P}, @var{kalman_tol}, @var{riccati_tol},@var{presample},@var{T},@var{Q},@var{R},@var{H},@var{Z},@var{mm},@var{pp},@var{rr},@var{Zflag},@var{diffuse_periods})
%! @anchor{kalman_filter}
%! @sp 1
%! Computes the likelihood of a stationary state space model, given initial condition for the states (mean and variance).
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item Y
%! Matrix (@var{pp}*T) of doubles, data.
%! @item start
%! Integer scalar, first period.
%! @item last
%! Integer scalar, last period (@var{last}-@var{first} has to be inferior to T).
%! @item a
%! Vector (@var{mm}*1) of doubles, initial mean of the state vector.
%! @item P
%! Matrix (@var{mm}*@var{mm}) of doubles, initial covariance matrix of the state vector.
%! @item kalman_tol
%! Double scalar, tolerance parameter (rcond, inversibility of the covariance matrix of the prediction errors).
%! @item riccati_tol
%! Double scalar, tolerance parameter (iteration over the Riccati equation).
%! @item presample
%! Integer scalar, presampling if strictly positive (number of initial iterations to be discarded when evaluating the likelihood).
%! @item T
%! Matrix (@var{mm}*@var{mm}) of doubles, transition matrix of the state equation.
%! @item Q
%! Matrix (@var{rr}*@var{rr}) of doubles, covariance matrix of the structural innovations (noise in the state equation).
%! @item R
%! Matrix (@var{mm}*@var{rr}) of doubles, second matrix of the state equation relating the structural innovations to the state variables.
%! @item H
%! Matrix (@var{pp}*@var{pp}) of doubles, covariance matrix of the measurement errors (if no measurement errors set H as a zero scalar).
%! @item Z
%! Matrix (@var{pp}*@var{mm}) of doubles or vector of integers, matrix relating the states to the observed variables or vector of indices (depending on the value of @var{Zflag}).
%! @item mm
%! Integer scalar, number of state variables.
%! @item pp
%! Integer scalar, number of observed variables.
%! @item rr
%! Integer scalar, number of structural innovations.
%! @item Zflag
%! Integer scalar, equal to 0 if Z is a vector of indices targeting the obseved variables in the state vector, equal to 1 if Z is a @var{pp}*@var{mm} matrix.
%! @item diffuse_periods
%! Integer scalar, number of diffuse filter periods in the initialization step.
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item LIK
%! Double scalar, value of (minus) the likelihood.
%! @item likk
%! Column vector of doubles, values of the density of each observation.
%! @item a
%! Vector (@var{mm}*1) of doubles, mean of the state vector at the end of the (sub)sample.
%! @item P
%! Matrix (@var{mm}*@var{mm}) of doubles, covariance of the state vector at the end of the (sub)sample.
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{DsgeLikelihood}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{kalman_filter_ss}
%! @end deftypefn
%@eod:
% Copyright (C) 2004-2013 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Set defaults.
if nargin<17
Zflag = 0;
end
if nargin<18
diffuse_periods = 0;
end
if nargin<19
analytic_derivation = 0;
end
if isempty(Zflag)
Zflag = 0;
end
if isempty(diffuse_periods)
diffuse_periods = 0;
end
% Get sample size.
smpl = last-start+1;
% Initialize some variables.
dF = 1;
QQ = R*Q*transpose(R); % Variance of R times the vector of structural innovations.
t = start; % Initialization of the time index.
likk = zeros(smpl,1); % Initialization of the vector gathering the densities.
LIK = Inf; % Default value of the log likelihood.
oldK = Inf;
notsteady = 1;
F_singular = 1;
asy_hess=0;
DLIK=[];
Hess=[];
LIKK=[];
if Zflag
K = T*P*Z';
F = Z*P*Z' + H;
else
K = T*P(:,Z);
F = P(Z,Z) + H;
end
W = K;
iF = inv(F);
Kg = K*iF;
M = -iF;
while notsteady && t<=last
s = t-start+1;
if Zflag
v = Y(:,t)-Z*a;
else
v = Y(:,t)-a(Z);
end
if rcond(F) < kalman_tol
if ~all(abs(F(:))<kalman_tol)
return
else
a = T*a;
P = T*P*transpose(T)+QQ;
end
else
F_singular = 0;
dF = det(F);
likk(s) = log(dF)+transpose(v)*iF*v;
a = T*a+Kg*v;
if Zflag
ZWM = Z*W*M;
ZWMWp = ZWM*W';
M = M + ZWM'*iF*ZWM;
F = F + ZWMWp*Z';
iF = inv(F);
K = K + T*ZWMWp';
Kg = K*iF;
W = (T - Kg*Z)*W;
else
ZWM = W(Z,:)*M;
ZWMWp = ZWM*W';
M = M + ZWM'*iF*ZWM;
F = F + ZWMWp(:,Z);
iF = inv(F);
K = K + T*ZWMWp';
Kg = K*iF;
W = T*W - Kg*W(Z,:);
end
% notsteady = max(abs(K(:)-oldK))>riccati_tol;
oldK = K(:);
end
t = t+1;
end
if F_singular
error('The variance of the forecast error remains singular until the end of the sample')
end
% Add observation's densities constants and divide by two.
likk(1:s) = .5*(likk(1:s) + pp*log(2*pi));
if analytic_derivation,
DLIK = DLIK/2;
dlikk = dlikk/2;
if analytic_derivation==2 || asy_hess,
if asy_hess==0,
Hess = Hess + tril(Hess,-1)';
end
Hess = -Hess/2;
end
end
% Call steady state Kalman filter if needed.
if t <= last
if analytic_derivation,
if analytic_derivation==2,
[tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,D2a,D2T,D2Yss);
else
[tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
analytic_derivation,Da,DT,DYss,asy_hess);
end
likk(s+1:end)=tmp2{1};
dlikk(s+1:end,:)=tmp2{2};
DLIK = DLIK + tmp{2};
if analytic_derivation==2 || asy_hess,
Hess = Hess + tmp{3};
end
else
[tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag);
end
end
% Compute minus the log-likelihood.
if presample>diffuse_periods,
LIK = sum(likk(1+(presample-diffuse_periods):end));
else
LIK = sum(likk);
end
if analytic_derivation,
if analytic_derivation==2 || asy_hess,
LIK={LIK, DLIK, Hess};
else
LIK={LIK, DLIK};
end
LIKK={likk, dlikk};
else
LIKK=likk;
end
// See fs2000.mod in the examples/ directory for details on the model
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
dA = exp(gam);
e = 1;
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
gp_obs = mst/dA;
gy_obs = dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
// this disable the switch to steady state filter
options_.riccati_tol = 0;
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8,plot_priors=0,mode_compute=0);
// See fs2000.mod in the examples/ directory for details on the model
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
dA = exp(gam);
e = 1;
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
gp_obs = mst/dA;
gy_obs = dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady;
check;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, beta_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
end;
varobs gp_obs gy_obs;
options_.solve_tolf = 1e-12;
options_.fast_kalman=1;
estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8,plot_priors=0,mode_compute=0);
// country_nbr must be passed on the command line: dynare swk_0 -Dcountry_nbr=3
@#define countries = [1:country_nbr]
var
@#for c in countries
MC_@{c} E_@{c} EF_@{c} R_KF_@{c} QF_@{c} CF_@{c} IF_@{c} YF_@{c} LF_@{c} PIEF_@{c} WF_@{c} RF_@{c} R_K_@{c} Q_@{c} C_@{c} I_@{c} Y_@{c} L_@{c} PIE_@{c} W_@{c} R_@{c} EE_A_@{c} PIE_BAR_@{c} EE_B_@{c} EE_G_@{c} EE_L_@{c} EE_I_@{c} KF_@{c} K_@{c}
@#endfor
;
varexo
@#for c in countries
E_A_@{c} E_B_@{c} E_G_@{c} E_L_@{c} E_I_@{c} ETA_R_@{c} E_PIE_BAR_@{c} ETA_Q_@{c} ETA_P_@{c} ETA_W_@{c}
@#endfor
;
parameters xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i ;
alpha=.30;
beta=.99;
tau=0.025;
ccs=0.6;
cinvs=.22;
lambda_w = 0.5;
phi_i= 6.771;
sig_c= 1.353;
hab= 0.573;
xi_w= 0.737;
sig_l= 2.400;
xi_p= 0.908;
xi_e= 0.599;
gamma_w= 0.763;
gamma_p= 0.469;
czcap= 0.169;
phi_y= 1.408;
r_pie= 1.684;
r_dpi= 0.14;
rho= 0.961;
r_y= 0.099;
r_dy= 0.159;
rho_a= 0.823;
rho_b= 0.855;
rho_g= 0.949;
rho_l= 0.889;
rho_i= 0.927;
rho_pb= 0.924;
model(linear);
@#for c in countries
CF_@{c} = (1/(1+hab))*(CF_@{c}(1)+hab*CF_@{c}(-1))-((1-hab)/((1+hab)*sig_c))*(RF_@{c}-PIEF_@{c}(1)-EE_B_@{c}) ;
0 = alpha*R_KF_@{c}+(1-alpha)*WF_@{c} -EE_A_@{c} ;
PIEF_@{c} = 0;
IF_@{c} = (1/(1+beta))* (( IF_@{c}(-1) + beta*(IF_@{c}(1)))+(1/phi_i)*QF_@{c})+0*ETA_Q_@{c}+EE_I_@{c} ;
QF_@{c} = -(RF_@{c}-PIEF_@{c}(1))+(1-beta*(1-tau))*((1+czcap)/czcap)*R_KF_@{c}(1)+beta*(1-tau)*QF_@{c}(1) +0*EE_I_@{c} ;
KF_@{c} = (1-tau)*KF_@{c}(-1)+tau*IF_@{c}(-1) ;
YF_@{c} = (ccs*CF_@{c}+cinvs*IF_@{c})+EE_G_@{c} ;
YF_@{c} = 1*phi_y*( alpha*KF_@{c}+alpha*(1/czcap)*R_KF_@{c}+(1-alpha)*LF_@{c}+EE_A_@{c} ) ;
WF_@{c} = (sig_c/(1-hab))*(CF_@{c}-hab*CF_@{c}(-1)) + sig_l*LF_@{c} - EE_L_@{c} ;
LF_@{c} = R_KF_@{c}*((1+czcap)/czcap)-WF_@{c}+KF_@{c} ;
EF_@{c} = EF_@{c}(-1)+EF_@{c}(1)-EF_@{c}+(LF_@{c}-EF_@{c})*((1-xi_e)*(1-xi_e*beta)/(xi_e));
C_@{c} = (hab/(1+hab))*C_@{c}(-1)+(1/(1+hab))*C_@{c}(1)-((1-hab)/((1+hab)*sig_c))*(R_@{c}-PIE_@{c}(1)-EE_B_@{c}) ;
I_@{c} = (1/(1+beta))* (( I_@{c}(-1) + beta*(I_@{c}(1)))+(1/phi_i)*Q_@{c} )+1*ETA_Q_@{c}+1*EE_I_@{c} ;
Q_@{c} = -(R_@{c}-PIE_@{c}(1))+(1-beta*(1-tau))*((1+czcap)/czcap)*R_K_@{c}(1)+beta*(1-tau)*Q_@{c}(1) +EE_I_@{c}*0+0*ETA_Q_@{c} ;
K_@{c} = (1-tau)*K_@{c}(-1)+tau*I_@{c}(-1) ;
Y_@{c} = (ccs*C_@{c}+cinvs*I_@{c})+ EE_G_@{c} ;
Y_@{c} = phi_y*( alpha*K_@{c}+alpha*(1/czcap)*R_K_@{c}+(1-alpha)*L_@{c} ) +phi_y*EE_A_@{c} ;
PIE_@{c} = (1/(1+beta*gamma_p))*((beta)*(PIE_@{c}(1)) +(gamma_p)*(PIE_@{c}(-1))+((1-xi_p)*(1-beta*xi_p)/(xi_p))*(MC_@{c})) + ETA_P_@{c} ;
MC_@{c} = alpha*R_K_@{c}+(1-alpha)*W_@{c} -EE_A_@{c};
W_@{c} = (1/(1+beta))*(beta*W_@{c}(+1)+W_@{c}(-1))+(beta/(1+beta))*(PIE_@{c}(+1))-((1+beta*gamma_w)/(1+beta))*(PIE_@{c})+(gamma_w/(1+beta))*(PIE_@{c}(-1))-(1/(1+beta))*(((1-beta*xi_w)*(1-xi_w))/(((1+(((1+lambda_w)*sig_l)/(lambda_w))))*xi_w))*(W_@{c}-sig_l*L_@{c}-(sig_c/(1-hab))*(C_@{c}-hab*C_@{c}(-1))+EE_L_@{c})+ETA_W_@{c};
L_@{c} = R_K_@{c}*((1+czcap)/czcap)-W_@{c}+K_@{c} ;
R_@{c} = r_dpi*(PIE_@{c}-PIE_@{c}(-1))+(1-rho)*(r_pie*(PIE_@{c}(-1)-PIE_BAR_@{c})+r_y*(Y_@{c}-YF_@{c}))+r_dy*(Y_@{c}-YF_@{c}-(Y_@{c}(-1)-YF_@{c}(-1)))+rho*(R_@{c}(-1)-PIE_BAR_@{c})+PIE_BAR_@{c}+ETA_R_@{c};
E_@{c} = E_@{c}(-1)+E_@{c}(1)-E_@{c}+(L_@{c}-E_@{c})*((1-xi_e)*(1-xi_e*beta)/(xi_e));
EE_A_@{c} = (rho_a)*EE_A_@{c}(-1) + E_A_@{c};
PIE_BAR_@{c} = rho_pb*PIE_BAR_@{c}(-1)+ E_PIE_BAR_@{c} ;
EE_B_@{c} = rho_b*EE_B_@{c}(-1) + E_B_@{c} ;
EE_G_@{c} = rho_g*EE_G_@{c}(-1) + E_G_@{c} ;
EE_L_@{c} = rho_l*EE_L_@{c}(-1) + E_L_@{c} ;
EE_I_@{c} = rho_i*EE_I_@{c}(-1) + E_I_@{c} ;
@#endfor
end;
estimated_params;
@#for c in countries
stderr E_A_@{c},0.543,0.01,4,INV_GAMMA_PDF,0.4,2;
stderr E_PIE_BAR_@{c},0.072,0.001,4,INV_GAMMA_PDF,0.02,10;
stderr E_B_@{c},0.2694,0.01,4,INV_GAMMA_PDF,0.2,2;
stderr E_G_@{c},0.3052,0.01,4,INV_GAMMA_PDF,0.3,2;
stderr E_L_@{c},1.4575,0.1,6,INV_GAMMA_PDF,1,2;
stderr E_I_@{c},0.1318,0.01,4,INV_GAMMA_PDF,0.1,2;
stderr ETA_R_@{c},0.1363,0.01,4,INV_GAMMA_PDF,0.1,2;
stderr ETA_Q_@{c},0.4842,0.01,4,INV_GAMMA_PDF,0.4,2;
stderr ETA_P_@{c},0.1731,0.01,4,INV_GAMMA_PDF,0.15,2;
stderr ETA_W_@{c},0.2462,0.1,4,INV_GAMMA_PDF,0.25,2;
@#endfor
rho_a,.9722,.1,.9999,BETA_PDF,0.85,0.1;
rho_pb,.85,.1,.999,BETA_PDF,0.85,0.1;
rho_b,.7647,.1,.99,BETA_PDF,0.85,0.1;
rho_g,.9502,.1,.9999,BETA_PDF,0.85,0.1;
rho_l,.9542,.1,.9999,BETA_PDF,0.85,0.1;
rho_i,.6705,.1,.99,BETA_PDF,0.85,0.1;
phi_i,5.2083,1,15,NORMAL_PDF,4,1.5;
sig_c,0.9817,0.25,3,NORMAL_PDF,1,0.375;
hab,0.5612,0.3,0.95,BETA_PDF,0.7,0.1;
xi_w,0.7661,0.3,0.9,BETA_PDF,0.75,0.05;
sig_l,1.7526,0.5,5,NORMAL_PDF,2,0.75;
xi_p,0.8684,0.3,0.95,BETA_PDF,0.75,0.05;
xi_e,0.5724,0.1,0.95,BETA_PDF,0.5,0.15;
gamma_w,0.6202,0.1,0.99,BETA_PDF,0.75,0.15;
gamma_p,0.6638,0.1,0.99,BETA_PDF,0.75,0.15;
czcap,0.2516,0.01,2,NORMAL_PDF,0.2,0.075;
phi_y,1.3011,1.001,2,NORMAL_PDF,1.45,0.125;
r_pie,1.4616,1.2,2,NORMAL_PDF,1.7,0.1;
r_dpi,0.1144,0.01,0.5,NORMAL_PDF,0.3,0.1;
rho,0.8865,0.5,0.99,BETA_PDF,0.8,0.10;
r_y,0.0571,0.01,0.2,NORMAL_PDF,0.125,0.05;
r_dy,0.2228,0.05,0.5,NORMAL_PDF,0.0625,0.05;
end;
varobs
@#for c in countries
Y_@{c} C_@{c} I_@{c} E_@{c} PIE_@{c} W_@{c} R_@{c}
@#endfor
;
// set estimation globals
estimation(datafile=pseudo_data,presample=40, first_obs=1, nobs=118, mh_replic=0, mode_compute=0,plot_priors=0);
// get variables
[dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list_, M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_);
// this disable the switch to steady state filter
options_.riccati_tol = 0;
// actual computation for timing
tic;
for i=1:100;
fval = dsge_likelihood(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_,[]);
end;
toc;
// country_nbr must be passed on the command line: dynare swk_1 -Dcountry_nbr=3
@#define countries = [1:country_nbr]
var
@#for c in countries
MC_@{c} E_@{c} EF_@{c} R_KF_@{c} QF_@{c} CF_@{c} IF_@{c} YF_@{c} LF_@{c} PIEF_@{c} WF_@{c} RF_@{c} R_K_@{c} Q_@{c} C_@{c} I_@{c} Y_@{c} L_@{c} PIE_@{c} W_@{c} R_@{c} EE_A_@{c} PIE_BAR_@{c} EE_B_@{c} EE_G_@{c} EE_L_@{c} EE_I_@{c} KF_@{c} K_@{c}
@#endfor
;
varexo
@#for c in countries
E_A_@{c} E_B_@{c} E_G_@{c} E_L_@{c} E_I_@{c} ETA_R_@{c} E_PIE_BAR_@{c} ETA_Q_@{c} ETA_P_@{c} ETA_W_@{c}
@#endfor
;
parameters xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i ;
alpha=.30;
beta=.99;
tau=0.025;
ccs=0.6;
cinvs=.22;
lambda_w = 0.5;
phi_i= 6.771;
sig_c= 1.353;
hab= 0.573;
xi_w= 0.737;
sig_l= 2.400;
xi_p= 0.908;
xi_e= 0.599;
gamma_w= 0.763;
gamma_p= 0.469;
czcap= 0.169;
phi_y= 1.408;
r_pie= 1.684;
r_dpi= 0.14;
rho= 0.961;
r_y= 0.099;
r_dy= 0.159;
rho_a= 0.823;
rho_b= 0.855;
rho_g= 0.949;
rho_l= 0.889;
rho_i= 0.927;
rho_pb= 0.924;
model(linear);
@#for c in countries
CF_@{c} = (1/(1+hab))*(CF_@{c}(1)+hab*CF_@{c}(-1))-((1-hab)/((1+hab)*sig_c))*(RF_@{c}-PIEF_@{c}(1)-EE_B_@{c}) ;
0 = alpha*R_KF_@{c}+(1-alpha)*WF_@{c} -EE_A_@{c} ;
PIEF_@{c} = 0;
IF_@{c} = (1/(1+beta))* (( IF_@{c}(-1) + beta*(IF_@{c}(1)))+(1/phi_i)*QF_@{c})+0*ETA_Q_@{c}+EE_I_@{c} ;
QF_@{c} = -(RF_@{c}-PIEF_@{c}(1))+(1-beta*(1-tau))*((1+czcap)/czcap)*R_KF_@{c}(1)+beta*(1-tau)*QF_@{c}(1) +0*EE_I_@{c} ;