PyMPDATA_examples.Arabas_and_Farhat_2020.Bjerksund_and_Stensland_1993

View Source

 1import numpy as np
 2import PyMPDATA_examples.Arabas_and_Farhat_2020.Black_Scholes_1973 as BS
 3
 4
 5def _phi(
 6    S: [np.ndarray, float],
 7    gamma: float,
 8    H: float,
 9    I: float,
10    r: float,
11    b: float,
12    var: float,
13    T: float,
14):
15    lmbd = (-r + gamma * b + 0.5 * gamma * (gamma - 1) * var) * T
16    d = -(np.log(S / H) + (b + (gamma - 0.5) * var) * T) / np.sqrt(var * T)
17    kappa = 2 * b / var + (2 * gamma - 1)
18    return (
19        np.exp(lmbd)
20        * np.power(S, gamma)
21        * (
22            BS.N(d)
23            - pow((I / S), kappa) * BS.N(d - 2 * np.log(I / S) / np.sqrt(var * T))
24        )
25    )
26
27
28def c_amer(
29    S: [np.ndarray, float],
30    K: [float, np.ndarray],
31    T: float,
32    r: float,
33    b: float,
34    sgma: float,
35):
36    if b >= r:
37        return BS.c_euro(S, K=K, T=T, r=r, b=b, sgma=sgma)
38
39    var = sgma * sgma
40    beta = (0.5 - b / var) + np.sqrt(pow((b / var - 0.5), 2) + 2 * r / var)
41    BInf = beta / (beta - 1) * K
42    B0 = np.maximum(K, r / (r - b) * K)
43    ht = -(b * T + 2 * sgma * np.sqrt(T)) * B0 / (BInf - B0)
44    I = B0 + (BInf - B0) * (1 - np.exp(ht))
45    alpha = (I - K) * pow(I, -beta)
46
47    return np.where(
48        S >= I,
49        S - K,
50        alpha * np.power(S, beta)
51        + (
52            -alpha * _phi(S, gamma=beta, H=I, I=I, r=r, b=b, var=var, T=T)
53            + _phi(S, gamma=1, H=I, I=I, r=r, b=b, var=var, T=T)
54            - _phi(S, gamma=1, H=K, I=I, r=r, b=b, var=var, T=T)
55            - K * _phi(S, gamma=0, H=I, I=I, r=r, b=b, var=var, T=T)
56            + K * _phi(S, gamma=0, H=K, I=I, r=r, b=b, var=var, T=T)
57        ),
58    )
59
60
61def p_amer(S: [np.ndarray, float], K: float, T: float, r: float, b: float, sgma: float):
62    return c_amer(K, K=S, T=T, r=r - b, b=-b, sgma=sgma)

def c_amer( S: [<class 'numpy.ndarray'>, <class 'float'>], K: [<class 'float'>, <class 'numpy.ndarray'>], T: float, r: float, b: float, sgma: float): View Source

29def c_amer(
30    S: [np.ndarray, float],
31    K: [float, np.ndarray],
32    T: float,
33    r: float,
34    b: float,
35    sgma: float,
36):
37    if b >= r:
38        return BS.c_euro(S, K=K, T=T, r=r, b=b, sgma=sgma)
39
40    var = sgma * sgma
41    beta = (0.5 - b / var) + np.sqrt(pow((b / var - 0.5), 2) + 2 * r / var)
42    BInf = beta / (beta - 1) * K
43    B0 = np.maximum(K, r / (r - b) * K)
44    ht = -(b * T + 2 * sgma * np.sqrt(T)) * B0 / (BInf - B0)
45    I = B0 + (BInf - B0) * (1 - np.exp(ht))
46    alpha = (I - K) * pow(I, -beta)
47
48    return np.where(
49        S >= I,
50        S - K,
51        alpha * np.power(S, beta)
52        + (
53            -alpha * _phi(S, gamma=beta, H=I, I=I, r=r, b=b, var=var, T=T)
54            + _phi(S, gamma=1, H=I, I=I, r=r, b=b, var=var, T=T)
55            - _phi(S, gamma=1, H=K, I=I, r=r, b=b, var=var, T=T)
56            - K * _phi(S, gamma=0, H=I, I=I, r=r, b=b, var=var, T=T)
57            + K * _phi(S, gamma=0, H=K, I=I, r=r, b=b, var=var, T=T)
58        ),
59    )

def p_amer( S: [<class 'numpy.ndarray'>, <class 'float'>], K: float, T: float, r: float, b: float, sgma: float): View Source

62def p_amer(S: [np.ndarray, float], K: float, T: float, r: float, b: float, sgma: float):
63    return c_amer(K, K=S, T=T, r=r - b, b=-b, sgma=sgma)