matlab/+occbin/map_regime.m

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

matlab/+occbin/match_function.m

 function [resids, grad, state_out, E, M_, out] = match_function(err_0, obs_list,current_obs, opts_simul,...
-                                                                            M_, oo_, options_)
+                                                                            M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_)
 % function [resids, grad, stateout, E, M_, out] = match_function(err_0, obs_list,current_obs, opts_simul,...
-%                                                                             M_, oo_, options_)
+%                                                                             M_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, options_)
 % Outputs:
 %  - resids         [double]        [n_exo by 1] vector of residuals
 %  - 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_).
-% - 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_).
+% - 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_).
 
 % Original authors: Pablo Cuba-Borda, Luca Guerrieri, Matteo Iacoviello, and Molin Zhong
 % Original file downloaded from:
@@ -36,7 +39,7 @@ opts_simul.SHOCKS = err_0';
 options_.occbin.simul=opts_simul;
 options_.occbin.simul.full_output=1;
 options_.noprint = 1;
-[~, out, ss] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+[~, out, ss] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
 
 nobs = size(obs_list,1);
 resids = zeros(nobs,1);

matlab/+occbin/mkdatap_anticipated_2constraints_dyn.m

@@ -98,11 +98,7 @@ if T_max > 0
         % check if last binding regime was already stored
         tmp = 0*binding_indicator;
         tmp(1:end-T_max+1,:) = binding_indicator(T_max:end,:);
-        if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
-            itmp = find(~any(bsxfun(@minus, dictionary.binding_indicator(1:length(tmp)*2,:), tmp(:))));
-        else
         itmp = find(~any(dictionary.binding_indicator(1:length(tmp)*2,:)-tmp(:)));
-        end
         if ~isempty(itmp)
             ireg(T_max) = itmp;
         else
@@ -134,11 +130,7 @@ if T_max > 0
     for i = T_max-1:-1:1        
         tmp = 0*binding_indicator;
         tmp(1:end-i+1,:) = binding_indicator(i:end,:);
-        if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
-            itmp = find(~any(bsxfun(@minus, dictionary.binding_indicator(1:length(tmp)*2,:), tmp(:))));
-        else
         itmp = find(~any(dictionary.binding_indicator(1:length(tmp)*2,:)-tmp(:)));
-        end
         if ~isempty(itmp)
             ireg(i) = itmp;
         else

matlab/+occbin/mkdatap_anticipated_dyn.m

@@ -89,11 +89,7 @@ if T_max > 0
         
         tmp = 0*binding_indicator;
         tmp(1:end-i+1) = binding_indicator(i:end);
-        if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
-            itmp = find(~any(bsxfun(@minus, dictionary.binding_indicator, tmp)));
-        else
         itmp = find(~any(dictionary.binding_indicator-tmp));
-        end
         if ~isempty(itmp)
             ireg(i) = itmp;
         else

matlab/+occbin/plot_irfs.m

