matlab/gsa/scatter_analysis.m → matlab/+gsa/scatter_analysis.m

-function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions)
+function scatter_analysis(lpmat, xdata, options_scatter, options_)
+% scatter_analysis(lpmat, xdata, options_scatter, options_)
+% Plot scatter plot analysis
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
@@ -6,7 +8,7 @@ function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions)
 %
 
 % Copyright © 2017 European Commission
-% Copyright © 2017 Dynare Team
+% Copyright © 2017-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -25,7 +27,7 @@ function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions)
 
 param_names = options_scatter.param_names;
 
-if DynareOptions.TeX
+if options_.TeX
     if ~isfield(options_scatter,'param_names_tex')
         param_names_tex = options_scatter.param_names;
     else
@@ -34,7 +36,6 @@ if DynareOptions.TeX
 end
 amcf_name = options_scatter.amcf_name;
 amcf_title = options_scatter.amcf_title;
-title = options_scatter.title;
 fname_ = options_scatter.fname_;
 xparam1=[];
 if isfield(options_scatter,'xparam1')
@@ -42,11 +43,15 @@ if isfield(options_scatter,'xparam1')
 end
 OutputDirectoryName = options_scatter.OutputDirectoryName;
 
-if ~DynareOptions.nograph
+if ~options_.nograph
     skipline()
     xx=[];
     if ~isempty(xparam1)
         xx=xparam1;
     end
-    scatter_plots(lpmat, xdata, param_names, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, DynareOptions)
+    if options_.TeX
+        gsa.scatter_plots(lpmat, xdata, param_names_tex, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
+    else 
+        gsa.scatter_plots(lpmat, xdata, param_names, '.', [fname_, '_', amcf_name], OutputDirectoryName, amcf_title, xx, options_)
+    end
 end

matlab/gsa/scatter_mcf.m → matlab/+gsa/scatter_mcf.m

-function  scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, DynareOptions, beha_name, non_beha_name)
+function  scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_, beha_name, non_beha_name, beha_name_latex, non_beha_name_latex)
+% scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_, beha_name, non_beha_name, beha_name_latex, non_beha_name_latex)
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
@@ -37,12 +38,6 @@ function  scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, D
 
 
 Z=[X;Y];
-[n,p] = size(X);
-% X = X - ones(n,1)*min(Z);
-% X = X ./ (ones(n,1)*max(Z));
-[n,p] = size(Y);
-% Y = Y - ones(n,1)*min(Z);
-% Y = Y ./ (ones(n,1)*max(Z));
 [n,p] = size(Z);
 clear Z;
 
@@ -52,8 +47,10 @@ if nargin >=3
 end
 
 if nargin<4 || isempty(plotsymbol)
-    if n*p<100, plotsymbol = 'o';
-    else plotsymbol = '.';
+    if n*p<100
+        plotsymbol = 'o';
+    else 
+        plotsymbol = '.';
     end
 end
 
@@ -82,9 +79,9 @@ end
 
 figtitle_tex=strrep(figtitle,'_','\_');
 
-fig_nam_=[fnam];
+fig_nam_=fnam;
 if ~nograph
-    hh_fig=dyn_figure(DynareOptions.nodisplay,'name',figtitle);
+    hh_fig=dyn_figure(options_.nodisplay,'name',figtitle);
 end
 
 bf = 0.1;
@@ -99,11 +96,10 @@ for i = 1:p
     for j = 1:p
         h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]);
         if i==j
-            h1=cumplot(X(:,j));
-            %             set(h1,'color',[0 0 1], 'linestyle','--','LineWidth',1.5)
+            h1=gsa.cumplot(X(:,j));
             set(h1,'color',[0 0 1],'LineWidth',1.5)
             hold on,
-            h2=cumplot(Y(:,j));
+            h2=gsa.cumplot(Y(:,j));
             set(h2,'color',[1 0 0],'LineWidth',1.5)
             if ~isempty(xparam1)
                 hold on, plot(xparam1([j j]),[0 1],'k--')
@@ -125,10 +121,10 @@ for i = 1:p
                 plot(X(:,i),X(:,j),[plotsymbol,'b'])
             end
             if ~isempty(xparam1)
-                hold on, plot(xparam1(i),xparam1(j),'s','MarkerFaceColor',[0 0.75 0],'MarkerEdgeColor',[0 0.75 0])
+                hold on 
+                plot(xparam1(i),xparam1(j),'s','MarkerFaceColor',[0 0.75 0],'MarkerEdgeColor',[0 0.75 0])
             end
             hold off;
-            %             axis([-0.1 1.1 -0.1 1.1])
             if i<p
                 set(gca,'YTickLabel',[],'YTick',[]);
             else
@@ -143,16 +139,26 @@ for i = 1:p
         end
         if i==1
             if nflag == 1
+                if options_.TeX
                     ylabel(vnames(j,:),'Rotation',45, ...
-                       'HorizontalAlignment','right','VerticalAlignment','middle');
+                        'HorizontalAlignment','right','VerticalAlignment','middle','Interpreter','latex');
+                else
+                    ylabel(vnames(j,:),'Rotation',45, ...
+                        'HorizontalAlignment','right','VerticalAlignment','middle','Interpreter','none');
+                end
             else
                 ylabel([num2str(j),' '],'Rotation',90)
             end
         end
         if j==1
             if nflag == 1
+                if options_.TeX
                     title(vnames(i,:),'Rotation',45, ...
-                      'HorizontalAlignment','left','VerticalAlignment','bottom')
+                      'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','latex')
+                else
+                    title(vnames(i,:),'Rotation',45, ...
+                      'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','none')
+                end
             else
                 title(num2str(i))
             end
@@ -161,13 +167,18 @@ for i = 1:p
     end
 end
 if ~isoctave
+    if options_.TeX
+        annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name_latex,'Color','Blue','horizontalalignment','center','interpreter','latex');
+        annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name_latex,'Color','Red','horizontalalignment','center','interpreter','latex');
+    else
         annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name,'Color','Blue','horizontalalignment','center','interpreter','none');
         annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name,'Color','Red','horizontalalignment','center','interpreter','none');
     end
+end
 
 if ~nograph
-    dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],DynareOptions.nodisplay,DynareOptions.graph_format);
-    if DynareOptions.TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
+    dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],options_.nodisplay,options_.graph_format);
+    if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fidTeX = fopen([dirname,'/',fig_nam_ '.tex'],'w');
         fprintf(fidTeX,'%% TeX eps-loader file generated by scatter_mcf.m (Dynare).\n');
         fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);

matlab/gsa/scatter_plots.m → matlab/+gsa/scatter_plots.m

-function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, DynareOptions)
+function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_)
+% scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, options_)
+% Pairwise scatter plots of the columns of x and y after Monte Carlo filtering
+% Inputs:
+%  - X              [double]    nxk matrix with columns containing behavioural sample
+%  - xp             [double]    mxk matrix with columns containing non-behavioural sample
+%  - vnames         [char]      vector of variable names (default = numeric labels 1,2,3 etc.)
+%  - plotsymbol     [char]      plt symbol (default = '.' for npts > 100, 'o' for npts < 100
+%  - fnam           [char]      figure name
+%  - dirname        [char]      directory name
+%  - figtitle       [char]      figure title
+%  - xparam1        [double]    parameter vector
+%  - options_       [struct]    option structure
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
-%
+
 
 % Copyright © 2017 European Commission
 % Copyright © 2017-2023 Dynare Team
@@ -23,23 +35,7 @@ function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1,
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-% PURPOSE: Pairwise scatter plots of the columns of x and y after
-% Monte Carlo filtering
-%---------------------------------------------------
-% USAGE:    scatter_mcf(x,y,vnames,pltsym,diagon)
-%        or scatter_mcf(x,y) which relies on defaults
-% where:
-%        x = an nxk matrix with columns containing behavioural sample
-%        y = an mxk matrix with columns containing non-behavioural sample
-%   vnames = a vector of variable names
-%            (default = numeric labels 1,2,3 etc.)
-%   pltsym = a plt symbol
-%            (default = '.' for npts > 100, 'o' for npts < 100
-
-
 [n,p] = size(X);
-% X = X - ones(n,1)*min(Z);
-% X = X ./ (ones(n,1)*max(Z));
 
 nflag = 0;
 if nargin >=3
@@ -47,8 +43,10 @@ if nargin >=3
 end
 
 if nargin<4 || isempty(plotsymbol)
-    if n*p<100, plotsymbol = 'o';
-    else plotsymbol = '.';
+    if n*p<100
+        plotsymbol = 'o';
+    else 
+        plotsymbol = '.';
     end
 end
 
@@ -58,7 +56,7 @@ end
 if nargin<6 || isempty(dirname)
     dirname='';
     nograph=1;
-    DynareOptions.nodisplay=0;
+    options_.nodisplay=0;
 else
     nograph=0;
 end
@@ -71,9 +69,9 @@ end
 
 figtitle_tex=strrep(figtitle,'_','\_');
 
-fig_nam_=[fnam];
+fig_nam_=fnam;
 
-hh_fig=dyn_figure(DynareOptions.nodisplay,'name',figtitle);
+hh_fig=dyn_figure(options_.nodisplay,'name',figtitle);
 set(hh_fig,'userdata',{X,xp})
 
 bf = 0.1;
@@ -88,9 +86,8 @@ for i = 1:p
     for j = 1:p
         h = axes('position',[fL(i),fL(p+1-j),ffl,ffl]);
         if i==j
-            h1=cumplot(X(:,j));
+            h1=gsa.cumplot(X(:,j));
             set(h,'Tag','cumplot')
-            %             set(h1,'color',[0 0 1], 'linestyle','--','LineWidth',1.5)
             set(h1,'color',[0 0 1],'LineWidth',1.5)
             if ~isempty(xparam1)
                 hold on, plot(xparam1([j j]),[0 1],'k--')
@@ -101,40 +98,14 @@ for i = 1:p
                 grid off
             end
             set(gca,'YTickLabel',[],'YTick',[]);
-        else
-            if j>i
-                plot(X(:,i),X(:,j),[plotsymbol,'b'])
         else
             plot(X(:,i),X(:,j),[plotsymbol,'b'])
-            end
             set(h,'Tag','scatter')
 
-            %%
-            if ~isoctave
-                % Define a context menu; it is not attached to anything
-                hcmenu = uicontextmenu('Callback','pick','Tag','Run viewer');
-                % Define callbacks for context menu
-                % items that change linestyle
-                hcb1 = 'scatter_callback';
-                % hcb2 = ['set(gco,''LineStyle'','':'')'];
-                % hcb3 = ['set(gco,''LineStyle'',''-'')'];
-                % % Define the context menu items and install their callbacks
-                item1 = uimenu(hcmenu,'Label','save','Callback',hcb1,'Tag','save params');
-                item2 = uimenu(hcmenu,'Label','eval','Callback',hcb1,'Tag','eval params');
-                % item3 = uimenu(hcmenu,'Label','solid','Callback',hcb3);
-                % Locate line objects
-                hlines = findall(h,'Type','line');
-                % Attach the context menu to each line
-                for line = 1:length(hlines)
-                    set(hlines(line),'uicontextmenu',hcmenu)
-                end
-            end
-            %%
             if ~isempty(xparam1)
                 hold on, plot(xparam1(i),xparam1(j),'s','MarkerFaceColor',[0 0.75 0],'MarkerEdgeColor',[0 0.75 0])
             end
             hold off;
-            %             axis([-0.1 1.1 -0.1 1.1])
             if i<p
                 set(gca,'YTickLabel',[],'YTick',[]);
             else
@@ -149,16 +120,26 @@ for i = 1:p
         end
         if i==1
             if nflag == 1
+                if options_.TeX
+                    ylabel(vnames(j,:),'Rotation',45,'interpreter','latex', ...
+                        'HorizontalAlignment','right','VerticalAlignment','middle');
+                else
                     ylabel(vnames(j,:),'Rotation',45,'interpreter','none', ...
                         'HorizontalAlignment','right','VerticalAlignment','middle');
+                end
             else
                 ylabel([num2str(j),' '],'Rotation',90)
             end
         end
         if j==1
             if nflag == 1
+                if options_.TeX
+                    title(vnames(i,:),'interpreter','latex','Rotation',45, ...
+                          'HorizontalAlignment','left','VerticalAlignment','bottom')
+                else
                     title(vnames(i,:),'interpreter','none','Rotation',45, ...
                           'HorizontalAlignment','left','VerticalAlignment','bottom')            
+                end
             else
                 title(num2str(i))
             end
@@ -166,14 +147,10 @@ for i = 1:p
         drawnow
     end
 end
-% if ~isoctave
-%     annotation('textbox', [0.1,0,0.35,0.05],'String', beha_name,'Color','Blue','horizontalalignment','center','interpreter','none');
-%     annotation('textbox', [0.55,0,0.35,0.05],'String', non_beha_name,'Color','Red','horizontalalignment','center','interpreter','none');
-% end
 
 if ~nograph
-    dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],DynareOptions.nodisplay,DynareOptions.graph_format);
-    if DynareOptions.TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
+    dyn_saveas(hh_fig,[dirname,filesep,fig_nam_],options_.nodisplay,options_.graph_format);
+    if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
         fidTeX = fopen([dirname,'/',fig_nam_ '.tex'],'w');
         fprintf(fidTeX,'%% TeX eps-loader file generated by scatter_plots.m (Dynare).\n');
         fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);

