Efficiency & Cosmetic changes related to the design of bayestopt_ and

the computation of the prior density.

bayestopt_.p1 is always the prior mean
bayestopt_.p2 is always the prior standard deviation
bayestopt_.p3 is unchanged
bayestopt_.p4 is unchanged
bayestopt_.p5 [new field] is the prior mode
bayestopt_.p6 [new field] is the first hyper-parameter of the prior density
bayestopt_.p7 [new field] is the second hyper-parameter of the prior density
 
These fields are defined in  set_prior and are never changed after. In
the previous version of Dynare,  the hyper parameters of the densities
were  updated at  each iteration  of the  optimization routine  or the
metropolis.

Removed fields pmean and pstdev.

Vectorized the code in priordens.

Fixed the bug mentionned by Gianni. If a (logged) density is evaluated
outside the  prior domain, the  output of priordens if  minus infinity
(instead of a complex number).


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2556 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
stepan 2009-04-06 14:38:37 +00:00
parent f02ea822ec
commit aa31417a05
23 changed files with 644 additions and 514 deletions

View File

@ -74,7 +74,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
end
offset = offset+estim_params_.nvn;
end
end
if estim_params_.ncx
for i=1:estim_params_.ncx
k1 =estim_params_.corrx(i,1);
@ -330,6 +330,6 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior);
options_.kalman_algo = kalman_algo;

View File

@ -322,7 +322,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,g
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior);
options_.kalman_algo = kalman_algo;
llik=[-lnprior; .5*lik(start:end)];

View File

@ -217,7 +217,7 @@ else
lik = .5*lik;% Minus likelihood
end
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (lik-lnprior);
if (nargout == 6)

View File

@ -15,7 +15,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2006-2008 Dynare Team
% Copyright (C) 2006-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -96,14 +96,14 @@ if np
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
end
disp(sprintf(pformat,name,bayestopt_.pmean(ip),...
disp(sprintf(pformat,name,bayestopt_.p1(ip),...
post_mean, ...
hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
bayestopt_.p2(ip)));
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
@ -144,11 +144,11 @@ if nvx
M_.Sigma_e(k,k) = post_mean*post_mean;
end
end
disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval,...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
disp(sprintf(pformat,name,bayestopt_.p1(ip),post_mean,hpd_interval,...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
@ -184,11 +184,11 @@ if nvn
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
end
disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
disp(sprintf(pformat,name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
@ -237,11 +237,11 @@ if ncx
M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
end
end
disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
disp(sprintf(pformat, name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end
@ -288,11 +288,11 @@ if ncn
post_median,post_var,post_deciles,density);
end
end
disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
disp(sprintf(pformat, name,bayestopt_.p1(ip),post_mean,hpd_interval, ...
pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.p2(ip)));
if TeX
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.p1(ip),...
bayestopt_.p2(ip),post_mean,sqrt(post_var),hpd_interval);
end
ip = ip+1;
end

View File

@ -0,0 +1,87 @@
function m = compute_prior_mode(hyperparameters,shape)
% This function computes the mode of the prior distribution given the (two, three or four) hyperparameters
% of the prior distribution.
%
% INPUTS
% hyperparameters [double] 1*n vector of hyper parameters.
% shape [integer] scalar specifying the prior shape:
% shape=1 => Beta distribution,
% shape=2 => Gamma distribution,
% shape=3 => Gaussian distribution,
% shape=4 => Inverse Gamma (type 1) distribution,
% shape=5 => Uniform distribution,
% shape=6 => Inverse Gamma (type 2) distribution.
%
% OUTPUTS
% m [double] scalar or 2*1 vector, the prior mode.
%
% REMARKS
% [1] The size of the vector of hyperparameters is 3 when the Gamma or Inverse Gamma is shifted and 4 when
% the support of the Beta distribution is not [0,1].
% [2] The hyperparameters of the uniform distribution are the lower and upper bounds.
% [3] The uniform distribution has an infinity of modes. In this case the function returns the prior mean.
% [4] For the beta distribution we can have 1, 2 or an infinity of modes.
% Copyright (C) 2009 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/>.
m = NaN ;
switch shape
case 1
if (hyperparameters(1)>1 && hyperparameters(2)>1)
m = (hyperparameters(1)-1)/(hyperparameters(1)+hyperparameters(2)-2) ;
elseif (hyperparameters(1)<1 && hyperparameters(2)<1)
m = [ 0 ; 1 ] ;
elseif ( hyperparameters(1)<1 && hyperparameters(2)>1-eps ) || ( abs(hyperparameters(1)-1)<2*eps && hyperparameters(2)>1 )
m = 0;
elseif ( hyperparameters(1)>1 && hyperparameters(2)<1+eps ) || ( abs(hyperparameters(1)-1)<2*eps && hyperparameters(2)<1 )
m = 1;
elseif ( abs(hyperparameters(1)-1)<2*eps && abs(hyperparameters(2)-1)<2*eps )% Uniform distribution!
m = .5 ;
end
if length(hyperparameters)==4
m = m*(hyperparameters(4)-hyperparameters(3)) + hyperparameters(3) ;
end
case 2
if hyperparameters(1)<1
m = 0;
else
m = (hyperparameters(1)-1)*hyperparameters(2);
end
if length(hyperparameters)>2
m = m + hyperparameters(3);
end
case 3
m = hyperparameters(1);
case 4
% s = hyperparameters(1)
% nu = hyperparameters(2)
m = sqrt((hyperparameters(2)-1)/hyperparameters(1));
if length(hyperparameters)>2
m = m + hyperparameters(3);
end
case 5
m = .5*(hyperparameters(2)-hyperparameters(1)) ;
case 6
% s = hyperparameters(1)
% nu = hyperparameters(2)
m = hyperparameters(1)/(hyperparameters(2)+2) ;
if length(hyperparameters)>2
m = m + hyperparameters(3) ;
end
otherwise
error('Unknown prior shape!')
end

View File

@ -1,6 +1,5 @@
function [x,f,abscissa,dens,binf,bsup] = draw_prior_density(indx,bayestopt_);
% function [x,f,abscissa,dens,binf,bsup] = draw_prior_density(indx)
% Computes values of prior density at many points (before plotting)
% Computes values of prior densities at many points (before plotting)
%
% INPUTS
% indx [integer] Parameter number.
@ -13,9 +12,7 @@ function [x,f,abscissa,dens,binf,bsup] = draw_prior_density(indx,bayestopt_);
% dens [double] Row vector, density
% binf: [double] Scalar, first element of x
% bsup: [double] Scalar, last element of x
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2004-2009 Dynare Team
%
@ -34,73 +31,58 @@ function [x,f,abscissa,dens,binf,bsup] = draw_prior_density(indx,bayestopt_);
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
pmean = bayestopt_.pmean;
pshape = bayestopt_.pshape;
p1 = bayestopt_.p1;
p2 = bayestopt_.p2;
pshape = bayestopt_.pshape;
p3 = bayestopt_.p3;
p4 = bayestopt_.p4;
p6 = bayestopt_.p6;
p7 = bayestopt_.p7;
truncprior = 1e-3;
steps = 200;
indx
switch pshape(indx)
case 1 % Beta prior
case 1% Beta prior
density = @(x,a,b,aa,bb) betapdf((x-aa)/(bb-aa), a, b)/(bb-aa);
mu = (p1(indx)-p3(indx))/(p4(indx)-p3(indx));
stdd = p2(indx)/(p4(indx)-p3(indx));
a = (1-mu)*mu^2/stdd^2 - mu;
b = a*(1/mu-1);
aa = p3(indx);
bb = p4(indx);
infbound = betainv(truncprior,a,b)*(bb-aa)+aa;
supbound = betainv(1-truncprior,a,b)*(bb-aa)+aa;
infbound = betainv(truncprior,p6(indx),p7(indx))*(p4(indx)-p3(indx))+p3(indx);
supbound = betainv(1-truncprior,p6(indx),p7(indx))*(p4(indx)-p3(indx))+p3(indx);
stepsize = (supbound-infbound)/steps;
abscissa = infbound:stepsize:supbound;
dens = density(abscissa,a,b,aa,bb);
case 2 % Generalized Gamma prior
dens = density(abscissa,p6(indx),p7(indx),p3(indx),p4(indx));
case 2% Generalized Gamma prior
density = @(x,a,b,c) gampdf(x-c,a,b);
mu = p1(indx)-p3(indx);
b = p2(indx)^2/mu;
a = mu/b;
c = p3(indx);
infbound = gaminv(truncprior,a,b)+c;
supbound = gaminv(1-truncprior,a,b)+c;
infbound = gaminv(truncprior,p6(indx),p7(indx))+p3(indx);
supbound = gaminv(1-truncprior,p6(indx),p7(indx))+p3(indx);
stepsize = (supbound-infbound)/steps;
abscissa = infbound:stepsize:supbound;
dens = density(abscissa,a,b,c);
case 3 % Gaussian prior
a = p1(indx);
b = p2(indx);
infbound = norminv(truncprior,a,b);
supbound = norminv(1-truncprior,a,b);
dens = density(abscissa,p6(indx),p7(indx),p3(indx));
case 3% Gaussian prior
infbound = norminv(truncprior,p6(indx),p7(indx));
supbound = norminv(1-truncprior,p6(indx),p7(indx));
stepsize = (supbound-infbound)/steps;
abscissa = infbound:stepsize:supbound;
dens = normpdf(abscissa,a,b);
case 4 % Inverse-gamma of type 1 prior
nu = p2(indx);
s = p1(indx);
infbound = 1/sqrt(gaminv(1-10*truncprior, nu/2, 2/s));
supbound = 1/sqrt(gaminv(10*truncprior, nu/2, 2/s));
dens = normpdf(abscissa,p6(indx),p7(indx));
case 4% Inverse-gamma of type 1 prior
infbound = 1/sqrt(gaminv(1-10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
supbound = 1/sqrt(gaminv(10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
stepsize = (supbound-infbound)/steps;
abscissa = infbound:stepsize:supbound;
dens = exp(lpdfig1(abscissa,s,nu));
case 5 % Uniform prior
infbound = p1(indx);
supbound = p2(indx);
dens = exp(lpdfig1(abscissa-p3(indx),p6(indx),p7(indx)));
case 5% Uniform prior
infbound = p6(indx);
supbound = p7(indx);
stepsize = (supbound-infbound)/steps;
abscissa = infbound:stepsize:supbound;
dens = ones(1, steps) / (supbound-infbound);
case 6 % Inverse-gamma of type 2 prior
nu = p2(indx);
s = p1(indx);
infbound = 1/(gaminv(1-10*truncprior, nu/2, 2/s));
supbound = 1/(gaminv(10*truncprior, nu/2, 2/s));
stepsize = (supbound-infbound)/steps;
case 6% Inverse-gamma of type 2 prior
infbound = 1/(gaminv(1-10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
supbound = 1/(gaminv(10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
stepsize = (supbound-infbound)/steps ;
abscissa = infbound:stepsize:supbound;
dens = exp(lpdfig2(abscissa,s,nu));
otherwise
error(sprintf('draw_prior_density: unknown distribution shape (index %d, type %d)', indx, pshape(indx)));
dens = exp(lpdfig2(abscissa-p3(indx),p6(indx),p7(indx)));
otherwise
error(sprintf('draw_prior_density: unknown distribution shape (index %d, type %d)', indx, pshape(indx)));
end
k = [1:length(dens)];

View File

@ -314,6 +314,6 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior);
options_.kalman_algo = kalman_algo;

View File

@ -118,6 +118,9 @@ else% If estim_params_ is empty...
bayestopt_.p2 = [];
bayestopt_.p3 = [];
bayestopt_.p4 = [];
bayestopt_.p5 = [];
bayestopt_.p6 = [];
bayestopt_.p7 = [];
estim_params_.nvx = 0;
estim_params_.nvn = 0;
estim_params_.ncx = 0;
@ -431,13 +434,13 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation
% covariance (a diagonal matrix) %%Except for infinite variances ;-)
varinit = 'prior';
if strcmpi(varinit,'prior')
stdev = bayestopt_.pstdev;
stdev = bayestopt_.p2;
indx = find(isinf(stdev));
stdev(indx) = ones(length(indx),1)*sqrt(10);
vars = stdev.^2;
CovJump = diag(vars);
elseif strcmpi(varinit,'eye')
vars = ones(length(bayestopt_.pstdev),1)*0.1;
vars = ones(length(bayestopt_.p2),1)*0.1;
CovJump = diag(vars);
else
disp('gmhmaxlik :: Error!')
@ -552,7 +555,7 @@ if options_.posterior_mode_estimation
invhess = inv(hh);
stdh = sqrt(diag(invhess));
else
variances = bayestopt_.pstdev.^2;
variances = bayestopt_.p2.^2;
invhess = 0.01*diag(variances);
%invhess = 0.01*eye(length(variances));
end
@ -575,9 +578,9 @@ if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
name = bayestopt_.name{ip};
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
name, ...
bayestopt_.pmean(ip),xparam1(ip),stdh(ip),tstath(ip), ...
bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
bayestopt_.p2(ip)));
eval(['oo_.posterior_mode.parameters.' name ' = xparam1(ip);']);
eval(['oo_.posterior_std.parameters.' name ' = stdh(ip);']);
ip = ip+1;
@ -591,9 +594,9 @@ if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
k = estim_params_.var_exo(i,1);
name = deblank(M_.exo_names(k,:));
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
name,bayestopt_.pmean(ip),xparam1(ip), ...
name,bayestopt_.p1(ip),xparam1(ip), ...
stdh(ip),tstath(ip),pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
bayestopt_.p2(ip)));
M_.Sigma_e(k,k) = xparam1(ip)*xparam1(ip);
eval(['oo_.posterior_mode.shocks_std.' name ' = xparam1(ip);']);
eval(['oo_.posterior_std.shocks_std.' name ' = stdh(ip);']);
@ -607,10 +610,10 @@ if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
for i=1:nvn
name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
name,bayestopt_.pmean(ip), ...
name,bayestopt_.p1(ip), ...
xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pstdev(ip)));
bayestopt_.p2(ip)));
eval(['oo_.posterior_mode.measurement_errors_std.' name ' = xparam1(ip);']);
eval(['oo_.posterior_std.measurement_errors_std.' name ' = stdh(ip);']);
ip = ip+1;
@ -626,8 +629,8 @@ if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', name, ...
bayestopt_.pmean(ip),xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.pstdev(ip)));
bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.p2(ip)));
M_.Sigma_e(k1,k2) = xparam1(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
eval(['oo_.posterior_mode.shocks_corr.' NAME ' = xparam1(ip);']);
@ -645,8 +648,8 @@ if any(bayestopt_.pshape > 0) & options_.posterior_mode_estimation
name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', name, ...
bayestopt_.pmean(ip),xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.pstdev(ip)));
bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.p2(ip)));
eval(['oo_.posterior_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.posterior_std.measurement_errors_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
@ -772,8 +775,8 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
M_.param_names_tex(estim_params_.param_vals(i,1),:),...
deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),...
bayestopt_.p1(ip),...
bayestopt_.p2(ip),...
xparam1(ip),...
stdh(ip));
ip = ip + 1;
@ -808,8 +811,8 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,[ '$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n'],...
deblank(M_.exo_names_tex(k,:)),...
deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
bayestopt_.pmean(ip),...
bayestopt_.pstdev(ip),...
bayestopt_.p1(ip),...
bayestopt_.p2(ip),...
xparam1(ip), ...
stdh(ip));
ip = ip+1;
@ -843,8 +846,8 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
deblank(M_.endo_names_tex(idx,:)), ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
bayestopt_.pmean(ip), ...
bayestopt_.pstdev(ip),...
bayestopt_.p1(ip), ...
bayestopt_.p2(ip),...
xparam1(ip),...
stdh(ip));
ip = ip+1;
@ -878,8 +881,8 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,[ '$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n'],...
[deblank(M_.exo_names_tex(k1,:)) ',' deblank(M_.exo_names_tex(k2,:))], ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
bayestopt_.pmean(ip), ...
bayestopt_.pstdev(ip), ...
bayestopt_.p1(ip), ...
bayestopt_.p2(ip), ...
xparam1(ip), ...
stdh(ip));
ip = ip+1;
@ -913,8 +916,8 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
[deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:))], ...
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.pmean(ip), ...
bayestopt_.pstdev(ip), ...
bayestopt_.p1(ip), ...
bayestopt_.p2(ip), ...
xparam1(ip), ...
stdh(ip));
ip = ip+1;

View File

@ -123,7 +123,7 @@ function forcst_unc(y0,var_list)
end
% compute shock uncertainty around forecast with mean prior
set_parameters(bayestopt_.pmean);
set_parameters(bayestopt_.p1);
[dr,info] = resol(oo_.steady_state,0);
[yf3,yf3_intv] = forcst(dr,y0,periods,var_list);
yf3_1 = yf3'-[zeros(maximum_lag,n); yf3_intv];

View File

@ -54,8 +54,8 @@ function get_prior_info(info)
for i=1:size(bayestopt_.name,1)
[tmp,TexName] = get_the_name(i,1);
PriorShape = PriorNames{ bayestopt_.pshape(i) };
PriorMean = bayestopt_.pmean(i);
PriorStandardDeviation = bayestopt_.pstdev(i);
PriorMean = bayestopt_.p1(i);
PriorStandardDeviation = bayestopt_.p2(i);
switch bayestopt_.pshape(i)
case { 1 , 5 }
LowerBound = bayestopt_.p3(i);
@ -102,10 +102,9 @@ function get_prior_info(info)
function format_string = build_format_string(bayestopt,i)
format_string = ['%s & %s & %6.4f &'];
if isinf(bayestopt.pstdev(i))
if isinf(bayestopt.p2(i))
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];

View File

@ -1,20 +1,19 @@
function ldens = lpdfgam(x,a,b);
% function ldens = lpdfgam(x,a,b)
% log GAMMA PDF
% Evaluates the logged GAMMA PDF at x.
%
% INPUTS
% x [double] m*n matrix of locations,
% a [double] m*n matrix or scalar, First GAMMA distribution parameters,
% b [double] m*n matrix or scalar, Second GAMMA distribution parameters,
%
% INPUTS
% x: density evatuated at x
% a: GAMMA distribution paramerer
% b: GAMMA distribution paramerer
%
% OUTPUTS
% ldens: the log GAMMA PDF
% ldens [double] m*n matrix of logged GAMMA densities evaluated at x.
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -31,8 +30,11 @@ function ldens = lpdfgam(x,a,b);
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
ldens = -gammaln(a) -a*log(b)+ (a-1)*log(x) -x/b ;
% 10/11/03 MJ adapted from an earlier GAUSS version by F. Schorfeide,
% translated to MATLAB by R. Wouters
% use MATLAB gammaln rather than lngam
ldens = -Inf( size(x) ) ;
idx = find( x>0 );
if length(a)==1
ldens(idx) = -gammaln(a) - a*log(b) + (a-1)*log(x(idx)) - x(idx)/b ;
else
ldens(idx) = -gammaln(a(idx)) - a(idx).*log(b(idx)) + (a(idx)-1).*log(x(idx)) - x(idx)./b(idx) ;
end

View File

@ -1,22 +1,20 @@
function ldens = lpdfgbeta(x,a,b,aa,bb);
% function ldens = lpdfgbeta(x,a,b,aa,bb);
% log (generalized) BETA PDF
% Evaluates the logged BETA PDF at x.
%
% INPUTS
% x: density evatuated at x
% a: BETA distribution paramerer
% b: BETA distribution paramerer
% aa: lower bound
% bb: upper bound
% OUTPUTS
% ldens: the log (generalized) BETA PDF
% INPUTS
% x [double] m*n matrix of loactions,
% a [double] m*n matrix of First BETA distribution parameters,
% b [double] m*n matrix of Second BETA distribution parameters,
% aa [double] m*n matrix of lower bounds,
% bb [double] m*n matrix of upper bounds.
%
% OUTPUTS
% ldens [double] m*n matrix of logged (generalized) BETA densities.
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -33,9 +31,11 @@ function ldens = lpdfgbeta(x,a,b,aa,bb);
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
ldens = -betaln(a,b) + (a-1)*log(x-aa) + (b-1)*log(bb-x) - (a+b-1)*log(bb-aa);
%gammaln(a+b) - gammaln(a) - gammaln(b)
%betaln(a,b)
%pause
% 02/16/04 SA Interval [aa,bb] is the support of the probability density function.
ldens = -Inf( size(x) ) ;
idx = find( (x-aa)>0 & (x-bb)<0 ) ;
if length(a)==1
ldens(idx) = -betaln(a,b) + (a-1)*log(x(idx)-aa) + (b-1)*log(bb-x(idx)) - (a+b-1)*log(bb-aa) ;
else
ldens(idx) = -betaln(a(idx),b(idx)) + (a(idx)-1).*log(x(idx)-aa(idx)) + (b(idx)-1).*log(bb(idx)-x(idx)) - (a(idx)+b(idx)-1).*log(bb(idx)-aa(idx));
end

View File

