Итеративный метод оценки американского опциона "пут"

Question

Я пытаюсь реализовать алгоритм определения цены американского опциона с использованием итеративного метода, как в статье «Простой итеративный метод оценки американских опционов»

Математика может быть сложной, но мы этого не делаем надо много беспокоиться об этом. Основная идея 1046 * заключается в том, что мы делим время до зрелости на n интервалов (n + 1 узел), в каждом узле у нас есть одна оптимальная граница упражнения, и мы пытаемся подобрать кривую, которая проходит через все эти точки (n- степень полинома). Название «итеративный метод» предполагает, что мы проводим это приближение несколько раз, каждый раз, когда у нас будет лучшая кривая, где она будет ближе к правильной кривой.

Достаточно сказано со следующими параметрами: S, T , К, п ,, г, сигма, к; алгоритм следующий:

1. Разделите время до погашения на n интервалов, что создает n + 1 границы упражнения. t0, t1, t2, ... ti, .., tn
2. При t = 0 (срок погашения) оптимальной границей упражнения является K (цена страйка).
3. Начальное предположение для B (ti) = K
4. Примените формулу (6), чтобы найти первый раунд оценки для B (i).

введите описание изображения здесь

Приблизительная полиномиальная функция для n + 1 точек данных B (ti)

Используйте уравнение (5) до go для следующего цикла итерации

Продолжайте повторять шаг 5. k раз или до схождения

введите описание изображения здесь

Ниже приведены мои коды для этих 2 шагов в Python:

dt=T/n
B = [K]*(n+1)
#first iteration - equation (6) in paper
for i in range(n):
    i+=1
    tt = i*dt
    A1 = K*(NCDF(d12(K,tt,K)[0])+(1/sig/np.sqrt(2*np.pi*tt))*np.exp(-0.5*d12(K,tt,K)[0]**2))**(-1)
    A2 = 1/sig/np.sqrt(2*np.pi*tt)*np.exp(-(r*tt+0.5*d12(K,tt,K)[1]**2))
    A3 = (2*sig*r)/(2*r+sig**2)*(2*NCDF((2*r+sig**2)*np.sqrt(tt)/(2*sig))-1)
    B[i]=A1*(A2+A3)

for i in range(k): #kth iteration as in equation (5)
    xvals=np.linspace(0,T,n+1)
    params=np.polyfit(xvals,B,deg=n)
    for j in range(n): #equation (5) at each node point
        j+=1
        tt = j*dt
        A1 = (NCDF(d12(B[j],tt,K)[0])+(1/sig/np.sqrt(2*np.pi*tt))*np.exp(-0.5*d12(B[j],tt,K)[0]**2))**(-1)
        A2 = (1/sig/np.sqrt(2*np.pi*tt)*K*np.exp(-(r*tt+0.5*d12(B[j],tt,K)[1]**2)))
        A3 = r*K*sci.fixed_quad(lambda x: (1/sig/np.sqrt(2*np.pi*(tt-x)))*r*np.exp(-(r*(tt-x)+\
        0.5*d12(B[j],tt,Polyf(params,x))[1]**2)),0,tt)[0]
        B[j]=A1*(A2+A3)

xvals=np.linspace(0,T,n+1)
params=np.polyfit(xvals,B,deg=n)

Price = BSPut(S0,K,r,sig,T) + sci.fixed_quad(lambda x: r*K*np.exp(-r*(T-x))*NCDF(-d12(S0,T-x,Polyf(params,x))[1]),0,T)[0]
return Price

Не беспокойтесь о функциях Polyf и BSPut, так как я тщательно их проверил. Проблема в том, что я никогда не смог бы воспроизвести правильный результат с этими кодами. Даже кривая B (ti) после первого раунда итерации не соответствует результату в статье при S = ​​45, K = 45, r = 0,05, сигма = 0,2, T = 0,5, n = 16. Если вы Запустите мои коды, B1 (ti) намного ниже, чем B1 (ti) в статье (я думаю, что минимальное значение составляет приблизительно 41, в то время как моя программа дала мне B (16) = 37 после первой итерации.

Пример результата

Я был бы очень признателен, если бы кто-нибудь показал мне, что может быть не так в моей реализации.

Большое спасибо!