Simulate Common Stochastic Process

import numpy as npimport matplotlib.pyplot as pltimport numpy.random as npr####################################### simulate stock price in one step######################################S0 = 100r = 0.05sigma = 0.25T = 2.0I = 10000   # number of pathsST1 = 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 stepsdt = T/MS = np.zeros((M+1, I))S[0] = S0for 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.05kappa = 3.0theta = 0.02sigma = 0.1I = 10000M = 50dt = T/Mxh = np.zeros((M+1, I))x1 = np.zeros_like(xh)xh[0] = x0x1[0] = x0for 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 = 100r = 0.05v0 = 0.1kappa = 3.0theta = 0.25sigma = 0.1rho = 0.6T = 1.0M = 50I = 10000z_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 volatilityxh = np.zeros((M+1, I))x1 = np.zeros_like(xh)xh[0] = x0x1[0] = x0for 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 priceS = np.zeros((M+1, I))S[0] = S0for 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 = 100r = 0.05sigma = 0.2lamb = 0.75mu = -0.6delta = 0.25T = 1.0M = 50I = 10000dt = T/Mrj = lamb*(np.exp(mu+0.5*delta**2)-1)S = np.zeros((M+1, I))S[0]  = S0rn1 = 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])