<jupyter_start><jupyter_text>BSCSIR AES<jupyter_code>def bates_SC_SIR_AES(numberPaths, N, s0, v0, T, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5,
                    xip, muJ, sigmaJ, r0, rho0):
    X = np.zeros([numberPaths, N + 1])
    S = np.zeros([numberPaths, N + 1])
    V = np.zeros([numberPaths, N + 1])
    R = np.zeros([numberPaths, N + 1])
    rho = np.zeros([numberPaths, N + 1])

    M_t = np.ones([numberPaths, N + 1])
    time = np.zeros(N + 1)
    dt = T/float(N)

    Zx = np.random.normal(0, 1, [numberPaths, N])
    Zrho = np.random.normal(0, 1, [numberPaths, N])
    ZP = np.random.poisson(xip*dt, [numberPaths, N])
    J = np.random.normal(muJ, sigmaJ, [numberPaths, N])

    X[:, 0] = np.log(s0)
    V[:, 0] = v0
    rho[:, 0] = rho0
    R[:, 0] = abs(r0)

    EeJ = np.exp(muJ + 0.5 * sigmaJ**2)

    for t in range(N):

        V[:, t+1] = CIR_exact(numberPaths, k, gamma, vb, 0, dt, V[:, t])
        R[:, t+1] = CIR_exact(numberPaths, kr, gammar, mur, 0, dt, R[:, t])
        
        if (R[:, t+1] == 0).any():
            R[np.where(R[:, t+1] == 0)[0], t+1] = 1e-4

        rho[:, t+1] = rho[:, t]*np.exp(-krho*dt) + murho*(1 - np.exp(-krho*dt)) +\
                    sigmarho*np.sqrt((1-np.exp(-2*krho*dt))/(2*krho)) * Zrho[:, t]

        if (rho[:, t+1] > 1).any():
            rho[np.where(rho[:, t+1] > 1)[0], t+1] = 0.9999

        if (rho[:, t+1] < -1).any():
            rho[np.where(rho[:, t+1] < -1)[0], t+1] = -0.9999

        X[:, t+1] = X[:, t] + (R[:, t] - 0.5* V[:, t] - xip*(EeJ-1))*dt + rho[:, t]/gamma * (V[:, t+1] - V[:, t] - k*(vb - V[:, t])*dt) + \
            rho4 * np.sqrt(V[:, t])/sigmarho * (rho[:, t+1] - rho[:, t] - krho*(murho - rho[:, t])*dt) + \
            rho5 * np.sqrt(V[:, t])/(gammar*np.sqrt(R[:, t])) * (R[:, t+1] - R[:, t] - kr*(mur-R[:, t])*dt) +\
            np.sqrt(V[:, t])* np.sqrt(abs(1 - (rho[:, t])**2 - rho4**2 - rho5**2)) * np.sqrt(dt) * Zx[:, t] + J[:, t] * ZP[:, t]

        M_t[:,t+1] = M_t[:,t] * np.exp(0.5*(R[:,t+1] + R[:,t])*dt)
        time[t+1] = time[t] + dt

    S = np.exp(X)

    return time, S, M_t
# v0, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5, muJ, sigmaJ, r0, rho0
numberPaths = 1000
N = 100
T = maturity
type_option = 'c'
# xip = 0.05

def error_BSCSIR_AES(x, prices, K, s0):
    time, S, M = mcm.bates_SC_SIR_AES(numberPaths, N, s0, x[0], T, x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11],
                    x[16], x[12], x[13], x[14], x[15])
    
    price_calib = mcm.optionPriceMC_Stoch(type_option, S[:, -1], K, T, M[:, -1])
    
    return np.mean((price_calib - prices)**2)

init_vals = np.array([0.04, 0.3, 0.04, 0.3, 0.1, 0.2, 0.5, 0.05, 0.4, 0.5, 0.1, 0.1, 0, 1, 0.001, -0.7, 0.05])

xmin = [0.001, 0.0001, 0.0001, 0.0005, 0.001, 0.001, 0.0005, 0.0001, 0.0001, 0.4, -0.89, -0.89, -10, 0.0001, 0.0001, -0.99, 1e-6]
xmax = [1, 0.7, 1.5, 2, 3, 2, 1, 1.7, 4, 2, 0.89, 0.89, 10, 3, 0.3, 0.99, 5]
bounds = [(low, high) for low, high in zip(xmin, xmax)]