+function plot_irfs(M_,irfs,options_,var_list)
+% plot_irfs(M_,irfs,options_,var_list)
+%
+% INPUTS
+% - M_                      [structure]     MATLAB's structure describing the model
+% - irfs                    [structure]     IRF results
+% - options_                [structure]     MATLAB's structure describing the current options
+% - var_list                [character array]  list of endogenous variables specified
+%
+% OUTPUTS
+%   none
+%
+% SPECIAL REQUIREMENTS
+%   none.
+
+% Copyright © 2022-2023 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+if (isfield(options_,'irf_shocks')==0)
+    shock_names  = M_.exo_names;
+else
+    shock_names  = options_.irf_shocks;
+end
+
+simul_name = options_.occbin.plot_irf.simulname;
+if isempty(simul_name)
+    save_name = simul_name;
+else
+    save_name = [ simul_name '_' ];
+end
+
+if isempty(var_list)
+    var_list = M_.endo_names(1:M_.orig_endo_nbr);
+end
+
+[i_var, ~, index_uniques] = varlist_indices(var_list, M_.endo_names);
+vars_irf=var_list(index_uniques);
+
+endo_scaling_factor = options_.occbin.plot_irf.endo_scaling_factor;
+length_irf = options_.irf;
+
+steps_irf = 1;
+
+DirectoryName = CheckPath('graphs',M_.dname);
+latexFolder = CheckPath('latex',M_.dname);
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fidTeX = fopen([latexFolder '/' M_.fname '_occbin_irfs.tex'],'w');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by occbin.plot_irfs.m (Dynare).\n');
+    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+    fprintf(fidTeX,' \n');
+end
+
+iexo=[];
+for var_iter=1:size(shock_names,1)
+    itemp = strmatch(shock_names{var_iter},M_.exo_names,'exact');
+    if isempty(itemp)
+        error(['Shock ',shock_names{var_iter},' is not defined!'])
+    else
+        iexo=[iexo, itemp];
+    end
+end
+
+ncols     = options_.occbin.plot_irf.ncols;
+nrows     = options_.occbin.plot_irf.nrows;
+npan      = ncols*nrows;
+
+shocksigns = options_.occbin.plot_irf.shocksigns;
+threshold = options_.impulse_responses.plot_threshold;
+
+for shock_sign_iter = 1:numel(shocksigns)
+    shocksign = shocksigns{shock_sign_iter};
+    if strcmp(shocksign,'pos')
+        plot_title_sign='positive';
+    elseif strcmp(shocksign,'neg')
+        plot_title_sign='negative';
+    else
+       error('Unknown shock sign %s',shocksign);
+    end
+        
+    for shock_iter=1:size(shock_names,1)
+        j1   = 0;
+        isub = 0;
+        ifig = 0;
+        % Variables
+        for var_iter = 1:length(vars_irf)
+            j1=j1+1;
+            if mod(j1,npan)==1
+                % vector corresponds to [left bottom width height]. 680 and 678 for the left and bottom elements correspond to the default values used by MATLAB while creating a figure and width, .
+                screensize = get( groot, 'Screensize' );
+                hfig = dyn_figure(options_.nodisplay,'name',['OccBin IRFs to ' plot_title_sign ' ' shock_names{shock_iter} ' shock ' simul_name],'OuterPosition' ,[50 50 min(1000,screensize(3)-50) min(750,screensize(4)-50)]);
+                ifig=ifig+1;
+                isub=0;
+            end
+            isub=isub+1;
+            if isempty(endo_scaling_factor)
+                exofactor = 1;
+            else
+                exofactor = endo_scaling_factor{var_iter};
+            end
+
+            subplot(nrows,ncols,isub)
+            irf_field   = strcat(vars_irf{var_iter},'_',shock_names{shock_iter},'_',shocksign);
+
+            %%linear IRFs
+            if ~isfield(irfs.linear,irf_field)
+                warning('occbin.plot_irfs: no linear IRF for %s to %s detected',vars_irf{var_iter,1},shock_names{shock_iter})
+            else
+                irfvalues   = irfs.linear.(irf_field);
+                irfvalues(abs(irfvalues) <threshold) = 0;
+                if options_.occbin.plot_irf.add_steadystate
+                    irfvalues = irfvalues + get_mean(vars_irf{var_iter,1});
+                end
+                max_irf_length_1=min(length_irf,length(irfvalues));
+                plot(irfvalues(1:steps_irf:max_irf_length_1)*exofactor,'linewidth',2);
+            end
+            hold on
+
+            if ~isfield(irfs.piecewise,irf_field)
+                warning('occbin.plot_irfs: no piecewise linear IRF for %s to %s detected',vars_irf{var_iter,1},shock_names{shock_iter})
+            else
+                irfvalues   = irfs.piecewise.(irf_field);
+                irfvalues(abs(irfvalues) <threshold) = 0;
+                if options_.occbin.plot_irf.add_steadystate
+                    irfvalues = irfvalues + get_mean(vars_irf{var_iter,1});
+                end
+                max_irf_length_2=min(length_irf,length(irfvalues));
+                plot(irfvalues(1:steps_irf:max_irf_length_2)*exofactor,'r--','linewidth',2);
+            end
+
+            plot(irfvalues(1:steps_irf:max(max_irf_length_1,max_irf_length_2))*0,'k-','linewidth',1.5);
+
+            if options_.occbin.plot_irf.grid
+                grid on
+            end     
+            xlim([1 max(max_irf_length_1,max_irf_length_2)]);
+            if options_.TeX
+                title(['$' M_.endo_names_tex{i_var(var_iter)}, '$'],'Interpreter','latex')
+            else
+                title(M_.endo_names{i_var(var_iter)},'Interpreter','none')
+            end
+            % Annotation Box + save figure
+            % ----------------------
+            if mod(j1,npan)==0 || (mod(j1,npan)~=0 && var_iter==length(vars_irf))
+                annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
+                annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
+                dyn_saveas(hfig,[DirectoryName,filesep,M_.fname,'_irf_occbin_',save_name,shock_names{shock_iter},'_',shocksign,'_',int2str(ifig)],options_.nodisplay,options_.graph_format);
+                if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+                    fprintf(fidTeX,'\\begin{figure}[H]\n');
+                    fprintf(fidTeX,'\\centering \n');
+                    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_irf_occbin_%s}\n',options_.figures.textwidth*min((j1-1)/ncols,1),...
+                        [DirectoryName '/' ,M_.fname],[save_name,shock_names{shock_iter},'_',shocksign,'_',int2str(ifig)]);
+                    fprintf(fidTeX,'\\caption{OccBin IRFs to  %s shock to %s}\n',plot_title_sign,shock_names{shock_iter});
+                    fprintf(fidTeX,'\\label{Fig:occbin_irfs:%s}\n',[save_name,shock_names{shock_iter},'_',shocksign,'_',int2str(ifig)]);
+                    fprintf(fidTeX,'\\end{figure}\n');
+                    fprintf(fidTeX,'\n');
+                end
+            end
+        end
+    end
+end
+
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fprintf(fidTeX,'%% End Of TeX file.');
+    fclose(fidTeX);
+end
\ No newline at end of file

