matlab/+occbin/forecast.m

-function [oo_, error_flag] = forecast(options_,M_,oo_,forecast) %,hist_period)
-%function oo_ = forecast(options_,M_,oo_,forecast)
-% forecast
+function [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
+% [forecast, error_flag] = forecast(options_,M_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state,forecast_horizon)
+% Occbin forecasts
+%
+% INPUTS
+% - options_                [structure]     Matlab's structure describing the current options
+% - M_                      [structure]     Matlab's structure describing the model
+% - dr_in                   [structure]     model information structure
+% - endo_steady_state       [double]        steady state value for endogenous variables
+% - exo_steady_state        [double]        steady state value for exogenous variables
+% - exo_det_steady_state    [double]        steady state value for exogenous deterministic variables                                    
+% - forecast_horizon        [integer]       forecast horizon
+%
+% OUTPUTS
+% - forecast                [structure]     forecast results
+% - error_flag              [integer]       error code
+%
+% 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/>.
 
 opts = options_.occbin.forecast;
 
 options_.occbin.simul.maxit = opts.maxit;
 options_.occbin.simul.check_ahead_periods = opts.check_ahead_periods;
-options_.occbin.simul.periods = forecast;
-SHOCKS0 = opts.SHOCKS0;
-if ~isempty(SHOCKS0)
-    if iscell(SHOCKS0)
-        for j=1:length(SHOCKS0)
-            sname = SHOCKS0{j}{1};
-            inds(j)=strmatch(sname,M_.exo_names);
-            SHOCKS1(j,:)=SHOCKS0{j}{2};
+options_.occbin.simul.periods = forecast_horizon;
+shocks_input = opts.SHOCKS0;
+
+if ~isempty(shocks_input)
+    n_shocks=size(shocks_input,1);
+    if iscell(shocks_input)
+        inds=NaN(n_shocks,1);
+        periods=length(shocks_input{1}{2});
+        shock_mat=NaN(n_shocks,periods);
+        for j=1:n_shocks
+            exo_pos=strmatch(shocks_input{j}{1},M_.exo_names,'exact');
+            if isempty(exo_pos)
+        	    error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
+            else
+                inds(j)=exo_pos;
+            end
+            if length(shocks_input{j}{2})~=periods
+        	    error('occbin.forecast: unknown exogenous shock %s',shocks_input{j}{1})
+            else
+                shock_mat(j,:)=shocks_input{j}{2};
             end
-    elseif isreal(SHOCKS0)
-        SHOCKS1=SHOCKS0;
-        inds = 1:M_.exo_nbr;
+        end
+    elseif isreal(shocks_input)
+        shock_mat=shocks_input;
+        inds = (1:M_.exo_nbr)';
     end
 end
 
+options_.occbin.simul.endo_init = M_.endo_histval(:,1)-endo_steady_state; %initial condition
+options_.occbin.simul.init_regime = opts.frcst_regimes;
+options_.occbin.simul.init_binding_indicator = [];
+
+shocks_base = zeros(forecast_horizon,M_.exo_nbr);
+if ~isempty(shocks_input)
+    for j=1:n_shocks
+        shocks_base(:,inds(j))=shock_mat(j,:);
+    end
+end
 
 if opts.replic
     h = dyn_waitbar(0,'Please wait occbin forecast replic ...');
     ishock = find(sqrt(diag((M_.Sigma_e))));
+    options_.occbin.simul.exo_pos=ishock;
     effective_exo_nbr=  length(ishock);
-    effective_exo_names = M_.exo_names(ishock);
-    effective_Sigma_e = M_.Sigma_e(ishock,ishock);
+    effective_Sigma_e = M_.Sigma_e(ishock,ishock);  % does not take heteroskedastic shocks into account
     [U,S] = svd(effective_Sigma_e);
+    % draw random shocks
     if opts.qmc
         opts.replic =2^(round(log2(opts.replic+1)))-1;
-        SHOCKS_ant =   qmc_sequence(forecast*effective_exo_nbr, int64(1), 1, opts.replic)';
+        SHOCKS_add = qmc_sequence(forecast_horizon*effective_exo_nbr, int64(1), 1, opts.replic);
     else
-        SHOCKS_ant = randn(forecast*effective_exo_nbr,opts.replic)';
+        SHOCKS_add = randn(forecast_horizon*effective_exo_nbr,opts.replic);
     end
-    zlin0=zeros(forecast,M_.endo_nbr,opts.replic);
-    zpiece0=zeros(forecast,M_.endo_nbr,opts.replic);
+    SHOCKS_add=reshape(SHOCKS_add,effective_exo_nbr,forecast_horizon,opts.replic);
+    z.linear=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
+    z.piecewise=NaN(forecast_horizon,M_.endo_nbr,opts.replic);
+    error_flag=true(opts.replic,1);
+    simul_SHOCKS=NaN(forecast_horizon,M_.exo_nbr,opts.replic);
     for iter=1:opts.replic
-
-        if ~isempty(SHOCKS0)
-            for j=1:length(SHOCKS0)
-                SHOCKS(:,inds(j))=SHOCKS1(j,:);
-            end
-        end
-
-        error_flagx=1;
-        % while error_flagx,
-        % SHOCKS=transpose(sqrt(diag(diag(effective_Sigma_e)))*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
-        SHOCKS=transpose(U*sqrt(S)*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
-        %             SHOCKS=transpose(U*sqrt(S)*randn(forecast,M_.exo_nbr)');   %realized shocks
-        options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
-        options_.occbin.simul.init_regime = opts.frcst_regimes;
-        options_.occbin.simul.init_binding_indicator = [];
-        options_.occbin.simul.exo_pos=ishock;
-        options_.occbin.simul.SHOCKS = SHOCKS;
+        options_.occbin.simul.SHOCKS = shocks_base+transpose(U*sqrt(S)*SHOCKS_add(:,:,iter));
         options_.occbin.simul.waitbar=0;
-        [~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-        zlin0(:,:,iter)=out.linear;
-        zpiece0(:,:,iter)=out.piecewise;
-        ys=out.ys;
+        [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
+        error_flag(iter)=out.error_flag;
+        if ~error_flag(iter)
+            z.linear(:,:,iter)=out.linear;
+            z.piecewise(:,:,iter)=out.piecewise;
             frcst_regime_history(iter,:)=out.regime_history;
             error_flag(iter)=out.error_flag;
-        error_flagx = error_flag(iter);
-        % end
-        simul_SHOCKS(:,:,iter) = SHOCKS;
-
-        if error_flag(iter) && debug_flag
-            %             display('no solution')
-
-%             keyboard;
-            save no_solution SHOCKS zlin0 zpiece0 iter frcst_regime_history
+            simul_SHOCKS(:,:,iter) = shocks_base;
+        else
+            if options_.debug
+                save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag','out','shocks_base')
+            end
         end
         dyn_waitbar(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
     end
     dyn_waitbar_close(h);
-    save temp zlin0 zpiece0 simul_SHOCKS error_flag
+    if options_.debug
+         save('Occbin_forecast_debug','simul_SHOCKS','z','iter','frcst_regime_history','error_flag')
+    end
     inx=find(error_flag==0);
-    zlin0=zlin0(:,:,inx);
-    zpiece0=zpiece0(:,:,inx);
-    zlin   = mean(zlin0,3);
-    zpiece = mean(zpiece0,3);
-    zpiecemin = min(zpiece0,[],3);
-    zpiecemax = max(zpiece0,[],3);
-    zlinmin = min(zlin0,[],3);
-    zlinmax = max(zlin0,[],3);
+    z.linear=z.linear(:,:,inx);
+    z.piecewise=z.piecewise(:,:,inx);
+    z.min.piecewise = min(z.piecewise,[],3);
+    z.max.piecewise = max(z.piecewise,[],3);
+    z.min.linear = min(z.linear,[],3);
+    z.max.linear = max(z.linear,[],3);
     
+    field_names={'linear','piecewise'};
+    post_mean=NaN(forecast_horizon,1);
+    post_median=NaN(forecast_horizon,1);
+    post_var=NaN(forecast_horizon,1);
+    hpd_interval=NaN(forecast_horizon,2);
+    post_deciles=NaN(forecast_horizon,9);
+    for field_iter=1:2
         for i=1:M_.endo_nbr
-        for j=1:forecast
-            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),options_.forecasts.conf_sig);
-        end
-        oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})=post_mean;
-        oo_.occbin.linear_forecast.Median.(M_.endo_names{i})=post_median;
-        oo_.occbin.linear_forecast.Var.(M_.endo_names{i})=post_var;
-        oo_.occbin.linear_forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
-        oo_.occbin.linear_forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
-        oo_.occbin.linear_forecast.Deciles.(M_.endo_names{i})=post_deciles;
-        oo_.occbin.linear_forecast.Min.(M_.endo_names{i})=zlinmin(:,i);
-        oo_.occbin.linear_forecast.Max.(M_.endo_names{i})=zlinmax(:,i);
-        for j=1:forecast
-            [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),options_.forecasts.conf_sig);
-        end
-        oo_.occbin.forecast.Mean.(M_.endo_names{i})=post_mean;
-        oo_.occbin.forecast.Median.(M_.endo_names{i})=post_median;
-        oo_.occbin.forecast.Var.(M_.endo_names{i})=post_var;
-        oo_.occbin.forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
-        oo_.occbin.forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
-        oo_.occbin.forecast.Deciles.(M_.endo_names{i})=post_deciles;
-        oo_.occbin.forecast.Min.(M_.endo_names{i})=zpiecemin(:,i);
-        oo_.occbin.forecast.Max.(M_.endo_names{i})=zpiecemax(:,i);
-        %   eval([M_.endo_names{i},'_ss=zdatass(i);']);
+            for j=1:forecast_horizon
+                [post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(z.(field_names{field_iter})(j,i,:)),options_.forecasts.conf_sig);
             end
-
-else
-    SHOCKS = zeros(forecast,M_.exo_nbr);
-    if ~isempty(SHOCKS0)
-        for j=1:length(SHOCKS0)
-            SHOCKS(:,inds(j))=SHOCKS1(j,:);
+            forecast.(field_names{field_iter}).Mean.(M_.endo_names{i})=post_mean;
+            forecast.(field_names{field_iter}).Median.(M_.endo_names{i})=post_median;
+            forecast.(field_names{field_iter}).Var.(M_.endo_names{i})=post_var;
+            forecast.(field_names{field_iter}).HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
+            forecast.(field_names{field_iter}).HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
+            forecast.(field_names{field_iter}).Deciles.(M_.endo_names{i})=post_deciles;
+            forecast.(field_names{field_iter}).Min.(M_.endo_names{i})=z.min.(field_names{field_iter})(:,i);
+            forecast.(field_names{field_iter}).Max.(M_.endo_names{i})=z.max.(field_names{field_iter})(:,i);
         end
     end
-    options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
-    options_.occbin.simul.init_regime = opts.frcst_regimes;
-    options_.occbin.simul.init_violvecbool = [];
+else
     options_.occbin.simul.irfshock = M_.exo_names;
-    options_.occbin.simul.SHOCKS = SHOCKS;
-    [~, out] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+    options_.occbin.simul.SHOCKS = shocks_base;
+    [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
+    error_flag=out.error_flag;
+    if out.error_flag
+        fprintf('occbin.forecast: forecast simulation failed.')
+        return;
+    end
 
-    zlin=out.linear;
-    zpiece=out.piecewise;
     frcst_regime_history=out.regime_history;
     error_flag=out.error_flag;
     for i=1:M_.endo_nbr
-        oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})= zlin(:,i);
-        oo_.occbin.forecast.Mean.(M_.endo_names{i})= zpiece(:,i);
-        oo_.occbin.forecast.HPDinf.(M_.endo_names{i})= nan;
-        oo_.occbin.forecast.HPDsup.(M_.endo_names{i})= nan;
+        forecast.linear.Mean.(M_.endo_names{i})= out.linear(:,i);
+        forecast.piecewise.Mean.(M_.endo_names{i})= out.piecewise(:,i);
     end
-
 end
 
-oo_.occbin.forecast.regimes=frcst_regime_history;
-
+forecast.regimes=frcst_regime_history;
\ No newline at end of file

matlab/+occbin/irf.m

-function [oo_] = irf(M_,oo_,options_)
+function irfs = irf(M_,oo_,options_)
+% irfs = irf(M_,oo_,options_)
+% Calls a minimizer
+%
+% INPUTS
+% - M_                  [structure]     Matlab's structure describing the model
+% - oo_                 [structure]     Matlab's structure containing the results
+% - options_            [structure]     Matlab's structure describing the current options
+%
+% OUTPUTS
+% - irfs                [structure]     IRF results
+%
+% 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/>.
 
 shocknames = options_.occbin.irf.exo_names;
-shocksigns = options_.occbin.irf.shocksigns;
+shocksigns = options_.occbin.irf.shocksigns; %'pos','neg'
 shocksize = options_.occbin.irf.shocksize;
-t0 = options_.occbin.irf.t0;
+t_0 = options_.occbin.irf.t0;
+
+%% set simulation options based on IRF options
 options_.occbin.simul.init_regime = options_.occbin.irf.init_regime;
 options_.occbin.simul.check_ahead_periods = options_.occbin.irf.check_ahead_periods;
 options_.occbin.simul.maxit = options_.occbin.irf.maxit;
 options_.occbin.simul.periods = options_.irf;
 
-% Run inital conditions + other shocks
-if t0 == 0
-    shocks0= zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
+%% Run initial conditions + other shocks
+if t_0 == 0
+    shocks_base = zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
     options_.occbin.simul.endo_init = [];
 else
-    % girf conditional to smoothed states in t0 and shocks in t0+1
-    shocks0= [oo_.occbin.smoother.etahat(:,t0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
-    options_.occbin.simul.SHOCKS=shocks0;
-    options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
+    if ~isfield(oo_.occbin,'smoother')
+        error('occbin.irfs: smoother must be run before requesting GIRFs based on smoothed results')
+    end
+    % GIRF conditional on smoothed states in t_0 and shocks in t_0+1
+    shocks_base= [oo_.occbin.smoother.etahat(:,t_0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
+    options_.occbin.simul.SHOCKS=shocks_base;
+    options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
+end
+options_.occbin.simul.SHOCKS=shocks_base;
+
+[~, out_base] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+if out_base.error_flag
+    error('occbin.irfs: could not compute the solution')
 end
-options_.occbin.simul.SHOCKS=shocks0;
-[~, out0] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-zlin0 = out0.linear;
-zpiece0 = out0.piecewise;
 
-% Select shocks of interest
-jexo_all =zeros(size(shocknames,1),1);
+irfs.linear = struct();
+irfs.piecewise = struct();
 
+% Get indices of shocks of interest
+exo_index =zeros(size(shocknames,1),1);
 for i=1:length(shocknames)
-    jexo_all(i) = strmatch(shocknames{i},M_.exo_names,'exact');
+    exo_index(i) = strmatch(shocknames{i},M_.exo_names,'exact');
 end
 
-oo_.occbin.linear_irfs = struct();
-oo_.occbin.irfs = struct();
+% cs=get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
+% irf_shocks_indx = getIrfShocksIndx(M_, options_);
+
 % Set shock size
 if isempty(shocksize)
-    %         if isfield(oo_.posterior_mode.shocks_std,M_.exo_names{jexo})
-    shocksize = sqrt(diag(M_.Sigma_e(jexo_all,jexo_all))); %oo_.posterior_mode.shocks_std.(M_.exo_names{jexo});
+    shocksize = sqrt(diag(M_.Sigma_e(exo_index,exo_index)));
     if any(shocksize < 1.e-9)
         shocksize(shocksize < 1.e-9) = 0.01;
     end
 end
+
 if numel(shocksize)==1
     shocksize=repmat(shocksize,[length(shocknames),1]);
 end
 
 % Run IRFs
-for counter = 1:length(jexo_all)
-  
-    jexo = jexo_all(counter);
-    
-    if ~options_.noprint
-        fprintf('Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},counter,size(jexo_all,1));
+for sign_iter=1:length(shocksigns)
+    for IRF_counter = 1:length(exo_index)
+        jexo = exo_index(IRF_counter);
+        if ~options_.noprint && options_.debug
+            fprintf('occbin.irf: Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},IRF_counter,size(exo_index,1));
+        end
+        shocks1=shocks_base;
+        if ismember('pos',shocksigns{sign_iter})
+            shocks1(1,jexo)=shocks_base(1,jexo)+shocksize(IRF_counter);
+        elseif ismember('neg',shocksigns{sign_iter})
+            shocks1(1,jexo)=shocks_base(1,jexo)-shocksize(IRF_counter);
         end
-    
-    if ismember('pos',shocksigns)
-        % (+) shock
-        shocks1=shocks0;
-        shocks1(1,jexo)=shocks0(1,jexo)+shocksize(counter);
-        if t0 == 0
         options_.occbin.simul.SHOCKS=shocks1;
+        if t_0 == 0
             options_.occbin.simul.endo_init = [];
-            [~, out_pos] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
         else
-            options_.occbin.simul.SHOCKS=shocks1;
-            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
-            [~, out_pos] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
         end
-        if out_pos.error_flag
-            warning('Occbin error.')
-            return
+        [~, out_sim] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+        if out_sim.error_flag
+            warning('occbin.irfs: simulation failed')
+            skip
         end
-        zlin_pos = out_pos.linear;
-        zpiece_pos = out_pos.piecewise;
         % Substract inital conditions + other shocks
-        zlin_pos_diff       = zlin_pos-zlin0;
-        zpiece_pos_diff      = zpiece_pos-zpiece0;
-    end
-    
-    if ismember('neg',shocksigns)
-        % (-) shock
-        shocks_1=shocks0;
-        shocks_1(1,jexo)=shocks0(1,jexo)-shocksize(counter);
-        if t0 == 0
-            options_.occbin.simul.SHOCKS=shocks_1;
-            options_.occbin.simul.endo_init = [];
-            [~, out_neg] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-        else
-            options_.occbin.simul.SHOCKS=shocks_1;
-            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
-            [~, out_neg] = occbin.solver(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
-        end
-        if out_neg.error_flag
-            warning('Occbin error.')
-            return
-        end
-        zlin_neg    = out_neg.linear;
-        zpiece_neg  = out_neg.piecewise;
-        zlin_neg_diff    = zlin_neg-zlin0;
-        zpiece_neg_diff  = zpiece_neg-zpiece0;
-    end
-
-    % Save
-    if ~isfield(oo_.occbin,'linear_irfs')
-        oo_.occbin.linear_irfs = struct();
-    end
-    if ~isfield(oo_.occbin,'irfs')
-        oo_.occbin.irfs = struct();
-    end
-
-    for jendo=1:M_.endo_nbr
-%         oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])   = zpiece_pos (:,jendo);
-%         oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])  = zpiece_neg (:,jendo);
-%         oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])   = zlin_pos (:,jendo);
-%         oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])  = zlin_neg(:,jendo);
+        zdiff.linear.(shocksigns{sign_iter}) = out_sim.linear-out_base.linear;
+        zdiff.piecewise.(shocksigns{sign_iter}) = out_sim.piecewise-out_base.piecewise;
 
