Running a decomposition

It is important to note that deseats decomposes a time series under the assumption that the series follows an additive component model. To implement a multiplicative component model, one may use a log-transformation, so that the log-data follows an additive component model.

In the following, we show the application of deseats() to the object NOLABORFORCE, which is already formatted as a "ts" object with frequency = 12, as it is a monthly time series with yearly repeating seasonal pattern. Moreover, the variation around its overall level seems stable enough to assume an additive component model.

yt <- NOLABORFORCE
decomp <- deseats(yt, smoothing_options = set_options())
decomp
#> 
#> ************************************************
#> *                                              *
#> *               Smoothing Results              *
#> *                                              *
#> ************************************************
#> 
#>  Series:    yt
#>  Estimated: Trend + Seasonality
#>  Kernel:    Epanechnikov
#>  Boundary:  Shorten
#>  Bandwidth: 0.2113
#> 
#>  Trend:
#>  ------
#>    Order of local polynomial: 3
#> 
#>  Seasonality:
#>  ------------
#>    Frequency: 12

The console output shows you the main results of the decomposition, including the selected bandwidth. Since deseats() returns an S4-class object, we may also obtain the selected bandwidth as follows:

decomp@bwidth
#> [1] 0.2113423

Aside from the bandwidth, we may obtain the three estimated components, i.e. the fitted trend, the estimated seasonal component, and the residuals, using designated methods, namely trend(), season() and residuals(). To obtain the seasonally adjusted time series, one may call deseasonalize().

trend_est <- trend(decomp)           # Fitted trend
season_est <- season(decomp)         # Fitted seasonal component
error_est <- residuals(decomp)       # Residuals
season_adj <- deseasonalize(decomp)  # Seasonally adjusted series

Illustrating decomposition results

The package contains designated plot() and autoplot() methods to create figures in the style of base R and ggplot2. The user may either select a figure from the options displayed in the R console, or the selection may already be passed to the argument which in the function call. A list of the available figures can be seen from the console after calling plot(decomp) or from the package manual. To select a decomposition plot with the original series and the estimated components, we may run:

plot(decomp, which = 1, xlab = "Year")

Note that the label of the x-axis was adjusted with xlab = "Year" in the function call. Similarly, other options to modify the plot are available. See the documentation for more information.

A more informative plot is available with:

plot(decomp, which = 5, xlab = "Year", ylab = "Millions of person", 
     main = "US persons not in the US labor force", s_around = 55)

It shows the original series together with the fitted trend, and it also displays the estimated seasonality. By default, the estimated seasonal component, even though it is fluctuating around zero, will be shifted and shown around the original series’ sample mean for improved visibility, but using s_around we may re-center the estimated seasonal component within the figure around any arbitrary value.

Forecasting

To produce point and interval forecasts, we employ s_semiarma() first. It runs the decomposition by deseats() and selectes the best ARMA model for the residuals in accordance with either the BIC (default) or the AIC.

full_model <- s_semiarma(yt, smoothing_options = set_options())
full_model
#> 
#> *******************************************
#> *                                         *
#> *        Fitted Seasonal Semi-ARMA        *
#> *                                         *
#> *******************************************
#> 
#> Series: yt
#> 
#> Nonparametric part (Trend + Seasonality):
#> -----------------------------------------
#> 
#>  Kernel:    Epanechnikov
#>  Boundary:  Shorten
#>  Bandwidth: 0.2113
#> 
#>   Trend:
#>   ------
#>     Order of local polynomial: 3
#> 
#>   Seasonality:
#>   ------------
#>     Frequency: 12
#> 
#> Parametric part (Rest):
#> -----------------------
#> 
#> Call:
#> stats::arima(x = res, order = c(ar, 0, ma), include.mean = arma_mean)
#> 
#> Coefficients:
#>          ar1
#>       0.6885
#> s.e.  0.0380
#> 
#> sigma^2 estimated as 0.07362:  log likelihood = -41.54,  aic = 87.08

An AR(1) model was selected for the residuals. Now we can produce the forecasts using predict(). If the residuals resulting from the fitted ARMA model are roughly normal, one may use the setting method = "norm" to obtain forecasting intervals under the assumption of normality. Otherwise, assuming that the errors are at least independent, we may use a bootstrap via method = "boot". If the latter is used, set.seed() should be called to allow for reproducibility. Forecasts are produced as follows: the estimated trend is extrapolated linearly using the last two trend values, the last estimated seasonal cycle is continued into the future, and the residual forecasts are obtained from the fitted ARMA. We condition on the trend and seasonal forecasts, so that the forecasting intervals are only obtained from the ARMA model fitted to the decomposition residuals.

set.seed(1)
fc <- predict(full_model, n.ahead = 12, method = "boot")

The returned forecasting object contains the point forecasts (fc@pred) as well as the lower and the upper bounds of the forecasting intervals (fc@interv). By default, those bounds will be obtained at the 95% and the 99% level; however, this can be adjusted with the argument alpha in the call to predict(). A simple forecasting plot may now be created as follows:

plot(fc, xlab = "Year", ylab = "Millions of person", main = "US persons not in the US labor force",
     xlim = c(2010, 2021 - 1 / 12), ylim = c(80, 97))

xlim and ylim were used to adjust the plot window, so that only the last observations will be shown together with the forecasts. If we had implemented an initial log-transformation of the data, we may want to retransform the forecasts. In such a case, we should set expo = TRUE in the call to the predict() method. By default, with adjust.bias = TRUE, the point forecasts are then adjusted for biases on the original scale. Under normality, this is done using the theoretical relation between the normal distribution (for the log-data forecasts) and the log-normal distribution (for the original data forecasts). When the bootstrap is applied, the sample mean of the future retransformed paths is considered.

Note: forecasts using this method should only be created for the next few future time points, as extrapolating estimated components can lead to highly unreliable results, if going too far into the future.