All forecasts are wrong, but some generate fewer complaints

In the engineering business, it’s common to be asked for time series forecasts.   Production forecasts, man-hour forecasts, sales forecasts, the list is endless.   In those cases where the only available information is the time series itself, the “go to” forecasting method is commonly a “Black Box”, or some form of the ARIMA (autoregressive integrated moving average) technique.  These techniques do a good job of torturing usable patterns out of the original data, however they don’t provide insights into why the patterns exist.          

If additional information is available, it can be worthwhile to explore its use in the forecast.  An underlying real world model structure can not only improve the forecast, but can boost the user’s confidence in the results.

In the following example, the Electric Reliability Council of Texas (ERCOT) monthly average load data, which is the total electric demand on the Texas grid, will be used to generate a “Black Box” forecast.  Those results will be compared to a regression model forecast.

As shown in Graph 1 below, the gray time series is the ERCOT monthly average load from 2003 to 2012.  It is a 31 day, end of month average, that exhibits seasonality and a long term trend.

Graph 1 – ARIMA 48 month forecast of the ERCOT Monthly Average Load

The solid blue line is a 48 month forecast, generated from standard ARIMA methods, allowing for a non-zero drift term.   The dashed lines are the upper and lower 1% tails, such that the actual ERCOT load is not expected to exceed either tail more than 1% of the time.

A common question about this graph is: Compared to the past, isn’t the forecast summer upper 1% range too high?  The hot 2011 summer was at least a 1% event, yet the dashed 1% forecast line, if extended backward, appears to be well above that level.  As you might expect, there isn’t an answer to that question that is acceptable to most people.   The results are simply the results, and a long drawn out explanation of the ARIMA details doesn’t help.

A better solution can sometimes be a model that is designed according to common sense, to be used for common sense.  The following video describes a simple regression model relating temperature, and economic activity to ERCOT’s monthly average load.

In Graph 2, the forecast in the video is compared with the original ARIMA forecast.  It shows a significant difference in both the summer 1% bands, and the 2013 winter 1% bands.  Which forecast do you believe, and why?

Graph 2 – Comparison of the ARIMA forecast with the simple regression forecast that was discussed in the video.

It turns out that the uncertainty in Texas temperatures determines most of the 1% spread, such that the faith in the video model is the direct result of the stability in the Texas climate.

The bottom line?   If the supporting data is available, it might be worthwhile to use it.


Links to the data sources referenced in the video:

http://ercot.com/gridinfo/load/load_hist/

http://research.stlouisfed.org/fred2/series/TXNA

http://www.srh.noaa.gov/images/ewx/sat/satmontemp.pdf


The following is the R script that was used to fit and forecast the ARIMA model, and the regression model that was discussed in the video.

#==========================================================================
# revgr.com
#
# The following R script is used to generate two models that fit the
# ERCOT Monthly Average Load Data.  The first model is a typical
# ARIMA model, allowing for a non-zero drift term.   The second
# model is a simple regression model that uses temperature and
# economic activity as inputs.
#==========================================================================

#Load the necessary libraries

library(zoo) #Used to represent time series data
library(DEoptim) #Used to determine the "Comfortable Temperature Zone"
library(forecast) #Used to generate the ARIMA model and forecast