+        for j_endo=1:M_.endo_nbr
             if ismember('pos',shocksigns)
-        oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos'])   = zpiece_pos_diff (:,jendo);
-        oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos'])   = zlin_pos_diff (:,jendo);
+                irfs.piecewise.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}])   = zdiff.piecewise.(shocksigns{sign_iter})(:,j_endo);
+                irfs.linear.([M_.endo_names{j_endo} '_' M_.exo_names{jexo} '_' shocksigns{sign_iter}])   = zdiff.linear.(shocksigns{sign_iter})(:,j_endo);
             end
-        
-        if ismember('neg',shocksigns)
-            oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg'])  = zpiece_neg_diff (:,jendo);
-            oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg'])  = zlin_neg_diff (:,jendo);
         end
-        
-% %         
-%         oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])   = zpiece0(:,jendo);
-%         oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1'])  = zlin0(:,jendo);
-%         oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])   = zpiece0(:,jendo);
-%         oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1'])  = zlin0(:,jendo);
-
     end
 end
\ No newline at end of file
-
-
-
-end
-

matlab/+occbin/plot_irfs.m

-function  plot_irfs(M_,oo_,options_,irf3,irf4)
+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
 
-shocknames = options_.occbin.plot_irf.exo_names;
-simulname = options_.occbin.plot_irf.simulname;
-if isempty(simulname)
-   simulname_ = simulname;
+simul_name = options_.occbin.plot_irf.simulname;
+if isempty(simul_name)
+    save_name = simul_name;
 else
