# Generated by `rjournal_pdf_article()` using `knitr::purl()`: do not edit by hand
# Please edit schulz.Rmd to modify this file

## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# ### R Script for Article "deseats: An R Package for Data-Driven Trend and
# ###                       Seasonality Estimation in Time Series"


## ----setup2, include=FALSE, purl = TRUE---------------------------------------

##### Setup and definition of helper functions for table formatting etc.

library(kableExtra)
library(magrittr)
library(purrr)

# Make landscape orientation for table in PDF
landscape_tab <- function(tab) {
  if(knitr::is_latex_output()) {
    return(kableExtra::landscape(tab))
  } else {
    return(tab)
  }
}

# Reduce vertical space in PDF
reduce_vspace <- function(space = "-0.9cm") {
  if(knitr::is_latex_output()) {
    return(knitr::asis_output(paste0("\\vspace{", space, "}")))
  }
  return(invisible(""))
}

# Change from bold in PDF to bold in HTML for table entries
adjust_bold_html <- function(df, cols = 3:20) {

  df[, cols] <- df[, cols] %>%
    purrr::map_dfc(
      function(.x) {
        w <- which(grepl("\\bftab", .x, fixed = TRUE))
        .x_sub <- gsub("\\bftab ", replacement = "",
                       x = .x[w], fixed = TRUE)
        .x[w] <- paste0(
          "<strong>", 
          .x_sub, 
          "</strong>"
        )
        .x
      }
    )
  df
}

adjust_html_tab_size <- function(x) {
  `%>%` <- magrittr::`%>%`
  if (knitr::is_html_output()) {
    x <- x %>% 
      kableExtra::kable_styling(full_width = FALSE) %>%
      kableExtra::scroll_box(width = "100%")
  }
  x
}










## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### DeSeaTS decomposition with deseats


## ----echo = TRUE, eval = TRUE, include = TRUE, purl = TRUE--------------------
library(deseats)      # Attach the deseats-package
LGDP <- log(GDP)      # Get data
NLF <- NOLABORFORCE
LCO <- log(COVID)
sm_options3 <- set_options()                # Local cubic trend
sm_options1 <- set_options(order_poly = 1)  # Local linear trend
LGDP_deseats <- deseats(LGDP, sm_options1)  # Run DeSeaTS
NLF_deseats <- deseats(NLF, sm_options3)
LCO_deseats <- deseats(LCO, sm_options3)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Plotting of results




## ----plotdec, fig.cap = "The estimated trend and seasonal components for the LGDP, NLF and LCO data following DeSeaTS. The estimated seasonal components are shifted for better visibility.", fig.height = 4.5, fig.width = 5.6, echo = FALSE, fig.pos = 'h', purl = TRUE----
par_old <- par(no.readonly = TRUE)
par(fig = c(0, 0.5, 0.5, 1), cex = 0.5)
plot(LGDP_deseats, which = 5, s_around = 3.9, main = "(a) LGDP", xlab = "Year")
par(fig = c(0.5, 1, 0.5, 1), new = TRUE)
plot(NLF_deseats, which = 5, s_around = 55, main = "(b) NLF", xlab = "Year")
par(fig = c(0.25, 0.75, 0, 0.5), new = TRUE)
plot(LCO_deseats, which = 5, s_around = -2.2, main = "(c) LCO", 
     xlab = "Calendar week of year 2021")
par(par_old)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Benchmark method


## ----echo = FALSE, eval = TRUE, include=TRUE, warning = FALSE, message = FALSE----
library(forecast)
stl_GDP <- mstl(log(GDP))
stl_NLF <- mstl(NLF)
stl_COVID <- mstl(log(COVID)) 
adj_stl_GDP <- as.numeric(exp(stl_GDP[, 1] - stl_GDP[, 3]))
adj_stl_NLF <- as.numeric(stl_NLF[, 1] - stl_NLF[, 3]) 
adj_stl_COVID <- as.numeric(exp(stl_COVID[, 1] - stl_COVID[, 3]))
t1 <- as.numeric(time(stl_GDP))
t2 <- as.numeric(time(stl_NLF))
t3 <- as.numeric(time(stl_COVID))

fc_GDP_stl <- forecast(stl_GDP, method = "arima", h = 4, lambda = 0, biasadj = TRUE)$mean
fc_NLF_stl <- forecast(stl_NLF, method = "arima", h = 12, lambda = NULL)$mean
fc_COVID_stl <- forecast(stl_COVID, method = "arima", h = 7, lambda = 0, biasadj = TRUE)$mean
t1_fc <- as.numeric(time(fc_GDP_stl))
t2_fc <- as.numeric(time(fc_NLF_stl))
t3_fc <- as.numeric(time(fc_COVID_stl))
fc_GDP_stl <- as.numeric(fc_GDP_stl)
fc_NLF_stl <- as.numeric(fc_NLF_stl)
fc_COVID_stl <- as.numeric(fc_COVID_stl)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Plots of seasonally adjusted series




