Option

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import scipy.stats as scs
from math import *
%matplotlib inline

S0 = 100 # initial value
r = 0.05 # constant short rate
sigma = 0.25 # constant volatility
T = 1.0 # in years
I = 10000 # number of random draws
ST1 = S0 * np.exp((r - 0.5 * sigma ** 2) * T + sigma * np.sqrt(T) * npr.standard_normal(I))
ST2 = S0 * npr.lognormal((r - 0.5 * sigma ** 2) * T, sigma * np.sqrt(T), size=I)

#plt.hist(ST1, bins=50)
#plt.xlabel(‘index level’)
#plt.grid(True)
#plt.ylabel(‘frequency’)

#plt.hist(ST2, bins=50)
#plt.xlabel(‘index level’)
#plt.ylabel(‘frequency’)
#plt.grid(True)

def print_statistics(a1, a2):
    sta1 = scs.describe(a1)
    sta2 = scs.describe(a2)
    print ('size', sta1[0], sta2[0])
    print ('min', sta1[1][0], sta2[1][0])
    print ('max', sta1[1][1], sta2[1][1])
    print ('mean', sta1[2], sta2[2])
    print ('std', np.sqrt(sta1[3]), np.sqrt(sta2[3]))
    print ('skew', sta1[4], sta2[4])
    print ('kurtosis', sta1[5], sta2[5])
print_statistics(ST1, ST2)

#### BSM ####

I = 10000
M = 50
dt = T / M
S = np.zeros((M + 1, I))
S[0] = S0
for t in range(1, M + 1):
    S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * npr.standard_normal(I))

print_statistics(S[-1], ST2)
#plt.hist(S[-1], bins=50)
#plt.xlabel(‘index level’)
#plt.ylabel(‘frequency’)
#plt.grid(True)


#plt.plot(S[:, :10], lw=1.5)
#plt.xlabel(‘time’)
#plt.ylabel(‘index level’)
#plt.grid(True)

#### CIR - Euler ####
x0 = 0.05
kappa = 3.0
theta = 0.02
sigma = 0.1
I = 10000
M = 50
T = 1
dt = T / M

def srd_euler():
    xh = np.zeros((M + 1, I))
    x1 = np.zeros_like(xh)
    xh[0] = x0
    x1[0] = x0
    for t in range(1, M + 1):
        xh[t] = (xh[t - 1] + kappa * (theta - np.maximum(xh[t - 1], 0)) * dt + sigma * np.sqrt(np.maximum(xh[t - 1], 0)) * np.sqrt(dt) * npr.standard_normal(I))
    x1 = np.maximum(xh, 0)
    return x1
x1 = srd_euler()

#plt.hist(x1[-1], bins=50)
#plt.xlabel(‘value’)
#plt.ylabel(‘frequency’)
#plt.grid(True)

#plt.plot(x1[:, :10], lw=1.5)
#plt.ylabel(‘index level’)
#plt.xlabel(‘time’)
#plt.grid(True)

def srd_exact():
    x2 = np.zeros((M + 1, I))
    x2[0] = x0
    for t in range(1, M + 1):
        df = 4 * theta * kappa / sigma ** 2
        c = (sigma ** 2 * (1 - np.exp(-kappa * dt))) / (4 * kappa)
        nc = np.exp(-kappa * dt) / c * x2[t - 1]
        x2[t] = c * npr.noncentral_chisquare(df, nc, size=I)
    return x2
x2 = srd_exact()

#plt.hist(x2[-1], bins=50)
#plt.xlabel(‘value’)
#plt.ylabel(‘frequency’)
#plt.grid(True)

#plt.xlabel(‘time’)
#plt.plot(x2[:, :10], lw=1.5)
#plt.ylabel(‘index level’)
#plt.grid(True)

print_statistics(x1[-1], x2[-1])

I = 250000
%time x1 = srd_euler()
%time x2 = srd_exact()

