spooky: a brief introduction

Expanding time features using Spectral analysis and Jack-Knife resampling

spooky uses Discrete Fast Fourier Transformation to extrapolate time features beyond their boundaries. As you can see below, moving inside out the picture, the main points of spooky are the following:

1. We can decompose a time feature signal in all its components using a Discrete Fast-Fourier Transformation: some components will be highly related to future trends and will be relevant for any purpose of future generalizations, while others will probably be no more than noises and errors without any representative power for extrapolation: we sift through all the possible frequencies using a jack-knife resampling to magnify the predicting power of best frequencies. Each time feature is expanded inside each validation window: after the spectral analysis, the signal of each time feature is extended into the future for the requested sequence of time steps (seq_len), then the expansion is being carried out on the resampled set from a jack-knife leaving N points out (lno).

2. Some basic transformation are directly managed in background. Differentiation and integration are automatically managed by spooky using the maximal p-value in a recursive F-test for de-trending each time-feature: this allows to easily determine the different dynamic characteristics of each time feature (random walk, trend, exponential). Maybe not cool like Augmented Dickey-Fuller or Ljung Box Test, but it works for our humble purpose. If you have limited missing values in your time features, spooky automatically proceeds with the imputation using the Kalman filter method¹. If you prefer to project into the future the smoothed version, you can set smoother = TRUE to use loess² function.

3. The test errors are cross-validated through an expanding validation n_windows: the default value is set to 10, meaning that the time features are divided into 10 + 1 segments guaranteeing at least ten validation sets to measure the error on unforeseen data. For each point in the prediction sequence, a thousand samples are collected for the calculation of quantiles, mean, mode, standard deviation, skewness and kurtosis, and other less common measures that are provided for each time step (see below).

4. Concluding this brief excursus inside out the mode, beyond the windowing phase, spooky manages the random search sampling n_samp models, randomly selected within the cartesian space delimited by seq_len and lno. The process is really fast.

   Easy, that’s all.

The process flow of spooky

A look at the price dynamics of tech giants’ stocks

The dataset time features included with spooky is a recent take on IBM and Microsoft stock prices (source: Yahoo Finance). The data is expected in a dataframe format, where each column represents a different time series (the date information is not mandatory and could be provided separately).

In the first example, we are predicting the close price for IBM and Microsoft. In this example we try to set seq_len= 100 (sequence length) and lno= 1 (leave-N-out), using a cross-validation scheme of 10 n_windows for error measurement.

example1 <- spooky(time_features, n_samp = 1, seq_len = 50, n_windows = 10,  dates = as.Date(rownames(time_features)), lno = 1)
  time: 10.11 sec elapsed

The result is a list of different components, as you can see below.

names(example1)
  [1] "exploration" "history"     "best_model"  "time_log"

The first variable, exploration, includes all the models evaluated during the random search within the boundaries decided by users: besides the predictions for each feature, exploration includes the test error statistics for each time sequence (length equals to seq_len) and model information. The second variable, history, summarizes the hyper-parameters selected by users, the relative max error metrics across different time features (me, mae, rmsse, mpe, mape, rmae, rrmse, rame, mase, smse, sce, gmrae ³) and the weighted average rank used to sort out the best model. The last variable, best_model, collects testing_errors, preds, plots for the selected model.

The prediction is a list including the predicted results for each time-feature (quantile, min, max, mean, mode, sd, skewness, kurtosis, iqr to range, risk ratio, upside probability, divergence for each time point in the seq_len sequence). The IQR to range is the interquartile range to the min-max range, the risk ratio is the range above median to the range below it, the upside probability is the probability of growth compared to the former point in the time sequence, the divergence is the maximum distance of cumulative normal curve of each point to the former point in the sequence.

For each time features included in the model, you get a plot of the median with the chosen confidence interval (ci default is 0.8). As in other packages⁴, we provide different stats to give a better hint on the different dynamics related to aleatoric and epistemic uncertainty.

  $IBM.Close

  
  $MSFT.Close

Fast model optimization in a compact space

Now, let’s try a fast random search to optimize our estimations. The following example show how to sample 30 different models for a sequence of 50 time steps. Setting lno = NULL starts a search between 1 and the square root of time feature length.

example2 <- spooky(time_features, n_samp = 30, seq_len = 50, n_windows = 10,  dates = as.Date(rownames(time_features)), lno = NULL)
  time: 489.36 sec elapsed

If we compare the error statistics from the best model in example2 with the model in example1, we see an interesting improvement in the testing errors for IBM and Microsoft (the default error measurement for relative and scaled error is naive). The error statistics from example1, evaluated as the mean on 10 validation windows, are the following:

example1$best_model$testing_errors
                me   mae     mse rmsse   mpe  mape  rmae rrmse  rame  mase   smse
  IBM.Close  4.283 8.575 122.207 7.580 0.032 0.067 0.068 0.080 0.062 4.793 67.670
  MSFT.Close 4.911 7.782 132.215 6.302 0.021 0.042 0.042 0.048 0.039 4.330 57.677
                 sce gmrae
  IBM.Close  117.787 0.048
  MSFT.Close 124.995 0.034

The error statistics from example2.

example2$best_model$testing_errors
                me   mae     mse rmsse   mpe  mape  rmae rrmse  rame  mase   smse
  IBM.Close  3.478 7.763 107.472 6.804 0.027 0.061 0.062 0.073 0.051 4.266 58.027
  MSFT.Close 4.229 7.439 114.796 6.092 0.018 0.041 0.041 0.046 0.037 4.189 51.030
                 sce gmrae
  IBM.Close   92.959 0.046
  MSFT.Close 110.937 0.033

The improvement is clear for both the time features, but we are still using a naive approach to measure scaled and relative errors. Let’s try to shift to deviation as scale, and average as benchmark, that are more challenging evaluation criteria.

example3 <- spooky(time_features, n_samp = 30, seq_len = 50, n_windows = 10,  dates = as.Date(rownames(time_features)), lno = NULL, error_scale = "deviation", error_benchmark = "average")
  time: 736.65 sec elapsed

As you can see, the relative and scaled measures change sensibly as we raise the bar of our expectations

example3$best_model$testing_errors
                me   mae     mse rmsse   mpe  mape  rmae rrmse rame  mase   smse
  IBM.Close  3.478 7.763 107.472  6.71 0.027 0.061 1.141 1.180    1 4.151 56.380
  MSFT.Close 4.229 7.439 114.796  6.02 0.018 0.041 0.837 0.843    1 4.143 49.313
                 sce gmrae
  IBM.Close   88.589 1.099
  MSFT.Close 109.244 0.840

Prospecting the forecasting horizon

You can also explore for the best possible forecasting horizon, choosing a range for seq_len besides lno. Let’s see the best foreseeability using the ranges lno = c(1, 50) and seq_len = c(10, 20).

example4 <- spooky(time_features, n_samp = 10, n_windows = 10,  dates = as.Date(rownames(time_features)), lno = c(1, 50), seq_len = c(10, 20))
  time: 117.86 sec elapsed

As we can see from history, 13 days and 47 left-out are the best parameters to reduce the max joint error on the time features.

Let’s see the plots.

  $IBM.Close

  
  $MSFT.Close

Some useful references