Simulando Movimiento Browniano con la expansión de Karhunen-Loeve (Python)

Con la expansión de Karhunen-Loeve, representamos un proceso estocástico

$\{X(t), t\in T\}$ por medio de una secuencia de variables aleatorias independientes

$\{\xi_i, i\in \mathbb{N}\}$ . Se supone que el proceso tiene esperanza

$0$ y segundo momento finito. Consideramos la función de covarianzas

$\mathbb{E}X(t)X(s)=K(t,s)$ . Esta función es importante ya que nos generará un espacio dado de cierta ecuación integral. Sabiendo esto podemos enunciar el siguiente teorema:

Teorema:

Sea

$\{X(t), t\in T\}$ un proceso estocástico con media 0 y segundo momento finito. Con función de covarianzas continua

$K(t,s)$ .

1. Supongamos que

$\lambda _k$ y

$\varphi_k$ satisfacen la siguiente ecuación integral:

$\int_T K(s,t)\varphi_k (t)dt = \lambda _k \varphi_k (s)$ Y tomamos

$\xi_k = \frac{1}{\sqrt{\lambda_k}}\int_{T} X(t)\varphi_k (t) dt$ Entonces:

$X(t)=\lim_{n\rightarrow \infty} \sum_{k=1}^{n}\sqrt{\lambda _k} \xi_k \varphi_k (t)$ uniformemente de manera que

$\mathbb{E}\left( X(t) -\sum_{k=1}^{n}\sqrt{\lambda _k} \xi_k \varphi_k (t)\right) ^2\rightarrow 0$

2. Inversamente si

$X(t)=\sum_{k=1}^{\infty}\sqrt{\lambda _k} \xi_k \varphi_k (t)$ donde

$\xi_{k}$ son variables aleatorias independientes con media 0 y varianza 1, entonces:

$\int_T K(s,t)\varphi_k (t)dt = \lambda _k \varphi_k (s)$

De esta manera podríamos tomar las sumas parciales

$X^{(n)}(t)=\sum_{k=1}^{n}\sqrt{\lambda _k} \xi_k \varphi_k (t)$ para aproximar nuestro proceso

$X(t)$ con un valor suficientemente grande de

$n$ , nuestro problema reside en nuestra capacidad para encontrar soluciones analíticas de la ecuación integral

$\int_T K(s,t)\varphi_k (t)dt = \lambda _k \varphi_k (s)$ lo cual en la mayoría de los casos no nos sera posible y tendremos que aproximar nuestras soluciones de

$\varphi_k$ y

$\lambda_k$ de manera computacional.

Nuestro objetivo es poder aproximar el movimiento browniano con la expansión de Karhunen-Loeve. El movimiento Browniano es un proceso gausiano que empieza en

$B_0=0$ y tiene función de covarianzas

$K(t,s)=\min\{t,s\}$ con

$s,t\in T$ por lo que debemos solucionar la ecuación integral

$\mathcal {K}f =\int_T \min\{s,t\}f(t)dt=\lambda f(s)$ Notamos que

$-\frac{d^2}{ds^2}\mathcal{K}f=f$ entonces es equivalente a solucionar la ecuación diferencial

$f''+\frac{1}{\lambda}f=0$ con

$f(0)=f'(1)=0$ , cuya solución es sencilla de conocer:

$\varphi_k(t) =\sqrt{2}\sin\left[ (k-\frac{1}{2})\pi t\right]$

$\lambda _k = \frac{4}{\pi ^2 (2k-1)^2}$ Así:

$B(t)=\sqrt{2}\sum_{k=1}^{\infty}\frac{2}{\pi (2k-1)}\sin\left[\left( k-\frac{1}{2}\right)\pi t\right]\xi_k$

En el caso de los puentes brownianos, que son procesos con la propiedad que

$W_0=W_1=0$ y solo están definidos en el intervalo

$[0,1]$ , tienen una matriz de covarianzas

$K(t,s)=\min\{t,s\} -st$ por lo que nuestro problema es muy similar, en este caso tambien logramos obtener una solución analitica dada por

$W(t)= \sqrt{2}\sum_{k=1}^{\infty}\frac{1}{k\pi}\sin(k\pi t)\xi_k$

Podemos ahora aplicar nuestros resultados para simular trayectorias de ambos procesos en el intervalo

$[0,1]$ , para esto utilizaremos Python y las expresiones

$X^{(n)}(t)$ explicitamente, con

