Getting Started with Time Series Analysis (1)

Overview

I would like to organize my knowledge on the theoretical and practical aspects of "Time Series Analysis," which is particularly challenging among statistical methods.

The motivation was that when I tried to study the theoretical parts using "Time Series Analysis" (Kyoritsu Shuppan, 2017), I didn't make much progress, perhaps because I wasn't used to the thinking patterns specific to time series analysis. That said, simply knowing how to use libraries and write code just to get the analysis done is meaningless.^[1]

Recently, I watched the Course on "Time Series Analysis" by Professor Kitagawa, who has authored many masterpieces, offered as a public lecture by the University of Tokyo, and my understanding improved significantly.

In this article, referring to the public lecture and other literature, I will organize "preprocessing and analysis/visualization of statistics," which are important in the analysis of time series data. In doing so, I will also provide programming examples while using datasets.

What is Time Series Data?

First, time series data is data that records phenomena that fluctuate over time (e.g., rainfall, stock prices, number of infected persons, seismic motion). And time series analysis is a method for understanding the complex phenomena behind it based on time series data, and for making predictions, control, and decisions.

In time series analysis, it is important to bring in time series models that match the characteristics of the data, which greatly affects the results and interpretation of the analysis. In particular, if time series data is assumed to be linear, stationary, and normal, very easy-to-handle and powerful models can be used.

However, it is rare for real-world time series data to have such ideal properties. While one could handle such complex data as is, one should also consider transforming the data itself into something simpler (as much as possible) to simplify the modeling.

Among these, the property called stationarity is particularly important in time series analysis. "Stationary" means that the properties of the time series data (mean, variance, covariance) do not change over time. Due to these properties, several time series models work very effectively (AR models, MA models, ARMA models, etc.).

Why is Stationarity Important?

A reason why the stationarity of time series data is important is the ease of modeling mentioned earlier. Related to that, there is also the aim of matching the premises necessary for analyzing time series data using statistical methods.

For general data that is not time series data, the Law of Large Numbers is effectively used to estimate the characteristics of the population from the samples at hand. This states that if the sample size is sufficient, the estimator for the sample will match (approach) the actual quantity of the population.

You might think that the Law of Large Numbers would work for time series data as well since a lot of data is collected over time, but that is actually not the case.

The following passage often seen in time series analysis books explains the difficulty of handling time series data:

What makes time series analysis most different from other statistical analyses is the order of observations, that is, the fact that it deals with values that appear along the flow of time. Since time only moves forward and cannot go back, it cannot be observed repeatedly in real time, and statistical inference must be performed from a single time series, that is, a single sample path.
(From "Time Series Analysis" (Kyoritsu Shuppan, 2017))

To explain in more detail, originally, to obtain a model or analytical quantity to predict the observed value $y_n$ at time point $n$, one would have to collect many observed values $y_n$ at this time point $n$... but this is a task that simply cannot be done with time series data.^[2]

Despite such troublesome nature, if stationarity is assumed, it becomes acceptable to use data from other time points (e.g., $y_{n-3}$, $y_{n+1}$) for the analysis of the observed value $y_n$ at time point $n$. This is because no matter which point in time the data is taken from, there is no change in the properties of the generative model of the time series data throughout time.

Therefore, treating stationary time series data as a time average (ergodicity) is important for analysis.

Visualizing Time Series Data

We would like time series data to show "clean" stationarity if possible, but it is rare for real-world data to have ideal stationarity. Therefore, it is necessary to convert it into data that is easy to analyze through some kind of preprocessing.

To perform preprocessing, it is important to know what properties the original time series data has. So, once you have time series data, I recommend visualizing it in a graph or similar format and checking it visually.

Exploring Various Datasets

The datasets for explanation will be taken from the R package TSSS.

install.package("TSSS")

library("TSSS") # Loading the TSSS package

Ship's directional angular velocity data: HAKUSAN

data(HAKUSAN)
hakusan1 <- HAKUSAN[:,1] # Obtain only directional angular velocity data
plot(hakusan1)

Looking at the data,

1. A uniform mean line could likely be drawn somewhere around $y \in [-2, 0]$ (Mean Stationarity)
2. It shows similar fluctuations throughout the entire time period (Variance/Covariance Stationarity)
3. The fluctuation patterns seem to have periodicity

Sunspot number data: Sunspot

data(Sunspot) # Reading sunspot number data
plot(Sunspot)

