Simulate Common Stochastic Process

import numpy as np
import matplotlib.pyplot as plt
import numpy.random as npr

######################################
# simulate stock price in one step
######################################
S0 = 100
r = 0.05
sigma = 0.25
T = 2.0
I = 10000   # number of paths

ST1 = S0*np.exp((r-0.5*sigma**2)*T+sigma*np.sqrt(T)*npr.standard_normal(I))
#plt.hist(ST1, bins=200)

######################################
# simulate stock price paths
######################################
M = 50     # number of steps
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))

#plt.hist(S[-1], bins=200)

#######################################
# simulate volatility
#######################################
x0 = 0.05
kappa = 3.0
theta = 0.02
sigma = 0.1
I = 10000
M = 50
dt = T/M

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)
#for i in range(100):
#    plt.plot(x1[:,i], lw=1.5)
    
########################################    
# simulation Heston volatility model
########################################
S0 = 100
r = 0.05
v0 = 0.1
kappa = 3.0
theta = 0.25
sigma = 0.1
rho = 0.6
T = 1.0

M = 50
I = 10000

z_ret = npr.standard_normal((M+1, I))
z_v = npr.standard_normal((M+1, I))
z_vcorr = rho*z_ret + np.sqrt(1-rho**2)*z_v

# simulate volatility
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)*z_vcorr[t])
x1 = np.maximum(xh, 0)

# simulate stock price
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*x1[t])*dt + np.sqrt(x1[t]*dt)*z_ret[t])
    
#for i in range(100):
#    plt.plot(x1[:,i])
    
###########################################
# simulate Merton jump diffusion model
###########################################
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
rn1 = npr.standard_normal((M+1, I))
rn2 = npr.standard_normal((M+1, I))
poi = npr.poisson(lamb*dt, (M+1, I))

for t in range(1, M+1):
    S[t] = S[t-1]*(np.exp((r-rj-0.5*sigma**2)*dt+sigma*np.sqrt(dt)*rn1[t]) + (np.exp(mu+delta*rn2[t])-1)*poi[t])

for i in range(500):
    plt.plot(S[:,i])