Me pregunto por qué nos gustaría modelar una serie temporal de temperatura como un proceso Ornstein-Uhlenbeck, y si esto tiene alguna implicación práctica. Bueno, lo hace. Si bien el proceso OU se usa en física para modelar la velocidad de las partículas bajo fricción, también se usa en finanzas para modelar las tasas de interés, la volatilidad o la reversión media de propagación. Esta aplicación de temperatura se origina en un área donde tengo experiencia personal. Mi tesis, “Enfoques de red neuronal para fijar el precio de la derivada del climaS ”, se centra en desarrollar Nuevos métodos para fijar el precio de estos derivados.
El proceso OU se usa comúnmente en el precio de los derivados del clima. Un derivado financiero diseñado para entregar un pago basado en un índice climático subyacente (temperatura, lluvia, velocidad de viento, nevadas). Estos derivados permiten a cualquier persona expuesta a riesgos climáticos para manejar su exposición al riesgo. Eso podría ser agricultores preocupados por la sequía, las compañías energéticas preocupadas por el aumento de los costos de calefacción o los minoristas estacionales preocupados por una temporada pobre. Puedes leer mucho más sobre esto en Derivados del climaAntonis Alexandridis K. (con quien estoy trabajando actualmente), y Achilleas D. Zapranis.
Los SDE aparecen donde sea la incertidumbre. El Scholes negros, Ecuaciones de Schrödingery Movimiento geométrico brownian son todas las ecuaciones diferenciales estocásticas. Estas ecuaciones describen cómo evolucionan los sistemas cuando está involucrado un componente aleatorio. Los SDE son particularmente útiles para pronosticar bajo incertidumbre, ya sea crecimiento de la población, propagación epidémicao precios de las acciones. Hoy veremos cómo podemos modelar la temperatura utilizando el Ornstein – Uhlenbeck proceso, que se comporta bastante como el Ecuación de calor, un PDE que describe la temperatura FLAy.
En mis dos historias anteriores:
Los datos del clima también están disponibles en mi GitHub, por lo que podemos omitir este paso.
Hoy lo haremos:
- Explique cómo el proceso OU puede modelar una serie de tiempo de temperatura.
- Hable sobre los procesos de reversión media, las ecuaciones diferenciales estocásticas y cómo el proceso OU se relaciona con la ecuación de calor.
- Use Python para adaptarse a un proceso Ornstein-Uhlenbeck (OU), una ecuación diferencial estocástica (SDE), a los datos climáticos de la NASA.
En la figura anterior, podemos ver cómo cambia la volatilidad con las estaciones. Esto enfatiza una razón importante por la cual el proceso OU es adecuado para modelar esta serie de tiempo, ya que no asume una volatilidad constante.
El proceso OU, la reversión media y la ecuación de calor
Utilizaremos el proceso OU para modelar cómo cambia la temperatura con el tiempo en un punto fijo. En nuestra ecuación, t
\[dT
Kappa (κ) es la velocidad de reversión media. Al igual que su título implicaría, esta constante controla la rapidez con que nuestra temperatura vuelve a sus medios estacionales. La reversión media en este contexto significa que si hoy es un día helado, más frío que el promedio, entonces esperaríamos que la temperatura de los mañana o la próxima semana sea más cálida a medida que fluctúa sobre esa media estacional.
\[
dT
Esta constante juega un papel similar en la ecuación de calor, donde describe qué tan rápido aumenta una caña de temperatura a medida que se calienta. En la ecuación de calor, κ es la difusividad térmica, que describe cómo fluye el calor rápido a través de un material.
\[
\frac{\partial u(\mathbf{x},t)}{\partial t} = {\color{red}{\kappa}} \, \nabla^2 u(\mathbf{x},t) + q(\mathbf{x},t)
\]
- κ = 0: Sin reversión media. El proceso se reduce al movimiento browniano con la deriva.
- 0 <κ <1: Reversión media débil. Los choques persisten, la lenta descomposición de la autocorrelación.
- κ> 1: Fuerte reversión media. Las desviaciones se corrigen rápidamente, y el proceso se agrupa estrechamente alrededor de μ.
Serie de Fourier
Estimación de media y volatilidad
\[s
\[\sigma^2
Ajuste del proceso Ornstein -Uhlenbeck a los datos de Mumbai
Encajamos en la media de nuestro proceso de temperatura
- Ajuste la media usando la serie Fourier
- Modelo de dependencia a corto plazo, y calcule el parámetro de reversión media κ usando el proceso AR (1)
- Ajuste la volatilidad utilizando la serie Fourier
Ajustando la media
Al ajustar nuestros datos al 80% de los datos, luego podemos usar nuestro SDE para pronosticar temperaturas. Nuestra media es un ajuste de regresión de OLS S (T) a nuestros datos.
# --------------------
# Config (parameters and paths for analysis)
# --------------------
CITY = "Mumbai, India" # Name of the city being analyzed (used in labels/plots)
SEED = 42 # Random seed for reproducibility of results
SPLIT_FRAC = 0.80 # Fraction of data to use for training (rest for testing)
MEAN_HARM = 2 # Number of harmonic terms to use for modeling the mean (seasonality)
VOL_HARM = 3 # Number of harmonic terms to use for modeling volatility (seasonality)
LJUNG_LAGS = 10 # Number of lags for Ljung-Box test (check autocorrelation in residuals)
EPS = 1e-12 # Small value to avoid division by zero or log(0) issues
MIN_TEST_N = 8 # Minimum number of test points required to keep a valid test set
# --------------------
# Paths (where input/output files are stored)
# --------------------
# Base directory in Google Drive where climate data and results are stored
BASE_DIR = Path("/content/drive/MyDrive/TDS/Climate")
# Subdirectory specifically for outputs of the Mumbai analysis
OUT_BASE = BASE_DIR / "Benth_Mumbai"
# Subfolder for saving plots generated during analysis
PLOTS_DIR = OUT_BASE / "plots"
# Ensure the output directories exist (create them if they don’t)
OUT_BASE.mkdir(parents=True, exist_ok=True)
PLOTS_DIR.mkdir(parents=True, exist_ok=True)
# Path to the climate data CSV file (input dataset)
CSV_PATH = BASE_DIR / "climate_data.csv"
Es mucho más difícil ver la señal en el ruido al mirar los residuos de nuestra regresión. Podríamos estar tentados a asumir que este es un ruido blanco, pero estaríamos equivocados. Todavía hay dependencia en estos datos, una dependencia en serie. Si deseamos modelar la volatilidad más tarde utilizando una serie de Fourier, debemos tener en cuenta esta dependencia en los datos.
# =================================================================
# 1) MEAN MODEL: residuals after fitting seasonal mean + diagnostics
# =================================================================
# Residuals after subtracting fitted seasonal mean (before AR step)
mu_train = mean_fit.predict(train)
x_train = train["DAT"] - mu_train
# Save mean model regression summary to text file
save_model_summary_txt(mean_fit, DIAG_DIR / f"{m_slug}_mean_OLS_summary.txt")
# --- Plot: observed DAT vs fitted seasonal mean (train + test) ---
fig = plt.figure(figsize=(12,5))
plt.plot(train["Date"], train["DAT"], lw=1, alpha=0.8, label="DAT (train)")
plt.plot(train["Date"], mu_train, lw=2, label="μ̂
# Predict and plot fitted mean for test set
mu_test = mean_fit.predict(test[["trend"] + [c for c in train.columns if c.startswith(("cos","sin"))][:2*MEAN_HARM]])
plt.plot(test["Date"], test["DAT"], lw=1, alpha=0.8, label="DAT (test)")
plt.plot(test["Date"], mu_test, lw=2, label="μ̂
plt.title("Mumbai — DAT and seasonal mean fit")
plt.xlabel("Date"); plt.ylabel("Temperature (DAT)")
plt.legend()
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_mean_fit_timeseries.png", dpi=160); plt.close(fig)
# --- Plot: mean residuals (time series) ---
fig = plt.figure(figsize=(12,4))
plt.plot(train["Date"], x_train, lw=1)
plt.axhline(0, color="k", lw=1)
plt.title("Mumbai — Residuals after mean fit (x_t = DAT - μ̂)")
plt.xlabel("Date"); plt.ylabel("x_t")
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_mean_residuals_timeseries.png", dpi=160); plt.close(fig)κ, ϕ y el proceso autorregresivo
Ajuste un proceso AR (1). Esto te dice la velocidad de la reversión media.
Un proceso AR (1) puede capturar la dependencia entre nuestro residual en el paso T de tiempo T en función del valor en T-1. El valor exacto para κ depende de la ubicación. Para Mumbai, κ = 0.861. Esto significa que las temperaturas vuelven a la media con bastante rapidez.
\[X_t = c + \phi X_{t-1} + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0,\sigma^2)\]
# =======================================================
# 2) AR(1) MODEL: fit and residual diagnostics
# =======================================================
# Extract AR(1) parameter φ and save summary
phi = float(ar_fit.params.iloc[0]) if len(ar_fit.params) else np.nan
save_model_summary_txt(ar_fit, DIAG_DIR / f"{m_slug}_ar1_summary.txt")
# --- Scatterplot: x_t vs x_{t-1} with fitted line φ x_{t-1} ---
x_t = x_train.iloc[1:].values
x_tm1 = x_train.iloc[:-1].values
x_pred = phi * x_tm1
fig = plt.figure(figsize=(5.8,5.2))
plt.scatter(x_tm1, x_t, s=10, alpha=0.6, label="Observed")
xline = np.linspace(np.min(x_tm1), np.max(x_tm1), 2)
plt.plot(xline, phi*xline, lw=2, label=f"Fitted: x_t = {phi:.3f} x_(t-1)")
plt.title("Mumbai — AR(1) scatter: x_t vs x_{t-1}")
plt.xlabel("x_{t-1}"); plt.ylabel("x_t"); plt.legend()
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_ar1_scatter.png", dpi=160); plt.close(fig)
# --- Time series: actual vs fitted AR(1) values ---
fig = plt.figure(figsize=(12,4))
plt.plot(train["Date"].iloc[1:], x_t, lw=1, label="x_t")
plt.plot(train["Date"].iloc[1:], x_pred, lw=2, label=r"$\hat{x}_t = \phi x_{t-1}$")
plt.axhline(0, color="k", lw=1)
plt.title("Mumbai — AR(1) fitted vs observed deviations")
plt.xlabel("Date"); plt.ylabel("x_t")
plt.legend()
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_ar1_timeseries_fit.png", dpi=160); plt.close(fig)
# --- Save AR diagnostics (φ + Ljung-Box test p-value) ---
from statsmodels.stats.diagnostic import acorr_ljungbox
try:
lb = acorr_ljungbox(e_train, lags=[10], return_df=True)
lb_p = float(lb["lb_pvalue"].iloc[0]) # test for autocorrelation
lb_stat = float(lb["lb_stat"].iloc[0])
except Exception:
lb_p = np.nan; lb_stat = np.nan
pd.DataFrame([{"phi": phi, "ljungbox_stat_lag10": lb_stat, "ljungbox_p_lag10": lb_p}])\
.to_csv(DIAG_DIR / f"{m_slug}_ar1_diagnostics.csv", index=False)En el gráfico anterior, podemos ver cómo nuestro proceso AR (1), en naranja, estima la temperatura del día siguiente de 1981 a 2016. Podemos justificar aún más el uso del proceso AR (1) y obtener una intuición adicional en él trazando Xₜ y Xₜ₋₁. Aquí podemos ver cómo estos valores están claramente correlacionados. Al enmarcarlo de esta manera, también podemos ver que el proceso AR (1) es simplemente una regresión lineal. Ni siquiera tenemos que agregar un término de deriva/intercepción, ya que hemos eliminado la media. Por lo tanto, nuestra intersección siempre es cero.
Ahora, mirando nuestros residuos AR (1), podemos comenzar a pensar en cómo podemos modelar la volatilidad de nuestra temperatura a lo largo del tiempo. Desde el gráfico a continuación, podemos ver cómo nuestros residuos pulsan a través del tiempo de una manera aparentemente periódica. Podemos postular que los puntos en el año con residuos más altos corresponden a períodos de tiempo con mayor volatilidad.
# --- Residuals from AR(1) model ---
e_train = ar_fit.resid.astype(float).values
fig = plt.figure(figsize=(12,4))
plt.plot(train["Date"].iloc[1:], e_train, lw=1)
plt.axhline(0, color="k", lw=1)
plt.title("Mumbai — AR(1) residuals ε_t")
plt.xlabel("Date"); plt.ylabel("ε_t")
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_ar1_residuals_timeseries.png", dpi=160); plt.close(fig)
Serie de volatilidad Fourier
Nuestra solución es modelar la volatilidad utilizando una serie de Fourier, muy parecido a lo que hicimos por nuestra media. Para hacer esto, tendremos que hacer algunas transformaciones en nuestro εₜ residual. Consideramos el εₜ² cuadrado porque la volatilidad no puede ser negativa por definición.
\[\varepsilon_t = \sigma_t \eta_t \qquad \varepsilon_t^{2} = \sigma_t^{2}\eta_t^{2}\]
Podemos pensar que estos residuos se componen de choques estandarizados de pieza ηₜ, con una media de 0 y una varianza de 1 (que representa aleatoriedad) y volatilidad σₜ. Queremos aislar para σₜ. Esto se puede lograr más fácilmente tomando el registro. Pero lo que es más importante, hacer esto hace que nuestros términos de error sea aditivo.
\[\displaystyle y_t := \log(\varepsilon_t^2) = \underbrace{\log\sigma_t^2}_{\text{deterministic, seasonal}} + \underbrace{\log\eta_t^2}_{\text{random noise}}\]
En resumen, el registro de modelado (εₜ²) ayuda a:
- Mantenga σ^t₂> 0
- Reduzca el efecto que los valores atípicos tienen en el ajuste, lo que significa que tendrán un impacto menos significativo en el ajuste.
Después de ajustar el registro (εₜ²), si alguna vez deseamos recuperar εₜ², lo exponemos más tarde.
\[\displaystyle \widehat{\sigma}_t^{\,2} \;=\; \exp\!\big(\widehat y_t\big)\]
# =======================================================
# 3) VOLATILITY REGRESSION: fit log(ε_t^2) on seasonal harmonics
# =======================================================
if vol_fit is not None:
# Compute log-squared residuals (proxy for variance)
log_eps2 = np.log(e_train**2 + 1e-12)
# Use cosine/sine harmonics as regressors for volatility
feats = train.iloc[1:][[c for c in train.columns if c.startswith(("cos","sin"))][:2*VOL_HARM]]
vol_terms = [f"{b}{k}" for k in range(1, VOL_HARM + 1) for b in ("cos","sin")]
Xbeta = np.asarray(vol_fit.predict(train.iloc[1:][vol_terms]), dtype=float)
# --- Time plot: observed vs fitted log-variance ---
fig = plt.figure(figsize=(12,4))
plt.plot(train["Date"].iloc[1:], log_eps2, lw=1, label="log(ε_t^2)")
plt.plot(train["Date"].iloc[1:], Xbeta, lw=2, label="Fitted log-variance")
plt.title("Mumbai — Volatility regression (log ε_t^2)")
plt.xlabel("Date"); plt.ylabel("log variance")
plt.legend()
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_volatility_logvar_timeseries.png", dpi=160); plt.close(fig)
# --- Plot: estimated volatility σ̂
fig = plt.figure(figsize=(12,4))
plt.plot(train["Date"].iloc[1:], sigma_tr, lw=1.5, label="σ̂ (train)")
if 'sigma_te' in globals():
plt.plot(test["Date"], sigma_te, lw=1.5, label="σ̂ (test)")
plt.title("Mumbai — Conditional volatility σ̂
plt.xlabel("Date"); plt.ylabel("σ̂")
plt.legend()
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_volatility_sigma_timeseries.png", dpi=160); plt.close(fig)
# --- Coefficients of volatility regression (with CI) ---
vol_coef = coef_ci_frame(vol_fit).sort_values("estimate")
fig = plt.figure(figsize=(8, max(4, 0.4*len(vol_coef))))
y = np.arange(len(vol_coef))
plt.errorbar(vol_coef["estimate"], y, xerr=1.96*vol_coef["stderr"], fmt="o", capsize=3)
plt.yticks(y, vol_coef["term"])
plt.axvline(0, color="k", lw=1)
plt.title("Mumbai — Volatility model coefficients (95% CI)")
plt.xlabel("Estimate")
fig.tight_layout(); fig.savefig(DIAG_DIR / f"{m_slug}_volatility_coefficients.png", dpi=160); plt.close(fig)
# Save volatility regression summary
save_model_summary_txt(vol_fit, DIAG_DIR / f"{m_slug}_volatility_summary.txt")
else:
print("Volatility model not available (too few points or regression failed). Skipping vol plots.")
print("Diagnostics saved to:", DIAG_DIR)\[dT
Demostración de estabilización de log
Los gráficos a continuación ilustran cómo el logaritmo ayuda a ajustar la volatilidad. El εₜ² residual contiene valores atípicos significativos, en parte debido a la operación cuadrada realizada para garantizar la positividad en la escala estabilizada. Sin la estabilización del registro, los grandes choques dominan el ajuste y el resultado final está sesgado.
# Fix the seasonal plot from the previous cell (indexing bug) and
# add a *second scenario* with rare large outliers to illustrate instability
# when regressing on skewed eps^2 directly.
# ------------------------------
# 1) Seasonal view (fixed indexing)
# ------------------------------
# Recreate the arrays from the prior simulation cell by re-executing the same seed & setup
np.random.seed(7) # deterministic reproducibility
n_days = 730 # two synthetic years (365 * 2)
t = np.arange(n_days) # day index 0..729
omega = 2 * np.pi / 365.0 # seasonal frequency (one-year period)
# True (data-generating) log-variance is seasonal via sin/cos
a_true, b_true, c_true = 0.2, 0.9, -0.35
log_var_true = a_true + b_true * np.sin(omega * t) + c_true * np.cos(omega * t)
# Convert log-variance to standard deviation: sigma = exp(0.5 * log var)
sigma_true = np.exp(0.5 * log_var_true)
# White noise innovations; variance scaled by sigma_true
z = np.random.normal(size=n_days)
eps = sigma_true * z # heteroskedastic residuals
eps2 = eps**2 # "variance-like" target if regressing raw eps^2
# Design matrix with intercept + annual sin/cos harmonics
X = np.column_stack([np.ones(n_days), np.sin(omega * t), np.cos(omega * t)])
# --- Fit A: OLS on raw eps^2 (sensitive to skew/outliers) ---
beta_A, *_ = np.linalg.lstsq(X, eps2, rcond=None)
var_hat_A = X @ beta_A # fitted variance (can be negative from OLS)
var_hat_A = np.clip(var_hat_A, 1e-8, None) # clip to avoid negatives
sigma_hat_A = np.sqrt(var_hat_A) # convert to sigma
# --- Fit B: OLS on log(eps^2) (stabilizes scale & reduces skew) ---
eps_safe = 1e-12 # small epsilon to avoid log(0)
y_log = np.log(eps2 + eps_safe) # stabilized target
beta_B, *_ = np.linalg.lstsq(X, y_log, rcond=None)
log_var_hat_B = X @ beta_B
sigma_hat_B = np.exp(0.5 * log_var_hat_B)
# Day-of-year index for seasonal averaging across years
doy = t % 365
# Correct ordering for the first year's DOY values
order365 = np.argsort(doy[:365])
def seasonal_mean(x):
"""
Average the two years day-by-day to get a single seasonal curve.
Assumes x has length 730 (two years); returns length-365 array.
"""
return 0.5 * (x[:365] + x[365:730])
# Plot one "synthetic year" view of the seasonal sigma pattern
plt.figure(figsize=(12, 5))
plt.plot(np.sort(doy[:365]), seasonal_mean(sigma_true)[order365], label="True sigma seasonality")
plt.plot(np.sort(doy[:365]), seasonal_mean(sigma_hat_A)[order365], label="Fitted from eps^2 regression", alpha=0.9)
plt.plot(np.sort(doy[:365]), seasonal_mean(sigma_hat_B)[order365], label="Fitted from log(eps^2) regression", alpha=0.9)
plt.title("Seasonal Volatility Pattern: True vs Fitted (one-year view) – Fixed")
plt.xlabel("Day of year")
plt.ylabel("Sigma")
plt.legend()
plt.tight_layout()
plt.show()
# ------------------------------
# 2) OUTLIER SCENARIO to illustrate instability
# ------------------------------
np.random.seed(21) # separate seed for the outlier experiment
def run_scenario(n_days=730, outlier_rate=0.05, outlier_scale=8.0):
"""
Generate two-year heteroskedastic residuals with occasional huge shocks
to mimic heavy tails. Compare:
- Fit A: regress raw eps^2 on sin/cos (can be unstable, negative fits)
- Fit B: regress log(eps^2) on sin/cos (more stable under heavy tails)
Return fitted sigmas and error metrics (MAE/MAPE), plus diagnostics.
"""
t = np.arange(n_days)
omega = 2 * np.pi / 365.0
# Same true seasonal log-variance as above
a_true, b_true, c_true = 0.2, 0.9, -0.35
log_var_true = a_true + b_true * np.sin(omega * t) + c_true * np.cos(omega * t)
sigma_true = np.exp(0.5 * log_var_true)
# Base normal innovations
z = np.random.normal(size=n_days)
# Inject rare, huge shocks to create heavy tails in eps^2
mask = np.random.rand(n_days) < outlier_rate
z[mask] *= outlier_scale
# Heteroskedastic residuals and their squares
eps = sigma_true * z
eps2 = eps**2
# Same sin/cos design
X = np.column_stack([np.ones(n_days), np.sin(omega * t), np.cos(omega * t)])
# --- Fit A: raw eps^2 on X (OLS) ---
beta_A, *_ = np.linalg.lstsq(X, eps2, rcond=None)
var_hat_A_raw = X @ beta_A
neg_frac = np.mean(var_hat_A_raw < 0.0) # fraction of negative variance predictions
var_hat_A = np.clip(var_hat_A_raw, 1e-8, None) # clip to ensure non-negative variance
sigma_hat_A = np.sqrt(var_hat_A)
# --- Fit B: log(eps^2) on X (OLS on log-scale) ---
y_log = np.log(eps2 + 1e-12)
beta_B, *_ = np.linalg.lstsq(X, y_log, rcond=None)
log_var_hat_B = X @ beta_B
sigma_hat_B = np.exp(0.5 * log_var_hat_B)
# Error metrics comparing fitted sigmas to the true sigma path
mae = lambda a, b: np.mean(np.abs(a - b))
mape = lambda a, b: np.mean(np.abs((a - b) / (a + 1e-12))) * 100
mae_A = mae(sigma_true, sigma_hat_A)
mae_B = mae(sigma_true, sigma_hat_B)
mape_A = mape(sigma_true, sigma_hat_A)
mape_B = mape(sigma_true, sigma_hat_B)
return {
"t": t,
"sigma_true": sigma_true,
"sigma_hat_A": sigma_hat_A,
"sigma_hat_B": sigma_hat_B,
"eps2": eps2,
"y_log": y_log,
"neg_frac": neg_frac,
"mae_A": mae_A, "mae_B": mae_B,
"mape_A": mape_A, "mape_B": mape_B
}
# Run with 5% outliers scaled 10x to make the point obvious
res = run_scenario(outlier_rate=0.05, outlier_scale=10.0)
print("\nOUTLIER SCENARIO (5% of days have 10x shocks) — illustrating instability when using eps^2 directly")
print(f" MAE (sigma): raw eps^2 regression = {res['mae_A']:.4f} | log(eps^2) regression = {res['mae_B']:.4f}")
print(f" MAPE (sigma): raw eps^2 regression = {res['mape_A']:.2f}% | log(eps^2) regression = {res['mape_B']:.2f}%")
print(f" Negative variance predictions before clipping (raw fit): {res['neg_frac']:.2%}")
# Visual comparison: true sigma vs two fitted approaches under outliers
plt.figure(figsize=(12, 5))
plt.plot(res["t"], res["sigma_true"], label="True sigma
plt.plot(res["t"], res["sigma_hat_A"], label="Fitted sigma from eps^2 regression", alpha=0.9)
plt.plot(res["t"], res["sigma_hat_B"], label="Fitted sigma from log(eps^2) regression", alpha=0.9)
plt.title("True vs Fitted Volatility with Rare Large Shocks")
plt.xlabel("Day")
plt.ylabel("Sigma")
plt.legend()
plt.tight_layout()
plt.show()
# Show how the targets behave when outliers are present
plt.figure(figsize=(12, 5))
plt.plot(res["t"], res["eps2"], label="eps^2 (now extremely heavy-tailed due to outliers)")
plt.title("eps^2 under outliers: unstable target for regression")
plt.xlabel("Day")
plt.ylabel("eps^2")
plt.legend()
plt.tight_layout()
plt.show()
plt.figure(figsize=(12, 5))
plt.plot(res["t"], res["y_log"], label="log(eps^2): compressed & stabilized")
plt.title("log(eps^2) under outliers: stabilized scale")
plt.xlabel("Day")
plt.ylabel("log(eps^2)")
plt.legend()
plt.tight_layout()
plt.show()
Conclusión
Ahora que hemos ajustado todos estos parámetros a este SDE, podemos considerar usarlo para pronosticar. Pero llegar aquí no ha sido fácil. Nuestra solución dependía en gran medida de dos conceptos clave.
- Las funciones se pueden aproximar a través de una serie de Fourier.
- El parámetro de reversión media κ es equivalente a ϕ de nuestro proceso AR (1).
Siempre que se ajuste a una ecuación diferencial estocástica, siempre que haga algo estocásticoMe resulta divertido pensar en cómo se deriva la palabra Stokhos, significado objetivo. Incluso más divertido es cómo se convirtió esa palabra en Stokhazesthai, significado objetivo/adivinar. Ese proceso debe haber involucrado algún error y un objetivo horrible.
Referencias
Alexandridis, AK y Zapranis, AD (2013). Derivados del clima: modelado y fijación de precios de riesgo relacionado con el clima. Saltador. https://doi.org/10.1007/978-1-4614-6071-8
Benth, Fe y Šaltytė Benth, J. (2007). La volatilidad de la temperatura y el precio de los derivados del clima. Finanzas cuantitativas, 7 (5), 553–561. https://doi.org/10.1080/14697680601155334