Classical Decomposition of a time series

A time series is a data set that has a time component. Yes, it is just what you think about, in the optimal case you have one value in a fixed time interval. It is a pretty actual topic to create massive datasets on behalf of sensors in an Internet of Things scenario, and you usually get too much values (>100 measurements per second), so to create a useful time series requires some preprocessing (filter or average). Missing values can be a problem, also it is usually not difficult to find a good estimate.

You can therefore build a time line which could look like this:

Here I used R's default dataset "AirPassengers" which reflects the monthly international airline passenger numbers in the years 1949 to 1960.

Natural questions are: 

• Can you see which of the years were good, which of them were bad?
• Can you split up the time series into more homogenious components?
• Can you predict the values for future months? And how good are these predictions?

Time series analysis seems not too complex, however in reality often there are combinations of models required to make good predictions. In this post I want to show you how standard approches work. In reality, especially in retail, the timelines are usually not that easy to handle, as it is widely influenced by customer reviews.

The central idea is, that a time series $Y = (Y_t)$ is a combination of a three independent sub - time series:
- A trend component T is a long-term tendence in the data, it does not have to be linear.
- A seasonal component S is a pattern that reoccurs regularily after a fixed period (like every summer, every january or every day at 10:30).
- A random component I, also called irregular or noise.

We want to try to find these three time series in the upper mentioned example. First we have to decide on the type of decomposition, we can choose from additive and multiplicative.
In an additive model we add the 3 sub time series up to get the original time series: $$Y_t = T_t + S_t + I_t$$ You should use it when the seasonal variance does not change that much.

In a multiplicative model we multiply the 3 sub time series: $$Y_t= T_t * S_t * I_t$$ Use it when you see the peeks growing with time, like in the earlier mentioned example of airplane passengers. Here we should go for a multiplicative model.

Tip: A multiplicative model often can be changed into an additive model using the log function.

How would we get the values for the trend, seasonal and random sub time series? We will go step by step, just to motivate, here the result calculated by R with function decompose:

Here the corresponding coding in R: $$plot(decompose(ts(AirPassengers, frequency = 12, start =1949), type = "mult"))$$

To go on:

1. Here is how to determine the trend component
2. Here is how to determine the seasonal and random component
3. Here is a summary on the classical decomposition of time series

Seasonality and Random Determination

In this post we saw how the three components trend, seasonality and random of a time series. How to extract the trend was shown there, now we focus on how the seasonal component and the random component is determined.

Assume we have a detrended time series (we take here the AirPassengers time series and remove the trend). We assume a seasonality of a fixed period. In reality the assumption to have a fixed seasonality is too strict, as the period could shorten or change its structure over time. But under this assumption the determination of the seasonality is easy: To get the seasonal value of January, we take all values of January and build the average. This is the pattern we use for all periods.

The last step is to determine the random component $I$, we get it by simply removing the trend $T$ and seasonal component $S$ from the original time series $Y$, in an additive model this would be $I_t = Y_t - T_t - S_t$ and in a multiplicative model $I_t = Y_t / (T_t * S_t)$.
In our example, this is how the random components looks like:

What can we get out of it?
The random component shows the noise in the data, the values that do not fit the model. It helps to get a feeling how well the data is explained by the assumption to have a trend and a seasonality. The classical decomposition also could help to find outliers, which will show up with a high peek.

For completeness here again the whole picture holding all the steps discussed:

 

Trend determination with moving averages

We already saw in the previous post that you can decompose a trend from a time series. Here a classical approach how to determine an underlying trend.
 
A trend $T$ is actually a smoothened version a time series, it helps to capture global tendencies. To get a trend line here is what you have to do:

•  Determine the seasonal periodicity of the time series (if there is one). These periodic patterns are usually visible, but if you cannot see them form the plotted chart, there are also methods using Fourier Transform Algorithms to determine them. In our example we see a yearly periodicity, as the values for the airplane passengers come in monthly, the periodic value is m=12.

• With this number m use methods like the moving average of order m to determine the values of the trend

Now what is the moving average?
Once you understand the concept, it is easy to remember: Imagine that your dataset consists of 5 values $y_1$, ..., $y_5$. To determine the value of the trend of order $m = 3$ you would take the original value, the value of the predecessor and the value of the successor and create an average. In the simplest approach you simply take the sum of the values devided by the number of values (the so called Simple Moving Average SMA).
In the example of 5 values this would look like:
           