$n=1000$ , se agrega como un extra el caso del Brownian Sheet donde en lugar de tener un tiempo se tiene 2, para este caso se uso la expansión de Karhunen-Loeve de varias variables, pero en escencia es un resultado muy similar.

Importamos todo lo necesario en este caso usaremos ployly para graficar nuestras simulaciones:

from __future__ import division
from math import *
import pandas as pd
import numpy as np
import plotly.graph_objs as go
import plotly.plotly as py
from plotly.tools import FigureFactory as FF
import plotly.plotly as py

En el caso del Movimiento Browniano tenemos la expansión:

$B(t)=\sqrt{2}\sum_{k=1}^{\infty}\frac{2}{\pi (2k-1)}\sin\left[\left( k-\frac{1}{2}\right)\pi t\right]\xi_k$ para la cual escribimos la siguiente función de

$dt$ :

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

Simulamos varias veces y graficamos el resultado

W = []
DAT = []
for k in range(8):
    X, Y = BM(0.001)
    W.append(Y)
    DAT.append(go.Scatter(x=X, y=Y,  line = dict(color = '#ff{}000'.format(k)),name="B{}".format(k)))    

Layout=go.Layout(title='Movimiento Browniano $B_t(\omega)$')
fig = go.Figure(data=DAT, layout=Layout)
py.iplot(fig,filename='MBM')

Hacemos lo mismo en el caso del Puente Browniano; ;

def BB(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+1)*pi)
        B = B + xi*np.array([sin((i+1)*pi*t) for t in T])
    return T, B

Y simulamos una trayectoria para probar:

W = []
DAT = []
for k in range(8):
    X, Y = BB(0.001)
    W.append(Y)
    DAT.append(go.Scatter(x=X, y=Y,  line = dict(color = '#ff{}580'.format(k)),name="W{}".format(k)))    

Layout=go.Layout(title='Puente Browniano $W_t(\omega)$')
fig = go.Figure(data=DAT, layout=Layout)
py.iplot(fig,filename='PB')

Finalmente como extra simularemos un Brownian Sheet, que es un proceso estocástico con tiempo bivariado cuya función de covarianzas es:

$\mathbb{E}X(t_1,s_1)X(t_2,s_2)=\min\{s_1,s_2\}\min\{t_1,t_2\}$ , se cuenta con una versión más general de la expansión de Karhunen-Loeve para más variables, esta se calcula de manera similar resolviendo ecuaciones integrales, la expansión del Brownian Sheet es la siguiente:

$W(t,s)=\sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{2}{\pi^2 (2i-1)(2j-1)}\sin \left[\left( i-\frac{1}{2}\right) \pi t\right] \sin \left[\left( j-\frac{1}{2}\right) \pi s\right] \xi_{ij}$
Programamos un par de funcíones, para calcular

$W(t,s)$ puntualmente y luego para una gradilla dada por

$dt$

def BSp(x,y,n, Nij):
    z = 0
    for i in range(n):
        for j in range(n):
            z += Nij[i][j]*(2/(pi**2*(i-0.5)*(j-0.5)))*sin((i-0.5)*pi*x)*sin((j-0.5)*pi*y)
    return z

def BS(dt):
    T = np.arange(0,1+dt,dt)
    Z = pd.DataFrame(columns=T)
    n = len(T)
    Nij = np.random.rand(n,n)
    for x in T:
        Z[x] = [BSp(x,y,n, Nij) for y in T]
    return  Z

En este caso se obtiene una grafica en 3-D:

Z = BS(0.01)
data = [
    go.Surface(
        z=Z.as_matrix()
    )
]
layout = go.Layout(
    title='Brownian Sheet',
    autosize=False,
    width=800,
    height=500,
    margin=dict(
        l=65,
        r=50,
        b=65,
        t=90
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='Brownian Sheet')

Pueden consultar el código en https://github.com/aleespa/Karhunen-Loeve-Expansion.git
Fuentes:
1. Karhunen-Loeve Expansions and their Applications, Limin Wang, The London School of Economics and Political Science, 2008.
2. Quasi-Monte Carlo simulation of Brownian sheet with application to option pricing, Xinyu Song and Yazhen Wang, Statistical Theory and Related Fields, pag 82-91, 2017, https://doi.org/10.1080/24754269.2017.1332965

Ciencia Actuarial

Buscar este blog

Simulando Movimiento Browniano con la expansión de Karhunen-Loeve (Python)

Comentarios

Publicar un comentario