matlab/+occbin/plot_regimes.m

+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
+
+% Copyright © 2021-2023 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+nperiods = size(regimes,2);
+nconstr = length(fieldnames(regimes(1)))/2;
+if nconstr ==1
+    regime(1) = {'regime'};
+    regimestart(1) = {'regimestart'};
+else
+    regime(1) = {'regime1'};
+    regimestart(1) = {'regimestart1'};
+    regime(2) = {'regime2'};
+    regimestart(2) = {'regimestart2'};
+end    
+
+GraphDirectoryName = CheckPath('graphs',M_.dname);
+fhandle = dyn_figure(options_.nodisplay,'Name',[M_.fname ': OccBin regimes']);
+
+latexFolder = CheckPath('latex',M_.dname);
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fidTeX = fopen([latexFolder '/' M_.fname '_occbin_regimes.tex'],'w');
+    fprintf(fidTeX,'%% TeX eps-loader file generated by occbin.plot_regimes.m (Dynare).\n');
+    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+    fprintf(fidTeX,' \n');
+end
+
+for k=1:nconstr
+    subplot(nconstr,1,k)
+    for t=1:nperiods
+        nregimes_in_t = length(regimes(t).(regime{k}));
+        start_periods = regimes(t).(regimestart{k});
+        start_periods = [start_periods max(start_periods)];
+        for r=1:nregimes_in_t
+            isconstrained = regimes(t).(regime{k})(r);
+            if isconstrained
+                plot(t,start_periods(r),'*r')
+                hold all,
+                plot([t t],start_periods(r:r+1),'-r')
+            else
+                plot(t,start_periods(r),'ob')
+                hold all,
+                plot([t t],start_periods(r:r+1),'-b')
+            end
+        end
+    end
+    xlim([1 nperiods])
+    title(['Regime ' int2str(k)])
+    xlabel('Historic period')
+    ylabel('Regime: expected start')
+end
+annotation('textbox',[.25,0,.15,.05],'String','Slack','Color','blue');
+annotation('textbox',[.65,0,.2,.05],'String','Binding','Color','red');
+
+dyn_saveas(fhandle,[GraphDirectoryName, filesep, M_.fname '_occbin_regimes'],options_.nodisplay,options_.graph_format);
+
+if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
+    fprintf(fidTeX,'\\begin{figure}[H]\n');
+    fprintf(fidTeX,'\\centering \n');
+    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_occbin_regimes}\n',options_.figures.textwidth,[GraphDirectoryName '/' M_.fname]);
+    fprintf(fidTeX,'\\caption{OccBin: regime evolution over time.}\n');
+    fprintf(fidTeX,'\\label{Fig:occbin_regimes}\n');
+    fprintf(fidTeX,'\\end{figure}\n');
+    fprintf(fidTeX,'\n');
+    fclose(fidTeX);
+end

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
 %
@@ -42,30 +42,36 @@ if ismember(flag,{'all'})
 end
 
 if ismember(flag,{'filter','all'})
+    options_occbin_.filter.state_covariance = false;
+    options_occbin_.filter.guess_regime = false;
+    options_occbin_.filter.periodic_solution = true;
     options_occbin_.filter.use_relaxation = false;