-   simulname_ = [ simulname '_' ];
+    save_name = [ simul_name '_' ];
 end
-vars_irf = options_.occbin.plot_irf.endo_names;
-endo_names_long = options_.occbin.plot_irf.endo_names_long;
-endo_scaling_factor = options_.occbin.plot_irf.endo_scaling_factor;
-length_irf = options_.occbin.plot_irf.tplot;
-if isempty(length_irf)
-    length_irf = options_.irf;
+
+if isempty(var_list)
+    var_list = M_.endo_names(1:M_.orig_endo_nbr);
 end
 
-irflocation_lin   = oo_.occbin.linear_irfs;
-irflocation_piece = oo_.occbin.irfs;
+[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;
-warning('off','all')
 
 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 i=1:size(shocknames,1)
-    itemp = strmatch(shocknames{i},M_.exo_names,'exact');
+for var_iter=1:size(shock_names,1)
+    itemp = strmatch(shock_names{var_iter},M_.exo_names,'exact');
     if isempty(itemp)
-        error(['Shock ',shocknames{i},' is not defined!'])
+        error(['Shock ',shock_names{var_iter},' is not defined!'])
     else
         iexo=[iexo, itemp];
     end
@@ -38,104 +78,102 @@ ncols     = options_.occbin.plot_irf.ncols;
 nrows     = options_.occbin.plot_irf.nrows;
 npan      = ncols*nrows;
 
-plot_grid = options_.occbin.plot_irf.grid;
 shocksigns = options_.occbin.plot_irf.shocksigns;
-threshold = options_.occbin.plot_irf.threshold; 
-
-% Add steady_state
-if options_.occbin.plot_irf.add_steadystate
-    add_stst = options_.occbin.plot_irf.add_steadystate;
+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
-    add_stst = 0;
+       error('Unknown shock sign %s',shocksign);
     end
-for sss = 1:numel(shocksigns)
-    
-    shocksign = shocksigns{sss};
-    
-    for j=1:size(shocknames,1)
-        %shocknames = M_.exo_names{j};
         
+    for shock_iter=1:size(shock_names,1)
         j1   = 0;
         isub = 0;
         ifig = 0;
-
         % Variables
-        % ----------------------
-        for i = 1:length(vars_irf)
-
+        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, .
-                hfig = dyn_figure(options_.nodisplay,'name',['OccbinIRFs ' shocknames{j} ' ' simulname ' ' shocksign],'PaperPositionMode', 'auto','PaperType','A4','PaperOrientation','portrait','renderermode','auto','position',[10 10 950 650]);
+                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{i};
+                exofactor = endo_scaling_factor{var_iter};
             end
 
             subplot(nrows,ncols,isub)
-            irf_field   = strcat(vars_irf{i,1},'_',shocknames{j},'_',shocksign);
-            
-            irfvalues   = irflocation_lin.(irf_field);
-            if add_stst
-                irfvalues = irfvalues + get_mean(vars_irf{i,1});
-            end          
+            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;
-            
-            plot(irfvalues(1:steps_irf:length_irf)*exofactor,'linewidth',2);
-            hold on
-
-            irfvalues   = irflocation_piece.(irf_field);  
-            if add_stst
-                irfvalues = irfvalues + get_mean(vars_irf{i,1});
+                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
-            irfvalues(abs(irfvalues) <threshold) = 0;
-            plot(irfvalues(1:steps_irf:length_irf)*exofactor,'r--','linewidth',2);
-            
             hold on
-            plot(irfvalues(1:steps_irf:length_irf)*0,'k-','linewidth',1.5);
-            % Optional additional IRFs
-            if nargin > 10
-                  irfvalues   = irf3.(irf_field) ;
+
+            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;
-                  plot(irfvalues(1:steps_irf:length_irf)*exofactor,'k:','linewidth',2);
+                if options_.occbin.plot_irf.add_steadystate
+                    irfvalues = irfvalues + get_mean(vars_irf{var_iter,1});
                 end
-            if nargin > 11
-                  irfvalues   = irf4.(irf_field)   ;   
-                  irfvalues(abs(irfvalues) <threshold) = 0;
-                  plot(irfvalues(1:steps_irf:length_irf)*exofactor,'g-.','linewidth',2);
+                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 plot_grid
+            if options_.occbin.plot_irf.grid
                 grid on
             end     
-
-            xlim([1 (length_irf/steps_irf)]);
-            
-            % title
-            if isempty(endo_names_long)
-                title(regexprep(vars_irf{i},'_',' '))
+            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(endo_names_long{i})
+                title(M_.endo_names{i_var(var_iter)},'Interpreter','none')
             end
-
             % Annotation Box + save figure
             % ----------------------
-            if mod(j1,npan)==0 || (mod(j1,npan)~=0 && i==length(vars_irf))
+            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_',simulname_,shocknames{j},'_',shocksign,'_',int2str(ifig)],options_.nodisplay,options_.graph_format);
+                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
 
-warning('on','all')
+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;
@@ -13,9 +35,15 @@ else
 end    
 
 GraphDirectoryName = CheckPath('graphs',M_.dname);
+fhandle = dyn_figure(options_.nodisplay,'Name',[M_.fname ': OccBin regimes']);
 
-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)
@@ -36,12 +64,23 @@ for k=1:nconstr
             end
         end
     end
