2010-01-15 10:57:05 +01:00
function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification ( options_ident, pdraws0)
2011-02-28 16:55:21 +01:00
%function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_ident, pdraws0)
%
% INPUTS
% o options_ident [structure] identification options
% o pdraws0 [matrix] optional: matrix of MC sample of model params.
%
% OUTPUTS
% o pdraws [matrix] matrix of MC sample of model params used
% o TAU, [matrix] MC sample of entries in the model solution (stacked vertically)
% o GAM, [matrix] MC sample of entries in the moments (stacked vertically)
% o LRE, [matrix] MC sample of entries in LRE model (stacked vertically)
% o gp, [matrix] derivatives of the Jacobian (LRE model)
% o H, [matrix] derivatives of the model solution
% o JJ [matrix] derivatives of the moments
%
% SPECIAL REQUIREMENTS
% None
2009-12-16 18:17:34 +01:00
% main
%
2013-06-12 16:42:09 +02:00
% Copyright (C) 2010-2013 Dynare Team
2009-12-16 18:17:34 +01:00
%
% 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 <http://www.gnu.org/licenses/>.
global M_ options_ oo_ bayestopt_ estim_params_
2013-11-04 10:54:45 +01:00
if isoctave
2011-03-18 11:05:40 +01:00
warning ( ' off' ) ,
2010-02-23 10:47:43 +01:00
else
2011-03-18 11:05:40 +01:00
warning off ,
2010-02-23 10:47:43 +01:00
end
2011-02-28 16:55:21 +01:00
fname_ = M_ . fname ;
2011-05-04 10:39:30 +02:00
options_ident = set_default_option ( options_ident , ' gsa_sample_file' , 0 ) ;
options_ident = set_default_option ( options_ident , ' parameter_set' , ' prior_mean' ) ;
2009-12-16 18:17:34 +01:00
options_ident = set_default_option ( options_ident , ' load_ident_files' , 0 ) ;
2010-03-08 16:12:01 +01:00
options_ident = set_default_option ( options_ident , ' useautocorr' , 0 ) ;
2011-06-15 23:47:27 +02:00
options_ident = set_default_option ( options_ident , ' ar' , 1 ) ;
2010-09-29 16:12:48 +02:00
options_ident = set_default_option ( options_ident , ' prior_mc' , 1 ) ;
2010-03-08 16:12:01 +01:00
options_ident = set_default_option ( options_ident , ' prior_range' , 0 ) ;
2010-09-29 16:12:48 +02:00
options_ident = set_default_option ( options_ident , ' periods' , 300 ) ;
options_ident = set_default_option ( options_ident , ' replic' , 100 ) ;
options_ident = set_default_option ( options_ident , ' advanced' , 0 ) ;
2011-09-01 11:18:24 +02:00
options_ident = set_default_option ( options_ident , ' normalize_jacobians' , 1 ) ;
2015-06-05 16:39:17 +02:00
%Deal with non-stationary cases
if isfield ( options_ident , ' diffuse_filter' ) %set lik_init and options_
options_ident . lik_init = 3 ; %overwrites any other lik_init set
options_ . diffuse_filter = options_ident . diffuse_filter ; %make options_ inherit diffuse filter; will overwrite any conflicting lik_init in dynare_estimation_init
else
if options_ . diffuse_filter == 1 %warning if estimation with diffuse filter was done, but not passed
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.' )
end
if isfield ( options_ident , ' lik_init' )
if options_ident . lik_init == 3 %user specified diffuse filter using the lik_init option
options_ident . analytic_derivation = 0 ; %diffuse filter not compatible with analytic derivation
end
end
end
2011-06-23 10:54:06 +02:00
options_ident = set_default_option ( options_ident , ' lik_init' , 1 ) ;
2012-07-05 10:14:10 +02:00
options_ident = set_default_option ( options_ident , ' analytic_derivation' , 1 ) ;
2015-06-05 16:39:17 +02:00
2012-10-02 08:59:49 +02:00
if isfield ( options_ident , ' nograph' ) ,
options_ . nograph = options_ident . nograph ;
end
if isfield ( options_ident , ' nodisplay' ) ,
options_ . nodisplay = options_ident . nodisplay ;
end
if isfield ( options_ident , ' graph_format' ) ,
options_ . graph_format = options_ident . graph_format ;
end
2015-06-05 16:39:17 +02:00
if isfield ( options_ident , ' prior_trunc' ) ,
options_ . prior_trunc = options_ident . prior_trunc ;
end
2012-10-02 08:59:49 +02:00
2011-05-04 10:39:30 +02:00
if options_ident . gsa_sample_file ,
2012-04-04 10:35:18 +02:00
GSAFolder = checkpath ( ' gsa' , M_ . dname ) ;
2011-05-04 10:39:30 +02:00
if options_ident . gsa_sample_file == 1 ,
load ( [ GSAFolder , filesep , fname_ , ' _prior' ] , ' lpmat' , ' lpmat0' , ' istable' ) ;
elseif options_ident . gsa_sample_file == 2 ,
load ( [ GSAFolder , filesep , fname_ , ' _mc' ] , ' lpmat' , ' lpmat0' , ' istable' ) ;
else
load ( [ GSAFolder , filesep , options_ident . gsa_sample_file ] , ' lpmat' , ' lpmat0' , ' istable' ) ;
end
if isempty ( istable ) ,
istable = 1 : size ( lpmat , 1 ) ;
end
if ~ isempty ( lpmat0 ) ,
lpmatx = lpmat0 ( istable , : ) ;
else
lpmatx = [ ] ;
end
pdraws0 = [ lpmatx lpmat ( istable , : ) ] ;
clear lpmat lpmat0 istable ;
2011-10-24 18:45:12 +02:00
elseif nargin == 1 ,
2011-05-04 10:39:30 +02:00
pdraws0 = [ ] ;
end
external_sample = 0 ;
if nargin == 2 || ~ isempty ( pdraws0 ) ,
2009-12-16 18:17:34 +01:00
options_ident . prior_mc = size ( pdraws0 , 1 ) ;
2011-05-04 10:39:30 +02:00
options_ident . load_ident_files = 0 ;
external_sample = 1 ;
2009-12-16 18:17:34 +01:00
end
2010-09-29 16:12:48 +02:00
if isempty ( estim_params_ ) ,
options_ident . prior_mc = 1 ;
options_ident . prior_range = 0 ;
prior_exist = 0 ;
else
prior_exist = 1 ;
2011-05-04 10:39:30 +02:00
parameters = options_ident . parameter_set ;
2011-02-28 16:55:21 +01:00
end
2009-12-16 18:17:34 +01:00
2011-05-02 11:14:29 +02:00
% options_ident.load_ident_files=1;
2009-12-16 18:17:34 +01:00
iload = options_ident . load_ident_files ;
2011-05-02 11:14:29 +02:00
%options_ident.advanced=1;
2010-09-29 16:12:48 +02:00
advanced = options_ident . advanced ;
2009-12-16 18:17:34 +01:00
nlags = options_ident . ar ;
2010-09-29 16:12:48 +02:00
periods = options_ident . periods ;
replic = options_ident . replic ;
2009-12-16 18:17:34 +01:00
useautocorr = options_ident . useautocorr ;
2012-04-04 10:35:18 +02:00
options_ . order = 1 ;
2009-12-16 18:17:34 +01:00
options_ . ar = nlags ;
options_ . prior_mc = options_ident . prior_mc ;
options_ . options_ident = options_ident ;
2010-02-23 10:47:43 +01:00
options_ . Schur_vec_tol = 1.e-8 ;
2011-06-06 15:28:46 +02:00
options_ . nomoments = 0 ;
2015-06-05 16:39:17 +02:00
options_ = set_default_option ( options_ , ' analytic_derivation' , 1 ) ;
2012-09-25 15:49:59 +02:00
options_ = set_default_option ( options_ , ' datafile' , ' ' ) ;
2009-12-16 18:17:34 +01:00
options_ . mode_compute = 0 ;
2011-02-28 16:55:21 +01:00
options_ . plot_priors = 0 ;
2013-03-19 16:56:10 +01:00
options_ . smoother = 1 ;
2014-06-17 12:04:58 +02:00
[ dataset_ , dataset_info , xparam1 , hh , M_ , options_ , oo_ , estim_params_ , bayestopt_ ] = dynare_estimation_init ( M_ . endo_names , fname_ , 1 , M_ , options_ , oo_ , estim_params_ , bayestopt_ ) ;
2012-07-05 10:14:10 +02:00
options_ident . analytic_derivation_mode = options_ . analytic_derivation_mode ;
2009-12-16 18:17:34 +01:00
2011-02-28 16:55:21 +01:00
if prior_exist
2013-03-19 16:56:10 +01:00
if any ( bayestopt_ . pshape > 0 )
2011-06-06 16:09:55 +02:00
if options_ident . prior_range
prior_draw ( 1 , 1 ) ;
else
prior_draw ( 1 ) ;
end
2011-02-28 16:55:21 +01:00
else
2011-06-06 16:09:55 +02:00
options_ident . prior_mc = 1 ;
2011-02-28 16:55:21 +01:00
end
2010-03-08 16:12:01 +01:00
end
2010-09-03 15:23:05 +02:00
2011-06-06 16:09:55 +02:00
SampleSize = options_ident . prior_mc ;
2012-04-04 16:43:44 +02:00
if ~ ( exist ( ' sylvester3' , ' file' ) == 2 ) ,
2009-12-16 18:17:34 +01:00
dynareroot = strrep ( which ( ' dynare' ) , ' dynare.m' , ' ' ) ;
addpath ( [ dynareroot ' gensylv' ] )
end
2011-12-15 11:56:23 +01:00
IdentifDirectoryName = CheckPath ( ' identification' , M_ . dname ) ;
2010-09-29 16:12:48 +02:00
if prior_exist ,
2009-12-16 18:17:34 +01:00
2011-03-18 11:05:40 +01:00
indx = [ ] ;
if ~ isempty ( estim_params_ . param_vals ) ,
indx = estim_params_ . param_vals ( : , 1 ) ;
end
2010-09-29 16:12:48 +02:00
indexo = [ ] ;
if ~ isempty ( estim_params_ . var_exo )
indexo = estim_params_ . var_exo ( : , 1 ) ;
end
nparam = length ( bayestopt_ . name ) ;
np = estim_params_ . np ;
name = bayestopt_ . name ;
2011-03-18 11:05:40 +01:00
name_tex = char ( M_ . exo_names_tex ( indexo , : ) , M_ . param_names_tex ( indx , : ) ) ;
2009-12-16 18:17:34 +01:00
2010-09-29 16:12:48 +02:00
offset = estim_params_ . nvx ;
offset = offset + estim_params_ . nvn ;
offset = offset + estim_params_ . ncx ;
offset = offset + estim_params_ . ncn ;
else
indx = [ 1 : M_ . param_nbr ] ;
indexo = [ 1 : M_ . exo_nbr ] ;
offset = M_ . exo_nbr ;
np = M_ . param_nbr ;
nparam = np + offset ;
name = [ cellstr ( M_ . exo_names ) ; cellstr ( M_ . param_names ) ] ;
2011-03-18 11:05:40 +01:00
name_tex = [ cellstr ( M_ . exo_names_tex ) ; cellstr ( M_ . param_names_tex ) ] ;
2010-09-29 16:12:48 +02:00
end
2009-12-16 18:17:34 +01:00
2013-12-06 08:01:01 +01:00
skipline ( )
disp ( [ ' ==== Identification analysis ====' ] ) ,
skipline ( )
if nparam < 2 ,
options_ident . advanced = 0 ;
advanced = options_ident . advanced ;
disp ( ' There is only one parameter to study for identitification.' )
disp ( ' The advanced option is re-set to 0.' )
skipline ( )
end
2011-05-02 11:14:29 +02:00
options_ident = set_default_option ( options_ident , ' max_dim_cova_group' , min ( [ 2 , nparam - 1 ] ) ) ;
2011-05-30 14:39:27 +02:00
options_ident . max_dim_cova_group = min ( [ options_ident . max_dim_cova_group , nparam - 1 ] ) ;
2011-04-12 12:08:14 +02:00
2009-12-16 18:17:34 +01:00
MaxNumberOfBytes = options_ . MaxNumberOfBytes ;
2011-05-04 10:39:30 +02:00
store_options_ident = options_ident ;
2009-12-16 18:17:34 +01:00
2011-02-28 16:55:21 +01:00
if iload < = 0 ,
2009-12-16 18:17:34 +01:00
2011-05-02 11:14:29 +02:00
[ I , J ] = find ( M_ . lead_lag_incidence ' ) ;
if prior_exist ,
2011-06-06 16:09:55 +02:00
% if exist([fname_,'_mean.mat'],'file'),
% % disp('Testing posterior mean')
% load([fname_,'_mean'],'xparam1')
% pmean = xparam1';
% clear xparam1
% end
% if exist([fname_,'_mode.mat'],'file'),
% % disp('Testing posterior mode')
% load([fname_,'_mode'],'xparam1')
% pmode = xparam1';
% clear xparam1
% end
2011-05-02 11:14:29 +02:00
params = set_prior ( estim_params_ , M_ , options_ ) ' ;
2013-03-19 16:56:10 +01:00
if all ( bayestopt_ . pshape == 0 )
parameters = ' ML_Starting_value' ;
disp ( ' Testing ML Starting value' )
2011-06-06 16:09:55 +02:00
else
2011-05-04 10:39:30 +02:00
switch parameters
case ' posterior_mode'
disp ( ' Testing posterior mode' )
params ( 1 , : ) = get_posterior_parameters ( ' mode' ) ;
case ' posterior_mean'
disp ( ' Testing posterior mean' )
params ( 1 , : ) = get_posterior_parameters ( ' mean' ) ;
case ' posterior_median'
disp ( ' Testing posterior median' )
params ( 1 , : ) = get_posterior_parameters ( ' median' ) ;
case ' prior_mode'
disp ( ' Testing prior mode' )
params ( 1 , : ) = bayestopt_ . p5 ( : ) ;
case ' prior_mean'
disp ( ' Testing prior mean' )
params ( 1 , : ) = bayestopt_ . p1 ;
otherwise
disp ( ' The option parameter_set has to be equal to:' )
disp ( ' ' ' posterior_mode' ' , ' )
disp ( ' ' ' posterior_mean' ' , ' )
disp ( ' ' ' posterior_median' ' , ' )
disp ( ' ' ' prior_mode' ' or' )
disp ( ' ' ' prior_mean' ' .' )
error ;
end
2011-06-06 16:09:55 +02:00
end
2011-05-02 11:14:29 +02:00
else
params = [ sqrt ( diag ( M_ . Sigma_e ) ) ' , M_ . params ' ] ;
2011-05-04 10:39:30 +02:00
parameters = ' Current_params' ;
2011-06-06 16:09:55 +02:00
disp ( ' Testing current parameter values' )
2011-05-02 11:14:29 +02:00
end
2012-08-24 17:09:13 +02:00
[ idehess_point , idemoments_point , idemodel_point , idelre_point , derivatives_info_point , info ] = ...
2014-06-17 12:04:58 +02:00
identification_analysis ( params , indx , indexo , options_ident , dataset_ , dataset_info , prior_exist , name_tex , 1 ) ;
2012-08-24 17:09:13 +02:00
if info ( 1 ) ~= 0 ,
2013-07-10 17:12:34 +02:00
skipline ( )
2012-08-24 17:09:13 +02:00
disp ( ' ----------- ' )
disp ( ' Parameter error:' )
disp ( [ ' The model does not solve for ' , parameters , ' with error code info = ' , int2str ( info ( 1 ) ) ] ) ,
2013-07-10 17:12:34 +02:00
skipline ( )
2013-03-19 16:56:10 +01:00
if info ( 1 ) == 1 ,
disp ( ' info==1 %! The model doesn' ' t determine the current variables uniquely.' )
elseif info ( 1 ) == 2 ,
disp ( ' info==2 %! MJDGGES returned an error code.' )
elseif info ( 1 ) == 3 ,
disp ( ' info==3 %! Blanchard & Kahn conditions are not satisfied: no stable equilibrium. ' )
elseif info ( 1 ) == 4 ,
disp ( ' info==4 %! Blanchard & Kahn conditions are not satisfied: indeterminacy. ' )
elseif info ( 1 ) == 5 ,
disp ( ' info==5 %! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure. ' )
elseif info ( 1 ) == 6 ,
disp ( ' info==6 %! The jacobian evaluated at the deterministic steady state is complex.' )
elseif info ( 1 ) == 19 ,
disp ( ' info==19 %! The steadystate routine thrown an exception (inconsistent deep parameters). ' )
elseif info ( 1 ) == 20 ,
disp ( ' info==20 %! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations). ' )
elseif info ( 1 ) == 21 ,
disp ( ' info==21 %! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.' )
elseif info ( 1 ) == 22 ,
disp ( ' info==22 %! The steady has NaNs. ' )
elseif info ( 1 ) == 23 ,
disp ( ' info==23 %! M_.params has been updated in the steadystate routine and has complex valued scalars. ' )
elseif info ( 1 ) == 24 ,
disp ( ' info==24 %! M_.params has been updated in the steadystate routine and has some NaNs. ' )
elseif info ( 1 ) == 30 ,
disp ( ' info==30 %! Ergodic variance can' ' t be computed. ' )
end
2012-08-24 17:09:13 +02:00
disp ( ' ----------- ' )
2013-07-10 17:12:34 +02:00
skipline ( )
2013-03-19 16:56:10 +01:00
if any ( bayestopt_ . pshape )
disp ( ' Try sampling up to 50 parameter sets from the prior.' )
kk = 0 ;
while kk < 50 && info ( 1 ) ,
kk = kk + 1 ;
params = prior_draw ( ) ;
[ idehess_point , idemoments_point , idemodel_point , idelre_point , derivatives_info_point , info ] = ...
2014-06-18 15:37:13 +02:00
identification_analysis ( params , indx , indexo , options_ident , dataset_ , dataset_info , prior_exist , name_tex , 1 ) ;
2013-03-19 16:56:10 +01:00
end
end
if info ( 1 )
2013-07-10 17:12:34 +02:00
skipline ( )
2013-03-19 16:56:10 +01:00
disp ( ' ----------- ' )
disp ( ' Identification stopped:' )
if any ( bayestopt_ . pshape )
disp ( ' The model did not solve for any of 50 attempts of random samples from the prior' )
end
disp ( ' ----------- ' )
2013-07-10 17:12:34 +02:00
skipline ( )
return
2013-12-10 11:05:54 +01:00
else
parameters = ' Random_prior_params' ;
2013-03-19 16:56:10 +01:00
end
2014-07-23 12:10:00 +02:00
end
2011-05-04 10:39:30 +02:00
idehess_point . params = params ;
% siH = idemodel_point.siH;
% siJ = idemoments_point.siJ;
% siLRE = idelre_point.siLRE;
2011-05-02 11:14:29 +02:00
% normH = max(abs(siH)')';
% normJ = max(abs(siJ)')';
% normLRE = max(abs(siLRE)')';
2011-05-04 10:39:30 +02:00
save ( [ IdentifDirectoryName ' /' M_ . fname ' _identif.mat' ] , ' idehess_point' , ' idemoments_point' , ' idemodel_point' , ' idelre_point' , ' store_options_ident' )
2013-12-10 11:05:54 +01:00
save ( [ IdentifDirectoryName ' /' M_ . fname ' _' parameters ' _identif.mat' ] , ' idehess_point' , ' idemoments_point' , ' idemodel_point' , ' idelre_point' , ' store_options_ident' )
2012-06-15 13:45:26 +02:00
disp_identification ( params , idemodel_point , idemoments_point , name , advanced ) ;
2012-06-15 13:07:26 +02:00
if ~ options_ . nograph ,
plot_identification ( params , idemoments_point , idehess_point , idemodel_point , idelre_point , advanced , parameters , name , IdentifDirectoryName ) ;
end
2009-12-16 18:17:34 +01:00
2010-09-29 16:12:48 +02:00
if SampleSize > 1 ,
2013-07-10 17:12:34 +02:00
skipline ( )
2011-05-02 11:14:29 +02:00
disp ( ' Monte Carlo Testing' )
2011-10-27 13:09:26 +02:00
h = dyn_waitbar ( 0 , ' Monte Carlo identification checks ...' ) ;
2011-05-02 11:14:29 +02:00
iteration = 0 ;
loop_indx = 0 ;
file_index = 0 ;
run_index = 0 ;
options_MC = options_ident ;
options_MC . advanced = 0 ;
else
iteration = 1 ;
pdraws = [ ] ;
2010-09-29 16:12:48 +02:00
end
2009-12-16 18:17:34 +01:00
while iteration < SampleSize ,
loop_indx = loop_indx + 1 ;
2011-05-04 10:39:30 +02:00
if external_sample ,
2011-05-02 11:14:29 +02:00
params = pdraws0 ( iteration + 1 , : ) ;
2010-09-29 16:12:48 +02:00
else
2011-05-02 11:14:29 +02:00
params = prior_draw ( ) ;
end
2011-06-21 09:43:24 +02:00
[ dum1 , ideJ , ideH , ideGP , dum2 , info ] = ...
2014-06-18 15:37:13 +02:00
identification_analysis ( params , indx , indexo , options_MC , dataset_ , dataset_info , prior_exist , name_tex , 0 ) ;
2012-07-15 22:15:06 +02:00
if iteration == 0 && info ( 1 ) == 0 ,
2011-05-02 11:14:29 +02:00
MAX_tau = min ( SampleSize , ceil ( MaxNumberOfBytes / ( size ( ideH . siH , 1 ) * nparam ) / 8 ) ) ;
stoH = zeros ( [ size ( ideH . siH , 1 ) , nparam , MAX_tau ] ) ;
stoJJ = zeros ( [ size ( ideJ . siJ , 1 ) , nparam , MAX_tau ] ) ;
stoLRE = zeros ( [ size ( ideGP . siLRE , 1 ) , np , MAX_tau ] ) ;
TAU = zeros ( size ( ideH . siH , 1 ) , SampleSize ) ;
GAM = zeros ( size ( ideJ . siJ , 1 ) , SampleSize ) ;
LRE = zeros ( size ( ideGP . siLRE , 1 ) , SampleSize ) ;
pdraws = zeros ( SampleSize , nparam ) ;
idemoments . indJJ = ideJ . indJJ ;
idemodel . indH = ideH . indH ;
idelre . indLRE = ideGP . indLRE ;
idemoments . ind0 = zeros ( SampleSize , nparam ) ;
idemodel . ind0 = zeros ( SampleSize , nparam ) ;
idelre . ind0 = zeros ( SampleSize , np ) ;
idemoments . jweak = zeros ( SampleSize , nparam ) ;
idemodel . jweak = zeros ( SampleSize , nparam ) ;
idelre . jweak = zeros ( SampleSize , np ) ;
idemoments . jweak_pair = zeros ( SampleSize , nparam * ( nparam + 1 ) / 2 ) ;
idemodel . jweak_pair = zeros ( SampleSize , nparam * ( nparam + 1 ) / 2 ) ;
idelre . jweak_pair = zeros ( SampleSize , np * ( np + 1 ) / 2 ) ;
idemoments . cond = zeros ( SampleSize , 1 ) ;
idemodel . cond = zeros ( SampleSize , 1 ) ;
idelre . cond = zeros ( SampleSize , 1 ) ;
idemoments . Mco = zeros ( SampleSize , nparam ) ;
idemodel . Mco = zeros ( SampleSize , nparam ) ;
idelre . Mco = zeros ( SampleSize , np ) ;
2011-05-10 10:10:42 +02:00
idemoments . S = zeros ( SampleSize , min ( 8 , nparam ) ) ;
idemoments . V = zeros ( SampleSize , nparam , min ( 8 , nparam ) ) ;
2011-05-02 11:14:29 +02:00
delete ( [ IdentifDirectoryName ' /' M_ . fname ' _identif_*.mat' ] )
2010-09-29 16:12:48 +02:00
end
2009-12-16 18:17:34 +01:00
if info ( 1 ) == 0 ,
2011-04-15 15:28:58 +02:00
iteration = iteration + 1 ;
run_index = run_index + 1 ;
2011-05-02 11:14:29 +02:00
TAU ( : , iteration ) = ideH . TAU ;
LRE ( : , iteration ) = ideGP . LRE ;
GAM ( : , iteration ) = ideJ . GAM ;
idemoments . cond ( iteration , 1 ) = ideJ . cond ;
idemodel . cond ( iteration , 1 ) = ideH . cond ;
idelre . cond ( iteration , 1 ) = ideGP . cond ;
idemoments . ino ( iteration , 1 ) = ideJ . ino ;
idemodel . ino ( iteration , 1 ) = ideH . ino ;
idelre . ino ( iteration , 1 ) = ideGP . ino ;
idemoments . ind0 ( iteration , : ) = ideJ . ind0 ;
idemodel . ind0 ( iteration , : ) = ideH . ind0 ;
idelre . ind0 ( iteration , : ) = ideGP . ind0 ;
idemoments . jweak ( iteration , : ) = ideJ . jweak ;
idemodel . jweak ( iteration , : ) = ideH . jweak ;
idelre . jweak ( iteration , : ) = ideGP . jweak ;
idemoments . jweak_pair ( iteration , : ) = ideJ . jweak_pair ;
idemodel . jweak_pair ( iteration , : ) = ideH . jweak_pair ;
idelre . jweak_pair ( iteration , : ) = ideGP . jweak_pair ;
idemoments . Mco ( iteration , : ) = ideJ . Mco ;
idemodel . Mco ( iteration , : ) = ideH . Mco ;
idelre . Mco ( iteration , : ) = ideGP . Mco ;
2011-05-10 10:10:42 +02:00
idemoments . S ( iteration , : ) = ideJ . S ;
idemoments . V ( iteration , : , : ) = ideJ . V ;
2011-05-02 11:14:29 +02:00
stoLRE ( : , : , run_index ) = ideGP . siLRE ;
stoH ( : , : , run_index ) = ideH . siH ;
stoJJ ( : , : , run_index ) = ideJ . siJ ;
pdraws ( iteration , : ) = params ;
if run_index == MAX_tau || iteration == SampleSize ,
file_index = file_index + 1 ;
if run_index < MAX_tau ,
stoH = stoH ( : , : , 1 : run_index ) ;
stoJJ = stoJJ ( : , : , 1 : run_index ) ;
stoLRE = stoLRE ( : , : , 1 : run_index ) ;
2009-12-16 18:17:34 +01:00
end
2011-05-02 11:14:29 +02:00
save ( [ IdentifDirectoryName ' /' M_ . fname ' _identif_' int2str ( file_index ) ' .mat' ] , ' stoH' , ' stoJJ' , ' stoLRE' )
run_index = 0 ;
stoH = zeros ( size ( stoH ) ) ;
stoJJ = zeros ( size ( stoJJ ) ) ;
stoLRE = zeros ( size ( stoLRE ) ) ;
2011-02-28 16:55:21 +01:00
2011-05-02 11:14:29 +02:00
end
if SampleSize > 1 ,
2013-11-04 10:54:45 +01:00
% if isoctave || options_.console_mode,
2011-10-27 13:09:26 +02:00
% console_waitbar(0,iteration/SampleSize);
% else
dyn_waitbar ( iteration / SampleSize , h , [ ' MC identification checks ' , int2str ( iteration ) , ' /' , int2str ( SampleSize ) ] )
% end
2009-12-16 18:17:34 +01:00
end
end
2011-05-02 11:14:29 +02:00
2009-12-16 18:17:34 +01:00
end
2011-02-28 16:55:21 +01:00
2010-09-29 16:12:48 +02:00
if SampleSize > 1 ,
2013-11-04 10:54:45 +01:00
if isoctave || options_ . console_mode ,
2011-10-27 13:09:26 +02:00
fprintf ( ' \n' ) ;
diary on ;
else
close ( h ) ,
end
2011-04-12 12:08:14 +02:00
normTAU = std ( TAU ' ) ' ;
normLRE = std ( LRE ' ) ' ;
normGAM = std ( GAM ' ) ' ;
normaliz1 = std ( pdraws ) ;
iter = 0 ;
for ifile_index = 1 : file_index ,
load ( [ IdentifDirectoryName ' /' M_ . fname ' _identif_' int2str ( ifile_index ) ' .mat' ] , ' stoH' , ' stoJJ' , ' stoLRE' )
for irun = 1 : size ( stoH , 3 ) ,
iter = iter + 1 ;
siJnorm ( iter , : ) = vnorm ( stoJJ ( : , : , irun ) ./ repmat ( normGAM , 1 , nparam ) ) .* normaliz1 ;
siHnorm ( iter , : ) = vnorm ( stoH ( : , : , irun ) ./ repmat ( normTAU , 1 , nparam ) ) .* normaliz1 ;
siLREnorm ( iter , : ) = vnorm ( stoLRE ( : , : , irun ) ./ repmat ( normLRE , 1 , nparam - offset ) ) .* normaliz1 ( offset + 1 : end ) ;
end
end
2011-05-02 11:14:29 +02:00
idemoments . siJnorm = siJnorm ;
idemodel . siHnorm = siHnorm ;
idelre . siLREnorm = siLREnorm ;
save ( [ IdentifDirectoryName ' /' M_ . fname ' _identif.mat' ] , ' pdraws' , ' idemodel' , ' idemoments' , ' idelre' , ... %'indJJ', 'indH', 'indLRE', ...
' TAU' , ' GAM' , ' LRE' , ' -append' )
else
2011-05-04 10:39:30 +02:00
siJnorm = idemoments_point . siJnorm ;
siHnorm = idemodel_point . siHnorm ;
siLREnorm = idelre_point . siLREnorm ;
2010-09-29 16:12:48 +02:00
end
2011-02-28 16:55:21 +01:00
2009-12-16 18:17:34 +01:00
else
2011-05-02 11:14:29 +02:00
load ( [ IdentifDirectoryName ' /' M_ . fname ' _identif' ] )
2011-05-04 10:39:30 +02:00
% identFiles = dir([IdentifDirectoryName '/' M_.fname '_identif_*']);
parameters = store_options_ident . parameter_set ;
options_ident . parameter_set = parameters ;
2010-03-08 16:12:01 +01:00
options_ident . prior_mc = size ( pdraws , 1 ) ;
2011-02-28 16:55:21 +01:00
SampleSize = options_ident . prior_mc ;
2010-03-08 16:12:01 +01:00
options_ . options_ident = options_ident ;
2010-01-15 10:57:05 +01:00
end
2009-12-16 18:17:34 +01:00
2011-02-10 17:23:22 +01:00
if nargout > 3 && iload ,
2009-12-16 18:17:34 +01:00
filnam = dir ( [ IdentifDirectoryName ' /' M_ . fname ' _identif_*.mat' ] ) ;
H = [ ] ;
JJ = [ ] ;
2011-02-28 16:55:21 +01:00
gp = [ ] ;
2009-12-16 18:17:34 +01:00
for j = 1 : length ( filnam ) ,
load ( [ IdentifDirectoryName ' /' M_ . fname ' _identif_' , int2str ( j ) , ' .mat' ] ) ;
H = cat ( 3 , H , stoH ( : , abs ( iload ) , : ) ) ;
JJ = cat ( 3 , JJ , stoJJ ( : , abs ( iload ) , : ) ) ;
2010-01-15 10:57:05 +01:00
gp = cat ( 3 , gp , stoLRE ( : , abs ( iload ) , : ) ) ;
2011-02-28 16:55:21 +01:00
2009-12-16 18:17:34 +01:00
end
end
2011-05-02 11:14:29 +02:00
if iload ,
2011-05-04 10:39:30 +02:00
disp ( [ ' Testing ' , parameters ] )
2012-06-15 13:45:26 +02:00
disp_identification ( idehess_point . params , idemodel_point , idemoments_point , name , advanced ) ;
2012-06-15 13:07:26 +02:00
if ~ options_ . nograph ,
plot_identification ( idehess_point . params , idemoments_point , idehess_point , idemodel_point , idelre_point , advanced , parameters , name , IdentifDirectoryName ) ;
end
2011-04-12 18:18:18 +02:00
end
2010-09-29 16:12:48 +02:00
if SampleSize > 1 ,
2011-05-04 10:39:30 +02:00
fprintf ( ' \n' )
2011-05-02 11:14:29 +02:00
disp ( ' Testing MC sample' )
disp_identification ( pdraws , idemodel , idemoments , name ) ;
2012-06-15 13:07:26 +02:00
if ~ options_ . nograph ,
2014-03-10 20:57:29 +01:00
plot_identification ( pdraws , idemoments , idehess_point , idemodel , idelre , advanced , ' MC sample ' , name , IdentifDirectoryName ) ;
2012-06-15 13:07:26 +02:00
end
2011-05-02 11:14:29 +02:00
if advanced ,
jcrit = find ( idemoments . ino ) ;
2011-05-10 10:10:42 +02:00
if length ( jcrit ) < SampleSize ,
if isempty ( jcrit ) ,
2011-06-21 09:43:24 +02:00
[ dum , jmax ] = max ( idemoments . cond ) ;
2011-05-10 10:10:42 +02:00
fprintf ( ' \n' )
tittxt = ' Draw with HIGHEST condition number' ;
fprintf ( ' \n' )
2011-09-28 20:48:47 +02:00
disp ( [ ' Testing ' , tittxt , ' . Press ENTER' ] ) , pause ( 5 ) ,
2011-05-10 10:10:42 +02:00
if ~ iload ,
[ idehess_max , idemoments_max , idemodel_max , idelre_max , derivatives_info_max ] = ...
2014-06-27 18:40:09 +02:00
identification_analysis ( pdraws ( jmax , : ) , indx , indexo , options_ident , dataset_ , dataset_info , prior_exist , name_tex , 1 ) ;
2011-05-10 10:10:42 +02:00
save ( [ IdentifDirectoryName ' /' M_ . fname ' _identif.mat' ] , ' idehess_max' , ' idemoments_max' , ' idemodel_max' , ' idelre_max' , ' jmax' , ' -append' ) ;
end
2012-06-15 13:45:26 +02:00
disp_identification ( pdraws ( jmax , : ) , idemodel_max , idemoments_max , name , 1 ) ;
2011-05-30 14:39:27 +02:00
close all ,
2012-06-15 13:07:26 +02:00
if ~ options_ . nograph ,
plot_identification ( pdraws ( jmax , : ) , idemoments_max , idehess_max , idemodel_max , idelre_max , 1 , tittxt , name , IdentifDirectoryName ) ;
end
2011-06-21 09:43:24 +02:00
[ dum , jmin ] = min ( idemoments . cond ) ;
2011-05-10 10:10:42 +02:00
fprintf ( ' \n' )
tittxt = ' Draw with SMALLEST condition number' ;
fprintf ( ' \n' )
2011-09-28 20:48:47 +02:00
disp ( [ ' Testing ' , tittxt , ' . Press ENTER' ] ) , pause ( 5 ) ,
2011-05-10 10:10:42 +02:00
if ~ iload ,
[ idehess_min , idemoments_min , idemodel_min , idelre_min , derivatives_info_min ] = ...
2014-06-18 15:37:13 +02:00
identification_analysis ( pdraws ( jmin , : ) , indx , indexo , options_ident , dataset_ , dataset_info , prior_exist , name_tex , 1 ) ;
2011-05-10 10:10:42 +02:00
save ( [ IdentifDirectoryName ' /' M_ . fname ' _identif.mat' ] , ' idehess_min' , ' idemoments_min' , ' idemodel_min' , ' idelre_min' , ' jmin' , ' -append' ) ;
end
2012-06-15 13:45:26 +02:00
disp_identification ( pdraws ( jmin , : ) , idemodel_min , idemoments_min , name , 1 ) ;
2011-05-30 14:39:27 +02:00
close all ,
2012-06-15 13:07:26 +02:00
if ~ options_ . nograph ,
plot_identification ( pdraws ( jmin , : ) , idemoments_min , idehess_min , idemodel_min , idelre_min , 1 , tittxt , name , IdentifDirectoryName ) ;
end
2011-05-10 10:10:42 +02:00
else
for j = 1 : length ( jcrit ) ,
tittxt = [ ' Rank deficient draw n. ' , int2str ( j ) ] ;
fprintf ( ' \n' )
2011-09-28 20:48:47 +02:00
disp ( [ ' Testing ' , tittxt , ' . Press ENTER' ] ) , pause ( 5 ) ,
2011-05-10 10:10:42 +02:00
if ~ iload ,
[ idehess_ ( j ) , idemoments_ ( j ) , idemodel_ ( j ) , idelre_ ( j ) , derivatives_info_ ( j ) ] = ...
2014-06-18 15:37:13 +02:00
identification_analysis ( pdraws ( jcrit ( j ) , : ) , indx , indexo , options_ident , dataset_ , dataset_info , prior_exist , name_tex , 1 ) ;
2011-05-10 10:10:42 +02:00
end
2012-06-15 13:45:26 +02:00
disp_identification ( pdraws ( jcrit ( j ) , : ) , idemodel_ ( j ) , idemoments_ ( j ) , name , 1 ) ;
2011-05-30 14:39:27 +02:00
close all ,
2012-06-15 13:07:26 +02:00
if ~ options_ . nograph ,
plot_identification ( pdraws ( jcrit ( j ) , : ) , idemoments_ ( j ) , idehess_ ( j ) , idemodel_ ( j ) , idelre_ ( j ) , 1 , tittxt , name , IdentifDirectoryName ) ;
end
2011-05-10 10:10:42 +02:00
end
if ~ iload ,
save ( [ IdentifDirectoryName ' /' M_ . fname ' _identif.mat' ] , ' idehess_' , ' idemoments_' , ' idemodel_' , ' idelre_' , ' jcrit' , ' -append' ) ;
end
2011-03-18 11:05:40 +01:00
end
2011-04-15 15:28:58 +02:00
end
2011-03-18 11:05:40 +01:00
end
2010-09-29 16:12:48 +02:00
end
2010-02-23 10:47:43 +01:00
2013-11-04 10:54:45 +01:00
if isoctave
2011-03-18 11:05:40 +01:00
warning ( ' on' ) ,
2010-02-23 10:47:43 +01:00
else
2011-03-18 11:05:40 +01:00
warning on ,
2010-02-23 10:47:43 +01:00
end
2010-09-29 16:12:48 +02:00
2013-07-10 17:12:34 +02:00
skipline ( )
2010-09-29 16:12:48 +02:00
disp ( [ ' ==== Identification analysis completed ====' ] ) ,
2013-07-10 17:12:34 +02:00
skipline ( 2 )