#San Antonio monthly average temperature starting January 1885
SATem1 <- c(52.3, 56.7, 59.4, 70, 71.5, 81.2, 80.5, 83.2, 79, 66.9, 61.6, 
           55.1, 44.4, 55.3, 59.3, 66.4, 75.4, 80.2, 83.2, 83.9, 78, 69.1, 
           59.3, 53.5, 50.7, 59.4, 65.6, 69.2, 75.3, 79.7, 84.7, 83.9, 78.3, 
           67.6, 59.9, 49.7, 48.4, 58, 57.5, 70.5, 73.9, 79.3, 82.5, 82.9, 
           75.3, 70.5, 57.2, 55.6, 51.6, 55.1, 59.4, 69.1, 72.7, 79.1, 81.7, 
           80.7, 73.6, 69.5, 55.8, 66.1, 59.4, 60.8, 62.6, 67.9, 74.4, 78.5, 
           83.4, 83.2, 76.1, 70.1, 60.9, 56.7, 51, 57.4, 57.3, 67.4, 71.8, 
           81.7, 84.4, 82.8, 78.6, 68.9, 61.1, 55.1, 47.1, 62.3, 58.9, 73.7, 
           78.3, 82, 84.2, 80.9, 77.9, 71.9, 61.9, 53.2, 53.1, 55.6, 62.3, 
           74.3, 75.3, 81, 85.2, 83.8, 81.8, 71.1, 58.9, 58.8, 55.5, 51.9, 
           65.1, 74.1, 77.1, 80.9, 84, 80.1, 79.4, 73.4, 62.1, 58.4, 54.3, 
           45.2, 62.6, 70, 72.9, 80.2, 84.7, 84.5, 80.8, 68.8, 58, 54.4, 
           53.9, 56.7, 60.9, 70.3, 78.9, 84, 83.8, 84.8, 78.7, 68.2, 62, 
           56.5, 49.3, 59.9, 67.6, 69.2, 74.8, 81.4, 85, 82.4, 77.5, 72.6, 
           62.9, 51.1, 56, 58.9, 62.4, 69, 76.4, 80.6, 83.1, 84, 79.6, 71.2, 
           56.9, 47.4, 50.9, 45.8, 65.6, 68.8, 78.4, 79.6, 82.8, 85.8, 78.1, 
           72.8, 61.3, 53.7, 53.3, 54.4, 61.9, 67.9, 74.2, 82.8, 82.4, 82.8, 
           82.5, 72.7, 62.4, 54.2, 56.4, 52.7, 62.9, 67.8, 76.5, 82.7, 84.1, 
           85.8, 78.4, 72.3, 63.3, 53.5, 51.4, 53.5, 64.4, 71.7, 77.4, 83.7, 
           83, 85.8, 78.6, 70.6, 64.8, 53.2, 51.4, 50.1, 59.2, 67.1, 72, 
           75.8, 79.8, 81.2, 76.7, 68.4, 60.3, 54.2, 51.6, 59.7, 69.1, 69.6, 
           73.6, 79.8, 81, 82.1, 79.8, 71, 60.5, 52.8, 49.9, 43.9, 65.3, 
           67.5, 77, 80.2, 81.3, 84.8, 80.5, 70.6, 63.9, 48.5, 52.4, 53.3, 
           58.5, 68.8, 75.6, 83.4, 82.5, 82, 80, 66.1, 60.7, 58.3, 60.5, 
           58.6, 71.1, 67.7, 71, 82.7, 83.3, 85, 81.5, 71.8, 56.5, 54.7, 
           54.7, 55.3, 67.8, 70.2, 76.3, 82.7, 82.7, 82.7, 77.2, 68.1, 62.3, 
           56.6, 57.2, 58.5, 65.5, 68.9, 75.2, 82.8, 85.8, 84.5, 79, 71.1, 
           68.2, 47.8, 54.9, 51.6, 67.8, 69.5, 74.7, 81.3, 84.8, 86.1, 81.5, 
           69.1, 63.4, 55.4, 59.1, 60.1, 66, 68.7, 74, 84.2, 84.1, 86.3, 
           84.6, 70.5, 56.1, 51.6, 48.9, 51.4, 56.7, 68.6, 76.8, 78.4, 85.1, 
           86, 81.6, 71.9, 59.7, 49.9, 52.4, 52.1, 59.4, 66.9, 75.6, 79.3, 
           84.5, 84, 75.6, 67, 66.2, 52.4, 56.4, 53.2, 58.9, 66.8, 74.4, 
           82, 85.6, 82.7, 79.6, 70.2, 61.2, 46.3, 50.6, 58.5, 53.2, 67.5, 
           75.7, 83.9, 84.8, 82.5, 79.8, 72.2, 63.8, 56.8, 56.2, 58.7, 68.7, 
           67.6, 76.1, 84.1, 82.8, 82, 78, 70.9, 59, 54.7, 55.2, 57.7, 63.6, 
           69.1, 71.6, 83.1, 84.8, 85.7, 79.3, 68.7, 63.5, 50.4, 47, 56.6, 
           66.6, 68.9, 75.9, 83.7, 85.3, 85.1, 76.4, 71.6, 57.4, 54.7, 49.7, 
           53, 61, 68.4, 73.4, 77.5, 80.7, 82.2, 77.9, 73.9, 60.8, 51.7, 
           49.4, 57.7, 60.4, 69.5, 76.8, 78.5, 83.8, 82.9, 82.1, 71.1, 57.5, 
           54.7, 58.4, 58.4, 67, 67.5, 75.4, 81, 83.6, 85.2, 81.7, 70.4, 
           65.6, 59, 49.5, 58.2, 61.7, 70.2, 77.1, 79.4, 84.1, 85.9, 79.7, 
           71.4, 61.9, 59.7, 62, 52, 58.7, 69.2, 77.2, 83.2, 83.6, 84.2, 
           79.3, 68, 58.7, 55.6, 47.4, 54.4, 58.6, 69.2, 71.6, 81.2, 82.5, 
           86, 78.5, 72, 66.5, 52.7, 50.5, 61.5, 67.7, 74.2, 77.3, 84.6, 
           86.1, 84, 80.8, 68.7, 60.2, 49.7, 50.2, 61.2, 58.7, 65.7, 74.6, 
           81.6, 82.3, 84.4, 82.5, 75.2, 59.6, 54.5, 55.8, 61.6, 63.3, 73.3, 
           80.5, 80, 84.9, 86, 80.5, 73.3, 69.8, 50.4, 53.1, 55.6, 65.6, 
           65.8, 74.7, 81.8, 85, 86.6, 77.4, 74.1, 59.4, 53.5, 53.8, 50, 
           65.9, 72.8, 74.7, 81.7, 81.9, 85.1, 80.3, 72, 55, 54.6, 43.3, 
           61.6, 58.9, 72.2, 75, 79.6, 84.3, 85.9, 80.7, 69.8, 59.7, 51, 
           52.6, 57.4, 57.9, 64.2, 72.5, 81.3, 82.6, 82.1, 83.3, 78, 66.6, 
           55.2, 55.4, 63, 57.4, 70.5, 75.5, 82.2, 85.3, 84.1, 77.3, 69.2, 
           57.8, 50.3, 60, 52.9, 65.9, 71.6, 79.8, 80.9, 85.1, 83.3, 82.9, 
           75.6, 65.3, 62.5, 55.2, 55.5, 61, 70.9, 76.7, 85.2, 84.8, 84.6, 
           80.2, 76.4, 65.6, 56, 57, 56.3, 67.3, 70.4, 72.6, 79.5, 83.7, 
           84.9, 77, 74, 60.3, 52, 52.1, 51.1, 66.3, 68.1, 74, 82.5, 81.6, 
           83.8, 79.3, 66.3, 57.9, 56.5, 50.4, 56.3, 58, 70, 76.9, 82.5, 
           83.6, 86.1, 82, 72.6, 58.5, 52.9, 54.5, 61.3, 68.3, 66.9, 75.2, 
           82.4, 85.1, 84.6, 79.7, 73.1, 58.5, 55.9, 54.9, 54.1, 65.4, 70.9, 
           79.2, 83.1, 85.2, 83.4, 81.2, 73.5, 58.5, 58.8, 43.9, 56, 64.1, 
           68.8, 75.5, 79.6, 83.7, 85.1, 79, 72.3, 59.5, 56.8, 55.7, 52.8, 
           56.1, 68.6, 76.3, 81, 84.4, 84.9, 81.6, 75.8, 59.1, 54.9, 50.2, 
           53.9, 61.5, 69.7, 76, 83.4, 81.3, 82.7, 75.6, 69.6, 63, 55.1, 
           50.3, 58.4, 58.3, 71.8, 77.2, 81.3, 83.4, 85.9, 76.4, 67.7, 57.1, 
           50.5, 50.3, 57.1, 60.5, 68.1, 71.8, 81.3, 84.5, 84, 78.1, 69.7, 
           60.3, 48.9, 51.5, 55.6, 67.2, 67.2, 75.2, 82.4, 84.2, 84.6, 80.3, 
           68.4, 63.9, 50.9, 49.2, 55.6, 63, 71.1, 75.4, 80, 84.4, 83.8, 
           78.2, 72.4, 60.4, 55.9, 48.5, 48.6, 56.1, 68.2, 75.3, 83.9, 84.2, 
           83.3, 80.4, 76.3, 57.4, 53.5, 45.4, 52.7, 60.4, 72.9, 78.2, 84.4, 
           85.4, 85.3, 77.7, 70.4, 57.7, 55.9, 46.5, 57.5, 62.3, 64.6, 78.4, 
           81.1, 83.6, 82.5, 80.8, 70.6, 62.7, 55.9, 58.4, 56.7, 60.6, 67.1, 
           76.9, 80.7, 84.2, 82.5, 79.5, 73.3, 59.4, 52.8, 50.5, 54.3, 62.6, 
           68.9, 75.4, 81.9, 86.7, 87.1, 80.5, 72.9, 58.1, 55.2, 59.5, 58.6, 
           60.7, 65.8, 73.3, 81.9, 83.6, 86.2, 77.5, 65.6, 58.4, 51.6, 55.9, 
           54.8, 67.8, 69.2, 76, 85, 86.1, 84.6, 78.3, 71.2, 59, 49.4, 54.9, 
           60.4, 62.3, 73.7, 74.9, 83.5, 85.8, 86.7, 82.6, 73.1, 60.6, 58.3, 
           53.1, 56.4, 63.4, 73.8, 78.9, 81.4, 84.5, 84.6, 81.9, 70.5, 59.8, 
           53.5, 52.1, 57.1, 63.3, 69.7, 78.9, 84.7, 85.5, 84.8, 80, 73, 
           57.8, 56.9, 53.3, 62.9, 61.9, 66.3, 73.4, 81, 85.7, 86, 77.8, 
           66.9, 57.1, 55.5, 50.1, 50, 55.8, 67.7, 75.1, 82.8, 84.6, 84.8, 
           79.2, 67.7, 60.8, 50.1, 48.9, 54, 60, 65.9, 77.1, 82.9, 84, 84.4, 
           81.4, 70.4, 53.3, 54, 50.1, 49.9, 56.1, 69.7, 74.1, 83.3, 84.2, 
           83.6, 78.6, 73.3, 62.2, 50.1, 47.9, 55.9, 65.7, 68.5, 78.5, 81.3, 
           82.6, 82.5, 80.5, 71.1, 58, 54.2, 45.9, 62.8, 59.1, 69.7, 77.9, 
           82.3, 86.9, 87.5, 80.9, 75.5, 60.4, 52.2, 46.2, 52.6, 65.6, 74.6, 
           77.7, 83.4, 85.4, 85.7, 81.1, 74.1, 62.4, 45.7, 51, 49.8, 61.5, 
           70.5, 77.6, 82.4, 86.3, 86.2, 80, 66.4, 62.6, 52.3, 54.4, 49.8, 
           54.9, 71.6, 75, 81.6, 84.9, 84, 80.7, 66.8, 64.5, 55.5, 45.4, 
           49.8, 60, 68.6, 73.5, 78.8, 84.2, 81.9, 77.5, 67, 63, 50.7, 50.2, 
           51.8, 66.9, 76.6, 76.6, 84.5, 85.3, 82.7, 75.5, 66.9, 60.5, 51, 
           49.8, 48.3, 58, 68.1, 75.3, 80.5, 82.7, 84.2, 76, 72.2, 56.4, 
           50.7, 52.5, 53.6, 54.9, 69, 73.5, 81.2, 86.8, 85.7, 79.6, 69.8, 
           58.1, 55.1, 45.6, 54.8, 56.8, 70.2, 72.9, 80.7, 84, 85.7, 81.1, 
           67.7, 58, 60.1, 56, 57.4, 64.6, 69.4, 78.1, 83.6, 85.9, 81.2, 
           80.1, 73.9, 63.2, 57.2, 52.8, 56.7, 66.3, 73.7, 72.8, 80.3, 82.2, 
           82.1, 82, 71.9, 54, 50.3, 47.2, 51.9, 66.1, 66, 74.7, 79.2, 83.2, 
           82.1, 79.3, 72.5, 65.8, 52.2, 51, 56.5, 67.9, 69.7, 77.3, 79.4, 
           83, 81.2, 72.3, 68.2, 57.3, 50.9, 53.2, 53.5, 61.4, 68.4, 73.5, 
           80, 80.9, 81.7, 76, 71.1, 60.3, 53.1, 49.6, 61.2, 63.8, 68.9, 
           71.3, 79.8, 79.8, 81.6, 77.5, 61.1, 52.1, 49.9, 44.1, 52.8, 61.8, 
           66.9, 74.8, 81.5, 84.9, 84.7, 82.3, 71.2, 61.4, 53.4, 43.4, 46.4, 
           59.6, 68.9, 77.1, 82.7, 86.1, 83.1, 78.5, 69.3, 62.4, 51.7, 43.7, 
           52.4, 63.3, 69.7, 73.9, 80.9, 84.7, 83.1, 78.7, 74.7, 58.2, 55.4, 
           52.6, 53.7, 61.5, 67.6, 76.1, 85.1, 88.1, 85.3, 83.7, 70.7, 58.3, 
           55, 50.8, 53.7, 60.7, 72.9, 75.3, 81.5, 84.2, 84.7, 78.9, 71.9, 
           62.4, 53, 50.8, 49.7, 63.1, 66.9, 74.5, 81.6, 85.5, 86, 80.1, 
           69.3, 59.4, 52.4, 48.9, 52.1, 58.7, 65.2, 73.6, 79.2, 82.9, 84.5, 
           78.5, 70.9, 62.5, 43, 46.7, 54.1, 64.2, 69.7, 77.1, 82.8, 85, 
           84.7, 77.6, 71.2, 58.8, 59.6, 44.2, 50.5, 64.1, 69.4, 76.7, 80.2, 
           82.2, 85.5, 79.4, 71.7, 64.4, 49.9, 53.4, 58, 62.9, 72.6, 74.6, 
           81.5, 85.8, 85.7, 83.7, 69.7, 59.4, 51.6, 50.7, 55.9, 57.8, 66.1, 
           75.8, 80.5, 83.8, 86, 79.2, 71.2, 60.6, 54.2, 47.6, 54.3, 61.3, 
           69.1, 76.1, 81.2, 84.6, 86.4, 80.7, 73.2, 65.1, 56, 56.2, 51.6, 
           61.9, 70.4, 81.7, 83.3, 86.6, 86, 79.1, 71.3, 61.8, 43.4, 56.4, 
           58.9, 61.5, 69.7, 79.3, 87.5, 83.4, 85.3, 80, 69.3, 63, 51.9, 
           48.9, 56.6, 64, 72.4, 77.7, 82.8, 84.5, 85.8, 77.8, 73.3, 57.4, 
           55.5, 50.8, 59.1, 63.3, 69, 73.7, 82.5, 84.7, 82.2, 81.7, 73.4, 
           57.3, 56.3, 51.2, 55.5, 61.5, 67.3, 73.9, 81.6, 86.1, 87.3, 81.5, 
           70.7, 56.3, 55.1, 52.3, 56.2, 63.9, 69.8, 76, 84.5, 87.9, 86.1, 
           78.4, 72.7, 64.8, 57, 53.5, 57.4, 61.9, 69.8, 78.6, 79.3, 84.3, 
           85.5, 80.1, 69.8, 59.5, 55.7, 51.1, 57.9, 57.6, 69.5, 81.9, 84.1, 
           87.3, 84.4, 78.4, 71.1, 61.3, 54.5, 49.2, 53.1, 63.3, 63.9, 74, 
           79.8, 85.1, 86.1, 82.2, 70.2, 57.3, 50.2, 56.4, 55.3, 59.8, 66.7, 
           79.8, 86.3, 88.1, 83.6, 80.5, 71.3, 62.4, 52.7, 54.6, 61.8, 62.7, 
           71.2, 76.2, 81.9, 82.9, 86.1, 80.3, 69.6, 63.1, 54, 55.2, 62.6, 
           67, 70.7, 78.6, 81, 85.9, 86.3, 81, 71.1, 56.9, 46.4, 49.2, 57.5, 
           56.6, 70.8, 76.3, 82.6, 85.4, 85.6, 76.9, 67.9, 62.9, 53.8, 53.7, 
           50.8, 60.3, 73.2, 76.8, 83.4, 82.5, 85.3, 78.7, 70.7, 57.8, 53.8, 
           50.1, 53.1, 60.6, 71.7, 80.3, 81.7, 81.9, 83.7, 76.7, 70.6, 63, 
           53.9, 54.5, 52.6, 65.9, 67.2, 76.1, 80.8, 82.9, 83.3, 80.5, 76.9, 
           61.1, 53.1, 55.9, 56.3, 61.3, 68.4, 75, 82.6, 85.3, 85.7, 84.3, 
           70.9, 64.9, 53, 58.2, 55.9, 67.5, 76.7, 79, 83.6, 85.7, 88.3, 
           79.7, 72.4, 63.8, 54.4, 48.3, 54.8, 65, 65.2, 75.3, 80.7, 80.4, 
           83.7, 80.2, 73.1, 62.7, 56.1, 51.8, 61.7, 64.5, 70.6, 80.1, 86.8, 
           84.1, 84.4, 79.5, 71.4, 63.7, 55, 54.5, 62.9, 65.1, 69.4, 79.5, 
           86.3, 88.7, 88.3, 78.4, 69.9, 60.7, 48.3, 49.7, 49.4, 59.3, 68.6, 
           77.5, 83.5, 84, 87.5, 80.1, 70.2, 62.1, 53.8, 50.5, 55.4, 66.8, 
           75.7, 78.6, 86.2, 87.9, 90, 82.9, 71, 62.9, 53.4, 56.2, 57.4, 
           66.4, 73.8, 78.1, 84.8, 85.4, 87.2, 79.6)

#Texas monthly Nonfarm Payrolls starting May 2003
TexasNFP1 <- c(9359.5, 9352.7, 9345.3, 9356.7, 9370.4, 9374.9, 9384.4, 9393.8, 
              9419, 9429.1, 9440, 9459.9, 9469.2, 9480.1, 9500.6, 9520.3, 9515.8, 
              9563.8, 9569.2, 9583.7, 9605.2, 9620.7, 9639.3, 9679.4, 9693.6, 
              9697.6, 9763.5, 9778.3, 9807.6, 9823.3, 9867.9, 9887.6, 9924.3, 
              9946.1, 9985.8, 9998, 10024.6, 10049.2, 10058.1, 10101.5, 10136.7, 
              10150.9, 10177.5, 10212.5, 10231.7, 10268.6, 10312.3, 10334.6, 
              10359.3, 10393.6, 10414.7, 10430, 10447.7, 10482.5, 10507.5, 
              10529.3, 10562.2, 10594.8, 10584.4, 10606.5, 10609.8, 10619.7, 
              10634.9, 10639.7, 10603.9, 10618.3, 10609, 10578.7, 10515.8, 
              10465, 10398.7, 10329.6, 10310.6, 10281.4, 10248, 10227.3, 10219.1, 
              10213, 10214.9, 10212.1, 10229.2, 10245.2, 10271.1, 10290.4, 
              10351.1, 10354.3, 10346.3, 10360.9, 10370.2, 10407.6, 10417.4, 
              10434.5, 10452.2, 10464.9, 10496, 10529.8, 10532.6, 10551.9, 
              10581.5, 10586.1, 10594.9, 10605.8, 10618.7, 10643.2, 10710.9, 
              10730.8, 10742.4, 10757.2, 10770.8, 10786.2, 10807.6)

