diff --git a/doc/dynare.texi b/doc/dynare.texi index d7b88d109..c58e95b79 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4389,6 +4389,7 @@ the total number of Metropolis draws available. Default: The fraction of initially generated parameter vectors to be dropped as a burnin before using posterior simulations. Default: @code{0.5} @item mh_jscale = @var{DOUBLE} +@anchor{mh_jscale} The scale to be used for the jumping distribution in Metropolis-Hastings algorithm. The default value is rarely satisfactory. This option must be tuned to obtain, ideally, an @@ -4485,6 +4486,26 @@ value of that function as the posterior mode. @noindent Default value is @code{4}. +@item MCMC_jumping_covariance = hessian|prior_variance|identity_matrix|@var{filename} +Tells Dynare which covariance to use for the proposal density of the MCMC sampler. @code{MCMC_jumping_covariance} can be one of the following: + +@table @code +@item hessian +Uses the Hessian matrix computed at the mode. + +@item prior_variance +Uses the prior variances. No infinite prior variances are allowed in this case. + +@item identity_matrix +Uses an identity matrix. + +@item filename +Loads an arbitrary user-specified covariance matrix from @code{filename.mat}. The covariance matrix must be saved in a variable named @code{jumping_covariance}, must be square, positive definite, and have the same dimension as the number of estimated parameters. + +@end table +@noindent +Note that the covariance matrices are stil scaled with @ref{mh_jscale}. Default value is @code{hessian}. + @item mode_check Tells Dynare to plot the posterior density for values around the computed mode for each estimated parameter in turn. This is helpful to diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 0a553921a..67de8c602 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -609,8 +609,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!']) end end - if ~isequal(options_.mode_compute,6) - if options_.cova_compute == 1 + if ~isequal(options_.mode_compute,6) %always already computes covariance matrix + if options_.cova_compute == 1 %user did not request covariance not to be computed if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'), ana_deriv = options_.analytic_derivation; options_.analytic_derivation = 2; @@ -636,6 +636,35 @@ if options_.cova_compute == 0 hh = [];%NaN(length(xparam1),length(xparam1)); end +switch options_.MCMC_jumping_covariance + case 'hessian' %Baseline + %do nothing and use hessian from mode_compute + case 'prior_variance' %Use prior variance + if any(isinf(bayestopt_.p2)) + error('Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.') + else + hh = diag(1./(bayestopt_.p2.^2)); + end + case 'identity_matrix' %Use identity + hh = eye(nx); + otherwise %user specified matrix in file + try + load(options_.MCMC_jumping_covariance,'jumping_covariance') + hh=jumping_covariance; + catch + error(['No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat']) + end + [nrow, ncol]=size(hh); + if ~isequal(nrow,ncol) && ~isequal(nrow,nx) %check if square and right size + error(['jumping_covariance matrix must be square and have ',num2str(nx),' rows and columns']) + end + try %check for positive definiteness + chol(hh); + catch + error(['Specified jumping_covariance is not positive definite']) + end +end + if ~options_.mh_posterior_mode_estimation && options_.cova_compute try chol(hh); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 36c86ef28..a663dba26 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -381,6 +381,7 @@ options_.mh_nblck = 2; options_.mh_recover = 0; options_.mh_replic = 20000; options_.recursive_estimation_restart = 0; +options_.MCMC_jumping_covariance='hessian'; options_.mode_compute = 4; options_.mode_file = ''; diff --git a/tests/estimation/fs2000_MCMC_jumping_covariance.mod b/tests/estimation/fs2000_MCMC_jumping_covariance.mod new file mode 100644 index 000000000..e7ff61da3 --- /dev/null +++ b/tests/estimation/fs2000_MCMC_jumping_covariance.mod @@ -0,0 +1,86 @@ +// See fs2000.mod in the examples/ directory for details on the model + +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +initval; +k = 6; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +gam, normal_pdf, 0.0085, 0.003; +mst, normal_pdf, 1.0002, 0.007; +rho, beta_pdf, 0.129, 0.223; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, .1; +stderr e_m, inv_gamma_pdf, 0.008862, .1; +end; + +varobs gp_obs gy_obs; + +options_.solve_tolf = 1e-12; +options_.mode_compute=4; +options_.plot_priors=0; +options_.MCMC_jumping_covariance='hessian'; +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=1000,mh_nblocks=1,mh_jscale=0.8); +load fs2000_MCMC_jumping_covariance_mode hh; +jumping_covariance=diag(diag(hh)); +save test_matrix jumping_covariance; +options_.MCMC_jumping_covariance='prior_variance'; +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=1000,mh_nblocks=1,mh_jscale=0.01); +options_.MCMC_jumping_covariance='identity_matrix'; +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=1000,mh_nblocks=1,mh_jscale=0.0001); +options_.MCMC_jumping_covariance='test_matrix'; +estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=1000,mh_nblocks=1,mh_jscale=0.8); diff --git a/tests/estimation/fs2000_MCMC_jumping_covariance_steadystate.m b/tests/estimation/fs2000_MCMC_jumping_covariance_steadystate.m new file mode 100644 index 000000000..2ac140da5 --- /dev/null +++ b/tests/estimation/fs2000_MCMC_jumping_covariance_steadystate.m @@ -0,0 +1,73 @@ +% computes the steady state of fs2000 analyticaly +% largely inspired by the program of F. Schorfheide + +% Copyright (C) 2004-2010 Dynare Team +% +% 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 . + +function [ys,check] = fs2000_steadystate(ys,exe) + global M_ + + alp = M_.params(1); + bet = M_.params(2); + gam = M_.params(3); + mst = M_.params(4); + rho = M_.params(5); + psi = M_.params(6); + del = M_.params(7); + + check = 0; + + dA = exp(gam); + gst = 1/dA; + m = mst; + + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; + + ys =[ +m +P +c +e +W +R +k +d +n +l +gy_obs +gp_obs +y +dA ];