params_BSCSIR = minimize(error_BSCSIR_AES, x0=init_vals, args=(marketPrice, K, s0), bounds=bounds, tol=1e-8, options={"maxiter": 100})
params_BSCSIR
v0, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5, muJ, sigmaJ, r0, rho0, xip = params_BSCSIR.x
paths_BSCSIR = mcm.bates_SC_SIR_AES(numberPaths, N, s0, v0, T, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5,
                    xip, muJ, sigmaJ, r0, rho0)
option_BSCSIR_AES = mcm.optionPriceMC_Stoch(type_option, paths_BSCSIR[1][:, -1], K, T, paths_BSCSIR[2][:, -1])
plt.subplots(figsize=(10, 5), dpi=100)

plt.plot(K, marketPrice, label='Initial Prices')
plt.plot(K, option_BSCSIR_AES, color='green', label='BSCSIR AES')
        
plt.title('Option Prices for AAPL with BSCSIR AES', fontsize=16)
plt.xlabel('Strikes', fontsize=14)
plt.ylabel('Values', fontsize=14)
plt.legend()
plt.show()
list_params = np.array([v0, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5, muJ, sigmaJ, r0, rho0, xip])
pd.DataFrame(list_params.reshape(1, -1), columns=['v0', 'k', 'gamma', 'vb', 'kr', 'gammar', 'mur', 'krho', 'murho', 'sigmarho', 'rho4', 'rho5', 'muJ', 'sigmaJ', 'r0', 'rho0', 'xip'])<jupyter_output><empty_output><jupyter_text>BSCSIR with stochastic correlation by DCL process<jupyter_code>def calibrationBates_SIR_SC_DCL(CP, xip, K, marketPrice, s0, T):
    K = np.array(K)
    marketPrice = np.array(marketPrice)
   
    # k, gamma, vb, kr, gammar, mur, krho, delta, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0
    f_obj = lambda x: errorFun(CP, tau, T, x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9],
                               x[10], x[11], x[12], x[13], x[14], x[15], K, marketPrice, s0)

                        # k, gamma, vb, kr, gammar, mur, krho, delta, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0
    initial = np.array([0.4, 0.3, 0.04, 0.3, 0.2,   0.2,   2,   1.5,  0.1, 0.5, 0.05,  0,    1, 0.04,   0.05, 0.7])

            # k,    gamma,   vb,     kr,   gammar, mur, krho, delta, rho4, rho5,  xip,    muJ, sigmaJ, v0, r0, rho0
    xmin = [0.0001, 0.0001, 0.0001, 0.0005, 0.001, 0.0001, 1, 1,     -0.99, -0.99, 0.0001, -10,   -10,   0.0001, 0.001, -0.9999]
    xmax = [5,       0.7,   1.5,     2,      3,     2,   5,  4,      0.98, 0.98,  5,      10, 10,      0.3, 0.3, 0.1] 
    bounds = [(low, high) for low, high in zip(xmin, xmax)]
    
    minimizer_kwargs = dict(method="L-BFGS-B", bounds=bounds)
    
    pars  = minimize(f_obj, x0=initial, tol=1e-8, bounds=bounds, options = {'maxiter': 200, 'disp':False})
    print(pars)
    
    k_est = pars.x[0]
    gamma_est = pars.x[1]
    vb_est = pars.x[2]
    kr_est = pars.x[3]
    gammar_est = pars.x[4]
    mur_est = pars.x[5]
    krho_est = pars.x[6]
    delta_est = pars.x[7]
    rho4_est = pars.x[8]
    rho5_est = pars.x[9]
    xip_est = pars.x[10]
    muJ_est = pars.x[11]
    sigmaJ_est = pars.x[12]
    v0_est = pars.x[13]
    r0_est = pars.x[14]
    rho0_est = pars.x[15]
   
   
    
    parmCalibr =  {"k":k_est, "gamma":gamma_est, "vb":vb_est, "gammar":gammar_est, "mur":mur_est,\
                   "rho4":rho4_est, "rho5": rho5_est, "v0":v0_est, "kr": kr_est, "krho": krho_est, \
                   "rho0": rho0_est, "r0": r0_est, "muJ": muJ_est, "sigmaJ": sigmaJ_est, "xip": xip_est, "delta": delta_est, 'ErrorFinal':pars.fun}
    return parmCalibr
def errorFun(CP, tau, T, k, gamma, vb, kr, gammar, mur, krho, delta, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0, K, marketPrice, s0):
    if K is not np.array:
        K = np.array(K).reshape([len(K),1])

    N = 500
    L = 10
    cf = chf.ChFBates_StochIR_StochCor_DCL(tau, T, k, gamma, vb, kr, gammar, mur, krho, delta, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0)
    valCOS = exm.optionPriceCOSMthd_StochIR(cf, CP, s0, T, K, N, L)

    errorVector = np.mean((valCOS.T[0] - marketPrice)**2)
    value = errorVector

    return value