## ----plotdeseas, fig.cap = "The original example series (grey and solid) together with the seasonally adjusted series following DeSeaTS (blue and solid). The seasonally adjusted series following automated STL decompositions (red and dashed) are displayed as a benchmark.", fig.height = 4.5, fig.width = 5.6, echo = FALSE, fig.pos = "h", purl = TRUE----
GDP_sa <- exp(deseasonalize(LGDP_deseats))
NLF_sa <- deseasonalize(NLF_deseats)
COVID_sa <- exp(deseasonalize(LCO_deseats))
par_old <- par(no.readonly = TRUE)
par(fig = c(0, 0.5, 0.5, 1), cex = 0.5)
ts.plot(GDP, GDP_sa, col = c("grey60", "blue"), main = "(a) Real GDP of the US", 
     xlab = "Year", ylab = "Billions of USD")
lines(t1, adj_stl_GDP, col = "red", lty = 2)
par(fig = c(0.5, 1, 0.5, 1), new = TRUE)
ts.plot(NLF, NLF_sa, col = c("grey60", "blue"), ylab = "Millions of persons",
        main = "(b) US persons not in the labor force", xlab = "Year")
lines(t2, adj_stl_NLF, col = "red", lty = 2)
par(fig = c(0.25, 0.75, 0, 0.5), new = TRUE)
ts.plot(COVID, COVID_sa, col = c("grey60", "blue"), 
        xlab = "Calendar week of year 2021", ylab = "Thousands of cases",
        main = "(c) Newly confirmed COVID-19 Cases in Germany")
lines(t3, adj_stl_COVID, col = "red", lty = 2)
par(par_old)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Fit S-Semi-ARMA models


## ----echo = TRUE, eval = TRUE, include = TRUE---------------------------------
nar_lim <- nma_lim <- c(0, 2)
model_LGDP <- s_semiarma(LGDP, sm_options1, nar_lim = nar_lim, nma_lim = nma_lim) 
model_NLF <- s_semiarma(NLF, sm_options3, nar_lim = nar_lim, nma_lim = nma_lim)
model_LCO <- s_semiarma(LCO, sm_options3, nar_lim = nar_lim, nma_lim = nma_lim) 


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Run seasonality tests on residuals


## ----tabtests, echo = FALSE, eval = TRUE, include = TRUE----------------------
library(knitr)
library(kableExtra)

# For seasonality tests
source("scripts/seasonality_tests.R") 

sprintf4 <- function(x) {sprintf("%.4f", x)}

# Apply seasonality tests (to residual data of full models)
res_LGDP <- model_LGDP@par_model$residuals
res_NLF <- model_NLF@par_model$residuals
res_LCO <- model_LCO@par_model$residuals

qs_LGDP <- qs_test_internal(res_LGDP, 2) 
qs_NLF <- qs_test_internal(res_NLF, 2)
qs_LCO <- qs_test_internal(res_LCO, 2)

qs_p_LGDP <- sprintf4(qs_LGDP$p_val)
qs_p_NLF <- sprintf4(qs_NLF$p_val)
qs_p_LCO <- sprintf4(qs_LCO$p_val)

fried_LGDP <- fried_test(res_LGDP)
fried_NLF <- fried_test(res_NLF)
fried_LCO <- fried_test(res_LCO)

fried_p_LGDP <- sprintf4(fried_LGDP$p_val)
fried_p_NLF <- sprintf4(fried_NLF$p_val)
fried_p_LCO <- sprintf4(fried_LCO$p_val)

bwidth_LGDP <- sprintf4(bwidth(LGDP_deseats))
bwidth_NLF <- sprintf4(bwidth(NLF_deseats))
bwidth_LCO <- sprintf4(bwidth(LCO_deseats))

ar1_LGDP <- sprintf4(model_LGDP@par_model$coef[["ar1"]])
ar2_LGDP <- sprintf4(model_LGDP@par_model$coef[["ar2"]])
ma1_LGDP <- sprintf4(model_LGDP@par_model$coef[["ma1"]])
se_ar1_LGDP <- sprintf4(sqrt(model_LGDP@par_model$var.coef[1, 1]))
se_ar2_LGDP <- sprintf4(sqrt(model_LGDP@par_model$var.coef[2, 2]))
se_ma1_LGDP <- sprintf4(sqrt(model_LGDP@par_model$var.coef[3, 3]))

ar1_NLF <- sprintf4(model_NLF@par_model$coef[["ar1"]])
se_ar1_NLF <- sprintf4(sqrt(model_NLF@par_model$var.coef[[1, 1]]))

ma1_LCO <- sprintf4(model_LCO@par_model$coef[["ma1"]])
se_ma1_LCO <- sprintf4(sqrt(model_LCO@par_model$var.coef[[1, 1]]))

df_tab <- data.frame(
  Series = c("LGDP", "NLF", "LCO"),
  Bandwidth = c(bwidth_LGDP, bwidth_NLF, bwidth_LCO),
  ar1 = c(ar1_LGDP, ar1_NLF, "---"),
  se_ar1 = c(se_ar1_LGDP, se_ar1_NLF, "---"),
  ar2 = c(ar2_LGDP, "---", "---"),
  se_ar2 = c(se_ar2_LGDP, "---", "---"),
  ma1 = c(ma1_LGDP, "---", ma1_LCO),
  se_ma1 = c(se_ma1_LGDP, "---", se_ma1_LCO),
  qs_p = c(qs_p_LGDP, qs_p_NLF, qs_p_LCO),
  fried_p = c(fried_p_LGDP, fried_p_NLF, fried_p_LCO)
)

fmt <- if(knitr::is_latex_output()) {
  "latex"
} else if (knitr::is_html_output()) {
  "html"
}

kable(df_tab, booktabs = TRUE, col.names = c("Series", "Bandwidth", "$\\hat\\phi_1$", "$se\\left(\\hat\\phi_1\\right)$", "$\\hat\\phi_2$", "$se\\left(\\hat\\phi_2\\right)$", "$\\hat\\theta_1$", "$se\\left(\\hat\\theta_1\\right)$", "$p$-value (QS)", "$p$-value (Friedman)"), escape = FALSE, format = fmt, caption = "Selected bandwidths for the three series LGDP, NLF and LCO, estimated ARMA parameters and corresponding standard errors for DeSeaTS residuals and $p$-values for QS tests and Friedman tests applied to S-Semi-ARMA residuals.") %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Forecasting


## ----eval = TRUE, echo = TRUE-------------------------------------------------
set.seed(1) 
fc_GDP_DS <- predict(model_LGDP, n.ahead = 4, method = "boot", expo = TRUE)
fc_NLF_DS <- predict(model_NLF, n.ahead = 12, method = "boot")
fc_COVID_DS <- predict(model_LCO, n.ahead = 7, method = "boot", expo = TRUE)




## ----plotfcast, fig.cap = "S-Semi-ARMA point and 95% / 99% interval forecasts (blue) together with automated ARIMA-based STL point forecasts (red) as a benchmark.", fig.height = 4.5, fig.width = 5.6, eval = TRUE, fig.pos = 't', echo = FALSE, purl = TRUE----
par_old <- par(no.readonly = TRUE)
par(fig = c(0, 0.5, 0.5, 1), cex = 0.5)
plot(fc_GDP_DS, main = "(a) Real GDP of the US", xlab = "Year", 
     ylab = "Billions of USD", xlim = c(2010, 2020.75), ylim = c(3500, 6000))
lines(t1_fc, fc_GDP_stl, col = "red", lty = 2)
par(fig = c(0.5, 1, 0.5, 1), new = TRUE)
plot(fc_NLF_DS, main = "(b) US persons not in the labor force", xlab = "Year", 
     ylab = "Millions of persons", xlim = c(2010, 2020.917), ylim = c(80, 97.2))
lines(t2_fc, fc_NLF_stl, col = "red", lty = 2)
par(fig = c(0.25, 0.75, 0, 0.5), new = TRUE)
plot(fc_COVID_DS, main = "(c) Newly confirmed COVID-19 Cases in Germany", 
     xlab = "Calendar week of year 2021", ylab = "Thousands of cases",
     xlim = c(45, 49.143), ylim = c(10, 85))
lines(t3_fc, fc_COVID_stl, col = "red", lty = 2)
par(par_old)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Examples for first simulation


## ----plotsimexamples, fig.cap = "Examples of quarterly simulated series (grey and solid) together with the true underlying trends (blue and dashed) for $n_1$ and the estimated trends (red and solid) using the DeSeaTS framework with $r=3$. All examples are with number of observations $n_1$ and the same simulated errors from $\\{e_{1,t}\\}$.", echo = FALSE, fig.height = 5, fig.pos = "t", fig.align='center'----
library(deseats)
library(magrittr)
library(tidyverse)

set.seed(123) 

### Number of observations

ny <- 15
n1 <- ny * 4

# Time points

t1 <- 1:n1
xt_n1 <- (t1) / n1

### Trend functions

# Trend 1

t1_n1 <- (5 * xt_n1 + 1) %>%
  ts(frequency = 4)

# Trend 2

t2_n1 <- (2.5 * tanh(5 * (xt_n1 - 0.4)) + 3.5) %>%
  ts(frequency = 4)

# Trend 3

t3_n1 <- (3.2 * xt_n1 + 0.92 * (sin(3.2 * pi * (xt_n1 - 0.2))) + 1.85) %>%
  ts(frequency = 4)

### Seasonalities

## Exactly periodic and with frequency 4