#ERCOT monthly load starting May 2003
ERCOT1 <- c(35302.697705803, 38306.4197031039, 40790.1848852901, 42058.1983805668, 
           34595.6680161943, 29973.338731444, 28641.6477732794, 29102.9770580297, 
           29770.6180836707, 30092.016194332, 27212.9527665317, 28358.0904183536, 
           33172.0566801619, 37295.036437247, 40477.5060728745, 40304.1902834008, 
           37456.9338731444, 33293.5425101215, 28436.9392712551, 31177.7854251012, 
           30726.5789473684, 30241.6720647773, 28351.9784075573, 28864.3414304993, 
           33719.1889338731, 41634.9419703104, 41865.3873144399, 43151.1443994602, 
           40773.5721997301, 31928.004048583, 29733.2145748988, 31483.2739541161, 
           29049.0458839406, 30937.669365722, 29450.3603238866, 33063.0823211876, 
           37119.9568151147, 41875.7395411606, 44004.6032388664, 46722.975708502, 
           38431.983805668, 32980.9649122807, 30039.624831309, 31890.2685560054, 
           35506.1538461538, 32992.7638326586, 29342.8798920378, 29775.6680161943, 
           34041.5897435897, 39308.7260458839, 39520.2429149798, 44312.7246963563, 
           39644.7719298246, 34195.8380566802, 30794.6099865047, 32199.632928475, 
           33962.7683099593, 31211.9402442968, 30430.910197913, 31094.7171275795, 
           37402.3486293477, 43990.5906599034, 44056.108160517, 43287.7098482178, 
           35784.830165608, 31924.0328127956, 29270.8523495385, 32507.8446495164, 
           32416.1110773279, 29973.6791068232, 29824.9746483131, 29999.4942956815, 
           35264.5936004318, 43189.956772529, 45379.1591635047, 45009.3181117328, 
           36612.6374397773, 31205.5939856572, 28569.2043987652, 34694.8886355587, 
           35058.6019490836, 35673.8250112645, 29671.5797991647, 29893.8308174075, 
           36682.9767228246, 44077.1310952496, 43589.0989841066, 48128.7481663199, 
           40326.9601261012, 32198.4185196505, 30156.4563538353, 32348.9193960594, 
           34784.5079397179, 34821.0702141767, 30508.9018582524, 34090.5066567746, 
           37378.7814707747, 46646.7495415263, 48900.5389791592, 51495.7376787787, 
           41967.9449820203, 33263.9903463873, 30949.573562475, 33778.3111987746, 
           32470.1768471512, 31937.0263402591, 31437.6482162038, 34068.7761726991, 
           38859.650306533, 44701.9764394535, 45822.4162995182)

