stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'occbin_bugfixes_marco' into 'master' 

Improvements to OccBin filtering/smoothing

See merge request Dynare/dynare!2237
covariance-quadratic-approximation
Sébastien Villemot 2023-12-18 18:55:12 +00:00
parent c363355dbf d7ae0d0039
commit 70b9d9277a
16 changed files with 948 additions and 337 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
matlab/+occbin/DSGE_smoother.m
View File

@ -65,8 +65,10 @@ regime_history=[];
if options_.occbin.smoother.linear_smoother && nargin==12
%% linear smoother
options_.occbin.smoother.status=false;
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
tmp_smoother=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty);
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = ...
DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
tmp_smoother=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,...
aK,P,PK,decomp,Trend,state_uncertainty);
for jf=1:length(smoother_field_list)
oo_.occbin.linear_smoother.(smoother_field_list{jf}) = tmp_smoother.(smoother_field_list{jf});
end
@ -80,7 +82,9 @@ if options_.occbin.smoother.linear_smoother && nargin==12
oo_.occbin.linear_smoother.T0=T0;
oo_.occbin.linear_smoother.R0=R0;
oo_.occbin.linear_smoother.decomp=decomp;
oo_.occbin.linear_smoother.alphahat0=alphahat0;
oo_.occbin.linear_smoother.state_uncertainty0=state_uncertainty0;
fprintf('\nOccbin: linear smoother done.\n')
options_.occbin.smoother.status=true;
end
@ -115,10 +119,43 @@ opts_simul.piecewise_only = options_.occbin.smoother.piecewise_only;
occbin_options = struct();
occbin_options.first_period_occbin_update = options_.occbin.smoother.first_period_occbin_update;
occbin_options.opts_regime = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.likelihood.init_regime_history;
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
occbin_options.opts_simul = opts_simul; % this builds the opts_simul options field needed by occbin.solver
occbin_options.opts_regime.binding_indicator = options_.occbin.smoother.init_binding_indicator;
occbin_options.opts_regime.regime_history=options_.occbin.smoother.init_regime_history;
error_indicator=false;
try
%blanket try-catch should be replaced be proper error handling, see https://git.dynare.org/Dynare/dynare/-/merge_requests/2226#note_20318
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);% T1=TT;
catch ME
error_indicator=true;
disp(ME.message)
for iter = 1:numel(ME.stack)
ME.stack(iter)
end
end
if error_indicator || isempty(alphahat0)
etahat= oo_.occbin.linear_smoother.etahat;
alphahat0= oo_.occbin.linear_smoother.alphahat0;
base_regime = struct();
if M_.occbin.constraint_nbr==1
base_regime.regime = 0;
base_regime.regimestart = 1;
else
base_regime.regime1 = 0;
base_regime.regimestart1 = 1;
base_regime.regime2 = 0;
base_regime.regimestart2 = 1;
end
oo_.occbin.smoother.regime_history = [];
for jper=1:size(alphahat,2)+1
if jper == 1
oo_.occbin.smoother.regime_history = base_regime;
else
oo_.occbin.smoother.regime_history(jper) = base_regime;
end
end
end
oo_.occbin.smoother.realtime_regime_history = oo_.occbin.smoother.regime_history;
regime_history = oo_.occbin.smoother.regime_history;
@ -139,6 +176,7 @@ opts_simul.SHOCKS = [etahat(:,1:end)'; zeros(1,M_.exo_nbr)];
opts_simul.exo_pos = 1:M_.exo_nbr;
opts_simul.endo_init = alphahat0(oo_.dr.inv_order_var,1);
opts_simul.init_regime=regime_history; % use realtime regime for guess, to avoid multiple solution issues!
opts_simul.periods = size(opts_simul.SHOCKS,1);
options_.occbin.simul=opts_simul;
options_.noprint = true;
[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
@ -261,13 +299,13 @@ while is_changed && maxiter>iter && ~is_periodic
eee(:,k) = eig(TT(:,:,k));
end
if options_.debug
err_eig(iter-1) = max(max(abs(sort(eee)-sort(sto_eee))));
err_alphahat(iter-1) = max(max(max(abs(alphahat-sto_alphahat))));
err_etahat(iter-1) = max(max(max(abs(etahat-sto_etahat{iter-1}))));
err_CC(iter-1) = max(max(max(abs(CC-sto_CC))));
err_RR(iter-1) = max(max(max(abs(RR-sto_RR))));
err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
end
err_eig(iter-1) = max(max(abs(sort(eee)-sort(sto_eee))));
err_alphahat(iter-1) = max(max(max(abs(alphahat-sto_alphahat))));
err_etahat(iter-1) = max(max(max(abs(etahat-sto_etahat{iter-1}))));
err_CC(iter-1) = max(max(max(abs(CC-sto_CC))));
err_RR(iter-1) = max(max(max(abs(RR-sto_RR))));
err_TT(iter-1) = max(max(max(abs(TT-sto_TT))));
end
end
if occbin_smoother_debug || is_periodic
@ -391,6 +429,10 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
oo_.occbin.smoother.T0=TT;
oo_.occbin.smoother.R0=RR;
oo_.occbin.smoother.C0=CC;
oo_.occbin.smoother.simul.piecewise = out.piecewise(1:end-1,:);
if ~options_.occbin.simul.piecewise_only
oo_.occbin.smoother.simul.linear = out.linear(1:end-1,:);
end
if options_.occbin.smoother.plot
GraphDirectoryName = CheckPath('graphs',M_.fname);
latexFolder = CheckPath('latex',M_.dname);
@ -404,7 +446,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
j1=0;
ifig=0;
for j=1:M_.exo_nbr
if M_.Sigma_e(j,j)
if max(abs(oo_.occbin.smoother.etahat(j,:)))>1.e-8
j1=j1+1;
if mod(j1,9)==1
hh_fig = dyn_figure(options_.nodisplay,'name','Occbin smoothed shocks');
@ -441,15 +483,15 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
fprintf(fidTeX,'\\label{Fig:smoothedshocks_occbin:%s}\n',int2str(ifig));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
end
end
end
if mod(j1,9)~=0 && j==M_.exo_nbr
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
dyn_saveas(hh_fig,[GraphDirectoryName filesep M_.fname,'_smoothedshocks_occbin',int2str(ifig)],options_.nodisplay,options_.graph_format);
if mod(j1,9)~=0 && j==M_.exo_nbr
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
dyn_saveas(hh_fig,[GraphDirectoryName filesep M_.fname,'_smoothedshocks_occbin',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
% TeX eps loader file
fprintf(fidTeX,'\\begin{figure}[H]\n');
@ -463,6 +505,6 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fclose(fidTeX);
end
end
end
end

87
matlab/+occbin/backward_map_regime.m Normal file
View File

@ -0,0 +1,87 @@
function [binding_indicator, A, regime_string] = backward_map_regime(regime, regime_start)
% [binding_indicator, A, regime_string] = backward_map_regime(regime, regime_start)
% Map regime information into regime indicator
%
% Inputs:
% - regime [integer] [1 by n_transitions] vector of regime number indices
% - regime_start [integer] [1 by n_transitions] vectors with period numbers in which regime starts
%
% Outputs:
% - binding_indicator [integer] [nperiods by 1] vector of regime indices
% - A [bin] binary representation of binding indicator
% - error_flag [boolean] 1 if regime never leaves 1 or is still there at the end of nperiods
% 0 otherwise
% Copyright © 2023 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
if nargin ==1
% polymorphism
if ~isstruct(regime)
disp('error::backward_map_regime')
disp('input arguments may be 1 structure with regime info')
disp('or two arrays: regime and regimestart ')
error('wrong input')
end
fnam = fieldnames(regime);
if length(fnam) == 2
[binding_indicator, A, regime_string] = occbin.backward_map_regime(regime.regime, regime.regimestart);
else
for k=1:2
nperiods(k) = regime.(['regimestart' int2str(k)])(end);
number_of_binary_tokens(k) = ceil((nperiods(k)-1)/50);
end
binding_indicator = false(max(nperiods),2);
A = int64(zeros(max(number_of_binary_tokens),2));
for k=1:2
[binding_indicator(1:nperiods(k),k), A(1:number_of_binary_tokens(k),k), tmp{k}] = ...
occbin.backward_map_regime(regime.(['regime' int2str(k)]), regime.(['regimestart' int2str(k)]));
end
regime_string = char(tmp{1},tmp{2});
end
return
else
if isstruct(regime)
disp('error::backward_map_regime')
disp('input arguments may be ONE structure with regime info')
disp('or TWO arrays: regime and regimestart ')
error('wrong input')
end
end
regime_string = char(mat2str(double(regime)),mat2str(regime_start));
nperiods_0 = regime_start(end);
number_of_binary_tokens = max(1,ceil((nperiods_0-1)/50));
A = int64(zeros(number_of_binary_tokens,1));
binding_indicator = false(nperiods_0,1);
if length(regime)>1
for ir=1:length(regime)-1
binding_indicator(regime_start(ir):regime_start(ir+1)-1,1) = regime(ir);
for k=regime_start(ir):regime_start(ir+1)-1
this_token = ceil(k/50);
A(this_token) = int64(bitset(A(this_token),k-50*(this_token-1),regime(ir)));
end
end
end
binding_indicator = logical(binding_indicator);
% to convert regime in a readable string array
% a = dec2bin(A);
% bindicator = [a(end:-1:1) '0'];

50
matlab/+occbin/cost_function.m Normal file
View File

@ -0,0 +1,50 @@
function [cost, out] = cost_function(err_0, current_obs, weights, opts_simul,...
M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% [cost, out] = cost_function(err_0, current_obs, opts_simul,...
% M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% Outputs:
% - cost [double] penalty
% - out [structure] Occbin's results structure
%
% Inputs
% - err_0 [double] value of shocks
% - current_obs [double] [1 by n_obs] current value of observables
% - weights [double] [1 by n_obs] variance of observables,
% - opts_simul [structure] Structure with simulation options
% used in cost function
% - M_ [structure] Matlab's structure describing the model (M_).
% - dr_ [structure] model information structure
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] Matlab's structure describing the current options (options_).
% Copyright © 2023 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
opts_simul.SHOCKS = err_0';
options_.occbin.simul=opts_simul;
options_.occbin.simul.full_output=1;
options_.noprint = 1;
[~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
cost = 0;
if ~out.error_flag
cost = mean((out.piecewise(1,opts_simul.varobs_id)'-current_obs').^2./weights);
else
cost = cost+1.e10;
end

63
matlab/+occbin/findmin.m Normal file
View File

@ -0,0 +1,63 @@
function [y, out, cost] = findmin(d_index, a0, P1, Qt, Y, ZZ, opts_simul,M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% [y, out, cost] = findmin(d_index, a0, P1, Qt, Y, ZZ, opts_simul,M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_)
% Outputs:
% - cost [double] penalty
% - out [structure] Occbin's results structure
%
% Inputs
% - opts_simul [structure] Structure with simulation options
% used in cost function
% - M_ [structure] Matlab's structure describing the model (M_).
% - dr_ [structure] model information structure
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% - options_ [structure] Matlab's structure describing the current options (options_).
% Copyright © 2023 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
current_obs = Y(d_index,2)'+dr.ys(options_.varobs_id(d_index))';
err_index = find(diag(Qt(:,:,2))~=0);
F = ZZ(d_index,:)*P1(:,:,2)*ZZ(d_index,:)' ;
weights=diag(F);
filtered_errs_init = zeros(1,length(err_index));
opts_simul.varobs_id=options_.varobs_id(d_index)';
opts_simul.exo_pos=err_index; %err_index is predefined mapping from observables to shocks
opts_simul.SHOCKS = filtered_errs_init;
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(dr.restrict_var_list,1)=a0(:,1); %updated state
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = a0(dr.inv_order_var,1);
end
[y] = fminsearch(@cost_function,filtered_errs_init');
[cost, out] = occbin.cost_function(y, current_obs, weights, opts_simul,...
M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_);
function cost = cost_function(x)
cost = occbin.cost_function(x, current_obs, weights, opts_simul,...
M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_);
end
end

233
matlab/+occbin/kalman_update_algo_1.m
View File

@ -1,4 +1,4 @@
function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat, alphahat, V] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
% function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat] = kalman_update_algo_1(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options)
% INPUTS
% - a [N by 1] t-1's state estimate
@ -86,6 +86,8 @@ else
base_regime.regimestart2 = 1;
end
regimes_ = [base_regime base_regime base_regime];
opts_simul = occbin_options.opts_simul;
options_.occbin.simul=opts_simul;
mm=size(a,1);
%% store info in t=1
@ -103,16 +105,46 @@ PZI = P1(:,:,t)*ZZ'*iF(di,di,t);
% L(:,:,t) = T-K(:,di,t)*ZZ;
L(:,:,t) = eye(mm)-PZI*ZZ;
if ~options_.occbin.filter.use_relaxation
[a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
if ~isempty(fieldnames(regimes0))
if options_.occbin.filter.guess_regime
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, regimes0(1),'reduced_state_space',T0,R0);
if M_.occbin.constraint_nbr==1
bindx = occbin.backward_map_regime(regimes0(1).regime, regimes0(1).regimestart);
bindx = bindx(2:end);
[regimes0(2).regime, regimes0(2).regimestart, error_flag]=occbin.map_regime(bindx);
bindx = bindx(2:end);
[regimes0(3).regime, regimes0(3).regimestart, error_flag]=occbin.map_regime(bindx);
else
bindx1 = occbin.backward_map_regime(regimes0(1).regime1, regimes0(1).regimestart1);
bindx2 = occbin.backward_map_regime(regimes0(1).regime2, regimes0(1).regimestart2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(2).regime1, regimes0(2).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(2).regime2, regimes0(2).regimestart2, error_flag]=occbin.map_regime(bindx2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(3).regime1, regimes0(3).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(3).regime2, regimes0(3).regimestart2, error_flag]=occbin.map_regime(bindx2);
end
% regimes0=[regimes0 base_regime base_regime];
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
end
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
else
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, base_regime,'reduced_state_space',T0,R0);
regimes0(1)=base_regime;
if isempty(fieldnames(regimes0))
regimes0 = regimes_;
else
regimes0(1)=base_regime;
end
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
[a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
end
if error_flag
etahat=NaN(size(QQQ,1),1);
@ -120,7 +152,6 @@ if error_flag
end
%% run here the occbin simul
opts_simul = occbin_options.opts_simul;
opts_simul.SHOCKS = zeros(3,M_.exo_nbr);
opts_simul.exo_pos=1:M_.exo_nbr;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
@ -134,14 +165,24 @@ else
my_order_var = dr.order_var;
end
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
if options_.occbin.filter.guess_regime
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
else
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
end
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
return;
end
@ -159,53 +200,18 @@ end
lik_hist=lik;
niter=1;
is_periodic=0;
if options_.occbin.filter.use_relaxation || isequal(regimes0(1),base_regime)
nguess=1;
else
nguess=0;
end
newguess=0;
if any(myregime) || ~isequal(regimes_(1),regimes0(1))
while ~isequal(regimes_(1),regimes0(1)) && ~is_periodic && ~out.error_flag && niter<=options_.occbin.likelihood.max_number_of_iterations
niter=niter+1;
oldstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes0(1).regimestart)>1
oldstart = regimes0(1).regimestart(end);
end
newstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes_(1).regimestart)>1
newstart = regimes_(1).regimestart(end);
end
if M_.occbin.constraint_nbr==1 && (newstart-oldstart)>2 && options_.occbin.filter.use_relaxation
regimestart = max(oldstart+2,round(0.5*(newstart+oldstart)));
regimestart = min(regimestart,oldstart+4);
if regimestart<=regimes_(1).regimestart(end-1)
if length(regimes_(1).regimestart)<=3
regimestart = max(regimestart, min(regimes_(1).regimestart(end-1)+2,newstart));
else
regimes_(1).regime = regimes_(1).regime(1:end-2);
regimes_(1).regimestart = regimes_(1).regimestart(1:end-2);
regimestart = max(regimestart, regimes_(1).regimestart(end-1)+1);
end
end
regimes_(1).regimestart(end)=regimestart;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space', T0, R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
elseif newguess==0
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
end
newguess=0;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
regime_hist(niter) = {regimes_(1)};
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
[a, a1, P, P1, v, alphahat, etahat, lik] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
[a, a1, P, P1, v, alphahat, etahat, lik, V] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
@ -216,9 +222,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
else
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
if not(options_.occbin.filter.use_relaxation)
opts_simul.init_regime=regimes_(1);
end
opts_simul.init_regime=regimes_(1);
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
@ -227,9 +231,14 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=[];
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
return;
end
regimes0=regimes_;
@ -238,36 +247,14 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
for kiter=1:niter-1
is_periodic(kiter) = isequal(regime_hist{kiter}, regimes_(1));
end
is_periodic_iter = find(is_periodic);
is_periodic = any(is_periodic);
if is_periodic
if nguess<3 && M_.occbin.constraint_nbr==1
newguess=1;
is_periodic=0;
nguess=nguess+1;
if nguess==1
% change starting regime
regimes_(1).regime=0;
regimes_(1).regimestart=1;
elseif nguess==2
% change starting regime
regimes_(1).regime=[0 1 0];
regimes_(1).regimestart=[1 2 3];
else
regimes_(1).regime=[1 0];
regimes_(1).regimestart=[1 2];
end
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
regime_hist = regime_hist(1);
niter=1;
else
% re-set to previous regime
regimes_ = regimes0;
% force projection conditional on previous regime
opts_simul.init_regime=regimes0(1);
% re-set to previous regime
if options_.occbin.filter.periodic_solution
% force projection conditional on most likely regime
[m, im]=min(lik_hist(is_periodic_iter:end));
opts_simul.init_regime=regime_hist{is_periodic_iter+im-1};
if M_.occbin.constraint_nbr==1
myregimestart = [regimes0.regimestart];
else
@ -280,8 +267,20 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
return;
else
regimes_ = out.regime_history;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
[a, a1, P, P1, v, alphahat, etahat, lik, V] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood, options_.occbin.filter.state_covariance);
end
else
error_flag = 330;
etahat=etahat(:,2);
lik=inf;
return;
end
end
end
@ -291,61 +290,6 @@ end
error_flag = out.error_flag;
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1))
error_flag = 1;
if M_.occbin.constraint_nbr==1 % try some other regime
[~, il]=sort(regime_end);
rr=regime_hist(il(2:3));
newstart=1;
if length(rr{1}.regimestart)>1
newstart = rr{1}.regimestart(end)-rr{1}.regimestart(end-1)+1;
end
oldstart=1;
if length(rr{2}.regimestart)>1
oldstart = rr{2}.regimestart(end)-rr{2}.regimestart(end-1)+1;
end
nstart=sort([newstart oldstart]);
regimes_=rr{1}(1);
for k=(nstart(1)+1):(nstart(2)-1)
niter=niter+1;
regimes_(1).regimestart(end)=k;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state, [base_regime regimes_(1)],'reduced_state_space',T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, alphahat, etahat, lik] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, options_.rescale_prediction_error_covariance, options_.occbin.likelihood.IF_likelihood);
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
regime_hist(niter) = {regimes_(1)};
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
% opts_simul.init_regime=regimes_(1);
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
myregimestart = [regimes_.regimestart1 regimes_.regimestart2];
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
if isequal(out.regime_history(1),regimes_(1))
error_flag=0;
break
end
end
regimes_ = out.regime_history;
end
end
if ~error_flag
@ -358,17 +302,26 @@ C = ss.C(my_order_var,1:2);
QQ = R(:,:,2)*QQQ(:,:,3)*transpose(R(:,:,2));
P(:,:,1) = P(:,:,2);
P(:,:,2) = T(:,:,2)*P(:,:,1)*transpose(T(:,:,2))+QQ;
etahat=etahat(:,2);
if nargout>=13
etahat=etahat(:,2);
end
warning_config;
end
function [a, a1, P, P1, v, alphahat, etahat, lik, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, rescale_prediction_error_covariance, IF_likelihood)
function [a, a1, P, P1, v, alphahat, etahat, lik, V, error_flag] = occbin_kalman_update0(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,iF,L,mm, rescale_prediction_error_covariance, IF_likelihood, state_uncertainty_flag)
alphahat=NaN(size(a));
etahat=NaN(size(QQQ,1),2);
lik=Inf;
error_flag=0;
if state_uncertainty_flag
V = zeros(mm,mm,2);
N = zeros(mm,mm,3);
else
V=[];
end
warning off
if nargin<18
IF_likelihood=0;
@ -444,9 +397,15 @@ while t > 1
if isempty(di)
% in this case, L is simply T due to Z=0, so that DK (2012), eq. 4.93 obtains
r(:,t) = L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38 with Z=0
if state_uncertainty_flag
N(:,:,t)=L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42 with Z=0
end
else
ZZ = Z(di,:);
r(:,t) = ZZ'*iF(di,di,t)*v(di,t) + L(:,:,t)'*r(:,t+1); %compute r_{t-1}, DK (2012), eq. 4.38
if state_uncertainty_flag
N(:,:,t)=ZZ'*iF(di,di,t)*ZZ+L(:,:,t)'*N(:,:,t+1)*L(:,:,t); %compute N_{t-1}, DK (2012), eq. 4.42
end
end
Q=QQQ(:,:,t);
QRt = Q*transpose(RR(:,:,t));
@ -454,6 +413,10 @@ while t > 1
alphahat(:,t) = a1(:,t) + P1(:,:,t)*r(:,t); %DK (2012), eq. 4.35
etahat(:,t) = QRt*r(:,t); %DK (2012), eq. 4.63
r(:,t) = T'*r(:,t); % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
if state_uncertainty_flag
V(:,:,t) = P1(:,:,t)-P1(:,:,t)*N(:,:,t)*P1(:,:,t); %DK (2012), eq. 4.43
N(:,:,t) = T'*N(:,:,t)*T; %DK (2012), eq. 4.43
end
if IF_likelihood && t==2 && not(isempty(di))
ishocks = any(ZZ*RR(:,:,t));

