3.2.1 Curve Fitting

This technique is not more than removing the trend and analyze the “residuals.” In general, particularly for yearly data, we can use a simple polynomial curve.

oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))

# Our old friend, AirPassengers dataset. See how it increases both the mean and the variance.
plot(AirPassengers, main = "Monthly totals of international airline passengers", ylab = "Passengers")
# Let's just fit a line on the data
fit <- lm((AirPassengers) ~ seq_along(AirPassengers))
pred <- predict(fit, data.frame(seq_along(AirPassengers)))
pred <- ts(pred, start = start(AirPassengers), frequency = frequency(AirPassengers))
# `pred` contains our line based on the simple linear model we did.
data <- AirPassengers / pred # why divide and not subtract? because the variance also increases (we have a multiplicative seasonality)
lines(pred, col = "red")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(data, ylab = "Passengers / fitted")

par(oldpar)

3.2.2 Filtering

Filtering is a little more complex operation than curve fitting. We are not just trying to find a polynomial that best fits our data, but we are transforming our original TS into another TS using a formula (that here we call filter). This filter can be one of several kinds of “Moving Averages,” locally weighted regressions (e.g., LOESS), or “Splines” (a piecewise polynomial). One caveat of these smoothing techniques is the end-effect problem (since in one end of the time series, we do not have all the values to compute, for example, the moving average).

Simple moving average ex: \(Sm(x_t) = \frac{1}{2q+1} \sum^{+q}_{r = -q}x_{t+r}\)

Example with splines:

oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))

plot(USAccDeaths, main = "Accidental Deaths in the US 1973-1978", ylab = "Accidental Deaths")

# here is another way to "fit a line". The number of knots is your choice, but an old Professor once told that 4-5 per year is sufficient.
pred <- smooth.spline(USAccDeaths, nknots = 24)$y
pred <- ts(pred, start = start(USAccDeaths), frequency = frequency(USAccDeaths))
lines(pred, col = "red")
par(mar = c(5.1, 4.1, 0.1, 1.1))
data <- USAccDeaths - pred # see how subtraction and division only affects the mean on this case.
plot(data, ylab = "Deaths - spline")
par(mar = c(5.1, 4.1, 0.1, 1.1))
data <- USAccDeaths / pred
plot(data, ylab = "Deaths / spline")

par(oldpar)

3.2.3 Differencing

It is a special kind of filtering, where we compute the difference between the current value and the next. It is helpful to remove trends, making a TS stationary. We can use differencing multiple times (we call “orders”), but usually, one (“first-order”) iteration is sufficient. The mathematical operator used to denote differencing is the “nabla” (\(\nabla\)). \(\nabla^2\) means second-order differencing.

oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))

plot(co2, main = "Mauna Loa Atmospheric CO2 Concentration", ylab = "CO2 concentration")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(diff(co2), ylab = "diff(co2)") # the first difference removes all the trend!

par(oldpar)

See that the example above removed the trend, but kept the seasonality.

The Slutzky-Yule effect[2]: They showed that by using operations like differencing and moving average, one could induce sinusoidal variation in the data that, in fact, is not real information.

3.4.1 The correlogram

The graphic representation of the autocorrelation coefficients is called a correlogram, in which the coefficient \(r\) is plotted against the lag \(k\).

Interpretation. Here we offer some general advice:

• Random series: for a large \(N\), \(r_k \simeq 0\) for all non-zero values of \(k\). Usually, 19 out of 20 of the values of \(r_k\) lie between \(\pm 2/\sqrt{N}\).

• Short-term correlation: Stationary series usually presents with a short-term correlation. What we see is a large value for \(r_1\), followed by a geometric decay.

• Alternating series: If the data tends to alternate sequentially around the overall mean, the correlogram will also show this behavior: the value of \(r_1\) will be negative, \(r_2\) will be positive, and so on.

• Non-stationary series: If the data has a trend, the values of \(r_k\) will not come to zero, only for large lag values. This kind of correlogram has little to offer since the trend muffles any other features we may be interested in. Here we can conclude that the correlogram is only helpful after removing any trend (in other words, turn the series stationary).

• Seasonal fluctuations: If the data contains seasonal fluctuations, the correlogram will display an oscillation of the same frequency.

3.4.2 Testing for randomness

