Diferenciación numérica

Tanto en ciencias puras como negocios e ingeniería, e incluso en otras áreas sociales como la economía, dinámica de poblaciones, etc., se estudian sistemas que cambian, y por tanto existe el interés de entender la manera en que ocurren esos cambios.

En general, los cambios de una funcion $f$ son ocasionados por los cambios en la variable, $x$; ambos están relacionados mediante la tasa de cambio \[\delta f(x_0)=\frac{f(x_0+\Delta x)-f(x_0)}{\Delta x}\] donde el limite cuando $\Delta x \rightarrow 0$ corresponde a la derivada de la funcion en $x_0$: $f^{\prime}(x_0)$. Muchas veces no se conoce la regla de la funcion, $f(x)$, pero solamente una tabla de datos, la derivada se puede aproximar \[f^{\prime}(x_i)\approx\frac{f(x_{i+1})-f(x_i)}{\Delta x}\] donde la cercania entre los puntos es $\Delta x$ pequenyo.

Se puede hacer una mejor aproximacion, utilizando una serie de Taylor y despejando la derivada \[f^{\prime}(x_i)=\frac{f(x_{i+1})-f(x_i)}{\Delta x}-\frac{1}{2!}f^{\prime \prime}(x_i)\Delta x-\ldots\] lo cual nos dice que la primera aproximacion tiene un error \[\varepsilon =\frac{1}{2!}f^{\prime \prime}(x_i)\Delta x+\ldots\] Si el iintervalo $\Delta x\ll1$, el prmer termino de la serie dominara su comportamento, se dice que $\varepsilon \sim \mathcal{O}(\Delta x)$.

ejemplo

Dada la funcion

x <- seq(0.2,0.7, by=0.1); fx <- c(0.45, 0.57,0.71, 0.88,1.08,1.27)
(cbind(x,fx))

##        x   fx
## [1,] 0.2 0.45
## [2,] 0.3 0.57
## [3,] 0.4 0.71
## [4,] 0.5 0.88
## [5,] 0.6 1.08
## [6,] 0.7 1.27

Calcule las tasas de cambio, usando (diferencias hacia adelante) \[f(x_i)^{\prime}\approx\frac{f(x_{i+1})-f(x_{i})}{h}\] o bien (diferencias hacia atras) \[f(x_i)^{\prime}\approx\frac{f(x_{i})-f(x_{i-1})}{h}\] y grafique, junto con los errores

plot(x,fx, type="b",col="red")

 x<- seq(0.2,0.7, by=0.1); fx <- c(0.45, 0.57,0.71, 0.88,1.08,1.27)
# diferencia hacia adelante
c(x[1], (dfa = (fx[2]-fx[1])/(x[2]-x[1])))

## [1] 0.2 1.2

c(x[2], (dfa = (fx[3]-fx[2])/(x[3]-x[2])))

## [1] 0.3 1.4

x <-seq(0.2,0.7, by=0.1); fx <- c(0.45, 0.57,0.71, 0.88,1.08,1.27)
n <- length(x)-1
dfadelante <- seq(1:n)
dfa<- function (x, fx) {
for (i in 1:n) {
dfadelante[i] <- (fx[i+1]-fx[i])/(x[i+1]-x[1])
}}
dfa(x,fx);
dfadelante

## [1] 1 2 3 4 5

#en forma vectorial
x <- seq(0.2,0.7, by =0.1); fx <- c(0.45, 0.57,0.71, 0.88,1.08,1.27) 
fu <- fx[-1]; fu

## [1] 0.57 0.71 0.88 1.08 1.27

xu<- x[-1]; xu

## [1] 0.3 0.4 0.5 0.6 0.7

(dfa <- (fu- fx[1:5])/(xu-x[1:5]))

## [1] 1.2 1.4 1.7 2.0 1.9

(cbind(x[1:5],fx[1:5],dfa))

##               dfa
## [1,] 0.2 0.45 1.2
## [2,] 0.3 0.57 1.4
## [3,] 0.4 0.71 1.7
## [4,] 0.5 0.88 2.0
## [5,] 0.6 1.08 1.9

El principal problema de este metodo es que el parametro $h\equiv \Delta x$ debe ser pequeño, pero en terminos prácticos no puede ser tan pequenyo, pues la division entre cero afectara la precision de la maquina.

Una mejor aproximacion es \[f(x_i)^{\prime}\approx\frac{f(x_{i+1})-f(x_{i-1})}{2h}\sim \mathcal{O}(h^2)\] por tanto con un error menor, en un orden de magnitud. Tambien podemos escribir \[f(x_i)^{\prime}\approx\frac{f(x_{i}+h)-f(x_{i}-h)}{2h}\] que se conoce como la primera derivada en diferencias centradas.

#diferencia hacia atras
#en forma vectorial


x<- seq (0.2,0.7, by=0.1); fx <- c(0.45, 0.57,0.71, 0.88,1.08,1.27)
c(x [2], (fx[2]-fx[1])/(x[2]-x[1]))

## [1] 0.3 1.2

(cbind (x, fx))

##        x   fx
## [1,] 0.2 0.45
## [2,] 0.3 0.57
## [3,] 0.4 0.71
## [4,] 0.5 0.88
## [5,] 0.6 1.08
## [6,] 0.7 1.27

