dynare_identification.m 23.6 KB
Newer Older
ratto's avatar
ratto committed
1
function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_ident, pdraws0)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%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
19
20
21

% main 
%
Sébastien Villemot's avatar
Sébastien Villemot committed
22
% Copyright (C) 2010-2013 Dynare Team
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
%
% 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_

41
if isoctave
42
    warning('off'),
Marco Ratto's avatar
Marco Ratto committed
43
else
44
    warning off,
Marco Ratto's avatar
Marco Ratto committed
45
46
end

47
fname_ = M_.fname;
48
49
50
51
52
if ~isfield(M_,'dname'),
    M_.dname = M_.fname;
end
options_ident = set_default_option(options_ident,'gsa_sample_file',0);
options_ident = set_default_option(options_ident,'parameter_set','prior_mean');
53
options_ident = set_default_option(options_ident,'load_ident_files',0);
54
options_ident = set_default_option(options_ident,'useautocorr',0);
55
options_ident = set_default_option(options_ident,'ar',1);
56
options_ident = set_default_option(options_ident,'prior_mc',1);
57
options_ident = set_default_option(options_ident,'prior_range',0);
58
59
60
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);
61
options_ident = set_default_option(options_ident,'normalize_jacobians',1);
62
options_ident = set_default_option(options_ident,'lik_init',1);
63
options_ident = set_default_option(options_ident,'analytic_derivation',1);
64
65
66
67
68
69
70
71
72
73
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

74
if options_ident.gsa_sample_file,
Marco Ratto's avatar
Marco Ratto committed
75
    GSAFolder = checkpath('gsa',M_.dname);
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    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;
93
elseif nargin==1,
94
95
96
97
    pdraws0=[];
end
external_sample=0;
if nargin==2 || ~isempty(pdraws0),
98
    options_ident.prior_mc=size(pdraws0,1);
99
100
    options_ident.load_ident_files = 0;
    external_sample=1;
101
end
102
103
104
105
106
107
if isempty(estim_params_),
    options_ident.prior_mc=1;
    options_ident.prior_range=0;
    prior_exist=0;
else
    prior_exist=1;
108
    parameters = options_ident.parameter_set;
109
end
110

Marco Ratto's avatar
Marco Ratto committed
111
% options_ident.load_ident_files=1;
112
iload = options_ident.load_ident_files;
Marco Ratto's avatar
Marco Ratto committed
113
%options_ident.advanced=1;
114
advanced = options_ident.advanced;
115
nlags = options_ident.ar;
116
117
periods = options_ident.periods;
replic = options_ident.replic;
118
useautocorr = options_ident.useautocorr;
Marco Ratto's avatar
Marco Ratto committed
119
options_.order=1;
120
121
122
options_.ar=nlags;
options_.prior_mc = options_ident.prior_mc;
options_.options_ident = options_ident;
Marco Ratto's avatar
Marco Ratto committed
123
options_.Schur_vec_tol = 1.e-8;
124
options_.nomoments=0;
Marco Ratto's avatar
Marco Ratto committed
125
options_.analytic_derivation=1;
126