@ -1,22 +1,23 @@
function ldens = lpdfig1(x,s,nu)
% function ldens = lpdfig1(x,s,nu)
% log INVERSE GAMMA (type 1)
% X ~ IG1(s,nu)
% X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution)
% Evaluates the logged INVERSE-GAMMA-1 PDF at x.
%
% INPUTS
% x: density evatuated at x
% s: shape parameter
% nu: scale parameter
% X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution)
%
% See L. Bauwens, M. Lubrano and J-F. Richard [1999, appendix A] for more details.
%
%
% INPUTS
% x [double] m*n matrix of locations,
% s [double] m*n matrix or scalar, First INVERSE-GAMMA-1 distribution parameters,
% nu [double] m*n matrix or scalar, Second INVERSE-GAMMA-1 distribution parameters.
%
% OUTPUTS
% ldens: the log INVERSE GAMMA density function (type 1)
%
% ldens [double] m*n matrix of logged INVERSE-GAMMA-1 densities evaluated at x.
%
% SPECIAL REQUIREMENTS
% See L. Bauwens, M. Lubrano and J-F. Richard [1999, appendix A] for more
% details.
% none
% Copyright (C) 2004-2008 Dynare Team
% Copyright (C) 2004-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -33,4 +34,11 @@ function ldens = lpdfig1(x,s,nu)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
ldens = log(2) - gammaln(nu/2) - (nu/2).*log(2/s) - (nu+1)*log(x) - .5*s./(x.^2);
ldens = -Inf( size(x) ) ;
idx = find( x>0 ) ;
if length(s)==1
ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ;
else
ldens(idx) = log(2) - gammaln(.5*nu(idx)) - .5*nu(idx).*(log(2)-log(s(idx))) - (nu(idx)+1).*log(x(idx)) - .5*s(idx)./(x(idx).*x(idx)) ;
end

View File

@ -1,22 +1,23 @@
function ldens = lpdfig2(x,s,nu)
% function ldens = lpdfig2(x,s,nu)
% log INVERSE GAMMA (type 2)
% X ~ IG2(s,nu)
% X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution)
% Evaluates the logged INVERSE-GAMMA-2 PDF at x.
%
% INPUTS
% x: density evatuated at x
% s: shape parameter
% nu: scale parameter
% X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution)
%
% See L. Bauwens, M. Lubrano and J-F. Richard [1999, appendix A] for more details.
%
%
% INPUTS
% x [double] m*n matrix of locations,
% s [double] m*n matrix or scalar, First INVERSE-GAMMA-2 distribution parameters,
% nu [double] m*n matrix or scalar, Second INVERSE-GAMMA-2 distribution parameters.
%
% OUTPUTS
% ldens: the log INVERSE GAMMA density function (type 2)
%
% ldens [double] m*n matrix of logged INVERSE-GAMMA-2 densities evaluated at x.
%
% SPECIAL REQUIREMENTS
% See L. Bauwens, M. Lubrano and J-F. Richard [1999, appendix A] for more
% details.
% none
% Copyright (C) 2004-2008 Dynare Team
% Copyright (C) 2004-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -33,4 +34,11 @@ function ldens = lpdfig2(x,s,nu)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
ldens = - gammaln(nu/2) - (nu/2)*log(2/s) - .5*(nu+2)*log(x) -.5*s./x;
ldens = -Inf( size(x) ) ;
idx = find( x>0 ) ;
if length(s)==1
ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx);
else
ldens(idx) = -gammaln(.5*nu(idx)) - (.5*nu(idx)).*(log(2)-log(s(idx))) - .5*(nu(idx)+2).*log(x(idx)) -.5*s(idx)./x(idx);
end

View File

@ -1,19 +1,19 @@
function f = lpdfnorm(x,m,s)
% function f = lpdfnorm(x,m,s)
% The log of the normal density function
function ldens = lpdfnorm(x,a,b)
% Evaluates the logged UNIVARIATE GAUSSIAN PDF at x.
%
% INPUTS
% x: density evatuated at x
% m: mean
% s: standard deviation
% INPUTS
% x [double] m*n matrix of locations,
% a [double] m*n matrix or scalar, First GAUSSIAN distribution parameters (expectation)
% b [double] m*n matrix or scalar, Second GAUSSIAN distribution parameters (standard deviation).
%
% OUTPUTS
% f: the log of the normal density function
% OUTPUTS
% ldens [double] m*n matrix of logged GAUSSIAN densities evaluated at x.
%
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -30,7 +30,6 @@ function f = lpdfnorm(x,m,s)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<3, s=1; end
if nargin<2, m=0; end
f = -log(s)-log(2*pi)/2-((x-m)./s).^2/2;
if nargin<3, b=1; end
if nargin<2, a=0; end
ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;

View File

@ -13,7 +13,7 @@ function plot_priors(bayestopt_,M_,options_)
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2004-2008 Dynare Team
% Copyright (C) 2004-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -33,7 +33,7 @@ function plot_priors(bayestopt_,M_,options_)
TeX = options_.TeX;
figurename = 'Priors';
npar = length(bayestopt_.pmean);
npar = length(bayestopt_.p1);
[nbplt,nr,nc,lr,lc,nstar] = pltorg(npar);
if TeX

View File

