## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(  
  collapse = TRUE,  
  comment = "#>"  
)  

## ----setup--------------------------------------------------------------------
library(stlpro)  

## ----additive-data, fig.width=8, fig.height=4---------------------------------
set.seed(2106)  

n <- 15 * 12  
tt <- seq_len(n)  
  
trend <- 100 + 0.20 * tt  
season <- 12 * sin(2 * pi * tt / 12) +  
  4 * cos(2 * pi * tt / 12)  
noise <- rnorm(n, sd = 3)  
  
y <- trend + season + noise  
y_ts <- ts(y, start = c(2010, 1), frequency = 12)  
  
plot(y_ts, main = "Simulated monthly time series", ylab = "y")  

## ----additive-fit-------------------------------------------------------------
fit_additive <- stl(  
  y_ts,  
  period = 12,  
  s_window = 13
)  
  
lapply(fit_additive, head)

## ----additive-components, fig.width=8, fig.height=8---------------------------
plot(fit_additive)  

## ----additive-manual-plot, eval = FALSE, fig.width=8, fig.height=8------------
# par(mfrow = c(4, 1), mar = c(3, 4, 2, 1))
# 
# plot(y_ts, main = "Observed series", ylab = "")
# plot(fit_additive$trend, main = "Trend", ylab = "")
# plot(fit_additive$seasonal, main = "Seasonal component", ylab = "")
# plot(fit_additive$remainder, main = "Remainder", ylab = "")
# 
# par(mfrow = c(1, 1))

## ----multiplicative-data, fig.width=8, fig.height=4---------------------------
n <- 20 * 12  
tt <- seq_len(n)  
  
trend_mult <- exp(4.5 + 0.005 * tt)  
season_mult <- exp(  
  0.20 * sin(2 * pi * tt / 12) +  
    0.08 * cos(2 * pi * tt / 12)  
)  
remainder_mult <- exp(rnorm(n, sd = 0.06))  
  
y_mult <- trend_mult * season_mult * remainder_mult  
y_mult_ts <- ts(y_mult, start = c(2000, 1), frequency = 12)  
  
plot(  
  y_mult_ts,  
  main = "Simulated multiplicative monthly series",  
  ylab = "y"  
)  

## ----multiplicative-fit-------------------------------------------------------
fit_mult <- stl(  
  y_mult_ts,  
  period = 12,  
  s_window = 13,  
  multiplicative = TRUE
)  

## ----multiplicative-plot, fig.width=8, fig.height=8---------------------------
plot(fit_mult)  

## ----weekly-data, fig.width=8, fig.height=4-----------------------------------
n_years <- 25  
period_weekly <- 52.18  
n <- round(n_years * period_weekly)  
  
tt <- seq_len(n)  
  
trend_weekly <- 200 + 0.12 * tt  
season_weekly <- 25 * sin(2 * pi * tt / period_weekly) +  
  8 * cos(2 * pi * tt / period_weekly)  
noise_weekly <- rnorm(n, sd = 3)  
  
y_weekly <- trend_weekly + season_weekly + noise_weekly  
  
weekly_start <- as.Date("2015-01-05")  
weekly_dates <- weekly_start + 7 * (tt - 1)  
  
plot(  
  weekly_dates,  
  y_weekly,  
  type = "l",  
  xlab = "",  
  ylab = "y",  
  main = "Simulated weekly series with a 52.18-week seasonal cycle"  
)  

## ----weekly-fit---------------------------------------------------------------
fit_weekly <- stl(  
  y_weekly,  
  period = 52.18,  
  s_window = 7 
)  

## ----weekly-components, eval = TRUE, fig.width=8, fig.height=8----------------
plot(fit_weekly)

## ----weekly-compare-----------------------------------------------------------
fit_weekly_52 <- stl(  
  y_weekly,  
  period = 52,  
  s_window = 7  
)  
  
fit_weekly_5218 <- stl(  
  y_weekly,  
  period = 52.18,  
  s_window = 7  
)  

## ----weekly-compare-plot, eval = TRUE, fig.width=8, fig.height=8--------------
plot(  
  weekly_dates,  
  fit_weekly_52$seasonal,  
  type = "l",  
  col = "orange",  
  lwd = 1.5,  
  xlab = "",  
  ylab = "Seasonal component",  
  main = "Seasonal estimates using 52 and 52.18 week periods",
  ylim = c(-30, 35)
)  
  
lines(  
  weekly_dates,  
  fit_weekly_5218$seasonal,  
  col = "darkblue",  
  lwd = 1.5  
)  
  
legend(  
  "topright",  
  legend = c("period = 52", "period = 52.18"),  
  col = c("orange", "darkblue"),  
  lwd = 1.5,  
  bty = "n"  
)  

## ----cv-grid------------------------------------------------------------------
candidate_ns <- seq(from = 7, to = 51, by = 2) 

## ----cv-fit-direct, eval = TRUE-----------------------------------------------
fit_weekly_cv <- stl(  
  y_weekly,  
  period = 52.18,  
  s_window_candidates = candidate_ns
)  
  
fit_weekly_cv$cv_results$cv_choice

knitr::kable(fit_weekly_cv$cv_results$error_metrics)

## ----cv-plot, eval = TRUE, fig.width=8, fig.height=4--------------------------
plot(  
  candidate_ns,  
  fit_weekly_cv$cv_results$error_metrics,  
  type = "b",  
  pch = 19,  
  xlab = "Seasonal window (ns)",  
  ylab = "Cross-validation error",  
  main = "Seasonal-window cross-validation"  
)  

## ----outlier-data, fig.width=8, fig.height=4----------------------------------
y_outlier <- y_weekly  
y_outlier[c(80, 180, 260)] <- y_outlier[c(80, 180, 260)] + c(80, -100, 120)  
  
plot(  
  weekly_dates,  
  y_outlier,  
  type = "l",  
  xlab = "",  
  ylab = "y",  
  main = "Weekly series with artificial outliers"  
)  

## ----robust-fit---------------------------------------------------------------
fit_nonrobust <- stl(  
  y_outlier,  
  period = 52.18,  
  s_window = 7  
)  
  
fit_robust <- stl(  
  y_outlier,  
  period = 52.18,  
  s_window = 7,
  robust = TRUE
)  

## ----robust-plot, eval = TRUE, fig.width=8, fig.height=4----------------------
plot(  
  weekly_dates,  
  fit_nonrobust$seasonal,  
  type = "l",  
  col = "grey40",  
  xlab = "",  
  ylab = "Seasonal component",  
  main = "Effect of robust fitting"  
)  
  
lines(  
  weekly_dates,  
  fit_robust$seasonal,  
  col = "steelblue",  
  lwd = 1.5  
)  
  
legend(  
  "topright",  
  legend = c("Non-robust", "Robust"),  
  col = c("grey40", "steelblue"),  
  lwd = c(1, 1.5),  
  bty = "n"  
)  

## ----comparison-data----------------------------------------------------------
n <- 12 * 13  
tt <- seq_len(n)  
  
y_benchmark <- 100 +  
  0.15 * tt +  
  12 * sin(2 * pi * tt / 12) +  
  rnorm(n, sd = 4)  
  
y_benchmark_ts <- ts(y_benchmark, frequency = 12)  

## ----comparison_results, eval = TRUE------------------------------------------

stlpro_result <- stl(
  y_benchmark_ts,
  period = 12,
  s_window = 13,
  s_jump = 1,
  l_jump = 1,
  t_jump = 1,
  robust = FALSE
)


stats_stl_result <- stats::stl(
  y_benchmark_ts,
  s.window = 13,
  s.jump = 1,
  t.jump = 1,
  l.jump = 1,
  robust = FALSE
)  

res <- data.frame(stlpro = stlpro_result$seasonal, 
           stats_stl = stats_stl_result$time.series[,1],
           diff = stlpro_result$seasonal - stats_stl_result$time.series[,1])

res

max(res[,3])


 

## ----benchmark-system-time, eval = TRUE---------------------------------------
n_rep <- 1000 

time_stlpro <- system.time({
  for (i in seq_len(n_rep)) {
    stl(
      y_benchmark_ts,
      period = 12,
      s_window = 13,
      s_jump = 1,
      l_jump = 1,
      t_jump = 1,
      robust = FALSE
    )
  }
})  
  
time_stats <- system.time({  
  for (i in seq_len(n_rep)) {  
    stats::stl(  
      y_benchmark_ts,  
      s.window = 13,
      s.jump = 1,
      t.jump = 1, 
      l.jump = 1,
      robust = FALSE  
    )  
  }  
})  
  
rbind(  
  stlpro = time_stlpro,  
  stats_stl = time_stats  
)  

## ----benchmark-daily, eval = TRUE---------------------------------------------
n <- 12 * 366 
tt <- seq_len(n)  
  
y_benchmark <- 100 +  
  0.15 * tt +  
  20 * sin(2 * pi * tt / 365.2425) +  
  rnorm(n, sd = 4)  
  
y_daily <- ts(y_benchmark, frequency = 365.2425) 

n_rep <- 1 

time_stlpro_daily <- system.time({
  for (i in seq_len(n_rep)) {
    stl(
      y_daily,
      period = 12
    )
  }
})

time_stlpro_daily