diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index 60ce7d1b3..e41bd4bdd 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -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;
\ No newline at end of file
diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m
index 82d1a3784..bc8360ece 100644
--- a/matlab/DsgeLikelihood_hh.m
+++ b/matlab/DsgeLikelihood_hh.m
@@ -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)];
diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index a31ab9dc4..b1916cae4 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -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)
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index a23363712..00664f3e1 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -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
diff --git a/matlab/distributions/compute_prior_mode.m b/matlab/distributions/compute_prior_mode.m
new file mode 100644
index 000000000..121c45261
--- /dev/null
+++ b/matlab/distributions/compute_prior_mode.m
@@ -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 .
+ 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
\ No newline at end of file
diff --git a/matlab/draw_prior_density.m b/matlab/draw_prior_density.m
index 085cc25f3..4d5154aca 100644
--- a/matlab/draw_prior_density.m
+++ b/matlab/draw_prior_density.m
@@ -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 .
-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)];
diff --git a/matlab/dsge_posterior_kernel.m b/matlab/dsge_posterior_kernel.m
index f6d52d2b5..8fd2c95fd 100644
--- a/matlab/dsge_posterior_kernel.m
+++ b/matlab/dsge_posterior_kernel.m
@@ -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;
\ No newline at end of file
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index ac872d32c..9b795ca52 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -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;
diff --git a/matlab/forcst_unc.m b/matlab/forcst_unc.m
index b4f243060..ffaf1295d 100644
--- a/matlab/forcst_unc.m
+++ b/matlab/forcst_unc.m
@@ -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];
diff --git a/matlab/get_prior_info.m b/matlab/get_prior_info.m
index b9e5df3c2..a2df20ad3 100644
--- a/matlab/get_prior_info.m
+++ b/matlab/get_prior_info.m
@@ -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 &'];
diff --git a/matlab/lpdfgam.m b/matlab/lpdfgam.m
index c3a4e778a..7f7690b2a 100644
--- a/matlab/lpdfgam.m
+++ b/matlab/lpdfgam.m
@@ -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 .
- 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
\ No newline at end of file
diff --git a/matlab/lpdfgbeta.m b/matlab/lpdfgbeta.m
index 7623bc44e..75a5e8fec 100644
--- a/matlab/lpdfgbeta.m
+++ b/matlab/lpdfgbeta.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/lpdfig1.m b/matlab/lpdfig1.m
index cf5446448..4b7b061a6 100644
--- a/matlab/lpdfig1.m
+++ b/matlab/lpdfig1.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/lpdfig2.m b/matlab/lpdfig2.m
index 1cc398204..59c4d07c7 100644
--- a/matlab/lpdfig2.m
+++ b/matlab/lpdfig2.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/lpdfnorm.m b/matlab/lpdfnorm.m
index 0279bb399..0e232884a 100644
--- a/matlab/lpdfnorm.m
+++ b/matlab/lpdfnorm.m
@@ -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 .
-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) ;
\ No newline at end of file
diff --git a/matlab/plot_priors.m b/matlab/plot_priors.m
index 2ed37a56d..1a10e2ddd 100644
--- a/matlab/plot_priors.m
+++ b/matlab/plot_priors.m
@@ -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
diff --git a/matlab/prior_bounds.m b/matlab/prior_bounds.m
index d4b817817..610c42e07 100644
--- a/matlab/prior_bounds.m
+++ b/matlab/prior_bounds.m
@@ -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
diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m
index 0059c6a2f..6614b0829 100644
--- a/matlab/prior_draw.m
+++ b/matlab/prior_draw.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/prior_sampler.m b/matlab/prior_sampler.m
index ac4626eb3..e79a1cd1f 100644
--- a/matlab/prior_sampler.m
+++ b/matlab/prior_sampler.m
@@ -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
diff --git a/matlab/priordens.m b/matlab/priordens.m
index a59196caf..1b8920220 100644
--- a/matlab/priordens.m
+++ b/matlab/priordens.m
@@ -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 .
-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
\ No newline at end of file
diff --git a/matlab/rndprior.m b/matlab/rndprior.m
index 5b45aa8aa..3ff7f6e0b 100644
--- a/matlab/rndprior.m
+++ b/matlab/rndprior.m
@@ -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 .
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
\ No newline at end of file
diff --git a/matlab/set_prior.m b/matlab/set_prior.m
index 5e10d9224..938592f33 100644
--- a/matlab/set_prior.m
+++ b/matlab/set_prior.m
@@ -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 .
- 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);
\ No newline at end of file
+ % 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);
\ No newline at end of file
diff --git a/matlab/uniform_specification.m b/matlab/uniform_specification.m
index 9d09b7bd5..c45f94243 100644
--- a/matlab/uniform_specification.m
+++ b/matlab/uniform_specification.m
@@ -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 .
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
\ No newline at end of file