stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

GSA: Remove globals and clean up map_ident_.m and friends

new-samplers
Johannes Pfeifer 2023-11-28 15:12:24 +01:00
parent 4199f57788
commit 5060f246ea
3 changed files with 61 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
matlab/gsa/map_ident_.m
View File

@ -60,7 +60,7 @@ if opt_gsa.load_ident_files==0
mss = teff(mss(:,istable),Nsam,istable);
yys = teff(yys(dr.order_var,istable),Nsam,istable);
if exist('T','var')
[vdec, cc, ac] = mc_moments(T, lpmatx, dr);
[vdec, cc, ac] = mc_moments(T, lpmatx, dr, M_.exo_nbr, options_);
else
return
end
@ -77,14 +77,14 @@ if opt_gsa.load_ident_files==0
ifig=0;
for j=1:M_.exo_nbr
if mod(j,6)==1
hh_fig=dyn_figure(options_.nodisplay,'name',['Variance decomposition shocks']);
hh_fig=dyn_figure(options_.nodisplay,'name','Variance decomposition shocks');
ifig=ifig+1;
iplo=0;
end
iplo=iplo+1;
subplot(2,3,iplo)
myboxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:size(options_.varobs,1)])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:size(options_.varobs,1))
set(gca,'xlim',[0.5 size(options_.varobs,1)+0.5])
set(gca,'ylim',[-2 102])
for ip=1:size(options_.varobs,1)
@ -95,7 +95,7 @@ if opt_gsa.load_ident_files==0
title(M_.exo_names{j},'interpreter','none')
if mod(j,6)==0 || j==M_.exo_nbr
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],ifig,['Variance decomposition shocks'],'vdec_exo',options_.figures.textwidth*min(iplo/3,1))
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],ifig,'Variance decomposition shocks','vdec_exo',options_.figures.textwidth*min(iplo/3,1))
end
end
end
@ -113,8 +113,8 @@ if opt_gsa.load_ident_files==0
iv = (1:endo_nbr)';
ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
dr.ghx = T(:, [1:(nc1-M_.exo_nbr)],1);
dr.ghu = T(:, [(nc1-M_.exo_nbr+1):end], 1);
dr.ghx = T(:, 1:(nc1-M_.exo_nbr),1);
dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, 1);
[Aa,Bb] = kalman_transition_matrix(dr,iv,ic);
A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
if ~isempty(lpmatx)
@ -122,8 +122,8 @@ if opt_gsa.load_ident_files==0
end
A(:,:,1)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
for j=2:length(istable)
dr.ghx = T(:, [1:(nc1-M_.exo_nbr)],j);
dr.ghu = T(:, [(nc1-M_.exo_nbr+1):end], j);
dr.ghx = T(:, 1:(nc1-M_.exo_nbr),j);
dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, j);
[Aa,Bb] = kalman_transition_matrix(dr, iv, ic);
if ~isempty(lpmatx)
set_shocks_param(lpmatx(j,:));
@ -155,12 +155,12 @@ if opt_gsa.morris==1
SAvdec = squeeze(SAMorris(:,1,:))';
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
load([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec')
end
hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition');
myboxplot(SAvdec,[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
ydum = get(gca,'ylim');
set(gca,'ylim',[0 ydum(2)])
@ -193,7 +193,7 @@ if opt_gsa.morris==1
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: theoretical moments');
myboxplot(SAcc,[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
@ -231,11 +231,11 @@ if opt_gsa.morris==1
end
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAnorm','SAmunorm','SAsignorm','-append')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm','SAmunorm','SAsignorm')
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model');
myboxplot(SAnorm',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
@ -261,7 +261,7 @@ else % main effects analysis
fsuffix = '_rank';
end
imap=[1:npT];
imap=1:npT;
if isempty(lpmat0)
x0=lpmat(istable,:);
@ -309,12 +309,12 @@ else % main effects analysis
save([OutputDirectoryName,'/map_cc',fsuffix,'.mat'],'gsa_')
save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'imap_cc','SAcc','ccac','-append')
else
load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_cc','SAcc','ccac')
load([OutputDirectoryName,'/',fname_,'_main_eff'],'SAcc')
end
hh_fig=dyn_figure(options_.nodisplay,'Name',['Identifiability indices in the ',fsuffix,' moments.']);
bar(sum(SAcc))
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
ydum = get(gca,'ylim');
set(gca,'ylim',[0 ydum(2)])

34
matlab/gsa/mc_moments.m
View File

@ -1,6 +1,18 @@
function [vdec, cc, ac] = mc_moments(mm, ss, dr)
function [vdec, cc, ac] = mc_moments(mm, ss, dr, exo_nbr, options_)
% [vdec, cc, ac] = mc_moments(mm, ss, dr, exo_nbr, options_)
% Conduct Monte Carlo simulation of second moments for GSA
% Inputs:
% - dr [structure] decision rules
% - exo_nbr [double] number of exogenous shocks
% - options_ [structure] Matlab's structure describing the current options
%
% Outputs:
% - vdec [double] variance decomposition matrix
% - cc [double] vector of unique elements of cross correlation matrix
% - ac [cell] autocorrelation matrix
% Copyright © 2012-2018 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
@ -17,27 +29,25 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global options_ M_ estim_params_ oo_
[nr1, nc1, nsam] = size(mm);
[~, nc1, nsam] = size(mm);
nobs=length(options_.varobs);
disp('Computing theoretical moments ...')
disp('mc_moments: Computing theoretical moments ...')
h = dyn_waitbar(0,'Theoretical moments ...');
vdec = zeros(nobs,M_.exo_nbr,nsam);
vdec = zeros(nobs,exo_nbr,nsam);
cc = zeros(nobs,nobs,nsam);
ac = zeros(nobs,nobs*options_.ar,nsam);
for j=1:nsam
oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
oo_.dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
dr.ghx = mm(:, 1:(nc1-exo_nbr),j);
dr.ghu = mm(:, (nc1-exo_nbr+1):end, j);
if ~isempty(ss)
set_shocks_param(ss(j,:));
end
[vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs);
[vdec(:,:,j), corr, autocorr] = th_moments(dr,options_,M_);
cc(:,:,j)=triu(corr);
dum=[];
dum=NaN(nobs,nobs*options_.ar);
for i=1:options_.ar
dum=[dum, autocorr{i}];
dum(:,(i-1)*nobs+1:i*nobs)=autocorr{i};
end
ac(:,:,j)=dum;
dyn_waitbar(j/nsam,h)

42
matlab/gsa/th_moments.m
View File

@ -1,7 +1,22 @@
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_.var_list);
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_.var_list{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