127
options_ = set_default_option(options_,'datafile','');
128
options_.mode_compute = 0;
129
options_.plot_priors = 0;
130
options_.smoother=1;
131
[dataset_,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
132
options_ident.analytic_derivation_mode = options_.analytic_derivation_mode;
133
134
135
136
137
138
if isempty(dataset_),
    dataset_.info.ntobs = periods;
    dataset_.info.nvobs = rows(options_.varobs);
    dataset_.info.varobs = options_.varobs;
    dataset_.rawdata = [];
    dataset_.missing.state = 0;
139
140
141
142
    for jdata=1:periods,
        temp1{jdata}=[1:dataset_.info.nvobs]';
    end
    dataset_.missing.aindex = temp1;
143
144
    dataset_.missing.vindex = [];
    dataset_.missing.number_of_observations = [];
145
    dataset_.missing.no_more_missing_observations = 1;
146
147
148
149
150
151
152
153
154
    dataset_.descriptive.mean = [];
    dataset_.data = [];

%     data_info.gend = periods;
%     data_info.data = [];
%     data_info.data_index = [];
%     data_info.number_of_observations = periods*size(options_.varobs,1);
%     data_info.no_more_missing_observations = 0;
%     data_info.missing_value = 0;
155
end
156
157
158

% results = prior_sampler(0,M_,bayestopt_,options_,oo_);

159
if prior_exist
160
    if any(bayestopt_.pshape > 0)
161
162
163
164
165
        if options_ident.prior_range
            prior_draw(1,1);
        else
            prior_draw(1);
        end
166
    else
167
        options_ident.prior_mc=1;
168
    end
169
end
170

171
172
SampleSize = options_ident.prior_mc;

Marco Ratto's avatar
Marco Ratto committed
173
if ~(exist('sylvester3','file')==2),
174
175
176
177
178

    dynareroot = strrep(which('dynare'),'dynare.m','');
    addpath([dynareroot 'gensylv'])
end

179
IdentifDirectoryName = CheckPath('identification',M_.dname);
180
if prior_exist,
181

182
183
184
185
    indx = [];
    if ~isempty(estim_params_.param_vals),
        indx = estim_params_.param_vals(:,1);
    end
186
187
188
189
190
191
192
193
    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;
194
    name_tex = char(M_.exo_names_tex(indexo,:),M_.param_names_tex(indx,:));
195

196
197
198
199
200
201
202
203
204
205
206
    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)];
207
    name_tex = [cellstr(M_.exo_names_tex); cellstr(M_.param_names_tex)];
208
end
209

210
211
212
213
214
215
216
217
218
219
220
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

Marco Ratto's avatar
Marco Ratto committed
221
options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,nparam-1]));
222
options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,nparam-1]);
Marco Ratto's avatar
Marco Ratto committed
223
224


225
MaxNumberOfBytes=options_.MaxNumberOfBytes;
226
store_options_ident = options_ident;
227

228
if iload <=0,
229
    