-    title(['regime ' int2str(k)])
-    xlabel('historic period')
-    ylabel('regime expected start')
+    xlim([1 nperiods])
+    title(['Regime ' int2str(k)])
+    xlabel('Historic period')
+    ylabel('Regime: expected start')
 end
-annotation('textbox',[.25,0,.15,.05],'String','Unbinding','Color','blue');
-annotation('textbox',[.65,0,.15,.05],'String','Binding','Color','red');
-
+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

@@ -55,9 +55,7 @@ if ismember(flag,{'forecast','all'})
     options_occbin_.forecast.maxit=30;
     options_occbin_.forecast.qmc=0;
     options_occbin_.forecast.replic=0;
-    options_occbin_.forecast.sepath=0;    
     options_occbin_.forecast.SHOCKS0=[];
-    options_occbin_.forecast.treepath=1; % number of branches
 end
 
 if ismember(flag,{'irf','all'})    
@@ -98,17 +96,12 @@ end
 
 if ismember(flag,{'plot_irf','all'})
     options_occbin_.plot_irf.add_steadystate = 0;
-    options_occbin_.plot_irf.exo_names = [];
-    options_occbin_.plot_irf.endo_names = M_.endo_names;
-    options_occbin_.plot_irf.endo_names_long = [];
     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='';
-    options_occbin_.plot_irf.threshold = 10^-6;
-    options_occbin_.plot_irf.tplot = [];
 end
 
 if ismember(flag,{'plot_shock_decomp','all'})

matlab/+pac/+update/parameters.m

@@ -13,7 +13,7 @@ function M_ = parameters(pacname, M_, oo_, verbose)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2018-2023 Dynare Team
+% Copyright © 2018-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -78,6 +78,11 @@ else
     numberofcomponents = 0;
 end
 
