matlab/+identification/simulated_moment_uncertainty.m

@@ -34,7 +34,7 @@ mm=zeros(length(indx),replic);
 disp('Evaluating simulated moment uncertainty ... please wait')
 disp(['Doing ',int2str(replic),' replicas of length ',int2str(periods),' periods.'])
 
-h = dyn_waitbar(0,'Simulated moment uncertainty ...');
+h = waitbar.run(0,'Simulated moment uncertainty ...');
 %Do check whether simulation is possible
 if options_.periods == 0
     error('simulated_moment_uncertainty: Periods must be bigger than 0')
@@ -96,9 +96,9 @@ for j=1:replic
         dum=[dum; vec(oo_.autocorr{i}.*(sd*sd'))];
     end
     mm(:,j)=dum(indx);
-    dyn_waitbar(j/replic,h,['Simulated moment uncertainty. Replic  ',int2str(j),'/',int2str(replic)])
+    waitbar.run(j/replic,h,['Simulated moment uncertainty. Replic  ',int2str(j),'/',int2str(replic)])
 end
-dyn_waitbar_close(h);
+waitbar.close(h);
 
 cmm = cov(mm');
 disp('Simulated moment uncertainty ... done!')

matlab/+identification/unfold_g3.m

 function g3_unfolded = unfold_g3(g3, ny)
 % Given the 3rd order derivatives stored in a sparse matrix and without
 % symmetric elements (as returned by the static/dynamic files) and the number
-% of (static or dynamic )variables in the jacobian, returns
+% of (static or dynamic )variables in the Jacobian, returns
 % an unfolded version of the same matrix (i.e. with symmetric elements).
 
 % Copyright © 2019 Dynare Team

matlab/+identification/unfold_g4.m

 function g4_unfolded = unfold_g4(g4, ny)
 % Given the 4th order derivatives stored in a sparse matrix and without
 % symmetric elements (as returned by the static/dynamic files) and the number
-% of (static or dynamic) variables in the jacobian, returns
+% of (static or dynamic) variables in the Jacobian, returns
 % an unfolded version of the same matrix (i.e. with symmetric elements).
 
 % Copyright © 2019 Dynare Team

matlab/+mom/mode_compute_gmm_smm.m

@@ -37,7 +37,7 @@ function [xparam1, weighting_info, mom_verbose] = mode_compute_gmm_smm(xparam0,
 % o prior_dist_names
 % -------------------------------------------------------------------------
 
-% Copyright © 2023 Dynare Team
+% Copyright © 2023-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -109,6 +109,7 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
         if options_mom_.optimizer_vec{optim_iter} == 0
             xparam1 = xparam0; % no minimization, evaluate objective at current values            
             fval = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
+            optimization_info = [];
         else
             if options_mom_.optimizer_vec{optim_iter} == 13
                 options_mom_.mom.vector_output = true;
@@ -120,7 +121,8 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
             else
                 options_mom_.mom.compute_derivs = false;
             end
-            [xparam1, fval] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, [],...
+            [xparam1, fval, ~, ~, ~, ~, optimization_info] = ...
+                dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, [], ...
                                           data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
             if options_mom_.mom.vector_output
                 fval = fval'*fval;
@@ -133,6 +135,7 @@ for stage_iter = 1:size(options_mom_.mom.weighting_matrix,1)
             tbl_title_iter = sprintf('FREQUENTIST %s (STAGE %d ITERATION %d) VERBOSE',options_mom_.mom.mom_method,stage_iter,optim_iter);
             field_name_iter = sprintf('%s_stage_%d_iter_%d',lower(options_mom_.mom.mom_method),stage_iter,optim_iter);
             mom_verbose.(field_name_iter) = display_estimation_results_table(xparam1,std_via_invhessian_xparam1_iter,M_,options_mom_,estim_params_,bayestopt_,[],prior_dist_names,tbl_title_iter,field_name_iter);
+            mom_verbose.(field_name_iter).optimization_info = optimization_info;
         end
         xparam0 = xparam1;
     end

matlab/+mom/mode_compute_irf_matching.m

@@ -9,7 +9,7 @@ function [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matchi
 % -------------------------------------------------------------------------
 % INPUTS
 % xparam0:               [vector]       initialized parameters
-% hessian_xparam0:       [matrix]       initialized hessian at xparam0
+% hessian_xparam0:       [matrix]       initialized Hessian at xparam0
 % objective_function:    [func handle]  name of the objective function
 % doBayesianEstimation:  [logical]      true if Bayesian estimation
 % weighting_info:        [structure]    information on weighting matrix
@@ -26,7 +26,7 @@ function [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matchi
 % -------------------------------------------------------------------------
 % OUTPUT
 % xparam1:               [vector]       mode of objective function
-% hessian_xparam1:       [matrix]       hessian at xparam1
+% hessian_xparam1:       [matrix]       Hessian at xparam1
 % fval:                  [double]       function value at mode
 % mom_verbose:           [structure]    information on intermediate estimation results
 % Also saves the computed mode and hessian to a file.
@@ -41,7 +41,7 @@ function [xparam1, hessian_xparam1, fval, mom_verbose] = mode_compute_irf_matchi
 %  o mom.objective_function
 % -------------------------------------------------------------------------
 
-% Copyright © 2023 Dynare Team
+% Copyright © 2023-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -67,8 +67,10 @@ for optim_iter = 1:length(options_mom_.optimizer_vec)
         xparam1 = xparam0;
         hessian_xparam1 = hessian_xparam0;
         fval = feval(objective_function, xparam1, data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
+        optimization_info = [];
     else
-        [xparam1, fval, exitflag, hessian_xparam1, options_mom_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, hessian_xparam0,...
+        [xparam1, fval, exitflag, hessian_xparam1, ~, ~, optimization_info] = ...
+            dynare_minimize_objective(objective_function, xparam0, options_mom_.optimizer_vec{optim_iter}, options_mom_, [BoundsInfo.lb BoundsInfo.ub], bayestopt_.name, bayestopt_, hessian_xparam0, ...
                                       data_moments, weighting_info, options_mom_, M_, estim_params_, bayestopt_, BoundsInfo, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
     end
     fprintf('\nMode Compute Iteration %d: Value of minimized moment distance objective function: %12.10f.\n',optim_iter,fval);
@@ -98,6 +100,7 @@ for optim_iter = 1:length(options_mom_.optimizer_vec)
             field_name_iter = sprintf('iter_%d',optim_iter);
         end
         mom_verbose.(field_name_iter) = display_estimation_results_table(xparam1,std_via_invhessian_xparam1_iter,M_,options_mom_,estim_params_,bayestopt_,[],prior_dist_names,tbl_title_iter,field_name_iter);
+        mom_verbose.(field_name_iter).optimization_info = optimization_info;
     end
     xparam0 = xparam1;
     hessian_xparam0 = hessian_xparam1;

matlab/+mom/objective_function.m

@@ -22,11 +22,11 @@ function [fval, info, exit_flag, df, junk_hessian, Q, model_moments, model_momen
 %  - fval:                         [double]  value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
 %  - info:                         [vector]  information on error codes and penalties
 %  - exit_flag:                    [double]  flag for exit status (0 if error, 1 if no error)
-%  - df:                           [matrix]  analytical jacobian of the moment difference (wrt paramters), currently for GMM only
+%  - df:                           [matrix]  analytical Jacobian of the moment difference (wrt parameters), currently for GMM only
 %  - junk_hessian:                 [matrix]  empty matrix required for optimizer interface (Hessian would typically go here)
 %  - Q:                            [double]  value of the quadratic form of the moment difference
 %  - model_moments:                [vector]  model moments
-%  - model_moments_params_derivs:  [matrix]  analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
+%  - model_moments_params_derivs:  [matrix]  analytical Jacobian of the model moments wrt estimated parameters (currently for GMM only)
 %  - irf_model_varobs:             [matrix]  model IRFs for observable variables (used for plotting matched IRfs in mom.run)
 % -------------------------------------------------------------------------
 % This function is called by
@@ -116,7 +116,11 @@ if info(1)
             info(1) == 81 || info(1) == 84 ||  info(1) == 85 ||  info(1) == 86
         % meaningful second entry of output that can be used
         fval = Inf;
+        if ~isfinite(info(2))
+            info(4) = 0.1;
+        else
             info(4) = info(2);
+        end
         exit_flag = 0;
         if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
             fval = ones(options_mom_.mom.mom_nbr,1)*options_mom_.huge_number;
@@ -147,7 +151,7 @@ if strcmp(options_mom_.mom.mom_method,'GMM')
         if ~isempty(estim_params_.var_exo)
             indpstderr = estim_params_.var_exo(:,1); % values correspond to varexo declaration order, row number corresponds to order in estimated_params
         end
-        indpcorr=[]; % initialize matrix for corr paramters
+        indpcorr=[]; % initialize matrix for corr parameters
         if ~isempty(estim_params_.corrx)
             indpcorr = estim_params_.corrx(:,1:2); % values correspond to varexo declaration order, row number corresponds to order in estimated_params
         end
@@ -293,7 +297,7 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING') && strcmp(options_mom_.mom
         [model_irf, check] = feval(str2func(options_mom_.mom.irf_matching_file.name), model_irf, M_, options_mom_, dr.ys);
         if check
             fval = Inf;
-            info(1) = 180;
+            info(1) = 182;
             info(4) = 0.1;
             exit_flag = 0;
             if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
@@ -330,11 +334,36 @@ if strcmp(options_mom_.mom.mom_method,'IRF_MATCHING')
     end
     % add log prior if necessary
     lnprior = priordens(xparam,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
+    if isinf(lnprior)
+        fval = Inf; info(1) = 40; info(4) = 0.1; exit_flag = 0;
+        return
+    end
+    if isnan(lnprior)
+        fval = Inf; info(1) = 47; info(4) = 0.1; exit_flag = 0;
+        return
+    end
+    if imag(lnprior)~=0
+        fval = Inf; info(1) = 48; info(4) = 0.1; exit_flag = 0;
+        return
+    end
+    % compute the posterior kernel
     fval = - (lnlik + lnprior);
 
 elseif strcmp(options_mom_.mom.mom_method,'GMM') || strcmp(options_mom_.mom.mom_method,'SMM')
     residuals = sqrt(options_mom_.mom.weighting_matrix_scaling_factor)*weighting_info.Sw*moments_difference;
     Q = residuals'*residuals;
+    if isinf(Q)
+        fval = Inf; info(1) = 177; info(4) = 0.1; exit_flag = 0;
+        return
+    end
+    if isnan(Q)
+        fval = Inf; info(1) = 178; info(4) = 0.1; exit_flag = 0;
+        return
+    end
+    if imag(Q)~=0
+        fval = Inf; info(1) = 179; info(4) = 0.1; exit_flag = 0;
+        return
+    end
     if options_mom_.mom.vector_output == 1 % lsqnonlin requires vector output
         fval = residuals;
         if options_mom_.mom.penalized_estimator

matlab/+mom/run.m

@@ -273,7 +273,7 @@ end
 % set priors and bounds over the estimated parameters
 [xparam0, estim_params_, bayestopt_, lb, ub, M_] = set_prior(estim_params_, M_, options_mom_);
 number_of_estimated_parameters = length(xparam0);
-hessian_xparam0 = []; % initialize hessian
+hessian_xparam0 = []; % initialize Hessian
 % check if enough moments for estimation
 if options_mom_.mom.mom_nbr < length(xparam0)
     skipline;
@@ -503,7 +503,7 @@ test_for_deep_parameters_calibration(M_);
 % -------------------------------------------------------------------------
 objective_function = str2func('mom.objective_function');
 try
-    % check for NaN or complex values of moment-distance-funtion evaluated at initial parameters
+    % check for NaN or complex values of moment-distance-function evaluated at initial parameters
     if strcmp(options_mom_.mom.mom_method,'SMM') || strcmp(options_mom_.mom.mom_method,'GMM')
         oo_.mom.weighting_info.Sw = eye(options_mom_.mom.mom_nbr); % initialize with identity weighting matrix
     end
@@ -733,7 +733,7 @@ if do_bayesian_estimation_mcmc
                 oo_.mom.posterior.metropolis = oo_load_mh.oo_.mom.posterior.metropolis;
             end
         end
-        [error_flag,~,options_mom_]= metropolis_draw(1,options_mom_,estim_params_,M_);
+        [options_mom_.sub_draws, error_flag]=set_number_of_subdraws(M_,options_mom_); %check whether number of sub_draws is feasible
         if ~(~isempty(options_mom_.sub_draws) && options_mom_.sub_draws==0)
             % THIS IS PROBABLY NOT USEFUL HERE AND CAN BE REMOVED (PREPROCESSOR: REMOVE bayesian_irf, moments_varendo)
             %if options_mom_.bayesian_irf

matlab/+mom/standard_errors.m

@@ -10,7 +10,7 @@ function [stderr_values, asympt_cov_mat] = standard_errors(xparam, objective_fun
 %  - xparam:               [vector]     value of estimated parameters as returned by set_prior()
 %  - objective_function    [func]       function handle with string of objective function
 %  - model_moments:        [vector]     model moments
-%  - model_moments_params_derivs:  [matrix]  analytical jacobian of the model moments wrt estimated parameters (currently for GMM only)
+%  - model_moments_params_derivs:  [matrix]  analytical Jacobian of the model moments wrt estimated parameters (currently for GMM only)
 %  - m_data                [matrix]     selected empirical moments at each point in time
 %  - data_moments:         [vector]     data with moments/IRFs to match
 %  - weighting_info:       [structure]  storing information on weighting matrices

matlab/+occbin/DSGE_smoother.m

@@ -8,13 +8,13 @@ function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,P
 % - Y             [double]        (n*T) matrix of data.
 % - data_index    [cell]          1*smpl cell of column vectors of indices.
 % - missing_value [boolean]       1 if missing values, 0 otherwise
-% - M_            [structure]     Matlab's structure describing the model (M_).
-% - oo_           [structure]     Matlab's structure containing the results (oo_).
-% - options_      [structure]     Matlab's structure describing the current options (options_).
+% - M_            [structure]     MATLAB's structure describing the model (M_).
+% - oo_           [structure]     MATLAB's structure containing the results (oo_).
+% - options_      [structure]     MATLAB's structure describing the current options (options_).
 % - bayestopt_    [structure]     describing the priors
 % - estim_params_ [structure]     characterizing parameters to be estimated
 % - dataset_      [structure]     the dataset after required transformation
-% - dataset_info  [structure]     Various informations about the dataset (descriptive statistics and missing observations)
+% - dataset_info  [structure]     Various information about the dataset (descriptive statistics and missing observations)
 %
 % OUTPUTS
 % - alphahat      [double]  (m*T) matrix, smoothed endogenous variables (a_{t|T})  (decision-rule order)
@@ -35,7 +35,7 @@ function [alphahat,etahat,epsilonhat,ahat0,SteadyState,trend_coeff,aKK,T0,R0,P,P
 % - Trend         [double] (n*T) pure trend component; stored in options_.varobs order
 % - state_uncertainty [double] (K,K,T) array, storing the uncertainty
 %                                   about the smoothed state (decision-rule order)
-% - M_            [structure] decribing the model
+% - M_            [structure] describing the model
 % - oo_           [structure] storing the results
 % - options_      [structure] describing the options
 % - bayestopt_    [structure] describing the priors
@@ -85,7 +85,7 @@ if  options_.occbin.smoother.linear_smoother && nargin==12
     oo_.occbin.linear_smoother.alphahat0=alphahat0;
     oo_.occbin.linear_smoother.state_uncertainty0=state_uncertainty0;
 
-    fprintf('\nOccbin: linear smoother done.\n')
+    disp_verbose('OccBin: linear smoother done.',options_.verbosity)
     options_.occbin.smoother.status=true;
 end
 % if init_mode
@@ -123,20 +123,20 @@ occbin_options.opts_simul = opts_simul; % this builds the opts_simul options fie
 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)
+options_.noprint = true;
+
+[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T0,R0,P,PK,decomp,Trend,state_uncertainty,oo_,bayestopt_,alphahat0,state_uncertainty0,~,error_indicator] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_,occbin_options);%     T1=TT;
+
+if error_indicator(1) || isempty(alphahat0)
+    if ~options_.occbin.smoother.linear_smoother || nargin~=12 %make sure linear smoother results are set before using them
+        options_.occbin.smoother.status=false;
+        [~,etahat,~,~,~,~,~,~,~,~,~,~,~,~,~,~,alphahat0] = ...
+            DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
+         options_.occbin.smoother.status=true;
+    else
         etahat= oo_.occbin.linear_smoother.etahat;
         alphahat0= oo_.occbin.linear_smoother.alphahat0;
+    end
     base_regime = struct();
     if M_.occbin.constraint_nbr==1
         base_regime.regime = 0;
@@ -171,18 +171,16 @@ occbin_options.first_period_occbin_update = inf;
 opts_regime.binding_indicator=[];
 regime_history0 = regime_history;
 
-fprintf('Occbin smoother iteration 1.\n')
+disp_verbose('OccBin smoother iteration 1.',options_.verbosity)
 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);
 if out.error_flag
-    fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
-    print_info(out.error_flag, options_.noprint, options_)
+    disp_verbose('OccBin smoother:: simulation within smoother did not converge.',options_.verbosity)    
     oo_.occbin.smoother.error_flag=321;
     return;
 end
@@ -225,7 +223,7 @@ end
 
 while is_changed && maxiter>iter && ~is_periodic
     iter=iter+1;
-    fprintf('Occbin smoother iteration %u.\n', iter)
+    disp_verbose(sprintf('OccBin smoother iteration %u.', iter),options_.verbosity)
     occbin_options.opts_regime.regime_history=regime_history;
     [alphahat,etahat,epsilonhat,~,SteadyState,trend_coeff,~,T0,R0,P,~,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,TT,RR,CC);
@@ -244,8 +242,7 @@ while is_changed && maxiter>iter && ~is_periodic
     options_.occbin.simul=opts_simul;
     [~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
     if out.error_flag
-        fprintf('Occbin smoother:: simulation within smoother did not converge.\n')
-        print_info(out.error_flag, false, options_)
+        disp_verbose('OccBin smoother:: simulation within smoother did not converge.',options_.verbosity)
         oo_.occbin.smoother.error_flag=321;
         return;
     end
@@ -390,22 +387,22 @@ if occbin_smoother_debug
 end
 
 if (maxiter==iter && is_changed) || is_periodic
-    disp('occbin.DSGE_smoother: smoother did not converge.')
-    fprintf('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.\n')
+    disp_verbose('occbin.DSGE_smoother: smoother did not converge.',options_.verbosity)
+    disp_verbose('occbin.DSGE_smoother: The algorithm did not reach a fixed point for the smoothed regimes.',options_.verbosity)
     if is_periodic
         oo_.occbin.smoother.error_flag=0;
-        fprintf('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.\n')
-        fprintf('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.\n')
+        disp_verbose('occbin.DSGE_smoother: For the periods indicated above, regimes loops between the "regime_" and the "regime_new_" pattern displayed above.',options_.verbosity)
+        disp_verbose('occbin.DSGE_smoother: We provide smoothed shocks consistent with "regime_" in oo_.',options_.verbosity)
     else
-        fprintf('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.\n')
+        disp_verbose('occbin.DSGE_smoother: The respective fields in oo_ will be left empty.',options_.verbosity)
         oo_.occbin.smoother=[];
         oo_.occbin.smoother.error_flag=322;
     end
 else
-    disp('occbin.DSGE_smoother: smoother converged.')
+    disp_verbose('occbin.DSGE_smoother: smoother converged.',options_.verbosity)
     oo_.occbin.smoother.error_flag=0;
     if occbin_smoother_fast && is_changed_start
-        disp('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge')
+        disp_verbose('occbin.DSGE_smoother: WARNING: fast algo is used, regime duration was not forced to converge',options_.verbosity)
     end
 end
 if (~is_changed || occbin_smoother_debug) && nargin==12
@@ -449,7 +446,7 @@ if (~is_changed || occbin_smoother_debug) && nargin==12
             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');
+                    hh_fig = dyn_figure(options_.nodisplay,'name','OccBin smoothed shocks');
                     ifig=ifig+1;
                     isub=0;
                 end

matlab/+occbin/IVF_core.m

-function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
+function [filtered_errs, resids, Emat, stateval, error_code, regime_history] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
 % [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_,err_index,filtered_errs_init,my_obs_list,obs,init_val)
 % Calls the solver to get the shocks explaining the data for the inversion filter
 %
@@ -10,12 +10,12 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, e
 %  - error_code             [4 by 1]            error code
 %
 % Inputs
-% - M_                      [structure]     Matlab's structure describing the model (M_).
+% - M_                      [structure]     MATLAB's structure describing the model (M_).
 % - dr                      [structure]     Reduced form model.
 % - 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_).
+% - options_                [structure]     MATLAB's structure describing the current options (options_).
 % - err_index               [double]        index of shocks with strictly positive variance in M_.exo_names
 % - filtered_errs_init      [T by N_obs]    initial values for the shocks
 % - my_obs_list             [cell]          names of observables
@@ -28,7 +28,7 @@ function [filtered_errs, resids, Emat, stateval, error_code] = IVF_core(M_,dr, e
 % Adapted for Dynare by Dynare Team.
 %
 % This code is in the public domain and may be used freely.
-% However the authors would appreciate acknowledgement of the source by
+% However the authors would appreciate acknowledgment of the source by
 % citation of any of the following papers:
 %
 % Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong (2019): "Likelihood evaluation of models
@@ -68,13 +68,13 @@ opts_simul.piecewise_only = 1;
 filtered_errs=zeros(sample_length,n_obs);
 
 if options_.occbin.likelihood.waitbar
-    hh_fig = dyn_waitbar(0,'IVF_core: Filtering the shocks');
+    hh_fig = waitbar.run(0,'IVF_core: Filtering the shocks');
     set(hh_fig,'Name','IVF_core: Filtering the shocks.');
 end
 
 for this_period=1:sample_length
     if options_.occbin.likelihood.waitbar
-        dyn_waitbar(this_period/sample_length, hh_fig, sprintf('Period %u of %u', this_period,sample_length));
+        waitbar.run(this_period/sample_length, hh_fig, sprintf('Period %u of %u', this_period,sample_length));
     end
     current_obs = obs(this_period,:);
     init_val_old = init_val;
@@ -92,21 +92,25 @@ for this_period=1:sample_length
         filtered_errs=NaN;
         error_code(1) = 304;
         error_code(4) = 1000;
-        if options_.occbin.likelihood.waitbar; dyn_waitbar_close(hh_fig); end
+        if this_period == 1
+            regime_history(this_period) = [];
+        end
+        if options_.occbin.likelihood.waitbar; waitbar.close(hh_fig); end
         return
     end
     filtered_errs(this_period,inan)=err_vals_out';
 
     opts_simul.SHOCKS = err_vals_out;
 
-    [ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_] = occbin.match_function(...
+    [ resids(this_period,inan), ~, stateval(this_period,:), Emat(:,inan,this_period), M_, out] = occbin.match_function(...
         err_vals_out,obs_list,current_obs,opts_simul, M_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state,options_);
     init_val = stateval(this_period,:); %update
+    regime_history(this_period) = out.regime_history(1);
     if max(abs(err_vals_out))>1e8
         error_code(1) = 306;
         error_code(4) = max(abs(err_vals_out))/1000;
         filtered_errs=NaN;
-        if options_.occbin.likelihood.waitbar; dyn_waitbar_close(hh_fig); end
+        if options_.occbin.likelihood.waitbar; waitbar.close(hh_fig); end
         return
     end
     if max(abs(resids(this_period,:)))>0.001
@@ -115,12 +119,12 @@ for this_period=1:sample_length
         filtered_errs=NaN;
         error_code(1) = 303;
         error_code(4) = max(abs(resids(this_period,:)))*100;
-        if options_.occbin.likelihood.waitbar; dyn_waitbar_close(hh_fig); end
+        if options_.occbin.likelihood.waitbar; waitbar.close(hh_fig); end
         return
     end
 end
 if options_.occbin.likelihood.waitbar
-    dyn_waitbar_close(hh_fig);
+    waitbar.close(hh_fig);
 end
 
 end
\ No newline at end of file

matlab/+occbin/IVF_posterior.m

-function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
+function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov, regime_history] = IVF_posterior(xparam1,...
     dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
 % [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,bayestopt_,dr, atT, innov] = IVF_posterior(xparam1,...
 %     dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,BoundsInfo,dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
@@ -7,10 +7,10 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,baye
 % INPUTS
 % - xparam1             [double]        current values for the estimated parameters.
 % - dataset_            [structure]     dataset after transformations
-% - dataset_info        [structure]     storing informations about the
+% - dataset_info        [structure]     storing information about the
 %                                       sample; not used but required for interface
-% - options_            [structure]     Matlab's structure describing the current options
-% - M_                  [structure]     Matlab's structure describing the model
+% - options_            [structure]     MATLAB's structure describing the current options
+% - M_                  [structure]     MATLAB's structure describing the model
 % - estim_params_       [structure]     characterizing parameters to be estimated
 % - bayestopt_          [structure]     describing the priors
 % - BoundsInfo          [structure]     containing prior bounds
@@ -21,7 +21,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,M_,options_,baye
 %
 % OUTPUTS
 % - fval                    [double]        scalar, value of the likelihood or posterior kernel.
-% - info                    [integer]       4×1 vector, informations resolution of the model and evaluation of the likelihood.
+% - info                    [integer]       4×1 vector, information resolution of the model and evaluation of the likelihood.
 % - exit_flag               [integer]       scalar, equal to 1 (no issues when evaluating the likelihood) or 0 (not able to evaluate the likelihood).
 % - DLIK                    [double]        Empty array.
 % - Hess                    [double]        Empty array.
@@ -58,6 +58,8 @@ trend_coeff = [];
 obs = dataset_.data;
 obs_list = options_.varobs(:);
 exit_flag   = 1;
+% initialize output argument in case of early return
+regime_history = [];
 
 
 if size(xparam1,1)<size(xparam1,2)
@@ -90,7 +92,11 @@ if info(1)
                 info(1) == 81 || info(1) == 84 ||  info(1) == 85 ||  info(1) == 86
         %meaningful second entry of output that can be used
         fval = Inf;
+        if ~isfinite(info(2))
+            info(4) = 0.1;
+        else
             info(4) = info(2);
+        end
         exit_flag = 0;
         return
     else
@@ -107,7 +113,7 @@ end
 sample_length = size(obs,1);
 filtered_errs_init = zeros(sample_length,sum(err_index));
 
-[filtered_errs, resids, Emat, stateval, info] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
+[filtered_errs, resids, Emat, stateval, info, regime_history] = occbin.IVF_core(M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,options_,err_index,filtered_errs_init,obs_list,obs);
 if info(1)
     fval = Inf;
     exit_flag = 0;
@@ -153,26 +159,20 @@ end
 like = 0.5*sum(likei(options_.presample+1:end));
 
 if isinf(like) 
-    fval = Inf;
-    info(1) = 301;
-    info(4) = 1000;
-    exit_flag = 0;
+    fval = Inf; info(1) = 301; info(4) = 1000; exit_flag = 0;
     return
 elseif isnan(like)
-    fval = Inf;
-    info(1) = 302;
-    info(4) = 1000;
-    exit_flag = 0;
+    fval = Inf; info(1) = 302; info(4) = 1000; exit_flag = 0;
+    return
+elseif imag(like)~=0
+    fval = Inf; info(1) = 300; info(4) = 1000; exit_flag = 0;
     return
 end
 
 maxresid = max(abs(resids(:)));
 if maxresid>1e-3
     disp_verbose('Penalize failure of residuals to be zero',options_.verbosity)    
-    fval = Inf;
-    info(1) = 303;
-    info(4) = sum(resids(:).^2);
-    exit_flag = 0;
+    fval = Inf; info(1) = 303; info(4) = sum(resids(:).^2); exit_flag = 0;
     return
 end
 
@@ -181,11 +181,17 @@ if ~isempty(xparam1)
 else
     prior = 0;
 end
-if prior == Inf
+if isinf(prior)
     % If parameters outside prior bound, minus prior density is very large
-    fval = Inf;
-    info(4) = 1000;
-    exit_flag = 0;
+    fval = Inf; info(1) = 40; info(4) = 1000; exit_flag = 0;
+    return
+end
+if isnan(prior)
+    fval = Inf; info(1) = 47; info(4) = 1000; exit_flag = 0;
+    return
+end
+if imag(prior)~=0
+    fval = Inf; info(1) = 48; info(4) = 1000; exit_flag = 0;
     return
 end
 

matlab/+occbin/check_regimes.m

@@ -6,8 +6,8 @@ function [TT, RR, CC, regime_history] = check_regimes(TT, RR, CC, opts_regime, M
 % - RR            [N by N_exo]      shock impact matrix of state space
 % - CC            [N by 1]          constant of state space
 % - opts_regime_  [structure]       structure describing the regime
-% - M_            [structure]       Matlab's structure describing the model
-% - options_      [structure]       Matlab's structure describing the current options
+% - M_            [structure]       MATLAB's structure describing the model
+% - options_      [structure]       MATLAB's structure describing the current options
 % - dr                  [structure]     Reduced form model.
 % - endo_steady_state   [vector]        steady state value for endogenous variables
 % - exo_steady_state    [vector]        steady state value for exogenous variables

matlab/+occbin/cost_function.m

@@ -4,7 +4,7 @@ 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_)
 % Outputs:
 %  - cost               [double]        penalty 
-%  - out                [structure]     Occbin's results structure
+%  - out                [structure]     OccBin's results structure
 %
 % Inputs
 % - err_0               [double]        value of shocks 
@@ -12,12 +12,12 @@ function [cost, out] = cost_function(err_0, current_obs, weights, opts_simul,...
 % - 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_).
+% - 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_).
+% - options_            [structure]     MATLAB's structure describing the current options (options_).
 
 % Copyright © 2023 Dynare Team
 %

matlab/+occbin/dynare_resolve.m

@@ -6,8 +6,8 @@ function [A,B,ys,info,dr,params,TT, RR, CC, A0, B0] ...
 % transition equation. Mirrors dynare_resolve
 %
 % Inputs:
-% - M_                  [structure]     Matlab's structure describing the model
-% - 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]     Reduced form model.
 % - endo_steady_state   [vector]        steady state value for endogenous variables
 % - exo_steady_state    [vector]        steady state value for exogenous variables

matlab/+occbin/findmin.m

@@ -2,17 +2,17 @@ function [y, out, cost] = findmin(d_index, a0, P1, Qt, Y, ZZ, opts_simul,M_, dr,
 % [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
+%  - 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_).
+% - M_                  [structure]     MATLAB's structure describing the model
 % - 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_).
+% - options_            [structure]     MATLAB's structure describing the current options
 
 % Copyright © 2023 Dynare Team
 %

matlab/+occbin/forecast.m

@@ -3,8 +3,8 @@ function [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_
 % Occbin forecasts
 %
 % INPUTS
-% - options_                [structure]     Matlab's structure describing the current options
-% - M_                      [structure]     Matlab's structure describing the model
+% - options_                [structure]     MATLAB's structure describing the current options
+% - M_                      [structure]     MATLAB's structure describing the model
 % - dr_in                   [structure]     model information structure
 % - endo_steady_state       [double]        steady state value for endogenous variables
 % - exo_steady_state        [double]        steady state value for exogenous variables
@@ -18,7 +18,7 @@ function [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_
 % SPECIAL REQUIREMENTS
 %   none.
 
-% Copyright © 2022-2023 Dynare Team
+% Copyright © 2022-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -79,7 +79,8 @@ if ~isempty(shocks_input)
 end
 
 if opts.replic
-    h = dyn_waitbar(0,'Please wait occbin forecast replic ...');
+    options_.noprint=true;
+    h = waitbar.run(0,'Please wait occbin forecast replic ...');
     ishock = find(sqrt(diag((M_.Sigma_e))));
     options_.occbin.simul.exo_pos=ishock;
     effective_exo_nbr=  length(ishock);
@@ -98,7 +99,7 @@ if opts.replic
     error_flag=true(opts.replic,1);
     simul_SHOCKS=NaN(forecast_horizon,M_.exo_nbr,opts.replic);
     for iter=1:opts.replic
-        options_.occbin.simul.SHOCKS = shocks_base+transpose(U*sqrt(S)*SHOCKS_add(:,:,iter));
+        options_.occbin.simul.SHOCKS = shocks_base(:,ishock)+transpose(U*sqrt(S)*SHOCKS_add(:,:,iter));
         options_.occbin.simul.waitbar=0;
         [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
         error_flag(iter)=out.error_flag;
@@ -113,13 +114,16 @@ if opts.replic
                 save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag','out','shocks_base')
             end
         end
-        dyn_waitbar(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
+        waitbar.run(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
     end
-    dyn_waitbar_close(h);
+    waitbar.close(h);
     if options_.debug
          save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag')
     end
     inx=find(error_flag==0);
+    if length(inx)<opts.replic
+        fprintf('\noccbin.forecast: forecast simulation was only successful in %u of %u simulations.\n',length(inx),opts.replic);
+    end
     z.linear=z.linear(:,:,inx);
     z.piecewise=z.piecewise(:,:,inx);
     z.min.piecewise = min(z.piecewise,[],3);

matlab/+occbin/get_deriv.m

@@ -13,7 +13,7 @@ function [h_minus_1, h, h_plus_1, h_exo, resid] = get_deriv(M_, ys_)
 % - h_exo      [N by N_exo] derivative matrix with respect to exogenous variables
 % - resid      [N by 1]     vector of residuals
 
-% Copyright © 2021 Dynare Team
+% Copyright © 2021-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,40 +30,16 @@ function [h_minus_1, h, h_plus_1, h_exo, resid] = get_deriv(M_, ys_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-x = zeros(M_.maximum_lag + M_.maximum_lead + 1,M_.exo_nbr);
+dyn_endo_ss = repmat(ys_, 3, 1);
+x = zeros(M_.exo_nbr);
 
-iyv = M_.lead_lag_incidence';
-iyr0 = find(iyv(:)) ;
-z=repmat(ys_,1,M_.maximum_lag + M_.maximum_lead + 1);
+[resid, T_order, T] = feval([M_.fname '.sparse.dynamic_resid'], dyn_endo_ss, x, M_.params, ys_);
 
-[resid,g1]=feval([M_.fname,'.dynamic'],z(iyr0),x, M_.params, ys_, M_.maximum_exo_lag+1);
+g1 = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, x, M_.params, ys_, ...
+           M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
+           M_.dynamic_g1_sparse_colptr, T_order, T);
 
-% Initialize matrices
-h_minus_1=zeros(M_.endo_nbr);
-h = h_minus_1;
-h_plus_1 = h_minus_1;
-
-% build h_minus_1
-if M_.maximum_lag
-    lag_columns=find(iyv(:,1));
-    n_lag_columns=length(lag_columns);
-    h_minus_1(:,lag_columns) = g1(:,1:n_lag_columns);
-else
-    n_lag_columns=0;
-end
-
-% build h
-contemporaneous_columns=find(iyv(:,M_.maximum_lag+1));
-n_contemporaneous_columns = length(contemporaneous_columns);
-h(:,contemporaneous_columns) = g1(:,1+n_lag_columns:n_lag_columns+n_contemporaneous_columns);
-
-%build h_plus_1
-if M_.maximum_lead
-    lead_columns=find(iyv(:,end));
-    n_lead_columns = length(lead_columns);
-    h_plus_1(:,lead_columns) = g1(:,n_lag_columns+n_contemporaneous_columns+1:n_lag_columns+n_contemporaneous_columns+n_lead_columns);
-else
-    n_lead_columns=0;
-end
-
-h_exo =g1(:,n_lag_columns+n_contemporaneous_columns+n_lead_columns+1:end);
+h_minus_1 = full(g1(:,1:M_.endo_nbr));
+h = full(g1(:,M_.endo_nbr + (1:M_.endo_nbr)));
+h_plus_1 = full(g1(:,2*M_.endo_nbr + (1:M_.endo_nbr)));
+h_exo = full(g1(:,3*M_.endo_nbr+1:end));

matlab/+occbin/graph.m

@@ -2,10 +2,10 @@ function graph(M_, options_, options_occbin_, oo_, var_list)
 % function graph(M_, options_, options_occbin_, oo_, var_list)
 %
 % Inputs:
-% - M_                  [structure]     Matlab's structure describing the model
-% - options_            [structure]     Matlab's structure containing the options
-% - options_occbin_     [structure]     Matlab's structure containing Occbin options
-% - oo_                 [structure]     Matlab's structure containing the results
+% - M_                  [structure]     MATLAB's structure describing the model
+% - options_            [structure]     MATLAB's structure containing the options
+% - options_occbin_     [structure]     MATLAB's structure containing OccBin options
+% - oo_                 [structure]     MATLAB's structure containing the results
 % - var_list            [char]          list of the variables to plot 
 
 % Copyright © 2021-2023 Dynare Team
@@ -69,7 +69,7 @@ if number_of_plots_to_draw_exo>0
         if ~isempty(temp_index)
             exo_index(ii)=temp_index;
         else
-            error('%s was not part of the shocks for Occbin.', var_list{i_var_exo(ii)});
+            error('%s was not part of the shocks for OccBin.', var_list{i_var_exo(ii)});
         end
     end
     data_to_plot(:,end+1:end+number_of_plots_to_draw_exo,1)=[oo_.occbin.simul.shocks_sequence(:,exo_index); zeros(nperiods-size(oo_.occbin.simul.shocks_sequence,1),number_of_plots_to_draw_exo)];
@@ -79,7 +79,7 @@ end
 [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw_endo+number_of_plots_to_draw_exo);
 
 for fig = 1:nbplt
-    hh_fig = dyn_figure(options_.nodisplay,'Name',['Occbin simulated paths, figure ' int2str(fig)]);
+    hh_fig = dyn_figure(options_.nodisplay,'Name',['OccBin simulated paths, figure ' int2str(fig)]);
     for plt = 1:nstar
         if fig==nbplt && ~lr==0
             subplot(lr,lc,plt);

matlab/+occbin/irf.m

@@ -3,9 +3,9 @@ function irfs = irf(M_,oo_,options_)
 % Calls a minimizer
 %
 % INPUTS
-% - M_                  [structure]     Matlab's structure describing the model
-% - oo_                 [structure]     Matlab's structure containing the results
-% - options_            [structure]     Matlab's structure describing the current options
+% - M_                  [structure]     MATLAB's structure describing the model
+% - oo_                 [structure]     MATLAB's structure containing the results
+% - options_            [structure]     MATLAB's structure describing the current options
 %
 % OUTPUTS
 % - irfs                [structure]     IRF results

matlab/+occbin/kalman_update_algo_1.m

@@ -17,13 +17,13 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat, alpha
 % - RR              [N by N_exo by 2]       shock impact matrix at t-1:t
 % - CC              [N by 2]                state space constant state transition matrix at t-1:t
 % - regimes0        [structure]             regime info at t-1:t
-% - M_              [structure]             Matlab's structure describing the model (M_).
+% - M_              [structure]             MATLAB's structure describing the model
 % - dr                   [structure]        Reduced form model.
 % - 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_).
-% - occbin_options_ [structure]             Matlab's structure describing the Occbin options.
+% - options_        [structure]             MATLAB's structure describing the current options
+% - occbin_options_ [structure]             MATLAB's structure describing the OccBin options.
 % - kalman_tol      [double]                tolerance for reciprocal condition number
 % 
 % Outputs
@@ -37,7 +37,7 @@ function [a, a1, P, P1, v, T, R, C, regimes_, error_flag, M_, lik, etahat, alpha
 % - C               [N by 2]                state space constant state transition matrix at t-1:t
 % - regimes_        [structure]             regime info at t-1:t
 % - error_flag      [integer]               error code
-% - M_              [structure]             Matlab's structure describing the model (M_).
+% - M_              [structure]             MATLAB's structure describing the model
 % - lik             [double]                likelihood
 % - etahat:                                 smoothed shocks
 %
@@ -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