Модифицированная функция Бесселя с бесконечной суммой в R

Question

Я пытаюсь реализовать следующую формулу в R, где r0, t, theta0 и alpha являются константами.Также, я - модифицированная функция Бесселя первого рода.Моя проблема, я полагаю, связана с суммой до конца формулы.Я установил n = 150, учитывая, что функция быстро сходится к нулю, поэтому нет необходимости выходить за пределы 150. Я использую пакет "Bessel".

Формула1

Формула2

Результаты для воспроизведения первая строка = t, вторая строка = Defaultcorrв%

Вот что у меня есть до сих пор.Я не могу найти свою ошибку.Значение по умолчанию должно составлять 0,04% при t = 1 (согласно изображению «Результаты для воспроизведения»).Чтобы получить этот результат, «м» должен быть равен 6,234611709.

V1 = 5
V2 = 5
K1 = 1
K2 = 1
sigma1 = 0.3
sigma2 = 0.3
Z1 = log((V1/K1)/sigma1)
Z2 = log((V2/K2)/sigma2)
t = 1
rho = 0.4

#One firm default -> Firm #1 when lambda = mu
PD_asset1 = 2 * pnorm(-(Z1/sqrt(t)))
PD_asset1
PD_asset2 = 2 * pnorm(-(Z2/sqrt(t)))
PD_asset2

#Results assuming that lambda = mu

#Conditions for alpha, theta0, r0


if (rho < 0) {                        #alpha
  alpha = atan(-(sqrt(1-rho^2)) / rho)
} else {
  alpha = pi + atan(-(sqrt(1-rho^2)) / rho)
}

if (rho > 0) {                        #theta0
  theta0 = atan((Z2 * sqrt(1 - rho^2)) / (Z1 - (rho * Z2)))
} else {
  theta0 = pi + atan((Z2 * sqrt(1 - rho^2)) / (Z1 - (rho * Z2)))
}

r0 = (Z2 / sin(theta0))                 #r0



#Simplified function

h = function(n) {
  (sin((n * pi * theta0)/alpha)/n)
}

n = seq(1, 150, 2)

Bessel1 = (besselI(((r0^2)/(4*t)), (0.5*(((n*pi)/alpha) + 1)), FALSE))
Bessel2 = (besselI(((r0^2)/(4*t)), (0.5*(((n*pi)/alpha) - 1)), FALSE))


l = matrix(data = n, ncol = n)
m = apply((h(l)*(Bessel1 + Bessel2)), 2, FUN = sum) 


PD_asset1_or_asset2 = 1 - (((2 * r0)/(sqrt(2*pi*t))) * (exp(-(r0^2)/(4*t))) * m)
PD_asset1_or_asset2


Var_asset1 = PD_asset1 * (1 - PD_asset1)
Var_asset1

Var_asset2 = PD_asset2 * (1 - PD_asset2)
Var_asset2

PD_asset1_and_asset2 = PD_asset1 + PD_asset2 - PD_asset1_or_asset2
PD_asset1_and_asset2

Defaultcorr = (PD_asset1_and_asset2 - (PD_asset1 * PD_asset2)) / (sqrt(Var_asset1 * Var_asset2))
Defaultcorr

Любая помощь будет принята с благодарностью.Спасибо