diff --git a/matlab/distributions/ig1fun.m b/matlab/distributions/ig1fun.m index 25aa4624c68cc06f9305e32e6a4f25d83a588c91..684ef30824c1ab8a0cd3914014c910732e0f2fd0 100644 --- a/matlab/distributions/ig1fun.m +++ b/matlab/distributions/ig1fun.m @@ -16,4 +16,4 @@ function err = ig1fun(nu,mu2,sigma2) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see <http://www.gnu.org/licenses/>. -err = 2*mu2*gamma(nu/2).^2-(sigma2+mu2)*(nu-2).*gamma((nu-1)/2).^2; \ No newline at end of file +err = log(2*mu2) - log((sigma2+mu2)*(nu-2)) + 2*( gammaln(nu/2)-gammaln((nu-1)/2) ); diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m index 867a42783037771da721ef28dfa95111783d76fd..ea6a9981cf8383dc9533e33697d420e8ab0684b8 100644 --- a/matlab/distributions/inverse_gamma_specification.m +++ b/matlab/distributions/inverse_gamma_specification.m @@ -35,7 +35,12 @@ function [s,nu] = inverse_gamma_specification(mu,sigma,type,use_fzero_flag) %! @ref{set_prior} %! @sp 2 %! @strong{This function calls:} -%! +%! @sp 2 +%! @strong{Remark:} +%! The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails +%! more often in finding an interval for nu containing a signe change because it expands the interval on both sides and eventually +%! violates the condition nu>2. +%! %! @end deftypefn %@eod: @@ -77,10 +82,9 @@ elseif type == 1; % Inverse Gamma 1 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. + if err2>0 % Too short interval. + while nu2<500 % Shift the interval containing the root. nu1 = nu2; nu2 = nu2*1.01; err2 = ig1fun(nu2,mu2,sigma2); @@ -88,6 +92,9 @@ elseif type == 1; % Inverse Gamma 1 break end end + if err2>0 + error('inverse_gamma_specification:: Failed in finding an interval containing a sign change! You should check that the prior variance is not too small compared to the prior mean...'); + end end % Sove for nu using the secant method. while abs(nu2-nu1) > 1e-14 @@ -125,14 +132,16 @@ else end %@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); -%$ +%$ [s0,nu0] = inverse_gamma_specification(.75,.2,1,0); +%$ [s1,nu1] = inverse_gamma_specification(.75,.2,1,1); +%$ [s3,nu3] = inverse_gamma_specification(.75,.1,1,0); +%$ [s4,nu4] = 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(3) = isnan(s4); +%$ t(4) = isnan(nu4); +%$ t(5) = dyn_assert(s3,16.240907971002265,1e-6);; +%$ t(6) = dyn_assert(nu3,30.368398202624046,1e-6);; %$ T = all(t); %@eof:1 \ No newline at end of file