better handling errors in occbin update step

(cherry picked from commit 5ad067f1)

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

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

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

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

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_));