Simulación Estocástica de Monte Carlo Multivariante en Python: Factorización de Cholesky, Salto-Difusión de Merton y Gestión de Carteras (2010–2036)

1. Línea Base del Estudio Original en RStudio (Julio 2020) y sus 4 Limitaciones Estructurales

En julio de 2020 publicamos el Trabajo #10 en RStudio con el objetivo de superar la rigidez de las proyecciones financieras determinísticas —aquellas que extrapolan una única tasa compuesta fija— incorporando la naturaleza aleatoria de los mercados de capitales. Partiendo de las ponderaciones defensivas obtenidas en el Trabajo #4 sobre el Portafolio Tangente de Markowitz para seis corporaciones líderes del índice S&P 500 —Microsoft (MSFT, 4.10%), Johnson & Johnson (JNJ, 27.83%), Mastercard (MA, 7.41%), PepsiCo (PEP, 14.21%), Walmart (WMT, 23.02%) y McDonald's (MCD, 23.42%)— entre enero de 2010 y abril de 2020 (N = 123 meses), el modelo en R estimó una media mensual logarítmica μ̂₂₀₂₀ = 1.0794% y una desviación estándar mensual σ̂₂₀₂₀ = 3.2407%.

Con dichos estadísticos, aquel primer diseño ejecutó M = 100 trayectorias aleatorias a un horizonte de H = 120 meses (10 años, abril 2020 – abril 2030) mediante la instrucción recursiva cumprod(1 + rnorm(120, mean_port_return, stddev_port_return)) sobre una unidad monetaria inicial (W₀ = $1.00 USD), arrojando en esa corrida un crecimiento terminal mínimo de $1.48 USD, una mediana de $3.47 USD y un máximo muestral de $9.00 USD. Sin embargo, al someter ese diseño fundacional a los estándares contemporáneos de ingeniería cuantitativa institucional en 2026, emergen cuatro limitaciones matemáticas y econométricas críticas que esta investigación resuelve de raíz:

2. Auditoría Empírica Fuera de Muestra (Abril 2020 – Noviembre 2026): Realidad vs. Cono Estocástico

A diferencia de un ejercicio puramente teórico, hoy contamos con 6.53 años exactos de historia bursátil fuera de muestra transcurridos desde el cierre del estudio original: 1,641 sesiones diarias de negociación entre el 30 de abril de 2020 y el 10 de noviembre de 2026 (t = 78.36 meses dentro del horizonte original de 120 meses). Esta ventana abarca cinco regímenes macrofinancieros de altísima exigencia: la recuperación post-COVID (2020–2021), el choque inflacionario global y el ciclo restrictivo de la Reserva Federal en 2022 (donde el S&P 500 cedió -25.43%), la expansión estructural de la Inteligencia Artificial liderada por MSFT (2023–2024) y la reconfiguración comercial y de consumo defensivo de 2025–2026 liderada por WMT y MA.

Al reconstruir con precisión exacta (M = 10,000 réplicas) el cono probabilístico formulado en abril de 2020 y evaluar sus cuantiles teóricos exactamente en el mes t = 78.36 (noviembre de 2026), los umbrales proyectados por cada $1.00 USD invertido eran: Percentil 5% = $1.4015 ($14,015 por $10k), Percentil 25% = $1.8464, Mediana (P50) = $2.2403 ($22,403 por $10k), Media Teórica = $2.3241 ($23,241 por $10k), Percentil 75% = $2.7014 y Percentil 95% = $3.5430. La Tabla 2.1 confronta estos umbrales ex-ante contra el desempeño auditado de mercado:

3. Formulación Matemática y Cálculo Estocástico Multivariante

Para modelar rigurosamente un portafolio institucional compuesto por d = 6 activos riesgosos cuyos precios evolucionan en tiempo continuo sobre un espacio de probabilidad filtrado (Ω, &mathcal;F;, {&mathcal;F;_t}_t≥0, ℙ), abandonamos el escalar univariado y formulamos el sistema vectorial acoplado de ecuaciones diferenciales estocásticas (SDE).

3.1. Movimiento Browniano Geométrico Multivariante y Corrección del Lema de Itô

Sea S_t = (S_1,t, …, S_d,t)^⊤ el vector de cotizaciones ajustadas. Bajo difusión continua correlacionada, cada activo i ∈ {1, …, 6} obedece la SDE dS_i,t / S_i,t = μ_i dt + σ_i dW_i,t, donde 𝔼[dW_i,t dW_j,t] = ρ_ij dt. Aplicando el Lema de Itô a la transformación f(S_i,t) = ln(S_i,t), cuya segunda derivada f''(S) = −1/S² es estrictamente no nula, el término cuadrático ½f''(S_i,t)(dS_i,t)² = −½σ_i²dt genera la solución exacta libre de sesgo de discretización en intervalos Δt:

3.2. Factorización Triangular de Cholesky de la Matriz de Covarianza (6 × 6)

Dado que la matriz empírica de varianzas-covarianzas Σ ∈ ℝ^6×6 es simétrica y definida positiva, existe una única matriz triangular inferior L con entradas diagonales estrictamente positivas tal que Σ = LL^⊤. Al premultiplicar un vector de innovaciones ortogonales estándar Z_t (𝔼[Z_tZ_t^⊤] = I₆) por L, la covarianza del choque resultante reproduce exactamente la estructura empírica intersectorial:

3.3. Colas Pesadas Conjuntas (t-Student Multivariante) y Co-Dependencia de Extremos

En crisis sistémicas, las correlaciones efectivas se disparan simultáneamente: un choque extremo en tecnología co-ocurre con ventas forzadas de liquidez en consumo básico. Una cópula gaussiana posee dependencia asintótica de cola nula (λ_L = 0). Para capturar la leptocurtosis conjunta, modelamos el vector diario r_t ∈ ℝ⁶ bajo una Distribución t-Student Multivariante con ν > 2 grados de libertad y matriz de dispersión escalada S = ((ν−2)/ν)Σ (garantizando Cov(r_t) = Σ), estimando ν̂_MLE = 3.937 por Máxima Verosimilitud exacta:

3.4. Proceso de Salto-Difusión de Merton (1976) y Reducción de Varianza con Variables Antitéticas

Para incorporar discontinuidades macroeconómicas súbitas (guerras tarifarias, cierres sanitarios o sorpresas de política monetaria), acoplamos la difusión continua multivariante con un Proceso Compuesto de Poisson N_t ∼ Poisson(λt) de intensidad anual λ, donde cada salto multiplica el valor de la cartera por e^Y_j con Y_j ∼ &mathcal;N;(μ_J, σ_J²) y compensador de deriva insesgado κ = 𝔼[e^Y_j − 1] = exp(μ_J + ½σ_J²) − 1:

Dado que la riqueza terminal compuesta g(Z) es estrictamente monótona respecto a las innovaciones gaussianas acumuladas, el Teorema de Asociación de Esary-Proschan-Walkup garantiza que Corr[g(Z), g(−Z)] < 0. Al construir M/2 = 5,000 matrices independientes Z^(m) y reflejar su signo para obtener las restantes 5,000 trayectorias espejo −Z^(m), nuestro motor reduce empíricamente la varianza del estimador de Monte Carlo en un 80.26% a igual costo computacional, equivalente a simular más de 50,000 trayectorias estándar i.i.d.

4. Arquitectura Reproducible y Vectorización Tensorial en Python

Todo el flujo econométrico y estocástico ha sido programado en un módulo de producción autocontenido (models/monte_carlo_portafolios_2026.py) utilizando operaciones tensoriales de álgebra lineal en NumPy (np.einsum y np.linalg.cholesky) y optimización no lineal en SciPy, prescindiendo por completo de bucles iterativos lentos sobre trayectorias. A continuación se documentan los cuatro bloques algorítmicos centrales:

Paso 4.1. Ingesta Auditada (2010–2026) y Verificación Fuera de Muestra

import numpy as np
import pandas as pd
from scipy import optimize, special, stats

ASSETS = ["MSFT", "JNJ", "MA", "PEP", "WMT", "MCD"]
WEIGHTS_2020 = np.array([0.0410, 0.2783, 0.0741, 0.1421, 0.2302, 0.2342])
WEIGHTS_2020 = WEIGHTS_2020 / WEIGHTS_2020.sum()
RF_ANNUAL = 0.025  # Tasa libre de riesgo estructural T-Bill (2.50%)

# Carga de las 4,240 sesiones diarias ajustadas (04-Ene-2010 al 10-Nov-2026)
prices = pd.read_csv("data/sp500_montecarlo_2010_2026.csv", index_col=0, parse_dates=True)
daily_log = np.log(prices[ASSETS] / prices[ASSETS].shift(1)).dropna()
daily_simp = prices[ASSETS].pct_change().dropna()

# Auditoría Fuera de Muestra (30-Abr-2020 al 10-Nov-2026: 1,641 sesiones = 78.36 meses)
post_prices = prices.loc["2020-04-30":, ASSETS]
rel_prices = post_prices / post_prices.iloc[0]
bh_curve = rel_prices.values @ WEIGHTS_2020  # Buy & Hold ($1.00 -> $2.3135 USD)

Paso 4.2. Factorización Triangular de Cholesky 6D y Estimación MLE Student-t Multivariante

