matlab/+occbin/kalman_update_algo_3.m

@@ -20,13 +20,13 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, lik, 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_).
-% - options_        [structure]             Matlab's structure describing the current options (options_).
+% - M_              [structure]             MATLAB's structure describing the model (M_).
+% - options_        [structure]             MATLAB's structure describing the current options (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
 % - exo_det_steady_state [vector]           steady state value for exogenous deterministic variables
-% - occbin_options_ [structure]             Matlab's structure describing the Occbin options.
+% - occbin_options_ [structure]             MATLAB's structure describing the OccBin options.
 % - kalman_tol      [double]                tolerance for reciprocal condition number
 % - nk              [double]                number of forecasting periods
 % 
@@ -43,7 +43,7 @@ function [a, a1, P, P1, v, Fi, Ki, T, R, C, regimes_, error_flag, M_, lik, alpha
 % - CC              [N by 2]                state space constant state transition matrix at t-1:t
 % - regimes_        [structure]             regime info at t-1:t
 % - error_flag      [structure]             error flag
-% - M_              [structure]             Matlab's structure describing the model (M_).
+% - M_              [structure]             MATLAB's structure describing the model (M_).
 % - alphahat:                               smoothed variables (a_{t|T})
 % - etahat:                                 smoothed shocks
 % - TT              [N by N by 2]           state transition matrix at t-1:t
@@ -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

@@ -32,13 +32,14 @@ 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
+likvec = likx;
+regvec = regx(1);
 info0=info;
 if info
     if ~isequal(regimes_(1:2),[base_regime base_regime])
@@ -48,6 +49,8 @@ if info
             [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
+    likvec = likx;
+    regvec = regx(1);
     info1=info;
 else
     if ~isequal(regimes_(1:2),[base_regime base_regime])
@@ -56,6 +59,10 @@ else
         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 && not(isequal(regx1,regx))
+            likvec = [likvec likx1];
+            regvec = [regvec; regx1(1)];
+        end
         if info1==0 && likx1<likx
             ax=ax1;
             a1x=a1x1;
@@ -101,23 +108,88 @@ if M_.occbin.constraint_nbr==1
         oldstart = regimes_(1).regimestart(end);
         newstart = regx(1).regimestart(end);
         diffstart = newstart-oldstart;
+        regname = 'regimestart';
+        reg_string = 'regime';
     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);
+        [diffstart, diffregime] = max([newstart1-oldstart1,newstart2-oldstart2]);
+        switch diffregime
+            case 1
+                regname = 'regimestart1';
+                reg_string = 'regime1';
+            case 2
+                regname = 'regimestart2';
+                reg_string = 'regime2';
+        end
     end
 end
+if options_.occbin.filter.use_relaxation && diffstart>options_.occbin.filter.use_relaxation
+    guess_regime = [base_regime base_regime];
+    options_.occbin.filter.guess_regime = true;
+    guess_regime(1) = regx(1);
+    regx2= regx;
+    while isequal(guess_regime(1),regx2(1))
+        % we reduce length until the converged regime does not change
 
-if options_.occbin.filter.use_relaxation && diffstart>2
-    if info0==0
-        % make sure we match criteria to enter further solution attempts
-        info1=1;
+        guess_regime(1).(regname)(end) = regx2(1).(regname)(end)-1;
+        if guess_regime(1).(regname)(end)==1
+            % make sure we enforce base regime
+            guess_regime(1).(regname)=guess_regime(1).(regname)(end);
+            guess_regime(1).(reg_string)=0;
+        else
+            if guess_regime(1).(regname)(end-1)==guess_regime(1).(regname)(end)
+                guess_regime(1).(regname)(end-1) = guess_regime(1).(regname)(end-1)-1;
             end
-    options_.occbin.likelihood.brute_force_regime_guess   = true;
-    options_.occbin.likelihood.loss_function_regime_guess = true;
-    use_relaxation = true;
+        end        
+        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
+        isnew=true;
+        for kr=1:length(likvec)
+            % make sure likelihood does not differ by rounding issue
+            % but due to different regimes
+            if isequal(regx2(1),regvec(kr))
+                isnew = false;
+            end
+        end
+        if isnew
+            likvec = [likvec likx2];
+            regvec = [regvec; regx2(1)];
+        end
+        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
+    options_.occbin.filter.guess_regime = false;
+    options_.occbin.likelihood.brute_force_regime_guess   = false;
+    options_.occbin.likelihood.loss_function_regime_guess = false;
 end
 
 if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (info==0 &&  ~isequal(regx(1),base_regime))
@@ -137,9 +209,8 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
             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
+                if not(info==0 && isequal(regx2{1},regx)) 
+                    % found a solution, different from previous 
                     break
                 end
             end
@@ -153,8 +224,6 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
             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
@@ -216,18 +285,16 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
                             end
                             if info2==0
                                 use_index= gindex;
-                                if not(info==0 && isequal(regx2{gindex},regx)) && not(use_relaxation && likx2{gindex}>=likx)
+                                if not(info==0 && isequal(regx2{gindex},regx)) 
                                     % 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)
+                            if not(info==0 && isequal(regx2{gindex},regx)) 
                                 % found a solution, different from previous one
-                                % use_relaxation and likelihood improves
                                 break
                             end
                         end
@@ -235,9 +302,8 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
                     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)
+                        if not(info==0 && isequal(regx2{gindex},regx)) 
                             % found a solution, different from previous one
-                            % use_relaxation and likelihood improves
                             break
                         end
                     end
