## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  cache = FALSE,
  fig.show = "hold",
  fig.width = 7,
  fig.height = 3
)
cache <- function(object, file) {
  f <- paste0("../inst/extdata/", file)
  if (file.exists(f)) {
    readRDS(f)
  } else {
    saveRDS(
      object, 
      paste0("../inst/extdata/", file),
      compress = "xz",
    )
    object
  }
}

## ----libraries, message = FALSE-----------------------------------------------
library(stR)
library(forecast)
library(seasonal)

## ----classical, fig.height = 4.5----------------------------------------------
m <- decompose(co2)
plot(m)

## ----stl, fig.height = 4.5----------------------------------------------------
plot(stl(log(co2), s.window = "periodic", t.window = 30))

## ----tbats, include = FALSE---------------------------------------------------
co2.fit <- tbats(co2) |> cache("co2.fit.tbats.RDS")

## ----tbats_plot, fig.height = 4.5---------------------------------------------
plot(co2.fit)

## ----seas---------------------------------------------------------------------
co2.fit <- seas(co2)
plot(co2.fit, trend = TRUE)

## ----autostr, include = FALSE-------------------------------------------------
co2.fit <- AutoSTR(co2) |> cache("co2.fit.stR.RDS")

## ----autostr_plot, fig.height = 4---------------------------------------------
plot(co2.fit)

## ----taylor-------------------------------------------------------------------
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")

## ----taylor_fit, include = FALSE----------------------------------------------
taylor.fit <- AutoSTR(taylor.msts, gapCV = 48, confidence = 0.95) |>
  cache("taylor.fit.stR.RDS")

## ----taylor_plot, fig.height = 4.5--------------------------------------------
plot(taylor.fit)

## ----grocery------------------------------------------------------------------
plot(grocery, ylab = "NSW Grocery Turnover, $ 10^6")

## ----grocery_plot-------------------------------------------------------------
logGr <- log(grocery)
plot(logGr, ylab = "NSW Grocery Turnover, log scale")

## ----trendSeasonalStructure---------------------------------------------------
trendSeasonalStructure <- list(
  segments = list(c(0, 1)),
  sKnots = list(c(1, 0))
)

## ----seasonalStructure--------------------------------------------------------
seasonalStructure <- list(
  segments = list(c(0, 12)),
  sKnots = list(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, c(12, 0))
)

## ----seasons------------------------------------------------------------------
seasons <- as.vector(cycle(logGr))

## ----trendSeasons-------------------------------------------------------------
trendSeasons <- rep(1, length(logGr))

## ----times--------------------------------------------------------------------
times <- as.vector(time(logGr))

## ----data---------------------------------------------------------------------
data <- as.vector(logGr)

## ----trendTimeKnots-----------------------------------------------------------
trendTimeKnots <- seq(
  from = head(times, 1),
  to = tail(times, 1),
  length.out = 175
)

## ----seasonTimeKnots----------------------------------------------------------
seasonTimeKnots <- seq(
  from = head(times, 1),
  to = tail(times, 1),
  length.out = 15
)

## ----trendData----------------------------------------------------------------
trendData <- rep(1, length(logGr))
seasonData <- rep(1, length(logGr))

## ----trend--------------------------------------------------------------------
trend <- list(
  name = "Trend",
  data = trendData,
  times = times,
  seasons = trendSeasons,
  timeKnots = trendTimeKnots,
  seasonalStructure = trendSeasonalStructure,
  lambdas = c(0.5, 0, 0)
)

## ----season-------------------------------------------------------------------
season <- list(
  name = "Yearly seasonality",
  data = seasonData,
  times = times,
  seasons = seasons,
  timeKnots = seasonTimeKnots,
  seasonalStructure = seasonalStructure,
  lambdas = c(10, 0, 0)
)

## ----predictors---------------------------------------------------------------
predictors <- list(trend, season)

## ----gr.fit, include = FALSE--------------------------------------------------
gr.fit <- STR(data, predictors, confidence = 0.95, gap = 1, reltol = 0.00001) |>
  cache("gr.fit.stR.RDS")

## ----gr.fit_plot, fig.height=4------------------------------------------------
plot(gr.fit, xTime = times, forecastPanels = NULL)

## ----season2, echo=TRUE, warning=FALSE, results='hide'------------------------
season <- list(
  name = "Yearly seasonality",
  data = seasonData,
  times = times,
  seasons = seasons,
  timeKnots = seasonTimeKnots,
  seasonalStructure = seasonalStructure,
  lambdas = c(1, 1, 1)
)
predictors <- list(trend, season)

## ----gr.fit2, include = FALSE-------------------------------------------------
gr.fit <- STR(data,
  predictors,
  confidence = 0.95,
  gap = 1,
  reltol = 0.00001
) |>
  cache("gr.fit.2.stR.RDS")