@ -11,7 +11,7 @@ function bounds = prior_bounds(bayestopt)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -31,41 +31,33 @@ function bounds = prior_bounds(bayestopt)
global options_
pshape = bayestopt.pshape;
pmean = bayestopt.pmean;
p1 = bayestopt.p1;
p2 = bayestopt.p2;
p3 = bayestopt.p3;
p4 = bayestopt.p4;
p6 = bayestopt.p6;
p7 = bayestopt.p7;
n = length(pmean);
bounds = zeros(n,2);
bounds = zeros(length(p6),2);
for i=1:n
for i=1:length(p6)
switch pshape(i)
case 1
mu = (pmean(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
A = (1-mu)*mu^2/stdd^2 - mu;
B = A*(1/mu - 1);
bounds(i,1) = betainv(options_.prior_trunc,A,B)*(p4(i)-p3(i))+p3(i);
bounds(i,2) = betainv(1-options_.prior_trunc,A,B)*(p4(i)-p3(i))+p3(i);
bounds(i,1) = betainv(options_.prior_trunc,p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
bounds(i,2) = betainv(1-options_.prior_trunc,p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
case 2
b = p2(i)^2/(pmean(i)-p3(i));
a = (pmean(i)-p3(i))/b;
bounds(i,1) = gaminv(options_.prior_trunc,a,b)+p3(i);
bounds(i,2) = gaminv(1-options_.prior_trunc,a,b)+p3(i);
bounds(i,1) = gaminv(options_.prior_trunc,p6(i),p7(i))+p3(i);
bounds(i,2) = gaminv(1-options_.prior_trunc,p6(i),p7(i))+p3(i);
case 3
bounds(i,1) = norminv(options_.prior_trunc,pmean(i),p2(i));
bounds(i,2) = norminv(1-options_.prior_trunc,pmean(i),p2(i));
bounds(i,1) = norminv(options_.prior_trunc,p6(i),p7(i));
bounds(i,2) = norminv(1-options_.prior_trunc,p6(i),p7(i));
case 4
bounds(i,1) = 1/sqrt(gaminv(1-options_.prior_trunc, p2(i)/2, 2/p1(i)));
bounds(i,2) = 1/sqrt(gaminv(options_.prior_trunc, p2(i)/2, 2/p1(i)));
bounds(i,1) = 1/sqrt(gaminv(1-options_.prior_trunc, p7(i)/2, 2/p6(i)))+p3(i);
bounds(i,2) = 1/sqrt(gaminv(options_.prior_trunc, p7(i)/2, 2/p6(i)))+p3(i);
case 5
bounds(i,1) = p1(i)+(p2(i)-p1(i))*options_.prior_trunc;
bounds(i,2) = p2(i)-(p2(i)-p1(i))*options_.prior_trunc;
bounds(i,1) = p6(i)+(p7(i)-p6(i))*options_.prior_trunc;
bounds(i,2) = p7(i)-(p7(i)-p6(i))*options_.prior_trunc;
case 6
bounds(i,1) = 1/gaminv(1-options_.prior_trunc, p2(i)/2, 2/p1(i));
bounds(i,2) = 1/gaminv(options_.prior_trunc, p2(i)/2, 2/p1(i));
bounds(i,1) = 1/gaminv(1-options_.prior_trunc, p7(i)/2, 2/p6(i))+p3(i);
bounds(i,2) = 1/gaminv(options_.prior_trunc, p7(i)/2, 2/p6(i))+p3(i);
otherwise
error(sprintf('prior_bounds: unknown distribution shape (index %d, type %d)', i, pshape(i)));
end

View File

@ -32,63 +32,73 @@ function pdraw = prior_draw(init, prior_structure)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent prior_mean prior_standard_deviation a b p1 p2 p3 p4
persistent p6 p7 p3 p4
persistent uniform_index gaussian_index gamma_index beta_index inverse_gamma_1_index inverse_gamma_2_index
persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws
if nargin>0 && init
prior_shape = prior_structure.pshape;
prior_mean = prior_structure.pmean;
prior_standard_deviation = prior_structure.pstdev;
p1 = prior_structure.p1;
p2 = prior_structure.p2;
p6 = prior_structure.p6;
p7 = prior_structure.p7;
p3 = prior_structure.p3;
p4 = prior_structure.p4;
number_of_estimated_parameters = length(p1);
a = NaN(number_of_estimated_parameters,1);
b = NaN(number_of_estimated_parameters,1);
number_of_estimated_parameters = length(p6);
beta_index = find(prior_shape==1);
beta_draws = 1;
if isempty(beta_index)
beta_draws = 0;
end
gamma_index = find(prior_shape==2);
gamma_draws = 1;
if isempty(gamma_index)
gamma_draws = 0;
end
gaussian_index = find(prior_shape==3);
gaussian_draws = 1;
if isempty(gaussian_index)
gaussian_draws = 0;
end
inverse_gamma_1_index = find(prior_shape==4);
inverse_gamma_1_draws = 1;
if isempty(inverse_gamma_1_index)
inverse_gamma_1_draws = 0;
end
uniform_index = find(prior_shape==5);
uniform_draws = 1;
if isempty(uniform_index)
uniform_draws = 0;
end
inverse_gamma_2_index = find(prior_shape==6);
% Set parameters for the beta prior
mu = (p1(beta_index)-p3(beta_index))./(p4(beta_index)-p3(beta_index));
stdd = p2(beta_index)./(p4(beta_index)-p3(beta_index));
a(beta_index) = (1-mu).*mu.^2./stdd.^2 - mu;
b(beta_index) = a(beta_index).*(1./mu - 1);
% Set parameters for the gamma prior
mu = p1(gamma_index)-p3(gamma_index);
b(gamma_index) = p2(gamma_index).^2./mu;
a(gamma_index) = mu./b(gamma_index);
% Initialization of the vector of prior draws.
inverse_gamma_2_draws = 1;
if isempty(inverse_gamma_2_index)
inverse_gamma_2_draws = 0;
end
pdraw = zeros(number_of_estimated_parameters,1);
return
end
% Uniform draws.
if ~isempty(uniform_index)
if uniform_draws
pdraw(uniform_index) = rand(length(uniform_index),1).*(p4(uniform_index)-p3(uniform_index)) + p3(uniform_index);
end
% Gaussian draws.
if ~isempty(gaussian_index)
pdraw(gaussian_index) = randn(length(gaussian_index),1).*prior_standard_deviation(gaussian_index) + prior_mean(gaussian_index);
if gaussian_draws
pdraw(gaussian_index) = randn(length(gaussian_index),1).*p7(gaussian_index) + p6(gaussian_index);
end
% Gamma draws.
if ~isempty(gamma_index)
pdraw(gamma_index) = gamrnd(a(gamma_index),b(gamma_index))+p3(gamma_index);
if gamma_draws
pdraw(gamma_index) = gamrnd(p6(gamma_index),p7(gamma_index))+p3(gamma_index);
end
% Beta draws.
if ~isempty(beta_index)
pdraw(beta_index) = (p4(beta_index)-p3(beta_index)).*betarnd(a(beta_index),b(beta_index))+p3(beta_index);
if beta_draws
pdraw(beta_index) = (p4(beta_index)-p3(beta_index)).*betarnd(p6(beta_index),p7(beta_index))+p3(beta_index);
end
% Inverted gamma (type 1) draws.
if ~isempty(inverse_gamma_1_index)
if inverse_gamma_1_draws
pdraw(inverse_gamma_1_index) = ...
sqrt(1./gamrnd(p2(inverse_gamma_1_index)/2,2./p1(inverse_gamma_1_index)));
sqrt(1./gamrnd(p7(inverse_gamma_1_index)/2,2./p6(inverse_gamma_1_index)))+p3(inverse_gamma_1_index);
end
% Inverted gamma (type 2) draws.
if ~isempty(inverse_gamma_2_index)
if inverse_gamma_2_draws
pdraw(inverse_gamma_2_index) = ...
1./gamrnd(p2(inverse_gamma_2_index)/2,2./p1(inverse_gamma_2_index));
1./gamrnd(p7(inverse_gamma_2_index)/2,2./p6(inverse_gamma_2_index))+p3(inverse_gamma_2_index);
end

View File

@ -45,7 +45,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
count_dll_problem = 0;
count_unknown_problem = 0;
NumberOfSimulations = options_.prior_mc;
NumberOfParameters = length(bayestopt_.pmean);
NumberOfParameters = length(bayestopt_.p1);
NumberOfEndogenousVariables = size(M_.endo_names,1);
NumberOfElementsPerFile = ceil(options_.MaxNumberOfBytes/NumberOfParameters/NumberOfEndogenousVariables/8)
if NumberOfSimulations <= NumberOfElementsPerFile

View File

@ -1,29 +1,19 @@
function lnprior = priordens(para, pshape, p1, p2, p3, p4)
% function lnprior = priordens(para, pshape, p1, p2, p3, p4)
% computes a prior density for the structural parameters of DSGE models
function logged_prior_density = priordens(x, pshape, p6, p7, p3, p4)
% Computes a prior density for the structural parameters of DSGE models
%
% INPUTS
% para: parameter value
% pshape: 0 is point mass, both para and p2 are ignored
% 1 is BETA
% 2 is GAMMA
% 3 is NORMAL
% 4 is INVERTED GAMMA TYPE I
% 5 is UNIFORM
% 6 is INVERTED GAMMA TYPE II
% p1: mean
% p2: standard deviation
% p3: lower bound
% p4: upper bound
% INPUTS
% x [double] vector with n elements.
% pshape [integer] vector with n elements (bayestopt_.pshape).
% p6: [double] vector with n elements, first parameter of the prior distribution (bayestopt_.p6).
% p7: [double] vector with n elements, second parameter of the prior distribution (bayestopt_.p7).
% p3: [double] vector with n elements, lower bounds.
% p4: [double] vector with n elements, upper bound.
%
% OUTPUTS
% lnprior: log of the prior density
% OUTPUTS
% logged_prior_density [double] scalar, log of the prior density evaluated at x.
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -40,39 +30,89 @@ function lnprior = priordens(para, pshape, p1, p2, p3, p4)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
lnprior = 0;
nprio = length(pshape);
persistent pflag
persistent id1 id2 id3 id4 id5 id6
persistent tt1 tt2 tt3 tt4 tt5 tt6
i = 1;
while i <= nprio;
a = 0;
b = 0;
if pshape(i) == 1; % (generalized) BETA Prior
mu = (p1(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
a = (1-mu)*mu^2/stdd^2 - mu;
b = a*(1/mu - 1);
lnprior = lnprior + lpdfgbeta(para(i),a,b,p3(i),p4(i)) ;
elseif pshape(i) == 2; % GAMMA PRIOR
b = p2(i)^2/(p1(i)-p3(i));
a = (p1(i)-p3(i))/b;
lnprior = lnprior + lpdfgam(para(i)-p3(i),a,b);
elseif pshape(i) == 3; % GAUSSIAN PRIOR
lnprior = lnprior + lpdfnorm(para(i),p1(i),p2(i));
elseif pshape(i) == 4; % INVGAMMA1 PRIOR
lnprior = lnprior + lpdfig1(para(i),p1(i),p2(i));
elseif pshape(i) == 5; % UNIFORM PRIOR
lnprior = lnprior + log(1/(p2(i)-p1(i)));
elseif pshape(i) == 6; % INVGAMMA2 PRIOR
lnprior = lnprior + lpdfig2(para(i),p1(i),p2(i));
end;
i = i+1;
end;
if isempty(pflag)
Number0fParameters = length(pshape);
% Beta indices.
tt1 = 1;
id1 = find(pshape==1);
if isempty(id1)
tt1 = 0;
end
% Gamma indices.
tt2 = 1;
id2 = find(pshape==2);
if isempty(id2)
tt2 = 0;
end
% Gaussian indices.
tt3 = 1;
id3 = find(pshape==3);
if isempty(id3)
tt3 = 0;
end
% Inverse-Gamma-1 indices.
tt4 = 1;
id4 = find(pshape==4);
if isempty(id4)
tt4 = 0;
end
% Uniform indices.
tt5 = 1;
id5 = find(pshape==5);
if isempty(id5)
tt5 = 0;
end
% Inverse-Gamma-2 indices.
tt6 = 1;
id6 = find(pshape==6);
if isempty(id6)
tt6 = 0;
end
pflag = 1;
end
% 10/11/03 MJ adapted from an earlier version in GAUSS by F. Schorfheide
% and translated to Matlab by R. Wouters
% 11/18/03 MJ adopted M.Ratto's suggestion for inverse gamma
% changed order of input parameters
% 01/16/04 MJ added invgamma2
% for invgamma p2 is now standard error
% 16/02/04 SA changed beta prior call
logged_prior_density = 0.0;
if tt1
logged_prior_density = logged_prior_density + sum(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1))) ;
if isinf(logged_prior_density)
return
end
end
if tt2
logged_prior_density = logged_prior_density + sum(lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2))) ;
if isinf(logged_prior_density)
return
end
end
if tt3
logged_prior_density = logged_prior_density + sum(lpdfnorm(x(id3),p6(id3),p7(id3))) ;
end
if tt4
logged_prior_density = logged_prior_density + sum(lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4))) ;
if isinf(logged_prior_density)
return
end
end
if tt5
if any(x(id5)-p3(id5)<0) || any(x(id5)-p4(id5)>0)
logged_prior_density = -Inf ;
return
end
logged_prior_density = logged_prior_density + sum(log(1./(p4(id5)-p3(id5)))) ;
end
if tt6
logged_prior_density = logged_prior_density + sum(lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6))) ;
if isinf(logged_prior_density)
return
end
end

