Cálculo Estocástico con R: de la teoría a la simulación
Movimiento Browniano, Lema de Itô, Ornstein-Uhlenbeck, Puente Browniano y valoración de opciones
Economista Jeel Cueva
28/09/2026
1 Introducción
1.1 ¿Qué es el cálculo estocástico?
El cálculo tradicional (el de Newton y Leibniz) estudia funciones que cambian de manera suave y predecible: si conozco la posición de un auto y su velocidad, puedo saber dónde estará dentro de un segundo.
Pero muchas variables económicas no se comportan así. El precio de una acción, el tipo de cambio, la tasa de interés o el precio del cobre cambian a cada instante y de forma aleatoria. No basta con conocer su tendencia: también hay que modelar su incertidumbre.
El cálculo estocástico es la rama de las matemáticas que permite trabajar con funciones que evolucionan en el tiempo con un componente aleatorio. Su objeto de estudio son los procesos estocásticos: colecciones de variables aleatorias ordenadas en el tiempo, \(\{X_t\}_{t \geq 0}\).
💡 Intuición: una ecuación diferencial ordinaria dice “el cambio de X es igual a su tendencia”. Una ecuación diferencial estocástica dice “el cambio de X es igual a su tendencia más un shock aleatorio”.
\[ \underbrace{dX_t}_{\text{cambio}} = \underbrace{a(X_t, t)\,dt}_{\text{tendencia (drift)}} + \underbrace{b(X_t, t)\,dW_t}_{\text{shock aleatorio (difusión)}} \]
1.2 ¿Por qué le importa a un economista?
- Finanzas: la fórmula de Black-Scholes, que valora opciones y otros derivados, se construye sobre el cálculo estocástico.
- Macroeconomía y política monetaria: los modelos de tasas de interés (Vasicek, CIR) usan procesos con reversión a la media.
- Gestión de riesgos: el Valor en Riesgo (VaR) y las pruebas de estrés se basan en simular miles de escenarios posibles.
- Economía de recursos naturales: los precios de commodities como el cobre o el petróleo se modelan con procesos estocásticos para evaluar proyectos mineros o energéticos.
- Opciones reales: permiten valorar la flexibilidad de postergar, ampliar o abandonar una inversión bajo incertidumbre.
1.3 ¿Qué aprenderás en este notebook?
| Sección | Tema | Pregunta que responde |
|---|---|---|
| 3 | Movimiento Browniano estándar | ¿Cómo se modela el “ruido” puro? |
| 4 | Movimiento Browniano Geométrico | ¿Cómo evoluciona el precio de un activo? |
| 5 | Lema de Itô | ¿Cómo se deriva una función de un proceso aleatorio? |
| 6 | Proceso de Ornstein-Uhlenbeck | ¿Cómo se modela una variable que regresa a su promedio? |
| 7 | Puente Browniano | ¿Qué pasa si conozco el punto de llegada? |
| 8 | Valoración de opciones | ¿Cuánto vale hoy un derecho a comprar en el futuro? |
| 9 | Convergencia de Monte Carlo | ¿Cuántas simulaciones necesito? |
Cada sección sigue la misma estructura: teoría → código → interpretación.
2 Preparación del entorno
2.1 Librerías
Usaremos R con las siguientes librerías:
ggplot2: gráficos de alta calidad.dplyrytidyr: manipulación de datos.gridExtraypatchwork: combinación de gráficos.scales: formato de ejes (por ejemplo, separador de miles).
El siguiente bloque revisa si cada paquete está instalado; si no lo está, lo instala y luego lo carga.
2.2 Tema gráfico y paleta de colores
Para que todos los gráficos tengan un estilo uniforme y profesional,
definimos una función de tema propia (theme_jeel) y una
paleta de colores. Así evitamos repetir el mismo formato en cada
gráfico.
theme_jeel <- function() {
theme_minimal(base_size = 13, base_family = "sans") +
theme(
plot.title = element_text(face = "bold", size = 15, color = "#1a1a1a"),
plot.subtitle = element_text(size = 11, color = "#555555"),
plot.caption = element_text(size = 9, color = "#888888", face = "italic"),
axis.title = element_text(face = "bold", color = "#333333"),
axis.text = element_text(color = "#444444"),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey90", linewidth = 0.3),
legend.position = "bottom",
legend.title = element_text(face = "bold"),
plot.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA)
)
}
paleta <- c("#0B3C5D", "#328CC1", "#D9B310", "#B82601",
"#1D2731", "#4A7C59", "#7B2D8B", "#E36414")2.3 Parámetros globales: cómo “discretizamos” el tiempo
Una computadora no puede trabajar con tiempo continuo, así que dividimos el intervalo \([0, T]\) en \(N\) pasos pequeños de tamaño:
\[ \Delta t = \frac{T}{N} \]
Con \(T = 1\) año y \(N = 1000\) pasos, cada paso equivale a \(\Delta t = 0.001\) años (aproximadamente un tercio de día). Mientras más pequeño sea \(\Delta t\), más se parece la simulación al proceso continuo.
La instrucción set.seed(2024) fija la “semilla” del
generador de números aleatorios. Esto garantiza que cualquier
persona que ejecute este código obtenga exactamente los mismos
resultados, lo que hace el análisis reproducible.
set.seed(2024)
T_final <- 1 # Horizonte (1 año)
N <- 1000 # Pasos temporales
dt <- T_final / N
t_grid <- seq(0, T_final, length.out = N + 1)
n_sim <- 50 # Número de trayectorias3 Movimiento Browniano estándar
3.1 Teoría
El Movimiento Browniano (o proceso de Wiener), denotado \(W_t\), es el bloque fundamental del cálculo estocástico. Su nombre proviene del botánico Robert Brown, quien en 1827 observó el movimiento errático de partículas de polen en el agua. Louis Bachelier lo usó en 1900 para modelar precios en la Bolsa de París, y Norbert Wiener le dio su formalización matemática.
Un proceso \(\{W_t\}_{t \ge 0}\) es un Movimiento Browniano estándar si cumple cuatro propiedades:
- Comienza en cero: \(W_0 = 0\).
- Incrementos independientes: lo que ocurre entre \(t\) y \(s\) no depende de lo ocurrido antes. El pasado no ayuda a predecir el shock futuro.
- Incrementos normales: \(W_t - W_s \sim \mathcal{N}(0,\; t - s)\) para \(s < t\). El cambio tiene media cero y una varianza igual al tiempo transcurrido.
- Trayectorias continuas: el proceso no da saltos.
De la propiedad 3 se deduce que \(W_t \sim \mathcal{N}(0, t)\): la incertidumbre crece con el tiempo, pero su desviación estándar crece como \(\sqrt{t}\), no como \(t\).
💡 Intuición: imagina a una persona que camina dando pasos aleatorios hacia adelante o hacia atrás. En promedio no se mueve (\(E[W_t] = 0\)), pero mientras más tiempo pase, más lejos puede estar de su punto de partida.
3.1.1 ¿Cómo se simula?
Usando la propiedad 3, cada incremento se genera como:
\[ \Delta W = W_{t+\Delta t} - W_t = \sqrt{\Delta t}\; Z, \qquad Z \sim \mathcal{N}(0,1) \]
y la trayectoria completa se obtiene sumando acumuladamente estos incrementos:
\[ W_{t_k} = \sum_{i=1}^{k} \Delta W_i \]
3.1.2 La banda de confianza del 95 %
Como \(W_t \sim \mathcal{N}(0, t)\), el 95 % de las trayectorias debería estar, en cada instante \(t\), dentro de:
\[ \pm 1.96\,\sqrt{t} \]
Esta banda tiene forma de “embudo” que se abre con el tiempo.
3.2 Código
simular_bm <- function(n_sim, N, dt) {
# Matriz de incrementos: N filas (tiempo) x n_sim columnas (trayectorias)
dW <- matrix(rnorm(n_sim * N, mean = 0, sd = sqrt(dt)),
nrow = N, ncol = n_sim)
# Suma acumulada de cada columna, comenzando en W_0 = 0
W <- apply(dW, 2, function(x) c(0, cumsum(x)))
return(W)
}
W <- simular_bm(n_sim, N, dt)
# Pasamos la matriz a formato "largo" para ggplot
df_W <- data.frame(
t = rep(t_grid, times = n_sim),
W = as.vector(W),
id = factor(rep(1:n_sim, each = N + 1))
)
# Banda teórica ±1.96*sqrt(t)
banda <- data.frame(
t = t_grid,
upper = 1.96 * sqrt(t_grid),
lower = -1.96 * sqrt(t_grid)
)p1 <- ggplot() +
geom_line(data = df_W, aes(x = t, y = W, group = id, color = id),
alpha = 0.35, linewidth = 0.5) +
geom_ribbon(data = banda, aes(x = t, ymin = lower, ymax = upper),
fill = "#D9B310", alpha = 0.15) +
geom_line(data = banda, aes(x = t, y = upper),
color = "#B82601", linetype = "dashed", linewidth = 0.8) +
geom_line(data = banda, aes(x = t, y = lower),
color = "#B82601", linetype = "dashed", linewidth = 0.8) +
geom_hline(yintercept = 0, color = "black", linewidth = 0.3) +
scale_color_manual(values = colorRampPalette(paleta)(n_sim),
guide = "none") +
labs(
title = "Movimiento Browniano Estándar",
subtitle = paste0(n_sim, " trayectorias simuladas | banda teórica ±1.96·√t al 95%"),
x = "Tiempo (t)", y = expression(W[t]),
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p1)
3.3 Verificación numérica
Comprobamos que, al final del horizonte, la media y la varianza simuladas se acercan a sus valores teóricos: \(E[W_T] = 0\) y \(Var(W_T) = T = 1\).
## Media simulada de W_T : 0.0471 (teórica: 0)## Varianza simulada de W_T : 0.7936 (teórica: 1.0)cat(sprintf("Trayectorias fuera de la banda al final: %d de %d\n",
sum(abs(W[N + 1, ]) > 1.96 * sqrt(T_final)), n_sim))## Trayectorias fuera de la banda al final: 1 de 503.4 Interpretación
- Las trayectorias son continuas pero muy irregulares: el Movimiento Browniano no es derivable en ningún punto. Por eso el cálculo tradicional no sirve y necesitamos el cálculo de Itô.
- La banda amarilla se abre como un embudo: la incertidumbre sobre el valor futuro crece con el horizonte.
- Aproximadamente el 95 % de las trayectorias termina dentro de la banda. Con solo 50 trayectorias las cifras simuladas se desvían algo de las teóricas; con más simulaciones se acercarían.
4 Movimiento Browniano Geométrico (GBM)
4.1 Teoría
El Movimiento Browniano estándar tiene dos problemas para modelar precios: puede ser negativo y sus cambios son absolutos, no porcentuales. Un precio de S/ 100 y uno de S/ 10 no deberían tener shocks del mismo tamaño en soles.
El Movimiento Browniano Geométrico resuelve ambos problemas modelando los rendimientos (cambios porcentuales):
\[ \frac{dS_t}{S_t} = \mu\,dt + \sigma\,dW_t \qquad\Longleftrightarrow\qquad dS_t = \mu S_t\,dt + \sigma S_t\,dW_t \]
donde:
| Parámetro | Significado | Valor usado |
|---|---|---|
| \(S_0\) | Precio inicial | 100 |
| \(\mu\) | Rendimiento esperado anual (drift) | 10 % |
| \(\sigma\) | Volatilidad anual | 25 % |
La solución de esta ecuación, que se obtiene con el Lema de Itô (sección 5), es:
\[ S_t = S_0 \exp\!\left[\left(\mu - \tfrac{1}{2}\sigma^2\right)t + \sigma W_t\right] \]
Como la función exponencial siempre es positiva, el precio nunca es negativo. Este es el modelo que usa Black-Scholes.
4.1.1 ¿Cómo se simula?
Aprovechando la solución exacta, cada paso se calcula como:
\[ S_{t+\Delta t} = S_t \exp\!\left[\left(\mu - \tfrac{1}{2}\sigma^2\right)\Delta t + \sigma\sqrt{\Delta t}\,Z\right] \]
Esta forma de simular es exacta: no introduce error de discretización, a diferencia de otros métodos aproximados.
4.1.2 Media y percentiles teóricos
- Valor esperado: \(E[S_t] = S_0\,e^{\mu t}\).
- Percentil \(p\): \(S_0 \exp\!\left[(\mu - \tfrac{1}{2}\sigma^2)t + \sigma\sqrt{t}\,z_p\right]\), donde \(z_p\) es el cuantil de la normal estándar (\(z_{0.05} = -1.645\) y \(z_{0.95} = 1.645\)).
4.2 Código
S0 <- 100 # Precio inicial
mu <- 0.10 # Drift anual
sigma <- 0.25 # Volatilidad anual
simular_gbm <- function(S0, mu, sigma, n_sim, N, dt) {
S <- matrix(NA, nrow = N + 1, ncol = n_sim)
S[1, ] <- S0
for (i in 2:(N + 1)) {
Z <- rnorm(n_sim)
S[i, ] <- S[i - 1, ] * exp((mu - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z)
}
return(S)
}
S <- simular_gbm(S0, mu, sigma, n_sim, N, dt)
df_S <- data.frame(
t = rep(t_grid, times = n_sim),
S = as.vector(S),
id = factor(rep(1:n_sim, each = N + 1))
)
# Media y percentiles teóricos
media_teo <- S0 * exp(mu * t_grid)
q05 <- S0 * exp((mu - 0.5 * sigma^2) * t_grid + sigma * sqrt(t_grid) * qnorm(0.05))
q95 <- S0 * exp((mu - 0.5 * sigma^2) * t_grid + sigma * sqrt(t_grid) * qnorm(0.95))
bandas_S <- data.frame(t = t_grid, media = media_teo, q05 = q05, q95 = q95)p2 <- ggplot() +
geom_line(data = df_S, aes(x = t, y = S, group = id),
color = "#328CC1", alpha = 0.3, linewidth = 0.5) +
geom_ribbon(data = bandas_S, aes(x = t, ymin = q05, ymax = q95),
fill = "#0B3C5D", alpha = 0.12) +
geom_line(data = bandas_S, aes(x = t, y = media),
color = "#B82601", linewidth = 1.3) +
geom_line(data = bandas_S, aes(x = t, y = q05),
color = "#D9B310", linetype = "dashed", linewidth = 0.8) +
geom_line(data = bandas_S, aes(x = t, y = q95),
color = "#D9B310", linetype = "dashed", linewidth = 0.8) +
labs(
title = "Movimiento Browniano Geométrico (GBM)",
subtitle = bquote(S[0] == .(S0) ~ "|" ~ mu == .(mu) ~ "|" ~ sigma == .(sigma) ~
"| línea roja: E[S_t] | bandas amarillas: p5 y p95"),
x = "Tiempo (t)", y = expression(S[t]),
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p2)
4.3 Interpretación
- Todas las trayectorias se mantienen positivas, como debe ocurrir con un precio.
- La banda entre los percentiles 5 y 95 es asimétrica: se abre más hacia arriba que hacia abajo. Esto refleja que un precio puede subir sin límite, pero no puede caer por debajo de cero (la distribución es lognormal).
- La media teórica (línea roja) crece al 10 % anual: \(E[S_1] = 100\,e^{0.10} \approx 110.5\).
- Dato importante: la media está por encima de la mediana. La mayoría de trayectorias termina por debajo del valor esperado, porque unas pocas trayectorias con subidas muy fuertes elevan el promedio.
5 Lema de Itô
5.1 Teoría: ¿por qué necesitamos una “regla de la cadena” distinta?
En el cálculo tradicional, si \(y = f(x)\), la regla de la cadena dice \(dy = f'(x)\,dx\), porque los términos de segundo orden, como \((dx)^2\), son despreciables.
Con el Movimiento Browniano esto ya no es cierto. Como \(\Delta W \sim \mathcal{N}(0, \Delta t)\), su cuadrado tiene valor esperado \(\Delta t\), del mismo orden que el tiempo. En el límite se cumple la famosa regla:
\[ (dW_t)^2 = dt, \qquad dt \cdot dW_t = 0, \qquad (dt)^2 = 0 \]
Por eso, al desarrollar una función de un proceso estocástico, el término de segundo orden no desaparece. Esto da lugar al Lema de Itô, la regla de la cadena del cálculo estocástico.
5.2 Enunciado
Si \(X_t\) sigue \(dX_t = a\,dt + b\,dW_t\) y \(f(X_t, t)\) es una función dos veces derivable, entonces:
\[ df = \left(\frac{\partial f}{\partial t} + a\frac{\partial f}{\partial x} + \frac{1}{2}b^2\frac{\partial^2 f}{\partial x^2}\right)dt + b\,\frac{\partial f}{\partial x}\,dW_t \]
El término \(\tfrac{1}{2}b^2 f''\) es la corrección de Itô: no existe en el cálculo tradicional.
5.3 Aplicación paso a paso: el logaritmo del precio
Aplicamos el lema a \(f(S) = \ln S\), con \(dS = \mu S\,dt + \sigma S\,dW\):
Paso 1. Derivadas: \[ \frac{\partial f}{\partial t} = 0, \qquad \frac{\partial f}{\partial S} = \frac{1}{S}, \qquad \frac{\partial^2 f}{\partial S^2} = -\frac{1}{S^2} \]
Paso 2. Reemplazamos con \(a = \mu S\) y \(b = \sigma S\): \[ d(\ln S) = \left(0 + \mu S\cdot\frac{1}{S} + \frac{1}{2}\sigma^2 S^2\cdot\left(-\frac{1}{S^2}\right)\right)dt + \sigma S\cdot\frac{1}{S}\,dW \]
Paso 3. Simplificamos: \[ d(\ln S_t) = \left(\mu - \tfrac{1}{2}\sigma^2\right)dt + \sigma\,dW_t \]
Paso 4. Integramos entre \(0\) y \(T\): \[ \ln S_T = \ln S_0 + \left(\mu - \tfrac{1}{2}\sigma^2\right)T + \sigma W_T \]
Conclusión: como \(W_T \sim \mathcal{N}(0, T)\),
\[ \ln S_T \sim \mathcal{N}\!\left(\ln S_0 + \left(\mu - \tfrac{1}{2}\sigma^2\right)T,\;\; \sigma^2 T\right) \]
El precio sigue una distribución lognormal. Además, aparece el término \(-\tfrac{1}{2}\sigma^2\), llamado arrastre de volatilidad (volatility drag): a mayor volatilidad, menor es el crecimiento típico de una inversión, aunque su rendimiento esperado sea el mismo.
💡 Ejemplo del arrastre de volatilidad: si una inversión sube 50 % y luego cae 50 %, el rendimiento promedio aritmético es 0 %, pero el capital pasa de 100 a 150 y luego a 75. La volatilidad “castiga” el crecimiento compuesto.
5.4 Código: verificación con las simulaciones del GBM
Tomamos los precios finales \(S_T\) de las 50 trayectorias del GBM, calculamos su logaritmo y comparamos el histograma con la densidad normal teórica.
log_St <- log(S[N + 1, ]) # log de los precios finales
log_S0 <- log(S0)
df_log <- data.frame(log_S = log_St)
# Parámetros de la distribución teórica
media_log <- log_S0 + (mu - 0.5 * sigma^2) * T_final
sd_log <- sigma * sqrt(T_final)
cat(sprintf("Media teórica de log(S_T): %.4f | simulada: %.4f\n", media_log, mean(log_St)))## Media teórica de log(S_T): 4.6739 | simulada: 4.6718## Desv. estándar teórica : 0.2500 | simulada: 0.2431p3 <- ggplot(df_log, aes(x = log_S)) +
geom_histogram(aes(y = after_stat(density)),
bins = 20, fill = "#0B3C5D", color = "white", alpha = 0.85) +
stat_function(fun = dnorm, args = list(mean = media_log, sd = sd_log),
color = "#B82601", linewidth = 1.2) +
labs(
title = "Distribución de log(S_T) — Verificación del Lema de Itô",
subtitle = expression(log(S[T]) %~% N * "(" * log(S[0]) + (mu - sigma^2/2)*T * "," ~ sigma^2*T * ")"),
x = expression(log(S[T])), y = "Densidad",
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p3)
5.4.1 Prueba formal de normalidad
La prueba de Shapiro-Wilk contrasta la hipótesis nula de que los datos provienen de una distribución normal. Un p-valor mayor a 0.05 indica que no hay evidencia para rechazar la normalidad.
##
## Shapiro-Wilk normality test
##
## data: log_St
## W = 0.97275, p-value = 0.29835.5 Interpretación
- El histograma de \(\ln S_T\) se ajusta a la curva normal teórica (línea roja), lo que confirma el resultado del Lema de Itô.
- Con 50 observaciones el histograma es irregular; si se aumenta
n_sima 5 000, el ajuste es casi perfecto.
6 Proceso de Ornstein-Uhlenbeck
6.1 Teoría
No todas las variables económicas “vagan” libremente como un precio. Algunas tienden a regresar a un nivel de equilibrio: la tasa de interés, la inflación alrededor de su meta, el tipo de cambio real o el diferencial de precios entre dos activos relacionados. Para ellas se usa el proceso de Ornstein-Uhlenbeck (OU):
\[ dX_t = \theta\,(\mu - X_t)\,dt + \sigma\,dW_t \]
| Parámetro | Significado | Valor usado |
|---|---|---|
| \(\theta\) | Velocidad de reversión a la media | 2.0 |
| \(\mu\) | Nivel de equilibrio de largo plazo | 0 |
| \(\sigma\) | Volatilidad | 0.5 |
| \(X_0\) | Valor inicial | 2 |
6.1.1 ¿Cómo funciona la reversión a la media?
Observa el término de tendencia \(\theta(\mu - X_t)\):
- Si \(X_t > \mu\), el término es negativo y empuja al proceso hacia abajo.
- Si \(X_t < \mu\), el término es positivo y lo empuja hacia arriba.
- Mientras más lejos esté del equilibrio, más fuerte es el empuje, como un resorte.
6.1.2 Resultados teóricos importantes
- Media: \(E[X_t] = \mu + (X_0 - \mu)\,e^{-\theta t}\). La distancia al equilibrio se reduce exponencialmente.
- Varianza: \(Var(X_t) = \dfrac{\sigma^2}{2\theta}\left(1 - e^{-2\theta t}\right)\). A diferencia del Browniano, no crece indefinidamente: se estabiliza en \(\sigma^2/(2\theta)\).
- Vida media: el tiempo que tarda en cerrarse la mitad de la distancia al equilibrio es \(\dfrac{\ln 2}{\theta}\).
💡 Aplicación: el modelo de Vasicek (1977) para tasas de interés es exactamente un proceso OU. Si la tasa está muy por encima de su promedio histórico, el modelo predice que tenderá a bajar.
6.1.3 ¿Cómo se simula? El método de Euler-Maruyama
Reemplazamos los diferenciales por incrementos discretos:
\[ X_{t+\Delta t} = X_t + \theta(\mu - X_t)\,\Delta t + \sigma\sqrt{\Delta t}\,Z \]
Este es el método de Euler-Maruyama, la forma más sencilla de simular cualquier ecuación diferencial estocástica. Es aproximado, pero con \(\Delta t = 0.001\) el error es despreciable.
6.2 Código
theta <- 2.0 # Velocidad de reversión
mu_ou <- 0.0 # Nivel de largo plazo
sigma_ou <- 0.5 # Volatilidad
simular_ou <- function(x0, theta, mu, sigma, n_sim, N, dt) {
X <- matrix(NA, nrow = N + 1, ncol = n_sim)
X[1, ] <- x0
for (i in 2:(N + 1)) {
Z <- rnorm(n_sim)
X[i, ] <- X[i - 1, ] + theta * (mu - X[i - 1, ]) * dt + sigma * sqrt(dt) * Z
}
return(X)
}
X <- simular_ou(x0 = 2, theta, mu_ou, sigma_ou, n_sim, N, dt)
df_X <- data.frame(
t = rep(t_grid, times = n_sim),
X = as.vector(X),
id = factor(rep(1:n_sim, each = N + 1))
)p4 <- ggplot(df_X, aes(x = t, y = X, group = id)) +
geom_line(color = "#4A7C59", alpha = 0.35, linewidth = 0.5) +
geom_hline(yintercept = mu_ou, color = "#B82601",
linetype = "dashed", linewidth = 1) +
annotate("label", x = 0.05, y = mu_ou + 0.15,
label = "Nivel~de~largo~plazo~mu", parse = TRUE,
hjust = 0, fill = "white", color = "#B82601", size = 3.5) +
labs(
title = "Proceso de Ornstein-Uhlenbeck (Reversión a la Media)",
subtitle = bquote(dX[t] == theta*(mu - X[t])*dt + sigma*"d"*W[t] ~ "|" ~
theta == .(theta) ~ "," ~ sigma == .(sigma_ou)),
x = "Tiempo (t)", y = expression(X[t]),
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p4)
6.3 Verificación numérica
x0 <- 2
media_teo_T <- mu_ou + (x0 - mu_ou) * exp(-theta * T_final)
var_teo_T <- sigma_ou^2 / (2 * theta) * (1 - exp(-2 * theta * T_final))
cat(sprintf("Vida media (ln2/θ) : %.3f años (~%.0f días)\n", log(2) / theta, 365 * log(2) / theta))## Vida media (ln2/θ) : 0.347 años (~126 días)## Media en T -> teórica: 0.2707 | simulada: 0.3276## Varianza en T -> teórica: 0.0614 | simulada: 0.0528## Varianza de largo plazo σ²/(2θ): 0.06256.4 Interpretación
- Todas las trayectorias parten de \(X_0 = 2\) y caen rápidamente hacia el nivel de equilibrio \(\mu = 0\).
- Con \(\theta = 2\), la vida media es de unos 0.35 años (aproximadamente 4 meses): en ese tiempo se cierra la mitad de la distancia al equilibrio.
- Una vez cerca del equilibrio, las trayectorias oscilan alrededor de él dentro de una banda estable, en lugar de alejarse como en el Movimiento Browniano.
7 Puente Browniano
7.1 Teoría
Un Puente Browniano es un Movimiento Browniano condicionado a terminar en un valor conocido. En su versión estándar, comienza en cero y también termina en cero:
\[ B_0 = 0, \qquad B_T = 0 \]
7.1.1 Construcción
Se obtiene a partir de un Browniano \(W_t\) restándole, en cada instante, la parte proporcional de su valor final:
\[ B_t = W_t - \frac{t}{T}\,W_T \]
- En \(t = 0\): \(B_0 = W_0 - 0 = 0\).
- En \(t = T\): \(B_T = W_T - W_T = 0\).
7.1.2 Propiedades
- Media: \(E[B_t] = 0\).
- Varianza: \(Var(B_t) = \dfrac{t\,(T - t)}{T}\). Es cero en los extremos y máxima a la mitad del intervalo (\(t = T/2\)), donde vale \(T/4\).
💡 Intuición: es como una cuerda atada en ambos extremos que se sacude al azar: puede moverse mucho en el centro, pero no en los puntos donde está sujeta.
7.1.3 ¿Para qué sirve?
- Simulación eficiente: permite generar primero el valor final de un precio y luego “rellenar” la trayectoria intermedia.
- Valoración de opciones que dependen de la trayectoria, como las opciones barrera, que se activan si el precio toca cierto nivel.
- Estadística: es la base de la prueba de Kolmogorov-Smirnov.
- Finanzas: modela activos con valor final conocido, como un bono que al vencimiento paga su valor nominal.
7.2 Código
set.seed(777)
n_puente <- 20
W_full <- simular_bm(n_puente, N, dt)
W_T <- W_full[N + 1, ]
# B_t = W_t - (t/T) W_T, aplicado a todas las trayectorias a la vez
B <- W_full - outer(t_grid, W_T, "*") / T_final
df_B <- data.frame(
t = rep(t_grid, times = n_puente),
B = as.vector(B),
id = factor(rep(1:n_puente, each = N + 1))
)La función outer(t_grid, W_T, "*") construye una matriz
donde cada elemento es \(t \cdot W_T\).
Así se aplica la fórmula del puente a las 20 trayectorias de una sola
vez, sin usar bucles.
p5 <- ggplot(df_B, aes(x = t, y = B, group = id, color = id)) +
geom_line(alpha = 0.75, linewidth = 0.7) +
geom_hline(yintercept = 0, color = "grey30", linewidth = 0.4) +
geom_point(data = data.frame(t = c(0, T_final), B = c(0, 0)),
aes(x = t, y = B), inherit.aes = FALSE,
color = "#B82601", size = 3) +
scale_color_manual(values = colorRampPalette(paleta)(n_puente),
guide = "none") +
labs(
title = "Puente Browniano",
subtitle = "Trayectorias condicionadas a comenzar y terminar en cero",
x = "Tiempo (t)", y = expression(B[t]),
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p5)
7.3 Verificación numérica
## Valor en T (debe ser 0) : 0.00e+00cat(sprintf("Varianza en t = T/2 -> teórica: %.4f | simulada: %.4f\n",
(T_final / 2) * (T_final / 2) / T_final, var(B[mitad, ])))## Varianza en t = T/2 -> teórica: 0.2500 | simulada: 0.16707.4 Interpretación
- Todas las trayectorias comienzan y terminan exactamente en cero (puntos rojos).
- La dispersión es máxima alrededor de \(t = 0.5\) y se reduce hacia los extremos, tal como indica la fórmula de la varianza.
8 Valoración de opciones: Monte Carlo vs. Black-Scholes
8.1 Teoría
8.1.1 ¿Qué es una opción call europea?
Es un contrato que da a su poseedor el derecho, pero no la obligación, de comprar un activo a un precio fijo \(K\) (precio de ejercicio o strike) en una fecha futura \(T\).
- Si al vencimiento \(S_T > K\), conviene ejercerla: se compra a \(K\) algo que vale \(S_T\) y se gana \(S_T - K\).
- Si \(S_T \le K\), no conviene ejercerla y se gana cero.
Por eso el pago (payoff) es:
\[ \text{Payoff} = \max(S_T - K,\; 0) \]
La pregunta es: ¿cuánto debería pagarse hoy por ese derecho?
8.1.2 La valoración neutral al riesgo
La teoría financiera demuestra que, si no existen oportunidades de arbitraje, el precio de la opción es el valor esperado de su pago, descontado a la tasa libre de riesgo, calculado bajo una medida de probabilidad especial llamada medida neutral al riesgo (\(\mathbb{Q}\)):
\[ C_0 = e^{-rT}\;E^{\mathbb{Q}}\!\left[\max(S_T - K, 0)\right] \]
Bajo esta medida, el rendimiento esperado \(\mu\) del activo se reemplaza por la tasa libre de riesgo \(r\):
\[ S_T = S_0 \exp\!\left[\left(r - \tfrac{1}{2}\sigma^2\right)T + \sigma\sqrt{T}\,Z\right] \]
💡 Intuición: el precio de la opción no depende de cuánto crea cada inversionista que subirá la acción. Esto se debe a que la opción puede replicarse combinando la acción y un bono, y esa réplica no depende de \(\mu\).
8.1.3 Método 1: simulación de Monte Carlo
- Simular \(n\) precios finales \(S_T^{(i)}\) bajo la medida neutral al riesgo.
- Calcular el pago de cada escenario: \(\max(S_T^{(i)} - K, 0)\).
- Promediar los pagos y descontarlos:
\[ \hat{C}_{MC} = e^{-rT}\,\frac{1}{n}\sum_{i=1}^{n}\max\!\left(S_T^{(i)} - K,\,0\right) \]
- Medir la precisión con el error estándar: \(\;EE = e^{-rT}\,\dfrac{s_{\text{payoff}}}{\sqrt{n}}\). El intervalo de confianza al 95 % es \(\hat{C}_{MC} \pm 1.96\,EE\).
8.1.4 Método 2: la fórmula cerrada de Black-Scholes (1973)
Para una call europea existe una solución exacta:
\[ C_{BS} = S_0\,\Phi(d_1) - K e^{-rT}\,\Phi(d_2) \]
\[ d_1 = \frac{\ln(S_0/K) + \left(r + \tfrac{1}{2}\sigma^2\right)T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T} \]
donde \(\Phi(\cdot)\) es la función de distribución acumulada de la normal estándar. El término \(\Phi(d_2)\) es la probabilidad (neutral al riesgo) de que la opción termine siendo ejercida.
¿Por qué comparar ambos métodos? Black-Scholes solo tiene fórmula cerrada para casos sencillos. Monte Carlo funciona para opciones mucho más complejas (asiáticas, barrera, con varios activos). Si Monte Carlo reproduce bien el caso conocido, podemos confiar en él para los casos sin fórmula.
8.2 Código
| Parámetro | Valor |
|---|---|
| Precio actual \(S_0\) | 100 |
| Precio de ejercicio \(K\) | 100 (opción “en el dinero”) |
| Tasa libre de riesgo \(r\) | 5 % anual |
| Volatilidad \(\sigma\) | 20 % anual |
| Vencimiento \(T\) | 1 año |
| Simulaciones | 10 000 |
# Parámetros
S0_op <- 100
K <- 100 # Strike
r <- 0.05 # Tasa libre de riesgo
sigma_op <- 0.20
T_op <- 1
n_paths <- 10000
# --- Simulación Monte Carlo ---
Z <- rnorm(n_paths)
ST <- S0_op * exp((r - 0.5 * sigma_op^2) * T_op + sigma_op * sqrt(T_op) * Z)
payoff_call <- pmax(ST - K, 0)
precio_MC <- exp(-r * T_op) * mean(payoff_call)
error_MC <- exp(-r * T_op) * sd(payoff_call) / sqrt(n_paths)
# --- Fórmula cerrada Black-Scholes ---
d1 <- (log(S0_op / K) + (r + 0.5 * sigma_op^2) * T_op) / (sigma_op * sqrt(T_op))
d2 <- d1 - sigma_op * sqrt(T_op)
precio_BS <- S0_op * pnorm(d1) - K * exp(-r * T_op) * pnorm(d2)
cat("\n=============================================\n")##
## =============================================## PRICING DE OPCIÓN CALL EUROPEA## =============================================## Precio Monte Carlo : 10.4169 (± 0.2873)## Precio Black-Scholes: 10.4506## Diferencia absoluta : 0.033687## =============================================df_MC <- data.frame(ST = ST, payoff = payoff_call)
p6 <- ggplot(df_MC, aes(x = ST)) +
geom_histogram(aes(y = after_stat(density)),
bins = 60, fill = "#328CC1", color = "white", alpha = 0.8) +
geom_vline(xintercept = K, color = "#B82601",
linetype = "dashed", linewidth = 1) +
annotate("label", x = K + 3, y = Inf,
label = paste0("Strike K = ", K),
vjust = 1.5, hjust = 0, fill = "white",
color = "#B82601", size = 3.5) +
labs(
title = "Distribución de S_T bajo medida neutral al riesgo",
subtitle = sprintf("Precio Call MC = %.4f | Precio BS = %.4f | N = %s",
precio_MC, precio_BS, format(n_paths, big.mark = ",")),
x = expression(S[T]), y = "Densidad",
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p6)
8.3 Interpretación
- El precio de Black-Scholes es aproximadamente 10.45: por el derecho a comprar dentro de un año, a 100, un activo que hoy vale 100, se paga hoy alrededor de 10.45.
- El precio de Monte Carlo cae muy cerca de ese valor y el precio exacto queda dentro de su intervalo de confianza, lo que valida el método.
- En el histograma, solo los escenarios a la derecha de la línea roja (\(S_T > K\)) generan pago. La opción vale porque ofrece la ganancia de las subidas sin la pérdida de las caídas.
- El histograma está sesgado a la derecha: es la distribución lognormal derivada con el Lema de Itô.
9 Convergencia de Monte Carlo
9.1 Teoría
¿Por qué funciona Monte Carlo? Por dos resultados fundamentales de la estadística:
- Ley de los Grandes Números: el promedio de muchas observaciones independientes converge a su valor esperado. Por eso \(\hat{C}_{MC} \to C_{BS}\) cuando \(n \to \infty\).
- Teorema del Límite Central: el error de estimación se distribuye aproximadamente normal, con una desviación estándar proporcional a \(1/\sqrt{n}\).
Esta última propiedad tiene una consecuencia práctica muy importante:
📌 Para reducir el error a la mitad, hay que multiplicar el número de simulaciones por cuatro. Pasar de 10 000 a 40 000 simulaciones solo duplica la precisión.
9.2 Código
Estimamos el precio de la call con tamaños de muestra crecientes, desde 100 hasta 20 000 simulaciones. Cada punto del gráfico es una estimación independiente con \(n\) simulaciones.
n_seq <- seq(100, 20000, by = 100)
precios_acum <- sapply(n_seq, function(n) {
z <- rnorm(n)
sT <- S0_op * exp((r - 0.5 * sigma_op^2) * T_op + sigma_op * sqrt(T_op) * z)
exp(-r * T_op) * mean(pmax(sT - K, 0))
})
df_conv <- data.frame(n = n_seq, precio = precios_acum)p7 <- ggplot(df_conv, aes(x = n, y = precio)) +
geom_line(color = "#0B3C5D", linewidth = 0.9) +
geom_hline(yintercept = precio_BS, color = "#B82601",
linetype = "dashed", linewidth = 1) +
annotate("label", x = max(n_seq), y = precio_BS,
label = sprintf("Black-Scholes = %.4f", precio_BS),
hjust = 1.05, vjust = -0.6, fill = "white",
color = "#B82601", size = 3.5) +
scale_x_continuous(labels = scales::comma) +
labs(
title = "Convergencia del Precio Monte Carlo",
subtitle = "Aproximación al precio analítico de Black-Scholes por Ley de Grandes Números",
x = "Número de simulaciones", y = "Precio estimado del Call",
caption = "Elaborado por: Economista Jeel Cueva"
) +
theme_jeel()
print(p7)
9.3 Interpretación
- Con pocas simulaciones, las estimaciones oscilan mucho alrededor del precio verdadero.
- A medida que \(n\) aumenta, las oscilaciones se reducen y la línea azul se acerca a la línea roja de Black-Scholes.
- La reducción es lenta, al ritmo de \(1/\sqrt{n}\). Por eso en la práctica se usan técnicas de reducción de varianza (variables antitéticas, variables de control) que mejoran la precisión sin aumentar tanto las simulaciones.
10 Exportación de gráficos
Finalmente, guardamos los siete gráficos en alta resolución (300 dpi), listos para presentaciones, informes o publicaciones en redes.
dir.create("simulaciones_estocasticas", showWarnings = FALSE)
ggsave("simulaciones_estocasticas/01_browniano_estandar.png", p1, width = 10, height = 6, dpi = 300, bg = "white")
ggsave("simulaciones_estocasticas/02_GBM.png", p2, width = 10, height = 6, dpi = 300, bg = "white")
ggsave("simulaciones_estocasticas/03_ito_lognormal.png", p3, width = 10, height = 6, dpi = 300, bg = "white")
ggsave("simulaciones_estocasticas/04_ornstein_uhlenbeck.png", p4, width = 10, height = 6, dpi = 300, bg = "white")
ggsave("simulaciones_estocasticas/05_puente_browniano.png", p5, width = 10, height = 6, dpi = 300, bg = "white")
ggsave("simulaciones_estocasticas/06_pricing_MC.png", p6, width = 10, height = 6, dpi = 300, bg = "white")
ggsave("simulaciones_estocasticas/07_convergencia_MC.png", p7, width = 10, height = 6, dpi = 300, bg = "white")
cat("Simulaciones completadas. Gráficos guardados en: simulaciones_estocasticas/\n")Este bloque tiene la opción
eval=FALSEpara que no se ejecute al generar el documento HTML. Para guardar los gráficos, ejecútalo manualmente en RStudio.
11 Conclusiones
| Proceso | Ecuación | Característica clave | Aplicación económica |
|---|---|---|---|
| Browniano estándar | \(dW_t\) | Ruido puro, varianza \(= t\) | Base de todos los modelos |
| Browniano Geométrico | \(dS = \mu S\,dt + \sigma S\,dW\) | Siempre positivo, lognormal | Precios de acciones, Black-Scholes |
| Ornstein-Uhlenbeck | \(dX = \theta(\mu - X)\,dt + \sigma\,dW\) | Reversión a la media | Tasas de interés, inflación, spreads |
| Puente Browniano | \(B_t = W_t - \tfrac{t}{T}W_T\) | Extremos fijos | Opciones barrera, bonos |
Los mensajes principales de este notebook son:
- La incertidumbre se puede modelar. El Movimiento Browniano permite representar matemáticamente el azar y construir sobre él modelos para variables económicas.
- El cálculo estocástico tiene sus propias reglas. Como \((dW)^2 = dt\), el Lema de Itô agrega un término de corrección que explica por qué la volatilidad reduce el crecimiento compuesto.
- No todas las variables se comportan igual. Los precios de activos se modelan mejor con un GBM; las tasas de interés y otras variables con equilibrio, con un proceso de Ornstein-Uhlenbeck.
- La simulación es una herramienta poderosa. Monte Carlo reproduce el precio exacto de Black-Scholes y puede extenderse a problemas donde no hay fórmula cerrada.
11.1 Referencias
- Black, F. y Scholes, M. (1973). The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81(3), 637–654.
- Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer.
- Hull, J. C. (2022). Options, Futures, and Other Derivatives (11.ª ed.). Pearson.
- Øksendal, B. (2003). Stochastic Differential Equations: An Introduction with Applications (6.ª ed.). Springer.
- Shreve, S. E. (2004). Stochastic Calculus for Finance II: Continuous-Time Models. Springer.
- Vasicek, O. (1977). An Equilibrium Characterization of the Term Structure. Journal of Financial Economics, 5(2), 177–188.
Si este material te fue útil, compártelo y déjame tus comentarios.