#Define the 4 Texas Nonfarm Payroll scenarios that will be used for forecasting
#From August 2012 for the next 48 months
TexasNFPs1 <- c(10807.8, 10824.0117, 10840.2234, 10856.4351, 10872.6468, 10888.8585, 
                10905.0702, 10921.2819, 10937.4936, 10953.7053, 10969.917, 10958.490003125, 
                10947.06300625, 10935.636009375, 10924.2090125, 10912.782015625, 
                10901.35501875, 10889.928021875, 10878.501025, 10867.074028125, 
                10867.45797522, 10885.55833827, 10903.65870132, 10921.75906437, 
                10939.85942742, 10957.95979047, 10976.06015352, 10994.16051657, 
                11012.26087962, 11030.36124267, 11048.46160572, 11066.56196877, 
                11084.66233182, 11102.76269487, 11120.86305792, 11138.96342097, 
                11157.06378402, 11175.16414707, 11193.26451012, 11211.36487317, 
                11229.46523622, 11247.56559927, 11265.66596232, 11283.76632537, 
                11301.86668842, 11319.96705147, 11338.06741452, 11356.16777757)

TexasNFPs2 <- c(10807.8, 10826.71365, 10845.6273, 10864.54095, 10883.4546, 
                10902.36825, 10921.2819, 10940.19555, 10959.1092, 10978.02285, 
                10996.9365, 11015.85015, 11034.7638, 11053.67745, 11072.5911, 
                11091.50475, 11110.4184, 11129.33205, 11148.2457, 11167.15935, 
                11186.073, 11174.4208406249, 11162.76868125, 11151.116521875, 
                11139.4643624999, 11127.8122031249, 11116.16004375, 11104.5078843749, 
                11092.8557249999, 11081.203565625, 11081.59507818, 11100.05209863, 
                11118.50911908, 11136.96613953, 11155.42315998, 11173.88018043, 
                11192.33720088, 11210.79422133, 11229.25124178, 11247.70826223, 
                11266.16528268, 11284.62230313, 11303.07932358, 11321.53634403, 
                11339.99336448, 11358.45038493, 11376.90740538, 11395.36442583)

