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
  • 4.3
  • 4.4
  • 4.5
  • 4.6
  • 5.x
  • 6.x
  • asm
  • aux_func
  • clang+openmp
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • dragonfly
  • dynare_minreal
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exo_steady_state
  • gpm-optimal-policy
  • julia
  • madysson
  • master
  • mex-GetPowerDeriv
  • penalty
  • separateM_
  • slice
  • sphinx-doc-experimental
  • static_aux_vars
  • time-varying-information-set
  • various_fixes
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
  • 4.5.7
  • 4.6-beta1
  • 4.6.0
  • 4.6.0-rc1
  • 4.6.0-rc2
  • 4.6.1
  • 4.6.2
  • 4.6.3
  • 4.6.4
  • 4.7-beta1
  • 4.7-beta2
  • 4.7-beta3
  • 5.0
  • 5.0-rc1
  • 5.1
  • 5.2
  • 5.3
  • 5.4
  • 5.5
  • 6-beta1
  • 6-beta2
  • 6.0
  • 6.1
  • 6.2
  • 6.3
  • 6.4
91 results

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
  • bayesian_irf_matching
  • cet-octave-one-more-time
  • docker-meson
  • identificationME
  • irf-matching-example
  • irf_match_perfect_foresight
  • main
  • mom_fname
  • new-desktop-matlab
  • optimizers
  • pskf
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
46 results
Show changes
Showing
with 151 additions and 2162 deletions
......@@ -68,12 +68,10 @@ else
end
%% Call Dynamic Function
[junk, g1] = feval([M_.fname '.dynamic'], ...
ones(max(max(M_.lead_lag_incidence)), 1), ...
ones(1, M_.exo_nbr), ...
M_.params, ...
zeros(M_.endo_nbr, 1), ...
1);
g1 = feval([M_.fname '.sparse.dynamic_g1'], ones(3*M_.endo_nbr, 1), ones(1, M_.exo_nbr), ...
M_.params, zeros(M_.endo_nbr, 1), ...
M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
M_.dynamic_g1_sparse_colptr);
% Choose rows of Jacobian based on equation tags
ntags = length(M_.(model_type).(model_name).eqtags);
......@@ -88,8 +86,7 @@ g1 = -1 * g1(g1rows, :);
% Check for leads
if rows(M_.lead_lag_incidence) == 3
idxs = M_.lead_lag_incidence(3, M_.lead_lag_incidence(3, :) ~= 0);
assert(~any(any(g1(g1rows, idxs))), ...
assert(~any(any(g1(g1rows, 2*M_.endo_nbr+(1:M_.endo_nbr)))), ...
['You cannot have leads in the equations specified by ' strjoin(M_.(model_type).(model_name).eqtags, ',')]);
end
......@@ -152,19 +149,19 @@ for i = 1:length(lhs)
for j = 1:length(rhsvars{i}.vars)
var = rhsvars{i}.vars(j);
if rhsvars{i}.lags(j) == -1
g1col = M_.lead_lag_incidence(1, var);
g1col = var;
else
g1col = M_.lead_lag_incidence(2, var);
g1col = var + M_.endo_nbr;
end
if g1col ~= 0 && any(g1(:, g1col))
if any(g1(:, g1col))
if rhsvars{i}.arRhsIdxs(j) > 0
% Fill AR
[lag, ndiffs] = findLagForVar(var, -rhsvars{i}.lags(j), 0, lhs);
lag = findLagForVar(var, -rhsvars{i}.lags(j), 0, lhs);
oo_.(model_type).(model_name).ar(i, rhsvars{i}.arRhsIdxs(j), lag) = ...
oo_.(model_type).(model_name).ar(i, rhsvars{i}.arRhsIdxs(j), lag) + g1(i, g1col);
elseif rhsvars{i}.ecRhsIdxs(j) > 0
% Fill EC
[lag, ndiffs] = findLagForVar(var, -rhsvars{i}.lags(j), 0, ecRhsVars);
lag = findLagForVar(var, -rhsvars{i}.lags(j), 0, ecRhsVars);
if lag==1
if size(oo_.(model_type).(model_name).ec, 3) < lag
oo_.(model_type).(model_name).ec(i, rhsvars{i}.ecRhsIdxs(j), lag) = 0;
......
function [state_u,state_n] = get_dynare_random_generator_state()
% Get state of Matlab/Octave random generator depending on matlab
% (octave) version.
% In older versions, Matlab kept one generator for uniformly distributed numbers and
% Get state of MATLAB/Octave random generator depending on MATLAB
% (Octave) version.
% In older versions, MATLAB kept one generator for uniformly distributed numbers and
% one for normally distributed numbers.
% For backward compatibility, we return two vectors, but, in recent
% versions of Matlab and in Octave, we return two identical vectors.
% versions of MATLAB and in Octave, we return two identical vectors.
% Copyright © 2010-2020 Dynare Team
%
......
function eqnumber = get_equation_number_by_tag(eqname, M_)
% eqnumber = get_equation_number_by_tag(eqname, M_)
% Translates an equation name into an equation number.
%
% INPUTS
% - eqname [char] 1×n array, name of the equation.
% - DynareModel [struct] Structure describing the model, aka M_.
% - M_ [struct] Structure describing the model
%
% OUTPUTS
% - eqnumber [integer] Equation number.
% Copyright © 2018-2022 Dynare Team
% Copyright © 2018-2023 Dynare Team
%
% This file is part of Dynare.
%
......
function message = get_error_message(info, DynareOptions)
function message = get_error_message(info, options_)
% Returns error messages
%
% INPUTS
% info [double] vector returned by resol.m
% DynareOptions [structure] --> options_
% options_ [structure] Dynare options
% OUTPUTS
% message [string] corresponding error message
%
% SPECIAL REQUIREMENTS
% none
% Copyright © 2005-2020 Dynare Team
% Copyright © 2005-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -43,7 +43,7 @@ switch info(1)
case 6
message = 'The Jacobian matrix evaluated at the steady state contains elements that are not real or are infinite.';
case 7
message = sprintf('One of the eigenvalues is close to 0/0 (the absolute value of numerator and denominator is smaller than %5.4f!\n If you believe that the model has a unique solution you can try to reduce the value of qz_zero_threshold.',DynareOptions.qz_zero_threshold);
message = sprintf('One of the eigenvalues is close to 0/0 (the absolute value of numerator and denominator is smaller than %5.4f!\n If you believe that the model has a unique solution you can try to reduce the value of qz_zero_threshold.',options_.qz_zero_threshold);
case 8
if size(info,2)>=2 && info(2)~=0
global M_;
......@@ -66,7 +66,7 @@ switch info(1)
case 19
message = 'The steadystate file did not compute the steady state';
case 20
if DynareOptions.linear
if options_.linear
message = sprintf('Impossible to find the steady state (the sum of squared residuals of the static equations is %5.4f). Either the model doesn''t have a steady state or there are an infinity of steady states. Check whether your model is truly linear or whether there is a mistake in linearization.', info(2));
else
message = sprintf('Impossible to find the steady state (the sum of squared residuals of the static equations is %5.4f). Either the model doesn''t have a steady state, there are an infinity of steady states, or the guess values are too far from the solution', info(2));
......@@ -85,21 +85,23 @@ switch info(1)
message = 'The loglinearization of the model cannot be performed, because the steady state is not strictly positive.';
case 30
message = 'Ergodic variance can''t be computed.';
case 40
message = 'Prior density is Inf';
case 41
message = 'one (many) parameter(s) do(es) not satisfy the lower bound';
case 42
message = 'one (many) parameter(s) do(es) not satisfy the upper bound';
case 43
message = 'Covariance matrix of structural shocks is not positive definite';
case 44 %DsgeLikelihood_hh / dsge_likelihood
case 44
message = 'The covariance matrix of the measurement errors is not positive definite.';
case 45 %DsgeLikelihood_hh / dsge_likelihood
message = 'Likelihood is not a number (NaN) or a complex number';
case 46 %DsgeLikelihood_hh / dsge_likelihood
case 45
message = 'Likelihood is not a number (NaN)';
case 46
message = 'Likelihood is a complex number';
case 47 %DsgeLikelihood_hh / dsge_likelihood
case 47
message = 'Prior density is not a number (NaN)';
case 48 %DsgeLikelihood_hh / dsge_likelihood
case 48
message = 'Prior density is a complex number';
case 49
message = 'The model violates one (many) endogenous prior restriction(s)';
......@@ -107,13 +109,13 @@ switch info(1)
message = 'Likelihood is Inf';
case 51
message = sprintf('\n The dsge_prior_weight is dsge_var=%5.4f, but must be at least %5.4f for the prior to be proper.\n You are estimating a DSGE-VAR model, but the value of the dsge prior weight is too low!', info(2), info(3));
case 52 %dsge_var_likelihood
case 52
message = 'You are estimating a DSGE-VAR model, but the implied covariance matrix of the VAR''s innovations, based on artificial and actual sample is not positive definite!';
case 53 %dsge_var_likelihood
case 53
message = 'You are estimating a DSGE-VAR model, but the implied covariance matrix of the VAR''s innovations, based on the artificial sample, is not positive definite!';
case 55
message = 'Fast Kalman filter only works with stationary models [lik_init=1] or stationary observables for non-stationary models [lik_init=3]';
case 61 %Discretionary policy
case 61
message = 'Discretionary policy: maximum number of iterations has been reached. Procedure failed.';
case 62
message = 'Discretionary policy: some eigenvalues greater than options_.qz_criterium. Model potentially unstable.';
......@@ -129,7 +131,6 @@ switch info(1)
message = 'Calibrated covariance of the structural errors implies correlation larger than +-1.';
case 72
message = 'Calibrated covariance of the measurement errors implies correlation larger than +-1.';
% Aim Code Conversions by convertAimCodeToInfo.m
case 81
message = ['Ramsey: The solution to the static first order conditions for optimal policy could not be found. Either the model' ...
' doesn''t have a steady state, there are an infinity of steady states, ' ...
......@@ -163,15 +164,33 @@ switch info(1)
case 162
message = 'Aim: too many numeric shiftrights';
case 163
message = 'Aim: A is NAN or INF.';
message = 'Aim: A is NaN or Inf.';
case 164
message = 'Aim: Problem in SPEIG.';
case 174
message = 'OSR: the loss is Inf';
case 175
message = 'OSR: the loss is NaN';
case 176
message = 'OSR: the loss is a complex number';
case 177
message = 'method_of_moments: the sum of squared residuals is Inf';
case 178
message = 'method_of_moments: the sum of squared residuals is NaN';
case 179
message = 'method_of_moments: the sum of squared residuals is a complex number';
case 180
message = 'SMM: simulation resulted in NaN/Inf. You may need to enable pruning.';
case 181
message = 'IRF Matching: simulated IRFs were explosive. Either reduce the shock size, use pruning, or set the approximation order to 1.';
case 182
message = 'IRF Matching: transformations in irf_matching_file returned with an error.';
case 201
message = 'Particle Filter: Initial covariance of the states is not positive definite. Try a different nonlinear_filter_initialization';
case 202
message = 'Particle Filter: Initial covariance of the states based on simulation resulted in NaN/Inf. Use pruning or try a different nonlinear_filter_initialization';
case 300
message = 'IVF: The likelihood is a complex number.';
case 301
message = 'IVF: The likelihood is Inf.';
case 302
......@@ -183,19 +202,29 @@ switch info(1)
case 305
message = 'IVF: The returned shocks are bigger than 1e8.';
case 310
message = 'Occbin: Simulation terminated with periodic solution (no convergence).';
message = 'OccBin: Simulation terminated with periodic solution (no convergence).';
case 311
message = 'Occbin: Simulation did not converge, increase maxit or check_ahead_periods.';
message = 'OccBin: Simulation did not converge, increase maxit or check_ahead_periods.';
case 312
message = 'Occbin: Constraint(s) are binding at the end of the sample.';
message = 'OccBin: Constraint(s) are binding at the end of the sample.';
case 313
message = 'Occbin: Simulation did not converge -- infinite loop of guess regimes';
message = 'OccBin: Simulation did not converge -- infinite loop of guess regimes';
case 320
message = 'Piecewise linear Kalman filter: There was a problem in obtaining the likelihood.';
case 321
message = 'Occbin: there was a problem in running the smoother. Simulation within smoother failed.';
message = 'OccBin: there was a problem in running the smoother. Simulation within smoother failed.';
case 322
message = 'Occbin: smoother did not converge.';
message = 'OccBin: smoother did not converge.';
case 323
message = 'Piecewise linear Kalman filter: the likelihood is NaN.';
case 324
message = 'Piecewise linear Kalman filter: updated state vector is NaN.';
case 325
message = 'Piecewise linear Kalman filter: filter covariance NaN.';
case 330
message = 'Piecewise linear Kalman filter: update step did not reach a fixed point (periodic loop).';
case 331
message = 'Piecewise linear Kalman filter: update step did not reach a fixed point (max number of iterations reached).';
case 401
message = 'Cycle reduction reached the iteration limit. Try increasing maxit.';
case 402
......
......@@ -28,7 +28,7 @@ function ext = get_file_extension(file)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
[dir, fname, ext] = fileparts(file);
[~, ~, ext] = fileparts(file);
if ~isempty(ext)
% Removes the leading dot.
......
function [lhs, rhs, json] = get_lhs_and_rhs(eqname, DynareModel, original, json)
function [lhs, rhs, json] = get_lhs_and_rhs(eqname, M_, original, json)
% [lhs, rhs, json] = get_lhs_and_rhs(eqname, M_, original, json)
% Returns the left and right handsides of an equation.
%
% INPUTS
% - eqname [char] Name of the equation.
% - DynareModel [struct] Structure describing the current model (M_).
% - M_ [struct] Structure describing the current model.
% - original [logical] fetch equation in modfile-original.json or modfile.json
% - json [char] content of the JSON file
%
......@@ -23,7 +23,7 @@ function [lhs, rhs, json] = get_lhs_and_rhs(eqname, DynareModel, original, json)
% [name='Phillips curve']
% pi = beta*pi(1) + slope*y + lam;
% Copyright © 2018-2020 Dynare Team
% Copyright © 2018-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -47,9 +47,9 @@ end
% Load JSON file if nargin<4
if nargin<4
if original
json = loadjson_([DynareModel.fname filesep() 'model' filesep() 'json' filesep() 'modfile-original.json']);
json = loadjson_([M_.fname filesep() 'model' filesep() 'json' filesep() 'modfile-original.json']);
else
json = loadjson_([DynareModel.fname filesep() 'model' filesep() 'json' filesep() 'modfile.json']);
json = loadjson_([M_.fname filesep() 'model' filesep() 'json' filesep() 'modfile.json']);
end
end
......
......@@ -36,24 +36,24 @@ function indx = get_new_or_existing_ei_index(substructure_name, name1, name2)
global estimation_info
if eval(['isempty(estimation_info.' substructure_name ')'])
if isempty(estimation_info.(substructure_name))
indx = 1;
return
end
if isempty(name2) % parameter or std() statement
indx = eval(['find(strcmp(name1, estimation_info.' substructure_name ') == 1)']);
indx = find(strcmp(name1, estimation_info.(substructure_name)));
else % corr statement
indx = eval(['find(strcmp([''' name1 ':' name2 '''], estimation_info.' substructure_name ') == 1)']);
indx = find(strcmp(sprintf('%s:%s', name1, name2), estimation_info.(substructure_name)));
if isempty(indx)
indx = eval(['find(strcmp([''' name2 ':' name1 '''], estimation_info.' substructure_name ') == 1)']);
indx = find(strcmp(sprintf('%s:%s', name2, name1), estimation_info.(substructure_name)));
end
end
if isempty(indx)
indx = eval(['size(estimation_info.' substructure_name ', 2) + 1']);
indx = size(estimation_info.(substructure_name), 2) + 1;
end
if size(indx, 2) > 1
error(['Error: ' name1 ' ' name2 'found more than once in estimation_info.subsamples_index']);
error('Error: %s %s found more than once in estimation_info.subsamples_index', name1, name2);
end
......@@ -2,7 +2,7 @@ function mexpath = get_path_to_mex_files(dynareroot)
% Returns a cell array containing one or several directory paths
% which should contain the MEX files.
% Copyright © 2015-2023 Dynare Team
% Copyright © 2015-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -47,34 +47,20 @@ else
end
% Add win64 specific paths for Dynare Windows package
if strcmp(computer, 'PCWIN64')
if matlab_ver_less_than('9.4')
tmp = [dynareroot '../mex/matlab/win64-8.3-9.3/'];
tmp = [dynareroot '../mex/matlab/win64-9.8-25.1/'];
if exist(tmp, 'dir')
mexpath = tmp;
end
else
tmp = [dynareroot '../mex/matlab/win64-9.4-23.2/'];
if exist(tmp, 'dir')
mexpath = tmp;
end
end
end
% Add macOS paths for Dynare Mac package
if strcmp(computer, 'MACI64')
if matlab_ver_less_than('9.4')
tmp = [dynareroot '../mex/matlab/maci64-8.3-9.3/'];
tmp = [dynareroot '../mex/matlab/maci64-9.8-25.1/'];
if exist(tmp, 'dir')
mexpath = tmp;
end
else
tmp = [dynareroot '../mex/matlab/maci64-9.4-23.2/'];
if exist(tmp, 'dir')
mexpath = tmp;
end
end
end
if strcmp(computer, 'MACA64')
tmp = [dynareroot '../mex/matlab/maca64-23.2/'];
tmp = [dynareroot '../mex/matlab/maca64-23.2-25.1/'];
if exist(tmp, 'dir')
mexpath = tmp;
end
......
function [nam, texnam] = get_the_name(k, TeX, M_, estim_params_, options_)
%@info:
%! @deftypefn {Function File} {[@var{nam},@var{texnam}] =} get_the_name (@var{k},@var{TeX},@var{M_},@var{estim_params_},@var{options_})
%! @anchor{get_the_name}
%! @sp 1
%! Returns the name of the estimated parameter number @var{k}, following the internal ordering of the estimated parameters.
%! @sp 2
%! @strong{Inputs}
%! @sp 1
%! @table @ @var
%! @item k
%! Integer scalar, parameter number.
%! @item TeX
%! Integer scalar, if @var{TeX}==0 then @var{texnam} is not returned (empty matrix).
%! @item M_
%! Matlab's structure describing the model (initialized by @code{dynare}).
%! @item estim_params_
%! Matlab's structure describing the estimated parameters (initialized by @code{dynare}).
%! @item options_
%! Matlab's structure describing the options (initialized by @code{dynare}).
%! @end table
%! @sp 2
%! @strong{Outputs}
%! @sp 1
%! @table @ @var
%! @item nam
%! String, internal name of the variable
%! @item texnam
%! String, TeX name of the same variable (if defined in the mod file).
%! @end table
%! @sp 2
%! @strong{This function is called by:}
%! @sp 1
%! @ref{get_prior_info}, @ref{mcmc_diagnostics}, @ref{mode_check}, @ref{PlotPosteriorDistributions}, @ref{plot_priors}
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! None.
%! @end deftypefn
%@eod:
function [nam, texnam] = get_the_name(k, TeX, M_, estim_params_, varobs)
% [nam, texnam] = get_the_name(k, TeX, M_, estim_params_, varobs)
% Returns name of estimated parameter number k, following the internal ordering of
% the estimated parameters.
% Inputs:
% - k [integer] parameter number.
% - TeX [bool] if false, texnam is not returned (empty matrix)
% - M_ [structure] model
% - estim_params_ [structure] describing the estimated parameters
% - varobs [cell] name of observed variables
%
% Outputs
% - nam [char] internal name of the variable
% - texnam [char] TeX name of the same variable (if defined in the mod file)
%
% This function is called by:
% get_prior_info, mcmc_diagnostics, mode_check, PlotPosteriorDistributions, plot_priors
%
% This function calls:
% None.
%
% Copyright © 2004-2018 Dynare Team
% Copyright © 2004-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -57,7 +37,6 @@ function [nam, texnam] = get_the_name(k, TeX, M_, estim_params_, options_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
nam = [];
texnam = [];
nvx = estim_params_.nvx;
......@@ -73,7 +52,7 @@ if k <= nvx
texnam = sprintf('$ \\sigma_{%s} $', tname);
end
elseif k <= (nvx+nvn)
vname = options_.varobs{estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1)};
vname = varobs{estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1)};
nam = sprintf('SE_EOBS_%s', vname);
if TeX
tname = M_.endo_names_tex{estim_params_.var_endo(k-estim_params_.nvx,1)};
......
function [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, DynareModel)
function [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_)
% [pnames, enames, xnames, pid, eid, xid] = get_variables_and_parameters_in_equation(lhs, rhs, M_)
% Returns the lists of parameters, endogenous variables and exogenous variables in an equation.
%
% INPUTS
% - lhs [string] Left hand side of an equation.
% - rhs [string] Right hand side of an equation.
% - DynareModel [struct] Structure describing the current model (M_).
% - M_ [struct] Structure describing the current model.
%
% OUTPUTS
% - pnames [cell] Cell of row char arrays (p elements), names of the parameters.
......@@ -39,33 +39,33 @@ rhs_ = get_variables_and_parameters_in_expression(rhs);
lhs_ = get_variables_and_parameters_in_expression(lhs);
% Get list of parameters.
pnames = DynareModel.param_names;
pnames = M_.param_names;
pnames = intersect([rhs_, lhs_], pnames);
if nargout>1
% Get list of endogenous variables.
enames = DynareModel.endo_names;
enames = M_.endo_names;
enames = intersect([rhs_, lhs_], enames);
if nargout>2
% Get list of exogenous variables
xnames = DynareModel.exo_names;
xnames = M_.exo_names;
xnames = intersect([rhs_,lhs_], xnames);
if nargout>3
% Returns vector of indices for parameters endogenous and exogenous variables if required.
p = length(pnames);
pid = zeros(p, 1);
for i = 1:p
pid(i) = find(strcmp(pnames{i}, DynareModel.param_names));
pid(i) = find(strcmp(pnames{i}, M_.param_names));
end
p = length(enames);
eid = zeros(p, 1);
for i = 1:p
eid(i) = find(strcmp(enames{i}, DynareModel.endo_names));
eid(i) = find(strcmp(enames{i}, M_.endo_names));
end
p = length(xnames);
xid = zeros(p, 1);
for i = 1:p
xid(i) = find(strcmp(xnames{i}, DynareModel.exo_names));
xid(i) = find(strcmp(xnames{i}, M_.exo_names));
end
end
end
......
......@@ -3,7 +3,7 @@ function objects = get_variables_and_parameters_in_expression(expr)
% Returns the variables and parameters appearing in an expression.
%
% INPUTS
% - expr [char] 1×m char array, dynare model expression (typically RHS or LHS of an equation).
% - expr [char] 1×m char array, Dynare model expression (typically RHS or LHS of an equation).
%
% OUTPUTS
% - objects [cell] cell of row char arrays, names of the variables and parameters in expr.
......
......@@ -11,7 +11,7 @@ function global_initialization()
% SPECIAL REQUIREMENTS
% none
% Copyright © 2003-2018 Dynare Team
% Copyright © 2003-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -28,7 +28,7 @@ function global_initialization()
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global oo_ M_ options_ estim_params_ bayestopt_ estimation_info ex0_ ys0_ dataset_ dataset_info
global oo_ M_ options_ estim_params_ bayestopt_ estimation_info dataset_ dataset_info
estim_params_ = [];
bayestopt_ = [];
dataset_=[];
......@@ -86,13 +86,13 @@ estimation_info.joint_parameter = {'index','domain','interval','mean','median','
oo_.exo_simul = [];
oo_.endo_simul = [];
ys0_ = [];
ex0_ = [];
oo_.dr = [];
oo_.exo_steady_state = [];
oo_.exo_det_steady_state = [];
oo_.exo_det_simul = [];
oo_.initval_series = dseries();
oo_.initial_steady_state = [];
oo_.initial_exo_steady_state = [];
oo_.gui.ran_estimation = false;
oo_.gui.ran_stoch_simul = false;
......@@ -136,13 +136,11 @@ else
options_ = default_option_values(M_);
end
% initialize persistent variables in priordens()
priordens([],[],[],[],[],[],1);
% initialize persistent variables in dyn_first_order_solver()
dyn_first_order_solver();
% Set dynare random generator and seed.
set_dynare_seed('default');
options_=set_dynare_seed_local_options(options_,'default');
% Load user configuration file.
if isfield(options_, 'global_init_file')
......
function gsa_plotmatrix(type,varargin)
% function gsa_plotmatrix(type,varargin)
% extended version of the standard MATLAB plotmatrix
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2011-2012 European Commission
% Copyright © 2011-2017 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/>.
global bayestopt_ options_ M_
RootDirectoryName = CheckPath('gsa',M_.dname);
if options_.opt_gsa.pprior
load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable','iunstable','iindeterm','iwrong')
else
load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable','iunstable','iindeterm','iwrong')
eval(['load ' options_.mode_file ' xparam1;']');
end
iexplosive = iunstable(~ismember(iunstable,[iindeterm;iwrong]));
switch type
case 'all'
x=[lpmat0 lpmat];
NumberOfDraws=size(x,1);
B=NumberOfDraws;
case 'stable'
x=[lpmat0(istable,:) lpmat(istable,:)];
NumberOfDraws=size(x,1);
B=NumberOfDraws;
case 'nosolution'
x=[lpmat0(iunstable,:) lpmat(iunstable,:)];
NumberOfDraws=size(x,1);
B=NumberOfDraws;
case 'unstable'
x=[lpmat0(iexplosive,:) lpmat(iexplosive,:)];
NumberOfDraws=size(x,1);
B=NumberOfDraws;
case 'indeterm'
x=[lpmat0(iindeterm,:) lpmat(iindeterm,:)];
NumberOfDraws=size(x,1);
B=NumberOfDraws;
case 'wrong'
x=[lpmat0(iwrong,:) lpmat(iwrong,:)];
NumberOfDraws=size(x,1);
B=NumberOfDraws;
end
if isempty(x)
disp('Empty parameter set!')
return
end
for j=1:length(varargin)
jcol(j)=strmatch(varargin{j},bayestopt_.name,'exact');
end
[H,AX,BigA,P,PAx]=plotmatrix(x(:,jcol));
for j=1:length(varargin)
% axes(AX(1,j)), title(varargin{j})
% axes(AX(j,1)), ylabel(varargin{j})
% set(AX(1,j),'title',varargin{j}),
set(get(AX(j,1),'ylabel'),'string',varargin{j})
set(get(AX(end,j),'xlabel'),'string',varargin{j})
end
if options_.opt_gsa.pprior==0
xparam1=xparam1(jcol);
for j=1:length(varargin)
for i=1:j-1
axes(AX(j,i))
hold on, plot(xparam1(i),xparam1(j),'*r')
end
for i=j+1:length(varargin)
axes(AX(j,i))
hold on, plot(xparam1(i),xparam1(j),'*r')
end
end
end
function map_ident_(OutputDirectoryName,opt_gsa)
% 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/>.
global bayestopt_ M_ options_ estim_params_ oo_
% opt_gsa = options_.opt_gsa;
fname_ = M_.fname;
nliv = opt_gsa.morris_nliv;
ntra = opt_gsa.morris_ntra;
itrans = opt_gsa.trans_ident;
np = estim_params_.np;
if opt_gsa.load_ident_files
gsa_flag=0;
else
gsa_flag=-2;
end
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','nspred','nboth','nfwrd')
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
% th moments
% options_.ar = min(3,options_.ar);
mss = yys(bayestopt_.mfys,:);
mss = teff(mss(:,istable),Nsam,istable);
yys = teff(yys(oo_.dr.order_var,istable),Nsam,istable);
if exist('T')
[vdec, cc, ac] = mc_moments(T, lpmatx, oo_.dr);
else
return
end
if opt_gsa.morris==2
pdraws = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
% [pdraws, TAU, GAM] = dynare_identification(options_.options_ident,[lpmatx lpmat(istable,:)]);
if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0
disp(['Sample check OK ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
clear pdraws;
end
% for j=1:length(istable), gas(:,j)=[vech(cc(:,:,j)); vec(ac(:,:,j))]; end
% if ~isempty(mss),
% gas = [mss(istable,:)'; gas];
% end
% if max(max(abs(GAM-gas)))<=1.e-8,
% disp(['Moments check OK ',num2str(max(max(abs(GAM-gas))))]),
clear GAM gas
% end
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)
myboxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
% boxplot(squeeze(vdec(:,j,:))','whis',10,'symbol','.r')
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)
text(ip,-4,deblank(options_.varobs(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
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,:)=stand_(squeeze(log(cc(j,j,:))))./2;
end
[vdec, j0, ir_vdec, ic_vdec] = teff(vdec,Nsam,istable);
[cc, j0, ir_cc, ic_cc] = teff(cc,Nsam,istable);
[ac, j0, ir_ac, ic_ac] = teff(ac,Nsam,istable);
[nr1, nc1, nnn] = size(T);
endo_nbr = M_.endo_nbr;
nstatic = M_.nstatic;
nspred = M_.nspred;
iv = (1:endo_nbr)';
ic = [ nstatic+(1:nspred) endo_nbr+(1:size(oo_.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);
% bayestopt_.restrict_var_list, ...
% bayestopt_.restrict_columns, ...
% bayestopt_.restrict_aux, M_.exo_nbr);
A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
% Sig(estim_params_.var_exo(:,1))=lpmatx(1,:).^2;
if ~isempty(lpmatx)
set_shocks_param(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);
% bayestopt_.restrict_var_list, ...
% bayestopt_.restrict_columns, ...
% bayestopt_.restrict_aux, M_.exo_nbr);
if ~isempty(lpmatx)
set_shocks_param(lpmatx(j,:));
end
A(:,:,j)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
end
clear T
clear lpmatx
[nr,nc,nn]=size(A);
io=bayestopt_.mf2;
% T1=A(io,1:nr,:);
% ino=find(~ismember([1:nr],io));
% T2=A(ino,1:nr,:);
R=A(:,nr+1:nc,:);
% [tadj, iff] = gsa_speed(A(1:nr,1:nr,:),R,io,0.5);
% [tadj, j0, ir_tadj, ic_tadj] = teff(tadj,Nsam,istable);
% [iff, j0, ir_if, ic_if] = teff(iff,Nsam,istable);
[yt, j0]=teff(A,Nsam,istable);
yt = [yys yt];
if opt_gsa.morris==2
% iii=find(std(yt(istable,:))>1.e-8);
% if max(max(abs(TAU-yt(istable,iii)')))<= 1.e-8,
% err = max(max(abs(TAU-yt(istable,iii)')));
% disp(['Model check OK ',num2str(err)]),
clear TAU A
% end
else
clear A
end
% [yt1, j01]=teff(T1,Nsam,istable);
% [yt2, j02]=teff(T2,Nsam,istable);
% [ytr, j0r]=teff(R,Nsam,istable);
%
% yt=[yt1 yt2 ytr];
save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'ac','cc','vdec','yt','mss')
else
if opt_gsa.morris==2
% [pdraws, TAU, GAM] = dynare_identification([1:npT]); %,[lpmatx lpmat(istable,:)]);
% [pdraws, TAU, GAM] = dynare_identification(options_.options_ident);
pdraws = dynare_identification(options_.options_ident);
end
load([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'ac','cc','vdec','yt','mss')
end
% for j=1:nr,
% for i=1:nc,
% y0=squeeze(A(j,i,:));
% if max(y0)-min(y0)>1.e-10,
% j0=j0+1;
% y1=ones(size(lpmat,1),1)*NaN;
% y1(istable,1)=y0;
% yt(:,j0)=y1;
% end
% end
% end
% yt = yt(:,j0);
if opt_gsa.morris==1
%OutputDir = CheckPath('gsa/screen');
if ~isempty(vdec)
if opt_gsa.load_ident_files==0
SAMorris = [];
for i=1:size(vdec,2)
[SAmeas, SAMorris(:,:,i)] = 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'],'SAvdec','vdec','ir_vdec','ic_vdec')
end
hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition');
% boxplot(SAvdec,'whis',10,'symbol','r.')
myboxplot(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
text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
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
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['EET variance decomposition observed variables']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_vdec==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAvdec(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAvdec(iv,:),[],'.',[],3)
% else
% plot(SAvdec(iv,:),'r.');
% end
% 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)])
% for ip=1:npT,
% text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_vdec_varobs_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_vdec_varobs_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_vdec_varobs_',int2str(ifig)]);
% close(gcf)
% end
% end
%
% ifig = 0;
% for j=1:M_.exo_nbr,
% if mod(j,6)==1
% figure('name',['EET variance decomposition shocks']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ic_vdec==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAvdec(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAvdec(iv,:),[],'.',[],3)
% else
% plot(SAvdec(iv,:),'r.');
% end
% 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)])
% for ip=1:npT,
% text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(M_.exo_names(j,:),'interpreter','none')
% if mod(j,6)==0 | j==M_.exo_nbr,
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_vdec_exo_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_vdec_exo_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_vdec_exo_',int2str(ifig)]);
% close(gcf),
% end
% end
if opt_gsa.load_ident_files==0
SAMorris = [];
ccac = [mss cc ac];
for i=1:size(ccac,2)
[SAmeas, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv);
end
SAcc = squeeze(SAMorris(:,1,:))';
SAcc = SAcc./(max(SAcc')'*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');
% boxplot(SAcc,'whis',10,'symbol','r.')
myboxplot(SAcc,[],'.',[],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 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:npT
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
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)
% close(gcf),
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MORRIS FOR DERIVATIVES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if opt_gsa.load_ident_files==0,
% for j=1:npT,
% SAMorris = [];
% ddd=NaN(size(lpmat,1),size(JJ,1));
% ddd(istable,:) = squeeze(JJ(:,j,:))';
% for i=1:size(ddd,2),
% [SAmeas, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], [ddd(:,i)],nliv);
% end
% SAddd(:,:,j) = squeeze(SAMorris(:,1,:))';
% SAddd(:,:,j) = SAddd(:,:,j)./(max(SAddd(:,:,j)')'*ones(1,npT));
% sad(:,j) = median(SAddd(find(~isnan(squeeze(SAddd(:,1,j)))),:,j))';
% end
% save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAddd','sad','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAddd','sad')
% end
% figure,
% contourf(sad,10), colorbar
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
% set(gca,'yticklabel',' ','fontsize',10,'ytick',[1:npT])
% for ip=1:npT,
% text(ip,0.9,['D(',bayestopt_.name{ip},')'],'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(0.9,ip,[bayestopt_.name{ip}],'rotation',0,'HorizontalAlignment','right','interpreter','none')
% end
% [m,im]=max(sad);
% iii = find((im-[1:npT])==0);
% disp('Most identified params')
% disp(bayestopt_.name(iii))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END OF MORRIS FOR DERIVATIVES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['EET cross-correlations']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_cc==j);
% iv = [iv; find(ic_cc==j)];
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAcc(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAcc(iv,:),[],'.',[],3)
% else
% plot(SAcc(iv,:),'r.');
% end
% 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)])
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_cc_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_cc_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_cc_',int2str(ifig)]);
% close(gcf),
% end
% end
% if opt_gsa.load_ident_files==0,
% SAMorris = [];
% for i=1:size(ac,2),
% [SAmeas, SAMorris(:,:,i)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], ac(:,i),nliv);
% end
% %end
% SAac = squeeze(SAMorris(:,1,:))';
% save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAac','ac','ir_ac','ic_ac','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAac','ac','ir_ac','ic_ac')
% end
% figure,
% % boxplot(SAac,'whis',10,'symbol','r.')
% myboxplot(SAac,[],'.',[],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,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title('EET All auto-correlations')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_ac'])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_ac']);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_ac']);
% close(gcf),
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['EET auto-correlations']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_ac==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAac(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAac(iv,:),[],'.',[],3)
% else
% plot(SAac(iv,:),'r.');
% end
% 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)])
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_ac_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_ac_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_ac_',int2str(ifig)]);
% close(gcf),
% end
% end
% if opt_gsa.load_ident_files==0,
% js=0;
% %for j=1:size(tadj,1),
% SAMorris = [];
% for i=1:size(tadj,2),
% js=js+1;
% [SAmeas, SAMorris(:,:,js)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], tadj(:,i),nliv);
% end
% %end
% SAM = squeeze(SAMorris(nshock+1:end,1,:));
% for j=1:js,
% SAtadj(:,j)=SAM(:,j)./(max(SAM(:,j))+eps);
% end
% SAtadj = SAtadj';
% save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAtadj','tadj','ir_tadj','ic_tadj','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAtadj','tadj','ir_tadj','ic_tadj')
% end
% if opt_gsa.load_ident_files==0,
% js=0;
% SAMorris = [];
% for i=1:size(iff,2),
% js=js+1;
% [SAmeas, SAMorriss(:,:,js)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], iff(:,i),nliv);
% end
% SAM = squeeze(SAMorriss(nshock+1:end,1,:));
% for j=1:js,
% SAIF(:,j)=SAM(:,j)./(max(SAM(:,j))+eps);
% end
% SAIF = SAIF';
% save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAIF','iff','ir_if','ic_if','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAIF','iff','ir_if','ic_if')
% end
% figure,
% %bar(SAtadj),
% % boxplot(SAtadj,'whis',10,'symbol','r.')
% myboxplot(SAtadj,[],'.',[],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')
% end
% xlabel(' ')
% title('All half-life')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_tadj'])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_tadj']);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_tadj']);
% close(gcf),
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['EET speed of adjustment observed variables']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_tadj==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAtadj(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAtadj(iv,:),[],'.',[],3)
% else
% plot(SAtadj(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_tadj_varobs_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_tadj_varobs_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_tadj_varobs_',int2str(ifig)]);
% close(gcf),
% end
% end
% ifig = 0;
% for j=1:M_.exo_nbr,
% if mod(j,6)==1
% figure('name',['EET speed of adjustment shocks']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ic_tadj==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAtadj(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAtadj(iv,:),[],'.',[],3)
% else
% plot(SAtadj(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(M_.exo_names(j,:),'interpreter','none')
% if mod(j,6)==0 | j==M_.exo_nbr,
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_tadj_exo_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_tadj_exo_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_tadj_exo_',int2str(ifig)]);
% close(gcf),
% end
% end
% figure,
% %bar(SAIF),
% % boxplot(SAIF,'whis',10,'symbol','r.')
% myboxplot(SAIF,[],'.',[],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')
% end
% xlabel(' ')
% ylabel('Elementary Effects')
% title('Steady state gains (impact factors)')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_gain'])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_gain']);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_gain']);
% close(gcf),
%figure, bar(SAIF'), title('All Gain Relationships')
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['EET steady state gain observed series']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_if==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAIF(iv,:),'whis',10,'symbol','r.');
% myboxplot(SAIF(iv,:),[],'.',[],10)
% else
% plot(SAIF(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_gain_varobs_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_gain_varobs_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_gain_varobs_',int2str(ifig)]);
% close(gcf),
% end
% end
%
% ifig = 0;
% for j=1:M_.exo_nbr,
% if mod(j,6)==1
% figure('name',['EET steady state gain shocks']);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ic_if==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAIF(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAIF(iv,:),[],'.',[],3)
% else
% plot(SAIF(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(M_.exo_names(j,:),'interpreter','none')
% if mod(j,6)==0 | j==M_.exo_nbr,
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_gain_exo_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_gain_exo_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_gain_exo_',int2str(ifig)]);
% close(gcf),
% end
% end
if opt_gsa.load_ident_files==0
SAMorris = [];
for j=1:j0
[SAmeas, SAMorris(:,:,j)] = Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv);
end
% SAM = squeeze(SAMorris(nshock+1:end,1,:));
SAM = squeeze(SAMorris(1:end,1,:));
for j=1:j0
SAnorm(:,j)=SAM(:,j)./max(SAM(:,j));
irex(j)=length(find(SAnorm(:,j)>0.01));
end
[dum, irel]=sort(irex);
% SAMmu = squeeze(SAMorris(nshock+1:end,2,:));
SAMmu = squeeze(SAMorris(1:end,2,:));
for j=1:j0
SAmunorm(:,j)=SAMmu(:,j)./max(SAM(:,j)); % normalised w.r.t. mu*
end
% SAMsig = squeeze(SAMorris(nshock+1:end,3,:));
SAMsig = squeeze(SAMorris(1:end,3,:));
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','SAmunorm','SAsignorm')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model'); %bar(SAnorm(:,irel))
% boxplot(SAnorm','whis',10,'symbol','r.')
myboxplot(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
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
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)
% hh_fig=dyn_figure(options_.nodisplay); %bar(SAmunorm(:,irel))
% % boxplot(SAmunorm','whis',10,'symbol','r.')
% myboxplot(SAmunorm',[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
% set(gca,'xlim',[0.5 npT+0.5])
% set(gca,'ylim',[-1 1])
% set(gca,'position',[0.13 0.2 0.775 0.7])
% xlabel(' ')
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title('\mu in the model')
% dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_.nodisplay,options_.graph_format);
%
% hh_fig=dyn_figure(options_.nodisplay); %bar(SAsignorm(:,irel))
% % boxplot(SAsignorm','whis',10,'symbol','r.')
% myboxplot(SAsignorm',[],'.',[],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,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title('\sigma in the model')
% dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_.nodisplay,options_.graph_format);
% figure, bar(SAnorm(:,irel)')
% set(gca,'xtick',[1:j0])
% set(gca,'xlim',[0.5 j0+0.5])
% title('Elementary effects relationships')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_morris_redform'])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_morris_redform']);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_morris_redform']);
elseif opt_gsa.morris==3
return
np=estim_params_.np;
na=(4*np+1)*opt_gsa.Nsam;
for j=1:j0
[idex(j,:), yd(j,:)] = spop_ide(lpmat, yt(:,j), opt_gsa.Nsam, 5-1);
end
iok=find(~isnan(yt(1:opt_gsa.Nsam,1)));
yr=NaN*ones(size(lpmat,1),j0);
for j=1:j0,
ys(j,:)=yd(j,:)./max(yd(j,:));
[dum, is]=sort(yt(iok,j));
yr(iok(is),j)=[1:length(iok)]'./length(iok);
yr(istable(length(iok)+1:end),j) = interp1(yt(iok,j),yr(iok,j),yt(istable(length(iok)+1:end),j),'','extrap');
ineg=find(yr(:,j)<0);
if any(ineg)
[dum, is]=sort(yr(ineg,j));
yr(ineg(is),j)=-[length(ineg):-1:1]./length(iok);
end
[idex_r(j,:), yd_r(j,:)] = spop_ide(lpmat, yr(:,j), opt_gsa.Nsam, 5-1);
ys_r(j,:)=yd_r(j,:)./max(yd_r(j,:));
end,
figure, bar((idex.*ys)./opt_gsa.Nsam), title('Relationships')
figure, bar((idex.*ys)'./opt_gsa.Nsam), title('Parameters')
figure, bar((idex_r.*ys_r)./opt_gsa.Nsam), title('Relationships rank')
figure, bar((idex_r.*ys_r)'./opt_gsa.Nsam), title('Parameters rank')
[v0,d0]=eig(corrcoef(yt(iok,:)));
ee=diag(d0);
ee=ee([end:-1:1])./j0;
i0=length(find(ee>0.01));
v0=v0(:,[end:-1:1]);
for j=1:i0
[idex_pc(j,:), yd_pc(j,:)] = spop_ide(lpmat, yt*v0(:,j), opt_gsa.Nsam, 5-1);
end
for j=1:i0
ys_pc(j,:)=yd_pc(j,:)./max(yd_pc(j,:));
end,
figure, bar((idex_pc.*ys_pc)./opt_gsa.Nsam), title('Relationships PCA')
figure, bar((idex_pc.*ys_pc)'./opt_gsa.Nsam), title('Parameters PCA')
[vr,dr]=eig(corrcoef(yr(iok,:)));
er=diag(dr);
er=er([end:-1:1])./j0;
ir0=length(find(er>0.01));
vr=vr(:,[end:-1:1]);
for j=1:ir0
[idex_pcr(j,:), yd_pcr(j,:)] = spop_ide(lpmat, yr*vr(:,j), opt_gsa.Nsam, 5-1);
end
for j=1:ir0
ys_pcr(j,:)=yd_pcr(j,:)./max(yd_pcr(j,:));
end
figure, bar((idex_pcr.*ys_pcr)./opt_gsa.Nsam), title('Relationships rank PCA')
figure, bar((idex_pcr.*ys_pcr)'./opt_gsa.Nsam), title('Parameters rank PCA')
elseif opt_gsa.morris==2 % ISKREV staff
return
else % main effects analysis
if itrans==0
fsuffix = '';
elseif itrans==1
fsuffix = '_log';
else
fsuffix = '_rank';
end
imap=[1:npT];
if isempty(lpmat0)
x0=lpmat(istable,:);
else
x0=[lpmat0(istable,:), lpmat(istable,:)];
end
nrun=length(istable);
nest=min(250,nrun);
nfit=min(1000,nrun);
% opt_gsa.load_ident_files=0;
% if opt_gsa.load_ident_files==0,
% try
% EET=load([OutputDirectoryName,'/SCREEN/',fname_,'_morris_IDE'],'SAvdec','vdec','ir_vdec','ic_vdec');
% catch
% EET=[];
% end
% SAvdec=zeros(size(vdec,2),npT);
%
% for j=1:size(vdec,2),
% if itrans==0,
% y0 = vdec(istable,j);
% elseif itrans==1,
% y0 = log_trans_(vdec(istable,j));
% else
% y0 = trank(vdec(istable,j));
% end
% if ~isempty(EET),
% % imap=find(EET.SAvdec(j,:));
% % [dum, isort]=sort(-EET.SAvdec(j,:));
% imap=find(EET.SAvdec(j,:) >= (0.1.*max(EET.SAvdec(j,:))) );
% end
% gsa_(j) = gsa_sdp(y0(1:nest), x0(1:nest,imap), ...
% 2, [],[],[],0,[OutputDirectoryName,'/map_vdec',fsuffix,int2str(j)], pnames);
% if nfit>nest,
% gsa_(j) = gsa_sdp(y0(1:nfit), x0(1:nfit,imap), ...
% -2, gsa_(j).nvr*nest^3/nfit^3,[],[],0,[OutputDirectoryName,'/map_vdec',fsuffix,int2str(j)], pnames);
% end
%
% SAvdec(j,imap)=gsa_(j).si;
% imap_vdec{j}=imap;
% end
% save([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_vdec','SAvdec','vdec','ir_vdec','ic_vdec','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_vdec','SAvdec','vdec','ir_vdec','ic_vdec')
% end
% figure,
% % boxplot(SAvdec,'whis',10,'symbol','r.')
% myboxplot(SAvdec,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% 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,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects variance decomposition ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_vdec',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_vdec',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_vdec',fsuffix]);
% close(gcf),
%
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['Main effects observed variance decomposition ',fsuffix]);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_vdec==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAvdec(iv,:),'whis',10,'symbol','r.');
% myboxplot(SAvdec(iv,:),[],'.',[],10)
% else
% plot(SAvdec(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 npT+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_vdec',fsuffix,'_varobs_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_vdec',fsuffix,'_varobs_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_vdec',fsuffix,'_varobs_',int2str(ifig)]);
% close(gcf),
% end
% end
%
% ifig = 0;
% for j=1:M_.exo_nbr,
% if mod(j,6)==1
% figure('name',['Main effects shocks variance decomposition ',fsuffix]);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ic_vdec==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAvdec(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAvdec(iv,:),[],'.',[],10)
% else
% plot(SAvdec(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',3,'xtick',[1:np])
% set(gca,'xlim',[0.5 npT+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% set(gca,'fontsize',10)
% end
% title(M_.exo_names(j,:),'interpreter','none','fontsize',10)
% if mod(j,6)==0 | j==M_.exo_nbr
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_vdec',fsuffix,'_exo_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_vdec',fsuffix,'_exo_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_vdec',fsuffix,'_exo_',int2str(ifig)]);
% close(gcf),
% end
% end
if opt_gsa.load_ident_files==0
try
EET=load([OutputDirectoryName,'/SCREEN/',fname_,'_morris_IDE'],'SAcc','ir_cc','ic_cc');
catch
EET=[];
end
ccac = stand_([mss cc ac]);
[pcc, dd] = eig(cov(ccac(istable,:)));
[latent, isort] = sort(-diag(dd));
latent = -latent;
figure, bar(latent)
title('Eigenvalues in PCA')
pcc=pcc(:,isort);
ccac = ccac*pcc;
% npca = min(40, max(find(cumsum(latent)./length(latent)<0.99))+1);
npca = max(find(cumsum(latent)./length(latent)<0.99))+1;
siPCA = (EET.SAcc'*abs(pcc'))';
% siPCA = siPCA./(max(siPCA')'*ones(1,npT)).*(latent*ones(1,npT));
siPCA = siPCA./(max(siPCA')'*ones(1,npT));
% siPCA = sum(siPCA,1);
% siPCA = siPCA./max(siPCA);
SAcc=zeros(size(ccac,2),npT);
for j=1:npca %size(ccac,2),
if itrans==0
y0 = ccac(istable,j);
elseif itrans==1
y0 = log_trans_(ccac(istable,j));
else
y0 = trank(ccac(istable,j));
end
if ~isempty(EET)
% imap=find(EET.SAvdec(j,:));
% [dum, isort]=sort(-EET.SAvdec(j,:));
imap=find(siPCA(j,:) >= (0.1.*max(siPCA(j,:))) );
% imap=find(EET.SAcc(j,:) >= (0.1.*max(EET.SAcc(j,:))) );
end
gsa_(j) = gsa_sdp(y0(1:nest), x0(1:nest,imap), ...
2, [],[],[],0,[OutputDirectoryName,'/map_cc',fsuffix,int2str(j)], pnames);
% if nfit>nest,
% gsa_(j) = gsa_sdp(y0(1:nfit), x0(1:nfit,imap), ...
% -2, gsa_(j).nvr*nest^3/nfit^3,[],[],0,[OutputDirectoryName,'/map_cc',fsuffix,int2str(j)], pnames);
% end
SAcc(j,imap)=gsa_(j).si;
imap_cc{j}=imap;
end
save([OutputDirectoryName,'/map_cc',fsuffix,'.mat'],'gsa_')
save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'imap_cc','SAcc','ccac','-append')
else
load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_cc','SAcc','ccac')
end
% figure,
% % boxplot(SAcc,'whis',10,'symbol','r.')
% myboxplot(SAcc,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% 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,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% ylabel(' ')
% title(['Main effects moments''s PCA ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_cc',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_moments',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_moments',fsuffix]);
% close(gcf),
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['Main effects cross-covariances ',fsuffix]);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_cc==j);
% iv = [iv; find(ic_cc==j)];
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAcc(iv,:),'whis',10,'symbol','r.');
% myboxplot(SAcc(iv,:),[],'.',[],10)
% else
% plot(SAcc(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 npT+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% set(gca,'fontsize',10)
% end
% title(options_.varobs(j,:),'interpreter','none','fontsize',10)
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_cc',fsuffix,'_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_cc',fsuffix,'_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_cc',fsuffix,'_',int2str(ifig)]);
% close(gcf),
% end
% end
%
% if opt_gsa.load_ident_files==0,
% try
% EET=load([OutputDirectoryName,'/SCREEN/',fname_,'_morris_IDE'],'SAac','ir_ac','ic_ac');
% catch
% EET=[];
% end
% SAac=zeros(size(ac,2),npT);
% for j=1:size(ac,2),
% if itrans==0,
% y0 = ac(istable,j);
% elseif itrans==1,
% y0 = log_trans_(ac(istable,j));
% else
% y0 = trank(ac(istable,j));
% end
% if ~isempty(EET),
% imap=find(EET.SAac(j,:) >= (0.1.*max(EET.SAac(j,:))) );
% end
% % gsa_(j) = gsa_sdp_dyn( y0, lpmat(istable,:), ...
% % gsa_flag, [],[],[],0,[OutputDirectoryName,'/map_ac',fsuffix,int2str(j)], pnames);
% gsa_(j) = gsa_sdp(y0(1:nest), x0(1:nest,imap), ...
% 2, [],[],[],0,[OutputDirectoryName,'/map_ac',fsuffix,int2str(j)], pnames);
% if nfit>nest,
% gsa_(j) = gsa_sdp(y0(1:nfit), x0(1:nfit,imap), ...
% -2, gsa_(j).nvr*nest^3/nfit^3,[],[],0,[OutputDirectoryName,'/map_ac',fsuffix,int2str(j)], pnames);
% end
% SAac(j,imap)=gsa_(j).si;
% imap_ac{j}=imap;
%
% end
% save([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_ac','SAac','ac','ir_ac','ic_ac','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_ac','SAac','ac','ir_ac','ic_ac')
% end
%
% figure,
% % boxplot(SAac,'whis',10,'symbol','r.')
% myboxplot(SAac,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% 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:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects 1 lag auto-covariances ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_ac',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_ac',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_ac',fsuffix]);
% close(gcf),
%
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['Main effects auto-covariances ',fsuffix]);
% ifig=ifig+1;
% iplo = 0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_ac==j);
% %iv = [iv; find(ic_ac==j)];
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAac(iv,:),'whis',10,'symbol','r.');
% myboxplot(SAac(iv,:),[],'.',[],10)
% else
% plot(SAac(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 npT+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% set(gca,'fontsize',10)
% end
% title(options_.varobs(j,:),'interpreter','none','fontsize',10)
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_ac',fsuffix,'_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_ac',fsuffix,'_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_ac',fsuffix,'_',int2str(ifig)]);
% close(gcf),
% end
% end
% x0=x0(:,nshock+1:end);
imap=[1:npT];
% if opt_gsa.load_ident_files==0,
% try
% EET=load([OutputDirectoryName,'/SCREEN/',fname_,'_morris_IDE'],'SAtadj','ir_tadj','ic_tadj');
% ny=size(EET.SAtadj,1);
% catch
% EET=[];
% end
% SAtadj=zeros(size(tadj,2),np);
% for j=1:size(tadj,2),
% if itrans==0,
% y0 = tadj(istable,j);
% elseif itrans==1,
% y0 = log_trans_(tadj(istable,j));
% else
% y0 = trank(tadj(istable,j));
% end
% if ~isempty(EET),
% if size(tadj,2)~=ny,
% jj=find(EET.ir_tadj==ir_tadj(j));
% jj=jj(find(EET.ic_tadj(jj)==ic_tadj(j)));
% if ~isempty(jj),
% imap=find(EET.SAtadj(jj,:) >= (0.1.*max(EET.SAtadj(jj,:))) );
% else
% imap=[1:np];
% end
% else
% imap=find(EET.SAtadj(j,:) >= (0.1.*max(EET.SAtadj(j,:))) );
% end
% end
% % gsa_(j) = gsa_sdp_dyn( y0, lpmat(istable,:), ...
% % gsa_flag, [],[],[],0,[OutputDirectoryName,'/map_tadj',fsuffix,int2str(j)], pnames);
% gsa_(j) = gsa_sdp(y0(1:nest), x0(1:nest,imap), ...
% 2, [],[],[],0,[OutputDirectoryName,'/map_tadj',fsuffix,int2str(j)], pnames);
% if nfit>nest,
% gsa_(j) = gsa_sdp(y0(1:nfit), x0(1:nfit,imap), ...
% -2, gsa_(j).nvr*nest^3/nfit^3,[],[],0,[OutputDirectoryName,'/map_tadj',fsuffix,int2str(j)], pnames);
% end
% SAtadj(j,imap)=gsa_(j).si;
% imap_tadj{j}=imap;
%
% end
% save([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_tadj','SAtadj','tadj','ir_tadj','ic_tadj','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_tadj','SAtadj','tadj','ir_tadj','ic_tadj')
% end
%
% figure,
% % boxplot(SAtadj,'whis',10,'symbol','r.')
% myboxplot(SAtadj,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+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:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects speed of adjustment ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_tadj',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_tadj',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_tadj',fsuffix]);
% close(gcf),
%
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['Main effects observed speed adjustment ',fsuffix]);
% ifig=ifig+1;
% iplo = 0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_tadj==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAtadj(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAtadj(iv,:),[],'.',[],10)
% else
% plot(SAtadj(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_tadj',fsuffix,'_varobs_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_tadj',fsuffix,'_varobs_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_tadj',fsuffix,'_varobs_',int2str(ifig)]);
% close(gcf),
% end
% end
%
% ifig = 0;
% for j=1:M_.exo_nbr,
% if mod(j,6)==1
% figure('name',['Main effects shocks speed of adjustment ',fsuffix]);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ic_tadj==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAtadj(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAtadj(iv,:),[],'.',[],10)
% else
% plot(SAtadj(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(M_.exo_names(j,:),'interpreter','none')
% if mod(j,6)==0 | j==M_.exo_nbr,
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_tadj',fsuffix,'_exo_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_tadj',fsuffix,'_exo_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_tadj',fsuffix,'_exo_',int2str(ifig)]);
% close(gcf),
% end
% end
%
%
% if opt_gsa.load_ident_files==0,
% try
% EET=load([OutputDirectoryName,'/SCREEN/',fname_,'_morris_IDE'],'SAIF','ir_if','ic_if');
% catch
% EET=[];
% end
% SAif=zeros(size(iff,2),np);
% for j=1:size(iff,2),
% if itrans==0,
% y0 = iff(istable,j);
% elseif itrans==1,
% y0 = log_trans_(iff(istable,j));
% else
% y0 = trank(iff(istable,j));
% end
% if ~isempty(EET),
% imap=find(EET.SAIF(j,:) >= (0.1.*max(EET.SAIF(j,:))) );
% end
% % gsa_(j) = gsa_sdp_dyn( y0, lpmat(istable,:), ...
% % gsa_flag, [],[],[],0,[OutputDirectoryName,'/map_if',fsuffix,int2str(j)], pnames);
% gsa_(j) = gsa_sdp(y0(1:nest), x0(1:nest,imap), ...
% 2, [],[],[],0,[OutputDirectoryName,'/map_if',fsuffix,int2str(j)], pnames);
% if nfit>nest,
% gsa_(j) = gsa_sdp(y0(1:nfit), x0(1:nfit,imap), ...
% -2, gsa_(j).nvr*nest^3/nfit^3,[],[],0,[OutputDirectoryName,'/map_if',fsuffix,int2str(j)], pnames);
% end
% SAif(j,imap)=gsa_(j).si;
% imap_if{j}=imap;
%
% end
% save([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_if','SAif','iff','ir_if','ic_if','-append')
% else
% load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_if','SAif','iff','ir_if','ic_if')
% end
%
% figure,
% % boxplot(SAif,'whis',10,'symbol','r.')
% myboxplot(SAif,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+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:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects impact factors ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_if',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_if',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_if',fsuffix]);
% close(gcf),
%
% ifig = 0;
% for j=1:size(options_.varobs,1)
% if mod(j,6)==1
% figure('name',['Main effects observed impact factors ',fsuffix]);
% ifig=ifig+1;
% iplo = 0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ir_if==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAif(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAif(iv,:),[],'.',[],10)
% else
% plot(SAif(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(options_.varobs(j,:),'interpreter','none')
% if mod(j,6)==0 | j==size(options_.varobs,1)
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_if',fsuffix,'_varobs_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_if',fsuffix,'_varobs_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_if',fsuffix,'_varobs_',int2str(ifig)]);
% close(gcf),
% end
% end
%
% ifig = 0;
% for j=1:M_.exo_nbr,
% if mod(j,6)==1
% figure('name',['Main effects shocks impact factors ',fsuffix]);
% ifig=ifig+1;
% iplo=0;
% end
% iplo=iplo+1;
% subplot(3,2,iplo)
% iv = find(ic_if==j);
% if ~isempty(iv)
% if length(iv)>1
% % boxplot(SAif(iv,:),'whis',3,'symbol','r.');
% myboxplot(SAif(iv,:),[],'.',[],10)
% else
% plot(SAif(iv,:),'r.');
% end
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+0.5])
% ydum = get(gca,'ylim');
% set(gca,'ylim',[0 ydum(2)])
% for ip=1:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% end
% title(M_.exo_names(j,:),'interpreter','none')
% if mod(j,6)==0 | j==M_.exo_nbr
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_if',fsuffix,'_exo_',int2str(ifig)])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_if',fsuffix,'_exo_',int2str(ifig)]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_if',fsuffix,'_exo_',int2str(ifig)]);
% close(gcf),
% end
% end
% SAmom = [SAvdec' SAcc' SAac']';
% SAdyn = [SAtadj' SAif']';
% SAall = [SAmom(:,nshock+1:end)' SAdyn']';
%
% figure,
% % boxplot(SAtadj,'whis',10,'symbol','r.')
% myboxplot(SAmom,[],'.',[],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,
% % text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects theoretical moments ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_moments',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_moments',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_moments',fsuffix]);
% % close(gcf),
%
% figure,
% % boxplot(SAtadj,'whis',10,'symbol','r.')
% myboxplot(SAdyn,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+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:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects short-long term dynamics ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_dynamics',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_dynamics',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_dynamics',fsuffix]);
% % close(gcf),
%
% figure,
% % boxplot(SAtadj,'whis',10,'symbol','r.')
% myboxplot(SAall,[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
% set(gca,'xlim',[0.5 np+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:np,
% text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% % text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Main effects all ',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_map_ALL',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_map_ALL',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_map_ALL',fsuffix]);
% % close(gcf),
% for j=1:size(SAall,1),
% SAallN(j,:)=SAall(j,:)./max(SAall(j,:));
% end
% SAmean=mean(SAallN);
% for j=1:size(SAmom,1),
% SAmomN(j,:)=SAmom(j,1:nshock)./max(SAmom(j,1:nshock));
% end
% SAmomN(find(isnan(SAmomN)))=0;
% SAmeanexo=mean(SAmomN(:,1:nshock));
% figure, bar(latent'*SAcc),
hh_fig=dyn_figure(options_.nodisplay,'Name',['Identifiability indices in the ',fsuffix,' moments.']);
bar(sum(SAcc))
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'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
text(ip,-0.02*(ydum(2)),bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
xlabel(' ')
title(['Identifiability indices in the ',fsuffix,' moments.'],'interpreter','none')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_ident_ALL',fsuffix],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_ident_ALL',fsuffix],1,['Identifiability indices in the ',fsuffix,' moments.'],['ident_ALL',fsuffix]',1)
% figure, bar(SAmeanexo),
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:nshock])
% set(gca,'xlim',[0.5 nshock+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:nshock,
% % text(ip,-0.02*(ydum(2)),deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
% text(ip,-0.02*(ydum(2)),bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title(['Identifiability indices for shocks',fsuffix],'interpreter','none')
% saveas(gcf,[OutputDirectoryName,'/',fname_,'_ident_SHOCKS',fsuffix])
% eval(['print -depsc2 ' OutputDirectoryName '/' fname_ '_ident_SHOCKS',fsuffix]);
% eval(['print -dpdf ' OutputDirectoryName '/' fname_ '_ident_SHOCKS',fsuffix]);
end
return
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 sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2012 European Commission
% Copyright © 2012-2017 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/>.
% % % % endif
if nargin < 5 | isempty(maxwhisker), maxwhisker = 1.5; end
if nargin < 4 | isempty(vertical), vertical = 1; end
if nargin < 3 | isempty(symbol), symbol = ['+','o']; end
if nargin < 2 | isempty(notched), notched = 0; end
if length(symbol)==1, symbol(2)=symbol(1); end
if notched==1, notched=0.25; end
a=1-notched;
% ## figure out how many data sets we have
if iscell(data)
nc = length(data);
else
% if isvector(data), data = data(:); end
nc = size(data,2);
end
% ## compute statistics
% ## s will contain
% ## 1,5 min and max
% ## 2,3,4 1st, 2nd and 3rd quartile
% ## 6,7 lower and upper confidence intervals for median
s = zeros(7,nc);
box = zeros(1,nc);
whisker_x = ones(2,1)*[1:nc,1:nc];
whisker_y = zeros(2,2*nc);
outliers_x = [];
outliers_y = [];
outliers2_x = [];
outliers2_y = [];
for i=1:nc
% ## Get the next data set from the array or cell array
if iscell(data)
col = data{i}(:);
else
col = data(:,i);
end
% ## Skip missing data
% % % % % % % col(isnan(col) | isna (col)) = [];
col(isnan(col)) = [];
% ## Remember the data length
nd = length(col);
box(i) = nd;
if (nd > 1)
% ## min,max and quartiles
% s(1:5,i) = statistics(col)(1:5);
s(1,i)=min(col);
s(5,i)=max(col);
s(2,i)=myprctilecol(col,25);
s(3,i)=myprctilecol(col,50);
s(4,i)=myprctilecol(col,75);
% ## confidence interval for the median
est = 1.57*(s(4,i)-s(2,i))/sqrt(nd);
s(6,i) = max([s(3,i)-est, s(2,i)]);
s(7,i) = min([s(3,i)+est, s(4,i)]);
% ## whiskers out to the last point within the desired inter-quartile range
IQR = maxwhisker*(s(4,i)-s(2,i));
whisker_y(:,i) = [min(col(col >= s(2,i)-IQR)); s(2,i)];
whisker_y(:,nc+i) = [max(col(col <= s(4,i)+IQR)); s(4,i)];
% ## outliers beyond 1 and 2 inter-quartile ranges
outliers = col((col < s(2,i)-IQR & col >= s(2,i)-2*IQR) | (col > s(4,i)+IQR & col <= s(4,i)+2*IQR));
outliers2 = col(col < s(2,i)-2*IQR | col > s(4,i)+2*IQR);
outliers_x = [outliers_x; i*ones(size(outliers))];
outliers_y = [outliers_y; outliers];
outliers2_x = [outliers2_x; i*ones(size(outliers2))];
outliers2_y = [outliers2_y; outliers2];
elseif (nd == 1)
% ## all statistics collapse to the value of the point
s(:,i) = col;
% ## single point data sets are plotted as outliers.
outliers_x = [outliers_x; i];
outliers_y = [outliers_y; col];
else
% ## no statistics if no points
s(:,i) = NaN;
end
end
% % % % if isempty(outliers2_y)
% % % % outliers2_y=
% ## Note which boxes don't have enough stats
chop = find(box <= 1);
% ## Draw a box around the quartiles, with width proportional to the number of
% ## items in the box. Draw notches if desired.
box = box*0.23/max(box);
quartile_x = ones(11,1)*[1:nc] + [-a;-1;-1;1;1;a;1;1;-1;-1;-a]*box;
quartile_y = s([3,7,4,4,7,3,6,2,2,6,3],:);
% ## Draw a line through the median
median_x = ones(2,1)*[1:nc] + [-a;+a]*box;
% median_x=median(col);
median_y = s([3,3],:);
% ## Chop all boxes which don't have enough stats
quartile_x(:,chop) = [];
quartile_y(:,chop) = [];
whisker_x(:,[chop,chop+nc]) = [];
whisker_y(:,[chop,chop+nc]) = [];
median_x(:,chop) = [];
median_y(:,chop) = [];
% % % %
% ## Add caps to the remaining whiskers
cap_x = whisker_x;
cap_x(1,:) =cap_x(1,:)- 0.05;
cap_x(2,:) =cap_x(2,:)+ 0.05;
cap_y = whisker_y([1,1],:);
% #quartile_x,quartile_y
% #whisker_x,whisker_y
% #median_x,median_y
% #cap_x,cap_y
%
% ## Do the plot
mm=min(min(data));
MM=max(max(data));
if vertical
plot (quartile_x, quartile_y, 'b', ...
whisker_x, whisker_y, 'b--', ...
cap_x, cap_y, 'k', ...
median_x, median_y, 'r', ...
outliers_x, outliers_y, [symbol(1),'r'], ...
outliers2_x, outliers2_y, [symbol(2),'r']);
set(gca,'XTick',1:nc);
set(gca, 'XLim', [0.5, nc+0.5]);
set(gca, 'YLim', [mm-(MM-mm)*0.05-eps, MM+(MM-mm)*0.05+eps]);
else
% % % % % plot (quartile_y, quartile_x, "b;;",
% % % % % whisker_y, whisker_x, "b;;",
% % % % % cap_y, cap_x, "b;;",
% % % % % median_y, median_x, "r;;",
% % % % % outliers_y, outliers_x, [symbol(1),"r;;"],
% % % % % outliers2_y, outliers2_x, [symbol(2),"r;;"]);
end
if nargout
sout=s;
end
% % % endfunction
\ No newline at end of file
function pick
%
% Copyright © 2001-2017 European Commission
% Copyright © 2017-2023 DynareTeam
% This file is part of GLUEWIN
% GLUEWIN is a MATLAB code designed for analysing the output
% of Monte Carlo runs when empirical observations of the model output are available
% and implements the GSA-GLUE methodology by Ratto et al. [1], based on a combination
% of GLUE (Generalised Likelihood Uncertainty Estimation) by K. Beven [2] and GSA
% Global Sensitivity Analysis) [3].']
% The program has been developed by M. Ratto, European Commission, Joint Research Centre,
% Institute for the Protection and Security of The Citizen, Technological and Economic Risk Management,
% Applied Statistics, as a deliverable of the IMPACT project
% (EC Fifth Framework Programme, SCA Project, IST-1999-11313, DG-INFSO).
%
% The graphical layout of the code is inspired by the freeware GLUE package by K. Beven,
% vailable at the Lancaster University web site on the page [4]:
% http://www.es.lancs.ac.uk/hfdg/glue.html
% to which the GLUEWIN code introduces several extensions and additional options.
% Thanks are due to R. Girardi, A. Rossi, A. Saltelli, S. Tarantola and U. Callies for numerous
% comments and suggestions.
% For more information, please contact marco.ratto@ec.europa.eu
%
% Disclaimer: This software has been developed at the Joint Research Centre of European Commission
% by officers in the course of their official duties. 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.
%
% [1] Ratto, M., Tarantola, S., A. Saltelli, Sensitivity analysis in model calibration: GSA-GLUE approach.
% 'Computer Physics Communications, 136, 2001, 212-224
% [2] Beven K.J., Binley A., The Future of Distributed Models: Model Calibration and Uncertainty
% 'Prediction, Hydrological Processes, 6, 279-298, 1992
% [3] Saltelli, A., K. Chan, M. Scott, Editors, (2000), Sensitivity analysis, John Wiley & Sons
% 'publishers, Probability and Statistics series.
% [4] Beven K., GLUE for Windows User manual, 1998.
pmenu=findobj(gcf,'type','uicontextmenu','Tag','Run viewer');
button1=findobj(gcf,'type','uimenu','Tag','save params');
button2=findobj(gcf,'type','uimenu','Tag','eval params');
%button=get(pmenu,'children');
gg=gco;
ax0=gca;
set(gg,'buttondownfcn',[]);
c=get(gca,'currentpoint');
x=c(1,1);
y=c(1,2);
X=get(gco,'xdata');
Y=get(gco,'ydata');
dx=get(gca,'xlim');
dy=get(gca,'ylim');
pos=get(gca,'position');
scalex=dx(2)-dx(1);
scaley=dy(2)-dy(1);
if length(X)>1
K = dsearchn([(Y./scaley)' (X./scalex)'],[y/scaley x/scalex]);
else
az=get(gca,'children');
T =get(az(end),'ydata');
[dum K]=max(T);
end
KK=K;
set(button1,'Label',['Save ',num2str(K)],'Callback',['scatter_callback(',num2str(KK),',''save'')']);
set(button2,'Label',['Eval ',num2str(K)],'Callback',['scatter_callback(',num2str(KK),',''eval'')']);
hh_obj=findobj(gcf,'type','axes','Tag','scatter');
for k=1:length(hh_obj)
axes(hh_obj(k));
dum=get(gca,'children');
dumx=get(dum(end),'xdata');
dumy=get(dum(end),'ydata');
xmid=min(dumx) + 0.5*(max(dumx)-min(dumx));
hold on
plot(dumx(KK),dumy(KK),'or');
if dumx(KK) < xmid
text(dumx(KK),dumy(KK),[' ',num2str(K)], ...
'FontWeight','Bold',...
'Color','r');
else
text(dumx(KK),dumy(KK),[num2str(K),' '], ...
'HorizontalAlignment','right', ...
'FontWeight','Bold',...
'Color','r');
end
hold off
end
\ No newline at end of file
function [gend, data] = read_data()
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2012-2015 European Commission
% Copyright © 2012-2017 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/>.
global options_
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
data = transpose(rawdata-ones(gend,1)* mean(rawdata,1));
else
data = transpose(rawdata);
end
if ~isreal(rawdata)
error(['There are complex values in the data. Probably a wrong' ...
' transformation'])
end
function d=hess_element(func,element1,element2,args)
% function d=hess_element(func,element1,element2,args)
% returns an entry of the finite differences approximation to the hessian of func
% returns an entry of the finite differences approximation to the Hessian of func
%
% INPUTS
% func [function name] string with name of the function
% element1 [int] the indices showing the element within the hessian that should be returned
% element1 [int] the indices showing the element within the Hessian that should be returned
% element2 [int]
% args [cell array] arguments provided to func
%
% OUTPUTS
% d [double] the (element1,element2) entry of the hessian
% d [double] the (element1,element2) entry of the Hessian
%
% SPECIAL REQUIREMENTS
% none
......
......@@ -78,7 +78,7 @@ hessian_mat = zeros(size(f0,1), n*n);
for i=1:n
if i > 1
%fill symmetric part of Hessian based on previously computed results
k = [i:n:n*(i-1)];
k = i:n:n*(i-1);
hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1) = hessian_mat(:,k);
end
hessian_mat(:,(i-1)*n+i) = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23
......
function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
%function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
function [endo_histval, exo_histval, exo_det_histval] = histvalf(M_, options_)
%function [endo_histval, exo_histval, exo_det_histval] = histvalf(M_, options_)
% Sets initial values for simulation using values contained in `fname`, a
% file possibly created by a call to `smoother2histval`
%
......@@ -13,7 +13,7 @@ function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
% none
% Copyright © 2014-2021 Dynare Team
% Copyright © 2014-2025 Dynare Team
%
% This file is part of Dynare.
%
......@@ -30,33 +30,21 @@ function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
series = histvalf_initvalf('HISTVAL', M, options);
k = M.orig_maximum_lag - M.maximum_lag + 1;
series = histvalf_initvalf('HISTVAL', M_, options_);
if ~isoctave && matlab_ver_less_than('9.7')
% Workaround for MATLAB bug described in dseries#45
% The solution is to avoid using the "end" keyword
myend = nobs(series);
if series.nobs ~= M_.orig_maximum_lag
error('histval_file: number of observations does not correspond to the number of lags in the model')
end
endo_histval = series{M.endo_names{:}}.data(k:myend, :)';
k = M_.orig_maximum_lag - M_.maximum_lag + 1;
exo_histval = [];
if M.exo_nbr
exo_histval = series{M.exo_names{:}}.data(k:myend, :)';
end
exo_det_histval = [];
if M.exo_det_nbr
exo_det_histval = series{M.exo_names{:}}.data(k:myend, :)';
end
else
endo_histval = series{M.endo_names{:}}.data(k:end, :)';
endo_histval = series{M_.endo_names{:}}.data(k:end, :)';
exo_histval = [];
if M.exo_nbr
exo_histval = series{M.exo_names{:}}.data(k:end, :)';
if M_.exo_nbr
exo_histval = series{M_.exo_names{:}}.data(k:end, :)';
end
exo_det_histval = [];
if M.exo_det_nbr
exo_det_histval = series{M.exo_names{:}}.data(k:end, :)';
end
if M_.exo_det_nbr
exo_det_histval = series{M_.exo_names{:}}.data(k:end, :)';
end