Marco Ratto's avatar
Marco Ratto committed
230
231
    [I,J]=find(M_.lead_lag_incidence');
    if prior_exist,
232
233
234
235
236
237
238
239
240
241
242
243
%         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
Marco Ratto's avatar
Marco Ratto committed
244
        params = set_prior(estim_params_,M_,options_)';
245
246
247
        if all(bayestopt_.pshape == 0)
            parameters = 'ML_Starting_value';
            disp('Testing ML Starting value')
248
        else
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
        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
274
        end
Marco Ratto's avatar
Marco Ratto committed
275
276
    else
        params = [sqrt(diag(M_.Sigma_e))', M_.params'];
277
        parameters = 'Current_params';
278
        disp('Testing current parameter values')
Marco Ratto's avatar
Marco Ratto committed
279
    end
280
    [idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
281
        identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
282
    if info(1)~=0,
283
        skipline()
284
285
286
        disp('----------- ')
        disp('Parameter error:')
        disp(['The model does not solve for ', parameters, ' with error code info = ', int2str(info(1))]),
287
        skipline()
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
        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
315
        disp('----------- ')
316
        skipline()
317
318
319
320
321
322
323
324
325
326
327
        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] = ...
                    identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
            end
        end
        if info(1)
328
            skipline()
329
330
331
332
333
334
            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('----------- ')
335
336
            skipline()
            return
337
        end
338
    else
339
340
341
342
    idehess_point.params=params;
%     siH = idemodel_point.siH;
%     siJ = idemoments_point.siJ;
%     siLRE = idelre_point.siLRE;
Marco Ratto's avatar
Marco Ratto committed
343
344
345
%     normH = max(abs(siH)')';
%     normJ = max(abs(siJ)')';
%     normLRE = max(abs(siLRE)')';
346
    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_point', 'idemoments_point','idemodel_point', 'idelre_point','store_options_ident')
347
    disp_identification(params, idemodel_point, idemoments_point, name, advanced);
348
349
350
    if ~options_.nograph,
        plot_identification(params,idemoments_point,idehess_point,idemodel_point,idelre_point,advanced,parameters,name,IdentifDirectoryName);
    end
351
    end
352

353
    if SampleSize > 1,
354
        skipline()
Marco Ratto's avatar
Marco Ratto committed
355
        disp('Monte Carlo Testing')
356
        h = dyn_waitbar(0,'Monte Carlo identification checks ...');
Marco Ratto's avatar
Marco Ratto committed
357
358
359
360
361
362
363
364
365
        iteration = 0;
        loop_indx = 0;
        file_index = 0;
        run_index = 0;
        options_MC=options_ident;
        options_MC.advanced=0;
    else
        iteration = 1;
        pdraws = [];
366
    end
367
368
    while iteration < SampleSize,
        loop_indx = loop_indx+1;
369
        if external_sample,
Marco Ratto's avatar
Marco Ratto committed
370
            params = pdraws0(iteration+1,:);
371
        else
Marco Ratto's avatar
Marco Ratto committed
372
373
            params = prior_draw();
        end
374
        [dum1, ideJ, ideH, ideGP, dum2 , info] = ...
375
            identification_analysis(params,indx,indexo,options_MC,dataset_, prior_exist, name_tex,0);
376
        if iteration==0 && info(1)==0,
Marco Ratto's avatar
Marco Ratto committed
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
            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);
403
404
            idemoments.S = zeros(SampleSize,min(8,nparam));
            idemoments.V = zeros(SampleSize,nparam,min(8,nparam));
Marco Ratto's avatar
Marco Ratto committed
405
            delete([IdentifDirectoryName '/' M_.fname '_identif_*.mat'])
406
        end
407
        if info(1)==0,
408
409
            iteration = iteration + 1;
            run_index = run_index + 1;
Marco Ratto's avatar
Marco Ratto committed
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
            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;
431
432
            idemoments.S(iteration,:)=ideJ.S;
            idemoments.V(iteration,:,:)=ideJ.V;
Marco Ratto's avatar
Marco Ratto committed
433
434
435
436
437
438
439
440
441
442
            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);
443
                end
Marco Ratto's avatar
Marco Ratto committed
444
445
446
447
448
                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));
449
                
Marco Ratto's avatar
Marco Ratto committed
450
451
452
            end
            
            if SampleSize > 1,
453
%                 if isoctave || options_.console_mode,
454
455
456
457
%                     console_waitbar(0,iteration/SampleSize);
%                 else
                    dyn_waitbar(iteration/SampleSize,h,['MC identification checks ',int2str(iteration),'/',int2str(SampleSize)])
%                 end
458
459
            end
        end
Marco Ratto's avatar
Marco Ratto committed
460
        
461
    end
462
463
    
    
464
    if SampleSize > 1,
465
        if isoctave || options_.console_mode,
466
467
468
469
470
            fprintf('\n');
            diary on;
        else
            close(h),
        end
Marco Ratto's avatar
Marco Ratto committed
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
        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
Marco Ratto's avatar
Marco Ratto committed
486
487
488
489
490
491
        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
492
493
494
        siJnorm = idemoments_point.siJnorm;
        siHnorm = idemodel_point.siHnorm;
        siLREnorm = idelre_point.siLREnorm;
495
    end
496
    
497
else
Marco Ratto's avatar
Marco Ratto committed
498
    load([IdentifDirectoryName '/' M_.fname '_identif'])
499
500
501
%     identFiles = dir([IdentifDirectoryName '/' M_.fname '_identif_*']);
    parameters = store_options_ident.parameter_set;
    options_ident.parameter_set = parameters;
502
    options_ident.prior_mc=size(pdraws,1);
503
    SampleSize = options_ident.prior_mc;
504
    options_.options_ident = options_ident;
ratto's avatar
ratto committed
505
end  
506

507
if nargout>3 && iload,
508
509
510
    filnam = dir([IdentifDirectoryName '/' M_.fname '_identif_*.mat']);
    H=[];
    JJ = [];
511
    gp = [];
512
513
514
515
    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),:));