247
matlab/+occbin/kalman_update_algo_3.m
View File

@ -1,4 +1,4 @@
function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, lik, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
% function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, alphahat, etahat, TT, RR, CC] = kalman_update_algo_3(a,a1,P,P1,data_index,Z,v,Fi,Ki,Y,H,QQQ,T0,R0,TT,RR,CC,regimes0,M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,occbin_options,kalman_tol,nk)
%
% INPUTS
@ -85,7 +85,8 @@ if isempty(nk)
end
nk=max(nk,1);
opts_simul = occbin_options.opts_regime;
opts_simul = occbin_options.opts_simul;
options_.occbin.simul=opts_simul;
base_regime = struct();
if M_.occbin.constraint_nbr==1
base_regime.regime = 0;
@ -96,6 +97,7 @@ else
base_regime.regime2 = 0;
base_regime.regimestart2 = 1;
end
regimes_ = [base_regime base_regime base_regime];
myrestrict=[];
if options_.smoother_redux
opts_simul.restrict_state_space =1;
@ -103,15 +105,46 @@ if options_.smoother_redux
end
mm=size(a,1);
if ~options_.occbin.filter.use_relaxation
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
if ~isempty(fieldnames(regimes0))
if options_.occbin.filter.guess_regime
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, regimes0(1),myrestrict,T0,R0);
if M_.occbin.constraint_nbr==1
bindx = occbin.backward_map_regime(regimes0(1).regime, regimes0(1).regimestart);
bindx = bindx(2:end);
[regimes0(2).regime, regimes0(2).regimestart, error_flag]=occbin.map_regime(bindx);
bindx = bindx(2:end);
[regimes0(3).regime, regimes0(3).regimestart, error_flag]=occbin.map_regime(bindx);
else
bindx1 = occbin.backward_map_regime(regimes0(1).regime1, regimes0(1).regimestart1);
bindx2 = occbin.backward_map_regime(regimes0(1).regime2, regimes0(1).regimestart2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(2).regime1, regimes0(2).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(2).regime2, regimes0(2).regimestart2, error_flag]=occbin.map_regime(bindx2);
bindx1 = bindx1(2:end);
bindx2 = bindx2(2:end);
[regimes0(3).regime1, regimes0(3).regimestart1, error_flag]=occbin.map_regime(bindx1);
[regimes0(3).regime2, regimes0(3).regimestart2, error_flag]=occbin.map_regime(bindx2);
end
% regimes0=[regimes0 base_regime base_regime];
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
end
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
else
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, base_regime,myrestrict,T0,R0);
if isempty(fieldnames(regimes0))
regimes0 = regimes_;
else
regimes0(1)=base_regime;
end
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
regimes0(1)=base_regime;
end
@ -132,9 +165,23 @@ else
end
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if options_.occbin.filter.guess_regime
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
else
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
options_.occbin.simul.init_regime=regimes0;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
end
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
return;
end
@ -148,65 +195,23 @@ regime_hist = {regimes0(1)};
if M_.occbin.constraint_nbr==1
regime_end = regimes0(1).regimestart(end);
end
lik_hist=lik;
niter=1;
is_periodic=0;
if options_.occbin.filter.use_relaxation || isequal(regimes0(1),base_regime)
nguess=1;
else
nguess=0;
end
newguess=0;
if any(myregime) || ~isequal(regimes_(1),regimes0(1))
while ~isequal(regimes_(1),regimes0(1)) && ~is_periodic && ~out.error_flag && niter<=options_.occbin.likelihood.max_number_of_iterations
niter=niter+1;
newstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes_(1).regimestart)>1
newstart = regimes_(1).regimestart(end);
end
oldstart=1;
if M_.occbin.constraint_nbr==1 && length(regimes0(1).regimestart)>1
oldstart = regimes0(1).regimestart(end);
end
if M_.occbin.constraint_nbr==1 && (newstart-oldstart)>2 && options_.occbin.filter.use_relaxation
regimestart = max(oldstart+2,round(0.5*(newstart+oldstart)));
regimestart = min(regimestart,oldstart+4);
if regimestart<=regimes_(1).regimestart(end-1)
if length(regimes_(1).regimestart)<=3
regimestart = max(regimestart, min(regimes_(1).regimestart(end-1)+2,newstart));
else
regimes_(1).regime = regimes_(1).regime(1:end-2);
regimes_(1).regimestart = regimes_(1).regimestart(1:end-2);
regimestart = max(regimestart, regimes_(1).regimestart(end-1)+1);
end
end
% % if (newstart-oldstart)>3
% % regimestart = regimes_(1).regimestart(end-1)+oldstart+2;
% % % regimestart = regimes_(1).regimestart(end-1)+round(0.5*(newstart+oldstart))-1;
regimes_(1).regimestart(end)=regimestart;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
elseif newguess==0
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
end
newguess=0;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
regime_hist(niter) = {regimes_(1)};
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
% if opts_simul.restrict_state_space
% tmp=zeros(M_.endo_nbr,1);
% tmp(dr.restrict_var_list,1)=alphahat(:,1);
% opts_simul.endo_init = tmp(dr.inv_order_var,1);
% else
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
@ -214,10 +219,7 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
else
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
% end
if not(options_.occbin.filter.use_relaxation)
opts_simul.init_regime=regimes_(1);
end
opts_simul.init_regime=regimes_(1);
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
@ -228,6 +230,8 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
return;
end
regimes0=regimes_;
@ -236,36 +240,14 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
for kiter=1:niter-1
is_periodic(kiter) = isequal(regime_hist{kiter}, regimes_(1));
end
is_periodic_iter = find(is_periodic);
is_periodic = any(is_periodic);
if is_periodic
if nguess<3 && M_.occbin.constraint_nbr==1
newguess=1;
is_periodic=0;
nguess=nguess+1;
if nguess==1
% change starting regime
regimes_(1).regime=0;
regimes_(1).regimestart=1;
elseif nguess==2
% change starting regime
regimes_(1).regime=[0 1 0];
regimes_(1).regimestart=[1 2 3];
else
regimes_(1).regime=[1 0];
regimes_(1).regimestart=[1 2];
end
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
regime_hist = regime_hist(1);
niter=1;
else
% re-set to previous regime
regimes_ = regimes0;
% force projection conditional on previous regime
opts_simul.init_regime=regimes0(1);
if options_.occbin.filter.periodic_solution
% force projection conditional on most likely regime
[m, im]=min(lik_hist(is_periodic_iter:end));
opts_simul.init_regime=regime_hist{is_periodic_iter+im-1};
if M_.occbin.constraint_nbr==1
myregimestart = [regimes0.regimestart];
else
@ -274,68 +256,39 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
opts_simul.maxit=1;
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
end
[~, out, ss] = occbin.solver(M_,oo_,options_);
if out.error_flag
error_flag = out.error_flag;
etahat=etahat(:,2);
lik=inf;
return;
else
regimes_ = out.regime_history;
TT(:,:,2)=ss.T(my_order_var,my_order_var,1);
RR(:,:,2)=ss.R(my_order_var,:,1);
CC(:,2)=ss.C(my_order_var,1);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
end
else
error_flag = 330;
etahat=etahat(:,2);
lik=inf;
return;
end
end
end
end
end
error_flag = out.error_flag;
if error_flag==0 && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1)) %fixed point algorithm did not converge
if ~error_flag && niter>options_.occbin.likelihood.max_number_of_iterations && ~isequal(regimes_(1),regimes0(1)) %fixed point algorithm did not converge
error_flag = 1;
if M_.occbin.constraint_nbr==1
% try some other regime before giving up
[~, il]=sort(regime_end);
rr=regime_hist(il(2:3));
newstart=1;
if length(rr{1}(1).regimestart)>1
newstart = rr{1}(1).regimestart(end)-rr{1}(1).regimestart(end-1)+1;
end
oldstart=1;
if length(rr{2}(1).regimestart)>1
oldstart = rr{2}(1).regimestart(end)-rr{2}(1).regimestart(end-1)+1;
end
nstart=sort([newstart oldstart]);
regimes_=rr{1}(1);
for k=(nstart(1)+1):(nstart(2)-1)
niter=niter+1;
regimes_(1).regimestart(end)=k;
[~,~,~,~,~,~, TTx, RRx, CCx] ...
= occbin.dynare_resolve(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state, [base_regime regimes_(1)],myrestrict,T0,R0);
TT(:,:,2) = TTx(:,:,end);
RR(:,:,2) = RRx(:,:,end);
CC(:,2) = CCx(:,end);
[a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol);
opts_simul.SHOCKS(1,:) = etahat(:,2)';
if opts_simul.restrict_state_space
tmp=zeros(M_.endo_nbr,1);
tmp(dr.restrict_var_list,1)=alphahat(:,1);
opts_simul.endo_init = tmp(dr.inv_order_var,1);
else
opts_simul.endo_init = alphahat(dr.inv_order_var,1);
end
if M_.occbin.constraint_nbr==1
myregimestart = [regimes_.regimestart];
else
myregimestart = [regimes_.regimestart1 regimes_.regimestart2];
end
opts_simul.periods = max(opts_simul.periods,max(myregimestart));
options_.occbin.simul=opts_simul;
[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
if isequal(out.regime_history(1),regimes_(1))
error_flag=0;
break
end
end
regimes_ = out.regime_history;
end
end
regimes_=regimes_(1:3);
a = out.piecewise(1:nk+1,my_order_var)' - repmat(out.ys(my_order_var),1,nk+1);
if ~error_flag
a = out.piecewise(1:nk+1,my_order_var)' - repmat(out.ys(my_order_var),1,nk+1);
regimes_=regimes_(1:3);
end
T = ss.T(my_order_var,my_order_var,:);
R = ss.R(my_order_var,:,:);
C = ss.C(my_order_var,:);
@ -352,7 +305,7 @@ end
end
function [a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol)
function [a, a1, P, P1, v, Fi, Ki, alphahat, etahat, lik] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol)
% [a, a1, P, P1, v, Fi, Ki, alphahat, etahat] = occbin_kalman_update(a,a1,P,P1,data_index,Z,v,Y,H,QQQ,TT,RR,CC,Ki,Fi,mm,kalman_tol)
% - a
% - a1
@ -377,9 +330,18 @@ a(:,t) = a1(:,t);
P1(:,:,t) = T*P(:,:,t-1)*T' + QQ; %transition according to (6.14) in DK (2012)
P(:,:,t) = P1(:,:,t);
di = data_index{t}';
% store info for lik
if not(isempty(di))
vv = Y(di,t) - Z(di,:)*a(:,t);
F = Z(di,:)*P(:,:,t)*Z(di,:)' + diag(H(di));
sig=sqrt(diag(F));
lik=inf;
end
if isempty(di)
Fi(:,t) = 0;
Ki(:,:,t) = 0;
lik =0;
end
for i=di
Zi = Z(i,:);
@ -394,6 +356,11 @@ for i=di
% p. 157, DK (2012)
end
end
if not(isempty(di))
log_dF = log(det(F./(sig*sig')))+2*sum(log(sig));
iF = inv(F./(sig*sig'))./(sig*sig');
lik = log_dF + transpose(vv)*iF*vv + length(di)*log(2*pi);
end
%% do backward pass
ri=zeros(mm,1);
@ -416,4 +383,4 @@ while t > 1
ri = T'*ri; % KD (2003), eq. (23), equation for r_{t-1,p_{t-1}}
end
end
end

331
matlab/+occbin/kalman_update_engine.m Normal file
View File

@ -0,0 +1,331 @@
function [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V, Fix, Kix, TTx,RRx,CCx] = ...
kalman_update_engine(a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,...
dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
% [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V, Fix, Kix, TTx,RRx,CCx] = kalman_update_engine(
% a0,a1,P0,P1,t,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_,base_regime,d_index,M_,
% dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options, Fi,Ki,kalman_tol,nk)
% Copyright © 2023 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
is_multivariate = true;
Fix=[];
Kix=[];
TTx=[];
RRx=[];
CCx=[];
V=[];
if nargin>26
is_multivariate = false;
end
use_relaxation = false;
if is_multivariate
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, regx, info, M_, likx, alphahat, etahat,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,struct(),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
info0=info;
if info
if ~isequal(regimes_(1:2),[base_regime base_regime])
if is_multivariate
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx, etahat, alphahat, V] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, regx, info, M_, likx, alphahat, etahat,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
end
info1=info;
else
if ~isequal(regimes_(1:2),[base_regime base_regime])
if is_multivariate
[ax1, a1x1, Px1, P1x1, vx1, Tx1, Rx1, Cx1, regx1, info1, M_1, likx1, etahat1, alphahat1, V1] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax1, a1x1, Px1, P1x1, vx1, Fix1, Kix1, Tx1, Rx1, Cx1, regx1, info1, M_1, likx1, alphahat1, etahat1,TTx1,RRx1,CCx1] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,regimes_(1:2),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info1==0 && likx1<likx
ax=ax1;
a1x=a1x1;
Px=Px1;
P1x=P1x1;
vx=vx1;
Tx=Tx1;
Rx=Rx1;
Cx=Cx1;
regx=regx1;
info=info1;
M_= M_1;
likx=likx1;
etahat=etahat1;
alphahat=alphahat1;
if is_multivariate
V=V1;
else
Fix = Fix1;
Kix = Kix1;
TTx = TTx1;
RRx = RRx1;
CCx = CCx1;
end
end
else
if t>options_.occbin.likelihood.number_of_initial_periods_with_extra_regime_guess
info1=0;
else
% may help in first 2 periods to try some other guess regime, due to
% larger state uncertainty
info1=1;
options_.occbin.likelihood.brute_force_regime_guess = true;
options_.occbin.likelihood.loss_function_regime_guess = true;
end
end
end
diffstart=0;
if info==0
if M_.occbin.constraint_nbr==1
oldstart = regimes_(1).regimestart(end);
newstart = regx(1).regimestart(end);
diffstart = newstart-oldstart;
else
newstart1 = regx(1).regimestart1(end);
newstart2 = regx(1).regimestart2(end);
oldstart1 = regimes_(1).regimestart1(end);
oldstart2 = regimes_(1).regimestart2(end);
diffstart = max(newstart1-oldstart1,newstart2-oldstart2);
end
end
if options_.occbin.filter.use_relaxation && diffstart>2
if info0==0
% make sure we match criteria to enter further solution attempts
info1=1;
end
options_.occbin.likelihood.brute_force_regime_guess = true;
options_.occbin.likelihood.loss_function_regime_guess = true;
use_relaxation = true;
end
if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (info==0 && ~isequal(regx(1),base_regime))
guess_regime = [base_regime base_regime];
options_.occbin.filter.guess_regime = true;
use_index = 0;
if M_.occbin.constraint_nbr==1
for k=1:5
guess_regime(1).regimestart=[1 5 5+4*k];
guess_regime(1).regime=[0 1 0];
if is_multivariate
[ax2{1}, a1x2{1}, Px2{1}, P1x2{1}, vx2{1}, Tx2{1}, Rx2{1}, Cx2{1}, regx2{1}, info2, M_2{1}, likx2{1}, etahat2{1}, alphahat2{1}, V2{1}] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2{1}, a1x2{1}, Px2{1}, P1x2{1}, vx2{1}, Fix2{1}, Kix2{1}, Tx2{1}, Rx2{1}, Cx2{1}, regx2{1}, info2, M_2{1}, likx2{1}, alphahat2{1}, etahat2{1},TTx2{1},RRx2{1},CCx2{1}] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info2==0
use_index= 1;
if not(info==0 && isequal(regx2{1},regx)) && not(use_relaxation && likx2{1}>=likx)
% found a solution, different from previous or
% use_relaxation and likelihood is better
break
end
end
guess_regime(1).regimestart=[1 1+4*k];
guess_regime(1).regime=[1 0];
if is_multivariate
[ax2{2}, a1x2{2}, Px2{2}, P1x2{2}, vx2{2}, Tx2{2}, Rx2{2}, Cx2{2}, regx2{2}, info2, M_2{2}, likx2{2}, etahat2{2}, alphahat2{2}, V2{2}] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2{2}, a1x2{2}, Px2{2}, P1x2{2}, vx2{2}, Fix2{2}, Kix2{2}, Tx2{2}, Rx2{2}, Cx2{2}, regx2{2}, info2, M_2{2}, likx2{2}, alphahat2{2}, etahat2{2},TTx2{2},RRx2{2},CCx2{2}] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info2==0
use_index = 2;
% if use_relaxation and we are here, previous guess did not
% improve solution, so we test for this one
end
if use_index
% in case the second guess does not find a solution!
info2=0;
% a solution was found
break
end
end
end
if M_.occbin.constraint_nbr==2
for jk=0:1 % loop over other regime duration. this loop is shorter for parsimony. one may add an option ...
for k=1:5 % loop over current regime duration
gindex = 0;
for jr=1:2 % loop over current regime 1 or 2
if jr==1
regstart1 = 'regimestart1';
reg1 = 'regime1';
regstart2 = 'regimestart2';
reg2 = 'regime2';
else
regstart1 = 'regimestart2';
reg1 = 'regime2';
regstart2 = 'regimestart1';
reg2 = 'regime1';
end
for kk=1:2 % loop over current regime binding in expectation vs binding in current period
if kk==1
guess_regime(1).(regstart1)=[1 5 5+4*k];
guess_regime(1).(reg1)=[0 1 0];
else
guess_regime(1).(regstart1)=[1 1+4*k];
guess_regime(1).(reg1)=[1 0];
end
for kj=1:1+1*(jk>0)
% loop over other regime slack or binding in current period or binding in
% expectation
if jk==0
% other regime is slack
guess_regime(1).(regstart2) = 1;
guess_regime(1).(reg2) = 0;
else % jk>0
if kj==1
% other regime binding in current period
guess_regime(1).(regstart2)=[1 1+4*jk];
guess_regime(1).(reg2) = [1 0];
else
% other regime binding in expectation
guess_regime(1).(regstart2)=[1 5 5+4*jk];
guess_regime(1).(reg2) = [0 1 0];
end
end
gindex = gindex+1;
if is_multivariate
[ax2{gindex}, a1x2{gindex}, Px2{gindex}, P1x2{gindex}, vx2{gindex}, Tx2{gindex}, Rx2{gindex}, Cx2{gindex}, regx2{gindex}, info2, M_2{gindex}, likx2{gindex}, etahat2{gindex}, alphahat2{gindex}, V2{gindex}] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2{gindex}, a1x2{gindex}, Px2{gindex}, P1x2{gindex}, vx2{gindex}, Fix2{gindex}, Kix2{gindex}, Tx2{gindex}, Rx2{gindex}, Cx2{gindex}, regx2{gindex}, info2, M_2{gindex}, likx2{gindex}, alphahat2{gindex}, etahat2{gindex},TTx2{gindex},RRx2{gindex},CCx2{gindex}] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
if info2==0
use_index= gindex;
if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
% found a solution, different from previous one
% use_relaxation and likelihood improves
break
end
end
end % loop over other regime slack, binding in expectation or binding in current period
if info2==0
if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
% found a solution, different from previous one
% use_relaxation and likelihood improves
break
end
end
end % loop over current regime binding in expectation vs binding in current period
if info2==0
if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
% found a solution, different from previous one
% use_relaxation and likelihood improves
break
end
end
end % loop over current regime 1 or 2
if use_index
info2=0;
break
end
end % loop over current regime duration
if use_index
break
end
end % loop over other regime duration
end % 2 constraints
if info2==0
% so that we DO NOT enter IVF step
info0=0;
info1=0;
end
if info2==0 && likx2{use_index}<likx
ax=ax2{use_index};
a1x=a1x2{use_index};
Px=Px2{use_index};
P1x=P1x2{use_index};
vx=vx2{use_index};
Tx=Tx2{use_index};
Rx=Rx2{use_index};
Cx=Cx2{use_index};
regx=regx2{use_index};
info=info2;
M_= M_2{use_index};
likx=likx2{use_index};
etahat=etahat2{use_index};
alphahat=alphahat2{use_index};
if is_multivariate
V=V2{use_index};
else
Fix = Fix2{use_index};
Kix = Kix2{use_index};
TTx = TTx2{use_index};
RRx = RRx2{use_index};
CCx = CCx2{use_index};
end
end
options_.occbin.filter.guess_regime = false;
end
if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %|| (info==0 && ~isequal(regx(1),base_regime))
[~, out] = occbin.findmin(d_index, a0, P1, Qt, Y, Z, occbin_options.opts_simul,M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_);
if out.error_flag==0
options_.occbin.filter.guess_regime = true;
guess_regime=out.regime_history;
guess_regime = [guess_regime base_regime];
if is_multivariate
[ax2, a1x2, Px2, P1x2, vx2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, etahat2, alphahat2, V2] = occbin.kalman_update_algo_1(a0,a1,P0,P1,data_index,Z,vv,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
else
[ax2, a1x2, Px2, P1x2, vx2, Fix2, Kix2, Tx2, Rx2, Cx2, regx2, info2, M_2, likx2, alphahat2, etahat2,TTx2,RRx2,CCx2] = occbin.kalman_update_algo_3(a0,a1,P0,P1,data_index,Z,vv,Fi,Ki,Y,H,Qt,T0,R0,TT,RR,CC,guess_regime,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
end
options_.occbin.filter.guess_regime = false;
if info2==0 && likx2<likx
ax=ax2;
a1x=a1x2;
Px=Px2;
P1x=P1x2;
vx=vx2;
Tx=Tx2;
Rx=Rx2;
Cx=Cx2;
regx=regx2;
info=info2;
likx=likx2;
M_= M_2;
etahat=etahat2;
alphahat=alphahat2;
if is_multivariate
V=V2;
else
Fix = Fix2;
Kix = Kix2;
TTx = TTx2;
RRx = RRx2;
CCx = CCx2;
end
end
end
end
end

3
matlab/+occbin/map_regime.m
View File

@ -26,6 +26,9 @@ function [regime, regime_start, error_flag]=map_regime(binding_indicator,debug_s
% Journal of Monetary Economics 70, 22-38
error_flag=0;
if isempty(binding_indicator)
binding_indicator = false;
end
% analyse violvec and isolate contiguous periods in the other regime.
regime(1) = binding_indicator(1);
regime_index = 1;

11
matlab/+occbin/set_default_options.m
View File

@ -42,6 +42,8 @@ if ismember(flag,{'all'})
end
if ismember(flag,{'filter','all'})
options_occbin_.filter.state_covariance = false;
options_occbin_.filter.guess_regime = false;
options_occbin_.filter.use_relaxation = false;
end
@ -69,6 +71,7 @@ if ismember(flag,{'irf','all'})
end
if ismember(flag,{'likelihood','all'})
options_occbin_.likelihood.brute_force_regime_guess = true;
options_occbin_.likelihood.curb_retrench = false;
options_occbin_.likelihood.first_period_occbin_update = 1;
options_occbin_.likelihood.full_output = false;
@ -77,10 +80,12 @@ if ismember(flag,{'likelihood','all'})
options_occbin_.likelihood.init_binding_indicator = false(0);
options_occbin_.likelihood.inversion_filter = false;
options_occbin_.likelihood.IVF_shock_observable_mapping = [];
options_occbin_.likelihood.loss_function_regime_guess = false;
options_occbin_.likelihood.maxit = 30; % this is for occbin solver algo
options_occbin_.likelihood.max_number_of_iterations = 10; % this is for occbin_kalman_update loop
options_occbin_.likelihood.max_check_ahead_periods=inf;
options_occbin_.likelihood.periods = 100;
options_occbin_.likelihood.number_of_initial_periods_with_extra_regime_guess=0;
options_occbin_.likelihood.periods = 3;
options_occbin_.likelihood.check_ahead_periods=200;
options_occbin_.likelihood.periodic_solution=false;
options_occbin_.likelihood.piecewise_only = true;
@ -193,6 +198,8 @@ if ismember(flag,{'simul','all'})
options_occbin_.simul.periods=100;
options_occbin_.simul.check_ahead_periods=200;
options_occbin_.simul.periodic_solution=false;
options_occbin_.simul.periodic_solution_threshold=1;
options_occbin_.simul.periodic_solution_strict=true;
options_occbin_.simul.piecewise_only = false;
options_occbin_.simul.reset_check_ahead_periods_in_new_period = false;
options_occbin_.simul.reset_regime_in_new_period = false;
@ -215,7 +222,7 @@ if ismember(flag,{'smoother','all'})
options_occbin_.smoother.maxit = 30; % this is for occbin solver algo
options_occbin_.smoother.max_check_ahead_periods=inf;
options_occbin_.smoother.max_number_of_iterations = 10; % this is for smoother loop
options_occbin_.smoother.periods = 100;
options_occbin_.smoother.periods = 3;
options_occbin_.smoother.check_ahead_periods=200;
options_occbin_.smoother.periodic_solution=false;
options_occbin_.smoother.piecewise_only = true;

26
matlab/+occbin/solve_one_constraint.m
View File

@ -136,6 +136,8 @@ for shock_period = 1:n_shocks_periods
binding_indicator_history={};
max_err = NaN(max_iter,1);
regime_violates_constraint_in_expectation = false(max_iter,1);
number_of_periods_with_violations = zeros(max_iter,1);
regime_exit_period = ones(max_iter,1);
while (regime_change_this_iteration && iter<max_iter && ~is_periodic && ~is_periodic_loop)
iter = iter +1;
@ -194,9 +196,11 @@ for shock_period = 1:n_shocks_periods
if ~isinf(opts_simul_.max_check_ahead_periods) && opts_simul_.max_check_ahead_periods<length(binding_indicator)
end_periods = opts_simul_.max_check_ahead_periods;
last_indicator = false(length(binding_indicator)-end_periods,1);
regime_violates_constraint_in_last_period(iter) = any(binding.constraint_1(end_periods+1:end));
else
end_periods = length(binding_indicator);
last_indicator = false(0);
regime_violates_constraint_in_last_period(iter) = binding.constraint_1(end_periods);
end
% check if changes to the hypothesis of the duration for each
% regime
@ -215,6 +219,9 @@ for shock_period = 1:n_shocks_periods
% if max_check_ahead_periods<inf, enforce unconstrained regime for period larger than max_check_ahead_periods
binding_constraint_new=[binding.constraint_1(1:end_periods);last_indicator];
relaxed_constraint_new = [relax.constraint_1(1:end_periods);not(last_indicator)];
number_of_periods_with_violations(iter) = sum(binding_indicator -((binding_indicator | binding_constraint_new) & ~(binding_indicator & relaxed_constraint_new)));
regime_exit_period(iter) = max(regime_history(shock_period).regimestart);
if curb_retrench % apply Gauss-Seidel idea of slowing down the change in the guess
% for the constraint -- only relax one
% period at a time starting from the last
@ -252,20 +259,35 @@ for shock_period = 1:n_shocks_periods
if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
% vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
% is_periodic(kiter) = isequal(vvv, binding_indicator);
is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}-binding_indicator))==1;
is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && sum(binding_indicator_history{iter}-binding_indicator)<=opts_simul_.periodic_solution_threshold;
else
is_periodic(kiter)=false;
end
end
is_periodic_all =is_periodic;
is_periodic = any(is_periodic);
is_periodic_strict = is_periodic;
if is_periodic_loop && ~is_periodic
if ~opts_simul_.periodic_solution_strict && any(number_of_periods_with_violations(~regime_violates_constraint_in_expectation(1:iter))<=opts_simul_.periodic_solution_threshold)
is_periodic=true;
is_periodic_all=false(size(is_periodic_loop_all));
is_periodic_all(1) = true; % here we consider all guesses and pick the best one according to the criteria below
else
do_nothing=true;
end
end
if is_periodic && periodic_solution
inx = find(is_periodic_all,1):iter;
inx1 = inx(find(~regime_violates_constraint_in_expectation(inx)));
if not(isempty(inx1))
inx=inx1;
end
[merr,imerr]=min(max_err(inx));
if is_periodic_strict || isempty(inx1)
[merr,imerr]=min(max_err(inx));
else
[merr,imerr]=min(number_of_periods_with_violations(inx));
end
inx = inx(imerr);
binding_indicator=binding_indicator_history{inx};
if inx<iter

28
matlab/+occbin/solve_two_constraints.m
View File

@ -228,7 +228,7 @@ for shock_period = 1:n_shocks_periods
err_binding_constraint_new = [err.binding_constraint_1(1:end_periods); err.binding_constraint_2(1:end_periods)];
err_relaxed_constraint_new = [err.relax_constraint_1(1:end_periods); err.relax_constraint_2(1:end_periods)];
% check if changes_
if any(binding_constraint_new & ~my_binding_indicator(:)) || any(relaxed_constraint_new & my_binding_indicator(:))
err_violation = err_binding_constraint_new(binding_constraint_new & ~my_binding_indicator(:));
@ -243,6 +243,11 @@ for shock_period = 1:n_shocks_periods
binding_constraint_new=[binding.constraint_1(1:end_periods);last_indicator; binding.constraint_2(1:end_periods);last_indicator];
relaxed_constraint_new = [relax.constraint_1(1:end_periods);not(last_indicator); relax.constraint_2(1:end_periods);not(last_indicator)];
tmp_nper(1) = sum(binding_indicator(:,1) - (binding_indicator(:,1) | [binding.constraint_1(1:end_periods);last_indicator]) & ~(binding_indicator(:,1) & [relax.constraint_1(1:end_periods);not(last_indicator)]));
tmp_nper(2) = sum(binding_indicator(:,2) - (binding_indicator(:,2) | [binding.constraint_2(1:end_periods);last_indicator]) & ~(binding_indicator(:,2) & [relax.constraint_2(1:end_periods);not(last_indicator)]));
number_of_periods_with_violations(iter) = max(tmp_nper);
regime_exit_period(iter,1) = max(regime_history(shock_period).regimestart1);
regime_exit_period(iter,2) = max(regime_history(shock_period).regimestart2);
if curb_retrench % apply Gauss-Seidel idea of slowing down the change in the guess
% for the constraint -- only relax one
% period at a time starting from the last
@ -276,7 +281,7 @@ for shock_period = 1:n_shocks_periods
is_periodic_loop(kiter) = false;
end
end
% is_periodic_loop_all =is_periodic_loop;
is_periodic_loop_all =is_periodic_loop;
is_periodic_loop = any(is_periodic_loop);
% only accept periodicity where regimes differ by one
% period!
@ -285,20 +290,35 @@ for shock_period = 1:n_shocks_periods
if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
% vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
% is_periodic(kiter) = isequal(vvv, binding_indicator);
is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}(:,1)-binding_indicator(:,1)))<=1 && length(find(binding_indicator_history{iter}(:,2)-binding_indicator(:,2)))<=1;
is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && sum(binding_indicator_history{iter}(:,1)-binding_indicator(:,1))<=opts_simul_.periodic_solution_threshold && sum(binding_indicator_history{iter}(:,2)-binding_indicator(:,2))<=opts_simul_.periodic_solution_threshold;
else
is_periodic(kiter)=false;
end
end
is_periodic_all = is_periodic;
is_periodic = any(is_periodic);
is_periodic_strict = is_periodic;
if is_periodic_loop && ~is_periodic
if ~opts_simul_.periodic_solution_strict && any(number_of_periods_with_violations(~regime_violates_constraint_in_expectation(1:iter))<opts_simul_.periodic_solution_threshold)
is_periodic=true;
is_periodic_all=false(size(is_periodic_loop_all));
is_periodic_all(1) = true; % here we consider all guesses and pick the best one according to the criteria below
else
do_nothing=true;
end
end
if is_periodic && periodic_solution
inx = find(is_periodic_all,1):iter;
inx1 = inx(find(~regime_violates_constraint_in_expectation(inx)));
if not(isempty(inx1))
inx=inx1;
end
[min_err,index_min_err]=min(max_err(inx));
if is_periodic_strict || isempty(inx1)
[min_err,index_min_err]=min(max_err(inx));
else
[min_err,index_min_err]=min(number_of_periods_with_violations(inx));
end
inx = inx(index_min_err);
binding_indicator=binding_indicator_history{inx}; %select regime history with same result, but smallest error
if inx<iter %update if needed

