@ -25,36 +25,332 @@ The [moving-average](https://en.wikipedia.org/wiki/Moving-average_model) aspect
Bottom line: ARIMA is used to make a model fit the special form of time series data as closely as possible.
### Preparation
Open the `/working` folder in this lesson and find the `notebook.ipynb` file. We have already loaded
Open the `/working` folder in this lesson and find the `notebook.ipynb` file. Run the notebook to load the `statsmodels` Python library; you will need this for ARIMA models.
---
## Load necessary libraries
[Step through content in blocks]
Now, load up several more libraries useful for plotting data:
## [Topic 1]
```python
import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt
import math
### Task:
from pandas.plotting import autocorrelation_plot
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.preprocessing import MinMaxScaler
from common.utils import load_data, mape
from IPython.display import Image
Work together to progressively enhance your codebase to build the project with shared code:
Now your data is loaded, so you can separate it into train and test sets. You'll train your model on the train set. As usual, after the model has finished training, you'll evaluate its accuracy using the test set. You need to ensure that the test set covers a later period in time from the training set to ensure that the model does not gain information from future time periods.
## [Topic 3]
Allocate a two-month period from September 1 to October 31, 2014 to the training set. The test set will include the two-month period of November 1 to December 31, 2014.
## 🚀Challenge
```python
train_start_dt = '2014-11-01 00:00:00'
test_start_dt = '2014-12-30 00:00:00'
```
Since this data reflects the daily consumption of energy, there is a strong seasonal pattern, but the consumption is most similar to the consumption in more recent days. You can visualize the differences:
Now that you have calibrated the scaled data, you can scale the test data:
```python
test['load'] = scaler.transform(test)
test.head()
```
## Implement ARIMA
It's time to implement ARIMA! You'll now use the `statsmodels` library that you installed earlier.
Now you need to follow several steps
1. Define the model by calling `SARIMAX()` and passing in the model parameters: p, d, and q parameters, and P, D, and Q parameters.
1. The model is prepared on the training data by calling the fit() function.
2. Predictions can be made by calling the `forecast()` function and specifying the number of steps (the `horizon`) to forecast
> 🎓 What are all these parameters for? In an ARIMA model there are 3 parameters that are used to help model the major aspects of a time series: seasonality, trend, and noise. These parameters are:
`p`: the parameter associated with the auto-regressive aspect of the model, which incorporates *past* values.
`d`: the parameter associated with the integrated part of the model, which affects the amount of *differencing* (🎓 remember differencing 👆?) to apply to a time series.
`q`: the parameter associated with the moving-average part of the model.
> Note: If your data has a seasonal aspect - which this one does - , we use a seasonal ARIMA model (SARIMA). In that case you need to use another set of parameters: `P`, `D`, and `Q` which describe the same associations as `p`, `d`, and `q`, but correspond to the seasonal components of the model.
Start by setting your preferred horizon value. Let's try 3 hours:
```python
# Specify the number of steps to forecast ahead
HORIZON = 3
print('Forecasting horizon:', HORIZON, 'hours')
```
Selecting the best values for an ARIMA model's parameters can be challenging as it's somewhat subjective and time intensive. You might consider using an `auto_arima()` function from the [`pyramid` library](https://alkaline-ml.com/pmdarima/0.9.0/modules/generated/pyramid.arima.auto_arima.html), but for now try some manual selections to find a good model.
```python
order = (4, 1, 0)
seasonal_order = (1, 1, 0, 24)
model = SARIMAX(endog=train, order=order, seasonal_order=seasonal_order)
results = model.fit()
print(results.summary())
```
TODO: Explain these results and show residuals
Add a challenge for students to work on collaboratively in class to enhance the project
You've built your first model! Now we need to find a way to evaluate it.
Optional: add a screenshot of the completed lesson's UI if appropriate
## Evaluate your model
To evaluate your model, you can perform the so-called `walk forward` validation. In practice, time series models are re-trained each time a new data becomes available. This allows the model to make the best forecast at each time step.
Starting at the beginning of the time series using this technique, train the model on the train data set. Then make a prediction on the next time step. The prediction is evaluated against the known value. The training set is then expanded to include the known value and the process is repeated.
> Note: You should keep the training set window fixed for more efficient training so that every time you add a new observation to the training set, you remove the observation from the beginning of the set.
This process provides a more robust estimation of how the model will perform in practice. However, it comes at the computation cost of creating so many models. This is acceptable if the data is small or if the model is simple, but could be an issue at scale.
Walk-forward validation is the gold standard of time series model evaluation and is recommended for your own projects.
First, create a test data point for each HORIZON step.
Observe the hourly data's prediction, compared to the actual load. How accurate is this?
Check the accuracy of your model by testing its mean absolute percentage error (MAPE) over all the predictions.
> **🧮 Show me the math**
>
> ![MAPE](images/mape.png)
>
> [MAPE](https://www.linkedin.com/pulse/what-mape-mad-msd-time-series-allameh-statistics/) is used to show prediction accuracy as a ratio defined by the above formula. The difference between Actual<sub>t</sub> and Predicted<sub>t</sub> is divided by the Actual<sub>t</sub>. "The absolute value in this calculation is summed for every forecasted point in time and divided by the number of fitted points n." [wikipedia](https://en.wikipedia.org/wiki/Mean_absolute_percentage_error)
"<sup>1</sup>Tao Hong, Pierre Pinson, Shu Fan, Hamidreza Zareipour, Alberto Troccoli and Rob J. Hyndman, \"Probabilistic energy forecasting: Global Energy Forecasting Competition 2014 and beyond\", International Journal of Forecasting, vol.32, no.3, pp 896-913, July-September, 2016."