matlab/gsa/set_shocks_param.m → matlab/+gsa/set_shocks_param.m

-function set_shocks_param(xparam1)
-% function set_shocks_param(xparam1)
+function M_=set_shocks_param(M_,estim_params_,xparam1)
+% function M_=set_shocks_param(M_,estim_params_,xparam1)
 % Set the structural and measurement error variances and covariances
+% Inputs
+%  - M_                     [structure]     Matlab's structure describing the model
+%  - estim_params_          [structure]     characterizing parameters to be estimated
+%  - xparam1                [double]        parameter vector
+% Outputs:
+%  - M_                     [structure]     Matlab's structure describing the model
+%
+% Notes: closely follows set_all_parameters.m
 
-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -19,8 +27,6 @@ function set_shocks_param(xparam1)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global estim_params_ M_
-
 nvx = estim_params_.nvx;
 ncx = estim_params_.ncx;
 nvn = estim_params_.nvn;

matlab/gsa/gsa_skewness.m → matlab/+gsa/skewness.m

-function s=gsa_skewness(y)
+function s=skewness(y)
+% s=skewness(y)
+% Compute normalized skewness of y
+% Inputs:
+%  - y  [double]  input vector
+% Outputs:
+%  - s  [double]  standardized skewness 
 
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
 
 % Copyright © 2012 European Commission
-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -22,8 +28,6 @@ function s=gsa_skewness(y)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-% y=stand_(y);
-% s=mean(y.^3);
 m2=mean((y-mean(y)).^2);
 m3=mean((y-mean(y)).^3);
 s=m3/m2^1.5;
\ No newline at end of file

matlab/gsa/smirnov.m → matlab/+gsa/smirnov_test.m

-function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
+function [H,prob,d] = smirnov_test(x1 , x2 , alpha, iflag )
+% [H,prob,d] = smirnov_test(x1 , x2 , alpha, iflag )
 % Smirnov test for 2 distributions
-%   [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
-%
+
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
@@ -34,11 +34,16 @@ end
 % empirical cdfs.
 xmix= [x1;x2];
 bin = [-inf ; sort(xmix) ; inf];
-
+if isoctave
     ncount1 = histc(x1 , bin);
-ncount1 = ncount1(:);
+else
+    ncount1 = histcounts(x1 , bin);
+end
+if isoctave
     ncount2 = histc(x2 , bin);
-ncount2 = ncount2(:);
+else
+    ncount2 = histcounts(x2 , bin);
+end
 
 cum1  =  cumsum(ncount1)./sum(ncount1);
 cum1  =  cum1(1:end-1);

matlab/gsa/stab_map_.m → matlab/+gsa/stability_mapping.m

-function x0 = stab_map_(OutputDirectoryName,opt_gsa)
-%
-% function x0 = stab_map_(OutputDirectoryName)
-%
+function x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
+% x0 = stability_mapping(OutputDirectoryName,opt_gsa,M_,oo_,options_,bayestopt_,estim_params_)
 % Mapping of stability regions in the prior ranges applying
 % Monte Carlo filtering techniques.
 %
-% INPUTS (from opt_gsa structure)
+% Inputs
+%  - OutputDirectoryName    [string]        name of the output directory
+%  - opt_gsa                [structure]     GSA options structure
+%  - M_                     [structure]     Matlab's structure describing the model
+%  - oo_                    [structure]     Matlab's structure describing the results
+%  - options_               [structure]     Matlab's structure describing the current options
+%  - bayestopt_             [structure]     describing the priors
+%  - estim_params_          [structure]     characterizing parameters to be estimated
+%
+% Outputs:
+%  - x0                                 one parameter vector for which the model is stable.
+%
+%
+% Inputs from opt_gsa structure
 % Nsam = MC sample size
 % fload = 0 to run new MC; 1 to load prevoiusly generated analysis
 % alpha2 =  significance level for bivariate sensitivity analysis
@@ -16,8 +27,6 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa)
 %            _prior.mat   file
 %        = 0: sample from posterior ranges: sample saved in
 %            _mc.mat file
-% OUTPUT:
-% x0: one parameter vector for which the model is stable.
 %
 % GRAPHS
 % 1) Pdf's of marginal distributions under the stability (dotted
@@ -28,14 +37,14 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa)
 % 3) Bivariate plots of significant correlation patterns
 %  ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets
 %
-% USES qmc_sequence, stab_map_1, stab_map_2
+% USES qmc_sequence, gsa.stability_mapping_univariate, gsa.stability_mapping_bivariate
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
 
 % Copyright © 2012-2016 European Commission
-% Copyright © 2012-2018 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -52,11 +61,6 @@ function x0 = stab_map_(OutputDirectoryName,opt_gsa)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
-global bayestopt_ estim_params_ options_ oo_ M_
-
-% opt_gsa=options_.opt_gsa;
-
 Nsam   = opt_gsa.Nsam;
 fload  = opt_gsa.load_stab;
 alpha2 = opt_gsa.alpha2_stab;
@@ -84,9 +88,14 @@ nshock = nshock + estim_params_.ncn;
 lpmat0=zeros(Nsam,0);
 xparam1=[];
 
-[~,~,~,lb,ub,~] = set_prior(estim_params_,M_,options_); %Prepare bounds
+%% prepare prior bounds
+[~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); %Prepare bounds
 if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
     % Set prior bounds
+    if options_.prior_trunc==0
+        fprintf('\nstability_mapping: GSA with priors requires bounded support. Setting options_.prior_trunc=1e-10.\n')
+        options_.prior_trunc=1e-10;
+    end
     bounds = prior_bounds(bayestopt_, options_.prior_trunc);
     bounds.lb = max(bounds.lb,lb);
     bounds.ub = min(bounds.ub,ub);
@@ -109,15 +118,16 @@ options_mcf.pvalue_ks = pvalue_ks;
 options_mcf.pvalue_corr = pvalue_corr;
 options_mcf.alpha2 = alpha2;
 
+%% get LaTeX names
 name=cell(np,1);
 name_tex=cell(np,1);
 for jj=1:np
     if options_.TeX
-        [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
-        name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
+        [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_.varobs);
+        name_tex{jj,1} = param_name_tex_temp;
         name{jj,1} = param_name_temp;
     else
-        param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
+        param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_.varobs);
         name{jj,1} = param_name_temp;
     end
 end
@@ -130,19 +140,18 @@ options_mcf.fname_ = fname_;
 options_mcf.OutputDirectoryName = OutputDirectoryName;
 options_mcf.xparam1 = [];
 
-opt=options_;
 options_.periods=0;
 options_.nomoments=1;
 options_.irf=0;
 options_.noprint=1;
-if fload==0
+if fload==0 %run new MC
     if isfield(dr_,'ghx')
         egg=zeros(length(dr_.eigval),Nsam);
     end
     yys=zeros(length(dr_.ys),Nsam);
 
     if opt_gsa.morris == 1
-        [lpmat, OutFact] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
+        [lpmat] = gsa.Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
         lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
         Nsam=size(lpmat,1);
         lpmat0 = lpmat(:,1:nshock);
@@ -160,10 +169,9 @@ if fload==0
             for j=1:np
                 lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
             end
-
         end
     end
-    dummy=prior_draw_gsa(1); 
+    gsa.prior_draw(M_,bayestopt_,options_,estim_params_,1); %initialize
     if pprior
         for j=1:nshock
             if opt_gsa.morris~=1
@@ -180,16 +188,16 @@ if fload==0
                 lpmat(:,j)=lpmat(:,j).*(upper_bound-lower_bound)+lower_bound;
             end
         else
-            xx=prior_draw_gsa(0,[lpmat0 lpmat]);
+            xx=gsa.prior_draw(M_,bayestopt_,options_,estim_params_,0,[lpmat0 lpmat]);
             lpmat0=xx(:,1:nshock);
             lpmat=xx(:,nshock+1:end);
             clear xx;
         end
-    else
+    else %posterior analysis
         if neighborhood_width>0 && isempty(options_.mode_file)
             xparam1 = get_all_parameters(estim_params_,M_);
         else
-            eval(['load ' options_.mode_file '.mat;']);
+            load([options_.mode_file '.mat'],'hh','xparam1');
         end
         if neighborhood_width>0
             for j=1:nshock
@@ -218,8 +226,8 @@ if fload==0
             for j=1:Nsam*2
                 lnprior(j) = any(lp(j,:)'<=bounds.lb | lp(j,:)'>=bounds.ub);
             end
-            ireal=[1:2*Nsam];
-            ireal=ireal(find(lnprior==0));
+            ireal=1:2*Nsam;
+            ireal=ireal(lnprior==0);
             lp=lp(ireal,:);
             Nsam=min(Nsam, length(ireal));
             lpmat0=lp(1:Nsam,1:nshock);
@@ -229,9 +237,9 @@ if fload==0
     end
     %
     h = dyn_waitbar(0,'Please wait...');
-    istable=[1:Nsam];
+    istable=1:Nsam;
     jstab=0;
-    iunstable=[1:Nsam];
+    iunstable=1:Nsam;
     iindeterm=zeros(1,Nsam);
     iwrong=zeros(1,Nsam);
     inorestriction=zeros(1,Nsam);
@@ -239,12 +247,11 @@ if fload==0
     infox=zeros(Nsam,1);
     for j=1:Nsam
         M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
-        %try stoch_simul([]);
         try
             if ~isempty(options_.endogenous_prior_restrictions.moment)
-                [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
+                [Tt,Rr,~,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
             else
-                [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
+                [Tt,Rr,~,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
             end
             infox(j,1)=info(1);
             if infox(j,1)==0 && ~exist('T','var')
@@ -252,8 +259,7 @@ if fload==0
                 if prepSA
                     try
                         T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
-                    catch
-                        ME = lasterror();
+                    catch ME                         
                         if strcmp('MATLAB:nomem',ME.identifier)
                             prepSA=0;
                             disp('The model is too large for storing state space matrices ...')
@@ -296,7 +302,7 @@ if fload==0
                 nboth = dr_.nboth;
                 nfwrd = dr_.nfwrd;
             end
-            info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_);
+            info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
             infox(j,1)=info(1);
             if info(1)
                 inorestriction(j)=j;
@@ -326,22 +332,23 @@ if fload==0
         end
         ys_=real(dr_.ys);
         yys(:,j) = ys_;
-        ys_=yys(:,1);
+        if mod(j,3)
             dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
         end
+    end
     dyn_waitbar_close(h);
     if prepSA && jstab
         T=T(:,:,1:jstab);
     else
         T=[];
     end
-    istable=istable(find(istable));  % stable params ignoring restrictions
-    irestriction=irestriction(find(irestriction));  % stable params & restrictions OK
-    inorestriction=inorestriction(find(inorestriction));  % stable params violating restrictions
-    iunstable=iunstable(find(iunstable));   % violation of BK & restrictions & solution could not be found (whatever goes wrong)
-    iindeterm=iindeterm(find(iindeterm));  % indeterminacy
-    iwrong=iwrong(find(iwrong));  % dynare could not find solution
-    ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong,inorestriction]))); % explosive roots
+    istable=istable(istable~=0);  % stable params ignoring restrictions
+    irestriction=irestriction(irestriction~=0);  % stable params & restrictions OK
+    inorestriction=inorestriction(inorestriction~=0);  % stable params violating restrictions
+    iunstable=iunstable(iunstable~=0);   % violation of BK & restrictions & solution could not be found (whatever goes wrong)
+    iindeterm=iindeterm(iindeterm~=0);  % indeterminacy
+    iwrong=iwrong(iwrong~=0);  % dynare could not find solution
+    ixun=iunstable(~ismember(iunstable,[iindeterm,iwrong,inorestriction])); % explosive roots
 
     bkpprior.pshape=bayestopt_.pshape;
     bkpprior.p1=bayestopt_.p1;
