

In [1]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Brownian Motion



In [2]:

    
T = 1.0
N = 1000
dt = T/N
dW = zeros(N)
W = zeros(N)

dW = sqrt(dt) * random.normal(size=N)
W = dW.cumsum()
plot(W)









    Out[2]:





[<matplotlib.lines.Line2D at 0x7f6da8fae990>]



In [3]:

    
T = 1.0
N = 1000
M = 10
dt = T/N
dW = zeros((N, M))
W = zeros((N, M))

random.seed(100)

for i in range(M):
    dW = sqrt(dt) * random.normal(size=N)
    W[:,i] = dW.cumsum()
    plot(W[:,i])

Function of W(t)



In [4]:

    
T = 1.0
N = 100
M = 500
dt = T/N
dW = zeros((N, M))
W = zeros((N, M))

t = linspace(0,T,N)
U = W.copy()
for i in range(M):
    dW = sqrt(dt) * random.normal(size=N)
    W[:,i] = dW.cumsum()
    U[:,i] = exp(t + 0.5 * W[:,i])

for i in random.choice(arange(M),size=5):
    plot(t, U[:,i])
plot(t, U.mean(axis=1), lw=5)









    Out[4]:





[<matplotlib.lines.Line2D at 0x7f6da8f89150>]



In [5]:

    
plot(t, U.mean(axis=1))
plot(t, exp(9.0/8.0 * t))
print norm((U.mean(axis=1)- exp(9.0/8.0 * t)),ord=inf)









    



0.0969578450577



In [6]:

    
plot(t, abs(U.mean(axis=1) - exp(9.0/8.0 * t)))









    Out[6]:





[<matplotlib.lines.Line2D at 0x7f6da842e450>]

Ito and Stratonovich Integrals



In [7]:

    
T = 1.0
N = 1000
dt = T/N
dW = zeros(N)
W = zeros(N)

dW = sqrt(dt) * random.normal(size=N)
W = dW.cumsum()
W = insert(W,0,0)

ito = sum(W[0:-1] * dW)

strat = sum((0.5 * (W[0:-1] + W[1:]) + 0.5 * sqrt(dt) * random.normal(size=N)) * dW)

itoerr = abs(ito - 0.5*(W[-1]**2-T))
straterr = abs(strat - 0.5 * W[-1]**2)

print "Ito:", ito
print "Strat:", strat
print "Error of Ito:",itoerr
print "Error of Strat:", straterr









    



Ito: 0.284702555399
Strat: 0.79074137579
Error of Ito: 0.00318051864766
Error of Strat: 0.00921933903828

Euler Maruyama para la Ec. Black-Scholes.



In [8]:

    
gamma=2
mu=1
Xzero=1

T=1
N=2**7
dt = float(T)/N
t = linspace(0,T,N+1)

dW = sqrt(dt) * random.randn(1,N)
W = cumsum(dW)

Xtrue = Xzero * exp((gamma-0.5*mu**2)*t[1:] + mu* W)
Xtrue = insert(Xtrue,0,Xzero)

figure(figsize=(8,8))
plot(t,Xtrue)

R=4
Dt = R*dt
L = float(N)/R
Xem = zeros(L+1) # Euler Maruyama Solution 
Xem[0] = Xzero

for j in range(1,int(L)+1):
    Winc = sum(dW[0][range(R*(j-1), R*j)]) # subsample
    Xem[j] = Xem[j-1] + Dt * gamma * Xem[j-1] + mu * Xem[j-1] * Winc

emerr=np.abs(Xem[-1]-Xtrue[-1])
print "Error at endpoint: ", emerr

plot(np.linspace(0,T,L+1),Xem,'r--*')









    



Error at endpoint:  0.0321733814685






    



/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:22: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future






    Out[8]:





[<matplotlib.lines.Line2D at 0x7f6da83f5c10>]