dymo: a brief introduction

Predicting non-linear systems using Dynamic Mode Decomposition

dymo uses Dynamic Mode Decomposition¹ (DMD) to predict multiple time features correlated by non-linear dynamics. dymo is a pretty simple implementation of SVD-based DMD with a couple of add-ons to make it a little more interesting:

1. An introduction to DMD is out-of-scope of this vignette. What you have to know, simplifying a bit, is that the DMD linearizes the non-linear dynamics of multiple time features by projecting them in an Hilbert space to approximates their value. What DMD does is understanding the non-linear dynamics of a system in time and approximates them as linear multiplicative system (you can think of it as a kind of regression, but no ordinary least squares here). Key points: you need at least two different time features (a single one is not a system by definition) and you need to choose the rank of the SVD approximation (dymo tests all the different ranks and selects the best one). You can also explore the space of all possible combinations of time features by selecting the minimum and the maximun number of features (from 2 to the total number of features).

2. Some basic transformation are directly managed in background. Differentiation and integration are automatically managed by dymo 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 (somehow more simple and practical compared to other formal approaches like Augmented Dickey-Fuller or Ljung Box Test). If you have limited missing values in your time features, dymo 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).

The process flow of dymo

A look at the price dynamics of tech giants’ stocks

The dataset time features included with dymo is a recent take on some Big Techs’ 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), using a cross-validation scheme of 10 n_windows for error measurement.

example1 <- dymo(time_features[, 2:6], seq_len = 100, n_windows = 10,  dates = time_features$date)
  time: 13.96 sec elapsed

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

names(example1)
  [1] "comb_metrics" "best_model"   "time_log"
names(example1$best_model)
  [1] "best_combination" "quant_preds"      "testing_errors"   "plots"

comb_metrics includes error metrics for all possible combinations of time features (beside prediction score, me, mae, mse, rmsse, mpe, mape, rmae, rrmse, rame, mase, smse, sce, gmrae ⁴, averaged across features, ranks and validation windows). best_model collects the best combination of features and rank (best_combination), the error metric for each time feature (testing_errors), the prediction intervals (quant_preds) and the visualizations (plots).

The quant_preds 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

  
  $AAPL.Close

  
  $AMZN.Close

  
  $GOOGL.Close

  
  $MSFT.Close

The best combination is not always what you may think

Now, let’s try a fast grid search for the best combination setting the minimum and the max number of time features to consider (obviously, the number can’t be lower than two and greater than the total number of features). The following example shows how to sample 75 different models at different rank values for a sequence of 100 time steps, setting min_feats and max_feats between 2 and 5.

example2 <- dymo(time_features[, 2:6], seq_len = 100, n_windows = 10,  min_feats = 2, max_feats = 5, dates = time_features$date)
  time: 146.4 sec elapsed

The best model includes time features number 1 and number 3 (respectively, IBM and Amazon). If we compare the error statistics from the best model in example2 with the model in example1, for IBM and Amazon we see consistent improvement. All the relative and scaled error metrics defaults to naive, but you can choose more challenging thresholds (like the deviation of the whole time feature or the average of the whole predicted sequence).

The error statistics from example1 (averaged across 10 expanding validation windows):

example1$best_model$testing_errors
              pred_scores       me      mae        mse   rmsse    mpe   mape
  IBM.Close     0.3960911   2.3043   7.5842    96.5223  8.4726 0.0155 0.0590
  AAPL.Close    0.2657489   7.8821   8.7555   158.3111 13.2152 0.0841 0.0992
  AMZN.Close    0.3512089 122.8044 185.2975 74715.1934 55.4967 0.0581 0.0820
  GOOGL.Close   0.2576022  69.0797 113.7313 25695.1495 40.1203 0.0433 0.0696
  MSFT.Close    0.2363667   9.4318  13.1982   300.0342 14.6394 0.0596 0.0763
                rmae  rrmse   rame    mase      smse       sce  gmrae
  IBM.Close   1.0135 1.0083 0.9825  6.3671   80.9444  162.0369 1.0220
  AAPL.Close  1.0079 1.0064 0.9982 15.0325  239.2960 1325.7503 1.0117
  AMZN.Close  1.0062 1.0062 1.0237 11.8917 4389.4121  926.2679 0.9998
  GOOGL.Close 1.0049 1.0037 0.9642  9.6565 1938.6584  614.0788 1.0023
  MSFT.Close  1.0113 1.0108 1.0154 12.7054  235.8639 1017.8694 1.0166

The error statistics from example2 (as above, averaged across 10 expanding validation windows):

example2$best_model$testing_errors
             pred_scores       me      mae        mse   rmsse    mpe   mape
  IBM.Close    0.4018044   2.1896   7.5408    96.1175  8.4489 0.0147 0.0585
  AMZN.Close   0.3514911 120.9444 184.5771 74493.9610 55.2893 0.0576 0.0818
               rmae  rrmse   rame    mase      smse      sce  gmrae
  IBM.Close  1.0076 1.0052 0.9964  6.3319   80.6405 153.6033 1.0118
  AMZN.Close 0.9997 0.9995 1.0022 11.8531 4375.5401 916.9243 1.0012

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 <- dymo(time_features[, 2:6], seq_len = 100, n_windows = 10,  min_feats = 2, max_feats = 5, dates = time_features$date, error_scale = "deviation", error_benchmark = "average")
  time: 133.97 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
             pred_scores       me      mae        mse   rmsse    mpe   mape rmae
  IBM.Close    0.4018044   2.1896   7.5408    96.1175  2.5574 0.0147 0.0585    1
  AMZN.Close   0.3514911 120.9444 184.5771 74493.9610 12.6031 0.0576 0.0818    1
             rrmse rame   mase     smse     sce  gmrae
  IBM.Close      1    1 0.5788   7.3931 13.3151 0.9995
  AMZN.Close     1    1 0.6886 242.9613 59.2111 1.0000

Some useful references