Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
Loading items
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 1151 additions and 676 deletions
matlab/+gsa/map_identification.m 0 → 100644
View file @ 4a32acc2
function map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
% map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
% 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
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fname_ = M_.fname;
dr=oo_.dr;
nliv = opt_gsa.morris_nliv;
itrans = opt_gsa.trans_ident;
np = size(estim_params_.param_vals,1);
pnames = M_.param_names(estim_params_.param_vals(:,1));
if opt_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior'];
else
filetoload=[OutputDirectoryName '/' fname_ '_mc'];
end
load(filetoload,'lpmat','lpmat0','istable','T','yys')
if ~isempty(lpmat0)
lpmatx=lpmat0(istable,:);
else
lpmatx=[];
end
Nsam = size(lpmat,1);
nshock = size(lpmat0,2);
npT = np+nshock;
fname_ = M_.fname;
if opt_gsa.load_ident_files==0
mss = yys(bayestopt_.mfys,:);
mss = gsa.teff(mss(:,istable),Nsam,istable);
yys = gsa.teff(yys(dr.order_var,istable),Nsam,istable);
if exist('T','var')
[vdec, cc, ac] = gsa.monte_carlo_moments(T, lpmatx, dr, M_, options_, estim_params_);
else
return
end
if opt_gsa.morris==2
pdraws = identification.run(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]);
if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0
disp(['Sample check OK. Largest difference: ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
clear pdraws;
end
clear GAM gas
end
if opt_gsa.morris~=1 && M_.exo_nbr>1
ifig=0;
for j=1:M_.exo_nbr
if mod(j,6)==1
hh_fig=dyn_figure(options_.nodisplay,'name','Variance decomposition shocks');
ifig=ifig+1;
iplo=0;
end
iplo=iplo+1;
subplot(2,3,iplo)
gsa.boxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
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)
if options_.TeX
text(ip,-4,deblank(opt_gsa.varobs_tex(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-4,deblank(options_.varobs(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
ylabel(' ')
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))
end
end
end
for j=1:size(cc,1)
cc(j,j,:)=gsa.standardize_columns(squeeze(log(cc(j,j,:))))./2;
end
[vdec, ~, ir_vdec, ic_vdec] = gsa.teff(vdec,Nsam,istable);
[cc, ~, ir_cc, ic_cc] = gsa.teff(cc,Nsam,istable);
[ac, ~, ir_ac, ic_ac] = gsa.teff(ac,Nsam,istable);
nc1= size(T,2);
endo_nbr = M_.endo_nbr;
nstatic = M_.nstatic;
nspred = M_.nspred;
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);
[Aa,Bb] = kalman_transition_matrix(dr,iv,ic);
A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
if ~isempty(lpmatx)
M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(1,:));
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);
[Aa,Bb] = kalman_transition_matrix(dr, iv, ic);
if ~isempty(lpmatx)
M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(j,:));
end
A(:,:,j)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
end
clear T
clear lpmatx
[yt, j0]=gsa.teff(A,Nsam,istable);
yt = [yys yt];
if opt_gsa.morris==2
clear TAU A
else
clear A
end
save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'ac','cc','vdec','yt','mss')
else %load identification files
load([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'ac','cc','vdec','yt','mss')
end
if opt_gsa.morris==1
if ~isempty(vdec)
if opt_gsa.load_ident_files==0
SAMorris=NaN(npT,3,size(vdec,2));
for i=1:size(vdec,2)
[~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], vdec(:,i),nliv);
end
SAvdec = squeeze(SAMorris(:,1,:))';
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec')
end
hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition');
gsa.boxplot(SAvdec,[],'.',[],10)
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)])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:npT
if options_.TeX
[~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_.varobs);
text(ip,-2,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
title('Elementary effects variance decomposition')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morris_vdec'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_morris_vdec'],1,'Screening identification: variance decomposition','morris_vdec',1)
else
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'vdec')
end
if opt_gsa.load_ident_files==0
ccac = [mss cc ac];
SAMorris=NaN(npT,3,size(ccac,2));
for i=1:size(ccac,2)
[~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv);
end
SAcc = squeeze(SAMorris(:,1,:))';
SAcc = SAcc./(max(SAcc,[],2)*ones(1,npT));
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAcc','cc','ir_cc','ic_cc','-append')
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'ac','ir_ac','ic_ac','-append')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAcc','cc','ir_cc','ic_cc')
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'ac','ir_ac','ic_ac')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: theoretical moments');
gsa.boxplot(SAcc,[],'.',[],10)
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])
for ip=1:npT
if options_.TeX
[~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_.varobs);
text(ip,-0.02,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
title('Elementary effects in the moments')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morris_moments'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_morris_moments'],1,'Screening identification: theoretical moments','morris_moments',1)
if opt_gsa.load_ident_files==0
SAMorris=NaN(npT,3,j0);
for j=1:j0
[~, SAMorris(:,:,j)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv);
end
SAM = squeeze(SAMorris(1:end,1,:));
SAnorm=NaN(npT,j0);
irex=NaN(j0);
for j=1:j0
SAnorm(:,j)=SAM(:,j)./max(SAM(:,j));
irex(j)=length(find(SAnorm(:,j)>0.01));
end
SAMmu = squeeze(SAMorris(1:end,2,:));
SAmunorm=NaN(npT,j0);
for j=1:j0
SAmunorm(:,j)=SAMmu(:,j)./max(SAM(:,j)); % normalised w.r.t. mu*
end
SAMsig = squeeze(SAMorris(1:end,3,:));
SAsignorm=NaN(npT,j0);
for j=1:j0
SAsignorm(:,j)=SAMsig(:,j)./max(SAMsig(:,j));
end
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAnorm','SAmunorm','SAsignorm','-append')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model');
gsa.boxplot(SAnorm',[],'.',[],10)
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])
xlabel(' ')
for ip=1:npT
if options_.TeX
[~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_.varobs);
text(ip,-0.02,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
title('Elementary effects in the model')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morris_par'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_morris_par'],1,'Screening identification: model','morris_par',1)
elseif opt_gsa.morris==3
return
elseif opt_gsa.morris==2 % ISKREV stuff
return
else
error('gsa/map_identification: unsupported option morris=%u',opt_gsa.morris)
end
function []=create_TeX_loader(options_,figpath,ifig_number,caption,label_name,scale_factor)
if nargin<6
scale_factor=1;
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([figpath '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by map_ident_.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',scale_factor,strrep(figpath,'\','/'));
fprintf(fidTeX,'\\caption{%s.}',caption);
fprintf(fidTeX,'\\label{Fig:%s:%u}\n',label_name,ifig_number);
fprintf(fidTeX,'\\end{figure}\n\n');
fprintf(fidTeX,'%% End Of TeX file. \n');
fclose(fidTeX);
end
function yr = trank(y)
% yr is the rank transformation of y
yr=NaN(size(y));
[nr, nc] = size(y);
for j=1:nc
[~, is]=sort(y(:,j));
yr(is,j)=[1:nr]'./nr;
end
matlab/gsa/filt_mc_.m → matlab/+gsa/monte_carlo_filtering.m
View file @ 4a32acc2
function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% function [rmse_MC, ixx] = filt_mc_(OutDir)
function [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_)
% [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_
% Inputs:
% - OutputDirectoryName [string] name of the output directory
% - options_gsa_ [structure] GSA options
% - dataset_ [dseries] object storing the dataset
% - dataset_info [structure] storing information about the sample.
% - M_ [structure] MATLAB's structure describing the model
% - oo_ [structure] storing 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:
% - rmse_MC [double] RMSE by nvar matrix of the RMSEs
% - ixx [double] RMSE by nvar matrix of sorting
% indices (descending order of RMSEs)
%
% Notes: the R^2 definition is 1-var(ymodel-ydata)/var(ydata). It ranges
% between (-inf, 1], with negative values indicating that the model is a worse
% predictor than the sample mean of the data
% inputs (from opt_gsa structure)
% vvarvecm = options_gsa_.var_rmse;
% loadSA = options_gsa_.load_rmse;
......@@ -7,14 +27,13 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% alpha = options_gsa_.alpha_rmse;
% alpha2 = options_gsa_.alpha2_rmse;
% istart = options_gsa_.istart_rmse;
% alphaPC = 0.5;
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright (C) 2012-2016 European Commission
% Copyright (C) 2012-2018 Dynare Team
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -29,11 +48,8 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% 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 <http://www.gnu.org/licenses/>.
global bayestopt_ estim_params_ M_ options_ oo_
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% options_gsa_=options_.opt_gsa;
vvarvecm = options_gsa_.var_rmse;
if options_.TeX
vvarvecm_tex = options_gsa_.var_rmse_tex;
......@@ -46,13 +62,10 @@ alpha = options_gsa_.alpha_rmse;
alpha2 = 0;
pvalue = options_gsa_.alpha2_rmse;
istart = max(2,options_gsa_.istart_rmse);
alphaPC = 0.5;
fname_ = M_.fname;
lgy_ = M_.endo_names;
dr_ = oo_.dr;
skipline(2)
skipline(1)
disp('Starting sensitivity analysis')
disp('for the fit of EACH observed series ...')
skipline()
......@@ -61,12 +74,12 @@ if ~options_.nograph
a=dir([OutDir,filesep,'*.*']);
tmp1='0';
if options_.opt_gsa.ppost
tmp=['_rmse_post'];
tmp='_rmse_post';
else
if options_.opt_gsa.pprior
tmp=['_rmse_prior'];
tmp='_rmse_prior';
else
tmp=['_rmse_mc'];
tmp='_rmse_mc';
end
if options_gsa_.lik_only
tmp1 = [tmp,'_post_SA'];
......@@ -75,29 +88,37 @@ if ~options_.nograph
end
for j=1:length(a)
if strmatch([fname_,tmp],a(j).name)
if options_.debug
disp(a(j).name)
end
delete([OutDir,filesep,a(j).name])
end
if strmatch([fname_,tmp1],a(j).name)
if options_.debug
disp(a(j).name)
end
delete([OutDir,filesep,a(j).name])
end
end
disp('done !')
end
[param_names,param_names_tex]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
nshock=estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
npar=estim_params_.np;
if ~isempty(options_.mode_file)
load(options_.mode_file,'xparam1')
end
if options_.opt_gsa.ppost
c=load([fname_,'_mean.mat'],'xparam1');
c=load([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'xparam1');
xparam1_mean=c.xparam1;
xparam1=c.xparam1;
clear c
elseif ~isempty(options_.mode_file) && exist([fname_,'_mean.mat'])==2
c=load([fname_,'_mean.mat'],'xparam1');
elseif ~isempty(options_.mode_file) && exist([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'file')==2
c=load([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'xparam1');
xparam1_mean=c.xparam1;
xparam1=c.xparam1;
clear c
end
......@@ -124,31 +145,11 @@ if loadSA
end
end
if ~loadSA
if exist('xparam1','var')
M_ = set_all_parameters(xparam1,estim_params_,M_);
ys_mode=steady_(M_,options_,oo_);
end
if exist('xparam1_mean','var')
M_ = set_all_parameters(xparam1_mean,estim_params_,M_);
ys_mean=steady_(M_,options_,oo_);
end
Y = transpose(dataset_.data);
gend = dataset_.nobs;
data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state;
for jx=1:gend
data_indx(jx,data_index{jx})=true;
end
load([DirectoryName filesep M_.fname '_data.mat']);
filfilt = dir([DirectoryName filesep M_.fname '_filter_step_ahead*.mat']);
temp_smooth_file_list = dir([DirectoryName filesep M_.fname '_smooth*.mat']);
jfile=0;
for j=1:length(temp_smooth_file_list)
if isempty(strfind(temp_smooth_file_list(j).name,'smoothed')),
jfile=jfile+1;
filsmooth(jfile)=temp_smooth_file_list(j);
end
end
filupdate = dir([DirectoryName filesep M_.fname '_update*.mat']);
filparam = dir([DirectoryName filesep M_.fname '_param*.mat']);
x=[];
......@@ -156,11 +157,11 @@ if ~loadSA
sto_ys=[];
for j=1:length(filparam)
if isempty(strmatch([M_.fname '_param_irf'],filparam(j).name))
load([DirectoryName filesep filparam(j).name]);
x=[x; stock];
logpo2=[logpo2; stock_logpo];
sto_ys=[sto_ys; stock_ys];
clear stock stock_logpo stock_ys;
temp=load([DirectoryName filesep filparam(j).name]); % from prior_posterior_statistics_core
x=[x; temp.stock];
logpo2=[logpo2; temp.stock_logpo];
sto_ys=[sto_ys; temp.stock_ys];
clear temp;
end
end
nruns=size(x,1);
......@@ -168,38 +169,41 @@ if ~loadSA
if options_.opt_gsa.ppost || (options_.opt_gsa.ppost==0 && options_.opt_gsa.lik_only==0)
skipline()
disp('Computing RMSE''s...')
jxj=NaN(length(vvarvecm),1);
js=NaN(length(vvarvecm),1);
yss=NaN(length(vvarvecm),gend,size(sto_ys,1));
for i = 1:length(vvarvecm)
vj = vvarvecm{i};
jxj(i) = strmatch(vj, lgy_(dr_.order_var), 'exact');
js(i) = strmatch(vj, lgy_, 'exact');
jxj(i) = strmatch(vj, M_.endo_names(oo_.dr.order_var), 'exact');
js(i) = strmatch(vj, M_.endo_names, 'exact');
yss(i,:,:)=repmat(sto_ys(:,js(i))',[gend,1]);
end
if exist('xparam1','var')
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);% + kron(ys_mode(js),ones(1,gend)));
yobs = transpose( ahat(jxj,:));% + kron(ys_mode(js),ones(1,gend)));
[~,~,~,ahat,~,~,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);
yobs = transpose( ahat(jxj,:));
rmse_mode = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2));
r2_mode = 1-sum((yobs(istart:end,:)-y0(istart:end,:)).^2)./sum(yobs(istart:end,:).^2);
end
y0=-yss;
y0=-yss; %demean everything using the theoretical mean, i.e. steady state
nbb=0;
for j=1:length(filfilt)
load([DirectoryName filesep M_.fname '_filter_step_ahead',num2str(j),'.mat']);
nb = size(stock,4);
y0(:,:,nbb+1:nbb+nb)=y0(:,:,nbb+1:nbb+nb)+reshape(stock(1,js,1:gend,:),[length(js) gend nb]);
temp=load([DirectoryName filesep M_.fname '_filter_step_ahead',num2str(j),'.mat']);
nb = size(temp.stock,4);
y0(:,:,nbb+1:nbb+nb)=y0(:,:,nbb+1:nbb+nb)+reshape(temp.stock(1,js,1:gend,:),[length(js) gend nb]);
nbb=nbb+nb;
clear stock;
clear temp;
end
yobs=-yss;
nbb=0;
for j=1:length(filupdate)
load([DirectoryName filesep M_.fname '_update',num2str(j),'.mat']);
nb = size(stock,3);
yobs(:,:,nbb+1:nbb+nb)=yobs(:,:,nbb+1:nbb+nb)+reshape(stock(js,1:gend,:),[length(js) gend nb]);
temp=load([DirectoryName filesep M_.fname '_update',num2str(j),'.mat']);
nb = size(temp.stock,3);
yobs(:,:,nbb+1:nbb+nb)=yobs(:,:,nbb+1:nbb+nb)+reshape(temp.stock(js,1:gend,:),[length(js) gend nb]);
nbb=nbb+nb;
clear stock;
clear temp;
end
y0M=mean(y0,2);
rmse_MC=zeros(nruns,length(js));
r2_MC=zeros(nruns,length(js));
for j=1:nruns
......@@ -207,14 +211,15 @@ if ~loadSA
r2_MC(j,:) = 1-mean((yobs(:,istart:end,j)'-y0(:,istart:end,j)').^2)./mean((yobs(:,istart:end,j)').^2);
end
if exist('xparam1_mean','var')
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);% + kron(ys_mean(js),ones(1,gend)));
yobs = transpose( ahat(jxj,:));% + kron(ys_mean(js),ones(1,gend)));
[~,~,~,ahat,~,~,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);
yobs = transpose( ahat(jxj,:));
rmse_pmean = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2));
r2_pmean = 1-mean((yobs(istart:end,:)-y0(istart:end,:)).^2)./mean(yobs(istart:end,:).^2);
end
clear stock_filter;
end
lnprior=NaN(nruns,1);
for j=1:nruns
lnprior(j,1) = priordens(x(j,:)',bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
end
......@@ -242,7 +247,7 @@ if ~loadSA
end
end
end
else
else % loadSA
if options_.opt_gsa.lik_only && options_.opt_gsa.ppost==0
load([OutDir,filesep,fnamtmp, '.mat'],'x','logpo2','likelihood');
else
......@@ -252,27 +257,27 @@ else
nruns=size(x,1);
nfilt=floor(pfilt*nruns);
end
% smirnov tests
% Smirnov tests
nfilt0 = nfilt*ones(length(vvarvecm), 1);
logpo2=logpo2(:);
if ~options_.opt_gsa.ppost
[dum, ipost]=sort(-logpo2);
[dum, ilik]=sort(-likelihood);
[~, ipost]=sort(-logpo2);
[~, ilik]=sort(-likelihood);
end
% visual scatter analysis!
if options_.opt_gsa.ppost
tmp_title='R2 Posterior:';
atitle='R2 Posterior:';
tmp_title='R2 Scatter plot: Posterior';
atitle='R2 Scatter plot: Posterior';
asname='r2_post';
else
if options_.opt_gsa.pprior
tmp_title='R2 Prior:';
atitle='R2 Prior:';
tmp_title='R2 Scatter plot: Prior';
atitle='R2 Scatter plot: Prior';
asname='r2_prior';
else
tmp_title='R2 MC:';
atitle='R2 MC:';
tmp_title='R2 Scatter plot: MC';
atitle='R2 Scatter plot: MC';
asname='r2_mc';
end
end
......@@ -283,7 +288,7 @@ options_scatter.OutputDirectoryName = OutDir;
options_scatter.amcf_name = asname;
options_scatter.amcf_title = atitle;
options_scatter.title = tmp_title;
scatter_analysis(r2_MC, x,options_scatter, options_);
gsa.scatter_analysis(r2_MC, x,options_scatter, options_);
% end of visual scatter analysis
if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
......@@ -297,13 +302,10 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.pvalue_ks = alpha;
options_mcf.pvalue_corr = pvalue;
options_mcf.alpha2 = alpha2;
options_mcf.param_names = param_names;
if options_.TeX
[pnames,pnames_tex]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
options_mcf.param_names = pnames;
options_mcf.param_names_tex = pnames_tex;
options_mcf.param_names_tex = param_names_tex;
else
[pnames]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
options_mcf.param_names = pnames;
options_mcf.param_names_tex = {};
end
options_mcf.fname_ = fname_;
......@@ -313,7 +315,12 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.title = atitle;
options_mcf.beha_title = 'better posterior kernel';
options_mcf.nobeha_title = 'worse posterior kernel';
mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, options_);
if options_.TeX
options_mcf.beha_title_latex = 'better posterior kernel';
options_mcf.nobeha_title_latex = 'worse posterior kernel';
end
gsa.monte_carlo_filtering_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
if options_.opt_gsa.pprior
anam = 'rmse_prior_lik';
atitle = 'RMSE prior: Log Likelihood Kernel';
......@@ -326,15 +333,20 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.title = atitle;
options_mcf.beha_title = 'better likelihood';
options_mcf.nobeha_title = 'worse likelihood';
mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, options_);
if options_.TeX
options_mcf.beha_title_latex = 'better likelihood';
options_mcf.nobeha_title_latex = 'worse likelihood';
end
gsa.monte_carlo_filtering_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
else
if options_.opt_gsa.ppost
rmse_txt=rmse_pmean;
r2_txt=r2_pmean;
else
if options_.opt_gsa.pprior || ~exist('rmse_pmean')
if exist('rmse_mode')
if options_.opt_gsa.pprior || ~exist('rmse_pmean','var')
if exist('rmse_mode','var')
rmse_txt=rmse_mode;
r2_txt=r2_mode;
else
......@@ -346,18 +358,19 @@ else
r2_txt=r2_pmean;
end
end
ixx=NaN(size(rmse_MC,1),length(vvarvecm));
for i = 1:length(vvarvecm)
[dum, ixx(:,i)] = sort(rmse_MC(:,i));
[~, ixx(:,i)] = sort(rmse_MC(:,i));
end
PP = ones(npar+nshock, length(vvarvecm));
PPV = ones(length(vvarvecm), length(vvarvecm), npar+nshock);
SS = zeros(npar+nshock, length(vvarvecm));
for j = 1:npar+nshock
for i = 1:length(vvarvecm)
[H, P, KSSTAT] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
[H1, P1, KSSTAT1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
[H2, P2, KSSTAT2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
if H1 & H2==0
[~, P] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
[H1] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
[H2] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
if H1==0 && H2==0
SS(j,i)=1;
elseif H1==0
SS(j,i)=-1;
......@@ -369,7 +382,7 @@ else
for i = 1:length(vvarvecm)
for l = 1:length(vvarvecm)
if l~=i && PP(j,i)<alpha && PP(j,l)<alpha
[H,P,KSSTAT] = smirnov(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
[~,P] = gsa.smirnov_test(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
PPV(i,l,j) = P;
elseif l==i
PPV(i,l,j) = PP(j,i);
......@@ -391,34 +404,38 @@ else
end
if mod(i,9)==1
ifig=ifig+1;
hh=dyn_figure(options_.nodisplay,'name',[temp_name,' ',int2str(ifig)]);
hh_fig=dyn_figure(options_.nodisplay,'name',[temp_name,' ',int2str(ifig)]);
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(lnprior(ixx(1:nfilt0(i),i)));
h=gsa.cumplot(lnprior(ixx(1:nfilt0(i),i)));
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(lnprior);
hold on, h=gsa.cumplot(lnprior);
set(h,'color','k','linewidth',1)
h=cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
h=gsa.cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
set(h,'color','red','linewidth',2)
if options_.TeX
title(vvarvecm_tex{i},'interpreter','latex')
else
title(vvarvecm{i},'interpreter','none')
end
if mod(i,9)==0 || i==length(vvarvecm)
if ~isoctave
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Log-prior for BETTER R2','Color','Blue','horizontalalignment','center');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Log-prior for WORSE R2', 'Color','Red','horizontalalignment','center');
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_post_lnprior',int2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_post_lnprior',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir '/' fname_ '_rmse_post_lnprior',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_post_lnprior',options_.figures.textwidth*min((i-9*(ifig-1))/3,1))
end
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_lnprior',int2str(ifig) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_prior_lnprior',int2str(ifig) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir '/' fname_ '_rmse_prior_lnprior',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_prior_lnprior',options_.figures.textwidth*min((i-9*(ifig-1))/3,1))
end
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_lnprior',int2str(ifig) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_mc_lnprior',int2str(ifig) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir '/' fname_ '_rmse_mc_lnprior',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_mc_lnprior',options_.figures.textwidth*min((i-9*(ifig-1))/3,1))
end
......@@ -439,16 +456,20 @@ else
end
if mod(i,9)==1
ifig=ifig+1;
hh = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
hh_fig = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(likelihood(ixx(1:nfilt0(i),i)));
h=gsa.cumplot(likelihood(ixx(1:nfilt0(i),i)));
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(likelihood);
hold on, h=gsa.cumplot(likelihood);
set(h,'color','k','linewidth',1)
h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
h=gsa.cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
set(h,'color','red','linewidth',2)
if options_.TeX
title(vvarvecm_tex{i},'interpreter','latex')
else
title(vvarvecm{i},'interpreter','none')
end
if options_.opt_gsa.ppost==0
set(gca,'xlim',[min( likelihood(ixx(1:nfilt0(i),i)) ) max( likelihood(ixx(1:nfilt0(i),i)) )])
end
......@@ -458,18 +479,18 @@ else
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Log-likelihood for WORSE R2', 'Color','Red','horizontalalignment','center');
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_post_lnlik',int2str(ifig) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_post_lnlik',int2str(ifig) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_post_lnprior',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_post_lnprior',options_.figures.textwidth*min((i-9*(ifig-1))/3,1));
end
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_lnlik',int2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_prior_lnlik',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_lnlik',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_prior_lnlik',options_.figures.textwidth*min((i-9*(ifig-1))/3,1));
end
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_lnlik',int2str(ifig) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_mc_lnlik',int2str(ifig) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_lnlik',int2str(ifig) ],ifig,[temp_name,' ',int2str(ifig)],'rmse_mc_lnlik',options_.figures.textwidth*min((i-9*(ifig-1))/3,1));
end
......@@ -490,16 +511,20 @@ else
end
if mod(i,9)==1
ifig=ifig+1;
hh = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
hh_fig = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(logpo2(ixx(1:nfilt0(i),i)));
h=gsa.cumplot(logpo2(ixx(1:nfilt0(i),i)));
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(logpo2);
hold on, h=gsa.cumplot(logpo2);
set(h,'color','k','linewidth',1)
h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
h=gsa.cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
set(h,'color','red','linewidth',2)
if options_.TeX
title(vvarvecm_tex{i},'interpreter','latex')
else
title(vvarvecm{i},'interpreter','none')
end
if options_.opt_gsa.ppost==0
set(gca,'xlim',[min( logpo2(ixx(1:nfilt0(i),i)) ) max( logpo2(ixx(1:nfilt0(i),i)) )])
end
......@@ -509,18 +534,18 @@ else
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Log-posterior for WORSE R2', 'Color','Red','horizontalalignment','center');
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_post_lnpost',int2str(ifig) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_post_lnpost',int2str(ifig) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_post_lnpost',int2str(ifig) ],ifig,[temp_name,' ',int2str(ifig)],'rmse_post_lnpost',options_.figures.textwidth*min((i-9*(ifig-1))/3,1));
end
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_lnpost',int2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_prior_lnpost',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_lnpost',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_prior_lnpost',options_.figures.textwidth*min((i-9*(ifig-1))/3,1));
end
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_lnpost',int2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_mc_lnpost',int2str(ifig)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_lnpost',int2str(ifig)],ifig,[temp_name,' ',int2str(ifig)],'rmse_mc_lnpost',options_.figures.textwidth*min((i-9*(ifig-1))/3,1));
end
......@@ -529,15 +554,6 @@ else
end
end
end
if options_.TeX
[pnames,pnames_tex]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
param_names = pnames;
param_names_tex = pnames_tex;
else
[pnames]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
param_names = pnames;
param_names_tex = {};
end
skipline()
title_string='RMSE over the MC sample:';
data_mat=[min(rmse_MC)' max(rmse_MC)'];
......@@ -549,7 +565,7 @@ else
end
invar = find( std(rmse_MC)./mean(rmse_MC)<=0.0001 );
if ~isempty(invar)
skipline(2)
skipline(1)
disp('RMSE is not varying significantly over the MC sample for the following variables:')
disp(vvarvecm{invar})
disp('These variables are excluded from SA')
......@@ -561,8 +577,7 @@ else
rmse_MC = rmse_MC(:,ivar);
skipline()
disp(['Sample filtered the ',num2str(pfilt*100),'% best RMSE''s for each observed series ...' ])
skipline(2)
disp('RMSE ranges after filtering:')
skipline(1)
title_string='RMSE ranges after filtering:';
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
headers = {'Variable'; 'min'; 'max'; 'min'; 'max'; 'posterior mode'};
......@@ -589,7 +604,7 @@ else
else
values_length = max(ceil(max(max(log10(abs(data_mat(isfinite(data_mat))))))),1)+val_precis+1;
end
if any(data_mat) < 0 %add one character for minus sign
if any(data_mat < 0) %add one character for minus sign
values_length = values_length+1;
end
headers_length = cellofchararraymaxlength(headers(2:end));
......@@ -598,7 +613,6 @@ else
else
val_width = max(headers_length, values_length)+2;
end
value_format = sprintf('%%%d.%df',val_width,val_precis);
header_string_format = sprintf('%%%ds',val_width);
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
optional_header=sprintf([label_format_leftbound,header_string_format,header_string_format,header_string_format,header_string_format],'','',['best ',num2str(pfilt*100),'% filtered'],'','remaining 90%');
......@@ -610,7 +624,7 @@ else
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
optional_header={[' & \multicolumn{2}{c}{best ',num2str(pfilt*100),' filtered} & \multicolumn{2}{c}{remaining 90\%}\\']};
else
optional_header={[' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\']};
optional_header={' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\'};
end
dyn_latex_table(M_, options_, title_string, 'RMSE_ranges_after_filtering', headers_tex, vvarvecm_tex, data_mat, 0, val_width, val_precis, optional_header);
end
......@@ -657,7 +671,7 @@ else
else
values_length = max(ceil(max(max(log10(abs(data_mat(isfinite(data_mat))))))),1)+val_precis+1;
end
if any(data_mat) < 0 %add one character for minus sign
if any(data_mat < 0) %add one character for minus sign
values_length = values_length+1;
end
headers_length = cellofchararraymaxlength(headers(2:end));
......@@ -666,7 +680,6 @@ else
else
val_width = max(headers_length, values_length)+2;
end
value_format = sprintf('%%%d.%df',val_width,val_precis);
header_string_format = sprintf('%%%ds',val_width);
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
......@@ -679,7 +692,7 @@ else
if ~options_.opt_gsa.ppost && options_.opt_gsa.pprior
optional_header = {[' & \multicolumn{2}{c}{best ',num2str(pfilt*100),' filtered} & \multicolumn{2}{c}{remaining 90\%}\\']};
else
optional_header = {[' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\']};
optional_header = {' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\'};
end
dyn_latex_table(M_, options_, title_string, 'R2_ranges_after_filtering', headers_tex, vvarvecm_tex, data_mat, 0, val_width, val_precis, optional_header);
end
......@@ -690,16 +703,15 @@ else
SP(ns,j)=ones(size(ns));
SS(:,j)=SS(:,j).*SP(:,j);
end
for j=1:npar+nshock %estim_params_.np,
nsp=NaN(npar+nshock,1);
for j=1:npar+nshock
nsp(j)=length(find(SP(j,:)));
end
snam0=param_names(find(nsp==0));
snam1=param_names(find(nsp==1));
snam2=param_names(find(nsp>1));
snam=param_names(find(nsp>0));
snam0=param_names(nsp==0);
snam1=param_names(nsp==1);
snam2=param_names(nsp>1);
nsnam=(find(nsp>1));
skipline(2)
skipline(1)
disp('These parameters do not affect significantly the fit of ANY observed series:')
disp(char(snam0))
skipline()
......@@ -708,7 +720,6 @@ else
skipline()
disp('These parameters affect MORE THAN ONE observed series: trade off exists!')
disp(char(snam2))
pnam=bayestopt_.name;
% plot trade-offs
if ~options_.nograph
a00=jet(length(vvarvecm));
......@@ -740,62 +751,74 @@ else
options_mcf.amcf_title = [atitle ' ' vvarvecm{iy}];
options_mcf.beha_title = ['better fit of ' vvarvecm{iy}];
options_mcf.nobeha_title = ['worse fit of ' vvarvecm{iy}];
if options_.TeX
options_mcf.beha_title_latex = ['better fit of ' vvarvecm_tex{iy}];
options_mcf.nobeha_title_latex = ['worse fit of ' vvarvecm_tex{iy}];
end
options_mcf.title = ['the fit of ' vvarvecm{iy}];
mcf_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, options_);
gsa.monte_carlo_filtering_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, M_, options_, bayestopt_, estim_params_);
end
for iy = 1:length(vvarvecm)
ipar = find(any(squeeze(PPV(iy,:,:))<alpha));
for ix=1:ceil(length(ipar)/5)
hh = dyn_figure(options_.nodisplay,'name',[temp_name,' observed variable ', vvarvecm{iy}]);
hh_fig = dyn_figure(options_.nodisplay,'name',[temp_name,' observed variable ', vvarvecm{iy}]);
for j=1+5*(ix-1):min(length(ipar),5*ix)
subplot(2,3,j-5*(ix-1))
h0=cumplot(x(:,ipar(j)));
h0=gsa.cumplot(x(:,ipar(j)));
set(h0,'color',[0 0 0])
hold on,
iobs=find(squeeze(PPV(iy,:,ipar(j)))<alpha);
for i = 1:length(vvarvecm)
if any(iobs==i) || i==iy
h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
if ~isoctave
hcmenu = uicontextmenu;
uimenu(hcmenu,'Label',vvarvecm{i});
set(h0,'uicontextmenu',hcmenu)
end
else
h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
end
set(h0,'color',a00(i,:),'linewidth',2)
end
ydum=get(gca,'ylim');
if exist('xparam1')
if exist('xparam1','var')
xdum=xparam1(ipar(j));
h1=plot([xdum xdum],ydum);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
end
xlabel('')
title([pnam{ipar(j)}],'interpreter','none')
if options_.TeX
title([param_names_tex{ipar(j)}],'interpreter','latex')
else
title([param_names{ipar(j)}],'interpreter','none')
end
end
if isoctave
legend(vertcat('base',vvarvecm),'location','eastoutside');
else
h0=legend(vertcat('base',vvarvecm));
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none');
if options_.TeX
h0=legend(vertcat('base',vvarvecm_tex),'interpreter','latex');
else
h0=legend(vertcat('base',vvarvecm),'interpreter','none');
end
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3]);
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable $',vvarvecm_tex{iy} '$'],['rmse_post_' vvarvecm{iy}],1)
create_TeX_loader(options_,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable ',vvarvecm_tex{iy} ],['rmse_post_' vvarvecm{iy}],1)
end
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],ix,[temp_name,' observed variable $',vvarvecm_tex{iy} '$'],['rmse_prior_' vvarvecm{iy}],1)
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],ix,[temp_name,' observed variable ',vvarvecm_tex{iy}],['rmse_prior_' vvarvecm{iy}],1)
end
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable $',vvarvecm_tex{iy} '$'],['rmse_mc_' vvarvecm{iy}],1)
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable ',vvarvecm_tex{iy}],['rmse_mc_' vvarvecm{iy}],1)
end
end
end
......@@ -803,18 +826,18 @@ else
end
% now I plot by individual parameters
for ix=1:ceil(length(nsnam)/5)
hh = dyn_figure(options_.nodisplay,'name',[temp_name,' estimated params and shocks ',int2str(ix)]);
hh_fig = dyn_figure(options_.nodisplay,'name',[temp_name,' estimated params and shocks ',int2str(ix)]);
for j=1+5*(ix-1):min(size(snam2,1),5*ix)
subplot(2,3,j-5*(ix-1))
h0=cumplot(x(:,nsnam(j)));
h0=gsa.cumplot(x(:,nsnam(j)));
set(h0,'color',[0 0 0])
hold on,
npx=find(SP(nsnam(j),:)==0);
for i = 1:length(vvarvecm)
if any(npx==i)
h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
else
h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
if ~isoctave
hcmenu = uicontextmenu;
uimenu(hcmenu,'Label', vvarvecm{i});
......@@ -824,34 +847,42 @@ else
set(h0,'color',a00(i,:),'linewidth',2)
end
ydum=get(gca,'ylim');
if exist('xparam1')
if exist('xparam1','var')
xdum=xparam1(nsnam(j));
h1=plot([xdum xdum],ydum);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
end
xlabel('')
title([pnam{nsnam(j)}],'interpreter','none')
if options_.TeX
title([param_names_tex{nsnam(j)}],'interpreter','latex')
else
title([param_names{nsnam(j)}],'interpreter','none')
end
end
%subplot(3,2,6)
if isoctave
legend(vertcat('base',vvarvecm),'location','eastoutside');
else
h0=legend(vertcat('base',vvarvecm));
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none');
if options_.TeX
h0=legend(vertcat('base',vvarvecm_tex),'interpreter','latex');
else
h0=legend(vertcat('base',vvarvecm),'interpreter','none');
end
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3]);
end
if options_.opt_gsa.ppost
dyn_saveas(hh,[ OutDir filesep fname_ '_rmse_post_params_' int2str(ix)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[ OutDir filesep fname_ '_rmse_post_params_' int2str(ix)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[ OutDir filesep fname_ '_rmse_post_params_' int2str(ix)],ix,[temp_name,' estimated params and shocks ',int2str(ix)],'rmse_post_params',1)
end
else
if options_.opt_gsa.pprior
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_prior_params_' int2str(ix) ],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_prior_params_' int2str(ix) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_params_' int2str(ix) ],ix,[temp_name,' estimated params and shocks ',int2str(ix)],'rmse_prior_params',1)
end
else
dyn_saveas(hh,[OutDir filesep fname_ '_rmse_mc_params_' int2str(ix)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_mc_params_' int2str(ix)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_params_' int2str(ix)],ix,[temp_name,' estimated params and shocks ',int2str(ix)],'rmse_mc_params',1)
end
......@@ -885,11 +916,11 @@ pnames=cell(np,1);
pnames_tex=cell(np,1);
for ii=1:length(bayestopt_.name)
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(ii,options_.TeX,M_,estim_params_,options_);
pnames_tex{ii,1} = strrep(param_name_tex_temp,'$','');
[param_name_temp, param_name_tex_temp]= get_the_name(ii,options_.TeX,M_,estim_params_,options_.varobs);
pnames_tex{ii,1} = param_name_tex_temp;
pnames{ii,1} = param_name_temp;
else
param_name_temp = get_the_name(ii,options_.TeX,M_,estim_params_,options_);
param_name_temp = get_the_name(ii,options_.TeX,M_,estim_params_,options_.varobs);
pnames{ii,1} = param_name_temp;
end
end
matlab/gsa/mcf_analysis.m → matlab/+gsa/monte_carlo_filtering_analysis.m
View file @ 4a32acc2
function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, DynareOptions)
function indmcf = monte_carlo_filtering_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
% indmcf = monte_carlo_filtering_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
% Inputs:
% - lpmat [double] Monte Carlo matrix
% - ibeha [integer] index of behavioural runs
% - inobeha [integer] index of non-behavioural runs
% - options_gsa_ [structure] GSA options_
% - M_ [structure] describing the model
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
% - estim_params_ [structure] characterizing parameters to be estimated
%
% Outputs:
% - indmcf [double] results of matrix
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
%
% Copyright (C) 2014 European Commission
% Copyright (C) 2016-2018 Dynare Team
% Copyright © 2014 European Commission
% Copyright © 2016-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -21,14 +34,14 @@ function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, DynareOptions
% 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 <http://www.gnu.org/licenses/>.
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
pvalue_ks = options_mcf.pvalue_ks;
pvalue_corr = options_mcf.pvalue_corr;
alpha2 = options_mcf.alpha2;
param_names = options_mcf.param_names;
if DynareOptions.TeX
if options_.TeX
if ~isfield(options_mcf,'param_names_tex')
param_names_tex = options_mcf.param_names;
else
......@@ -41,6 +54,10 @@ amcf_name = options_mcf.amcf_name;
amcf_title = options_mcf.amcf_title;
beha_title = options_mcf.beha_title;
nobeha_title = options_mcf.nobeha_title;
if options_.TeX
beha_title_latex = options_mcf.beha_title_latex;
nobeha_title_latex = options_mcf.nobeha_title_latex;
end
title = options_mcf.title;
fname_ = options_mcf.fname_;
xparam1=[];
......@@ -49,18 +66,18 @@ if isfield(options_mcf,'xparam1')
end
OutputDirectoryName = options_mcf.OutputDirectoryName;
[proba, dproba] = stab_map_1(lpmat, ibeha, inobeha, [],0);
[proba, dproba] = gsa.stability_mapping_univariate(lpmat, ibeha, inobeha, [],fname_, options_, bayestopt_.name, estim_params_,0);
indmcf=find(proba<pvalue_ks);
[tmp,jtmp] = sort(proba(indmcf),2,'ascend');
[~,jtmp] = sort(proba(indmcf),1,'ascend');
indmcf = indmcf(jtmp);
if ~isempty(indmcf)
skipline()
headers = {'Parameter','d-stat','p-value'};
labels = param_names(indmcf);
data_mat=[dproba(indmcf)' proba(indmcf)'];
data_mat=[dproba(indmcf) proba(indmcf)];
options_temp.noprint=0;
dyntable(options_temp,['Smirnov statistics in driving ', title],headers,labels,data_mat,size(labels,2)+2,16,3);
if DynareOptions.TeX
if options_.TeX
labels_TeX=param_names_tex(indmcf);
M_temp.dname=OutputDirectoryName ;
M_temp.fname=fname_;
......@@ -68,19 +85,31 @@ if ~isempty(indmcf)
end
end
if length(ibeha)>10 && length(inobeha)>10
indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, beha_title);
indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, nobeha_title);
if options_.TeX
indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title, beha_title_latex);
indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title, nobeha_title_latex);
else
indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title);
indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title);
end
indcorr = union(indcorr1(:), indcorr2(:));
indcorr = indcorr(~ismember(indcorr(:),indmcf));
indmcf = [indmcf(:); indcorr(:)];
end
if ~isempty(indmcf) && ~DynareOptions.nograph
if ~isempty(indmcf) && ~options_.nograph
skipline()
xx=[];
if ~ isempty(xparam1), xx=xparam1(indmcf); end
scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, DynareOptions, ...
if ~ isempty(xparam1)
xx=xparam1(indmcf);
end
if options_.TeX
gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
beha_title, nobeha_title, beha_title_latex, nobeha_title_latex)
else
gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
beha_title, nobeha_title)
end
end
matlab/gsa/mc_moments.m → matlab/+gsa/monte_carlo_moments.m
View file @ 4a32acc2
function [vdec, cc, ac] = mc_moments(mm, ss, dr)
function [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_, estim_params_)
% [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_,estim_params_)
% Conduct Monte Carlo simulation of second moments for GSA
% Inputs:
% - dr [structure] decision rules
% - M_ [structure] model structure
% - options_ [structure] MATLAB's structure describing the current options
% - estim_params_ [structure] characterizing parameters to be estimated
%
% Outputs:
% - vdec [double] variance decomposition matrix
% - cc [double] vector of unique elements of cross correlation matrix
% - ac [cell] autocorrelation matrix
% Copyright (C) 2012-2018 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -15,33 +28,33 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
% 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 <http://www.gnu.org/licenses/>.
% 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('monte_carlo_moments: Computing theoretical moments ...')
h = dyn_waitbar(0,'Theoretical moments ...');
vdec = zeros(nobs,M_.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-M_.exo_nbr),j);
dr.ghu = mm(:, (nc1-M_.exo_nbr+1):end, j);
if ~isempty(ss)
set_shocks_param(ss(j,:));
M_=gsa.set_shocks_param(M_,estim_params_,ss(j,:));
end
[vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs);
[vdec(:,:,j), corr, autocorr] = gsa.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;
if mod(j,3)==0
dyn_waitbar(j/nsam,h)
end
end
dyn_waitbar_close(h)
skipline()
disp('... done !')
matlab/gsa/prior_draw_gsa.m → matlab/+gsa/prior_draw.m
View file @ 4a32acc2
function pdraw = prior_draw_gsa(init,rdraw)
function pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw)
% pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw)
% Draws from the prior distributions for use with Sensitivity Toolbox for DYNARE
%
% INPUTS
% o init [integer] scalar equal to 1 (first call) or 0.
% o rdraw
% - M_ [structure] describing the model
% - bayestopt_ [structure] describing the priors
% - options_ [structure] describing the options
% - estim_params_ [structure] characterizing parameters to be estimated
% - init [integer] scalar equal to 1 (first call) or 0.
% - rdraw
%
% OUTPUTS
% o pdraw [double] draw from the joint prior density.
......@@ -18,8 +23,8 @@ function pdraw = prior_draw_gsa(init,rdraw)
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright (C) 2012-2015 European Commission
% Copyright (C) 2012-2017 Dynare Team
% Copyright © 2012-2015 European Commission
% Copyright © 2012-2017 Dynare Team
%
% This file is part of Dynare.
%
......@@ -34,9 +39,8 @@ function pdraw = prior_draw_gsa(init,rdraw)
% 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 <http://www.gnu.org/licenses/>.
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global bayestopt_ options_ estim_params_ M_
persistent npar pshape p6 p7 p3 p4 lbcum ubcum
if init
......@@ -49,7 +53,7 @@ if init
pdraw = zeros(npar,1);
lbcum = zeros(npar,1);
ubcum = ones(npar,1);
[~,~,~,lb,ub,~] = set_prior(estim_params_,M_,options_); %Prepare bounds
[~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Set prior bounds
bounds = prior_bounds(bayestopt_, options_.prior_trunc);
......@@ -64,29 +68,29 @@ if init
% set bounds for cumulative probabilities
for i = 1:npar
switch pshape(i)
case 5% Uniform prior.
p4(i) = min(p4(i),bounds.ub(i));
p3(i) = max(p3(i),bounds.lb(i));
case 3% Gaussian prior.
lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 2% Gamma prior.
lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
case 1% Beta distribution (TODO: generalized beta distribution)
lbcum(i) = betainc((bounds.lb(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
ubcum(i) = betainc((bounds.ub(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
case 2% Gamma prior.
lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
case 3% Gaussian prior.
lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 4% INV-GAMMA1 distribution
% TO BE CHECKED
lbcum(i) = gamcdf(1/(bounds.ub(i)-p3(i))^2,p7(i)/2,2/p6(i));
ubcum(i) = gamcdf(1/(bounds.lb(i)-p3(i))^2,p7(i)/2,2/p6(i));
case 5% Uniform prior.
p4(i) = min(p4(i),bounds.ub(i));
p3(i) = max(p3(i),bounds.lb(i));
case 6% INV-GAMMA2 distribution
% TO BE CHECKED
lbcum(i) = gamcdf(1/(bounds.ub(i)-p3(i)),p7(i)/2,2/p6(i));
ubcum(i) = gamcdf(1/(bounds.lb(i)-p3(i)),p7(i)/2,2/p6(i));
case 8
lbcum(i) = weibcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = weibcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
lbcum(i) = wblcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = wblcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
otherwise
% Nothing to do here.
end
......@@ -94,7 +98,7 @@ if init
return
end
pdraw=NaN(size(rdraw,1),npar);
for i = 1:npar
rdraw(:,i) = rdraw(:,i).*(ubcum(i)-lbcum(i))+lbcum(i);
switch pshape(i)
......
matlab/gsa/priorcdf.m → matlab/+gsa/priorcdf.m
View file @ 4a32acc2
function xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% This procedure transforms x vectors into cumulative values
% pshape: 0 is point mass, both para and p2 are ignored
% 1 is BETA(mean,stdd)
......@@ -11,7 +11,7 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% 8 is WEIBULL(s, k)
% Adapted by M. Ratto from MJ priordens.m
% Copyright (C) 2012-2015 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -26,8 +26,9 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% 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 <http://www.gnu.org/licenses/>.
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
xcum=NaN(size(para));
for i=1:length(pshape)
switch pshape(i)
case 1 % (generalized) BETA Prior
......
matlab/gsa/redform_map.m → matlab/+gsa/reduced_form_mapping.m
View file @ 4a32acc2
function redform_map(dirname,options_gsa_)
%function redform_map(dirname)
% inputs (from opt_gsa structure
% anamendo = options_gsa_.namendo;
% anamlagendo = options_gsa_.namlagendo;
% anamexo = options_gsa_.namexo;
% iload = options_gsa_.load_redform;
% pprior = options_gsa_.pprior;
% ilog = options_gsa_.logtrans_redform;
% threshold = options_gsa_.threshold_redform;
% ksstat = options_gsa_.ksstat_redform;
% alpha2 = options_gsa_.alpha2_redform;
function reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
% reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
% Inputs:
% - dirname [string] name of the output directory
% - options_gsa_ [structure] GSA options_
% - M_ [structure] describing the model
% - estim_params_ [structure] characterizing parameters to be estimated
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
% - oo_ [structure] storing the results
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright (C) 2012-2016 European Commission
% Copyright (C) 2012-2018 Dynare Team
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -31,25 +29,18 @@ function redform_map(dirname,options_gsa_)
% 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 <http://www.gnu.org/licenses/>.
global M_ oo_ estim_params_ options_ bayestopt_
% options_gsa_ = options_.opt_gsa;
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
anamexo = options_gsa_.namexo;
anamendo_tex = options_gsa_.namendo_tex;
anamlagendo_tex = options_gsa_.namlagendo_tex;
anamexo_tex = options_gsa_.namexo_tex;
iload = options_gsa_.load_redform;
pprior = options_gsa_.pprior;
ilog = options_gsa_.logtrans_redform;
threshold = options_gsa_.threshold_redform;
% ksstat = options_gsa_.ksstat_redform;
alpha2 = options_gsa_.alpha2_redform;
alpha2=0;
pvalue_ks = options_gsa_.ksstat_redform;
pvalue_corr = options_gsa_.alpha2_redform;
np = estim_params_.np;
nshock = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
......@@ -57,11 +48,11 @@ pnames=cell(np,1);
pnames_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_);
pnames_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);
pnames_tex{jj,1} = param_name_tex_temp;
pnames{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);
pnames{jj,1} = param_name_temp;
end
end
......@@ -93,14 +84,14 @@ end
options_mcf.fname_ = M_.fname;
options_mcf.OutputDirectoryName = adir;
if ~exist('T')
stab_map_(dirname,options_gsa_);
if ~exist('T','var')
gsa.stability_mapping(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_);
if pprior
load([dirname,filesep,M_.fname,'_prior'],'T');
else
load([dirname,filesep,M_.fname,'_mc'],'T');
end
if ~exist('T')
if ~exist('T','var')
disp('The model is too large!')
disp('Reduced form mapping stopped!')
return
......@@ -109,8 +100,6 @@ end
if isempty(dir(adir))
mkdir(adir)
end
adir0=pwd;
%cd(adir)
nspred=size(T,2)-M_.exo_nbr;
x0=lpmat(istable,:);
......@@ -121,7 +110,7 @@ else
xx0=lpmat0(istable,:);
nshocks=size(xx0,2);
end
[kn, np]=size(x0);
[~, np]=size(x0);
offset = length(bayestopt_.pshape)-np;
if options_gsa_.prior_range
pshape=5*(ones(np,1));
......@@ -144,28 +133,27 @@ options_map.pshape = pshape;
options_map.pd = pd;
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
lpmat=[];
lpmat0=[];
js=0;
for j = 1:length(anamendo)
namendo = anamendo{j};
namendo_tex = anamendo_tex{j};
iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact');
ifig = 0;
iplo = 0;
for jx = 1:length(anamexo)
namexo = anamexo{jx};
namexo_tex = anamexo_tex{jx};
iexo=strmatch(namexo, M_.exo_names, 'exact');
skipline()
disp(['[', namendo,' vs ',namexo,']'])
if ~isempty(iexo)
%y0=squeeze(T(iendo,iexo+nspred,istable));
y0=squeeze(T(iendo,iexo+nspred,:));
if (max(y0)-min(y0))>1.e-10
if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph
ifig=ifig+1;
hfig = dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping: ', namendo,' vs shocks ',int2str(ifig)]);
hh_fig = dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping: ', namendo,' vs shocks ',int2str(ifig)]);
iplo=0;
end
iplo=iplo+1;
......@@ -194,20 +182,23 @@ for j = 1:length(anamendo)
end
if ~options_.nograph
hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs ', namexo]);
hc = cumplot(y0);
hc = gsa.cumplot(y0);
a=axis; delete(hc);
% hist(mat_moment{ij}),
x1val=max(threshold(1),a(1));
x2val=min(threshold(2),a(2));
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
set(hp,'FaceColor', [0.7 0.8 1])
hold all,
hc = cumplot(y0);
hc = gsa.cumplot(y0);
set(hc,'color','k','linewidth',2)
hold off,
if options_.TeX
title([namendo_tex,' vs ', namexo_tex ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','latex')
else
title([namendo,' vs ', namexo ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','none')
end
dyn_saveas(hf,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],['Reduced Form Mapping (Monte Carlo Filtering): ',strrep(namendo,'_','\_'),' vs ', strrep(namexo,'_','\_')],[type '_' namendo,'_vs_', namexo])
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],['Reduced Form Mapping (Monte Carlo Filtering): ',namendo_tex,' vs ', namexo_tex],[type '_' namendo,'_vs_', namexo])
end
si(:,js) = NaN(np,1);
delete([xdir, '/*threshold*.*'])
......@@ -219,19 +210,23 @@ for j = 1:length(anamendo)
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'inside threshold';
options_mcf.nobeha_title = 'outside threshold';
if options_.TeX
options_mcf.beha_title_latex = 'inside threshold';
options_mcf.nobeha_title_latex = 'outside threshold';
end
options_mcf.title = atitle0;
options_mcf.OutputDirectoryName = xdir;
if ~isempty(iy) && ~isempty(iyc)
fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_);
icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
lpmat=x0(iy,:);
if nshocks
lpmat0=xx0(iy,:);
end
istable=[1:length(iy)];
istable=1:length(iy);
save([xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo '_threshold' ],'lpmat','lpmat0','istable','y0','x0','xx0','iy','iyc')
lpmat=[]; lpmat0=[]; istable=[];
lpmat0=[];
if length(iy)<=10 || length(iyc)<=10
icheck = []; % do the generic plot in any case
end
......@@ -255,12 +250,10 @@ for j = 1:length(anamendo)
end
atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namexo];
options_mcf.title = atitle0;
indmcf = redform_mcf(y0, x0, options_mcf, options_);
redform_mcf(y0, x0, options_mcf, options_, M_.fname, bayestopt_.name, estim_params_);
end
end
else
[yy, xdir] = log_trans_(y0,xdir0);
atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namexo];
aname=[type '_' namendo '_vs_' namexo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namexo];
......@@ -273,27 +266,34 @@ for j = 1:length(anamendo)
end
if isempty(threshold) && ~options_.nograph
figure(hfig)
figure(hh_fig)
subplot(3,3,iplo),
if ilog
[saso, iso] = sort(-silog(:,js));
[~, iso] = sort(-silog(:,js));
bar([silog(iso(1:min(np,10)),js)])
logflag='log';
else
[saso, iso] = sort(-si(:,js));
[~, iso] = sort(-si(:,js));
bar(si(iso(1:min(np,10)),js))
logflag='';
end
%set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
if options_.TeX
text(ip,-0.02,deblank(pnames_tex(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if options_.TeX
title([logflag,' ',namendo_tex,' vs ',namexo_tex],'interpreter','none')
else
title([logflag,' ',namendo,' vs ',namexo],'interpreter','none')
end
if iplo==9
dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namexo,'_','\_')],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],1)
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',namendo_tex,' vs ',namexo_tex],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],1)
end
end
......@@ -303,24 +303,24 @@ for j = 1:length(anamendo)
end
end
if iplo<9 && iplo>0 && ifig && ~options_.nograph
dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namexo,'_','\_')],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.figures.textwidth*min(iplo/3,1))
end
ifig=0;
iplo=0;
for je=1:length(anamlagendo)
namlagendo = anamlagendo{je};
namlagendo_tex = anamlagendo_tex{je};
ilagendo=strmatch(namlagendo, M_.endo_names(oo_.dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
skipline()
disp(['[', namendo,' vs lagged ',namlagendo,']'])
if ~isempty(ilagendo)
%y0=squeeze(T(iendo,ilagendo,istable));
y0=squeeze(T(iendo,ilagendo,:));
if (max(y0)-min(y0))>1.e-10
if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph
ifig=ifig+1;
hfig = dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping: ' namendo,' vs lags ',int2str(ifig)]);
hh_fig = dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping: ' namendo,' vs lags ',int2str(ifig)]);
iplo=0;
end
iplo=iplo+1;
......@@ -331,9 +331,9 @@ for j = 1:length(anamendo)
if isempty(dir(xdir0))
mkdir(xdir0)
end
atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs ', namlagendo];
aname=[type '_' namendo '_vs_' namlagendo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo];
atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs lagged', namlagendo];
aname=[type '_' namendo '_vs_lag_' namlagendo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs lagged',namlagendo];
options_map.amap_name = aname;
options_map.amap_title = atitle;
options_map.figtitle = atitle0;
......@@ -349,20 +349,23 @@ for j = 1:length(anamendo)
end
if ~options_.nograph
hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs lagged ', namlagendo]);
hc = cumplot(y0);
hc = gsa.cumplot(y0);
a=axis; delete(hc);
% hist(mat_moment{ij}),
x1val=max(threshold(1),a(1));
x2val=min(threshold(2),a(2));
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
set(hp,'FaceColor', [0.7 0.8 1])
hold all,
hc = cumplot(y0);
hc = gsa.cumplot(y0);
set(hc,'color','k','linewidth',2)
hold off,
hold off
if options_.TeX
title([namendo_tex,' vs lagged ', namlagendo_tex ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','latex')
else
title([namendo,' vs lagged ', namlagendo ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','none')
end
dyn_saveas(hf,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],['Reduced Form Mapping (Monte Carlo Filtering): ',strrep(namendo,'_','\_'),' vs lagged ', strrep(namlagendo,'_','\_')],[type '_' namendo,'_vs_', namlagendo],1)
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],['Reduced Form Mapping (Monte Carlo Filtering): ',namendo_tex,' vs lagged ', namlagendo_tex],[type '_' namendo,'_vs_', namlagendo],1)
end
delete([xdir, '/*threshold*.*'])
......@@ -374,24 +377,28 @@ for j = 1:length(anamendo)
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'inside threshold';
options_mcf.nobeha_title = 'outside threshold';
if options_.TeX
options_mcf.beha_title_latex = 'inside threshold';
options_mcf.nobeha_title_latex = 'outside threshold';
end
options_mcf.title = atitle0;
options_mcf.OutputDirectoryName = xdir;
if ~isempty(iy) && ~isempty(iyc)
fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_);
icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
lpmat=x0(iy,:);
if nshocks
lpmat0=xx0(iy,:);
end
istable=[1:length(iy)];
istable=1:length(iy);
save([xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo '_threshold' ],'lpmat','lpmat0','istable','y0','x0','xx0','iy','iyc')
lpmat=[]; lpmat0=[]; istable=[];
if length(iy)<=10 || length(iyc)<=10,
lpmat0=[];
if length(iy)<=10 || length(iyc)<=10
icheck = []; % do the generic plot in any case
end
else
icheck = [];
end
......@@ -412,11 +419,10 @@ for j = 1:length(anamendo)
end
atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namlagendo];
options_mcf.title = atitle0;
indmcf = redform_mcf(y0, x0, options_mcf, options_);
redform_mcf(y0, x0, options_mcf, options_, M_.fname, bayestopt_.name, estim_params_);
end
end
else
[yy, xdir] = log_trans_(y0,xdir0);
atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namlagendo];
aname=[type '_' namendo '_vs_' namlagendo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo];
......@@ -429,27 +435,30 @@ for j = 1:length(anamendo)
end
if isempty(threshold) && ~options_.nograph
figure(hfig),
figure(hh_fig),
subplot(3,3,iplo),
if ilog
[saso, iso] = sort(-silog(:,js));
[~, iso] = sort(-silog(:,js));
bar([silog(iso(1:min(np,10)),js)])
logflag='log';
else
[saso, iso] = sort(-si(:,js));
[~, iso] = sort(-si(:,js));
bar(si(iso(1:min(np,10)),js))
logflag='';
end
%set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
if options_.TeX
text(ip,-0.02,deblank(pnames_tex(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
title([logflag,' ',namendo,' vs ',namlagendo,'(-1)'],'interpreter','none')
if iplo==9
dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namlagendo,'_','\_'),'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],1)
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',namendo_tex,' vs ',namlagendo_tex,'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],1)
end
end
......@@ -459,49 +468,38 @@ for j = 1:length(anamendo)
end
end
if iplo<9 && iplo>0 && ifig && ~options_.nograph
dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namlagendo,'_','\_'),'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],options_.figures.textwidth*min(iplo/3,1));
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',namendo_tex,' vs ',namlagendo_tex,'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],options_.figures.textwidth*min(iplo/3,1));
end
end
if isempty(threshold) && ~options_.nograph
hh_fig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA');
if ilog==0
hfig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(si)
% boxplot(si','whis',10,'symbol','r.')
myboxplot(si',[],'.',[],10)
gsa.boxplot(si',[],'.',[],10)
else
gsa.boxplot(silog',[],'.',[],10)
end
xlabel(' ')
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np)
set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
if ilog==0
title('Reduced form GSA')
dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_gsa'],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_gsa'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa'],'Reduced Form GSA','redform_gsa')
else
hfig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(silog)
% boxplot(silog','whis',10,'symbol','r.')
myboxplot(silog',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
xlabel(' ')
set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Reduced form GSA - Log-transformed elements')
dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_gsa_log'],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_gsa_log'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa_log'],'Reduced form GSA - Log-transformed elements','redform_gsa_log')
end
end
function si = redform_private(x0, y0, options_map, options_)
np=size(x0,2);
x00=x0;
ilog = options_map.log_trans;
......@@ -515,7 +513,7 @@ if options_map.prior_range
x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3));
end
else
x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
x0=gsa.priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
end
if ilog
......@@ -531,14 +529,12 @@ if iload==0
nest=max(50,nrun/2);
nest=min(250,nest);
nfit=min(1000,nrun);
% dotheplots = (nfit<=nest);
% gsa_ = gsa_sdp(y0(1:nest), x0(1:nest,:), 2, [],[-1 -1 -1 -1 -1 0],[],0,[fname,'_est'], pnames);
[ys,is] = sort(y0);
[~,is] = sort(y0);
istep = ceil(nrun/nest);
if istep>1
iest = is(floor(istep/2):istep:end);
nest = length(iest);
irest = is(setdiff([1:nrun],[floor(istep/2):istep:nrun]));
irest = is(setdiff(1:nrun,floor(istep/2):istep:nrun));
istep = ceil(length(irest)/(nfit-nest));
ifit = union(iest, irest(1:istep:end));
else
......@@ -550,29 +546,34 @@ if iload==0
ifit = union(ifit, irest(end));
end
nfit=length(ifit);
% ifit = union(iest, irest(randperm(nrun-nest,nfit-nest)));
% ifit = iest;
% nfit=nest;
ipred = setdiff([1:nrun],ifit);
ipred = setdiff(1:nrun,ifit);
if ilog
[y1, tmp, isig, lam] = log_trans_(y0(iest));
[~, ~, isig, lam] = gsa.log_transform(y0(iest));
y1 = log(y0*isig+lam);
end
if ~options_.nograph
hfig=dyn_figure(options_.nodisplay,'name',options_map.figtitle);
hh_fig=dyn_figure(options_.nodisplay,'name',options_map.figtitle);
subplot(221)
if ilog
if isoctave
hist(y1,30)
else
histogram(y1,30)
end
else
if isoctave
hist(y0,30)
else
histogram(y0,30)
end
end
title(options_map.title,'interpreter','none')
subplot(222)
if ilog
hc = cumplot(y1);
hc = gsa.cumplot(y1);
else
hc = cumplot(y0);
hc = gsa.cumplot(y0);
end
set(hc,'color','k','linewidth',2)
title([options_map.title ' CDF'],'interpreter','none')
......@@ -582,15 +583,7 @@ if iload==0
if ilog
[gsa22, gsa1, gsax] = ss_anova_log(y1(iest), x0(iest,:), isig, lam, gsa0);
end
% if (gsa1.out.bic-gsa0.out.bic) < 10,
% y00=y0;
% gsa00=gsa0;
% gsa0=gsa1;
% y0=y1;
% ilog=1;
% end
if nfit>nest
% gsa_ = gsa_sdp(y0(1:nfit), x0(1:nfit,:), -2, gsa_.nvr*nest^3/nfit^3,[-1 -1 -1 -1 -1 0],[],0,fname, pnames);
nvr = gsa0.nvr*nest^3/nfit^3;
nvr(gsa0.stat<2) = gsa0.nvr(gsa0.stat<2)*nest^5/nfit^5;
gsa_ = ss_anova(y0(ifit), x0(ifit,:), 1, 0, 2, nvr);
......@@ -601,33 +594,6 @@ if iload==0
nvrx = gsax.nvr*nest^3/nfit^3;
nvrx(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5;
[gsa22, gsa1, gsax] = ss_anova_log(y1(ifit), x0(ifit,:), isig, lam, gsa0, [nvr1' nvrx']);
% gsa1 = ss_anova(y1(ifit), x0(ifit,:), 1, 0, 2, nvr);
% gsa2=gsa1;
% gsa2.y = gsa0.y;
% gsa2.fit = (exp(gsa1.fit)-lam)*isig;
% gsa2.f0 = mean(gsa2.fit);
% gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
% gsa2.out.bic = gsa2.out.bic-nest*log(gsa1.out.SSE)+nest*log(gsa2.out.SSE);
% gsa2.r2 = 1-cov(gsa2.fit-gsa2.y)/cov(gsa2.y);
% for j=1:np,
% gsa2.fs(:,j) = exp(gsa1.fs(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
% gsa2.f(:,j) = exp(gsa1.f(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
% gsa2.si(j) = var(gsa2.f(:,j))/var(gsa2.y);
% end
% nvr = gsax.nvr*nest^3/nfit^3;
% nvr(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5;
% gsax = ss_anova([gsa2.y-gsa2.fit], x0(ifit,:), 1, 0, 2, nvr);
% gsa22=gsa2;
% gsa22.fit = gsa2.fit+gsax.fit;
% gsa22.f0 = mean(gsa22.fit);
% gsa22.out.SSE = sum((gsa22.fit-gsa22.y).^2);
% gsa22.out.bic = nest*log(gsa22.out.SSE/nest) + (gsax.out.df+gsa2.out.df-1)*log(nest);
% gsa22.r2 = 1-sum((gsa22.fit-gsa22.y).^2)/sum((gsa22.y-mean(gsa22.y)).^2);
% for j=1:np,
% gsa22.fs(:,j) = gsa2.fs(:,j)+gsax.fs(:,j);
% gsa22.f(:,j) = gsa2.f(:,j)+gsax.f(:,j);
% gsa22.si(j) = var(gsa22.f(:,j))/var(gsa22.y);
% end
gsa_ = gsa22;
end
else
......@@ -638,31 +604,23 @@ if iload==0
end
end
save([fname,'_map.mat'],'gsa_')
[sidum, iii]=sort(-gsa_.si);
[~, iii]=sort(-gsa_.si);
gsa_.x0=x00(ifit,:);
if ~options_.nograph
hmap=gsa_sdp_plot(gsa_,[fname '_map'],pnames,iii(1:min(12,np)));
set(hmap,'name',options_map.amap_title);
end
gsa_.x0=x0(ifit,:);
% copyfile([fname,'_est.mat'],[fname,'.mat'])
if ~options_.nograph
figure(hfig);
figure(hh_fig);
subplot(223),
plot(y0(ifit),[gsa_.fit y0(ifit)],'.'),
r2 = gsa_.r2;
% if ilog,
% plot(y00(ifit),[log_trans_(gsa_.fit,'',isig,lam) y00(ifit)],'.'),
% r2 = 1 - cov(log_trans_(gsa_.fit,'',isig,lam)-y00(ifit))/cov(y00(ifit));
% else
% plot(y0(ifit),[gsa_.fit y0(ifit)],'.'),
% r2 = gsa_.r2;
% end
title(['Learning sample fit - R2=' num2str(r2,2)],'interpreter','none')
if nfit<nrun
if ilog
yf = ss_anova_fcast(x0(ipred,:), gsa1);
yf = log_trans_(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
yf = gsa.log_transform(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
else
yf = ss_anova_fcast(x0(ipred,:), gsa_);
end
......@@ -672,7 +630,7 @@ if iload==0
plot(yn,[yf yn],'.'),
title(['Out-of-sample prediction - R2=' num2str(r2,2)],'interpreter','none')
end
dyn_saveas(hfig,fname,options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,fname,options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,fname,['Out-of-sample prediction - R2=' num2str(r2,2)],'redform_gsa_log')
if options_.nodisplay
......@@ -680,20 +638,16 @@ if iload==0
end
end
else
% gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames);
% gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames);
load([fname,'_map.mat'],'gsa_')
if ~options_.nograph
yf = ss_anova_fcast(x0, gsa_);
hfig=dyn_figure(options_.nodisplay,'name',options_map.title);
hh_fig=dyn_figure(options_.nodisplay,'name',options_map.title);
plot(y0,[yf y0],'.'),
title([namy,' vs ', namx,' pred'],'interpreter','none')
dyn_saveas(hfig,[fname '_pred'],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[fname '_pred'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[fname '_pred'],options_map.title,[namy,' vs ', namx,' pred'])
end
end
% si = gsa_.multivariate.si;
si = gsa_.si;
return
......@@ -703,7 +657,7 @@ function gsa2 = log2level_map(gsa1, isig, lam)
nest=length(gsa1.y);
np = size(gsa1.x0,2);
gsa2=gsa1;
gsa2.y = log_trans_(gsa1.y,'',isig,lam);
gsa2.y = gsa.log_transform(gsa1.y,'',isig,lam);
gsa2.fit = (exp(gsa1.fit)-lam)*isig;
gsa2.f0 = mean(gsa2.fit);
gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
......@@ -730,6 +684,8 @@ else
end
gsa2 = log2level_map(gsa1, isig, lam);
if nargin >=5 && ~isempty(gsa0)
nvr2=NaN(np,1);
nvr0=NaN(np,1);
for j=1:np
nvr2(j) = var(diff(gsa2.fs(:,j),2));
nvr0(j) = var(diff(gsa0.fs(:,j),2));
......@@ -760,26 +716,24 @@ end
return
function indmcf = redform_mcf(y0, x0, options_mcf, options_)
function indmcf = redform_mcf(y0, x0, options_mcf, options_, fname, parnames, estim_params_)
hfig=dyn_figure(options_.nodisplay,'name',options_mcf.amcf_title);
hh_fig=dyn_figure(options_.nodisplay,'name',options_mcf.amcf_title);
[post_mean, post_median, post_var, hpd_interval, post_deciles, ...
density] = posterior_moments(y0,1,0.9);
[~, ~, ~, ~, post_deciles] = posterior_moments(y0,0.9);
post_deciles = [-inf; post_deciles; inf];
for jt=1:10
indy{jt}=find( (y0>post_deciles(jt)) & (y0<=post_deciles(jt+1)));
leg{jt}=[int2str(jt) '-dec'];
end
[proba, dproba] = stab_map_1(x0, indy{1}, indy{end}, [],0);
[proba] = gsa.stability_mapping_univariate(x0, indy{1}, indy{end}, [], fname, options_, parnames, estim_params_,0);
indmcf=find(proba<options_mcf.pvalue_ks);
if isempty(indmcf)
[tmp,jtmp] = sort(proba,2,'ascend');
[~,jtmp] = sort(proba,1,'ascend');
indmcf = jtmp(1);
% indmcf = jtmp(1:min(2,length(proba)));
end
[tmp,jtmp] = sort(proba(indmcf),2,'ascend');
[~,jtmp] = sort(proba(indmcf),1,'ascend');
indmcf = indmcf(jtmp);
nbr_par = length(indmcf);
nrow=ceil(sqrt(nbr_par+1));
......@@ -793,12 +747,16 @@ for jx=1:nbr_par
subplot(nrow,ncol,jx)
hold off
for jt=1:10
h=cumplot(x0(indy{jt},indmcf(jx)));
h=gsa.cumplot(x0(indy{jt},indmcf(jx)));
set(h,'color', cmap(jt,:), 'linewidth', 2)
hold all
end
if options_.TeX
title(options_mcf.param_names_tex(indmcf(jx),:),'interpreter','latex')
else
title(options_mcf.param_names(indmcf(jx),:),'interpreter','none')
end
end
hleg = legend(leg);
aa=get(hleg,'Position');
aa(1)=1-aa(3)-0.02;
......@@ -813,7 +771,7 @@ if ~isoctave
'horizontalalignment','center');
end
dyn_saveas(hfig,[options_mcf.OutputDirectoryName filesep options_mcf.fname_,'_',options_mcf.amcf_name],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[options_mcf.OutputDirectoryName filesep options_mcf.fname_,'_',options_mcf.amcf_name],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[options_mcf.OutputDirectoryName filesep options_mcf.fname_,'_',options_mcf.amcf_name],strrep(options_mcf.amcf_title,'_','\_'),[options_mcf.fname_,'_',options_mcf.amcf_name])
return
......@@ -824,7 +782,7 @@ if nargin<5
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([figpath '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by redform_map.m (Dynare).\n');
fprintf(fidTeX,'%% TeX eps-loader file generated by reduced_form_mapping.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
......
matlab/gsa/redform_screen.m → matlab/+gsa/reduced_form_screening.m
View file @ 4a32acc2
function redform_screen(dirname, options_gsa_)
%function redform_map(dirname, options_gsa_)
% inputs (from opt_gsa structure
% anamendo = options_gsa_.namendo;
% anamlagendo = options_gsa_.namlagendo;
% anamexo = options_gsa_.namexo;
% iload = options_gsa_.load_redform;
function reduced_form_screening(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
% reduced_form_screening(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
% Conduct reduced form screening
% Inputs:
% - dirname [string] name of the output directory
% - options_gsa_ [structure] GSA options_
% - estim_params [structure] describing the estimated parameters
% - M_ [structure] describing the model
% - dr [structure] decision rules
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright (C) 2012-2016 European Commission
% Copyright (C) 2012-2018 Dynare Team
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -26,19 +30,23 @@ function redform_screen(dirname, options_gsa_)
% 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 <http://www.gnu.org/licenses/>.
global M_ oo_ estim_params_ options_ bayestopt_
% options_gsa_ = options_.opt_gsa;
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
anamexo = options_gsa_.namexo;
iload = options_gsa_.load_redform;
anamendo_tex = options_gsa_.namendo_tex;
anamlagendo_tex = options_gsa_.namlagendo_tex;
anamexo_tex = options_gsa_.namexo_tex;
nliv = options_gsa_.morris_nliv;
pnames = M_.param_names(estim_params_.param_vals(:,1));
if options_.TeX
for par_iter=1:size(estim_params_.param_vals(:,1),1)
[~,tex_names{par_iter,1}]=get_the_name(estim_params_.param_vals(par_iter,1),options_.TeX, M_, estim_params_, options_.varobs);
end
end
if nargin==0
dirname='';
end
......@@ -54,37 +62,47 @@ nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
js=0;
for j=1:size(anamendo,1)
namendo = deblank(anamendo(j,:));
iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact');
namendo = anamendo{j,:};
namendo_tex = anamendo_tex{j,:};
iendo = strmatch(namendo, M_.endo_names(dr.order_var), 'exact');
iplo=0;
ifig=0;
for jx=1:size(anamexo,1)
namexo = deblank(anamexo(jx,:));
namexo = anamexo{jx};
namexo_tex = anamexo_tex{jx};
iexo = strmatch(namexo, M_.exo_names, 'exact');
if ~isempty(iexo)
y0=teff(T(iendo,iexo+nspred,:), kn, istable);
y0=gsa.teff(T(iendo,iexo+nspred,:), kn, istable);
if ~isempty(y0)
if mod(iplo,9)==0
ifig = ifig+1;
hh = dyn_figure(options_.nodisplay, 'name', [namendo,[' vs. shocks '], int2str(ifig)]);
hh_fig = dyn_figure(options_.nodisplay, 'name', [namendo,' vs. shocks ', int2str(ifig)]);
iplo = 0;
end
iplo = iplo+1;
js = js+1;
subplot(3, 3, iplo)
[SAmeas, SAMorris] = Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0, nliv);
[~, SAMorris] = gsa.Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0, nliv);
SAM = squeeze(SAMorris(nshock+1:end,1));
SA(:,js) = SAM./(max(SAM)+eps);
[saso, iso] = sort(-SA(:,js));
[~, iso] = sort(-SA(:,js));
bar(SA(iso(1:min(np,10)),js))
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
if options_.TeX
text(ip,-0.02,tex_names(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,pnames(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if options_.TeX
title([namendo_tex,' vs. ',namexo_tex],'interpreter','latex')
else
title([namendo,' vs. ',namexo],'interpreter','none')
end
if iplo==9
dyn_saveas(hh,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)],ifig,[namendo,' vs. shocks ',int2str(ifig)],[namendo,'_vs_shock'],1)
end
......@@ -92,68 +110,78 @@ for j=1:size(anamendo,1)
end
end
if iplo<9 && iplo>0 && ifig
dyn_saveas(hh,[dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)],ifig,[namendo,' vs. shocks ',int2str(ifig)],[namendo,'_vs_shock'],options_.figures.textwidth*min(iplo/3))
end
iplo=0;
ifig=0;
for je=1:size(anamlagendo,1)
namlagendo=deblank(anamlagendo(je,:));
ilagendo=strmatch(namlagendo, M_.endo_names(oo_.dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
namlagendo=anamlagendo{je};
namlagendo_tex=anamlagendo_tex{je};
ilagendo=strmatch(namlagendo, M_.endo_names(dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
if ~isempty(ilagendo)
y0=teff(T(iendo,ilagendo,:),kn,istable);
y0=gsa.teff(T(iendo,ilagendo,:),kn,istable);
if ~isempty(y0)
if mod(iplo,9)==0
ifig=ifig+1;
hh=dyn_figure(options_.nodisplay,'name',[namendo,' vs. lagged endogenous ',int2str(ifig)]);
hh_fig=dyn_figure(options_.nodisplay,'name',[namendo,' vs. lagged endogenous ',int2str(ifig)]);
iplo=0;
end
iplo=iplo+1;
js=js+1;
subplot(3,3,iplo),
[SAmeas, SAMorris] = Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0,nliv);
[~, SAMorris] = gsa.Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0,nliv);
SAM = squeeze(SAMorris(nshock+1:end,1));
SA(:,js)=SAM./(max(SAM)+eps);
[saso, iso] = sort(-SA(:,js));
[~, iso] = sort(-SA(:,js));
bar(SA(iso(1:min(np,10)),js))
%set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
if options_.TeX
text(ip,-0.02,tex_names(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,pnames{iso(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if options_.TeX
title([namendo_tex,' vs. ',namlagendo_tex,'(-1)'],'interpreter','latex')
else
title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
end
if iplo==9
dyn_saveas(hh,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],ifig,[namendo,' vs. lagged endogenous ',int2str(ifig)],[namendo,'_vs_lags'],1)
end
end
end
end
if iplo<9 && iplo>0 && ifig
dyn_saveas(hh,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],ifig,[namendo,' vs. lagged endogenous ',int2str(ifig)],[namendo,'_vs_lags'],options_.figures.textwidth*min(iplo/3))
end
end
hh=dyn_figure(options_.nodisplay,'Name','Reduced form screening');
%bar(SA)
% boxplot(SA','whis',10,'symbol','r.')
myboxplot(SA',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
hh_fig=dyn_figure(options_.nodisplay,'Name','Reduced form screening');
gsa.boxplot(SA',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np)
set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
if options_.TeX
text(ip,-0.02,tex_names(ip),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,pnames{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
ylabel('Elementary Effects')
title('Reduced form screening')
dyn_saveas(hh,[dirname,'/',M_.fname,'_redform_screen'],options_.nodisplay,options_.graph_format);
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_redform_screen'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_redform_screen'],1,'Reduced form screening','redform_screen',1)
......@@ -163,7 +191,7 @@ if nargin<6
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([figpath '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by redform_screen.m (Dynare).\n');
fprintf(fidTeX,'%% TeX eps-loader file generated by reduced_form_screening.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
......
matlab/dynare_sensitivity.m → matlab/+gsa/run.m
View file @ 4a32acc2
function x0=dynare_sensitivity(options_gsa)
function x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
% x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
% Frontend to the Sensitivity Analysis Toolbox for DYNARE
% Inputs:
% - 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
% - options_gsa [structure] MATLAB's structure describing the GSA options
%
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% M. Ratto (2008), Analysing DSGE Models with Global Sensitivity Analysis,
% Computational Economics (2008), 31, pp. 115–139
% Copyright (C) 2008-2018 Dynare Team
% Copyright © 2008-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -19,16 +28,13 @@ function x0=dynare_sensitivity(options_gsa)
% 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 <http://www.gnu.org/licenses/>.
global M_ options_ oo_ bayestopt_ estim_params_
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if options_.dsge_var
error('Identification does not support DSGE-VARs at the current stage')
end
fname_ = M_.fname;
lgy_ = M_.endo_names;
x0=[];
% check user defined options
......@@ -43,7 +49,6 @@ end
if isfield(options_gsa,'morris') && options_gsa.morris==1
if isfield(options_gsa,'identification') && options_gsa.identification==0
% options_gsa.redform=1;
end
if isfield(options_gsa,'ppost') && options_gsa.ppost
error('sensitivity:: Morris is incompatible with posterior sampling')
......@@ -84,7 +89,10 @@ elseif isfield(options_gsa,'neighborhood_width') && options_gsa.neighborhood_wid
options_.mode_file='';
end
if options_.order~=1
warning('dynare_sensitivity: dynare_sensitivity does only support order=1, resetting to order=1.')
options_.order = 1;
end
if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
if isempty(options_gsa.datafile) && options_gsa.rmse
......@@ -93,6 +101,9 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
disp('must be specified for RMSE analysis!');
error('Sensitivity anaysis error!')
end
if isfield(options_gsa,'nobs')
options_.nobs=options_gsa.nobs;
end
if ~isempty(options_.nobs) && length(options_.nobs)~=1
error('dynare_sensitivity does not support recursive estimation. Please specify nobs as a scalar, not a vector.')
end
......@@ -100,9 +111,6 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
if isfield(options_gsa,'first_obs')
options_.first_obs=options_gsa.first_obs;
end
if isfield(options_gsa,'nobs')
options_.nobs=options_gsa.nobs;
end
if isfield(options_gsa,'presample')
options_.presample=options_gsa.presample;
end
......@@ -124,15 +132,28 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
options_.mode_compute = 0;
options_.filtered_vars = 1;
options_.plot_priors = 0;
[dataset_,dataset_info,xparam1,hh, M_, options_, oo_, estim_params_, bayestopt_] = ...
[dataset_,dataset_info,~,~, M_, options_, oo_, estim_params_, bayestopt_] = ...
dynare_estimation_init(M_.endo_names, fname_, 1, M_, options_, oo_, estim_params_, bayestopt_);
% computes a first linear solution to set up various variables
else
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
if options_.prior_trunc==0
options_.prior_trunc=1.e-10;
end
end
if M_.exo_nbr==0
error('dynare_sensitivity does not support having no varexo in the model. As a workaround you could define a dummy exogenous variable.')
end
[~,~,~,~,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if isfield(oo_.dr,'eigval') && any(abs(oo_.dr.eigval-1)<abs(1-options_.qz_criterium)) && options_.qz_criterium<1
fprintf('\ngsa: The model features a unit root, but qz_criterium<1. Check whether that is intended.')
fprintf('\ngsa: If not, use the diffuse_filter option.\n')
end
[make,my,day,punk,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
options_gsa = set_default_option(options_gsa,'identification',0);
if options_gsa.identification
......@@ -152,7 +173,6 @@ if options_gsa.identification
options_.options_ident.load_ident_files = options_gsa.load_ident_files;
options_.options_ident.useautocorr = options_gsa.useautocorr;
options_.options_ident.ar = options_gsa.ar;
options_ident=options_.options_ident;
else
options_ident=[];
options_ident = set_default_option(options_ident,'load_ident_files',options_gsa.load_ident_files);
......@@ -179,7 +199,6 @@ options_gsa = set_default_option(options_gsa,'load_stab',0);
options_gsa = set_default_option(options_gsa,'alpha2_stab',0);
options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
options_gsa = set_default_option(options_gsa,'pvalue_corr',1.e-5);
%options_gsa = set_default_option(options_gsa,'load_mh',0);
% REDFORM mapping
options_gsa = set_default_option(options_gsa,'redform',0);
options_gsa = set_default_option(options_gsa,'load_redform',0);
......@@ -188,8 +207,29 @@ options_gsa = set_default_option(options_gsa,'threshold_redform',[]);
options_gsa = set_default_option(options_gsa,'ksstat_redform',0.001);
options_gsa = set_default_option(options_gsa,'alpha2_redform',1.e-5);
options_gsa = set_default_option(options_gsa,'namendo',{});
options_gsa = set_default_option(options_gsa,'namlagendo',[]);
options_gsa = set_default_option(options_gsa,'namlagendo',{});
options_gsa = set_default_option(options_gsa,'namexo',{});
options_gsa = set_default_option(options_gsa,'namendo_tex',{});
options_gsa = set_default_option(options_gsa,'namlagendo_tex',{});
options_gsa = set_default_option(options_gsa,'namexo_tex',{});
if strmatch(':',options_gsa.namendo,'exact')
options_gsa.namendo = M_.endo_names(1:M_.orig_endo_nbr);
end
if strmatch(':',options_gsa.namexo,'exact')
options_gsa.namexo = M_.exo_names;
end
if strmatch(':',options_gsa.namlagendo,'exact')
options_gsa.namlagendo = M_.endo_names(1:M_.orig_endo_nbr);
end
if options_.TeX
[~,Locb]=ismember(options_gsa.namendo,M_.endo_names);
options_gsa.namendo_tex=cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(Locb), 'UniformOutput', false);
[~,Locb]=ismember(options_gsa.namlagendo,M_.endo_names);
options_gsa.namlagendo_tex=cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(Locb), 'UniformOutput', false);
[~,Locb]=ismember(options_gsa.namexo,M_.exo_names);
options_gsa.namexo_tex=cellfun(@(x) horzcat('$', x, '$'), M_.exo_names_tex(Locb), 'UniformOutput', false);
end
% RMSE mapping
options_gsa = set_default_option(options_gsa,'load_rmse',0);
options_gsa = set_default_option(options_gsa,'lik_only',0);
......@@ -198,8 +238,9 @@ options_gsa = set_default_option(options_gsa,'var_rmse', options_.varobs);
options_gsa.var_rmse_tex={};
for ii=1:length(options_gsa.var_rmse)
temp_name = M_.endo_names_tex{strmatch(options_gsa.var_rmse{ii}, M_.endo_names, 'exact')};
options_gsa.var_rmse_tex = vertcat(options_gsa.var_rmse_tex, temp_name);
options_gsa.var_rmse_tex = vertcat(options_gsa.var_rmse_tex, ['$' temp_name '$']);
end
options_gsa.varobs_tex = cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(options_.varobs_id), 'UniformOutput', false);
options_gsa = set_default_option(options_gsa,'pfilt_rmse', 0.1);
options_gsa = set_default_option(options_gsa,'istart_rmse', options_.presample+1);
options_gsa = set_default_option(options_gsa,'alpha_rmse', 0.001);
......@@ -231,31 +272,22 @@ if options_gsa.morris==1
options_gsa.pprior=1;
end
options_gsa.ppost=0;
%options_gsa.stab=1;
options_gsa.glue=0;
options_gsa.rmse=0;
options_gsa.load_rmse=0;
options_gsa.alpha2_stab=1;
options_gsa.pvalue_ks=0;
options_gsa.pvalue_corr=0;
% if options_gsa.morris==3,
% options_gsa = set_default_option(options_gsa,'Nsam',256);
% OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
% else
OutputDirectoryName = CheckPath('gsa/screen',M_.dname);
% end
else
OutputDirectoryName = CheckPath('gsa',M_.dname);
end
% options_.opt_gsa = options_gsa;
if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform) && options_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior.mat'];
if ~exist(filetoload)
if ~exist(filetoload,'file')
disp([filetoload,' not found!'])
disp(['You asked to load a non existent analysis'])
%options_gsa.load_stab=0;
disp('You asked to load a non existent analysis')
return
else
if isempty(strmatch('bkpprior',who('-file', filetoload),'exact'))
......@@ -279,7 +311,7 @@ if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform)
end
if options_gsa.stab && ~options_gsa.ppost
x0 = stab_map_(OutputDirectoryName,options_gsa);
x0 = gsa.stability_mapping(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
if isempty(x0)
skipline()
disp('Sensitivity computations stopped: no parameter set provided a unique solution')
......@@ -287,16 +319,13 @@ if options_gsa.stab && ~options_gsa.ppost
end
end
% reduced form
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
options_.opt_gsa = options_gsa;
if ~isempty(options_gsa.moment_calibration) || ~isempty(options_gsa.irf_calibration)
map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_,bayestopt_);
gsa.map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_,bayestopt_);
end
if options_gsa.identification
map_ident_(OutputDirectoryName,options_gsa);
gsa.map_identification(OutputDirectoryName,options_gsa,M_,oo_,options_,estim_params_,bayestopt_);
end
if options_gsa.redform && ~isempty(options_gsa.namendo)
......@@ -304,7 +333,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
filnam = dir([M_.dname filesep 'metropolis' filesep '*param_irf*.mat']);
lpmat=[];
for j=1:length(filnam)
load ([M_.dname filesep 'metropolis' filesep M_.fname '_param_irf' int2str(j) '.mat'])
load ([M_.dname filesep 'metropolis' filesep M_.fname '_param_irf' int2str(j) '.mat'],'stock')
lpmat=[lpmat; stock];
end
clear stock
......@@ -322,35 +351,24 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
save([OutputDirectoryName filesep M_.fname '_mc.mat'],'lpmat','lpmat0','istable','iunstable','iwrong','iindeterm')
options_gsa.load_stab=1;
x0 = stab_map_(OutputDirectoryName,options_gsa);
end
if strmatch(':',options_gsa.namendo,'exact')
options_gsa.namendo = M_.endo_names(1:M_.orig_endo_nbr);
end
if strmatch(':',options_gsa.namexo,'exact')
options_gsa.namexo = M_.exo_names;
end
if strmatch(':',options_gsa.namlagendo,'exact')
options_gsa.namlagendo = M_.endo_names(1:M_.orig_endo_nbr);
x0 = gsa.stability_mapping(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
end
% options_.opt_gsa = options_gsa;
if options_gsa.morris==1
redform_screen(OutputDirectoryName,options_gsa);
gsa.reduced_form_screening(OutputDirectoryName,options_gsa, estim_params_, M_, oo_.dr, options_, bayestopt_);
else
% check existence of the SS_ANOVA toolbox
if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
fprintf('\nThe "SS-ANOVA-R: MATLAB Toolbox for the estimation of Smoothing Spline ANOVA models with Recursive algorithms" is missing.\n')
fprintf('To obtain it, go to:\n\n')
fprintf('http://ipsc.jrc.ec.europa.eu/?id=790 \n\n')
fprintf('https://joint-research-centre.ec.europa.eu/system/files/2025-01/ss_anova_recurs.zip \n\n')
fprintf('and follow the instructions there.\n')
fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
fprintf('After obtaining the files, you need to unpack them and set a MATLAB Path to those files.\n')
error('SS-ANOVA-R Toolbox missing!')
end
redform_map(OutputDirectoryName,options_gsa);
gsa.reduced_form_mapping(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_);
end
end
% RMSE mapping
% function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
options_.opt_gsa = options_gsa;
if options_gsa.rmse
if ~options_gsa.ppost
......@@ -377,7 +395,6 @@ if options_gsa.rmse
options_.forecast=0;
options_.filtered_vars=0;
end
% dynare_MC([],OutputDirectoryName,data,rawdata,data_info);
if options_gsa.pprior
TmpDirectoryName = ([M_.dname filesep 'gsa' filesep 'prior']);
else
......@@ -394,36 +411,19 @@ if options_gsa.rmse
delete([TmpDirectoryName filesep filparam(j).name]);
end
end
end
prior_posterior_statistics('gsa',dataset_, dataset_info);
oo_=prior_posterior_statistics('gsa',dataset_, dataset_info,M_,oo_,options_,estim_params_,bayestopt_,'gsa::mcmc');
if options_.bayesian_irf
PosteriorIRF('gsa');
oo_=PosteriorIRF('gsa',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,'gsa::mcmc');
end
options_gsa.load_rmse=0;
% else
% if options_gsa.load_rmse==0,
% disp('You already saved a MC filter/smoother analysis ')
% disp('Do you want to overwrite ?')
% pause;
% if options_gsa.pprior
% delete([OutputDirectoryName,'/',fname_,'_prior_*.mat'])
% else
% delete([OutputDirectoryName,'/',fname_,'_mc_*.mat'])
% end
% dynare_MC([],OutputDirectoryName);
% options_gsa.load_rmse=0;
% end
end
end
clear a;
% filt_mc_(OutputDirectoryName,data_info);
filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info);
gsa.monte_carlo_filtering(OutputDirectoryName,options_gsa,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_);
end
options_.opt_gsa = options_gsa;
if options_gsa.glue
dr_ = oo_.dr;
if options_gsa.ppost
......@@ -436,12 +436,11 @@ if options_gsa.glue
load([OutputDirectoryName,'/',fname_,'_mc']);
end
end
if ~exist('x')
disp(['No RMSE analysis is available for current options'])
disp(['No GLUE file prepared'])
if ~exist('x','var')
disp('No RMSE analysis is available for current options')
disp('No GLUE file prepared')
return,
end
nruns=size(x,1);
gend = options_.nobs;
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
......@@ -449,28 +448,20 @@ if options_gsa.glue
rawdata = log(rawdata);
end
if options_.prefilter == 1
%data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
data = transpose(rawdata-ones(gend,1)*mean(rawdata,1));
else
data = transpose(rawdata);
end
Obs.data = data;
Obs.time = [1:gend];
Obs.time = 1:gend;
Obs.num = gend;
for j=1:length(options_.varobs)
Obs.name{j} = options_.varobs{j};
vj = options_.varobs{j};
jxj = strmatch(vj,lgy_(dr_.order_var),'exact');
js = strmatch(vj,lgy_,'exact');
jxj = strmatch(vj,M_.endo_names(dr_.order_var),'exact');
if ~options_gsa.ppost
% y0=zeros(gend+1,nruns);
% nb = size(stock_filter,3);
% y0 = squeeze(stock_filter(:,jxj,:)) + ...
% kron(stock_ys(js,:),ones(size(stock_filter,1),1));
% Out(j).data = y0';
% Out(j).time = [1:size(y0,1)];
Out(j).data = jxj;
Out(j).time = [pwd,'/',OutputDirectoryName];
else
......@@ -484,17 +475,7 @@ if options_gsa.glue
Lik(j).isam = 1;
Lik(j).data = rmse_MC(:,j)';
if ~options_gsa.ppost
% y0 = squeeze(stock_smooth(:,jxj,:)) + ...
% kron(stock_ys(js,:),ones(size(stock_smooth,1),1));
% Out1(j).name = vj;
% Out1(j).ini = 'yes';
% Out1(j).time = [1:size(y0,1)];
% Out1(j).data = y0';
Out1=Out;
else
Out1=Out;
end
ismoo(j)=jxj;
end
......@@ -504,10 +485,6 @@ if options_gsa.glue
jsmoo=jsmoo+1;
vj = M_.endo_names{dr_.order_var(j)};
if ~options_gsa.ppost
% y0 = squeeze(stock_smooth(:,j,:)) + ...
% kron(stock_ys(j,:),ones(size(stock_smooth,1),1));
% Out1(jsmoo).time = [1:size(y0,1)];
% Out1(jsmoo).data = y0';
Out1(jsmoo).data = j;
Out1(jsmoo).time = [pwd,'/',OutputDirectoryName];
else
......@@ -518,7 +495,7 @@ if options_gsa.glue
Out1(jsmoo).ini = 'yes';
end
end
tit(M_.exo_names_orig_ord) = M_.exo_names;
tit = M_.exo_names;
for j=1:M_.exo_nbr
Exo(j).name = tit{j};
end
......@@ -530,36 +507,24 @@ if options_gsa.glue
end
Sam.name = bayestopt_.name;
Sam.dim = [size(x) 0];
Sam.data = [x];
Sam.data = x;
Rem.id = 'Original';
Rem.ind= [1:size(x,1)];
Rem.ind= 1:size(x,1);
Info.dynare=M_.fname;
Info.order_var=dr_.order_var;
Out=Out1;
if options_gsa.ppost
% Info.dynare=M_.fname;
% Info.order_var=dr_.order_var;
% Out=Out1;
Info.TypeofSample='post';
save([OutputDirectoryName,'/',fname_,'_glue_post.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
%save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info')
else
if options_gsa.pprior
Info.TypeofSample='prior';
save([OutputDirectoryName,'/',fname_,'_glue_prior.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
% save([OutputDirectoryName,'/',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
% Out=Out1;
% save([OutputDirectoryName,'/',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
else
Info.TypeofSample='mc';
save([OutputDirectoryName,'/',fname_,'_glue_mc.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
% save([OutputDirectoryName,'/',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
% Out=Out1;
% save([OutputDirectoryName,'/',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
end
end
end
matlab/gsa/scatter_analysis.m → matlab/+gsa/scatter_analysis.m
View file @ 4a32acc2
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,
% marco.ratto@ec.europa.eu
%
% Copyright (C) 2017 European Commission
% Copyright (C) 2017 Dynare Team
% Copyright © 2017 European Commission
% Copyright © 2017-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -21,11 +23,11 @@ function indmcf = scatter_analysis(lpmat, xdata, options_scatter, DynareOptions)
% 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 <http://www.gnu.org/licenses/>.
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
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
View file @ 4a32acc2
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,
% marco.ratto@ec.europa.eu
% Copyright (C) 2014-2016 European Commission
% Copyright (C) 2014-2017 Dynare Team
% Copyright © 2014-2016 European Commission
% Copyright © 2014-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -20,7 +21,7 @@ function scatter_mcf(X,Y,vnames,plotsymbol, fnam, dirname, figtitle, xparam1, D
% 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 <http://www.gnu.org/licenses/>.
% 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
......@@ -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=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,[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
View file @ 4a32acc2
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 (C) 2017 European Commission
% Copyright (C) 2017 Dynare Team
% Copyright © 2017 European Commission
% Copyright © 2017-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -21,25 +33,9 @@ function scatter_plots(X,xp,vnames,plotsymbol, fnam, dirname, figtitle, xparam1,
% 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 <http://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
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[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,10 +69,10 @@ end
figtitle_tex=strrep(figtitle,'_','\_');
fig_nam_=[fnam];
fig_nam_=fnam;
hh=dyn_figure(DynareOptions.nodisplay,'name',figtitle);
set(hh,'userdata',{X,xp})
hh_fig=dyn_figure(options_.nodisplay,'name',figtitle);
set(hh_fig,'userdata',{X,xp})
bf = 0.1;
ffs = 0.05/(p-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,[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
View file @ 4a32acc2
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 (C) 2012-2017 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -17,9 +25,7 @@ function set_shocks_param(xparam1)
% 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 <http://www.gnu.org/licenses/>.
global estim_params_ M_
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
......
matlab/gsa/gsa_skewness.m → matlab/+gsa/skewness.m
View file @ 4a32acc2
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 (C) 2012 European Commission
% Copyright (C) 2012-2017 Dynare Team
% Copyright © 2012 European Commission
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -20,10 +26,8 @@ function s=gsa_skewness(y)
% 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 <http://www.gnu.org/licenses/>.
% 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
View file @ 4a32acc2
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
% Copyright (C) 2012 European Commission
% Copyright (C) 2012-2017 Dynare Team
% Copyright © 2012 European Commission
% Copyright © 2012-2017 Dynare Team
%
% This file is part of Dynare.
%
......@@ -22,7 +22,7 @@ function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
% 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 <http://www.gnu.org/licenses/>.
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
if nargin<3
alpha = 0.05;
......@@ -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);
......