# Time Series Forecasting in Python and R

A friend recently made a prediction about the price of oil for the next three months. I thought I would perform some time series forecasting on the West Texas Intermediate prices and see if his numbers were reasonable from a dumb-numbers canned-forecasting perspective. I’m not making the claim that one can reasonably and accurately forecast oil prices with traditional time series techniques. (That’s bogus.) I’m simply doing this to learn more about forecasting.

Monthly petroleum prices can be found at the Energy Information Administration. Ever relevant, Wikipedia has a great write-up on recent trends in oil prices. Also, there is this Times article on the spike and drop in 2008 which had this apt summary,

[Oil prices are] the product of an extremely volatile mixture of speculation, oil production, weather, government policies, the global economy, the number of miles the average American is driving in any given week and so on. But the daily price is actually set — or discovered, in economic parlance — on the futures exchange.

# Data

If you go to the EIA site, you can download some monthly and weekly prices in an Excel spreadsheet. The WTI and Brent oil prices are on the second sheet of the Excel file. We can use Python to extract this data from Excel and put it into a Pandas `DataFrame` for easier processing.

```import pandas, seaborn

# create an Excel file object

# parse the first sheet
df = excel.parse( excel.sheet_names[1] )

# rename the columns
df = df.rename( columns=dict( zip( df.columns, ['Date','WTI','Brent'] ) ) )

# cut off the first 18 rows because these rows
# contain NaN values for the Brent prices
df = df[18:]

# index the data set by date
df.index = df['Date']

# remove the date column after re-indexing
df = df[['WTI','Brent']]
```

# Visualization

We can then visualize the data with `seaborn`. So, the takeaway form the big picture here is that oil prices have gone totally bonkers over the last ten years.

```df[['WTI','Brent']][::4].plot()
title('Crude Oil Prices')
xlabel('Year')
ylabel('Price [USD]')
savefig(&quot;wti_and_brent_all_data.png&quot;,dpi=200)
```

If we focus on the last three years, we get the following picture. The interesting thing about this perspective, apart from the recent marked drop, is that the WTI and Brent prices seem to start getting closer in value over time.

We can drop the last three years of weekly data, the smoothed data, and the de-trended data with the following line,

```df[-36*4:].to_csv('prices_last_three_years.csv')
```

# Seasonality

Rob Hyndman is the creator of the R forecast package, and you should really check out his blog. I found this post on R-bloggers where Hyndman proposed an interesting technique for detecting seasonality. Essentially, you fit one model with a seasonal component, and another model without a seasonal component, and then you compare their log-likelihood and their degrees of freedom using a chi-squared distribution.

For this data, I got errors attempting the `decompose()` and seasonal `HoltWinters()` functions, suggesting that there is not a strong seasonal component.

```data <- read.csv('prices_last_three_years.csv')
decompose( data\$WTI )
```
```Error in decompose(data\$WTI) : time series has no or less than 2 periods
```
```HoltWinters( data\$WTI )
```
```Error in decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) :
time series has no or less than 2 periods
```

# Holt-Winters Forecasting

I used this resource for learning about Holt-Winters forecasting. We can fit a model and make predictions using the Holt-Winters method with the following lines of code:

```# load Hyndman's forecast package
library(forecast)

# create a model with a trend
WTI.hwm <- HoltWinters( data\$WTI, gamma=FALSE )

# predict the values for the next 4 weeks #
WTI.hwf <- forecast.HoltWinters( WTI.hwm, h=4 )

# print a summary of the forecast
summary( WTI.hwf )
```
```Forecast method: HoltWinters

Model Information:
Holt-Winters exponential smoothing with trend and without seasonal component.

Call:
HoltWinters(x = data\$WTI, gamma = FALSE)

Smoothing parameters:
alpha: 1
beta : 0.3259544
gamma: FALSE

Coefficients:
[,1]
a 76.500000
b -2.011559

Error measures:
ME     RMSE      MAE        MPE     MAPE     MASE     ACF1
Training set -0.1612072 2.271556 1.785524 -0.1085958 1.855251 1.068586 0.104531

Forecasts:
Pt Forecast Lo 80 Hi 80 Lo 95 Hi 95
145       74.48 71.57 77.40 70.03 78.94
146       72.47 67.63 77.31 65.07 79.87
147       70.46 63.63 77.29 60.02 80.90
148       68.45 59.51 77.38 54.79 82.11
```

