diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index c51e0a2ba..7ecbd7b12 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -1,4 +1,4 @@
-function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_)
+function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames)
% This function prints and saves posterior estimates after the mcmc
% (+updates of oo_ & TeX output).
%
@@ -8,6 +8,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% options_ [structure]
% bayestopt_ [structure]
% oo_ [structure]
+% pnames [char] Array of char, names of the prior shapes available
%
% OUTPUTS
% oo_ [structure]
@@ -15,7 +16,7 @@ function oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bay
% SPECIAL REQUIREMENTS
% None.
-% Copyright (C) 2006-2013 Dynare Team
+% Copyright (C) 2006-2015 Dynare Team
%
% This file is part of Dynare.
%
@@ -59,7 +60,6 @@ FirstMhFile = record.KeepedDraws.FirstMhFile;
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
clear record;
-pnames=[' ';'beta ';'gamma';'norm ';'invg ';'unif ';'invg2'];
header_width = row_header_width(M_,estim_params_,bayestopt_);
hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval'];
tit2 = sprintf('%-*s %12s %12s %23s %8s %12s\n',header_width,' ','prior mean','post. mean',hpd_interval,'prior','pstdev');
diff --git a/matlab/distributions/compute_prior_mode.m b/matlab/distributions/compute_prior_mode.m
index e6d35920c..9270fc659 100644
--- a/matlab/distributions/compute_prior_mode.m
+++ b/matlab/distributions/compute_prior_mode.m
@@ -69,7 +69,7 @@ switch shape
case 4
% s = hyperparameters(1)
% nu = hyperparameters(2)
- m = 1/sqrt((hyperparameters(2)+1)/hyperparameters(1));%sqrt((hyperparameters(2)-1)/hyperparameters(1))
+ m = 1/sqrt((hyperparameters(2)+1)/hyperparameters(1));
if length(hyperparameters)>2
m = m + hyperparameters(3);
end
@@ -82,6 +82,18 @@ switch shape
if length(hyperparameters)>2
m = m + hyperparameters(3) ;
end
+ case 8
+ % s = hyperparameters(1) [scale parameter]
+ % k = hyperparameters(2) [shape parameter]
+ if hyperparameters(2)<=1
+ m = 0;
+ else
+ m = hyperparameters(1)*((hyperparameters(2)-1)/hyperparameters(2))^(1/hyperparameters(2))
+ end
+ if length(hyperparameters)>2
+ % Add location parameter
+ m = m + hyperparameters(3) ;
+ end
otherwise
error('Unknown prior shape!')
end
\ No newline at end of file
diff --git a/matlab/distributions/icdfweibull.m b/matlab/distributions/icdfweibull.m
new file mode 100644
index 000000000..9286d504e
--- /dev/null
+++ b/matlab/distributions/icdfweibull.m
@@ -0,0 +1,145 @@
+function t = icdfweibull(proba, scale, shape) % --*-- Unitary tests --*--
+
+% Inverse cumulative distribution function.
+%
+% INPUTS
+% - proba [double] Probability, scalar between 0 and 1.
+% - scale [double] Positive hyperparameter.
+% - shape [double] Positive hyperparameter.
+%
+% OUTPUTS
+% - t [double] scalar such that P(X<=t)=proba
+
+% Copyright (C) 2015 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 .
+
+if proba<2*eps()
+ t = 0;
+ return
+end
+
+if proba>1-2*eps()
+ t = Inf;
+ return
+end
+
+t = exp(log(scale)+log(-log(1-proba))/shape);
+
+%@test:1
+%$ try
+%$ x = icdfweibull(0, 1, 2);
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ t(2) = isequal(x, 0);
+%$ end
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ try
+%$ x = icdfweibull(1, 1, 2);
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ t(2) = isinf(x);
+%$ end
+%$ T = all(t);
+%@eof:2
+
+%@test:3
+%$ scales = [.5, 1, 5];
+%$ shapes = [.1, 1, 2];
+%$ x = NaN(9,1);
+%$
+%$ try
+%$ k = 0;
+%$ for i=1:3
+%$ for j=1:3
+%$ k = k+1;
+%$ x(k) = icdfweibull(.5, scales(i), shapes(j));
+%$ end
+%$ end
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ k = 1;
+%$ for i=1:3
+%$ for j=1:3
+%$ k = k+1;
+%$ t(k) = abs(x(k-1)-scales(i)*log(2)^(1/shapes(j)))<1e-12;
+%$ end
+%$ end
+%$ end
+%$ T = all(t);
+%@eof:3
+
+%@test:4
+%$ debug = false;
+%$ scales = [ .5, 1, 5];
+%$ shapes = [ 1, 2, 3];
+%$ x = NaN(9,1);
+%$ p = 1e-1;
+%$
+%$ try
+%$ k = 0;
+%$ for i=1:3
+%$ for j=1:3
+%$ k = k+1;
+%$ x(k) = icdfweibull(p, scales(i), shapes(j));
+%$ end
+%$ end
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ k = 1;
+%$ for i=1:3
+%$ for j=1:3
+%$ k = k+1;
+%$ shape = shapes(j);
+%$ scale = scales(i);
+%$ density = @(z) exp(lpdfgweibull(z,shape,scale));
+%$ if debug
+%$ [shape, scale, x(k-1)]
+%$ end
+%$ if isoctave
+%$ s = quadv(density, 0, x(k-1));
+%$ else
+%$ s = integral(density, 0, x(k-1));
+%$ end
+%$ if debug
+%$ [s, abs(p-s)]
+%$ end
+%$ t(k) = abs(p-s)<1e-12;
+%$ end
+%$ end
+%$ end
+%$ T = all(t);
+%@eof:4
+
diff --git a/matlab/distributions/lpdfgweibull.m b/matlab/distributions/lpdfgweibull.m
new file mode 100644
index 000000000..cdfc4ef7e
--- /dev/null
+++ b/matlab/distributions/lpdfgweibull.m
@@ -0,0 +1,431 @@
+function [ldens,Dldens,D2ldens] = lpdfgweibull(x,a,b,c) % --*-- Unitary tests --*--
+
+% Evaluates the logged Weibull PDF at x.
+%
+% INPUTS
+% - x [double] m*n matrix of points where the (logged) density will be evaluated,
+% - a [double] m*n matrix of First Weibull distribution parameters,
+% - b [double] m*n matrix of Second Weibull distribution parameters,
+% - c [double] m*n matrix of Third Weibull distribution parameters (location parameter, default is 0).
+%
+% OUTPUTS
+% - ldens [double] m*n matrix of logged (generalized) Weibull densities.
+% - Dldens [double] m*n matrix (first order derivatives w.r.t. x)
+% - D2ldens [double] m*n matrix (second order derivatives w.r.t. x)
+%
+% SPECIAL REQUIREMENTS
+% none
+
+% Copyright (C) 2015 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 .
+
+% Initialize output arguments
+ldens = -Inf(size(x));
+if nargout>1
+ Dldens = NaN(size(x));
+ D2ldens = NaN(size(x));
+end
+
+% Check the number of input arguments
+if nargin<3
+ error('CheckInputArgs','At least three input arguments required!')
+end
+
+% Set default value for location parameter(s).
+if nargin<4
+ c = zeros(size(x));
+end
+
+% Reshape inputs if needed (and possible)
+if ~isscalar(x)
+ if isscalar(a)
+ a = repmat(a, size(x));
+ end
+ if isscalar(b)
+ b = repmat(b, size(x));
+ end
+ if isscalar(c)
+ c = repmat(c, size(x));
+ end
+end
+
+% Get indices for which the densty is defined
+idx = find((x-c)>=0);
+
+% Check size of the inputs
+if (~isequal(size(x), size(a)) || ~isequal(size(x), size(b)) || ~isequal(size(x), size(c)))
+ error('CheckInputArgs','All input arguments must have the same dimensions!')
+end
+
+if isempty(idx), return, end
+
+% Compute the logged density
+
+jdx = find( abs(a-1)<1e-12 & x>=c & (x-c)<1e-12) ;
+ldens(jdx) = 1.0;
+
+if ~isempty(idx)
+ x0 = x(idx)-c(idx);
+ x1 = x0./b(idx);
+ x2 = x1.^a(idx);
+ idx = setdiff(idx, jdx);
+ ldens(idx) = log(a(idx)) - a(idx).*log(b(idx)) + (a(idx)-1).*log(x0) - x2 ;
+end
+
+% Compute the first and second derivatives.
+if nargout>1
+ x3 = (a(idx)-1)./x0;
+ x4 = a(idx).*x2./x1./b(idx);
+ Dldens(idx) = x3 - x4;
+ if nargout>2
+ D2ldens(idx) = -x3./x0 - (a(idx)-1).*x4./x1./b(idx);
+ end
+end
+
+%@test:1
+%$ try
+%$ lpdfgweibull(1.0);
+%$ t(1) = false;
+%$ catch
+%$ t(1) = true;
+%$ end
+%$
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ try
+%$ lpdfgweibull(1.0, .5);
+%$ t(1) = false;
+%$ catch
+%$ t(1) = true;
+%$ end
+%$
+%$ T = all(t);
+%@eof:2
+
+%@test:3
+%$ try
+%$ lpdfgweibull(ones(2,2), .5*ones(2,1), .1*ones(3,1));
+%$ t(1) = false;
+%$ catch
+%$ t(1) = true;
+%$ end
+%$
+%$ T = all(t);
+%@eof:3
+
+%@test:4
+%$ try
+%$ a = lpdfgweibull(-1, .5, .1);
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ t(2) = isinf(a);
+%$ end
+%$
+%$ T = all(t);
+%@eof:4
+
+%@test:5
+%$ try
+%$ a = lpdfgweibull(0, 1, 1);
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ t(2) = abs(a-1.0)<1e-10;
+%$ end
+%$
+%$ T = all(t);
+%@eof:5
+
+%@test:6
+%$ try
+%$ a = lpdfgweibull([0, -1], [1 1], [1 1]);
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$ t(2) = abs(a(1)-1.0)<1e-10;
+%$ t(3) = isinf(a(2));
+%$ end
+%$
+%$ T = all(t);
+%@eof:6
+
+%@test:7
+%$ scale = 1;
+%$ shape = 2;
+%$ mode = scale*((shape-1)/shape)^(1/shape);
+%$
+%$ try
+%$ [a, b, c] = lpdfgweibull(mode, shape, scale);
+%$ p = rand(1000,1)*4;
+%$ am = lpdfgweibull(p, shape*ones(size(p)), scale*ones(size(p)));
+%$ t(1) = true;
+%$ catch
+%$ t(1) = false;
+%$ end
+%$
+%$ if t(1)
+%$
+%$ t(2) = abs(b)<1e-8;
+%$ t(3) = c<0;
+%$ t(4) = all(am.
+
+scale = NaN;
+shape = NaN;
+
+mu2 = mu*mu;
+
+eqn = @(k) gammaln(1+2./k) - 2*gammaln(1+1./k) - log(1+sigma2/mu2);
+eqn2 = @(k) eqn(k).*eqn(k);
+
+kstar = fminbnd(eqn2, 1e-9, 100);
+[shape, fval, exitflag] = fzero(eqn, kstar);
+
+if exitflag<1
+ shape = NaN;
+ return
+end
+
+scale = mu/gamma(1+1/shape);
+
+%@test:1
+%$ debug = false;
+%$ scales = 1:.01:5;
+%$ shapes = .5:.01:2;
+%$ n_scales = length(scales);
+%$ n_shapes = length(shapes);
+%$ mu = NaN(n_scales, n_shapes);
+%$ s2 = NaN(n_scales, n_shapes);
+%$ for i=1:n_shapes
+%$ g1 = gamma(1+1/shapes(i));
+%$ g2 = gamma(1+2/shapes(i));
+%$ g3 = g1*g1;
+%$ for j=1:n_scales
+%$ mu(j, i) = scales(j)*g1;
+%$ s2(j, i) = scales(j)*scales(j)*(g2-g3);
+%$ end
+%$ end
+%$ if debug
+%$ success = [];
+%$ failed1 = [];
+%$ failed1_ = [];
+%$ failed2 = [];
+%$ end
+%$ try
+%$ for i=1:n_shapes
+%$ for j=1:n_scales
+%$ if debug
+%$ disp(sprintf('... mu=%s and s2=%s', num2str(mu(j,i)),num2str(s2(j,i))))
+%$ end
+%$ if ~isnan(mu(j,i)) && ~isnan(s2(j,i)) && ~isinf(mu(j,i)) && ~isinf(s2(j,i))
+%$ [scale, shape] = weibull_specification(mu(j,i), s2(j,i));
+%$ if isnan(scale)
+%$ t = false;
+%$ else
+%$ if abs(scales(j)-scale)<1e-9 && abs(shapes(i)-shape)<1e-9
+%$ t = true;
+%$ else
+%$ t = false;
+%$ end
+%$ end
+%$ if ~t
+%$ failed1 = [failed1; mu(j,i) s2(j,i)];
+%$ failed1_ = [failed1_; shapes(i) scales(j)];
+%$ error('UnitTest','Cannot compute scale and shape hyperparameters for mu=%s and s2=%s', num2str(mu(j,i)), num2str(s2(j,i)))
+%$ end
+%$ if debug
+%$ success = [success; mu(j,i) s2(j,i)];
+%$ end
+%$ else
+%$ failed2 = [failed2; shapes(i) scales(j)];
+%$ continue % Pass this test
+%$ end
+%$ end
+%$ end
+%$ catch
+%$ t = false;
+%$ end
+%$
+%$ if debug
+%$ figure(1)
+%$ plot(success(:,1),success(:,2),'ok');
+%$ if ~isempty(failed1)
+%$ hold on
+%$ plot(failed1(:,1),failed1(:,2),'or');
+%$ hold off
+%$ figure(2)
+%$ plot(failed1_(:,1),failed1_(:,2),'or')
+%$ end
+%$ if ~isempty(failed2)
+%$ figure(2)
+%$ plot(failed2(:,1),failed2(:,2),'or');
+%$ end
+%$ end
+%$ T = all(t);
+%@eof:1
+
diff --git a/matlab/draw_prior_density.m b/matlab/draw_prior_density.m
index 0f7b70450..cf587d61f 100644
--- a/matlab/draw_prior_density.m
+++ b/matlab/draw_prior_density.m
@@ -100,6 +100,12 @@ switch pshape(indx)
end
abscissa = linspace(infbound,supbound,steps);
dens = exp(lpdfig2(abscissa-p3(indx),p6(indx),p7(indx)));
+ case 8
+ density = @(x,a,b,c) exp(lpdfgweibull(x, a, b, c));
+ infbound = p3(indx)+icdfweibull(truncprior,p6(indx),p7(indx));
+ supbound = p3(indx)+icdfweibull(1-truncprior,p6(indx),p7(indx));
+ abscissa = linspace(infbound,supbound,steps);
+ dens = density(abscissa,p6(indx),p7(indx),p3(indx));
otherwise
error(sprintf('draw_prior_density: unknown distribution shape (index %d, type %d)', indx, pshape(indx)));
end
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index b764f0d52..55b40654c 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -130,7 +130,7 @@ ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of
np = estim_params_.np ; % Number of deep parameters.
nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
%% Set the names of the priors.
-pnames = [' ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
+pnames = [' '; 'beta '; 'gamm '; 'norm '; 'invg '; 'unif '; 'invg2'; ' '; 'weibl'];
dr = oo_.dr;
@@ -470,7 +470,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if options_.mh_replic
[marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
% Store posterior statistics by parameter name
- oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_);
+ oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames);
if ~options_.nograph
oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
end
diff --git a/matlab/missing/stats/wblrnd.m b/matlab/missing/stats/wblrnd.m
new file mode 100644
index 000000000..2f55bc8b7
--- /dev/null
+++ b/matlab/missing/stats/wblrnd.m
@@ -0,0 +1,43 @@
+function rnd = wblrnd(a, b)
+% This function produces independent random variates from the Weibull distribution.
+%
+% INPUTS
+% a [double] m*n matrix of positive parameters (scale).
+% b [double] m*n matrix of positive parameters (shape).
+%
+% OUTPUT
+% rnd [double] m*n matrix of independent variates from the beta(a,b) distribution.
+
+% Copyright (C) 2015 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 .
+
+if (nargin ~= 2)
+ error('Two input arguments required!');
+end
+
+if (any(a(:)<0)) || (any(b(:)<0)) || (any(a(:)==Inf)) || (any(b(:)==Inf))
+ error('Input arguments must be finite and positive!');
+end
+
+[ma,na] = size(a);
+[mb,nb] = size(b);
+
+if ma~=mb || na~=nb
+ error('Input arguments must have the same size!');
+end
+
+rnd = a.*(-log(rand(ma, na))).^(1./b);
\ No newline at end of file
diff --git a/matlab/prior_bounds.m b/matlab/prior_bounds.m
index a224c69a0..66baf8494 100644
--- a/matlab/prior_bounds.m
+++ b/matlab/prior_bounds.m
@@ -151,6 +151,14 @@ for i=1:length(p6)
end
end
end
+ case 8
+ if prior_trunc == 0
+ bounds.lb(i) = p3(i);
+ bounds.ub(i) = Inf;
+ else
+ bounds.lb(i) = p3(i)+icdfweibull(prior_trunc,p6(i),p7(i));
+ bounds.ub(i) = p3(i)+icdfweibull(1-prior_trunc,p6(i),p7(i));
+ end
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 7bf27e260..5f89ed0b9 100644
--- a/matlab/prior_draw.m
+++ b/matlab/prior_draw.m
@@ -43,9 +43,8 @@ function pdraw = prior_draw(init,uniform) % --*-- Unitary tests --*--
% along with Dynare. If not, see .
persistent p6 p7 p3 p4 lb ub
-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
-
+persistent uniform_index gaussian_index gamma_index beta_index inverse_gamma_1_index inverse_gamma_2_index weibull_index
+persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws weibull_draws
if nargin>0 && init
p6 = evalin('base', 'bayestopt_.p6');
@@ -97,6 +96,12 @@ if nargin>0 && init
else
inverse_gamma_2_draws = 1;
end
+ weibull_index = find(prior_shape==8);
+ if isempty(weibull_index)
+ weibull_draws = 0;
+ else
+ weibull_draws = 1;
+ end
pdraw = NaN(number_of_estimated_parameters,1);
return
end
@@ -159,6 +164,15 @@ if inverse_gamma_2_draws
end
end
+if weibull_draws
+ pdraw(weibull_index) = weibrnd(p6(weibull_index), p7(weibull_index)) + p3(weibull_index);
+ out_of_bound = find( (pdraw(weibull_index)'>ub(weibull_index)) | (pdraw(weibull_index)'ub(weibull_index)) | (pdraw(weibull_index)'.
-persistent id1 id2 id3 id4 id5 id6
-persistent tt1 tt2 tt3 tt4 tt5 tt6
+persistent id1 id2 id3 id4 id5 id6 id8
+persistent tt1 tt2 tt3 tt4 tt5 tt6 tt8
info=0;
@@ -74,6 +74,12 @@ if nargin > 6 && initialization == 1
if isempty(id6)
tt6 = 0;
end
+ % Weibull indices.
+ tt8 = 1;
+ id8 = find(pshape==8);
+ if isempty(id8)
+ tt8 = 0;
+ end
pflag = 1;
end
@@ -167,6 +173,21 @@ if tt6
end
end
+if tt8
+ logged_prior_density = logged_prior_density + sum(lpdfgweibull(x(id8),p6(id8),p7(id8)));
+ if isinf(logged_prior_density)
+ if nargout ==4
+ info=id8(isinf(log(lpdfgweibull(x(id8),p6(id8),p7(id8)))));
+ end
+ return
+ end
+ if nargout==2
+ [tmp, dlprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8))
+ elseif nargout==3
+ [tmp, dlprior(id8), ds2lprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8))
+ end
+end
+
if nargout==3,
d2lprior = diag(d2lprior);
end
\ No newline at end of file
diff --git a/matlab/set_prior.m b/matlab/set_prior.m
index f94b25962..901e6b5c4 100644
--- a/matlab/set_prior.m
+++ b/matlab/set_prior.m
@@ -260,6 +260,20 @@ for i=1:length(k)
bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ;
end
+% Weibull distribution
+k = find(bayestopt_.pshape == 8);
+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)
+ if (bayestopt_.p1(k(i))bayestopt_.p4(k(i)))
+ error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be above the lower (' num2str(bayestopt_.p3(k(i))) ') bound of the Weibull prior density!']);
+ end
+ [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = weibull_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)), bayestopt_.p2(k(i)));
+ bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 8) ;
+end
+
k = find(isnan(xparam1));
if ~isempty(k)
xparam1(k) = bayestopt_.p1(k);
@@ -273,6 +287,8 @@ if options_.initialize_estimated_parameters_with_the_prior_mode
end
end
+
+
% I create subfolder M_.dname/prior if needed.
CheckPath('prior',M_.dname);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 1e71110ac..00f587964 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,5 +1,6 @@
MODFILES = \
estimation/fs2000.mod \
+ estimation/fs2000_with_weibull_prior.mod \
estimation/fs2000_MCMC_jumping_covariance.mod \
estimation/fs2000_initialize_from_calib.mod \
estimation/fs2000_calibrated_covariance.mod \
diff --git a/tests/estimation/fs2000_with_weibull_prior.mod b/tests/estimation/fs2000_with_weibull_prior.mod
new file mode 100644
index 000000000..54117b221
--- /dev/null
+++ b/tests/estimation/fs2000_with_weibull_prior.mod
@@ -0,0 +1,84 @@
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+ W = l/n;
+ ist = y-c;
+ q = 1 - d;
+
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, weibull_pdf, 0.035449, 25;
+stderr e_m, inv_gamma_pdf, 0.008862, Inf;
+end;
+
+varobs gp_obs gy_obs;
+
+options_.solve_tolf = 1e-12;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8,moments_varendo,consider_only_observed);