Monte Carlo Simulation of Temperature for Weather Derivative Pricing

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\}\)


 – 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('')
min_temp = pd.read_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.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()


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 =[0]
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):
    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[-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.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.legend(loc='lower right')

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
    - dataframe row with columns: T, Tbar, dTbar and vol
    - kappa: rate of mean reversion
    - temp: simulated next day temperatures
    if row['Tbar_shift'] != np.nan:
        T_i = row['Tbar']
        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
    - 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
    - mc_temps: DataFrame of all components individual components
    - mc_sims: DataFrame of all simulated temerpature paths
    if isinstance(trading_dates, pd.DatetimeIndex):
    # 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]], 
    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_Sims DataFrame

Simulated temperature paths for DatetimeIndex for \(\large M\) simulated paths.


Plot an Individual Simulation

mc_sims[1].plot(alpha=0.6,linewidth=2, marker='*')
plt.legend(loc='lower right')

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.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')