For reference, I’ll replace the forecasts with dates. You’ll note that today is 11/23/2014, and that the first predicted data point is for 11/21/2014, but there’s a lag between production and reporting, so the 11/21/2014 numbers won’t be published for a few days yet.

```        Pt Forecast Lo 80 Hi 80 Lo 95 Hi 95
11/21/2014    74.48 71.57 77.40 70.03 78.94
11/28/2014    72.47 67.63 77.31 65.07 79.87
12/05/2014    70.46 63.63 77.29 60.02 80.90
12/12/2014    68.45 59.51 77.38 54.79 82.11
```

We can then plot the forecast as follows:

```plot.forecast( WTI.hwf )
```

# ARIMA Forecasting

I prefer the outlook of the Holt-Winters forecast, but I decided to do an ARIMA(1,0,0) and ARIMA(2,0,0) forecast also, for completeness. I’m still trying to understand how to interpret the ACF and PACF plots of time series. The canonical data sets make a ton of sense, but things get murkier with more typical data sets. You end up scratching your head and saying, “Hell, let’s just try an AR(1) process.” (This is where the `forecast::auto.arima()` function comes in handy.)

```WTI.arima <- arima( data\$WTI, order=c(1,0,0) )
WTI.arimaf <- forecast.Arima( WTI.arima, h=4 )
summary( WTI.arimaf )
```
```Forecast method: ARIMA(1,0,0) with non-zero mean

Model Information:
Series: data\$WTI
ARIMA(1,0,0) with non-zero mean

Coefficients:
ar1  intercept
0.9700    94.0842
s.e.  0.0218     5.1415

sigma^2 estimated as 4.727:  log likelihood=-317.58
AIC=641.16   AICc=641.34   BIC=650.07

Error measures:
ME     RMSE      MAE        MPE     MAPE     MASE     ACF1
Training set -0.09167768 2.174233 1.676698 -0.1544393 1.771276 0.997412 0.298963

Forecasts:
Pt Forecast Lo 80 Hi 80 Lo 95 Hi 95
11/21       77.02 74.24 79.81 72.76 81.28
11/28       77.54 73.65 81.42 71.60 83.47
12/05       78.03 73.35 82.72 70.87 85.20
12/12       78.52 73.19 83.84 70.36 86.67
```

And for the AR(2) process, we have the following:

```Forecast method: ARIMA(2,0,0) with non-zero mean

Model Information:
Series: data\$WTI
ARIMA(2,0,0) with non-zero mean

Coefficients:
ar1      ar2  intercept
1.2725  -0.3220    94.9581
s.e.  0.0806   0.0829     3.2659

sigma^2 estimated as 4.278:  log likelihood=-310.4
AIC=628.79   AICc=629.08   BIC=640.67

Error measures:
ME     RMSE      MAE         MPE     MAPE      MASE
Training set -0.03991496 2.068382 1.599582 -0.09327853 1.687477 0.9515379
ACF1
Training set -0.01486163

Forecasts:
Pt Forecast Lo 80 Hi 80 Lo 95 Hi 95
11/21       76.85 74.20 79.50 72.79 80.90
11/28       77.86 73.57 82.15 71.30 84.42
12/05       79.03 73.53 84.53 70.62 87.44
12/12       80.20 73.79 86.60 70.40 89.99
```

# Conclusions

I will revisit this post in two weeks with some more details on forecasting, and an update on these forecasts. For future reference, the `ets()` and `auto.arima()` functions from the `forecast` package are a great starting point for forecasting, they will automatically fit a model and return model performance metrics.