[R-SIG-Finance] bessel functions

[Diethelm Wuertz]

You can just  use for nu > 10 and x > 1000

x = 1000
nu = 10
mu = 4*nu^2

    A1 =  1
    A2 = A1 * (mu-  1) / (1 * (8*x))
    A3 = A2 * (mu-  9) / (2 * (8*x))
    A4 = A3 * (mu- 25) / (3 * (8*x))
    A5 = A4 * (mu- 49) / (4 * (8*x))
    A6 = A5 * (mu- 81) / (5 * (8*x))
    A7 = A6 * (mu-121) / (6 * (8*x))
   
    print(1/sqrt(2*pi*x) * (A1 - A2 + A3 - A4 + A5 - A6 + A7), digits = 12)
    # [1] 0.0120015950241
    print(besselI(x, nu, TRUE), digits = 12)
    # [1] 0.0120015950241
   
 
x = 5539.947191
nu = 11
mu = 4*nu^2

    A1 =  1
    A2 = A1 * (mu-  1) / (1 * (8*x))
    A3 = A2 * (mu-  9) / (2 * (8*x))
    A4 = A3 * (mu- 25) / (3 * (8*x))
    A5 = A4 * (mu- 49) / (4 * (8*x))
    A6 = A5 * (mu- 81) / (5 * (8*x))
    A7 = A6 * (mu-121) / (6 * (8*x))
   
    print(1/sqrt(2*pi*x) * (A1 - A2 + A3 - A4 + A5 - A6 + A7), digits = 12)
    # 0.00530180603357
    # BesselI fails
    # Mathematica:
    # 0.00530180603357
    # Maple:
    # 0.005301806035

regards
Diethelm Wuertz




Martin Maechler wrote:

>Hi Stefano,  et al.
>
>I have been too busy to get this earlier,
>
>  
>
>>>>>>"Stefano" == stefano iacus <stefano.iacus at unimi.it>
>>>>>>    on Thu, 14 Sep 2006 19:28:15 +0900 writes:
>>>>>>            
>>>>>>
>
>    Stefano> Hi, does anyone know about a reliable algorithm to
>    Stefano> evaluate the exponentially rescaled bessel function
>    Stefano> of first find, any order, real argument?
>
>    Stefano> what I need is the equivalent of R function
>
>    Stefano> besselI(x, nu, expon.scaled = TRUE)
>
>    Stefano> for large real x (say 5000) and large nu (say 10 to
>    Stefano> 30)
>
>    Stefano> The R implementation does not satisfy my needs
>    Stefano> because, for example,
>
>
>    >> besselI(5539.947191,11,exp=TRUE)
>    Stefano> [1] NaN
>
>    Stefano> but if you ask Mathematica you get this number:
>    Stefano> 0.00530180603357
>
>    Stefano> To be more specific I need this to evaluate the
>    Stefano> likelihood of the CIR model. I know there are many
>    Stefano> different estimators for the CIR but I need to
>    Stefano> evaluate the likelihood.
>
>    Stefano> Googling around did not help me.
>
>But  RSiteSearch("bessel")  does immediately help
>
>The answer is *not* to use numerical recipes 
>(they have ``source'' but have heavily restricted copyrights;
> and also very often, their code is far suboptimal) as 'mel'
>indicated, but  for "special math functions" like these,
>rather look at the GSL = GNU Scientific Library.
>
>It has all kinds of Bessel functions and uses newer versions of
>the algorithm than R, since
>--- as we have the 'gsl' package on CRAN - thanks to Robin Hankin ---
>you get the desired result from
>
> > library(gsl)
> > bessel_In_scaled(x= 5539.947191, n= 11)
> [1] 0.005301806
>
>
>Martin Maechler
>
>_______________________________________________
>R-SIG-Finance at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>
>  
>