dynare_identification.m 23.8 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
338
        else
            parameters = 'Random_prior_params';
339
        end
340
    else
341
342
343
344
    idehess_point.params=params;
%     siH = idemodel_point.siH;
%     siJ = idemoments_point.siJ;
%     siLRE = idelre_point.siLRE;
Marco Ratto's avatar
Marco Ratto committed
345
346
347
%     normH = max(abs(siH)')';
%     normJ = max(abs(siJ)')';
%     normLRE = max(abs(siLRE)')';
348
    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_point', 'idemoments_point','idemodel_point', 'idelre_point','store_options_ident')
349
    save([IdentifDirectoryName '/' M_.fname '_' parameters '_identif.mat'], 'idehess_point', 'idemoments_point','idemodel_point', 'idelre_point','store_options_ident')
350
    disp_identification(params, idemodel_point, idemoments_point, name, advanced);
351
352
353
    if ~options_.nograph,
        plot_identification(params,idemoments_point,idehess_point,idemodel_point,idelre_point,advanced,parameters,name,IdentifDirectoryName);
    end
354
    end
355

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

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

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

595
if isoctave
596
    warning('on'),
Marco Ratto's avatar
Marco Ratto committed
597
else
598
    warning on,
Marco Ratto's avatar
Marco Ratto committed
599
end
600

601
skipline()
602
disp(['==== Identification analysis completed ====' ]),
603
skipline(2)