Commit 4e1122f0 authored by michel's avatar michel
Browse files

v4: merge changes to dynare_estimation.m

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1183 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 5fa4cb9e
......@@ -113,7 +113,7 @@ while newRank & t < smpl
% options_.diffuse_d=i; %this is buggy
% end
% end new terminiation criteria by M. Ratto
else
elseif Fstar(i,t) > crit
%% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
%% rank(Pinf)=0. [stphane,11-03-2004].
......
......@@ -55,8 +55,7 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
iFinf(:,:,t) = inv(Pinf(mf,mf,t));
Kinf(:,:,t) = T*Pinf(:,mf,t)*iFinf(:,:,t);
a(:,t+1) = T*a(:,t) + Kinf(:,:,t)*v(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
......@@ -85,8 +84,7 @@ while notsteady & t<smpl
K(:,:,t) = T*P(:,mf,t)*iF(:,:,t);
L(:,:,t) = T-K(:,:,t)*Z;
a(:,t+1) = T*a(:,t) + K(:,:,t)*v(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ;
......@@ -106,8 +104,7 @@ while t<smpl
t=t+1;
v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
a(:,t+1) = T*a(:,t) + K_s*v(:,t);
aK(1,:,t+1) = a(:,t+1);
for jnk=2:nk,
for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
end
end
......
......@@ -26,17 +26,24 @@ v = zeros(pp,smpl);
a = zeros(mm,smpl+1);
a1 = a;
aK = zeros(nk,mm,smpl+nk);
Fstar = zeros(pp,smpl);
Finf = zeros(pp,smpl);
if isempty(options_.diffuse_d),
smpl_diff = 1;
else
smpl_diff=rank(Pinf1);
end
Fstar = zeros(pp,smpl_diff);
Finf = zeros(pp,smpl_diff);
Ki = zeros(mm,pp,smpl);
Li = zeros(mm,mm,pp,smpl);
Linf = zeros(mm,mm,pp,smpl);
L0 = zeros(mm,mm,pp,smpl);
Kstar = zeros(mm,pp,smpl);
Linf = zeros(mm,mm,pp,smpl_diff);
L0 = zeros(mm,mm,pp,smpl_diff);
Kstar = zeros(mm,pp,smpl_diff);
P = zeros(mm,mm,smpl+1);
P1 = P;
Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1;
Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1;
Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1;
Pstar1 = Pstar;
Pinf1 = Pinf;
crit = options_.kalman_tol;
......@@ -106,7 +113,7 @@ while newRank & t < smpl
% options_.diffuse_d=i; % this line is buggy!
% end
% end new terminiation criteria by M. Ratto
else
elseif Fstar(i,t) > crit
%% Note that : (1) rank(Pinf)=0 implies that Finf = 0, (2) outside this loop (when for some i and t the condition
%% rank(Pinf)=0 is satisfied we have P = Pstar and F = Fstar and (3) Finf = 0 does not imply that
%% rank(Pinf)=0. [stphane,11-03-2004].
......@@ -149,8 +156,8 @@ while notsteady & t<smpl
P1(:,:,t) = P(:,:,t);
for i=1:pp
v(i,t) = Y(i,t) - a(mf(i),t) - trend(i,t);
Fi(i,t) = P(mf(i),mf(i),t);
Ki(:,i,t) = P(:,mf(i),t) + H(i,i);
Fi(i,t) = P(mf(i),mf(i),t)+H(i,i);
Ki(:,i,t) = P(:,mf(i),t);
if Fi(i,t) > crit
Li(:,:,i,t) = eye(mm)-Ki(:,i,t)*Z(i,:)/Fi(i,t);
a(:,t) = a(:,t) + Ki(:,i,t)*v(i,t)/Fi(i,t);
......@@ -253,3 +260,4 @@ else
end
epsilonhat = Y-alphahat(mf,:)-trend;
a=a(:,1:end-1);
\ No newline at end of file
......@@ -28,13 +28,13 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
return;
end
Q = M_.Sigma_e;
H = M_.H;
for i=1:estim_params_.nvx
k =estim_params_.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i);
end
offset = estim_params_.nvx;
if estim_params_.nvn
H = zeros(nobs,nobs);
for i=1:estim_params_.nvn
k = estim_params_.var_endo(i,1);
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
......@@ -87,6 +87,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
% M_.params(estim_params_.param_vals(i,1)) = xparam1(i+offset);
%end
M_.Sigma_e = Q;
M_.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
......@@ -150,7 +151,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if estim_params_.nvn
if any(any(H ~= 0)) % should be replaced by a flag
if options_.kalman_algo == 1
LIK = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start);
if isinf(LIK) & ~estim_params_.ncn %% The univariate approach considered here doesn't
......
......@@ -87,6 +87,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,g
% M_.params(estim_params_.param_vals(i,1)) = xparam1(i+offset);
%end
M_.Sigma_e = Q;
M_.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
......@@ -152,7 +153,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,g
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if estim_params_.nvn
if any(any(H ~= 0))
if options_.kalman_algo == 1
[LIK, lik] = DiffuseLikelihoodH1(T,R,Q,H,Pinf,Pstar,data,trend,start);
if isinf(LIK) & ~estim_params_.ncn %% The univariate approach considered here doesn't
......
......@@ -70,7 +70,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoo
% -----------------------------------------------------------------------------
% 4. Kalman smoother
% -----------------------------------------------------------------------------
if estim_params_.nvn
if any(any(H ~= 0)) % should be replaced by a flag
if options_.kalman_algo == 1
[alphahat,epsilonhat,etahat,ahat,aK] = DiffuseKalmanSmootherH1(T,R,Q,H,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
if all(alphahat(:)==0)
......
......@@ -3,8 +3,6 @@ function dynare_estimation(var_list_)
global M_ options_ oo_ estim_params_
global bayestopt_ dsge_prior_weight
% temporary fix until M_.H is initialized by the parser
M_.H = [];
options_.varlist = var_list_;
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
......
......@@ -13,7 +13,7 @@ function global_initialization()
options_.slowc = 1;
options_.timing = 0;
options_.gstep = 1e-2;
options_.debug = 0
options_.debug = 0;
options_.initval_file = 0;
options_.Schur_vec_tol = 1e-8; % used to find nonstationary variables
% in Schur decomposition of the
......@@ -134,4 +134,7 @@ function global_initialization()
oo_.endo_simul = [];
oo_.dr = [];
oo_.exo_det_steady_state = [];
oo_.exo_det_simul = [];
\ No newline at end of file
oo_.exo_det_simul = [];
% Variance matrix for measurement errors
M_.H = 0;
\ No newline at end of file
......@@ -38,6 +38,10 @@ function [xparam1,estim_params_,bayestopt_,lb,ub]=set_prior(estim_params_)
bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
end
if nvn
if M_.H == 0
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
end
for i=1:nvn
estim_params_.var_endo(i,1) = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact');
end
......@@ -69,6 +73,10 @@ function [xparam1,estim_params_,bayestopt_,lb,ub]=set_prior(estim_params_)
cellstr(M_.exo_names(estim_params_.corrx(:,2),:))))));
end
if ncn
if M_.H == 0
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
end
xparam1 = [xparam1; estim_params_.corrn(:,3)];
ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
lb = [lb; max(min(estim_params_.corrn(:,4),1),-1)];
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment