diff --git a/matlab/+occbin/DSGE_smoother.m b/matlab/+occbin/DSGE_smoother.m
index 265f631cb94878447a8309a742e76b5bdf9c88ca..2a5d919d53a64070ed381babd76c51f5ec127b3f 100644
--- a/matlab/+occbin/DSGE_smoother.m
+++ b/matlab/+occbin/DSGE_smoother.m
@@ -185,7 +185,9 @@ while is_changed && maxiter>iter && ~is_periodic
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,M_,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options,TT,RR,CC);% T1=TT;
sto_etahat(iter)={etahat};
regime_history0(iter,:) = regime_history;
- save('info1','regime_history0');
+ if occbin_smoother_debug
+ save('info1','regime_history0');
+ end
sto_CC = CC;
sto_RR = RR;
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 8922a2b5a37b01aad94bd924a3ae5893ab9fbdfc..18d0f286c08fba86940de23fdf9c5aa18e8be786 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -319,7 +319,7 @@ if kalman_algo == 2 || kalman_algo == 4
a_initial = zeros(np,1);
a_initial=set_Kalman_smoother_starting_values(a_initial,M_,oo_,options_);
a_initial=ST*a_initial; %set state prediction for first Kalman step;
- [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
+ [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp,state_uncertainty, aahat, eehat, d, regimes_,TT,RR,CC,TTx,RRx,CCx] = missing_DiffuseKalmanSmootherH3_Z(a_initial,ST, ...
Z,R1,Q,diag(H), ...
Pinf,Pstar,data1,vobs,np,smpl,data_index, ...
options_.nk,kalman_tol,diffuse_kalman_tol, ...
@@ -386,9 +386,11 @@ else
else
opts_simul = options_.occbin.simul;
end
+ % reconstruct smoothed variables
aaa=zeros(M_.endo_nbr,gend);
aaa(oo_.dr.restrict_var_list,:)=alphahat;
- for k=2:gend
+ iTx = zeros(size(TTx));
+ for k=1:gend
if isoccbin
A = TT(:,:,k);
B = RR(:,:,k);
@@ -396,20 +398,53 @@ else
else
C=0;
end
- aaa(:,k) = C+A*aaa(:,k-1)+B*etahat(:,k);
+ iT = pinv(TTx(:,:,k));
+ % store pinv
+ iTx(:,:,k) = iT;
+ Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
+ Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
+ Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
+ AS = Tstar*iT;
+ BS = Rstar-AS*RRx(:,:,k);
+ CS = Cstar-AS*CCx(:,k);
+ static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
+ ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
+ static_var_list0 = static_var_list;
+ static_var_list0(static_var_list) = ilagged;
+ static_var_list(static_var_list) = ~ilagged;
+ aaa(static_var_list,k) = AS(~ilagged,:)*alphahat(:,k)+BS(~ilagged,:)*etahat(:,k)+CS(~ilagged);
+ if any(ilagged) && k>1
+ aaa(static_var_list0,k) = Tstar(ilagged,:)*alphahat(:,k-1)+Rstar(ilagged,:)*etahat(:,k)+Cstar(ilagged);
+ end
+
end
alphahat=aaa;
- aaa=zeros(M_.endo_nbr,gend);
+
+ % reconstruct updated variables
bbb=zeros(M_.endo_nbr,gend);
- bbb(oo_.dr.restrict_var_list,:)=ahat;
- aaa(oo_.dr.restrict_var_list,:)=aahat;
- for k=d+2:gend
+ bbb(oo_.dr.restrict_var_list,:)=ahat; % this is t|t
+ for k=1:gend
if isoccbin
A = TT(:,:,k);
B = RR(:,:,k);
C = CC(:,k);
- bbb(:,k) = C+A*aaa(:,k-1)+B*eehat(:,k);
- else
+ iT = iTx(:,:,k);
+ Tstar = A(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),oo_.dr.restrict_var_list);
+ Rstar = B(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list),:);
+ Cstar = C(~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list));
+ AS = Tstar*iT;
+ BS = Rstar-AS*RRx(:,:,k);
+ CS = Cstar-AS*CCx(:,k);
+ static_var_list = ~ismember(1:M_.endo_nbr,oo_.dr.restrict_var_list);
+ ilagged = any(abs(AS*TTx(:,:,k)-Tstar)'>1.e-12);
+ static_var_list0 = static_var_list;
+ static_var_list0(static_var_list) = ilagged;
+ static_var_list(static_var_list) = ~ilagged;
+ bbb(static_var_list,k) = AS(~ilagged,:)*ahat(:,k)+BS(~ilagged,:)*eehat(:,k)+CS(~ilagged);
+ if any(ilagged) && k>d+1
+ bbb(static_var_list0,k) = Tstar(ilagged,:)*aahat(:,k-1)+Rstar(ilagged,:)*eehat(:,k)+Cstar(ilagged);
+ end
+ elseif k>d+1
opts_simul.curb_retrench = options_.occbin.smoother.curb_retrench;
opts_simul.waitbar = options_.occbin.smoother.waitbar;
opts_simul.maxit = options_.occbin.smoother.maxit;
@@ -430,11 +465,9 @@ else
% period ahead (so if regimestart was [1 5] it should be [1 4]
% in out
% end
- bbb(oo_.dr.inv_order_var,k) = out.piecewise(1,:);
+ bbb(oo_.dr.inv_order_var,k) = out.piecewise(1,:) - out.ys';
end
end
- % do not overwrite accurate computations using reduced st. space
- bbb(oo_.dr.restrict_var_list,:)=ahat;
ahat0=ahat;
ahat=bbb;
if ~isempty(P)
@@ -475,7 +508,7 @@ else
% in out
% end
for jnk=1:options_.nk
- aaa(jnk,oo_.dr.inv_order_var,k+jnk-1) = out.piecewise(jnk,:);
+ aaa(jnk,oo_.dr.inv_order_var,k+jnk-1) = out.piecewise(jnk,:) - out.ys';
end
end
aK=aaa;
diff --git a/matlab/compute_Pinf_Pstar.m b/matlab/compute_Pinf_Pstar.m
index 9158f27b5bf5f4c7bc2d886e7e9814e6356b77f6..9a05188a5e6b385377bae8812d787d676eb570df 100644
--- a/matlab/compute_Pinf_Pstar.m
+++ b/matlab/compute_Pinf_Pstar.m
@@ -48,6 +48,16 @@ function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_colum
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
np = size(T,1);
+if iszero(Q)
+ % this may happen if users set Q=0 and use heteroskedastic shocks to set
+ % variances period by period
+ % this in practice triggers a form of conditional filter where states
+ % are initialized at st. state with zero variances
+ Pstar=T*0;
+ Pinf=T*0;
+ return
+end
+
if nargin == 6
indx = restrict_columns;
indx0=find(~ismember([1:np],indx));
diff --git a/matlab/expand_group.m b/matlab/expand_group.m
index 9c2e087944dabea2eef6e8ba58a3480d6973f90c..8e54f995a35a8a5c9851b959e1b3ee68684402ca 100644
--- a/matlab/expand_group.m
+++ b/matlab/expand_group.m
@@ -46,6 +46,9 @@ options.nobs=mydata.nobs;
label = mydata.shock_group.label;
label = strrep(label,' ','_');
label = strrep(label,'-','_');
+label = strrep(label,'(','');
+label = strrep(label,')','');
+label = strrep(label,'.','');
shocks = mydata.shock_group.shocks;
options.plot_shock_decomp.fig_name = [mydata.fig_name '. Expand'];
options.plot_shock_decomp.use_shock_groups = label; %[use_shock_groups_old int2str(ic)];
diff --git a/matlab/marginal_density.m b/matlab/marginal_density.m
index 489fb8cf61e5b05770886326bfc38b67a3bc2ccb..8bca13e6681c7204709d7dcc4d12e4983356b9f8 100644
--- a/matlab/marginal_density.m
+++ b/matlab/marginal_density.m
@@ -69,7 +69,14 @@ parameter_names = bayestopt_.name;
save([M_.dname filesep 'Output' filesep M_.fname '_mean.mat'],'xparam1','hh','parameter_names','SIGMA');
fprintf('Estimation::marginal density: I''m computing the posterior log marginal density (modified harmonic mean)... ');
-logdetSIGMA = log(det(SIGMA));
+try
+ % use this robust option to avoid inf/nan
+ logdetSIGMA = 2*sum(log(diag(chol(SIGMA))));
+catch
+ % in case SIGMA is not positive definite
+ logdetSIGMA = nan;
+ fprintf('Estimation::marginal density: the covariance of MCMC draws is not positive definite. You may have too few MCMC draws.');
+end
invSIGMA = hh;
marginal = zeros(9,2);
linee = 0;
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index aa06d67dc76187cbad7826bb38cb45e026506a0e..6ba116349471c14480f39fe4375d9a36e07bfd53 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -433,7 +433,7 @@ while notsteady && t<smpl
Pf = P(:,:,t+1);
end
aK(1,:,t+1) = a1(:,t+1);
- if ~isempty(nk) && nk>1 && isoccbin
+ if ~isempty(nk) && nk>1 && isoccbin && (t>=first_period_occbin_update || isinf(first_period_occbin_update))
opts_simul = occbin_options.opts_regime;
opts_simul.SHOCKS = zeros(nk,M_.exo_nbr);
if smoother_redux
@@ -455,11 +455,11 @@ while notsteady && t<smpl
end
PK(jnk,:,:,t+jnk) = Pf;
end
- if isoccbin
+ if isoccbin && (t>=first_period_occbin_update || isinf(first_period_occbin_update))
if smoother_redux
- aK(jnk,:,t+jnk) = out.piecewise(jnk,oo_.dr.order_var(oo_.dr.restrict_var_list));
+ aK(jnk,:,t+jnk) = out.piecewise(jnk,oo_.dr.order_var(oo_.dr.restrict_var_list)) - out.ys(oo_.dr.order_var(oo_.dr.restrict_var_list))';
else
- aK(jnk,oo_.dr.inv_order_var,t+jnk) = out.piecewise(jnk,:);
+ aK(jnk,oo_.dr.inv_order_var,t+jnk) = out.piecewise(jnk,:) - out.ys';
end
elseif jnk>1
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
@@ -486,6 +486,9 @@ varargout{1} = regimes_;
varargout{2} = TTT;
varargout{3} = RRR;
varargout{4} = CCC;
+varargout{5} = TT;
+varargout{6} = RR;
+varargout{7} = CC;
% $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
% $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
% $$$ Fi_s = Fi(:,t);