Commit 8a6e2384 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Merge remote-tracking branch 'jpfeifer/mode_compute'

parents 27f525e1 0d669a73
...@@ -4389,6 +4389,7 @@ the total number of Metropolis draws available. Default: ...@@ -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} 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} @item mh_jscale = @var{DOUBLE}
@anchor{mh_jscale}
The scale to be used for the jumping distribution in The scale to be used for the jumping distribution in
Metropolis-Hastings algorithm. The default value is rarely Metropolis-Hastings algorithm. The default value is rarely
satisfactory. This option must be tuned to obtain, ideally, an satisfactory. This option must be tuned to obtain, ideally, an
...@@ -4485,6 +4486,26 @@ value of that function as the posterior mode. ...@@ -4485,6 +4486,26 @@ value of that function as the posterior mode.
@noindent @noindent
Default value is @code{4}. 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 @item mode_check
Tells Dynare to plot the posterior density for values around the Tells Dynare to plot the posterior density for values around the
computed mode for each estimated parameter in turn. This is helpful to computed mode for each estimated parameter in turn. This is helpful to
......
...@@ -609,8 +609,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation ...@@ -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!']) error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!'])
end end
end end
if ~isequal(options_.mode_compute,6) if ~isequal(options_.mode_compute,6) %always already computes covariance matrix
if options_.cova_compute == 1 if options_.cova_compute == 1 %user did not request covariance not to be computed
if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'), if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'),
ana_deriv = options_.analytic_derivation; ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 2; options_.analytic_derivation = 2;
...@@ -636,6 +636,35 @@ if options_.cova_compute == 0 ...@@ -636,6 +636,35 @@ if options_.cova_compute == 0
hh = [];%NaN(length(xparam1),length(xparam1)); hh = [];%NaN(length(xparam1),length(xparam1));
end 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 if ~options_.mh_posterior_mode_estimation && options_.cova_compute
try try
chol(hh); chol(hh);
......
...@@ -381,6 +381,7 @@ options_.mh_nblck = 2; ...@@ -381,6 +381,7 @@ options_.mh_nblck = 2;
options_.mh_recover = 0; options_.mh_recover = 0;
options_.mh_replic = 20000; options_.mh_replic = 20000;
options_.recursive_estimation_restart = 0; options_.recursive_estimation_restart = 0;
options_.MCMC_jumping_covariance='hessian';
options_.mode_compute = 4; options_.mode_compute = 4;
options_.mode_file = ''; options_.mode_file = '';
......
// 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);
% 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 <http://www.gnu.org/licenses/>.
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 ];
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment