8.1. Stock Market Prediction

A recent post on Towards Data Science (TDS) demonstrated the use of ARIMA models to predict stock market data with raw statsmodels. This post addresses the same problem using pmdarima’s auto-ARIMA, and ends up achieving a different result with an even lower error rate.

8.1.1. Imports & data loading

For this example, all we’ll need is Numpy, Pandas and pmdarima. Matplotlib is optional, but highly encouraged in order to qualitatively validate the results of the model fit. To run this example, you’ll need pmdarima version 1.3.0 or greater. If you’re running this in a notebook, make sure to include %matplotlib inline, or the plots will not show up under your cells!

This example was also designed for Python 3.6+. If you’re running python 3.5, simply remove the F-string print statements and it will work.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# %matplotlib inline

import pmdarima as pm
print(f"Using pmdarima {pm.__version__}")
# Using pmdarima 1.3.0

The pmdarima module conveniently includes the dataset we’ll be using as an internal utility. Rather than carting around .csv files, you can simply load the data from the package:

from pmdarima.datasets.stocks import load_msft

df = load_msft()
df.head()

Date	Open	High	Low	Close	Volume	OpenInt
1986-03-13	0.06720	0.07533	0.06720	0.07533	1371330506	0
1986-03-14	0.07533	0.07533	0.07533	0.07533	409569463	0
1986-03-17	0.07533	0.07533	0.07533	0.07533	176995245	0
1986-03-18	0.07533	0.07533	0.07533	0.07533	90067008	0
1986-03-19	0.07533	0.07533	0.07533	0.07533	63655515	0

8.1.2. Data splitting

As with all statistical and ML modeling, we need to make sure we’ve split our data so we can evaluate model performance on a hold-out set. However, unlike other traditional supervised learning, time series models intrinsically introduce endogenous temporality, meaning that the values at any given point \(y_{t}\) in our time series likely have some effect on some future value, \(y_{t+n}\). Therefore, we cannot simply split our data randomly; we must make a clean split in our time series (and exogenous variables, if present).

As in the TDS example, we’ll use \(0.8 * dataSize\) as our training sample.

train_len = int(df.shape[0] * 0.8)
train_data, test_data = df[:train_len], df[train_len:]

y_train = train_data['Open'].values
y_test = test_data['Open'].values

print(f"{train_len} train samples")
print(f"{df.shape[0] - train_len} test samples")
# 6386 train samples
# 1597 test samples

8.1.3. Pre-modeling analysis

As you may know (if not, venture over to Tips to using auto_arima before continuing), an ARIMA model has 3 core hyper-parameters, known as “order”:

• \(p\): The order of the auto-regressive (AR) model (i.e., the number of lag observations)
• \(d\): The degree of differencing.
• \(q\): The order of the moving average (MA) model. This is essentially the size of the “window” function over your time series data.

Part of the science behind the auto-arima approach is intelligently finding the proper combination of p, d, q such that you achieve the best fit. The TDS article took the approach of fixing the \(p\) parameter at 5 after examining auto-correlations with lag plots. A lag plot can provide clues about the underlying structure of your data [1]:

• A linear shape to the plot suggests that an autoregressive model is probably a better choice.
• An elliptical plot suggests that the data comes from a single-cycle sinusoidal model.

from pandas.plotting import lag_plot

fig, axes = plt.subplots(3, 2, figsize=(8, 12))
plt.title('MSFT Autocorrelation plot')

# The axis coordinates for the plots
ax_idcs = [
    (0, 0),
    (0, 1),
    (1, 0),
    (1, 1),
    (2, 0),
    (2, 1)
]

for lag, ax_coords in enumerate(ax_idcs, 1):
    ax_row, ax_col = ax_coords
    axis = axes[ax_row][ax_col]
    lag_plot(df['Open'], lag=lag, ax=axis)
    axis.set_title(f"Lag={lag}")

plt.show()

As you can see, all the lags look fairly linear, so it’s a good indicator that an auto-regressive model is a good choice. But since we don’t want to allow simple visual bias to impact our decision here, we’ll allow the auto_arima to select the proper lag term for us.

8.1.3.1. Estimating the differencing term

The TDS article selecting \(d=1\) as the differencing term. But how did they make that choice? With pmdarima, we can run several differencing tests against the time series to select the best number of differences such that the time series will be stationary.

Here, we’ll use the KPSS test and ADF test, selecting the maximum value between the two to be conservative. Fortunately, in this case, both tests indicated that \(d=1\) was the best answer, but in the case where they disagreed, we could try both or allow auto_arima to auto-select the d term.

from pmdarima.arima import ndiffs

kpss_diffs = ndiffs(y_train, alpha=0.05, test='kpss', max_d=6)
adf_diffs = ndiffs(y_train, alpha=0.05, test='adf', max_d=6)
n_diffs = max(adf_diffs, kpss_diffs)

print(f"Estimated differencing term: {n_diffs}")
# Estimated differencing term: 1

Therefore, we will use \(d=1\).

8.1.4. Fitting our model

Now it’s time to let the auto_arima method do its magic:

auto = pm.auto_arima(y_train, d=n_diffs, seasonal=False, stepwise=True,
                     suppress_warnings=True, error_action="ignore", max_p=6,
                     max_order=None, trace=True)

Notice that we preset d=n_diffs, since we’ve already settled on a value for d. However, we’re allowing our ARIMA models explore various values of p and q. After a few seconds, we arrive at the following solution:

print(auto.order)
# (1, 1, 1)

Where the TDS model was of order (5, 1, 0), we ended up selecting a significantly more simple model. But how does it perform?

8.1.5. Updating the model

Now that the heavy lifting of selecting model hyper-parameters has been performed, we can update our model by simulating days passing with our test set. For each new observation, we’ll let our model progress for several more iterations, allowing MLE to update its discovered parameters and shifting the latest observed value. Then we can measure the error on the forecasts:

from sklearn.metrics import mean_squared_error
from pmdarima.metrics import smape

model = auto  # seeded from the model we've already fit

def forecast_one_step():
    fc, conf_int = model.predict(n_periods=1, return_conf_int=True)
    return (
        fc.tolist()[0],
        np.asarray(conf_int).tolist()[0])

forecasts = []
confidence_intervals = []

for new_ob in y_test:
    fc, conf = forecast_one_step()
    forecasts.append(fc)
    confidence_intervals.append(conf)

    # Updates the existing model with a small number of MLE steps
    model.update(new_ob)

print(f"Mean squared error: {mean_squared_error(y_test, forecasts)}")
print(f"SMAPE: {smape(y_test, forecasts)}")
# Mean squared error: 0.3416473178248818
# SMAPE: 0.981464018635346

In the end, our model ended up way out-performing the TDS model!

Source	MSE	SMAPE
pmdarima	0.342	0.981 (!!)
TDS article	0.343	40.776

8.1.6. Viewing forecasts

Let’s take a look at the forecasts our model produces overlaid on the actuals (in the first plot), and the confidence intervals of the forecasts (in the second plot):

fig, axes = plt.subplots(2, 1, figsize=(12, 12))

# --------------------- Actual vs. Predicted --------------------------
axes[0].plot(y_train, color='blue', label='Training Data')
axes[0].plot(test_data.index, forecasts, color='green', marker='o',
             label='Predicted Price')

axes[0].plot(test_data.index, y_test, color='red', label='Actual Price')
axes[0].set_title('Microsoft Prices Prediction')
axes[0].set_xlabel('Dates')
axes[0].set_ylabel('Prices')

axes[0].set_xticks(np.arange(0, 7982, 1300).tolist(), df['Date'][0:7982:1300].tolist())
axes[0].legend()


# ------------------ Predicted with confidence intervals ----------------
axes[1].plot(y_train, color='blue', label='Training Data')
axes[1].plot(test_data.index, forecasts, color='green',
             label='Predicted Price')

axes[1].set_title('Prices Predictions & Confidence Intervals')
axes[1].set_xlabel('Dates')
axes[1].set_ylabel('Prices')

conf_int = np.asarray(confidence_intervals)
axes[1].fill_between(test_data.index,
                     conf_int[:, 0], conf_int[:, 1],
                     alpha=0.9, color='orange',
                     label="Confidence Intervals")

axes[1].set_xticks(np.arange(0, 7982, 1300).tolist(), df['Date'][0:7982:1300].tolist())
axes[1].legend()

8.1.7. Conclusion

The TDS article provided an awesome example of how to use ARIMAs to predict stocks. Our hope in this example was to show how using pmdarima can simplify and enhance the models you build. If you’d like to run the already-setup notebook for yourself, head on over to the project’s Git page and grab the example notebook.