diff --git a/matlab/+occbin/forecast.m b/matlab/+occbin/forecast.m
index 3cc855223134928d0fe1a29f177ee47ea8f873f4..480fc2067ed9c3930f2eaa8fdd5a77626db76944 100644
--- a/matlab/+occbin/forecast.m
+++ b/matlab/+occbin/forecast.m
@@ -1,143 +1,169 @@
-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
         end
-    elseif isreal(SHOCKS0)
-        SHOCKS1=SHOCKS0;
-        inds = 1:M_.exo_nbr;
+    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;
-        frcst_regime_history(iter,:)=out.regime_history;
+        [~, out] = occbin.solver(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
         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
+        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;
+            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);
-
-    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);
+    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_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
+            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
-        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);']);
     end
-
 else
-    SHOCKS = zeros(forecast,M_.exo_nbr);
-    if ~isempty(SHOCKS0)
-        for j=1:length(SHOCKS0)
-            SHOCKS(:,inds(j))=SHOCKS1(j,:);
-        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 = [];
     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
diff --git a/matlab/+occbin/irf.m b/matlab/+occbin/irf.m
index a94b42b9d957ce4b56b9161d5bd98337556890f5..6266f1b69eed24f8ad6432f27a091989a61f897b 100644
--- a/matlab/+occbin/irf.m
+++ b/matlab/+occbin/irf.m
@@ -1,140 +1,124 @@
-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));
-    end 
-    
-    if ismember('pos',shocksigns)
-        % (+) shock
-        shocks1=shocks0;
-        shocks1(1,jexo)=shocks0(1,jexo)+shocksize(counter);
-        if t0 == 0
-            options_.occbin.simul.SHOCKS=shocks1;
-            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);
+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
-        if out_pos.error_flag
-            warning('Occbin error.')
-            return
+        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
-        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.SHOCKS=shocks1;
+        if t_0 == 0
             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);
+            options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t_0);
         end
-        if out_neg.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_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
+        % Substract inital conditions + other shocks
+        zdiff.linear.(shocksigns{sign_iter}) = out_sim.linear-out_base.linear;
+        zdiff.piecewise.(shocksigns{sign_iter}) = out_sim.piecewise-out_base.piecewise;
 
-    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);
-   
-        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);
+        for j_endo=1:M_.endo_nbr
+            if ismember('pos',shocksigns)
+                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
         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
-
-
-
-end
-
+end
\ No newline at end of file
diff --git a/matlab/+occbin/plot_irfs.m b/matlab/+occbin/plot_irfs.m
index efe92854789c66490e625d6be5762b811807cef5..29c9c63c39039d49f8c6648e4d3aff828902591b 100644
--- a/matlab/+occbin/plot_irfs.m
+++ b/matlab/+occbin/plot_irfs.m
@@ -1,34 +1,74 @@
-function  plot_irfs(M_,oo_,options_,irf3,irf4)
-
-shocknames = options_.occbin.plot_irf.exo_names;
-simulname = options_.occbin.plot_irf.simulname;
-if isempty(simulname)
-   simulname_ = simulname;
-else 
-   simulname_ = [ simulname '_' ];
+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
-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;
+
+simul_name = options_.occbin.plot_irf.simulname;
+if isempty(simul_name)
+    save_name = simul_name;
+else
+    save_name = [ simul_name '_' ];
 end
 
-irflocation_lin   = oo_.occbin.linear_irfs;
-irflocation_piece = oo_.occbin.irfs;
+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; 
-warning('off','all')
+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 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;
-else
-    add_stst = 0;
-end 
-for sss = 1:numel(shocksigns)
-    
-    shocksign = shocksigns{sss};
-    
-    for j=1:size(shocknames,1)
-        %shocknames = M_.exo_names{j};
+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 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]);
+                % 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;      
-            
+            isub=isub+1;
             if isempty(endo_scaling_factor)
                 exofactor = 1;
-            else               
-                exofactor = endo_scaling_factor{i};
-            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          
-            
-            irfvalues(abs(irfvalues) <threshold) = 0;
-            
-            plot(irfvalues(1:steps_irf:length_irf)*exofactor,'linewidth',2);
-            hold on
+            else
+                exofactor = endo_scaling_factor{var_iter};
+            end
 
-            irfvalues   = irflocation_piece.(irf_field);  
-            if add_stst
-                irfvalues = irfvalues + get_mean(vars_irf{i,1});
+            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
-            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) ;
-                  irfvalues(abs(irfvalues) <threshold) = 0;
-                  plot(irfvalues(1:steps_irf:length_irf)*exofactor,'k:','linewidth',2);
-            end
-            if nargin > 11
-                  irfvalues   = irf4.(irf_field)   ;   
-                  irfvalues(abs(irfvalues) <threshold) = 0;
-                  plot(irfvalues(1:steps_irf:length_irf)*exofactor,'g-.','linewidth',2);
-            end
-      
 
-            if plot_grid
-                grid 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
 