fx[-1]; fx[1:5]

## [1] 0.57 0.71 0.88 1.08 1.27

## [1] 0.45 0.57 0.71 0.88 1.08

x [-1]; x [1:5]

## [1] 0.3 0.4 0.5 0.6 0.7

## [1] 0.2 0.3 0.4 0.5 0.6

dfb <- (fx[-1]-fx[1:5])/(x[-1]-x[1:5]); dfb

## [1] 1.2 1.4 1.7 2.0 1.9

cbind(x[2:6],fx[2:6],dfb)

##               dfb
## [1,] 0.3 0.57 1.2
## [2,] 0.4 0.71 1.4
## [3,] 0.5 0.88 1.7
## [4,] 0.6 1.08 2.0
## [5,] 0.7 1.27 1.9

u<- cbind(x[1:5], fx[1:5], dfa); v <- cbind(x[2:6],fx[2:6], dfb)
(cbind(u,v))

##               dfa          dfb
## [1,] 0.2 0.45 1.2 0.3 0.57 1.2
## [2,] 0.3 0.57 1.4 0.4 0.71 1.4
## [3,] 0.4 0.71 1.7 0.5 0.88 1.7
## [4,] 0.5 0.88 2.0 0.6 1.08 2.0
## [5,] 0.6 1.08 1.9 0.7 1.27 1.9

Extrapolación de Richardson

Dado que el cálculo de la derivada numérica depende en general del parámetro $h^k;~k\in \mathbb{N}$, entre mayor sea el valor del exponente $k$ menos será la necesidad de contar con un parámetro $h\ll1$, y con ello se evitan los riesgos inherentes a dividir entre cero. Dado que, en el caso de la derivada numérica, no se conoce el valor real, $D_r$, sino una aproximación, $D_a$, el error bruto en el cálculo se puede escribir \[\varepsilon_r=D_r-D_a\] En el caso de la derivada centrada, el error es del orden de $h^2$, por tanto \[D_r(h)=(D_a)_{h}+Ch^2\] donde $(D_a)_{h}$ significa aproximar la derivada con un paso $h$, y $C$ es una constante desconocida. Si ahora se reduce el paso a la mitad, \[D_r(h/2)=(D_a)_{h/2}+C(h^2/4)\]

multiplicando por 4, $4D_r(h/2)=4(D_a)_{h/2}+C\cdot h^2$ y restando para eliminar $C$, tenemos \[D_r=\frac{4}{3}(D_a)_{h/2}-\frac{1}{3}(D_a)_{h}\] En otras palabras \[f^{\prime}(x_i)=\frac{4}{3}\cdot f^{\prime}_{\text{cent}}(x_i)_{h/2}-\frac{1}{3}\cdot f^{\prime}_{\text{cent}}(x_i)_{h}\] Esta ecuación se conoce como la ecuación de interpolación de Richardson. Aquí $f^{\prime}_{\text{cent}}(x_i)_{\Delta}$ es la derivada de diferencias centradas, valuada en $x_i$, con parámetro $\Delta$.

#ejercicio TI hecho en clase obtener las derivadas con diferencias centradas, para los punto que sea posible). para $x_OS$ tomar diferencia hacia adelante; para $x_n$ calcular diferencia hacia atrás

x<- seq (0.2,0.7, by=0.1); fx <- c(0.45, 0.57,0.71, 0.88,1.08,1.27)
(cbind (x, fx))

##        x   fx
## [1,] 0.2 0.45
## [2,] 0.3 0.57
## [3,] 0.4 0.71
## [4,] 0.5 0.88
## [5,] 0.6 1.08
## [6,] 0.7 1.27

n <- length(x)-1; dta <- seq(1:n+1)
dfa[1] <- (fx[2]-fx[1])/(x[2]-x[1])
for (i in 2:n) {
            dfa[i] <- 0.5*(fx[i+1]-fx[i-1])/(x[i+1]-x[i-1])
 }

dfa [n+1] <- (fx[n+1]-fx [n])/(x[n+1]-x [n])
(cbind(x,fx, dfa))

##        x   fx   dfa
## [1,] 0.2 0.45 1.200
## [2,] 0.3 0.57 0.650
## [3,] 0.4 0.71 0.775
## [4,] 0.5 0.88 0.925
## [5,] 0.6 1.08 0.975
## [6,] 0.7 1.27 1.900

#getwd()
#write.csv(cbind(x,fx,dfa), file #="C:/Users/jony2/OneDrive/Documentos/maestria/MATERIAS/ALGORITMOS/derivada numerica.csv")

ejercicio (Diferencias)

Completar el cuadro a continuación i) usando diferenecias adelante/atrás; ii) usando diferencias centradas

tabla 1

$i~x_i~f(x_i)~f^{\prime}(x_i)$ 0 0.4 1.45 1 0.5 2.57 2 0.6 3.71 3 0.7 4.88 4 0.8 6.12

tabla 2

0 −0.2 −1.23 1 0.0 −2.34 2 0.2 −0.34 3 0.4 2.46 4 0.6 4.35

