2020HW-Options1.ipynb

Kernel: Python 3 (system-wide)

В задачах 1 и 2 есть замечания по их расширению/уличшению решения. Эти требования дополнительные. Решение будет полностью защитано и если если задача без учета замечания полностью решена.

Задача 1

Напишите код, который используя метод Монте-Карло расчитывал бы цену опциона EKI/KO.

EKI/KO опцион это Европейский опцион, для которого добавлено два барьера. Непрерывный барьер - действует точно так же как в KO опционе. EKI барьер действует только в момент экспирации опциона : для call опциона выплота происходит только если курс актива на момент экспирации больше уровня EKI, для опциона put выплата происходит только если курс ниже уровня EKI на момент экприрации.

Сделайте расчет для опциона call с параметрами указанными ниже

from __future__ import print_function

import numpy as np
import scipy
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import scipy.linalg as linalg
%matplotlib inline


import math as m
from time import time
from scipy import stats

S0 = 80.; K = 85.; T = 1.0; r = 0.05; q = 0.0; sigma = 0.2
B = 110 # continues barrier level
EKI = 90 # EKI level
M = 360 # number of steps in time
I = 50000 # number of paths to generate

def mc_call_eki_ko_pv_with_paths(S0, K, T, r, sigma, M, I,EKI,B):
    # Simulating I paths with M time steps
    S = np.zeros((2, I))#будем запоминать только текущие два временных слоя, тем самым по памяти получим 2*I = O(I)
    S[0] = S0
    dt = float(T) / M
    for t in range(1, M + 1):
        z = standard_normal(I)
        S[1] = S[0] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * m.sqrt(dt) * z)
        S[1]=np.array([s if s<B else 0 for s in S[1]])
        S[0]=S[1]

    S[0]=np.array([s if s>EKI else 0 for s in S[0]])
    # PV is extected discounted payoff
    C = np.sum(m.exp(-r * T) * np.maximum(S[0]-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


np.random.seed(12345)
t0 = time()
C, SPaths = mc_call_eki_ko_pv_with_paths(S0, K, T, r, sigma, M, I,EKI,B)
calcTime = time() - t0

print( "PV: %.5f" % C )
print( "Calculation time   %.5f" % calcTime )

PV: 2.27471
Calculation time   10.26824

замечание

Код в задачах из лекции использует $O( M \cdot I )$ памяти. Попробуйте сделать код, который использует $O( I )$ памяти.

Задача 2

Модифицируйте функцию, разобранную на лекции, для расчета цены опциона через решения уравнения в частных производных таким образом, что бы можно было задавать верхний и/или нижний барьер.

Возвращаемые значения функции оставить как в примере на лекции.

Сигнатура модифицированной функции и прамеры и примеры расчетов, которые нужно сделать:

замечание

Можно ли вычислить PV и Delta опциона в одной функции не увеличив существенно колличечество вычислений? Если получится, то сделайте функцию, которая возвращает не только цену опциона, но и его дельту.

def priceWithPDE1(S0, K, r, T, sigma, M, N, DB=None, UB=None, Smax=None, is_call=True):
    # TODO

    if DB and UB:
        S = np.linspace(DB, UB, M+1)
        Smax=UB
    elif DB:
        S = np.linspace(DB, Smax, M+1)
    elif UB:
        S = np.linspace(0, UB, M+1)
        Smax=UB
    else:
        S = np.linspace(0, Smax, M+1)

    M, N = int(M), int(N)
    dS = float(Smax) / M
    dt = float(T) / N
    k_values = np.arange(M)
    j_values = np.arange(N)
    grid = np.zeros(shape=(M+1, N+1))

    # setup boundary conditions
    if is_call:
        grid[:, -1] = np.maximum(S - K, 0)
        grid[-1, :-1] = Smax - K* np.exp(-r *dt*(N-j_values))
    else:
        grid[:, -1] = np.maximum(K-S, 0)

        grid[0, :-1] = K * np.exp(-r *dt *(N-j_values))

    # calculate coefficients and matrix to solve
    alpha = 0.25*dt*((sigma**2)*(k_values**2) -r*k_values)

    beta = -dt*0.5*((sigma**2)*(k_values**2) + r)         

    gamma = 0.25*dt*((sigma**2)*(k_values**2) +r*k_values)
    

    A1 = -np.diag(alpha[2:M], -1) + \
        np.diag(1-beta[1:M]) - np.diag(gamma[1:M-1], 1)

    A2 = np.diag(alpha[2:M], -1) + \
        np.diag(1+beta[1:M]) + np.diag(gamma[1:M-1], 1)
        
    #print(A1, "A1")
    #print(A2, "A2")

    print(f"{alpha.shape} {beta.shape} {gamma.shape}")
    print(f"{A1.shape} {A2.shape}")
    
    # solve systems of linear equstions 
    for j in reversed(range(N)):
        B = A2 @ grid[1:M, j+1]
        B[0] += alpha[0]*(grid[0,j+1] + grid[0,j])
        B[-1] += gamma[M-1]*(grid[-1,j+1] + grid[-1,j])
        x2 = linalg.solve(A1, B)
        grid[1:M, j] = x2
        

    
    if DB:
        #print("DB")
        for i in range(M+1):
            grid[i,:]=0 if S[i]<=float(DB) else grid[i,:]
                
    if UB:
        #print("UB")
        for i in range(M+1):
            grid[i,:]=0 if S[i]>=float(UB) else grid[i,:]          
            
    print(np.where(grid==0)[0].size)
    plt.figure(figsize=(12, 12))
    plt.imshow(grid)
        
        

    # use linear interpolation to get result value 
    pv_0 = np.interp(S0, S, grid[:, 0])
    return pv_0, grid, S 

S0 = 80.; K = 85.; T = 1.0; r = 0.05; sigma = 0.2
UB = 110
DB = 70
Smax = 160.
M=360; N=100

print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, Smax=Smax, is_call=True)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, Smax=Smax, is_call=False)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, UB=UB, is_call=True)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, UB=UB, is_call=False)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, DB=DB, Smax=Smax, is_call=True)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, DB=DB, Smax=Smax, is_call=False)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, DB=DB, UB=UB, is_call=False)[0] )
print( priceWithPDE1(S0=S0, K=K, r=r, T=T, sigma=sigma,M=M, N=N, DB=DB, UB=UB, is_call=True)[0] )

