metropolis_hastings_initialization.m 15.3 KB
Newer Older
1
function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
2
    metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin)
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% Metropolis-Hastings initialization.
% 
% INPUTS 
%   o TargetFun  [char]     string specifying the name of the objective
%                           function (posterior kernel).
%   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
%   o vv         [double]   (p*p) matrix, posterior covariance matrix (at the mode).
%   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters. 
%   o varargin              list of argument following mh_bounds
%  
% OUTPUTS 
%   None  
%
% SPECIAL REQUIREMENTS
%   None.
18

19
% Copyright (C) 2006-2011 Dynare Team
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
%
% 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/>.

36
37
global M_ options_ bayestopt_

38
39
40
41
42
43
44
45
46
47
48
49
50
ix2 = [];
ilogpo2 = [];
ModelName = []; 
MhDirectoryName = [];
fblck = [];
fline = [];
npar = [];
nblck = [];
nruns = [];
NewFile = [];
MAX_nruns = [];
d = [];

51
52
53
54
55
56
57
58
ModelName = M_.fname;

if ~isempty(M_.bvar)
    ModelName = [M_.fname '_bvar'];
end

bayestopt_.penalty = 1e8;

59
MhDirectoryName = CheckPath('metropolis',M_.dname);
60
61
62
63

nblck = options_.mh_nblck;
nruns = ones(nblck,1)*options_.mh_replic;
npar  = length(xparam1);
64
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
65
66
d = chol(vv);

67
if ~options_.load_mh_file && ~options_.mh_recover
68
69
70
71
72
73
74
75
76
77
78
    %% Here we start a new metropolis-hastings, previous draws are not
    %% considered.
    if nblck > 1
        disp('MH: Multiple chains mode.')
    else
        disp('MH: One Chain mode.')
    end
    %% Delete old mh files...
    files = dir([ MhDirectoryName '/' ModelName '_mh*_blck*.mat']);
    if length(files)
        delete([ MhDirectoryName '/' ModelName '_mh*_blck*.mat']);
79
        disp('MH: Old _mh files successfully erased!')
80
81
82
83
    end
    file = dir([ MhDirectoryName '/metropolis.log']);
    if length(file)
        delete([ MhDirectoryName '/metropolis.log']);
84
        disp('MH: Old metropolis.log file successfully erased!')
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
        disp('MH: Creation of a new metropolis.log file.')
    end
    fidlog = fopen([MhDirectoryName '/metropolis.log'],'w');
    fprintf(fidlog,'%% MH log file (Dynare).\n');
    fprintf(fidlog,['%% ' datestr(now,0) '.\n']);
    fprintf(fidlog,' \n\n');
    fprintf(fidlog,'%% Session 1.\n');
    fprintf(fidlog,' \n');
    fprintf(fidlog,['  Number of blocks...............: ' int2str(nblck) '\n']);
    fprintf(fidlog,['  Number of simulations per block: ' int2str(nruns(1)) '\n']);
    fprintf(fidlog,' \n');
    %% Initial values...
    if nblck > 1% Case 1: multiple chains
        fprintf(fidlog,['  Initial values of the parameters:\n']);
        disp('MH: Searching for initial values...')
        ix2 = zeros(nblck,npar);
        ilogpo2 = zeros(nblck,1);
        for j=1:nblck
103
104
            validate    = 0;
            init_iter   = 0;
105
            trial = 1;
106
            while validate == 0 && trial <= 10
107
                candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
108
                if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2)) 
