stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'cmr-prior'

time-shift
Stéphane Adjemian (Charybdis) 2015-12-04 15:06:28 +01:00
parent 954e57c421 a0119404c5
commit 5ef86416bf
14 changed files with 914 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
matlab/GetPosteriorParametersStatistics.m
View File

@ -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');

14
matlab/distributions/compute_prior_mode.m
View File

@ -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

145
matlab/distributions/icdfweibull.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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

431
matlab/distributions/lpdfgweibull.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
% 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<a);
%$ end
%$
%$ T = all(t);
%@eof:7
%@test:8
%$ scale = 1;
%$ shape = 2;
%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
%$
%$ try
%$ if isoctave
%$ s = quadv(density, .0000000001, 100000);
%$ else
%$ s = integral(density, 0, 100000);
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = abs(s-1)<1e-6;
%$ end
%$
%$ T = all(t);
%@eof:8
%@test:9
%$ scale = 1;
%$ shape = 1;
%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
%$
%$ try
%$ if isoctave
%$ s = quadv(density, .0000000001, 100000);
%$ else
%$ s = integral(density, 0, 100000);
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = abs(s-1)<1e-6;
%$ end
%$
%$ T = all(t);
%@eof:9
%@test:10
%$ scale = 1;
%$ shape = .5;
%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
%$
%$ try
%$ if isoctave
%$ s = quadv(density, .0000000001, 100000);
%$ else
%$ s = integral(density, 0, 100000);
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = abs(s-1)<1e-6;
%$ end
%$
%$ T = all(t);
%@eof:10
%@test:11
%$ scale = 1;
%$ shape = 2;
%$ xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
%$
%$ try
%$ if isoctave
%$ s = quadv(xdens, .0000000001, 100000);
%$ else
%$ s = integral(xdens, 0, 100000);
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
%$ end
%$
%$ T = all(t);
%@eof:11
%@test:12
%$ scale = 1;
%$ shape = 1;
%$ xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
%$
%$ try
%$ if isoctave
%$ s = quadv(xdens, .0000000001, 100000);
%$ else
%$ s = integral(xdens, 0, 100000);
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
%$ end
%$
%$ T = all(t);
%@eof:12
%@test:13
%$ scale = 1;
%$ shape = .5;
%$ xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
%$
%$ try
%$ if isoctave
%$ s = quadv(xdens, .0000000001, 100000);
%$ else
%$ s = integral(xdens, 0, 100000);
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
%$ end
%$
%$ T = all(t);
%@eof:13
%@test:14
%$ scale = 1;
%$ shape = 2;
%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
%$ n = 200;
%$
%$ try
%$ s = NaN(n, 1);
%$ for i=1:n
%$ if isoctave()
%$ s(i) = quadv(density, .0000000001, .1*i);
%$ else
%$ s(i) = integral(density, 0, .1*i);
%$ end
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ for i=1:n
%$ x = .1*i;
%$ q = 1-exp(-(x/scale)^shape);
%$ t(i+1) = abs(s(i)-q)<1e-6;
%$ end
%$ end
%$
%$ T = all(t);
%@eof:14
%@test:15
%$ scale = 1;
%$ shape = 1;
%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
%$ n = 200;
%$
%$ try
%$ s = NaN(n, 1);
%$ for i=1:n
%$ if isoctave()
%$ s(i) = quadv(density, .0000000001, .1*i);
%$ else
%$ s(i) = integral(density, 0, .1*i);
%$ end
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ for i=1:n
%$ x = .1*i;
%$ q = 1-exp(-(x/scale)^shape);
%$ t(i+1) = abs(s(i)-q)<1e-6;
%$ end
%$ end
%$
%$ T = all(t);
%@eof:15
%@test:16
%$ scale = 1;
%$ shape = .5;
%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
%$ n = 200;
%$
%$ try
%$ s = NaN(n, 1);
%$ for i=1:n
%$ if isoctave()
%$ s(i) = quadv(density, .0000000001, .1*i);
%$ else
%$ s(i) = integral(density, 0, .1*i);
%$ end
%$ end
%$ t(1) = true;
%$ catch
%$ t(1) = false;
%$ end
%$
%$ if t(1)
%$ for i=1:n
%$ x = .1*i;
%$ q = 1-exp(-(x/scale)^shape);
%$ t(i+1) = abs(s(i)-q)<1e-6;
%$ end
%$ end
%$
%$ T = all(t);
%@eof:16

122
matlab/distributions/weibull_specification.m Normal file
View File

@ -0,0 +1,122 @@
function [scale, shape] = weibull_specification(mu, sigma2) % --*-- Unitary tests --*--
% Returns the hyperparameters of the Weibull distribution given the expectation and variance.
%
% INPUTS
%
%
% OUTPUTS
%
%
% 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 <http://www.gnu.org/licenses/>.
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

6
matlab/draw_prior_density.m
View File

@ -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

4
matlab/dynare_estimation_1.m
View File

@ -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

43
matlab/missing/stats/wblrnd.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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);

8
matlab/prior_bounds.m
View File

@ -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

20
matlab/prior_draw.m
View File

@ -43,9 +43,8 @@ function pdraw = prior_draw(init,uniform) % --*-- Unitary tests --*--
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
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)'<lb(weibull_index)));
while ~isempty(out_of_bound),
pdraw(weibull_index(out_of_bound)) = weibrnd(p6(weibull_index(out_of_bound)),p7(weibull_index(out_of_bound)))+p3(weibull_index(out_of_bound));
out_of_bound = find( (pdraw(weibull_index)'>ub(weibull_index)) | (pdraw(weibull_index)'<lb(weibull_index)));
end
end
%@test:1
%$ %% Initialize required structures
%$ options_.prior_trunc=0;

25
matlab/priordens.m
View File

@ -32,8 +32,8 @@ function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
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

16
matlab/set_prior.m
View File

@ -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_.p3(k(i))) || (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);

1
tests/Makefile.am
View File

@ -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 \

84
tests/estimation/fs2000_with_weibull_prior.mod Normal file
View File

@ -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);