In this online tutorial series dedicated to weather derivatives we have estimated the parameters of our modified mean-reverting Ornstein-Uhlenbeck process which defines our Temperature dynamics, and have now implemented different models for our time varying volatility. Now we move on to simulating temperature paths using Monte Carlo simulation method under the physical probability measure.
Once we have completed these simulation functions we can move onto implementing the Monte Carlo method under the risk-neutral pricing methodology for valuation of our temperature options.
Quick Summary on Weather Derivative Contract Features
Aim: we want to price temperature options.
Underlying: HDD/CDD index over given period.
The underlying over a temperature option is the heating/cooling degree days (HDD/CDD) index based on ‘approximation’ of average temperature and reference (/base) temperature.
1. Temperature Underlying
\(\large T_n = \frac{T^{max}+T^{min}}{2}\)
2. Degree Days
For a day \(n \in N\):
– \(\large HDD_n = (T_{ref}-T_n)^+\)
– \(\large CDD_n = (T_n – T_{ref})^+\)
3. Payoff Functions
Here the buyer of an option with receive an amount:
\(\large \xi = f(DD)\)
Payoff function \(f\) is computed on the cumulative index over a period \(P\):
– heating degree seasons \(DD = H_n = HDD^{N} = \sum^N_n HDD_n\)
– cooling degree seasons \(DD = C_n = CDD^{N} = \sum^N_n CDD_n\)
4. Typical Seasons OTC
– CDD season: 15-May to 15-Sep
– HDD season: 15-Dec(Nov) to 15-Mar
5. Popular Payoff Functions
Call with Cap
\(\large \xi = min\{\alpha(DD – K)^+, C\}\)
where:
– payoff rates \(\large \alpha\) is commonly US\(\$2,500\) or US\(\$5,000\)
– while caps \(\large C\) is commonly US\(\$500,000\) or US\(\$1,000,000\)
Quick Summary on Temperature Modelling
NOTE: Everything we have done so far in this series has been under the Physical probability measure \(\large \mathbb{P}\), including for our parameter estimation of the following models.
import os
import numpy as np
import pandas as pd
import datetime as dt
from scipy import interpolate
import matplotlib.pyplot as plt
max_temp = pd.read_csv('https://raw.githubusercontent.com/ASXPortfolio/jupyter-notebooks-data/main/maximum_temperature.csv')
min_temp = pd.read_csv('https://raw.githubusercontent.com/ASXPortfolio/jupyter-notebooks-data/main/minimum_temperature.csv')
def datetime(row):
return dt.datetime(row.Year,row.Month,row.Day)
max_temp['Date'] = max_temp.apply(datetime,axis=1)
min_temp['Date'] = min_temp.apply(datetime,axis=1)
max_temp.set_index('Date', inplace=True)
min_temp.set_index('Date', inplace=True)
drop_cols = [0,1,2,3,4,6,7]
max_temp.drop(max_temp.columns[drop_cols],axis=1,inplace=True)
min_temp.drop(min_temp.columns[drop_cols],axis=1,inplace=True)
max_temp.rename(columns={'Maximum temperature (Degree C)':'Tmax'}, inplace=True)
min_temp.rename(columns={'Minimum temperature (Degree C)':'Tmin'}, inplace=True)
temps = max_temp.merge(min_temp,how='inner',left_on=['Date'],right_on=['Date'])
def avg_temp(row):
return (row.Tmax+row.Tmin)/2
temps['T'] = temps.apply(avg_temp,axis=1)
# drop na values here
temps = temps.dropna()
temp_t = temps['T'].copy(deep=True)
temp_t = temp_t.to_frame()
temp_t.head()
Modified mean-reverting Ornstein-Uhlenbeck (OU) process
Our model for our mean-reverting Daily Average Temperature (DAT) is described in the following OU process SDE with estimated parameters.
\(\large dT_t = \left[\frac{d\bar{T_t}}{dt} + 0.438 (\bar{T_t} – T_t)\right]dt + \sigma_t dW_t\)
Where our changing average of DAT \(\large \bar{T_t}\) is:
\(\large \bar{T_t} = 16.8 + (3.32e-05)t + 5.05 sin((\frac{2\pi}{365.25})t + 1.27)\)
First derivative does not need finite difference approximation here because the function is differentiable \(\large \bar{T’_t}\)
\(\large \bar{T’_t} = (3.32e-05) + 5.05 (\frac{2\pi}{365.25}) cos((\frac{2\pi}{365.25})t + 1.27)\)
Where the date 01-Jan 1859 corresponds with the first ordinal number 0.
if isinstance(temp_t.index , pd.DatetimeIndex):
first_ord = temp_t.index.map(dt.datetime.toordinal)[0]
temp_t.index=temp_t.index.map(dt.datetime.toordinal)
def T_model(x, a, b, alpha, theta):
omega = 2*np.pi/365.25
T = a + b*x + alpha*np.sin(omega*x + theta)
return T
def dT_model(x, a, b, alpha, theta):
omega=2*np.pi/365.25
dT = b + alpha*omega*np.cos(omega*x + theta)
return dT
Tbar_params = [16.8, 3.32e-05, 5.05, 1.27]
temp_t['model_fit'] = T_model(temp_t.index-first_ord, *Tbar_params)
if not isinstance(temp_t.index , pd.DatetimeIndex):
temp_t.index=temp_t.index.map(dt.datetime.fromordinal)
temp_t[-2000:].plot(figsize=(12,4), style=['s','^-','k--'] , markersize=4, linewidth=2 )
Volatility of Temperature Process
The volatility estimator is based on the the quadratic variation \(\Large \sigma_t^2\) of of the temperature process \(\Large T_t\)
\(\Large \hat{\sigma}^2_t = \frac{1}{N_t} \sum^{N-1}_{i=0} (T_{i+1} – T_i)^2 \)
\(\Large \sigma_t$ is the dynamic volatility of the Temperature process. This could be both daily (as with our temperature dynamics) or seasonal (for example monthly)
temp_vol = temps['T'].copy(deep=True)
temp_vol = temp_vol.to_frame()
temp_vol['day'] = temp_vol.index.dayofyear
temp_vol['month'] = temp_vol.index.month
vol = temp_vol.groupby(['day'])['T'].agg(['mean','std'])
days = np.array(vol['std'].index)
T_std = np.array(vol['std'].values)
def spline(knots, x, y):
x_new = np.linspace(0, 1, knots+2)[1:-1]
t, c, k = interpolate.splrep(x, y, t=np.quantile(x, x_new), s=3)
yfit = interpolate.BSpline(t,c, k)(x)
return yfit
volatility = spline(5, days, T_std)
plt.figure(figsize=(8,4))
plt.plot(days, T_std, marker='*',label='Volatility')
plt.plot(days, volatility, linewidth=4, label='Spline Fit')
plt.ylabel('Std Dev (deg C)')
plt.xlabel('Day in Year')
plt.xlim(0,366)
plt.legend(loc='lower right')
plt.show()
Monte Carlo Simulaton of Temperature Paths
Now with our Euler scheme of approximation:
\(\large T_{i+1} = T_{i} + \bar{T’}_{i} + \kappa(\bar{T}_{i} – T_{i}) + \sigma_i z_i\)
def euler_step(row, kappa, M):
"""Function for euler scheme approximation step in
modified OH dynamics for temperature simulations
Inputs:
- dataframe row with columns: T, Tbar, dTbar and vol
- kappa: rate of mean reversion
Output:
- temp: simulated next day temperatures
"""
if row['Tbar_shift'] != np.nan:
T_i = row['Tbar']
else:
T_i = row['Tbar_shift']
T_det = T_i + row['dTbar']
T_mrev = kappa*(row['Tbar'] - T_i)
sigma = row['vol']*np.random.randn(M)
return T_det + T_mrev + sigma
def monte_carlo_temp(trading_dates, Tbar_params, vol_model, first_ord, M=1, kappa=0.438):
"""Monte Carlo simulation of temperature
Inputs:
- trading_dates: pandas DatetimeIndex from start to end dates
- M: number of simulations
- Tbar_params: parameters used for Tbar model
- vol_model: fitted volatility model with days in year index
- first_ord: first ordinal of fitted Tbar model
Outputs:
- mc_temps: DataFrame of all components individual components
- mc_sims: DataFrame of all simulated temerpature paths
"""
if isinstance(trading_dates, pd.DatetimeIndex):
trading_date=trading_dates.map(dt.datetime.toordinal)
# Use Modified Ornstein-Uhlenbeck process with estimated parameters to simulate Tbar DAT
Tbars = T_model(trading_date-first_ord, *Tbar_params)
# Use derivative of modified OH process SDE to calculate change of Tbar
dTbars = dT_model(trading_date-first_ord, *Tbar_params)
# Create DateFrame with thi
mc_temps = pd.DataFrame(data=np.array([Tbars, dTbars]).T,
index=trading_dates, columns=['Tbar','dTbar'])
# Create columns for day in year
mc_temps['day'] = mc_temps.index.dayofyear
# Apply BSpline volatility model depending on day of year
mc_temps['vol'] = vol_model[mc_temps['day']-1]
# Shift Tbar by one day (lagged Tbar series)
mc_temps['Tbar_shift'] = mc_temps['Tbar'].shift(1)
# Apply Euler Step Pandas Function
data = mc_temps.apply(euler_step, args=[kappa, M], axis=1)
# Create final DataFrame of all simulations
mc_sims = pd.DataFrame(data=[x for x in [y for y in data.values]],
index=trading_dates,columns=range(1,M+1))
return mc_temps, mc_sims
Simulation of Temperature
Here we finally simulate temperature under the Physical probability measure \(\large \mathbb{P}\)
# define trading date range no_sims = 5 trading_dates = pd.date_range(start='2022-09-01', end='2025-08-31', freq='D') mc_temps, mc_sims = monte_carlo_temp(trading_dates, Tbar_params, volatility, first_ord, no_sims)
MC_Temps DataFrame
Individual components used to simulate temperature path simulations
mc_temps
MC_Sims DataFrame
Simulated temperature paths for DatetimeIndex for \(\large M\) simulated paths.
mc_sims
Plot an Individual Simulation
plt.figure(figsize=(12,6)) mc_sims[1].plot(alpha=0.6,linewidth=2, marker='*') mc_temps["Tbar"].plot(linewidth=4) plt.legend(loc='lower right') plt.show()
Temperature distributions
Let’s finally observe the difference between our monte carlo temperature simulations distribution between the peaks of summer and winter days.
Here we simulate 10,000 temperatures on 1st July, 2023 for Winter and 1st Jan, 2023 for a representation of peak Summer temps and volatility.
no_sims = 10000
trading_dates_winter = pd.date_range(start='2023-07-01', end='2023-07-01', freq='D')
mc_temps_winter, mc_sims_winter = monte_carlo_temp(trading_dates_winter, Tbar_params, volatility, first_ord, no_sims)
trading_dates_summer = pd.date_range(start='2023-01-01', end='2023-01-01', freq='D')
mc_temps_summer, mc_sims_summer = monte_carlo_temp(trading_dates_summer, Tbar_params, volatility, first_ord, no_sims)
plt.figure(figsize=(12,6))
plt.title('Winter vs Summer Temperature MC Sims')
Tbar_summer = mc_temps_summer.iloc[-1,:]['Tbar']
Tbar_winter = mc_temps_winter.iloc[-1,:]['Tbar']
plt.hist(mc_sims_summer.iloc[-1,:],bins=20, alpha=0.5, label='Summer', color='tab:orange')
plt.plot([Tbar_summer,Tbar_summer],[0,1600], linewidth=4, label='Tbar_summer', color='tab:orange')
plt.title('Summer vs Winter Temperature MC Sims')
plt.hist(mc_sims_winter.iloc[-1,:],bins=20, alpha=0.5, label='Winter', color='tab:blue')
plt.plot([Tbar_winter,Tbar_winter],[0,1600], linewidth=4, label='Tbar_winter', color='tab:blue')
plt.legend()
plt.show()
