# coding: utf-8

from matplotlib import pyplot as plt
import numpy as np

def levyStep(oldBt, T):
    """
    Given n+1 points oldBt[i] of a Brownian motion at equidistant times 
    `t_i = i/n T`, i=0...n, simulate the values of that Borwnian motion at 
    times `(t_i + t_(i+1))/2` and return all points on the finer time grid.
    """
    assert len(oldBt) >= 2
    n = len(oldBt) - 1
    left = oldBt[:-1]
    right = oldBt[1:]
    # Sei dt = T/n, Bl:= B_(t_i), Br:= B_(t_(i+1)). 
    # Also Br-Bl ~ N(0, dt) und Br-Bl ist unabhängig von Bl.
    #
    # Bmitte:= (Bl+Br)/2 + sZ  mit unabh. sZ ~ N(0,s^2)
    # Bmitte = B_(t_i + dt/2), also Bmitte-Bl ~ N(0, dt/2)
    # Somit ist
    #    Bmitte-Bl = (Br-Bl)/2  + sZ           ~ N(0, dt/4 + s^2)
    #                \____/      \_/
    #                 ~N(0,dt)     ~N(0,s^2)
    #                \_______/
    #                 ~N(0, dt/4)
    #
    # Also muss s^2 = dt/4 = T/(4n) sein.
    scaledZ = np.random.normal(loc=0, scale=np.sqrt(T/(4.0*n)), size=n)
    newBt = np.zeros(2*n + 1)
    newBt[0::2] = oldBt
    newBt[1::2] = (left + right) / 2.0 + scaledZ
    return newBt

def initBrownianMotion(n, T):
    assert n > 0
    dBt = np.random.normal(loc=0, scale=np.sqrt(float(T)/n), size=(n+1))
    dBt[0] = 0.0
    return np.cumsum(dBt)

fig = plt.figure()
fig.canvas.set_window_title("Lévy-Konstruktion der brownschen Bewegung")
#plt.axes().set_aspect('equal','box')
T = 1.0
Bt = initBrownianMotion(1, T)

#ts = np.linspace(0,T,100)
#xs = np.sqrt(ts)
#plt.plot(ts, xs, 'r-', ts, -xs, 'r-')

for step in range(0,10):
    plt.title("Schritt %d (warte auf Mausklick...)" % step)
    plt.plot(np.linspace(0,T,len(Bt)), Bt, linewidth=0.5)
    plt.tight_layout()
    plt.show()
    plt.waitforbuttonpress()

    Bt = levyStep(Bt, T)

plt.title("Lévy-Konstruktion der brownschen Bewegung")
plt.plot(np.linspace(0,T,len(Bt)), Bt, 'k-', linewidth=1)
plt.tight_layout()
plt.show()