T = maturity
tau = T
type_option = 'c'

# P0T = lambda T: np.exp(0.0005*T)

# dt = 0.0001 
# # P0T = lambda T: np.exp(0.0012*T+0.0007) 
# f0T = lambda t: - (np.log(P0T(t+dt))-np.log(P0T(t-dt)))/(2.0*dt)
# r0 = f0T(0.00001)

K = apple_strikes

marketPrice = apple_market_price

calibratedParms =  calibrationBates_SIR_SC_DCL(type_option, xip, K, marketPrice, s0, T)

k = calibratedParms.get('k')
gamma = calibratedParms.get('gamma')
vb = calibratedParms.get('vb')
gammar = calibratedParms.get('gammar')
mur = calibratedParms.get('mur')
rho4 = calibratedParms.get('rho4')
rho5 = calibratedParms.get('rho5')
v0 = calibratedParms.get('v0')
kr = calibratedParms.get('kr')
krho = calibratedParms.get('krho')
rho0 = calibratedParms.get('rho0')
r0 = calibratedParms.get('r0')
muJ = calibratedParms.get('muJ')
sigmaJ = calibratedParms.get('sigmaJ')
xip = calibratedParms.get('xip')
delta = calibratedParms.get('delta')

errorFinal = calibratedParms.get('ErrorFinal') 
np.sqrt(errorFinal)
round(errorFinal, 9)
L = 10
N = 500

cf = chf.ChFBates_StochIR_StochCor_DCL(tau, T, k, gamma, vb, kr, gammar, mur, krho, delta, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0)
valCOS_DCL = exm.optionPriceCOSMthd_StochIR(cf, type_option, s0, T, K, N, L)
plt.subplots(figsize=(10, 5), dpi=100)

plt.plot(K, marketPrice, label='Initial Prices')
plt.plot(K, valCOS_DCL, '--', color='green', label='Bates Stoch Corr Stoch IR')
        
plt.title('Option Pricing for AAPL with Bates Stoch. Corr. DCL Stoch. IR ', fontsize=16)
plt.xlabel('Strikes', fontsize=14)
plt.ylabel('Values', fontsize=14)
plt.legend()
plt.show()
list_params = np.array([v0, k, gamma, vb, kr, gammar, mur, krho, delta, rho4, rho5, muJ, sigmaJ, r0, rho0, xip])
pd.DataFrame(list_params.reshape(1, -1), columns=['v0', 'k', 'gamma', 'vb', 'kr', 'gammar', 'mur', 'krho', 'delta', 'rho4', 'rho5', 'muJ', 'sigmaJ', 'r0', 'rho0', 'xip'])
# v0, k, gamma, vb, kr, gammar, mur, theta, delta, rho4, rho5, muJ, sigmaJ, r0, rho0, xip
numberPaths = 10000
N = 100
T = maturity
type_option = 'c'

def error_BSCSIR_DCL_AES(x, prices, K, s0):
    time, S, M = mcm.bates_SC_SIR_DCL_AES(numberPaths, N, s0, x[0], T, x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[15],
                    x[11], x[12], x[13], x[14])
    
    price_calib = mcm.optionPriceMC_Stoch(type_option, S[:, -1], K, T, M[:, -1])
    
    return np.mean((price_calib - prices)**2)

init_vals = np.array([0.04, 0.3, 0.04, 0.3, 0.1, 0.2, 0.5, 0.1, 1, 0.1, 0.1, 0, 1, 0.001, -0.7, 0.05])

xmin = [0.001, 0.0001, 0.0001, 0.0005, 0.001, 0.001, 0.0005, 0.1, 1, -0.89, -0.89, -10, 0.0001, 0.0001, -0.99, 1e-5]
xmax = [1, 0.7, 1.5, 2, 3, 2, 1, 5, 3, 0.89, 0.89, 10, 3, 0.3, 0.99, 5]
bounds = [(low, high) for low, high in zip(xmin, xmax)]

