The following code provides examples on how to decompose a time series comprised of a deterministic trend component, a deterministic seasonal component and a (non-deterministic) stationary component, which could be either a WN, an ARMA or a similar kind of process.

Looking back on the White Noise process

Often it is important to establish whether a stationary process under consideration is WN. To achieve this, we can use the graphs of the ACF and PACF functions produced by tsdisplay:

suppressPackageStartupMessages({library(forecast)})
set.seed(50)
rn1 = ts(rnorm(100))
rn2 = ts(rnorm(100))
tsdisplay(rn1)

tsdisplay(rn2)

while both processes were generated from the normal distribution, we can see that the correlogram of the second sample is ambiguous - the bars of ACF and PACF are close or slightly outside the confidence band.

We can create a simple model:

\[\text{rn1}_t = \beta_0 + \epsilon_t\] and test if the residuals are WN. If \(\epsilon_t\) is a zero-mean WN, then \(\beta_0\) would correspond to the mean of the WN process. We can carry out the Ljung-Box test graphically:

mod1 = arima(rn1, c(0,0,0))
tsdiag(mod1)

All the bubbles (p-values) in the bottom graphs are above the critical 0.05 blue line, thus rn1 is white noise.

Trend Decomposition via Moving Average

We will provide an example of the moving-average method for decomposition.

If the time series contains a seasonal component and we want to average it out, the length of the moving average must be equal to the seasonal frequency (for monthly series, we would take \(l=12\)). However, there is a slight hurdle. Suppose, our time series begins in January (\(t=1\)) and we average up to December (\(t=12\)). This averages corresponds to a time \(t = 6.5\) (time between June and July). When we come to estimate seasonal effects, we need a moving average at integer times. This can be achieved by averaging the average of January to December and the average of February (\(t= 2\)) up to January (\(t=13\)). This average if the two moving averages corresponds to \(t=7\) and the process is called centering.

Thus, the trend at time \(t\) can be estimated by the centered moving average.

Two-sided averaging:

\[ \begin{aligned} \widehat{T}_t &= \dfrac{(Y_{t-6} +...+ Y_{t+5})/12 + (Y_{t-5} +...+ Y_{t+6})/12}{2} = \dfrac{(1/2)Y_{t-6} + Y_{t-5}...+ Y_{t+5}+ (1/2)Y_{t+6}}{12} \end{aligned} \]

One-sided averaging:

\[ \begin{aligned} \widehat{T}_t &= \dfrac{(Y_{t-12} +...+ Y_{t-1})/12 + (Y_{t-11} +...+ Y_{t})/12}{2} = \dfrac{(1/2)Y_{t-12} + Y_{t-11}...+ Y_{t-1}+ (1/2)Y_{t}}{12} \end{aligned} \]

Note that the weights are different for the first and last elements used in the averaging. One way to perform the moving-average decomposition in R is to use the filter function:

suppressPackageStartupMessages({library(fma)})
data(airpass)
lat <- log(airpass)
#one-sided averaging:
lat.trend.s1 <- filter(lat, filter = c(1/24, rep(1/12, 11), 1/24), sides = 1)
#two-sided averaging:
lat.trend.s2 <- filter(lat, filter = c(1/24, rep(1/12, 11), 1/24), sides = 2)

Alternatively, the same can be achieved by computing the formulas from above:

#Two-sided:
T.hat.s1 <- rep(NA, length(lat))
#one-sided:
T.hat.s2 <- rep(NA, length(lat))
for(j in 1:length(lat)){
  if(j <= 6 | j > (length(lat)-6)){
    T.hat.s2[j] <- NA
  }else{
    T.hat.s2[j] <- sum(lat[(j+(-6:6))] * c(1/24, rep(1/12, 11), 1/24))
  }
  if(j <= 12){
    T.hat.s1[j] <- NA
  }else{
    T.hat.s1[j] <- sum(lat[(j+(-12:0))] * c(1/24, rep(1/12, 11), 1/24))
  }
}
T.hat.s1 <- ts(T.hat.s1, start = start(lat), frequency = frequency(lat))
T.hat.s2 <- ts(T.hat.s2, start = start(lat), frequency = frequency(lat))

The total sum of the absolute differences between the values estimated via the build-in function and the ones calculated manually are very small:

#one-sided average calculation comparisons
sum(abs(T.hat.s1 - lat.trend.s1), na.rm = TRUE)

## [1] 6.039613e-14

