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()

def srd_eulerale():
    xh = np.zeros((M+1, I))
    x1 = 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 * sqrt(np.maximum(xh[t-1], 0)) * sqrt(dt)
                * npr.standard_normal(I))
    x1 = np.maximum(xh, 0)
    return x1
x1 = srd_eulerale()


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

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