knitr::opts_chunk$set(cache = FALSE, fig.show = 'hold', fig.width = 7, fig.height = 3) library(stR) library(forecast) library(seasonal)
This vignette describes some functionality and provides a few examples of how to use package stR. stR implements method STR, where STR stands for Seasonal-Trend decomposition by Regression and capital R emphasizes R, the name of the popular statistical software. The whole name is also reminiscent of STL, the method which inspired me to introduce STR.
There are many packages and methods which work with seasonal data. For example the oldest method for decomposition -- classical additive decomposition -- is implemented in package stats
. The method splits the data into trend, seasonal and random components:
m <- decompose(co2) plot(m)
Another well known method is STL, implemented in packages stats
and stlplus
:
plot(stl(log(co2), s.window = "periodic", t.window = 30))
Other R packages which implement various versions of seasonal decomposition and seasonal adjustment include forecast
, x12
, seasonal
, season
, seas
and deseasonalize
.
A few more examples are provided below:
f = "co2.fit.tbats.RDS" if(file.exists(f)) { co2.fit = readRDS(f) } else { co2.fit <- tbats(co2) saveRDS(co2.fit, f, compress = "xz") }
library(forecast) co2.fit <- tbats(co2)
plot(co2.fit)
library(seasonal) co2.fit <- seas(co2) plot(co2.fit, trend = TRUE)
After looking at the above examples, a reader might ask "Why do we need another method of seasonal decomposition?"
A short answer is that
This vignette provides some details on the first claim.
For time series decomposition with objects of class ts
or class msts
, and with no regressors or complex seasonality, it is simple to do an STR decomposition using the AutoSTR
function.
For example, the co2
time series can be decomposed as follows:
f = "co2.fit.stR.RDS" if(file.exists(f)) { co2.fit = readRDS(f) } else { co2.fit <- AutoSTR(co2) saveRDS(co2.fit, f, compress = "xz") }
co2.fit <- AutoSTR(co2)
plot(co2.fit)
The time series taylor
from package forecast
provides us with half-hourly electricity demand in England and Wales. It exhibits (at least) two seasonalities -- daily and weekly. They can be observed in a 4 weeks' subset of the data:
taylor.msts <- msts(log(head(as.vector(taylor), 336*4)), seasonal.periods=c(48,48*7,48*7*52.25), start=2000+22/52) plot(taylor.msts, ylab = "Electricity demand")
Since the data is at half hour granularity, the daily seasonality has a period of 48 observations and weekly has a period of 336.
The data is of class msts
(multiple seasonal time series), which can also be handled automatically with AutoSTR
:
f = "taylor.fit.stR.RDS" if(file.exists(f)) { taylor.fit = readRDS(f) } else { taylor.fit <- AutoSTR(taylor.msts, gapCV = 48, confidence = 0.95) saveRDS(taylor.fit, f, compress = "xz") }
taylor.fit <- AutoSTR(taylor.msts, gapCV = 48, confidence = 0.95)
plot(taylor.fit)
The parameters supplied to AutoSTR
are:
gapCV = 48
-- gaps of 48 observations are used for cross validationreltol = 0.001
-- this parameter is passed directly to optim
function to control how well (and for how long) the model parameters are optimizedconfidence = 0.95
-- 95% confidence intervals are calculated (with assumptions: errors are uncorrelated, model parameters are estimated exactly)This example shows how to tune an STR decomposition, rather than use the automated AutoSTR
function. STR is a flexible method, which can be adjusted to the data in multiple ways, making the interface rather complex.
Let us consider the dataset grocery
which contains monthly data of supermarkets and grocery stores turnover in New South Wales:
plot(grocery, ylab = "NSW Grocery Turnover, $ 10^6")
We will use a log transformation to stabilize seasonal variance in the data:
logGr <- log(grocery) plot(logGr, ylab = "NSW Grocery Turnover, log scale")
At the next step we define trend and seasonal structures.
Trend does not have seasonality, therefore its seasonal structure contains only a single knot: c(1,0)
. Here a knot is defined as a pair of numbers, both referring to the same point. The segments
component contains only one segment c(0,1)
.
trendSeasonalStructure <- list( segments = list(c(0,1)), sKnots = list(c(1,0)))
The seasonal structure of the data is defined in the seasonalStructure
variable. The segments
component contains one pair of numbers c(0,12)
which defines ends the segment, while sKnots
variable contains seasonal knots from 1 to 12 (months). The last knot c(0,12)
also defines that ends of the segment c(0,12)
are connected and 0 and 12 represent the same knot (month December).
seasonalStructure <- list( segments = list(c(0,12)), sKnots = list(1,2,3,4,5,6,7,8,9,10,11,c(12,0)))
Variable seasons
contains months corresponding to every data point:
seasons <- as.vector(cycle(logGr))
Since trend does not have seasonality (the seasonal structure for the trend contains only one season), the trend seasons are all ones:
trendSeasons <- rep(1, length(logGr))
The times
vector contains times corresponding to data points:
times <- as.vector(time(logGr))
The data
vector contains observations (in this case log turnover):
data <- as.vector(logGr)
trendTimeKnots
vector contains times where time knots for the trend are positioned:
trendTimeKnots <- seq( from = head(times, 1), to = tail(times, 1), length.out = 175)
The seasonTimeKnots
vector contains times where time knots for the seasonal component are positioned:
seasonTimeKnots <- seq( from = head(times, 1), to = tail(times, 1), length.out = 15)
In stR package every component of a decomposition is a regressor. In the case of trend and seasonal components values for such regressors are constants (vectors of ones):
trendData <- rep(1, length(logGr)) seasonData <- rep(1, length(logGr))
The complete trend
structure contains all components relevant to the trend.
Component lambdas
is always a vector with three elements. Trend, since it does not have seasonality, has only the first element different from zero. This element defines smoothness of the trend at the starting point of the optimization procedure.
trend <- list( name = "Trend", data = trendData, times = times, seasons = trendSeasons, timeKnots = trendTimeKnots, seasonalStructure = trendSeasonalStructure, lambdas = c(0.5,0,0))
The complete season
structure contains all components relevant to the trend. Component lambdas
is a vector with three elements. The elements define smoothness of the seasonal component at starting point of the optimization procedure.
According to STR approach every non-static seasonal component has two-dimensional structure (see the corresponding article on STR method). In this particular case it has topology of a tube. The "observed" seasonal component is a spiral subset of that "tube organised" data ("observed" means here that the seasonal component is observed as part of the data together with the trend and other components).
The first element of the vector defines smoothness of data along time dimension of the tube. The second component defines smoothness along seasonal dimension. And the third component defines smoothness in some way in both dimensions (by restricting partial discrete derivative in both directions).
The last two zeros in lambdas
component mean that those two components will not be optimized (and effectively two dimensional structure of the seasonal component will not be used).
season <- list( name = "Yearly seasonality", data = seasonData, times = times, seasons = seasons, timeKnots = seasonTimeKnots, seasonalStructure = seasonalStructure, lambdas = c(10,0,0))
All components of STR decomposition are considered to be predictors. For example trend is a predictor with no seasonality and independent variable which is constant (and equal to one). The seasonal component is a predictor with some predefined seasonality and a constant independent variable (also equal to one).
predictors <- list(trend, season)
To calculate STR decomposition we supply data points, predictors, required confidence intervals, the gap, to perform cross validation, and reltol parameter which was described earlier.
f = "gr.fit.stR.RDS" if(file.exists(f)) { gr.fit = readRDS(f) } else { gr.fit <- STR(data, predictors, confidence = 0.95, gap = 1, reltol = 0.00001) saveRDS(gr.fit, f, compress = "xz") }
gr.fit <- STR(data, predictors, confidence = 0.95, gap = 1, reltol = 0.00001)
plot(gr.fit, xTime = times, forecastPanels = NULL)
In plot
function forecastPanels = NULL
means that fit/forecast are not displayed.
If we decide to use two-dimensional structure for the seasonal component we need to redefine lambdas
component:
season <- list(name = "Yearly seasonality", data = seasonData, times = times, seasons = seasons, timeKnots = seasonTimeKnots, seasonalStructure = seasonalStructure, lambdas = c(1,1,1)) predictors = list(trend, season)
f = "gr.fit.2.stR.RDS" if(file.exists(f)) { gr.fit = readRDS(f) } else { gr.fit <- STR(data, predictors, confidence = 0.95, gap = 1, reltol = 0.00001) saveRDS(gr.fit, f, compress = "xz") }
gr.fit <- STR(data, predictors, confidence = 0.95, gap = 1, reltol = 0.00001)
This allows to find a set of parameters with a smaller cross validated error, and therefore potentially more insightful decomposition:
plot(gr.fit, xTime = times, forecastPanels = NULL)
Since STR is based on Ordinary Least Squares (OLS), it does not tolerate outliers well. In the examples below I compare STR with its robust version -- Robust STR. The latter is based on quantile regression approach (only 0.5 quantile is used), and therefore is robust to various types of outliers.
Let us create a time series with two "spikes" to model two isolated outliers:
outl <- rep(0,length(grocery)) outl[14] <- 900 outl[113] <- -700 tsOutl <- ts(outl, start = c(2000,1), frequency = 12)
and combine it with grocery
time series
logGrOutl <- log(grocery + tsOutl) plot(logGrOutl, ylab = "Log turnover with outliers")
Decomposition of this time series with STR shows considerable distortions in all components:
trendSeasonalStructure <- list(segments = list(c(0,1)), sKnots = list(c(1,0))) seasonalStructure <- list(segments = list(c(0,12)), sKnots = list(1,2,3,4,5,6,7,8,9,10,11,c(12,0))) seasons <- as.vector(cycle(logGrOutl)) trendSeasons <- rep(1, length(logGrOutl)) times <- as.vector(time(logGrOutl)) data <- as.vector(logGrOutl) timeKnots <- times trendData <- rep(1, length(logGrOutl)) seasonData <- rep(1, length(logGrOutl)) trend <- list(data = trendData, times = times, seasons = trendSeasons, timeKnots = timeKnots, seasonalStructure = trendSeasonalStructure, lambdas = c(0.1,0,0)) season <- list(data = seasonData, times = times, seasons = seasons, timeKnots = timeKnots, seasonalStructure = seasonalStructure, lambdas = c(10,0,0)) predictors <- list(trend, season)
f = "logGrOutl.stR.RDS" if(file.exists(f)) { fit.str = readRDS(f) } else { fit.str <- STR(as.vector(logGrOutl), predictors, confidence = 0.95, gapCV = 1, reltol = 0.001) saveRDS(fit.str, f, compress = "xz") }
fit.str <- STR(as.vector(logGrOutl), predictors, confidence = 0.95, gapCV = 1, reltol = 0.001)
plot(fit.str, xTime = times, forecastPanels = NULL)
On the other hand Robust STR decomposition results in much cleaner decomposition. The outliers appear only as residuals (component with name "Random") and trend and seasonal components are not distorted:
f = "logGrOutl.stR.robust.RDS" if(file.exists(f)) { fit.rstr = readRDS(f) } else { fit.rstr <- STR(as.vector(logGrOutl), predictors, confidence = 0.95, gapCV = 1, reltol = 0.001, nMCIter = 200, robust = TRUE) saveRDS(fit.rstr, f, compress = "xz") }
fit.rstr <- STR(as.vector(logGrOutl), predictors, confidence = 0.95, gapCV = 1, reltol = 0.001, nMCIter = 200, robust = TRUE)
plot(fit.rstr, xTime = times, forecastPanels = NULL)
Data set calls
provides data about number of call arrivals per 5-minute interval handled on weekdays between 7:00 am and 9:05 pm from March 3, 2003 in a large North American commercial bank.
Below is an example of decomposition of calls
data using STR:
times <- as.vector(time(calls)) timeKnots <- seq(min(times), max(times), length.out=25) trendSeasonalStructure <- list(segments = list(c(0,1)), sKnots = list(c(1,0))) trendSeasons <- rep(1, length(calls)) sKnotsDays <- as.list(seq(1, 169, length.out = 169)) seasonalStructureDays <- list(segments = list(c(1, 169)), sKnots = sKnotsDays) seasonsDays <- seq_along(calls) %% 169 + 1 sKnotsWeeks <- as.list(seq(0,169*5, length.out=13*5)) seasonalStructureWeeks <- list(segments = list(c(0, 169*5)), sKnots = sKnotsWeeks) seasonsWeeks <- seq_along(calls) %% (169*5) + 1 data <- as.vector(calls) trendData <- rep(1, length(calls)) seasonData <- rep(1, length(calls)) trend <- list(data = trendData, times = times, seasons = trendSeasons, timeKnots = timeKnots, seasonalStructure = trendSeasonalStructure, lambdas = c(0.02, 0, 0)) seasonDays <- list(data = seasonData, times = times, seasons = seasonsDays, timeKnots = seq(min(times), max(times), length.out=25), seasonalStructure = seasonalStructureDays, lambdas = c(0, 11, 30)) seasonWeeks <- list(data = seasonData, times = times, seasons = seasonsWeeks, timeKnots = seq(min(times), max(times), length.out=25), seasonalStructure = seasonalStructureWeeks, lambdas = c(30, 500, 0.02)) predictors <- list(trend, seasonDays, seasonWeeks)
f = "calls.fit.RDS" if(file.exists(f)) { calls.fit = readRDS(f) } else { calls.fit <- STR(data = data, predictors = predictors, confidence = 0.95, reltol = 0.003, nFold = 4, gap = 169) saveRDS(calls.fit, f, compress = "xz") }
calls.fit <- STR(data = data, predictors = predictors, confidence = 0.95, reltol = 0.003, nFold = 4, gap = 169)
plot(calls.fit, xTime = as.Date("2003-03-03") + ((seq_along(data)-1)/169) + (((seq_along(data)-1)/169) / 5)*2, forecastPanels = NULL)
Electricity consumption dataset electricity
provides information about electricity consumption in Victoria, Australia during the 115 days starting on 10th of January, 2000, and comprises the maximum electricity demand in Victoria during 30-minute periods (48 observations per day). For each 30-minute period, the dataset also provides the air temperature in Melbourne.
In the example below the data is decomposed using weekly seasonal pattern, daily seasonal pattern which takes into account weekends and holidays and transition periods between them, and two flexible predictors:
TrendSeasonalStructure <- list(segments = list(c(0,1)), sKnots = list(c(1,0))) DailySeasonalStructure <- list(segments = list(c(0,48)), sKnots = c(as.list(1:47), list(c(48,0)))) WeeklySeasonalStructure <- list(segments = list(c(0,336)), sKnots = c(as.list(seq(4,332,4)), list(c(336,0)))) WDSeasonalStructure <- list(segments = list(c(0,48), c(100,148)), sKnots = c(as.list(c(1:47,101:147)), list(c(0,48,100,148)))) TrendSeasons <- rep(1, nrow(electricity)) DailySeasons <- as.vector(electricity[,"DailySeasonality"]) WeeklySeasons <- as.vector(electricity[,"WeeklySeasonality"]) WDSeasons <- as.vector(electricity[,"WorkingDaySeasonality"]) Data <- as.vector(electricity[,"Consumption"]) Times <- as.vector(electricity[,"Time"]) TempM <- as.vector(electricity[,"Temperature"]) TempM2 <- TempM^2 TrendTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 116) SeasonTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 24) SeasonTimeKnots2 <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 12) TrendData <- rep(1, length(Times)) SeasonData <- rep(1, length(Times)) Trend <- list(name = "Trend", data = TrendData, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(1500,0,0)) WSeason <- list(name = "Weekly seas", data = SeasonData, times = Times, seasons = WeeklySeasons, timeKnots = SeasonTimeKnots2, seasonalStructure = WeeklySeasonalStructure, lambdas = c(0.8,0.6,100)) WDSeason <- list(name = "Daily seas", data = SeasonData, times = Times, seasons = WDSeasons, timeKnots = SeasonTimeKnots, seasonalStructure = WDSeasonalStructure, lambdas = c(0.003,0,240)) TrendTempM <- list(name = "Trend temp Mel", data = TempM, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(1e7,0,0)) TrendTempM2 <- list(name = "Trend temp Mel^2", data = TempM2, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(1e7,0,0)) Predictors <- list(Trend, WSeason, WDSeason, TrendTempM, TrendTempM2)
f = "elec.fit.RDS" if(file.exists(f)) { elec.fit = readRDS(f) } else { elec.fit <- STR(data = Data, predictors = Predictors, confidence = 0.95, gapCV = 48*7) saveRDS(elec.fit, f, compress = "xz") }
elec.fit <- STR(data = Data, predictors = Predictors, confidence = 0.95, gapCV = 48*7)
plot(elec.fit, xTime = as.Date("2000-01-11")+((Times-1)/48-10), forecastPanels = NULL)
The example below shows a simple way of forecasting seasonal data using STR.
TrendSeasonalStructure <- list(segments = list(c(0,1)), sKnots = list(c(1,0))) DailySeasonalStructure <- list(segments = list(c(0,48)), sKnots = c(as.list(1:47), list(c(48,0)))) WeeklySeasonalStructure <- list(segments = list(c(0,336)), sKnots = c(as.list(seq(4,332,4)), list(c(336,0)))) WDSeasonalStructure <- list(segments = list(c(0,48), c(100,148)), sKnots = c(as.list(c(1:47,101:147)), list(c(0,48,100,148)))) TrendSeasons <- rep(1, nrow(electricity)) DailySeasons <- as.vector(electricity[,"DailySeasonality"]) WeeklySeasons <- as.vector(electricity[,"WeeklySeasonality"]) WDSeasons <- as.vector(electricity[,"WorkingDaySeasonality"]) Data <- as.vector(electricity[,"Consumption"]) Times <- as.vector(electricity[,"Time"]) TempM <- as.vector(electricity[,"Temperature"]) TempM2 <- TempM^2 TrendTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 116) SeasonTimeKnots <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 24) SeasonTimeKnots2 <- seq(from = head(Times, 1), to = tail(Times, 1), length.out = 12) TrendData <- rep(1, length(Times)) SeasonData <- rep(1, length(Times)) Trend <- list(name = "Trend", data = TrendData, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(1500,0,0)) WSeason <- list(name = "Weekly seas", data = SeasonData, times = Times, seasons = WeeklySeasons, timeKnots = SeasonTimeKnots2, seasonalStructure = WeeklySeasonalStructure, lambdas = c(0.8,0.6,100)) WDSeason <- list(name = "Daily seas", data = SeasonData, times = Times, seasons = WDSeasons, timeKnots = SeasonTimeKnots, seasonalStructure = WDSeasonalStructure, lambdas = c(0.003,0,240)) TrendTempM <- list(name = "Trend temp Mel", data = TempM, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(1e7,0,0)) TrendTempM2 <- list(name = "Trend temp Mel^2", data = TempM2, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(1e7,0,0)) Predictors <- list(Trend, WSeason, WDSeason, TrendTempM, TrendTempM2)
The values, which need to be forecast, are supplied to STR as NAs. In our case we are going to forecast the last week of the original data.
Data[(length(Data)-7*48):length(Data)] <- NA
The forecasting is performed at the same time when the model is fitted.
f = "elec.fit.forecasting.RDS" if(file.exists(f)) { elec.fit = readRDS(f) } else { elec.fit <- STR(data = Data, predictors = Predictors, confidence = 0.95, gapCV = 48*7) saveRDS(elec.fit, f, compress = "xz") }
elec.fit <- STR(data = Data, predictors = Predictors, confidence = 0.95, gapCV = 48*7)
The result of the decomposition and forecasting is depicted below.
plot(elec.fit, xTime = as.Date("2000-01-11")+((Times-1)/48-10), forecastPanels = 7)
To check meaningfulness of the forecast it is advisable to plot beta coefficients of the decomposition. If the coefficients are too "wiggly" the forecast can be suboptimal.
Beta coefficients of the trend look a bit too "wiggly". It explains an uptrend in the forecast.
plotBeta(elec.fit, predictorN = 1)
Beta coefficients of weekly and daily seasonalities look smooth at the end of the time series.
plotBeta(elec.fit, predictorN = 2) plotBeta(elec.fit, predictorN = 3)
Beta coefficients for temperature and squared temperature predictors look smooth.
plotBeta(elec.fit, predictorN = 4) plotBeta(elec.fit, predictorN = 5)
Probably, the model can be re-estimated with a higher lambda coefficient for the trend to provide a better forecast.
Trend <- list(name = "Trend", data = TrendData, times = Times, seasons = TrendSeasons, timeKnots = TrendTimeKnots, seasonalStructure = TrendSeasonalStructure, lambdas = c(150000,0,0)) Predictors <- list(Trend, WSeason, WDSeason, TrendTempM, TrendTempM2)
f = "elec.fit.2.forecasting.RDS" if(file.exists(f)) { elec.fit.2 = readRDS(f) } else { elec.fit.2 <- STR(data = Data, predictors = Predictors, confidence = 0.95, gapCV = 48*7) saveRDS(elec.fit.2, f, compress = "xz") }
elec.fit.2 <- STR(data = Data, predictors = Predictors, confidence = 0.95, gapCV = 48*7)
The result gets much lower cross validated mean squared error.
plot(elec.fit.2, xTime = as.Date("2000-01-11")+((Times-1)/48-10), forecastPanels = 7)
Beta coefficients of the trend, seasonal components and predictors look smooth.
for(i in 1:5) plotBeta(elec.fit.2, predictorN = i)
To achieve higher performance, it is recommended to use Intel MKL for matrix operations. The simplest way to switch to Intel MKL is to install R distribution from Microsoft.
Note, that registering a parallel backend while using Intel MKL can reduce performance significanly since Intel MKL already tries to utilise all available cores. In case if a user has options to use Intel MKL or register a parallel backend (while using an ordinary computer, for example 4 phisical cores, 8 virtual cores), it is recommended to use Intel MKL and do not register a parallel backend.
For testing and exploration purposes it is recommended to avoid calculation of the confidence/forecasting intervals. Such calculation currently involves inversion of a big matrix and can take long time.
To monitor progress of the parameters' estimation it is recommended to set parameter trace to TRUE in STR call.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.