diff --git a/matlab/distributions/ig1fun.m b/matlab/distributions/ig1fun.m
new file mode 100644
index 0000000000000000000000000000000000000000..a695bb2887ba6da16c6f5e3b9ef3a1e4a12c8b37
--- /dev/null
+++ b/matlab/distributions/ig1fun.m
@@ -0,0 +1,2 @@
+function err = ig1fun(nu,mu2,sigma2)
+ err = 2*mu2*gamma(nu/2).^2-(sigma2+mu2)*(nu-2).*gamma((nu-1)/2).^2;
\ No newline at end of file
diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m
index 98ba04c145caae1a3c905a47247c9430990471d3..867a42783037771da721ef28dfa95111783d76fd 100644
--- a/matlab/distributions/inverse_gamma_specification.m
+++ b/matlab/distributions/inverse_gamma_specification.m
@@ -1,23 +1,45 @@
-function [s,nu] = inverse_gamma_specification(mu,sigma,type)
+function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag)
+% Computes the inverse Gamma hyperparameters from the prior mean and standard deviation.
-% function [s,nu] = inverse_gamma_specification(mu,sigma,type)
-% Specification of the inverse Gamma function parameters
-% X ~ IG(s,nu)
-%
-% INPUTS
-% mu: expectation
-% sigma: standard deviation
-% type=1: inverse Gamma 1
-% type=2: inverse Gamma 2
-
-% OUTPUTS
-% s: shape parameter
-% nu: scale parameter
-%
-% SPECIAL REQUIREMENTS
-% none
+%@info:
+%! @deftypefn {Function File} {[@var{s}, @var{nu} ]=} colon (@var{mu}, @var{sigma}, @var{type}, @var{use_fzero_flag})
+%! @anchor{distributions/inverse_gamma_specification}
+%! @sp 1
+%! Computes the inverse Gamma (type 1 or 2) hyperparameters from the prior mean (@var{mu}) and standard deviation (@var{sigma}).
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item mu
+%! Double scalar, prior mean.
+%! @item sigma
+%! Positive double scalar, prior standard deviation.
+%! @item type
+%! Integer scalar equal to one or two, type of the Inverse-Gamma distribution.
+%! @item use_fzero_flag
+%! Integer scalar equal to 0 (default) or 1. Use (matlab/octave's implementation of) fzero to solve for @var{nu} if equal to 1, use
+%! dynare's implementation of the secant method otherwise.
+%! @end table
+%! @sp 1
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item s
+%! Positive double scalar (greater than two), first hypermarameter of the Inverse-Gamma prior.
+%! @item nu
+%! Positive double scala, second hypermarameter of the Inverse-Gamma prior.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 1
+%! @ref{set_prior}
+%! @sp 2
+%! @strong{This function calls:}
+%!
+%! @end deftypefn
+%@eod:
-% Copyright (C) 2003-2010 Dynare Team
+% Copyright (C) 2003-2011 Dynare Team
%
% This file is part of Dynare.
%
@@ -34,43 +56,83 @@ function [s,nu] = inverse_gamma_specification(mu,sigma,type)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+check_solution_flag = 1;
+
+if nargin==3
+ use_fzero_flag = 0;
+end
+
sigma2 = sigma^2;
mu2 = mu^2;
-if type == 2; % Inverse Gamma 2
+if type == 2; % Inverse Gamma 2
nu = 2*(2+mu2/sigma2);
s = 2*mu*(1+mu2/sigma2);
-elseif type == 1; % Inverse Gamma 1
- if sigma2 < Inf;
+elseif type == 1; % Inverse Gamma 1
+ if sigma2 < Inf
nu = sqrt(2*(2+mu2/sigma2));
- nu2 = 2*nu;
- nu1 = 2;
- err = 2*mu2*gamma(nu/2)^2-(sigma2+mu2)*(nu-2)*gamma((nu-1)/2)^2;
- while abs(nu2-nu1) > 1e-12
- if err > 0
- nu1 = nu;
- if nu < nu2
- nu = nu2;
+ if use_fzero_flag
+ nu = fzero(@(nu)ig1fun(nu,mu2,sigma2),nu);
+ else
+ nu2 = 2*nu;
+ nu1 = 2;
+ err = ig1fun(nu,mu2,sigma2);
+ err1 = ig1fun(nu1,mu2,sigma2);
+ err2 = ig1fun(nu2,mu2,sigma2);
+ if sign(err2) % Too short interval.
+ while nu2/nu<1000 % Shift the interval contaioning the root.
+ nu1 = nu2;
+ nu2 = nu2*1.01;
+ err2 = ig1fun(nu2,mu2,sigma2);
+ if err2<0
+ break
+ end
+ end
+ end
+ % Sove for nu using the secant method.
+ while abs(nu2-nu1) > 1e-14
+ if err > 0
+ nu1 = nu;
+ if nu < nu2
+ nu = nu2;
+ else
+ nu = 2*nu;
+ nu2 = nu;
+ end
else
- nu = 2*nu;
nu2 = nu;
end
- else
- nu2 = nu;
+ nu = (nu1+nu2)/2;
+ err = ig1fun(nu,mu2,sigma2);
end
- nu = (nu1+nu2)/2;
- err = 2*mu2*gamma(nu/2)^2-(sigma2+mu2)*(nu-2)*gamma((nu-1)/2)^2;
end
s = (sigma2+mu2)*(nu-2);
- else;
+ if check_solution_flag
+ if abs(mu-sqrt(s/2)*gamma((nu-1)/2)/gamma(nu/2))>1e-9
+ error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
+ end
+ if abs(sigma-sqrt(s/(nu-2)-mu^2))>1e-9
+ error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
+ end
+ end
+ else
nu = 2;
s = 2*mu2/pi;
- end;
-else;
+ end
+else
s = -1;
nu = -1;
-end;
+end
-% 01/18/2004 MJ replaced fsolve with secant
-% suppressed chck
-% changed order of output parameters
\ No newline at end of file
+%@test:1
+%$ %addpath ../matlab/distributions
+%$
+%$ % Define two dates
+%$ [s0,nu0] = inverse_gamma_specification(.75,.1,1,0);
+%$ [s1,nu1] = inverse_gamma_specification(.75,.1,1,1);
+%$
+%$ % Check the results.
+%$ t(1) = dyn_assert(s0,s1,1e-6);
+%$ t(2) = dyn_assert(nu0,nu1,1e-6);
+%$ T = all(t);
+%@eof:1
\ No newline at end of file
diff --git a/matlab/set_prior.m b/matlab/set_prior.m
index b5b720b5784b590d8e1cee64849d211f7dadd22d..579dc3a76705e971126477259791a26a2c7fcd19 100644
--- a/matlab/set_prior.m
+++ b/matlab/set_prior.m
@@ -222,7 +222,7 @@ bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
[bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
- inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),1) ;
+ inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),1,0) ;
bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4) ;
end
@@ -244,7 +244,7 @@ bayestopt_.p3(k(k1)) = zeros(length(k1),1);
bayestopt_.p4(k(k2)) = Inf(length(k2),1);
for i=1:length(k)
[bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
- inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),2);
+ inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),2,0);
bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ;
end