metropolis_hastings_initialization

 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.

0001 function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
0002     metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin)
0003 % Metropolis-Hastings initialization.
0004 %
0005 % INPUTS
0006 %   o TargetFun  [char]     string specifying the name of the objective
0007 %                           function (posterior kernel).
0008 %   o xparam1    [double]   (p*1) vector of parameters to be estimated (initial values).
0009 %   o vv         [double]   (p*p) matrix, posterior covariance matrix (at the mode).
0010 %   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters.
0011 %   o varargin              list of argument following mh_bounds
0012 %
0013 % OUTPUTS
0014 %   None
0015 %
0016 % SPECIAL REQUIREMENTS
0017 %   None.
0018 
0019 % Copyright (C) 2006-2011 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 global M_ options_ bayestopt_
0037 
0038 ix2 = [];
0039 ilogpo2 = [];
0040 ModelName = []; 
0041 MhDirectoryName = [];
0042 fblck = [];
0043 fline = [];
0044 npar = [];
0045 nblck = [];
0046 nruns = [];
0047 NewFile = [];
0048 MAX_nruns = [];
0049 d = [];
0050 
0051 ModelName = M_.fname;
0052 
0053 if ~isempty(M_.bvar)
0054     ModelName = [M_.fname '_bvar'];
0055 end
0056 
0057 bayestopt_.penalty = 1e8;
0058 
0059 MhDirectoryName = CheckPath('metropolis',M_.dname);
0060 
0061 nblck = options_.mh_nblck;
0062 nruns = ones(nblck,1)*options_.mh_replic;
0063 npar  = length(xparam1);
0064 MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
0065 d = chol(vv);
0066 
0067 if ~options_.load_mh_file && ~options_.mh_recover
0068     %% Here we start a new metropolis-hastings, previous draws are not
0069     %% considered.
0070     if nblck > 1
0071         disp('MH: Multiple chains mode.')
0072     else
0073         disp('MH: One Chain mode.')
0074     end
0075     %% Delete old mh files...
0076     files = dir([ MhDirectoryName '/' ModelName '_mh*_blck*.mat']);
0077     if length(files)
0078         delete([ MhDirectoryName '/' ModelName '_mh*_blck*.mat']);
0079         disp('MH: Old _mh files successfully erased!')
0080     end
0081     file = dir([ MhDirectoryName '/metropolis.log']);
0082     if length(file)
0083         delete([ MhDirectoryName '/metropolis.log']);
0084         disp('MH: Old metropolis.log file successfully erased!')
0085         disp('MH: Creation of a new metropolis.log file.')
0086     end
0087     fidlog = fopen([MhDirectoryName '/metropolis.log'],'w');
0088     fprintf(fidlog,'%% MH log file (Dynare).\n');
0089     fprintf(fidlog,['%% ' datestr(now,0) '.\n']);
0090     fprintf(fidlog,' \n\n');
0091     fprintf(fidlog,'%% Session 1.\n');
0092     fprintf(fidlog,' \n');
0093     fprintf(fidlog,['  Number of blocks...............: ' int2str(nblck) '\n']);
0094     fprintf(fidlog,['  Number of simulations per block: ' int2str(nruns(1)) '\n']);
0095     fprintf(fidlog,' \n');
0096     %% Initial values...
0097     if nblck > 1% Case 1: multiple chains
0098         fprintf(fidlog,['  Initial values of the parameters:\n']);
0099         disp('MH: Searching for initial values...')
0100         ix2 = zeros(nblck,npar);
0101         ilogpo2 = zeros(nblck,1);
0102         for j=1:nblck
0103             validate    = 0;
0104             init_iter   = 0;
0105             trial = 1;
0106             while validate == 0 && trial <= 10
0107                 candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
0108                 if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2)) 
0109                     ix2(j,:) = candidate;
0110                     ilogpo2(j) = - feval(TargetFun,ix2(j,:)',varargin{:});
0111                     if ilogpo2(j) <= - bayestopt_.penalty+1e-6
0112                         validate = 0;
0113                     else
0114                         fprintf(fidlog,['    Blck ' int2str(j) ':\n']);
0115                         for i=1:length(ix2(1,:))
0116                             fprintf(fidlog,['      params:' int2str(i) ': ' num2str(ix2(j,i)) '\n']);
0117                         end
0118                         fprintf(fidlog,['      logpo2: ' num2str(ilogpo2(j)) '\n']);
0119                         j = j+1;
0120                         validate = 1;
0121                     end
0122                 end
0123                 init_iter = init_iter + 1;
0124                 if init_iter > 100 && validate == 0
0125                     disp(['MH: I couldn''t get a valid initial value in 100 trials.'])
0126                     disp(['MH: You should Reduce mh_init_scale...'])
0127                     disp(sprintf('MH: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale))
0128                     options_.mh_init_scale = input('MH: Enter a new value...  ');
0129                     trial = trial+1;
0130                 end
0131             end
0132             if trial > 10 && ~validate
0133                 disp(['MH: I''m unable to find a starting value for block ' int2str(j)])
0134                 return
0135             end
0136         end
0137         fprintf(fidlog,' \n');
0138         disp('MH: Initial values found!')
0139         disp(' ')
0140     else% Case 2: one chain (we start from the posterior mode)
0141         fprintf(fidlog,['  Initial values of the parameters:\n']);
0142         candidate = transpose(xparam1);
0143         if all(candidate' > mh_bounds(:,1)) && all(candidate' < mh_bounds(:,2)) 
0144             ix2 = candidate;
0145             ilogpo2 = - feval(TargetFun,ix2',varargin{:});
0146             disp('MH: Initialization at the posterior mode.')
0147             disp(' ')
0148             fprintf(fidlog,['    Blck ' int2str(1) 'params:\n']);
0149             for i=1:length(ix2(1,:))
0150                 fprintf(fidlog,['      ' int2str(i)  ':' num2str(ix2(1,i)) '\n']);
0151             end
0152             fprintf(fidlog,['    Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']);
0153         else
0154             disp('MH: Initialization failed...')
0155             disp('MH: The posterior mode lies outside the prior bounds.')
0156             return
0157         end
0158         fprintf(fidlog,' \n');
0159     end
0160     fprintf(fidlog,' \n');
0161     fblck = 1;
0162     fline = ones(nblck,1);
0163     NewFile = ones(nblck,1);
0164     %%
0165     %% Creation of the mh-history file:
0166     %%
0167     file = dir([MhDirectoryName '/'  ModelName '_mh_history.mat']);
0168     if length(files)
0169         delete([ MhDirectoryName '/' ModelName '_mh_history.mat']);
0170         disp('MH: Old mh_history file successfully erased!')
0171     end
0172     AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
0173     AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
0174     record.Nblck = nblck;
0175     record.MhDraws = zeros(1,3);
0176     record.MhDraws(1,1) = nruns(1);
0177     record.MhDraws(1,2) = AnticipatedNumberOfFiles;
0178     record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
0179     record.AcceptationRates = zeros(1,nblck);
0180     % separate initializaton for each chain
0181     JSUM = 0;
0182     for j=1:nblck,
0183         JSUM  = JSUM + sum(100*clock);
0184         randn('state',JSUM);
0185         rand('state',JSUM);
0186         record.Seeds(j).Normal = randn('state');
0187         record.Seeds(j).Unifor = rand('state');
0188     end
0189     record.InitialParameters = ix2;
0190     record.InitialLogLiK = ilogpo2;
0191     record.LastParameters = zeros(nblck,npar);
0192     record.LastLogLiK = zeros(nblck,1);
0193     record.LastFileNumber = AnticipatedNumberOfFiles ;
0194     record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
0195     save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');  
0196     fprintf(fidlog,['  CREATION OF THE MH HISTORY FILE!\n\n']);
0197     fprintf(fidlog,['    Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']);
0198     fprintf(fidlog,['    Expected number of lines in the last file: ' ...
0199                     int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
0200     fprintf(fidlog,['\n']);
0201     for j = 1:nblck,
0202         fprintf(fidlog,['    Initial seed (randn) for chain number ',int2str(j),':\n']);
0203         for i=1:length(record.Seeds(j).Normal)
0204             fprintf(fidlog,['      ' num2str(record.Seeds(j).Normal(i)') '\n']);
0205         end
0206         fprintf(fidlog,['    Initial seed (rand) for chain number ',int2str(j),':\n']);
0207         for i=1:length(record.Seeds(j).Unifor)
0208             fprintf(fidlog,['      ' num2str(record.Seeds(j).Unifor(i)') '\n']);
0209         end
0210     end,
0211     fprintf(fidlog,' \n');
0212     fclose(fidlog);
0213 elseif options_.load_mh_file && ~options_.mh_recover
0214     %% Here we consider previous mh files (previous mh did not crash).
0215     disp('MH: I''m loading past metropolis-hastings simulations...')
0216     file = dir([ MhDirectoryName '/'  ModelName '_mh_history.mat' ]);
0217     files = dir([ MhDirectoryName filesep ModelName '_mh*.mat']);
0218     if ~length(files)
0219         disp('MH:: FAILURE! there is no MH file to load here!')
0220         return
0221     end
0222     if ~length(file)
0223         disp('MH:: FAILURE! there is no MH-history file!')
0224         return
0225     else
0226         load([ MhDirectoryName '/'  ModelName '_mh_history.mat'])
0227     end
0228     fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
0229     fprintf(fidlog,'\n');
0230     fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']);
0231     fprintf(fidlog,' \n');
0232     fprintf(fidlog,['  Number of blocks...............: ' int2str(nblck) '\n']);
0233     fprintf(fidlog,['  Number of simulations per block: ' int2str(nruns(1)) '\n']);
0234     fprintf(fidlog,' \n');
0235     past_number_of_blocks = record.Nblck;
0236     if past_number_of_blocks ~= nblck
0237         disp('MH:: The specified number of blocks doesn''t match with the previous number of blocks!')
0238         disp(['MH:: You declared ' int2str(nblck) ' blocks, but the previous number of blocks was ' ...
0239               int2str(past_number_of_blocks) '.'])
0240         disp(['MH:: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
0241         nblck = past_number_of_blocks;
0242         options_.mh_nblck = nblck;
0243     end
0244     % I read the last line of the last mh-file for initialization
0245     % of the new metropolis-hastings simulations:
0246     LastFileNumber = record.LastFileNumber;
0247     LastLineNumber = record.LastLineNumber;
0248     if LastLineNumber < MAX_nruns
0249         NewFile = ones(nblck,1)*LastFileNumber;
0250     else
0251         NewFile = ones(nblck,1)*(LastFileNumber+1);
0252     end
0253     ilogpo2 = record.LastLogLiK;
0254     ix2 = record.LastParameters;
0255     fblck = 1;
0256     fline = ones(nblck,1)*(LastLineNumber+1);
0257     NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
0258     record.MhDraws = [record.MhDraws;zeros(1,3)];
0259     NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
0260     NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile;
0261     AnticipatedNumberOfFiles = ceil(NumberOfDrawsToBeSaved/MAX_nruns);
0262     AnticipatedNumberOfLinesInTheLastFile = NumberOfDrawsToBeSaved - (AnticipatedNumberOfFiles-1)*MAX_nruns;  
0263     record.LastFileNumber = LastFileNumber + AnticipatedNumberOfFiles;
0264     record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
0265     record.MhDraws(end,1) = nruns(1);
0266     record.MhDraws(end,2) = AnticipatedNumberOfFiles;
0267     record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
0268     save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
0269     disp(['MH: ... It''s done. I''ve loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
0270     disp(' ')
0271     fclose(fidlog);
0272 elseif options_.mh_recover
0273     %% The previous metropolis-hastings crashed before the end! I try to
0274     %% recover the saved draws...
0275     disp('MH: Recover mode!')
0276     disp(' ')
0277     file = dir([MhDirectoryName '/'  ModelName '_mh_history.mat']);
0278     if ~length(file)
0279         disp('MH:: FAILURE! there is no MH-history file!')
0280         return
0281     else
0282         load([ MhDirectoryName '/'  ModelName '_mh_history.mat'])
0283     end
0284     nblck = record.Nblck;% Number of "parallel" mcmc chains.
0285     options_.mh_nblck = nblck;
0286     if size(record.MhDraws,1) == 1
0287         OldMh = 0;% The crashed metropolis was the first session.
0288     else
0289         OldMh = 1;% The crashed metropolis wasn't the first session.
0290     end
0291     %% Default initialization:
0292     if OldMh
0293         ilogpo2 = record.LastLogLiK;
0294         ix2 = record.LastParameters;
0295     else
0296         ilogpo2 = record.InitialLogLiK;
0297         ix2 = record.InitialParameters;
0298     end
0299     %% Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers.
0300     if OldMh
0301         LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session.
0302         LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions.
0303         if LastLineNumberInThePreviousMh < MAX_nruns% Test if the last mh files of the previous session are not complete (yes)
0304             NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh;
0305             fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1);
0306         else% The last mh files of the previous session are complete.
0307             NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1);
0308             fline = ones(nblck,1);
0309         end
0310     else
0311         LastLineNumberInThePreviousMh = 0;
0312         LastFileNumberInThePreviousMh = 0;
0313         NewFile = ones(nblck,1);
0314         fline = ones(nblck,1);
0315     end
0316     %% Set fblck (First block), an integer targeting the crashed mcmc chain.
0317     fblck = 1;
0318     %% How many mh files should we have ?
0319     ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
0320     ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
0321     %% I count the total number of saved mh files...
0322     AllMhFiles = dir([MhDirectoryName '/' ModelName '_mh*_blck*.mat']);
0323     TotalNumberOfMhFiles = length(AllMhFiles);
0324     %% And I quit if I can't find a crashed mcmc chain.
0325     if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
0326         disp('MH: It appears that you don''t need to use the mh_recover option!')
0327         disp('    You have to edit the mod file and remove the mh_recover') 
0328         disp('    in the estimation command')
0329         error
0330     end
0331     %% I count the number of saved mh files per block.
0332     NumberOfMhFilesPerBlock = zeros(nblck,1);
0333     for b = 1:nblck
0334         BlckMhFiles = dir([ MhDirectoryName '/' ModelName '_mh*_blck' int2str(b) '.mat']);
0335         NumberOfMhFilesPerBlock(b) = length(BlckMhFiles);
0336     end
0337     %% Is there a chain with less mh files than expected ?
0338     while fblck <= nblck
0339         if  NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock
0340             disp(['MH: Chain ' int2str(fblck) ' is not complete!'])
0341             break
0342             % The mh_recover session will start from chain fblck.
0343         else
0344             disp(['MH: Chain ' int2str(fblck) ' is complete!'])
0345         end
0346         fblck = fblck+1;
0347     end
0348     %% How many mh-files are saved in this block?
0349     NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck);
0350     %% Correct the number of saved mh files if the crashed metropolis was not the first session (so
0351     %% that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
0352     %% of the current session).
0353     if OldMh
0354         NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
0355     end
0356     NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
0357     %% Correct initial conditions.
0358     if NumberOfSavedMhFiles
0359         load([MhDirectoryName '/' ModelName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']);
0360         ilogpo2(fblck) = logpo2(end);
0361         ix2(fblck,:) = x2(end,:);
0362     end
0363 end