print_statistics(x1[-1], x2[-1])
x1 = 0.0; x2 = 0.0

#### HESTON MODEL ####
S0 = 100.
r = 0.05
v0 = 0.1
kappa = 3.0
theta = 0.25
sigma = 0.1
rho = 0.6
T = 1.0

corr_mat = np.zeros((2, 2))
corr_mat[0, :] = [1.0, rho]
corr_mat[1, :] = [rho, 1.0]
cho_mat = np.linalg.cholesky(corr_mat)

cho_mat

M = 50
I = 10000
ran_num = npr.standard_normal((2, M + 1, I))
dt = T / M
v = np.zeros_like(ran_num[0])
vh = np.zeros_like(v)
v[0] = v0
vh[0] = v0
for t in range(1, M + 1):
    ran = np.dot(cho_mat, ran_num[:, t, :])
    vh[t] = (vh[t - 1] + kappa * (theta - np.maximum(vh[t - 1], 0)) * dt
            + sigma * np.sqrt(np.maximum(vh[t - 1], 0)) * np.sqrt(dt)
            * ran[1])
v = np.maximum(vh, 0)

S = np.zeros_like(ran_num[0])
S[0] = S0
for t in range(1, M + 1):
    ran = np.dot(cho_mat, ran_num[:, t, :])
    S[t] = S[t - 1] * np.exp((r - 0.5 * v[t]) * dt +
    np.sqrt(v[t]) * ran[0] * np.sqrt(dt))
print_statistics(S[-1], v[-1])


'''
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 5))
ax1.hist(S[-1], bins=50)
ax1.set_xlabel(‘index level’)
ax1.set_ylabel(‘frequency’)
ax1.grid(True)
ax2.hist(v[-1], bins=50)
ax2.set_xlabel(‘volatility’)
ax2.grid(True)

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(7, 6))
ax1.plot(S[:, :10], lw=1.5)
ax1.set_ylabel(‘index level’)
ax1.grid(True)
ax2.plot(v[:, :10], lw=1.5)
ax2.set_xlabel(‘time’)
ax2.set_ylabel(‘volatility’)
ax2.grid(True)
'''

#### JUMP ####
S0 = 100.
r = 0.05
sigma = 0.2
lamb = 0.75
mu = -0.6
delta = 0.25
T = 1.0

M = 50
I = 10000
dt = T / M
rj = lamb * (np.exp(mu + 0.5 * delta ** 2) - 1)
S = np.zeros((M + 1, I))
S[0] = S0
sn1 = npr.standard_normal((M + 1, I))
sn2 = npr.standard_normal((M + 1, I))
poi = npr.poisson(lamb * dt, (M + 1, I))
for t in range(1, M + 1, 1):
    S[t] = S[t - 1] * (np.exp((r - rj - 0.5 * sigma ** 2) * dt
            + sigma * np.sqrt(dt) * sn1[t])
            + (np.exp(mu + delta * sn2[t]) - 1)
            * poi[t])
S[t] = np.maximum(S[t], 0)

plt.hist(S[-1], bins=50)
plt.xlabel(‘value’)
plt.ylabel(‘frequency’)
plt.grid(True)

plt.plot(S[:, :10], lw=1.5)
plt.xlabel(‘time’)
plt.ylabel(‘index level’)
plt.grid(True)

#### Variance Reduction Function ####
sn = npr.standard_normal(10000 / 2)
sn = np.concatenate((sn, -sn))
np.shape(sn)

sn_new = (sn - sn.mean()) / sn.std()

def gen_sn(M, I, anti_paths=True, mo_match=True):
    if anti_paths is True:
        sn = npr.standard_normal((M + 1, I / 2))
        sn = np.concatenate((sn, -sn), axis=1)
    else:
        sn = npr.standard_normal((M + 1, I))
    if mo_match is True:
        sn = (sn - sn.mean()) / sn.std()
    return sn