#two-sided average calculation comparisons
sum(abs(T.hat.s2 - lat.trend.s2), na.rm = TRUE)

## [1] 6.039613e-14

#plot the data
plot(lat, type = "l", main = "Moving-average trend decomposition")
lines(lat.trend.s1, col = 2)
lines(lat.trend.s2, col = 4, lty = 2, lwd = 2)
legend(x = 1949, y = 6.4,
       legend = c("Actual", "One-sided", "Two-sided "),
       col = c(1, 2, 4),
       lty = c(1, 1, 2))

Seasonal Effect Decomposition

Having estimated the trend, we can move on to estimating the seasonal effect.

We will decompose the seasonal effect using both trend averaging methods and examine how the decomposition fits the data:

Estimate the trend component \(T_t\) using the moving average method.
Estimate \(\widehat{S}_t = Y_t - \widehat{T_t}\), \(t = 1,...,T\);
Decide on the seasonal period \(d\): \(S_t = S_{t+d}\) and \(\sum_{t = 1}^d S_t = 0\). Then, calculate \(\overline{S}_k = \text{mean}\left( \widehat{S}_k, \widehat{S}_{2k}, \widehat{S}_{3k}, ... \right)\), \(k = 1,...,d\);
Calculate the total sum: \(S_{tot} = \sum_{k = 1}^d \overline{S}_k\);
Calculate \(w = \dfrac{1}{d} \sum_{k = 1}^d S_{tot}\)
Adjust each seasonal factor: \(\tilde{S}_t = \widehat{S}_t - w\)
We now have estimated the seasonal component \(\tilde{S}_t\) such that \(\tilde{S}_t = \tilde{S}_{t+d}\) and \(\sum_{t= 1}^d \tilde{S}_t = 0\).

#seasonal component of one-sided averaging
lat.sez1   <- lat - lat.trend.s1
lat.se1    <- matrix(lat.sez1, ncol = 12, byrow = TRUE)
lat.s1     <- apply(lat.se1, 2, mean, na.rm = TRUE)
#seasonal component of two-sided averaging
lat.sez2   <- lat - lat.trend.s2
lat.se2    <- matrix(lat.sez2, ncol = 12, byrow = TRUE)
lat.s2     <- apply(lat.se2, 2, mean, na.rm = TRUE)

plot(lat, lwd = 2, main = "One-sided averaging")
lines(lat.trend.s1, col = 2)
lines(lat.trend.s1 + lat.s1, col = 3)
legend(1949, 6.4, c("actual", "trend", "trend + seas "), lty = 1, col = c(1, 2, 3))

plot(lat, lwd = 2, main = "Two-sided averaging")
lines(lat.trend.s2, col = 2)
lines(lat.trend.s2 + lat.s2, col = 3)
legend(1949, 6.4, c("actual", "trend", "trend + seas "), lty = 1, col = c(1, 2, 3))

The decomposed trend and seasonality fit the data much nicer compared to only the trend component. We can also look at what the seasonality component looks like on its own:

par(mfrow = c(1,2))
plot(1:12, lat.s1, type = "l", main = "Seasonal (one-sided decomp.)")
plot(1:12, lat.s2, type = "l", main = "Seasonal (two-sided decomp.)")

If needed, we can also adjust the seasonal components, so that they sum up to zero:

adj   <- sum(lat.s2)/length(lat.s2)
adj

## [1] -0.00086568

The adjustment is very small. In this case, the adjustment wasn’t needed. The plots of the adjusted seasonal component are identical to the unadjusted decomposition:

lat.adj.s2 <- lat.s2 - adj
par(bg = "grey60")
plot(lat, lwd = 2, col = 1)
lines(lat.trend.s2 + lat.s2, col = 2)
lines(lat.trend.s2 + lat.adj.s2, col = 3, lty = 2)
legend(x = 1949, y= 6.4, 
       legend = c("actual", "trend + seas.", "trend + adj. seas."),
       lwd = c(2, 1, 1),
       lty = c(1, 1, 2), 
       col = c(1, 2, 3))

Finally, we can check if the remainder is stationary:

resid.s1 <- lat - lat.trend.s1 - lat.s1
resid.s2 <- lat - lat.trend.s2 - lat.s2
par(mfrow = c(1,2)); plot.ts(resid.s1); plot.ts(resid.s2)

par(mfrow = c(1, 2))
Acf(resid.s1, main = "Residuals of one-sided decomp.")
Acf(resid.s2, main = "Residuals of two-sided decomp.")

Pacf(resid.s1, main = "Residuals of one-sided decomp.")
Pacf(resid.s2, main = "Residuals of two-sided decomp.")

The series is not white noise, so we could estimate an ARMA model on the residuals. Although, it seems that the two-sided decomposition still leaves a seasonality effect as shown by the oscillating residual ACF plot.

Note that this method does not allow to produce forecasts. However, this is possible when using the methods introduced below.

Note that it is also possible to decompose a time series into seasonal, trend and irregular components using moving-averages via the decompose function from forecast package.

Decomposition using `decompose` function

Say we want to decompose the GDP of Lithuania into trend and seasonal components and also a random component (which is hopefully stationary). While many functions in R can decompose the data, only some allow to forecast the series.

suppressPackageStartupMessages({require("forecast")})
txt1 <- "http://uosis.mif.vu.lt/~rlapinskas/(data%20R&GRETL/"
txt2 <- "L.gdp.txt"
l.gdp <- read.csv(url(paste0(txt1, txt2)), 
                     header = TRUE, as.is = TRUE, skip = 2)
l.gdp  <- ts(l.gdp, start = c(1993, 1), frequency = 4)
ll.gdp <- log(l.gdp)
ll <- window(ll.gdp, start = 2000)
par(mfrow = c(1,2))
plot(l.gdp, main = "Lithuanian GDP")
plot(ll, main = "log of Lithuanian GDP")

We restrict ourselves to 2000:01-2012:04 period, because of the transition to the free market economy (irregular up to 1996) and the Russian financial crisis in 1998. We choose to analyse logarithms because its seasonal part appears to be of constant amplitude. thus, we want to find the additive decomposition: \(log(l.gdp_t) = T_t + S_t + \epsilon_t\).

Looking at the graphs, it would be difficult to describe the trend in analytic terms, therefore, we shall apply the decompose function to use non-parametric decomposition:

dec.ll <- decompose(ll)
plot(dec.ll)

tsdisplay(dec.ll$rand)

The random component of dec.ll seems to be stationary, but not a white noise. However, the main problem with our procedure is that it does not allow to forecast ll.

The Global Method of OLS

The OLS method estimates coefficients of an additive model. If necessary, we can estimate a polynomial trend. We will use the monthly airline passenger numbers 1949-1960.

plot(airpass)

#time sequence of dates:
at = time(airpass)
#time sequence + 36 months for forecast:
time.f <- seq(1949, 1964-1/12, by = 1/12)
#seasonal dummy variables:
month   <- seasonaldummy(airpass)
#seasonal dummy forecast:
month.f <- seasonaldummy(airpass, h = 36)

We will estimate three different models via OLS:

Model 1: \(AP_t = \beta_0 + \beta_1 t + \beta_2 t^2 + \epsilon_t\);
Model 2: \(AP_t = \beta_0 + \beta_1 t + \beta_2 t^2 + \gamma_1 dm1_t +...+\gamma_{11} dm11_t + \epsilon_t\);
Model 3: \(log \left( AP_t \right) = \beta_0 + \beta_1 t + \beta_2 t^2 + \gamma_1 dm1_t +...+\gamma_{11} dm11_t + \epsilon_t\);

#only quadratic trend:
air1.lm <- lm(airpass ~ at + I(at^2))
summary(air1.lm)

## 
## Call:
## lm(formula = airpass ~ at + I(at^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -100.353  -27.339   -7.442   21.603  146.116 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  3.795e+06  1.333e+06   2.848  0.00506 **
## at          -3.914e+03  1.363e+03  -2.871  0.00473 **
## I(at^2)      1.009e+00  3.487e-01   2.894  0.00441 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 44.91 on 141 degrees of freedom
## Multiple R-squared:  0.8618, Adjusted R-squared:  0.8599 
## F-statistic: 439.8 on 2 and 141 DF,  p-value: < 2.2e-16

air1.f <- predict(air1.lm, data.frame(at = time.f))

#quadratic trend and monthly seasonality:
air2.lm <- lm(airpass ~ at + I(at^2) + month)
summary(air2.lm)