Tómese un intervalo $[a,b]$, función $f(x)$, y entero $n$ dados. Haga na partición uniforme del intervalo en nodos $x_i$, con $h=\Delta x=\frac{b-a}{n}$
(es decir $x_i=a+i\Delta x$)
Calcule las derivadas $f^{\prime}_{\text{cent}}(x_i)_{h}$ y $f^{\prime}_{\text{cent}}(x_i)_{h/2}$ en los nodos $x_i$
Calcule la derivada $f^{\prime}(x_i)$ usando la fórmula de estrapolación de Richardson

\[ f(x)=e^x;~ n=10;~ [0,1]\\ f(x)=\sin x-x;~ n=12;~ [0,\pi] \\ f(x)=x^4-x^3 -1;~ n=20;~ [2,6] \\ f(x)=\ln x+2x;~ n=30;~ [3,7] \\ f(x)=\tan x-x^2;~ n=20;~ [-1,0]\]

En los mismos casos de 3. Calcular la derivada rwal, tabular/graficar ambas para comparar.

#!/usr/bin/env python3 # -- coding: utf-8 -- # ——————————————————————— # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program. If not, see http://www.gnu.org/licenses/ # ——————————————————————— # Implementación extrapolación de Richardson y algunos casos de salida.

from math import *

def pol(x): return x**3 + 4*x**2 - 10

def trig(x): “““Función de prueba”“” return x*cos(x-1) - sin(x)

def dercentrada(f, x, h): “““Retorna diferencias centradas”“” return (f(x + h) - f(x - h))/(2*h)

def richardson(f, x, h): ““” Implementación extrapolación de Richardson Entradas: f – función x – punto h – paso
```
  Salida:
  d -- aproximación a la derivada
 """

  d = (4/3)*dercentrada(f, x, h/2)-(1/3)*dercentrada(f, x, h)
  return d


  #pol(x), x = 1.5, h = 0.1
  print("Derivada función pol(x):")
  print("{0:.12f}".format(richardson(pol, 1.5, 0.1)))

  # trig(x), x = 4, h = 0.2
  print("Derivada función trig(x):")
  print("{0:.12f}".format(richardson(trig, 4, 0.2)))
```

Integración numérica

En cálculo integral, la integral de Riemann se construye sobre na partición del intervalo $[a,b]=[x_0,x_n]$, que se ha dividido en subintervalos $[x_{i-1},x_i]$ entonces \[R=\sum_{i=0}^nf(x^\star_i)\cdot \Delta x_i;\quad \Delta x_i=x_i-x_{i-1}\]

para $x^\star_i\in [x_{i-1},x_i]$. Esta aproximación llegará al verdadero valor de la itegral cuando $\Delta x\rightarrow 0$.

ejemplo (Integral de Gauss)

Integraremos $e^{-x^2}$ en un intervalo arbitrario $[a,b]$

##Regla de los trapecios

Normalmente se aproxima el área bajo la curva de $f_x$ mediante el cálculo del área de un trapecio.

\[\int_a^bf(x)dx\approx \frac{f(a)+f(b)}{2}\cdot (b-a)\] Normalmente, el intervalo $[a,b]$ estará dividido en subintervalos $[x_{i-1},x_i)$ de modo que \[\int_a^bf(x)dx\approx \sum_{i=1}^n\frac{f(x_{i-1})+f(x_i)}{2}\cdot \Delta x_i\] conocida como la regla de los trapecios. SI la partición del intervalo tiene longitudes uniforme \[\int_a^bf(x)dx\approx \frac{h}{2}\left( f(x_0)+2f(x_1)+2f(x_2)+\ldots+2f(x_{n-1})+f(x_n)\right)\] la cual es más usada en la práctica.

El término de error \[\exists \xi\in (a,b):\quad \int_a^bf(x)dx= \frac{h}{2}\left( f(x_0)+2\sum_{i=1}^{n-1}f(x_i)+ f(x_n) \right)-\frac{b-a}{12}h^2f^{\prime \prime}(\xi)\]

ejemplo (estimación de error)

Estimar el valor de $n$ para que usando la regla del trapecio se pueda tener una aproximación del valor exacto de $\int_{-1}^{1}e^{-x^2}dx$ con una precisión menor a $10^{-6}$

Regla de Simpson

En la regla del trapecio, para aproximar el área bajo la curva, se usa interpolación lineal a trozos. Un mejora natural sería usar interpolación de grado superior, y de hecho, si se usan polinomios cuadráticos, se llega a un método conocido como la regla de Simpson. Esta aproximación utiliza un poinomio de Lagrange $P_3(x$)

\[\int^{x_2}_{x_0} f(x) \space dx = \frac{h}{3}[f(x_0) + 4f(x_1) + f(x_2)] – \frac{h^2}{90}f^{(4)}(\epsilon)\] donde $x_1 = x_0 + h$, $h = (x_2 – x_0)/2$ y $\frac{h^2}{90}f^{(4)}(\epsilon)$ es el t´rmino de error.

ejemplo

Analizar \[a)~ \int _0^{\pi/4}x\sin xdx \\ b)~\int _0^1 x^2 e^{-x}dx\]

En R se puede realizar la integración numérica direcamente con la función integrate.

Ejemplo 1 (integración numérca)