@@ -358,8 +365,7 @@ if fload==0
                  'bkpprior','lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
                  'egg','yys','T','nspred','nboth','nfwrd','infox')
         end
-
-    else
+    else %~pprior
         if ~prepSA
             save([OutputDirectoryName filesep fname_ '_mc.mat'], ...
                  'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
@@ -370,19 +376,18 @@ if fload==0
                  'egg','yys','T','nspred','nboth','nfwrd','infox')
         end
     end
-else
+else %load old run
     if pprior
         filetoload=[OutputDirectoryName filesep fname_ '_prior.mat'];
     else
         filetoload=[OutputDirectoryName filesep fname_ '_mc.mat'];
     end
-    load(filetoload,'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun','egg','yys','nspred','nboth','nfwrd','infox')
+    load(filetoload,'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun','infox')
     Nsam = size(lpmat,1);
     if pprior==0 && ~isempty(options_.mode_file)
-        eval(['load ' options_.mode_file '.mat;']);
+        load([options_.mode_file '.mat'],'xparam1');
     end
 
-
     if prepSA && isempty(strmatch('T',who('-file', filetoload),'exact'))
         h = dyn_waitbar(0,'Please wait...');
         options_.periods=0;
@@ -394,30 +399,24 @@ else
         yys=NaN(length(ys_),ntrans);
         for j=1:ntrans
             M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
-            %stoch_simul([]);
-            [Tt,Rr,SteadyState,info,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
+            [~,~,~,~,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state,'restrict');
             if ~exist('T','var')
                 T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
             end
             dr_ = oo_.dr;
             T(:,:,j) = [dr_.ghx dr_.ghu];
-            if ~exist('nspred','var')
-                nspred = dr_.nspred; %size(dr_.ghx,2);
-                nboth = dr_.nboth;
-                nfwrd = dr_.nfwrd;
-            end
             ys_=real(dr_.ys);
             yys(:,j) = ys_;
-            ys_=yys(:,1);
+            if mod(j,3)
                 dyn_waitbar(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
             end
+        end
         dyn_waitbar_close(h);
         save(filetoload,'T','-append')
-    elseif prepSA
-        load(filetoload,'T')
     end
 end
 
+%% display and save output
 if pprior
     aunstname='prior_unstable'; aunsttitle='Prior StabMap: explosiveness of solution';
     aindname='prior_indeterm'; aindtitle='Prior StabMap: Indeterminacy';
@@ -437,17 +436,18 @@ delete([OutputDirectoryName,filesep,fname_,'_',aindname,'.*']);
 delete([OutputDirectoryName,filesep,fname_,'_',aunstname,'.*']);
 delete([OutputDirectoryName,filesep,fname_,'_',awrongname,'.*']);
 
-if length(iunstable)>0 || length(iwrong)>0
-    fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
-    fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(ixun) )/Nsam*100)
+fprintf('\nSensitivity Analysis: Stability mapping:\n')
+if ~isempty(iunstable) || ~isempty(iwrong)
+    fprintf('%4.1f%% of the prior support gives unique saddle-path solution.\n',length(istable)/Nsam*100)
+    fprintf('%4.1f%% of the prior support gives explosive dynamics.\n',(length(ixun) )/Nsam*100)
     if ~isempty(iindeterm)
-        fprintf(['%4.1f%% of the prior support gives indeterminacy.\n'],length(iindeterm)/Nsam*100)
+        fprintf('%4.1f%% of the prior support gives indeterminacy.\n',length(iindeterm)/Nsam*100)
     end
-    inorestriction = istable(find(~ismember(istable,irestriction))); % violation of prior restrictions
+    inorestriction = istable(~ismember(istable,irestriction)); % violation of prior restrictions
     if ~isempty(iwrong) || ~isempty(inorestriction)
         skipline()
         if any(infox==49)
-            fprintf(['%4.1f%% of the prior support violates prior restrictions.\n'],(length(inorestriction) )/Nsam*100)
+            fprintf('%4.1f%% of the prior support violates prior restrictions.\n',(length(inorestriction) )/Nsam*100)
         end
         if ~isempty(iwrong)
             skipline()
@@ -488,50 +488,66 @@ if length(iunstable)>0 || length(iwrong)>0
     end
     skipline()
     if length(iunstable)<Nsam || length(istable)>1
-        itot = [1:Nsam];
-        isolve = itot(find(~ismember(itot,iwrong))); % dynare could find a solution
+        itot = 1:Nsam;
+        isolve = itot(~ismember(itot,iwrong)); % dynare could find a solution
                                                      % Blanchard Kahn
         if neighborhood_width
             options_mcf.xparam1 = xparam1(nshock+1:end);
         end
-        itmp = itot(find(~ismember(itot,istable)));
+        itmp = itot(~ismember(itot,istable));
         options_mcf.amcf_name = asname;
         options_mcf.amcf_title = atitle;
         options_mcf.beha_title = 'unique Stable Saddle-Path';
         options_mcf.nobeha_title = 'NO unique Stable Saddle-Path';
+        if options_.TeX
+            options_mcf.beha_title_latex = 'unique Stable Saddle-Path';
+            options_mcf.nobeha_title_latex = 'NO unique Stable Saddle-Path';            
+        end
         options_mcf.title = 'unique solution';
-        mcf_analysis(lpmat, istable, itmp, options_mcf, options_);
+        gsa.monte_carlo_filtering_analysis(lpmat, istable, itmp, options_mcf, M_, options_, bayestopt_, estim_params_);
 
         if ~isempty(iindeterm)
-            itmp = isolve(find(~ismember(isolve,iindeterm)));
+            itmp = isolve(~ismember(isolve,iindeterm));
             options_mcf.amcf_name = aindname;
             options_mcf.amcf_title = aindtitle;
             options_mcf.beha_title = 'NO indeterminacy';
             options_mcf.nobeha_title = 'indeterminacy';
+            if options_.TeX
+                options_mcf.beha_title_latex = 'NO indeterminacy';
+                options_mcf.nobeha_title_latex = 'indeterminacy';
+            end
             options_mcf.title = 'indeterminacy';
-            mcf_analysis(lpmat, itmp, iindeterm, options_mcf, options_);
+            gsa.monte_carlo_filtering_analysis(lpmat, itmp, iindeterm, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
 
         if ~isempty(ixun)
-            itmp = isolve(find(~ismember(isolve,ixun)));
+            itmp = isolve(~ismember(isolve,ixun));
             options_mcf.amcf_name = aunstname;
             options_mcf.amcf_title = aunsttitle;
             options_mcf.beha_title = 'NO explosive solution';
             options_mcf.nobeha_title = 'explosive solution';
+            if options_.TeX
+                options_mcf.beha_title_latex = 'NO explosive solution';
+                options_mcf.nobeha_title_latex = 'explosive solution';
+            end
             options_mcf.title = 'instability';
-            mcf_analysis(lpmat, itmp, ixun, options_mcf, options_);
+            gsa.monte_carlo_filtering_analysis(lpmat, itmp, ixun, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
 
-        inorestriction = istable(find(~ismember(istable,irestriction))); % violation of prior restrictions
-        iwrong = iwrong(find(~ismember(iwrong,inorestriction))); % what went wrong beyond prior restrictions
+        inorestriction = istable(~ismember(istable,irestriction)); % violation of prior restrictions
+        iwrong = iwrong(~ismember(iwrong,inorestriction)); % what went wrong beyond prior restrictions
         if ~isempty(iwrong)
-            itmp = itot(find(~ismember(itot,iwrong)));
+            itmp = itot(~ismember(itot,iwrong));
             options_mcf.amcf_name = awrongname;
             options_mcf.amcf_title = awrongtitle;
             options_mcf.beha_title = 'NO inability to find a solution';
             options_mcf.nobeha_title = 'inability to find a solution';
+            if options_.TeX
+                options_mcf.beha_title_latex = 'NO inability to find a solution';
+                options_mcf.nobeha_title_latex = 'inability to find a solution';
+            end
             options_mcf.title = 'inability to find a solution';
-            mcf_analysis(lpmat, itmp, iwrong, options_mcf, options_);
+            gsa.monte_carlo_filtering_analysis(lpmat, itmp, iwrong, options_mcf, M_, options_, bayestopt_, estim_params_);
         end
 
         if ~isempty(irestriction)
@@ -543,11 +559,11 @@ if length(iunstable)>0 || length(iwrong)>0
             name_tex=cell(np,1);
             for jj=1:np
                 if options_.TeX
-                    [param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_);
-                    name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
+                    [param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_.varobs);
+                    name_tex{jj,1} = param_name_tex_temp;
                     name{jj,1} = param_name_temp;
                 else
-                    param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_);
+                    param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_.varobs);
                     name{jj,1} = param_name_temp;
                 end
             end
@@ -559,8 +575,12 @@ if length(iunstable)>0 || length(iwrong)>0
             options_mcf.amcf_title = acalibtitle;
             options_mcf.beha_title = 'prior IRF/moment calibration';
             options_mcf.nobeha_title = 'NO prior IRF/moment calibration';
+            if options_.TeX
+                options_mcf.beha_title_latex = 'prior IRF/moment calibration';
+                options_mcf.nobeha_title_latex = 'NO prior IRF/moment calibration';
+            end
             options_mcf.title = 'prior restrictions';
-            mcf_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, options_);
+            gsa.monte_carlo_filtering_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, M_, options_, bayestopt_, estim_params_);
             iok = irestriction(1);
             x0 = [lpmat0(iok,:)'; lpmat(iok,:)'];
         else
@@ -570,7 +590,7 @@ if length(iunstable)>0 || length(iwrong)>0
         end
 
         M_ = set_all_parameters(x0,estim_params_,M_);
-        [oo_.dr,info,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
+        [oo_.dr,~,M_.params] = resol(0,M_,options_,oo_.dr ,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
     else
         disp('All parameter values in the specified ranges are not acceptable!')
         x0=[];
@@ -580,15 +600,8 @@ else
     disp('and match prior IRF/moment restriction(s) if any!')
     x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
     x0 = [x0; lpmat(istable(1),:)'];
-
 end
+skipline(1);
 
 xparam1=x0;
 save([OutputDirectoryName filesep 'prior_ok.mat'],'xparam1');
-
-options_.periods=opt.periods;
-if isfield(opt,'nomoments')
-    options_.nomoments=opt.nomoments;
-end
-options_.irf=opt.irf;
-options_.noprint=opt.noprint;

matlab/gsa/stab_map_2.m → matlab/+gsa/stability_mapping_bivariate.m

-function indcorr = stab_map_2(x,alpha2, pvalue_crit, fnam, dirname,xparam1,figtitle)
-% function stab_map_2(x, alpha2, pvalue, fnam, dirname,xparam1)
+function indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, case_name_plain, case_name_latex, dirname,xparam1,figtitle,fig_caption_latex)
+% indcorr = stability_mapping_bivariate(x,alpha2, pvalue_crit, M_,options_,bayestopt_,estim_params_, fnam, fnam_latex, dirname,xparam1,figtitle,fig_caption_latex)
+% Inputs:
+%  - x
+%  - alpha2
+%  - pvalue_crit
+%  - M_                     [structure]     Matlab's structure describing the model
+%  - options_               [structure]     Matlab's structure describing the current options
+%  - bayestopt_             [structure]     describing the priors
+%  - estim_params_          [structure]     characterizing parameters to be estimated
+%  - fnam                   [string]        file name
+%  - dirnam                 [string]        directory name
+%  - xparam1                [double]        parameter vector
+%  - figtitle               [string]        figure title
 %
+% Output:
+%  - indcorr
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
@@ -22,29 +36,30 @@ function indcorr = stab_map_2(x,alpha2, pvalue_crit, fnam, dirname,xparam1,figti
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
-global bayestopt_ estim_params_ options_ oo_ M_
-
 npar=size(x,2);
-nsam=size(x,1);
 ishock= npar>estim_params_.np;
 nograph = options_.nograph;
-if nargin<4
-    fnam='';
+
+if nargin<8
+    case_name_plain='';
 end
-if nargin<5
+if nargin<9
+    case_name_latex=case_name_plain;
+end
+if nargin<10
     dirname='';
     nograph=1;
 end
-if nargin<6
+if nargin<11
     xparam1=[];
 end
-if nargin<7
-    figtitle=fnam;
+if nargin<12
+    figtitle=case_name_plain;
+end
+if nargin<13
+    fig_caption_latex=case_name_latex;
 end
 
-ys_ = oo_.dr.ys;
-dr_ = oo_.dr;
 fname_ = M_.fname;
 nshock = estim_params_.nvx;
 nshock = nshock + estim_params_.nvn;
@@ -53,9 +68,9 @@ nshock = nshock + estim_params_.ncn;
 
 [c0, pvalue] = corrcoef(x);
 c00=tril(c0,-1);
-fig_nam_=[fname_,'_',fnam,'_corr_'];
-fig_nam_tex_table=strrep([fnam,'_corr'],' ','_');
-fig_nam_=strrep(fig_nam_,' ','_');
+fig_nam_save=[fname_,'_',case_name_plain,'_corr_'];
+fig_nam_save=strrep(fig_nam_save,' ','_');
+fig_nam_tex_table_save=strrep([case_name_plain,'_corr'],' ','_');
 
 ifig=0;
 j2=0;
@@ -67,47 +82,43 @@ if ishock==0
 else
     npar=estim_params_.np+nshock;
 end
-title_string=['Correlation analysis for ',fnam];
-title_string_tex=['Correlation analysis for ',strrep(fnam,'_','\\_')];
+title_string=['Correlation analysis for ',case_name_plain];
+title_string_tex=['Correlation analysis for ',case_name_latex];
 
 indcorr = [];
 entry_iter=1;
 for j=1:npar
     i2=find(abs(c00(:,j))>alpha2);
-    if length(i2)>0
+    if ~isempty(i2)
         for jx=1:length(i2)
             if pvalue(j,i2(jx))<pvalue_crit
                 indcorr = [indcorr; [j i2(jx)]];
                 j2=j2+1;
                 if ishock
                     if options_.TeX
-                        [param_name_temp1, param_name_tex_temp1]= get_the_name(j,options_.TeX,M_,estim_params_,options_);
-                        param_name_tex_temp1 = strrep(param_name_tex_temp1,'$','');
-                        [param_name_temp2, param_name_tex_temp2]= get_the_name(i2(jx),options_.TeX,M_,estim_params_,options_);
-                        param_name_tex_temp2 = strrep(param_name_tex_temp2,'$','');
+                        [param_name_temp1, param_name_tex_temp1]= get_the_name(j,options_.TeX,M_,estim_params_,options_.varobs);
+                        [param_name_temp2, param_name_tex_temp2]= get_the_name(i2(jx),options_.TeX,M_,estim_params_,options_.varobs);
                         tmp_name=(['[',param_name_temp1,',',param_name_temp2,']']);
                         tmp_name_tex=(['[',param_name_tex_temp1,',',param_name_tex_temp2,']']);
                         name{entry_iter,1}=tmp_name;
-                        name_tex{entry_iter,1}=tmp_name_tex;
+                        name_tex{entry_iter,1}=strrep(tmp_name_tex,'$',''); %prevent $ inside of expression for table
                     else
-                        [param_name_temp1]= get_the_name(j,options_.TeX,M_,estim_params_,options_);
-                        [param_name_temp2]= get_the_name(i2(jx),options_.TeX,M_,estim_params_,options_);
+                        [param_name_temp1]= get_the_name(j,options_.TeX,M_,estim_params_,options_.varobs);
+                        [param_name_temp2]= get_the_name(i2(jx),options_.TeX,M_,estim_params_,options_.varobs);
                         tmp_name=(['[',param_name_temp1,',',param_name_temp2,']']);
                         name{entry_iter,1}=tmp_name;
                     end
                 else
                     if options_.TeX
-                        [param_name_temp1, param_name_tex_temp1]= get_the_name(j+nshock,options_.TeX,M_,estim_params_,options_);
-                        param_name_tex_temp1 = strrep(param_name_tex_temp1,'$','');
-                        [param_name_temp2, param_name_tex_temp2]= get_the_name(i2(jx)+nshock,options_.TeX,M_,estim_params_,options_);
-                        param_name_tex_temp2 = strrep(param_name_tex_temp2,'$','');
+                        [param_name_temp1, param_name_tex_temp1]= get_the_name(j+nshock,options_.TeX,M_,estim_params_,options_.varobs);
+                        [param_name_temp2, param_name_tex_temp2]= get_the_name(i2(jx)+nshock,options_.TeX,M_,estim_params_,options_.varobs);
                         tmp_name=(['[',param_name_temp1,',',param_name_temp2,']']);
                         tmp_name_tex=(['[',param_name_tex_temp1,',',param_name_tex_temp2,']']);
                         name{entry_iter,1}=tmp_name;
-                        name_tex{entry_iter,1}=tmp_name_tex;
+                        name_tex{entry_iter,1}=strrep(tmp_name_tex,'$',''); %prevent $ inside of expression for table
                     else
-                        [param_name_temp1]= get_the_name(j+nshock,options_.TeX,M_,estim_params_,options_);
-                        [param_name_temp2]= get_the_name(i2(jx)+nshock,options_.TeX,M_,estim_params_,options_);
+                        [param_name_temp1]= get_the_name(j+nshock,options_.TeX,M_,estim_params_,options_.varobs);
+                        [param_name_temp2]= get_the_name(i2(jx)+nshock,options_.TeX,M_,estim_params_,options_.varobs);
                         tmp_name=(['[',param_name_temp1,',',param_name_temp2,']']);
                         name{entry_iter,1}=tmp_name;
                     end
@@ -121,17 +132,10 @@ for j=1:npar
                         hh_fig=dyn_figure(options_.nodisplay,'name',[figtitle,' sample bivariate projection ', num2str(ifig)]);
                     end
                     subplot(3,4,j2-(ifig-1)*12)
-                    %             bar(c0(i2,j)),
-                    %             set(gca,'xticklabel',bayestopt_.name(i2)),
-                    %             set(gca,'xtick',[1:length(i2)])
-                    %plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
-                    %hold on,
                     plot(x(:,j),x(:,i2(jx)),'.')
                     if ~isempty(xparam1)
                         hold on, plot(xparam1(j),xparam1(i2(jx)),'ro')
                     end
-                    %             xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
-                    %             ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
                     if ishock
                         xlabel(bayestopt_.name{j},'interpreter','none'),
                         ylabel(bayestopt_.name{i2(jx)},'interpreter','none'),
@@ -141,16 +145,16 @@ for j=1:npar
                     end
                     title(['cc = ',num2str(c0(i2(jx),j))])
                     if (mod(j2,12)==0) && j2>0
-                        dyn_saveas(hh_fig,[dirname,filesep,fig_nam_,int2str(ifig)],options_.nodisplay,options_.graph_format);
+                        dyn_saveas(hh_fig,[dirname,filesep,fig_nam_save,int2str(ifig)],options_.nodisplay,options_.graph_format);
                         if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-                            fidTeX = fopen([dirname,filesep,fig_nam_,int2str(ifig),'.tex'],'w');
+                            fidTeX = fopen([dirname,filesep,fig_nam_save,int2str(ifig),'.tex'],'w');
                             fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_2.m (Dynare).\n');
                             fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                             fprintf(fidTeX,'\\begin{figure}[H]\n');
                             fprintf(fidTeX,'\\centering \n');
-                            fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s}\n',strrep([dirname,'/',fig_nam_,int2str(ifig)],'\','/'));
-                            fprintf(fidTeX,'\\caption{%s.}',[figtitle,' sample bivariate projection ', num2str(ifig)]);
-                            fprintf(fidTeX,'\\label{Fig:%s:%u}\n',fig_nam_,ifig);
+                            fprintf(fidTeX,'\\includegraphics[width=0.8\\textwidth]{%s}\n',strrep([dirname,'/',fig_nam_save,int2str(ifig)],'\','/'));
+                            fprintf(fidTeX,'\\caption{%s.}',[fig_caption_latex,' sample bivariate projection ', num2str(ifig)]);
+                            fprintf(fidTeX,'\\label{Fig:%s:%u}\n',fig_nam_save,ifig);
                             fprintf(fidTeX,'\\end{figure}\n\n');
                             fprintf(fidTeX,'%% End Of TeX file. \n');
                             fclose(fidTeX);
@@ -162,16 +166,16 @@ for j=1:npar
         end
     end
     if ~nograph && (j==(npar)) && j2>0 && (mod(j2,12)~=0)
-        dyn_saveas(hh_fig,[dirname,filesep,fig_nam_,int2str(ifig)],options_.nodisplay,options_.graph_format);
+        dyn_saveas(hh_fig,[dirname,filesep,fig_nam_save,int2str(ifig)],options_.nodisplay,options_.graph_format);
         if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
-            fidTeX = fopen([dirname,filesep,fig_nam_,int2str(ifig),'.tex'],'w');
+            fidTeX = fopen([dirname,filesep,fig_nam_save,int2str(ifig),'.tex'],'w');
             fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_2.m (Dynare).\n');
             fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
             fprintf(fidTeX,'\\begin{figure}[H]\n');
             fprintf(fidTeX,'\\centering \n');
-            fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s}\n',options_.figures.textwidth*min((j2-(ifig-1)*12)/3,1),strrep([dirname,'/',fig_nam_,int2str(ifig)],'\','/'));
-            fprintf(fidTeX,'\\caption{%s.}',[figtitle,' sample bivariate projection ', num2str(ifig)]);
-            fprintf(fidTeX,'\\label{Fig:%s:%u}\n',fig_nam_,ifig);
+            fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s}\n',options_.figures.textwidth*min((j2-(ifig-1)*12)/3,1),strrep([dirname,'/',fig_nam_save,int2str(ifig)],'\','/'));
+            fprintf(fidTeX,'\\caption{%s.}',[fig_caption_latex,' sample bivariate projection ', num2str(ifig)]);
+            fprintf(fidTeX,'\\label{Fig:%s:%u}\n',fig_nam_save,ifig);
             fprintf(fidTeX,'\\end{figure}\n\n');
             fprintf(fidTeX,'%% End Of TeX file. \n');
             fclose(fidTeX);
@@ -181,7 +185,7 @@ end
 
 if j2==0
     skipline();
-    disp(['No correlation term with pvalue <', num2str(pvalue_crit),' and |corr. coef.| >',num2str(alpha2),' found for ',fnam])
+    disp(['No correlation term with pvalue <', num2str(pvalue_crit),' and |corr. coef.| >',num2str(alpha2),' found for ',case_name_plain])
 else
     headers={'Parameters'; 'corrcoef'};
     if ~options_.noprint
@@ -189,7 +193,6 @@ else
     end
     dyntable(options_,title_string,headers, name, data_mat, 0, 7, 3);
     if options_.TeX
-        dyn_latex_table(M_, options_, title_string_tex, fig_nam_tex_table, headers, name_tex, data_mat, 0, 7, 3);
+        dyn_latex_table(M_, options_, title_string_tex, fig_nam_tex_table_save, headers, name_tex, data_mat, 0, 7, 3);
     end
 end
\ No newline at end of file
-%close all

matlab/gsa/stab_map_1.m → matlab/+gsa/stability_mapping_univariate.m

-function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname, pcrit, atitle)
-%function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname, pcrit)
-%
-% lpmat =  Monte Carlo matrix
-% ibehaviour = index of behavioural runs
-% inonbehaviour = index of non-behavioural runs
-% aname = label of the analysis
-% iplot = 1 plot cumulative distributions (default)
-% iplot = 0 no plots
-% ipar = index array of parameters to plot
-% dirname = (OPTIONAL) path of the directory where to save
+function [proba, dproba] = stability_mapping_univariate(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle)
+% [proba, dproba] = stability_mapping_univariate(lpmat, ibehaviour, inonbehaviour, aname, fname_, options_, parnames, estim_params_, iplot, ipar, dirname, pcrit, atitle)
+% Inputs:
+% - lpmat               [double]        Monte Carlo matrix
+% - ibehaviour          [integer]       index of behavioural runs
+% - inonbehaviour       [integer]       index of non-behavioural runs
+% - aname               [string]        label of the analysis
+% - fname_              [string]        file name
+% - options_            [structure]     options structure
+% - parnames            [char]          parameter name vector
+% - estim_params_       [structure]     characterizing parameters to be estimated
+% - iplot               [boolean]       1 plot cumulative distributions (default)
+%                                       0 no plots
+% - ipar                [integer]       index array of parameters to plot
+% - dirname             [string]        (OPTIONAL) path of the directory where to save
 %                                       (default: current directory)
