.py

import matplotlib.pyplot as plt
import seaborn
from tqdm import tqdm
import numpy as np
plt.ion()
T_mat = 5  # years
vol0 = .3
var_mean = vol0**2
var_var=.2**2
kappa = .0
n_steps = 50
h = T_mat / n_steps

# def dz(dt):
#     return np.random.choice([-1, 1]) * np.sqrt(dt)

def dv(vt, h, var_var):
    return kappa*(var_mean - vt)*h + var_var * np.sqrt(vt) * np.sqrt(h) * np.random.normal()

def simulate_vols(v0, h, var_var, n_steps):
    '''
    Simulate volatilities according to that kappa model.
    returns the volatilities
    '''
    L = [v0]
    for t in range(n_steps):
        d = dv(v0, h, var_var)
        v0 = v0
        if v0 + d > 0:
            v0 = v0 + d
        else:
            while True:
                N = simulate_vols(v0, h**2/100, var_var, n_steps=100)[-1]
                if not np.isnan(N) and v0 + N > 0:
                    v0 = N
                    break

        L.append(v0)
    return np.sqrt(np.array(L))


X = []
ax = plt.gca()
ax.cla()

n_runs = 100
for col, var_var in zip('bg', [.2**2, .5**2]):
    X = []
    for t in tqdm(range(n_runs)):
        L = simulate_vols(var_mean, h, var_var, n_runs)
        X.append(L)
        label = None if t != n_runs-1 else r'$\sigma_v =$ ' + f'{np.sqrt(var_var):.2f}^2'
        ax.plot(range(n_runs+1), L, col, alpha=.07, label=label)

    means = np.array(X).mean(axis=1)
    ax.plot(means, col, alpha=1, label='mean with ' + r'$\sigma_v = $' + f'{np.sqrt(var_var):.2f}^2')

ax.set_yticklabels([f'{np.sqrt(e):.3f}^2' for e in ax.get_yticks()])
ax.legend()
ax.xaxis.grid(0), ax.yaxis.grid(0)
seaborn.despine()