Hallar el valor de la integral \[I=\int _0^{\infty} \frac{1}{\sqrt{x}(1+x)}dx\]

## define the integrated function
integrand <- function(x) {1/((x+1)*sqrt(x))}
## integrate the function from 0 to infinity
integrate(integrand, lower = 0, upper = Inf)

## 3.141593 with absolute error < 2.7e-05

Podemos decir que el resultado numérico es muy cercano a $\pi$, \[I\approx \pi \]

Ejemplo 2

Un ejemplo similar

Hallar el valor de \[I_k = \frac{k!}{\sqrt{2\pi}}\int _{-0.96}^{0.96}e^{-x^2/4}\cdot e^{-kx}dx; \qquad k\in[1,3]\] donde queremos integrar más de una vez, pues $k$ es un parámetro que toma valores en un intervalo dado. Nos sirve tener una función que puede integrar numéricamente.

## define the integrated function
Pi <-3.141579; a=1; b=3; n=10; 
  dk <- (b-a)/n 
#Intervalo donde está el parámtero:
  KI <- seq(n); kvar <- seq(n)
for (i in 1:(n+1)){
  k <- a + (i-1)*dk 
  integrand <- function(x) {(factorial(k)/(sqrt(2*Pi)))*exp(-x^2/4)*exp(-k*x)}
  ## integrate the function 
  KI[i] <- integrate(integrand, lower = -0.96, upper = 0.96)$value
  kvar[i] <- k
}
#  Tenemos los valores de las integrales en el intervalo [1,3]
 plot(kvar,KI, type = "b", col="blue", xlab = "k", ylab = "I")

Integración numérica (derivados asiáticos)

Observemos que en un modelo continuo, el promedio del activo de una opción asiática se debe calcular con una integral \[\bar{S}=\frac{1}{T}\int _0 ^TdS_t\] donde $S$ es de hecho la variable aleatoria, con riesgo, cuya distribución piede ser normal, u otra. Pero ¿cómo se pueden realizar las integrales numéricamente?

Pensemos en alguna función real, $f(x)=e^{-x^2} +1$, la cual se quiere integrar en un intervalo $x\in[a, b]$. La función $f$ no es aleatoria, es decir, para cada valor de $x=a$ dado, hay un valor $y=e^{-a^2} +1$ con $100$% de probabilidad. Pero la función integral $F(y)=\int ^y_{-\infty } df(x)$ no se conoce.

Observamos el integrando

f <- function(x){
  x*exp(-x^2) +1
}
curve(f, from = 0, to= 5, ylab="función", col="red")

De hecho, ya sabemos cómo hacerlo:

#definir el integrando
f <- function(x){
  x*exp(-x^2) +1
}
#Definir nuestra propia función Riemann-integradora
INtf <- function(f,a,b,N){
  dx = (b-a)/N
  # damos valores dentro del intervalo a integrar:
  x <- seq(a,b, by = dx)
  #los valores desde a + 1/N
  y <- x[-1]
  x <- x[-length(x)]
  #sumas inferiores y superiores
  Lft <- dx*f(x); Lsum <- sum(Lft)
  Rgt <- dx*f(y); Rsum <- sum(Rgt)
  #c("min/max sums",Lsum, Rsum)
  (Lsum+Rsum)/2
  #f(x)-f(y)
}

INtf(f,0,3,N=100)

## [1] 3.499863

# verificamos, usando el default de R
("verificación")

## [1] "verificación"

integrate(f,lower= 0,upper= 3)

## 3.499938 with absolute error < 1.6e-11

CHECAR? OK

Ejercicio TAREA 6

Hallar un valor para la integral de la función logística \[y=\frac{b}{1+e^{-ax}}\] donde $a=0.5$ y $b=10$ para $x\in[0,10]$ Observamos el integrando.

a <- 0.5; b<-10
logst <- function(x,a,b){
  b/(1+exp(-a*x))
}
curve(logst(x,0.5,10), from = 0, to = 15, col="red")

Integral aleatoria (I)

Supongamos ahora que la variable $X$ no es más que una función aleatoria \[X: \Omega \rightarrow \mathbb{R}\] Intentemos hacer el mismo cálculo de arriba, es decir integrando $f(X)=e^{-X^2} +1$. Supongamos aún que la variable es de Laplace, $f_X \sim e^{-X/s}$.

Usaremos el paquete LaplacesDemon

#install.packages("LaplacesDemon")
library(LaplacesDemon)

## Warning: package 'LaplacesDemon' was built under R version 4.2.2

#x <- dlaplace(1,0,1)
#x <- plaplace(1,0,1)
#x <- qlaplace(0.5,0,1)
#x <- rlaplace(100,0,1)

#Plot Probability Functions
x <- seq(from=-5, to=5, by=0.1)
plot(x, dlaplace(x,0,0.5), ylim=c(0,1), type="l", main="Probability Function",
     ylab="density", col="red")
lines(x, dlaplace(x,0,1), type="l", col="green")
lines(x, dlaplace(x,0,2), type="l", col="blue")
legend(2, 0.9, expression(paste(mu==0, ", ", lambda==0.5),
     paste(mu==0, ", ", lambda==1), paste(mu==0, ", ", lambda==2)),
     lty=c(1,1,1), col=c("red","green","blue"))