TexasNFPs3 <- c(10807.8, 10829.4156, 10851.0312, 10872.6468, 10894.2624, 10915.878, 
                10937.4936, 10959.1092, 10980.7248, 11002.3404, 11023.956, 11045.5716, 
                11067.1872, 11088.8028, 11110.4184, 11132.034, 11153.6496, 11175.2652, 
                11196.8808, 11218.4964, 11240.112, 11261.7276, 11283.3432, 11304.9588, 
                11326.5744, 11348.19, 11369.8056, 11391.4212, 11413.0368, 11434.6524, 
                11456.268, 11477.8836, 11499.4992, 11487.5205549999, 11475.5419099999, 
                11463.563265, 11451.5846199999, 11439.6059749999, 11427.62733, 
                11415.6486849999, 11403.6700399999, 11391.691395, 11392.093877472, 
                11411.068051152, 11430.042224832, 11449.016398512, 11467.990572192, 
                11486.964745872)

TexasNFPs4 <- c(10807.8, 10832.11755, 10856.4351, 10880.75265, 10905.0702, 
                10929.38775, 10953.7053, 10978.02285, 11002.3404, 11026.65795, 
                11050.9755, 11075.29305, 11099.6106, 11123.92815, 11148.2457, 
                11172.56325, 11196.8808, 11221.19835, 11245.5159, 11269.83345, 
                11294.151, 11318.46855, 11342.7861, 11367.10365, 11391.4212, 
                11415.73875, 11440.0563, 11464.37385, 11488.6914, 11513.00895, 
                11537.3265, 11561.64405, 11585.9616, 11610.27915, 11634.5967, 
                11658.91425, 11683.2318, 11707.54935, 11731.8669, 11756.18445, 
                11780.502, 11804.81955, 11829.1371, 11853.45465, 11877.7722, 
                11902.08975, 11926.4073, 11950.72485)