+    options_occbin_.filter.use_relaxation_tol_period = 1;
 end
 
 if ismember(flag,{'forecast','all'})
+    options_occbin_.forecast.check_ahead_periods=30;
     options_occbin_.forecast.debug_flag=false;
     options_occbin_.forecast.frcst_regimes=[];
     options_occbin_.forecast.maxit=30;
-    options_occbin_.forecast.periods=30;
     options_occbin_.forecast.qmc=0;
     options_occbin_.forecast.replic=0;
     options_occbin_.forecast.SHOCKS0=[];
-    options_occbin_.forecast.treepath=1; % number of branches
 end
 
 if ismember(flag,{'irf','all'})    
+    options_occbin_.irf.check_ahead_periods=30;
+    options_occbin_.irf.exo_names=M_.exo_names;
     options_occbin_.irf.init_regime=[];
     options_occbin_.irf.maxit=30;
-    options_occbin_.irf.threshold = 10^-6;
 %     options_occbin_.irf.periods=options_.irf;
     options_occbin_.irf.shocksize=[];
-    options_occbin_.irf.shocksigns = {'1','_1'}; 
+    options_occbin_.irf.shocksigns = {'pos','neg'}; 
+    options_occbin_.irf.t0=0;
 end
 
 if ismember(flag,{'likelihood','all'})
+    options_occbin_.likelihood.brute_force_regime_guess = true;
     options_occbin_.likelihood.curb_retrench = false;
     options_occbin_.likelihood.first_period_occbin_update = 1;
     options_occbin_.likelihood.full_output = false;
@@ -74,10 +80,12 @@ if ismember(flag,{'likelihood','all'})
     options_occbin_.likelihood.init_binding_indicator = false(0);
     options_occbin_.likelihood.inversion_filter = false;
     options_occbin_.likelihood.IVF_shock_observable_mapping = [];
+    options_occbin_.likelihood.loss_function_regime_guess = false;
     options_occbin_.likelihood.maxit = 30; % this is for occbin solver algo
     options_occbin_.likelihood.max_number_of_iterations = 10; % this is for occbin_kalman_update loop
     options_occbin_.likelihood.max_check_ahead_periods=inf;
-    options_occbin_.likelihood.periods = 100;
+    options_occbin_.likelihood.number_of_initial_periods_with_extra_regime_guess=0;
+    options_occbin_.likelihood.periods = 3;
     options_occbin_.likelihood.check_ahead_periods=200;
     options_occbin_.likelihood.periodic_solution=false;
     options_occbin_.likelihood.piecewise_only = true;
@@ -87,6 +95,16 @@ if ismember(flag,{'likelihood','all'})
     options_occbin_.likelihood.waitbar=false;
 end
 
+if ismember(flag,{'plot_irf','all'})
+    options_occbin_.plot_irf.add_steadystate = 0;
+    options_occbin_.plot_irf.endo_scaling_factor = [];
+    options_occbin_.plot_irf.grid            = true;
+    options_occbin_.plot_irf.ncols            = 3;
+    options_occbin_.plot_irf.nrows            = 3;
+    options_occbin_.plot_irf.shocksigns = {'pos','neg'}; 
+    options_occbin_.plot_irf.simulname='';
+end
+
 if ismember(flag,{'plot_shock_decomp','all'})
     options_occbin_.plot_shock_decomp.add_steadystate = false;
     options_occbin_.plot_shock_decomp.add_zero_line = false;
@@ -175,6 +193,8 @@ if ismember(flag,{'simul','all'})
     options_occbin_.simul.periods=100;
     options_occbin_.simul.check_ahead_periods=200;
     options_occbin_.simul.periodic_solution=false;
+    options_occbin_.simul.periodic_solution_threshold=1;
+    options_occbin_.simul.periodic_solution_strict=true;
     options_occbin_.simul.piecewise_only = false;
     options_occbin_.simul.reset_check_ahead_periods_in_new_period = false;
     options_occbin_.simul.reset_regime_in_new_period = false;
@@ -197,7 +217,7 @@ if ismember(flag,{'smoother','all'})
     options_occbin_.smoother.maxit = 30; % this is for occbin solver algo
     options_occbin_.smoother.max_check_ahead_periods=inf;
     options_occbin_.smoother.max_number_of_iterations = 10; % this is for smoother loop