fr <- function(x){
  x*exp(-x^2) +1
}
s = 0.2
INtXLfr <- function(fr,a,b,N){
  dx = (b-a)/N
  x <- rlaplace(N,0,s)
  xl <- sort(x)
  y <- xl[-1]
  x <- xl[-length(x)]
  Lft <- dx*fr(x)*exp(-x*s); Lsum <- sum(Lft)
  Rgt <- dx*fr(y)*exp(-x*s); Rsum <- sum(Rgt)
  #c("min/max sums",Lsum, Rsum) 
    c("Sampling INtegral", (Lsum+Rsum)/2)
}
INtXLfr(fr,a=0,b=3,N=100)

## [1] "Sampling INtegral" "2.90416257432609"

("verificación")

## [1] "verificación"

integrate(fr,0,3)

## 3.499938 with absolute error < 1.6e-11

Pero ¿qué función hemos integrado numéricamente?

fr <- function(x){
  (x*exp(-x^2) +1)
}
s = 0.2
a=0; b=3; N=100
x <- rlaplace(N,0,s)
xl <- sort(x)
  
#    fr(xl)
#    plot(fr(xl), type="p", main="Random Function",
#     ylab="function", col="blue")
    plot(xl,fr(xl), type="b", main="Evaluated Function",
     ylab="function", col="red")

Como observamos, al tomar simplemente $N$ valores de esta distribución, la variable aleatoria $X$ cae incluso fuera del rango de integración.

Cada integral que calculemos será una variable aleatoria
Cada integral que calculemos tendrá un rango aleatorio.

NOTA Ahora restringimos la variable aleatoria $X$ al intervalo de integración localizado:

fr <- function(x){
  x*exp(-x^2) +1
}
s = 0.2
INtXLfr <- function(fr,a,b,N){
  dx = (b-a)/N
  xl <- seq(N); i <- 1;
  while(i <= N){
    x <- rlaplace(1,0,s)
    # selecting in the integration domain
    if((a <= x) & (x <=b)){ xl[i] <- x; i = i +1 }
  }
    xlr <- sort(xl)
    
  y <- xlr[-1]
  x <- xlr[-length(x)]
  Lft <- dx*fr(x); Lsum <- sum(Lft)
  Rgt <- dx*fr(y); Rsum <- sum(Rgt)
  #c("min/max sums",Lsum, Rsum) 
    c("INtegral", (Lsum+Rsum)/2)
}
INtXLfr(fr,a=0,b=3,N=10000)

## [1] "INtegral"         "3.49287062106691"

("verificación")

## [1] "verificación"

integrate(fr,0,3)

## 3.499938 with absolute error < 1.6e-11

Esta vez obtendremos un valor similar al de arriba, dado que a final de cuentas sólo hemos escogido $N$ valores aleatorios en el dominio de integración. Esta integral es una integral de Riemann

Observemos ahora la función que hemos integrado numéricamente:

fr <- function(x){
  x*exp(-x^2) +1
}
s = 0.2
a=0; b=3; N=1000
xl <- seq(N); i <- 1;
# Condition: in the Integration Domain
  while(i <= N){
    x <- rlaplace(1,0,s)
    if((a <= x) & (x <=b)){ xl[i] <- x; i = i +1 }
  }
    xlr <- sort(xl)
  
#    fr(xl)
#    plot(fr(xl), type="p", main="Random Function",
#     ylab="function", col="blue")
    plot(xlr,fr(xlr), type="p", main="Evaluated Function",
     ylab="function", col="red")

Conclusión

El método de integración por muestreo que hemos visto se puede emplear como sigue:

Si queremos hallar la integral \[I=\int _a^bF(x)dx\] donde $F(x)$ es una función desconocida, podemos muestrear esta función con una secuencia de números aleatorios $X_1, X_2, \ldots , X_N\in [a,b]$, con $N >>1$, y calcular una suma de Riemann \[\sum ^N F(X_i)\Delta X_i \underset{N\rightarrow \infty }{\longrightarrow}I\]

Entonces, entre más puntos muestreemos $X\in [a,b]$ mejor será la aproximación. ¿Cómo sabremos cuátos puntos son suficientes? Esto puede ser complicado para funciones con muchas variaciones en $[a,b]$

Hay varios métodos de integración estocástica que resultan mejores.

Integración: con Monte Carlo

El método de Monte Carlo abarca cualquier técnica de muestreo estadístico empleada para aproximar soluciones a problemas cuantitativos. Esencialmente, el método Monte Carlo resuelve un problema simulando directamente el proceso subyacente y luego calculando el resultado promedio del proceso. Este enfoque general es válido en áreas como física, química, informática, etc.

En finanzas, el método de Monte Carlo se utiliza para simular las diversas fuentes de incertidumbre que afectan el valor del instrumento, cartera o portafolio de inversión en cuestión. Posteriormente, se calcula un valor representativo (e.g. el valor esperado) dado en contexto de posibles valores de los datos subyacentes. En otras palabras, cubriendo todas las contingencias imaginables del mundo real en proporción a su probabilidad. En términos de teoría financiera, esto es, esencialmente, una aplicación de la valoración neutral al riesgo.

