3349.sagews

import numpy as np
iBar = 5
sigma = 1
T = 252 * 2
dt = 1/252

ep = np.zeros((T, iBar))
z = np.zeros((T, iBar))
dz = np.zeros((T, iBar))
mu = np.zeros((T, iBar)) 

for c in range(0,iBar):
    ep[:,c] = np.random.normal(0, sqrt(dt) * 1, T)
    z[0,c]  = 0.05 * (c / iBar)

for t in range(1,T):
    for c in range(iBar,1):
        if (z[t-1,c] > z[t-1,c-1]):
            z1 = z[t-1,c-1]
            z2 = z[t-1,c]
            z[t-1,c-1] = z2
            z[t-1,c] = z1      
    for c1 in range(0,iBar):        
        for c2 in range(0,iBar):
            if c2 != c1:
                mu[t,c1] = mu[t,c1] +  1 / ( z[t-1,c1] - z[t-1,c2] )
    for c in range(0,iBar):        
        dz[t,c] = (2 * sigma^2 * mu[t,c] / iBar - z[t-1,c]) * dt + (2 * sigma / sqrt(iBar)) * ep[t,c]
        z[t,c] = z[t-1,c] + dz[t,c]

p1 = points([((t/T - 0.5)*2, z[t,0] ) for t in range(1/dt, T)], alpha = 0.50, rgbcolor = (1,0,0))
p2 = points([((t/T - 0.5)*2, z[t,1] ) for t in range(1/dt, T)], alpha = 0.50, rgbcolor = (0,1,0))
p3 = points([((t/T - 0.5)*2, z[t,2] ) for t in range(1/dt, T)], alpha = 0.50, rgbcolor = (0,0,1))
p4 = points([((t/T - 0.5)*2, z[t,3] ) for t in range(1/dt, T)], alpha = 0.50, rgbcolor = (0,1,1))
p5 = points([((t/T - 0.5)*2, z[t,4] ) for t in range(1/dt, T)], alpha = 0.50, rgbcolor = (1,0,1))

eqn = text("$dx_{i,t} = \\left( \\sum_{i \\neq j} \\frac{2 \\cdot \\sigma^2}{\\bar{i} \\cdot \\left( x_{i,t} - x_{j,t} \\right)} - x_{i,t} \\right) \\cdot dt + \\left( \\frac{2 \\cdot \\sigma}{\\sqrt{\\bar{i}}} \\right) \\cdot dB_{i,t}$", (0.60, 2.60), rgbcolor=(0,0,0), fontsize = 20)
yLab = text("$x_{i,t}$", (-0.025, 0.5), rgbcolor=(0,0,0), fontsize = 20, rotation = 90)
param = text("$\\sigma = 1, \\bar{i} = 5$", (0.25, 0.25), rgbcolor=(0,0,0), fontsize = 16)

show(p1 + p2 + p3 + p4 + p5 + eqn + yLab + param, figsize = [9,5], ymax = 3.0, ymin = -2.5, xmin = -0.03)