params_BSCSIR_DCL = minimize(error_BSCSIR_DCL_AES, x0=init_vals, args=(marketPrice, K, s0), bounds=bounds, tol=1e-4, options={"maxiter": 1000})
params_BSCSIR_DCL
v0, k, gamma, vb, kr, gammar, mur, theta, delta, rho4, rho5, muJ, sigmaJ, r0, rho0, xip = params_BSCSIR_DCL.x
paths_BSCSIR_DCL = mcm.bates_SC_SIR_DCL_AES(numberPaths, N, s0, v0, T, k, gamma, vb, kr, gammar, mur, theta, delta, rho4, rho5,
                    xip, muJ, sigmaJ, r0, rho0)
option_BSCSIR_DCL_AES = mcm.optionPriceMC_Stoch(type_option, paths_BSCSIR_DCL[1][:, -1], K, T, paths_BSCSIR_DCL[2][:, -1])
plt.subplots(figsize=(10, 5), dpi=100)

plt.plot(K, marketPrice, label='Initial Prices')
plt.plot(K, option_BSCSIR_DCL_AES, color='green', label='BSCSIR with DCL AES')
        
plt.title('Option Prices for AAPL by BSCSIR with DCL AES', fontsize=16)
plt.xlabel('Strikes', fontsize=14)
plt.ylabel('Values', fontsize=14)
plt.legend()
plt.show()
tsla = ops.get_calls("TSLA", 'May 13, 2022')
# tsla.to_csv('tsla_option_call_may_13_2022.csv')
# s0_tsla = 766.37 #15.03.2022
tslaticker = yf.Ticker('TSLA')
s0_tsla = tslaticker.history().iloc[-1, 3]
s0_tsla
tsla_strikes = tsla.loc[:, 'Strike']
tsla_market_price = tsla.loc[:, 'Last Price']
tsla.head()
plt.plot(tsla_strikes, tsla_market_price)
plt.show()
from datetime import datetime
expr_date = datetime.strptime('May 13, 2022', '%B %d, %Y')
expr_date
today_date = datetime.now()
today_date
maturity = expr_date - today_date
maturity = maturity.days/365.25
maturity
(expr_date - today_date).days
N = 500
T = maturity
r = 0.01
rho1 = 0.2
s0 = s0_tsla

type_option = 'c'
# type_option = 'p'
# r, kappa, gamma, vbar, v0, rho

def error_fHM(x, prices, strikes, s0):
    cf = chf.ChFHestonModel(x[0], 1, x[1], x[2], x[3], x[4], x[5])
    
    price_calib = exm.CallPutOptionPriceCOS(cf, type_option, s0, x[0], T, strikes, N, 8).T[0]
#     errorVector = (price_calib - prices)
    errorVector = np.mean(((price_calib - prices))**2)
    
#     value = np.linalg.norm(errorVector)
    value = errorVector
    return value
init_vals = [0.05, 2, 0.3, 0.4, 0.3, -0.7]
bounds = ((0.001, 0.3), (0.001, 4), (0.0001, 2), (0.0001, 1.5), (0.0001, 1.5), (-0.9, 0.9))

# P0T = lambda T: np.exp(0.0012*T+0.0007) 
# s0 = 174.62

# referenceIV = np.array([48.53, 24.49, 19.23])/100.0
        
# T = 5/365.25
# referenceIV = np.array([36.76, 29.18, 27.66])/100.0
    
# frwd = s0/P0T(T)
# K = [0.4*frwd, frwd, 1.2*frwd]
# K = np.array(K).reshape([len(K),1])
      
# referenceIV = np.array(referenceIV).reshape([len(referenceIV),1])
# marketPrice = P0T(T)* exm.BS_Call_Option_Price('c', s0 / P0T(T), K, referenceIV, T, 0)


# params_heston_COS = minimize(error_fHM, x0=init_vals, args=(apple_market_price, apple_strikes, s0_apple, apple_spread),
#                              bounds=bounds, tol=1e-4, options={"maxiter": 1000})

K = tsla_strikes
marketPrice = tsla_market_price

params_heston_COS = minimize(error_fHM, x0=init_vals, args=(marketPrice, K, s0),
                             bounds=bounds, tol=1e-4, options={"maxiter": 1000})

params_heston_COS
np.sqrt(params_heston_COS.fun)
r, kappa, gamma, vbar, v0, rho = params_heston_COS.x
cf = chf.ChFHestonModel(r, 1, kappa, gamma, vbar, v0, rho)
option_heston_cos = exm.CallPutOptionPriceCOS(cf, type_option, s0, r, T, K, N, 10).T[0]
plt.subplots(figsize=(10, 5), dpi=100)

plt.plot(K, marketPrice, label='Initial Prices')
plt.plot(K, option_heston_cos, color='green', label='Heston COS')
        