En finanzas los casos donde puede llegar a utilizarse son

  Finanzas corporativas, proyecciones financieras;
  Análiss de opciones reales
  Valuación de opciones (asiáticas)
  Valuación de Bonos; CMO's bajo tasa de interés
  Valor en Riesgo 
  Simulación de portafolios e Inversiones

Aunque los métodos de Monte Carlo brindan flexibilidad, y pueden manejar múltiples fuentes de incertidumbre, el uso de estas técnicas no siempre es apropiado. En general, los métodos de simulación MC se preferirán a otras técnicas de valuación solo cuando hay varias variables de estado (es decir, varias fuentes de incertidumbre).

Veamos casos simples. En muchos casos, lo común el cálculo de un valor esperado \[\mathbb{E}[f(X)]=\int _a^bf(X)P(X)dx\] deonde $P(X)$ es la distribución de $X$.

Paquete SI (stochastic Integration)

El paquete SI (Stochastic Integration) proporciona varios métodos de integración Monte Carlo (MC) que incluyen

Método de punto estocástico
Método de valor medio
Método de muestreo relevante (importante)
Método de muestreo estratificado

Tengamos en cuenta que Stochastic Point Methodes es simplemente un método de punto estocástico para estimar la integración. Sin embargo, se proporciona éste en el paquete con el propósito de ser fácilmente comparado con los otros métodos de MC. Estos cuatro métodos pertenecen a métodos estocásticos MC.

Integración estocástica

Supongamos que $f(x)\in \mathbb{C}^1$ es una función de derivada continua, la cual está definida en el intervalo $D=[a,b]$ y es acotada (s.p.g.) \[0\leq f(x)\leq M\] donde $M$ es una cota superior.

Necesitamos calcular la integral \[I=\int_a^bf(x)\mathrm{d}x\] En otras palabras, queremos obtener el área bajo $f(x)$, \[A=\{(x,y):0\leq y\leq f(x),\ x\in D\}\]

Método de punto estocástico

Hacemos un muestreo aleatorio de $N$ puntos en el rectángulo $R=D\times(0,M)$, obteniendo variables aleatorias $Z_1,\cdots,Z_N$, donde $Z_i=(X_i,Y_i),\ i=1,2,\ldots,N$.

Para $i=1,2,\ldots, N$, tomamos la función indicadora \[\xi_i=\begin{cases} 1&,Z_i\in A\\ 0&,\text{otra forma} \end{cases}\]

Entonces, $\{ \xi _i \}$ son ensayos de Bernoulli $\xi_1,\cdots,\xi_N\overset{i.i.d.}{\sim}Bernoulli(1,p)$ donde el parámetro de probabilidad \[p=\mathbb{P}(Z_i\in D)=\dfrac{I}{M(b-a)}\] Al muestrear $N$ veces, obtenemos el estimador $\hat{p}=(\sum ^N\xi _i)/N \underset{N\rightarrow \infty }{\longrightarrow} p$ (de acuerdo con la Ley de los Grandes Números), que se calcula como la proporción de puntos $\{Z_i\}$ que quedan debajo de la curva de $f(x)$. De estamanera, obtendremos también un estimador de la integral \[\hat{I}=\hat{p}\cdot M(b-a) \underset{N\rightarrow \infty }{\longrightarrow} I\] Usando el Teorma del Límite Central podemos estimar la varianza de la aproximación $\hat{I}$ \[ \dfrac{\hat{p}-p}{Var(p)}=\dfrac{\hat{p}-p}{\sqrt{\dfrac{p(1-p)}{N}}}\xrightarrow{D}N(0,1)\\ \dfrac{\hat{I}_1-I}{\sqrt{\frac{1}{N}}}\xrightarrow{D}N(0,[M(b-a)]^2p(1-p)) \] de modo que \[Var(\hat{I})\approx \dfrac{1}{N}[M(b-a)]^2\hat{p}(1-\hat{p})\] ### Ejemplo (Stochastic Point Method) en R En R esto se puede implementar de la siguiente manera

#install.packages("SI")
library(SI)

## Warning: package 'SI' was built under R version 4.2.2

Entonces usaremos

## To integrate (exp(x)+1)^exp(x) from -1 to 1
set.seed(0)
h <- function(x){
    (exp(x)+1)^exp(x)
}
N <- 1000000
# Upper bound
M <- h(1)
SPMresult <- SI.SPM (h,-1,1,M,N)
I1 <- SPMresult[[1]]
VarI1 <- SPMresult[[2]]
(c("integral approximation", I1, "st. dev.", sqrt(VarI1)))

## [1] "integral approximation" "9.875270193367"         "st. dev."              
## [4] "0.0245727787011359"

Método del valor medio