View File

@ -11,7 +11,7 @@ function y = rndprior(bayestopt_)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -29,52 +29,29 @@ function y = rndprior(bayestopt_)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
pshape=bayestopt_.pshape;
pmean=bayestopt_.pmean;
p1=bayestopt_.p1;
p2=bayestopt_.p2;
p3=bayestopt_.p3;
p4=bayestopt_.p4;
for i=1:length(pmean),
p6=bayestopt_.p6;
p7=bayestopt_.p7;
y = NaN(1,length(pshape));
for i=1:length(pshape)
switch pshape(i)
case 1 % Beta
mu = (pmean(i)-p3(i))/(p4(i)-p3(i));
stdd = p2(i)/(p4(i)-p3(i));
A = (1-mu)*mu^2/stdd^2 - mu;
B = A*(1/mu - 1);
y(1,i) = betarnd(A, B);
y(1,i) = y(1,i) * (p4(i)-p3(i)) + p3(i);
case 2 % Generalized gamma
mu = pmean(i)-p3(i);
B = p2(i)^2/mu;
A = mu/B;
y(1,i) = gamrnd(A, B) + p3(i);
case 3 % Gaussian
MU = pmean(i);
SIGMA = p2(i);
y(1,i) = randn*SIGMA+ MU;
case 4 % Inverse-gamma type 1
nu = p2(i);
s = p1(i);
y(1,i) = 1/sqrt(gamrnd(nu/2, 2/s));
case 5 % Uniform
y(1,i) = rand*(p2(i)-p1(i)) + p1(i);
case 6 % Inverse-gamma type 2
nu = p2(i);
s = p1(i);
y(1,i) = 1/gamrnd(nu/2, 2/s);
otherwise
error(sprintf('rndprior: unknown distribution shape (index %d, type %d)', i, pshape(i)));
case 1 % Beta
y(i) = betarnd(p6(i),p7(i));
y(i) = y(1,i) * (p4(i)-p3(i)) + p3(i);
case 2 % Generalized gamma
y(i) = gamrnd(p6(i),p7(i)) + p3(i);
case 3 % Gaussian
y(i) = randn*p7(i) + p6(i) ;
case 4 % Inverse-gamma type 1
y(i) = 1/sqrt(gamrnd(p7(i)/2, 2/p6(i))) + p3(i);
case 5 % Uniform
y(i) = rand*(p4(i)-p3(i)) + p3(i);
case 6 % Inverse-gamma type 2
y(i) = 1/gamrnd(p7(i)/2, 2/p6(i)) + p3(i);
otherwise
error(sprintf('rndprior: unknown distribution shape (index %d, type %d)', i, pshape(i)));
end
end
% initial version by Marco Ratto
end

View File

