Risk-Neutral Pricing of Weather Derivatives

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^{max}}{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. 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 numpy as np
import pandas as pd
import datetime as dt
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy import stats, interpolate

max_temp = pd.read_csv('max_temp.csv')
min_temp = pd.read_csv('min_temp.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()
first_ord = temp_t.index.map(dt.datetime.toordinal)[0]

temp_vol = temps['T'].copy(deep=True).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)

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.

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)

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

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

Tbar_params = [16.8, 3.32e-05, 5.05, 1.27]

Simulating 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, lamda):
    """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['T'] != np.nan:
        T_i = row['Tbar']
    else:
        T_i = row['T']  
    T_det = T_i + row['dTbar']
    T_mrev =  kappa*(row['Tbar'] - T_i)
    sigma = row['vol']*np.random.randn(M)
    riskn = lamda*row['vol']
    return T_det + T_mrev + sigma - riskn

def monte_carlo_temp(trading_dates, Tbar_params, vol_model, first_ord, M=1, lamda=0):
    """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 and simulated temperatures
    """
    kappa=0.438
    if isinstance(trading_dates, pd.DatetimeIndex):
        trading_date=trading_dates.map(dt.datetime.toordinal)
  
    Tbars = T_model(trading_date-first_ord, *Tbar_params) 
    dTbars = dT_model(trading_date-first_ord, *Tbar_params) 
    mc_temps = pd.DataFrame(data=np.array([Tbars, dTbars]).T,
                            index=trading_dates, columns=['Tbar','dTbar'])
    mc_temps['day'] = mc_temps.index.dayofyear
    mc_temps['vol'] = vol_model[mc_temps['day']-1]
    
    mc_temps['T'] = mc_temps['Tbar'].shift(1)  
    data = mc_temps.apply(euler_step, args=[kappa, M, lamda], axis=1)
    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

Risk-Neutral Pricing of Weather Derivatives

A summary of everything we have done in this series in only three points:

  1. Aim: Trying to calculate fair value of a temperature option
  2. Temperature SDE: Modified mean-reverting Ornstein-Uhlenbeck (OU) process
  3. Volatility: B-Spline interpolation as our non-constant volatility model

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

Temperature Simulation

Here is an example of a single sample path of temperature simulated under the Physical probability measure \(\large \mathbb{P}\)

# define trading date range
trading_dates = pd.date_range(start='2022-09-01', end='2023-08-31', freq='D')
volatility = spline(5, days, T_std)
mc_temps, mc_sims = monte_carlo_temp(trading_dates, Tbar_params, volatility, first_ord)

plt.figure(figsize=(10,6))
mc_sims[1].plot(alpha=0.5,linewidth=1, marker='*')
mc_temps["Tbar"].plot(linewidth=3)
plt.legend(loc='lower right')
plt.show()

Risk-Neutral Pricing Methodology

The market for weather derivatives is a typical example of an incomplete market,
because the underlying variable, the temperature, is not tradable. Therefore we
have to consider the market price of risk \(\large \lambda\), in order to obtain unique prices for such contracts.

Since there is not yet a real market from which we can
obtain prices, we assume for simplicity that the market price of risk is constant.
Furthermore, we assume that we are given a risk free asset with constant interest rate \(\large r\) and a contract that for each degree Celsius pays one unit of currency.
Thus, under a martingale measure \(\large \mathbb{Q}\), characterized by the market price of risk
\(\large \lambda\), our price process also denoted by \(\large T_t\) satisfies the following dynamics:

\(\large dT_t^\mathbb{Q} = \left[\frac{d\bar{T_t}}{dt} + \kappa (\bar{T_t} – T_t) – \lambda \sigma_t \right]dt + \sigma_t dW_t^\mathbb{Q}\)

Compared to the following dynamics under the physical probability measure \(\large \mathbb{P}\):

\(\large dT_t^\mathbb{P} = \left[\frac{d\bar{T_t}}{dt} + \kappa (\bar{T_t} – T_t)\right]dt + \sigma_t dW_t^\mathbb{P}\)

We showed in a previous video we showed by applying Ito-Doeblin to our SDE where our ito process is \(T_t\) with dynamics \(dT_t\) as described by our mean-reverting SDE, we can arrive at the following expectation and variance.

\(\large T_t = \bar{T_t} + (T_s – \bar{T_s}) e^{-\kappa (t-s)} + \int_s^t e^{-\kappa (t-u)} \sigma_u dW_u\)

\(\large \mathbb{E}^\mathbb{P}[T_t | F_s] = \bar{T_t} + (T_s – \bar{T_s}) e^{-\kappa (t-s)}\)

\(\large \mathbb{Var}^\mathbb{P}[T_t | F_s] = \int_s^t e^{-2 \kappa (t-u)} \sigma^2_u du\)

Remember: If the \(\large \frac{d\bar{T_t}}{dt}\) is included in the temperature dynamics, under a change of measure (using Girasanov’s theorem) the expectation can be shown to exactly equal to \(\large \bar{T_t}\)

Under the martingale measure \(\large \mathbb{Q}\) the equation now becomes:

\(\large T_t = \bar{T_t} + (T_s – \bar{T_s}) e^{-\kappa (t-s)} – \lambda \int_s^t e^{-\kappa (t-u)} \sigma_u du + \int_s^t e^{-\kappa (t-u)} \sigma_u dW_u\)

\(\large \mathbb{E}^\mathbb{Q}[T_t | F_s] = \mathbb{E}^\mathbb{P}[T_t | F_s] – \lambda \int_s^t e^{-\kappa (t-u)} \sigma_u du\)

\(\large \mathbb{Var}^\mathbb{Q}[T_t | F_s] = \int_s^t e^{-2 \kappa (t-u)} \sigma^2_u du\)

Evaluating these integrals in one interval

As long as volatility is constant throughout the interval than the following equations for risk-neutral mean, variance and covariance hold true.

Constant Volatility Along Interval: \(\large \sigma_n, n \in [s,t]\)

\(\large \mathbb{E}^\mathbb{Q}[T_t | F_s] = \mathbb{E}^\mathbb{P}[T_t | F_s] – \frac{\lambda \sigma_n}{\kappa} (1-e^{-\kappa (t-s)})\)

\(\large \mathbb{Var}^\mathbb{Q}[T_t | F_s] = \frac{\sigma^2_n}{2 \kappa} (1 – e^{-2 \kappa (t-s)})\)

Where: \(\large 0 \leq s \leq t \leq u\)

\(\large \mathbb{Cov}^\mathbb{Q}[T_t T_u| F_s] = e^{-\kappa(u-t)} \mathbb{Var}^\mathbb{Q}[T_t | F_s]\)

Note: Variance does not change.

Temperature Options Valuation

Method 1: Alternative Black Scholes Approach

The reference temperature for weather derivatives contracts traded at the CME is 18 degC.

\(\large \xi = \alpha(DD – K)^+\)

\(\large DD = H_n = HDD^{N} = \sum^N_n (18-T_n)^+\)

Now unfortunately in Sydney, Australia the temperature often rises well above the 18 degree threshold, however in the US during Heating Degree periods, you may often observe that for the winter months there is a very low frequency of temperatures rising above the reference temperature of 18 degC. In mathematical terms [1]:

\(\large \mathbb{P}{(18-T_n)^+ = 0} \approx 0\)

This means we can use the Gaussian OU process that under the risk-neutral measure \(\large \mathbb{Q}\) is normally distributed with \(\large \mu_t = \mathbb{E}^\mathbb{Q}[T_t | F_s]\) and \(\large \sigma^2_t = \mathbb{Var}^\mathbb{Q}[T_t | F_s]\).

\(\large T_t \sim N(\mu_t, \sigma_t)\)

Therefore \(\large DD\) simplifies to:

\(\large DD = H_n = HDD^{N} = \sum^N_n (18-T_n)^{+} \approx 18N – \sum^N_n T_n\)

\(\large \mu_t = \mathbb{E}^\mathbb{Q}[H_n | F_s] \approx 18N – \sum^N_n \mathbb{E}^\mathbb{Q}[T_n | F_s]\)

\(\large \sigma^2_t = \mathbb{Var}^\mathbb{Q}[H_n | F_s] = Var \left( \sum^N_n T_n \right)\)

\(\large \sigma^2_t = \sum^N_n Var(T_n) + 2 \sum^N_{n=1} \sum^N_{n<m} Cov(T_n,T_m)\)

References:

[1] Djehiche, Boualem, Alaton, Peter, and Stillberger, David. On modelling and pricing weather
derivatives. Applied Mathematical Finance, 9:1–20, 2002

trading_dates = pd.date_range(start='2023-06-01', end='2023-08-31', freq='D')
volatility = spline(5, days, T_std)
mc_temps, mc_sims = monte_carlo_temp(trading_dates, Tbar_params, volatility, first_ord)

print('Probability P(max(18-Tn, 0) = 0): {0:1.1f}%'.format(len(mc_sims[mc_sims[1] >= 18])/len(mc_sims)*100))

Black Scholes Approach

Risk-Neutral \(\mu\)

\(\large \mu_t = \mathbb{E}^\mathbb{Q}[H_n | F_s] \approx 18N – \sum^N_n \mathbb{E}^\mathbb{Q}[T_n | F_s]\)

\(\large \mathbb{E}^\mathbb{Q}[T_t | F_s] = \mathbb{E}^\mathbb{P}[T_t | F_s] – \frac{\lambda \sigma_n}{\kappa} (1-e^{-\kappa (t-s)})\)

Risk-Neutral \(\sigma^2_t\)

\(\large \sigma^2_t = \sum^N_n Var(T_n) + 2 \sum^N_{n=1} \sum^N_{n<m} Cov(T_n,T_m)\)

\(\large \mathbb{Var}^\mathbb{Q}[T_t | F_s] = \frac{\sigma^2_n}{2 \kappa} (1 – e^{-2 \kappa (t-s)})\)

Where: \(\large 0 \leq s \leq t \leq u\)

\(\large \mathbb{Cov}^\mathbb{Q}[T_t T_u| F_s] = e^{-\kappa(u-t)} \mathbb{Var}^\mathbb{Q}[T_t | F_s]\)

def rn_mean(time_arr, vol_arr, Tbars, lamda, kappa):
    """Evaluate the risk neutral integral above for each segment of constant volatility
    Rectangular integration with step size of days
    """
    dt = 1/365.25
    N = len(time_arr)
    mean_intervals = -vol_arr*(1 - np.exp(-kappa*dt))/kappa
    return 18*N - (np.sum(Tbars) - lamda*np.sum(mean_intervals))

def rn_var(time_arr, vol_arr, kappa):
    """Evaluate the risk neutral integral above for each segment of constant volatility
    Rectangular integration with step size of days
    """
    dt = 1/365.25 
    var_arr = np.power(vol_arr,2) 
    var_intervals = var_arr/(2*kappa)*(1-np.exp(-2*kappa*dt))
    cov_sum = 0
    for i, ti in enumerate(time_arr):
        for j, tj in enumerate(time_arr):
            if j > i:
                cov_sum += np.exp(-kappa*(tj-ti)) * var_intervals[i]
    return np.sum(var_intervals) + 2*cov_sum

def risk_neutral(trading_dates, Tbar_params, vol_model, first_ord, lamda, kappa=0.438):
    if isinstance(trading_dates, pd.DatetimeIndex):
        trading_date=trading_dates.map(dt.datetime.toordinal)
  
    Tbars = T_model(trading_date-first_ord, *Tbar_params) 
    dTbars = dT_model(trading_date-first_ord, *Tbar_params) 
    mc_temps = pd.DataFrame(data=np.array([Tbars, dTbars]).T,
                            index=trading_dates, columns=['Tbar','dTbar'])
    mc_temps['day'] = mc_temps.index.dayofyear
    mc_temps['vol'] = vol_model[mc_temps['day']-1]
    time_arr = np.array([i/365.25 for i in range(1,len(trading_dates)+1)])
    vol_arr = mc_temps['vol'].values
    mu_rn = rn_mean(time_arr, vol_arr, Tbars, lamda, kappa)
    var_rn = rn_var(time_arr, vol_arr, kappa)
    return mu_rn, var_rn

Closed-form Solution

Using the fundamental theorem of asset pricing, we have the following Risk-Neutral Pricing formulation.

\(\begin{equation}\large
\frac{C_t}{B_t} = \mathbb{E}_{\mathbb{Q}}[\frac{C_T}{B_T}\mid F_t]
\end{equation}\)

\(\begin{equation}\large
C_t = B_t \mathbb{E}_{\mathbb{Q}}[\frac{C_T}{B_T}\mid F_t]
\end{equation}\)

Call Option: \(\large \xi = \alpha(DD – K)^+\)

\(\large C_t = B_t \mathbb{E}_{\mathbb{Q}}[\alpha(DD – K)^+]\)

\(\large C_t = \alpha e^{-r(T-t)} \int^\infty_K (x-K)f_{DD}(x)dx\)

\(\large C_t = \alpha e^{-r(T-t)} \left( (\mu_t – K) \Phi(-Z_t) + \frac{\sigma_t}{\sqrt{2\pi}} e^{-\frac{Z_t^2}{2}}\right)\)

where: \(\large Z_t = \frac{K – \mu_t}{\sigma_t}\) and \(\large \Phi\) is the CDF of the Normal Distribution.

Put Option: \(\large \xi = \alpha(K – DD)^+\)

\(\large P_t = B_t \mathbb{E}_{\mathbb{Q}}[\alpha(K – DD)^+]\)

\(\large P_t = \alpha e^{-r(T-t)} \int^K_0 (K-x)f_{DD}(x)dx\)

\(\large P_t = \alpha e^{-r(T-t)} \left( (K – \mu_t) \left(\Phi(Z_t) – \Phi(-\frac{\mu_t}{\sigma_t})\right)+ \frac{\sigma_t}{\sqrt{2\pi}} \left( e^{-\frac{Z_t^2}{2}} – e^{-\frac{u_t^2}{2 \sigma^2_t}} \right) \right)\)

def winter_option(trading_dates, r, alpha, K, tau, opt='c', lamda=0.0):
    """Evaluate the fair value of temperature option in winter
    Based on heating degree days only if the physical probability that 
    the average temperature exceeds the Tref (18 degC) is close to 0 
    """
    mu_rn, var_rn = risk_neutral(trading_dates, Tbar_params, volatility, first_ord, lamda)
    disc = np.exp(-r*tau)
    vol_rn = np.sqrt(var_rn)
    zt = (K - mu_rn)/vol_rn
    exp = np.exp(-zt**2/2)
    if opt == 'c':
        return alpha*disc*((mu_rn - K)*stats.norm.cdf(-zt) + vol_rn*exp/np.sqrt(2*np.pi))
    else:
        exp2 = np.exp(-mu_rn**2/(2*vol_rn**2))
        return alpha*disc*((K - mu_rn)*(stats.norm.cdf(zt) - stats.norm.cdf(-mu_rn/vol_rn)) +
                           vol_rn/np.sqrt(2*np.pi)*(exp-exp2))


trading_dates = pd.date_range(start='2023-06-01', end='2023-08-31', freq='D')
r=0.05
K=300
alpha=2500

def years_between(d1, d2):
    d1 = dt.datetime.strptime(d1, "%Y-%m-%d")
    d2 = dt.datetime.strptime(d2, "%Y-%m-%d")
    return abs((d2 - d1).days)/365.25

start = dt.datetime.today().strftime('%Y-%m-%d')
end = '2023-08-31'

tau = years_between(start, end)

print('Start Valuation Date:', start, 
      '\nEnd of Contract Date:', end, 
      '\nYears between Dates :', round(tau,3))

print('Call Price: ', round(winter_option(trading_dates, r, alpha, K, tau, 'c'),2))
print('Put Price : ', round(winter_option(trading_dates, r, alpha, K, tau, 'p'),2))

Method 2: Monte Carlo Valuation

\(\begin{equation}\large
\frac{C_t}{B_t} = \mathbb{E}_{\mathbb{Q}}[\frac{C_T}{B_T}\mid F_t]
\end{equation}\)

\(\begin{equation}\large
C_t = B_t \mathbb{E}_{\mathbb{Q}}[\frac{C_T}{B_T}\mid F_t]
\end{equation}\)

Call Option: \(\large \xi = \alpha(DD – K)^+\)

\(\large C_t = e^{-r\tau} \mathbb{E}_{\mathbb{Q}}[\alpha(DD – K)^+ \mid F_t]\)

For each simulation \(\Large m \in M\), there is a resulting number of Degree Days in the valuation period.

  • \(\large DD = H_n = HDD^{N} = \sum^N_n (18-T_n)^{+}\)

\(\large C_t = e^{-r\tau} \frac{1}{M} \sum^M_{m=1} [\alpha(DD_m – K)^+]\)

# define trading date range
trading_dates = pd.date_range(start='2023-06-01', end='2023-08-31', freq='D')
no_sims = 10000
vol_model = spline(5, days, T_std)

def temperature_option(trading_dates, no_sims, Tbar_params, vol_model, r, alpha, K, tau, first_ord, opt='c'):
    "Evaluates the price of a temperature call option"
    mc_temps, mc_sims = monte_carlo_temp(trading_dates, Tbar_params, volatility, first_ord, no_sims)
    N, M = np.shape(mc_sims)
    mc_arr = mc_sims.values
    DD = np.sum(np.maximum(18-mc_arr,0), axis=0)
    if opt == 'c':
        CT = alpha*np.maximum(DD-K,0)
    else:
        CT = alpha*np.maximum(K-DD,0)
    C0 = np.exp(-r*tau)*np.sum(CT)/M
    sigma = np.sqrt( np.sum( (np.exp(-r*tau)*CT - C0)**2) / (M-1) )
    SE = sigma/np.sqrt(M)
    return C0, SE

call = np.round(temperature_option(trading_dates, no_sims, Tbar_params, vol_model, r, alpha, K, tau, first_ord, 'c'),2)
put = np.round(temperature_option(trading_dates, no_sims, Tbar_params, vol_model, r, alpha, K, tau, first_ord, 'p'),2)
print('Call Price: {0} +/- {1} (2se)'.format(call[0], call[1]*2))
print('Put Price : {0} +/- {1} (2se)'.format(put[0], put[1]*2))

Comparing the two methods

strikes = np.arange(180,520,20)
data = np.zeros(shape=(len(strikes),4))
for i, strike in enumerate(strikes):
    data[i,0] = temperature_option(trading_dates, no_sims, Tbar_params, vol_model, r, alpha, strike, tau, first_ord, 'c')[0]
    data[i,1] = winter_option(trading_dates, r, alpha, strike, tau, 'c')
    data[i,2] = temperature_option(trading_dates, no_sims, Tbar_params, vol_model, r, alpha, strike, tau, first_ord, 'p')[0]
    data[i,3] = winter_option(trading_dates, r, alpha, strike, tau, 'p')

df = pd.DataFrame({'MC Call': data[:, 0], 'BSA Call': data[:, 1], 'MC Put': data[:, 2], 'BSA Put': data[:, 3]})
df.index = strikes

plt.plot(df)
plt.title('Winter Temperature Options')
plt.ylabel('Option Price $USD')
plt.xlabel('Strikes (HDD)')
plt.legend(df.columns, loc=4)
plt.show()