@@ -261,6 +327,18 @@ if options_.occbin.likelihood.brute_force_regime_guess && (info0 || info1) %|| (
         % so that we DO NOT enter IVF step
         info0=0;
         info1=0;
+        isnew=true;
+        for kr=1:length(likvec)
+            % make sure likelihood does not differ by rounding issue
+            % but truly for different regimes
+            if isequal(regx2{use_index}(1),regvec(kr))
+                isnew = false;
+            end
+        end
+        if isnew
+            likvec = [likvec likx2{use_index}];
+            regvec = [regvec; regx2{use_index}(1)];
+        end
     end
     if info2==0 && likx2{use_index}<likx
         ax=ax2{use_index};
@@ -301,6 +379,18 @@ if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %||
             [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;
+        isnew=true;
+        for kr=1:length(likvec)
+            % make sure likelihood does not differ by rounding issue
+            % but due to different regimes
+            if isequal(regx2(1),regvec(kr))
+                isnew = false;
+            end
+        end
+        if isnew
+            likvec = [likvec likx2];
+            regvec = [regvec; regx2(1)];
+        end
         if info2==0 && likx2<likx
             ax=ax2;
             a1x=a1x2;
@@ -328,4 +418,17 @@ if options_.occbin.likelihood.loss_function_regime_guess && (info0 || info1) %||
         end
     end
 end
+if length(likvec)>1
+    % sum the likelihood of multiple solutions
+    likx = -2*log(sum(exp(-likvec./2)));
+end
+
+if info(1)==0
+    if isnan(likx)
+        info = 323;
+    elseif any(any(isnan(ax)))
+        info = 324;
+    end
+end
+
 end

matlab/+occbin/match_function.m

@@ -7,20 +7,20 @@ function [resids, grad, state_out, E, M_, out] = match_function(err_0, obs_list,
 %  - grad           [double]        [n by n_exo] gradient (response of observables to shocks)
 %  - state_out      [double]        [ny by 1] value of endogenous variables
 %  - E              [double]        response of endogenous variables to shocks
-%  - M_             [structure]     Matlab's structure describing the model (M_).
-%  - out            [structure]     Occbin's results structure
+%  - M_             [structure]     MATLAB's structure describing the model (M_).
+%  - out            [structure]     OccBin's results structure
 %
 % Inputs
 % - err_            [double]        value of shocks 
 % - obs_list        [cell]          names of observables
 % - current_obs     [double]        [1 by n_obs] current value of observables
 % - opts_simul      [structure]     Structure with simulation options
-% - 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_).
 
 % Original authors: Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong
 % Original file downloaded from:

matlab/+occbin/plot_irfs.m

@@ -2,9 +2,9 @@ function plot_irfs(M_,irfs,options_,var_list)
 % plot_irfs(M_,irfs,options_,var_list)
 %
 % INPUTS
-% - M_                      [structure]     Matlab's structure describing the model
+% - M_                      [structure]     MATLAB's structure describing the model
 % - irfs                    [structure]     IRF results
-% - options_                [structure]     Matlab's structure describing the current options
+% - options_                [structure]     MATLAB's structure describing the current options
 % - var_list                [character array]  list of endogenous variables specified
 %
 % OUTPUTS

matlab/+occbin/plot_regimes.m

@@ -2,8 +2,8 @@ function plot_regimes(regimes,M_,options_)
 % plot_regimes(regimes,M_,options_)
 % Inputs: 
 % - regimes     	[structure]     OccBin regime information
-% - 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
 
 % Copyright © 2021-2023 Dynare Team
 %

matlab/+occbin/set_default_options.m

 function options_occbin_ = set_default_options(options_occbin_,M_,flag)
 % function options_occbin_ = set_default_options(options_occbin_,M_,flag)
-% Sets default options for Occbin
+% Sets default options for OccBin
 %
 % INPUTS
-% - options_occbin_ [structure]     Matlab's structure describing the current options
-% - M_              [structure]     Matlab's structure describing the model
-% - flag            [cell]          govern what/how much to initialize
+% - options_occbin_ [structure]     MATLAB's structure describing the current options
+% - M_              [structure]     MATLAB's structure describing the model
+% - flag            [cell]          governs what/how much to initialize
 %
 % OUTPUTS
-% - options_occbin_ [structure]     Matlab's structure describing the current options
+% - options_occbin_ [structure]     MATLAB's structure describing the current options
 
 % Copyright © 2021 Dynare Team
 %
@@ -46,6 +46,7 @@ if ismember(flag,{'filter','all'})
     options_occbin_.filter.guess_regime = false;
     options_occbin_.filter.periodic_solution = true;
     options_occbin_.filter.use_relaxation = false;
+    options_occbin_.filter.use_relaxation_tol_period = 1;
 end
 
 if ismember(flag,{'forecast','all'})

matlab/+occbin/set_option.m

 function options_=set_option(options_,options_occbin_,fieldname)
 % function options_=set_option(options_,options_occbin_,fieldname)
-% Set local option for Occbin
+% Set local option for OccBin
 %
 % Inputs:
-% - options_            [structure]     Matlab's structure containing the options
-% - options_occbin_     [structure]     Matlab's structure containing Occbin options
+% - options_            [structure]     MATLAB's structure containing the options
+% - options_occbin_     [structure]     MATLAB's structure containing OccBin options
 % - fieldname           [string]        name of the options field to set 
 %
 % Outputs:
-% - options_            [structure]     Matlab's structure containing the options
+% - options_            [structure]     MATLAB's structure containing the options
 
 % Copyright © 2021 Dynare Team
 %

matlab/+occbin/setup.m

 function [M_, options_] = setup(M_,options_, options_occbin_)
 % function [M_, options_] = setup(M_, options_, options_occbin_)
-% Sets up run of Occbin: creates shock matrix, sets options
+% Sets up run of OccBin: creates shock matrix, sets options
 %
 % INPUT:
-% - 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
+% - 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
 %
 % OUTPUT:
-% - M_                  [structure]     Matlab's structure describing the model
-% - options_occbin_     [structure]     Matlab's structure containing Occbin options
+% - M_                  [structure]     MATLAB's structure describing the model
+% - options_occbin_     [structure]     MATLAB's structure containing OccBin options
 
 % Copyright © 2021 Dynare Team
 %

matlab/+occbin/shock_decomposition.m

@@ -2,12 +2,12 @@ function oo_ = shock_decomposition(oo_, M_, options_, vname)
 % function oo_ = shock_decomposition(oo_, M_, options_, vname)
 %
 % INPUTS
-% - oo_           [structure]     Matlab's structure containing the results (oo_).
-% - M_            [structure]     Matlab's structure describing the model (M_).
-% - options_      [structure]     Matlab's structure describing the current options (options_).
+% - oo_           [structure]     MATLAB's structure containing the results
+% - M_            [structure]     MATLAB's structure describing the model
+% - options_      [structure]     MATLAB's structure describing the current options
 % - vname         [cell]          array of variable names
 % OUTPUT
-% - oo_           [structure]     Matlab's structure containing the results (oo_).
+% - oo_           [structure]     MATLAB's structure containing the results
 
 % Copyright © 2021 Dynare Team
 %
@@ -45,7 +45,7 @@ use_shock_groups = shock_decomp_options.use_shock_groups;
 init_names_ = shock_decomp_options.init_names_;
 nfrcst = shock_decomp_options.nfrcst;
 
-%% new dynare grouping
+%% new Dynare grouping
 if isempty(use_shock_groups)
     use_shock_groups = 'ALL';
     ngroups = M_.exo_nbr;

matlab/+occbin/solve_one_constraint.m

@@ -2,9 +2,9 @@ function [data, SS_out, error_flag ] = solve_one_constraint(M_,dr, opts_simul_,
 %function [data, SS_out, error_flag ] = solve_one_constraint(M_,dr, opts_simul_, solve_DM)
 %
 % INPUT: 
-% - M_                  [structure]     Matlab's structure describing the model (M_).
+% - M_                  [structure]     MATLAB's structure describing the model (M_).
 % - dr                  [structure]     decision rules for the model
-% - opts_simul          [structure]     Matlab's structure containing the Occbin options (opts_simul).
+% - opts_simul          [structure]     MATLAB's structure containing the OccBin options (opts_simul).
 % - solve_DM            [double]        indicator on whether to recompute decision rules
 %
 % OUTPUT:
@@ -17,7 +17,7 @@ function [data, SS_out, error_flag ] = solve_one_constraint(M_,dr, opts_simul_,
 %                                           - 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
-% - error_flag          [integer]       1 if a problem was encoutered, 0 otherwise
+% - error_flag          [integer]       1 if a problem was encountered, 0 otherwise
 
 % Original authors: Luca Guerrieri and Matteo Iacoviello 
 % Original file downloaded from:
@@ -25,7 +25,7 @@ function [data, SS_out, error_flag ] = solve_one_constraint(M_,dr, opts_simul_,
 % 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:
 %
 % Luca Guerrieri and Matteo Iacoviello (2015): "OccBin: A toolkit for solving
@@ -114,13 +114,13 @@ else
 end
 
 if opts_simul_.waitbar
-    hh_fig = dyn_waitbar(0,'Occbin: Solving the model');
-    set(hh_fig,'Name','Occbin: Solving the model.');
+    hh_fig = waitbar.run(0,'OccBin: Solving the model');
+    set(hh_fig,'Name','OccBin: Solving the model.');
 end
 
 for shock_period = 1:n_shocks_periods
     if opts_simul_.waitbar
-        dyn_waitbar(shock_period/n_shocks_periods, hh_fig, sprintf('Period %u of %u', shock_period,n_shocks_periods));
+        waitbar.run(shock_period/n_shocks_periods, hh_fig, sprintf('Period %u of %u', shock_period,n_shocks_periods));
     end
     
     regime_change_this_iteration=true;
@@ -330,11 +330,11 @@ for shock_period = 1:n_shocks_periods
         if max_iter>opts_simul_.algo_truncation
             disp_verbose(['occbin solver:: period ' int2str(shock_period) '::'],opts_simul_.debug)
             if is_periodic
-                disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
+                disp_verbose('OccBin solver loops between two regimes.',opts_simul_.debug)
                 if periodic_solution
                     disp_verbose(['Max error:' num2str(merr) '.'],opts_simul_.debug)
                 else
-                    if opts_simul_.waitbar; dyn_waitbar_close(hh_fig); end
+                    if opts_simul_.waitbar; waitbar.close(hh_fig); end
                     error_flag = 310;
                     return
                 end
@@ -346,7 +346,7 @@ for shock_period = 1:n_shocks_periods
                 disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
                     error_flag = 311;
                 end
-                if opts_simul_.waitbar; dyn_waitbar_close(hh_fig); end
+                if opts_simul_.waitbar; waitbar.close(hh_fig); end
                 return
             end
         else
@@ -355,7 +355,7 @@ for shock_period = 1:n_shocks_periods
     end
     if any(error_code_period)
         disp_verbose('Increase nperiods.',opts_simul_.debug)
-        if opts_simul_.waitbar; dyn_waitbar_close(hh_fig); end
+        if opts_simul_.waitbar; waitbar.close(hh_fig); end
         error_flag = 312;
         return
     end
@@ -385,5 +385,5 @@ if ~opts_simul_.piecewise_only
 end
 
 if opts_simul_.waitbar
-    dyn_waitbar_close(hh_fig); 
+    waitbar.close(hh_fig); 
 end

matlab/+occbin/solve_two_constraints.m

@@ -2,9 +2,9 @@ function [ data, SS_out, error_flag] = solve_two_constraints(M_,dr, opts_simul_,
 % function [ data, SS_out, error_flag] = solve_two_constraints(M_,dr, opts_simul_, solve_DM)
 %
 % INPUT: 
-% - M_                  [structure]     Matlab's structure describing the model (M_).
+% - M_                  [structure]     MATLAB's structure describing the model
 % - dr                  [structure]     decision rules for the model
-% - opts_simul          [structure]     Matlab's structure containing the Occbin options (opts_simul).
+% - opts_simul          [structure]     MATLAB's structure containing the OccBin options (opts_simul).
 % - solve_DM            [double]        indicator on whether to recompute decision rules
 %
 % OUTPUT:
@@ -17,7 +17,7 @@ function [ data, SS_out, error_flag] = solve_two_constraints(M_,dr, opts_simul_,
 %                                           - 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
-% - error_flag          [integer]       1 if a problem was encoutered, 0 otherwise
+% - error_flag          [integer]       1 if a problem was encountered, 0 otherwise
 
 % Original authors: Luca Guerrieri and Matteo Iacoviello 
 % Original file downloaded from:
@@ -25,7 +25,7 @@ function [ data, SS_out, error_flag] = solve_two_constraints(M_,dr, opts_simul_,
 % 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:
 %
 % Luca Guerrieri and Matteo Iacoviello (2015): "OccBin: A toolkit for solving
@@ -124,13 +124,13 @@ else
 end
 
 if opts_simul_.waitbar
-    hh_fig = dyn_waitbar(0,'Occbin: Solving the model');
-    set(hh_fig,'Name','Occbin: Solving the model.');
+    hh_fig = waitbar.run(0,'OccBin: Solving the model');
+    set(hh_fig,'Name','OccBin: Solving the model.');
 end
 
 for shock_period = 1:n_shocks_periods
     if opts_simul_.waitbar
-        dyn_waitbar(shock_period/n_shocks_periods, hh_fig, sprintf('Period %u of %u', shock_period,n_shocks_periods));
+        waitbar.run(shock_period/n_shocks_periods, hh_fig, sprintf('Period %u of %u', shock_period,n_shocks_periods));
     end
     regime_change_this_iteration=true;
     nperiods_endogenously_increased = false;
@@ -363,12 +363,12 @@ for shock_period = 1:n_shocks_periods
         if max_iter>opts_simul_.algo_truncation
             disp_verbose(['occbin solver: period ' int2str(shock_period) ':'],opts_simul_.debug)
             if is_periodic
-                disp_verbose('Occbin solver loops between two regimes.',opts_simul_.debug)
+                disp_verbose('OccBin solver loops between two regimes.',opts_simul_.debug)
                 if periodic_solution
                     disp_verbose(['Max error:' num2str(min_err) '.'],opts_simul_.debug)
                 else
                     error_flag = 310;
-                    if opts_simul_.waitbar; dyn_waitbar_close(hh_fig); end
+                    if opts_simul_.waitbar; waitbar.close(hh_fig); end
                     return;
                 end
             else
@@ -379,7 +379,7 @@ for shock_period = 1:n_shocks_periods
                     disp_verbose('Did not converge -- increase maxit.',opts_simul_.debug)
                     error_flag = 311;
                 end
-                if opts_simul_.waitbar; dyn_waitbar_close(hh_fig); end
+                if opts_simul_.waitbar; waitbar.close(hh_fig); end
                 return;
             end
         else
@@ -391,7 +391,7 @@ for shock_period = 1:n_shocks_periods
     if any(error_code_period)
         disp_verbose('Increase nperiods.',opts_simul_.debug)
         error_flag = 312;
-        if opts_simul_.waitbar; dyn_waitbar_close(hh_fig); end
+        if opts_simul_.waitbar; waitbar.close(hh_fig); end
         return;
     end
     
@@ -416,5 +416,5 @@ if ~opts_simul_.piecewise_only
 end
 
 if opts_simul_.waitbar
-    dyn_waitbar_close(hh_fig); 
+    waitbar.close(hh_fig); 
 end

matlab/+occbin/solver.m

@@ -3,8 +3,8 @@ function [dr, out, ss] = solver(M_, options_, dr ,steady_state, exo_steady_state
 % Solves the model with an OBC and produces simulations/IRFs
 %
 % INPUT: 
-% - 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]     model information structure
 % - endo_steady_state       [vector]        steady state value for endogenous variables
 % - exo_steady_state        [vector]        steady state value for exogenous variables

matlab/+occbin/unpack_simulations.m

 function oo_=unpack_simulations(M_,oo_,options_)
 % function oo_=unpack_simulations(M_,oo_,options_)
-% Writes Occbin simulations from matrix to structure
+% Writes OccBin simulations from matrix to structure
 % 
 % Inputs
-% - 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
+% - oo_                 [structure]     MATLAB's structure containing the results
+% - options_            [structure]     MATLAB's structure containing the options
 %
 % Outputs
-% - oo_                 [structure]     Matlab's structure containing the results
+% - oo_                 [structure]     MATLAB's structure containing the results
 
 % Copyright © 2021 Dynare Team
 %

matlab/+occbin/write_regimes_to_xls.m

@@ -4,8 +4,8 @@ function write_regimes_to_xls(occbin_struct,M_,options_)
 %
 % INPUTS
 % - occbin_struct   [struct]    occbin structure containing information on the regimes
-% - M_              [struct]    Matlab's structure describing the model
-% - options_        [struct]    Matlab's structure describing the current options
+% - M_              [struct]    MATLAB's structure describing the model
+% - options_        [struct]    MATLAB's structure describing the current options
 
 % Copyright © 2021-2023 Dynare Team
 %

matlab/+osr/objective.m

@@ -50,7 +50,11 @@ if info(1)
             info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
             info(1) == 81 || info(1) == 84 ||  info(1) == 85
         loss = 1e8;
+        if ~isfinite(info(2))
+            info(4) = 0.1;
+        else
             info(4) = info(2);
+        end
         return
     else
         loss = 1e8;
@@ -82,3 +86,16 @@ else
         df(jp,1) = sum(weights(:).*model_moments_params_derivs(:,jp));
     end
 end
+
+if isinf(loss)
+    loss = 1e8; info(1) = 174; info(4) = 0.1;
+    return
+end
+if isnan(loss)
+    loss = 1e8; info(1) = 175; info(4) = 0.1;
+    return
+end
+if imag(loss)~=0
+    loss = 1e8; info(1) = 176; info(4) = 0.1;
+    return
+end
\ No newline at end of file

matlab/+osr/run.m

@@ -26,7 +26,7 @@ function [info, oo_, options_, M_] = run(M_, options_, oo_, var_list, params, i_
 % SPECIAL REQUIREMENTS
 %   none.
 
-% Copyright © 2001-2023 Dynare Team
+% Copyright © 2001-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -140,11 +140,13 @@ else
         error('OSR: OSR with bounds on parameters requires a constrained optimizer, i.e. opt_algo= 1,2 or 9.')
     end
     %%do actual optimization
-    [p, f] = dynare_minimize_objective(str2func('osr.objective'),par_0,options_.osr.opt_algo,options_,M_.osr.param_bounds,M_.param_names(i_params),[],[], M_,oo_,options_,i_params,...
+    [p, f, ~, ~, ~, ~, optimization_info] = ...
+        dynare_minimize_objective(str2func('osr.objective'),par_0,options_.osr.opt_algo,options_,M_.osr.param_bounds,M_.param_names(i_params),[],[], M_,oo_,options_,i_params,...
                                   inv_order_var(i_var),W(i_var,i_var));
 end
 
 osr_res.objective_function = f;
+osr_res.optimization_info = optimization_info;
 M_.params(i_params)=p; %make sure optimal parameters are set (and not the last draw used in csminwel)
 for i=1:length(i_params)
     osr_res.optim_params.(M_.param_names{i_params(i)}) = p(i);

matlab/+pruned_SS/prodmom.m

 %
 % prodmom.m		Date: 4/29/2006
-% This Matlab program computes the product moment of X_{i_1}^{nu_1}X_{i_2}^{nu_2}...X_{i_m}^{nu_m},
+% This MATLAB program computes the product moment of X_{i_1}^{nu_1}X_{i_2}^{nu_2}...X_{i_m}^{nu_m},
 % where X_{i_j} are elements from X ~ N(0_n,V).  
 % V only needs to be positive semidefinite.
 % V: variance-covariance matrix of X

matlab/+var_expectation/print.m

 function print(varexpectationmodelname, withcalibration)
 
-% Prints the exansion of the VAR_EXPECTATION term in files.
+% Prints the expansion of the VAR_EXPECTATION term in files.
 %
 % INPUTS
 % - varepxpectationmodelname       [string]    Name of the expectation model.
@@ -16,7 +16,7 @@ function print(varexpectationmodelname, withcalibration)
 % - {varexpectationmodelname}-expression.inc     which contains the expanded version of the expectation model kind term.
 %
 % These routines are saved under the {modfilename}/model/varexpectationmodel subfolder, and can be used
-% after in another mod file (ie included with the macro directive @#include). A matlab routine is also
+% after in another mod file (ie included with the macro directive @#include). A MATLAB routine is also
 % created for evaluating (dseries) the var-expectations.
 
 % Copyright © 2018-2021 Dynare Team

matlab/dyn_waitbar_close.m → matlab/+waitbar/close.m

-function dyn_waitbar_close(h)
-% h = dyn_waitbar_close(h)
+function close(h)
+% h = close(h)
 % adaptive close waitbar, compatible with
 % octave and when console_mode=1
 
@@ -23,11 +23,11 @@ function dyn_waitbar_close(h)
 global options_
 
 if options_.console_mode
-    clear dyn_waitbar;
+    clear waitbar.run;
     diary on
     fprintf('\n');
 else
     close(h)
 end
 
-clear dyn_waitbar;
+clear waitbar.run;

matlab/+waitbar/multi.m

+function [cancel,figh] = multi( label, varargin )
+%waitbar.multi: add, remove or update an entry on the multi waitbar
+%
+%   waitbar.multi(TAG,VALUE) adds a waitbar for the specified label, or
+%   if it already exists updates the value. TAG must be a string and
+%   VALUE a number between zero and one or the string 'Close' to remove the
+%   entry Setting value equal to 0 or 'Reset' will cause the progress bar
+%   to reset and the time estimate to be re-initialized.
+%
+%   waitbar.multi(TAG,COMMAND,VALUE,...)  or
+%   waitbar.multi(TAG,VALUE,COMMAND,VALUE,...)
+%   passes one or more command/value pairs for changing the named waitbar
+%   entry. Possible commands include:
+%   'Value'       Set the value of the named waitbar entry. The
+%                 corresponding value must be a number between 0 and 1.
+%   'Increment'   Increment the value of the named waitbar entry. The
+%                 corresponding value must be a number between 0 and 1.
+%   'Color'       Change the color of the named waitbar entry. The
+%                 value must be an RGB triple, e.g. [0.1 0.2 0.3], or a
+%                 single-character color name, e.g. 'r', 'b', 'm'.
+%   'Relabel'     Change the label of the named waitbar entry. The
+%                 value must be the new name.
+%   'Reset'       Set the named waitbar entry back to zero and reset its
+%                 timer. No value need be specified.
+%   'CanCancel'   [on|off] should a "cancel" button be shown for this bar
+%                 (default 'off').
+%   'CancelFcn'   Function to call in the event that the user cancels.
+%   'ResetCancel' Reset the "cancelled" flag for an entry (ie. if you
+%                 decide not to cancel).
+%   'Close'       Remove the named waitbar entry.
+%   'Busy'        Puts this waitbar in "busy mode" where a small bar
+%                 bounces back and forth. Return to normal progress display
+%                 using the 'Reset' command.
+%
+%   cancel = waitbar.multi(LABEL,VALUE) also returns whether the user has
+%   clicked the "cancel" button for this entry (true or false). Two
+%   mechanisms are provided for cancelling an entry if the 'CanCancel'
+%   setting is 'on'. The first is just to check the return argument and if
+%   it is true abort the task. The second is to set a 'CancelFcn' that is
+%   called when the user clicks the cancel button, much as is done for
+%   MATLAB's built-in WAITBAR. In either case, you can use the
+%   'ResetCancel' command if you don't want to cancel after all.
+%
+%   waitbar.multi('CLOSEALL') closes the waitbar window.
+%
+%   Example:
+%   waitbar.multi( 'CloseAll' );
+%   waitbar.multi( 'Task 1', 0 );
+%   waitbar.multi( 'Task 2', 0.5, 'Color', 'b' );
+%   waitbar.multi( 'Task 3', 'Busy');
+%   waitbar.multi( 'Task 1', 'Value', 0.1 );
+%   waitbar.multi( 'Task 2', 'Increment', 0.2 );
+%   waitbar.multi( 'Task 3', 'Reset' ); % Disables "busy" mode
+%   waitbar.multi( 'Task 3', 'Value', 0.3 );
+%   waitbar.multi( 'Task 2', 'Close' );
+%   waitbar.multi( 'Task 3', 'Close' );
+%   waitbar.multi( 'Task 1', 'Close' );
+%
+%   Example:
+%   waitbar.multi( 'Task 1', 0, 'CancelFcn', @(a,b) disp( ['Cancel ',a] ) );
+%   for ii=1:100
+%      abort = waitbar.multi( 'Task 1', ii/100 );
+%      if abort
+%         % Here we would normally ask the user if they're sure
+%         break
+%      else
+%         pause( 1 )
+%      end
+%   end
+%   waitbar.multi( 'Task 1', 'Close' )
+%
+%   Example:
+%   waitbar.multi( 'CloseAll' );
+%   waitbar.multi( 'Red...',    7/7, 'Color', [0.8 0.0 0.1] );
+%   waitbar.multi( 'Orange...', 6/7, 'Color', [1.0 0.4 0.0] );
+%   waitbar.multi( 'Yellow...', 5/7, 'Color', [0.9 0.8 0.2] );
+%   waitbar.multi( 'Green...',  4/7, 'Color', [0.2 0.9 0.3] );
+%   waitbar.multi( 'Blue...',   3/7, 'Color', [0.1 0.5 0.8] );
+%   waitbar.multi( 'Indigo...', 2/7, 'Color', [0.4 0.1 0.5] );
+%   waitbar.multi( 'Violet...', 1/7, 'Color', [0.8 0.4 0.9] );
+
+%   Based on the work by Ben Tordoff
+%
+%   Copyright 2007-2025 The MathWorks, Inc.
+%   All rights reserved.
+% 
+%   Redistribution and use in source and binary forms, with or without
+%   modification, are permitted provided that the following conditions are met:
+% 
+%   * Redistributions of source code must retain the above copyright notice, this
+%     list of conditions and the following disclaimer.
+% 
+%   * Redistributions in binary form must reproduce the above copyright notice,
+%     this list of conditions and the following disclaimer in the documentation
+%     and/or other materials provided with the distribution
+% 
+%   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+%   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+%   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+%   FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+%   DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+%   SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+%   CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+%   OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+%   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+persistent FIGH;
+cancel = false;
+
+% Check basic inputs
+error( nargchk( 1, inf, nargin ) ); %#ok<NCHKN> - kept for backwards compatibility
+if ~ischar( label ) && ~isa( label, 'string' )
+    error( 'multiWaitbar:BadArg', 'LABEL must be the name of the progress entry (i.e. a string)' );
+end
+
+% Try to get hold of the figure
+if isempty( FIGH ) || ~ishandle( FIGH )
+    FIGH = findall( 0, 'Type', 'figure', 'Tag', 'multiWaitbar:Figure' );
+    if isempty(FIGH)
+        FIGH = iCreateFig();
+    else
+        FIGH = handle( FIGH(1) );
+    end
+end
+
+% Check for close all and stop early
+if any( strcmpi( label, {'CLOSEALL','CLOSE ALL'} ) )
+    iDeleteFigure(FIGH);
+    figh = [];
+    return;
+end
+if nargout>1
+    figh = FIGH;
+end
+
+% Make sure we're on-screen
+if ~strcmpi( FIGH.Visible, 'on' )
+    FIGH.Visible = 'on';
+end
+
+% Make sure the timer is still valid - it can be found and deleted
+% externally.
+if ~isvalid( getappdata( FIGH, 'BusyTimer' ) )
+    setappdata( FIGH, 'BusyTimer', iCreateTimer(FIGH) );
+end
+
+% Get the list of entries and see if this one already exists
+entries = getappdata( FIGH, 'ProgressEntries' );
+if isempty(entries)
+    idx = [];
+else
+    idx = find( strcmp( label, {entries.Tag} ), 1, 'first' );
+end
+bgcol = getappdata( FIGH, 'DefaultProgressBarBackgroundColor' );
+
+% If it doesn't exist, create it
+needs_redraw = false;
+entry_added = isempty(idx);
+if entry_added
+    % Create a new entry
+    defbarcolor = getappdata( FIGH, 'DefaultProgressBarColor' );
+    entries = iAddEntry( FIGH, entries, label, 0, defbarcolor, bgcol );
+    idx = numel( entries );
+end
+
+% Check if the user requested a cancel
+if nargout
+    cancel = entries(idx).Cancel;
+end
+
+% Parse the inputs. We shortcut the most common case as an efficiency
+force_update = false;
+if nargin==2 && isnumeric( varargin{1} )
+    entries(idx).LastValue = entries(idx).Value;
+    entries(idx).Value = max( 0, min( 1, varargin{1} ) );
+    entries(idx).Busy = false;
+    needs_update = true;
+else
+    [params,values] = iParseInputs( varargin{:} );
+
+    needs_update = false;
+    for ii=1:numel( params )
+        switch upper( params{ii} )
+            case 'BUSY'
+                entries(idx).Busy = true;
+                needs_update = true;
+
+            case 'VALUE'
+                entries(idx).LastValue = entries(idx).Value;
+                entries(idx).Value = max( 0, min( 1, values{ii} ) );
+                entries(idx).Busy = false;
+                needs_update = true;
+
+            case {'INC','INCREMENT'}
+                entries(idx).LastValue = entries(idx).Value;
+                entries(idx).Value = max( 0, min( 1, entries(idx).Value + values{ii} ) );
+                entries(idx).Busy = false;
+                needs_update = true;
+
+            case {'COLOR','COLOUR'}
+                entries(idx).CData = iMakeColors( values{ii}, 16 );
+                needs_update = true;
+                force_update = true;
+
+            case {'RELABEL', 'UPDATELABEL'}
+                % Make sure we have a string as the value and that it
+                % doesn't already appear
+                if ~ischar( values{ii} ) && ~isStringScalar( values{ii} )
+                    error( 'multiWaitbar:BadString', 'Value for ''Relabel'' must be a string.' );
+                end
+
+                entries(idx).Label = values{ii};
+                if ~isempty(entries(idx).CancelButton)
+                    set( entries(idx).CancelButton, 'Callback', @(src,evt) iCancelEntry(src, entries(idx).Tag) );
+                end
+                needs_update = true;
+                force_update = true;
+
+            case {'CANCANCEL'}
+                if ~ischar( values{ii} ) || ~any( strcmpi( values{ii}, {'on','off'} ) )
+                    error( 'multiWaitbar:BadString', 'Parameter ''CanCancel'' must be a ''on'' or ''off''.' );
+                end
+                entries(idx).CanCancel = strcmpi( values{ii}, 'on' );
+                entries(idx).Cancel = false;
+                needs_redraw = true;
+
+            case {'RESETCANCEL'}
+                entries(idx).Cancel = false;
+                needs_redraw = true;
+
+            case {'CANCELFCN'}
+                if ~isa( values{ii}, 'function_handle' )
+                    error( 'multiWaitbar:BadFunction', 'Parameter ''CancelFcn'' must be a valid function handle.' );
+                end
+                entries(idx).CancelFcn = values{ii};
+                if ~entries(idx).CanCancel
+                    entries(idx).CanCancel = true;
+                end
+                needs_redraw = true;
+
+            case {'CLOSE','DONE'}
+                if ~isempty(idx)
+                    % Remove the selected entry
+                    entries = iDeleteEntry( entries, idx );
+                end
+                if isempty( entries )
+                    iDeleteFigure( FIGH );
+                    % With the window closed, there's nothing else to do
+                    return;
+                else
+                    needs_redraw = true;
+                end
+                % We can't continue after clearing the entry, so jump out
+                break;
+
+            otherwise
+                error( 'multiWaitbar:BadArg', 'Unrecognized command: ''%s''', params{ii} );
+
+        end
+    end
+end
+
+% Now work out what to update/redraw
+if needs_redraw
+    setappdata( FIGH, 'ProgressEntries', entries );
+    iRedraw( FIGH );
+    % NB: Redraw includes updating all bars, so never need to do both
+elseif needs_update
+    [entries(idx),needs_redraw] = iUpdateEntry( entries(idx), force_update );
+    setappdata( FIGH, 'ProgressEntries', entries );
+    % NB: if anything was updated onscreen, "needs_redraw" is now true.
+end
+if entry_added || needs_redraw
+    % If the shape or size has changed, do a full redraw, including events
+    drawnow();
+end
+
+% If we have any "busy" entries, start the timer, otherwise stop it.
+myTimer = getappdata( FIGH, 'BusyTimer' );
+if any([entries.Busy])
+    if strcmpi(myTimer.Running,'off')
+        start(myTimer);
+    end
+else
+    if strcmpi(myTimer.Running,'on')
+        stop(myTimer);
+    end
+end
+
+end % multiWaitbar
+
+
+%-------------------------------------------------------------------------%
+function [params, values] = iParseInputs( varargin )
+% Parse the input arguments, extracting a list of commands and values
+idx = 1;
+params = {};
+values = {};
+if nargin==0
+    return;
+end
+if isnumeric( varargin{1} )
+    params{idx} = 'Value';
+    values{idx} = varargin{1};
+    idx = idx + 1;
+end
+
+while idx <= nargin
+    param = varargin{idx};
+    if ~ischar( param )
+        error( 'multiWaitbar:BadSyntax', 'Additional properties must be supplied as property-value pairs' );
+    end
+    params{end+1,1} = param; %#ok<AGROW>
+    values{end+1,1} = []; %#ok<AGROW>
+    switch upper( param )
+        case {'DONE','CLOSE','RESETCANCEL'}
+            % No value needed, and stop
+            break;
+        case {'BUSY'}
+            % No value needed but parsing should continue
+            idx = idx + 1;
+        case {'RESET','ZERO','SHOW'}
+            % All equivalent to saying ('Value', 0)
+            params{end} = 'Value';
+            values{end} = 0;
+            idx = idx + 1;
+        otherwise
+            if idx==nargin
+                error( 'multiWaitbar:BadSyntax', 'Additional properties must be supplied as property-value pairs' );
+            end
+            values{end,1} = varargin{idx+1};
+            idx = idx + 2;
+    end
+end
+if isempty( params )
+    error( 'multiWaitbar:BadSyntax', 'Must specify a value or a command' );
+end
+end % iParseInputs
+
+%-------------------------------------------------------------------------%
+function fobj = iCreateFig()
+% Create the progress bar group window
+bgcol = get(0,'DefaultUIControlBackgroundColor');
+f = figure( ...
+    'Name', 'Progress', ...
+    'Tag', 'multiWaitbar:Figure', ...
+    'Color', bgcol, ...
+    'MenuBar', 'none', ...
+    'ToolBar', 'none', ...
+    'WindowStyle', 'normal', ... % We don't want to be docked!
+    'HandleVisibility', 'off', ...
+    'IntegerHandle', 'off', ...
+    'Visible', 'off', ...
+    'NumberTitle', 'off' );
+% Resize and centre on the first screen
+screenSize = get(0,'ScreenSize');
+figSz = [360 42];
+figPos = ceil((screenSize(1,3:4)-figSz)/2);
+fobj = handle( f );
+fobj.Position = [figPos, figSz];
+setappdata( fobj, 'ProgressEntries', [] );
+% Make sure we have the image
+defbarcolor = [0.8 0.0 0.1];
+barbgcol = uint8( 255*0.75*bgcol );
+setappdata( fobj, 'DefaultProgressBarBackgroundColor', barbgcol );
+setappdata( fobj, 'DefaultProgressBarColor', defbarcolor );
+setappdata( fobj, 'DefaultProgressBarSize', [350 16] );
+setappdata( fobj, 'MaxEntryRows', 10 );
+% Create the timer to use for "Busy" mode, being sure to delete any
+% existing ones
+delete( timerfind('Tag', 'MultiWaitbarTimer') );
+myTimer = iCreateTimer(f);
+setappdata( fobj, 'BusyTimer', myTimer );
+
+% Setup the resize function after we've finished setting up the figure to
+% avoid excessive redraws
+fobj.ResizeFcn = @iRedraw;
+fobj.CloseRequestFcn = @iCloseFigure;
+end % iCreateFig
+
+%-------------------------------------------------------------------------%
+function t = iCreateTimer(fig)
+t = timer( ...
+    'TimerFcn', @(src,evt) iTimerFcn(fig), ...
+    'Period', 0.02, ...
+    'ExecutionMode', 'FixedRate', ...
+    'Tag', 'MultiWaitbarTimer' );
+end
+
+%-------------------------------------------------------------------------%
+function cdata = iMakeColors( baseColor, height )
+% Creates a shiny bar from a single base color
+lightColor = [1 1 1];
+badColorErrorID = 'multiWaitbar:BadColor';
+badColorErrorMsg = 'Colors must be a three element vector [R G B] or a single character (''r'', ''g'' etc.)';
+
+if ischar(baseColor)
+    switch upper(baseColor)
+        case 'K'
+            baseColor = [0.1 0.1 0.1];
+        case 'R'
+            baseColor = [0.8 0 0];
+        case 'G'
+            baseColor = [0 0.6 0];
+        case 'B'
+            baseColor = [0 0 0.8];
+        case 'C'
+            baseColor = [0.2 0.8 0.9];
+        case 'M'
+            baseColor = [0.6 0 0.6];
+        case 'Y'
+            baseColor = [0.9 0.8 0.2];
+        case 'W'
+            baseColor = [0.9 0.9 0.9];
+        otherwise
+            error( badColorErrorID, badColorErrorMsg );
+    end
+else
+    if numel(baseColor) ~= 3
+        error( badColorErrorID, badColorErrorMsg );
+    end
+    if isa( baseColor, 'uint8' )
+        baseColor = double( baseColor ) / 255;
+    elseif isa( baseColor, 'double' )
+        if any(baseColor>1) || any(baseColor<0)
+            error( 'multiWaitbar:BadColorValue', 'Color values must be in the range 0 to 1 inclusive.' );
+        end
+    else
+        error( badColorErrorID, badColorErrorMsg );
+    end
+end
+
+% By this point we should have a double precision 3-element vector.
+cols = repmat( baseColor, [height, 1] );
+
+breaks = max( 1, round( height * [1 25 50 75 88 100] / 100 ) );
+cols(breaks(1),:) = 0.6*baseColor;
+cols(breaks(2),:) = lightColor - 0.4*(lightColor-baseColor);
+cols(breaks(3),:) = baseColor;
+cols(breaks(4),:) = min( baseColor*1.2, 1.0 );
+cols(breaks(5),:) = min( baseColor*1.4, 0.95 ) + 0.05;
+cols(breaks(6),:) = min( baseColor*1.6, 0.9 ) + 0.1;
+
+y = 1:height;
+cols(:,1) = max( 0, min( 1, interp1( breaks, cols(breaks,1), y, 'pchip' ) ) );
+cols(:,2) = max( 0, min( 1, interp1( breaks, cols(breaks,2), y, 'pchip' ) ) );
+cols(:,3) = max( 0, min( 1, interp1( breaks, cols(breaks,3), y, 'pchip' ) ) );
+cdata = uint8( 255 * cat( 3, cols(:,1), cols(:,2), cols(:,3) ) );
+end % iMakeColors
+
+
+%-------------------------------------------------------------------------%
+function cdata = iMakeBackground( baseColor, height )
+% Creates a shaded background
+if isa( baseColor, 'uint8' )
+    baseColor = double( baseColor ) / 255;
+end
+
+ratio = 1 - exp( -0.5-2*(1:height)/height )';
+cdata = uint8( 255 * cat( 3, baseColor(1)*ratio, baseColor(2)*ratio, baseColor(3)*ratio ) );
+end % iMakeBackground
+
+%-------------------------------------------------------------------------%
+function entries = iAddEntry( parent, entries, label, value, color, bgcolor )
+% Add a new entry to the progress bar
+
+% Create bar coloring
+psize = getappdata( parent, 'DefaultProgressBarSize' );
+cdata = iMakeColors( color, 16 );
+% Create background image
+barcdata = iMakeBackground( bgcolor, psize(2) );
+
+% Work out the size in advance
+mypanel = uipanel( 'Parent', parent, 'Units', 'Pixels', 'BorderType', 'line' );
+labeltext = uicontrol( 'Style', 'Text', ...
+    'String', label, ...
+    'Parent', parent, ...
+    'HorizontalAlignment', 'Left' );
+etatext = uicontrol( 'Style', 'Text', ...
+    'String', '', ...
+    'Parent', parent, ...
+    'HorizontalAlignment', 'Right' );
+% From R2025a onwards we must use UIImage rather than a checkbox.
+if matlab_ver_less_than('25.1')
+    progresswidget = uicontrol( 'Style', 'Checkbox', ...
+        'String', '', ...
+        'Parent', parent, ...
+        'Position', [5 5 psize], ...
+        'CData', barcdata );
+else
+    progresswidget = uiimage( ...
+        'Parent', parent, ...
+        'Position', [5 5 psize], ...
+        'ImageSource', barcdata );
+end
+cancelwidget = uicontrol( 'Style', 'PushButton', ...
+    'String', '', ...
+    'FontWeight', 'Bold', ...
+    'Parent', parent, ...
+    'Position', [5 5 16 16], ...
+    'CData', iMakeCross( 8 ), ...
+    'Callback', @(src,evt) iCancelEntry( src, label ), ...
+    'Visible', 'off' );
+
+newentry = struct( ...
+    'Tag', label, ...
+    'Label', label, ...
+    'Value', value, ...
+    'LastValue', inf, ...
+    'Created', tic(), ...
+    'LabelText', labeltext, ...
+    'ETAText', etatext, ...
+    'ETAString', '', ...
+    'Progress', progresswidget, ...
+    'ProgressSize', psize, ...
+    'Panel', mypanel, ...
+    'BarCData', barcdata, ...
+    'CData', cdata, ...
+    'BackgroundCData', barcdata, ...
+    'CanCancel', false, ...
+    'CancelFcn', [], ...
+    'CancelButton', cancelwidget, ...
+    'Cancel', false, ...
+    'Busy', false );
+if isempty( entries )
+    entries = newentry;
+else
+    entries = [entries;newentry];
+end
+% Store in figure before the redraw
+setappdata( parent, 'ProgressEntries', entries );
+if strcmpi( get( parent, 'Visible' ), 'on' )
+    iRedraw( parent, [] );
+    % Redraw may have modified the entries
+    entries = getappdata( parent, 'ProgressEntries');
+else
+    set( parent, 'Visible', 'on' );
+end
+end % iAddEntry
+
+%-------------------------------------------------------------------------%
+function entries = iDeleteEntry( entries, idx )
+delete( entries(idx).LabelText );
+delete( entries(idx).ETAText );
+delete( entries(idx).CancelButton );
+delete( entries(idx).Progress );
+delete( entries(idx).Panel );
+entries(idx,:) = [];
+end % iDeleteEntry
+
+%-------------------------------------------------------------------------%
+function entries = iCancelEntry( src, name )
+figh = ancestor( src, 'figure' );
+entries = getappdata( figh, 'ProgressEntries' );
+if isempty(entries)
+    % The entries have been lost - nothing can be done.
+    return
+end
+idx = find( strcmp( name, {entries.Tag} ), 1, 'first' );
+
+% Set the cancel flag so that the user is told on next update
+entries(idx).Cancel = true;
+setappdata( figh, 'ProgressEntries', entries );
+
+% If a user function is supplied, call it
+if ~isempty( entries(idx).CancelFcn )
+    feval( entries(idx).CancelFcn, name, 'Cancelled' );
+end
+
+end % iCancelEntry
+
+
+%-------------------------------------------------------------------------%
+function [entry,updated] = iUpdateEntry( entry, force )
+% Update one progress bar
+
+% Deal with busy entries separately
+if entry.Busy
+    entry = iUpdateBusyEntry(entry);
+    updated = true;
+    return;
+end
+
+% Some constants
+marker_weight = 0.8;
+
+% Check if the label needs updating
+updated = force;
+val = entry.Value;
+lastval = entry.LastValue;
+
+% Now update the bar
+psize = entry.ProgressSize;
+filled = max( 1, round( val*psize(1) ) );
+lastfilled = max( 1, round( lastval*psize(1) ) );
+
+% We do some careful checking so that we only redraw what we have to. This
+% makes a small speed difference, but every little helps!
+if force || (filled<lastfilled)
+    % Create the bar background
+    startIdx = 1;
+    bgim = entry.BackgroundCData(:,ones( 1, ceil(psize(1)-filled) ),:);
+    barim = iMakeBarImage(entry.CData, startIdx, filled);
+    progresscdata = [barim,bgim];
+
+    % Add light/shadow around the markers
+    markers = round( (0.1:0.1:val)*psize(1) );
+    markers(markers<=startIdx | markers>(filled-2)) = [];
+    highlight = [marker_weight*entry.CData, 255 - marker_weight*(255-entry.CData)];
+    for ii=1:numel( markers )
+        progresscdata(:,markers(ii)+[-1,0],:) = highlight;
+    end
+    updated = true;
+elseif filled > lastfilled
+    % Just need to update the existing data
+    progresscdata = entry.BarCData;
+    startIdx = max(1,lastfilled-1);
+    % Repmat is the obvious way to fill the bar, but BSXFUN is often
+    % faster. Indexing is obscure but faster still.
+    progresscdata(:,startIdx:filled,:) = iMakeBarImage(entry.CData, startIdx, filled);
+
+    % Add light/shadow around the markers
+    markers = round( (0.1:0.1:val)*psize(1) );
+    markers(markers<startIdx | markers>(filled-2)) = [];
+    highlight = [marker_weight*entry.CData, 255 - marker_weight*(255-entry.CData)];
+    for ii=1:numel( markers )
+        progresscdata(:,markers(ii)+[-1,0],:) = highlight;
+    end
+    updated = true;
+end
+
+if updated
+    % Set the image into the entry and in the widget
+    entry.BarCData = progresscdata;
+    if isa(entry.Progress, "matlab.ui.control.Image")
+        entry.Progress.ImageSource = progresscdata;
+    else
+        set( entry.Progress, 'cdata', progresscdata );
+    end
+end
+
+% As an optimization, don't update any text if the bar didn't move and the
+% percentage hasn't changed
+decval = round( val*100 );
+lastdecval = round( lastval*100 );
+
+if ~updated && (decval == lastdecval)
+    return
+end
+
+% Now work out the remaining time
+minTime = 3; % secs
+if val <= 0
+    % Zero value, so clear the eta
+    entry.Created = tic();
+    elapsedtime = 0;
+    etaString = '';
+else
+    elapsedtime = round(toc( entry.Created )); % in seconds
+
+    % Only show the remaining time if we've had time to estimate
+    if elapsedtime < minTime
+        % Not enough time has passed since starting, so leave blank
+        etaString = '';
+    else
+        % Calculate a rough ETA
+        eta = elapsedtime * (1-val) / val;
+        etaString = iGetTimeString( eta );
+    end
+end
+
+if ~isequal( etaString, entry.ETAString )
+    set( entry.ETAText, 'String', etaString );
+    entry.ETAString = etaString;
+    updated = true;
+end
+
+% Update the label too
+if force || elapsedtime > minTime
+    if force || (decval ~= lastdecval)
+        labelstr = [entry.Label, sprintf( ' (%d%%)', decval )];
+        set( entry.LabelText, 'String', labelstr );
+        updated = true;
+    end
+end
+
+end % iUpdateEntry
+
+function eta = iGetTimeString( remainingtime )
+if remainingtime > 172800 % 2 days
+    eta = sprintf( '%d days', round(remainingtime/86400) );
+else
+    if remainingtime > 7200 % 2 hours
+        eta = sprintf( '%d hours', round(remainingtime/3600) );
+    else
+        if remainingtime > 120 % 2 mins
+            eta = sprintf( '%d mins', round(remainingtime/60) );
+        else
+            % Seconds
+            remainingtime = round( remainingtime );
+            if remainingtime > 1
+                eta = sprintf( '%d secs', remainingtime );
+            elseif remainingtime == 1
+                eta = '1 sec';
+            else
+                eta = ''; % Nearly done (<1sec)
+            end
+        end
+    end
+end
+end % iGetTimeString
+
+
+%-------------------------------------------------------------------------%
+function entry = iUpdateBusyEntry( entry )
+% Update a "busy" progress bar
+% Make sure the widget is still OK
+if ~ishandle(entry.Progress)
+    return
+end
+% Work out the new position. Since the bar is 0.1 long and needs to bounce,
+% the position varies from 0 up to 0.9 then back down again. We achieve
+% this with judicious use of "mod" with 1.8.
+entry.Value = mod(entry.Value+0.01,1.8);
+val = entry.Value;
+if val>0.9
+    % Moving backwards
+    val = 1.8-val;
+end
+psize = entry.ProgressSize;
+startIdx = max( 1, round( val*psize(1) ) );
+endIdx = max( 1, round( (val+0.1)*psize(1) ) );
+barLength = endIdx - startIdx + 1;
+
+% Create the image
+bgim = entry.BackgroundCData(:,ones( 1, psize(1) ),:);
+barim = iMakeBarImage(entry.CData, 1, barLength);
+bgim(:,startIdx:endIdx,:) = barim;
+
+% Put it into the widget
+entry.BarCData = bgim;
+if isa(entry.Progress,"matlab.ui.control.Image")
+    entry.Progress.ImageSource = bgim;
+else
+    set(entry.Progress, 'cdata', bgim);
+end
+end % iUpdateBusyEntry
+
+
+%-------------------------------------------------------------------------%
+function barim = iMakeBarImage(strip, startIdx, endIdx)
+shadow1_weight = 0.4;
+shadow2_weight = 0.7;
+barLength = endIdx - startIdx + 1;
+% Repmat is the obvious way to fill the bar, but BSXFUN is often
+% faster. Indexing is obscure but faster still.
+barim = strip(:,ones(1, barLength),:);
+% Add highlight to the start of the bar
+if startIdx <= 2 && barLength>=2
+    barim(:,1,:) = 255 - shadow1_weight*(255-strip);
+    barim(:,2,:) = 255 - shadow2_weight*(255-strip);
+end
+% Add shadow to the end of the bar
+if endIdx>=4 && barLength>=2
+    barim(:,end,:) = shadow1_weight*strip;
+    barim(:,end-1,:) = shadow2_weight*strip;
+end
+end % iMakeBarImage
+
+%-------------------------------------------------------------------------%
+function iCloseFigure( fig, evt ) %#ok<INUSD>
+% Closing the figure just makes it invisible
+set( fig, 'Visible', 'off' );
+end % iCloseFigure
+
+%-------------------------------------------------------------------------%
+function iDeleteFigure( fig )
+% Actually destroy the figure
+busyTimer = getappdata( fig, 'BusyTimer' );
+if isvalid( busyTimer )
+    stop( busyTimer );
+end
+delete( busyTimer );
+delete( fig );
+end % iDeleteFigure
+
+%-------------------------------------------------------------------------%
+function iRedraw( fig, evt ) %#ok<INUSD>
+entries = getappdata( fig, 'ProgressEntries' );
+fobj = handle( fig );
+p = fobj.Position;
+% p = get( fig, 'Position' );
+border = 5;
+textheight = 16;
+barheight = 16;
+panelheight = 10;
+maxRows = getappdata( fig, 'MaxEntryRows' );
+N = max( 1, numel( entries ) );
+Nrows = min(maxRows, N);
+Ncols = ceil(N ./ Nrows);
+
+% Check the height is correct
+heightperentry = textheight+barheight+panelheight;
+requiredheight = 2*border + Nrows*heightperentry - panelheight;
+if ~isequal( p(4), requiredheight )
+    p(2) = p(2) + p(4) - requiredheight;
+    p(4) = requiredheight;
+    % In theory setting the position should re-trigger this callback, but
+    % in practice it doesn't, probably because we aren't calling "drawnow".
+    set( fig, 'Position', p )
+end
+width = floor((p(3) - Ncols*2*border) ./ Ncols);
+setappdata( fig, 'DefaultProgressBarSize', [width barheight] );
+
+for ii=1:numel( entries )
+    col = ceil(ii./Nrows)-1;
+    row = ii - col.*Nrows - 1;
+    xpos = border + (width+2*border).*col;
+    ypos = p(4) - border - heightperentry.*row;
+
+    set( entries(ii).Panel, 'Position', [xpos-border ypos+panelheight/2-heightperentry width+2*border+1 heightperentry] );
+    set( entries(ii).LabelText, 'Position', [xpos ypos-textheight width*0.75 textheight] );
+    set( entries(ii).ETAText, 'Position', [xpos+width*0.75 ypos-textheight width*0.25 textheight] );
+    ypos = ypos - textheight;
+
+    if entries(ii).CanCancel
+        set( entries(ii).Progress, 'Position', [xpos ypos-barheight max(1,width-barheight+1) barheight] );
+        entries(ii).ProgressSize = max(1,[width-barheight barheight]);
+        set( entries(ii).CancelButton, 'Visible', 'on', 'Position', [p(3)-border-barheight ypos-barheight barheight barheight] );
+    else
+        set( entries(ii).Progress, 'Position', [xpos ypos-barheight width+1 barheight] );
+        entries(ii).ProgressSize = [width barheight];
+        set( entries(ii).CancelButton, 'Visible', 'off' );
+    end
+    entries(ii) = iUpdateEntry( entries(ii), true );
+end
+setappdata( fig, 'ProgressEntries', entries );
+end % iRedraw
+
+function cdata = iMakeCross( sz )
+% Create a cross-shape icon of size sz*sz*3
+
+cdata = diag(ones(sz,1),0) + diag(ones(sz-1,1),1) + diag(ones(sz-1,1),-1);
+cdata = cdata + flip(cdata,2);
+
+% Convert zeros to nans (transparent) and non-zeros to zero (black)
+cdata(cdata == 0) = nan;
+cdata(~isnan(cdata)) = 0;
+
+% Convert to RGB
+cdata = cat( 3, cdata, cdata, cdata );
+end % iMakeCross
+
+
+function iTimerFcn(fig)
+% Timer callback for updating stuff every so often
+entries = getappdata( fig, 'ProgressEntries' );
+for ii=1:numel(entries)
+    if entries(ii).Busy
+        entries(ii) = iUpdateBusyEntry(entries(ii));
+    end
+end
+setappdata( fig, 'ProgressEntries', entries );
+end % iTimerFcn
\ No newline at end of file