2020HW-Options3.ipynb

Kernel: Python 3 (system-wide)

Задача 1

Реализовать функции для расчета цены Европейкого или Американкого опциона используя биномиальную модель.

def gen_lattice(S0, u, d, N):
    """generates a binomial lattice
    
    for a given up, down, start value and number of steps (N).
    Resulting lattice has N+1 levels.
    """
    S = [float(S0)]

    for i in range(1, N+1):
        for j in range(0, i+1):
            S.append( S0 * d**j * u**(i-j) )

    return S

#def lattice_levels(S):
#    return int( round( (m.sqrt(8*len(S)+1)-1)/2 ) )

import numpy as np

import sys
print (sys.version)
#if sys.version_info >= (3,4):
#    print( "with enums" )

from enum import Enum
class CallPut(Enum):
    call = 1
    put = 2

class ExerciseStyle(Enum):
    euro = 1
    amer = 2

3.6.9 (default, Nov  7 2019, 10:44:02) 
[GCC 8.3.0]

def pv_crr(amerEuro, callPut, S0, K, T, r, sigma, N):
    import math as m
    dt = T / N
    df = m.exp(-r * dt)

    u = m.exp(sigma * m.sqrt(dt))
    d = 1 / u
    p = ( m.exp(r * dt) - d ) / (u - d)


    S = gen_lattice(S0=S0, N=N, u=u, d=d)

    L = N+1 #lattice_levels(S)

    payoff = lambda x: max( 0, x - K ) if CallPut.call == callPut else max( 0, K - x)

    for i in range(1,N+2):

        S[int((N+1)*(N+2)/2) - i] = payoff(S[int((N+1)*(N+2)/2) - i])

    # Calculate payoff at the last lattice level
    ## TODO ##h

    for j in reversed(range(0,N+1)):
        for i in range(0,j):

            if ExerciseStyle.amer == amerEuro: 
                S[int(j*(j-1)/2) + i] = np.maximum((p*S[int((j+1)*j/2) + i] + (1-p)*S[int((j+1)*j/2) + i+1])*np.exp(-r*dt), payoff(S[int(j*(j-1)/2) + i]))

            else:
                S[int(j*(j-1)/2) + i] = (p*S[int((j+1)*j/2) + i] + (1-p)*S[int((j+1)*j/2) + i+1])*np.exp(-r*dt)

    # Go backwards, calculate extected value for prev node, based on known nodes
    # Calculate payoff at node k, based on expected value of S
    ## TODO ##

    return S[0], S

# parameters
S0 = 100.
T = 1.
r = 0.05
sigma = 0.20
K = 100.
N = 1000

es = ExerciseStyle.euro
pvC, _ = pv_crr(es, CallPut.call, S0=S0, K=K, T=T, r=r, sigma=sigma, N=N)
print( pvC )
pvP, _ = pv_crr(es, CallPut.put, S0=S0, K=K, T=T, r=r, sigma=sigma, N=N)
print( pvP )

es = ExerciseStyle.amer
pvC, _ = pv_crr(es, CallPut.call, S0=S0, K=K, T=T, r=r, sigma=sigma, N=N)
print( pvC )
pvP, _ = pv_crr(es, CallPut.put, S0=S0, K=K, T=T, r=r, sigma=sigma, N=N)
print( pvP )

10.448584103764572
5.571526553833685
10.448584103764572
6.089595282977953

Задача 2

Напишите код для расчета цены европейского и барьерного опционов для модели Хестона

Сделайте расчет

* implied волотильности для европейского опциона
* цены барьерного опциона knock-out 

# Параметры оционов
S0 = 80.; K = 85.; T = 1.0
B=110
#
r = 0.05
# Параметры модели
sigma0 = 0.2  # значение волатальности в начальный момент времени
long_term_variance = 0.25  # долгосрочный средний уровень дисперсии
mean_revertion_rate = 0.02 # скорость возврата к среднему
vol_of_vol=0.05  # волатильность волатильности
corr=0.07 # корреляция случайных факторов в модели