Looking at the data,

1. Since it is count data, waves are only in the upper half (Positive values, Vertical asymmetry)
2. Counts reach local maxima periodically, but it is unclear if the period is stable (Pseudo-periodicity)
3. The increase and decrease in counts do not seem to have clean symmetry... if anything, the decrease seems slightly slower (Time-reversal asymmetry)

Daily maximum temperature data: Temperature

data(Temperature) # Daily maximum temperature data
maxtemp <- Temperature
plot(maxtemp, ylim = c(0, 35)) # Set the y-axis scale from 0 to 35

Looking at the data (considering it is collected daily),

1. There seems to be a clear trend in temperature changes within the year (Long-term periodicity)
2. Around the trend, other components seem mostly stationary

MYE1F (East-West Seismic Wave) Data: MYE1F

data(MYE1F)
plot(MYE1F)

Looking at the data,

1. There seems to be no trend.
2. The mean is stationary (constant over time?), but the variance and covariance are non-stationary.
3. However, it feels stationary if restricted to certain intervals (Piecewise Stationary).

Haibara (Groundwater Level and Air Pressure) Data: Haibara

data(Haibara)
plot(Haibara)

Looking at the data,

1. It is bivariate time series data observed at the same time points (groundwater level and air pressure data).
2. The two time series fluctuate in opposite directions (Inverse Correlation).
3. There are missing or jumping data points here and there (Abnormal values, Missing values).

Classification of Time Series Data Properties

Time series data has various properties.
In general, time series can be classified as follows.

Continuous Time Series and Discrete Time Series

Data recorded continuously in time is called a continuous time series, and data collected at certain time intervals (e.g., daily) is called a discrete time series.

In practice, since data collection usually happens at specific timings, we mostly deal with discrete data as time series data. Furthermore, discrete time series can be divided into equispaced time series (e.g., every 10 minutes) and unequally spaced time series where the collection timing is irregular.

Univariate Time Series and Multivariate Time Series

Data where only one type of information is obtained at each observation point is a univariate time series, and data where two or more types of information are acquired simultaneously is a multivariate time series.

In particular, for multivariate time series, considering the correlation between time series data can improve prediction accuracy compared to analyzing a univariate time series alone.

Stationary Time Series and Non-stationary Time Series

When looking at time series data, all of them show irregular fluctuations, but in time series analysis, we express these "irregular fluctuations" using stochastic models.
At this time, even for data that appears to fluctuate irregularly at first glance, if the properties of the underlying stochastic model do not change over time, it is called a stationary time series.
On the other hand, if the properties of the stochastic model change as time progresses, it is called a non-stationary time series.

Care must be taken when handling non-stationary time series data. However, some non-stationary time series data can be considered stationary by applying some kind of preprocessing (e.g., log transformation, differencing, etc.).

Linear Time Series and Non-linear Time Series

When a time series model is expressed as a linear sum (e.g., an AR model determined by a weighted sum of past data, or when the original data consists of the sum of several components), it is called a linear time series. On the other hand, those requiring more complex non-linear models are called non-linear time series.

Gaussian Time Series and Non-Gaussian Time Series

Data where the distribution of values observed as time series data follows a normal distribution is called a Gaussian time series, and other cases are called non-Gaussian time series.

If it is difficult to imagine how to determine if a time series is Gaussian... it is easier to understand if you count the time series data regarding observed values—in other words, create a histogram.

data(HAKUSAN)
hakusan1<- HAKUSAN[, 1]
plot(hakusan1, main = "")
hist(hakusan1, col = "#539952")

This is a histogram created for the ship's directional angular velocity data. Looking closely at the histogram, we can see that the observed values $y_n$ spread symmetrically. This is the appearance of a Gaussian time series.

In that case, if there is a positive trend (the overall amount increases) or if the data is vertically asymmetrical, the histogram will not form a clean normal distribution. This is a non-Gaussian time series.

Even for non-Gaussian time series, it is sometimes possible to bring it closer to a normal distribution by performing some kind of preprocessing (e.g., log transformation).

Missing Values and Outliers

In some time series data, observed values may not be recorded for some reason (missing values), or clearly abnormal data (abnormal values, outliers) may be obtained due to measurement errors or problems during data transmission.

By using state-space models, one can perform time series analysis without being overly conscious of such abnormalities.

Preprocessing of Time Series Data