-            xlim([1 (length_irf/steps_irf)]);
-            
-            % title
-            if isempty(endo_names_long)
-                title(regexprep(vars_irf{i},'_',' '))
+            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(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
 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
diff --git a/matlab/+occbin/plot_regimes.m b/matlab/+occbin/plot_regimes.m
index 24e4aeb0c34c778cb95301af23995084d42ece93..a9fd4932be6577a146ba5aa0225e6a6c084b6bfd 100644
--- a/matlab/+occbin/plot_regimes.m
+++ b/matlab/+occbin/plot_regimes.m
@@ -1,4 +1,26 @@
 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
diff --git a/matlab/+occbin/set_default_options.m b/matlab/+occbin/set_default_options.m
index 99a1e2e1f8af76bff9413dfba5c886efe9f0e820..b720160c5392e6bcbba83ba096921b5141cf4d5d 100644
--- a/matlab/+occbin/set_default_options.m
+++ b/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'})
diff --git a/tests/occbin/filter/NKM.mod b/tests/occbin/filter/NKM.mod
index 325e7740262a022936d1cee23b9f4c297b77e296..1b552290f4dce06214ae6430c3b87fdea665d224 100644
--- a/tests/occbin/filter/NKM.mod
+++ b/tests/occbin/filter/NKM.mod
@@ -13,66 +13,66 @@
 
 // ----------------- Defintions -----------------------------------------//
 var        
-    c          //1  Consumption 
-    n          //2  Labor
-    y          //5  Output
-    yf         //6  Final goods       
-    yg         //11 Output growth gap
-    w          //12 Real wage rate
-    wf         //13 Flexible real wage
-    pigap      //15 Inflation rate -> pi(t)/pibar = pigap
-    inom    ${i^{nom}}$       //16 Nominal interest rate
-    inomnot    //17 Notional interest rate
-    mc         //19 Real marginal cost
-    lam  ${\lambda}$      //20 Inverse marginal utility of wealth  
-    g          //21 Growth shock       
-    s          //22 Risk premium shock
-    mp         //23 Monetary policy shock    
-    pi   ${\pi}$   //24 Observed inflation    
+    c       $c$         (long_name='Consumption')
+    n                   (long_name='Labor')
+    y       $y$         (long_name='Output')
+    yf                  (long_name='Final goods')
+    yg                  (long_name='Output growth gap')
+    w                   (long_name='Real wage rate')
+    wf                  (long_name='Flexible real wage')
+    pigap               (long_name='Inflation rate -> pi(t)/pibar = pigap')
+    inom    ${i^{nom}}$ (long_name='Nominal interest rate')
+    inomnot             (long_name='Notional interest rate')
+    mc                  (long_name='Real marginal cost')
+    lam  ${\lambda}$    (long_name='Inverse marginal utility of wealth')
+    g                   (long_name='Growth shock')
+    s                   (long_name='Risk premium shock')
+    mp                  (long_name='Monetary policy shock')
+    pi   ${\pi}$        (long_name='Observed inflation')
     @#if !(small_model)
-        x          //3  Investment
-        k          //4  Capital    
-        u          //7  Utilization cost
-        ups        //8  Utilization choice
-        wg         //9  Real wage growth gap    
-        xg         //10 Investment growth
-        rk         //14 Real rental rate
-        q          //18 Tobins q   
+        x               (long_name='Investment')
+        k               (long_name='Capital')
+        u               (long_name='Utilization cost')
+        ups             (long_name='Utilization choice')
+        wg              (long_name='Real wage growth gap')    
+        xg              (long_name='Investment growth')
+        rk              (long_name='Real rental rate')
+        q               (long_name='Tobins q')
     @#endif
 ;
 varexo          
-    epsg    ${\varepsilon_g}$   // Productivity growth shock
-    epsi       // Notional interest rate shock
-    epss       // Risk premium shock
+    epsg    ${\varepsilon_g}$   (long_name='Productivity growth shock')
+    epsi                        (long_name='Notional interest rate shock')
+    epss                        (long_name='Risk premium shock')
 ;        
 parameters
     // Calibrated Parameters    
-    beta    $\beta$    // Discount factor
-    chi         // Labor disutility scale
-    thetap      // Elasticity of subs. between intermediate goods
-    thetaw      // Elasticity of subs. between labor types
-    nbar        // Steady state labor
-    eta         // Inverse frish elasticity of labor supply
-    delta       // Depreciation
-    alpha       // Capital share
-    gbar        // Mean growth rate
-    pibar       // Inflation target
-    inombar     // Steady gross nom interest rate    
-    inomlb      // Effective lower bound on gross nominal interest rate    
-    sbar        // Average risk premium
+    beta    $\beta$    (long_name='Discount factor')
+    chi                 (long_name='Labor disutility scale')
+    thetap              (long_name='Elasticity of subs. between intermediate goods')
+    thetaw              (long_name='Elasticity of subs. between labor types')
+    nbar                (long_name='Steady state labor')
+    eta                 (long_name='Inverse frish elasticity of labor supply')
+    delta               (long_name='Depreciation')
+    alpha               (long_name='Capital share')
+    gbar                (long_name='Mean growth rate')
+    pibar               (long_name='Inflation target')
+    inombar             (long_name='Steady gross nom interest rate')
+    inomlb              (long_name='Effective lower bound on gross nominal interest rate')
+    sbar                (long_name='Average risk premium')
     // Parameters for DGP and Estimated parameters
-    varphip     // Rotemberg price adjustment cost
-    varphiw     // Rotemberg wage adjustment cost
-    h          	// Habit persistence
-    rhos        // Persistence
-    rhoi    	// Persistence
-    sigz        // Standard deviation technology
-    sigs        // Standard deviation risk premia
-    sigi        // Standard deviation mon pol
-    phipi       // Inflation responsiveness
-    phiy        // Output responsiveness
-    nu          // Investment adjustment cost
-    sigups      // Utilization    
+    varphip             (long_name='Rotemberg price adjustment cost')
+    varphiw             (long_name='Rotemberg wage adjustment cost')
+    h          	        (long_name='Habit persistence')
+    rhos                (long_name='Persistence')
+    rhoi    	        (long_name='Persistence')
+    sigz                (long_name='Standard deviation technology')
+    sigs                (long_name='Standard deviation risk premia')
+    sigi                (long_name='Standard deviation mon pol')
+    phipi               (long_name='Inflation responsiveness')
+    phiy                (long_name='Output responsiveness')
+    nu                  (long_name='Investment adjustment cost')
+    sigups              (long_name='Utilization')
  ;
 
 
@@ -316,36 +316,34 @@ varobs yg inom pi;
 
     // forecast starting from period 42, zero shocks (default)
     smoother2histval(period=42);
-    [oo, error_flag] = occbin.forecast(options_,M_,oo_,8);
+    [forecast, error_flag] = occbin.forecast(options_,M_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,8);
     // forecast with stochastic shocks
     options_.occbin.forecast.qmc=true;
     options_.occbin.forecast.replic=127;
-    [oo1, error_flag] = occbin.forecast(options_,M_,oo_,8);
+    [forecast1, error_flag] = occbin.forecast(options_,M_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,8);
+
+    figure('Name','OccBin: Forecasts')
+    subplot(2,1,1)
+    plot(1:8,forecast.piecewise.Mean.inom,'b-',1:8,forecast1.piecewise.Mean.inom,'r--')
+    subplot(2,1,2)
+    plot(1:8,forecast.piecewise.Mean.y,'b-',1:8,forecast1.piecewise.Mean.y,'r--')
 
     // GIRF given states in 42 and shocks in 43
     t0=42;
     options_.occbin.irf.exo_names=M_.exo_names;
     options_.occbin.irf.t0=t0;
-    oo_ = occbin.irf(M_,oo_,options_);
+    oo_.occbin.irfs = occbin.irf(M_,oo_,options_);
 
-    vars_irf = {
-    'c', 'consumption'    
-    'n', 'labor'    
-    'y', 'output'    
-    'pigap', 'inflation rate'
-    'inom', 'interest rate'  
-    'inomnot', 'shadow rate'
-    };
+    var_list_ = {'c','n','y','pigap','inom','inomnot'};
 
-    options_.occbin.plot_irf.exo_names = M_.exo_names;
-    options_.occbin.plot_irf.endo_names = vars_irf(:,1);
-    options_.occbin.plot_irf.endo_names_long = vars_irf(:,2);
     // if you want to scale ...
     // options_occbin_.plot_irf.endo_scaling_factor = vars_irf(:,3);
     options_.occbin.plot_irf.simulname = ['t0_' int2str(t0)];
-    options_.occbin.plot_irf.tplot = min(40,options_.irf);
-    occbin.plot_irfs(M_,oo_,options_);
-
+    options_.irf=40;
+    occbin.plot_irfs(M_,oo_.occbin.irfs,options_,var_list_);
+    var_list_={};
+    options_.occbin.plot_irf.simulname = ['t0_' int2str(t0) '_full'];
+    occbin.plot_irfs(M_,oo_.occbin.irfs,options_,var_list_);
     oo0=oo_;
     // use smoother_redux
     estimation(
@@ -374,7 +372,7 @@ varobs yg inom pi;
             consider_all_endogenous,heteroskedastic_filter,filter_step_ahead=[1],smoothed_state_uncertainty);
             
     // show initial condition effect of IF
-    figure,
+    figure('Name','OccBin: Smoothed shocks')
     subplot(221)
     plot([oo0.SmoothedShocks.epsg oo_.SmoothedShocks.epsg]), title('epsg')
     subplot(222)
@@ -382,7 +380,7 @@ varobs yg inom pi;
     subplot(223)
     plot([oo0.SmoothedShocks.epss oo_.SmoothedShocks.epss]), title('epss')
     legend('PKF','IF')
-    figure,
+    figure('Name','OccBin: Smoothed Variables')
     subplot(221)
     plot([oo0.SmoothedVariables.inom oo_.SmoothedVariables.inom]), title('inom')
     subplot(222)