Introducción

Las ecuaciones diferenciales estocásticas son de suma importancia en áreas como biología, química o Finanzas.

Estas ecuaciones tienen como base el movimiento browniano, el primero en describir matemáticamente al movimiento browniano fue N. Thiele en un documento sobre mínimos cuadrados. Seguido independientemente de L.Bachelier (1900) en su tesis doctoral Théorie de la spéculation en la que presenta un análisis estocástico de acciones y opciones financieras.

A mediados del siglo XX Kiyoshi Itô extiende los métodos del cálculo a procesos estocásticos tales como el MB, teoría llamada como cálculo de Itô que tiene importantes aplicaciones en Finanzas y en ecuaciones diferenciales estocásticas. Su concepto central es la integral estocástica de Itô, que es una generalización de la integral de Riemann-Stieltjes.Itô construyó una ecuación diferencial estocástica (EDE) de la forma

\[ dX_t = a(X_t)dt + b(X_t)dW_t \]

A Kiyoshi Itô se le conoce como el padre de la integración estocástica. Tiempo después, el físico Ruslan Stratonovich construyó una integral alternativa a la de Itô, conocida como Integral de Stratonovich (o de Fisk-Stratonovich), integrales que son, en cierto sentido, fáciles de manipular. Su trabajo llega a manos de Kiyoshi Itô, y éste decide perfeccionar esta alternativa. Es desde ese entonces donde se sientan las bases de la teoría del cálculo estocástico y el de las ecuaciones diferenciales estocásticas, presentes en las aplicaciones que involucran sistemas con el in‡ujo de perturbaciones aleatorias: finanzas, dinámica poblacional, etc.

Movimiento Browniano

El movimiento browniano lleva su nombre gracias al botánico Robert Brown quien, en 1820, fue el primero observar el movimiento irregular de granos de polen inmersos en agua. En 1905 Einstein, a pesar de ignorar el descubrimiento, predijo su existencia de forma puramente teórica. Cinco años antes, L.Bachelier lo había empleado para modelar precios de acciones en el tiempo, a través de la fórmula :

\[ S(t) = S(0) + \sigma W(t) \]

Definición

$W(0) = 0 \ \ \ \ c.s$
Para $ 0 < s < t < T $ la variable aleatoria dada por el incremento\ $W(t)-W(s)$ tiene una distribucion normal con media cero y varianza $t-s$ o lo que es lo mismo que: \[ W(t)-W(s) \backsim \sqrt{t-s}N(0,1) \]
Para $s,t,u,v$ tal que $0 \leq s < t < u < v \leq T$ los incrementos $W(t)-W(s) $ y $ W(v)-W(u) $ son independientes.

Aproximación Por Euler-Maruyama.

Para presentar las principales cuestiones relativas a la aproximación discreta en el tiempo de los procesos Ito, se examina el esquema estocástico de Euler con algún detalle en este capítulo. Consideraró algunos ejemplos de problemas típicos que pueden manejarse mediante la simulación de trayectorias discretas de tiempo aproximado. Además, se darán definiciones generales para discretizaciones de tiempo y aproximaciones discretas de tiempo, y se introducirán los criterios de convergencia fuertes y débiles para aproximaciones discretas de tiempo.

Una de las aproximaciones discretas de tiempo más simples de un proceso Ito es la aproximación Euler, o la aproximación Euler-Maruyama, como a veces se la llama. Consideraremos un proceso Ito $X=\{ X_t, t_0\leq t\leq T \}$ que satisface la ecuación diferencial estocástica escalar.

\[ dX_t = a(t,X_t)dt + b(t,X_t)dW_t \]

sobre $t_0\leq t\leq T$ con la condición inicial

\[ X_{t_0} = X_0 \]

Para una discretización dada $t_0 = \tau_0<...<\tau_N = T$ del intervalo $[t_0,T]$ una aproximación de Euler es un proceso estocástico de tiempo continuo $Y = \{Y(t),t_o<t<T\}$ satisfaciendo el esquema iterativo.

\[ Y_{n+1} = Y_n + a(\tau_n,Y_n)(\tau_{n+1}-\tau_n) + b(\tau_n,Y_n)(W_{\tau_{n+1}}-W_{\tau_n})\]

Para $n= 0,1,2...N-1$ con valor inicial \[ Y_0 = X_0\]

Para el valor de la aproximación en el tiempo de discretización $\tau_n$.Escribiremos

\[ \Delta _n = \tau_{n+1} - \tau_n \]

para el enésimo incremento de tiempo y llamaremos

\[ \delta = max_n \{ \Delta _n \}\]

el paso de tiempo mínimo. Para gran parte de este capítulo consideraremos tiempos de discretización equidistantes

\[ \tau_n = t_o + n \delta \]

Cuando el coeficiente de difusión es idénticamente cero, es decir, cuando $b \equiv 0$, el esquema iterativo estocástico se reduce al esquema determinista de Euler para la ecuación diferencial ordinaria