I explained that time series data can be classified in various ways. Among various properties of data, time series such as normality, stationarity, and linearity are very easy to handle. However, the fact remains that such ideal time series data is scarce.

Nevertheless, even data that at first glance is non-stationary, non-Gaussian, or non-linear can be converted into easy-to-handle data—such as having normality, stationarity, and linearity—by applying some kind of preprocessing.

Variable Transformation

First, as a basic element, I will introduce how to perform variable transformations.
As an example, assume that the data $X$ at hand is generated from some distribution $f(x)$ (assume this $f(x)$ is unknown).

And assume that the data $x$ can be transformed into $y$ using a transformation function $h$. At this time, assume that the inverse transformation is also defined ($y = h(x), x = h(y)^{-1}$). It is known that $y$ after transformation is generated from some known distribution $g(y)$.

Then, the distribution $f(x)$ that generates data $X$ can be obtained as follows:

$$f(x) = g(h(x))\Big|\frac{dh}{dx}\Big|$$

Here, $\frac{df}{dx}$ is called the Jacobian matrix (Jacobian) and is obtained by the following formula:

$$\frac{\partial {h}}{\partial {x}} = \begin{pmatrix} \frac{\partial {h_1}}{\partial {x_1}} & \cdots & \frac{\partial {h_1}}{\partial {x_m}} \\ \vdots & \ddots & \vdots \\ \frac{\partial {h_m}}{\partial {x_1}} & \cdots & \frac{\partial {h_m}}{\partial {x_m}} \end{pmatrix}$$

By variable transformation, the probability distribution is also deformed, and at the same time, the domain may change.

Log Transformation

As a simple example, I will explain the log transformation. Since the transformation function $h(x)$ is defined as $h(x) = \log(x)$, the Jacobian is obtained as $\frac{dh}{dx} = \frac{1}{x}$. If we assume that the data $y$ follows a normal distribution $g(y)$, then $g(y) = \frac{1}{\sqrt{2\pi\sigma^{2}}}\exp^{-\frac{(y-\mu)^2}{2\sigma^{2}}}$.

Then, the shape of $f(x)$ that generates data $x$ becomes the following log-normal distribution.

$$f(x) = g(h(x))\Big|\frac{dh}{dx}\Big| = \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp^{-\frac{(\log{x} - \mu)^2}{2\sigma^{2}}} \frac{1}{x}$$

As an example, we apply log transformation to the WHARD data in the TSSS package.

data(WHARD)
plot(WHARD) # WHARD data before transformation
plot(log(WHARD)) # WHARD data after log transformation

The figure above shows the WHARD data (data recording the monthly wholesale sales of certain hardware) plotted before (original) and after log transformation.

The non-stationarity of variance seen in the original data is mitigated, and the trend is maintained after transformation.

Box-Cox Transformation

I will introduce the convenient Box-Cox transformation, which can create various power-type transformations, including log transformation.

First, for the definition:

$$\begin{equation} y_{\lambda} = \begin{cases} \frac{x^{\lambda}-1}{\lambda} & \lambda \neq 0 \\ \log{x} & \lambda = 0 \\ \end{cases} \end{equation}$$

As you can see, various transformations can be expressed with the value of $\lambda$. For example, the aforementioned log transformation corresponds to the Box-Cox transformation where $\lambda = 0$.

Then, the question arises: how is $\lambda$ determined? This can be done by performing Box-Cox transformations with various values of $\lambda$ and evaluating whether the transformed data is close to a normal distribution to find the best value for $\lambda$.

Parameter Selection for Box-Cox Transformation

I will explain how to determine the parameter $\lambda$ of the Box-Cox transformation using AIC (Akaike Information Criterion).

Broadly speaking, this involves performing Box-Cox transformations on the original data for each candidate $\lambda$ and evaluating with AIC whether the distribution of the transformed data is close to a normal distribution.

The TSSS package includes a boxcox function that automatically adjusts $\lambda$ and performs the Box-Cox transformation using that value.

data(Sunspot)
boxcox(Sunspot) # Box-Cox transformation (automatic parameter adjustment)

Plot of Sunspot data before conversion and after Box-Cox transformation ($\lambda: 0.4$)

For the Sunspot data, the optimal parameter $\lambda$ was determined to be 0.4.

Looking at the results, it can be seen that the original data, which was vertically asymmetrical, has become vertically symmetrical after transformation. In other words, the time series data has been converted to a Gaussian type.

Actually incorrect AIC comparison after transformation