## 
## Call:
## lm(formula = airpass ~ at + I(at^2) + month)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.122 -13.394   0.825  12.733  75.773 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.866e+06  7.047e+05   5.485 2.08e-07 ***
## at          -3.986e+03  7.210e+02  -5.529 1.70e-07 ***
## I(at^2)      1.028e+00  1.844e-01   5.573 1.38e-07 ***
## monthJan     9.180e+00  9.709e+00   0.946  0.34612    
## monthFeb    -1.587e-01  9.706e+00  -0.016  0.98698    
## monthMar     3.240e+01  9.704e+00   3.339  0.00110 ** 
## monthApr     2.670e+01  9.702e+00   2.752  0.00676 ** 
## monthMay     2.882e+01  9.700e+00   2.971  0.00353 ** 
## monthJun     6.601e+01  9.699e+00   6.806 3.33e-10 ***
## monthJul     1.030e+02  9.697e+00  10.623  < 2e-16 ***
## monthAug     1.001e+02  9.696e+00  10.323  < 2e-16 ***
## monthSep     4.874e+01  9.695e+00   5.027 1.62e-06 ***
## monthOct     1.020e+01  9.695e+00   1.052  0.29475    
## monthNov    -2.627e+01  9.694e+00  -2.710  0.00764 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.75 on 130 degrees of freedom
## Multiple R-squared:  0.9644, Adjusted R-squared:  0.9608 
## F-statistic: 270.8 on 13 and 130 DF,  p-value: < 2.2e-16

air2.f <- predict(air2.lm, data.frame(at = time.f,
                                      month = I(rbind(month, month.f))))

#model for logs:
air3.lm <- lm(log(airpass) ~ at + I(at^2) + month)
summary(air3.lm)

## 
## Call:
## lm(formula = log(airpass) ~ at + I(at^2) + month)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12748 -0.03709  0.00418  0.03197  0.11529 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.205e+04  1.430e+03  -8.426 5.79e-14 ***
## at           1.222e+01  1.463e+00   8.347 8.95e-14 ***
## I(at^2)     -3.093e-03  3.743e-04  -8.265 1.41e-13 ***
## monthJan     2.132e-02  1.971e-02   1.082    0.281    
## monthFeb    -9.486e-04  1.970e-02  -0.048    0.962    
## monthMar     1.291e-01  1.970e-02   6.555 1.19e-09 ***
## monthApr     9.771e-02  1.969e-02   4.962 2.15e-06 ***
## monthMay     9.525e-02  1.969e-02   4.838 3.66e-06 ***
## monthJun     2.174e-01  1.969e-02  11.041  < 2e-16 ***
## monthJul     3.213e-01  1.968e-02  16.323  < 2e-16 ***
## monthAug     3.120e-01  1.968e-02  15.855  < 2e-16 ***
## monthSep     1.675e-01  1.968e-02   8.511 3.62e-14 ***
## monthOct     2.947e-02  1.968e-02   1.497    0.137    
## monthNov    -1.141e-01  1.968e-02  -5.797 4.84e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0482 on 130 degrees of freedom
## Multiple R-squared:  0.9892, Adjusted R-squared:  0.9881 
## F-statistic: 912.7 on 13 and 130 DF,  p-value: < 2.2e-16

air3.f <- predict(air3.lm, data.frame(at = time.f,
                                      month = I(rbind(month, month.f))))

We can view the estimated model forecasts:

par(mfrow = c(1,3))
plot(time.f, air1.f, type = "l", lty = 2, 
     main = "quadratic")
lines(airpass, col = 2)
plot(time.f, air2.f, type = "l", lty = 2, 
     main = "quadratic + seasonal")
lines(airpass, col = 2)
plot(time.f, air3.f, type = "l", lty = 2,
     main = "quadratic + seasonal for \n log(airpass)")
lines(log(airpass), col = 2)

While the central graph attempts to capture the seasonality effect and it does so better than the left-side graph, we can see that the actual data fluctuations increase in size together with the level of airpass. so an additive model is unable to catch this effect. The common approach to this problem is to take the logarithms and estimate the additive model (right-side graph).

The Global Methods of Nonlinear Least Squares

If an economic indicator increases every year by the same percent, then it grows exponentially. To model this kind of time series, an exponential model is needed. however, modeling is rather complicated because to extract the trend here we need a nonlinear regression (nls function).

suppressPackageStartupMessages({library(fma)})
data(shampoo)
TIME = as.numeric(time(shampoo))
#nonlinear least squares with starting values for a,b , c
shampoo.nls <- nls(shampoo ~ a + b*exp(c*TIME), start = list(a = 100, b = 10, c = 1))
plot(shampoo, main = "Shampoo sales data")
lines(TIME, predict(shampoo.nls))

Let us also check the residuals:

#The residuals are not normal:
par(mfrow = c(1,2))
plot(summary(shampoo.nls)$res, main = "Residuals")
abline(0, 0, col = "red")
hist(summary(shampoo.nls)$res, main = "Histogram of Residuals")