quart <- (t1 - 1) %% 4
s1 <- rep(-0.39, length(t1))
s1[quart == 1] <- 0.01
s1[quart == 2] <- -0.04
s1[quart == 3] <- 0.42

s1_n1 <- s1 %>%
  ts(frequency = 4)

#### Slowly changing and with frequency 4
f_sc <- 0.83 * (1  +0.1 * sin(50 * pi * xt_n1 / 9) + 0.2 * cos(3 * pi * xt_n1))
s2_n1 <-  f_sc * s1_n1

### Errors

arma1 <- list(ar = 0.75, sd_e = 0.083)

# Example simulations
e1 <- arima.sim(model = list(ar = arma1$ar), sd = arma1$sd_e, n = n1) %>%
  ts(frequency = 4)

series1s1 <- t1_n1 + s1_n1 + e1
series1s2 <- t1_n1 + s2_n1 + e1
series2s1 <- t2_n1 + s1_n1 + e1
series2s2 <- t2_n1 + s2_n1 + e1
series3s1 <- t3_n1 + s1_n1 + e1
series3s2 <- t3_n1 + s2_n1 + e1

est1s1 <- deseats(series1s1, set_options(order_poly = 3))
est1s2 <- deseats(series1s2, set_options(order_poly = 3))
est2s1 <- deseats(series2s1, set_options(order_poly = 3))
est2s2 <- deseats(series2s2, set_options(order_poly = 3))
est3s1 <- deseats(series3s1, set_options(order_poly = 3))
est3s2 <- deseats(series3s2, set_options(order_poly = 3))

par_old <- par(no.readonly = TRUE)
par(mfrow = c(3, 2), cex.axis = 0.9, cex.lab = 0.9, cex.main = 0.9,
    mar = c(2, 4, 3, 2) + 0.1, oma = c(2, 0, 0, 0))
ts.plot(series1s1, t1_n1, trend(est1s1),
     col = c("grey50", "blue", "red"),
     main = expression(paste(bold("(a) "), italic(g[1]), ", ", italic(s[1]))),
     lty = c(1, 2, 1))
ts.plot(series1s2, t1_n1, trend(est1s2),
     col = c("grey50", "blue", "red"),
     main = expression(paste(bold("(b) "), italic(g[1]), ", ", italic(s[2]))),
     lty = c(1, 2, 1))
ts.plot(series2s1, t2_n1, trend(est2s1),
     col = c("grey50", "blue", "red"),
     main = expression(paste(bold("(c) "), italic(g[2]), ", ", italic(s[1]))),
     lty = c(1, 2, 1))
ts.plot(series2s2, t2_n1, trend(est2s2),
     col = c("grey50", "blue", "red"),
     main = expression(paste(bold("(d) "), italic(g[2]), ", ", italic(s[2]))),
     lty = c(1, 2, 1))
par(xpd = NA)
ts.plot(series3s1, t3_n1, trend(est3s1),
     col = c("grey50", "blue", "red"),
     main = expression(paste(bold("(e) "), italic(g[3]), ", ", italic(s[1]))),
     lty = c(1, 2, 1))
par(xpd = NA)
ts.plot(series3s2, t3_n1, trend(est3s2),
     col = c("grey50", "blue", "red"),
     main = expression(paste(bold("(f) "), italic(g[3]), ", ", italic(s[2]))),
     lty = c(1, 2, 1))
par(par_old)


## ----plotsimexamples2, fig.cap = "True underlying seasonal components (grey and solid) of the quarterly simulated examples and the estimated seasonal components (red and solid) using the DeSeaTS framework with $r=3$. On the time axis, the true values are slightly shifted to the left and the estimates slightly to the right for better comparability. All examples are with number of observations $n_1$ and the same simulated errors from $\\{e_{1,t}\\}$.", echo = FALSE, fig.height = 5, fig.pos = "t", fig.align='center'----
par_old <- par(no.readonly = TRUE)
par(mfrow = c(3, 2), cex.axis = 0.9, cex.lab = 0.9, cex.main = 0.9,
    mar = c(2, 4, 3, 2) + 0.1, oma = c(2, 0, 0, 0))
tp <- c(time(s1_n1)) 
matplot(cbind(tp - 0.07, tp + 0.07), cbind(c(s1_n1), c(season(est1s1))),
    type = "hh",
     col = c("grey50", "red"),
     main = expression(paste(bold("(a) "), italic(g[1]), ", ", italic(s[1]))),
    xlab = "Time", ylab = "", lty = c(1, 1))
matplot(cbind(tp - 0.07, tp + 0.07), cbind(c(s2_n1), c(season(est1s2))),
    type = "hh",
     col = c("grey50", "red"),
     main = expression(paste(bold("(b) "), italic(g[1]), ", ", italic(s[2]))),
    xlab = "Time", ylab = "", lty = c(1, 1))