#Build the zoo objects
SATem <- zooreg(SATem1, frequency=12, start=as.yearmon("Jan 1885"))
TexasNFP <- zooreg(TexasNFP1, frequency=12, start=as.yearmon("May 2003"))
ERCOT <- zooreg(ERCOT1, frequency=12, start=as.yearmon("May 2003"))

#Forecast the data using the forecast() "black box"
for_arima <- forecast(auto.arima(log(as.ts(ERCOT))), level=0.98, h=48)
timind <- index(for_arima$fitted)
plot(timind, ERCOT, type="l", xlim=c(2003, 2017), ylim=c(10000, 60000), log="y", xlab="", xaxs="i", yaxs="i", lwd=2, col="darkgray")
mtext("ARIMA Forecast - ERCOT Monthly Average Load (MW)", side=3, line=1, font=2, cex=1.3)
lines(index(for_arima$mean), exp(for_arima$mean), lwd=2, col="blue")
lines(index(for_arima$mean), exp(for_arima$lower), lty=2, col="blue")
lines(index(for_arima$mean), exp(for_arima$upper), lty=2, col="blue")
grid()

#Build a zoo object with the important input data
datmod <- merge(ERCOT, TexasNFP, SATem, all=FALSE)

#Plot the data
plot(datmod)

#Define the model equation as a formula so it can be used in mutiple places
modequ <- formula(log(ERCOT) ~ log(TexasNFP) + pmax(0.0, SATem-TemSum) + pmin(0.0, SATem-TemWin))

#Build a function that DEoptim will monitor that fits the model with the
#assumed temperature range for the "Comfortable Temperature Zone".
minadjrsqu <- function(x){

  #Extract the upper and lower temperatures for the "Comfortable Temperature Zone"
  TemSum <- x[1]
  TemWin <- x[2]

  #Assignment is required because modequ is defined outside of the function's
  #scope and the two local variables are needed in modequ.
  assign("TemSum", TemSum, pos=.GlobalEnv)
  assign("TemWin", TemWin, pos=.GlobalEnv)

  #Build a model to estimate ERCOT loads
  mod <- lm(formula=modequ, data=datmod)

  #Because DEoptim searches for a minimum, subtract the model Adjusted R Square
  #value from 1.0 and return it.
  1.0 - summary(mod)$adj.r.squared

}

