1. Introduction
2. Approach & Methodology
3. Results
4. Future Work

Despite great advances in technology over the years, we humans remain unable to control many forces of Mother Nature, weather being one of them. Weather affects everyone and has not only a prominent role from an environmental standpoint, but also from a business one. In fact, it is estimated that nearly 20% of the U.S. economy is affected by the weather, including industries ranging from obvious ones like agriculture and energy to less transparents ones like entertainment and travel. As former commerce secretary William Daley stated in a 1998 testimony to Congress, "Weather is not just an environmental issue, it is a major economic factor. At least $1 trillion of our economy is weather-sensitive."

Weather risk is incredibly unique. It's localized, can't be modified or controlled and still cannot be precisely and consistently predicted. Until recently, there weren't many financial instruments corporations could use to hedge against weather risk. Typically, insurance has been the go-to tool for protection against weather conditions, but it has one big problem: it only provides protection against calamitous disasters. It doesn't have the capability to protect against slight changes in weather that can still have adverse profitability consequences. And thus, in the late 1990s, the weather derivative was born.

Weather derivatives come in multiple flavors; common underlyings include rainfall or temperature (we'll be pricing these). The general structure of a weather derivative with a temperature (˚F) underlying is relatively simple: weather is quantified and indexed in terms of monthly or seasonal average temperatures. Following that, a dollar amount is then attached to each index value, and this package is traded on an exchange (we'll be looking at weather futures trading on the Chicago Mercantile Exchange, or CME).

Weather contracts in U.S. cities are tied to an index of heating degree day (HDD) or cooling degree day (CDD) values. Both are calculated according to how many degrees a day’s average temperature varies from a baseline of 65° F, where the average temperature is based on the maximum and minimum temperature of the given day.

An HDD value equals the number of degrees the day's average temperature is lower than 65° F. For example, a day's average temperature of 40° F would give you an HDD value of 25 (65 - 40). If the temperature exceeded 65° F, the value of the HDD would be zero. This is because in theory there would be no need for heating on a day warmer than 65°. As a result, this index is often used for winter months.

A CDD value equals the number of degrees the average daily temperature exceeds 65° F. For example, a day's average temperature of 80° F would give you a daily CDD value of 15 (80 - 65). If the temperature were lower than 65° F, the value of the CDD would be zero. Again, remember that in theory there would be no need for air conditioning if the temperature were less than 65°F, and thus this index is typically used for summer months.

These index values can then be aggregated over a given time period. Thus, a monthly HDD or CDD index value is simply the sum of all daily HDD or CDD values recorded that month. Let it be noted that HDD and CDD values must be positive, as all negatives are simply recorded as 0. For weather derivatives, an expected total HDD or CDD is established for the relevant time period. The party that will benefit more will be the one that successfully predicts whether the actual value will be above or below the index. Therefore, the accuracy of the temperature model is the key component to achieving profit.

The contract price of the future can be determined by multiplying the expected HDD or CDD value by a specified tick size. The tick size is simply a dollar amount per HDD or CDD and can vary based on how much money the parties are looking to exchange. For example, the value of a CME weather futures contract is determined by multiplying the monthly HDD or CDD value by a tick size of $20.

According to the former executive director of the Weather Risk Management Association, Valerie Cooper, an $8 billion weather derivatives industry developed within a few years of its inception. While current users of the derivative are primarily energy-related companies, there has been increased interest from agricultural firms, restaurants and tourism/travel companies. This can be seen in the numbers: The volume of traded CME weather futures in 2003 more than quadrupled from 2002, totaling roughly $1.6 billion in notional value. The momentum of this volume continues to increase as well, meaning we should expect to see continued growth and interest in this exotic derivative.

Shiny app: https://shalini-s.shinyapps.io/weather-app/

• Model futures by setting a specific contract range and tick size

Alternatively, if you want to use the code yourself, run in the following order:

• $ python scrape_data.py
• $ python clean_data.py
• $ python avg_temperature_model.py or avg_temperature_model.R in RStudio

We used a service called Wunderground - a "commercial weather service providing real-time weather information via the Internet" - to gather data. Using a web scraper (please see scrape_data.py) that was fortunately already developed by someone else specifically for this purpose, we harvested 8 years worth of temperature data ranging from July 30th, 2010 to July 30th, 2018. The data was derived from the "KCASANFR49" weather station, a small station situated in the Mission District within San Francisco (this goes back to the previous point in the Introduction stating that weather data is localized). We could have passed in other stations to use as well, but for the purposes of this demonstration, we just stuck to one station within San Francisco. Some undesirable characteristics of the data were that it:

• included extraneous information like humidity, rainfall and much more
• was recorded in 15 minute intervals throughout the day (in other words, the temp_data_raw.csv file generated after scrape_data.py was run has ~500k+ observations)

The first point was combated by just appending temperature data in Fahrenheit as a column to the pandas DataFrame used, and the second point was fixed in the next step of the process.

Using pandas, aggregation was done by day reducing the total number of observation from ~500k+ to 2923. Three columns were created through this aggregation process: min, max, and mean. min recorded the minimum temperature value from each day, max recorded the maximum temperature value from each day and mean recorded the average temperature value from each day. The data was now grouped together and easier to work with, but no weather station is perfect and "KCASANFR49" is no exception to that rule. There were a couple days where the min and mean columns displayed either extreme negative values (-999˚F) or NaN values. Using the forward fill pandas method, negative values were turned into NaNs and NaNs were turned into the previous valid data point. The resulting DataFrame was then exported as temp_data_cleaned.csv.