matplot(cbind(tp - 0.07, tp + 0.07), cbind(c(s1_n1), c(season(est2s1))),
    type = "hh",
     col = c("grey50", "red"),
     main = expression(paste(bold("(c) "), italic(g[2]), ", ", italic(s[1]))),
    xlab = "Time", ylab = "", lty = c(1, 1))
matplot(cbind(tp - 0.07, tp + 0.07), cbind(c(s2_n1), c(season(est2s2))),
    type = "hh",
     col = c("grey50", "red"),
     main = expression(paste(bold("(d) "), italic(g[2]), ", ", italic(s[2]))),
    xlab = "Time", ylab = "", lty = c(1, 1))
par(xpd = NA)
matplot(cbind(tp - 0.07, tp + 0.07), cbind(c(s1_n1), c(season(est3s1))),
    type = "hh",
     col = c("grey50", "red"),
     main = expression(paste(bold("(e) "), italic(g[3]), ", ", italic(s[1]))),
    xlab = "Time", ylab = "", lty = c(1, 1))
par(xpd = NA)
matplot(cbind(tp - 0.07, tp + 0.07), cbind(c(s2_n1), c(season(est3s2))),
    type = "hh",
     col = c("grey50", "red"),
     main = expression(paste(bold("(f) "), italic(g[3]), ", ", italic(s[2]))),
    xlab = "Time", ylab = "", lty = c(1, 1))
par(par_old)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Examples for stochastic components for second simulation study


## ----stochexample, echo = FALSE, fig.cap = "A realization of the stochastic trend $G_{2,t}$ and of the stochastic seasonal component $S_{3,t}$ for 15 years.", fig.height = 2.5, fig.width = 5.6, fig.pos = 'h'----
t2_fun <- function(xt) {
  n <- length(xt)
  et <- rnorm(n = n + 200) + 0.05

  g_out <- tail(cumsum(et), n)
  (g_out - mean(g_out)) / sd(g_out)

}

s3_fun <- function(t, sig2) {
  sig <- sqrt(sig2)
  ps <- 12

  n <- length(t)
  m <- 300 * 12

  sx <- rep(0, n + m)
  sx[1:11] <- rnorm(n = ps - 1, mean = 0, sd = sig)

  for (i in ps:(n + m)) {
    sx[[i]] <- -sum(sx[(i - 11):(i - 1)]) + rnorm(n = 1, mean = 0, sd = sig)
  }
  sx <- tail(sx, n)

  (sx - mean(sx)) / sd(sx)

}

set.seed(1)

n <- 15 * 12
t <- 1:n
xt <- t / n

sig2 <- 0.025

example_trend2 <- ts(t2_fun(xt), frequency = 12)
example_season3 <- ts(s3_fun(t, sig2), frequency = 12)

main1 <- expression(paste(bold("(a) "), italic(G[2*","*t])))
main2 <- expression(paste(bold("(b) "), italic(S[3*","*t])))

par_old <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), cex = 0.6, mar = c(5, 3, 4, 1) + 0.1)
plot(example_trend2, main = main1, ylab = "")
plot(example_season3, main = main2, ylab = "")
par(par_old)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### X13 and TRAMO-SEATS applied to real examples


## ----figuresx13, eval = TRUE, echo = FALSE, include=TRUE, fig.cap = "Seasonally adjusted series for the US-GDP series and of the series about US persons not in the labor force following the automated X13 (blue and solid) and TRAMO-SEATS (red and dashed) subroutines of RJDemetra+. X13 used a 5-terms Henderson trend filter and a 3x5 seasonal filter for the series in (a), whereas a 13-terms Henderson trend filter and a 3x5 seasonal filter were considered for the series in (b).", fig.height = 2.6, warning = FALSE, fig.width = 5.6, fig.pos='h'----
library(RJDemetra)
library(deseats)
X13_GDP <- x13(GDP)
X13_NLF <- x13(NOLABORFORCE)
TS_GDP <- tramoseats(GDP)
TS_NLF <- tramoseats(NOLABORFORCE)
X13_sa_GDP <- X13_GDP$final$series[, "sa"]
X13_sa_NLF <- X13_NLF$final$series[, "sa"]
TS_sa_GDP <- TS_GDP$final$series[, "sa"]
TS_sa_NLF <- TS_NLF$final$series[, "sa"]

old_par <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), cex = 0.5)
ts.plot(GDP, X13_sa_GDP, TS_sa_GDP, col = c("grey60", "blue", "red"), lty = c(1, 1, 2),
        xlab = "Year", ylab = "Billions of Australian dollars", main = "(a) Real consumption expenditure for Australia")
ts.plot(NOLABORFORCE, X13_sa_NLF, TS_sa_NLF, col = c("grey60", "blue", "red"), lty = c(1, 1, 2),
        xlab = "Year", main = "(b) US persons not in the labor force",
        ylab = "Millions of persons")
par(old_par)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Tables for first simulation study