\[ \frac{dx}{dt} = a(t,x) \]

La secuencia $\{Yn, n = 0, 1, ..., N\}$ de valores de la aproximación de Euler en los instantes de la discretización de tiempo $()_= { _n, n = 0,1,…,N } $ se puede calcular de manera similar a las del caso determinista. La principal diferencia es que ahora necesitamos generar los incrementos aleatorios

\[ \Delta W_n = W_{\tau_{n+1}}-W_{\tau_n} \]

Para $n = 0, 1, ..., N -1$, del proceso Wiener $W = \{Wt, t \geq 0\}$ sabemos que estos incrementos son variables aleatorias gaussianas independientes con media

\[ \mathbb{E}(\Delta W_n) = 0\]

y varianza

\[ \mathbb{E}(\Delta W_n)^2 = \Delta_n \]

Podemos usar una secuencia independiente de números aleatorios seudo-gaussianos.

Aqui se expone para el caso particular:

\[ S_t = 254.1221 \cdot e ^{ -0.003697016 \cdot t + 0.04323157 W(t) } \]

Euler_Maruyama = function(mu,sigma,Xo,n){
  # Ecuacion explicita 
     # n es el numero de nodos en la aproximacion 
  mu=mu ; sigma=sigma; Xzero=Xo  # Parametros y condicion inicial. 
  
  T=1 # Fin del intervalo 
  N=2^9 # Buscamos una particion fina para la simulacion. 
  dt=1/N # numero de pasos queremos muchos porque este es continua  
  dw=sqrt(dt)*rnorm(N) # Derivada del movimiento Browniano. 
  w=cumsum(dw) # Aqui acumulamos. 
  t=seq(dt,T,by=dt)
  t1=seq(0,T,by=dt)
  # Esta es la solucion explicita obtenida en el capitulo 9. 
  Xtrue=Xzero*exp((mu-0.5*sigma^2)*t+sigma*w)
  
  # Aqui comienza la aproximacion E-M
  
  R=N/n ; Dt=R*dt ; L=N/R
  Xem=0
  Xtemp=Xzero
  for (j in 1:L) {
    Winc=sum(dw[(R*(j-1)+1):(R*j)])
    Xtemp=Xtemp + Dt*mu*Xtemp + sigma*Xtemp*Winc
    Xem[j]=Xtemp
  } 
  t2=seq(0,T,by=Dt)
  # Error de aproximacion. 
    # Proximanete. 
  real = c(Xzero,Xtrue)
  indices = c()
  
  for(i in 1:length(t2)){
    indices = append(indices,which(t1==t2[i]))
  }
  
  # Aqui comienzan los ploteos. 
  plot(t1,c(Xzero,Xtrue),type="l",col="blue",xlab="t",ylab="X")
  grid()
  lines(t2,c(Xzero,Xem),col="red",lty=2)
  points(t2,c(Xzero,Xem),pch=20)
  
 
 
}

# En esta parte llamamos a la funcion que creamos. 

Euler_Maruyama(1.5,1,1,4)

Error en Euler Maruyama

# M es el número de caminos
# N es el número de pasos
# k es la semilla 
E_Maruyama = function(mu,sigma,Xo,n,M,k){
  set.seed(k)
  alma = c()
  
  for(i in 1:M){
  
  # Ecuacion explicita 
  mu=mu ; sigma=sigma; Xzero=Xo  # Parametros y condicion inicial. 
  T=1 # Fin del intervalo 
  N=2^9 # Buscamos una particion fina para la simulacion. 
  dt=1/N # numero de pasos  
  dw=sqrt(dt)*rnorm(N) # Browniano. 
  w=cumsum(dw) # Aqui acumulamos. 
  t=seq(dt,T,by=dt)
  t1=seq(0,T,by=dt)
  # Esta es la solucion explicita obtenida en el capitulo 9. 
  Xtrue=Xzero*exp((mu-0.5*sigma^2)*t+sigma*w)
  
  # Aqui comienza la aproximacion E-M
    # L aqui 
  
  R=N/n ; Dt=R*dt ; L=N/R
  Xem=0
  Xtemp=Xzero
  for (j in 1:L) {
    Winc=sum(dw[(R*(j-1)+1):(R*j)])
    Xtemp=Xtemp + Dt*mu*Xtemp + sigma*Xtemp*Winc
    Xem[j]=Xtemp
  } 
  t2=seq(0,T,by=Dt)
  # Error de aproximacion. 
  # Proximanete. 
  real = c(Xzero,Xtrue)
  indices = c()
  
  for(i in 1:length(t2)){
    indices = append(indices,which(t1==t2[i]))
  }
  emerr=abs(Xem[length(Xem)]-Xtrue[length(Xtrue)])
  alma = append(alma,emerr)
  }
  print(mean(alma))
}


# En esta parte llamamos a la funcion que creamos. 

# M es el número de caminos
# N es el número de pasos
# k es la semilla 
# E_Maruyama = function(mu,sigma,Xo,n,M,k)