The boxcox function automatically provides the optimal value of $\lambda$ based on AIC, but caution is required if you implement the process from scratch or calculate AIC using other software.

When considering the optimal $\lambda$ based on AIC, the $AIC_z$ calculated using the transformed data $z$ should actually not be compared. If you want to compare, you need to find the $AIC'_{z}$ in terms of the original data.

This value is obtained by correcting $AIC_z$ with the Jacobian of the transformation function.

$$AIC'_{z} = AIC_{z} - 2\log{|\frac{df}{dx}\Big|}$$

Below are the $AIC_z$ for the transformed data and the $AIC'_{z}$ for the original data. Note that $AIC_z$ is not minimized at $\lambda = 0.4$.

AIC Results

lambda	aic'	LL'	aic	LL	mean	variance
1.00	2360.37	-1178.19	2360.37	-1178.19	4.910866e+01	1.576345e+03
0.90	2335.48	-1165.74	2174.59	-1085.29	3.546734e+01	7.052906e+02
0.80	2313.88	-1154.94	1992.09	-994.04	2.591715e+01	3.200762e+02
0.70	2295.85	-1145.93	1813.16	-904.58	1.917797e+01	1.475255e+02
0.60	2281.75	-1138.87	1638.16	-817.08	1.438208e+01	6.915929e+01
0.50	2271.96	-1133.98	1467.48	-731.74	1.093834e+01	3.303407e+01
0.40	2267.00	-1131.50	1301.62	-648.81	8.441855e+00	1.611155e+01
0.30	2267.51	-1131.76	1141.24	-568.62	6.613792e+00	8.046373e+00
0.20	2274.40	-1135.20	987.23	-491.61	5.260971e+00	4.131024e+00
0.10	2289.01	-1142.51	840.94	-418.47	4.248735e+00	2.192978e+00
0.00	2313.40	-1154.70	704.44	-350.22	3.482608e+00	1.214490e+00
-	-	-	(omitted below)	-	-	-

lambda = 0.40 AIC' minimum = 2267.00

:::

Logit Transformation

A transformation that converts ratio data $x (0 < x < 1)$ into the natural logarithm of the odds $\frac{x}{1-x}$ is called a logit transformation. Because the original data is a probability (ratio), its domain is $[0, 1]$, but the characteristic feature is that after transformation, it becomes $[-\infty, \infty]$.

$$y = \log{\frac{x}{1-x}}$$

Differencing

If the mean or variance is not constant due to some influence, extracting fluctuations from a certain point in time may actually reveal stationary changes. Therefore, preprocessing to calculate differences can also be effective.

Let the $k$-th order difference of the time series $y_n$ be defined as $\Delta^{k}y_{n}$:

$$\begin{align} \Delta^{1}y_{n} &= y_n - y_{n-1} (\text{for first-order difference}) \\ \Delta^{2}y_{n} &= \Delta^{1}y_{n} - \Delta^{1}y_{n-1} = y_n - 2y_{n-1} + y_{n-2} (\text{for second-order difference}) \\ \end{align}$$

In fact, taking the (log) difference of the Nikkei225 data in the TSSS package:

data(Nikkei225)
plot(Nikkei225) # Nikkei225 data before differencing

plot(diff(log(Nikkei225)))  # Nikkei225 data after log differencing

The top figure is the original Nikkei225 (Nikkei 225 Stock Average) data, and the bottom figure shows the result after taking the log difference. Data that was originally mean non-stationary has been transformed into clean mean-stationary data centered around zero.

Here, instead of the data from one time point prior, you can also take the difference from even earlier data. Especially for data with periodicity like seasonality, seasonal fluctuations can be removed by comparing it with data from $s$ periods ago, where $s$ is the length of the season.

$$\Delta_{s} y_{n} = y_n - y_{n-s}$$

Seasonal differencing with s = 12 for WHARD data

data(WHARD)
y <- log(WHARD)
n <- length(WHARD)
period <- 12 # Periodicity every 12 months
z <- rep(NA, n)
for (i in period + 1:n) {
    z[i] <- y[i] - y[i-period] # Seasonal differencing
}
plot(as.ts(z))

Other ways to take differences include the year-over-year ratio $x_n = y_n/y_{n-1}$ or the year-on-year ratio for the same period $x_n = y_n / y_{n-p}$.

Moving Average Filter and Moving Median Filter