$y_1$      3
$y_2$      5     -> (3+5+4)/3 = 4
$y_3$      4     -> (5+4+1)/3 = 3.33
$y_4$      1     -> (4+1+3)/3 = 2.67
$y_5$      3

You already see the problem here: for the beginning and the ending there are no values available, the missing tail depends on $m$.

What if we choose $m=4$? We would have to decide to take more points from past than from future. In these cases the algorithm is not symmetric anymore, usually you therefore either change to the next odd number ($m=13$) or you choose a so called centered moving average. In the centered moving average you use a simple moving average of order 2 in the first step to determine values like $y_{1.5}$. Then you use the SMA of order $m$.

$y_1$      3
                           -> $y_{1.5}$ = (3+5)/2 = 4
$y_2$      5
                           -> $y_{2.5}$ = (5+4)/2 = 4.5
$y_3$      4                                                                           ->  (4+4.5+2.5+2) / 4 = 3.25
                           -> $y_{3.5}$ = (4+1)/2 = 2.5
$y_4$      1
                           -> $y_{4.5}$ = (1+3)/2 = 2
$y_5$      3

In our earlier example of air passengers we determined $m=12$ (even) and the trend values in yellow are determined by a centered moving average of order 12.

Here the command in R (ma stands for Moving Average):

lines(ma(x, order = 3, centre = TRUE))

After successfully determining the trend, we can remove it from the original data. In an additive model we get the de-trended time series by substracting it ($Y_t - T_t$), in a multiplicative model by deviding $Y_t/T_T$. This detrended time series of our AirPassenger example looks like this:

The next step is now to get the seasonality component and the random component from the de-trened time series.

Classical Decomposition - Summary

The classical decomposition of a time series can help to get an overview on the tendencies (trend component), periodic patterns (seasonal component) and quality of the model (random component). In addition it helps to identify outliers in a time series.

To forecast a time series it is often useful to have a decomposition and to forecast each of the components in the decomposition seperately. A seasonal component would just be repeated constantly (naive forecast), meanwhile you could use exponential smoothing methods to forecast the trend and random component.

On the other hand the classical decomposition shows some disadvantages: We saw in this post that the trend and therefore also the random component cannot be determined at the beginning and at the end of a time series. Also we saw in that post that it relies on the assumption that we have a stable period with a pretty constant pattern. In reality this is often not the case: e.g. 100 years ago the energy consumption was high in winter was high due to heaters, now in summer it is equally high due to air condition.

To overcome these bounderies other decomposition methods have been developed, see for instance the Seasonal and Trend Determination using LOESS (1990). I will describe it in a new post.

I hope you liked and got a picture on the classical decomposition, I really enjoyed building up this example and encourage you to comment and extend.

How is the weather in the Black Forest? - Nearest Neighbour Approaches

Planning a city or hiking trip on the weekend, a barbequeue or a romantic picnic in the mountains requires a pretty good weather forecast. What would you do if you wanted to know the weather? Eventually you would just check the weather forecast in the desired city, these weather forcasts are availyble online. But what if your route is lying in the mountains, without any reference city or village? You would then choose a location near your route, for which a weather forecast is available and assume you will get the same weather on your route.
Assume you have a sunny weather forecast for one side of your route and rainy weather forecast on the other side, you would probably rely on the forecast closer to the route, but also consider the forcast farer away.
This approach is quite reasonable assuming that the weather is similar on geographic similar locations. The same reasoning is also applied in the business and research world, to forecast key figures like preferences, behaviours and prizes. One approach is the Nearest Neighbour Approach which starts with the assumption that objects "close" to each other behave in the "same" way, so if you have to predict an unknown value, just ask the "nearest" neighbours and use their results.

However there are some problems in this approach: it might be easy to find neighbours considering numeric values like age, income or location, but how would you find neighbours for, lets say, prefered music?

As often there is no general answer and different ways to find neighbours exists. Usually the tasks to find a good neighbour consist of two main tasks:
1. Find a measure for distance between two datasets
2. Combine the datasets for the closest neighbours to make a prediction

Usually you have a dataset consisting of a collection of features, for a customer this could be a so called Customer ID hodling all the relevant data that could possibly influence the target variable. Often age, income, location are part of it. For these features is it easy to determine a distance, usually the difference (direct difference or euklidian difference) is used. For better comparison normalization is recommended, as the features have different ranges.

You could then go for the k nearest neigbours, check the known values of the target variables and would the have to combine these results to a target value. Here business knowledge is required to determine to which extend e.g. the age is important or if location is more relevant than income.

The advantages of this approach is that they are usually pretty easy to understand, often show additional insights and change whenever the data changes.
However they could be time consuming as for finding good neighbours you would have to check every single known dataset and measure the distance. Also the prediction could depend considerably on the value of k, further analysis is required. And the algorithm is discontinuous, meaning that a new dataset could have a huge impact on the existing predictions.

To overcome the bounderies different approaches have been established (e.g. reducing the number of datasets by choosing "important" neigbours), and the nearest neigbour approaches ars widely used in different classification or regression problems (e.g. supporting breast cancer detection, estimations on house prizes).
And defining the right distance measure you can even find neighbours to music files or detect song titles.

As a closing example here some simplified excercise. Consider the following data:

Target: Money spent online for hobbies per year
Gender Age    Income Nof Children Target
Male     35      65.000            1            6.500
Female 24      25.000           0           2.000
Female   60      60.000            4              380
Male      48      45.000            2           2.100
Female   39     60.000            0            6.000
Male      49       75.000            2  7.500
Male      18           800             0               670

Which are your closest neigbours? And how good does the estimation of the 1-nearest neighbour model fit to you?
For the distance choose the absolute difference value in age, income and number of children, take the weights age_weight = 7, income_weight = 10, nof_children_weight = 1.
(Note that this is a very simplified example, the money spent on hobbies would in reality rely on way more factors and the variables are not independent of each other).

Hopefully you get a good approximation, so your next romantic outdoor picnic in the Black Forest does not look like this:

Association Analysis

 The typical question behind Association Analysis or often also called Basket Analysis is:Which products are bought together? This question is important as based on the result different measures can be taken: You could place those products together, increase the price of one of the products and lower the price of the other one, advertice only one of them or create combo offers.

To find out about depending products, create rules R like

R: If product A is bought, then also product B is bought

Here parameters A is called Antecedents and B is called Concequent. To determine the importance of such a rule, three statisical key figures are defined:

$SUPPORT(R) := \frac{\text{number of baskets the support the rule}}{\text{number of overall baskets}}$

$CONFIDENCE(A, B) := \frac{\text{number of baskets that support the rule}}{\text{nof of baskets that contain B}}$

In lots of examples, both of the these key figures can be high, but the result is not a very useful rule (e.g. in case product A is bought by 95% of the customer). Therefore the lift, or also called improvement is introduced:
$$LIFT(A,B) := \frac{CONFIDENCE(A, B)}{SUPPORT(A)}$$
Meanwhile the support and the lift are symmetric respect A and B, the confidence is not.

Now the lift decides, if our rule is valid:


If the lift is < 1, the rule does not describe an association. For a lift of 1 the antecedents and concequents are independent of each other, and a lift > 1 describe to which degree the products depend on each other.


A typical example for an association algorithm is the so called APRIORI algorithm which creates rules for all possible subsets having a minimal support. The big advantage of it is that it produces clear, easily understandable results, which can be directly used, however the performance grows exponentially with the set of products, also very rare data is not included into the analysis.

Anomaly detection

A classical Data Science problem is to identify outliers, meaning anomalous behaviours or unexpected high or low values. These unusual values should always be analyzed, they could mean errors in the test data (to be corrected or removed from the data), they could occur naturally or they could actually be the target of an analysis. Typical applications for those analysis are e.g. fraud detection , in which a company wants to detect misuses of their products, fault detection, in which quality or security problems can be identified, but also monitoring of server and computer landscapes in order to reduce or even avoid downtimes. In these problems, the challenge is to identifying the outliers. Characteristics of such problems are, that there are only few negativ examples, but large sets of positive examples.

There are several approaches to adress this kind of challenges. Apart from the a recommended visual analysis there are a lot of algorithms that adress this problem:

A simple algorithm called Interquartile Range Algorithm calculates the Interqartile Range (IQR or also called midspread) on a set of values to find anomalous data points. It splits the number of values into 4 (equal) parts and takes the three borders between them as variables $Q_1$, $Q_2$ and $Q_3$ (ascending). The $IQR$ is then defined as $IQR = Q_3-Q_1$.

Outliers can then be defined in different ways, e.g. via a definition of the american statist John Wilder Tukey, there are two kinds of outliers:
- suspected outliers that lie $1.5 * IQR$ or more above $Q_3$ or below $Q_1$
- outliers that lie $3*IQR$ or more above $Q_3$ or below $Q_1$.
Visually they are represented by a so called "Boxplot", that is easy to understand:

The "whisker" represent the still accepted data sets, what lies outside is considered an outlier. The length is not symmetrical, it is defined by the smallest or largest values that are not yet considered an outlier. Suspected outliers are drawn in transparent circles, real outliers in filled circles.
So the middle 50% lie inside the box, the median is just another word for $Q2$.

But there are other definitions for outliers as well...


This is a simple algorithm that works on one-dimensional values. There are other algorithms that focus on distances to neigbours to find anomalous data points.

A more sophisticated approach is the following anomaly detection algorithm using the Gaussian Distribution:
To adress the issue, the idea is to start choosing "normal" behaviour for all the features that might be indicators of anomalous examples. Use those normal examples as training data for an anomaly detection algorithm:
Assume that all the features $x_1, ... x_n$ are normally distributed (Gaussian), therefore the means $\mu_i$ and the variances $\sigma_i$ of all the features in the training data is needed.

For new examples use than $p(x) = \prod_{i=1}^n (p(x_i; \mu_i, \sigma_i))$ to calculate the probability to be "normal", choose a decision boundary $\epsilon$ (e.g. $\epsilon = 0,02$) and predict anomaly examples to be the the ones with $p(x) <= \epsilon$. Here $$p(x; \mu, \sigma) = \frac{1}{\sqrt{(2\pi\sigma)}} \exp(-\frac{(x-\mu)}{2\sigma^2})$$ is the Gaussian Distribution for an $n$-dimensional value $x$ under the means $\mu = \frac{1}{m}\sum_{i = 1}^{m}x^{(i)}$ and standard deviation $\sigma = \frac{1}{m}\sum_{i = 1}^{m}(x^{(i)}-\mu)^2$.

How to choose $\epsilon$?
You can use the cross validation set which should also hold examples of anomalies to find an appropriate value. Please also make sure, that you have anomalies left for your test data.

How to measure the algorithm?
As we work here with skewed classes (there are much more positive than negative examples), accuracy is not the right way to measure the quality of the algorithm. Instead use the number of True/False Positives/Negatives or other statistical values like Precision, Recall or the $F_1$-Score.

How to come up with features:
If anomaly is not clearly distinguishable from the normal examples, try to find a new property that in compination with the existing features distinguish the normal and the anomalous examples. Start features that take on very small or very large values.


A generalization of the Gaussean Distribution Algorithm is the Multiple Gaussian Distribution Algorithm. As the upper mentioned algorithm usually yields to concentric acceptance areas (which is usually not wanted), they cannot not detect all anomalies sufficiently well. Here multiple Gaussian distribution is a useful tool to further improve the upper mentioned algorithm. However the price is paid with mor expensive calsulation:
$$p(x; \mu, \sigma) = \frac{1}{(2\pi)^{n/2}|\Sigma|^{(1/2)}}  * \exp(-\frac{1}{2} \frac{(x-\mu)^T }{\Sigma^{-1} (x-\mu)} )$$ where $|\Sigma|$ is the determinante of the matrix $\Sigma$
Here the acceptance areas depend on the values of the covariant matrix $\Sigma$ and are in general not concentric and not axis aligned.

Principal Component Analysis

Typical Data Science problems have a huge amout of input features, they are often deciding factors for model performance and make it difficult to visualize and structure data. Therefore reduction of the input features is often helpful or even necessary. The central question here is: Which are the important features? Which ones can be dropped without loosing too much data?
There are quite a few techniques to adress these questions.

Classical Techniques adressed the problem in two ways:
Backword elimination creates a model for all features, then deletes the feature that less raises the error when dropping. Forward selection goes the other way round, it chooses features stepwise choosing always the feature that improves the model the most.
As a combination of the upper mentioned techniques there is the stepwise regression: after adding a feature, a check is executed to decide which features can be deleted without increasing the error.
In all mentioned approaches the main problem is to decide, when the algorithm should stop. Also a problem is that all variables are considered independent, which in general is not the case.