There are valuable tools to test if the data is random or not. For the sake of simplicity, this subject will not be covered in this article. However, testing residuals for randomness is a different problem[1] and will be discussed later.

4.1.1 Second-order Stationarity

A more helpful definition, for practical reasons, is a less restricted definition of stationarity where the mean is constant, and ACVF only depends on the lag. This process is called second-order stationary.

This simplified definition of stationarity will be generally as long as the properties of the processes depend only on its structure as specified by its first and second moments.

4.2.1 Property 1

The ACF is an even function of the lag in that

\[ \rho(\tau)=\rho(-\tau) \]

This property just states that the correlation between \(X(t)\) and \(X(t + \tau)\) is the same as that between \(X(t)\) and \(X(t-\tau)\). The result is easily proved using \(\gamma(\tau)=\rho(\tau)\sigma^2\) by

\[ \begin{aligned} \gamma(\tau)&=Cov[X(t), X(t+\tau)] \\ &=Cov[X(t-\tau), X(t)] \qquad \text{since} \; X(t) \; \text{stationary} \\ &=\gamma(-\tau) \end{aligned} \]

4.2.2 Property 2

\(|\rho(\tau)| \leq 1\). This is the “usual” property of a correlation. It is proved by noting that

\[ Var[\lambda_1 X(t) + \lambda_2 X(t+\tau)] \geq 0 \]

for any constants \(\lambda_1, \lambda_2\) since variance is always non-negative. This variance is equal to

\[ \lambda_1^2 Var[X(t)] + \lambda_2^2 Var[X(t + \tau)] + 2 \lambda_1 \lambda_2 Cov[X(t), X(t + \tau)] = (\lambda_1^2 + \lambda_2^2) \sigma^2 + 2 \lambda_1 \lambda_2 \gamma(\tau) \]

When \(\lambda_1 = \lambda_2 = 1\), we find

\[ \gamma(\tau) \ge -\sigma^2 \]

so that \(\rho(\tau) \ge -1\). When \(\lambda_1 = 1, \lambda_2 = -1\), we find

\[ \sigma^2 \ge \gamma(\tau) \]

so that \(\rho(\tau) \le +1\)

4.2.3 Property 3

Lack of uniqueness. A stochastic process has a unique covariance structure. However, the opposite is not valid. We can find other processes that produce the same ACF, adding another level of difficulty to the sample ACF interpretation. To overcome this problem, we have the invertibility condition that will be described later on Moving average processes.

4.3.1 Purely random process

A process can be called purely random if it is composed of a sequence of random variables \(\{Z_t\}\) that are independent and identically distributed. By definition, it has constant mean and variance, given that

\[ \gamma(k)=Cov(Z_t, Z_{t+k}) = 0 \qquad \text{for}\;k=\pm 1, 2, \dots \]

Since the mean and ACVF do not depend on time, the process is second-order stationary. In fact, it also satisfies the condition for a strictly stationary process. The ACF is given by

\[ \rho(k)=\left\{ \begin{array}{ll} 1 & k=0 \\ 0 & k=\pm 1, \pm 2, \ldots \end{array}\right. \]

set.seed(2021)
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))

normal_random <- rnorm(500)
plot.ts(normal_random, main = "Purely random series - N(0, 1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(normal_random, lag.max = 1000)

par(oldpar)

This type of process is of particular importance as building blocks of more complicated processes such as moving average processes.

4.3.2 Random walk

The Random Walk is a process very similar to the previous process. The difference lies in that the current observation sums the current random variable to the previous observation instead of being independent. The definition is given by

\[ X_t=X_{t-1}+Z_t \]

Usually, the process starts at zero for \(t=0\), so that

\[ X_1=Z_1 \]

and

\[ X_t= \sum_{i=1}^t Z_i \]

Then we find that \(E(X_t)=t\mu\) and that \(Var(X_t)=t\sigma^2_Z\). As the mean and variance change with \(t\), the process is non-stationary.

set.seed(2021)

random_walk <- cumsum(rnorm(500))
plot.ts(random_walk, main = "Random Walk - N(0, 1)", ylab = "Value")

Meanwhile, the interesting feature is that the first difference of a random walk forms a purely random process, which is therefore stationary.

set.seed(2021)

random_norm <- rnorm(500)
random_walk <- cumsum(random_norm) # same from the last plot
diff_walk <- diff(random_walk)
# below we use 2:500 because with diff, we lose the first observation
all.equal(diff_walk, random_norm[2:500])
## [1] TRUE

4.3.3 Moving average processes

First, not to be confused with the moving average algorithm. The moving average process is a common approach to model a univariate time series. The concept is that the current value \(X_t\) depends linearly on \(q\) past values of a stochastic process. Another practical way to see this process is imagining the process as a finite impulse applied to a white noise. This impulse “has” affected the \(q\) previous values and the current.

The moving average process only remembers the \(q\) previous components of the random process³, so it is also limited to \(q\) steps in the future. After that, one cannot predict any value without new random values being generated[4].

Here we will say that \(\{Z_t\}\) is a process that only generates purely random values with mean zero and variance \(\sigma^2_Z\). Then the process \(\{X_t\}\) can be said to be a moving average process of order \(q\) (abbreviated to a \(\text{MA}(q)\) process) if

\[ X_t=\beta_0Z_t+\beta_1Z_{t-1}+\cdots+\beta_qZ_{t-q} \tag{4.4.3.1} \]

where \(\{\beta_i\}\) are constants.

This equation is very similar to a linear regression \(y = a + bx\) where the dependent process \(\{X_t\}\) is modeled by an independent process \(\{Z_t\}\) (a purely random process). Here we omit, for simplicity, the “intercept,” that may be a constant \(\mu\) added to the end of the right side of the equation. The “intercept” is the overall mean of this process.

Usually, the random process is scaled so that \(\beta_0=1\). Then we find that

\[ E(X_t)=0 \\ Var(X_t)=\sigma^2_Z \sum^q_{i=0} \beta_i^2 \]

since the \(Z\)s are independent. We also have

\[ \begin{aligned} \gamma(k) &= Cov(X_t, X_{t+k}) \\ &= Cov(\beta_0Z_t+\cdots+\beta_qZ_{t-q}, \beta_0Z_{t-k}+\cdots+\beta_qZ_{t+k-q}) \\ &= \left\{\begin{array}{cc} 0 & k>q \\ \sigma_{Z}^{2} \sum_{i=0}^{q-k} \beta_{i} \beta_{i+k} & k=0,1, \ldots, q \\ \gamma(-k) & k<0 \end{array}\right. \end{aligned} \]

since

\[ Cov(Z_{s}, Z_{t})=\left\{\begin{array}{ll} \sigma_{Z}^{2} & s=t \\ 0 & s \neq t \end{array}\right. \]

As \(\gamma(k)\) does not depend on \(t\), and the mean is constant, the process is second-order stationary for all values of the \(\{\beta_i\}\). Furthermore, if the \(Z\)s are normally distributed, so are the \(X\)s, and we have a strictly stationary normal process.

The ACF of the \(\text{MA}(q)\) process is given by

\[ \rho(k)=\left\{\begin{array}{cl} 1 & k=0 \\ \sum_{i=0}^{q-k} \beta_{i} \beta_{i+k} / \sum_{i=0}^{q} \beta_{i}^{2} & k=1, \ldots, q \\ 0 & k>q \\ \rho(-k) & k<0 \end{array}\right. \]

Note that the ACF is “clipped” to zero at lag \(q\). This is a feature of \(\text{MA}\) processes that we can spot on ACF plots. In practice, the “zero” will be a value below the significance line.

oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
set.seed(2021)

mov_avg <- arima.sim(list(order = c(0,0,1), ma = 0.8), n = 200)
plot.ts(mov_avg, main = "Moving Average MA(1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(mov_avg)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(mov_avg)

par(oldpar)

There is no restriction on \(\{\beta_i\}\) values to produce a stationary \(\text{MA}\) process. However, as we briefly mentioned at the ACF Property 3, it is desirable that the process is invertible (e.g., Box and Jenkins, 1970, p. 50).

4.3.4 Autoregressive processes

In the previous section about the moving average process, we imagined the process as an experiment where you had an impulse applied to a random process with a finite time span influence. Here we can think of an experiment where this impulse persists in time, but its influence is not visible immediately, but we see it as a repeatable pattern over and over. This also implies that the \(\text{AR}\) process may not be stationary, in contrast with \(\text{MA}\) process. Moreover, \(\text{MA}\) process only “remembers” the previous components of the underlying random process, where the \(\text{AR}\) process depends directly on the previous observations, hence the prefix “auto” regressive.

Again, we will say that \(\{Z_t\}\) is a process that only generates purely random values with mean zero and variance \(\sigma^2_Z\). Then the process \(\{X_t\}\) can be said to be an autoregressive process of order \(p\) if

\[ X_t=\alpha_1 X_{t-1} + \cdots + \alpha_p X_{t-p} + Z_t \tag{4.3.4.1} \]

In contrast with the \(\text{MA}\) process, the \(\text{AR}\) process looks like a multiple regression model since \(X_t\) is regressed not on independent variables but past values of \(X_t\). An autoregressive process of order \(p\) will be abbreviated to an \(\text{AR}(p)\) process.

oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
set.seed(2021)

auto_reg <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 200)
plot.ts(auto_reg, main = "Autoregressive AR(1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(auto_reg)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(auto_reg)

par(oldpar)

4.3.5 Mixed \(\text{ARMA}\) models

Using the previous knowledge of \(\text{MA}\) and \(\text{AR}\) processes, and their relations, we can combine both into a mixed autoregressive/moving-average process containing \(p\) \(\text{AR}\) terms and \(q\) \(\text{MA}\) terms. This is the \(\text{ARMA}\) process of order \((p, q)\), and it is given by

\[ X_t=\alpha_1X_{t-1}+ \cdots + \alpha_pX_{t-p} + Z_t + \beta_1Z_{t-1}+ \cdots + \beta_qZ_{t-q} \tag{4.3.5.1} \]

Using the backward shift operator \(B\), equation (4.3.5.1) may be written in the form

\[ \phi(B)X_t = \theta(B)Z_t \tag{4.3.5.1a} \]

where \(\phi(B)\), \(\theta(B)\) are polynomials of order \(p\), \(q\) respectively, such that

\[ \phi(B) = 1 – \alpha_1B-\cdots-\alpha_pB^p \]

and

\[ \theta(B) = 1 + \beta_1B+\cdots+\beta_qB^q \]

set.seed(2021)
oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))

stoch_arma <- arima.sim(list(order = c(1,0,1), ma = 0.8, ar = 0.8), n = 200)
plot.ts(stoch_arma, main = "ARMA AR(1) MA(1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(stoch_arma)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(stoch_arma)

par(oldpar)

4.3.6 Integrated \(\text{ARIMA}\) models

In real life, most of the time series we have are non-stationary. This means we have to first remove this source of variation before working with the models we have seen until now, or, we use another composition that already takes in account the non-stationarity. As suggested in section 3.2.3, we can difference the time series to turn it stationary.

Formally we replace \(X_t\) by \(\nabla^d X_t\) where \(d\) is how many times we take the difference (\(\nabla\)) of \(X_t\). This model is called an “integrated” model because the fitted model on the differenced data needs to be summed (or “integrated”) to fit the original data.

Here we define the \(\text{ARIMA}\) model as

\[ W_t = \nabla^d X_t = (1-B)^d X_t \qquad d \in \mathbb{N}_0 \]

the general autoregressive integrated moving average process (abbreviated \(\text{ARIMA}\) process) is of the form

\[ W_t = \alpha_1 W_{t-1} + \cdots + \alpha_p W_{t-p} + Z_t + \cdots + \beta_qZ_{t-q} \tag{4.3.6.1} \]

By analogy with equation (4.3.5.1a), we may write equation (4.3.6.1) in the form

\[ \phi(B)W_t = \theta(B)Z_t \tag{4.3.6.1a} \]

or

\[ \phi(B)(1-B)^dX_t=\theta(B)Z_t \tag{4.3.6.1b} \]

oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
set.seed(2021)

stoch_arma <- arima.sim(list(order = c(2,1,2), ma = c(-0.2279, 0.2488), ar = c(0.8897, -0.4858)), n = 200)
plot.ts(stoch_arma, main = "ARIMA(2, 1, 2)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(stoch_arma)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(stoch_arma)

par(oldpar)

Thus we have an \(\text{ARIMA}\) process of order \((p,d,q)\). The model for \(X_t\) is clearly non-stationary, as the \(\text{AR}\) operator \(\phi(B)(1 – B)^d\) has \(d\) roots on the unit circle. Just for curiosity, see that the random walk process can be modeled in the form of an \(\text{ARIMA}(0, 1, 0)\) process.