+% Makes no sense to have a composite PAC target with a trend component auxiliary model (where all the variables are non stationnary)
+if numberofcomponents>0 && strcmp(pacmodel.auxiliary_model_type, 'trend_component')
+    error('Composite PAC target not allowed with trend component model.')
+end
+
 % Build the vector of PAC parameters (ECM parameter + autoregressive parameters).
 pacvalues = M_.params([pacmodel.ec.params; pacmodel.ar.params(1:pacmodel.max_lag)']);
 
@@ -90,7 +95,7 @@ if numberofcomponents
             % Find the auxiliary variables if any
             ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
             if isempty(ad)
-                error('Cannot find the trend variable in the Companion VAR/VECM model.')
+                error('Cannot find the trend variable in the Companion VAR model.')
             else
                 for j=1:length(ad)
                     auxinfo = M_.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(j)}));
@@ -105,7 +110,7 @@ if numberofcomponents
                 end
             end
             if isempty(id{i})
-                error('Cannot find the trend variable in the Companion VAR/VECM model.')
+                error('Cannot find the trend variable in the Companion VAR model.')
             end
         end
     end
@@ -115,7 +120,7 @@ else
         % Find the auxiliary variables if any
         ad = find(cell2mat(cellfun(@(x) isauxiliary(x, [8 10]), varmodel.list_of_variables_in_companion_var, 'UniformOutput', false)));
         if isempty(ad)
-            error('Cannot find the trend variable in the Companion VAR/VECM model.')
+            error('Cannot find the trend variable in the auxiliary VAR / Trend component model.')
         else
             for i=1:length(ad)
                 auxinfo = M_.aux_vars(get_aux_variable_id(varmodel.list_of_variables_in_companion_var{ad(i)}));
@@ -130,7 +135,7 @@ else
             end
         end
         if isempty(id{1})
-            error('Cannot find the trend variable in the Companion VAR/VECM model.')
+            error('Cannot find the trend variable in the auxiliary VAR / Trend component model.')
         end
     end
 end

matlab/collapse_figures_in_tabgroup.m

+function collapse_figures_in_tabgroup
+
+% Copyright © 2023 Eduard Benet Cerda
+% Copyright © 2024 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/>.
+
+% Create a new figure with results
+fig = uifigure(Name = 'Dynare Results');
+
+% Add a grid layout to make sure it spans the entire width
+g = uigridlayout(fig, [1,1], Padding = 0);
+
+% Add a tabgroup
+tg = uitabgroup(g);
+
+% Find all figures with Dynare Tag 
+f = findobj('-regexp','tag','dynare-figure');
+
+% Loop over all figures and reparent them to a tab. Avoid legends, they are
+% automatically tied.
+for j = 1 : numel(f)
+    t = uitab(tg);
+    types = arrayfun(@class, f(j).Children, 'UniformOutput', false);
+    idx = ismember(types, 'matlab.graphics.illustration.Legend'); % no need to reparent legends
+    set(f(j).Children(~idx),'Parent',t)
+    t.Title = f(j).Name;
+    delete(f(j))
+end
\ No newline at end of file

matlab/dyn_figure.m

@@ -12,7 +12,7 @@ function h = dyn_figure(nodisplay, varargin)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2012-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,7 +30,7 @@ function h = dyn_figure(nodisplay, varargin)
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
 if nodisplay
-    h = figure(varargin{:},'visible','off');
+    h = figure(varargin{:},'visible','off','tag','dynare-figure');
 else
-    h = figure(varargin{:});
+    h = figure(varargin{:},'tag','dynare-figure');
 end

matlab/estimation/dynare_estimation_1.m

@@ -493,15 +493,22 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
         end
         if ~issmc(options_)
             [error_flag, ~, options_]= metropolis_draw(1, options_, estim_params_, M_);
+        else
+            error_flag=false;
         end
         if ~(~isempty(options_.sub_draws) && options_.sub_draws==0)
             if options_.bayesian_irf
+                if ~issmc(options_)
                     if error_flag
                         error('%s: I cannot compute the posterior IRFs!',dispString)
                     end
                     oo_=PosteriorIRF('posterior',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,dispString);
+                else
+                    fprintf('%s: SMC does not yet support the bayesian_irf option. Skipping computation.\n',dispString);
+                end
             end
             if options_.moments_varendo
+                if ~issmc(options_)
                     if error_flag
                         error('%s: I cannot compute the posterior moments for the endogenous variables!',dispString)
                     end
@@ -525,8 +532,12 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
                         end
                     end
                     oo_ = compute_moments_varendo('posterior',options_,M_,oo_,estim_params_,var_list_);
+                else
+                    fprintf('%s: SMC does not yet support the moments_varendo option. Skipping computation.\n',dispString);
+                end
             end
             if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
+                if ~ishssmc(options_)
                     if error_flag
                         error('%s: I cannot compute the posterior statistics!',dispString)
                     end