The spread of the residuals is more or less constant, however, they do not appear normal.

One of the drawbacks of the OLS is its globality - if the time series departs for some reason from the trend at the final moment \(T\), this will change all the coefficients of the model and therefore the predicted value at moment \(t = 1\). The local methods are free of this defect - their forecast at moment \(t\) is (mostly) influenced by the values at close time moments.

The Local Method of Exponential Smoothing

We know three variants of exponential smoothing - they are for a series without trend, with a linear trend and with seasonality.

R allows to use an automated procedure to choose the smoothing parameters - the function ets selects the best variant automatically (by the minimum of AIC) and estimates its parameters:

suppressPackageStartupMessages({
  library(forecast)
  library(datasets)
})
data("AirPassengers")
l.ap.ets <- ets(log(AirPassengers))
summary(l.ap.ets)

## ETS(A,A,A) 
## 
## Call:
##  ets(y = log(AirPassengers)) 
## 
##   Smoothing parameters:
##     alpha = 0.6534 
##     beta  = 1e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 4.8022 
##     b = 0.01 
##     s=-0.1047 -0.2186 -0.0761 0.0636 0.2083 0.217
##            0.1145 -0.011 -0.0111 0.0196 -0.1111 -0.0905
## 
##   sigma:  0.0359
## 
##       AIC      AICc       BIC 
## -208.3619 -203.5047 -157.8750 
## 
## Training set error measures:
##                         ME       RMSE        MAE         MPE     MAPE
## Training set -0.0006710596 0.03592072 0.02773886 -0.01250262 0.508256
##                   MASE       ACF1
## Training set 0.2291672 0.09431354

plot(forecast(l.ap.ets), include = 36)

tsdisplay(l.ap.ets$residuals)

It seems that the seasonality is removed.

Looking at the log of GDP of Lithuania

ets.ll = ets(ll)
plot(forecast(ets.ll), include = 16, 
     main = "Exponential smoothing forecast of log of Lithuanian GDP")

tsdisplay(ets.ll$resid)

It seems that this decomposition method worked quite well, although it doesn’t provide much in terms of coefficient interpretation and the point-forecast (the blue line) seems to be less affected by the trend, compared to the historical data.

Local Linear Forecast Using Cubic Splines

The cubic smoothing spline model is equivalent to an ARIMA(0, 2, 2) (this model will be presented later) but with a restricted parameter space. The advantage of the spline model over the full ARIMA model is that it provides a smooth historical trend as well as a linear forecast function.

We shall use the splinef from the forecast package to extract the linear forecast from shampoo time series and forecast it for 24 months:

suppressPackageStartupMessages({
  library(forecast)
  library(fma)
})
data(shampoo)
par(mfrow=c(1,2))
fcast <- splinef(shampoo, h = 24)
plot(fcast, 
     main = "Cubic smoothing spline of sales of shampoo over a three year period")
fcast.l <- splinef(log(shampoo), h = 24)
plot(fcast.l, 
     main = "Cubic smoothing spline of sales of log of shampoo over a three year period")

STL decomposition

STL - Seasonal and Trend decomposition using Loess.

It is a very versatile and robust decomposition method - the seasonal component is allowed to change over time (the rate of change can be controlled by the user), the smoothness of the trend cycle can also be controlled by the user.

This method is robust to outliers, however it is only available for additive decompositions.

suppressPackageStartupMessages(library(fpp))
data(elecequip)
plot.ts(elecequip)

par(mfrow = c(1,2))
Acf(elecequip)
Pacf(elecequip)

par(mfrow = c(1,1))
fit <- stl(elecequip, t.window = 15, s.window = "periodic", robust=TRUE)
plot(fit)

here t.window controls the ‘wiggliness’ of the trend components and s.window controls the variation on the seasonal component.

resid <- fit$time.series[, "remainder"]
tsdisplay(resid)

Box.test(resid, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  resid
## X-squared = 41.109, df = 20, p-value = 0.003607

Box.test(resid, lag = 20, type = "Box-Pierce")

## 
##  Box-Pierce test
## 
## data:  resid
## X-squared = 38.758, df = 20, p-value = 0.007144

The remainder seems to not have a trend and seasonal component anymore, though it does not appear to be a white noise process - using Ljung-Box and Box-Pierce tests, we reject the null hypothesis (p-value < 0.05) of no residual correlation in the first 20 lags.

Task 02