We have time series data, and a popular model used to forecast time series data is the autoregressive integrated moving average model, or ARIMA for short. At a high level, an ARIMA model makes use of past data to model the existing data as well as to make predictions of future behavior.

Breaking down the process of determining an accurate ARIMA model to fit the data, we have to:

• Visualize the time series data
• Test the time series data for stationarity
• Plot the autocorrelation and partial autocorrelation charts
• Determine if seasonality or differencing would improve the model
• Construct the ARIMA model
• Use the model to make predictions

Looking at our plot, we can see that the temperature data almost resembles a sinusoid. This is due to the four seasons taking on a cyclical pattern. From an initial glance, the data appears to have little to no trend component and be stationary, but this still has to be tested.

Augmented Dickey-Fuller (ADF) Unit Root Tests are used to test the stationarity of data. Passing in df["mean"], we get the following output:

The key statistic here is the p-value. Roughly speaking, if the p-value is less than 0.05, we can make the assumption that our data is stationary. In this case, it is. In the avg_temperature_model.py code, we do indeed do some additional differencing and run the ADF test on df["Seasonal First Difference"] but we don't end up using this when fitting our model so it can be ignored.

Using a maximum lag of 30, we see that we have a gradual, downward-sloping autocorrelation and a sharp drop-off in the partial autocorrelation. This indicates that both AR and MA terms will be required in the model. An eacf chart was examined in order to determine the exact orders. This plot showed significant lags at (1,2), (1,1), and (4,0). Each of these models were tested, and the model with the smallest AIC value was selected.

Next, the model was examined and other methods of making the model more accurate were analyzed. The data was differenced once and was determined to still be stationary (by ADF testing). The eacf chart of the differenced data was computed and each of the resulting models was compared through AIC values to determine the optimal model. To make sure that we were not overfitting the model, the standard deviation of both models was calculated. The sd of the differenced model was less than the sd of the original model, which does not indicate overfitting. Furthermore, the lag-1 autocorrelation of the model was greater than -0.5, which also does not indicate overfitting.

To account for global warming, a seasonal difference was tested. The seasonal difference was determined to be stationary, and the same procedure from above was applied to avoid overfitting. The frequency of the seasonality was determined to be 365, which indicates a comparison between the temperature at a specific date and the temperature at that same date during each year in the data. The eacf of the seasonally differenced data was computed and no terms were determined to be significant.

Given the seasonality of our data, we used the SARIMAX functionality from the statsmodels library in Python and the standard ARIMA in our R script. From the analysis above, the p, d and q values for the order were 1, 1 and 2 respectively and P, D, Q and frequency values for the seasonal order were 0, 1, 0 and 365 respectively. These values were picked using the eacf functionality in R. After deriving them with R, they were hard-coded into the Python model.

Furthermore, R was used for a forecast over all the existing data as well as one year into the future. While this functionality was added to the Python model in the form of imported csv files (forecasted_existing.csv and forecasted_unknown_1y.csv) retrieved from the R model, SARIMAX can't handle daily period lags very well and does better with monthly or quarterly data. As a result, the .predict methods cause the program to crash because there isn't enough RAM to do computations on so many dense arrays, and there is currently no seasonal ARIMA version in Python that uses sparse arrays for the same purpose (statespace models are optimized for smaller arrays using dense LAPACK functions). Perhaps it would work on a computer with 32GB or 64GB of RAM, but this hasn't been tested at the time of writing.

Residual plot from our model on existing data:

Residual plot showing kernel density estimation:

Forecast across existing time values:

Forecast for the next year (July 31st, 2018 - July 31st, 2019)

Using a hypothetical scenario, let's say the month is currently November in 2018 and we think the temperature in December is going to be particularly low. We are a farmer and are worried that low temperatures will cause more of our crops to die on average than usual, resulting in a decrease in profitability for the month. We would like to hedge this risk by purchasing a future with heating degree days (HDDs) set as the underlying.

The contract range for this future is December 1st, 2018 - December 31th, 2018. It's trading on the CME which means the tick size is $20 ($20 per index value). Our model predicts that there will be 365.264087215738 HDDs during this time period. As a result, we predict that a single futures contract will be worth $7305.29. The next step is to check what one of these contracts is currently going for on the CME. If there are any that are listed for a price less than our predicted one, we can recognize an arbitrage opportunity and purchase one or more contracts.

This project is ongoing as opportunities for future work remain. Currently, we're only pricing futures. The obvious next step is to use an options pricing model like the Black-Scholes model or a Monte Carlo simulation alongside our ARIMA model to price options. Furthermore, we may look at alternative modeling approaches other than our current ARIMA model. A cubic spline interpolation on three dimensional data (years, days, temperature) prior to an ARMA process may be more accurate as it can smooth the global warming effect across years. Additionally, using more than our chosen 8 years of data will help improve model robustness and yield a more accurate result. Finally, a distribution of potential payouts could be constructed to help the investor determine whether the contract fits in the desired risk profile. Finding the probability of payouts given our predicted CDD and HDD values can be more useful than simply returning the fair price of the derivative.