## ----sim1tab1, echo = FALSE---------------------------------------------------
library(knitr)
library(kableExtra)
data0 <- read.csv("sim-data/Sim1Table-Quarterly-Error1.csv") 
data0$X.r.[is.na(data0$X.r.)] <- ""
kable(data0, col.names = c("$r$", "Inflation", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)),
      booktabs = TRUE, escape = FALSE, format = fmt,
      linesep = c("\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5), "\\addlinespace", "\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5)),
    caption = "DeSeaTS simulation results under consideration of quarterly data and error process $\\{e_{1,t}\\}$. $\\text{SD}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat g\\right)$, $\\text{MSE}\\left(\\hat s\\right)$ and $\\text{MSE}\\left(\\hat m\\right)$ were multiplied by 1000.") %>%
  add_header_above(
    c(" ", " ", " ", "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3, "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3), escape = FALSE
  ) %>%
  add_header_above(
    c(" ", " ", " ", "$s_1$" = 9, "$s_2$" = 9), escape = FALSE
  ) %>%  
  kable_styling(font_size = 7.5) %>%
  adjust_html_tab_size() %>%
  landscape_tab() 


## ----sim1tab2, echo = FALSE---------------------------------------------------
library(knitr)
library(kableExtra)
data0 <- read.csv("sim-data/Sim1Table-Quarterly-Error2.csv") 
data0$X.r.[is.na(data0$X.r.)] <- ""
kable(data0, col.names = c("$r$", "Inflation", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)),
      booktabs = TRUE, escape = FALSE, format = fmt,
      linesep = c("\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5), "\\addlinespace", "\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5)),
    caption = "DeSeaTS simulation results under consideration of quarterly data and error process $\\{e_{2,t}\\}$. $\\text{SD}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat g\\right)$, $\\text{MSE}\\left(\\hat s\\right)$ and $\\text{MSE}\\left(\\hat m\\right)$ were multiplied by 1000.") %>%
  add_header_above(
    c(" ", " ", " ", "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3, "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3), escape = FALSE
  ) %>%
  add_header_above(
    c(" ", " ", " ", "$s_1$" = 9, "$s_2$" = 9), escape = FALSE
  ) %>%  
  kable_styling(font_size = 7.5) %>%
  adjust_html_tab_size() %>%
  landscape_tab() 


## ----sim1tab3, echo = FALSE---------------------------------------------------
data0 <- read.csv("sim-data/Sim1Table-Monthly-Error1.csv") 
data0$X.r.[is.na(data0$X.r.)] <- ""
kable(data0, col.names = c("$r$", "Inflation", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)),
      booktabs = TRUE, escape = FALSE, format = fmt,
      linesep = c("\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5), "\\addlinespace", "\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5)),
    caption = "DeSeaTS simulation results under consideration of monthly data and error process $\\{e_{1,t}\\}$. $\\text{SD}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat g\\right)$, $\\text{MSE}\\left(\\hat s\\right)$ and $\\text{MSE}\\left(\\hat m\\right)$ were multiplied by 1000.") %>%
  add_header_above(
    c(" ", " ", " ", "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3, "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3), escape = FALSE
  ) %>%
  add_header_above(
    c(" ", " ", " ", "$s_1$" = 9, "$s_2$" = 9), escape = FALSE
  ) %>%  
  kable_styling(font_size = 7.5) %>%
  adjust_html_tab_size() %>%
  landscape_tab() 