-% pcrit = (OPTIONAL) critical value of the pvalue below which show the plots
+% - pcrit               [double]        (OPTIONAL) critical value of the pvalue below which show the plots
 %
 % Plots: dotted lines for BEHAVIOURAL
 %        solid lines for NON BEHAVIOURAL
-% USES smirnov
+% USES gsa.smirnov_test.m
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
@@ -38,16 +42,13 @@ function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, i
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global estim_params_ bayestopt_ M_ options_
-
-if nargin<5
+if nargin<9
     iplot=1;
 end
-fname_ = M_.fname;
-if nargin<7
+if nargin<11
     dirname='';
 end
-if nargin<9,
+if nargin<13
     atitle=aname;
 end
 
@@ -59,49 +60,49 @@ nshock = nshock + estim_params_.ncn;
 npar=size(lpmat,2);
 ishock= npar>estim_params_.np;
 
-if nargin<6
+if nargin<10
     ipar=[];
 end
-if nargin<8 || isempty(pcrit)
+if nargin<12 || isempty(pcrit)
     pcrit=1;
 end
 
-% Smirnov test for Blanchard;
+% Smirnov test for Blanchard
+proba=NaN(npar,1);
+dproba=NaN(npar,1);
 for j=1:npar
-    [H,P,KSSTAT] = smirnov(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
+    [~,P,KSSTAT] = gsa.smirnov_test(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
     proba(j)=P;
     dproba(j)=KSSTAT;
 end
 if isempty(ipar)
-    %     ipar=find(dproba>dcrit);
     ipar=find(proba<pcrit);
 end
 nparplot=length(ipar);
 if iplot && ~options_.nograph
     lpmat=lpmat(:,ipar);
-    ftit=bayestopt_.name(ipar+nshock*(1-ishock));
+    ftit=parnames(ipar+nshock*(1-ishock));
 
     for i=1:ceil(nparplot/12)
         hh_fig=dyn_figure(options_.nodisplay,'name',atitle);
         for j=1+12*(i-1):min(nparplot,12*i)
             subplot(3,4,j-12*(i-1))
             if ~isempty(ibehaviour)
-                h=cumplot(lpmat(ibehaviour,j));
+                h=gsa.cumplot(lpmat(ibehaviour,j));
                 set(h,'color',[0 0 1], 'linestyle',':','LineWidth',1.5)
             end
             hold on
             if ~isempty(inonbehaviour)
-                h=cumplot(lpmat(inonbehaviour,j));
+                h=gsa.cumplot(lpmat(inonbehaviour,j));
                 set(h,'color',[0 0 0],'LineWidth',1.5)
             end
-            %     title([ftit{j},'. D-stat ', num2str(dproba(ipar(j)),2)],'interpreter','none')
             title([ftit{j},'. p-value ', num2str(proba(ipar(j)),2)],'interpreter','none')
         end
         if nparplot>12
             dyn_saveas(hh_fig,[dirname,filesep,fname_,'_',aname,'_SA_',int2str(i)],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([dirname,filesep,fname_,'_',aname,'_SA_',int2str(i) '.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_1.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by gsa.stability_mapping_univariate.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');
@@ -116,7 +117,7 @@ if iplot && ~options_.nograph
             dyn_saveas(hh_fig,[dirname,filesep,fname_,'_',aname,'_SA'],options_.nodisplay,options_.graph_format);
             if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                 fidTeX = fopen([dirname,filesep,fname_,'_',aname,'_SA.tex'],'w');
-                fprintf(fidTeX,'%% TeX eps-loader file generated by stab_map_1.m (Dynare).\n');
+                fprintf(fidTeX,'%% TeX eps-loader file generated by gsa.stability_mapping_univariate.m (Dynare).\n');
                 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
                 fprintf(fidTeX,'\\begin{figure}[H]\n');
                 fprintf(fidTeX,'\\centering \n');

matlab/gsa/stand_.m → matlab/+gsa/standardize_columns.m

-function [y, meany, stdy] = stand_(x)
-% STAND_  Standardise a matrix by columns
+function [y, meany, stdy] = standardize_columns(x)
+% [y, meany, stdy] = standardize_columns(x)
+% Standardise a matrix by columns
 %
 % [x,my,sy]=stand(y)
 %
-% y: Time series (column matrix)
+% Inputs:
+%  - x: Time series (column matrix)
 %
-% x: standardised equivalent of y
-% my: Vector of mean values for each column of y
-% sy: Vector of standard deviations for each column of y
+%  - y:         standardised equivalent of x
+%  - meany:     Vector of mean values for each column of x
+%  - stdy:      Vector of standard deviations for each column of x
 %
 % Written by Marco Ratto
 % Joint Research Centre, The European Commission,
 % marco.ratto@ec.europa.eu
 
 % Copyright © 2012 European Commission
-% Copyright © 2012-2017 Dynare Team%
+% Copyright © 2012-2023 Dynare Team
 % This file is part of Dynare.
 %
 % Dynare is free software: you can redistribute it and/or modify
@@ -34,9 +36,11 @@ if nargin==0
     return
 end
 
+meany=NaN(size(x,2),1);
+stdy=NaN(size(x,2),1);
+y=NaN(size(x));
 for j=1:size(x,2)
-    meany(j)=mean(x(find(~isnan(x(:,j))),j));
-    stdy(j)=std(x(find(~isnan(x(:,j))),j));
+    meany(j)=mean(x(~isnan(x(:,j)),j));
+    stdy(j)=std(x(~isnan(x(:,j)),j));
     y(:,j)=(x(:,j)-meany(j))./stdy(j);
 end
\ No newline at end of file
-% end of m-file
\ No newline at end of file

matlab/gsa/teff.m → matlab/+gsa/teff.m

@@ -37,15 +37,12 @@ if ndim==3
     [ir, ic]=(find( (tmax-tmin)>1.e-8));
     j0 = length(ir);
     yt=zeros(Nsam, j0);
-
     for j=1:j0
         y0=squeeze(T(ir(j),ic(j),:));
-        %y1=ones(size(lpmat,1),1)*NaN;
         y1=ones(Nsam,1)*NaN;
         y1(istable,1)=y0;
         yt(:,j)=y1;
     end
-
 else
     tmax=max(T,[],2);
     tmin=min(T,[],2);
@@ -53,7 +50,4 @@ else
     j0 = length(ir);
     yt=NaN(Nsam, j0);
     yt(istable,:)=T(ir,:)';
-
-
 end
-%clear y0 y1;

matlab/gsa/th_moments.m → matlab/+gsa/th_moments.m

-function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
-% [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
+function [vdec, corr, autocorr, z, zz] = th_moments(dr,options_,M_)
+% [vdec, corr, autocorr, z, zz] = th_moments(dr,options_,M_)
+% Computes theoretical moments for GSA
+%
+% INPUTS
+% - dr            [structure]     model information structure
+% - options_      [structure]     Matlab's structure describing the current options
+% - M_            [structure]     Matlab's structure describing the model
+%
+% OUTPUTS
+% - vdec          [double]        variance decomposition matrix
+% - corr          [double]        correlation matrix
+% - autocorr      [cell]          contains autocorrelation or
+%                                 auto- and cross-covariance matrices
+% - z             [double]        matrix containing mean, standard
+%                                 deviation, and variance vector
+% - zz            [double]        autocorrelation matrix
 
-% Copyright © 2012-2018 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -18,18 +33,16 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global M_ oo_ options_
-
-nvar = length(var_list);
+nvar = length(options_.varobs);
 if nvar == 0
     nvar = length(dr.order_var);
     ivar = [1:nvar]';
 else
     ivar=zeros(nvar,1);
     for i=1:nvar
-        i_tmp = strmatch(var_list{i}, M_.endo_names, 'exact');
+        i_tmp = strmatch(options_.varobs{i}, M_.endo_names, 'exact');
         if isempty(i_tmp)
-            error(['One of the variables specified does not exist']) ;
+            error('th_moments: One of the variables specified does not exist');
         else
             ivar(i) = i_tmp;
         end
@@ -39,21 +52,11 @@ end
 [gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
 m = dr.ys(ivar(stationary_vars));
 
-
-%   i1 = find(abs(diag(gamma_y{1})) > 1e-12);
-i1 = [1:length(ivar)];
+i1 = 1:length(ivar);
 s2 = diag(gamma_y{1});
 sd = sqrt(s2);
 
-
 z = [ m sd s2 ];
-mean = m;
-var = gamma_y{1};
-
-
-%'THEORETICAL MOMENTS';
-%'MEAN','STD. DEV.','VARIANCE');
-z;
 
 %'VARIANCE DECOMPOSITION (in percent)';
 if M_.exo_nbr>1
@@ -69,6 +72,7 @@ else
     corr = gamma_y{1}(i1,i1);
 end
 if options_.ar > 0
+    zz=NaN(length(ivar),options_.ar);
     %'COEFFICIENTS OF AUTOCORRELATION';
     for i=1:options_.ar
         if options_.opt_gsa.useautocorr

matlab/+gui/+perfect_foresight/run.m

@@ -13,7 +13,7 @@ function run(json)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2019-2022 Dynare Team
+% Copyright © 2019-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,7 +30,7 @@ function run(json)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
 
-global M_ options_ oo_ ys0_ ex0_
+global M_ options_ oo_
 
 %% Check Inputs
 if nargin ~= 1 || ~ischar(json)
@@ -64,12 +64,12 @@ end
 
 %% ENDVAL instructions
 % initialize exogenous shocks to zero and compute final ss unless there is a permanent shock
-ys0_ = [];
-ex0_ = [];
+oo_.initial_steady_state = [];
+oo_.initial_exo_steady_state = [];
 M_.det_shocks = [];
 if ~isempty(jm.anticipated_permanent_shocks) || ~isempty(jm.endval_endo)
-    ys0_= oo_.steady_state;
-    ex0_ = oo_.exo_steady_state;
+    oo_.initial_steady_state= oo_.steady_state;
+    oo_.initial_exo_steady_state = oo_.exo_steady_state;
     for i = 1:length(jm.endval_endo)
         oo_.steady_state(jm.endval_endo(i).id) = jm.endval_endo(i).value;
     end
@@ -84,7 +84,7 @@ if ~isempty(jm.anticipated_permanent_shocks) || ~isempty(jm.endval_endo)
                 struct(...
                 'exo_det', false, ...
                 'exo_id', s.exo_id, ...
-                'multiplicative', false, ...
+                'type', 'level', ...
                 'periods', 1:s.start_date, ...
                 'value', 0)];
         end
@@ -99,11 +99,10 @@ if ~isempty(jm.anticipated_transitory_shocks)
             M_.det_shocks; ...
             struct('exo_det', false, ...
             'exo_id', s.exo_id, ...
-            'multiplicative', false, ...
+            'type', 'level', ...
             'periods', s.start_date:s.end_date, ...
             'value', s.value)];
     end
-    M_.exo_det_length = 0;
 end
 
 %% Make unanticipated shock map
@@ -142,26 +141,26 @@ mapkeys = unique(cell2mat([keys(unanticipated_p_shocks) keys(unanticipated_t_sho
 
 %% Simulation
 options_.periods = jm.periods;
-perfect_foresight_setup;
+oo_=perfect_foresight_setup(M_, options_, oo_);
 
 % no surprise shocks present
 if isempty(mapkeys)
-    perfect_foresight_solver;
+    oo_=perfect_foresight_solver(M_, options_, oo_);
     return
 end
 
 % surprise shocks present
 % in case there are unanticipated shocks...
-if isempty(ys0_)
+if isempty(oo_.initial_steady_state)
     yy = oo_.steady_state;
 else
-    yy = ys0_;
+    yy = oo_.initial_steady_state;
 end
 
 if mapkeys(1) ~= 1
     % if first unanticipated shock is not in period 1
     % simulate until first unanticipated shock and save
-    perfect_foresight_solver;
+    oo_=perfect_foresight_solver(M_, options_, oo_);
     yy = [yy oo_.endo_simul(:, 2:mapkeys(1)+1)];
 end
 
@@ -203,7 +202,7 @@ for i = 1:length(mapkeys)
     last_period = this_period;
     assert(rows(oo_.exo_simul) == oo_exo_simul_rows, 'error encountered setting oo_.exo_simul');
     oo_.endo_simul(:, 1) = yy(:, end);
-    perfect_foresight_solver;
+    oo_=perfect_foresight_solver(M_, options_, oo_);
     if next_period > 0
         yy = [yy oo_.endo_simul(:, 2:next_period-this_period+1)];
     else

matlab/+gui/+stochastic-simulation/read.m

@@ -13,7 +13,7 @@ function read(json)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright © 2019-2020 Dynare Team
+% Copyright © 2019-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,7 +36,6 @@ global M_ options_ oo_
 jm = loadjson_(json, 'SimplifyCell', 1);
 data2json=struct();
 
-M_.exo_det_length = 0;
 for nshocks = 1:length(jm.stochasticshocksdescription)
 	covartype=jm.stochasticshocksdescription{nshocks}.shockattributevalue;
 	thisshock=(jm.stochasticshocksdescription{nshocks}.shockindex)+1;

matlab/+hank/check_steady_state_input.m

+% Copyright © 2025 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/>.
+%
+% Compute the stochastic simulations of heterogeneus-agent models
+%
+% INPUTS
+% - M_       [structure] Matlab's structure describing the model
+% - options_ [structure] Matlab's structure describing the current options
+% - ss       [structure] Matlab's structure with the model steady-state
+%                       information. It contains the following fields:
+%    - pol [structure] Matlab's structure containing the policy functions
+%                      discretization
+%       - pol.grids [structure]: Matlab's structure containing the nodes of the
+%                                state grids as column vectors.
+%       - pol.values [structure]: Matlab's structure containing the policy
+%                                 function as matrices. Row indices and column
+%                                 indices follow the lexicographic order
+%                                 specified in pol.shocks and pol.states
+%       - pol.shocks [array]: a list containing the order of shock/innovation
+%                             variables used for the rows of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement (resp. varexo(heterogeneity=)
+%                             statement) for discretizes AR(1) processes (resp.
+%                             discretized gaussian i.i.d innovations).
+%       - pol.states [array]: a list containing the order of state variables
+%                             used for the columns of pol.values. If not
+%                             specified, it follows the declaration order in the
+%                             var(heterogeneity=) statement.
+%    - shocks [structure]: Matlab's stucture describing the discretization of
+%                          individual shocks or innovations:
+%       - shocks.grids [structure]: Matlab's structure containing the nodes of
+%                                   the shock grids
+%       - shocks.Pi [structure]: Matlab's structure containing the Markov
+%                                matrices if the shock processes are discretized
+%                                AR(1) processes. The field should be absent
+%                                otherwise.
+%       - shocks.w [structure]: Matlab's structure containing the Gauss-Hermite
+%                               weights if the i.i.d gaussian innovation
+%                               processes are discretized
+%    - d [structure]: Matlab's structure describing the steady-state
+%                     distribution
+%       - d.grids [structure]: structure containing the states grids as column
+%                              vectors. If one of the states grid is not
+%                              specified, it is assumed to be identical to
+%                              its policy-function value stored in pol.grids
+%       - d.hist [matrix]: the histogram of the distribution as a matrix with
+%                          shocks/innovation indices as rows and state indices
+%                          as columns. Row and column indices follow the
+%                          lexicographic order specified in d.shocks and
+%                          d.states
+%       - d.shocks [array]: a list containing the order of shock/innovation
+%                           variables used for the rows of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.shocks.
+%       - d.states [array]: a list containing the order of state variables used
+%                           for the columns of d.hist. If it is not
+%                           specified, it falls back to the value induced by
+%                           pol.states.
+%    - agg [structure]: Matlab's tructure containing the steady-state values of
+%                       aggregate variables
+% OUTPUTS
+% - out_ss [structure]: validated and normalized steady-state input structure.
+%                       It has a similar structure as the ss structure in inputs.
+% - sizes [structure]: structure containing sizes information
+%    - n_e [integer]: number of shocks
+%    - N_e [integer]: total number of nodes in the shocks grids
+%    - shocks [structure]: number of nodes for each shock grid
+%    - n_a [integer]: number of states
+%    - pol [structure]: 
+%       - pol.N_a [integer]: total number of nodes in the policy state grids
+%       - pol.states [structure]: number of nodes for each state grid
+%    - d [structure]:
+%       - d.states [structure]: number of nodes for each distribution state grid
+%       - d.N_a [integer]: total number of nodes in the distribution grids
+%    - n_pol [integer]: number of policy variables
+%    - agg [integer]: number of aggregate variables
+function [out_ss, sizes] = check_steady_state_input(M_, options_, ss)
+   %% Checks
+   % Retrieve variables information from M_
+   H_ = M_.heterogeneity(1);
+   state_symbs = H_.endo_names(H_.state_var);
+   inn_symbs = H_.exo_names;
+   agg_symbs = M_.endo_names(1:M_.orig_endo_nbr);
+   required_pol_symbs = H_.endo_names(1:H_.orig_endo_nbr);
+   mult_symbs = {};
+   for i = 1:numel(H_.aux_vars)
+      if (H_.aux_vars(i).type == 15)
+         mult_symbs{end+1} = H_.endo_names{H_.aux_vars(i).endo_index};
+      end
+   end
+   % Initialize output variables
+   sizes = struct;
+   out_ss = struct;
+   if ~isstruct(ss)
+      error('Misspecified steady-state input `ss`: the steady-state input `ss` is not a structure.')
+   end
+   %% Shocks
+   % Check that the field `ss.shocks` exists
+   check_isfield('shocks', ss, 'ss.shocks', ' Models without idiosyncratic shocks but with ex-post heterogeneity over individual state variables cannot be solved yet.')
+   shocks = ss.shocks;
+   % Check that `ss.shocks` is a structure
+   check_isstruct(shocks, 'ss.shocks');
+   % Check that the field `ss.shocks.grids` exists
+   check_isfield('grids', ss.shocks, 'ss.shocks.grids');
+   % Check that `ss.shocks.grids` is a structure
+   check_isstruct(shocks.grids, 'ss.shocks.grids');
+   shock_symbs = fieldnames(shocks.grids);
+   % Check the types of `ss.shocks.grids` values
+   check_fun(shocks.grids, 'ss.shocks.grids', shock_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(shock_symbs)
+      s = shock_symbs{i};
+      if isrow(ss.shocks.grids.(s))
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s)';
+      else
+         out_ss.shocks.grids.(s) = ss.shocks.grids.(s);
+      end
+   end
+   % Check shock discretization
+   flag_ar1 = isfield(shocks, 'Pi');
+   flag_gh = isfield(shocks, 'w');
+   if ~flag_ar1 && ~flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both non-specified, while exactly one of them should be set.');
+   end
+   % Check that either individual AR(1) processes or i.i.d gaussion innovation
+   % processes are discretized
+   if flag_ar1 && flag_gh
+      error('Misspecified steady-state input `ss`: the `ss.shocks.Pi` and `ss.shocks.w` fields are both specified, while only one of them should be.');
+   end
+   % Case of discretized AR(1) exogenous processes
+   if flag_ar1
+      % Check that the grids specification for AR(1) shocks and the
+      % var(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, required_pol_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).endo_names');
+      % Check that shocks.Pi is a structure
+      check_isstruct(shocks.Pi, 'ss.shocks.Pi');
+      % Check that the grids specification for individual AR(1) shocks, the
+      % information in the M_ structure and the Markov transition matrices
+      % specification are compatible and display a warning if redundancies are
+      % detected
+      check_missingredundant(shocks.Pi, 'ss.shocks.Pi', shock_symbs, options_.hank.nowarningredundant);
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the type of shocks.Pi field values
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) ismatrix(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real matrices');
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the size of their counterparts in `ss.shocks.grids`', grid_nb);
+      % Check that the matrices provided in shocks.Pi are Markov matrices:
+      %     - all elements are non-negative
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(x >= 0, 'all'), 'contain negative elements')
+      %     - the sum of each row equals 1.
+      check_fun(shocks.Pi, 'ss.shocks.Pi', shock_symbs, @(x) all(abs(sum(x,2)-1) < options_.hank.tol_check_sum), 'have row sums different from 1.')
+      % Remove AR(1) shock processes from state and policy variables
+      state_symbs = setdiff(state_symbs, shock_symbs);
+      required_pol_symbs = setdiff(required_pol_symbs, shock_symbs);
+      % Copy relevant shock-related elements in out_ss 
+      out_ss.shocks.Pi = ss.shocks.Pi;
+   end
+   relevant_pol_symbs = [required_pol_symbs; mult_symbs];
+   % Case of discretized i.i.d gaussian innovations
+   if flag_gh
+      % Check that shocks.w is a structure
+      check_isstruct(shocks.w, 'shocks.w');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the varexo(heterogeneity=) statement are compatible
+      check_consistency(shock_symbs, inn_symbs, 'ss.shocks.grids', 'M_.heterogeneity(1).exo_names');
+      % Verify that the grids specification for individual i.i.d innovations and
+      % the Gauss-Hermite weights specification are compatible and display a
+      % warning if redundancies are detected
+      check_missingredundant(shocks.w, 'ss.shocks.w', shock_symbs, options_.hank.nowarningredundant);
+      % Check the type of `ss.shocks.w` elements
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) isvector(x) && isnumeric(x) && isreal(x) && ~issparse(x), 'are not dense real vectors');
+      % Store:
+      %    - the number of shocks
+      sizes.n_e = numel(shock_symbs);
+      %    - the size of the tensor product of the shock grids 
+      grid_nb = cellfun(@(s) numel(shocks.grids.(s)), shock_symbs); 
+      sizes.N_e = prod(grid_nb);
+      %    - the size of each shock grid
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         sizes.shocks.(s) = numel(shocks.grids.(s));
+      end
+      % Check the internal size compatibility of discretized shocks processes
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) numel(x), 'have incompatible sizes with those of `ss.shocks.grids`', grid_nb);
+      % Check that the arrays provided in shocks.w have the properties of
+      % Gauss-Hermite weights:
+      %     - all elements are non-negative
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(x >= 0.), 'contain negative elements');
+      %     - the sum of Gauss-Hermite weights is 1
+      check_fun(shocks.w, 'ss.shocks.w', shock_symbs, @(x) all(abs(sum(x)-1) < options_.hank.tol_check_sum), 'have sums different from 1.');
+      % Make `ss.shocks.w` values column vectors in the output steady-state
+      % structure
+      for i=1:numel(shock_symbs)
+         s = shock_symbs{i};
+         if isrow(ss.shocks.w.(s))
+            out_ss.shocks.w.(s) = ss.shocks.w.(s)';
+         else
+            out_ss.shocks.w.(s) = ss.shocks.w.(s);
+         end
+      end
+   end
+   %% Policy functions
+   % Check that the `ss.pol` field exists
+   check_isfield('pol', ss, 'ss.pol');
+   pol = ss.pol;
+   % Check that `ss.pol` is a structure
+   check_isstruct(ss.pol, 'ss.pol');
+   % Check that the `ss.pol`.grids field exists
+   check_isfield('grids', ss.pol, 'ss.pol.grids');
+   % Check that `ss.pol.grids` is a structure
+   check_isstruct(ss.pol.grids, 'ss.pol.grids');
+   % Check that the state grids specification and the var(heterogeneity=)
+   % statement are compatible
+   check_missingredundant(pol.grids, 'ss.pol.grids', state_symbs, options_.hank.nowarningredundant);
+   % Check that `ss.pol.grids` values are dense real vectors
+   check_fun(pol.grids, 'ss.pol.grids', state_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+   % Store:
+   %    - the number of states
+   sizes.n_a = numel(state_symbs);
+   %    - the size of the tensor product of the states grids 
+   grid_nb = cellfun(@(s) numel(pol.grids.(s)), state_symbs);
+   sizes.pol.N_a = prod(grid_nb);
+   %    - the size of each state grid
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      sizes.pol.states.(s) = numel(pol.grids.(s));
+   end
+   % Make `ss.shocks.grids` values column vector in the output steady-state
+   % structure
+   for i=1:numel(state_symbs)
+      s = state_symbs{i};
+      if isrow(ss.pol.grids.(s))
+         out_ss.pol.grids.(s) = ss.pol.grids.(s)';
+      else
+         out_ss.pol.grids.(s) = ss.pol.grids.(s);
+      end
+   end
+   % Check that the field `ss.pol.values` exists 
+   check_isfield('values', pol, 'ss.pol.values');
+   % Check that `ss.pol.values` is a struct
+   check_isstruct(pol.values, 'ss.pol.values');
+   % Check the missing and redundant variables in `ss.pol.values`
+   check_missingredundant(pol.values, 'ss.pol.values', required_pol_symbs, true);
+   pol_values = fieldnames(pol.values);
+   pol_values_in_relevant_pol_symbs = ismember(pol_values, relevant_pol_symbs);
+   if ~options_.hank.nowarningredundant
+      if ~all(pol_values_in_relevant_pol_symbs)
+         warning('Steady-state input `ss`. The following fields are redundant in `%s`: %s.', 'ss.pol.values', strjoin(pol_values(~pol_values_in_relevant_pol_symbs)));
+      end
+   end
+   provided_relevant_pol_symbs = pol_values(pol_values_in_relevant_pol_symbs);
+   % Check that `ss.pol.values` values are dense real matrices
+   check_fun(pol.values, 'ss.pol.values', provided_relevant_pol_symbs, @(x) isnumeric(x) && isreal(x) && ismatrix(x) && ~issparse(x), 'are not dense real matrices');
+   % Check the internal size compatibility of `ss.pol.values`
+   sizes.n_pol = H_.endo_nbr;
+   check_fun(pol.values, 'ss.pol.values', provided_relevant_pol_symbs, @(x) size(x,1), 'have a number of rows that is not consistent with the sizes of `ss.shocks.grids` elements', sizes.N_e);
+   check_fun(pol.values, 'ss.pol.values', provided_relevant_pol_symbs, @(x) size(x,2), 'have a number of columns that is not consistent with the sizes of `ss.pol.grids` elements', sizes.pol.N_a);
+   % Copy `ss.pol.values` in `out_ss`
+   out_ss.pol.values = ss.pol.values;
+   % Check the permutation of state variables for policy functions
+   out_ss.pol.states = check_permutation(pol, 'states', 'ss.pol', state_symbs);
+   % Check the permutation of shock variables for policy functions
+   out_ss.pol.shocks = check_permutation(pol, 'shocks', 'ss.pol', shock_symbs);
+   %% Distribution
+   % Check that the field `ss.d` exists 
+   check_isfield('d', ss, 'ss.d');
+   d = ss.d;
+   % Check that the field `ss.d.hist` exists 
+   check_isfield('hist', ss.d, 'ss.d.hist');
+   % Check the type of `ss.d.hist`
+   if ~(ismatrix(d.hist) && isnumeric(d.hist) && isreal(d.hist) && ~issparse(d.hist))
+      error('Misspecified steady-state input `ss`: `ss.d.hist` is not a dense real matrix.');
+   end
+   % Check the consistency of `ss.d.grids` 
+   if ~isfield(d, 'grids')
+      if ~options_.hank.nowarningdgrids
+         warning('In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states %s. The policy grids in `ss.pol.grids` shall be used.' , strjoin(state_symbs));
+      end
+      % Copy the relevant sizes from the pol field
+      sizes.d = sizes.pol;
+      out_ss.d.grids = out_ss.pol.grids;
+   else
+      % Check that `ss.d.grids` is a struct
+      check_isstruct(d.grids, 'ss.d.grids');
+      % Check redundant variables in `ss.d.grids`
+      d_grids_symbs = fieldnames(d.grids);
+      if isempty(d_grids_symbs)
+         if ~options_.hank.nowarningredundant
+            warning('In the steady-state input `ss.d.grids`, no distribution-specific grid is set for states %s. The policy grids in `ss.pol.grids` shall be used.' , strjoin(state_symbs));
+         end
+         % Copy the relevant sizes from the pol field
+         sizes.d = sizes.pol;
+         out_ss.d.grids = out_ss.pol.grids;
+      else
+         d_grids_symbs_in_states = ismember(d_grids_symbs, state_symbs);
+         if ~all(d_grids_symbs_in_states)
+            if ~options_.hank.nowarningredundant
+               warning('In the steady-state input `ss.d.states`, the following specification for the states grids in the distribution structure are not useful: %s.', strjoin(d_grids_symbs(~d_grids_symbs_in_states)));
+            end
+            d_grids_symbs = d_grids_symbs(d_grids_symbs_in_states);
+         end
+         % Check the types of `ss.d.grids` elements
+         check_fun(pol.grids, 'ss.d.grids', d_grids_symbs, @(x) isnumeric(x) && isreal(x) && isvector(x) && ~issparse(x), 'are not dense real vectors');
+         % Store the size of the states grids for the distribution
+         for i=1:numel(d_grids_symbs)
+            s = d_grids_symbs{i};
+            sizes.d.states.(s) = numel(d.grids.(s));
+            out_ss.d.grids.(s) = d.grids.(s);
+         end
+         states_out_of_d = setdiff(state_symbs, d_grids_symbs);
+         if ~isempty(states_out_of_d)
+            if ~options_.hank.nowarningdgrids
+               warning('hank.bbeg.het_stoch_simul: in the steady-state structure, no distribution-specific grid set for states %s. The policy grids specified in pol.grids shall be used.', strjoin(states_out_of_d));
+            end
+            % Store the sizes of the unspecified states grids from the pol field
+            for i=1:numel(states_out_of_d)
+               s = states_out_of_d{i};
+               sizes.d.states.(s) = sizes.pol.states.(s);
+               out_ss.d.grids.(s) = pol.grids.(s);
+            end
+         end
+      end
+   end
+   % Check the internal size compatibility of the distribution histogram
+   if size(d.hist,1) ~= sizes.N_e
+      error('Misspecified steady-state input `ss`: the number of rows of the histogram matrix `ss.d.hist` is not consistent with the sizes of shocks grids `ss.shocks.grids`.');
+   end
+   sizes.d.N_a = prod(structfun(@(x) x, sizes.d.states));
+   if size(d.hist,2) ~= sizes.d.N_a
+      error('Misspecified steady-state input `ss`: the number of columns of the histogram matrix `ss.d.hist` is not consistent with the sizes of states grids `ss.d.grids`/`ss.pol.grids`.');
+   end
+   % Copy `ss.d.hist` in `out_ss` 
+   out_ss.d.hist = ss.d.hist;
+   % Check the permutation of state variables in the distribution
+   out_ss.d.states = check_permutation(d, 'states', 'ss.d', state_symbs);
+   % Check the permutation of shock variables in the distribution
+   out_ss.d.shocks = check_permutation(d, 'shocks', 'ss.d', shock_symbs);
+   %% Aggregate variables
+   % Check that the field `ss.agg` exists
+   check_isfield('agg', ss, 'ss.agg');
+   agg = ss.agg;
+   % Check that `ss.agg` is a structure
+   check_isstruct(agg, 'ss.agg');
+   % Check that the aggregate variables specification and the var statement are
+   % compatible
+   check_missingredundant(agg, 'ss.agg', agg_symbs, options_.hank.nowarningredundant);
+   % Store the number of aggregate variables
+   sizes.agg = numel(fieldnames(agg));
+   % Check the types of `ss.agg` values
+   check_fun(agg, 'ss.agg', agg_symbs, @(x) isreal(x) && isscalar(x), 'are not real scalars');
+   % Copy `ss.agg` into `out_ss`
+   out_ss.agg = agg;
+end
\ No newline at end of file

matlab/+hank/compute_steady_state_residuals.m

+% Copyright © 2025 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/>.
+%
+% COMPUTE_STEADY_STATE_RESIDUALS Computes residuals for a heterogeneous-agent model at the steady state
+%
+% This function evaluates the residuals of both aggregate and individual (heterogeneous-agent) 
+% equilibrium conditions given a steady-state object stored in `oo_.hank`. It uses previously 
+% computed policy functions and basis matrices to evaluate the model’s residuals pointwise.
+%
+% INPUTS
+%   M_   [struct] : Dynare model structure
+%   oo_  [struct] : Dynare output structure. Must contain the fields:
+%                    - oo_.hank.ss     : steady-state structure (see initialize_steady_state)
+%                    - oo_.hank.sizes  : tensor-product grid sizes
+%                    - oo_.hank.mat    : basis and interpolation matrices
+%
+% OUTPUTS
+%   F [matrix] : Residuals of the heterogeneous-agent equations (size: H_.endo_nbr × N_sp),
+%                where N_sp is the number of points in the full state space.
+%   G [vector] : Residuals of the aggregate equilibrium equations (size: M_.endo_nbr × 1)
+%
+% NOTES
+% - For the aggregate block:
+%     - The individual state distribution `ss.d.hist` is reordered to match the policy grid.
+%     - Aggregate variables defined as expectations over the distribution are recomputed
+%       using `mat.pol.x_bar_dash` and the `Phi` basis matrix.
+%     - The residuals are evaluated using the sparse dynamic function `.sparse.dynamic_resid`.
+%
+% - For the heterogeneous block:
+%     - Constructs the full stacked endogenous (`yh`) and exogenous (`xh`) state vectors
+%       at times t-1, t, and t+1 using the stored policy functions and interpolation structure.
+%     - Calls `.sparse.dynamic_het1_resid` at each point in the tensor-product state space.
+%
+%           dynamic_resid, dynamic_het1_resid
+function [F,G] = compute_steady_state_residuals(M_, oo_)
+   ss = oo_.hank.ss;
+   sizes = oo_.hank.sizes;
+   mat = oo_.hank.mat;
+   % Rearrange ss.d.hist so that the state and shock orders is similar to the one
+   % used for policy functions. order_shocks(i) contains the position of
+   % ss.pol.shocks(i) in ss.d.shocks. In other words, ss.d.shocks(order_shocks)
+   % equals ss.pol.shocks. A similar logic applies to order_states.
+   order_shocks = cellfun(@(x) find(strcmp(x,ss.d.shocks)), ss.pol.shocks);
+   order_states = cellfun(@(x) find(strcmp(x,ss.d.states)), ss.pol.states);
+   dim_states = cellfun(@(s) sizes.d.states.(s), ss.d.states);
+   dim_shocks = cellfun(@(s) sizes.shocks.(s), ss.d.shocks);
+   hist = reshape(ss.d.hist, [dim_shocks dim_states]);
+   hist = permute(hist, [order_shocks sizes.n_e+order_states]);
+   % Compute aggregated heterogeneous variables
+   Ix = mat.pol.x_bar_dash*mat.Phi*reshape(hist,[],1);
+   % Computing aggregate residuals
+   % - Extended aggregate endogenous variable vector at the steady state - %
+   y = NaN(3*M_.endo_nbr,1);
+   for i=1:M_.orig_endo_nbr
+      var = M_.endo_names{i};
+      y(i) = ss.agg.(var); 
+      y(M_.endo_nbr+i) = ss.agg.(var); 
+      y(2*M_.endo_nbr+i) = ss.agg.(var); 
+   end
+   % - Variables resulting from an aggregation operator - %
+   iIx = zeros(size(M_.heterogeneity_aggregates,1),1);
+   ix = zeros(size(M_.heterogeneity_aggregates,1),1);
+   for i=1:size(M_.heterogeneity_aggregates,1)
+      if (M_.heterogeneity_aggregates{i,1} == 'sum')
+         iIx(i) = M_.aux_vars(M_.heterogeneity_aggregates{i,2}).endo_index;
+         ix(i) = M_.heterogeneity_aggregates{i,3};
+      end
+   end
+   yagg = Ix(ix);
+   y(iIx) = yagg;
+   y(M_.endo_nbr+iIx) = yagg;
+   y(2*M_.endo_nbr+iIx) = yagg;
+   % - Exogenous variables - %
+   x = zeros(M_.exo_nbr,1);
+   % - Call to the residual function - %
+   G = feval([M_.fname '.sparse.dynamic_resid'], y, x, M_.params, [], yagg);
+   % Compute heterogeneous residuals
+   N_sp = sizes.N_e*sizes.pol.N_a;
+   H_ = M_.heterogeneity(1);
+   % Heterogeneous endogenous variable vector
+   yh = NaN(3*H_.endo_nbr, N_sp);
+   % Heterogeneous exogenous variable vector
+   xh = NaN(H_.exo_nbr, N_sp);
+   % Get the H_.endo_names indices for the steady-state ordering ss.pol.shocks and ss.pol.states 
+   if isfield(ss.shocks, 'Pi')
+      order_shocks = cellfun(@(x) find(strcmp(x,H_.endo_names)), ss.pol.shocks);
+   else
+      order_shocks = cellfun(@(x) find(strcmp(x,H_.exo_names)), ss.pol.shocks);
+   end
+   % Get the H_.endo_names indices for the steady-state ordering ss.pol.shocks and ss.pol.states 
+   order_states = cellfun(@(x) find(strcmp(x,H_.endo_names)), ss.pol.states);
+   % Get the indices of originally declared pol values in H_.endo_names
+   pol_names = H_.endo_names;
+   pol_names(ismember(H_.endo_names, ss.pol.shocks)) = [];
+   pol_ind = cellfun(@(x) find(strcmp(x,H_.endo_names)), pol_names);
+   % t-1
+   yh(order_states,:) = mat.pol.sm(1:sizes.n_a,:);
+   % t
+   yh(H_.endo_nbr+pol_ind,:) = mat.pol.x_bar;
+   if isfield(ss.shocks, 'Pi')
+      yh(H_.endo_nbr+order_shocks,:) = mat.pol.sm(sizes.n_a+1:end,:);
+   else
+      xh(order_shocks,:) = mat.pol.sm(sizes.n_a+1:end,:);
+   end
+   % t+1
+   yh(2*H_.endo_nbr+pol_ind,:) = mat.pol.x_bar_dash*mat.Phi_tilde_e;
+   % Call the residual function
+   F = NaN(H_.endo_nbr, N_sp);
+   for j=1:N_sp
+      F(:,j) = feval([M_.fname '.sparse.dynamic_het1_resid'], y, x, M_.params, [], yh(:,j), xh(:,j), []);
+   end
+end
\ No newline at end of file

matlab/+hank/initialize_steady_state.m

+% Copyright © 2025 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/>.
+%
+% INITIALIZE_STEADY_STATE Wrapper to validate, complete, and store steady-state input for HANK models.
+%
+% This function performs validation of the user-provided steady-state structure `ss`, constructs
+% interpolation objects and basis matrices, computes missing policy functions for complementarity 
+% multipliers (if not supplied), and populates the global `oo_.hank` structure with the results.
+%
+% INPUTS
+%   M_       [struct] : Dynare model structure
+%   options_ [struct] : Dynare options structure
+%   oo_      [struct] : Dynare output structure to which results will be written
+%   ss       [struct] : User-supplied steady-state structure (partial or complete)
+% OUTPUTS
+%   oo_      [struct] : Updated Dynare output structure containing:
+%                        - oo_.hank.ss    : validated and completed steady-state structure
+%                        - oo_.hank.sizes : grid size structure
+%                        - oo_.hank.mat   : interpolation and policy function matrices
+function oo_ = initialize_steady_state(M_, options_, oo_, ss)
+   assert(isscalar(M_.heterogeneity), 'Heterogeneous-agent models with more than one heterogeneity dimension are not allowed yet!');
+   % Check steady-state input
+   [out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
+   % Build useful basis and interpolation matrices
+   H_ = M_.heterogeneity(1);
+   mat = build_grids_and_interpolation_matrices(out_ss, sizes);
+   % Compute policy matrices without the complementarity conditions multipliers
+   [mat.pol.x_bar, mat.pol.x_bar_dash] = compute_pol_matrices(H_, out_ss, sizes, mat);
+   % Compute the policy values of complementarity conditions multipliers
+   mult_values = compute_pol_mcp_mul(M_, out_ss, sizes, mat);
+   % Update ss.pol.values accordingly
+   f = fieldnames(mult_values);
+   for i=1:numel(f)
+     out_ss.pol.values.(f{i}) = mult_values.(f{i});
+   end
+   % Update the policy matrices
+   [mat.pol.x_bar, mat.pol.x_bar_dash] = compute_pol_matrices(H_, out_ss, sizes, mat);
+   % Store the results in oo_
+   oo_.hank.ss = out_ss;
+   oo_.hank.sizes = sizes;
+   oo_.hank.mat = mat;
+end
\ No newline at end of file

matlab/+hank/private/build_grids_and_interpolation_matrices.m

+% Copyright © 2025 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/>.
+%
+% BUILD_GRIDS_AND_INTERPOLATION_MATRICES Constructs interpolation and basis matrices
+% for evaluating policy functions and distributions in heterogeneous-agent models.
+%
+% This function builds all numerical structures needed to evaluate policy residuals,
+% compute expectations, and perform interpolation over the state space. It constructs:
+% - tensor-product basis matrices (`Phi`)
+% - expectation operators (`Phi_tilde_e`)
+% - grid point matrices (`pol.sm`, `d.sm`)
+% - LU decomposition of the interpolation basis for efficient linear solves.
+%
+% INPUTS
+%   ss    [struct] : Steady-state structure validated by `check_steady_state_input`
+%   sizes [struct] : Structure containing grid dimensions and counts
+%
+% OUTPUT
+%   mat   [struct] : Structure containing matrices used in solution and simulation:
+%       - Phi            : basis matrix for evaluating interpolated functions over distribution grid
+%       - Phi_tilde      : expanded identity matrix for state variables (used in solving)
+%       - Phi_tilde_e    : expectation operator incorporating both interpolation and transition shocks
+%       - Mu             : Kronecker product of exogenous transition or quadrature weights
+%       - pol.sm         : full tensor grid of state and shock combinations for policy evaluation
+%       - d.sm           : full tensor grid of state and shock combinations for distribution evaluation
+%       - L_Phi_tilde,
+%         U_Phi_tilde,
+%         P_Phi_tilde    : LU decomposition of `Phi_tilde` (used to project policy values)
+%
+% Notes:
+% - The structure `ss` must contain valid grids and values in `ss.pol.values`, `ss.pol.states`,
+%   `ss.pol.shocks`, `ss.d.grids`, and `ss.shocks`.
+% - Assumes linear interpolation and separable basis construction across state variables.
+% - `Phi_tilde_e` applies the expected operator `E_t[f(x')]` for the simulation step.
+function mat = build_grids_and_interpolation_matrices(ss, sizes)
+   mat = struct;
+   % Compute Φ and ̃Φ with useful basis matrices
+   Phi = speye(sizes.N_e);
+   Phi_tilde = speye(sizes.N_e);
+   B = struct;
+   B_tilde_bar = struct;
+   for i=sizes.n_a:-1:1
+      s = ss.pol.states{i};
+      [ind, w] = find_bracket_linear_weight(ss.pol.grids.(s), ss.d.grids.(s));
+      B.(s) = lin_spli(ind, w, sizes.pol.states.(s));
+      [ind, w] = find_bracket_linear_weight(ss.pol.grids.(s), reshape(ss.pol.values.(s), 1, []));
+      B_tilde_bar.(s) = lin_spli(ind, w, sizes.pol.states.(s));
+      Phi = kron(B.(s), Phi);
+      Phi_tilde = kron(speye(sizes.pol.states.(s)), Phi_tilde);
+   end
+   mat.Phi = Phi;
+   mat.Phi_tilde = Phi_tilde;
+   % Compute ̃Φᵉ and Mu
+   N_sp = sizes.pol.N_a*sizes.N_e;
+   hadamard_prod = ones(sizes.pol.N_a,N_sp);
+   prod_a_1_to_jm1 = 1;
+   prod_a_jp1_to_n_a = sizes.pol.N_a;
+   for j=1:sizes.n_a
+      s = ss.pol.states{j};
+      prod_a_jp1_to_n_a = prod_a_jp1_to_n_a/sizes.pol.states.(s);
+      hadamard_prod = hadamard_prod .* kron(kron(ones(prod_a_1_to_jm1,1), B_tilde_bar.(s)), ones(prod_a_jp1_to_n_a,1));
+      prod_a_1_to_jm1 = prod_a_1_to_jm1*sizes.pol.states.(s);
+   end
+   Mu = eye(1);
+   if isfield(ss.shocks, 'Pi')
+      for i=sizes.n_e:-1:1
+         e = ss.pol.shocks{i};
+         Mu = kron(ss.shocks.Pi.(e), Mu);
+      end
+      Phi_tilde_e = kron(hadamard_prod, ones(sizes.N_e,1)).*kron(ones(sizes.pol.N_a), Mu');
+   else
+      for i=sizes.n_e:-1:1
+         e = ss.pol.shocks{i};
+         Mu = kron(ss.shocks.w.(e), Mu);
+      end
+      Phi_tilde_e = kron(hadamard_prod, Mu);
+   end
+   mat.Mu = Mu;
+   mat.Phi_tilde_e = Phi_tilde_e;
+   % Compute state and shock grids for the policy and distribution state grids
+   mat.pol.sm = set_state_matrix(ss, sizes, "pol");
+   mat.d.sm = set_state_matrix(ss, sizes, "d");
+   % LU decomposition of Phi_tilde 
+   [L_Phi_tilde,U_Phi_tilde,P_Phi_tilde] = lu(Phi_tilde);
+   mat.L_Phi_tilde = L_Phi_tilde;
+   mat.U_Phi_tilde = U_Phi_tilde;
+   mat.P_Phi_tilde = P_Phi_tilde;
+end
\ No newline at end of file

matlab/+hank/private/check_consistency.m

+% Copyright © 2025 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/>.
+%
+% Check consistency between two sets of variable names.
+%
+% INPUTS
+% - f1 [cell array of strings]: list of variable names to be checked for consistency
+% - f2 [cell array of strings]: reference list of valid variable names
+% - f1_name [char]: name of the first set (for error messages)
+% - f2_name [char]: name of the second set (for error messages)
+%
+% OUTPUTS
+% - (none) This function throws an error if there are elements in `f1` that are not present in `f2`.
+%
+% DESCRIPTION
+% Checks that all elements of `f1` are included in `f2`. If not, throws a descriptive error 
+% mentioning the missing variables and the corresponding structure names for easier debugging.
+function check_consistency(f1, f2, f1_name, f2_name)
+   f1_in_f2 = ismember(f1,f2);
+   if ~all(f1_in_f2)
+      error('Misspecified steady-state input `ss`: the following variables are mentioned in `%s`, but are missing in `%s`: %s', f1_name, f2_name, strjoin(f1(~f1_in_f2)));
+   end
+end
\ No newline at end of file