Supongamos una variable aleatoria uniforme $U\sim \text{Unif}(a,b)$ Queremos calcular \[\mathbb{E}[f(U)]=\int_a^bf(x)\dfrac{1}{b-a}\mathrm{d}x;\\ \mathbb{E}[f(U)]=\dfrac{I}{b-a}\\ I=(b-a)\mathbb{E}[f(U)].\] Supóngase que se muestrean $N$ valores uniformes $U_1,\cdots,U_N\overset{i.i.d.}{\sim}Unif(a,b)$; entonces las variables aleatorias $Y_i=f(U_i),\ i=1,2,\cdots,N$, cuya media cumple, por Ley de Grandes Números, \[\overline{Y}=\dfrac{1}{N}\sum\limits_{i=1}^Nf(U_i)\xrightarrow{a.s.}\mathbb{E}[f(U)]=\dfrac{I}{b-a}\] De modo que se tiene un estimador de la integral, calculando el valor promedio \[\hat{I}_2=\dfrac{b-a}{N}\sum\limits_{i=1}^Nf(U_i)\] Para estimar la varianza, usando el Teorma del Límite Central \[\dfrac{\hat{I}-I}{\sqrt{1/N}}\xrightarrow{Dist.}N(0,(b-a)^2Var[f(U)])\\ \hat{I}\xrightarrow{Dist.}N(I,\dfrac{(b-a)^2}{N}Var[f(U)])\] donde \[Var[f(U)]=\int_a^b(f(u)-\mathbb{E}[f(U)])^2\dfrac{1}{b-a}\mathrm{d}u\approx \dfrac{1}{N}\sum\limits_{i=1}^N(Y_i-\overline{Y})^2\] dado quese puede estimar con la varianza muestral. Entonces tenemos, finalmente \[Var(\hat{I}_2)\approx \dfrac{(b-a)^2}{N^2}\sum\limits_{i=1}^N(Y_i-\overline{Y})^2\] ### Ejemplo (MVM) en R En R esto se puede implementar de la siguiente manera

## To integrate exp(x) from -1 to 1
set.seed(0)
h <- function(x){
    (exp(x)+1)^exp(x)
}
N <- 100000
MVMresult <- SI.MVM(h,-1,1,N)
I2 <- MVMresult[[1]]
VarI2 <- MVMresult[[2]]
(c("integral approximation", I2, "st. dev.", sqrt(VarI2)))

## [1] "integral approximation" "9.89444080651512"       "st. dev."              
## [4] "0.0420044576408042"

Verifiquemos los ejemplos anteriores usando la integral de Riemann. Para obtener \[I=\int _{-1}^1dx(1+e^x)^{e^x}\]

## define the integrated function
integrand <- function(x) {(1+exp(x))^exp(x)}
## integrate the function from 0 to infinity
integrate(integrand, lower = -1, upper = 1)

## 9.873631 with absolute error < 5.1e-07

Método del mapeo relevante

Supongamos que $g(x)$ es una densidad similar a $|f(x)|$; en el siguiente sentido, $g(x)=0\Rightarrow f(x)=0$, y además $f(x)=o(g(x)) $ cuando $x\rightarrow \infty$.

Supongamos además que tomammos $N$ variables aleatorias con la misma densidad $X_i\overset{i.i.d.}{\sim}g(x),\ i=1,2,\cdots,N$ y sea \[\eta_i=\dfrac{f(X_i)}{g(X_i)}\] entonces, \[\mathbb{E}[\eta_i]=\int_D\dfrac{f(x)}{g(x)}\cdot g(x)\mathrm{d}x=I\] Tenemos entonces la aproximación de la integral mediante el estimador de muestreo relevante \[\hat{I}_3=\overline{\eta}=\dfrac{1}{N}\sum\limits_{i=1}^N\dfrac{f(X_i)}{g(X_i)}\] cuya varianza puede ser estimada mediante \[Var(\hat{I}_3)=Var(\overline{\eta})=\dfrac{1}{N}Var\left(\dfrac{f(X)}{g(X)}\right)\\ Var(\hat{I}_3)\approx \dfrac{1}{N^2}\sum\limits_{i=1}^N \left[\dfrac{f(X_i)}{g(X_i)}-\hat{I}_3\right]^2\]

Ejemplo (Important Sampling Method) en `R`

En R este método se puede implementar de la siguiente manera

## To integrate exp(x) from -1 to 0.5
## Use the sampling density exp(x) ~ (3/2+x)/3
set.seed(0)
h <- function(x){
    (exp(x)+1)^2
}
N <- 1000000    
g <- function(x){return( ( 1+(3+2*x)/6)^4 )}
# Valores max/min de g(x)
mg <- g(-1); Mg <- g(0.5)
G_inv <- function(y){return(3*y^(1/4)-9/2)}
ISMresult <- SI.ISM(h,g,G_inv,N, min_G = mg, max_G = Mg )
I3 <- ISMresult[[1]]; I3

## [1] 4.896199

VarI3 <- ISMresult[[2]]; sqrt(VarI3)

## [1] 0.0002973158

Se observan las funciones $f(x)=(1+e^x)^2$ (azul) y $g(x)=(1+(3+2x)/6)^4$ (rojo)del ejemplo ISM (Important Sampling Method)

Verificación

## define the integrated function
integrand <- function(x) {(1+exp(x))^2}
## integrate the function from 0 to infinity
integrate(integrand, lower = -1, upper = 0.5)

## 5.353157 with absolute error < 5.9e-14

Método de muestreo estratificado

