## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(cache = FALSE, fig.show = "hold", fig.width = 7, fig.height = 3)
library(stR)
library(forecast)
library(seasonal)

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

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

## ----include = FALSE----------------------------------------------------------
f <- "co2.fit.tbats.RDS"
if (file.exists(f)) {
  co2.fit <- readRDS(f)
} else {
  co2.fit <- tbats(co2)
  saveRDS(co2.fit, f, compress = "xz")
}

## ----eval = FALSE-------------------------------------------------------------
#  library(forecast)
#  co2.fit <- tbats(co2)

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

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

## ----include = FALSE----------------------------------------------------------
f <- "co2.fit.stR.RDS"
if (file.exists(f)) {
  co2.fit <- readRDS(f)
} else {
  co2.fit <- AutoSTR(co2)
  saveRDS(co2.fit, f, compress = "xz")
}

## ----eval = FALSE-------------------------------------------------------------
#  co2.fit <- AutoSTR(co2)

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  taylor.fit <- AutoSTR(taylor.msts, gapCV = 48, confidence = 0.95)

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

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

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

## -----------------------------------------------------------------------------
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(logGr))

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

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

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

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

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

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

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

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  gr.fit <- STR(data, predictors, confidence = 0.95, gap = 1, reltol = 0.00001)

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  gr.fit <- STR(data,
#    predictors,
#    confidence = 0.95,
#    gap = 1,
#    reltol = 0.00001
#  )

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

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

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  fit.str <- STR(as.vector(logGrOutl), predictors, confidence = 0.95, gapCV = 1, reltol = 0.001)

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  fit.rstr <- STR(as.vector(logGrOutl), predictors, confidence = 0.95, gapCV = 1, reltol = 0.001, nMCIter = 200, robust = TRUE)

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  calls.fit <- STR(
#    data = data,
#    predictors = predictors,
#    confidence = 0.95,
#    reltol = 0.003,
#    nFold = 4,
#    gap = 169
#  )

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  elec.fit <- STR(
#    data = Data,
#    predictors = Predictors,
#    confidence = 0.95,
#    gapCV = 48 * 7
#  )

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

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  elec.fit <- STR(
#    data = Data,
#    predictors = Predictors,
#    confidence = 0.95,
#    gapCV = 48 * 7
#  )

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

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

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

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

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

## ----include = FALSE----------------------------------------------------------
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")
}

## ----eval = FALSE-------------------------------------------------------------
#  elec.fit.2 <- STR(
#    data = Data,
#    predictors = Predictors,
#    confidence = 0.95,
#    gapCV = 48 * 7
#  )

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

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