Time series data with irregular fluctuations contains error fluctuations (noise) in addition to components meaningful for interpretation, such as trends and seasonal fluctuations. If this error fluctuation (noise) is too large, it is difficult to read meaningful patterns from the original time series data.

Therefore, smoothing is effective as a method for removing error fluctuations (noise) and cleanly extracting patterns of change.

Well-known methods include the moving average filter and the moving median filter. These are smoothing methods that replace data at a certain time point with the average or median of the data before and after that point, respectively.

The moving average filter can be expressed as:

$$t_n = \frac{1}{2k+1}(y_{n-k} + \dots + y_{n} + \dots + y_{n-k})$$

The moving median filter can be expressed as:

$$t_n = median(y_{n-k},\dots,y_{n},\dots,y_{n-k})$$

data(Temperature)
maxtemp <- Temperature

# In the case of a moving average filter (using the SMA function)
plot(maxtemp, ylim = c(0,40))
x <- SMA(maxtemp, 34) # Use the SMA function from the TTR package; window size is 34
lines(x, col = 2, lwd = 2) # col: color, lwd: thickness

# In the case of a moving average filter (implemented from scratch)
plot(maxtemp, ylim = c(0,40))
y <- maxtemp
ndata <- length(maxtemp)
y[1: ndata] <- NA
kfilter <- 17
n0 <- kfilter + 1
n1 <- ndata - kfilter
for (i in n0:n1) {
    i0 <- i - kfilter
    i1 <- i + kfilter
    y[i] <- mean(maxtemp[i0:i1])
}
lines(y, col=2, ylim = c(0, 40), lwd = 2)

# In the case of a moving median filter (using the runMedian function from TTR)
data(Temperature)
maxtemp <- Temperature
plot(maxtemp, ylim = c(0,40))
x <- runMedian(maxtemp, 34) # Use the runMedian function from the TTR package; window size is 34 (17 + 17)
lines(x, col = 3, lwd = 2) # col: color, lwd: thickness

# In the case of a moving median filter (implemented from scratch)
plot(maxtemp, ylim = c(0,40))
y <- maxtemp
ndata <- length(maxtemp)
y[1: ndata] <- NA
kfilter <- 17
n0 <- kfilter + 1
n1 <- ndata - kfilter
for (i in n0:n1) {
    i0 <- i - kfilter
    i1 <- i + kfilter
    y[i] <- median(maxtemp[i0:i1])
}
lines(y, col=3, ylim = c(0, 40), lwd = 2)

For the moving average filter

For the moving median filter

Since the results of both methods change depending on the window size (the number of time series data points to consider), it is necessary to consider smoothing with an appropriate window size.

Also, a comparison between the moving average filter and the moving median filter is as follows.

	Moving Average Filter	Moving Median Filter
Pros	Smooth estimates	Accurately detects structural changes. Robust to outliers.
Cons	Cannot accurately detect structural changes (ignores them to smooth over). Sensitive to outliers.	May conversely increase fluctuations.

Summary

In this post, I organized the nature of time series data by visualizing actual datasets and explaining what time series data is. I also reviewed techniques such as differencing, variable transformation, and smoothing as preprocessing steps to facilitate modeling. In the next post, I will introduce important quantities of time series data.

Miscellaneous Notes

Although there are already many articles on methods of time series analysis, few explain the quirks and procedures of time series analysis from a practical perspective. I remember having a hard time resolving the question of why we transform data into stationary data.

I feel fortunate to have learned not only the theory of time series models but also practical analysis procedures and how to interpret results through Professor Kitagawa's course on time series analysis.

However, it's terrifying that the volume of this post covers only about the first lecture of the course... 😱

References

https://ocw.u-tokyo.ac.jp/course_11477/

Course on time series analysis by Professor Kitagawa. Very easy to understand, and includes advanced topics not found in common articles.

https://www.kyoritsu-pub.co.jp/book/b10003870.html

Introduces theoretical explanations of stationary time series models. I felt it was difficult as a beginner (like myself), but after learning various things, I realized the theory is summarized beautifully and concisely.

https://www.kyoritsu-pub.co.jp/book/b10003204.html

Introduces code examples in Python.

脚注

1. As for data science, I want to be able to deeply understand the nature of the data in front of me, and then choose appropriate models and interpret the results!

2. For example, the price of ice cream in 1980 can only be obtained once at that single point in time. To force the Law of Large Numbers to work, you would have to collect data on the price of ice cream in 1980 in a "different world line."