Suponga que hacemos una m-partición del dominio $D=\bigcup\limits_{j=1}^m D_j$, donde $D_i\cap D_j=\varnothing$. Para cada elemento de la partición $D_j$ usamos el método del valro medio para aproximar la integral \[\hat{I}_{D_j}=\int_{D_j}f(x)\mathrm{d}x\] y obtenemos \[\hat{I}_4=\sum\limits_{j=1}^m\hat{I}_{D_j}\] Se puede mostrar quela varianza de esta estimación es menor que $Var(I_2)$

Ejemplo (Stratified Sampling Method)

## To integrate exp(x) from -1 to 1
set.seed(0)
h <- function(x){
   exp(x)
    #(1+exp(x))^2
}
N <- 10000
m <- 10
SSMresult <- SI.SSM(h,-1,0.5,m,N)

##       to                                           
## [1,] 0.5 0.65 0.8 0.95 1.1 1.25 1.4 1.55 1.7 1.85 2

I4 <- SSMresult[[1]]
VarI4 <- SSMresult[[2]]
(c("integral approximation", I4, "st. dev.", sqrt(VarI4)))

## [1] "integral approximation" "5.74069597245384"       "st. dev."              
## [4] "0.00271224065505527"

(p: HAY UN ERROR EN L PARTICIÓN?) NO obstante, usando la integral de Riemann

## define the integrated function
integrand <- function(x) {(1+exp(x))^2}
## integrate the function from 0 to infinity
integrate(integrand, lower = -1, upper = 0.5)

## 5.353157 with absolute error < 5.9e-14

(VER) Introducing Monte Carlo Methods with R, 2009, Springer-Verlag•Data and R programs for the course available at http://www.stat.ufl.edu/ casella/IntroMonte/

Ejercicio TAREA 7

Evaluar la integral \[I=\int _0^{\infty}\frac{3}{4}xe^{-x^{3/4}}dx\]

El resultado exacto es una función Gamma $I = \Gamma \left(\frac{20}{3}\right) \approx 389.035$

Ejemplo (VaR)

Calcularemos un Valor en Riesgo para un portafolio hipotético sobre dos activos con riesgo, \[\pi (x,y) =w_1S_1+w_2S_2\] donde $w_1=w_2=\frac{1}{2}$

# constructing hypothetical portfolio consisting of x1 and x2
n <- 10000
pi <- matrix(c(rnorm(n,50,4),rnorm(n,5,0.5)),ncol=2);
colnames(pi)<-c("x1","x2"); #pi 

#-x1 and x2 weights

weights<-c(0.5,0.5); c("weights",weights)

## [1] "weights" "0.5"     "0.5"

# calculated means and covariance matrix

mu <-apply(pi,2,mean); sigma<-cov(pi)
c("medias",mu);

##                                    x1                 x2 
##           "medias" "50.0247078182965"  "4.9990596105608"

("El modelo de riesgo")

## [1] "El modelo de riesgo"

sigma

##              x1           x2
## x1 15.992756270 -0.008218783
## x2 -0.008218783  0.249513496

# generate n scenarios for x1 and x2 with given covariance matrix sigma
# `Modern Applied Statistics with S`

library(MASS);

# MASS::mvnorm Simulate from a Multivariate Normal Distribution
# Simulación

MC <- mvrnorm(n,mu,sigma)

#calculate portfolio value for simulated x1 and x2

MCportfolio <- MC%*%as.matrix(weights);

str(MCportfolio)

##  num [1:10000, 1] 27.1 26.3 29.6 30.5 28 ...

# find VaR for 95%

("VaR portfolio")

## [1] "VaR portfolio"

quantile(MCportfolio,p=0.05)

##       5% 
## 24.18799

Rnotebook30012020

A Salazar

14/02/2023

Diferenciación numérica

ejemplo

Extrapolación de Richardson

ejercicio (Diferencias)

tabla 1

tabla 2

Integración numérica

ejemplo (Integral de Gauss)

ejemplo (estimación de error)

Regla de Simpson

ejemplo

Ejemplo 1 (integración numérca)

Ejemplo 2

Integración numérica (derivados asiáticos)

Ejercicio TAREA 6

Integral aleatoria (I)

Conclusión

Integración: con Monte Carlo

Paquete SI (stochastic Integration)

Integración estocástica

Método de punto estocástico

Método del valor medio

Método del mapeo relevante

Ejemplo (Important Sampling Method) en `R`

Método de muestreo estratificado

Ejemplo (Stratified Sampling Method)

Ejercicio TAREA 7

Ejemplo (VaR)

Rnotebook30012020

A Salazar

14/02/2023

Diferenciación numérica

ejemplo

Extrapolación de Richardson

ejercicio (Diferencias)

tabla 1

tabla 2

Integración numérica

ejemplo (Integral de Gauss)

ejemplo (estimación de error)

Regla de Simpson

ejemplo

Ejemplo 1 (integración numérca)

Ejemplo 2

Integración numérica (derivados asiáticos)

Ejercicio TAREA 6

Integral aleatoria (I)

Conclusión

Integración: con Monte Carlo

Paquete SI (stochastic Integration)

Integración estocástica

Método de punto estocástico

Método del valor medio

Método del mapeo relevante

Ejemplo (Important Sampling Method) en R

Método de muestreo estratificado

Ejemplo (Stratified Sampling Method)

Ejercicio TAREA 7

Ejemplo (VaR)

Ejemplo (Important Sampling Method) en `R`