@ -18,7 +18,7 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2003-2008 Dynare Team
% Copyright (C) 2003-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -35,165 +35,190 @@ function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nvx = size(estim_params_.var_exo,1);
nvn = size(estim_params_.var_endo,1);
ncx = size(estim_params_.corrx,1);
ncn = size(estim_params_.corrn,1);
np = size(estim_params_.param_vals,1);
nvx = size(estim_params_.var_exo,1);
nvn = size(estim_params_.var_endo,1);
ncx = size(estim_params_.corrx,1);
ncn = size(estim_params_.corrn,1);
np = size(estim_params_.param_vals,1);
estim_params_.nvx = nvx;
estim_params_.nvn = nvn;
estim_params_.ncx = ncx;
estim_params_.ncn = ncn;
estim_params_.np = np;
estim_params_.nvx = nvx;
estim_params_.nvn = nvn;
estim_params_.ncx = ncx;
estim_params_.ncn = ncn;
estim_params_.np = np;
xparam1 = [];
ub = [];
lb = [];
bayestopt_.pshape = [];
bayestopt_.pmean = [];
bayestopt_.pstdev = [];
bayestopt_.p1 = [];
bayestopt_.p2 = [];
bayestopt_.p3 = [];
bayestopt_.p4 = [];
bayestopt_.jscale = [];
bayestopt_.name = [];
if nvx
xparam1 = estim_params_.var_exo(:,2);
ub = estim_params_.var_exo(:,4);
lb = estim_params_.var_exo(:,3);
bayestopt_.pshape = estim_params_.var_exo(:,5);
bayestopt_.pmean = estim_params_.var_exo(:,6);
bayestopt_.pstdev = estim_params_.var_exo(:,7);
bayestopt_.p3 = estim_params_.var_exo(:,8);
bayestopt_.p4 = estim_params_.var_exo(:,9);
bayestopt_.jscale = estim_params_.var_exo(:,10);
bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
end
if nvn
if M_.H == 0
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
xparam1 = [];
ub = [];
lb = [];
bayestopt_.pshape = [];
bayestopt_.p1 = []; % prior mean
bayestopt_.p2 = []; % prior standard deviation
bayestopt_.p3 = []; % lower bound
bayestopt_.p4 = []; % upper bound
bayestopt_.p5 = []; % prior mode
bayestopt_.p6 = []; % first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
bayestopt_.p7 = []; % second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
bayestopt_.jscale = [];
bayestopt_.name = [];
if nvx
xparam1 = estim_params_.var_exo(:,2);
ub = estim_params_.var_exo(:,4);
lb = estim_params_.var_exo(:,3);
bayestopt_.pshape = estim_params_.var_exo(:,5);
bayestopt_.p1 = estim_params_.var_exo(:,6);
bayestopt_.p2 = estim_params_.var_exo(:,7);
bayestopt_.p3 = estim_params_.var_exo(:,8);
bayestopt_.p4 = estim_params_.var_exo(:,9);
bayestopt_.jscale = estim_params_.var_exo(:,10);
bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
end
for i=1:nvn
estim_params_.var_endo(i,1) = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact');
if nvn
if M_.H == 0
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
end
for i=1:nvn
estim_params_.var_endo(i,1) = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact');
end
xparam1 = [xparam1; estim_params_.var_endo(:,2)];
ub = [ub; estim_params_.var_endo(:,4)];
lb = [lb; estim_params_.var_endo(:,3)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.var_endo(:,5)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.var_endo(:,6)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.var_endo(:,7)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name), M_.endo_names(estim_params_.var_endo(:,1),:)));
end
xparam1 = [xparam1; estim_params_.var_endo(:,2)];
ub = [ub; estim_params_.var_endo(:,4)];
lb = [lb; estim_params_.var_endo(:,3)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.var_endo(:,5)];
bayestopt_.pmean = [ bayestopt_.pmean; estim_params_.var_endo(:,6)];
bayestopt_.pstdev = [ bayestopt_.pstdev; estim_params_.var_endo(:,7)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name),...
M_.endo_names(estim_params_.var_endo(:,1),:)));
end
if ncx
xparam1 = [xparam1; estim_params_.corrx(:,3)];
ub = [ub; max(min(estim_params_.corrx(:,5),1),-1)];
lb = [lb; max(min(estim_params_.corrx(:,4),1),-1)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrx(:,6)];
bayestopt_.pmean = [ bayestopt_.pmean; estim_params_.corrx(:,7)];
bayestopt_.pstdev = [ bayestopt_.pstdev; estim_params_.corrx(:,8)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrx(:,9)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrx(:,10)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrx(:,11)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name),...
char(strcat(cellstr(M_.exo_names(estim_params_.corrx(:,1),:)),...
',',...
cellstr(M_.exo_names(estim_params_.corrx(:,2),:))))));
end
if ncn
if M_.H == 0
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
if ncx
xparam1 = [xparam1; estim_params_.corrx(:,3)];
ub = [ub; max(min(estim_params_.corrx(:,5),1),-1)];
lb = [lb; max(min(estim_params_.corrx(:,4),1),-1)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrx(:,6)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrx(:,7)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrx(:,8)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrx(:,9)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrx(:,10)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrx(:,11)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name), char(strcat(cellstr(M_.exo_names(estim_params_.corrx(:,1),:)), ...
',' , cellstr(M_.exo_names(estim_params_.corrx(:,2),:))))));
end
if ncn
if M_.H == 0
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
end
xparam1 = [xparam1; estim_params_.corrn(:,3)];
ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
lb = [lb; max(min(estim_params_.corrn(:,4),1),-1)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrn(:,6)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrn(:,7)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrn(:,8)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrn(:,9)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name), char(strcat(cellstr(M_.endo_names(estim_params_.corrn(:,1),:)),...
',' , cellstr(M_.endo_names(estim_params_.corrn(:,2),:))))));
end
if np
xparam1 = [xparam1; estim_params_.param_vals(:,2)];
ub = [ub; estim_params_.param_vals(:,4)];
lb = [lb; estim_params_.param_vals(:,3)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.param_vals(:,5)];
bayestopt_.p1 = [ bayestopt_.p1; estim_params_.param_vals(:,6)];
bayestopt_.p2 = [ bayestopt_.p2; estim_params_.param_vals(:,7)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.param_vals(:,8)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.param_vals(:,9)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.param_vals(:,10)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name),M_.param_names(estim_params_.param_vals(:,1),:)));
end
xparam1 = [xparam1; estim_params_.corrn(:,3)];
ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
lb = [lb; max(min(estim_params_.corrn(:,4),1),-1)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrn(:,6)];
bayestopt_.pmean = [ bayestopt_.pmean; estim_params_.corrn(:,7)];
bayestopt_.pstdev = [ bayestopt_.pstdev; estim_params_.corrn(:,8)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrn(:,9)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name),...
char(strcat(cellstr(M_.endo_names(estim_params_.corrn(:,1),:)),...
',',...
cellstr(M_.endo_names(estim_params_.corrn(:,2),:))))));
end
if np
xparam1 = [xparam1; estim_params_.param_vals(:,2)];
ub = [ub; estim_params_.param_vals(:,4)];
lb = [lb; estim_params_.param_vals(:,3)];
bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.param_vals(:,5)];
bayestopt_.pmean = [ bayestopt_.pmean; estim_params_.param_vals(:,6)];
bayestopt_.pstdev = [ bayestopt_.pstdev; estim_params_.param_vals(:,7)];
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.param_vals(:,8)];
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.param_vals(:,9)];
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.param_vals(:, ...
10)];
bayestopt_.name = cellstr(strvcat(char(bayestopt_.name),M_.param_names(estim_params_.param_vals(:,1),:)));
end
bayestopt_.ub = ub;
bayestopt_.lb = lb;
bayestopt_.ub = ub;
bayestopt_.lb = lb;
bayestopt_.p6 = NaN(size(bayestopt_.p1)) ;
bayestopt_.p7 = bayestopt_.p6 ;
bayestopt_.p1 = bayestopt_.pmean;
bayestopt_.p2 = bayestopt_.pstdev;
% generalized location parameters by default for beta distribution
k = find(bayestopt_.pshape == 1);
k1 = find(isnan(bayestopt_.p3(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
k1 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p4(k(k1)) = ones(length(k1),1);
% generalized location parameter by default for gamma distribution
k = find(bayestopt_.pshape == 2);
k1 = find(isnan(bayestopt_.p3(k)));
k2 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
% truncation parameters by default for normal distribution
k = find(bayestopt_.pshape == 3);
k1 = find(isnan(bayestopt_.p3(k)));
bayestopt_.p3(k(k1)) = -Inf*ones(length(k1),1);
k1 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p4(k(k1)) = Inf*ones(length(k1),1);
% generalized location parameters by default for beta distribution
k = find(bayestopt_.pshape == 1);
k1 = find(isnan(bayestopt_.p3(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
k1 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p4(k(k1)) = ones(length(k1),1);
for i=1:length(k)
mu = (bayestopt_.p1(k(i))-bayestopt_.p3(k(i)))/(bayestopt_.p4(k(i))-bayestopt_.p3(k(i)));
stdd = bayestopt_.p2(k(i))/(bayestopt_.p4(k(i))-bayestopt_.p3(k(i)));
bayestopt_.p6(k(i)) = (1-mu)*mu^2/stdd^2 - mu ;
bayestopt_.p7(k(i)) = bayestopt_.p6(k(i))*(1/mu-1) ;
m = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) ],1);
if length(m)==1
bayestopt_.p5(k(i)) = m;
else
disp(['Prior distribution for parameter ' int2str(k(i)) ' has two modes!'])
bayestopt_.p5(k(i)) = bayestopt_.p1(k(i)) ;
end
end
% generalized location parameter by default for gamma distribution
k = find(bayestopt_.pshape == 2);
k1 = find(isnan(bayestopt_.p3(k)));
k2 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
mu = bayestopt_.p1(k(i))-bayestopt_.p3(k(i));
bayestopt_.p7(k(i)) = bayestopt_.p2(k(i))^2/mu ;
bayestopt_.p6(k(i)) = mu/bayestopt_.p7(k(i)) ;
bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) ], 2) ;
end
% truncation parameters by default for normal distribution
k = find(bayestopt_.pshape == 3);
k1 = find(isnan(bayestopt_.p3(k)));
bayestopt_.p3(k(k1)) = -Inf*ones(length(k1),1);
k1 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p4(k(k1)) = Inf*ones(length(k1),1);
for i=1:length(k)
bayestopt_.p6(k(i)) = bayestopt_.p1(k(i)) ;
bayestopt_.p7(k(i)) = bayestopt_.p2(k(i)) ;
bayestopt_.p5(k(i)) = bayestopt_.p1(k(i)) ;
end
% inverse gamma distribution
k = find(bayestopt_.pshape == 4);
for i=1:length(k)
[bayestopt_.p1(k(i)),bayestopt_.p2(k(i))] = ...
inverse_gamma_specification(bayestopt_.pmean(k(i)),bayestopt_.pstdev(k(i)),1);
end
k1 = find(isnan(bayestopt_.p3(k)));
k2 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
% uniform distribution
k = find(bayestopt_.pshape == 5);
for i=1:length(k)
[bayestopt_.pmean(k(i)),bayestopt_.pstdev(k(i)),bayestopt_.p1(k(i)),bayestopt_.p2(k(i))] = ...
uniform_specification(bayestopt_.pmean(k(i)),bayestopt_.pstdev(k(i)),bayestopt_.p3(k(i)),bayestopt_.p4(k(i)));
end
% inverse gamma distribution (type 2)
k = find(bayestopt_.pshape == 6);
for i=1:length(k)
[bayestopt_.p1(k(i)),bayestopt_.p2(k(i))] = ...
inverse_gamma_specification(bayestopt_.pmean(k(i)),bayestopt_.pstdev(k(i)),2);
end
k1 = find(isnan(bayestopt_.p3(k)));
k2 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
k = find(isnan(xparam1));
xparam1(k) = bayestopt_.pmean(k);
% inverse gamma distribution
k = find(bayestopt_.pshape == 4);
k1 = find(isnan(bayestopt_.p3(k)));
k2 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
[bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),1) ;
bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4) ;
end
% uniform distribution
k = find(bayestopt_.pshape == 5);
for i=1:length(k)
[bayestopt_.p1(k(i)),bayestopt_.p2(k(i)),bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
uniform_specification(bayestopt_.p1(k(i)),bayestopt_.p2(k(i)),bayestopt_.p3(k(i)),bayestopt_.p4(k(i)));
bayestopt_.p3(k(i)) = bayestopt_.p6(k(i)) ;
bayestopt_.p4(k(i)) = bayestopt_.p7(k(i)) ;
bayestopt_.p5(k(i)) = NaN ;
end
% inverse gamma distribution (type 2)
k = find(bayestopt_.pshape == 6);
k1 = find(isnan(bayestopt_.p3(k)));
k2 = find(isnan(bayestopt_.p4(k)));
bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
[bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),2);
bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ;
end
k = find(isnan(xparam1));
xparam1(k) = bayestopt_.p1(k);

View File

@ -1,6 +1,4 @@
function [m,s,p1,p2] = uniform_specification(m,s,p3,p4)
% function [m,s,p1,p2] = uniform_specification(m,s,p3,p4)
function [m,s,p6,p7] = uniform_specification(m,s,p3,p4)
% Specification of the uniform density function parameters
%
% INPUTS
@ -18,7 +16,7 @@ function [m,s,p1,p2] = uniform_specification(m,s,p3,p4)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2004-2008 Dynare Team
% Copyright (C) 2004-2009 Dynare Team
%
% This file is part of Dynare.
%
@ -36,11 +34,11 @@ function [m,s,p1,p2] = uniform_specification(m,s,p3,p4)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if ~(isnan(p3) | isnan(p4))
p1 = p3;
p2 = p4;
m = (p3+p4)/2;
s = (p4-p3)/(sqrt(12));
p6 = p3;
p7 = p4;
m = (p3+p4)/2;
s = (p4-p3)/(sqrt(12));
else
p1 = m-s*sqrt(3);
p2 = m+s*sqrt(3);
end
p6 = m-s*sqrt(3);
p7 = m+s*sqrt(3);
end