ratto's avatar
ratto committed
516
        gp = cat(3,gp, stoLRE(:,abs(iload),:));
517
        
518
519
520
    end
end

Marco Ratto's avatar
Marco Ratto committed
521
if iload,
522
    disp(['Testing ',parameters])
523
    disp_identification(idehess_point.params, idemodel_point, idemoments_point, name,advanced);
524
525
526
    if ~options_.nograph,
        plot_identification(idehess_point.params,idemoments_point,idehess_point,idemodel_point,idelre_point,advanced,parameters,name,IdentifDirectoryName);
    end
527
end
528
if SampleSize > 1,
529
    fprintf('\n')
Marco Ratto's avatar
Marco Ratto committed
530
531
    disp('Testing MC sample')
    disp_identification(pdraws, idemodel, idemoments, name);
532
533
534
    if ~options_.nograph,
        plot_identification(pdraws,idemoments,idehess_point,idemodel,idelre,advanced,'MC sample - ',name, IdentifDirectoryName);
    end
Marco Ratto's avatar
Marco Ratto committed
535
536
    if advanced,
        jcrit=find(idemoments.ino);
537
538
        if length(jcrit)<SampleSize,
            if isempty(jcrit),
539
                [dum,jmax]=max(idemoments.cond);
540
541
542
                fprintf('\n')
                tittxt = 'Draw with HIGHEST condition number';
                fprintf('\n')
543
                disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
544
545
                if ~iload,
                    [idehess_max, idemoments_max, idemodel_max, idelre_max, derivatives_info_max] = ...
546
                        identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
547
548
                    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_max', 'idemoments_max','idemodel_max', 'idelre_max', 'jmax', '-append');
                end
549
                disp_identification(pdraws(jmax,:), idemodel_max, idemoments_max, name,1);
550
                close all,
551
552
553
                if ~options_.nograph,
                    plot_identification(pdraws(jmax,:),idemoments_max,idehess_max,idemodel_max,idelre_max,1,tittxt,name,IdentifDirectoryName);
                end
554
                [dum,jmin]=min(idemoments.cond);
555
556
557
                fprintf('\n')
                tittxt = 'Draw with SMALLEST condition number';
                fprintf('\n')
558
                disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
559
560
                if ~iload,
                    [idehess_min, idemoments_min, idemodel_min, idelre_min, derivatives_info_min] = ...
561
                        identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
562
563
                    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_min', 'idemoments_min','idemodel_min', 'idelre_min', 'jmin', '-append');
                end
564
                disp_identification(pdraws(jmin,:), idemodel_min, idemoments_min, name,1);
565
                close all,
566
567
568
                if ~options_.nograph,
                    plot_identification(pdraws(jmin,:),idemoments_min,idehess_min,idemodel_min,idelre_min,1,tittxt,name,IdentifDirectoryName);
                end
569
570
571
572
            else
                for j=1:length(jcrit),
                    tittxt = ['Rank deficient draw n. ',int2str(j)];
                    fprintf('\n')
573
                    disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
574
575
                    if ~iload,
                        [idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j)] = ...
576
                            identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
577
                    end
578
                    disp_identification(pdraws(jcrit(j),:), idemodel_(j), idemoments_(j), name,1);
579
                    close all,
580
581
582
                    if ~options_.nograph,
                        plot_identification(pdraws(jcrit(j),:),idemoments_(j),idehess_(j),idemodel_(j),idelre_(j),1,tittxt,name,IdentifDirectoryName);
                    end
583
584
585
586
                end
                if ~iload,
                    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_', 'idemoments_','idemodel_', 'idelre_', 'jcrit', '-append');
                end
587
            end
588
        end
589
    end
590
end
Marco Ratto's avatar
Marco Ratto committed
591

592
if isoctave
593
    warning('on'),
Marco Ratto's avatar
Marco Ratto committed
594
else
595
    warning on,
Marco Ratto's avatar
Marco Ratto committed
596
end
597

598
skipline()
599
disp(['==== Identification analysis completed ====' ]),
600
skipline(2)