Блог

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

Формула 1

Формула 2

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

Вот что у меня есть на данный момент. Кажется, я не могу найти свою ошибку. Исправление по умолчанию должно составлять 0,04% при t = 1 (согласно изображению «Результаты для воспроизведения»). Для получения этого результата «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
  

Любая помощь была бы оценена. Спасибо