109
110
                    ix2(j,:) = candidate;
                    ilogpo2(j) = - feval(TargetFun,ix2(j,:)',varargin{:});
111
112
113
114
115
116
117
118
119
120
                    if ilogpo2(j) <= - bayestopt_.penalty+1e-6
                        validate = 0;
                    else
                        fprintf(fidlog,['    Blck ' int2str(j) ':\n']);
                        for i=1:length(ix2(1,:))
                            fprintf(fidlog,['      params:' int2str(i) ': ' num2str(ix2(j,i)) '\n']);
                        end
                        fprintf(fidlog,['      logpo2: ' num2str(ilogpo2(j)) '\n']);
                        j = j+1;
                        validate = 1;
121
122
123
                    end
                end
                init_iter = init_iter + 1;
124
                if init_iter > 100 && validate == 0
125
126
127
128
129
130
131
                    disp(['MH: I couldn''t get a valid initial value in 100 trials.'])
                    disp(['MH: You should Reduce mh_init_scale...'])
                    disp(sprintf('MH: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale))
                    options_.mh_init_scale = input('MH: Enter a new value...  ');
                    trial = trial+1;
                end
            end
132
            if trial > 10 && ~validate
133
134
135
136
137
138
139
140
141
142
                disp(['MH: I''m unable to find a starting value for block ' int2str(j)])
                return
            end
        end
        fprintf(fidlog,' \n');
        disp('MH: Initial values found!')
        disp(' ')
    else% Case 2: one chain (we start from the posterior mode)
        fprintf(fidlog,['  Initial values of the parameters:\n']);
        candidate = transpose(xparam1);
143
        if all(candidate' > mh_bounds(:,1)) && all(candidate' < mh_bounds(:,2)) 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
            ix2 = candidate;
            ilogpo2 = - feval(TargetFun,ix2',varargin{:});
            disp('MH: Initialization at the posterior mode.')
            disp(' ')
            fprintf(fidlog,['    Blck ' int2str(1) 'params:\n']);
            for i=1:length(ix2(1,:))
                fprintf(fidlog,['      ' int2str(i)  ':' num2str(ix2(1,i)) '\n']);
            end
            fprintf(fidlog,['    Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']);
        else
            disp('MH: Initialization failed...')
            disp('MH: The posterior mode lies outside the prior bounds.')
            return
        end
        fprintf(fidlog,' \n');
    end
    fprintf(fidlog,' \n');
    fblck = 1;
    fline = ones(nblck,1);
    NewFile = ones(nblck,1);
    %%
    %% Creation of the mh-history file:
    %%
    file = dir([MhDirectoryName '/'  ModelName '_mh_history.mat']);
    if length(files)
        delete([ MhDirectoryName '/' ModelName '_mh_history.mat']);
170
        disp('MH: Old mh_history file successfully erased!')
171
172
173
174
175
176
177
178
179
    end
    AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
    AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
    record.Nblck = nblck;
    record.MhDraws = zeros(1,3);
    record.MhDraws(1,1) = nruns(1);
    record.MhDraws(1,2) = AnticipatedNumberOfFiles;
    record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
    record.AcceptationRates = zeros(1,nblck);
180
181
182
    % separate initializaton for each chain
    JSUM = 0;
    for j=1:nblck,
183
184
        record.Seeds(j).Normal = randn('state');
        record.Seeds(j).Unifor = rand('state');
185
    end
186
187
188
189
190
191
    record.InitialParameters = ix2;
    record.InitialLogLiK = ilogpo2;
    record.LastParameters = zeros(nblck,npar);
    record.LastLogLiK = zeros(nblck,1);
    record.LastFileNumber = AnticipatedNumberOfFiles ;
    record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
192
    save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');  
193
194
195
196
197
    fprintf(fidlog,['  CREATION OF THE MH HISTORY FILE!\n\n']);
    fprintf(fidlog,['    Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']);
    fprintf(fidlog,['    Expected number of lines in the last file: ' ...
                    int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
    fprintf(fidlog,['\n']);
198
    for j = 1:nblck,
199
200
201
202
203
204
205
206
        fprintf(fidlog,['    Initial seed (randn) for chain number ',int2str(j),':\n']);
        for i=1:length(record.Seeds(j).Normal)
            fprintf(fidlog,['      ' num2str(record.Seeds(j).Normal(i)') '\n']);
        end
        fprintf(fidlog,['    Initial seed (rand) for chain number ',int2str(j),':\n']);
        for i=1:length(record.Seeds(j).Unifor)
            fprintf(fidlog,['      ' num2str(record.Seeds(j).Unifor(i)') '\n']);
        end
207
    end,
208
209
    fprintf(fidlog,' \n');
    fclose(fidlog);
210
elseif options_.load_mh_file && ~options_.mh_recover
211
212
213
    %% Here we consider previous mh files (previous mh did not crash).
    disp('MH: I''m loading past metropolis-hastings simulations...')
    file = dir([ MhDirectoryName '/'  ModelName '_mh_history.mat' ]);
214
    files = dir([ MhDirectoryName filesep ModelName '_mh*.mat']);
215
216
217
218
219
220
221
222
    if ~length(files)
        disp('MH:: FAILURE! there is no MH file to load here!')
        return
    end
    if ~length(file)
        disp('MH:: FAILURE! there is no MH-history file!')
        return
    else
223
        load([ MhDirectoryName '/'  ModelName '_mh_history.mat'])
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
    end
    fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
    fprintf(fidlog,'\n');
    fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']);
    fprintf(fidlog,' \n');
    fprintf(fidlog,['  Number of blocks...............: ' int2str(nblck) '\n']);
    fprintf(fidlog,['  Number of simulations per block: ' int2str(nruns(1)) '\n']);
    fprintf(fidlog,' \n');
    past_number_of_blocks = record.Nblck;
    if past_number_of_blocks ~= nblck
        disp('MH:: The specified number of blocks doesn''t match with the previous number of blocks!')
        disp(['MH:: You declared ' int2str(nblck) ' blocks, but the previous number of blocks was ' ...
              int2str(past_number_of_blocks) '.'])
        disp(['MH:: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
        nblck = past_number_of_blocks;
        options_.mh_nblck = nblck;
    end
    % I read the last line of the last mh-file for initialization 
    % of the new metropolis-hastings simulations:
    LastFileNumber = record.LastFileNumber;
    LastLineNumber = record.LastLineNumber;
    if LastLineNumber < MAX_nruns
        NewFile = ones(nblck,1)*LastFileNumber;
    else
        NewFile = ones(nblck,1)*(LastFileNumber+1);
    end
    ilogpo2 = record.LastLogLiK;
    ix2 = record.LastParameters;
    fblck = 1;
    fline = ones(nblck,1)*(LastLineNumber+1);
    NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
    record.MhDraws = [record.MhDraws;zeros(1,3)];
    NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
    NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile;
    AnticipatedNumberOfFiles = ceil(NumberOfDrawsToBeSaved/MAX_nruns);
    AnticipatedNumberOfLinesInTheLastFile = NumberOfDrawsToBeSaved - (AnticipatedNumberOfFiles-1)*MAX_nruns;  
    record.LastFileNumber = LastFileNumber + AnticipatedNumberOfFiles;
    record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
    record.MhDraws(end,1) = nruns(1);
    record.MhDraws(end,2) = AnticipatedNumberOfFiles;
    record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
265
    save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
266
267
268
269
270
271
272
273
274
275
276
277
278
    disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
    disp(' ')
    fclose(fidlog);
elseif options_.mh_recover
    %% The previous metropolis-hastings crashed before the end! I try to
    %% recover the saved draws...
    disp('MH: Recover mode!')
    disp(' ')
    file = dir([MhDirectoryName '/'  ModelName '_mh_history.mat']);
    if ~length(file)
        disp('MH:: FAILURE! there is no MH-history file!')
        return
    else
279
        load([ MhDirectoryName '/'  ModelName '_mh_history.mat'])
280
    end
stepan's avatar
stepan committed
281
    nblck = record.Nblck;% Number of "parallel" mcmc chains.
282
283
284
285
286
287
288
289
290
291
292
293
294
295
    options_.mh_nblck = nblck;
    if size(record.MhDraws,1) == 1
        OldMh = 0;% The crashed metropolis was the first session.
    else
        OldMh = 1;% The crashed metropolis wasn't the first session.
    end
    %% Default initialization:
    if OldMh
        ilogpo2 = record.LastLogLiK;
        ix2 = record.LastParameters;
    else
        ilogpo2 = record.InitialLogLiK;
        ix2 = record.InitialParameters;
    end
stepan's avatar
stepan committed
296
    %% Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers.
297
    if OldMh
stepan's avatar
stepan committed
298
299
300
        LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session.
        LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions.
        if LastLineNumberInThePreviousMh < MAX_nruns% Test if the last mh files of the previous session are not complete (yes)
301
            NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh;
stepan's avatar
stepan committed
302
303
            fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1);
        else% The last mh files of the previous session are complete.
304
            NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1);
stepan's avatar
stepan committed
305
            fline = ones(nblck,1);
306
307
        end
    else
stepan's avatar
stepan committed
308
309
        LastLineNumberInThePreviousMh = 0;
        LastFileNumberInThePreviousMh = 0;
310
311
312
        NewFile = ones(nblck,1);
        fline = ones(nblck,1);
    end
stepan's avatar
stepan committed
313
    %% Set fblck (First block), an integer targeting the crashed mcmc chain.
314
315
316
317
318
319
320
    fblck = 1;
    %% How many mh files should we have ?
    ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
    ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
    %% I count the total number of saved mh files...
    AllMhFiles = dir([MhDirectoryName '/' ModelName '_mh*_blck*.mat']);
    TotalNumberOfMhFiles = length(AllMhFiles);
stepan's avatar
stepan committed
321
322
323
324
325
326
327
328
    %% And I quit if I can't find a crashed mcmc chain. 
    if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
        disp('MH: It appears that you don''t need to use the mh_recover option!')
        disp('    You have to edit the mod file and remove the mh_recover') 
        disp('    in the estimation command')
        error
    end
    %% I count the number of saved mh files per block.
329
    NumberOfMhFilesPerBlock = zeros(nblck,1);
stepan's avatar
stepan committed
330
331
332
    for b = 1:nblck
        BlckMhFiles = dir([ MhDirectoryName '/' ModelName '_mh*_blck' int2str(b) '.mat']);
        NumberOfMhFilesPerBlock(b) = length(BlckMhFiles);
333
334
    end
    %% Is there a chain with less mh files than expected ? 
stepan's avatar
stepan committed
335
336
337
    while fblck <= nblck
        if  NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock
            disp(['MH: Chain ' int2str(fblck) ' is not complete!'])
338
            break
stepan's avatar
stepan committed
339
            % The mh_recover session will start from chain fblck.
340
        else
stepan's avatar
stepan committed
341
            disp(['MH: Chain ' int2str(fblck) ' is complete!'])
342
        end
stepan's avatar
stepan committed
343
        fblck = fblck+1;
344
    end
345
    %% How many mh-files are saved in this block?
stepan's avatar
stepan committed
346
347
348
349
    NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck);
    %% Correct the number of saved mh files if the crashed metropolis was not the first session (so
    %% that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain 
    %% of the current session).  
350
    if OldMh
stepan's avatar
stepan committed
351
        NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
352
    end
stepan's avatar
stepan committed
353
354
355
356
357
358
    NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
    %% Correct initial conditions.
    if NumberOfSavedMhFiles
        load([MhDirectoryName '/' ModelName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']);
        ilogpo2(fblck) = logpo2(end);
        ix2(fblck,:) = x2(end,:);
359
    end
stepan's avatar
stepan committed
360
end