More modern techniques collect subsets of features by filters or wrappers, these techniques allow insights on relations of variables, however these techniques are quite performance intense (note that the number of subsets is essentially bigger than the number of features!).


I want to present another algorithm using the Singular value decomposition of the covariant matrix:
Assume that you have an $n$-dimensial training data set of size m that is mean normalized and feature scaled, every data point $x = {x_1, ..., x_n}$ has $n$ components.
In the first step you calculate the covariance Matrix $$\Sigma =  \frac{1}{m} * \sum_{i=1}^{m} (x^{(i)} (x^{(i)})^T) $$ which is an ($n$ x $n$)-Matrix. Then determine the decomposition $$\Sigma = USV*$$ with a unitary matrices $U$ and $V$ and a diagonal matrix $S$ (use a predefined algorithm for that).
From the returning matrix $U$ choose only the first $k < n$ colums to denote a matrix $U_{k}$.

Use then $$z = U_{k}^T * x$$ with $k < n$ new features $z = {z_1, ..., z_k}$ instead of the existing $n$ ones.

To get back to the original parameters from $z$ to $x$ use the following approximation:
$$x \sim x_{approx} = U_{k}^T * z$$ which only works, if $k$ is properly chosen.

How should $k$ be chosen?
You can choose $k$ = smallest value so that the average squared projection error devided by the total variation in the data is below a boundary $\epsilon$, e.g. $\epsilon = 0.01$ to say that "99% of variance is retained":

$$\frac{\frac{1}{m} \sum_{i=1}^{m} || x^{(i)} - x_{approx}^{(i)} ||^2}{\frac{1}{m} \sum_{i=1}{m} ||x(i)||^2 } <= \epsilon$$
There is even a faster, direkt way to determine $k$, using the diagonal matrix $S$:
Chosse the smallest value for $k$ so that
$$1 -  \frac{\sum_{i=1}^{k} S_{ii}}{\sum_{i=1}^n S_{ii}} <= \epsilon$$ or $$ \frac{\sum_{i=1}^{k} S_{ii}}{\sum_{i=1}^n S_{ii}} >=1-\epsilon$$


Notes: 
- PCA should be applied to the training set to determine the parameters of the mapping from $x$ to $z$. However the reduced matrix $U_{reduce}$ can later still be used for the data in the cross and test set.
- PCA should not be used to adress overfitting (as data is getting lost, not using the exising results), therefore the better solution is using regulatization parameters $\lambda$.
- PCA is often not needed, so the recommended approach is to first try to run an algorithm without PCA and only apply it, if it is really needed (visualization, performance).

Machine Learning

What does "machine learning" mean, how can a computer "learn"?

The term "learning" is used in this context to refer to the fact that a machine uses new data in order to better up previous results. For sure the machine will not develop a brain or a physically similar organ. But for optimization or minimizing tasks, especially tasks including complex calculations and transformations, a machine controlled adaption is a useful and often indispensable tool. The basics for such tasks are often machine learning algorithms like artificial neural networks or clustering algorithms which run iterative trying to minimize errors in each calculation step.

Machine learning algorithms usually can be devided in two main groups:

• Supervised learning algorithms, in which the output is known and the rules are being trained by using the input and trying to minimize the error between the predicted output and the given result. Neuronal networks are here the most promising examples, they are called like that as they reflect the way our brains work: given some input the human brain learns by trying out and correct until it finds the perfect rule to explain the result.

• Unsupervised learning algorithms in which data is given to an algorithm which then tries to find pattern in the data. As an example think about astronomical data: if you can cluster the stellar data you could find a structure in it and learn about the past and future.

In addition modern algorithms can also be active learning algorithms, in which the input can be completed by additional requests for input trying to minimize the numbers of these additional inputs.

What are these algorithms needed for and why are they considered promising?

Machines have the capacity to calculate fast and storage is cheap nowadays, also in presicion they are unbeatable and often parallelize their tasks. As nearly unlimited data is available and often more data leads to better results (not always, an intelligent way to sort out data is one of the reasons for a data scientist!) machines and computer visualizations are essential reasons for machine supported analysis. In addition cloud computing and distributed systems improve the way data is collected, loaded and analyzed.
Internet of Things scenarios are considered the modern way to improve business processes and drive Industry 4.0 by collecting huge amounts of (sensor) data. To analyze those datasets machines are not only helpful, but necessary.