def fit_multivariate_tail_and_cholesky(daily_log_df: pd.DataFrame) -> dict:
    """Factorización triangular de Cholesky 6D y MLE t-Student Multivariante."""
    x = daily_log_df.values
    n_obs, d = x.shape
    mu_d = np.mean(x, axis=0)
    cov_d = np.cov(x, rowvar=False, ddof=1)
    chol_d = np.linalg.cholesky(cov_d)  # Sigma = L @ L.T

    # Distancias cuadráticas de Mahalanobis D_t^2 = (r_t - mu)' Sigma^{-1} (r_t - mu)
    inv_cov = np.linalg.inv(cov_d)
    _, logdet_cov = np.linalg.slogdet(cov_d)
    diff = x - mu_d
    maha_sq_raw = np.einsum("ti,ij,tj->t", diff, inv_cov, diff)

    def neg_loglik_mvt(params):
        nu = params[0]
        if nu <= 2.05:
            return 1e12
        scale_factor = (nu - 2.0) / nu
        logdet_s = logdet_cov + d * np.log(scale_factor)
        maha_s = maha_sq_raw / scale_factor
        ll = (
            n_obs * special.gammaln((nu + d) / 2.0)
            - n_obs * special.gammaln(nu / 2.0)
            - 0.5 * n_obs * (d * np.log(nu * np.pi) + logdet_s)
            - 0.5 * (nu + d) * np.sum(np.log1p(maha_s / nu))
        )
        return -ll

    opt = optimize.minimize(neg_loglik_mvt, x0=[4.5], bounds=[(2.1, 80.0)], method="L-BFGS-B")
    nu_hat = float(opt.x[0])
    ll_mvt = float(-opt.fun)
    ll_mvn = float(-0.5 * n_obs * (d * np.log(2.0 * np.pi) + logdet_cov) - 0.5 * np.sum(maha_sq_raw))
    return {"chol_d": chol_d, "nu_mvt_mle": nu_hat, "lr_stat": 2.0 * (ll_mvt - ll_mvn)}

Paso 4.3. Simulación Monte Carlo Transversal (15,000 Carteras Dirichlet) + Solver SLSQP

def optimize_cross_sectional_markowitz(daily_simp_df: pd.DataFrame, n_portfolios: int = 15_000):
    """Simulación Monte Carlo Transversal de Dirichlet + Optimización Cuadrática SLSQP."""
    rng = np.random.default_rng(20261111)
    mu_ann = daily_simp_df.mean().values * 252.0
    cov_ann = daily_simp_df.cov().values * 252.0
    n_assets = len(mu_ann)

    # 1. Nube estocástica de 15,000 vectores sobre el símplex w_i >= 0, sum(w_i) = 1
    w_rand = rng.dirichlet(np.ones(n_assets), size=n_portfolios)
    ret_cloud = w_rand @ mu_ann
    vol_cloud = np.sqrt(np.einsum("mi,ij,mj->m", w_rand, cov_ann, w_rand))
    sr_cloud = (ret_cloud - RF_ANNUAL) / vol_cloud

    # 2. Solución exacta del Portafolio Tangente 2026 mediante programación cuadrática SLSQP
    bounds = tuple((0.0, 1.0) for _ in range(n_assets))
    cons_sum = {"type": "eq", "fun": lambda w: np.sum(w) - 1.0}
    res_tan = optimize.minimize(
        lambda w: -(float(w @ mu_ann) - RF_ANNUAL) / np.sqrt(float(w @ cov_ann @ w)),
        x0=np.repeat(1.0 / n_assets, n_assets),
        method="SLSQP", bounds=bounds, constraints=[cons_sum],
    )
    return res_tan.x, ret_cloud, vol_cloud, sr_cloud

Paso 4.4. Motor Tensorial 3D de Salto-Difusión de Merton con Variables Antitéticas

def simulate_merton_jump_antithetic(mu_d, chol_d, weights, lam_d, mu_j, sd_j, n_sims=10_000, n_days=2520):
    """Motor acoplado GBM Multivariante (Cholesky + Antitéticas) + Saltos Compuestos de Poisson."""
    rng = np.random.default_rng(20261111)
    half = n_sims // 2
    # Pares antitéticos (Z, -Z) -> Reducción de varianza verificada del 80.26%
    z_half = rng.standard_normal((half, n_days, len(weights)))
    z_all = np.concatenate([z_half, -z_half], axis=0)

    # Difusión continua correlacionada + compensación del salto esperado kappa = E[e^Y - 1]
    r_diff_log = mu_d[None, None, :] + np.einsum("mtd,kd->mtk", z_all, chol_d)
    r_port_diff = np.expm1(r_diff_log) @ weights

    kappa = np.exp(mu_j + 0.5 * sd_j ** 2) - 1.0
    n_jumps = rng.poisson(lam_d, size=(n_sims, n_days))
    j_shock = np.where(n_jumps > 0, n_jumps * mu_j + np.sqrt(n_jumps) * sd_j * rng.standard_normal((n_sims, n_days)), 0.0)
    r_total = np.maximum(r_port_diff - lam_d * kappa + np.expm1(j_shock), -0.95)

    # Trayectorias diarias acumuladas (M=10,000 x T=2,520) y Maximum Drawdown exacto
    wealth_daily = np.hstack([np.ones((n_sims, 1)), np.cumprod(1.0 + r_total, axis=1)])
    running_peaks = np.maximum.accumulate(wealth_daily, axis=1)
    mdd_dist = np.min(wealth_daily / running_peaks - 1.0, axis=1)
    return wealth_daily, mdd_dist

