function [PostMod,PostVar,Scale,PostMean] = ...
gmhmaxlik(ObjFun,xparam1,mh_bounds,num,iScale,info,MeanPar,VarCov,varargin)
(Dirty) Global minimization routine of (minus) a likelihood (or posterior density) function.
INPUTS
o ObjFun [char] string specifying the name of the objective function.
o xparam1 [double] (p*1) vector of parameters to be estimated.
o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
o num [integer] scalar specifying the number of MH iterations in step 2.
o iScale [double] scalar specifying the initial of the jumping distribution's scale parameter.
o info [char] string, empty or equal to 'LastCall'.
o MeanPar [double] (p*1) vector specifying the initial posterior mean.
o VarCov [double] (p*p) matrix specifying the initial posterior covariance matrix.
o gend [integer] scalar specifying the number of observations ==> varargin{1}.
o data [double] (T*n) matrix of data ==> varargin{2}.
OUTPUTS
o PostMod [double] (p*1) vector, evaluation of the posterior mode.
o PostVar [double] (p*p) matrix, evaluation of the posterior covariance matrix.
o Scale [double] scalar specifying the scale parameter that should be used in
an eventual metropolis-hastings algorithm.
o PostMean [double] (p*1) vector, evaluation of the posterior mean.
ALGORITHM
Metropolis-Hastings with an constantly updated covariance matrix for
the jump distribution. The posterior mean, variance and mode are
updated (in step 2) with the following rules:
\[
\mu_t = \mu_{t-1} + \frac{1}{t}\left(\theta_t-\mu_{t-1}\right)
\]
\[
\Sigma_t = \Sigma_{t-1} + \mu_{t-1}\mu_{t-1}'-\mu_{t}\mu_{t}' +
\frac{1}{t}\left(\theta_t\theta_t'-\Sigma_{t-1}-\mu_{t-1}\mu_{t-1}'\right)
\]
and
\[
\mathrm{mode}_t = \left\{
\begin{array}{ll}
\theta_t, & \hbox{if } p(\theta_t|\mathcal Y) > p(\mathrm{mode}_{t-1}|\mathcal Y) \\
\mathrm{mode}_{t-1}, & \hbox{otherwise.}
\end{array}
\right.
\]
where $t$ is the iteration, $\mu_t$ the estimate of the posterior mean
after $t$ iterations, $\Sigma_t$ the estimate of the posterior
covariance matrix after $t$ iterations, $\mathrm{mode}_t$ is the
evaluation of the posterior mode after $t$ iterations and
$p(\theta_t|\mathcal Y)$ is the posterior density of parameters
(specified by the user supplied function "fun").
SPECIAL REQUIREMENTS
None.0001 function [PostMod,PostVar,Scale,PostMean] = ... 0002 gmhmaxlik(ObjFun,xparam1,mh_bounds,num,iScale,info,MeanPar,VarCov,varargin) 0003 0004 %function [PostMod,PostVar,Scale,PostMean] = ... 0005 %gmhmaxlik(ObjFun,xparam1,mh_bounds,num,iScale,info,MeanPar,VarCov,varargin) 0006 % (Dirty) Global minimization routine of (minus) a likelihood (or posterior density) function. 0007 % 0008 % INPUTS 0009 % o ObjFun [char] string specifying the name of the objective function. 0010 % o xparam1 [double] (p*1) vector of parameters to be estimated. 0011 % o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. 0012 % o num [integer] scalar specifying the number of MH iterations in step 2. 0013 % o iScale [double] scalar specifying the initial of the jumping distribution's scale parameter. 0014 % o info [char] string, empty or equal to 'LastCall'. 0015 % o MeanPar [double] (p*1) vector specifying the initial posterior mean. 0016 % o VarCov [double] (p*p) matrix specifying the initial posterior covariance matrix. 0017 % o gend [integer] scalar specifying the number of observations ==> varargin{1}. 0018 % o data [double] (T*n) matrix of data ==> varargin{2}. 0019 % 0020 % OUTPUTS 0021 % o PostMod [double] (p*1) vector, evaluation of the posterior mode. 0022 % o PostVar [double] (p*p) matrix, evaluation of the posterior covariance matrix. 0023 % o Scale [double] scalar specifying the scale parameter that should be used in 0024 % an eventual metropolis-hastings algorithm. 0025 % o PostMean [double] (p*1) vector, evaluation of the posterior mean. 0026 % 0027 % ALGORITHM 0028 % Metropolis-Hastings with an constantly updated covariance matrix for 0029 % the jump distribution. The posterior mean, variance and mode are 0030 % updated (in step 2) with the following rules: 0031 % 0032 % \[ 0033 % \mu_t = \mu_{t-1} + \frac{1}{t}\left(\theta_t-\mu_{t-1}\right) 0034 % \] 0035 % 0036 % \[ 0037 % \Sigma_t = \Sigma_{t-1} + \mu_{t-1}\mu_{t-1}'-\mu_{t}\mu_{t}' + 0038 % \frac{1}{t}\left(\theta_t\theta_t'-\Sigma_{t-1}-\mu_{t-1}\mu_{t-1}'\right) 0039 % \] 0040 % 0041 % and 0042 % 0043 % \[ 0044 % \mathrm{mode}_t = \left\{ 0045 % \begin{array}{ll} 0046 % \theta_t, & \hbox{if } p(\theta_t|\mathcal Y) > p(\mathrm{mode}_{t-1}|\mathcal Y) \\ 0047 % \mathrm{mode}_{t-1}, & \hbox{otherwise.} 0048 % \end{array} 0049 % \right. 0050 % \] 0051 % 0052 % where $t$ is the iteration, $\mu_t$ the estimate of the posterior mean 0053 % after $t$ iterations, $\Sigma_t$ the estimate of the posterior 0054 % covariance matrix after $t$ iterations, $\mathrm{mode}_t$ is the 0055 % evaluation of the posterior mode after $t$ iterations and 0056 % $p(\theta_t|\mathcal Y)$ is the posterior density of parameters 0057 % (specified by the user supplied function "fun"). 0058 % 0059 % SPECIAL REQUIREMENTS 0060 % None. 0061 0062 % Copyright (C) 2006-2011 Dynare Team 0063 % 0064 % This file is part of Dynare. 0065 % 0066 % Dynare is free software: you can redistribute it and/or modify 0067 % it under the terms of the GNU General Public License as published by 0068 % the Free Software Foundation, either version 3 of the License, or 0069 % (at your option) any later version. 0070 % 0071 % Dynare is distributed in the hope that it will be useful, 0072 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0073 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0074 % GNU General Public License for more details. 0075 % 0076 % You should have received a copy of the GNU General Public License 0077 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0078 0079 global bayestopt_ estim_params_ options_ 0080 0081 options_.lik_algo = 1; 0082 npar = length(xparam1); 0083 0084 NumberOfIterations = num; 0085 MaxNumberOfTuningSimulations = 200000; 0086 MaxNumberOfClimbingSimulations = 200000; 0087 AcceptanceTarget = 1/3; 0088 0089 CovJump = VarCov; 0090 ModePar = xparam1; 0091 0092 %% [1] I tune the scale parameter. 0093 hh = dyn_waitbar(0,'Tuning of the scale parameter...'); 0094 set(hh,'Name','Tuning of the scale parameter.'); 0095 j = 1; jj = 1; 0096 isux = 0; jsux = 0; test = 0; 0097 ix2 = ModePar;% initial condition! 0098 ilogpo2 = - feval(ObjFun,ix2,varargin{:});% initial posterior density 0099 mlogpo2 = ilogpo2; 0100 try 0101 dd = transpose(chol(CovJump)); 0102 catch 0103 dd = eye(length(CovJump)); 0104 end 0105 while j<=MaxNumberOfTuningSimulations 0106 proposal = iScale*dd*randn(npar,1) + ix2; 0107 if all(proposal > mh_bounds(:,1)) && all(proposal < mh_bounds(:,2)) 0108 logpo2 = - feval(ObjFun,proposal,varargin{:}); 0109 else 0110 logpo2 = -inf; 0111 end 0112 % I move if the proposal is enough likely... 0113 if logpo2 > -inf && log(rand) < logpo2 - ilogpo2 0114 ix2 = proposal; 0115 if logpo2 > mlogpo2 0116 ModePar = proposal; 0117 mlogpo2 = logpo2; 0118 end 0119 ilogpo2 = logpo2; 0120 isux = isux + 1; 0121 jsux = jsux + 1; 0122 end% ... otherwise I don't move. 0123 prtfrc = j/MaxNumberOfTuningSimulations; 0124 if mod(j, 10)==0 0125 dyn_waitbar(prtfrc,hh,sprintf('Acceptance rates: %f [%f]',isux/j,jsux/jj)); 0126 end 0127 if j/500 == round(j/500) 0128 test1 = jsux/jj; 0129 cfactor = test1/AcceptanceTarget; 0130 if cfactor>0 0131 iScale = iScale*cfactor; 0132 else 0133 iScale = iScale/10; 0134 end 0135 jsux = 0; jj = 0; 0136 if cfactor>0.90 && cfactor<1.10 0137 test = test+1; 0138 end 0139 if test>4 0140 break 0141 end 0142 end 0143 j = j+1; 0144 jj = jj + 1; 0145 end 0146 0147 dyn_waitbar_close(hh); 0148 %% [2] One block metropolis, I update the covariance matrix of the jumping distribution 0149 hh = dyn_waitbar(0,'Metropolis-Hastings...'); 0150 set(hh,'Name','Estimation of the posterior covariance...'), 0151 j = 1; 0152 isux = 0; 0153 ilogpo2 = - feval(ObjFun,ix2,varargin{:}); 0154 while j<= NumberOfIterations 0155 proposal = iScale*dd*randn(npar,1) + ix2; 0156 if all(proposal > mh_bounds(:,1)) && all(proposal < mh_bounds(:,2)) 0157 logpo2 = - feval(ObjFun,proposal,varargin{:}); 0158 else 0159 logpo2 = -inf; 0160 end 0161 % I move if the proposal is enough likely... 0162 if logpo2 > -inf && log(rand) < logpo2 - ilogpo2 0163 ix2 = proposal; 0164 if logpo2 > mlogpo2 0165 ModePar = proposal; 0166 mlogpo2 = logpo2; 0167 end 0168 ilogpo2 = logpo2; 0169 isux = isux + 1; 0170 jsux = jsux + 1; 0171 end% ... otherwise I don't move. 0172 prtfrc = j/NumberOfIterations; 0173 if mod(j, 10)==0 0174 dyn_waitbar(prtfrc,hh,sprintf('Acceptance rate: %f',isux/j)); 0175 end 0176 % I update the covariance matrix and the mean: 0177 oldMeanPar = MeanPar; 0178 MeanPar = oldMeanPar + (1/j)*(ix2-oldMeanPar); 0179 CovJump = CovJump + oldMeanPar*oldMeanPar' - MeanPar*MeanPar' + ... 0180 (1/j)*(ix2*ix2' - CovJump - oldMeanPar*oldMeanPar'); 0181 j = j+1; 0182 end 0183 dyn_waitbar_close(hh); 0184 PostVar = CovJump; 0185 PostMean = MeanPar; 0186 %% [3 & 4] I tune the scale parameter (with the new covariance matrix) if 0187 %% this is the last call to the routine, and I climb the hill (without 0188 %% updating the covariance matrix)... 0189 if strcmpi(info,'LastCall') 0190 hh = dyn_waitbar(0,'Tuning of the scale parameter...'); 0191 set(hh,'Name','Tuning of the scale parameter.'), 0192 j = 1; jj = 1; 0193 isux = 0; jsux = 0; 0194 test = 0; 0195 ilogpo2 = - feval(ObjFun,ix2,varargin{:});% initial posterior density 0196 dd = transpose(chol(CovJump)); 0197 while j<=MaxNumberOfTuningSimulations 0198 proposal = iScale*dd*randn(npar,1) + ix2; 0199 if all(proposal > mh_bounds(:,1)) && all(proposal < mh_bounds(:,2)) 0200 logpo2 = - feval(ObjFun,proposal,varargin{:}); 0201 else 0202 logpo2 = -inf; 0203 end 0204 % I move if the proposal is enough likely... 0205 if logpo2 > -inf && log(rand) < logpo2 - ilogpo2 0206 ix2 = proposal; 0207 if logpo2 > mlogpo2 0208 ModePar = proposal; 0209 mlogpo2 = logpo2; 0210 end 0211 ilogpo2 = logpo2; 0212 isux = isux + 1; 0213 jsux = jsux + 1; 0214 end% ... otherwise I don't move. 0215 prtfrc = j/MaxNumberOfTuningSimulations; 0216 if mod(j, 10)==0 0217 dyn_waitbar(prtfrc,hh,sprintf('Acceptance rates: %f [%f]',isux/j,jsux/jj)); 0218 end 0219 if j/1000 == round(j/1000) 0220 test1 = jsux/jj; 0221 cfactor = test1/AcceptanceTarget; 0222 iScale = iScale*cfactor; 0223 jsux = 0; jj = 0; 0224 if cfactor>0.90 && cfactor<1.10 0225 test = test+1; 0226 end 0227 if test>4 0228 break 0229 end 0230 end 0231 j = j+1; 0232 jj = jj + 1; 0233 end 0234 dyn_waitbar_close(hh); 0235 Scale = iScale; 0236 %% 0237 %% Now I climb the hill 0238 %% 0239 climb = 1; 0240 if climb 0241 hh = dyn_waitbar(0,' '); 0242 set(hh,'Name','Now I am climbing the hill...'), 0243 disp('Now I am climbing the hill...'), 0244 j = 1; jj = 1; 0245 jsux = 0; 0246 test = 0; 0247 while j<=MaxNumberOfClimbingSimulations 0248 proposal = iScale*dd*randn(npar,1) + ModePar; 0249 if all(proposal > mh_bounds(:,1)) && all(proposal < mh_bounds(:,2)) 0250 logpo2 = - feval(ObjFun,proposal,varargin{:}); 0251 else 0252 logpo2 = -inf; 0253 end 0254 if logpo2 > mlogpo2% I move if the proposal is higher... 0255 ModePar = proposal; 0256 mlogpo2 = logpo2; 0257 jsux = jsux + 1; 0258 end% otherwise I don't move... 0259 prtfrc = j/MaxNumberOfClimbingSimulations; 0260 if mod(j, 10)==0 0261 dyn_waitbar(prtfrc,hh,sprintf('%f Jumps / MaxStepSize %f',jsux,sqrt(max(diag(iScale*CovJump))))); 0262 end 0263 if j/200 == round(j/200) 0264 if jsux<=1 0265 test = test+1; 0266 else 0267 test = 0; 0268 end 0269 jsux = 0; 0270 jj = 0; 0271 if test>4% If I do not progress enough I reduce the scale parameter 0272 % of the jumping distribution (cooling down the system). 0273 iScale = iScale/1.10; 0274 end 0275 if sqrt(max(diag(iScale*CovJump)))<10^(-4) 0276 break% Steps are too small! 0277 end 0278 end 0279 j = j+1; 0280 jj = jj + 1; 0281 end 0282 dyn_waitbar_close(hh); 0283 end%climb 0284 else 0285 Scale = iScale; 0286 end 0287 PostMod = ModePar;