lower <- c(65, 50) #Lower search limits of "Comfortable Temperature Zone". The format is c(TemSum, TemWin)
upper <- c(80, 65) #Upper search limits of "Comfortable Temperature Zone". The format is c(TemSum, TemWin)

#Fit the model, searching for the best "Comfortable Temperature Zone".
outDEoptim <- DEoptim(minadjrsqu, lower, upper, DEoptim.control(itermax = 100))

#Using the results from DEoptim, define the "Comfortable Temperature Zone".
TemSum <- outDEoptim$optim$bestmem[1] #Temperature at top of "Comfortable Temperature Zone"
TemWin <- outDEoptim$optim$bestmem[2] #Temperature at bottom of "Comfortable Temperature Zone"
cat(paste("The comfortable Temperature Zone is between ", format(TemWin, digits=3) , " and ", format(TemSum, digits=3), " degF.", sep=""))

#Fit the model
mod <- lm(formula=modequ, data=datmod)
summary(mod)
(ressig <- summary(mod)$sigma)

#Plot the original data and the model fit
plot(datmod$ERCOT)
modfit <- zooreg(exp(predict(mod)), frequency=12, start=start(datmod))
lines(modfit, col="green")

#Combine the 4 Texas Nonfarm Payroll Scenarios into a single zoo object
s1 <- zooreg(TexasNFPs1, frequency=12, start=end(datmod) + (1/12))
s2 <- zooreg(TexasNFPs2, frequency=12, start=end(datmod) + (1/12))
s3 <- zooreg(TexasNFPs3, frequency=12, start=end(datmod) + (1/12))
s4 <- zooreg(TexasNFPs4, frequency=12, start=end(datmod) + (1/12))
TexasNFPs <- merge(s1, s2, s3, s4)
plot(TexasNFPs, las=2)

#Build a temperature matrix to speed up the temperature sampling process
numyea <- floor(length(SATem)/12)
SATemmat <- matrix(SATem[1:(12*numyea)], nrow=numyea, ncol=12, byrow=TRUE)

#Setup the forecast
alpha <- 0.98 #alpha value that is used for predict.lm(...level=alpha)
onetai <- (1 - alpha)/2 #confidence value to be used with the quantile() function.
forsam <- 1000 #Number of samples used each month in the forecast
formon <- 48 #Number of months to forecast

#Build the data structure that will be used for the forecast
fordatmod <- zooreg(matrix(NA, nrow=formon, ncol=3, dimnames=list(c(), c("ERCOT", "TexasNFP", "SATem"))), frequency=12, start=end(datmod) + (1/12))

#Build a vector with the desired month to speed up the determination
#of the month number
mon <- as.numeric(cycle(fordatmod))

#Setup the matrix to hold the predictions
premodfit <- zooreg(matrix(0, nrow=formon, ncol = forsam), frequency=12,  start=end(datmod) + (1/12))
premodlwr <- zooreg(matrix(0, nrow=formon, ncol = forsam), frequency=12,  start=end(datmod) + (1/12))
premodupr <- zooreg(matrix(0, nrow=formon, ncol = forsam), frequency=12,  start=end(datmod) + (1/12))

#Generate a prediction, with its "prediction interval" for each sample and
#store the results so the upper and lower tails can be calculated.
for (j in 1:forsam) {
  #Load the sampled scenario for the Texas Nonfarm Payroll data
  fordatmod[, 2] <- TexasNFPs[, sample(x=1:4, size=1, prob=c(0.2, 0.2, 0.4, 0.2))]

  #Load the monthly temperature forecast data into fordatmod
  for (i in 1:formon) {
    fordatmod[i, 3] <- sample(SATemmat[,mon[i]], 1)
  }

  #Forecast the ERCOT load for the desired number of months
  #premod[, j] <- exp(predict(mod, newdata=fordatmod, interval="none") + rnorm(formon, 0, sd=ressig))
  pre <- predict(mod, newdata=fordatmod, interval="prediction", level=alpha)
  premodfit[, j] <- exp(pre[, 1])
  premodlwr[, j] <- exp(pre[, 2])
  premodupr[, j] <- exp(pre[, 3])

}

#Plot the forecast for the regression model.
plot(timind, ERCOT, type="l", xlim=c(2003, 2017), ylim=c(10000, 60000), log="y", xlab="", xaxs="i", yaxs="i", lwd=2, col="darkgray")
mtext("Regression Forecast - ERCOT Monthly Average Load (MW)", side=3, line=1, font=2, cex=1.3)
lines(index(premodfit), apply(premodfit, 1, quantile, 0.50), lwd=2, col="red")
lines(index(premodfit), apply(premodlwr, 1, quantile, onetai), lty=2, col="red")
lines(index(premodfit), apply(premodupr, 1, quantile, (1 - onetai)), lty=2, col="red")
grid()

Web Host Advertising:

Advertisements
This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s