Compare commits
7 Commits
master
...
remove-pri
Author | SHA1 | Date |
---|---|---|
Stéphane Adjemian (Ryûk) | a4152295d7 | |
Stéphane Adjemian (Ryûk) | 3229919762 | |
Stéphane Adjemian (Ryûk) | 302d991a25 | |
Stéphane Adjemian (Ryûk) | b448780c34 | |
Stéphane Adjemian (Ryûk) | 394c8e21b8 | |
Stéphane Adjemian (Ryûk) | 7abea05d33 | |
Stéphane Adjemian (Ryûk) | 3e1fb6dd2f |
|
@ -0,0 +1,35 @@
|
||||||
|
function lpd = densities(o, X)
|
||||||
|
|
||||||
|
% Evaluate the logged prior densities at X (for each column).
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
% - X [double] m×n matrix, n points where the prior density is evaluated.
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - lpd [double] 1×n, values of the logged prior density at X.
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
n = columns(X);
|
||||||
|
|
||||||
|
lpd = NaN(1, n);
|
||||||
|
|
||||||
|
parfor i=1:n
|
||||||
|
lpd(i) = density(o, X(:,i));
|
||||||
|
end
|
|
@ -0,0 +1,384 @@
|
||||||
|
function [lpd, dlpd, d2lpd, info] = density(o, x)
|
||||||
|
|
||||||
|
% Evaluate the logged prior density at x.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
% - x [double] m×1 vector, point where the prior density is evaluated.
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - lpd [double] scalar, value of the logged prior density at x.
|
||||||
|
% - dlpd [double] m×1 vector, first order derivatives.
|
||||||
|
% - d2lpd [double] m×1 vector, second order derivatives.
|
||||||
|
%
|
||||||
|
% REMARKS
|
||||||
|
% Second order derivatives holder, d2lpd, has the same rank and shape than dlpd because the priors are
|
||||||
|
% independent (we would have to use a matrix if non orthogonal priors were allowed in Dynare).
|
||||||
|
%
|
||||||
|
% EXAMPLE
|
||||||
|
%
|
||||||
|
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
|
||||||
|
% >> lpd = Prior.dsensity(x)
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
lpd = 0.0;
|
||||||
|
if nargout>1
|
||||||
|
dlpd = zeros(1, length(x));
|
||||||
|
if nargout>2
|
||||||
|
d2lpd = dlpd;
|
||||||
|
if nargout>3
|
||||||
|
info = [];
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isuniform
|
||||||
|
if any(x(o.iduniform)-o.p3(o.iduniform)<0) || any(x(o.iduniform)-o.p4(o.iduniform)>0)
|
||||||
|
lpd = -Inf ;
|
||||||
|
if nargout==4
|
||||||
|
info = o.iduniform((x(o.iduniform)-o.p3(o.iduniform)<0) || (x(o.iduniform)-o.p4(o.iduniform)>0));
|
||||||
|
end
|
||||||
|
return
|
||||||
|
end
|
||||||
|
lpd = lpd - sum(log(o.p4(o.iduniform)-o.p3(o.iduniform))) ;
|
||||||
|
if nargout>1
|
||||||
|
dlpd(o.iduniform) = zeros(length(o.iduniform), 1);
|
||||||
|
if nargout>2
|
||||||
|
d2lpd(o.iduniform) = zeros(length(o.iduniform), 1);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isgaussian
|
||||||
|
switch nargout
|
||||||
|
case 1
|
||||||
|
lpd = lpd + sum(lpdfnorm(x(o.idgaussian), o.p6(o.idgaussian), o.p7(o.idgaussian)));
|
||||||
|
case 2
|
||||||
|
[tmp, dlpd(o.idgaussian)] = lpdfnorm(x(o.idgaussian), o.p6(o.idgaussian), o.p7(o.idgaussian));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
case {3,4}
|
||||||
|
[tmp, dlpd(o.idgaussian), d2lpd(o.idgaussian)] = lpdfnorm(x(o.idgaussian), o.p6(o.idgaussian), o.p7(o.idgaussian));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isgamma
|
||||||
|
switch nargout
|
||||||
|
case 1
|
||||||
|
lpd = lpd + sum(lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma)));
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 2
|
||||||
|
[tmp, dlpd(o.idgamma)] = lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 3
|
||||||
|
[tmp, dlpd(o.idgamma), d2lpd(o.idgamma)] = lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 4
|
||||||
|
[tmp, dlpd(o.idgamma), d2lpd(o.idgamma)] = lpdfgam(x(o.idgamma)-o.p3(o.idgamma), o.p6(o.idgamma), o.p7(o.idgamma));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd)
|
||||||
|
info = o.idgamma(isinf(tmp));
|
||||||
|
return
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isbeta
|
||||||
|
switch nargout
|
||||||
|
case 1
|
||||||
|
lpd = lpd + sum(lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta)));
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 2
|
||||||
|
[tmp, dlpd(o.idbeta)] = lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 3
|
||||||
|
[tmp, dlpd(o.idbeta), d2lpd(o.idbeta)] = lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 4
|
||||||
|
[tmp, dlpd(o.idbeta), d2lpd(o.idbeta)] = lpdfgbeta(x(o.idbeta), o.p6(o.idbeta), o.p7(o.idbeta), o.p3(o.idbeta), o.p4(o.idbeta));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd)
|
||||||
|
info = o.idbeta(isinf(tmp));
|
||||||
|
return
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isinvgamma1
|
||||||
|
switch nargout
|
||||||
|
case 1
|
||||||
|
lpd = lpd + sum(lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1)));
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 2
|
||||||
|
[tmp, dlpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 3
|
||||||
|
[tmp, dlpd(o.idinvgamma1), d2lpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 4
|
||||||
|
[tmp, dlpd(o.idinvgamma1), d2lpd(o.idinvgamma1)] = lpdfig1(x(o.idinvgamma1)-o.p3(o.idinvgamma1), o.p6(o.idinvgamma1), o.p7(o.idinvgamma1));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd)
|
||||||
|
info = o.idinvgamma1(isinf(tmp));
|
||||||
|
return
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isinvgamma2
|
||||||
|
switch nargout
|
||||||
|
case 1
|
||||||
|
lpd = lpd + sum(lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2)));
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 2
|
||||||
|
[tmp, dlpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 3
|
||||||
|
[tmp, dlpd(o.idinvgamma2), d2lpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 4
|
||||||
|
[tmp, dlpd(o.idinvgamma2), d2lpd(o.idinvgamma2)] = lpdfig2(x(o.idinvgamma2)-o.p3(o.idinvgamma2), o.p6(o.idinvgamma2), o.p7(o.idinvgamma2));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd)
|
||||||
|
info = o.idinvgamma2(isinf(tmp));
|
||||||
|
return
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isweibull
|
||||||
|
switch nargout
|
||||||
|
case 1
|
||||||
|
lpd = lpd + sum(lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull)));
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 2
|
||||||
|
[tmp, dlpd(o.idweibull)] = lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 3
|
||||||
|
[tmp, dlpd(o.idweibull), d2lpd(o.idweibull)] = lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd), return, end
|
||||||
|
case 4
|
||||||
|
[tmp, dlpd(o.idweibull), d2lpd(o.idweibull)] = lpdfgweibull(x(o.idweibull), o.p6(o.idweibull), o.p7(o.idweibull));
|
||||||
|
lpd = lpd + sum(tmp);
|
||||||
|
if isinf(lpd)
|
||||||
|
info = o.idweibull(isinf(tmp));
|
||||||
|
return
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
return % --*-- Unit tests --*--
|
||||||
|
|
||||||
|
%@test:1
|
||||||
|
% Fill global structures with required fields...
|
||||||
|
prior_trunc = 1e-10;
|
||||||
|
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
|
||||||
|
p1 = .4*ones(14,1); % Prior mean
|
||||||
|
p2 = .2*ones(14,1); % Prior std.
|
||||||
|
p3 = NaN(14,1);
|
||||||
|
p4 = NaN(14,1);
|
||||||
|
p5 = NaN(14,1);
|
||||||
|
p6 = NaN(14,1);
|
||||||
|
p7 = NaN(14,1);
|
||||||
|
|
||||||
|
for i=1:14
|
||||||
|
switch p0(i)
|
||||||
|
case 1
|
||||||
|
% Beta distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = 1;
|
||||||
|
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
|
||||||
|
case 2
|
||||||
|
% Gamma distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
|
||||||
|
case 3
|
||||||
|
% Normal distribution
|
||||||
|
p3(i) = -Inf;
|
||||||
|
p4(i) = Inf;
|
||||||
|
p6(i) = p1(i);
|
||||||
|
p7(i) = p2(i);
|
||||||
|
p5(i) = p1(i);
|
||||||
|
case 4
|
||||||
|
% Inverse Gamma (type I) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
|
||||||
|
case 5
|
||||||
|
% Uniform distribution
|
||||||
|
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
|
||||||
|
p3(i) = p6(i);
|
||||||
|
p4(i) = p7(i);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
|
||||||
|
case 6
|
||||||
|
% Inverse Gamma (type II) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
|
||||||
|
case 8
|
||||||
|
% Weibull distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
|
||||||
|
otherwise
|
||||||
|
error('This density is not implemented!')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
BayesInfo.pshape = p0;
|
||||||
|
BayesInfo.p1 = p1;
|
||||||
|
BayesInfo.p2 = p2;
|
||||||
|
BayesInfo.p3 = p3;
|
||||||
|
BayesInfo.p4 = p4;
|
||||||
|
BayesInfo.p5 = p5;
|
||||||
|
BayesInfo.p6 = p6;
|
||||||
|
BayesInfo.p7 = p7;
|
||||||
|
|
||||||
|
% Call the tested routine
|
||||||
|
try
|
||||||
|
Prior = dprior(BayesInfo, prior_trunc, false);
|
||||||
|
|
||||||
|
% Compute density at the prior mode
|
||||||
|
lpdstar = Prior.density(p5);
|
||||||
|
|
||||||
|
% Draw random deviates in a loop and evaluate the density.
|
||||||
|
LPD = NaN(10000,1);
|
||||||
|
parfor i = 1:10000
|
||||||
|
x = Prior.draw();
|
||||||
|
LPD(i) = Prior.density(x);
|
||||||
|
end
|
||||||
|
t(1) = true;
|
||||||
|
catch
|
||||||
|
t(1) = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if t(1)
|
||||||
|
t(2) = all(LPD<=lpdstar);
|
||||||
|
end
|
||||||
|
T = all(t);
|
||||||
|
%@eof:1
|
||||||
|
|
||||||
|
%@test:2
|
||||||
|
% Fill global structures with required fields...
|
||||||
|
prior_trunc = 1e-10;
|
||||||
|
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
|
||||||
|
p1 = .4*ones(14,1); % Prior mean
|
||||||
|
p2 = .2*ones(14,1); % Prior std.
|
||||||
|
p3 = NaN(14,1);
|
||||||
|
p4 = NaN(14,1);
|
||||||
|
p5 = NaN(14,1);
|
||||||
|
p6 = NaN(14,1);
|
||||||
|
p7 = NaN(14,1);
|
||||||
|
|
||||||
|
for i=1:14
|
||||||
|
switch p0(i)
|
||||||
|
case 1
|
||||||
|
% Beta distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = 1;
|
||||||
|
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
|
||||||
|
case 2
|
||||||
|
% Gamma distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
|
||||||
|
case 3
|
||||||
|
% Normal distribution
|
||||||
|
p3(i) = -Inf;
|
||||||
|
p4(i) = Inf;
|
||||||
|
p6(i) = p1(i);
|
||||||
|
p7(i) = p2(i);
|
||||||
|
p5(i) = p1(i);
|
||||||
|
case 4
|
||||||
|
% Inverse Gamma (type I) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
|
||||||
|
case 5
|
||||||
|
% Uniform distribution
|
||||||
|
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
|
||||||
|
p3(i) = p6(i);
|
||||||
|
p4(i) = p7(i);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
|
||||||
|
case 6
|
||||||
|
% Inverse Gamma (type II) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
|
||||||
|
case 8
|
||||||
|
% Weibull distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
|
||||||
|
otherwise
|
||||||
|
error('This density is not implemented!')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
BayesInfo.pshape = p0;
|
||||||
|
BayesInfo.p1 = p1;
|
||||||
|
BayesInfo.p2 = p2;
|
||||||
|
BayesInfo.p3 = p3;
|
||||||
|
BayesInfo.p4 = p4;
|
||||||
|
BayesInfo.p5 = p5;
|
||||||
|
BayesInfo.p6 = p6;
|
||||||
|
BayesInfo.p7 = p7;
|
||||||
|
|
||||||
|
% Call the tested routine
|
||||||
|
try
|
||||||
|
Prior = dprior(BayesInfo, prior_trunc, false);
|
||||||
|
mu = NaN(14,1);
|
||||||
|
std = NaN(14,1);
|
||||||
|
|
||||||
|
for i=1:14
|
||||||
|
% Define conditional density (it's also a marginal since priors are orthogonal)
|
||||||
|
f = @(x) exp(Prior.densities(substitute(p5, i, x)));
|
||||||
|
% TODO: Check the version of Octave we use (integral is available as a compatibility wrapper in latest Octave version)
|
||||||
|
m = integral(f, p3(i), p4(i));
|
||||||
|
density = @(x) f(x)/m; % rescaling is required since the probability mass depends on the conditioning.
|
||||||
|
% Compute the conditional expectation
|
||||||
|
mu(i) = integral(@(x) x.*density(x), p3(i), p4(i));
|
||||||
|
std(i) = sqrt(integral(@(x) ((x-mu(i)).^2).*density(x), p3(i), p4(i)));
|
||||||
|
end
|
||||||
|
|
||||||
|
t(1) = true;
|
||||||
|
catch
|
||||||
|
t(1) = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if t(1)
|
||||||
|
t(2) = all(abs(mu-.4)<1e-6);
|
||||||
|
t(3) = all(abs(std-.2)<1e-6);
|
||||||
|
end
|
||||||
|
T = all(t);
|
||||||
|
%@eof:2
|
|
@ -1,26 +1,48 @@
|
||||||
classdef dprior
|
classdef dprior < handle
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
properties
|
properties
|
||||||
p6 = []; % Prior first hyperparameter.
|
p1 = []; % Prior mean.
|
||||||
p7 = []; % Prior second hyperparameter.
|
p2 = []; % Prior stddev.
|
||||||
p3 = []; % Lower bound of the prior support.
|
p3 = []; % Lower bound of the prior support.
|
||||||
p4 = []; % Upper bound of the prior support.
|
p4 = []; % Upper bound of the prior support.
|
||||||
|
p5 = []; % Prior mode.
|
||||||
|
p6 = []; % Prior first hyperparameter.
|
||||||
|
p7 = []; % Prior second hyperparameter.
|
||||||
|
p11 = []; % Prior median
|
||||||
lb = []; % Truncated prior lower bound.
|
lb = []; % Truncated prior lower bound.
|
||||||
ub = []; % Truncated prior upper bound.
|
ub = []; % Truncated prior upper bound.
|
||||||
uniform_index = []; % Index for the uniform priors.
|
name = {}; % Name of the parameter
|
||||||
gaussian_index = []; % Index for the gaussian priors.
|
iduniform = []; % Index for the uniform priors.
|
||||||
gamma_index = []; % Index for the gamma priors.
|
idgaussian = []; % Index for the gaussian priors.
|
||||||
beta_index = []; % Index for the beta priors.
|
idgamma = []; % Index for the gamma priors.
|
||||||
inverse_gamma_1_index = []; % Index for the inverse gamma type 1 priors.
|
idbeta = []; % Index for the beta priors.
|
||||||
inverse_gamma_2_index = []; % Index for the inverse gamma type 2 priors.
|
idinvgamma1 = []; % Index for the inverse gamma type 1 priors.
|
||||||
weibull_index = []; % Index for the weibull priors.
|
idinvgamma2 = []; % Index for the inverse gamma type 2 priors.
|
||||||
uniform_draws = false;
|
idweibull = []; % Index for the weibull priors.
|
||||||
gaussian_draws = false;
|
isuniform = false;
|
||||||
gamma_draws = false;
|
isgaussian = false;
|
||||||
beta_draws = false;
|
isgamma = false;
|
||||||
inverse_gamma_1_draws = false;
|
isbeta = false;
|
||||||
inverse_gamma_2_draws = false;
|
isinvgamma1 = false;
|
||||||
weibull_draws = false;
|
isinvgamma2 = false;
|
||||||
|
isweibull = false;
|
||||||
end
|
end
|
||||||
|
|
||||||
methods
|
methods
|
||||||
|
@ -38,10 +60,17 @@ classdef dprior
|
||||||
%
|
%
|
||||||
% REQUIREMENTS
|
% REQUIREMENTS
|
||||||
% None.
|
% None.
|
||||||
o.p6 = bayestopt_.p6;
|
if ~nargin % Empty object
|
||||||
o.p7 = bayestopt_.p7;
|
return
|
||||||
o.p3 = bayestopt_.p3;
|
end
|
||||||
o.p4 = bayestopt_.p4;
|
if isfield(bayestopt_, 'p1'), o.p1 = bayestopt_.p1; end
|
||||||
|
if isfield(bayestopt_, 'p2'), o.p2 = bayestopt_.p2; end
|
||||||
|
if isfield(bayestopt_, 'p3'), o.p3 = bayestopt_.p3; end
|
||||||
|
if isfield(bayestopt_, 'p4'), o.p4 = bayestopt_.p4; end
|
||||||
|
if isfield(bayestopt_, 'p5'), o.p5 = bayestopt_.p5; end
|
||||||
|
if isfield(bayestopt_, 'p6'), o.p6 = bayestopt_.p6; end
|
||||||
|
if isfield(bayestopt_, 'p7'), o.p7 = bayestopt_.p7; end
|
||||||
|
if isfield(bayestopt_, 'p11'), o.p11 = bayestopt_.p11; end
|
||||||
bounds = prior_bounds(bayestopt_, PriorTrunc);
|
bounds = prior_bounds(bayestopt_, PriorTrunc);
|
||||||
o.lb = bounds.lb;
|
o.lb = bounds.lb;
|
||||||
o.ub = bounds.ub;
|
o.ub = bounds.ub;
|
||||||
|
@ -50,138 +79,38 @@ classdef dprior
|
||||||
else
|
else
|
||||||
prior_shape = bayestopt_.pshape;
|
prior_shape = bayestopt_.pshape;
|
||||||
end
|
end
|
||||||
o.beta_index = find(prior_shape==1);
|
o.idbeta = find(prior_shape==1);
|
||||||
if ~isempty(o.beta_index)
|
if ~isempty(o.idbeta)
|
||||||
o.beta_draws = true;
|
o.isbeta = true;
|
||||||
end
|
end
|
||||||
o.gamma_index = find(prior_shape==2);
|
o.idgamma = find(prior_shape==2);
|
||||||
if ~isempty(o.gamma_index)
|
if ~isempty(o.idgamma)
|
||||||
o.gamma_draws = true;
|
o.isgamma = true;
|
||||||
end
|
end
|
||||||
o.gaussian_index = find(prior_shape==3);
|
o.idgaussian = find(prior_shape==3);
|
||||||
if ~isempty(o.gaussian_index)
|
if ~isempty(o.idgaussian)
|
||||||
o.gaussian_draws = true;
|
o.isgaussian = true;
|
||||||
end
|
end
|
||||||
o.inverse_gamma_1_index = find(prior_shape==4);
|
o.idinvgamma1 = find(prior_shape==4);
|
||||||
if ~isempty(o.inverse_gamma_1_index)
|
if ~isempty(o.idinvgamma1)
|
||||||
o.inverse_gamma_1_draws = true;
|
o.isinvgamma1 = true;
|
||||||
end
|
end
|
||||||
o.uniform_index = find(prior_shape==5);
|
o.iduniform = find(prior_shape==5);
|
||||||
if ~isempty(o.uniform_index)
|
if ~isempty(o.iduniform)
|
||||||
o.uniform_draws = true;
|
o.isuniform = true;
|
||||||
end
|
end
|
||||||
o.inverse_gamma_2_index = find(prior_shape==6);
|
o.idinvgamma2 = find(prior_shape==6);
|
||||||
if ~isempty(o.inverse_gamma_2_index)
|
if ~isempty(o.idinvgamma2)
|
||||||
o.inverse_gamma_2_draws = true;
|
o.isinvgamma2 = true;
|
||||||
end
|
end
|
||||||
o.weibull_index = find(prior_shape==8);
|
o.idweibull = find(prior_shape==8);
|
||||||
if ~isempty(o.weibull_index)
|
if ~isempty(o.idweibull)
|
||||||
o.weibull_draws = true;
|
o.isweibull = true;
|
||||||
end
|
end
|
||||||
end
|
end % dprior (constructor)
|
||||||
|
|
||||||
function p = draw(o)
|
|
||||||
% Return a random draw from the prior distribution.
|
|
||||||
%
|
|
||||||
% INPUTS
|
|
||||||
% - o [dprior]
|
|
||||||
%
|
|
||||||
% OUTPUTS
|
|
||||||
% - p [double] m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
|
|
||||||
%
|
|
||||||
% REMARKS
|
|
||||||
% None.
|
|
||||||
%
|
|
||||||
% EXAMPLE
|
|
||||||
%
|
|
||||||
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
|
|
||||||
% >> d = Prior.draw()
|
|
||||||
p = NaN(rows(o.lb), 1);
|
|
||||||
if o.uniform_draws
|
|
||||||
p(o.uniform_index) = rand(length(o.uniform_index),1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
|
|
||||||
out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.uniform_index) = rand(length(o.uniform_index), 1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index);
|
|
||||||
out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)<o.lb(o.uniform_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
if o.gaussian_draws
|
|
||||||
p(o.gaussian_index) = randn(length(o.gaussian_index), 1).*o.p7(o.gaussian_index) + o.p6(o.gaussian_index);
|
|
||||||
out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.gaussian_index(out_of_bound)) = randn(length(o.gaussian_index(out_of_bound)), 1).*o.p7(o.gaussian_index(out_of_bound)) + o.p6(o.gaussian_index(out_of_bound));
|
|
||||||
out_of_bound = find( (p(o.gaussian_index)>o.ub(o.gaussian_index)) | (p(o.gaussian_index)<o.lb(o.gaussian_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
if o.gamma_draws
|
|
||||||
p(o.gamma_index) = gamrnd(o.p6(o.gamma_index), o.p7(o.gamma_index))+o.p3(o.gamma_index);
|
|
||||||
out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.gamma_index(out_of_bound)) = gamrnd(o.p6(o.gamma_index(out_of_bound)), o.p7(o.gamma_index(out_of_bound)))+o.p3(o.gamma_index(out_of_bound));
|
|
||||||
out_of_bound = find( (p(o.gamma_index)>o.ub(o.gamma_index)) | (p(o.gamma_index)<o.lb(o.gamma_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
if o.beta_draws
|
|
||||||
p(o.beta_index) = (o.p4(o.beta_index)-o.p3(o.beta_index)).*betarnd(o.p6(o.beta_index), o.p7(o.beta_index))+o.p3(o.beta_index);
|
|
||||||
out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.beta_index(out_of_bound)) = (o.p4(o.beta_index(out_of_bound))-o.p3(o.beta_index(out_of_bound))).*betarnd(o.p6(o.beta_index(out_of_bound)), o.p7(o.beta_index(out_of_bound)))+o.p3(o.beta_index(out_of_bound));
|
|
||||||
out_of_bound = find( (p(o.beta_index)>o.ub(o.beta_index)) | (p(o.beta_index)<o.lb(o.beta_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
if o.inverse_gamma_1_draws
|
|
||||||
p(o.inverse_gamma_1_index) = ...
|
|
||||||
sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index)/2, 2./o.p6(o.inverse_gamma_1_index)))+o.p3(o.inverse_gamma_1_index);
|
|
||||||
out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.inverse_gamma_1_index(out_of_bound)) = ...
|
|
||||||
sqrt(1./gamrnd(o.p7(o.inverse_gamma_1_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_1_index(out_of_bound))))+o.p3(o.inverse_gamma_1_index(out_of_bound));
|
|
||||||
out_of_bound = find( (p(o.inverse_gamma_1_index)>o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)<o.lb(o.inverse_gamma_1_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
if o.inverse_gamma_2_draws
|
|
||||||
p(o.inverse_gamma_2_index) = ...
|
|
||||||
1./gamrnd(o.p7(o.inverse_gamma_2_index)/2, 2./o.p6(o.inverse_gamma_2_index))+o.p3(o.inverse_gamma_2_index);
|
|
||||||
out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.inverse_gamma_2_index(out_of_bound)) = ...
|
|
||||||
1./gamrnd(o.p7(o.inverse_gamma_2_index(out_of_bound))/2, 2./o.p6(o.inverse_gamma_2_index(out_of_bound)))+o.p3(o.inverse_gamma_2_index(out_of_bound));
|
|
||||||
out_of_bound = find( (p(o.inverse_gamma_2_index)>o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)<o.lb(o.inverse_gamma_2_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
if o.weibull_draws
|
|
||||||
p(o.weibull_index) = wblrnd(o.p7(o.weibull_index), o.p6(o.weibull_index)) + o.p3(o.weibull_index);
|
|
||||||
out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
|
|
||||||
while ~isempty(out_of_bound)
|
|
||||||
p(o.weibull_index(out_of_bound)) = wblrnd(o.p7(o.weibull_index(out_of_bound)), o.p6(o.weibull_index(out_of_bound)))+o.p3(o.weibull_index(out_of_bound));
|
|
||||||
out_of_bound = find( (p(o.weibull_index)>o.ub(o.weibull_index)) | (p(o.weibull_index)<o.lb(o.weibull_index)));
|
|
||||||
end
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
function P = draws(o, n)
|
|
||||||
% Return n independent random draws from the prior distribution.
|
|
||||||
%
|
|
||||||
% INPUTS
|
|
||||||
% - o [dprior]
|
|
||||||
%
|
|
||||||
% OUTPUTS
|
|
||||||
% - P [double] m×n matrix, random draw from the prior distribution.
|
|
||||||
%
|
|
||||||
% REMARKS
|
|
||||||
% If the Parallel Computing Toolbox is available, the main loop is run in parallel.
|
|
||||||
%
|
|
||||||
% EXAMPLE
|
|
||||||
%
|
|
||||||
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
|
|
||||||
% >> Prior.draws(1e6)
|
|
||||||
P = NaN(rows(o.lb), 1);
|
|
||||||
parfor i=1:n
|
|
||||||
P(:,i) = draw(o);
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
end % methods
|
end % methods
|
||||||
|
|
||||||
end % classdef --*-- Unit tests --*--
|
end % classdef --*-- Unit tests --*--
|
||||||
|
|
||||||
%@test:1
|
%@test:1
|
||||||
|
@ -263,114 +192,22 @@ end % classdef --*-- Unit tests --*--
|
||||||
%$ try
|
%$ try
|
||||||
%$ % Instantiate dprior object
|
%$ % Instantiate dprior object
|
||||||
%$ o = dprior(bayestopt_, prior_trunc, false);
|
%$ o = dprior(bayestopt_, prior_trunc, false);
|
||||||
%$ % Do simulations in a loop and estimate recursively the mean and the variance.
|
|
||||||
%$ for i = 1:ndraws
|
|
||||||
%$ draw = o.draw();
|
|
||||||
%$ m1 = m0 + (draw-m0)/i;
|
|
||||||
%$ m2 = m1*m1';
|
|
||||||
%$ v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
|
|
||||||
%$ m0 = m1;
|
|
||||||
%$ end
|
|
||||||
%$ t(1) = true;
|
%$ t(1) = true;
|
||||||
%$ catch
|
%$ catch
|
||||||
%$ t(1) = false;
|
%$ t(1) = false;
|
||||||
%$ end
|
%$ end
|
||||||
%$
|
%$
|
||||||
%$ if t(1)
|
|
||||||
%$ t(2) = all(abs(m0-bayestopt_.p1)<3e-3);
|
|
||||||
%$ t(3) = all(all(abs(v0-diag(bayestopt_.p2.^2))<5e-3));
|
|
||||||
%$ end
|
|
||||||
%$ T = all(t);
|
%$ T = all(t);
|
||||||
%@eof:1
|
%@eof:1
|
||||||
|
|
||||||
%@test:2
|
%@test:2
|
||||||
%$ % Fill global structures with required fields...
|
|
||||||
%$ prior_trunc = 1e-10;
|
|
||||||
%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
|
|
||||||
%$ p1 = .4*ones(14,1); % Prior mean
|
|
||||||
%$ p2 = .2*ones(14,1); % Prior std.
|
|
||||||
%$ p3 = NaN(14,1);
|
|
||||||
%$ p4 = NaN(14,1);
|
|
||||||
%$ p5 = NaN(14,1);
|
|
||||||
%$ p6 = NaN(14,1);
|
|
||||||
%$ p7 = NaN(14,1);
|
|
||||||
%$
|
|
||||||
%$ for i=1:14
|
|
||||||
%$ switch p0(i)
|
|
||||||
%$ case 1
|
|
||||||
%$ % Beta distribution
|
|
||||||
%$ p3(i) = 0;
|
|
||||||
%$ p4(i) = 1;
|
|
||||||
%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
|
||||||
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
|
|
||||||
%$ case 2
|
|
||||||
%$ % Gamma distribution
|
|
||||||
%$ p3(i) = 0;
|
|
||||||
%$ p4(i) = Inf;
|
|
||||||
%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
|
||||||
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
|
|
||||||
%$ case 3
|
|
||||||
%$ % Normal distribution
|
|
||||||
%$ p3(i) = -Inf;
|
|
||||||
%$ p4(i) = Inf;
|
|
||||||
%$ p6(i) = p1(i);
|
|
||||||
%$ p7(i) = p2(i);
|
|
||||||
%$ p5(i) = p1(i);
|
|
||||||
%$ case 4
|
|
||||||
%$ % Inverse Gamma (type I) distribution
|
|
||||||
%$ p3(i) = 0;
|
|
||||||
%$ p4(i) = Inf;
|
|
||||||
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
|
|
||||||
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
|
|
||||||
%$ case 5
|
|
||||||
%$ % Uniform distribution
|
|
||||||
%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
|
|
||||||
%$ p3(i) = p6(i);
|
|
||||||
%$ p4(i) = p7(i);
|
|
||||||
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
|
|
||||||
%$ case 6
|
|
||||||
%$ % Inverse Gamma (type II) distribution
|
|
||||||
%$ p3(i) = 0;
|
|
||||||
%$ p4(i) = Inf;
|
|
||||||
%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
|
|
||||||
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
|
|
||||||
%$ case 8
|
|
||||||
%$ % Weibull distribution
|
|
||||||
%$ p3(i) = 0;
|
|
||||||
%$ p4(i) = Inf;
|
|
||||||
%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
|
|
||||||
%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
|
|
||||||
%$ otherwise
|
|
||||||
%$ error('This density is not implemented!')
|
|
||||||
%$ end
|
|
||||||
%$ end
|
|
||||||
%$
|
|
||||||
%$ bayestopt_.pshape = p0;
|
|
||||||
%$ bayestopt_.p1 = p1;
|
|
||||||
%$ bayestopt_.p2 = p2;
|
|
||||||
%$ bayestopt_.p3 = p3;
|
|
||||||
%$ bayestopt_.p4 = p4;
|
|
||||||
%$ bayestopt_.p5 = p5;
|
|
||||||
%$ bayestopt_.p6 = p6;
|
|
||||||
%$ bayestopt_.p7 = p7;
|
|
||||||
%$
|
|
||||||
%$ ndraws = 1e5;
|
|
||||||
%$
|
|
||||||
%$ % Call the tested routine
|
|
||||||
%$ try
|
%$ try
|
||||||
%$ % Instantiate dprior object.
|
%$ % Instantiate dprior object
|
||||||
%$ o = dprior(bayestopt_, prior_trunc, false);
|
%$ o = dprior();
|
||||||
%$ X = o.draws(ndraws);
|
|
||||||
%$ m = mean(X, 2);
|
|
||||||
%$ v = var(X, 0, 2);
|
|
||||||
%$ t(1) = true;
|
%$ t(1) = true;
|
||||||
%$ catch
|
%$ catch
|
||||||
%$ t(1) = false;
|
%$ t(1) = false;
|
||||||
%$ end
|
%$ end
|
||||||
%$
|
%$
|
||||||
%$ if t(1)
|
|
||||||
%$ t(2) = all(abs(m-bayestopt_.p1)<3e-3);
|
|
||||||
%$ t(3) = all(all(abs(v-bayestopt_.p2.^2)<5e-3));
|
|
||||||
%$ end
|
|
||||||
%$ T = all(t);
|
%$ T = all(t);
|
||||||
%@eof:2
|
%@eof:2
|
||||||
|
|
|
@ -0,0 +1,197 @@
|
||||||
|
function p = draw(o)
|
||||||
|
|
||||||
|
% Return a random draw from the prior distribution.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - p [double] m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
|
||||||
|
%
|
||||||
|
% REMARKS
|
||||||
|
% None.
|
||||||
|
%
|
||||||
|
% EXAMPLE
|
||||||
|
%
|
||||||
|
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
|
||||||
|
% >> d = Prior.draw()
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
p = NaN(rows(o.lb), 1);
|
||||||
|
if o.isuniform
|
||||||
|
p(o.iduniform) = rand(length(o.iduniform),1).*(o.p4(o.iduniform)-o.p3(o.iduniform)) + o.p3(o.iduniform);
|
||||||
|
oob = find( (p(o.iduniform)>o.ub(o.iduniform)) | (p(o.iduniform)<o.lb(o.iduniform)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.iduniform) = rand(length(o.iduniform), 1).*(o.p4(o.iduniform)-o.p3(o.iduniform)) + o.p3(o.iduniform);
|
||||||
|
oob = find( (p(o.iduniform)>o.ub(o.iduniform)) | (p(o.iduniform)<o.lb(o.iduniform)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isgaussian
|
||||||
|
p(o.idgaussian) = randn(length(o.idgaussian), 1).*o.p7(o.idgaussian) + o.p6(o.idgaussian);
|
||||||
|
oob = find( (p(o.idgaussian)>o.ub(o.idgaussian)) | (p(o.idgaussian)<o.lb(o.idgaussian)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.idgaussian(oob)) = randn(length(o.idgaussian(oob)), 1).*o.p7(o.idgaussian(oob)) + o.p6(o.idgaussian(oob));
|
||||||
|
oob = find( (p(o.idgaussian)>o.ub(o.idgaussian)) | (p(o.idgaussian)<o.lb(o.idgaussian)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isgamma
|
||||||
|
p(o.idgamma) = gamrnd(o.p6(o.idgamma), o.p7(o.idgamma))+o.p3(o.idgamma);
|
||||||
|
oob = find( (p(o.idgamma)>o.ub(o.idgamma)) | (p(o.idgamma)<o.lb(o.idgamma)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.idgamma(oob)) = gamrnd(o.p6(o.idgamma(oob)), o.p7(o.idgamma(oob)))+o.p3(o.idgamma(oob));
|
||||||
|
oob = find( (p(o.idgamma)>o.ub(o.idgamma)) | (p(o.idgamma)<o.lb(o.idgamma)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isbeta
|
||||||
|
p(o.idbeta) = (o.p4(o.idbeta)-o.p3(o.idbeta)).*betarnd(o.p6(o.idbeta), o.p7(o.idbeta))+o.p3(o.idbeta);
|
||||||
|
oob = find( (p(o.idbeta)>o.ub(o.idbeta)) | (p(o.idbeta)<o.lb(o.idbeta)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.idbeta(oob)) = (o.p4(o.idbeta(oob))-o.p3(o.idbeta(oob))).*betarnd(o.p6(o.idbeta(oob)), o.p7(o.idbeta(oob)))+o.p3(o.idbeta(oob));
|
||||||
|
oob = find( (p(o.idbeta)>o.ub(o.idbeta)) | (p(o.idbeta)<o.lb(o.idbeta)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isinvgamma1
|
||||||
|
p(o.idinvgamma1) = ...
|
||||||
|
sqrt(1./gamrnd(o.p7(o.idinvgamma1)/2, 2./o.p6(o.idinvgamma1)))+o.p3(o.idinvgamma1);
|
||||||
|
oob = find( (p(o.idinvgamma1)>o.ub(o.idinvgamma1)) | (p(o.idinvgamma1)<o.lb(o.idinvgamma1)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.idinvgamma1(oob)) = ...
|
||||||
|
sqrt(1./gamrnd(o.p7(o.idinvgamma1(oob))/2, 2./o.p6(o.idinvgamma1(oob))))+o.p3(o.idinvgamma1(oob));
|
||||||
|
oob = find( (p(o.idinvgamma1)>o.ub(o.idinvgamma1)) | (p(o.idinvgamma1)<o.lb(o.idinvgamma1)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isinvgamma2
|
||||||
|
p(o.idinvgamma2) = ...
|
||||||
|
1./gamrnd(o.p7(o.idinvgamma2)/2, 2./o.p6(o.idinvgamma2))+o.p3(o.idinvgamma2);
|
||||||
|
oob = find( (p(o.idinvgamma2)>o.ub(o.idinvgamma2)) | (p(o.idinvgamma2)<o.lb(o.idinvgamma2)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.idinvgamma2(oob)) = ...
|
||||||
|
1./gamrnd(o.p7(o.idinvgamma2(oob))/2, 2./o.p6(o.idinvgamma2(oob)))+o.p3(o.idinvgamma2(oob));
|
||||||
|
oob = find( (p(o.idinvgamma2)>o.ub(o.idinvgamma2)) | (p(o.idinvgamma2)<o.lb(o.idinvgamma2)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isweibull
|
||||||
|
p(o.idweibull) = wblrnd(o.p7(o.idweibull), o.p6(o.idweibull)) + o.p3(o.idweibull);
|
||||||
|
oob = find( (p(o.idweibull)>o.ub(o.idweibull)) | (p(o.idweibull)<o.lb(o.idweibull)));
|
||||||
|
while ~isempty(oob)
|
||||||
|
p(o.idweibull(oob)) = wblrnd(o.p7(o.idweibull(oob)), o.p6(o.idweibull(oob)))+o.p3(o.idweibull(oob));
|
||||||
|
oob = find( (p(o.idweibull)>o.ub(o.idweibull)) | (p(o.idweibull)<o.lb(o.idweibull)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
return % --*-- Unit tests --*--
|
||||||
|
|
||||||
|
%@test:1
|
||||||
|
% Fill global structures with required fields...
|
||||||
|
prior_trunc = 1e-10;
|
||||||
|
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
|
||||||
|
p1 = .4*ones(14,1); % Prior mean
|
||||||
|
p2 = .2*ones(14,1); % Prior std.
|
||||||
|
p3 = NaN(14,1);
|
||||||
|
p4 = NaN(14,1);
|
||||||
|
p5 = NaN(14,1);
|
||||||
|
p6 = NaN(14,1);
|
||||||
|
p7 = NaN(14,1);
|
||||||
|
|
||||||
|
for i=1:14
|
||||||
|
switch p0(i)
|
||||||
|
case 1
|
||||||
|
% Beta distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = 1;
|
||||||
|
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
|
||||||
|
case 2
|
||||||
|
% Gamma distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
|
||||||
|
case 3
|
||||||
|
% Normal distribution
|
||||||
|
p3(i) = -Inf;
|
||||||
|
p4(i) = Inf;
|
||||||
|
p6(i) = p1(i);
|
||||||
|
p7(i) = p2(i);
|
||||||
|
p5(i) = p1(i);
|
||||||
|
case 4
|
||||||
|
% Inverse Gamma (type I) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
|
||||||
|
case 5
|
||||||
|
% Uniform distribution
|
||||||
|
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
|
||||||
|
p3(i) = p6(i);
|
||||||
|
p4(i) = p7(i);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
|
||||||
|
case 6
|
||||||
|
% Inverse Gamma (type II) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
|
||||||
|
case 8
|
||||||
|
% Weibull distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
|
||||||
|
otherwise
|
||||||
|
error('This density is not implemented!')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
BayesInfo.pshape = p0;
|
||||||
|
BayesInfo.p1 = p1;
|
||||||
|
BayesInfo.p2 = p2;
|
||||||
|
BayesInfo.p3 = p3;
|
||||||
|
BayesInfo.p4 = p4;
|
||||||
|
BayesInfo.p5 = p5;
|
||||||
|
BayesInfo.p6 = p6;
|
||||||
|
BayesInfo.p7 = p7;
|
||||||
|
|
||||||
|
ndraws = 1e5;
|
||||||
|
m0 = BayesInfo.p1; %zeros(14,1);
|
||||||
|
v0 = diag(BayesInfo.p2.^2); %zeros(14);
|
||||||
|
|
||||||
|
% Call the tested routine
|
||||||
|
try
|
||||||
|
% Instantiate dprior object
|
||||||
|
o = dprior(BayesInfo, prior_trunc, false);
|
||||||
|
% Do simulations in a loop and estimate recursively the mean and the variance.
|
||||||
|
for i = 1:ndraws
|
||||||
|
draw = o.draw();
|
||||||
|
m1 = m0 + (draw-m0)/i;
|
||||||
|
m2 = m1*m1';
|
||||||
|
v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
|
||||||
|
m0 = m1;
|
||||||
|
end
|
||||||
|
t(1) = true;
|
||||||
|
catch
|
||||||
|
t(1) = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if t(1)
|
||||||
|
t(2) = all(abs(m0-BayesInfo.p1)<3e-3);
|
||||||
|
t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3));
|
||||||
|
end
|
||||||
|
T = all(t);
|
||||||
|
%@eof:1
|
|
@ -0,0 +1,133 @@
|
||||||
|
function P = draws(o, n)
|
||||||
|
|
||||||
|
% Return n independent random draws from the prior distribution.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - P [double] m×n matrix, random draw from the prior distribution.
|
||||||
|
%
|
||||||
|
% REMARKS
|
||||||
|
% If the Parallel Computing Toolbox is available, the main loop is run in parallel.
|
||||||
|
%
|
||||||
|
% EXAMPLE
|
||||||
|
%
|
||||||
|
% >> Prior = dprior(bayestopt_, options_.prior_trunc);
|
||||||
|
% >> Prior.draws(1e6)
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
P = NaN(rows(o.lb), 1);
|
||||||
|
parfor i=1:n
|
||||||
|
P(:,i) = draw(o);
|
||||||
|
end
|
||||||
|
|
||||||
|
return % --*-- Unit tests --*--
|
||||||
|
|
||||||
|
%@test:1
|
||||||
|
% Fill global structures with required fields...
|
||||||
|
prior_trunc = 1e-10;
|
||||||
|
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
|
||||||
|
p1 = .4*ones(14,1); % Prior mean
|
||||||
|
p2 = .2*ones(14,1); % Prior std.
|
||||||
|
p3 = NaN(14,1);
|
||||||
|
p4 = NaN(14,1);
|
||||||
|
p5 = NaN(14,1);
|
||||||
|
p6 = NaN(14,1);
|
||||||
|
p7 = NaN(14,1);
|
||||||
|
|
||||||
|
for i=1:14
|
||||||
|
switch p0(i)
|
||||||
|
case 1
|
||||||
|
% Beta distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = 1;
|
||||||
|
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
|
||||||
|
case 2
|
||||||
|
% Gamma distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
|
||||||
|
case 3
|
||||||
|
% Normal distribution
|
||||||
|
p3(i) = -Inf;
|
||||||
|
p4(i) = Inf;
|
||||||
|
p6(i) = p1(i);
|
||||||
|
p7(i) = p2(i);
|
||||||
|
p5(i) = p1(i);
|
||||||
|
case 4
|
||||||
|
% Inverse Gamma (type I) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
|
||||||
|
case 5
|
||||||
|
% Uniform distribution
|
||||||
|
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
|
||||||
|
p3(i) = p6(i);
|
||||||
|
p4(i) = p7(i);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
|
||||||
|
case 6
|
||||||
|
% Inverse Gamma (type II) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
|
||||||
|
case 8
|
||||||
|
% Weibull distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
|
||||||
|
otherwise
|
||||||
|
error('This density is not implemented!')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
BayesInfo.pshape = p0;
|
||||||
|
BayesInfo.p1 = p1;
|
||||||
|
BayesInfo.p2 = p2;
|
||||||
|
BayesInfo.p3 = p3;
|
||||||
|
BayesInfo.p4 = p4;
|
||||||
|
BayesInfo.p5 = p5;
|
||||||
|
BayesInfo.p6 = p6;
|
||||||
|
BayesInfo.p7 = p7;
|
||||||
|
|
||||||
|
ndraws = 1e5;
|
||||||
|
|
||||||
|
% Call the tested routine
|
||||||
|
try
|
||||||
|
% Instantiate dprior object.
|
||||||
|
o = dprior(BayesInfo, prior_trunc, false);
|
||||||
|
X = o.draws(ndraws);
|
||||||
|
m = mean(X, 2);
|
||||||
|
v = var(X, 0, 2);
|
||||||
|
t(1) = true;
|
||||||
|
catch
|
||||||
|
t(1) = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if t(1)
|
||||||
|
t(2) = all(abs(m-BayesInfo.p1)<3e-3);
|
||||||
|
t(3) = all(all(abs(v-BayesInfo.p2.^2)<5e-3));
|
||||||
|
end
|
||||||
|
T = all(t);
|
||||||
|
%@eof:1
|
|
@ -0,0 +1,42 @@
|
||||||
|
function m = mean(o, resetmoments)
|
||||||
|
|
||||||
|
% Return the prior mean.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
% - resetmoments [logical] Force the computation of the prior mean
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - m [double] n×1 vector, prior mean
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
if nargin<2
|
||||||
|
resetmoments = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if any(isnan(o.p1))
|
||||||
|
resetmoments = true;
|
||||||
|
end
|
||||||
|
|
||||||
|
if resetmoments
|
||||||
|
o.p1 = NaN(size(o.p1));
|
||||||
|
o.moments('mean');
|
||||||
|
end
|
||||||
|
|
||||||
|
m = o.p1;
|
|
@ -0,0 +1,42 @@
|
||||||
|
function m = median(o, resetmoments)
|
||||||
|
|
||||||
|
% Return the prior median.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
% - resetmoments [logical] Force the computation of the prior median
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - m [double] n×1 vector, prior median
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
if nargin<2
|
||||||
|
resetmoments = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if any(isnan(o.p11))
|
||||||
|
resetmoments = true;
|
||||||
|
end
|
||||||
|
|
||||||
|
if resetmoments
|
||||||
|
o.p11 = NaN(size(o.p11));
|
||||||
|
o.moments('median');
|
||||||
|
end
|
||||||
|
|
||||||
|
m = o.p11;
|
|
@ -0,0 +1,42 @@
|
||||||
|
function m = mode(o, resetmoments)
|
||||||
|
|
||||||
|
% Return the prior mode.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
% - resetmoments [logical] Force the computation of the prior mode
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - m [double] n×1 vector, prior mode
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
if nargin<2
|
||||||
|
resetmoments = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if any(isnan(o.p5))
|
||||||
|
resetmoments = true;
|
||||||
|
end
|
||||||
|
|
||||||
|
if resetmoments
|
||||||
|
o.p5 = NaN(size(o.p5));
|
||||||
|
o.moments('mode');
|
||||||
|
end
|
||||||
|
|
||||||
|
m = o.p5;
|
|
@ -0,0 +1,291 @@
|
||||||
|
function o = moments(o, name)
|
||||||
|
|
||||||
|
% Compute the prior moments.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m = o.p1;
|
||||||
|
case 'median'
|
||||||
|
m = o.p11;
|
||||||
|
case 'std'
|
||||||
|
m = o.p2;
|
||||||
|
case 'mode'
|
||||||
|
m = o.p5;
|
||||||
|
otherwise
|
||||||
|
error('%s is not an implemented moemnt.', name)
|
||||||
|
end
|
||||||
|
id = isnan(m);
|
||||||
|
if any(id)
|
||||||
|
% For some parameters the prior mean is not defined.
|
||||||
|
% We compute the first order moment from the
|
||||||
|
% hyperparameters, if the hyperparameters are not
|
||||||
|
% available an error is thrown.
|
||||||
|
if o.isuniform
|
||||||
|
jd = intersect(o.iduniform, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p3(jd))) || any(isnan(o.p4(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (uniform distribution).')
|
||||||
|
end
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p3(jd) + .5*(o.p4(jd)-o.p3(jd));
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p3(jd) + .5*(o.p4(jd)-o.p3(jd));
|
||||||
|
case 'std'
|
||||||
|
m(jd) = (o.p4(jd)-o.p3(jd))/sqrt(12);
|
||||||
|
case 'mode' % Actually we have a continuum of modes with the uniform distribution.
|
||||||
|
m(jd) = o.p3(jd) + .5*(o.p4(jd)-o.p3(jd));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isgaussian
|
||||||
|
jd = intersect(o.idgaussian, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (gaussian distribution).')
|
||||||
|
end
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p6(jd);
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p6(jd);
|
||||||
|
case 'std'
|
||||||
|
m(jd) = o.p7(jd);
|
||||||
|
case 'mode' % Actually we have a continuum of modes with the uniform distribution.
|
||||||
|
m(jd) = o.p6(jd);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isgamma
|
||||||
|
jd = intersect(o.idgamma, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (gamma distribution).')
|
||||||
|
end
|
||||||
|
% α → o.p6, β → o.p7
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p3(jd) + o.p6(jd).*o.p7(jd);
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p3(jd) + gaminv(.5, o.p6(jd), o.p7(jda));
|
||||||
|
case 'std'
|
||||||
|
m(jd) = sqrt(o.p6(jd)).*o.p7(jd);
|
||||||
|
case 'mode'
|
||||||
|
m(jd) = 0;
|
||||||
|
hd = o.p6(jd)>1;
|
||||||
|
m(jd(hd)) = (o.p6(jd(hd))-1).*o.p7(jd(hd));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isbeta
|
||||||
|
jd = intersect(o.idbeta, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd))) || any(isnan(o.p4(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (beta distribution).')
|
||||||
|
end
|
||||||
|
% α → o.p6, β → o.p7
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p3(jd) + (o.p6(jd)./(o.p6(jd)+o.p7(jd))).*(o.p4(jd)-o.p3(jd));
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p3(jd) + betainv(.5, o.p6(jd), o.p7(jd)).*(o.p4(jd)-o.p3(jd));
|
||||||
|
case 'std'
|
||||||
|
m(jd) = (o.p4(jd)-o.p3(jd)).*sqrt(o.p6(jd).*o.p7(jd)./((o.p6(jd)+o.p7(jd)).^2.*(o.p6(jd)+o.p7(jd)+1)));
|
||||||
|
case 'mode'
|
||||||
|
h0 = true(jd, 1);
|
||||||
|
h1 = o.p6(jd)<=1 & o.p7(jd)>1; h0 = h0 & ~h1;
|
||||||
|
h2 = o.p7(jd)<=1 & o.p6(jd)>1; h0 = h0 & ~h2;
|
||||||
|
h3 = o.p6(jd)<1 & o.p7(jd)<1; h0 = h0 & ~h3;
|
||||||
|
h4 = ismembertol(o.p6(jd), 1) & ismembertol(o.p7(jd),1); h0 = h0 & ~h4;
|
||||||
|
m(jd(h1)) = o.p3(jd(h1)); % Standard β has a mode at 0
|
||||||
|
m(jd(h2)) = o.p4(jd(h2)); % Standard β has a mode at 1
|
||||||
|
m(jd(h3)) = o.p3(jd(h3)); % Standard β is bimodal, we pick the lowest mode (0)
|
||||||
|
m(jd(h4)) = o.p3(jd(h4)) + .5*(o.p4(jd(h4))-o.p3(jd(h4))); % Standard β is the uniform distribution (continuum of modes), we pick the mean as the mode
|
||||||
|
m(jd(h0)) = o.p3(jd(h0))+(o.p4(jd(h0))-o.p3(jd(h0))).*((o.p6(jd(h0))-1)./(o.p6(jd(h0))+o.p7(jd(h0))-2)); % β distribution is concave and has a unique interior mode.
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isinvgamma1
|
||||||
|
jd = intersect(o.idinvgamma1, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (inverse gamma type 1 distribution).')
|
||||||
|
end
|
||||||
|
% s → o.p6, ν → o.p7
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p3(jd) + sqrt(.5*o.p6(jd)) .*(gamma(.5*(o.p7(jd)-1))./gamma(.5*o.p7(jd)));
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p3(jd) + 1.0/sqrt(gaminv(.5, o.p7(jd)/2.0, 2.0/o.p6(jd)));
|
||||||
|
case 'std'
|
||||||
|
m(jd) = sqrt( o.p6(jd)./(o.p7(jd)-2)-(.5*o.p6(jd)).*(gamma(.5*(o.p7(jd)-1))./gamma(.5*o.p7(jd))).^2);
|
||||||
|
case 'mode'
|
||||||
|
m(jd) = sqrt((o.p7(jd)-1)./o.p6(jd));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isinvgamma2
|
||||||
|
jd = intersect(o.idinvgamma2, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (inverse gamma type 2 distribution).')
|
||||||
|
end
|
||||||
|
% s → o.p6, ν → o.p7
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p3(jd) + o.p6(jd)./(o.p7(jd)-2);
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p3(jd) + 1.0/gaminv(.5, o.p7(jd)/2.0, 2.0/o.p6(jd));
|
||||||
|
case 'std'
|
||||||
|
m(jd) = sqrt(2./(o.p7(jd)-4)).*o.p6(jd)./(o.p7(jd)-2);
|
||||||
|
case 'mode'
|
||||||
|
m(jd) = o.p6(jd)./(o.p7(jd)+2);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if o.isweibull
|
||||||
|
jd = intersect(o.idweibull, find(id));
|
||||||
|
if ~isempty(jd)
|
||||||
|
if any(isnan(o.p6(jd))) || any(isnan(o.p7(jd))) || any(isnan(o.p3(jd)))
|
||||||
|
error('dprior::mean: Some hyperparameters are missing (weibull distribution).')
|
||||||
|
end
|
||||||
|
% k → o.p6 (shape parameter), λ → o.p7 (scale parameter)
|
||||||
|
% See https://en.wikipedia.org/wiki/Weibull_distribution
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
m(jd) = o.p3(jd) + o.p7(jd).*gamma(1+1./o.p6(jd));
|
||||||
|
case 'median'
|
||||||
|
m(jd) = o.p3(jd) + o.p7(jd).*log(2).^(1./o.p6(jd));
|
||||||
|
case 'std'
|
||||||
|
m(jd) = o.p7(jd).*sqrt(gamma(1+2./o.p6(jd))-gamma(1+1./o.p6(jd)).^2);
|
||||||
|
case 'mode'
|
||||||
|
m(jd) = 0;
|
||||||
|
hd = o.p6(jd)>1;
|
||||||
|
m(jd(hd)) = o.p3(jd(hd)) + o.p7(jd(hd)).*((o.p6(jd(hd))-1)./o.p6(jd(hd))).^(1./o.p6(jd(hd)));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
switch name
|
||||||
|
case 'mean'
|
||||||
|
o.p1 = m;
|
||||||
|
case 'median'
|
||||||
|
o.p11 = m;
|
||||||
|
case 'std'
|
||||||
|
o.p2 = m;
|
||||||
|
case 'mode'
|
||||||
|
o.p5 = m;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
return % --*-- Unit tests --*--
|
||||||
|
|
||||||
|
%@test:5
|
||||||
|
% Fill global structures with required fields...
|
||||||
|
prior_trunc = 1e-10;
|
||||||
|
p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
|
||||||
|
p1 = .4*ones(14,1); % Prior mean
|
||||||
|
p2 = .2*ones(14,1); % Prior std.
|
||||||
|
p3 = NaN(14,1);
|
||||||
|
p4 = NaN(14,1);
|
||||||
|
p5 = NaN(14,1);
|
||||||
|
p6 = NaN(14,1);
|
||||||
|
p7 = NaN(14,1);
|
||||||
|
|
||||||
|
for i=1:14
|
||||||
|
switch p0(i)
|
||||||
|
case 1
|
||||||
|
% Beta distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = 1;
|
||||||
|
[p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
|
||||||
|
case 2
|
||||||
|
% Gamma distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
|
||||||
|
case 3
|
||||||
|
% Normal distribution
|
||||||
|
p3(i) = -Inf;
|
||||||
|
p4(i) = Inf;
|
||||||
|
p6(i) = p1(i);
|
||||||
|
p7(i) = p2(i);
|
||||||
|
p5(i) = p1(i);
|
||||||
|
case 4
|
||||||
|
% Inverse Gamma (type I) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
|
||||||
|
case 5
|
||||||
|
% Uniform distribution
|
||||||
|
[p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
|
||||||
|
p3(i) = p6(i);
|
||||||
|
p4(i) = p7(i);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
|
||||||
|
case 6
|
||||||
|
% Inverse Gamma (type II) distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
|
||||||
|
case 8
|
||||||
|
% Weibull distribution
|
||||||
|
p3(i) = 0;
|
||||||
|
p4(i) = Inf;
|
||||||
|
[p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
|
||||||
|
p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
|
||||||
|
otherwise
|
||||||
|
error('This density is not implemented!')
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
BayesInfo.pshape = p0;
|
||||||
|
BayesInfo.p1 = p1;
|
||||||
|
BayesInfo.p2 = p2;
|
||||||
|
BayesInfo.p3 = p3;
|
||||||
|
BayesInfo.p4 = p4;
|
||||||
|
BayesInfo.p5 = p5;
|
||||||
|
BayesInfo.p6 = p6;
|
||||||
|
BayesInfo.p7 = p7;
|
||||||
|
|
||||||
|
% Call the tested routine
|
||||||
|
try
|
||||||
|
Prior = dprior(BayesInfo, prior_trunc, false);
|
||||||
|
t(1) = true;
|
||||||
|
catch
|
||||||
|
t(1) = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if t(1)
|
||||||
|
t(2) = all(Prior.mean()==.4);
|
||||||
|
t(3) = all(ismembertol(Prior.mean(true),.4));
|
||||||
|
t(4) = all(ismembertol(Prior.variance(),.04));
|
||||||
|
t(5) = all(ismembertol(Prior.variance(true),.04));
|
||||||
|
end
|
||||||
|
T = all(t);
|
||||||
|
%@eof:5
|
|
@ -0,0 +1,41 @@
|
||||||
|
function p = subsref(o, S)
|
||||||
|
|
||||||
|
% Overload subsref method.
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
switch S(1).type
|
||||||
|
case '.'
|
||||||
|
if ismember(S(1).subs, {'p1','p2','p3','p4','p5','p6','p7','lb','ub'})
|
||||||
|
p = builtin('subsref', o, S(1));
|
||||||
|
elseif ismember(S(1).subs, {'draw'})
|
||||||
|
p = feval(S(1).subs, o);
|
||||||
|
elseif ismember(S(1).subs, {'draws', 'density', 'densities', 'moments'})
|
||||||
|
p = feval(S(1).subs, o , S(2).subs{:});
|
||||||
|
elseif ismember(S(1).subs, {'mean', 'median', 'variance', 'mode'})
|
||||||
|
if (length(S)==2 && isempty(S(2).subs)) || length(S)==1
|
||||||
|
p = feval(S(1).subs, o);
|
||||||
|
else
|
||||||
|
p = feval(S(1).subs, o , S(2).subs{:});
|
||||||
|
end
|
||||||
|
else
|
||||||
|
error('dprior::subsref: unknown method (%s).', S(1).subs)
|
||||||
|
end
|
||||||
|
otherwise
|
||||||
|
error('dprior::subsref: %s indexing not implemented.', S(1).type)
|
||||||
|
end
|
|
@ -0,0 +1,42 @@
|
||||||
|
function m = variance(o, resetmoments)
|
||||||
|
|
||||||
|
% Return the prior variance.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - o [dprior]
|
||||||
|
% - resetmoments [logical] Force the computation of the prior variance
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - m [double] n×1 vector, prior variance
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
if nargin<2
|
||||||
|
resetmoments = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
if any(isnan(o.p2))
|
||||||
|
resetmoments = true;
|
||||||
|
end
|
||||||
|
|
||||||
|
if resetmoments
|
||||||
|
o.p2 = NaN(size(o.p2));
|
||||||
|
o.moments('std');
|
||||||
|
end
|
||||||
|
|
||||||
|
m = o.p2.^2;
|
|
@ -1,9 +1,18 @@
|
||||||
function [xparams,lpd,hessian_mat] = ...
|
function [xparams, lpd, hessian_mat] = ...
|
||||||
maximize_prior_density(iparams, prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound,options_,M_,bayestopt_,estim_params_,oo_)
|
maximize_prior_density(iparams, names, options_, M_, Prior, estim_params_, oo_)
|
||||||
|
|
||||||
% Maximizes the logged prior density using Chris Sims' optimization routine.
|
% Maximizes the logged prior density using Chris Sims' optimization routine.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% iparams [double] vector of initial parameters.
|
% - iparams [double] vector of initial parameters.
|
||||||
|
% - Prior [dprior] vector specifying prior densities shapes.
|
||||||
|
% - DynareOptions [struct] Options, AKA options_
|
||||||
|
% - DynareModel [struct] Model description, AKA M_
|
||||||
|
% - EstimatedParams [struct] Info about estimated parameters, AKA estimated_params_
|
||||||
|
% - DynareResults [struct] Results, AKA oo_
|
||||||
|
|
||||||
|
%
|
||||||
|
%
|
||||||
% prior_shape [integer] vector specifying prior densities shapes.
|
% prior_shape [integer] vector specifying prior densities shapes.
|
||||||
% prior_hyperparameter_1 [double] vector, first hyperparameter.
|
% prior_hyperparameter_1 [double] vector, first hyperparameter.
|
||||||
% prior_hyperparameter_2 [double] vector, second hyperparameter.
|
% prior_hyperparameter_2 [double] vector, second hyperparameter.
|
||||||
|
@ -37,10 +46,18 @@ function [xparams,lpd,hessian_mat] = ...
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
[xparams, lpd, exitflag, hessian_mat]=dynare_minimize_objective('minus_logged_prior_density', ...
|
[xparams, lpd, exitflag, hessian_mat] = dynare_minimize_objective('minus_logged_prior_density', ...
|
||||||
iparams, options_.mode_compute, options_, [prior_inf_bound, prior_sup_bound], ...
|
iparams, ...
|
||||||
bayestopt_.name, bayestopt_, [], ...
|
options_.mode_compute, ...
|
||||||
prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound, ...
|
options_, ...
|
||||||
options_,M_,estim_params_,oo_);
|
[Prior.p3, Prior.p4], ...
|
||||||
|
names, ...
|
||||||
|
[], ...
|
||||||
|
[], ...
|
||||||
|
Prior,
|
||||||
|
options_, ...
|
||||||
|
M_, ...
|
||||||
|
estim_params_, ...
|
||||||
|
oo_);
|
||||||
|
|
||||||
lpd = -lpd;
|
lpd = -lpd;
|
||||||
|
|
|
@ -1,24 +1,19 @@
|
||||||
function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,options_,M_,estim_params_,oo_)
|
function [fval, info, exitflag, ~, ~] = minus_logged_prior_density(xparams, Prior, options_, M_, estim_params_, oo_)
|
||||||
% [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparams,pshape,p6,p7,p3,p4,options_,M_,estim_params_,oo_)
|
|
||||||
% Evaluates minus the logged prior density.
|
% Evaluates minus the logged prior density.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% xparams [double] vector of parameters.
|
% - xparams [double] vector of parameters.
|
||||||
% pshape [integer] vector specifying prior densities shapes.
|
% - Prior [dprior] vector specifying prior densities shapes.
|
||||||
% p6 [double] vector, first hyperparameter.
|
% - DynareOptions [struct] Options, AKA options_
|
||||||
% p7 [double] vector, second hyperparameter.
|
% - DynareModel [struct] Model description, AKA M_
|
||||||
% p3 [double] vector, prior's lower bound.
|
% - EstimatedParams [struct] Info about estimated parameters, AKA estimated_params_
|
||||||
% p4 [double] vector, prior's upper bound.
|
% - DynareResults [struct] Results, AKA oo_
|
||||||
% prior_sup_bound [double] vector, prior's upper bound.
|
|
||||||
% options_ [structure] describing the options
|
|
||||||
% M_ [structure] describing the model
|
|
||||||
% estim_params_ [structure] characterizing parameters to be estimated
|
|
||||||
% oo_ [structure] storing the results
|
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% f [double] value of minus the logged prior density.
|
% - fval [double] value of minus the logged prior density.
|
||||||
% info [double] vector: second entry stores penalty, first entry the error code.
|
% - info [double] 4×1 vector, second entry stores penalty, first entry the error code, last entry a penalty (used for optimization).
|
||||||
%
|
|
||||||
% Copyright © 2009-2023 Dynare Team
|
% Copyright © 2009-2023 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
|
@ -36,10 +31,7 @@ function [fval,info,exit_flag,fake_1,fake_2] = minus_logged_prior_density(xparam
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
fake_1 = [];
|
exitflag = true;
|
||||||
fake_2 = [];
|
|
||||||
|
|
||||||
exit_flag = 1;
|
|
||||||
info = zeros(4,1);
|
info = zeros(4,1);
|
||||||
|
|
||||||
%------------------------------------------------------------------------------
|
%------------------------------------------------------------------------------
|
||||||
|
@ -47,74 +39,75 @@ info = zeros(4,1);
|
||||||
%------------------------------------------------------------------------------
|
%------------------------------------------------------------------------------
|
||||||
|
|
||||||
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
|
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
|
||||||
if ~isequal(options_.mode_compute,1) && any(xparams<p3)
|
if ~isequal(options_.mode_compute, 1) && any(xparams<Prior.p3)
|
||||||
k = find(xparams<p3);
|
k = find(xparams<Prior.p3);
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
info(1) = 41;
|
info(1) = 41;
|
||||||
info(4) = sum((p3(k)-xparams(k)).^2);
|
info(4) = sum((Prior.p3(k)-xparams(k)).^2);
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
|
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
|
||||||
if ~isequal(options_.mode_compute,1) && any(xparams>p4)
|
if ~isequal(options_.mode_compute, 1) && any(xparams>Prior.p4)
|
||||||
k = find(xparams>p4);
|
k = find(xparams>Prior.p4);
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
info(1) = 42;
|
info(1) = 42;
|
||||||
info(4) = sum((xparams(k)-p4(k)).^2);
|
info(4) = sum((xparams(k)-Prior.p4(k)).^2);
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
|
% Get the diagonal elements of the covariance matrices for the structural innovations (Q) and the measurement error (H).
|
||||||
M_ = set_all_parameters(xparams,estim_params_,M_);
|
M_ = set_all_parameters(xparams, estim_params_, M_);
|
||||||
|
|
||||||
Q = M_.Sigma_e;
|
Q = M_.Sigma_e;
|
||||||
H = M_.H;
|
H = M_.H;
|
||||||
|
|
||||||
% Test if Q is positive definite.
|
% Test if Q is positive definite.
|
||||||
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_,'calibrated_covariances')
|
if ~issquare(Q) || estim_params_.ncx || isfield(estim_params_, 'calibrated_covariances')
|
||||||
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
|
% Try to compute the cholesky decomposition of Q (possible iff Q is positive definite)
|
||||||
[Q_is_positive_definite, penalty] = ispd(Q);
|
[Q_is_positive_definite, penalty] = ispd(Q);
|
||||||
if ~Q_is_positive_definite
|
if ~Q_is_positive_definite
|
||||||
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
|
% The variance-covariance matrix of the structural innovations is not definite positive. We have to compute the
|
||||||
|
% eigenvalues of this matrix in order to build the endogenous penalty.
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
info(1) = 43;
|
info(1) = 43;
|
||||||
info(4) = penalty;
|
info(4) = penalty;
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
if isfield(estim_params_,'calibrated_covariances')
|
if isfield(estim_params_, 'calibrated_covariances')
|
||||||
correct_flag=check_consistency_covariances(Q);
|
correct_flag = check_consistency_covariances(Q);
|
||||||
if ~correct_flag
|
if ~correct_flag
|
||||||
penalty = sum(Q(estim_params_.calibrated_covariances.position).^2);
|
penalty = sum(Q(estim_params_.calibrated_covariances.position).^2);
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
info(1) = 71;
|
info(1) = 71;
|
||||||
info(4) = penalty;
|
info(4) = penalty;
|
||||||
return4
|
return4
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
% Test if H is positive definite.
|
% Test if H is positive definite.
|
||||||
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_,'calibrated_covariances_ME')
|
if ~issquare(H) || estim_params_.ncn || isfield(estim_params_, 'calibrated_covariances_ME')
|
||||||
[H_is_positive_definite, penalty] = ispd(H);
|
[H_is_positive_definite, penalty] = ispd(H);
|
||||||
if ~H_is_positive_definite
|
if ~H_is_positive_definite
|
||||||
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues of this matrix in order to build the endogenous penalty.
|
% The variance-covariance matrix of the measurement errors is not definite positive. We have to compute the eigenvalues
|
||||||
|
% of this matrix in order to build the endogenous penalty.
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
info(1) = 44;
|
info(1) = 44;
|
||||||
info(4) = penalty;
|
info(4) = penalty;
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
if isfield(estim_params_,'calibrated_covariances_ME')
|
if isfield(estim_params_, 'calibrated_covariances_ME')
|
||||||
correct_flag=check_consistency_covariances(H);
|
correct_flag = check_consistency_covariances(H);
|
||||||
if ~correct_flag
|
if ~correct_flag
|
||||||
penalty = sum(H(estim_params_.calibrated_covariances_ME.position).^2);
|
penalty = sum(H(estim_params_.calibrated_covariances_ME.position).^2);
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
info(1) = 72;
|
info(1) = 72;
|
||||||
info(4) = penalty;
|
info(4) = penalty;
|
||||||
return
|
return
|
||||||
|
@ -127,7 +120,7 @@ end
|
||||||
% 2. Check BK and steady state
|
% 2. Check BK and steady state
|
||||||
%-----------------------------
|
%-----------------------------
|
||||||
|
|
||||||
[~,info] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
|
[~, info] = resol(0, M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
|
||||||
|
|
||||||
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
|
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
|
||||||
if info(1)
|
if info(1)
|
||||||
|
@ -137,14 +130,14 @@ if info(1)
|
||||||
%meaningful second entry of output that can be used
|
%meaningful second entry of output that can be used
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
info(4) = info(2);
|
info(4) = info(2);
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
return
|
return
|
||||||
else
|
else
|
||||||
fval = Inf;
|
fval = Inf;
|
||||||
info(4) = 0.1;
|
info(4) = 0.1;
|
||||||
exit_flag = 0;
|
exitflag = false;
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
fval = - priordens(xparams,pshape,p6,p7,p3,p4);
|
fval = - Prior.density(xparams);
|
||||||
|
|
|
@ -1,5 +1,5 @@
|
||||||
function optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
|
function optimize_prior(options_, M_, oo_, Prior, estim_params_, pnames)
|
||||||
% optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
|
|
||||||
% This routine computes the mode of the prior density using an optimization algorithm.
|
% This routine computes the mode of the prior density using an optimization algorithm.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
|
@ -26,24 +26,25 @@ function optimize_prior(options_, M_, oo_, bayestopt_, estim_params_)
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
oo_.dr = set_state_space(oo_.dr, M_, options_);
|
||||||
|
|
||||||
% Initialize to the prior mean
|
% Initialize to the prior mean
|
||||||
oo_.dr = set_state_space(oo_.dr,M_);
|
xparam1 = Prior.p1;
|
||||||
xparam1 = bayestopt_.p1;
|
|
||||||
|
|
||||||
% Pertubation of the initial condition.
|
% Pertubation of the initial condition.
|
||||||
look_for_admissible_initial_condition = 1; scale = 1.0; iter = 0;
|
look_for_admissible_initial_condition = true; scale = 1.0; iter = 0;
|
||||||
while look_for_admissible_initial_condition
|
while look_for_admissible_initial_condition
|
||||||
xinit = xparam1+scale*randn(size(xparam1));
|
xinit = xparam1+scale*randn(size(xparam1));
|
||||||
if all(xinit(:)>bayestopt_.p3) && all(xinit(:)<bayestopt_.p4)
|
if all(xinit>Prior.p3) && all(xinit<Prior.p4)
|
||||||
M_ = set_all_parameters(xinit,estim_params_,M_);
|
M_ = set_all_parameters(xinit, estim_params_, M_);
|
||||||
[oo_.dr,INFO,M_.params] = resol(0,M_,options_,oo_.dr,oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
|
[dr, INFO, M_, oo_] = resol(0, M_, options_, oo_);
|
||||||
if ~INFO(1)
|
if ~INFO(1)
|
||||||
look_for_admissible_initial_condition = 0;
|
look_for_admissible_initial_condition = false;
|
||||||
end
|
end
|
||||||
else
|
else
|
||||||
if iter == 2000
|
if iter==2000
|
||||||
scale = scale/1.1;
|
scale = scale/1.1;
|
||||||
iter = 0;
|
iter = 0;
|
||||||
else
|
else
|
||||||
iter = iter+1;
|
iter = iter+1;
|
||||||
end
|
end
|
||||||
|
@ -52,23 +53,19 @@ end
|
||||||
|
|
||||||
% Maximization of the prior density
|
% Maximization of the prior density
|
||||||
[xparams, lpd, hessian_mat] = ...
|
[xparams, lpd, hessian_mat] = ...
|
||||||
maximize_prior_density(xinit, bayestopt_.pshape, ...
|
maximize_prior_density(xinit, pnames, options_, M_, Prior, estim_params_, oo_);
|
||||||
bayestopt_.p6, ...
|
|
||||||
bayestopt_.p7, ...
|
|
||||||
bayestopt_.p3, ...
|
|
||||||
bayestopt_.p4,options_,M_,bayestopt_,estim_params_,oo_);
|
|
||||||
|
|
||||||
% Display the results.
|
% Display results.
|
||||||
skipline(2)
|
skipline(2)
|
||||||
disp('------------------')
|
disp('------------------')
|
||||||
disp('PRIOR OPTIMIZATION')
|
disp('PRIOR OPTIMIZATION')
|
||||||
disp('------------------')
|
disp('------------------')
|
||||||
skipline()
|
skipline()
|
||||||
for i = 1:length(xparams)
|
for i = 1:length(xparams)
|
||||||
disp(['deep parameter ' int2str(i) ': ' get_the_name(i,0,M_,estim_params_,options_) '.'])
|
dprintf('deep parameter %u: %s.', i, get_the_name(i, 0, M_, estim_params_, options_))
|
||||||
disp([' Initial condition ....... ' num2str(xinit(i)) '.'])
|
dprintf(' Initial condition ........ %s.', num2str(xinit(i)))
|
||||||
disp([' Prior mode .............. ' num2str(bayestopt_.p5(i)) '.'])
|
dprintf(' Prior mode ............... %s.', num2str(Prior.p5(i)))
|
||||||
disp([' Optimized prior mode .... ' num2str(xparams(i)) '.'])
|
dprintf(' Optimized prior mode ..... %s.', num2str(xparams(i)))
|
||||||
skipline()
|
skipline()
|
||||||
end
|
end
|
||||||
skipline()
|
skipline()
|
||||||
|
|
|
@ -1,51 +1,14 @@
|
||||||
function bounds = prior_bounds(bayestopt, prior_trunc)
|
function bounds = prior_bounds(bayestopt, priortrunc)
|
||||||
|
|
||||||
%@info:
|
|
||||||
%! @deftypefn {Function File} {@var{bounds} =} prior_bounds (@var{bayesopt},@var{option})
|
|
||||||
%! @anchor{prior_bounds}
|
|
||||||
%! @sp 1
|
|
||||||
%! Returns bounds for the prior densities. For each estimated parameter the lower and upper bounds
|
|
||||||
%! are such that the defined intervals contains a probability mass equal to 1-2*@var{option}.prior_trunc. The
|
|
||||||
%! default value for @var{option}.prior_trunc is 1e-10 (set in @ref{global_initialization}).
|
|
||||||
%! @sp 2
|
|
||||||
%! @strong{Inputs}
|
|
||||||
%! @sp 1
|
|
||||||
%! @table @ @var
|
|
||||||
%! @item bayestopt
|
|
||||||
%! Matlab's structure describing the prior distribution (initialized by @code{dynare}).
|
|
||||||
%! @item option
|
|
||||||
%! Matlab's structure describing the options (initialized by @code{dynare}).
|
|
||||||
%! @end table
|
|
||||||
%! @sp 2
|
|
||||||
%! @strong{Outputs}
|
|
||||||
%! @sp 1
|
|
||||||
%! @table @ @var
|
|
||||||
%! @item bounds
|
|
||||||
%! A structure with two fields lb and up (p*1 vectors of doubles, where p is the number of estimated parameters) for the lower and upper bounds.
|
|
||||||
%! @end table
|
|
||||||
%! @sp 2
|
|
||||||
%! @strong{This function is called by:}
|
|
||||||
%! @sp 1
|
|
||||||
%! @ref{get_prior_info}, @ref{dynare_estimation_1}, @ref{dynare_estimation_init}
|
|
||||||
%! @sp 2
|
|
||||||
%! @strong{This function calls:}
|
|
||||||
%! @sp 1
|
|
||||||
%! None.
|
|
||||||
%! @end deftypefn
|
|
||||||
%@eod:
|
|
||||||
|
|
||||||
|
|
||||||
% function bounds = prior_bounds(bayestopt)
|
% function bounds = prior_bounds(bayestopt)
|
||||||
% computes bounds for prior density.
|
% computes bounds for prior density.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% bayestopt [structure] characterizing priors (shape, mean, p1..p4)
|
% - bayestopt [struct] characterizing priors (shape, mean, p1..p4)
|
||||||
|
% - priortrunc [double] scalar, probability mass in the tails to be removed
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
% bounds [double] structure specifying prior bounds (lb and ub fields)
|
% - bounds [struct] prior bounds (lb, lower bounds, and ub, upper bounds, fields are n×1 vectors)
|
||||||
%
|
|
||||||
% SPECIAL REQUIREMENTS
|
|
||||||
% none
|
|
||||||
|
|
||||||
% Copyright © 2003-2023 Dynare Team
|
% Copyright © 2003-2023 Dynare Team
|
||||||
%
|
%
|
||||||
|
@ -64,74 +27,78 @@ function bounds = prior_bounds(bayestopt, prior_trunc)
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
if nargin<2, priortrunc = 0.0; end
|
||||||
|
|
||||||
|
assert(priortrunc>=0 && priortrunc<=1, 'Second input argument must be non negative and not larger than one.')
|
||||||
|
|
||||||
pshape = bayestopt.pshape;
|
pshape = bayestopt.pshape;
|
||||||
p3 = bayestopt.p3;
|
p3 = bayestopt.p3;
|
||||||
p4 = bayestopt.p4;
|
p4 = bayestopt.p4;
|
||||||
p6 = bayestopt.p6;
|
p6 = bayestopt.p6;
|
||||||
p7 = bayestopt.p7;
|
p7 = bayestopt.p7;
|
||||||
|
|
||||||
bounds.lb = zeros(length(p6),1);
|
bounds.lb = zeros(size(p6));
|
||||||
bounds.ub = zeros(length(p6),1);
|
bounds.ub = zeros(size(p6));
|
||||||
|
|
||||||
for i=1:length(p6)
|
for i=1:length(p6)
|
||||||
switch pshape(i)
|
switch pshape(i)
|
||||||
case 1
|
case 1
|
||||||
if prior_trunc == 0
|
if priortrunc==0
|
||||||
bounds.lb(i) = p3(i);
|
bounds.lb(i) = p3(i);
|
||||||
bounds.ub(i) = p4(i);
|
bounds.ub(i) = p4(i);
|
||||||
else
|
else
|
||||||
bounds.lb(i) = betainv(prior_trunc,p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
|
bounds.lb(i) = betainv(priortrunc, p6(i), p7(i))*(p4(i)-p3(i))+p3(i);
|
||||||
bounds.ub(i) = betainv(1-prior_trunc,p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
|
bounds.ub(i) = betainv(1.0-priortrunc, p6(i), p7(i))*(p4(i)-p3(i))+p3(i);
|
||||||
end
|
end
|
||||||
case 2
|
case 2
|
||||||
if prior_trunc == 0
|
if priortrunc==0
|
||||||
bounds.lb(i) = p3(i);
|
bounds.lb(i) = p3(i);
|
||||||
bounds.ub(i) = Inf;
|
bounds.ub(i) = Inf;
|
||||||
else
|
else
|
||||||
bounds.lb(i) = gaminv(prior_trunc,p6(i),p7(i))+p3(i);
|
bounds.lb(i) = gaminv(priortrunc, p6(i), p7(i))+p3(i);
|
||||||
bounds.ub(i) = gaminv(1-prior_trunc,p6(i),p7(i))+p3(i);
|
bounds.ub(i) = gaminv(1.0-priortrunc, p6(i), p7(i))+p3(i);
|
||||||
end
|
end
|
||||||
case 3
|
case 3
|
||||||
if prior_trunc == 0
|
if prior_trunc == 0
|
||||||
bounds.lb(i) = max(-Inf,p3(i));
|
bounds.lb(i) = max(-Inf, p3(i));
|
||||||
bounds.ub(i) = min(Inf,p4(i));
|
bounds.ub(i) = min(Inf, p4(i));
|
||||||
else
|
else
|
||||||
bounds.lb(i) = max(norminv(prior_trunc,p6(i),p7(i)),p3(i));
|
bounds.lb(i) = max(norminv(prior_trunc, p6(i), p7(i)), p3(i));
|
||||||
bounds.ub(i) = min(norminv(1-prior_trunc,p6(i),p7(i)),p4(i));
|
bounds.ub(i) = min(norminv(1-prior_trunc, p6(i), p7(i)), p4(i));
|
||||||
end
|
end
|
||||||
case 4
|
case 4
|
||||||
if prior_trunc == 0
|
if priortrunc==0
|
||||||
bounds.lb(i) = p3(i);
|
bounds.lb(i) = p3(i);
|
||||||
bounds.ub(i) = Inf;
|
bounds.ub(i) = Inf;
|
||||||
else
|
else
|
||||||
bounds.lb(i) = 1/sqrt(gaminv(1-prior_trunc, p7(i)/2, 2/p6(i)))+p3(i);
|
bounds.lb(i) = 1.0/sqrt(gaminv(1.0-priortrunc, p7(i)/2.0, 2.0/p6(i)))+p3(i);
|
||||||
bounds.ub(i) = 1/sqrt(gaminv(prior_trunc, p7(i)/2, 2/p6(i)))+p3(i);
|
bounds.ub(i) = 1.0/sqrt(gaminv(priortrunc, p7(i)/2.0, 2.0/p6(i)))+p3(i);
|
||||||
end
|
end
|
||||||
case 5
|
case 5
|
||||||
if prior_trunc == 0
|
if priortrunc == 0
|
||||||
bounds.lb(i) = p6(i);
|
bounds.lb(i) = p6(i);
|
||||||
bounds.ub(i) = p7(i);
|
bounds.ub(i) = p7(i);
|
||||||
else
|
else
|
||||||
bounds.lb(i) = p6(i)+(p7(i)-p6(i))*prior_trunc;
|
bounds.lb(i) = p6(i)+(p7(i)-p6(i))*priortrunc;
|
||||||
bounds.ub(i) = p7(i)-(p7(i)-p6(i))*prior_trunc;
|
bounds.ub(i) = p7(i)-(p7(i)-p6(i))*priortrunc;
|
||||||
end
|
end
|
||||||
case 6
|
case 6
|
||||||
if prior_trunc == 0
|
if priortrunc == 0
|
||||||
bounds.lb(i) = p3(i);
|
bounds.lb(i) = p3(i);
|
||||||
bounds.ub(i) = Inf;
|
bounds.ub(i) = Inf;
|
||||||
else
|
else
|
||||||
bounds.lb(i) = 1/gaminv(1-prior_trunc, p7(i)/2, 2/p6(i))+p3(i);
|
bounds.lb(i) = 1.0/gaminv(1.0-priortrunc, p7(i)/2.0, 2.0/p6(i))+p3(i);
|
||||||
bounds.ub(i) = 1/gaminv(prior_trunc, p7(i)/2, 2/p6(i))+ p3(i);
|
bounds.ub(i) = 1.0/gaminv(priortrunc, p7(i)/2.0, 2.0/p6(i))+ p3(i);
|
||||||
end
|
end
|
||||||
case 8
|
case 8
|
||||||
if prior_trunc == 0
|
if priortrunc == 0
|
||||||
bounds.lb(i) = p3(i);
|
bounds.lb(i) = p3(i);
|
||||||
bounds.ub(i) = Inf;
|
bounds.ub(i) = Inf;
|
||||||
else
|
else
|
||||||
bounds.lb(i) = p3(i)+wblinv(prior_trunc,p6(i),p7(i));
|
bounds.lb(i) = p3(i)+wblinv(priortrunc, p6(i), p7(i));
|
||||||
bounds.ub(i) = p3(i)+wblinv(1-prior_trunc,p6(i),p7(i));
|
bounds.ub(i) = p3(i)+wblinv(1.0-priortrunc, p6(i), p7(i));
|
||||||
end
|
end
|
||||||
otherwise
|
otherwise
|
||||||
error(sprintf('prior_bounds: unknown distribution shape (index %d, type %d)', i, pshape(i)));
|
error(sprintf('prior_bounds: unknown distribution shape (index %d, type %d)', i, pshape(i)));
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
|
@ -0,0 +1,59 @@
|
||||||
|
function v = substitute(v, i, x)
|
||||||
|
|
||||||
|
% Substitute a scalar in a vector.
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - v [double] m×1 vector
|
||||||
|
% - i [integer] scalar, index for the scalar to be replaced
|
||||||
|
% - x [double] scalar or 1×n vector.
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - v [double] m×1 vector or m×n matrix (with substituted value(s))
|
||||||
|
%
|
||||||
|
% REMARKS
|
||||||
|
% If x is a vector with n elements, then n substitutions are performed (returning n updated vectors in a matrix with n columns)
|
||||||
|
%
|
||||||
|
% EXAMPLES
|
||||||
|
% >> v = ones(2,1);
|
||||||
|
% >> substitude(v, 1, 0)
|
||||||
|
%
|
||||||
|
% ans = %
|
||||||
|
%
|
||||||
|
% 0
|
||||||
|
%
|
||||||
|
% 1
|
||||||
|
%
|
||||||
|
% >> substitute(v, 1, [3 4])
|
||||||
|
%
|
||||||
|
% ans =
|
||||||
|
%
|
||||||
|
% 3 4
|
||||||
|
% 1 1
|
||||||
|
|
||||||
|
% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
assert(isvector(v), 'First input argument must be a vector.')
|
||||||
|
assert(isvector(x), 'Last input argument must be a scalar or a vector.')
|
||||||
|
assert(isscalar(i) && isint(i) && i>0 && i<=length(v), 'Second input argument must be a scalar integer')
|
||||||
|
|
||||||
|
if length(x)==1
|
||||||
|
v(i) = x;
|
||||||
|
else
|
||||||
|
v = repmat(v, 1, length(x));
|
||||||
|
v(i,:) = x;
|
||||||
|
end
|
Loading…
Reference in New Issue