plt.title('Option Pricing for AAPL with Heston COS', fontsize=16)
plt.xlabel('Strikes', fontsize=14)
plt.ylabel('Values', fontsize=14)
plt.legend()
plt.show()
T = maturity
tau = T

xip = 1
type_option = 'c'

P0T = lambda T: np.exp(0.05*T)

# dt = 0.0001 
# # P0T = lambda T: np.exp(0.0012*T+0.0007) 
# f0T = lambda t: - (np.log(P0T(t+dt))-np.log(P0T(t-dt)))/(2.0*dt)
# r0 = f0T(0.00001)


K = tsla_strikes

marketPrice = tsla_market_price

calibratedParms =  calibrationBates_SIR_SC(type_option, xip, K, marketPrice, s0, T, P0T)

k = calibratedParms.get('k')
gamma = calibratedParms.get('gamma')
vb = calibratedParms.get('vb')
gammar = calibratedParms.get('gammar')
mur = calibratedParms.get('mur')
murho = calibratedParms.get('murho')
sigmarho = calibratedParms.get('sigmarho')
rho4 = calibratedParms.get('rho4')
rho5 = calibratedParms.get('rho5')
v0 = calibratedParms.get('v0')
kr = calibratedParms.get('kr')
krho = calibratedParms.get('krho')
rho0 = calibratedParms.get('rho0')
r0 = calibratedParms.get('r0')
muJ = calibratedParms.get('muJ')
sigmaJ = calibratedParms.get('sigmaJ')

errorFinal = calibratedParms.get('ErrorFinal') 
np.sqrt(errorFinal)
round(errorFinal, 9)
L = 10
N = 500

cf = chf.ChFBates_StochIR_StochCor(tau, T, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0)
valCOS = exm.optionPriceCOSMthd_StochIR(cf, type_option, s0, T, K, N, L, P0T(T))
plt.subplots(figsize=(10, 5), dpi=100)

plt.plot(K, marketPrice, label='Initial Prices')
plt.plot(K, valCOS, '--', color='green', label='Bates Stoch Corr Stoch IR')
        
plt.title('Option Pricing for AAPL with Bates Stoch Corr Stoch IR', fontsize=16)
plt.xlabel('Strikes', fontsize=14)
plt.ylabel('Values', fontsize=14)
plt.legend()
plt.show()
zcb0 = zcb_curve(tau, kr, mur, gammar, r0)
zcb0
K*zcb0 + s0 
#Zero-coupon parity
def optionPriceCOSMthd_StochIR(cf, CP, s0,tau,K,N,L,P0T):
    if K is not np.array:
        K = np.array(K).reshape([len(K),1])

    i = complex(0.0,1.0)
    x0 = np.log(s0 / K)

    # Truncation domain
    a = 0.0 - L * np.sqrt(tau)
    b = 0.0 + L * np.sqrt(tau)

    # Summation from k = 0 to k=N-1
    k = np.linspace(0,N-1,N).reshape([N,1])
    u = k * np.pi / (b - a)

    H_k = exm.Hk_Coefficients('p',a,b,k)
    mat = np.exp(i * np.outer((x0 - a) , u))
    temp = cf(u) * H_k
    temp[0] = 0.5 * temp[0]
    value = K * np.real(mat.dot(temp))
        
    if CP == 'c':
        value = value + s0 - K * P0T  

    return value
L = 10
N = 500

cf = chf.ChFBates_StochIR_StochCor(tau, T, k, gamma, vb, kr, gammar, mur, krho, murho, sigmarho, rho4, rho5, xip, muJ, sigmaJ, v0, r0, rho0)
valCOS = optionPriceCOSMthd_StochIR(cf, type_option, s0, T, K, N, L, zcb0)
plt.subplots(figsize=(10, 5), dpi=100)

plt.plot(K, marketPrice, label='Initial Prices')
plt.plot(K, valCOS, '--', color='green', label='Bates Stoch Corr Stoch IR')
        
plt.title('Option Pricing for AAPL with Bates Stoch Corr Stoch IR', fontsize=16)
plt.xlabel('Strikes', fontsize=14)
plt.ylabel('Values', fontsize=14)
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1980706202461196, 'fraction_numerical': 0.07092944159297508, 'mean_word_length': 3.502227171492205, 'pattern_counts': {'":': 20, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '294592', 'license': 'MIT', 'n_tokens_mistral': 8214, 'n_tokens_neox': 7417, 'n_words': 1648, 'path': 'Calibrations/Calibration_Main_Calls.ipynb', 'repo_name': 'kirillzx/Stochastic-Differential-Equations', 'size': '16171'}