wonder why we would want to model a temperature time series as an Ornstein-Uhlenbeck process, and if this has any practical implications. Well, it does. While the O-U process is used in Physics to model the velocity of particles under friction, it is also used in finance to model interest rates, volatility, or spread mean reversion. This temperature application originates from an area where I have personal experience. My thesis, “Neural Network Approaches to Pricing Weather Derivatives,” focuses on developing new methods for pricing these derivatives.
The O-U process is commonly used in pricing weather derivatives. A financial derivative designed to deliver a payout based on an underlying climate index (temperature, rainfall, windspeed, snowfall). These derivatives enable anyone exposed to climate risks to manage their risk exposure. That could be farmers concerned about drought, energy companies worried about rising heating costs, or seasonal retailers concerned about a poor season. You can read much more about this in Weather Derivatives, Antonis Alexandridis K. (with whom I am currently working), and Achilleas D. Zapranis.
SDEs show up wherever uncertainty lies. The Black Scholes, Schrödinger Equations, and Geometric Brownian Motion are all Stochastic Differential Equations. These equations describe how systems evolve when a random component is involved. SDEs are particularly useful for forecasting under uncertainty, whether that be population growth, epidemic spread, or stock prices. Today, we will be looking at how we can model temperature using the Ornstein–Uhlenbeck process, which behaves quite a bit like the Heat Equation, a PDE describing temperature flow.
In my past two stories:
The Climate data is also available on my GitHub, so we can skip this step.
Today we will:
- Explain how the O-U process can model a temperature time-series.
- Talk about mean-reverting processes, stochastic differential equations, and how the O-U process relates to the heat equation.
- Use Python to fit an Ornstein-Uhlenbeck (O-U) process, a Stochastic Differential Equation (SDE), to NASA climate data.
In the figure above, we can see how volatility changes with the seasons. This emphasizes an important reason why the O-U process is well-suited for modelling this time series, as it does not assume constant volatility.
The O-U process, mean-reversion, and the heat equation
We will use the O-U process to model how temperature changes over time at a fixed point. In our equation, T(t) represents the temperature as a function of time (t). As you intuitively understand, the Earth’s temperatures cycle around a 365-day mean. This means that seasonal changes are modelled using s(t). We can think of these means as the blue and orange lines in the temperature graphs for Barrow, Alaska, and Niamey, Niger.
\[dT(t) = ds(t) – \kappa \left( T(t) – s(t) \right) \, dt + \sigma(t) \, dB(t)\]
Kappa (κ) is the speed of mean reversion. Much like its title would imply, this constant controls how fast our temperature reverts back to its seasonal means. Mean reversion in this context means that if today is an icy day, colder than the average, then we would expect tomorrow’s or next week’s temperature to be warmer as it fluctuates around that seasonal mean.
\[
dT(t) = ds(t) – {\color{red}{\kappa}} \big( T(t) – s(t) \big)\, dt + \sigma(t)\, dB(t)
\]
This constant plays a similar role in the heat equation, where it describes how fast a rod increases in temperature as it is being heated. In the heat equation, κ is the thermal diffusivity, which describes how fast heat flows through a material.
\[
\frac{\partial u(\mathbf{x},t)}{\partial t} = {\color{red}{\kappa}} \, \nabla^2 u(\mathbf{x},t) + q(\mathbf{x},t)
\]
- κ = 0: No mean reversion. The process reduces to Brownian motion with drift.
- 0 < κ < 1: Weak mean reversion. The shocks persist, slow decay of autocorrelation.
- κ > 1: Strong mean reversion. Deviations are corrected quickly, and the process is tightly clustered around μ.
Fourier Series
Estimating Mean and Volatility
\[s(t) = a + b t + \sum_{k=1}^{K_1} \left( \alpha_k \sin\left( \frac{2\pi k t}{365} \right) + \beta_k \cos\left( \frac{2\pi k t}{365} \right) \right)\]
\[\sigma^2(t) = c + \sum_{k=1}^{K_2} \left( \gamma_k \sin\left( \frac{2\pi k t}{365} \right) + \delta_k \cos\left( \frac{2\pi k t}{365} \right) \right)\]
Fitting the Ornstein–Uhlenbeck process to Mumbai data
We fit the mean of our temperature process
- Fit the mean using the Fourier Series
- Model short-term dependence — and calculate the mean reversion parameter κ using the AR(1) process
- Fit Volatility using Fourier Series
Fitting the mean
By fitting our data to 80% of the data, we can later use our SDE to forecast temperatures. Our mean is an OLS regression fitting s(t) to our data.
# --------------------
# 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"
It’s far harder to see the signal in the noise when looking at the residuals of our regression. We might be tempted to assume that this is white-noise, but we would be mistaken. There is still dependence in this data, a serial dependence. If we wish to model volatility later using a Fourier series, we must account for this dependence in the data.
# =================================================================
# 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="μ̂(t) (train)")
# 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="μ̂(t) (test)")
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)κ, ϕ, and the autoregressive process
Fit an AR(1) process. This tells you the speed of mean reversion.
An AR(1) process can capture the dependence between our residual at time step t based on the value at t-1. The exact value for κ depends on location. For Mumbai, κ = 0.861. This means temperatures return to the mean fairly quickly.
\[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)
In the graph above, we can see how our AR(1) process, in orange, estimates the next day’s temperature from 1981 to 2016. We can further justify the use of the AR(1) process, and gain further intuition into it by plotting Xₜ and Xₜ₋₁. Here we can see how these values are clearly correlated. By framing it this way, we can also see that the AR(1) process is simply linear regression. We do not even have to add a drift/intercept term as we have removed the mean. Hence, our intercept is always zero.
Now, looking at our AR(1) residuals, we can begin to think about how we can model the volatility of our temperature across time. From the graph below, we can see how our residuals pulsate across time in a seemingly periodic way. We can postulate that the points in the year with higher residuals correspond to time periods with higher volatility.
# --- 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)
Volatility Fourier Series
Our solution is to model the volatility using a Fourier series, much like what we did for our mean. To do this, we will have to make some transformations to our residual εₜ. We consider the squared εₜ² because volatility cannot be negative by definition.
\[\varepsilon_t = \sigma_t \eta_t \qquad \varepsilon_t^{2} = \sigma_t^{2}\eta_t^{2}\]
We can think of these residuals as being comprised of part standardized shocks ηₜ, with a mean of 0 and a variance of 1 (representing randomness), and volatility σₜ. We want to isolate for σₜ. This can be accomplished more easily by taking the log. But more importantly, doing this makes our error terms additive.
\[\displaystyle y_t := \log(\varepsilon_t^2) = \underbrace{\log\sigma_t^2}_{\text{deterministic, seasonal}} + \underbrace{\log\eta_t^2}_{\text{random noise}}\]
In summary, modeling log(εₜ²) helps to:
- keep σ^t₂ > 0
- Reduce the effect that outliers have on the fit, meaning they will have a less significant impact on the fit.
After we fit log(εₜ²), if we ever wish to recover εₜ², we exponentiate it later.
\[\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 σ̂(t) over time ---
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 σ̂(t)")
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(t) = ds(t) – \kappa (T(t) – s(t)),dt + \exp\left(\frac{1}{2}(\alpha + \beta^T z(t))\right), dB(t)\]
Demo of log stabilization
The graphs below illustrate how the logarithm aids in fitting volatility. The residual εₜ² contains significant outliers, partly due to the square operation performed to ensure positivity on the stabilized scale. Without the log stabilization, the large shocks dominate the fit, and the final result is biased.
# 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(t)")
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()
Conclusion
Now that we have fitted all these parameters to this SDE, we can consider using it to forecast. But getting here hasn’t been easy. Our solution heavily relied on two key concepts.
- Functions can be approximated through a Fourier series.
- The mean reversion parameter κ is equivalent to ϕ from our AR(1) process.
Whenever fitting a stochastic differential equation, whenever doing anything stochastic, I find it funny to think of how the word derives from stokhos, meaning aim. Even funnier is how that word turned into stokhazesthai, meaning aim/guess. That process must’ve involved some error and horrible aim.
References
Alexandridis, A. K., & Zapranis, A. D. (2013). Weather Derivatives: Modelling and Pricing Weather-Related Risk. Springer. https://doi.org/10.1007/978-1-4614-6071-8
Benth, F. E., & Šaltytė Benth, J. (2007). The volatility of temperature and pricing of weather derivatives. Quantitative Finance, 7(5), 553–561. https://doi.org/10.1080/14697680601155334
