From 3eedbbf1ae859aaa96805daf6fe3ad3f251b3751 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@ec.europa.eu>
Date: Mon, 26 Aug 2024 15:57:35 +0200
Subject: [PATCH] better handling errors in occbin update step
(cherry picked from commit 5ad067f1eebd6647e7976cdff9944594df0a4e9a)
---
matlab/+occbin/kalman_update_algo_1.m | 18 ++++++++++++++----
matlab/+occbin/kalman_update_algo_3.m | 2 +-
matlab/+occbin/kalman_update_engine.m | 9 +++++++++
matlab/get_error_message.m | 10 ++++++++++
.../missing_observations_kalman_filter.m | 6 ------
5 files changed, 34 insertions(+), 11 deletions(-)
diff --git a/matlab/+occbin/kalman_update_algo_1.m b/matlab/+occbin/kalman_update_algo_1.m
index 11205e63e0..83cbd587ba 100644
--- a/matlab/+occbin/kalman_update_algo_1.m
+++ b/matlab/+occbin/kalman_update_algo_1.m
@@ -213,7 +213,12 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
if M_.occbin.constraint_nbr==1
regime_end(niter) = regimes_(1).regimestart(end);
end
- [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);
+ [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);
+ if error_flag
+ etahat=NaN(size(QQQ,1),1);
+ warning(orig_warning_state);
+ return;
+ end
etahat_hist(niter) = {etahat};
lik_hist(niter) = lik;
opts_simul.SHOCKS(1,:) = etahat(:,2)';
@@ -278,7 +283,12 @@ if any(myregime) || ~isequal(regimes_(1),regimes0(1))
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);
+ [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);
+ if error_flag
+ etahat=NaN(size(QQQ,1),1);
+ warning(orig_warning_state);
+ return;
+ end
end
else
error_flag = 330;
@@ -294,7 +304,7 @@ 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;
+ error_flag = 331;
end
if ~error_flag
@@ -357,7 +367,7 @@ else
F = ZZ*P(:,:,t)*ZZ' + H(di,di);
sig=sqrt(diag(F));
if any(any(isnan(F)))
- error_flag=1;
+ error_flag=325;
warning(orig_warning_state);
return;
end
diff --git a/matlab/+occbin/kalman_update_algo_3.m b/matlab/+occbin/kalman_update_algo_3.m
index 457c474806..efffc9badb 100644
--- a/matlab/+occbin/kalman_update_algo_3.m
+++ b/matlab/+occbin/kalman_update_algo_3.m
@@ -282,7 +282,7 @@ end
error_flag = out.error_flag;
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;
+ error_flag = 331;
end
if ~error_flag
diff --git a/matlab/+occbin/kalman_update_engine.m b/matlab/+occbin/kalman_update_engine.m
index fe2346436a..92c13d72e1 100644
--- a/matlab/+occbin/kalman_update_engine.m
+++ b/matlab/+occbin/kalman_update_engine.m
@@ -328,4 +328,13 @@ if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %||
end
end
end
+
+if info(1)==0
+ if isnan(likx)
+ info = 323;
+ elseif any(any(isnan(ax)))
+ info = 324;
+ end
+end
+
end
diff --git a/matlab/get_error_message.m b/matlab/get_error_message.m
index 1694dde4e2..f15f8c34a7 100644
--- a/matlab/get_error_message.m
+++ b/matlab/get_error_message.m
@@ -198,6 +198,16 @@ switch info(1)
message = 'Occbin: there was a problem in running the smoother. Simulation within smoother failed.';
case 322
message = 'Occbin: smoother did not converge.';
+ case 323
+ message = 'Piecewise linear Kalman filter: the likelihood is NaN.';
+ case 324
+ message = 'Piecewise linear Kalman filter: updated state vector is NaN.';
+ case 325
+ message = 'Piecewise linear Kalman filter: filter covariance NaN.';
+ case 330
+ message = 'Piecewise linear Kalman filter: update step did not reach a fixed point (periodic loop).';
+ case 331
+ message = 'Piecewise linear Kalman filter: update step did not reach a fixed point (max number of iterations reached).';
case 401
message = 'Cycle reduction reached the iteration limit. Try increasing maxit.';
case 402
diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter.m b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
index adfdda7652..67b3adf614 100644
--- a/matlab/kalman/likelihood/missing_observations_kalman_filter.m
+++ b/matlab/kalman/likelihood/missing_observations_kalman_filter.m
@@ -276,12 +276,6 @@ while notsteady && t<=last
[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 isnan(likx) || any(any(isnan(ax)))
- if options_.debug
- fprintf('\nmissing_observations_kalman_filter:PKF failed in period %u with: NaN in likelihood value\n', t);
- end
- return
- end
if info
if options_.debug
fprintf('\nmissing_observations_kalman_filter:PKF failed in period %u with: %s\n', t, get_error_message(info,options_));
--
GitLab