## ----sim1tab4, echo = FALSE---------------------------------------------------
data0 <- read.csv("sim-data/Sim1Table-Monthly-Error2.csv") 
data0$X.r.[is.na(data0$X.r.)] <- ""
kable(data0, col.names = c("$r$", "Inflation", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)),
      booktabs = TRUE, escape = FALSE, format = fmt,
      linesep = c("\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5), "\\addlinespace", "\\addlinespace", rep("", 5), "\\addlinespace", rep("", 5)),
    caption = "DeSeaTS simulation results under consideration of monthly data and error process $\\{e_{2,t}\\}$. $\\text{SD}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat h\\right)$, $\\text{MSE}\\left(\\hat g\\right)$, $\\text{MSE}\\left(\\hat s\\right)$ and $\\text{MSE}\\left(\\hat m\\right)$ were multiplied by 1000.") %>%
  add_header_above(
    c(" ", " ", " ", "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3, "$g_1$" = 3, "$g_2$" = 3, "$g_3$" = 3), escape = FALSE
  ) %>%
  add_header_above(
    c(" ", " ", " ", "$s_1$" = 9, "$s_2$" = 9), escape = FALSE
  ) %>%  
  kable_styling(font_size = 7.5) %>%
  adjust_html_tab_size() %>%
  landscape_tab() 


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Boxplot figures for first simulation study


## -----------------------------------------------------------------------------
unwrap <- function() {
  `%>%` <- magrittr::`%>%`
  data_bw <- read.csv("sim-data/sim1_all_bwidths.csv", check.names = FALSE) %>%
    tidyr::pivot_longer(
      cols = 1:288, 
      names_to = "Setting",
      values_to = "bwidth"
    ) %>%
    separate_wider_delim(
      cols = "Setting",
      delim = "_",
      names = c("N_Year", "Seasonal_Period", "Trend", "Season", "Error_Process",
                "r", "infr")
    ) %>%
    mutate(
      N_Year = as.numeric(N_Year),
      Seasonal_Period = as.numeric(Seasonal_Period),
      infr = factor(infr, levels = c("opt", "nai"))
    ) %>%
    arrange(Seasonal_Period, Error_Process, Season, Trend, r, N_Year, infr) %>%
    group_by(Seasonal_Period, Error_Process, r) %>%
    group_split()
  nams <- map_chr(data_bw, function(.x) {

    paste0("sp", .x$Seasonal_Period[[1]], "_r", .x$r[[1]], "_", .x$Error_Process[[1]])
  })
  names(data_bw) <- nams
  data_bw <- data_bw %>%
    map(function(.x) {
      .x %>%
        select(-Error_Process, -r)
    })
}

boxplot_fun <- function(df, ep, poly_ord) {
  parold <- par(no.readonly = TRUE)
  par(mfrow = c(6, 2), mar = c(4, 4, 3, 2) + 0.1)
  on.exit(expr = {par(parold)}, add = TRUE, after = TRUE)

  season_p <- df$Seasonal_Period[[1]]

  df %>%
    mutate(N_Year = factor(N_Year)) %>%
    mutate(infr = factor(infr, levels = c("opt", "nai"))) %>%
    arrange(Trend, Season, N_Year, infr) %>%
    group_by(Season, Trend, infr) %>%
    group_split() %>%
    walk(function(.x, ep, poly_ord, season_p) {
      tn <- switch(
        .x$Trend[[1]],
        "t1" = 1,
        "t2" = 2,
        "t3" = 3
      )
      sn <- switch(
        .x$Season[[1]],
        "s1" = 1,
        "s2" = 2
      )
      ifr <- switch(
        as.character(.x$infr[[1]]),
        "opt" = "optimal",
        "nai" = "naive"
      )

      arma_list <- if (ep == 1) {
        list(ar = 0.75, ma = NULL, sd_e = 0.083)
      } else if (ep == 2) {
        list(ar = c(1.05, -0.3), ma = NULL, sd_e = 0.06)
      }

      cb <- if (poly_ord == 1) {0.05} else if (poly_ord == 3) {0.1}

      if (tn == 1) {
        bw1 <- bw2 <- bw3 <- 0.49
      } else if (tn == 2) {
        fun <- expression(2.5 * tanh(5 * (x - 0.4)) + 3.5)
        bw1 <- hA_calc(fun, arma = arma_list, p = poly_ord, mu = 1, frequ = season_p, n = season_p * 15, cb = cb)$hA
        bw2 <- hA_calc(fun, arma = arma_list, p = poly_ord, mu = 1, frequ = season_p, n = season_p * 30, cb = cb)$hA
        bw3 <- hA_calc(fun, arma = arma_list, p = poly_ord, mu = 1, frequ = season_p, n = season_p * 60, cb = cb)$hA
      } else if (tn == 3) {
        fun <- expression(3.2 * x + 0.92 * (sin(3.2 * pi * (x - 0.2))) + 1.85)
        bw1 <- hA_calc(fun, arma = arma_list, p = poly_ord, mu = 1, frequ = season_p, n = season_p * 15, cb = cb)$hA
        bw2 <- hA_calc(fun, arma = arma_list, p = poly_ord, mu = 1, frequ = season_p, n = season_p * 30, cb = cb)$hA
        bw3 <- hA_calc(fun, arma = arma_list, p = poly_ord, mu = 1, frequ = season_p, n = season_p * 60, cb = cb)$hA
      }

      y_min <- min(bw1, bw2, bw3, .x$bwidth)
      y_max <- max(bw1, bw2, bw3, .x$bwidth)
      ylim <- c(y_min, y_max)

      boxplot(bwidth ~ N_Year, data = .x, xlab = "Number of years",
              ylab = "Estimated bandwidth", ylim = ylim, col = "white",
              main = substitute(italic(g[trn])*", "*italic(s[srn])*", "*infl*" inflation rate", list(trn = tn, srn = sn, infl = ifr)))
      points(c(1, 2, 3), c(bw1, bw2, bw3), pch = 4, col = "red", cex = 1.5)
    }, ep = ep, poly_ord = poly_ord, season_p = season_p)
}

# Get and unwrap data from the file 
# "sim1_all_bwidths.csv"
data_bw <- unwrap()


## ----sim1fig1, fig.cap = "Boxplots of estimated bandwidths using $r=1$ for simulated quarterly data with $\\{e_{1,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp4_r1_e1, ep = 1, poly_ord = 1)


## ----sim1fig2, fig.cap = "Boxplots of estimated bandwidths using $r=1$ for simulated quarterly data with $\\{e_{2,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp4_r1_e2, ep = 2, poly_ord = 1)


## ----sim1fig3, fig.cap = "Boxplots of estimated bandwidths using $r=3$ for simulated quarterly data with $\\{e_{1,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp4_r3_e1, ep = 1, poly_ord = 3)


## ----sim1fig4, fig.cap = "Boxplots of estimated bandwidths using $r=3$ for simulated quarterly data with $\\{e_{2,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp4_r3_e2, ep = 2, poly_ord = 3)


## ----sim1fig5, fig.cap = "Boxplots of estimated bandwidths using $r=1$ for simulated monthly data with $\\{e_{1,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp12_r1_e1, ep = 1, poly_ord = 1)


## ----sim1fig6, fig.cap = "Boxplots of estimated bandwidths using $r=1$ for simulated monthly data with $\\{e_{2,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp12_r1_e2, ep = 2, poly_ord = 1)


## ----sim1fig7, fig.cap = "Boxplots of estimated bandwidths using $r=3$ for simulated monthly data with $\\{e_{1,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp12_r3_e1, ep = 1, poly_ord = 3)


## ----sim1fig8, fig.cap = "Boxplots of estimated bandwidths using $r=3$ for simulated monthly data with $\\{e_{2,t}\\}$. Red crosses mark the true $h_A$.", fig.align='center', fig.height = 8.8, fig.width = 5.8, echo = FALSE, fig.pos="H"----
boxplot_fun(data_bw$sp12_r3_e2, ep = 2, poly_ord = 3)


## ----echo = FALSE, eval = FALSE, incude = FALSE, purl = TRUE------------------
# 
# ##### Tables for second simulation study


## ----sim2tab1, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)

df <- read.csv("sim-data/sim2_error1_gamma01_quarterly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for quarterly data with error $\\{E_{1,t}\\}$ and $\\gamma = 0.1$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%
  landscape_tab()



## ----sim2tab2, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error2_gamma01_quarterly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for quarterly data with error $\\{E_{2,t}\\}$ and $\\gamma = 0.1$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab3, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error3_gamma01_quarterly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for quarterly data with error $\\{E_{3,t}\\}$ and $\\gamma = 0.1$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab4, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error1_gamma02_quarterly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for quarterly data with error $\\{E_{1,t}\\}$ and $\\gamma = 0.2$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()



## ----sim2tab5, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error2_gamma02_quarterly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for quarterly data with error $\\{E_{2,t}\\}$ and $\\gamma = 0.2$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab6, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error3_gamma02_quarterly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for quarterly data with error $\\{E_{3,t}\\}$ and $\\gamma = 0.2$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab7, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error1_gamma01_monthly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for monthly data with error $\\{E_{1,t}\\}$ and $\\gamma = 0.1$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()



## ----sim2tab8, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error2_gamma01_monthly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for monthly data with error $\\{E_{2,t}\\}$ and $\\gamma = 0.1$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab9, echo = FALSE---------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error3_gamma01_monthly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for monthly data with error $\\{E_{3,t}\\}$ and $\\gamma = 0.1$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab10, echo = FALSE--------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error1_gamma02_monthly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for monthly data with error $\\{E_{1,t}\\}$ and $\\gamma = 0.2$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab11, echo = FALSE--------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error2_gamma02_monthly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for monthly data with error $\\{E_{2,t}\\}$ and $\\gamma = 0.2$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()


## ----sim2tab12, echo = FALSE--------------------------------------------------
library(knitr) 
library(kableExtra)
df <- read.csv("sim-data/sim2_error3_gamma02_monthly.csv")  
if (knitr::is_html_output()) {
  df <- df %>% adjust_bold_html()
}

kable(df, col.names = c("Method", "Quantity", rep(c("$n_1$", "$n_2$", "$n_3$"), 6)), escape = FALSE, format = fmt, booktabs = TRUE, 
  linesep = c(rep(c("", "\\addlinespace"), 8), ""),
  align = "llrrrrrrrrrrrrrrrrrr", caption = "Results of the seasonality estimation simulation study for monthly data with error $\\{E_{3,t}\\}$ and $\\gamma = 0.2$. The simulated mean MSE (Mean) and the MSE standard deviation (SD), both multiplied by 1000, are displayed for various methods. The smallest MSE values in each column are highlighted.") %>%
  add_header_above(c(" ", " ", rep(c("$G_1$" = 3, "$G_2$" = 3), 3)), escape = FALSE) %>%
  add_header_above(c(" ", " ", "$S_1$" = 6, "$S_2$" = 6, "$S_3$" = 6), escape = FALSE) %>%
  kable_styling(latex_options = c("scale_down")) %>%
  adjust_html_tab_size() %>%  
  landscape_tab()