5. Laboratorio Visual Interactivo del Modelo Monte Carlo (6 Gráficos Dinámicos)

En sustitución de las vistas estáticas de alta varianza muestral del post original de 2020, desplegamos a continuación una suite interactiva de seis visualizaciones dinámicas en Plotly.js alimentadas con las 4,240 sesiones empíricas y las 10,000 trayectorias estocásticas auditadas. Utiliza los botones superiores de cada lienzo para conmutar escalas monetarias, procesos generadores de datos y regímenes de correlación.

6. Interpretación Econométrica: Re-Optimización de Markowitz 2026, Matriz de Cholesky y Test LR

La extensión de la muestra hasta las 4,240 sesiones diarias (enero 2010 – noviembre 2026) revela dos fenómenos estructurales determinantes para la asignación estratégica de activos. En primer lugar, la correlación lineal media por pares se comprimió desde 0.4669 (2010–2020) hasta 0.2982 en el ciclo post-pandémico (2020–2026) (0.3976 en la muestra completa), ampliando sustancialmente el beneficio de diversificación entre crecimiento tecnológico/pagos digitales (MSFT, MA) y consumo defensivo/salud (WMT, JNJ, MCD). En segundo lugar, todas las series exhiben leptocurtosis individual extrema (curtosis diaria desde 11.59 en JNJ hasta 29.58 en MCD):

Lectura Económica de la Rotación de Markowitz y del Operador de Cholesky: Al resolver el problema cuadrático con restricción long-only sobre el historial completo 2010–2026, el Nuevo Portafolio Tangente Óptimo 2026 asigna WMT (28.10%) como ancla omnicanal antiinflacionaria, equilibra los dos motores de alto retorno compuesto MA (23.21%) y MSFT (22.86%), preserva a JNJ (19.79%) como estabilizador descorrelacionado (correlación post-2020 de apenas 0.078 frente a MSFT), reduce MCD a 6.04% y desplaza por completo a PEP (0.00%) debido a su compresión de márgenes reales en 2024–2026 (Sharpe individual de 0.402). Esta reasignación eleva el Ratio de Sharpe esperado desde 0.853 hasta 0.978 (retorno esperado de 18.26% anual con volatilidad contenida en 16.12%, inferior a la volatilidad de cualquiera de los seis activos por separado salvo la frontera mínima de GMV 2026 en 13.75%).

Por su parte, en la estimación conjunta de colas sobre los 4,239 vectores diarios en ℝ⁶, la log-verosimilitud asciende desde +77,607.79 bajo la Normal Multivariante hasta +80,929.86 bajo la t-Student Multivariante (ν̂_MLE = 3.937). El estadístico de Razón de Verosimilitud alcanza LR = 6,644.14 (p < 10^-100), rechazando de forma contundente la hipótesis gaussiana multivariante y demostrando que el cuarto momento conjunto diverge (2 < ν̂ < 4).

7. Playbook Operativo para el Inversor Institucional y Patrimonial (Mandato de $1,000,000 USD)

Una simulación estocástica institucional no tiene como fin contemplar curvas dispersas, sino transformar la distribución completa de probabilidades en reglas fiduciarias accionables de capitalización, presupuesto de riesgo y retiro sostenible. Tomando como motor rector el modelo más conservador en colas extremas —Salto-Difusión de Merton acoplado vía Cholesky (λ̂ = 1.96 saltos adversos/año con severidad media μ̂_J = −4.25%)— estructuramos cuatro directrices cuantitativas directas:

8. Torneo de los 5 Modelos Estocásticos a 10 Años (2026–2036) y Síntesis Ejecutiva

La Tabla 8.1 consolida el torneo comparativo completo ejecutado sobre un capital base estandarizado de W₀ = $10,000.00 USD a un horizonte de 10 años (H = 2,520 sesiones bursátiles / 120 meses, noviembre 2026 – noviembre 2036) con M = 10,000 trayectorias por especificación, contrastando tanto el Nuevo Portafolio Tangente Óptimo 2026 como la asignación histórica del Portafolio Defensivo 2020: