I have $M$ stacks of multivariate time series ($N$ variables sampled over $T$ time points), and each of the $M$ multivariate time series is associated with one of $K$ known classes. My goals are:

(1) being able to predict to which of the $K$ classes a given multivariate time series belongs;

(2) determine what is the subset of $N$ variables that contributes the most to each of the $K$ classes over the $M$ multivariate time series.


I have $M$ different trials, each of them storing the recorded activity of $N$ neurons ($N$>1000) while the animal is either moving one finger (class $K_1$) or another (class $K_2$). As such, each of the $M$ trials is associated with either one or the other class. My goals are:

(1) to predict which one of the two fingers the animal will move given the recorded activity of the $N$ neurons.

(2) identifying which are the subsets of neurons that are responsible for the movement of either one or the other finger.

  • $\begingroup$ Welcome to Cross Validated! This seems an interesting question, but it's hard to understand what you want, exactly. On one hand, it looks like you're interested in time series classification. You have $K$ classes, and $M$ time series, and you would like to know to which of the $K$ classes, each of the $M$ time series belong. But then you cite also $N$ time series, and this is where it becomes confusing. How many time series do you have? $N$ or $M$? Or do you mean that you have $M$ multivariate time series, i.e., each time series contains $N$ different variables? 1/ $\endgroup$
    – DeltaIV
    Commented Jan 20, 2018 at 10:56
  • $\begingroup$ 2/ Can you please edit your question by including an example of your data? If you don't clarify the question better, it may get closed. $\endgroup$
    – DeltaIV
    Commented Jan 20, 2018 at 10:57
  • 1
    $\begingroup$ Hi DeltaIV, and thanks for your comment. I am not familiar with the field of time series classification, thus the confusion. I am going to edit my question in order to make it clearer, and I'll also add an example of what I would like to achieve. Please let me know if you believe the question remains obscure. All the best $\endgroup$
    – Rugby
    Commented Jan 21, 2018 at 15:21
  • $\begingroup$ You're welcome. Can you add a data sample? 1000 variables is a lot, but maybe you could show just a few. $\endgroup$
    – DeltaIV
    Commented Jan 21, 2018 at 20:49
  • $\begingroup$ How long is each of your time series? How many time samples? $\endgroup$
    – DeltaIV
    Commented Jan 22, 2018 at 17:53

2 Answers 2


First of all, let's make sure I understood your question. You have $M$ time series (say, $M=100$ just to fix ideas), of probably different length. You don't say anything about how long these time series are, but let's assume they typically contain $\approx N_T=1000$ samples. Each time series has more than $N>1000$ features.

Suppose all your time series are stored as a tidy data frame, i.e., one feature for column, and one observation of each time series for row. You have a data frame bigger than $10^5\times 1000$. This is hefty. Assumption: you have enough memory to handle this baby, or you can afford to pay for some AWS/Google Cloud ML/whatever cloud Machine Learning instances. If you don't, my answer probably won't work for you, and you should look into Vopal Wabbit, though I have no idea if it supports time series classification.

I would try two different approaches, and use $k$-fold cross-validation (with $k=5$ or $10$) to choose among them. The problem seems complicated to me, so I would reserve a few of your $M$ time series as a test set, but you don't strictly have to do this and with $k$-fold CV you could just use all of them to choose your favorite classification procedure.

Approach 1: extract time series-related features and use a Random Forest or GBM classifier

In this approach, you try to extract features which are "typical" of a time series from your time series, using dedicated packages, such as for example tsfresh in Python.

Step 0: normalization

This is not strictly necessary, but it's definitely good practice. In practice, consider activations $X_i(t_0),\dots,X_i(t_{Nj})$ of neuron $i$. You sum together the values of neuron activation at all time instants and across all time series instances, and you get a sample mean activation $\bar{X}_i$. Similarly, you compute a sample standard deviation $S_i$. You then standardize the neuron activation as

$$X'_i(t_j) =\frac{X_i(t_j)-\bar{X}_i}{S_i}$$

Repeat for all neurons.

Step 1: feature extraction & filtering

"Typical" features may be the absolute energy, the spectral density, the maximum, minimum, median, mean, number of peaks, etc.

enter image description here

Note that in this approach each neuron will count as a separate time series. You can (you must!) include the fact that the $N$ neurons belong to the same instance (to the same multivariate time series) by using an ID (grouping) variable. For example, if you look here, the id column would contain the values $1,\dots,M=100$, while the $N$ neurons would correspond to the 6 columns $F_x, \dots, T_z$ in the example.

You now have an even bigger list of features: you may need to filter some features out to keep the training phase of your classifier manageable. The fresh (featuRe Extraction based on Scalable Hypothesis tests) algorithm should do that. It is described here:

Christ, M., Kempa-Liehr, A.W. and Feindt, M. (2016). Distributed and parallel time series feature extraction for industrial big data applications. ArXiv e-prints: 1610.07717 URL: http://adsabs.harvard.edu/abs/2016arXiv161007717C

tsfresh (as the name suggests) implements the algorithm.

Step 2: classification

Suppose your set of "relevant" features contains $f$ features $v_1,\dots,v_{f}$. You then have a training set containing $M$ samples, where each sample is $P_i=(v_{1i},\dots,v_{fi},y_i)$, with $y\in{1,\dots,K}$ being the class label of sample $P_i$. You can now train a random forest classifier or a GBM (Gradient Boosted Machine) on your training set.

Note: once trained, each time you use your classifier for a prediction, you'll need to extract the features $v_1,\dots\,v_f$ first. This means that if you want to classify a new unseen time series, which doesn't belong to the training set but of course has to come from the same probability distribution, first of all you compute the features $v_1,\dots\,v_f$ starting from the time series of your $N$ neurons.

Approach 2: Dynamic Time Warping + 1-NN

Another approach, loosely related to the first one, is to define some metric in the time series space which allows us to quantify "how close" two time series are. However, in this approach, which is commonly used for time series classification, you won't be able to say which neuron contributes the most to classification, which goes against one of your requirements.

Step 0: normalization

Same as before.

Step 1: Dynamic Time Warping

Dynamic Time Warping finds an optimal match between two sequences by allowing a non-linear mapping of one sequence to another, and minimizing the distance between two sequences. The sequences are "warped" non-linearly to determine their similarity independent of any nonlinear variations in the time dimension.

Suppose $$L_k=\{\mathbf{X}_k(t_{0})=(X_{k1}(t_{0}),\dots,X_{kN}(t_{0})),\dots,\mathbf{X}_k(t_{m})=(X_{k1}(t_{m}),\dots,X_{kN}(t_{m})\}$$ is the $N-$variate time series $k$, from your training set of $M$ time series. Each time sample is a vector of $N$ scalars, the $N$ neurons activation at time $t$. Let's denote $\mathbf{X}_k(t_{i})=\mathbf{a}_i$ for simplicty. Now, let


denote another time series, and let's denote its elements $\mathbf{X}_h(t_{j})=\mathbf{b}_j$. We thus have two sequences of vectors

$$\mathbf{a}_1,\dots,\mathbf{a}_m \quad \text{and} \quad \mathbf{b}_1,\dots,\mathbf{b}_n $$

having possibly different lenghts. For DTW, we build an $m\times n$ matrix $M$(called path matrix) where $m_{ij}=||\mathbf{a}_i-\mathbf{b}_j||$. We want to find a path $W=w_1,\dots,w_k$ through $M$, which minimizes the sum of the distances $m_{ij}$, and which is subject to various constraints such as for example endpoint matching, continuity and monotonicity (see here for details).

Computing the DTW in the most naive implementation requires $O(nm)$ operations, which is considerable in your case. There are numerous tricks you can employ to speed-up the computations: locality constraint or windowing approximates the exact DTW computation, while the Keogh lower bound is a rigorous lower bound or DTW, whose computation has only linear complexity. This is discussed in detail here, where you can also find a Python implementation of DTW.

Step 2: classify based on 1-NN

At this point, classification is pretty simple. As mentioned in the linked references, DTW + 1-NN seems to do pretty well for time series classification. In practice this means that each time you get a new time series $L$, you compute the DTW between $L$ and $L_1,\dots,L_M$ (your training set), you find the time series $L_i$ such that its DTW distance from $L$ is minimum, and you assign $L$ the same class as $L_i$. Note that using the Keogh lower bound, you usually don't need to compute all $M$ DTW distances. Suppose your DTW distance between $L_1$ and $L$ is $d$: if you then compute the Keogh l.b. for $L_2$ and you get a value $k>d$, you can safely skip the computation of the distance between $L$ and $L_2$, and go to $L_3$.

  • $\begingroup$ Excellent ! Could you have a look at a similar question here ? stats.stackexchange.com/questions/368100/… $\endgroup$ Commented Sep 22, 2018 at 18:04
  • $\begingroup$ Thank you for that @DeltaIV! Is it possible to apply DTW between two 2d matrices (in his case, each time series has 1000 time points x 1000 features)? $\endgroup$
    – noobalert
    Commented Feb 14, 2019 at 1:13

[This is an edited response, based on the updated project description and feedback from other people]

Main point: as stated earlier, the project description is a classification problem. We work with a collection of independent N-variate time series. Each time series to is to be classified into one of K classes. The nature of the data is

1) irregular,

2) multidimensional,

3) with potentially different importance between extreme movements and regular diffusive movements,

4) not amenable to good old continuous and "regular" basis functions,

5) etc, etc, etc.

For that reason, random forest has been suggested as a classification tool. The suggestion has been scrutinized and explored in detail in a well-thought-through response by DeltaIV. Now....... when a machine learning method is proposed in a complex setting, almost never does it mean that the data must be fed to the algorithm in the raw form. The data must be preprocessed and what is done to the data depends on the choice of the machine learning method. In fact, there was a recent thread raising a question: are the decisions made during the data preporcessing stage more important than the choice of the machine learning algorithm?

The answer is: both are important. Of course, how can it be any different?... To complement the well-written overview of DeltaIV, let me talk about one particular data preprocessing approach popular in world of systematic trading. Every N-variate time series is compressed into a set of features before being submitted to a random forest. The features can be

1) mean reversion as quantified by ADF statistics,

2) momentum as quantified by autocorrelation in the log-returns,

3) codependence of extreme values of log-returns (tail dependence),

4) n4-dimensional summary of time trends,

5) n5-dimensional summary of GARCH-implied volatilies,

6) Spearman's rho or Kendall's tau calculated off the N components of the time series,

7) etc, etc, etc.

The choice of features is specific to the objective of your project. Good luck!

  • $\begingroup$ First of all, many thanks for your answer. What you suggested may indeed be useful in a later phase, i.e. when I have tried out different statistical methods to try answer my question. I am going to edit the original question so to make it clearer! Please don't hesitate to let me know if you may have further doubts. $\endgroup$
    – Rugby
    Commented Jan 18, 2018 at 9:55
  • $\begingroup$ Thank you. Your edit is just a classification problem. You do not want to choose just a subset of features. You use all the features to classify the current state to "moving finger 1" or "moving finger 2". Random forests are one classification method. $\endgroup$
    – stans
    Commented Jan 18, 2018 at 11:50
  • $\begingroup$ Yes, it is indeed a classification problem. I am not familiar with random forests: do you know of any application of them on time series data? thanks again! $\endgroup$
    – Rugby
    Commented Jan 18, 2018 at 14:34
  • $\begingroup$ They are applied everywhere, across industries and data types. Random forests are just one option. Any machine learning (data mining) book will introduce you to a number of classification methods. However, in your case, random forest is likely to be superior to LDA or logistic regression because it will be able to identify some (not all) patterns in the time series. $\endgroup$
    – stans
    Commented Jan 18, 2018 at 14:38
  • 1
    $\begingroup$ @DeltaIV, I guess we have had slight miscommunication due to my response being unnecessarily schematic. Each N-variate time series may exhibit serial correlation but the created features are assumed to be independent over many different time series objects which are used to train the classification engine (e.g. random forest). That is what I meant... I've noticed you have posted pretty much the same approach as "approach 1". That is the approach that I referred to as popular in the world of finance. $\endgroup$
    – stans
    Commented Jan 24, 2018 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.