## ----gr.fit2_plot, fig.height=4-----------------------------------------------
plot(gr.fit, xTime = times, forecastPanels = NULL)

## ----spikes, fig.height = 4---------------------------------------------------
outl <- rep(0, length(grocery))
outl[14] <- 900
outl[113] <- -700
tsOutl <- ts(outl, start = c(2000, 1), frequency = 12)

## ----logGROutl----------------------------------------------------------------
logGrOutl <- log(grocery + tsOutl)
plot(logGrOutl, ylab = "Log turnover with outliers")

## ----strspikes, echo=TRUE, warning=FALSE, results='hide'----------------------
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)

## ----fit.str, include = FALSE-------------------------------------------------
fit.str <- STR(
  as.vector(logGrOutl),
  predictors,
  confidence = 0.95,
  gapCV = 1,
  reltol = 0.001
) |>
  cache("logGrOutl.stR.RDS")

## ----fit.str_plot, fig.height = 4---------------------------------------------
plot(fit.str, xTime = times, forecastPanels = NULL)

## ----fit.rstr, include = FALSE------------------------------------------------
fit.rstr <- STR(
  as.vector(logGrOutl),
  predictors,
  confidence = 0.95, gapCV = 1, reltol = 0.001, nMCIter = 200, robust = TRUE
) |>
  cache("logGrOutl.stR.robust.RDS")

## ----fit.rstr_plot, fig.height = 4--------------------------------------------
plot(fit.rstr, xTime = times, forecastPanels = NULL)

## ----calls, echo=TRUE, warning=FALSE, results='hide'--------------------------
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)

## ----calls.fit, include = FALSE-----------------------------------------------
calls.fit <- STR(
  data = data,
  predictors = predictors,
  confidence = 0.95,
  reltol = 0.003,
  nFold = 4,
  gap = 169
) |>
  cache("calls.fit.RDS")

## ----calls.fit_plot, fig.height = 4-------------------------------------------
plot(calls.fit,
  xTime = as.Date("2003-03-03") +
    ((seq_along(data) - 1) / 169) +
    (((seq_along(data) - 1) / 169) / 5) * 2,
  forecastPanels = NULL
)

## ----electricity, echo=TRUE, warning=FALSE, results='hide'--------------------
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)

## ----elec.fit, include = FALSE------------------------------------------------
elec.fit <- STR(
  data = Data,
  predictors = Predictors,
  confidence = 0.95,
  gapCV = 48 * 7
) |>
  cache("elec.fit.RDS")

## ----elec.fit_plot, fig.height = 6--------------------------------------------
plot(elec.fit,
  xTime = as.Date("2000-01-11") + ((Times - 1) / 48 - 10),
  forecastPanels = NULL
)

## ----forecasting, echo=TRUE, warning=FALSE, results='hide'--------------------
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)

## ----nas, echo=TRUE, warning=FALSE, results='hide'----------------------------
Data[(length(Data) - 7 * 48):length(Data)] <- NA

## ----elec.fit2, include = FALSE-----------------------------------------------
elec.fit <- STR(
  data = Data,
  predictors = Predictors,
  confidence = 0.95,
  gapCV = 48 * 7
) |>
  cache("elec.fit.forecasting.RDS")

## ----elec.fit2_plot, fig.height = 7.5-----------------------------------------
plot(elec.fit,
  xTime = as.Date("2000-01-11") + ((Times - 1) / 48 - 10),
  forecastPanels = 7
)

## ----plotBeta, fig.height = 4.5-----------------------------------------------
plotBeta(elec.fit, predictorN = 1)

## ----plotBeta2, fig.height = 4.5----------------------------------------------
plotBeta(elec.fit, predictorN = 2)
plotBeta(elec.fit, predictorN = 3)

## ----plotBeta3, fig.height = 4.5----------------------------------------------
plotBeta(elec.fit, predictorN = 4)
plotBeta(elec.fit, predictorN = 5)

## ----higher_lambda, echo=TRUE, warning=FALSE, results='hide'------------------
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)

## ----elec.fit3, include = FALSE-----------------------------------------------
elec.fit.2 <- STR(
  data = Data,
  predictors = Predictors,
  confidence = 0.95,
  gapCV = 48 * 7
) |>
  cache("elec.fit.2.forecasting.RDS")

## ----elec.fig3_plot, fig.height = 7.5-----------------------------------------
plot(elec.fit.2,
  xTime = as.Date("2000-01-11") + ((Times - 1) / 48 - 10),
  forecastPanels = 7
)

## ----betas, fig.height = 4.5--------------------------------------------------
for (i in 1:5) {
  plotBeta(elec.fit.2, predictorN = i)
}