(360,) (360,) (360,)
(359, 359) (359, 359)
292
5.98810225112211
(360,) (360,) (360,)
(359, 359) (359, 359)
269
6.84260324946887
(360,) (360,) (360,)
(359, 359) (359, 359)
480
5.979843283425456
(360,) (360,) (360,)
(359, 359) (359, 359)
182
6.834344285584508
(360,) (360,) (360,)
(359, 359) (359, 359)
161
0.0361192858548146
(360,) (360,) (360,)
(359, 359) (359, 359)
502
4.304560638502548
(360,) (360,) (360,)
(359, 359) (359, 359)
427
4.304434469271366
(360,) (360,) (360,)
(359, 359) (359, 359)
337
0.03599311662546681

Задача 3

Рассмотрим такую стратегию динамического хеджирования: в начальный момент времени мы покупаем базовый актив, в колличестве равном дельте call опциона. Далее мы перебалансируем наш портфель с каким-то шагом по времени. На каждом шаге мы покупаем или продаем базовый актив из такого расчета, что бы у нас на балансе опять был базовый актив, в колличестве равном дельте call опциона, с учетом изменившегося курса актива.

Используя симмуляции возможных реализаций процесса цены актива как в методе Монте-Карло сделайте код, который нарисует график для суммарной стоимости владания такой стратегией для некоторого кол-ва симмуляций.

График должен представлят из себя диаграмму рассеяния (scatter plot) в координатах $x = S_T$, $y = Total Profit_T$

Параметры для расчета:

S0 = 80.; K = 85.; T = 1.0; r = 0.05; q = 0.0; sigma = 0.2
M = 360 # number of steps in time
I = 1000 # number of paths to generate

Знаем, что для опциона Call $\Delta = e^{-q\tau}N(d_+)$, где $d_+ = \frac{log{\frac{S}{K}} + (r-q+\frac{\sigma^2}{2})\tau}{\sigma \sqrt\tau}$, а $N(x)$ - функция стандратного нормального распределения

В этом случае $t=0, T=1 => \tau = 1, S_t = S_0 = 80, r = 0.05, \sigma = 0.2, q = 0.01$

def Delta(S_0, t, T, q, r, sigma, K):
    return np.exp(-q*(T-t))*N( bsm_d1(S_0, K, T-t, r, q, sigma))

def N(x):
    return stats.norm.cdf(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 * m.sqrt(T)) 

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

# Значение spot на текущем шаге симуляции
S = np.full(I, S0)
# Значение payoff на текущем шаге симуляции
P = np.zeros(I)

delta = Delta(S,0, T, q, r, sigma,K)
temp = -S*delta

# TODO: set initial portfolio
#
P = np.full(I,-S0*delta)

dt = float(T) / M
for t in range(1, M):
    z = standard_normal(I)
    delta_prev = -temp/S
    S = S * np.exp((r - q - 0.5 * sigma ** 2) * dt + sigma * m.sqrt(dt) * z)
    # TODO: rebalance portfolio for the current step
    delta_new = Delta(S,t*dt, T,q,r,sigma,K)
    dif = (delta_new-delta_prev)*S
    P = P - dif
    temp = -S*delta_new
    P = P - temp

import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline

lw = 2.5 # line width
fs = 12  # fontsize
picWidth=12; picHeight=8;

matplotlib.rcParams.update({'font.size': 22})
plt.figure(figsize=(picWidth, picHeight))

plt.plot(S, P, 'bo', linewidth=lw)
plt.grid(True)
plt.xlabel('Spot at expiry')
plt.ylabel('Total payoff')

Text(0, 0.5, 'Total payoff')

Задача 4

Straddle это структура из опциона call и опционна put, у которых равны страйки.

Напишите код, который находит такой страйк для структуры straddle, при котором delta структуры равна нулю. (Delta neutral straddle). Для численного нахождения корней функции можно использовать любой метод из стантартных библиотек, например https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html

Параметры для расчета:

S0 = 80.; K = 85.; T = 1.0; r = 0.05; q = 0.01; sigma = 0.2
#Не очень понятно, зачем нам дан страйк

Put-Call Parity: $C - P = Se^{-q\tau}-Ke^{-r\tau}$

Тогда $\frac{\partial C}{\partial S} - \frac{\partial P}{\partial S} = e^{-q\tau}$ => $\frac{\partial Str}{\partial S} = 2 \frac{\partial C}{\partial S} - 1$ => $\frac{\partial Str}{\partial S} = e^{-q\tau}(2N(d_+)-1)$,

где $d_+ = \frac{log{\frac{S}{K}} + (r-q+\frac{\sigma^2}{2})\tau}{\sigma \sqrt\tau}$, а $N(x)$ - функция стандратного нормального распределения

То есть если мы значем $S_t$ и остальные параметры $r, \sigma, q, \tau$, то можем с помошью численных методов найти приближенно тот страйк $K$, которые обесчпечит $\Delta_{Str} = \frac{\partial Str}{\partial S} = 0$

В этом случае $t=0, T=1 => \tau = 1, S_t = S_0 = 80, r = 0.05, \sigma = 0.2, q = 0.01$

def Delta(S_0, t, T, q, r,sigma, K):
    return 2*N( bsm_d1(S_0, K, T, r, q, sigma))-1

def N(x):
    return stats.norm.cdf(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 * m.sqrt(T)) 

def Delta_K(K):
    return Delta(80,0,1,0.01,0.05,0.2,K)

scipy.optimize.brentq(Delta_K,1,1000)

84.94692372362877