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

Updated sensitivity package !!! 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1235 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ratto 2007-04-04 08:19:51 +00:00
parent 663fdccdec
commit 28da4a09c5
16 changed files with 1426 additions and 467 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
matlab/dat_fil_.m Normal file
View File

@ -0,0 +1,26 @@
function c = dat_fil_(data_file);
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
eval(data_file);
clear data_file;
a=who;
for j=1:length(a)
eval(['c.',a{j},'=',a{j},';']);
end

249
matlab/dynare_MC.m
View File

@ -1,179 +1,44 @@
function dynare_MC(var_list_,OutDir)
%
% Adapted by M. Ratto from dynare_estimation.m and posteriorsmoother.m
% (dynare_estimation.m and posteriorsmoother.m are part of DYNARE,
% copyright M. Juillard)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ options_ oo_ estim_params_
global bayestopt_
% temporary fix until M_.H is initialized by the parser
M_.H = [];
options_.varlist = var_list_;
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
for i = 1:size(M_.endo_names,1)
tmp = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
if ~isempty(tmp)
options_.lgyidx2varobs(i,1) = tmp;
end
end
options_ = set_default_option(options_,'first_obs',1);
options_ = set_default_option(options_,'prefilter',0);
options_ = set_default_option(options_,'presample',0);
options_ = set_default_option(options_,'lik_algo',1);
options_ = set_default_option(options_,'lik_init',1);
options_ = set_default_option(options_,'nograph',0);
options_ = set_default_option(options_,'mh_conf_sig',0.90);
options_ = set_default_option(options_,'mh_replic',20000);
options_ = set_default_option(options_,'mh_drop',0.5);
options_ = set_default_option(options_,'mh_jscale',0.2);
options_ = set_default_option(options_,'mh_init_scale',2*options_.mh_jscale);
options_ = set_default_option(options_,'mode_file','');
options_ = set_default_option(options_,'mode_compute',4);
options_ = set_default_option(options_,'mode_check',0);
options_ = set_default_option(options_,'prior_trunc',1e-10);
options_ = set_default_option(options_,'mh_mode',1);
options_ = set_default_option(options_,'mh_nblck',2);
options_ = set_default_option(options_,'load_mh_file',0);
options_ = set_default_option(options_,'nodiagnostic',0);
options_ = set_default_option(options_,'loglinear',0);
options_ = set_default_option(options_,'unit_root_vars',[]);
options_ = set_default_option(options_,'XTick',[]);
options_ = set_default_option(options_,'XTickLabel',[]);
options_ = set_default_option(options_,'bayesian_irf',0);
options_ = set_default_option(options_,'bayesian_th_moments',0);
options_ = set_default_option(options_,'TeX',0);
options_ = set_default_option(options_,'irf',40);
options_ = set_default_option(options_,'relative_irf',0);
options_ = set_default_option(options_,'order',1);
options_ = set_default_option(options_,'ar',5);
options_ = set_default_option(options_,'dr_algo',0);
options_ = set_default_option(options_,'linear',0);
options_ = set_default_option(options_,'drop',0);
options_ = set_default_option(options_,'replic',1);
options_ = set_default_option(options_,'hp_filter',0);
options_ = set_default_option(options_,'forecast',0);
options_ = set_default_option(options_,'smoother',0);
options_ = set_default_option(options_,'moments_varendo',0);
options_ = set_default_option(options_,'filtered_vars',0);
options_ = set_default_option(options_,'kalman_algo',1);
options_ = set_default_option(options_,'kalman_tol',10^(-12));
options_ = set_default_option(options_,'posterior_mode_estimation',1);
options_ = set_default_option(options_,'MaxNumberOfBytes',1e6);
options_ = set_default_option(options_,'xls_sheet','');
options_ = set_default_option(options_,'xls_range','');
options_ = set_default_option(options_,'filter_step_ahead',0);
options_ = set_default_option(options_,'diffuse_d',[]);
options_ = set_default_option(options_,'Opt6Iter',3);
options_ = set_default_option(options_,'Opt6Numb',100000);
options_ = set_default_option(options_,'steadystate_flag',0);
options_ = set_default_option(options_,'logdata',0);
options_ = set_default_option(options_,'use_mh_covariance_matrix',0);
if exist([M_.fname '_steadystate.m'])
options_.steadystate_flag = 1;
end
if options_.filtered_vars ~= 0 & options_.filter_step_ahead == 0
options_.filter_step_ahead = 1;
end
if options_.bayesian_irf & options_.irf == 0
options_.irf = 40;
end
if options_.filter_step_ahead ~= 0
options_.nk = max(options_.filter_step_ahead);
else
options_.nk = 0;
end
%% Add something to the parser ++>
M_.dname = M_.fname; % The user should be able to choose another name
% for the directory...
pnames = [' ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
n_varobs = size(options_.varobs,1);
[xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_);
bounds = prior_bounds(bayestopt_);
bounds(:,1)=max(bounds(:,1),lb);
bounds(:,2)=min(bounds(:,2),ub);
if any(xparam1 < bounds(:,1)) | any(xparam1 > bounds(:,2))
find(xparam1 < bounds(:,1))
find(xparam1 > bounds(:,2))
error('Initial parameter values are outside parameter bounds')
end
lb = bounds(:,1);
ub = bounds(:,2);
bayestopt_.lb = lb;
bayestopt_.ub = ub;
if ~isfield(options_,'trend_coeffs')
bayestopt_.with_trend = 0;
else
bayestopt_.with_trend = 1;
bayestopt_.trend_coeff = {};
trend_coeffs = options_.trend_coeffs;
nt = length(trend_coeffs);
for i=1:n_varobs
if i > length(trend_coeffs)
bayestopt_.trend_coeff{i} = '0';
else
bayestopt_.trend_coeff{i} = trend_coeffs{i};
end
end
end
bayestopt_.penalty = 1e8; % penalty
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
dr = set_state_space([]);
%% Initialization with unit-root variables
if ~isempty(options_.unit_root_vars)
n_ur = size(options_.unit_root_vars,1);
i_ur = zeros(n_ur,1);
for i=1:n_ur
i1 = strmatch(deblank(options_.unit_root_vars(i,:)),M_.endo_names(dr.order_var,:),'exact');
if isempty(i1)
error('Undeclared variable in unit_root_vars statement')
end
i_ur(i) = i1;
end
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
bayestopt_.Pinf = zeros(n1,n1);
l2 = find(l1');
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_ur,:) = l1(:,i_ur)';
l3 = l3(:);
i_ur1 = find(l3(l2));
i_stable = ones(M_.endo_nbr,1);
i_stable(i_ur) = zeros(n_ur,1);
i_stable = find(i_stable);
bayestopt_.Pinf(i_ur1,i_ur1) = diag(ones(1,length(i_ur1)));
bayestopt_.i_var_stable = i_stable;
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_stable,:) = l1(:,i_stable)';
l3 = l3(:);
bayestopt_.i_T_var_stable = find(l3(l2));
else
n1 = M_.endo_nbr;
bayestopt_.Pinf = zeros(n1,n1);
bayestopt_.Pinf(i_ur,i_ur) = diag(ones(1,length(i_ur)));
l1 = ones(M_.endo_nbr,1);
l1(i_ur,:) = zeros(length(i_ur),1);
bayestopt_.i_T_var_stable = find(l1);
end
options_.lik_init = 3;
end % if ~isempty(options_.unit_root_vars)
npar = nvx+nvn+ncx+ncn+np;
if isempty(options_.datafile)
error('ESTIMATION: datafile option is missing')
@ -183,36 +48,9 @@ if isempty(options_.varobs)
error('ESTIMATION: VAROBS is missing')
end
%% If jscale isn't specified for an estimated parameter, use
%% global option options_.jscale, set to 0.2, by default
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
%% Read and demean data
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
k = [];
k1 = [];
for i=1:n_varobs
k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')];
k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
end
% union of observed and state variables
k2 = union(k',[dr.nstatic+1:dr.nstatic+dr.npred]');
% including variables in t-2 and earlier, if any
k2 = [k2;[M_.endo_nbr+(1:dr.nspred-dr.npred)]'];
% set restrict_state to postion of observed + state variables
% in expanded state vector
bayestopt_.restrict_state = k2;
% set mf1 to positions of observed variables in restricted state vector
% for likelihood computation
[junk,bayestopt_.mf1] = ismember(k,k2);
% set mf2 to positions of observed variables in expanded state vector
% for filtering and smoothing
bayestopt_.mf2 = k;
bayestopt_.mfys = k1;
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;
@ -232,12 +70,6 @@ if ~isreal(rawdata)
' transformation'])
end
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
fname_=M_.fname;
@ -257,39 +89,40 @@ load(options_.mode_file)
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat lpmat0 istable %iunstable egg yys T
B = size(x,1);
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
else
n1 = M_.endo_nbr;
end
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(x(1,:)',gend,data);
n1=size(atT,1);
nfil=B/40;
stock_smooth = zeros(gend,n1,40);
stock_filter = zeros(gend+1,n1,40);
stock_ys = zeros(M_.endo_nbr,40);
stock_smooth = zeros(M_.endo_nbr,gend,40);
stock_filter = zeros(M_.endo_nbr,gend+1,40);
stock_ys = zeros(40, M_.endo_nbr);
logpo2=zeros(B,1);
%%
h = waitbar(0,'MC smoother ...');
delete([OutDir,'\',namfile,'_*.mat'])
ib=0;
ifil=0;
opt_gsa=options_.opt_gsa;
for b=1:B
ib=ib+1;
deep = x(b,:);
deep = x(b,:)';
set_all_parameters(deep);
dr = resol(oo_.steady_state,0);
%deep(1:offset) = xparam1(1:offset);
logpo2(b,1) = DsgeLikelihood(deep',gend,data);
logpo2(b,1) = DsgeLikelihood(deep,gend,data);
if opt_gsa.lik_only==0,
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(deep,gend,data);
stock_smooth(:,:,ib)=atT';
stock_filter(:,:,ib)=filtered_state_vector';
stock_ys(:,ib)=ys;
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(deep,gend,data);
stock_smooth(:,:,ib)=atT(1:M_.endo_nbr,:);
stock_filter(:,:,ib)=filtered_state_vector(1:M_.endo_nbr,:);
%stock_filter(:,:,ib)=aK(options_.filter_step_ahead,1:M_.endo_nbr,:);
stock_ys(ib,:)=ys';
if ib==40,
ib=0;
ifil=ifil+1;
save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
stock_smooth = zeros(gend,n1,40);
stock_filter = zeros(gend+1,n1,40);
stock_ys = zeros(M_.endo_nbr,40);
stock_smooth = zeros(M_.endo_nbr,gend,40);
stock_filter = zeros(M_.endo_nbr,gend+1,40);
stock_ys = zeros(40, M_.endo_nbr);
end
end
waitbar(b/B,h,['MC smoother ...',num2str(b),'/',num2str(B)]);
@ -300,7 +133,7 @@ if ib>0,
ifil=ifil+1;
stock_smooth = stock_smooth(:,:,1:ib);
stock_filter = stock_filter(:,:,1:ib);
stock_ys = stock_ys(:,1:ib);
stock_ys = stock_ys(1:ib,:);
save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
end
end

200
matlab/dynare_sensitivity.m
View File

@ -1,5 +1,19 @@
function x0=dynare_sensitivity()
% copyright Marco Ratto 2006
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ options_ oo_ bayestopt_
@ -10,29 +24,100 @@ x0=[];
options_ = set_default_option(options_,'opt_gsa',1);
options_gsa_ = options_.opt_gsa;
options_gsa_ = set_default_option(options_gsa_,'identification',0);
if options_gsa_.identification,
options_gsa_.redform=0;
options_gsa_ = set_default_option(options_gsa_,'morris',1);
end
% map stability
options_gsa_ = set_default_option(options_gsa_,'stab',1);
options_gsa_ = set_default_option(options_gsa_,'redform',0);
options_gsa_ = set_default_option(options_gsa_,'pprior',1);
options_gsa_ = set_default_option(options_gsa_,'prior_range',1);
options_gsa_ = set_default_option(options_gsa_,'ppost',0);
options_gsa_ = set_default_option(options_gsa_,'ilptau',1);
options_gsa_ = set_default_option(options_gsa_,'morris',0);
options_gsa_ = set_default_option(options_gsa_,'glue',0);
options_gsa_ = set_default_option(options_gsa_,'morris_nliv',6);
options_gsa_ = set_default_option(options_gsa_,'morris_ntra',20);
options_gsa_ = set_default_option(options_gsa_,'Nsam',2048);
options_gsa_ = set_default_option(options_gsa_,'load_redform',0);
options_gsa_ = set_default_option(options_gsa_,'load_rmse',0);
options_gsa_ = set_default_option(options_gsa_,'load_stab',0);
options_gsa_ = set_default_option(options_gsa_,'alpha2_stab',0.3);
options_gsa_ = set_default_option(options_gsa_,'ksstat',0.1);
%options_gsa_ = set_default_option(options_gsa_,'load_mh',0);
OutputDirectoryName = CheckPath('GSA');
if options_gsa_.redform,
options_gsa_.pprior=1;
options_gsa_.ppost=0;
end
if options_gsa_.morris,
if ~(exist('Sampling_Function_2','file')==6 | exist('Sampling_Function_2','file')==2),
disp('Download pre-parsed mapping routines at:')
disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
disp(' ' )
error('Mapping routines missing!')
end
options_gsa_.pprior=1;
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_.ksstat=1;
if options_gsa_.morris==2,
options_gsa_ = set_default_option(options_gsa_,'Nsam',256);
OutputDirectoryName = CheckPath('GSA/IDENTIF');
else
OutputDirectoryName = CheckPath('GSA/SCREEN');
end
else
OutputDirectoryName = CheckPath('GSA');
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 isempty(ls(filetoload)),
disp([filetoload,' not found!'])
disp(['You asked to load a non existent analysis'])
%options_gsa_.load_stab=0;
return,
else
if isempty(strmatch('bkpprior',who('-file', filetoload),'exact')),
disp('Warning! Missing prior info for saved sample') % trap for files previous
disp('The saved files are generated with previous version of GSA package') % trap for files previous
else
load(filetoload,'bkpprior'),
if any(bayestopt_.pshape~=bkpprior.pshape) | ...
any(bayestopt_.p1~=bkpprior.p1) | ...
any(bayestopt_.p2~=bkpprior.p2) | ...
any(bayestopt_.p3(~isnan(bayestopt_.p3))~=bkpprior.p3(~isnan(bkpprior.p3))) | ...
any(bayestopt_.p4(~isnan(bayestopt_.p4))~=bkpprior.p4(~isnan(bkpprior.p4))),
disp('WARNING!')
disp('The saved sample has different priors w.r.t. to current ones!!')
disp('')
disp('Press ENTER to continue')
pause;
end
end
end
end
if options_gsa_.stab & ~options_gsa_.ppost,
x0 = stab_map_(options_gsa_.Nsam, options_gsa_.load_stab, options_gsa_.ksstat, options_gsa_.alpha2_stab, ...
options_gsa_.redform, options_gsa_.pprior, options_gsa_.ilptau, OutputDirectoryName);
x0 = stab_map_(OutputDirectoryName);
end
% reduced form
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
options_gsa_ = set_default_option(options_gsa_,'load_redform',0);
options_gsa_ = set_default_option(options_gsa_,'logtrans_redform',0);
options_gsa_ = set_default_option(options_gsa_,'threshold_redform',[]);
options_gsa_ = set_default_option(options_gsa_,'ksstat_redform',0.1);
@ -41,20 +126,22 @@ options_gsa_ = set_default_option(options_gsa_,'namendo',[]);
options_gsa_ = set_default_option(options_gsa_,'namlagendo',[]);
options_gsa_ = set_default_option(options_gsa_,'namexo',[]);
if options_gsa_.redform & ~isempty(options_gsa_.namendo) & ~options_gsa_.ppost,
redform_map(options_gsa_.namendo, options_gsa_.namlagendo, options_gsa_.namexo, ...
options_gsa_.load_redform, options_gsa_.pprior, options_gsa_.logtrans_redform, ...
options_gsa_.threshold_redform, options_gsa_.ksstat_redform, ...
options_gsa_.alpha2_redform, OutputDirectoryName);
options_.opt_gsa = options_gsa_;
if options_gsa_.redform & ~isempty(options_gsa_.namendo) ...
& ~options_gsa_.ppost,
if options_gsa_.morris,
redform_screen(OutputDirectoryName);
else
redform_map(OutputDirectoryName);
end
end
% RMSE mapping
% function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
options_gsa_ = set_default_option(options_gsa_,'rmse',0);
options_gsa_ = set_default_option(options_gsa_,'lik_only',0);
options_gsa_ = set_default_option(options_gsa_,'var_rmse',options_.varobs);
options_gsa_ = set_default_option(options_gsa_,'load_rmse',0);
options_gsa_ = set_default_option(options_gsa_,'pfilt_rmse',0.1);
options_gsa_ = set_default_option(options_gsa_,'istart_rmse',1);
options_gsa_ = set_default_option(options_gsa_,'istart_rmse',options_.presample+1);
options_gsa_ = set_default_option(options_gsa_,'alpha_rmse',0.002);
options_gsa_ = set_default_option(options_gsa_,'alpha2_rmse',1);
options_.opt_gsa = options_gsa_;
@ -68,16 +155,27 @@ if options_gsa_.rmse,
if isempty(a),
dynare_MC([],OutputDirectoryName);
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_(options_gsa_.var_rmse, options_gsa_.load_rmse, options_gsa_.pfilt_rmse, ...
options_gsa_.alpha_rmse, options_gsa_.alpha2_rmse, OutputDirectoryName, ...
options_gsa_.istart_rmse);
filt_mc_(OutputDirectoryName);
end
options_gsa_ = set_default_option(options_gsa_,'glue',0);
if options_gsa_.glue,
dr_ = oo_.dr;
if options_gsa_.ppost
@ -90,6 +188,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'])
return,
end
nruns=size(x,1);
gend = options_.nobs;
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
@ -114,12 +217,14 @@ if options_gsa_.glue,
jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
js = strmatch(vj,lgy_,'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)];
% 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
Out(j).data = jxj;
Out(j).time = [pwd,'\',DirectoryName];
@ -132,12 +237,13 @@ if options_gsa_.glue,
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';
% 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
@ -150,10 +256,12 @@ if options_gsa_.glue,
jsmoo=jsmoo+1;
vj=deblank(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';
% 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
Out1(jsmoo).data = j;
Out1(jsmoo).time = [pwd,'\',DirectoryName];
@ -179,22 +287,30 @@ if options_gsa_.glue,
Rem.id = 'Original';
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;
save([OutputDirectoryName,'\',fname_,'_post_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
% Info.dynare=M_.fname;
% Info.order_var=dr_.order_var;
% Out=Out1;
Info.TypeofSample='post';
save([OutputDirectoryName,'\',fname_,'_glue_post'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
%save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info')
else
if options_gsa_.pprior
save([OutputDirectoryName,'\',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
Out=Out1;
save([OutputDirectoryName,'\',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
Info.TypeofSample='prior';
save([OutputDirectoryName,'\',fname_,'_glue_prior'], '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
save([OutputDirectoryName,'\',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
Out=Out1;
save([OutputDirectoryName,'\',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
Info.TypeofSample='mc';
save([OutputDirectoryName,'\',fname_,'_glue_mc'], '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

269
matlab/filt_mc_.m
View File

@ -1,29 +1,42 @@
function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2, OutDir, istart, alphaPC)
function [rmse_MC, ixx] = filt_mc_(OutDir)
% function [rmse_MC, ixx] = filt_mc_(OutDir)
% copyright Marco Ratto 2006
% inputs (from opt_gsa structure)
% vvarvecm = options_gsa_.var_rmse;
% loadSA = options_gsa_.load_rmse;
% pfilt = options_gsa_.pfilt_rmse;
% alpha = options_gsa_.alpha_rmse;
% alpha2 = options_gsa_.alpha2_rmse;
% istart = options_gsa_.istart_rmse;
% alphaPC = 0.5;
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global bayestopt_ estim_params_ M_ options_ oo_
if nargin<1 | isempty(vvarvecm),
vvarvecm = options_.varobs;
end
if nargin<2,
loadSA=0;
end
if nargin<3 | isempty(pfilt),
pfilt=0.1; % cut the best 10% of runs
end
if nargin<4 | isempty(alpha),
alpha=0.002;
end
if nargin<5 | isempty(alpha2),
alpha2=0.5;
end
if nargin<7 | isempty(istart),
istart=1;
end
if nargin<8,
alphaPC=0.5;
end
options_gsa_=options_.opt_gsa;
vvarvecm = options_gsa_.var_rmse;
loadSA = options_gsa_.load_rmse;
pfilt = options_gsa_.pfilt_rmse;
alpha = options_gsa_.alpha_rmse;
alpha2 = options_gsa_.alpha2_rmse;
istart = options_gsa_.istart_rmse;
alphaPC = 0.5;
fname_ = M_.fname;
lgy_ = M_.endo_names;
dr_ = oo_.dr;
@ -35,13 +48,18 @@ disp('for the fit of EACH observed series ...')
disp(' ')
disp('Deleting old SA figures...')
a=dir([OutDir,'\*.*']);
tmp1='0';
if options_.opt_gsa.ppost,
tmp=['_SA_fit_post'];
tmp=['_rmse_post'];
else
if options_.opt_gsa.pprior
tmp=['_SA_fit_prior'];
tmp=['_rmse_prior'];
else
tmp=['_SA_fit_mc'];
tmp=['_rmse_mc'];
end
if options_gsa_.lik_only,
tmp1 = [tmp,'_post_SA'];
tmp = [tmp,'_lik_SA'];
end
end
for j=1:length(a),
@ -49,28 +67,25 @@ for j=1:length(a),
disp(a(j).name)
delete([OutDir,'\',a(j).name])
end,
if strmatch([fname_,tmp1],a(j).name),
disp(a(j).name)
delete([OutDir,'\',a(j).name])
end,
end
disp('done !')
nshock=estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
npar=estim_params_.np;
for j=1:npar+nshock,
if j>nshock
if isfield(oo_,'posterior_mode'),
xparam1(j)=oo_.posterior_mode.parameters.(bayestopt_.name{j});
end
if isfield(oo_,'posterior_mean'),
xparam1_mean(j)=oo_.posterior_mean.parameters.(bayestopt_.name{j});
end
else
if isfield(oo_,'posterior_mode'),
xparam1(j)=oo_.posterior_mode.shocks_std.(bayestopt_.name{j});
end
if isfield(oo_,'posterior_mean'),
xparam1_mean(j)=oo_.posterior_mean.shocks_std.(bayestopt_.name{j});
end
end
load(options_.mode_file,'xparam1')
if options_.opt_gsa.ppost,
c=load([fname_,'_mean'],'xparam1');
xparam1_mean=c.xparam1;
clear c
elseif ~isempty(ls([fname_,'_mean.mat']))
c=load([fname_,'_mean'],'xparam1');
xparam1_mean=c.xparam1;
clear c
end
if options_.opt_gsa.ppost,
@ -94,11 +109,14 @@ if ~loadSA,
steady_;
ys_mean=oo_.steady_state;
end
eval(options_.datafile)
% eval(options_.datafile)
obs = dat_fil_(options_.datafile);
if ~options_.opt_gsa.ppost
load([OutDir,'\',fnamtmp]);
load([OutDir,'\',fnamtmp],'x','logpo2','stock_gend','stock_data');
logpo2=-logpo2;
else
load([DirectoryName '/' M_.fname '_data.mat']);
%load([DirectoryName '/' M_.fname '_data.mat']);
[stock_gend, stock_data] = read_data;
filfilt = dir([DirectoryName '/' M_.fname '_filter*.mat']);
filparam = dir([DirectoryName '/' M_.fname '_param*.mat']);
x=[];
@ -114,7 +132,6 @@ if ~loadSA,
clear stock stock_logpo stock_ys;
end
end
logpo2=-logpo2;
end
nruns=size(x,1);
nfilt=floor(pfilt*nruns);
@ -125,15 +142,17 @@ if ~loadSA,
nobs=options_.nobs;
for i=1:size(vvarvecm,1),
vj=deblank(vvarvecm(i,:));
eval(['vobs =obs.',vj,'(fobs:fobs-1+nobs);'])
if options_.prefilter == 1
%eval([vj,'=',vj,'-bayestopt_.mean_varobs(i);'])
eval([vj,'=',vj,'-mean(',vj,',1);'])
%eval([vj,'=',vj,'-mean(',vj,',1);'])
vobs = vobs-mean(vobs,1);
end
jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
js = strmatch(vj,lgy_,'exact');
if exist('xparam1','var')
eval(['rmse_mode(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));'])
eval(['rmse_mode(i) = sqrt(mean((vobs(istart:end)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));'])
end
y0=zeros(nobs+1,nruns);
if options_.opt_gsa.ppost
@ -151,19 +170,31 @@ if ~loadSA,
clear stock;
end
else
y0 = squeeze(stock_filter(:,jxj,:)) + ...
kron(stock_ys(js,:),ones(size(stock_filter,1),1));
filfilt=ls([OutDir,'\',fnamtmp,'_*.mat']);
nbb=0;
for j=1:size(filfilt,1),
load([OutDir,'\',fnamtmp,'_',num2str(j),'.mat'],'stock_filter','stock_ys');
nb = size(stock_filter,3);
y0(:,nbb+1:nbb+nb) = squeeze(stock_filter(jxj,:,:)) + ...
kron(stock_ys(:,js)',ones(nobs+1,1));
%y0(:,:,size(y0,3):size(y0,3)+size(stock,3))=stock;
nbb=nbb+nb;
clear stock_filter;
end
% y0 = squeeze(stock_filter(:,jxj,:)) + ...
% kron(stock_ys(js,:),ones(size(stock_filter,1),1));
end
y0M=mean(y0,2);
for j=1:nruns,
eval(['rmse_MC(j,i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0(istart:end-1,j)).^2));'])
rmse_MC(j,i) = sqrt(mean((vobs(istart:end)-y0(istart:end-1,j)).^2));
end
if exist('xparam1_mean','var')
%eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0M(istart:end-1)).^2));'])
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,stock_gend,stock_data);
y0 = ahat(jxj,:)' + ...
kron(ys_mean(js,:),ones(size(ahat,2),1));
eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0(istart:end-1)).^2));'])
rmse_pmean(i) = sqrt(mean((vobs(istart:end)-y0(istart:end-1)).^2));
end
end
clear stock_filter;
@ -171,7 +202,7 @@ if ~loadSA,
for j=1:nruns,
lnprior(j,1) = priordens(x(j,:),bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
end
likelihood=logpo2(:)+lnprior(:);
likelihood=logpo2(:)-lnprior(:);
disp('... done!')
if options_.opt_gsa.ppost
@ -180,7 +211,11 @@ if ~loadSA,
if options_.opt_gsa.lik_only
save([OutDir,'\',fnamtmp], 'likelihood', '-append')
else
save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean','-append')
if exist('xparam1_mean','var')
save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean','-append')
else
save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','-append')
end
end
end
else
@ -189,7 +224,7 @@ else
else
load([OutDir,'\',fnamtmp],'x','logpo2','likelihood','rmse_MC','rmse_mode','rmse_pmean');
end
lnprior=likelihood(:)-logpo2(:);
lnprior=logpo2(:)-likelihood(:);
nruns=size(x,1);
nfilt=floor(pfilt*nruns);
end
@ -197,21 +232,21 @@ end
nfilt0=nfilt*ones(size(vvarvecm,1),1);
logpo2=logpo2(:);
if ~options_.opt_gsa.ppost
[dum, ipost]=sort(logpo2);
[dum, ilik]=sort(likelihood);
[dum, ipost]=sort(-logpo2);
[dum, ilik]=sort(-likelihood);
end
if ~options_.opt_gsa.ppost & options_.opt_gsa.lik_only
if options_.opt_gsa.pprior
anam='SA_fit_prior_post';
anam='rmse_prior_post';
else
anam='SA_fit_mc_post';
anam='rmse_mc_post';
end
stab_map_1(x, ipost(1:nfilt), ipost(nfilt+1:end), anam, 1,[],OutDir);
stab_map_2(x(ipost(1:nfilt),:),alpha2,anam, OutDir);
if options_.opt_gsa.pprior
anam='SA_fit_prior_lik';
anam='rmse_prior_lik';
else
anam='SA_fit_mc_lik';
anam='rmse_mc_lik';
end
stab_map_1(x, ilik(1:nfilt), ilik(nfilt+1:end), anam, 1,[],OutDir);
stab_map_2(x(ilik(1:nfilt),:),alpha2,anam, OutDir);
@ -258,18 +293,18 @@ for i=1:size(vvarvecm,1),
title(vvarvecm(i,:))
if mod(i,9)==0 | i==size(vvarvecm,1)
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnprior',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_post_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_lnprior',int2str(ifig)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnprior',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_lnprior',int2str(ifig)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnprior',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_lnprior',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_lnprior',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_lnprior',int2str(ifig)]);
end
end
close(gcf)
@ -288,20 +323,23 @@ for i=1:size(vvarvecm,1),
h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
title(vvarvecm(i,:))
if options_.opt_gsa.ppost==0,
set(gca,'xlim',[min( likelihood(ixx(1:nfilt0(i),i)) ) max( likelihood(ixx(1:nfilt0(i),i)) )])
end
if mod(i,9)==0 | i==size(vvarvecm,1)
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnlik',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_post_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_lnlik',int2str(ifig)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnlik',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_lnlik',int2str(ifig)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnlik',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_lnlik',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_lnlik',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_lnlik',int2str(ifig)]);
end
end
close(gcf)
@ -320,20 +358,23 @@ for i=1:size(vvarvecm,1),
h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
title(vvarvecm(i,:))
if options_.opt_gsa.ppost==0,
set(gca,'xlim',[min( logpo2(ixx(1:nfilt0(i),i)) ) max( logpo2(ixx(1:nfilt0(i),i)) )])
end
if mod(i,9)==0 | i==size(vvarvecm,1)
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnpost',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_post_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_lnpost',int2str(ifig)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnpost',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_lnpost',int2str(ifig)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnpost',int2str(ifig)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_lnpost',int2str(ifig)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_lnpost',int2str(ifig)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_lnpost',int2str(ifig)]);
end
end
close(gcf)
@ -366,11 +407,11 @@ vvarvecm=vvarvecm(ivar,:);
rmse_MC=rmse_MC(:,ivar);
disp(' ')
if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior,
% if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior,
disp(['Sample filtered the ',num2str(pfilt*100),'% best RMSE''s for each observed series ...' ])
else
disp(['Sample filtered the best RMSE''s smaller than RMSE at the posterior mean ...' ])
end
% else
% disp(['Sample filtered the best RMSE''s smaller than RMSE at the posterior mean ...' ])
% end
% figure, boxplot(rmse_MC)
% set(gca,'xticklabel',vvarvecm)
% saveas(gcf,[fname_,'_SA_RMSE'])
@ -462,18 +503,18 @@ for ix=1:ceil(length(nsnam)/6),
set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
title([pnam{nsnam(j)}],'interpreter','none')
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_' int2str(ix)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_post_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_' int2str(ix)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_' int2str(ix)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_' int2str(ix)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_' int2str(ix)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_' int2str(ix)]);
end
end
end
@ -521,18 +562,18 @@ for ix=1:ceil(length(nsnam)/6),
set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
title([pnam{nsnam(j)}],'interpreter','none')
if options_.opt_gsa.ppost
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_dens_' int2str(ix)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_post_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_dens_' int2str(ix)]);
else
if options_.opt_gsa.pprior
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_dens_' int2str(ix)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_dens_' int2str(ix)]);
else
saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_dens_' int2str(ix)]);
saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_dens_',num2str(ix)])
eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_dens_' int2str(ix)]);
eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_dens_' int2str(ix)]);
end
end
end
@ -576,7 +617,7 @@ close all
% title([pnam{np(i)},'. K-S prob ', num2str(PP(np(i),j))],'interpreter','none')
% xlabel('')
% if mod(i,12)==0 | i==nsx(j),
% saveas(gcf,[fname_,'_SA_fit_',deblank(vvarvecm(j,:)),'_',int2str(nfig)])
% saveas(gcf,[fname_,'_rmse_',deblank(vvarvecm(j,:)),'_',int2str(nfig)])
% close(gcf)
% end
% end
@ -605,12 +646,12 @@ disp('Starting bivariate analysis:')
for i=1:size(vvarvecm,1)
if options_.opt_gsa.ppost
fnam = ['SA_fit_post_',deblank(vvarvecm(i,:))];
fnam = ['rmse_post_',deblank(vvarvecm(i,:))];
else
if options_.opt_gsa.pprior
fnam = ['SA_fit_prior_',deblank(vvarvecm(i,:))];
fnam = ['rmse_prior_',deblank(vvarvecm(i,:))];
else
fnam = ['SA_fit_mc_',deblank(vvarvecm(i,:))];
fnam = ['rmse_mc_',deblank(vvarvecm(i,:))];
end
end
stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,fnam, OutDir);

41
matlab/ghx2transition.m Normal file
View File

@ -0,0 +1,41 @@
function [A,B] = ghx2transition(mm,iv,ic,aux)
% [A,B] = ghx2transition(mm,iv,ic,aux)
%
% Adapted by M. Ratto from kalman_transition_matrix.m
% (kalman_transition_matrix.m is part of DYNARE, copyright M. Juillard)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_
[nr1, nc1] = size(mm);
ghx = mm(:, [1:(nc1-M_.exo_nbr)]);
ghu = mm(:, [(nc1-M_.exo_nbr+1):end] );
n_iv = length(iv);
n_ir1 = size(aux,1);
nr = n_iv + n_ir1;
A = zeros(nr,nr);
B = zeros(nr,M_.exo_nbr);
i_n_iv = 1:n_iv;
A(i_n_iv,ic) = ghx(iv,:);
if n_ir1 > 0
A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr);
end
B(i_n_iv,:) = ghu(iv,:);

74
matlab/gsa_sdp_dyn.m Normal file
View File

@ -0,0 +1,74 @@
function gsa_ = gsa_sdp_dyn(varargin)
% function gsa_ =
% gsa_sdp_dyn(y, x, tvp, nvr, ialias, nvrT, ifig, fname, pnames, yi_anal)
% 1 2 3 4 5 6 7 8 9 10
% INPUTS
% y, the output
% x, the inputs
% tvp: tipe of estimation
% 0 = load previous estimation
% 1 = RW for all parameters
% 2 = IRW for all parameters
% -1 = RW for all parameters, fixed nvr's (no estimation)
% -2 = IRW for all parameters, fixed nvr's (no estimation)
% an array of length(x) specifies the process for each parameter
% nvr initial nvr value for estimation or fixed nvr is no estimation
% ialias to be used with LPTAU samples < 1024
% nvrT maximum nvr allowable
% ifig: 1 plot figures, 0 no plots
% fname, file name to save the analysis
% pnames: names of the params
% yi_anal, analytical value of the first order HDMR terms (optional)
%
% OUTPUT
% gsa_.univariate.f
% gsa_.univariate.fs
% gsa_.univariate.fses
% gsa_.univariate.si
% gsa_.univariate.si_std
% gsa_.multivariate.f
% gsa_.multivariate.fs
% gsa_.multivariate.fses
% gsa_.multivariate.si
% gsa_.multivariate.si_std
% gsa_.multivariate.stat
% gsa_.x0
% gsa_.xx
% gsa_.y
%
% f, function estimates
% fs, sorted function estimates
% fses, sorted standard error of function estimates
% si, sensitivity indices
% si_std, st. error of sensitivity indices
% stat, euristic tstat for significance of f
% xx, transformed inputs used (rank-transformed)
% x0, original inputs
% y, output
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
if exist('sdr','file')==6 | exist('sdr','file')==2,
gsa_ = gsa_sdp_fn(varargin{:});
else
disp('Download the SDP mapping routines at:')
disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
disp(' ' )
error('SDP mapping routines missing!')
end

125
matlab/prior_draw_gsa.m Normal file
View File

@ -0,0 +1,125 @@
function pdraw = prior_draw_gsa(init,rdraw,cc)
% Draws from the prior distributions
% Adapted by M. Ratto from prior_draw (of DYNARE, copyright M. Juillard),
% for use with Sensitivity Toolbox for DYNARE
%
%
% INPUTS
% o init [integer] scalar equal to 1 (first call) or 0.
% o rdraw
% o cc [double] two columns matrix (same as in
% metropolis.m), constraints over the
% parameter space (upper and lower bounds).
%
% OUTPUTS
% o pdraw [double] draw from the joint prior density.
%
% ALGORITHM
% ...
%
% SPECIAL REQUIREMENTS
% MATLAB Statistics Toolbox
%
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ options_ estim_params_ bayestopt_
persistent fname npar bounds pshape pmean pstd a b p1 p2 p3 p4 condition
if init
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
MhDirectoryName = CheckPath('metropolis');
fname = [ MhDirectoryName '/' M_.fname];
pshape = bayestopt_.pshape;
pmean = bayestopt_.pmean;
pstd = bayestopt_.pstdev;
p1 = bayestopt_.p1;
p2 = bayestopt_.p2;
p3 = bayestopt_.p3;
p4 = bayestopt_.p4;
a = zeros(npar,1);
b = zeros(npar,1);
if nargin == 2
bounds = cc;
else
bounds = kron(ones(npar,1),[-Inf Inf]);
end
for i = 1:npar
switch pshape(i)
case 3% Gaussian prior
b(i) = pstd(i)^2/(pmean(i)-p3(i));
a(i) = (pmean(i)-p3(i))/b(i);
case 1% Beta prior
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
a(i) = (1-mu)*mu^2/stdd^2 - mu;
b(i) = a(i)*(1/mu - 1);
case 2;%Gamma prior
mu = p1(i)-p3(i);
b(i) = p2(i)^2/mu;
a(i) = mu/b(i);
case {5,4,6}
% Nothing to do here
%
% 4: Inverse gamma, type 1, prior
% p2(i) = nu
% p1(i) = s
% 6: Inverse gamma, type 2, prior
% p2(i) = nu
% p1(i) = s
% 5: Uniform prior
% p3(i) and p4(i) are used.
otherwise
disp('prior_draw :: Error!')
disp('Unknown prior shape.')
return
end
pdraw = zeros(npar,1);
end
condition = 1;
pdraw = zeros(npar,1);
return
end
for i = 1:npar
switch pshape(i)
case 5% Uniform prior.
pdraw(:,i) = rdraw(:,i)*(p4(i)-p3(i)) + p3(i);
case 3% Gaussian prior.
pdraw(:,i) = norminv(rdraw(:,i),pmean(i),pstd(i));
case 2% Gamma prior.
pdraw(:,i) = gaminv(rdraw(:,i),a(i),b(i))+p3(i);
case 1% Beta distribution (TODO: generalized beta distribution)
pdraw(:,i) = betainv(rdraw(:,i),a(i),b(i))*(p4(i)-p3(i))+p3(i);
case 4% INV-GAMMA1 distribution
% TO BE CHECKED
pdraw(:,i) = sqrt(1./gaminv(rdraw(:,i),p2(i)/2,1/p1(i)));
case 6% INV-GAMMA2 distribution
% TO BE CHECKED
pdraw(:,i) = 1./gaminv(rdraw(:,i),p2(i)/2,1/p1(i));
otherwise
% Nothing to do here.
end
end

39
matlab/read_data.m Normal file
View File

@ -0,0 +1,39 @@
function [gend, data] = read_data
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global options_ bayestopt_
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
if options_.loglinear == 1 & ~options_.logdata
rawdata = log(rawdata);
end
if options_.prefilter == 1
bayestopt_.mean_varobs = mean(rawdata,1);
data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
else
data = transpose(rawdata);
end
if ~isreal(rawdata)
error(['There are complex values in the data. Probably a wrong' ...
' transformation'])
end

314
matlab/redform_map.m Normal file
View File

@ -0,0 +1,314 @@
function redform_map(dirname)
%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;
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ oo_ estim_params_ options_ bayestopt_
options_gsa_ = options_.opt_gsa;
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;
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
if nargin==0,
dirname='';
end
if pprior
load([dirname,'\',M_.fname,'_prior']);
adir=[dirname '\redform_stab'];
else
load([dirname,'\',M_.fname,'_mc']);
adir=[dirname '\redform_mc'];
end
if isempty(dir(adir))
mkdir(adir)
end
adir0=pwd;
%cd(adir)
nspred=size(T,2)-M_.exo_nbr;
x0=lpmat(istable,:);
[kn, np]=size(x0);
offset = length(bayestopt_.pshape)-np;
pshape = bayestopt_.pshape(offset+1:end);
pd = [bayestopt_.p1(offset+1:end) bayestopt_.p2(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)];
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
clear lpmat lpmat0 egg iunstable yys
js=0;
for j=1:size(anamendo,1)
namendo=deblank(anamendo(j,:));
iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact');
ifig=0;
iplo=0;
for jx=1:size(anamexo,1)
namexo=deblank(anamexo(jx,:));
iexo=strmatch(namexo,M_.exo_names,'exact');
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,
ifig=ifig+1;
figure('name',[namendo,' vs. shocks ',int2str(ifig)]),
iplo=0;
end
iplo=iplo+1;
js=js+1;
xdir0 = [adir,'\',namendo,'_vs_', namexo];
if ilog==0,
if isempty(threshold)
si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0);
else
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
xdir = [xdir0,'_cut'];
si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir);
delete([xdir, '\*cut*.*'])
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
indsmirnov = find(dproba>ksstat);
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
end
else
[yy, xdir] = log_trans_(y0,xdir0);
silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir);
end
subplot(3,3,iplo),
if ilog,
[saso, iso] = sort(-silog(:,js));
bar([silog(iso(1:10),js)])
logflag='log';
else
[saso, iso] = sort(-si(:,js));
bar(si(iso(1:10),js))
logflag='';
end
%set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:10,
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none')
if iplo==9,
saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
close(gcf)
end
end
end
end
if iplo<9 & iplo>0 & ifig,
saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
close(gcf)
end
ifig=0;
iplo=0;
for je=1:size(anamlagendo,1)
namlagendo=deblank(anamlagendo(je,:));
ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
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,
ifig=ifig+1;
figure('name',[namendo,' vs. lags ',int2str(ifig)]),
iplo=0;
end
iplo=iplo+1;
js=js+1;
xdir0 = [adir,'\',namendo,'_vs_', namlagendo];
if ilog==0,
if isempty(threshold)
si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0);
else
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
xdir = [xdir0,'_cut'];
si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir);
delete([xdir, '\*cut*.*'])
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
indsmirnov = find(dproba>ksstat);
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
end
else
[yy, xdir] = log_trans_(y0,xdir0);
silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir);
end
subplot(3,3,iplo),
if ilog,
[saso, iso] = sort(-silog(:,js));
bar([silog(iso(1:10),js)])
logflag='log';
else
[saso, iso] = sort(-si(:,js));
bar(si(iso(1:10),js))
logflag='';
end
%set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:10,
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
if iplo==9,
saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
close(gcf)
end
end
end
end
if iplo<9 & iplo>0 & ifig,
saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
close(gcf)
end
end
if ilog==0,
figure, bar(si)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
set(gca,'xlim',[0.5 np+0.5])
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')
saveas(gcf,[dirname,'\',M_.fname,'_redform_gsa'])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_gsa']);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_gsa']);
else
figure, bar(silog)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
set(gca,'xlim',[0.5 np+0.5])
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')
saveas(gcf,[dirname,'\',M_.fname,'_redform_gsa_log'])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_gsa_log']);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_gsa_log']);
end
function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir)
global bayestopt_
fname=[xdir,'\map'];
if iload==0,
figure, hist(y0,30), title([namy,' vs. ', namx])
if isempty(dir(xdir))
mkdir(xdir)
end
saveas(gcf,[xdir,'\', namy,'_vs_', namx])
eval(['print -depsc2 ' xdir,'\', namy,'_vs_', namx]);
eval(['print -dpdf ' xdir,'\', namy,'_vs_', namx]);
close(gcf)
gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames);
else
gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames);
end
si = gsa_.multivariate.si;
function [yy, xdir]=log_trans_(y0,xdir0)
f=inline('skewness(log(y+lam))','lam','y');
if ~(max(y0)<0 | min(y0)>0)
isig=1;
if skewness(y0)<0,
isig=-1;
y0=-y0;
end
n=hist(y0,10);
if n(1)>20*n(end),
try lam=fzero(f,[-min(y0)+10*eps -min(y0)+abs(median(y0))],[],y0);
catch
yl(1)=f(-min(y0)+10*eps,y0);
yl(2)=f(-min(y0)+abs(median(y0)),y0);
if abs(yl(1))<abs(yl(2))
lam=-min(y0)+eps;
else
lam = -min(y0)+abs(median(y0)); %abs(100*(1+min(y0)));
end
end
yy = log(y0+lam);
xdir=[xdir0,'_logskew'];
else
yy = log(y0.^2);
xdir=[xdir0,'_logsquared'];
end
else
if max(y0)<0
y0=-y0;
%yy=log(-y0);
xdir=[xdir0,'_minuslog'];
elseif min(y0)>0
%yy=log(y0);
xdir=[xdir0,'_log'];
end
try lam=fzero(f,[-min(y0)+10*eps -min(y0)+median(y0)],[],y0);
catch
yl(1)=f(-min(y0)+10*eps,y0);
yl(2)=f(-min(y0)+abs(median(y0)),y0);
if abs(yl(1))<abs(yl(2))
lam=-min(y0)+eps;
else
lam = -min(y0)+abs(median(y0)); %abs(100*(1+min(y0)));
end
end
yy = log(y0+lam);
end

166
matlab/redform_screen.m Normal file
View File

@ -0,0 +1,166 @@
function redform_screen(dirname)
%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;
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ oo_ estim_params_ options_ bayestopt_
options_gsa_ = options_.opt_gsa;
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
anamexo = options_gsa_.namexo;
iload = options_gsa_.load_redform;
nliv = options_gsa_.morris_nliv;
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
if nargin==0,
dirname='';
end
load([dirname,'\',M_.fname,'_prior'],'lpmat','istable','T');
nspred=oo_.dr.nspred;
[kn, np]=size(lpmat);
offset = length(bayestopt_.pshape)-np;
pshape = bayestopt_.pshape(offset+1:end);
pd = [bayestopt_.p1(offset+1:end) bayestopt_.p2(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)];
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');
iplo=0;
ifig=0;
for jx=1:size(anamexo,1)
namexo=deblank(anamexo(jx,:));
iexo=strmatch(namexo,M_.exo_names,'exact');
if ~isempty(iexo),
y0=teff(T(iendo,iexo+nspred,:),kn,istable);
if ~isempty(y0),
if mod(iplo,9)==0,
ifig=ifig+1;
figure('name',[namendo,' vs. shocks ',int2str(ifig)]),
iplo=0;
end
iplo=iplo+1;
js=js+1;
subplot(3,3,iplo),
[SAmeas, SAMorris] = Morris_Measure_Groups(estim_params_.np, lpmat, y0,nliv);
SAM = squeeze(SAMorris(:,1));
SA(:,js)=SAM./(max(SAM)+eps);
[saso, iso] = sort(-SA(:,js));
bar(SA(iso(1:10),js))
%set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:10,
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title([namendo,' vs. ',namexo],'interpreter','none')
if iplo==9,
saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)]);
close(gcf)
end
end
end
end
if iplo<9 & iplo>0 & ifig,
saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]);
close(gcf)
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(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
if ~isempty(ilagendo),
y0=teff(T(iendo,ilagendo,:),kn,istable);
if ~isempty(y0),
if mod(iplo,9)==0,
ifig=ifig+1;
figure('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(estim_params_.np, lpmat, y0,nliv);
SAM = squeeze(SAMorris(:,1));
SA(:,js)=SAM./(max(SAM)+eps);
[saso, iso] = sort(-SA(:,js));
bar(SA(iso(1:10),js))
%set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:10,
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
if iplo==9,
saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
close(gcf)
end
end
end
end
if iplo<9 & iplo>0 & ifig,
saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
close(gcf)
end
end
figure,
%bar(SA)
boxplot(SA','whis',10,'symbol','r.')
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
xlabel(' ')
ylabel('Elementary Effects')
title('Reduced form screening')
saveas(gcf,[dirname,'\',M_.fname,'_redform_screen'])
eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_screen']);
eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_screen']);

15
matlab/smirnov.m
View File

@ -2,7 +2,20 @@ function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
% Smirnov test for 2 distributions
% [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
%
% Copyright (C) 2005 Marco Ratto
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%

273
matlab/stab_map_.m
View File

@ -1,14 +1,11 @@
function x0 = stab_map_(Nsam, fload, ksstat, alpha2, prepSA, pprior, ilptau, OutputDirectoryName)
function x0 = stab_map_(OutputDirectoryName)
%
% function x0 = stab_map_(Nsam, fload, alpha2, prepSA, pprior)
% function x0 = stab_map_(OutputDirectoryName)
%
% Mapping of stability regions in the prior ranges applying
% Monte Carlo filtering techniques.
%
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models
% I. Mapping stability, MIMEO, 2005.
%
% INPUTS
% INPUTS (from opt_gsa structure)
% Nsam = MC sample size
% fload = 0 to run new MC; 1 to load prevoiusly generated analysis
% alpha2 = significance level for bivariate sensitivity analysis
@ -34,100 +31,120 @@ function x0 = stab_map_(Nsam, fload, ksstat, alpha2, prepSA, pprior, ilptau, Out
% USES lptauSEQ,
% stab_map_1, stab_map_2
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% ALL COPIES MUST BE PROVIDED FREE OF CHARGE AND MUST INCLUDE THIS COPYRIGHT
% NOTICE.
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
opt_gsa=options_.opt_gsa;
Nsam = opt_gsa.Nsam;
fload = opt_gsa.load_stab;
ksstat = opt_gsa.ksstat;
alpha2 = opt_gsa.alpha2_stab;
prepSA = (opt_gsa.redform | opt_gsa.identification);
pprior = opt_gsa.pprior;
ilptau = opt_gsa.ilptau;
nliv = opt_gsa.morris_nliv;
ntra = opt_gsa.morris_ntra;
dr_ = oo_.dr;
if isfield(dr_,'ghx'),
%if isfield(dr_,'ghx'),
ys_ = oo_.dr.ys;
nspred = size(dr_.ghx,2);
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
%end
fname_ = M_.fname;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
lpmat0=[];
pshape = bayestopt_.pshape(nshock+1:end);
p1 = bayestopt_.p1(nshock+1:end);
p2 = bayestopt_.p2(nshock+1:end);
p3 = bayestopt_.p3(nshock+1:end);
p4 = bayestopt_.p4(nshock+1:end);
if nargin==0,
Nsam=2000; %2^13; %256;
end
if nargin<2,
fload=0;
end
if nargin<3,
ksstat=0.1;
end
if nargin<4,
alpha2=0.3;
end
if nargin<5,
prepSA=0;
end
if nargin<6,
pprior=1;
end
if nargin<7,
ilptau=1;
end
if nargin<8,
OutputDirectoryName='';
end
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
if fload==0 | nargin<2 | isempty(fload),
if prepSA
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2);
end
opt=options_;
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
options_.simul=0;
if fload==0,
% if prepSA
% T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2);
% end
if isfield(dr_,'ghx'),
egg=zeros(length(dr_.eigval),Nsam);
end
yys=zeros(length(dr_.ys),Nsam);
if estim_params_.np<52 & ilptau,
[lpmat] = lptauSEQ(Nsam,estim_params_.np);
if estim_params_.np>30
if opt_gsa.morris
if opt_gsa.morris == 1
[lpmat, OutFact] = Sampling_Function_2(nliv, estim_params_.np, ntra, ones(estim_params_.np, 1), zeros(estim_params_.np,1), []);
lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
Nsam=size(lpmat,1);
elseif opt_gsa.morris==2
lpmat = prep_ide(Nsam,estim_params_.np,5);
Nsam=size(lpmat,1);
end
else
if estim_params_.np<52 & ilptau>0,
[lpmat] = lptauSEQ(Nsam,estim_params_.np); % lptau
if estim_params_.np>30 | ilptau==2, % scrambled lptau
for j=1:estim_params_.np,
lpmat(:,j)=lpmat(randperm(Nsam),j);
end
end
else
else ilptau==0
%[lpmat] = rand(Nsam,estim_params_.np);
for j=1:estim_params_.np,
lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
end
end
end
prior_draw_gsa(1);
if pprior,
for j=1:nshock,
lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
lpmat0(:,j)=lpmat0(:,j).*(bayestopt_.ub(j)-bayestopt_.lb(j))+bayestopt_.lb(j);
if opt_gsa.prior_range
lpmat0(:,j)=lpmat0(:,j).*(bayestopt_.ub(j)-bayestopt_.lb(j))+bayestopt_.lb(j);
end
end
for j=1:estim_params_.np,
lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
if opt_gsa.prior_range
for j=1:estim_params_.np,
lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
end
else
xx=prior_draw_gsa(0,[lpmat0 lpmat]);
lpmat0=xx(:,1:nshock);
lpmat=xx(:,nshock+1:end);
clear xx;
end
else
% for j=1:nshock,
@ -183,7 +200,24 @@ if fload==0 | nargin<2 | isempty(fload),
iwrong=zeros(1,Nsam);
for j=1:Nsam,
M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
stoch_simul([]);
%try stoch_simul([]);
try
[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
bayestopt_.restrict_columns,...
bayestopt_.restrict_aux);
if ~exist('T')
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
end
catch
if isfield(oo_.dr,'eigval'),
oo_.dr=rmfield(oo_.dr,'eigval');
end
if isfield(oo_.dr,'ghx'),
oo_.dr=rmfield(oo_.dr,'ghx');
end
disp('No solution could be found'),
end
dr_ = oo_.dr;
if isfield(dr_,'ghx'),
egg(:,j) = sort(dr_.eigval);
@ -191,9 +225,13 @@ if fload==0 | nargin<2 | isempty(fload),
if prepSA
jstab=jstab+1;
T(:,:,jstab) = [dr_.ghx dr_.ghu];
[A,B] = ghx2transition(squeeze(T(:,:,jstab)), ...
bayestopt_.restrict_var_list, ...
bayestopt_.restrict_columns, ...
bayestopt_.restrict_aux);
end
if ~exist('nspred'),
nspred = size(dr_.ghx,2);
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
@ -212,7 +250,7 @@ if fload==0 | nargin<2 | isempty(fload),
end
else
if exist('egg'),
egg(:,j)=ones(size(egg,1),1).*1.1;
egg(:,j)=ones(size(egg,1),1).*NaN;
end
iwrong(j)=j;
end
@ -265,20 +303,31 @@ if fload==0 | nargin<2 | isempty(fload),
% end
% end
% iunstable=iunstable(find(iunstable)); % unstable params
bkpprior.pshape=bayestopt_.pshape;
bkpprior.p1=bayestopt_.p1;
bkpprior.p2=bayestopt_.p2;
bkpprior.p3=bayestopt_.p3;
bkpprior.p4=bayestopt_.p4;
if pprior,
if ~prepSA
save([OutputDirectoryName '\' fname_ '_prior'],'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
save([OutputDirectoryName '\' fname_ '_prior'], ...
'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
'egg','yys','nspred','nboth','nfwrd')
else
save([OutputDirectoryName '\' fname_ '_prior'],'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','T','nspred','nboth','nfwrd')
save([OutputDirectoryName '\' fname_ '_prior'], ...
'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
'egg','yys','T','nspred','nboth','nfwrd')
end
else
if ~prepSA
save([OutputDirectoryName '\' fname_ '_mc'], ...
'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
'egg','yys','nspred','nboth','nfwrd')
else
save([OutputDirectoryName '\' fname_ '_mc'], ...
'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','T','nspred','nboth','nfwrd')
'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
'egg','yys','T','nspred','nboth','nfwrd')
end
end
else
@ -289,6 +338,7 @@ else
end
load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
Nsam = size(lpmat,1);
if prepSA & isempty(strmatch('T',who('-file', filetoload),'exact')),
h = waitbar(0,'Please wait...');
@ -297,15 +347,21 @@ else
options_.irf=0;
options_.noprint=1;
stoch_simul([]);
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
%T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
ntrans=length(istable);
for j=1:ntrans,
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
stoch_simul([]);
%stoch_simul([]);
[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
bayestopt_.restrict_columns,...
bayestopt_.restrict_aux);
if ~exist('T')
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
end
dr_ = oo_.dr;
T(:,:,j) = [dr_.ghx dr_.ghu];
T(:,:,j) = [dr_.ghx dr_.ghu];
if ~exist('nspred')
nspred = size(dr_.ghx,2);
nspred = dr_.nspred; %size(dr_.ghx,2);
nboth = dr_.nboth;
nfwrd = dr_.nfwrd;
end
@ -316,6 +372,8 @@ else
end
close(h)
save(filetoload,'T','-append')
elseif prepSA
load(filetoload,'T')
end
end
@ -340,37 +398,54 @@ delete([OutputDirectoryName,'\',fname_,'_',aunstname,'_corr_*.*']);
delete([OutputDirectoryName,'\',fname_,'_',aindname,'_corr_*.*']);
if length(iunstable)>0 & length(iunstable)<Nsam,
disp([num2str(length(istable)/Nsam*100),'\% of the prior support is stable.'])
disp([num2str(length(istable)/Nsam*100,3),'\% of the prior support is stable.'])
disp([num2str( (length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100),'\% of the prior support is unstable.'])
if ~isempty(iindeterm),
disp([num2str(length(iindeterm)/Nsam*100),'\% of the prior support gives indeterminacy.'])
disp([num2str(length(iindeterm)/Nsam*100,3),'\% of the prior support gives indeterminacy.'])
end
if ~isempty(iwrong),
disp(' ');
disp(['For ',num2str(length(iwrong)/Nsam*100),'\% of the prior support dynare could not find a solution.'])
disp(['For ',num2str(length(iwrong)/Nsam*100,3),'\% of the prior support dynare could not find a solution.'])
end
disp(' ');
% Blanchard Kahn
[proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
indstab=find(dproba>ksstat);
disp('The following parameters mostly drive acceptable behaviour')
disp(M_.param_names(estim_params_.param_vals(indstab,1),:))
disp('Smirnov statistics in driving acceptable behaviour')
for j=1:estim_params_.np,
disp([M_.param_names(estim_params_.param_vals(j,1),:),' d-stat = ', num2str(dproba(j),3)])
end
disp(' ');
if ~isempty(indstab)
stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName);
end
ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
if ~isempty(iindeterm),
ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'],0);
indindet=find(dproba>ksstat);
disp('The following parameters mostly drive indeterminacy')
disp(M_.param_names(estim_params_.param_vals(indindet,1),:))
stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'], 1, indindet, OutputDirectoryName);
if ~isempty(ixun),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'],0);
indunst=find(dproba>ksstat);
disp('The following parameters mostly drive instability')
disp(M_.param_names(estim_params_.param_vals(indunst,1),:))
stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'], 1, indunst, OutputDirectoryName);
disp('Smirnov statistics in driving indeterminacy')
for j=1:estim_params_.np,
disp([M_.param_names(estim_params_.param_vals(j,1),:),' d-stat = ', num2str(dproba(j),3)])
end
disp(' ');
if ~isempty(indindet)
stab_map_1(lpmat, istable, iindeterm, [aname, '_indet'], 1, indindet, OutputDirectoryName);
end
end
if ~isempty(ixun),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'],0);
indunst=find(dproba>ksstat);
disp('Smirnov statistics in driving instability')
for j=1:estim_params_.np,
disp([M_.param_names(estim_params_.param_vals(j,1),:),' d-stat = ', num2str(dproba(j),3)])
end
disp(' ');
if ~isempty(indunst)
stab_map_1(lpmat, istable, ixun, [aname, '_unst'], 1, indunst, OutputDirectoryName);
end
end
disp(' ')
disp('Starting bivariate analysis:')
@ -378,14 +453,16 @@ if length(iunstable)>0 & length(iunstable)<Nsam,
c00=tril(c0,-1);
stab_map_2(lpmat(istable,:),alpha2, asname, OutputDirectoryName);
stab_map_2(lpmat(iunstable,:),alpha2, auname, OutputDirectoryName);
if ~isempty(iindeterm),
stab_map_2(lpmat(iindeterm,:),alpha2, aindname, OutputDirectoryName);
if ~isempty(ixun),
stab_map_2(lpmat(ixun,:),alpha2, aunstname, OutputDirectoryName);
end
if length(iunstable)>3,
stab_map_2(lpmat(iunstable,:),alpha2, auname, OutputDirectoryName);
end
if length(iindeterm)>3,
stab_map_2(lpmat(iindeterm,:),alpha2, aindname, OutputDirectoryName);
end
if length(ixun)>3,
stab_map_2(lpmat(ixun,:),alpha2, aunstname, OutputDirectoryName);
end
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
if istable(end)~=Nsam
@ -403,3 +480,17 @@ else
end
end
options_.periods=opt.periods;
if isfield(opt,'nomoments'),
options_.nomoments=opt.nomoments;
end
options_.irf=opt.irf;
options_.noprint=opt.noprint;
if isfield(opt,'simul'),
options_.simul=opt.simul;
end

20
matlab/stab_map_1.m
View File

@ -14,6 +14,22 @@ function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, i
% Plots: dotted lines for BEHAVIOURAL
% solid lines for NON BEHAVIOURAL
% USES smirnov
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global estim_params_ bayestopt_ M_ options_
@ -33,7 +49,7 @@ nshock = nshock + estim_params_.ncn;
npar=size(lpmat,2);
ishock= npar>estim_params_.np;
if nargin<6,
if nargin<6 | isempty(ipar),
ipar=[1:npar];
end
nparplot=length(ipar);
@ -68,4 +84,4 @@ for i=1:ceil(nparplot/12),
eval(['print -dpdf ' dirname '\' fname_ '_' aname '_SA_' int2str(i)]);
if options_.nograph, close(gcf), end
end
end
end

23
matlab/stab_map_2.m
View File

@ -1,17 +1,22 @@
%function stab_map_2(x,alpha2,istab,fnam)
function stab_map_2(x,alpha2,fnam, dirname)
% function stab_map_2(x,alpha2,fnam)
% function stab_map_2(x,alpha2,fnam, dirname)
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_

34
matlab/teff.m Normal file
View File

@ -0,0 +1,34 @@
function [yt, j0]=teff(T,Nsam,istable)
%
% [yt, j0]=teff(T,Nsam,istable)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
tmax=max(T,[],3);
tmin=min(T,[],3);
[ir, ic]=(find( (tmax-tmin)>1.e-8));
j0 = length(ir);
yt=zeros(Nsam, j0);
for j=1:j0,
y0=squeeze(T(ir(j),ic(j),:));
%y1=ones(size(lpmat,1),1)*NaN;
y1=ones(Nsam,1)*NaN;
y1(istable,1)=y0;
yt(:,j)=y1;
end
%clear y0 y1;

25
matlab/trank.m Normal file
View File

@ -0,0 +1,25 @@
function yr = trank(y);
% yr = trank(y);
% yr is the rank transformation of y
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
[nr, nc] = size(y);
for j=1:nc,
[dum, is]=sort(y(:,j));
yr(is,j)=[1:nr]'./nr;
end