1function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0) 2%function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0) 3% ------------------------------------------------------------------------- 4% This function is called, when the user specifies identification(...); in the mod file. It prepares all identification analysis: 5% (1) set options, local and persistent variables for a new identification 6% analysis either for a single point or a MC Sample. It also displays and plots the results 7% (2) load, display and plot a previously saved identification analysis 8% 9% Note 1: This function does not output the arguments to the workspace, but saves results to the folder identification 10% Note 2: If you want to use this function directly in the mod file and workspace, you still have 11% to put identification in your mod file, otherwise the preprocessor won't provide all necessary objects 12% ========================================================================= 13% INPUTS 14% * options_ident [structure] identification options 15% * pdraws0 [SampleSize by totparam_nbr] optional: matrix of MC sample of model parameters 16% ------------------------------------------------------------------------- 17% OUTPUTS 18% * pdraws [matrix] MC sample of model params used 19% * STO_REDUCEDFORM, [matrix] MC sample of entries in steady state and reduced form model solution (stacked vertically) 20% * STO_MOMENTS, [matrix] MC sample of entries in theoretical first two moments (stacked vertically) 21% * STO_DYNAMIC, [matrix] MC sample of entries in steady state and dynamic model derivatives (stacked vertically) 22% * STO_si_dDYNAMIC, [matrix] MC sample of derivatives of steady state and dynamic derivatives 23% * STO_si_dREDUCEDFORM, [matrix] MC sample of derivatives of steady state and reduced form model solution 24% * STO_si_dMOMENTS [matrix] MC sample of Iskrev (2010)'s J matrix 25% * STO_dSPECTRUM [matrix] MC sample of Qu and Tkachenko (2012)'s \bar{G} matrix 26% * STO_dMINIMAL [matrix] MC sample of Komunjer and Ng (2011)'s \bar{\Delta} matrix 27% ------------------------------------------------------------------------- 28% This function is called by 29% * driver.m 30% * map_ident_.m 31% ------------------------------------------------------------------------- 32% This function calls 33% * checkpath 34% * disp_identification 35% * dyn_waitbar 36% * dyn_waitbar_close 37% * get_all_parameters 38% * get_posterior_parameters 39% * get_the_name 40% * identification_analysis 41% * isoctave 42% * plot_identification 43% * prior_draw 44% * set_default_option 45% * set_prior 46% * skipline 47% * vnorm 48% ========================================================================= 49% Copyright (C) 2010-2020 Dynare Team 50% 51% This file is part of Dynare. 52% 53% Dynare is free software: you can redistribute it and/or modify 54% it under the terms of the GNU General Public License as published by 55% the Free Software Foundation, either version 3 of the License, or 56% (at your option) any later version. 57% 58% Dynare is distributed in the hope that it will be useful, 59% but WITHOUT ANY WARRANTY; without even the implied warranty of 60% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 61% GNU General Public License for more details. 62% 63% You should have received a copy of the GNU General Public License 64% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 65% ========================================================================= 66 67global M_ options_ oo_ bayestopt_ estim_params_ 68 69store_options_ = options_; % store options to restore them at the end 70fname = M_.fname; %model name 71dname = M_.dname; %model name 72 73%turn warnings off, either globally or only relevant ids 74if isoctave 75 %warning('off'), 76 warning('off','Octave:singular-matrix'); 77 warning('off','Octave:nearly-singular-matrix'); 78 warning('off','Octave:neg-dim-as-zero'); 79 warning('off','Octave:array-as-logical'); 80else 81 %warning off; 82 warning('off','MATLAB:rankDeficientMatrix'); 83 warning('off','MATLAB:singularMatrix'); 84 warning('off','MATLAB:nearlySingularMatrix'); 85 warning('off','MATLAB:plot:IgnoreImaginaryXYPart'); 86 warning('off','MATLAB:specgraph:private:specgraph:UsingOnlyRealComponentOfComplexData'); 87 warning('off','MATLAB:handle_graphics:exceptions:SceneNode'); 88 warning('off','MATLAB:divideByZero'); 89 warning('off','MATLAB:log:logOfZero'); 90end 91 92%% Set all options in options_ident and create objects 93options_ident = set_default_option(options_ident,'gsa_sample_file',0); 94 % 0: do not use sample file 95 % 1: triggers gsa prior sample 96 % 2: triggers gsa Monte-Carlo sample (i.e. loads a sample corresponding to pprior=0 and ppost=0 in dynare_sensitivity options) 97 % FILENAME: use sample file in provided path 98options_ident = set_default_option(options_ident,'parameter_set','prior_mean'); 99 % 'calibration': use values in M_.params and M_.Sigma_e to update estimated stderr, corr and model parameters (get_all_parameters) 100 % 'prior_mode': use values in bayestopt_.p5 prior to update estimated stderr, corr and model parameters 101 % 'prior_mean': use values in bayestopt_.p1 prior to update estimated stderr, corr and model parameters 102 % 'posterior_mode': use posterior mode values in estim_params_ to update estimated stderr, corr and model parameters (get_posterior_parameters('mode',...)) 103 % 'posterior_mean': use posterior mean values in estim_params_ to update estimated stderr, corr and model parameters (get_posterior_parameters('mean',...)) 104 % 'posterior_median': use posterior median values in estim_params_ to update estimated stderr, corr and model parameters (get_posterior_parameters('median',...)) 105options_ident = set_default_option(options_ident,'load_ident_files',0); 106 % 1: load previously computed analysis from identification/fname_identif.mat 107options_ident = set_default_option(options_ident,'useautocorr',0); 108 % 1: use autocorrelations in Iskrev (2010)'s theoretical second moments criteria 109 % 0: use autocovariances in Iskrev (2010)'s theoretical second moments criteria 110options_ident = set_default_option(options_ident,'ar',1); 111 % number of lags to consider for autocovariances/autocorrelations in Iskrev (2010)'s criteria 112options_ident = set_default_option(options_ident,'prior_mc',1); 113 % size of Monte-Carlo sample of parameter draws 114options_ident = set_default_option(options_ident,'prior_range',0); 115 % 1: sample uniformly from prior ranges implied by the prior specifications (overwrites prior shape when prior_mc > 1) 116 % 0: sample from specified prior distributions (when prior_mc > 1) 117options_ident = set_default_option(options_ident,'periods',300); 118 % length of stochastic simulation to compute simulated moment uncertainty (when analytic Hessian is not available) 119options_ident = set_default_option(options_ident,'replic',100); 120 % number of replicas to compute simulated moment uncertainty (when analytic Hessian is not available) 121options_ident = set_default_option(options_ident,'advanced',0); 122 % 1: show a more detailed analysis based on reduced-form model solution and dynamic model derivatives. Further, performs a brute force 123 % search of the groups of parameters best reproducing the behavior of each single parameter. 124options_ident = set_default_option(options_ident,'normalize_jacobians',1); 125 % 1: normalize Jacobians by either rescaling each row by its largest element in absolute value or for Gram (or Hessian-type) matrices by transforming into correlation-type matrices 126options_ident = set_default_option(options_ident,'grid_nbr',5000); 127 % number of grid points in [-pi;pi] to approximate the integral to compute Qu and Tkachenko (2012)'s criteria 128 % note that grid_nbr needs to be even and actually we use (grid_nbr+1) points, as we add the 0 frequency and use symmetry 129 if mod(options_ident.grid_nbr,2) ~= 0 130 options_ident.grid_nbr = options_ident.grid_nbr+1; 131 fprintf('IDENTIFICATION: ''grid_nbr'' needs to be even, hence it is reset to %d\n',options_ident.grid_nbr) 132 if mod(options_ident.grid_nbr,2) ~= 0 133 error('IDENTIFICATION: You need to set an even value for ''grid_nbr'''); 134 end 135 end 136options_ident = set_default_option(options_ident,'tol_rank','robust'); 137 % tolerance level used for rank computations 138options_ident = set_default_option(options_ident,'tol_deriv',1.e-8); 139 % tolerance level for selecting columns of non-zero derivatives 140options_ident = set_default_option(options_ident,'tol_sv',1.e-3); 141 % tolerance level for selecting non-zero singular values in identification_checks.m 142 143%check whether to compute identification strength based on information matrix 144if ~isfield(options_ident,'no_identification_strength') 145 options_ident.no_identification_strength = 0; 146else 147 options_ident.no_identification_strength = 1; 148end 149%check whether to compute and display identification criteria based on steady state and reduced-form model solution 150if ~isfield(options_ident,'no_identification_reducedform') 151 options_ident.no_identification_reducedform = 0; 152else 153 options_ident.no_identification_reducedform = 1; 154end 155%check whether to compute and display identification criteria based on Iskrev (2010), i.e. derivative of first two moments 156if ~isfield(options_ident,'no_identification_moments') 157 options_ident.no_identification_moments = 0; 158else 159 options_ident.no_identification_moments = 1; 160end 161%check whether to compute and display identification criteria based on Komunjer and Ng (2011), i.e. derivative of first moment, minimal state space system and observational equivalent spectral density transformation 162if ~isfield(options_ident,'no_identification_minimal') 163 options_ident.no_identification_minimal = 0; 164else 165 options_ident.no_identification_minimal = 1; 166end 167%Check whether to compute and display identification criteria based on Qu and Tkachenko (2012), i.e. Gram matrix of derivatives of first moment plus outer product of derivatives of spectral density 168if ~isfield(options_ident,'no_identification_spectrum') 169 options_ident.no_identification_spectrum = 0; 170else 171 options_ident.no_identification_spectrum = 1; 172end 173 174%overwrite setting, as identification strength and advanced need both reducedform and moments criteria 175if (isfield(options_ident,'no_identification_strength') && options_ident.no_identification_strength == 0) || (options_ident.advanced == 1) 176 options_ident.no_identification_reducedform = 0; 177 options_ident.no_identification_moments = 0; 178end 179 180%overwrite setting, as dynare_sensitivity does not make use of spectrum and minimal system 181if isfield(options_,'opt_gsa') && isfield(options_.opt_gsa,'identification') && options_.opt_gsa.identification == 1 182 options_ident.no_identification_minimal = 1; 183 options_ident.no_identification_spectrum = 1; 184end 185 186%Deal with non-stationary cases 187if isfield(options_ident,'diffuse_filter') 188 options_ident.lik_init=3; %overwrites any other lik_init set 189 options_.diffuse_filter=options_ident.diffuse_filter; %make options_ inherit diffuse filter; will overwrite any conflicting lik_init in dynare_estimation_init 190else 191 if options_.diffuse_filter==1 %warning if estimation with diffuse filter was done, but not passed 192 fprintf('WARNING IDENTIFICATION: Previously the diffuse_filter option was used, but it was not passed to the identification command. This may result in problems if your model contains unit roots.\n'); 193 end 194end 195options_ident = set_default_option(options_ident,'lik_init',1); 196options_.lik_init=options_ident.lik_init; %make options_ inherit lik_init 197if options_ident.lik_init==3 %user specified diffuse filter using the lik_init option 198 options_ident.analytic_derivation=0; %diffuse filter not compatible with analytic derivation 199end 200 % Type of initialization of Kalman filter: 201 % 1: stationary models: initial matrix of variance of error of forecast is set equal to the unconditional variance of the state variables 202 % 2: nonstationary models: wide prior is used with an initial matrix of variance of the error of forecast diagonal with 10 on the diagonal (follows the suggestion of Harvey and Phillips(1979)) 203 % 3: nonstationary models: use a diffuse filter (use rather the diffuse_filter option) 204 % 4: filter is initialized with the fixed point of the Riccati equation 205 % 5: i) option 2 for non-stationary elements by setting their initial variance in the forecast error matrix to 10 on the diagonal and all co-variances to 0 and 206 % ii) option 1 for the stationary elements 207options_ident = set_default_option(options_ident,'analytic_derivation',1); 208 % 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3 209options_ident = set_default_option(options_ident,'order',1); 210 % 1: first-order perturbation approximation, identification is based on linear state space system 211 % 2: second-order perturbation approximation, identification is based on second-order pruned state space system 212 % 3: third-order perturbation approximation, identification is based on third-order pruned state space system 213 214%overwrite values in options_, note that options_ is restored at the end of the function 215if isfield(options_ident,'TeX') 216 options_.TeX=options_ident.TeX; 217 % requests printing of results and graphs in LaTeX 218end 219if isfield(options_ident,'nograph') 220 options_.nograph=options_ident.nograph; 221 % do not display and do not save graphs 222end 223if isfield(options_ident,'nodisplay') 224 options_.nodisplay=options_ident.nodisplay; 225 % do not display, but save graphs 226end 227if isfield(options_ident,'graph_format') 228 options_.graph_format=options_ident.graph_format; 229 % specify file formats to save graphs: eps, pdf, fig, none (do not save but only display) 230end 231 232% check for external draws, i.e. set pdraws0 for a gsa analysis 233if options_ident.gsa_sample_file 234 GSAFolder = checkpath('gsa',dname); 235 if options_ident.gsa_sample_file==1 236 load([GSAFolder,filesep,fname,'_prior'],'lpmat','lpmat0','istable'); 237 elseif options_ident.gsa_sample_file==2 238 load([GSAFolder,filesep,fname,'_mc'],'lpmat','lpmat0','istable'); 239 else 240 load([GSAFolder,filesep,options_ident.gsa_sample_file],'lpmat','lpmat0','istable'); 241 end 242 if isempty(istable) 243 istable=1:size(lpmat,1); 244 end 245 if ~isempty(lpmat0) 246 lpmatx=lpmat0(istable,:); 247 else 248 lpmatx=[]; 249 end 250 pdraws0 = [lpmatx lpmat(istable,:)]; 251 clear lpmat lpmat0 istable; 252elseif nargin==1 253 pdraws0=[]; 254end 255external_sample=0; 256if nargin==2 || ~isempty(pdraws0) 257 % change settings if there is an external sample provided as input argument 258 options_ident.prior_mc = size(pdraws0,1); 259 options_ident.load_ident_files = 0; 260 external_sample = 1; 261end 262 263% Check if estimated_params block is provided 264if isempty(estim_params_) 265 prior_exist = 0; 266 %reset some options 267 options_ident.prior_mc = 1; 268 options_ident.prior_range = 0; 269 options_.identification_check_endogenous_params_with_no_prior = true; %needed to trigger endogenous steady state parameter check in dynare_estimation_init 270else 271 prior_exist = 1; 272 parameters = options_ident.parameter_set; 273end 274 275% overwrite settings in options_ and prepare to call dynare_estimation_init 276options_.order = options_ident.order; 277if options_ident.order > 1 278 %order>1 is not compatible with analytic derivation in dsge_likelihood.m 279 options_ident.analytic_derivation=0; 280 %order>1 is based on pruned state space system 281 options_.pruning = true; 282end 283if options_ident.order == 3 284 options_.k_order_solver = 1; 285end 286options_.ar = options_ident.ar; 287options_.prior_mc = options_ident.prior_mc; 288options_.Schur_vec_tol = 1.e-8; 289options_.nomoments = 0; 290options_.analytic_derivation=options_ident.analytic_derivation; 291 % 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3 292options_ = set_default_option(options_,'datafile',''); 293options_.mode_compute = 0; 294options_.plot_priors = 0; 295options_.smoother = 1; 296options_.options_ident = []; 297[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = dynare_estimation_init(M_.endo_names, fname, 1, M_, options_, oo_, estim_params_, bayestopt_); 298 299% set method to compute identification Jacobians (kronflag). Default:0 300options_ident = set_default_option(options_ident,'analytic_derivation_mode', options_.analytic_derivation_mode); % if not set by user, inherit default global one 301 % 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012) 302 % 1: kronecker products method to compute analytical derivatives as in Iskrev (2010) (only for order=1) 303 % -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification_numerical_objective.m (previously thet2tau.m) 304 % -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the identification Jacobians are then computed analytically as with 0 305options_.analytic_derivation_mode = options_ident.analytic_derivation_mode; %overwrite setting in options_ 306 307% initialize persistent variables in prior_draw 308if prior_exist 309 if any(bayestopt_.pshape > 0) 310 if options_ident.prior_range 311 %sample uniformly from prior ranges (overwrite prior specification) 312 prior_draw(bayestopt_, options_.prior_trunc, true); 313 else 314 %sample from prior distributions 315 prior_draw(bayestopt_, options_.prior_trunc, false); 316 end 317 else 318 options_ident.prior_mc = 1; %only one single point 319 end 320end 321 322% check if identification directory is already created 323IdentifDirectoryName = CheckPath('identification',dname); 324 325% Create indices for stderr, corr and model parameters 326if prior_exist % use estimated_params block 327 indpmodel = []; %initialize index for model parameters 328 if ~isempty(estim_params_.param_vals) 329 indpmodel = estim_params_.param_vals(:,1); %values correspond to parameters declaration order, row number corresponds to order in estimated_params 330 end 331 indpstderr=[]; %initialize index for stderr parameters 332 if ~isempty(estim_params_.var_exo) 333 indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params 334 end 335 indpcorr=[]; %initialize matrix for corr paramters 336 if ~isempty(estim_params_.corrx) 337 indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params 338 end 339 totparam_nbr = length(bayestopt_.name); % number of estimated stderr, corr and model parameters as declared in estimated_params 340 modparam_nbr = estim_params_.np; % number of model parameters as declared in estimated_params 341 stderrparam_nbr = estim_params_.nvx; % number of stderr parameters 342 corrparam_nbr = estim_params_.ncx; % number of corr parameters 343 if estim_params_.nvn || estim_params_.ncn %nvn is number of stderr parameters and ncn is number of corr parameters of measurement innovations as declared in estimated_params 344 error('Identification does not (yet) support measurement errors. Instead, define them explicitly as varexo and provide measurement equations in the model definition.') 345 end 346 name = cell(totparam_nbr,1); %initialize cell for parameter names 347 name_tex = cell(totparam_nbr,1); %initialize cell for TeX parameter names 348 for jj=1:totparam_nbr 349 if options_.TeX 350 [param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_); 351 name_tex{jj,1} = strrep(param_name_tex_temp,'$',''); %ordering corresponds to estimated_params 352 name{jj,1} = param_name_temp; %ordering corresponds to estimated_params 353 else 354 param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_); 355 name{jj,1} = param_name_temp; %ordering corresponds to estimated_params 356 end 357 end 358else % no estimated_params block, choose all model parameters and all stderr parameters, but no corr parameters 359 indpmodel = 1:M_.param_nbr; %all model parameters 360 indpstderr = 1:M_.exo_nbr; %all stderr parameters 361 indpcorr = []; %no corr parameters 362 stderrparam_nbr = M_.exo_nbr; 363 corrparam_nbr = 0; 364 modparam_nbr = M_.param_nbr; 365 totparam_nbr = modparam_nbr+stderrparam_nbr; 366 name = cellfun(@(x) horzcat('SE_', x), M_.exo_names, 'UniformOutput', false); %names for stderr parameters 367 name = vertcat(name, M_.param_names); 368 name_tex = cellfun(@(x) horzcat('$ SE_{', x, '} $'), M_.exo_names, 'UniformOutput', false); 369 name_tex = vertcat(name_tex, M_.param_names_tex); 370 if ~isequal(M_.H,0) 371 fprintf('\ndynare_identification:: Identification does not support measurement errors (yet) and will ignore them in the following. To test their identifiability, instead define them explicitly as varexo and provide measurement equations in the model definition.\n') 372 end 373end 374options_ident.name_tex = name_tex; 375 376fprintf('\n======== Identification Analysis ========\n') 377if options_ident.order > 1 378 fprintf('Based on Pruned State Space System (order=%d)\n',options_ident.order); 379end 380skipline() 381if totparam_nbr < 2 382 options_ident.advanced = 0; 383 fprintf('There is only one parameter to study for identitification. The advanced option is re-set to 0.\n') 384end 385 386% settings dependent on number of parameters 387options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,totparam_nbr-1])); 388options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,totparam_nbr-1]); 389 % In brute force search (ident_bruteforce.m) when advanced=1 this option sets the maximum dimension of groups of parameters that best reproduce the behavior of each single model parameter 390 391options_ident = set_default_option(options_ident,'checks_via_subsets',0); 392 % 1: uses identification_checks_via_subsets.m to compute problematic parameter combinations 393 % 0: uses identification_checks.m to compute problematic parameter combinations [default] 394options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4,totparam_nbr-1])); 395 % In identification_checks_via_subsets.m, when checks_via_subsets=1, this option sets the maximum dimension of groups of parameters for which the corresponding rank criteria is checked 396 397 398% store identification options 399options_.options_ident = options_ident; 400store_options_ident = options_ident; 401 402% get some options for quick reference 403iload = options_ident.load_ident_files; 404SampleSize = options_ident.prior_mc; 405 406if iload <=0 407 %% Perform new identification analysis, i.e. do not load previous analysis 408 if prior_exist 409 % use information from estimated_params block 410 params = set_prior(estim_params_,M_,options_)'; 411 if all(bayestopt_.pshape == 0) 412 % only bounds are specified in estimated_params 413 parameters = 'ML_Starting_value'; 414 parameters_TeX = 'ML starting value'; 415 fprintf('Testing ML Starting value\n'); 416 else 417 % use user-defined option 418 switch parameters 419 case 'calibration' 420 parameters_TeX = 'Calibration'; 421 fprintf('Testing calibration\n'); 422 params(1,:) = get_all_parameters(estim_params_,M_); 423 case 'posterior_mode' 424 parameters_TeX = 'Posterior mode'; 425 fprintf('Testing posterior mode\n'); 426 params(1,:) = get_posterior_parameters('mode',M_,estim_params_,oo_,options_); 427 case 'posterior_mean' 428 parameters_TeX = 'Posterior mean'; 429 fprintf('Testing posterior mean\n'); 430 params(1,:) = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); 431 case 'posterior_median' 432 parameters_TeX = 'Posterior median'; 433 fprintf('Testing posterior median\n'); 434 params(1,:) = get_posterior_parameters('median',M_,estim_params_,oo_,options_); 435 case 'prior_mode' 436 parameters_TeX = 'Prior mode'; 437 fprintf('Testing prior mode\n'); 438 params(1,:) = bayestopt_.p5(:); 439 case 'prior_mean' 440 parameters_TeX = 'Prior mean'; 441 fprintf('Testing prior mean\n'); 442 params(1,:) = bayestopt_.p1; 443 otherwise 444 fprintf('The option parameter_set has to be equal to: ''calibration'', ''posterior_mode'', ''posterior_mean'', ''posterior_median'', ''prior_mode'' or ''prior_mean''.\n'); 445 error('IDENTIFICATION: The option ''parameter_set'' has an invalid value'); 446 end 447 end 448 else 449 % no estimated_params block is available, all stderr and model parameters, but no corr parameters are chosen 450 params = [sqrt(diag(M_.Sigma_e))', M_.params']; 451 parameters = 'Current_params'; 452 parameters_TeX = 'Current parameter values'; 453 fprintf('Testing all current stderr and model parameter values\n'); 454 end 455 options_ident.tittxt = parameters; %title text for graphs and figures 456 % perform identification analysis for single point 457 [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, options_ident] = ... 458 identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables 459 if info(1)~=0 460 % there are errors in the solution algorithm 461 message = get_error_message(info,options_); 462 fprintf('-----------\n'); 463 fprintf('The model does not solve for %s (info = %d: %s)\n', parameters, info(1), message); 464 fprintf('-----------\n'); 465 if any(bayestopt_.pshape) 466 % if there are errors in the solution algorithm, try to sample a different point from the prior 467 fprintf('Try sampling up to 50 parameter sets from the prior.\n'); 468 kk=0; 469 while kk<50 && info(1) 470 kk=kk+1; 471 params = prior_draw(); 472 options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures 473 % perform identification analysis 474 [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, options_ident] = ... 475 identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); 476 end 477 end 478 if info(1) 479 fprintf('\n-----------\n'); 480 fprintf('Identification stopped:\n'); 481 if any(bayestopt_.pshape) 482 fprintf('The model did not solve for any of 50 attempts of random samples from the prior\n'); 483 end 484 fprintf('-----------\n'); 485 return 486 else 487 % found a (random) point that solves the model 488 fprintf('Found a random draw from the priors that solves the model:\n'); 489 disp(params); 490 fprintf('Identification now continues for this draw.'); 491 parameters = 'Random_prior_params'; 492 parameters_TeX = 'Random prior parameter draw'; 493 end 494 end 495 ide_hess_point.params = params; 496 % save all output into identification folder 497 save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident'); 498 save([IdentifDirectoryName '/' fname '_' parameters '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point', 'store_options_ident'); 499 % display results of identification analysis 500 disp_identification(params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident); 501 if ~options_ident.no_identification_strength && ~options_.nograph 502 % plot (i) identification strength and sensitivity measure based on the moment information matrix and (ii) plot advanced analysis graphs 503 plot_identification(params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, IdentifDirectoryName, parameters_TeX, name_tex); 504 end 505 506 if SampleSize > 1 507 % initializations for Monte Carlo Analysis 508 fprintf('\nMonte Carlo Testing\n'); 509 h = dyn_waitbar(0,'Monte Carlo identification checks ...'); 510 iteration = 0; % initialize counter for admissable draws 511 run_index = 0; % initialize counter for admissable draws after saving previous draws to file(s) 512 file_index = 0; % initialize counter for files (if options_.MaxNumberOfBytes is reached, we store results in files) 513 options_MC = options_ident; %store options structure for Monte Carlo analysis 514 options_MC.advanced = 0; %do not run advanced checking in a Monte Carlo analysis 515 options_ident.checks_via_subsets = 0; % for Monte Carlo analysis currently only identification_checks and not identification_checks_via_subsets is supported 516 else 517 iteration = 1; % iteration equals SampleSize and we are finished 518 pdraws = []; % to have output object otherwise map_ident.m may crash 519 end 520 while iteration < SampleSize 521 if external_sample 522 params = pdraws0(iteration+1,:); % loaded draws 523 else 524 params = prior_draw(); % new random draw from prior 525 end 526 options_ident.tittxt = []; % clear title text for graphs and figures 527 % run identification analysis 528 [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, ide_derivatives_info, info, options_MC] = ... 529 identification_analysis(params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore 530 531 if iteration==0 && info(1)==0 % preallocate storage in the first admissable run 532 delete([IdentifDirectoryName '/' fname '_identif_*.mat']) % delete previously saved results 533 MAX_RUNS_BEFORE_SAVE_TO_FILE = min(SampleSize,ceil(options_.MaxNumberOfBytes/(size(ide_reducedform.si_dREDUCEDFORM,1)*totparam_nbr)/8)); % set how many runs can be stored before we save to files 534 pdraws = zeros(SampleSize,totparam_nbr); % preallocate storage for draws in each row 535 536 % preallocate storage for dynamic model 537 STO_si_dDYNAMIC = zeros([size(ide_dynamic.si_dDYNAMIC, 1), modparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]); 538 STO_DYNAMIC = zeros(size(ide_dynamic.DYNAMIC, 1), SampleSize); 539 IDE_DYNAMIC.ind_dDYNAMIC = ide_dynamic.ind_dDYNAMIC; 540 IDE_DYNAMIC.in0 = zeros(SampleSize, modparam_nbr); 541 IDE_DYNAMIC.ind0 = zeros(SampleSize, modparam_nbr); 542 IDE_DYNAMIC.jweak = zeros(SampleSize, modparam_nbr); 543 IDE_DYNAMIC.jweak_pair = zeros(SampleSize, modparam_nbr*(modparam_nbr+1)/2); 544 IDE_DYNAMIC.cond = zeros(SampleSize, 1); 545 IDE_DYNAMIC.Mco = zeros(SampleSize, modparam_nbr); 546 547 % preallocate storage for reduced form 548 if ~options_MC.no_identification_reducedform 549 STO_si_dREDUCEDFORM = zeros([size(ide_reducedform.si_dREDUCEDFORM, 1), totparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]); 550 STO_REDUCEDFORM = zeros(size(ide_reducedform.REDUCEDFORM, 1), SampleSize); 551 IDE_REDUCEDFORM.ind_dREDUCEDFORM = ide_reducedform.ind_dREDUCEDFORM; 552 IDE_REDUCEDFORM.in0 = zeros(SampleSize, 1); 553 IDE_REDUCEDFORM.ind0 = zeros(SampleSize, totparam_nbr); 554 IDE_REDUCEDFORM.jweak = zeros(SampleSize, totparam_nbr); 555 IDE_REDUCEDFORM.jweak_pair = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2); 556 IDE_REDUCEDFORM.cond = zeros(SampleSize, 1); 557 IDE_REDUCEDFORM.Mco = zeros(SampleSize, totparam_nbr); 558 else 559 STO_si_dREDUCEDFORM = {}; 560 STO_REDUCEDFORM = {}; 561 IDE_REDUCEDFORM = {}; 562 end 563 564 % preallocate storage for moments 565 if ~options_MC.no_identification_moments 566 STO_si_dMOMENTS = zeros([size(ide_moments.si_dMOMENTS, 1), totparam_nbr, MAX_RUNS_BEFORE_SAVE_TO_FILE]); 567 STO_MOMENTS = zeros(size(ide_moments.MOMENTS, 1), SampleSize); 568 IDE_MOMENTS.ind_dMOMENTS = ide_moments.ind_dMOMENTS; 569 IDE_MOMENTS.in0 = zeros(SampleSize, 1); 570 IDE_MOMENTS.ind0 = zeros(SampleSize, totparam_nbr); 571 IDE_MOMENTS.jweak = zeros(SampleSize, totparam_nbr); 572 IDE_MOMENTS.jweak_pair = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2); 573 IDE_MOMENTS.cond = zeros(SampleSize, 1); 574 IDE_MOMENTS.Mco = zeros(SampleSize, totparam_nbr); 575 IDE_MOMENTS.S = zeros(SampleSize, min(8, totparam_nbr)); 576 IDE_MOMENTS.V = zeros(SampleSize, totparam_nbr, min(8, totparam_nbr)); 577 else 578 STO_si_dMOMENTS = {}; 579 STO_MOMENTS = {}; 580 IDE_MOMENTS = {}; 581 end 582 583 % preallocate storage for spectrum 584 if ~options_MC.no_identification_spectrum 585 STO_dSPECTRUM = zeros([size(ide_spectrum.dSPECTRUM, 1), size(ide_spectrum.dSPECTRUM, 2), MAX_RUNS_BEFORE_SAVE_TO_FILE]); 586 IDE_SPECTRUM.ind_dSPECTRUM = ide_spectrum.ind_dSPECTRUM; 587 IDE_SPECTRUM.in0 = zeros(SampleSize, 1); 588 IDE_SPECTRUM.ind0 = zeros(SampleSize, totparam_nbr); 589 IDE_SPECTRUM.jweak = zeros(SampleSize, totparam_nbr); 590 IDE_SPECTRUM.jweak_pair = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2); 591 IDE_SPECTRUM.cond = zeros(SampleSize, 1); 592 IDE_SPECTRUM.Mco = zeros(SampleSize, totparam_nbr); 593 else 594 STO_dSPECTRUM = {}; 595 IDE_SPECTRUM = {}; 596 end 597 598 % preallocate storage for minimal system 599 if ~options_MC.no_identification_minimal 600 STO_dMINIMAL = zeros([size(ide_minimal.dMINIMAL, 1), size(ide_minimal.dMINIMAL, 2), MAX_RUNS_BEFORE_SAVE_TO_FILE]); 601 IDE_MINIMAL.ind_dMINIMAL = ide_minimal.ind_dMINIMAL; 602 IDE_MINIMAL.in0 = zeros(SampleSize, 1); 603 IDE_MINIMAL.ind0 = zeros(SampleSize, totparam_nbr); 604 IDE_MINIMAL.jweak = zeros(SampleSize, totparam_nbr); 605 IDE_MINIMAL.jweak_pair = zeros(SampleSize, totparam_nbr*(totparam_nbr+1)/2); 606 IDE_MINIMAL.cond = zeros(SampleSize, 1); 607 IDE_MINIMAL.Mco = zeros(SampleSize, totparam_nbr); 608 else 609 STO_dMINIMAL = {}; 610 IDE_MINIMAL = {}; 611 end 612 end 613 614 if info(1)==0 % if admissable draw 615 iteration = iteration + 1; %increase total index of admissable draws 616 run_index = run_index + 1; %increase index of admissable draws after saving to files 617 pdraws(iteration,:) = params; % store draw 618 619 % store results for steady state and dynamic model derivatives 620 STO_DYNAMIC(:,iteration) = ide_dynamic.DYNAMIC; 621 STO_si_dDYNAMIC(:,:,run_index) = ide_dynamic.si_dDYNAMIC; 622 IDE_DYNAMIC.cond(iteration,1) = ide_dynamic.cond; 623 IDE_DYNAMIC.ino(iteration,1) = ide_dynamic.ino; 624 IDE_DYNAMIC.ind0(iteration,:) = ide_dynamic.ind0; 625 IDE_DYNAMIC.jweak(iteration,:) = ide_dynamic.jweak; 626 IDE_DYNAMIC.jweak_pair(iteration,:) = ide_dynamic.jweak_pair; 627 IDE_DYNAMIC.Mco(iteration,:) = ide_dynamic.Mco; 628 629 % store results for reduced form 630 if ~options_MC.no_identification_reducedform 631 STO_REDUCEDFORM(:,iteration) = ide_reducedform.REDUCEDFORM; 632 STO_si_dREDUCEDFORM(:,:,run_index) = ide_reducedform.si_dREDUCEDFORM; 633 IDE_REDUCEDFORM.cond(iteration,1) = ide_reducedform.cond; 634 IDE_REDUCEDFORM.ino(iteration,1) = ide_reducedform.ino; 635 IDE_REDUCEDFORM.ind0(iteration,:) = ide_reducedform.ind0; 636 IDE_REDUCEDFORM.jweak(iteration,:) = ide_reducedform.jweak; 637 IDE_REDUCEDFORM.jweak_pair(iteration,:) = ide_reducedform.jweak_pair; 638 IDE_REDUCEDFORM.Mco(iteration,:) = ide_reducedform.Mco; 639 end 640 641 % store results for moments 642 if ~options_MC.no_identification_moments 643 STO_MOMENTS(:,iteration) = ide_moments.MOMENTS; 644 STO_si_dMOMENTS(:,:,run_index) = ide_moments.si_dMOMENTS; 645 IDE_MOMENTS.cond(iteration,1) = ide_moments.cond; 646 IDE_MOMENTS.ino(iteration,1) = ide_moments.ino; 647 IDE_MOMENTS.ind0(iteration,:) = ide_moments.ind0; 648 IDE_MOMENTS.jweak(iteration,:) = ide_moments.jweak; 649 IDE_MOMENTS.jweak_pair(iteration,:) = ide_moments.jweak_pair; 650 IDE_MOMENTS.Mco(iteration,:) = ide_moments.Mco; 651 IDE_MOMENTS.S(iteration,:) = ide_moments.S; 652 IDE_MOMENTS.V(iteration,:,:) = ide_moments.V; 653 end 654 655 % store results for spectrum 656 if ~options_MC.no_identification_spectrum 657 STO_dSPECTRUM(:,:,run_index) = ide_spectrum.dSPECTRUM; 658 IDE_SPECTRUM.cond(iteration,1) = ide_spectrum.cond; 659 IDE_SPECTRUM.ino(iteration,1) = ide_spectrum.ino; 660 IDE_SPECTRUM.ind0(iteration,:) = ide_spectrum.ind0; 661 IDE_SPECTRUM.jweak(iteration,:) = ide_spectrum.jweak; 662 IDE_SPECTRUM.jweak_pair(iteration,:) = ide_spectrum.jweak_pair; 663 IDE_SPECTRUM.Mco(iteration,:) = ide_spectrum.Mco; 664 end 665 666 % store results for minimal system 667 if ~options_MC.no_identification_minimal 668 STO_dMINIMAL(:,:,run_index) = ide_minimal.dMINIMAL; 669 IDE_MINIMAL.cond(iteration,1) = ide_minimal.cond; 670 IDE_MINIMAL.ino(iteration,1) = ide_minimal.ino; 671 IDE_MINIMAL.ind0(iteration,:) = ide_minimal.ind0; 672 IDE_MINIMAL.jweak(iteration,:) = ide_minimal.jweak; 673 IDE_MINIMAL.jweak_pair(iteration,:) = ide_minimal.jweak_pair; 674 IDE_MINIMAL.Mco(iteration,:) = ide_minimal.Mco; 675 end 676 677 % save results to file: either to avoid running into memory issues, i.e. (run_index==MAX_RUNS_BEFORE_SAVE_TO_FILE) or if finished (iteration==SampleSize) 678 if run_index==MAX_RUNS_BEFORE_SAVE_TO_FILE || iteration==SampleSize 679 file_index = file_index + 1; 680 if run_index < MAX_RUNS_BEFORE_SAVE_TO_FILE 681 %we are finished (iteration == SampleSize), so get rid of additional storage 682 STO_si_dDYNAMIC = STO_si_dDYNAMIC(:,:,1:run_index); 683 if ~options_MC.no_identification_reducedform 684 STO_si_dREDUCEDFORM = STO_si_dREDUCEDFORM(:,:,1:run_index); 685 end 686 if ~options_MC.no_identification_moments 687 STO_si_dMOMENTS = STO_si_dMOMENTS(:,:,1:run_index); 688 end 689 if ~options_MC.no_identification_spectrum 690 STO_dSPECTRUM = STO_dSPECTRUM(:,:,1:run_index); 691 end 692 if ~options_MC.no_identification_minimal 693 STO_dMINIMAL = STO_dMINIMAL(:,:,1:run_index); 694 end 695 end 696 save([IdentifDirectoryName '/' fname '_identif_' int2str(file_index) '.mat'], 'STO_si_dDYNAMIC'); 697 STO_si_dDYNAMIC = zeros(size(STO_si_dDYNAMIC)); % reset storage 698 if ~options_MC.no_identification_reducedform 699 save([IdentifDirectoryName '/' fname '_identif_' int2str(file_index) '.mat'], 'STO_si_dREDUCEDFORM', '-append'); 700 STO_si_dREDUCEDFORM = zeros(size(STO_si_dREDUCEDFORM)); % reset storage 701 end 702 if ~options_MC.no_identification_moments 703 save([IdentifDirectoryName '/' fname '_identif_' int2str(file_index) '.mat'], 'STO_si_dMOMENTS','-append'); 704 STO_si_dMOMENTS = zeros(size(STO_si_dMOMENTS)); % reset storage 705 end 706 if ~options_MC.no_identification_spectrum 707 save([IdentifDirectoryName '/' fname '_identif_' int2str(file_index) '.mat'], 'STO_dSPECTRUM','-append'); 708 STO_dSPECTRUM = zeros(size(STO_dSPECTRUM)); % reset storage 709 end 710 if ~options_MC.no_identification_minimal 711 save([IdentifDirectoryName '/' fname '_identif_' int2str(file_index) '.mat'], 'STO_dMINIMAL','-append'); 712 STO_dMINIMAL = zeros(size(STO_dMINIMAL)); % reset storage 713 end 714 run_index = 0; % reset index 715 end 716 if SampleSize > 1 717 dyn_waitbar(iteration/SampleSize, h, ['MC identification checks ', int2str(iteration), '/', int2str(SampleSize)]); 718 end 719 end 720 end 721 722 if SampleSize > 1 723 dyn_waitbar_close(h); 724 normalize_STO_DYNAMIC = std(STO_DYNAMIC,0,2); 725 if ~options_MC.no_identification_reducedform 726 normalize_STO_REDUCEDFORM = std(STO_REDUCEDFORM,0,2); 727 end 728 if ~options_MC.no_identification_moments 729 normalize_STO_MOMENTS = std(STO_MOMENTS,0,2); 730 end 731 if ~options_MC.no_identification_minimal 732 normalize_STO_MINIMAL = 1; %not used (yet) 733 end 734 if ~options_MC.no_identification_spectrum 735 normalize_STO_SPECTRUM = 1; %not used (yet) 736 end 737 normaliz1 = std(pdraws); 738 iter = 0; 739 for ifile_index = 1:file_index 740 load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dDYNAMIC'); 741 maxrun_dDYNAMIC = size(STO_si_dDYNAMIC,3); 742 if ~options_MC.no_identification_reducedform 743 load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dREDUCEDFORM'); 744 maxrun_dREDUCEDFORM = size(STO_si_dREDUCEDFORM,3); 745 else 746 maxrun_dREDUCEDFORM = 0; 747 end 748 if ~options_MC.no_identification_moments 749 load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dMOMENTS'); 750 maxrun_dMOMENTS = size(STO_si_dMOMENTS,3); 751 else 752 maxrun_dMOMENTS = 0; 753 end 754 if ~options_MC.no_identification_spectrum 755 load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_dSPECTRUM'); 756 maxrun_dSPECTRUM = size(STO_dSPECTRUM,3); 757 else 758 maxrun_dSPECTRUM = 0; 759 end 760 if ~options_MC.no_identification_minimal 761 load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_dMINIMAL'); 762 maxrun_dMINIMAL = size(STO_dMINIMAL,3); 763 else 764 maxrun_dMINIMAL = 0; 765 end 766 for irun=1:max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]) 767 iter=iter+1; 768 % note that this is not the same si_dDYNAMICnorm as computed in identification_analysis 769 % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure 770 si_dDYNAMICnorm(iter,:) = vnorm(STO_si_dDYNAMIC(:,:,irun)./repmat(normalize_STO_DYNAMIC,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end); 771 if ~options_MC.no_identification_reducedform 772 % note that this is not the same si_dREDUCEDFORMnorm as computed in identification_analysis 773 % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure 774 si_dREDUCEDFORMnorm(iter,:) = vnorm(STO_si_dREDUCEDFORM(:,:,irun)./repmat(normalize_STO_REDUCEDFORM,1,totparam_nbr)).*normaliz1; 775 end 776 if ~options_MC.no_identification_moments 777 % note that this is not the same si_dMOMENTSnorm as computed in identification_analysis 778 % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure 779 si_dMOMENTSnorm(iter,:) = vnorm(STO_si_dMOMENTS(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1; 780 end 781 if ~options_MC.no_identification_spectrum 782 % note that this is not the same dSPECTRUMnorm as computed in identification_analysis 783 dSPECTRUMnorm(iter,:) = vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used 784 end 785 if ~options_MC.no_identification_minimal 786 % note that this is not the same dMINIMALnorm as computed in identification_analysis 787 dMINIMALnorm(iter,:) = vnorm(STO_dMINIMAL(:,:,irun)); %not yet used 788 end 789 end 790 end 791 IDE_DYNAMIC.si_dDYNAMICnorm = si_dDYNAMICnorm; 792 save([IdentifDirectoryName '/' fname '_identif.mat'], 'pdraws', 'IDE_DYNAMIC','STO_DYNAMIC','-append'); 793 if ~options_MC.no_identification_reducedform 794 IDE_REDUCEDFORM.si_dREDUCEDFORMnorm = si_dREDUCEDFORMnorm; 795 save([IdentifDirectoryName '/' fname '_identif.mat'], 'IDE_REDUCEDFORM', 'STO_REDUCEDFORM','-append'); 796 end 797 if ~options_MC.no_identification_moments 798 IDE_MOMENTS.si_dMOMENTSnorm = si_dMOMENTSnorm; 799 save([IdentifDirectoryName '/' fname '_identif.mat'], 'IDE_MOMENTS', 'STO_MOMENTS','-append'); 800 end 801 802 end 803 804else 805 %% load previous analysis 806 load([IdentifDirectoryName '/' fname '_identif']); 807 parameters = store_options_ident.parameter_set; 808 options_ident.parameter_set = parameters; 809 options_ident.prior_mc = size(pdraws,1); 810 SampleSize = options_ident.prior_mc; 811 options_.options_ident = options_ident; 812end 813 814%% if dynare_identification is called as it own function (not through identification command) and if we load files 815if nargout>3 && iload 816 filnam = dir([IdentifDirectoryName '/' fname '_identif_*.mat']); 817 STO_si_dDYNAMIC = []; 818 STO_si_dREDUCEDFORM = []; 819 STO_si_dMOMENTS = []; 820 STO_dSPECTRUM = []; 821 STO_dMINIMAL = []; 822 for j=1:length(filnam) 823 load([IdentifDirectoryName '/' fname '_identif_',int2str(j),'.mat']); 824 STO_si_dDYNAMIC = cat(3,STO_si_dDYNAMIC, STO_si_dDYNAMIC(:,abs(iload),:)); 825 if ~options_ident.no_identification_reducedform 826 STO_si_dREDUCEDFORM = cat(3,STO_si_dREDUCEDFORM, STO_si_dREDUCEDFORM(:,abs(iload),:)); 827 end 828 if ~options_ident.no_identification_moments 829 STO_si_dMOMENTS = cat(3,STO_si_dMOMENTS, STO_si_dMOMENTS(:,abs(iload),:)); 830 end 831 if ~options_ident.no_identification_spectrum 832 STO_dSPECTRUM = cat(3,STO_dSPECTRUM, STO_dSPECTRUM(:,abs(iload),:)); 833 end 834 if ~options_ident.no_identification_minimal 835 STO_dMINIMAL = cat(3,STO_dMINIMAL, STO_dMINIMAL(:,abs(iload),:)); 836 end 837 end 838end 839 840if iload 841 %if previous analysis is loaded 842 fprintf(['Testing %s\n',parameters]); 843 disp_identification(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident); 844 if ~options_.nograph 845 % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs 846 plot_identification(ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, IdentifDirectoryName, [], name_tex); 847 end 848end 849 850%displaying and plotting of results for MC sample 851if SampleSize > 1 852 fprintf('\nTesting MC sample\n'); 853 %print results to console but make sure advanced=0 854 advanced0 = options_ident.advanced; 855 options_ident.advanced = 0; 856 disp_identification(pdraws, IDE_REDUCEDFORM, IDE_MOMENTS, IDE_SPECTRUM, IDE_MINIMAL, name, options_ident); 857 options_ident.advanced = advanced0; % reset advanced setting 858 if ~options_.nograph 859 % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs 860 plot_identification(pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, IdentifDirectoryName, [], name_tex); 861 end 862 %advanced display and plots for MC Sample, i.e. look at draws with highest/lowest condition number 863 if options_ident.advanced 864 jcrit = find(IDE_MOMENTS.ino); 865 if length(jcrit) < SampleSize 866 if isempty(jcrit) 867 % Make sure there is no overflow of plots produced (these are saved to the disk anyways) 868 store_nodisplay = options_.nodisplay; 869 options_.nodisplay = 1; 870 % HIGHEST CONDITION NUMBER 871 [~, jmax] = max(IDE_MOMENTS.cond); 872 tittxt = 'Draw with HIGHEST condition number'; 873 fprintf('\nTesting %s.\n',tittxt); 874 if ~iload 875 options_ident.tittxt = tittxt; %title text for graphs and figures 876 [ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max, options_ident] = ... 877 identification_analysis(pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables 878 save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_dynamic_max', 'jmax', '-append'); 879 end 880 advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on 881 disp_identification(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident); 882 options_ident.advanced = advanced0; %reset advanced setting 883 if ~options_.nograph 884 % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs 885 plot_identification(pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex); 886 end 887 888 % SMALLEST condition number 889 [~, jmin] = min(IDE_MOMENTS.cond); 890 tittxt = 'Draw with SMALLEST condition number'; 891 fprintf('Testing %s.\n',tittxt); 892 if ~iload 893 options_ident.tittxt = tittxt; %title text for graphs and figures 894 [ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min, options_ident] = ... 895 identification_analysis(pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables 896 save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append'); 897 end 898 advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on 899 disp_identification(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident); 900 options_ident.advanced = advanced0; %reset advanced setting 901 if ~options_.nograph 902 % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs 903 plot_identification(pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,IdentifDirectoryName,tittxt,name_tex); 904 end 905 % reset nodisplay option 906 options_.nodisplay = store_nodisplay; 907 else 908 % Make sure there is no overflow of plots produced (these are saved to the disk anyways) 909 store_nodisplay = options_.nodisplay; 910 options_.nodisplay = 1; 911 for j=1:length(jcrit) 912 tittxt = ['Rank deficient draw nr. ',int2str(j)]; 913 fprintf('\nTesting %s.\n',tittxt); 914 if ~iload 915 options_ident.tittxt = tittxt; %title text for graphs and figures 916 [ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve, options_ident] = ... 917 identification_analysis(pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); 918 end 919 advanced0 = options_ident.advanced; options_ident.advanced = 1; %make sure advanced setting is on 920 disp_identification(pdraws(jcrit(j),:), ide_reducedform_(j), ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), name, options_ident); 921 options_ident.advanced = advanced0; % reset advanced 922 if ~options_.nograph 923 % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs 924 plot_identification(pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex); 925 end 926 end 927 if ~iload 928 save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_', 'ide_moments_', 'ide_reducedform_', 'ide_dynamic_', 'ide_spectrum_', 'ide_minimal_', 'jcrit', '-append'); 929 end 930 % reset nodisplay option 931 options_.nodisplay = store_nodisplay; 932 end 933 end 934 end 935end 936 937%reset warning state 938if isoctave 939 warning('on') 940else 941 warning on 942end 943 944fprintf('\n==== Identification analysis completed ====\n\n') 945 946options_ = store_options_; %restore options set