@@ -535,9 +546,12 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) ||  (any(
                     else
                         error('%s: Particle Smoothers are not yet implemented.',dispString)
                     end
+                else
+                    fprintf('%s: SMC does not yet support the smoother and forecast options. Skipping computation.\n',dispString);
+                end
             end
         else
-                fprintf('%s: sub_draws was set to 0. Skipping posterior computations.',dispString);
+            fprintf('%s: sub_draws was set to 0. Skipping posterior computations.\n',dispString);
         end
         xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
         M_ = set_all_parameters(xparam1,estim_params_,M_);

matlab/estimation/selec_posterior_draws.m

@@ -26,7 +26,7 @@ function SampleAddress = selec_posterior_draws(M_,options_,dr,endo_steady_state,
 %   None.
 %
 
-% Copyright © 2006-2022 Dynare Team
+% Copyright © 2006-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -146,7 +146,7 @@ if info
             pdraws(linee,1) = {x2(SampleAddress(i,4),:)};
             if info==2
                 M_ = set_parameters_locally(M_,pdraws{i,1});
-                [dr,~,M_.params] = compute_decision_rules(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+                [dr,~,M_.params] = compute_decision_rules(M_,options_,dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
                 pdraws(linee,2) = { dr };
             end
             old_mhfile = mhfile;

matlab/estimation/smc/dsmh.m

-function dsmh(TargetFun, xparam1, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_)
-
-% Dynamic Striated Metropolis-Hastings algorithm.
-%
-% INPUTS
-%   o TargetFun  [char]     string specifying the name of the objective
-%                           function (posterior kernel).
-%   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
-%   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters.
-%   o dataset_              data structure
-%   o dataset_info          dataset info structure
-%   o options_              options structure
-%   o M_                    model structure
-%   o estim_params_         estimated parameters structure
-%   o bayestopt_            estimation options structure
-%   o oo_                   outputs structure
-%
-% SPECIAL REQUIREMENTS
-%   None.
-%
-% PARALLEL CONTEXT
-% The most computationally intensive part of this function may be executed
-% in parallel. The code suitable to be executed in
-% parallel on multi core or cluster machine (in general a 'for' cycle)
-% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
-%
-% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
-% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
-% to manage the parallel computations.
-%
-% This function was the first function to be parallelized. Later, other
-% functions have been parallelized using the same methodology.
-% Then the comments write here can be used for all the other pairs of
-% parallel functions and also for management functions.
-
-% 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/>.
-
-opts = options_.posterior_sampler_options.dsmh;
-
-lambda = exp(bsxfun(@minus,options_.posterior_sampler_options.dsmh.H,1:1:options_.posterior_sampler_options.dsmh.H)/(options_.posterior_sampler_options.dsmh.H-1)*log(options_.posterior_sampler_options.dsmh.lambda1));
-c = 0.055 ;
-MM = int64(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/10) ;
-
-% Step 0: Initialization of the sampler
-[param, tlogpost_iminus1, loglik, bayestopt_] = ...
-    smc_samplers_initialization(TargetFun, 'dsmh', opts.particles, mh_bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
-
-ESS = zeros(options_.posterior_sampler_options.dsmh.H,1) ;
-zhat = 1 ;
-
-% The DSMH starts here
-for i=2:options_.posterior_sampler_options.dsmh.H
-    disp('');
-    disp('Tempered iteration');
-    disp(i) ;
-    % Step 1: sort the densities and compute IS weigths
-    [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param) ;
-    [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS) ;
-    % Step 2: tune c_i
-    c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-    % Step 3: Metropolis step
-    [param,tlogpost_iminus1,loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-end
-
-weights = exp(loglik*(lambda(end)-lambda(end-1)));
-weights = weights/sum(weights);
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.particles);
-distrib_param = param(:,indx_resmpl);
-
-mean_xparam = mean(distrib_param,2);
-npar  = length(xparam1);
-lb95_xparam = zeros(npar,1) ;
-ub95_xparam = zeros(npar,1) ;
-for i=1:npar
-    temp = sortrows(distrib_param(i,:)') ;
-    lb95_xparam(i) = temp(0.025*options_.posterior_sampler_options.dsmh.particles) ;
-    ub95_xparam(i) = temp(0.975*options_.posterior_sampler_options.dsmh.particles) ;
-end
-
-TeX = options_.TeX;
-
-str = sprintf(' Param. \t Lower Bound (95%%) \t Mean \t Upper Bound (95%%)');
-for l=1:npar
-    name = get_the_name(l,TeX,M_,estim_params_,options_.varobs);
-    str = sprintf('%s\n %s \t\t %5.4f \t\t %7.5f \t\t %5.4f', str, name, lb95_xparam(l), mean_xparam(l), ub95_xparam(l));
-end
-disp(str)
-disp('')
-
-%% Plot parameters densities
-
-if TeX
-    fidTeX = fopen([M_.fname '_param_density.tex'],'w');
-    fprintf(fidTeX,'%% TeX eps-loader file generated by DSMH.m (Dynare).\n');
-    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
-    fprintf(fidTeX,' \n');
-end
-
-number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.
-bandwidth = 0;                    % Rule of thumb optimal bandwidth parameter.
-kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourier Transform approximation.
-
-plt = 1 ;
-%for plt = 1:nbplt,
-if TeX
-    NAMES = [];
-    TeXNAMES = [];
-end
-hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
-for k=1:npar %min(nstar,npar-(plt-1)*nstar)
-    subplot(ceil(sqrt(npar)),floor(sqrt(npar)),k)
-    %kk = (plt-1)*nstar+k;
-    [name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
-    optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',options_.posterior_sampler_options.dsmh.particles,bandwidth,kernel_function);
-    [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
-        options_.posterior_sampler_options.dsmh.particles,optimal_bandwidth,kernel_function);
-    plot(density(:,1),density(:,2));
-    hold on
-    if TeX
-        title(texname,'interpreter','latex')
-    else
-        title(name,'interpreter','none')
-    end
-    hold off
-    axis tight
-    drawnow
-end
-dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
-if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-    % TeX eps loader file
-    fprintf(fidTeX,'\\begin{figure}[H]\n');
-    fprintf(fidTeX,'\\centering \n');
-    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/floor(sqrt(npar)),1),M_.fname,int2str(plt));
-    fprintf(fidTeX,'\\caption{Parameter densities based on the Dynamic Striated Metropolis-Hastings algorithm.}');
-    fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
-    fprintf(fidTeX,'\\end{figure}\n');
-    fprintf(fidTeX,' \n');
-end
-%end
-
-function [tlogpost_iminus1,loglik,param] = sort_matrices(tlogpost_iminus1,loglik,param)
-[~,indx_ord] = sortrows(tlogpost_iminus1);
-tlogpost_iminus1 = tlogpost_iminus1(indx_ord);
-param = param(:,indx_ord);
-loglik = loglik(indx_ord);
-
-function [tlogpost_i,weights,zhat,ESS,Omegachol] = compute_IS_weights_and_moments(param,tlogpost_iminus1,loglik,lambda,i,zhat,ESS)
-if i==1
-    tlogpost_i = tlogpost_iminus1 + loglik*lambda(i);
-else
-    tlogpost_i = tlogpost_iminus1 + loglik*(lambda(i)-lambda(i-1));
-end
-weights = exp(tlogpost_i-tlogpost_iminus1);
-zhat = (mean(weights))*zhat ;
-weights = weights/sum(weights);
-ESS(i) = 1/sum(weights.^2);
-% estimates of mean and variance
-mu = param*weights;
-z = bsxfun(@minus,param,mu);
-Omega = z*diag(weights)*z';
-Omegachol = chol(Omega)';
-
-function c = tune_c(TargetFun,param,tlogpost_i,lambda,i,c,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-disp('tuning c_i...');
-disp('Initial value =');
-disp(c) ;
-npar = size(param,1);
-lower_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^5;
-upper_prob = (.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))^(1/5);
-stop=0 ;
-while stop==0
-    acpt = 0.0;
-    indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
-    param0 = param(:,indx_resmpl);
-    tlogpost0 = tlogpost_i(indx_resmpl);
-    for j=1:options_.posterior_sampler_options.dsmh.G
-        for l=1:options_.posterior_sampler_options.dsmh.K
-            validate = 0;
-            while validate == 0
-                candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
-                if all(candidate >= mh_bounds.lb) && all(candidate <= mh_bounds.ub)
-                    [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-                    if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
-                        validate = 1;
-                        if rand(1,1)<exp(tlogpostx-tlogpost0(j)) % accept
-                            acpt = acpt + 1/(options_.posterior_sampler_options.dsmh.G*options_.posterior_sampler_options.dsmh.K);
-                            param0(:,j)= candidate;
-                            tlogpost0(j) = tlogpostx;
-                        end
-                    end
-                end
-            end
-        end
-    end
-    disp('Acceptation rate =') ;
-    disp(acpt) ;
-    if options_.posterior_sampler_options.dsmh.alpha0<=acpt && acpt<=options_.posterior_sampler_options.dsmh.alpha1
-        disp('done!');
-        stop=1;
-    else
-        if acpt<lower_prob
-            c = c/5;
-        elseif lower_prob<=acpt && acpt<=upper_prob
-            c = c*log(.5*(options_.posterior_sampler_options.dsmh.alpha0+options_.posterior_sampler_options.dsmh.alpha1))/log(acpt);
-        else
-            c = 5*c;
-        end
-        disp('Trying with c= ') ;
-        disp(c)
-    end
-end
-
-function [out_param,out_tlogpost_iminus1,out_loglik] = mutation_DSMH(TargetFun,param,tlogpost_i,tlogpost_iminus1,loglik,lambda,i,c,MM,Omegachol,weights,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
-indx_levels = (1:1:MM-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM;
-npar = size(param,1) ;
-p = 1/(10*options_.posterior_sampler_options.dsmh.tau);
-disp('Metropolis step...');
-% build the dynamic grid of levels
-levels = [0.0;tlogpost_iminus1(indx_levels)];
-% initialize the outputs
-out_param = param;
-out_tlogpost_iminus1 = tlogpost_i;
-out_loglik = loglik;
-% resample and initialize the starting groups
-indx_resmpl = smc_resampling(weights,rand(1,1),options_.posterior_sampler_options.dsmh.G);
-param0 = param(:,indx_resmpl);
-tlogpost_iminus10 = tlogpost_iminus1(indx_resmpl);
-tlogpost_i0 = tlogpost_i(indx_resmpl);
-loglik0 = loglik(indx_resmpl);
-% Start the Metropolis
-for l=1:options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.tau
-    for j=1:options_.posterior_sampler_options.dsmh.G
-        u1 = rand(1,1);
-        u2 = rand(1,1);
-        if u1<p
-            k=1 ;
-            for m=1:MM-1
-                if levels(m)<=tlogpost_iminus10(j) && tlogpost_iminus10(j)<levels(m+1)
-                    k = m+1;
-                    break
-                end
-            end
-            indx = floor( (k-1)*options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM+1 + u2*(options_.posterior_sampler_options.dsmh.N*options_.posterior_sampler_options.dsmh.G/MM-1) );
-            if i==1
-                alp = (loglik(indx)-loglik0(j))*lambda(i);
-            else
-                alp = (loglik(indx)-loglik0(j))*(lambda(i)-lambda(i-1));
-            end
-            if u2<exp(alp)
-                param0(:,j) = param(:,indx);
-                tlogpost_i0(j) = tlogpost_i(indx);
-                loglik0(j) = loglik(indx);
-                tlogpost_iminus10(j) = tlogpost_iminus1(indx);
-            end
-        else
-            validate= 0;
-            while validate==0
-                candidate = param0(:,j) + sqrt(c)*Omegachol*randn(npar,1);
-                if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
-                    [tlogpostx,loglikx] = tempered_likelihood(TargetFun,candidate,lambda(i),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
-                    if isfinite(loglikx) % if returned log-density is not Inf or Nan (penalized value)
-                        validate = 1;
-                        if u2<exp(tlogpostx-tlogpost_i0(j)) % accept
-                            param0(:,j) = candidate;
-                            tlogpost_i0(j) = tlogpostx;
-                            loglik0(j) = loglikx;
-                            if i==1
-                                tlogpost_iminus10(j) = tlogpostx-loglikx*lambda(i);
-                            else
-                                tlogpost_iminus10(j) = tlogpostx-loglikx*(lambda(i)-lambda(i-1));
-                            end
-                        end
-                    end
-                end
-            end
-        end
-    end
-    if mod(l,options_.posterior_sampler_options.dsmh.tau)==0
-        rang = (l/options_.posterior_sampler_options.dsmh.tau-1)*options_.posterior_sampler_options.dsmh.G+1:l*options_.posterior_sampler_options.dsmh.G/options_.posterior_sampler_options.dsmh.tau;
-        out_param(:,rang) = param0;
-        out_tlogpost_iminus1(rang) = tlogpost_i0;
-        out_loglik(rang) = loglik0;
-    end
-end

matlab/forecast_graphs.m

@@ -76,7 +76,7 @@ for j= 1:nvar
             fprintf(fidTeX,' \n');
         end
         n_fig =n_fig+1;
-        eval(['hh_fig=dyn_figure(options_.nodisplay,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
+        hh_fig=dyn_figure(options_.nodisplay,'Name',['Forecasts (' int2str(n_fig) ')']);
         m = 1;
     end
     subplot(nr, nc, m);
@@ -138,7 +138,7 @@ if isfield(oo_.forecast,'HPDinf_ME')
                 fprintf(fidTeX,' \n');
             end
             n_fig =n_fig+1;
-            eval(['hh_fig=dyn_figure(options_.nodisplay,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
+            hh_fig=dyn_figure(options_.nodisplay,'Name',['Forecasts including ME (' int2str(n_fig) ')']);
             m = 1;
         end
         subplot(nr,nc,m);

matlab/kalman/likelihood/univariate_kalman_filter.m

@@ -84,7 +84,7 @@ function [LIK, lik,a,P] = univariate_kalman_filter(data_index,number_of_observat
 %   Second Edition, Ch. 6.4 + 7.2.5
 
 
-% Copyright © 2004-2021 Dynare Team
+% Copyright © 2004-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -211,11 +211,21 @@ while notsteady && t<=last %loop over t
             end
             a = a + Ki*prediction_error; %filtering according to (6.13) in DK (2012)
             P = P - PZ*Ki';              %filtering according to (6.13) in DK (2012)
+        else
+            if Fi<0
+                %pathological numerical case where variance is negative
+                if analytic_derivation
+                    LIK={NaN,DLIK,Hess};
+                else
+                    LIK = NaN;
+                end
+                return
             else
                 % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
                 % p. 157, DK (2012)
             end
         end
+    end
     if analytic_derivation
         if analytic_derivation==2
             [Da,DP,D2a,D2P] = univariate_computeDstate(k,a,P,T,Da,DP,DT,DOm,notsteady,D2a,D2P,D2T,D2Om);

matlab/kalman/likelihood/univariate_kalman_filter_d.m

@@ -87,7 +87,7 @@ function [dLIK, dlikk, a, Pstar, llik] = univariate_kalman_filter_d(data_index,
 %   Series Analysis by State Space Methods", Oxford University Press,
 %   Second Edition, Ch. 5, 6.4 + 7.2.5
 
-% Copyright © 2004-2021 Dynare Team
+% Copyright © 2004-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -154,11 +154,17 @@ while newRank && (t<=last)
             dlikk(s) = dlikk(s) + llik(s,d_index(i));
             a = a+Kstar*(prediction_error/Fstar);
             Pstar = Pstar-Kstar*(Kstar'/Fstar);
+        else
+            if Fstar<0 || Finf<0
+                %pathological numerical case where variance is negative
+                dLIK = NaN;
+                return
             else
                 % do nothing as a_{t,i+1}=a_{t,i} and P_{t,i+1}=P_{t,i}, see
                 % p. 157, DK (2012)
             end
         end
+    end
     if newRank
         oldRank = rank(Z*Pinf*Z',diffuse_kalman_tol);
     else

matlab/model_info.m

@@ -280,21 +280,11 @@ function print_line(names,var_index,lead_lag,M_)
     else
         aux_index=find([M_.aux_vars(:).endo_index]==var_index);
         aux_type=M_.aux_vars(aux_index).type;
-        if ~isfield(M_.aux_vars(aux_index),'orig_lead_lag') || isempty(M_.aux_vars(aux_index).orig_lead_lag)
-            if ismember(aux_type,[1,3])
-                str = subst_auxvar(var_index, -1, M_);
-            elseif ismember(aux_type,[0,2])
-                str = subst_auxvar(var_index, 1, M_);
-            else
         if lead_lag==0
             str = subst_auxvar(var_index, [], M_);
         else
             str = subst_auxvar(var_index, lead_lag, M_);
         end
-            end
-        else
-            str = subst_auxvar(var_index, M_.aux_vars(aux_index).orig_lead_lag, M_);
-        end
         aux_orig_expression=M_.aux_vars(aux_index).orig_expr;
         if isempty(aux_orig_expression)
             fprintf('%s\n',str);

matlab/moments/compute_variance_decomposition.m

@@ -56,7 +56,7 @@ else
         var_decomp(stationary_vars,i) = vx2;
         variance_sum_loop = variance_sum_loop +vx2; %track overall variance over shocks
     end
-    if ~options_.pruning && max(abs(variance_sum_loop-var_stationary)./var_stationary) > 1e-4
+    if ~options_.pruning && max(abs(variance_sum_loop-var_stationary)./var_stationary) > 1e-4 && max(abs(variance_sum_loop-var_stationary))>1e-7
         warning(['Aggregate variance and sum of variances by shocks ' ...
             'differ by more than 0.01 %'])
     end

matlab/optimal_policy/dyn_ramsey_static.m

@@ -25,7 +25,7 @@ function [steady_state, params, check] = dyn_ramsey_static(ys_init, exo_ss, M_,
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2003-2023 Dynare Team
+% Copyright © 2003-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -137,7 +137,7 @@ end
 % Compute the value of the Lagrange multipliers that minimizes the norm of the
 % residuals, given the other endogenous
 if options_.bytecode
-    res = bytecode('static', M_, options, xx, exo_ss, M_.params, 'evaluate');
+    res = bytecode('static', M_, options_, xx, exo_ss, M_.params, 'evaluate');
 else
     res = feval([M_.fname '.sparse.static_resid'], xx, exo_ss, M_.params);
 end
@@ -167,7 +167,7 @@ end
 function result = check_static_model(ys,exo_ss,M_,options_)
 result = false;
 if (options_.bytecode)
-    res = bytecode('static', M_, options, ys, exo_ss, M_.params, 'evaluate');
+    res = bytecode('static', M_, options_, ys, exo_ss, M_.params, 'evaluate');
 else
     res = feval([M_.fname '.sparse.static_resid'], ys, exo_ss, M_.params);
 end

matlab/pac-tools/hVectors.m

@@ -18,7 +18,7 @@ function [h, lrcp] = hVectors(params, H, auxmodel, kind, id)
 %    params(2:end-1) ⟶ Autoregressive parameters.
 %    params(end)     ⟶ Discount factor.
 
-% Copyright © 2018-2021 Dynare Team
+% Copyright © 2018-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -52,17 +52,17 @@ n = length(H);
 tmp = eye(n*m)-kron(G, transpose(H)); % inv(W2)
 
 switch kind
-  case 'll'
+  case 'll' % (A.84), page 28 in Brayton, Davis and Tulip (2000) ⟹ The target is stationary (level-level).
     h = A_1*A_b*((kron(iota(m, m), H))'*(tmp\kron(iota(m, m), iota(n, id))));
-  case 'dd'
+  case 'dd' % (A.79), page 26 in Brayton, Davis and Tulip (2000) ⟹ The target appears in first difference as a dependent variable in the auxiliary model.
     h = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), H')*(tmp\kron(iota(m, m), iota(n, id))));
-  case 'dl'
+  case 'dl' % (A.74), page 24 in Brayton, Davis and Tulip (2000) ⟹ The target appears in level as a dependent variable in the auxiliary model.
     h = A_1*A_b*(kron(iota(m, m)'*inv(eye(m)-G), (H'-eye(length(H))))*(tmp\kron(iota(m, m), iota(n, id))));
   otherwise
     error('Unknown kind value in PAC model.')
 end
 
-if nargin>1
+if nargout>1
     if isequal(kind, 'll')
         lrcp = NaN;
     else

matlab/perfect-foresight-models/perfect_foresight_solver.m

@@ -460,12 +460,12 @@ function [steady_success, endo_simul, exo_simul, steady_state, exo_steady_state]
 %   exo_steady_state [vector]    steady state of exogenous corresponding to the scenario (equal to the input if terminal steady state not recomputed)
 
 % Compute convex combination for the path of exogenous
-exo_simul = exoorig*share/shareorig + exobase*(1-share);
+exo_simul = exoorig*share/shareorig + exobase*(1-share/shareorig);
 
 % Compute convex combination for the initial condition
 % In most cases, the initial condition is a steady state and this does nothing
 % This is for cases when the initial condition is out of equilibrium
-endo_simul(:, initperiods) = share/shareorig*endoorig(:, initperiods)+(1-share)*endobase(:, initperiods);
+endo_simul(:, initperiods) = share/shareorig*endoorig(:, initperiods)+(1-share/shareorig)*endobase(:, initperiods);
 
 % If there is a permanent shock, ensure that the rescaled terminal condition is
 % a steady state (if the user asked for this recomputation, or if the original
@@ -518,7 +518,7 @@ if recompute_final_steady_state
     options_.markowitz = saved_steady_markowitz;
 else
     % The terminal condition is not a steady state, compute a convex combination
-    endo_simul(:, lastperiods) = share/shareorig*endoorig(:, lastperiods)+(1-share)*endobase(:, lastperiods);
+    endo_simul(:, lastperiods) = share/shareorig*endoorig(:, lastperiods)+(1-share/shareorig)*endobase(:, lastperiods);
 end
 
 

matlab/perfect-foresight-models/solve_block_decomposed_problem.m

@@ -14,7 +14,7 @@ function [y, success, maxerror, per_block_status] = solve_block_decomposed_probl
 %   maxerror         [double]    ∞-norm of the residual
 %   per_block_status [struct]    vector structure with per-block information about convergence
 
-% Copyright © 2020-2023 Dynare Team
+% Copyright © 2020-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,18 +33,21 @@ function [y, success, maxerror, per_block_status] = solve_block_decomposed_probl
 
 cutoff = 1e-15;
 
-if options_.stack_solve_algo==0
-    mthd='Sparse LU';
-elseif options_.stack_solve_algo==1 || options_.stack_solve_algo==6
-    mthd='LBJ';
-elseif options_.stack_solve_algo==2
-    mthd='GMRES';
-elseif options_.stack_solve_algo==3
-    mthd='BICGSTAB';
-elseif options_.stack_solve_algo==4
-    mthd='OPTIMPATH';
-else
-    mthd='UNKNOWN';
+switch options_.stack_solve_algo
+    case 0
+        mthd='Sparse LU on stacked system';
+    case {1,6}
+        mthd='LBJ with LU solver';
+    case 2
+        mthd='GMRES on stacked system';
+    case 3
+        mthd='BiCGStab on stacked system';
+    case 4
+        mthd='Sparse LU solver with optimal path length on stacked system';
+    case 7
+        mthd='Solver from solve_algo option on stacked system';
+    otherwise
+        error('Unsupported stack_solve_algo value')
 end
 if options_.verbosity
     printline(41)