adaptive_metropolis_hastings

function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
 Random walk Metropolis-Hastings algorithm. 
 
 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  

 ALGORITHM 
   Metropolis-Hastings.       

 SPECIAL REQUIREMENTS
   None.

 PARALLEL CONTEXT
 The most computationally intensive part of this function may be executed
 in parallel. The code sutable to be executed in
 parallel on multi core or cluster machine (in general a 'for' cycle),
 is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
 Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in
 parallel and called name_function.m and name_function_core.m. 
 In addition in parallel package we have second set of functions used
 to manage the parallel computation.

 This function was the first function to be parallelized, later other
 functions have been parallelized using the same methodology.
 Then the comments write here can be used for all the other pairs of
 parallel functions and also for management funtions.

0001 function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
0002 %function adaptive_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
0003 % Random walk Metropolis-Hastings algorithm.
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 % ALGORITHM
0017 %   Metropolis-Hastings.
0018 %
0019 % SPECIAL REQUIREMENTS
0020 %   None.
0021 %
0022 % PARALLEL CONTEXT
0023 % The most computationally intensive part of this function may be executed
0024 % in parallel. The code sutable to be executed in
0025 % parallel on multi core or cluster machine (in general a 'for' cycle),
0026 % is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
0027 % Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in
0028 % parallel and called name_function.m and name_function_core.m.
0029 % In addition in parallel package we have second set of functions used
0030 % to manage the parallel computation.
0031 %
0032 % This function was the first function to be parallelized, later other
0033 % functions have been parallelized using the same methodology.
0034 % Then the comments write here can be used for all the other pairs of
0035 % parallel functions and also for management funtions.
0036 
0037 % Copyright (C) 2006-2011 Dynare Team
0038 %
0039 % This file is part of Dynare.
0040 %
0041 % Dynare is free software: you can redistribute it and/or modify
0042 % it under the terms of the GNU General Public License as published by
0043 % the Free Software Foundation, either version 3 of the License, or
0044 % (at your option) any later version.
0045 %
0046 % Dynare is distributed in the hope that it will be useful,
0047 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0048 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0049 % GNU General Public License for more details.
0050 %
0051 % You should have received a copy of the GNU General Public License
0052 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0053 
0054 global M_ options_ bayestopt_ estim_params_ oo_
0055 
0056 
0057 old_options = options_;
0058 
0059 accept_target = options_.amh.accept_target;
0060 m_directory = [M_.fname '/metropolis/']; 
0061 
0062 if options_.load_mh_file == 0
0063     delete([m_directory 'adaptive_metropolis_proposal_*.mat']);
0064     nP = 0;
0065 else
0066     D = dir([m_directory 'adaptive_metropolis_proposal_*.mat']);
0067     nP = size(D,1);
0068 end;
0069 
0070 if nP == 0
0071     jscale = options_.mh_jscale;
0072     bayestopt_.jscale = jscale;
0073     save([m_directory 'adaptive_metropolis_proposal_0'],'vv','jscale');
0074     nP = 1;
0075 else
0076     tmp = load([m_directory 'adaptive_metropolis_proposal_' ...
0077                 int2str(nP-1)],'vv','jscale');
0078     vv = tmp.vv;
0079     bayestopt_.jscale = tmp.jscale;
0080 end
0081 
0082 if options_.amh.cova_steps
0083     bayestopt_.jscale = tune_scale_parameter(TargetFun, ...
0084                                               ProposalFun,xparam1,vv,mh_bounds,varargin{:});
0085 end
0086 
0087 for i=1:options_.amh.cova_steps
0088     options_.mh_replic = options_.amh.cova_replic;
0089     random_walk_metropolis_hastings(TargetFun,ProposalFun, ...
0090                                     xparam1,vv,mh_bounds,varargin{:});
0091     tot_draws = total_draws(M_);
0092     options_.mh_drop = (tot_draws-options_.amh.cova_replic)/tot_draws;
0093     CutSample(M_,options_,estim_params_);
0094     [junk,vv] = compute_mh_covariance_matrix();
0095     jscale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin{:});
0096     bayestopt_.jscale = jscale;
0097     save([m_directory 'adaptive_metropolis_proposal_' ...
0098           int2str(nP)],'vv','jscale');
0099     nP = nP + 1;
0100 end
0101 
0102 options_.mh_replic = old_options.mh_replic;
0103 options_.mh_drop = old_options.mh_drop;
0104 record = random_walk_metropolis_hastings(TargetFun,ProposalFun, ...
0105                                          xparam1,vv,mh_bounds,varargin{:});
0106                 
0107 
0108 
0109 function selected_scale = tune_scale_parameter(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
0110 global options_ bayestopt_
0111 
0112 selected_scale = [];
0113 
0114 maxit = options_.amh.scale_tuning_maxit;
0115 accept_target = options_.amh.accept_target;
0116 test_runs = options_.amh.scale_tuning_test_runs;
0117 tolerance = options_.amh.scale_tuning_tolerance;
0118 Scales = zeros(maxit,1);
0119 AvRates = zeros(maxit,1);
0120 Scales(1) = bayestopt_.jscale;
0121 
0122 for i=1:maxit
0123     options_.mh_replic = options_.amh.scale_tuning_blocksize;
0124     bayestopt_.jscale = Scales(i);
0125     record = random_walk_metropolis_hastings(TargetFun,ProposalFun, ...
0126                                              xparam1,vv, ...
0127                                              mh_bounds,varargin{:});
0128     AvRates(i) = mean(record.AcceptationRates);    
0129 
0130     if i < test_runs
0131         i_kept_runs = 1:i;
0132     else
0133         i_kept_runs = i_kept_runs+1;
0134     end
0135     
0136     kept_runs_s = Scales(i_kept_runs);
0137     kept_runs_a = AvRates(i_kept_runs);
0138     
0139     if i > test_runs
0140         a_mean = mean(kept_runs_a);
0141         if (a_mean > (1-tolerance)*accept_target) && ...
0142                 (a_mean < (1+tolerance)*accept_target) && ...
0143                 all(kept_runs_a > (1-test_runs*tolerance)*accept_target) && ...
0144                 all(kept_runs_a < (1+test_runs*tolerance)*accept_target)
0145             selected_scale = mean(Scales((i-test_runs+1):i));
0146             disp(['Selected scale: ' num2str(selected_scale)])
0147             return
0148         end
0149     end
0150     
0151     mean_kept_runs_a = mean(kept_runs_a);
0152     if (mean_kept_runs_a/accept_target < 1-test_runs*tolerance) ...
0153             || (mean_kept_runs_a/accept_target > 1+test_runs*tolerance)
0154         damping_factor = 1
0155     else
0156         damping_factor = 1/3
0157     end
0158     Scales(i+1) = mean(kept_runs_s)*(mean(kept_runs_a)/ ...
0159                                      accept_target)^damping_factor;
0160 
0161 
0162     options_.load_mh_file = 1;
0163     
0164     disp(100*kept_runs_s')
0165     disp(100*kept_runs_a')
0166     disp(['Selected scale ' num2str(Scales(i+1))])    
0167 end
0168 
0169 error('AMH scale tuning: tuning didn''t converge')
0170 
0171 function y = total_draws(M_)
0172 load([M_.fname '/metropolis/' M_.fname '_mh_history'])
0173 y = sum(record.MhDraws(:,1));