32
matlab/+occbin/solver.m
View File

@ -1,23 +1,26 @@
function [dr, out, ss] = solver(M_, options_, dr ,steady_state, exo_steady_state, exo_det_steady_state)
% function [oo_, out, ss] = solver(M_,oo_,options_, dr ,steady_state, exo_steady_state, exo_det_steady_state
% [dr, out, ss] = solver(M_,oo_,options_, dr ,steady_state, exo_steady_state, exo_det_steady_state
% Solves the model with an OBC and produces simulations/IRFs
%
% INPUT:
% - M_ [structure] Matlab's structure describing the model
% - oo_ [structure] Matlab's structure containing the results
% - options_ [structure] Matlab's structure containing the options
% - M_ [structure] Matlab's structure describing the model
% - options_ [structure] Matlab's structure containing the options
% - dr [structure] model information structure
% - endo_steady_state [vector] steady state value for endogenous variables
% - exo_steady_state [vector] steady state value for exogenous variables
% - exo_det_steady_state [vector] steady state value for exogenous deterministic variables
%
% OUTPUT:
% - oo_ [structure] Matlab's structure containing the results
% - out [structure] simulation result containing fields:
% - linear: paths for endogenous variables ignoring OBC (linear solution)
% - piecewise: paths for endogenous variables satisfying the OBC (occbin/piecewise solution)
% - ys: vector of steady state values
% - regime_history: information on number and time of regime transitions
% - ss [structure] State space solution
% - T: [n_vars by n_vars by n_shock_period] array of transition matrices
% - R: [n_vars by n_exo by n_shock_period] array of shock response matrices
% - C: [n_vars by n_shock_period] array of constants
% - dr [structure] decision rules
% - out [structure] simulation result containing fields:
% - linear: paths for endogenous variables ignoring OBC (linear solution)
% - piecewise: paths for endogenous variables satisfying the OBC (occbin/piecewise solution)
% - ys: vector of steady state values
% - regime_history: information on number and time of regime transitions
% - ss [structure] State space solution
% - T: [n_vars by n_vars by n_shock_period] array of transition matrices
% - R: [n_vars by n_exo by n_shock_period] array of shock response matrices
% - C: [n_vars by n_shock_period] array of constants
% Copyright © 2021-2023 Dynare Team
%
@ -50,6 +53,7 @@ else
end
end
ss=struct();
if solve_dr
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;

6
matlab/estimation/dsge_likelihood.m
View File

@ -906,7 +906,5 @@ occbin_options.opts_simul.restrict_state_space = options_.occbin.likelihood.rest
occbin_options.opts_simul.full_output = options_.occbin.likelihood.full_output;
occbin_options.opts_simul.piecewise_only = options_.occbin.likelihood.piecewise_only;
if ~isempty(options_.occbin.smoother.init_binding_indicator)
occbin_options.opts_simul.init_binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_simul.init_regime_history=options_.occbin.likelihood.init_regime_history;
end
occbin_options.opts_regime.init_binding_indicator = options_.occbin.likelihood.init_binding_indicator;
occbin_options.opts_regime.init_regime_history=options_.occbin.likelihood.init_regime_history;

24
matlab/kalman/likelihood/missing_observations_kalman_filter.m
View File

@ -93,6 +93,7 @@ F_singular = true;
s = 0;
rescale_prediction_error_covariance0=rescale_prediction_error_covariance;
if occbin_.status
base_regime = struct();
Qt = repmat(Q,[1 1 3]);
a0 = zeros(mm,last);
a1 = zeros(mm,last);
@ -107,8 +108,9 @@ if occbin_.status
exo_det_steady_state=occbin_.info{5};
M_=occbin_.info{6};
occbin_options=occbin_.info{7};
opts_regime.regime_history = occbin_options.opts_simul.init_regime;
opts_regime.binding_indicator = occbin_options.opts_simul.init_binding_indicator;
occbin_options.opts_simul.SHOCKS = [];
opts_regime.regime_history = occbin_options.opts_regime.init_regime_history;
opts_regime.binding_indicator = occbin_options.opts_regime.init_binding_indicator;
if t>1
first_period_occbin_update = max(t+1,options_.occbin.likelihood.first_period_occbin_update);
else
@ -134,6 +136,15 @@ if occbin_.status
end
end
if M_.occbin.constraint_nbr==1
base_regime.regime = 0;
base_regime.regimestart = 1;
else
base_regime.regime1 = 0;
base_regime.regimestart1 = 1;
base_regime.regime2 = 0;
base_regime.regimestart2 = 1;
end
else
first_period_occbin_update = inf;
C=0;
@ -254,12 +265,16 @@ while notsteady && t<=last
RR01 = cat(3,R,RR(:,:,1));
CC01 = zeros(size(CC,1),2);
CC01(:,2) = CC(:,1);
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a00, a10, P00, P10, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), M_, dr, endo_steady_state,exo_steady_state,exo_det_steady_state, options_, occbin_options);
% insert here kalman update engine
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx] = occbin.kalman_update_engine(a00, a10, P00, P10, t, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), base_regime, d_index, M_, dr, endo_steady_state,exo_steady_state,exo_det_steady_state, options_, occbin_options);
% [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a00, a10, P00, P10, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), M_, dr, endo_steady_state,exo_steady_state,exo_det_steady_state, options_, occbin_options);
else
if isqvec
Qt = Qvec(:,:,t-1:t+1);
end
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
% insert here kalman update engine
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regx, info, M_, likx] = occbin.kalman_update_engine(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),t,data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),base_regime,d_index,M_,dr, endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
% [ax, a1x, Px, P1x, vx, Tx, Rx, Cx, regimes_(t:t+2), info, M_, likx] = occbin.kalman_update_algo_1(a0(:,t-1),a1(:,t-1:t),P0(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,vv(:,t-1:t),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options);
end
if info
if options_.debug
@ -270,6 +285,7 @@ while notsteady && t<=last
if options_.occbin.likelihood.use_updated_regime
lik(s) = likx;
end
regimes_(t:t+2) = regx;
a0(:,t) = ax(:,1);
a1(:,t) = a1x(:,2);
a = ax(:,2);

57
matlab/kalman/missing_DiffuseKalmanSmootherH3_Z.m
View File

@ -153,9 +153,9 @@ alphahat0=[];
aalphahat0=[];
V0=[];
C=0;
if ~occbin_.status
isoccbin = 0;
C=0;
TT=[];
RR=[];
CC=[];
@ -171,6 +171,15 @@ else
occbin_options=occbin_.info{7};
opts_regime = occbin_options.opts_regime;
% first_period_occbin_update = inf;
if M_.occbin.constraint_nbr==1
base_regime.regime = 0;
base_regime.regimestart = 1;
else
base_regime.regime1 = 0;
base_regime.regimestart1 = 1;
base_regime.regime2 = 0;
base_regime.regimestart2 = 1;
end
if isfield(opts_regime,'regime_history') && ~isempty(opts_regime.regime_history)
opts_regime.regime_history=[opts_regime.regime_history(1) opts_regime.regime_history];
else
@ -239,8 +248,19 @@ if newRank
% add this to get smoothed states in period 0
Pinf_init = Pinf(:,:,1);
Pstar_init = Pstar(:,:,1);
Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ;
ainit = a1(:,1);
if isoccbin && length(occbin_.info)>6
if isqvec
QQ = RR(:,:,1)*Qvec(:,:,1)*transpose(RR(:,:,1));
else
QQ = RR(:,:,1)*Q*transpose(RR(:,:,1));
end
T = TT(:,:,1);
C = CC(:,1);
a1(:,1) = T*a(:,1)+C; %transition according to (6.14) in DK (2012)
% Pinf(:,:,1) = T*Pinf(:,:,1)*T';
end
Pstar(:,:,1) = T*Pstar(:,:,1)*T' + QQ;
end
while newRank && t < smpl
@ -285,10 +305,20 @@ while newRank && t < smpl
oldRank = 0;
end
if isoccbin
TT(:,:,t+1)= T;
RR(:,:,t+1)= R;
if length(occbin_.info)>9
if isqvec
QQ = RR(:,:,t+1)*Qvec(:,:,t+1)*transpose(RR(:,:,t+1));
else
QQ = RR(:,:,t+1)*Q*transpose(RR(:,:,t+1));
end
T = TT(:,:,t+1);
C = CC(:,t+1);
else
TT(:,:,t+1)= T;
RR(:,:,t+1)= R;
end
end
a1(:,t+1) = T*a(:,t);
a1(:,t+1) = T*a(:,t)+C;
aK(1,:,t+1) = a1(:,t+1);
for jnk=2:nk
aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
@ -311,13 +341,14 @@ while newRank && t < smpl
end
end
d = t;
if isoccbin
first_period_occbin_update = occbin_options.first_period_occbin_update;
if d>0
first_period_occbin_update = max(t+2,occbin_options.first_period_occbin_update);
% kalman update is not yet robust to accommodate diffuse steps
end
if occbin_options.opts_regime.waitbar
if occbin_options.opts_simul.waitbar && first_period_occbin_update<smpl
hh_fig = dyn_waitbar(0,'Occbin: Piecewise Kalman Filter');
set(hh_fig,'Name','Occbin: Piecewise Kalman Filter.');
waitbar_indicator=1;
@ -328,7 +359,6 @@ else
first_period_occbin_update = inf;
waitbar_indicator=0;
end
d = t;
P(:,:,d+1) = Pstar(:,:,d+1);
Fstar = Fstar(:,1:d);
Finf = Finf(:,1:d);
@ -351,7 +381,7 @@ while notsteady && t<smpl
if waitbar_indicator
dyn_waitbar(t/smpl, hh_fig, sprintf('Period %u of %u', t,smpl));
end
occbin_options.opts_regime.waitbar=0;
occbin_options.opts_simul.waitbar=0;
if t==1
if isqvec
Qt = cat(3,Q,Qvec(:,:,t:t+1));
@ -371,12 +401,16 @@ while notsteady && t<smpl
RR01 = cat(3,R,RR(:,:,1));
CC01 = zeros(size(CC,1),2);
CC01(:,2) = CC(:,1);
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a10,P0,P10,data_index0,Z,v0,Fi0,Ki0,Y0,H,Qt,T0,R0,TT01,RR01,CC01,regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
% [ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a0,a10,P0,P10,data_index0,Z,v0,Fi0,Ki0,Y0,H,Qt,T0,R0,TT01,RR01,CC01,regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, tmp, error_flag, M_, likx, etaha, aha, V, Fix, Kix, TTx,RRx,CCx] = occbin.kalman_update_engine(a0, a10, P0, P10, t, data_index0, Z, v0, Y0, H, Qt, T0, R0, TT01, RR01, CC01, regimes_(t:t+1), base_regime, di, M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_, occbin_options, ...
Fi0,Ki0,kalman_tol,nk);
else
if isqvec
Qt = Qvec(:,:,t-1:t+1);
end
[ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,v(:,t-1:t),Fi(:,t-1),Ki(:,:,t-1),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
% [ax, a1x, Px, P1x, vx, Fix, Kix, Tx, Rx, Cx, tmp, error_flag, M_, aha, etaha,TTx,RRx,CCx] = occbin.kalman_update_algo_3(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),data_index(t-1:t),Z,v(:,t-1:t),Fi(:,t-1),Ki(:,:,t-1),Y(:,t-1:t),H,Qt,T0,R0,TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t),regimes_(t:t+1),M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,occbin_options,kalman_tol,nk);
[ax, a1x, Px, P1x, vx, Tx, Rx, Cx, tmp, error_flag, M_, likx, etaha, aha, V, Fix, Kix,TTx,RRx,CCx] = occbin.kalman_update_engine(a(:,t-1),a1(:,t-1:t),P(:,:,t-1),P1(:,:,t-1:t),t,data_index(t-1:t),Z,v(:,t-1:t), Y(:,t-1:t), H, Qt, T0, R0, TT(:,:,t-1:t),RR(:,:,t-1:t),CC(:,t-1:t), regimes_(t:t+1), base_regime, di, M_, dr,endo_steady_state,exo_steady_state,exo_det_steady_state, options_, occbin_options, ...
Fi(:,t-1),Ki(:,:,t-1),kalman_tol,nk);
end
if ~error_flag
regimes_(t:t+2)=tmp;
@ -425,6 +459,7 @@ while notsteady && t<smpl
T = TT(:,:,t);
C = CC(:,t);
a1(:,t) = T*a(:,t)+C; %transition according to (6.14) in DK (2012)
a(:,t) = a1(:,t);
P(:,:,t) = T*P(:,:,t)*T' + QQ; %transition according to (6.14) in DK (2012)
P1(:,:,t) = P(:,:,t);
end
@ -497,7 +532,7 @@ while notsteady && t<smpl
end
aK(1,:,t+1) = a1(:,t+1);
if ~isempty(nk) && nk>1 && isoccbin && (t>=first_period_occbin_update || isinf(first_period_occbin_update))
opts_simul = occbin_options.opts_regime;
opts_simul = occbin_options.opts_simul;
opts_simul.SHOCKS = zeros(nk,M_.exo_nbr);
if smoother_redux
tmp=zeros(M_.endo_nbr,1);

3
tests/occbin/model_irrcap_twoconstraints/dynrbc_common.inc
View File

@ -89,6 +89,9 @@ occbin_graph(noconstant) c erra lambdak k i a k;
occbin_solver(simul_periods=200,simul_maxit=200,simul_curb_retrench,simul_check_ahead_periods=200);
occbin_setup(smoother_periods=200,smoother_maxit=200,smoother_curb_retrench,smoother_check_ahead_periods=200);
calib_smoother(datafile=datasim);
oo0=oo_;
occbin_setup(filter_use_relaxation);
calib_smoother(datafile=datasim);
oo_= occbin.unpack_simulations(M_,oo_,options_);
titlelist = char('c','lambdak','k','i','a','erra');