def mc_call_ko_pv_with_paths(S0, K, T, r, sigma0, long_term_variance, mean_revertion_rate, vol_of_vol,corr, M, I,B=None, is_call=True):
    # Simulating I paths with M time steps
    S = np.zeros((M+1, I))
    v = np.zeros((M+1, I))
    S[0] = S0
    v[0]=sigma0**2
    dt = float(T) / M
    for t in range(1, M + 1):
        x1 = standard_normal(I)
        x2 = standard_normal(I)
        z1=x1
        z2 = corr * x1 + np.sqrt( 1 - corr**2 ) *x2
        
        S[t] = S[t-1] * np.exp((r - 0.5 * v[t-1]) * dt + np.sqrt(v[t-1]*dt) * z1)
        v[t] = (np.sqrt(v[t-1])+0.5*vol_of_vol*np.sqrt(dt)*z2)**2 + mean_revertion_rate*(long_term_variance-v[t-1])*dt - 0.25*(vol_of_vol)**2*dt
        
        if B:
            if is_call:
                S[t]=np.array([s if s<B else 0 for s in S[t]])
            else:
                S[t]=np.array([s if s>=B else 0 for s in S[t]])

    # PV is extected discounted payoff
    if is_call:
        C = np.sum(np.exp(-r * T) * np.maximum(S[-1]-K,0)) / I
    else:
        C = np.sum(np.exp(-r * T) * np.maximum(-S[-1]+K,0)) / I
        
    return C, S

def standard_normal(I):
    z = np.random.standard_normal(I)
    mean = np.mean(z)
    std = np.std(z)
    return (z - mean)/std
    

C1,S1 = mc_call_ko_pv_with_paths(S0, K, T, r, sigma0, long_term_variance, mean_revertion_rate, vol_of_vol,corr, 360, 50000,B)

print(C1)#цена барьерного опциона call с барьером B

2.446035919361902

C2,S2 = mc_call_ko_pv_with_paths(S0, K, T, r, sigma0, long_term_variance, mean_revertion_rate, vol_of_vol,corr, 360, 50000)

print(C2)#цена обычного опциона call

6.099997795375171

from scipy import stats

def N(x):
    return stats.norm.cdf(x, 0.0, 1.0)

def NPrime(x):
    return stats.norm.pdf(x, 0.0, 1.0)

def bsm_d1(S, K, T, r, q, sigma):
    return (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

def bsm_d2(S, K, T, r, q, sigma):
    return (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

def bsm_pv(isCall, S, K, T, r, q, sigma):
    d1 = bsm_d1(S, K, T, r, q, sigma)
    d2 = bsm_d2(S, K, T, r, q, sigma)
    if isCall:
        return S * N(d1) * np.exp(-q * T) - K * np.exp(-r * T) * N(d2)
    else:
        return K * N(-d2) * np.exp(-r * T)  - S * np.exp(-q * T) * N(-d1)

from scipy.optimize import brentq

S0 = 80.; K = 85.; T = 1.0; r = 0.05

q = brentq(lambda impvol: bsm_pv(1, S0, K, T, r, 0, impvol)-C2, -1, 1)

print(q)#волатильность для обычного опциона call

0.20350564483379144

Модель задается уравнениями:

$$dS_t = \mu S_t\,dt + \sqrt{\nu_t} S_t\,dW_1 \,$$$$d\nu_t = \theta(\omega - \nu_t)dt + \xi \sqrt{\nu_t}\,dW_2$$

где:

• $\omega$ - долгосрочный средний уровень (long-term mean)

• $\theta$ - скорость возврата к среднему

• $\xi$ - волатильность волатильности

• $W_1 \,$ и $W_2 \,$ винеровские процессы, с корреляцией $\rho \,$.

$$\left< dW_1\, dW_2 \right> = \rho dt$$

Есть разные схемы дискретизации для процесса $\nu$ в модели Хестона.

Для процесса $\nu$ в модели Хестона будем использовать такой вариант схемы дискретизации:

$$\nu_{t+dt} = \left( \sqrt{\nu_t} + \frac{1}{2}\xi\sqrt{dt}z_2 \right)^2 + \theta (\omega - \nu_t) dt - \frac{1}{4}\xi^2 dt$$

где $z_2$ случайная величена распределенная по нормальному закону (соответствует процессу $W_2$ )

См. Rouah F D. Euler and Milstein discretization. http://www.frouah.com/finance notes/Euler and Milstein Discretization.pdf раздел "2 Milstein Scheme"

Обратите внимание, что для генерации двух процессов нужно сгенерировать две случайные величены $z_1$ и $z_2$, которые бы были скорелированны с коэффицентом $\rho$