Integral de Wiener (Simulación en Python)

Un caso particular de la integral de Itò, es la integral de Wiener, es el caso donde se integra una función determinista $f$ con respecto a un Movimiento Browniano. Esta integral no puede ser definida de la manera tradicional (Integral de Riemann-Stieltjes) debido a que no cumple con tener variación acotada, lo que quiere decir que tomando sobre todas las particiones $\tau_n$: $$\sup_{\tau _n}\sum_{i=1}^{n}|B_{t_i}-B_{t_{i-1}}|=\infty$$ Esta propiedad es necesario para poder usar la integral de Riemann-Stieltjes, por lo tanto se tiene que desarrollar otra técnica.

Como primer paso consideraremos una clase de funciones escalonadas o definidas por intervalos finitos, esta clase la denotaremos por $\mathcal{S}$, la cual es un subconjunto de $L^2(\mathbb{R}^+,\lambda)$ de todas las funciones positivas tales que $\int f^2d\lambda<\infty$. Todas las funciones de $\mathcal{S}$ se pueden escribir como: $$h_t=\sum_{i=1}^{k}a_i 1_{(t_{j-1},t_j]}$$ para algún $n\in \mathbb{N}$, y $0<t_1<...<t_n$ y coeficientes $a_i\in \mathbb{R}$. Nuestro interés es encontrar un operador $\eta_t:L^{2}(\mathbb{R}^+,\lambda)\rightarrow L^2(\Omega,\mathbb{P})$, como es usual en probabilidad primero se define para $\mathcal{S}\subset L^{2}(\mathbb{R}^+,\lambda)$ y luego se extiende a cualquier función.
Se define $\eta_t$ para $h\in\mathcal{S}$ de la manera obvia: $$\eta_t h=\int_{0}^{t} h_s dB_t := \sum_{i=1}^{n}a_i(B_{t_{i}\wedge t}-B_{t_{i-1}\wedge t})$$
Notamos varias cosas importantes, la primera es notar que $\eta_t$ es un proceso gaussiano ya que $B_{t_{i}\wedge t}-B_{t_{i-1}\wedge t}$ son normales independientes por lo tanto cualquier combinación lineal no nula de estas se  distribuye normal, por otro lado, como  $B_{t_{i}\wedge t}-B_{t_{i-1}\wedge t}\sim\mathcal{N}(0, t_{i}\wedge t-t_{i-1}\wedge t)$ podemos calcular que $\mathbb{E}\eta_t h=0$ para cualquier $t$. La varianza la podemos calcular: $$\mathbb{E}(\eta_t h)^2=\sum_{i=1}^{n}a_i^2Var(B_{t_{i}\wedge t}-B_{t_{i-1}\wedge t})=\sum_{i=1}^{n}a_i^2(t_{i}\wedge t-t_{i-1}\wedge t))=\int_0^{t}h^2_sds$$
 Por lo tanto $\eta_t h\sim \mathcal{N}(0,\int_0^{t}h^2_sds)$ para cada $t>0$.

Tambien es facil notar que $\eta_t(\alpha h+g)=\alpha\eta_t h+\eta_t g$, de esta manera $h\mapsto \eta_t h$ denota una isometria lineal de $\mathcal{S}$ a $L^2(\Omega,\mathbb{P})$. Se puede probar que $\mathcal{S}$ es denso en $L^2(\mathbb{R}^+,\lambda)$ por lo que se puede extender la integral gracias a la continuidad.

Usando la función definida en este post para simular trayectorias del movimiento browniano:

def BM(dt):
    T = np.arange(0,1+dt,dt)
    n = len(T)
    B = np.ones(n)*0
    for i in range(n):
        xi = sqrt(2)*(np.random.randn())/((i+0.5)*pi)
        B = B + xi*np.array([sin((i+0.5)*pi*t) for t in T])
    return T, B

Entonces podemos simular trayectorias para el proceso $\eta _t f$ para alguna $f:[0,1]\rightarrow\mathbb{R}$ que pueda ser aproximada por funciones de la clase $\mathcal{S}$, podemos usar la aproximación $$f^{n}_t=\sum_{i=1}^{n}f\left(\frac{i}{n}\right)1_{[\frac{i-1}{n},\frac{i}{n})}(t)$$
Programamos un par de funciones, la primera para calcular la integral de Riemann-Stieltjes de la manera tradicional y la segunda para la integral de Wiener simulando una trayectoria de un movimiento Browniano:

def IntS(f, dt):
    T = np.arange(0,1+dt,dt)
    n = len(T)
    y = 0
    Y = [0]
    for i in range(1,n):
        y += f(T[i])*(T[i]-T[i-1])
        Y.append(y)
    return T, Y

def IntE(f, dt):
    T, B = BM(dt)
    n = len(T)
    y = 0
    Y = [0]
    for i in range(1,n):
        y += f(T[i])*(B[i]-B[i-1])
        Y.append(y)
    return T, Y

Ambas funciones son muy similares, IntS nos servirá para poder comparar $Var(\eta_t f) =\int_0^{t}f^2_sds$ cuando ya hayamos simulado una gran cantidad de trayectorias.
Por ejemplo para la función $f(t)= \sin(\pi t)$ simulamos 15 veces con 0.001 de aproximación y graficamos el resultado:

Esta gráfica no es realmente informativa acerca de nuestras simulaciones, si queremos probar el poder  nuestras simulaciones podemos comparar las varianzas con la integral de la función al cuadrado, para esto simulamos 100 veces y vemos la aproximación: 

Y aumentando el numero de simulaciones se aproximan mas ambas curvas. Este proceso lo simplificamos en la siguiente función para graficar $n$ simulaciones con $dt$ de aporximación, para una función arbitraria $f$:

def IntegralEst(f,dt,n):
    Y = []
    DAT = []
    W = []
    fig = tools.make_subplots(rows=1, cols=2,subplot_titles=('$$\eta_t f = \int_{0}^{t}f_sdW_s$$',
                                                             '$$\int_{0}^{1}f_s^2ds=Var(\eta_t f)$$'))
    for i in range(n):
        X, y = IntE(f,dt)
        Y.append(y)
        DAT.append(go.Scatter(x=X, y=y,showlegend=False))
        fig.append_trace(DAT[i], 1, 1)
    X, V = IntS(lambda x: f(x)**2,dt)   
    for i in range(len(X)):
        W.append(np.var([y[i] for y in Y]))
    data1 = go.Scatter(x=X, y=V,  line = dict(color = 'red'),name='Varianza Real')
    data2 = go.Scatter(x=X, y=W,  line = dict(color = 'blue'),name='Varianza Simulada')
    fig.append_trace(data1, 1, 2)
    fig.append_trace(data2, 1, 2)
    fig['layout'].update(height=500, width=1000, title='Integral de Wiener para funciones deterministas')
    return py.iplot(fig, filename='IntegralEstocastica4')

La usamos para funciones interesantes:

 

Fuente:
- Kallenberg, O. (2002). Foundations of modern probability. New York: Springer.
- Karatzas, I., & Shreve, S. E. (1988). Brownian motion and stochastic calculus. New York: Springer.