-    options_occbin_.smoother.periods = 100;
+    options_occbin_.smoother.periods = 3;
     options_occbin_.smoother.check_ahead_periods=200;
     options_occbin_.smoother.periodic_solution=false;
     options_occbin_.smoother.piecewise_only = true;

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

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

matlab/+occbin/solver.m

 function [dr, out, ss] = solver(M_, options_, dr ,steady_state, exo_steady_state, exo_det_steady_state)
-% function [oo_, out, ss] = solver(M_,oo_,options_, dr ,steady_state, exo_steady_state, exo_det_steady_state
+% [dr, out, ss] = solver(M_,oo_,options_, dr ,steady_state, exo_steady_state, exo_det_steady_state
 % Solves the model with an OBC and produces simulations/IRFs
 %
 % INPUT: 
-% - M_                  [structure]     Matlab's structure describing the model
-% - oo_                 [structure]     Matlab's structure containing the results
-% - options_            [structure]     Matlab's structure containing the options
+% - M_                      [structure]     MATLAB's structure describing the model
+% - options_                [structure]     MATLAB's structure containing the options
+% - dr                      [structure]     model information structure
+% - endo_steady_state       [vector]        steady state value for endogenous variables
+% - exo_steady_state        [vector]        steady state value for exogenous variables
+% - exo_det_steady_state    [vector]        steady state value for exogenous deterministic variables                                    
 %
 % OUTPUT: 
-% - oo_                 [structure]     Matlab's structure containing the results
+% - dr                      [structure]     decision rules
 % - out                     [structure]     simulation result containing fields:
 %                                               - linear: paths for endogenous variables ignoring OBC (linear solution)
 %                                               - piecewise: paths for endogenous variables satisfying the OBC (occbin/piecewise solution)
@@ -19,7 +22,7 @@ function [dr, out, ss] = solver(M_, options_, dr ,steady_state, exo_steady_state
 %                                               - R: [n_vars by n_exo by n_shock_period] array of shock response matrices
 %                                               - C: [n_vars by n_shock_period] array of constants
 
-% Copyright © 2021 Dynare Team
+% Copyright © 2021-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -50,6 +53,7 @@ else
     end
 end
 
+ss=struct();
 if solve_dr
     if isempty(options_.qz_criterium)
         options_.qz_criterium = 1+1e-6;
@@ -90,15 +94,7 @@ end
 
 % add back steady state
 if ~options_.occbin.simul.piecewise_only
-    if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
-        out.linear = bsxfun(@plus, out.linear, out.ys');
-    else
     out.linear = out.linear + out.ys';
 end
-end
-if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
-    out.piecewise = bsxfun(@plus, out.piecewise, out.ys');
-else
 out.piecewise = out.piecewise + out.ys';
-end
 out.exo_pos = options_.occbin.simul.exo_pos;

matlab/+occbin/squeeze_shock_decomposition.m

+function [oo_,options_] = squeeze_shock_decomposition(M_,oo_,options_, sd_vlist)
+
+if isstruct(options_.plot_shock_decomp.q2a)
+    
+    avname=char({options_.plot_shock_decomp.q2a.qname});
+    sda = options_.plot_shock_decomp.q2a(ismember(avname,sd_vlist,'rows'));
+    for k=1:length(sda)
+        if isstruct(sda(k).aux)
+             sd_vlist = [sd_vlist; cellstr(sda(k).aux.y)];
+        end
+    end
+end
+i_var = varlist_indices(sd_vlist,M_.endo_names);
+sd_vlist = M_.endo_names(i_var);
+% first we squeeze usual fields 
+oo_ = squeeze_shock_decomposition(M_,oo_,options_,sd_vlist);
+i_var = oo_.shock_decomposition_info.i_var;
+sd_vlist = M_.endo_names(i_var);
+
+% now we check for occbin SDs
+options_.occbin.shock_decomp.i_var = i_var;
+if isfield (oo_.occbin.smoother,'decomp')
+    oo_.occbin.smoother.decomp = oo_.occbin.smoother.decomp(i_var,:,:);
+    oo_.occbin.smoother.wdecomp = oo_.occbin.smoother.wdecomp(i_var,:,:);
+end
+
+if isfield(oo_.occbin,'shock_decomp')
+    fnames = fieldnames(oo_.occbin.shock_decomp);
+    for k=1:length(fnames)
+        nendo = numel(oo_.occbin.shock_decomp.(fnames{k}).vname);
+        tmp_i_var = varlist_indices(sd_vlist,char(oo_.occbin.shock_decomp.(fnames{k}).vname));
+        oo_.occbin.shock_decomp.(fnames{k}).vname = cellstr(sd_vlist);
+        tmpnames = fieldnames(oo_.occbin.shock_decomp.(fnames{k}));
+        for t=1:length(tmpnames)
+            if size(oo_.occbin.shock_decomp.(fnames{k}).(tmpnames{t}),3)==nendo
+                oo_.occbin.shock_decomp.(fnames{k}).(tmpnames{t})= oo_.occbin.shock_decomp.(fnames{k}).(tmpnames{t})(:,:,tmp_i_var);
+            end
+        end
+    end
+end
+
+end
+
+

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,10 +4,10 @@ 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 Dynare Team
+% Copyright © 2021-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -62,14 +62,7 @@ else
     end
 end
 
-if ~ispc && ~isoctave && matlab_ver_less_than('9.0')
-    % On GNU/Linux and macOS, with MATLAB < R2016a, “writeable” can’t write Excel files
-    % (and “xlswrite” can’t either)
-    warning('This version of MATLAB is too old and cannot create Excel files. The Occbin regimes will rather be written to a CSV file.')
-    filename=[OutputDirectoryName filesep xls_filename '.csv'];
-else
 filename=[OutputDirectoryName filesep xls_filename '.xls'];
-end
 
 if isfile(filename)
     delete(filename)
@@ -81,10 +74,6 @@ if isoctave
         error('The io package is required to write XLS files from Octave')
     end
     xlswrite(filename, vertcat(Header, xlsmat));
-elseif  ~ispc && matlab_ver_less_than('9.0')
-    % Use a CSV file. See the comment above about filename.
-    % We don’t use xlswrite because its CSV fallback does not support cell-arrays.
-    writetable(array2table(xlsmat,'VariableNames',Header), filename);
 else
     writetable(array2table(xlsmat,'VariableNames',Header), filename, 'Sheet', 'Regimes');
 end

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;
@@ -64,9 +68,9 @@ if ~options_.analytic_derivation
     loss = full(weights(:)'*vx(:));
 else
     totparam_nbr=length(i_params);
-    oo_.dr.derivs = get_perturbation_params_derivs(M_, options_, [], oo_, i_params, [], [], 0); %analytic derivatives of perturbation matrices
+    oo_.dr.derivs = identification.get_perturbation_params_derivs(M_, options_, [], oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, i_params, [], [], 0); %analytic derivatives of perturbation matrices
 
-    pruned_state_space = pruned_state_space_system(M_, options_, oo_.dr, i_var, 0, 0, 1);
+    pruned_state_space = pruned_SS.pruned_state_space_system(M_, options_, oo_.dr, i_var, 0, 0, 1);
     vx = pruned_state_space.Var_y + pruned_state_space.E_y*pruned_state_space.E_y';
     dE_yy = pruned_state_space.dVar_y;
     for jp=1:length(i_params)
@@ -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

@@ -105,7 +105,7 @@ if isfield(options_.osr,'maxit') || isfield(options_.osr,'tolf')
     end
 end
 
-oo_.dr = set_state_space(oo_.dr,M_,options_);
+oo_.dr = set_state_space(oo_.dr,M_);
 par_0 = M_.params(i_params);
 inv_order_var = oo_.dr.inv_order_var;
 

matlab/+pac/+estimate/init.m

@@ -18,7 +18,7 @@ function [pacmodl, lhs, rhs, pnames, enames, xnames, rname, pid, eid, xid, pname
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 % Get the original equation to be estimated
-[LHS, RHS] = get_lhs_and_rhs(eqname, M_, true);
+[~, RHS] = get_lhs_and_rhs(eqname, M_, true);
 
 % Check that the equation has a PAC expectation term.
 if ~contains(RHS, 'pac_expectation', 'IgnoreCase', true)

matlab/+pac/+estimate/iterative_ols.m

@@ -41,7 +41,7 @@ function iterative_ols(eqname, params, data, range)
 
 global M_ oo_ options_
 
-[pacmodl, ~, rhs, ~, ~, ~, rname, ~, ~, ~, ~, ipnames_, params, data, ~] = ...
+[pacmodl, ~, rhs, ~, ~, ~, rname, ~, ~, ~, ~, ipnames_, params, data] = ...
     pac.estimate.init(M_, oo_, eqname, params, data, range);
 
 % Set initial condition.