Computing logarithmic derivatives of modified Bessel function of the second kind

Question

I want to know how I can compute the first and second derivatives of logarithm of modified Bessel function of the second kind.

For instance, I'm interested in finding the following derivatives with respect to x:

• ∂/∂x log(K_(x))
• ∂^2/∂x^2 log(K_(x))

The first derivative is equivalent to (∂/∂x(K_(x))) / K_(x). I could use both R functions 'BesselK' (can be found in Bessel and fOptions packages) and 'BesselDK' (fOptions package) to calculate this. Is there a better alternative?

Particularly, how to calculate the second derivative above in R? I couldn't find anything specific anywhere.

Answer 1

Breaking out Python's sympy library (because R doesn't have native symbolic computation facilities):

from sympy import *
nu, x = symbols('nu, x')
b = log(besselk(nu, x))
b1 = b.diff(x)
print(b1)
b2 = b1.diff(x)
print(b2)

Translating the results into R code (for convenience I'm defining a besselk that swaps the order of the arguments to match Python)

besselk <- function(nu, x) besselK(x, nu)
blogd <- function(x, nu) {
   (-besselk(nu - 1, x)/2 - besselk(nu + 1, x)/2)/besselk(nu, x)
}
blogd2 <- function(x, nu) {
   (-besselk(nu - 1, x)/2 - besselk(nu + 1, x)/2)*(besselk(nu - 1, x)/2 + 
     besselk(nu + 1, x)/2)/besselk(nu, x)**2 + (besselk(nu, x)/2 + 
     besselk(nu - 2, x)/4 + besselk(nu + 2, x)/4)/besselk(nu, x)
}

Check results:

blogd(1.5, 2.5)
## [1] -2.051282
blogd2(1.5, 2.5)
## [1] 0.9375411

library(numDeriv)
lb <- function(x, nu) log(besselK(x, nu))
grad(lb, x = 1.5, nu = 2.5)
## [1] -2.051282
drop(hessian(lb, x = 1.5, nu = 2.5))
##  0.9375411

I was lazy and cut-and-pasted the output from sympy. These functions would be made more efficient by vectorizing the calls to besselK (e.g. bvec <- besselK(x, nu = nu + (-2:2))) and plugging the values in to the formula (as it is, besselK is called more times than necessary, especially in the second-derivative calculation; vectorization might not matter much, but not calling besselK more times than necessary would make a big difference).

PS BesselK, BesselDK are in the now-removed-from-CRAN fAsianOptions package