7. Evaluating models: persistence and residuals

Level 8: Time Series - Stationarity

The goal of time series forecasting is to make accurate predictions about the future. The fast and powerful methods we usually rely on in machine learning, such as random train-test splits or k-fold cross-validation, do not work with temporal data, because they ignore the temporal component inherent to the problem. Training on 2022 data to predict 2018 is indefensible.

In time series forecasting, this evaluation of models on historical data is called backtesting. There are three main methods: a train-test split that respects the temporal ordering; multiple successive splits that respect the temporal ordering; and walk-forward validation, in which the model can be updated every time a new observation arrives.

X = Re
train_size = int(len(X) * 0.66)
train, test = X[0:train_size], X[train_size:len(X)]

print('Observations: %d' % (len(X)))
print('Training Observations: %d' % (len(train)))
print('Testing Observations: %d' % (len(test)))

The persistence model. Establishing a benchmark is essential in any time series forecasting problem, because it tells you how well every other model is really performing. Here we develop a persistence forecast to use as a baseline. The steps are: turn the univariate dataset into a supervised problem, set up the training and test sets, define the persistence model, make the forecast and establish the baseline performance, and review the complete example with its plot.

# lagged dataset
values = pd.DataFrame(Re.values)
Re2 = pd.concat([values.shift(1), values], axis=1)
Re2.columns = ['t', 't+1']
Re2 = Re2[1:]

X = Re2.values
train_size = int(len(X) * 0.66)
train, test = X[1:train_size], X[train_size:]
train_X, train_y = train[:, 0], train[:, 1]
test_X, test_y = test[:, 0], test[:, 1]

def model_persistence(x):
    return x

predictions = []
for x in test_X:
    yhat = model_persistence(x)
    predictions.append(yhat)

rmse = np.sqrt(mean_squared_error(test_y, predictions))
print('Test RMSE: %.3f' % rmse)

The model literally returns what it receives: the forecast for tomorrow is today's value. Any sophisticated model that fails to beat this benchmark is worthless.

Residual forecast errors. Forecast errors in time series are called residual errors or residuals. A residual is computed as the expected outcome minus the forecast. Exploring the residuals carefully can tell you a great deal about your model and even suggest improvements. The first plot consists of looking at the residuals over time as a line plot: we would expect it to be random around the value 0, with no trend and no cyclical structure.

residuals = [test_y[i] - predictions[i] for i in range(len(predictions))]
residuals = pd.DataFrame(residuals)

residuals.plot()
residuals.describe()
residuals.hist()
residuals.plot(kind='kde')
qqplot(residuals, line='r'); plt.show()

We can compute descriptive statistics on the residuals. What matters most is their mean: a value close to zero suggests the forecasts are unbiased, whereas positive or negative values indicate bias in one direction or the other. Knowing that bias is useful because it can be corrected directly before using or evaluating the forecasts.

Autocorrelation of the residuals. We would not expect any correlation between the residuals, which would show up as autocorrelations below the significance threshold (the dashed horizontal lines on the plot). Significant autocorrelation in the residual plot suggests that the model could do a better job by incorporating the relationship between the observations and their lagged observations, that is, autoregression.

autocorrelation_plot(residuals)
plot_acf(residuals, lags=20); plt.show()

Random walk. A random walk is a model in which the value of the variable goes up or down with equal probability at each step: the expected value of the current observation equals the previous observation. Put in code, today's value is yesterday's plus a white noise term. If the series we are modeling is a random walk, it is unpredictable, because it is literally a walk at random. Its importance lies in the fact that it is not stationary: although the mean of the random walk remains zero, the covariance does depend on time. Specifically, the covariance equals the variance multiplied by time, so as time advances so does the variance.

np.random.seed(1)
n_samples = 1000
x = w = np.random.normal(size=n_samples)
for t in range(n_samples):
    x[t] = x[t-1] + w[t]

_ = tsplot(x, lags=30)          # clearly non-stationary
_ = tsplot(np.diff(x), lags=30) # the first difference should be white noise

If the random walk model is a good fit, today's value minus yesterday's must be exactly white noise: the first differences of a random walk are a white noise process. Applying it to real financial data:

_ = tsplot(np.diff(df), lags=30)   # first difference of the S&P 500

It looks fairly close to white noise, but some significant serial correlation appears in the ACF and PACF plots. In addition, the shape of the QQ and probability plots indicates that the process is close to normal but with fat tails. That means better models should exist to explain the price of the S&P 500.