r - R中具有无限和的修正贝塞尔函数

Question

我试图在 R 中实现以下公式，其中 r0、t、theta0 和 alpha 是常数。另外，I 是第一类修正贝塞尔函数。我想，我的问题是从 Sum 项到公式的末尾。我设置了 n = 150，因为函数很快收敛到零，所以不需要超过 150。我使用的是“Bessel”包。

公式1

公式2

结果重现 第一行 = t，第二行 = Defaultcorr in %

这是我到目前为止所拥有的。我似乎找不到我的错误。当 t = 1 时，Defaultcorr 应为 0.04 %（根据图像“重现的结果”）。要获得这个结果，“m”应该等于 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

任何帮助，将不胜感激。谢谢

标签： rsumbessel-functions