E_Maruyama(1.5,1,1,64,25,123)

## [1] 0.3605027

Intervalos para el error de aproximación.

La propagación de los errores iniciales y los errores de redondeo deben mantenerse bajo control. Este problema también surge para las ecuaciones diferenciales estocásticas rígidas, que se caracterizan en el caso lineal de la dimensión:

\[ dZ_t = AZ_tdt + \sum_{k=1}^{m}B^{k} Z_t \circ dW_t \]

La propagación de un error inicial quedará así limitada en cualquier intervalo acotado para un esquema numéricamente estable. Hacemos hincapié en que el criterio de estabilidad numérica se aplica solo a los tamaños de paso $\delta > 0$ que son menores que algún valor crítico $\Delta_0$ que generalmente dependerá del intervalo de tiempo $[t_0,T]$ y la ecuación diferencial en consideración. Este valor crítico puede ser extremadamente pequeño en algunos casos.y la ecuación diferencial en consideración. Este valor crítico puede ser extremadamente pequeño en algunos casos.

Podemos ver que a medida que el intervalo de tiempo $[t_0, T]$ se vuelve relativamente grande, el error propagado de un esquema numéricamente estable, que en teoría todavía está bajo control, puede, de hecho, llegar a ser tan irrealmente grande como para hacer la aproximación inútil para lgunos propósitos prácticos. Por ejemplo, al simular los primeros tiempos de salida, no conocemos el intervalo de tiempo apropiado de antemano y, por lo tanto, debemos permitir un intervalo arbitrariamente grande.

En este caso se calcula el error de aproximacion por trayectorias mediante el siguiente codigo.

# M es el número de caminos.
# N es el número de pasos.
# k es la semilla. 
# a es el intervalo de confianza.
E_Mar_I = function(mu,sigma,Xo,n,M,k,a,l){
  set.seed(k)
  alma = c()
  for(i in 1:M){
  
  # Ecuacion explicita 
  mu=mu ; sigma=sigma; Xzero=Xo  # Parametros y condicion inicial. 
  T=1 # Fin del intervalo 
  N=2^9 # Buscamos una particion fina para la simulacion. 
  dt=1/N # numero de pasos  
  dw=sqrt(dt)*rnorm(N) # Browniano. 
  w=cumsum(dw) # Aqui acumulamos. 
  t=seq(dt,T,by=dt)
  t1=seq(0,T,by=dt)
  # Esta es la solucion explicita obtenida en el capitulo 9. 
  Xtrue=Xzero*exp((mu-0.5*sigma^2)*t+sigma*w)
  
  # Aqui comienza la aproximacion E-M
  
  R=N/n ; Dt=R*dt ; L=N/R
  Xem=0
  Xtemp=Xzero
  for (j in 1:L) {
    Winc=sum(dw[(R*(j-1)+1):(R*j)])
    Xtemp=Xtemp + Dt*mu*Xtemp + sigma*Xtemp*Winc
    Xem[j]=Xtemp
  } 
  t2=seq(0,T,by=Dt)
  # Error de aproximacion. 
  # Proximanete. 
  real = c(Xzero,Xtrue)
  indices = c()
  
  for(i in 1:length(t2)){
    indices = append(indices,which(t1==t2[i]))
  }
  
  
  emerr=abs(Xem[length(Xem)]-Xtrue[length(Xtrue)])
  alma = append(alma,emerr)
  }
 e = mean(alma)
 v_e = var(alma)
 delta = qt(1-a,l-1)*(v_e/l)**0.5
 a = e - delta
 b = e + delta

 c = data.frame(a,b)
 print(c)
}

Intervalos de Confianza.

# Creacion de los distintos intervalos 
   
I1 = E_Mar_I(1.5,0.1,1,16,100,312,0.1,10)

##           a         b
## 1 0.2507475 0.3035163

I2 = E_Mar_I(1.5,0.1,1,16,100,312,0.1,20)

##           a         b
## 1 0.2592213 0.2950425

I3 = E_Mar_I(1.5,0.1,1,16,100,312,0.1,40)

##           a         b
## 1 0.2646969 0.2895668

I4 = E_Mar_I(1.5,0.1,1,16,100,312,0.1,100)

##           a         b
## 1 0.2693486 0.2849151

Invervalos de confianza variando delta

Se toman los intervalos de confianza considerando la tranformacion logaritmica para vizualizar mejor los cambios que ocurren.

x = c(log2(0.25),log2(0.125),log2(0.0625),log2(0.03125))
y = c( log2(0.8561853), log2(0.4958418) , log2(0.2685741), log2(0.1411855)  )

plot(x,y, main = "Intervalos de confianza variando log-Delta", xlab = "log2(Delta)", ylab = "log2(Eps)",type = "l",col="purple")
grid()

points(x,y,pch=19)

Payoff usando Euler-Maruyama

Ricardo Rodriguez Ovilla

19/6/2020