---
title: "Replication Materials for hdMTD"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replication Materials for hdMTD}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library("hdMTD")
library("dplyr")
library("ggplot2")
library("lubridate")
library("purrr")
library("tidyr")

replication_dir <- system.file(
  "extdata",
  "replication",
  package = "hdMTD"
)

# Allows the vignette to be knitted directly from the package source.
if (!nzchar(replication_dir)) {
  possible_dirs <- c(
    file.path("inst", "extdata", "replication"),
    file.path("..", "inst", "extdata", "replication")
  )
  replication_dir <- possible_dirs[dir.exists(possible_dirs)][1]
}

if (is.na(replication_dir) || !dir.exists(replication_dir)) {
  stop("Pre-computed replication files were not found.")
}

precomputed <- readRDS(
  file.path(replication_dir, "hdMTD_outputs.rds")
)

simulated_data <- readRDS(
  file.path(replication_dir, "simulated_data.rds")
)

results_sequential <- readRDS(
  file.path(replication_dir, "results_sequential_selection.rds")
)
```



This vignette provides replication materials for the examples and analyses presented in Section 5 of the article describing the `hdMTD` package, available as a preprint on arXiv at [https://arxiv.org/abs/2509.01808](https://arxiv.org/abs/2509.01808). Computations that take more than approximately two minutes are shown but not evaluated when the vignette is built. Their outputs are loaded from pre-computed files distributed with the package.

## Section 5: Using hdMTD

### 5.1 Data generation

1. Generate MTD model:

Parameters: $\Lambda = \{-30,-15,-1\}$, $\mathcal{A} = \{0,1\}$,
$\lambda_0= \{0.01\}$, $\lambda_{-1} = 0.39$, $\lambda_{-15} = \lambda_{-30} = 0.3$,
$p_0(0)=p_0(1)=0.5$, and transition matrices $p_j$, $j\in\Lambda$, sampled uniformly.

```{r binary-model}
set.seed(11)
Lambda <- c(1, 15, 30)
A <- c(0, 1)
lam0 <- 0.01
lamj <- c(0.39, 0.3, 0.3)
p0 <- c(0.5, 0.5)
MTD <- MTDmodel(Lambda = Lambda, A = A, lam0 = lam0, lamj = lamj, p0 = p0)
summary(MTD)
```


2. Sample from the invariant distribution

```{r binary-sample}
X <- perfectSample(MTD, N = 1000)
```


### 5.2 Estimation

3. Estimate relevant lags using FS method

```{r fs-specific}
hdMTD_FS(X, d = 40, l = 4)
```


```{r fs-general}
FS <- hdMTD(X, d = 40, method = "FS", l = 4)
S(FS); summary(FS)
```


4. Estimate relevant lags using BIC method

The unrestricted BIC computation below takes approximately 30 minutes and is therefore displayed without being evaluated. The stored result is printed immediately afterwards.

```{r bic-unrestricted, eval=FALSE}
BIC_d40_l4 <- hdMTD_BIC(X, d = 40, minl = 4, maxl = 4)
BIC_d40_l4
```

```{r bic-unrestricted-output, echo=FALSE}
BIC_d40_l4 <- precomputed$BIC_d40_l4
print(BIC_d40_l4)
```


Custom subset S
```{r bic-custom-set}
hdMTD_BIC(X, d = 40,
          S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40),
          minl = 4, maxl = 4)
```


Varying number of lags to be selected
```{r bic-variable-size}
hdMTD_BIC(X, d = 40,
          S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40),
          minl = 1, maxl = 4)
```


With BIC values by number of lags
```{r bic-by-size}
hdMTD_BIC(X, d = 40,
          S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40),
          minl = 1, maxl = 4,
          byl = TRUE, BICvalue = TRUE)
```


Setting $\xi=0.4$
```{r bic-xi}
hdMTD_BIC(X, d = 40,
          S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40),
          minl = 1, maxl = 4,
          byl = TRUE, BICvalue = TRUE,
          xi = 0.4)
```


All matrices $p_j$ are equal and $\lambda_0=0$
```{r bic-single-matrix}
hdMTD_BIC(X, d = 40,
          S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40),
          minl = 1, maxl = 4,
          byl = TRUE, BICvalue = TRUE,
          single_matrix = TRUE, indep_part = FALSE)
```


5. Estimate relevant lags using CUT method

The first two CUT computations below take more than two minutes and are therefore displayed without being evaluated. Their stored results are printed immediately afterwards.

```{r cut-full-candidate-set, eval=FALSE}
CUT_d40 <- hdMTD_CUT(
  X,
  d = 40,
  S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40)
)
CUT_d40
```

```{r cut-full-candidate-set-output, echo=FALSE}
CUT_d40 <- precomputed$CUT_d40
print(CUT_d40)
```

Setting $\alpha = 0.13$

```{r cut-alpha, eval=FALSE}
CUT_d40_alpha <- hdMTD_CUT(
  X,
  d = 40,
  S = c(1, 5, 10, 15, 17, 20, 27, 30, 35, 40),
  alpha = 0.13
)
CUT_d40_alpha
```

```{r cut-alpha-output, echo=FALSE}
CUT_d40_alpha <- precomputed$CUT_d40_alpha
print(CUT_d40_alpha)
```


Custom subset S
```{r cut-custom-set}
hdMTD_CUT(X, d = 40, S = c(1, 5, 17, 27, 30, 35), alpha = 0.13)
```


6. Estimate relevant lags using FSC method

```{r fsc}
hdMTD_FSC(X, d = 40, l = 4, alpha = 0.1, cut_fraction = 0.5)
```


FS method with halved sample
```{r fs-half-sample}
hdMTD_FS(X[501:1000], d = 40, l = 4)
```


7. Estimating transition probabilities

```{r empirical-probabilities-data-frame}
head(empirical_probs(X, S = c(1, 15, 30)), 6)
```


```{r empirical-probabilities-matrix}
empirical_probs(X, S = c(1, 15, 30), matrixform = TRUE)
```


8. Oscillations

Computing from MTD
```{r true-oscillations}
oscillation(MTD)
```

Estimating from sample
```{r estimated-oscillations}
oscillation(X, S = c(1, 15, 30))
```


9. Estimating MTD parameters through the EM algorithm

Initial parameters for EM method
```{r em-initial-values}
init <- list(
  'lambdas'= c(0.01, 0.33, 0.33, 0.33),
  'p0' = c(0.5, 0.5),
  'pj' = rep(list(matrix(c(0.5, 0.5, 0.5, 0.5), ncol = 2, nrow = 2)), 3)
)
```

Run EM
```{r em-fit}
emMTD <- MTDest(X, S = c(1, 15, 30), init = init, iter = TRUE)
summary(emMTD)
```

Stops after $9$ iterations
```{r em-fixed-iterations}
emMTD <- MTDest(X, S = c(1, 15, 30), M = NULL, nIter = 9, init = init, oscillations = TRUE)
summary(emMTD)
```


Compute estimated global transition matrix

```{r unnamed-chunk-26}
transitP(emMTD)
```

10. A three-state MTD example

All examples presented so far have used the binary state space $\mathcal{A}=\{0,1\}$. We now illustrate the use of the package with a three-state MTD model.

```{r three-state-model}
set.seed(11)
Lambda3 <- c(2, 5)
A3 <- c(1, 2, 3)
lam0.3 <- 0.05
lamj.3 <- c(0.20, 0.75)
p0.3 <- c(0.30, 0.30, 0.40)
MTD3 <- MTDmodel(
  Lambda = Lambda3,
  A = A3,
  lam0 = lam0.3,
  lamj = lamj.3,
  p0 = p0.3
)
summary(MTD3)
```

```{r three-state-sample}
Y <- perfectSample(MTD3, N = 2000)
oscillation(MTD3)
```

```{r three-state-fs}
Sfs3 <- hdMTD_FS(Y, d = 20, l = 3)
Sfs3
```

```{r three-state-bic}
Sbic3 <- hdMTD_BIC(Y, S = Sfs3, d = 20)
Sbic3
```

```{r three-state-probabilities}
P3 <- empirical_probs(Y, S = Sbic3, matrixform = TRUE)
P3
```


### 5.3 Testing hdMTD

1. MTD model specification:

```{r simulation-model}
set.seed(123)
Lambda <- c(1, 5)
A <- c(0, 1)
lam0 <- 0.01
p0 <- c(0.5, 0.5)
MTD <- MTDmodel(Lambda, A, lam0, p0 = p0) # Generates an MTD model
```


2. Simulation settings and results

Simulation parameters:
```{r simulation-settings}
n <- 100      # Number of replications
N <- 10000    # Full sample size
m <- c(1000, 1500, 2000, 2500, 3000, 5000, 10000) # Subsample sizes
d <- 100      # Max order for FS and Oracle
dNaive <- 5   # Max order for Naive
pairList <- t(combn(d, 2)) # All possible pairs with digits from 1 to 100
npairs <- nrow(pairList)
minP11_P12 <- min(MTD$P[1, 1], MTD$P[1, 2])
```



The following code reproduces the complete simulation study. It is not evaluated when the vignette is built because the computation may take approximately one day, depending on the number of workers.

```{r simulation-recompute, eval=FALSE}
FS  <- matrix(0, ncol = length(m), nrow = n)
Naive  <- matrix(0, ncol = length(m), nrow = n)
Oracle  <- matrix(0, ncol = length(m), nrow = n)
SFS     <- matrix(0, ncol = length(m) * 2, nrow = n)
ZOracle <- matrix(0, ncol = length(m) * 2, nrow = n)

X_list <- vector("list", n)
for (i in seq_len(n)) {
  X_list[[i]] <- perfectSample(MTD, N = N) # Generates n samples of size N
}

# Parallel
# Limit internal BLAS/OpenMP threads to avoid oversubscription across worker
Sys.setenv(OMP_NUM_THREADS = "1",
           MKL_NUM_THREADS = "1",
           OPENBLAS_NUM_THREADS = "1",
           BLIS_NUM_THREADS = "1")

# Detect cores and choose a safe number of workers
get_smart_workers <- function() {
  physical_cores <- parallel::detectCores(logical = FALSE)

  if (physical_cores >= 10) {
    return(6)
  } else if (physical_cores >= 6) {
    return(4)
  } else if (physical_cores >= 4) {
    return(2)
  } else {
    return(1)
  }
}

# Manually assign number of workers (n_workers) or leave NULL for automatic choice.
n_workers <- NULL
workers <- ifelse(is.null(n_workers), get_smart_workers(), n_workers)
message(sprintf("Using %d workers", workers))
future::plan(future::multisession, workers = workers)
op <- options(future.scheduling = 1)

# Parallelize over replications
one_rep <- function(X) {
  FS_i <- Naive_i <- Oracle_i <- numeric(length(m))
  SFS_row <- ZOracle_row <- integer(length(m) * 2)

  for (k in seq_along(m)) {
    Y <- X[seq_len(m[k])]
    ct <- countsTab(Y, d = d)

    # FS
    S <- hdMTD_FS(Y, d = d, l = 2)
    SFS_row[(k * 2 - 1):(k * 2)] <- S
    p_FS <- freqTab(S = S, A = A, countsTab = ct)$qax_Sj[1]
    FS_i[k]  <- abs(p_FS - MTD$P[1, 1])

    # Naive
    ct_dNaive <- countsTab(Y, dNaive)
    p_Naive <- freqTab(S = seq_len(dNaive), A = A, countsTab = ct_dNaive)$qax_Sj[1]
    Naive_i[k]  <- abs(p_Naive - MTD$P[1, 1])

    # Oracle
    p_pairs <- numeric(npairs)
    for (s in seq_len(npairs)) {
      p_pairs[s] <- freqTab(S = pairList[s, ], A = A, countsTab = ct)$qax_Sj[1]
    }
    minpos <- which.min(abs(p_pairs - MTD$P[1, 1]))
    ZOracle_row[(k * 2 - 1):(k * 2)] <- pairList[minpos, ]
    p_Oracle <- p_pairs[minpos]
    Oracle_i[k]  <- abs(p_Oracle - MTD$P[1, 1])
  }

  list(FS = FS_i, Naive = Naive_i, Oracle = Oracle_i,
       SFS = SFS_row, ZOracle = ZOracle_row)
}

indices_to_run <- seq_len(n)
chunks <- split(indices_to_run, ceiling(seq_along(indices_to_run) / workers))

# Note that, given X_list the following loop is deterministic
for (bi in seq_along(chunks)) { #takes ~ceiling(100/workers)*90 minutes to run)
  rows <- chunks[[bi]]
  res_block <- future.apply::future_lapply(X_list[rows], one_rep, future.seed = TRUE)

  #Fill block lines
  FS[rows, ]      <- do.call(rbind, lapply(res_block, `[[`, "FS"))
  Naive[rows, ]   <- do.call(rbind, lapply(res_block, `[[`, "Naive"))
  Oracle[rows, ]  <- do.call(rbind, lapply(res_block, `[[`, "Oracle"))
  SFS[rows, ]     <- do.call(rbind, lapply(res_block, `[[`, "SFS"))
  ZOracle[rows, ] <- do.call(rbind, lapply(res_block, `[[`, "ZOracle"))
}

future::plan(future::sequential)
options(op)

FSP <- FS/minP11_P12
OracleP <- Oracle/minP11_P12
NaiveP <- Naive/minP11_P12

# Store all simulated data in a list
simulated_data <- list(
  FS = FS, FSP = FSP,
  Oracle = Oracle, OracleP = OracleP,
  Naive = Naive, NaiveP = NaiveP,
  SFS = SFS, ZOracle = ZOracle
)
```

```{r simulation-load, echo=FALSE}
FS <- simulated_data$FS
FSP <- simulated_data$FSP
Oracle <- simulated_data$Oracle
OracleP <- simulated_data$OracleP
Naive <- simulated_data$Naive
NaiveP <- simulated_data$NaiveP
SFS <- simulated_data$SFS
ZOracle <- simulated_data$ZOracle
```


3. Generate Table 1
```{r unnamed-chunk-31}
means_table <- rbind(
  round(apply(FS, 2, mean), 5),
  round(apply(Oracle, 2, mean), 5),
  round(apply(Naive, 2, mean), 5),
  round(apply(FSP, 2, mean), 5),
  round(apply(OracleP, 2, mean), 5),
  round(apply(NaiveP, 2, mean), 5)
)
rownames(means_table) <- c(
  "Delta_FS(m)",
  "Delta_Oracle(m)",
  "Delta_Naive,5(m)",
  "std Delta_FS(m)",
  "std Delta_Oracle(m)",
  "std Delta_Naive,5(m)"
)
colnames(means_table) <- m
```


Table 1: Mean error of estimators
```{r mean-error-table, results='asis', message=FALSE, warning=FALSE}
suppressWarnings(
  tbl_1 <- knitr::kable(
    means_table,
    caption = "Mean error of estimators",
    format = "markdown",
    booktabs = TRUE,
    escape = FALSE
  )
)
suppressWarnings(print(tbl_1))
```


4. Compute how often the FS output differs from Oracle by subsample size
```{r unnamed-chunk-33}
m_index <- seq(1, length(m)*2, by = 2)
names(m_index) <- as.character(m)
SFS_vs_ZOracle_diff <- sapply(m_index, function(idx) {
  SFS_set <- apply(SFS[, c(idx, idx + 1)], 1, function(x) paste(sort(x), collapse = "-"))
  ZOracle_set <- apply(ZOracle[, c(idx, idx + 1)], 1, function(x) paste(sort(x), collapse = "-"))
  sum(SFS_set != ZOracle_set)
})
SFS_vs_ZOracle_diff
```


5. Generate Figure 1:

Data arrangement

FS
```{r unnamed-chunk-34}
tab <- FS
FStab <- rbind(apply(tab, 2, summary),'sd'=apply(tab, 2, sd))
FStab <- rbind(FStab,'sdLo'=FStab[4,]-FStab[7,],'sdUp'=FStab[4,]+FStab[7,])
Fmean <- FStab[4,]
FsdUp <- FStab[9,]
FsdLo <- FStab[8,]
Fq1 <- FStab[2,]
Fq2 <- FStab[3,]
Fq3 <- FStab[5,]
```


NAIVE
```{r unnamed-chunk-35}
tab <- Naive
Naivetab <- rbind(apply(tab, 2, summary),'sd'=apply(tab, 2, sd))
Naivetab <- rbind(Naivetab,'sdLo'=Naivetab[4,]-Naivetab[7,],'sdUp'=Naivetab[4,]+Naivetab[7,])
Nmean <- Naivetab[4,]
NsdUp <- Naivetab[9,]
NsdLo <- Naivetab[8,]
Nq1 <- Naivetab[2,]
Nq2 <- Naivetab[3,]
Nq3 <- Naivetab[5,]
```


ORACLE
```{r unnamed-chunk-36}
tab <- Oracle
Oracletab <- rbind(apply(tab, 2, summary),'sd'=apply(tab, 2, sd))
Oracletab <- rbind(Oracletab,'sdLo'=Oracletab[4,]-Oracletab[7,],'sdUp'=Oracletab[4,]+Oracletab[7,])
Omean <- Oracletab[4,]
OsdUp <- Oracletab[9,]
OsdLo <- Oracletab[8,]
Oq1 <- Oracletab[2,]
Oq2 <- Oracletab[3,]
Oq3 <- Oracletab[5,]
```


### Plot Figure 1: Estimators mean error across $N_{rep}=100$ replications.
```{r fig1-plot, fig.width=10, fig.height=6, message=FALSE, warning=FALSE}
 par(mfrow = c(1,2), oma = c(0,0,0,0))

 ## --- Left panel: Mean error with standard deviation bands ---
 par(mar = c(5,5,3,4), xpd = NA)

 plot(m/100, Fmean, type = "l", col = "#377EB8",
      xlab = "m (x100)", ylab = "Mean error", ylim = c(0, 0.12), lwd = 3,
      frame.plot = FALSE, xaxt = "n", yaxt = "n", xlim = c(10,100),
      cex.axis = 1.4, cex.lab = 1.6)
 lines(m/100, Omean, col = "#E41A1C", lwd = 3)
 lines(m/100, Nmean, col = "#4DAF4A", lwd = 3)
 points(m/100, Fmean, col = "#377EB8", pch = 19, cex = 0.7)
 points(m/100, Omean, col = "#E41A1C", pch = 19, cex = 0.7)
 points(m/100, Nmean, col = "#4DAF4A", pch = 19, cex = 0.7)
 lines(m/100, FsdUp, col = "#377EB8", lty = 2)
 lines(m/100, FsdLo, col = "#377EB8", lty = 2)
 lines(m/100, OsdUp, col = "#E41A1C", lty = 2)
 lines(m/100, OsdLo, col = "#E41A1C", lty = 2)
 cap <- 0.12
 NsdUp_cut <- ifelse(NsdUp > cap, NA, NsdUp)
 lines(m/100, NsdUp_cut, col = "#4DAF4A", lty = 2)
 lines(m/100, NsdLo, col = "#4DAF4A", lty = 2)
 axis(side = 1, at = m/100, labels = m/100, cex.axis = 1.4)
 axis(side = 2, cex.axis = 1.4)

 legend(
   "topright",
   inset = c(0.05, 0),
   legend = c(expression(bar(Delta) ~ "FS"),
              expression(bar(Delta) ~ "FS" %+-% "sd"),
              expression(bar(Delta) ~ "Oracle"),
              expression(bar(Delta) ~ "Oracle" %+-% "sd"),
              expression(bar(Delta) ~ "Naive"),
              expression(bar(Delta) ~ "Naive" %+-% "sd")),
   col = c("#377EB8","#377EB8","#E41A1C","#E41A1C","#4DAF4A","#4DAF4A"),
   lty = c(1,2,1,2,1,2),
   lwd = c(3,1.5,3,1.5,3,1.5), bty = "n",
   y.intersp = 0.95, x.intersp = 0.5, seg.len = 3, cex = 1.4
 )

 ## --- Right panel: Median and quartiles ---
 par(mar = c(5,5,3,4), xpd = NA)

 plot(m/100, Fq2, type = "l", col = "#377EB8",
      xlab = "m (x100)", ylab = "Quartiles of mean error",
      ylim = c(0, 0.12), lwd = 3, frame.plot = FALSE, xaxt = "n", yaxt = "n",
      xlim = c(10,100), cex.axis = 1.4, cex.lab = 1.6)
 lines(m/100, Oq2, col = "#E41A1C", lwd = 3)
 lines(m/100, Nq2, col = "#4DAF4A", lwd = 3)
 points(m/100, Fq2, col = "#377EB8", pch = 19, cex = 0.7)
 points(m/100, Oq2, col = "#E41A1C", pch = 19, cex = 0.7)
 points(m/100, Nq2, col = "#4DAF4A", pch = 19, cex = 0.7)
 lines(m/100, Fq1, col = "#377EB8", lty = 2)
 lines(m/100, Fq3, col = "#377EB8", lty = 2)
 lines(m/100, Oq1, col = "#E41A1C", lty = 2)
 lines(m/100, Oq3, col = "#E41A1C", lty = 2)
 lines(m/100, Nq1, col = "#4DAF4A", lty = 2)
 lines(m/100, Nq3, col = "#4DAF4A", lty = 2)
 axis(side = 1, at = m/100, labels = m/100, cex.axis = 1.4)
 axis(side = 2, cex.axis = 1.4)

 legend(
   "topright",
   inset = c(0.09, 0),
   legend = c(expression("Med  " ~ bar(Delta) ~ "FS"),
              expression("q1,q3" ~ bar(Delta) ~ "FS"),
              expression("Med  " ~ bar(Delta) ~ "Oracle"),
              expression("q1,q3" ~ bar(Delta) ~ "Oracle"),
              expression("Med  " ~ bar(Delta) ~ "Naive"),
              expression("q1,q3" ~ bar(Delta) ~ "Naive")),
   col = c("#377EB8","#377EB8","#E41A1C","#E41A1C","#4DAF4A","#4DAF4A"),
   lty = c(1,2,1,2,1,2),
   lwd = c(3,1.5,3,1.5,3,1.5), bty = "n",
   y.intersp = 0.95, x.intersp = 0.5, seg.len = 3, cex = 1.4
 )
```


### 5.4 Analysis of Real-World Data

```{r unnamed-chunk-38}
data("tempdata")
```


1. Treat NA data:

Removing days before "2010-08-05"
```{r unnamed-chunk-39}
tempdata <- hdMTD::tempdata %>% filter(DATE >= "2010-08-05")
```

Identify remaining $155$ NA positions
```{r unnamed-chunk-40}
posNA <- which(is.na(tempdata$MAXTEMP))
```

Fill short sequences of NAs ( $≤6$ NAs) using nearest neighbors mean
```{r unnamed-chunk-41}
for (i in posNA) {
  if(!is.na(tempdata$MAXTEMP[i - 1]) && !all(is.na(tempdata$MAXTEMP[(i + 1):(i + 6)]))) {
    aux <- which(!is.na(tempdata$MAXTEMP[(i + 1):(i + 6)]))[1]
    tempdata$MAXTEMP[i] <- mean(tempdata$MAXTEMP[i - 1],tempdata$MAXTEMP[i + aux])
  }
}
```

Identify remaining $88$ NA positions
```{r unnamed-chunk-42}
posNA <- which(is.na(tempdata$MAXTEMP))
```

Fill remaining NA with the mean of the previous hour temperature, next hour
temperature, and same hour of previous day temperature.
```{r unnamed-chunk-43}
for (i in posNA) {
  tempdata$MAXTEMP[i] <- mean(c(tempdata$MAXTEMP[i - 1], tempdata$MAXTEMP[i + 1],
                                tempdata$MAXTEMP[i - 24]),na.rm = TRUE)
}
```


2. Compute mean daily maximum temperatures:

```{r unnamed-chunk-44}
temp <- tempdata %>%
  group_by(DATE) %>%
  summarize(MAXTEMP = mean(MAXTEMP), .groups = 'drop')
head(temp, 4)
```


### Plot Figure 2: Time series with quarterly mean of daily maximum temperatures
```{r temp-plot, fig.width=9, fig.height=6, message=FALSE, warning=FALSE}
 TRIM_DATA <- temp %>%
   mutate(
     Y_TRIMESTER = paste0(year(DATE), "-T", quarter(DATE))
   ) %>%
   group_by(Y_TRIMESTER) %>%
   summarise(
     MEAN_TEMP = mean(MAXTEMP),
     DATA_REF = min(DATE)
   ) %>%
   ungroup() %>%
   arrange(DATA_REF)
 TRIM_DATA <- TRIM_DATA[-c(1, nrow(TRIM_DATA)),]

 ggplot(TRIM_DATA, aes(x = DATA_REF, y = MEAN_TEMP)) +
   geom_line(color = "steelblue", linewidth = 0.5) +
   geom_point(color = "steelblue", size = 1.2) +
   scale_x_date(
     date_breaks = "1 year",
     date_labels = "%Y",
     minor_breaks = NULL
   ) +
   theme_minimal() +
   theme(
     axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
     axis.text.y = element_text(size = 12),
     panel.grid.major = element_line(color = "gray90"),
     axis.title.x = element_text(size = 15),
     axis.title.y = element_text(size = 15),
     plot.title = element_text(hjust = 0.5, size = 17, face = "bold")
   ) +
   labs(
     title = "Quarterly mean of daily maximum temperatures across the years",
     x = "Year",
     y = "Mean Temperature (°C)"
   )
```


3. Create categories of temperature:

```{r unnamed-chunk-46}
xn <- max(temp$MAXTEMP)
x1 <- min(temp$MAXTEMP)
maxAmp <- xn - x1
temp$MAXTEMP1  <- ifelse(temp$MAXTEMP < x1 + maxAmp/2, 1, 2)
head(temp, 4)
prop.table(table(temp$MAXTEMP1)) # frequency of thermal regimes
```


4. Run FS for temp:

The lag-selection calls in this subsection take several minutes. Their code is displayed, while the vignette uses the corresponding pre-computed results.

```{r unnamed-chunk-47}
Temp12 <- rev(temp$MAXTEMP1)
```

hdMTD functions assume the sample is sorted from the latest observations to
oldest.

```{r temperature-fs-d400, eval=FALSE}
FS_Temp12_d400 <- hdMTD_FS(Temp12, d = 400, l = 3)
FS_Temp12_d400
```

```{r temperature-fs-d400-output, echo=FALSE}
FS_Temp12_d400 <- precomputed$FS_Temp12_d400
print(FS_Temp12_d400)
```


Note: The next code line is mentioned in the article but without a CodeChunk <br>
Reduce maximum order to improve estimation <br>
```{r temperature-fs-d364, eval=FALSE}
FS_Temp12_d364 <- hdMTD_FS(Temp12, d = 364, l = 3)
FS_Temp12_d364
```

```{r temperature-fs-d364-output, echo=FALSE}
FS_Temp12_d364 <- precomputed$FS_Temp12_d364
print(FS_Temp12_d364)
```


5. Split sample in Train and Test data:

```{r unnamed-chunk-52}
ndays <- nrow(temp %>%
                filter(DATE >= "2023-09-01")) # 366 days in the latest year of the sample
Temp12_Train <- Temp12[-seq_len(ndays)] # Training data with 4775 days
Temp12_Test <- Temp12[seq_len(ndays)] # Test data
```


6. Rerun FS for Train data:

```{r temperature-train-fs, eval=FALSE}
FS_Temp12Train <- hdMTD_FS(Temp12_Train, d = 364, l = 3)
FS_Temp12Train
```

```{r temperature-train-fs-output, echo=FALSE}
FS_Temp12Train <- precomputed$FS_Temp12Train
print(FS_Temp12Train)
```

7. Trim out irrelevant lags:

With CUT method
```{r unnamed-chunk-55}
hdMTD_CUT(Temp12_Train, d = 364, S = c(1, 364, 6))
```

With BIC method
```{r unnamed-chunk-56}
hdMTD_BIC(Temp12_Train, d = 364, S = c(1, 364, 6), minl = 1, maxl = 3,
          byl = TRUE, BICvalue = TRUE )
```

8. Estimated transition matrix for FS method output:

```{r unnamed-chunk-59}
P_FS <- empirical_probs(Temp12_Train, S = c(1, 6, 364), matrixform = TRUE)
P_FS
```


### Classic method for choosing relevant lag set:

9. Compute models:

```{r unnamed-chunk-60}
ct <- countsTab(Temp12_Train, d = 6) # Table with size 6 sequence counts
head(ct,4)
```


MC1 (does not explicitly appear in the article)
```{r unnamed-chunk-61}
ft <- freqTab(S = 1, A = c(1, 2), countsTab = ct)
LL <- sum(log(ft$qax_Sj) * ft$Nxa_Sj)
freeParam <- 2 * 1
BICMC1 <- -LL + 0.5 * log(length(Temp12_Train)) * freeParam
BICMC1
```

Comparable BIC if the model is a Markov chain of order $1$: $1869.162$

MC2 (does not explicitly appear in the article)
```{r unnamed-chunk-62}
ft <- freqTab(S = c(1, 2), A = c(1, 2), countsTab = ct)
LL <- sum(log(ft$qax_Sj) * ft$Nxa_Sj)
freeParam <- 2^2 * 1
BICMC2 <- -LL + 0.5 * log(length(Temp12_Train)) * freeParam
BICMC2
```

Comparable BIC if the model is a Markov chain of order $2$: $1850.598$

MC3. Used as example in the article
```{r unnamed-chunk-63}
ft <- freqTab(S = c(1, 2, 3), A = c(1, 2), countsTab = ct)
head(ft, 4)
LL <- sum(log(ft$qax_Sj) * ft$Nxa_Sj)
freeParam <- 2^3 * 1
BICMC3 <- -LL + 0.5 * log(length(Temp12_Train)) * freeParam
BICMC3
```

Comparable BIC if the model is a Markov chain of order $3$: $1854.029$

MC4 (does not explicitly appear in the article)
```{r unnamed-chunk-64}
ft <- freqTab(S = c(1, 2, 3, 4), A = c(1, 2), countsTab = ct)
LL <- sum(log(ft$qax_Sj) * ft$Nxa_Sj)
freeParam <- 2^4 * 1
BICMC4 <- -LL + 0.5 * log(length(Temp12_Train)) * freeParam
BICMC4
```

Comparable BIC if the model is a Markov chain of order $4$: $1877.888$

MC5 (does not explicitly appear in the article)
```{r unnamed-chunk-65}
ft <- freqTab(S = c(1, 2, 3, 4, 5), A = c(1, 2), countsTab = ct)
pos <- which(ft$Nxa_Sj > 0)
LL <- sum(log(ft$qax_Sj[pos]) * ft$Nxa_Sj[pos])
freeParam <- 2^5 * 1
BICMC5 <- -LL + 0.5 * log(length(Temp12_Train)) * freeParam
BICMC5
```

Comparable BIC if the model is a Markov chain of order $5$: $1925.962$

MC6 (does not explicitly appear in the article)
```{r unnamed-chunk-66}
ft <- freqTab(S = c(1, 2, 3, 4, 5, 6), A = c(1, 2), countsTab = ct)
pos <- which(ft$Nxa_Sj > 0)
LL <- sum(log(ft$qax_Sj[pos]) * ft$Nxa_Sj[pos])
freeParam <- 2^6 * 1
BICMC6 <- -LL + 0.5 * log(length(Temp12_Train)) * freeParam
BICMC6
```

Comparable BIC if the model is a Markov chain of order $6$: $2031.679$

10. Comparing models:

```{r unnamed-chunk-67}
BIC_vals <-  c(BICMC1, BICMC2, BICMC3, BICMC4, BICMC5, BICMC6)
model_names <- paste0("MC", 1:6)
minBIC_idx <- which.min(BIC_vals)
BIC_fmt <- sprintf("%.3f", BIC_vals)
BIC_fmt[minBIC_idx] <- paste0("**", BIC_fmt[minBIC_idx], "**")
minBIC_idx
```

The classic method chooses order $2$ ($S=\{-2,-1\}$).

Generating BIC values Table
```{r unnamed-chunk-68}
bic_matrix <- data.frame(matrix(BIC_fmt, nrow = 1))
colnames(bic_matrix) <- model_names
rownames(bic_matrix) <- "BIC"
```

Table 2: BIC values
```{r bic-table, results='asis', message=FALSE, warning=FALSE}
tbl_2 <- suppressWarnings(
  knitr::kable(
    bic_matrix,
    caption = "BIC values computed for classical Markov chain models of different orders."
  )
)
suppressWarnings(print(tbl_2))
```


Estimated matrix for $S=\{-2,-1\}$
```{r unnamed-chunk-70}
P_MC2 <- empirical_probs(Temp12_Train, S = c(1, 2), matrixform = TRUE)
P_MC2
```

Independent model distribution:
```{r unnamed-chunk-71}
P_Ind <- prop.table(table(Temp12_Train))
P_Ind
```


### Comparing methods

11. Computing values for Table 3.

```{r unnamed-chunk-72}
Days1 <- which(Temp12_Test == 1)
lenDays1 <- length(Days1)
lenDays1 / ndays # frequency of low temperature days
Temp12_Test <- c(Temp12_Test, Temp12_Train[seq_len(364)])
```


```{r unnamed-chunk-73}
set.seed(1)
hitInd <- numeric(1000)
hitMC2 <- numeric(1000)
hitFS <- numeric(1000)
T1Ind <- numeric(1000)
T1MC2 <- numeric(1000)
T1FS <- numeric(1000)
F1Ind <- numeric(1000)
F1MC2 <- numeric(1000)
F1FS <- numeric(1000)
pasts2 <- rownames(P_MC2)
pastsFS <- rownames(P_FS)
```


```{r unnamed-chunk-74}
for (j in seq_len(1000)){
  u <- runif(ndays)
  predInd <- numeric(ndays)
  predMC2 <- numeric(ndays)
  predFS <- numeric(ndays)

  for (i in ndays:1) {
    predInd[i] <- ifelse(u[i] <= P_Ind[1], 1, 2)
    pastRow <- which(pasts2 == paste0(rev(Temp12_Test[c(i+1, i+2)]), collapse = ""))
    predMC2[i] <- ifelse(u[i] <= P_MC2[pastRow, 1], 1, 2)
    pastRow <- which(pastsFS == paste0(rev(Temp12_Test[c(i+1, i+6, i+364)]), collapse = ""))
    predFS[i] <- ifelse(u[i] <= P_FS[pastRow, 1], 1, 2)
  }

  hitInd[j] <- sum(predInd == Temp12_Test[seq_len(ndays)])
  hitMC2[j] <- sum(predMC2 == Temp12_Test[seq_len(ndays)])
  hitFS[j] <- sum(predFS == Temp12_Test[seq_len(ndays)])
  T1Ind[j] <- sum(predInd[Days1] == 1)
  T1MC2[j] <- sum(predMC2[Days1] == 1)
  T1FS[j] <- sum(predFS[Days1] == 1)
  F1Ind[j] <- sum(predInd[-Days1] == 1)
  F1MC2[j] <- sum(predMC2[-Days1] == 1)
  F1FS[j] <- sum(predFS[-Days1] == 1)
}
```

### Accuracy
```{r unnamed-chunk-75}
AccInd <- mean(hitInd)/ndays
AccMC2 <- mean(hitMC2)/ndays
AccFS <- mean(hitFS)/ndays
AccInd; AccMC2; AccFS
```

### Precision
```{r unnamed-chunk-76}
PrecInd <- mean(T1Ind/(T1Ind + F1Ind))
PrecMC2 <- mean(T1MC2/(T1MC2 + F1MC2))
PrecFS <- mean(T1FS/(T1FS + F1FS))
PrecInd; PrecMC2; PrecFS
```

### Sensitivity (Recall)
```{r unnamed-chunk-77}
SensInd <- mean(T1Ind)/lenDays1
SensMC2 <- mean(T1MC2)/lenDays1
SensFS <- mean(T1FS)/lenDays1
SensInd; SensMC2; SensFS
```

### Specificity
```{r unnamed-chunk-78}
SpecInd <- 1 - mean(F1Ind)/(ndays - lenDays1)
SpecMC2 <- 1 - mean(F1MC2)/(ndays - lenDays1)
SpecFS <- 1 - mean(F1FS)/(ndays - lenDays1)
SpecInd; SpecMC2; SpecFS
```

### F1-Score
```{r unnamed-chunk-79}
F1ScoreInd <- 2 * (PrecInd * SensInd) / (PrecInd + SensInd)
F1ScoreMC2 <- 2 * (PrecMC2 * SensMC2) / (PrecMC2 + SensMC2)
F1ScoreFS <- 2 * (PrecFS * SensFS) / (PrecFS + SensFS)
F1ScoreInd; F1ScoreMC2; F1ScoreFS

metric <- c("Accuracy", "Precision", "Sensitivity (Recall)", "Specificity", "F1-Score")
formula <- c("(TP+TN)/(TP+TN+FP+FN)",
             "TP/(TP+FP)",
             "TP/(TP+FN)",
             "TN/(TN+FP)",
             "2(PPV*Recall)/(PPV+Recall)")
performance_table <- data.frame(
  Metric = metric,
  Formula = formula,
  indc = round(c(AccInd, PrecInd, SensInd, SpecInd, F1ScoreInd) * 100, 2),
  mc2c = round(c(AccMC2, PrecMC2, SensMC2, SpecMC2, F1ScoreMC2) * 100, 2),
  fsc = round(c(AccFS, PrecFS, SensFS, SpecFS, F1ScoreFS) * 100, 2),
  check.names = FALSE
)
names(performance_table) <- c("Metric", "Formula", "Ind (\\%)", "MC2 (\\%)", "FS (\\%)")
```

Table 3: Model performance metrics
```{r performance-table, results='asis', message=FALSE, warning=FALSE}
suppressWarnings(
  tbl_3 <- knitr::kable(performance_table, align = "l",
                        caption = "Model performance metrics.")
)
suppressWarnings(print(tbl_3))
```


### Plot Figure 3: Exploratory analysis of accuracies

```{r accuracy-plot, fig.width=9, fig.height=6, message=FALSE, warning=FALSE}
 accuracy_data <- data.frame(
   MC2 = hitMC2 / ndays,
   FS = hitFS / ndays
 ) %>%
   pivot_longer(
     everything(),
     names_to = "Model",
     values_to = "Accuracy"
   )

 ggplot(accuracy_data, aes(x = Model, y = Accuracy, fill = Model)) +
   geom_boxplot() +
   labs(
     title = "Accuracy distribution (1000 replications)",
     x = "Model",
     y = "Accuracy"
   ) +
   theme_minimal() +
   scale_fill_brewer(palette = "Paired") +
   theme(
     plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
     axis.title = element_text(size = 16),
     axis.text = element_text(size = 14),
     legend.title = element_text(size = 16),
     legend.text = element_text(size = 14)
   )
```


### Empirical $\nu$ Analysis

12. FS sequential selection based on $\hat{\nu}_{n,j,S}$ values:

```{r sequential-selection-function, eval=FALSE}
run_sequential_lag_selection <- function(Temp12_Train, d = 364) {
  # Initialization
  A <-  sort(unique(Temp12_Train))
  lenA <- length(A)
  lenX <- length(Temp12_Train)
  A_pairs <- matrix(A, ncol = 2) # All unique state pairs
  ct <- countsTab(X = Temp12_Train, d = 364) # Sequence counts table

  # Initialize storage
  results <- list(
    nuj1 = numeric(d),
    nuj2 = numeric(d-1),
    nuj3 = numeric(d-2),
    selected_lags = numeric(3)
  )

  # Helper function for empirical distribution calculation
  PI <- function(S, groupTab, x_S, lenX, d) {
    if (length(S) > 0) {
      filtr_S <- paste0("x", S)
      groupTab <- groupTab %>%
        dplyr::mutate(match = purrr::pmap_lgl(dplyr::pick(dplyr::all_of(filtr_S)),
                                              ~ all(c(...) == x_S))) %>%
        dplyr::filter(match) %>%
        dplyr::select(-match)
    }
    PI <- matrix(groupTab$Nx_Sj/(lenX - d),ncol = 1)
    PI
  }

  # Sequential Selection Process
  cat("=== Starting Sequential Lag Selection ===\n")

  # Step 1: Initial selection (S = ∅)
  cat("\n[Step 1] Selecting first lag (S = ∅)...\n")
  S <- NULL
  Sc <- sort(setdiff(seq_len(d), S), decreasing = TRUE) # Complement of S in 1:d

  for (z in seq_along(Sc)) { # Runs across all available lags
    j <- Sc[z]
    # Frequency tables
    b_Sja <- freqTab(S = S, j = j, A = A, countsTab = ct)
    b_Sj <- b_Sja %>%
      dplyr::group_by(dplyr::across(dplyr::all_of(paste0("x", j)))) %>%
      dplyr::summarise(Nx_Sj = sum(Nxa_Sj), .groups = "drop")

    # Compute νj
    PIs <- PI(S = S, groupTab = b_Sj, x_S = S, lenX = lenX, d = 364)
    dTVs <- dTV_sample(S = S, j = j, lenA = lenA, base = b_Sja,
                       A_pairs = A_pairs, x_S = S)
    results$nuj1[z] <- prod(PIs) * dTVs
  }
  results$selected_lags[1] <- Sc[which.max(results$nuj1)]
  cat(sprintf("Selected: j = %d (ν = %.4f)\n",
              results$selected_lags[1], max(results$nuj1)))

  # Step 2: Second selection (S = {1})
  cat(sprintf("\n[Step 2] Selecting second lag (S = {%d})...\n",
              results$selected_lags[1]))
  S <- results$selected_lags[1]
  Sc <- sort(setdiff(seq_len(d), S), decreasing = TRUE)

  for (z in seq_along(Sc)) {
    j <- Sc[z]
    Sj <- sort(c(S, j), decreasing = TRUE)

    # Frequency tables
    b_Sja <- freqTab(S = S, j = j, A = A, countsTab = ct)
    b_Sj <- b_Sja %>%
      dplyr::group_by(dplyr::across(dplyr::all_of(paste0("x", Sj)))) %>%
      dplyr::summarise(Nx_Sj = sum(Nxa_Sj), .groups = "drop")
    b_S <- b_Sja %>%
      dplyr::group_by(dplyr::across(dplyr::all_of(paste0("x", S)))) %>%
      dplyr::summarise(Nx_Sj = sum(Nxa_Sj), .groups = "drop")

    subx <- b_S[, -ncol(b_S)]

    # Compute ν_j
    for (t in which(b_S$Nx_Sj > 0)) {
      PIs <- PI(S = S, groupTab = b_Sj, x_S = subx[t, ],
                lenX = lenX, d = d)
      dTVs <- dTV_sample(S = S, j = j, lenA = lenA, base = b_Sja,
                         A_pairs = A_pairs, x_S = subx[t, ])
      PI_xS <- as.numeric(b_S[t, ncol(b_S)]/(lenX - d))
      results$nuj2[z] <- results$nuj2[z] + prod(PIs) * dTVs/PI_xS
    }
  }
  results$selected_lags[2] <- Sc[which.max(results$nuj2)]
  cat(sprintf("Selected: j = %d (ν = %.4f)\n",
              results$selected_lags[2], max(results$nuj2)))

  # Step 3: Third selection (S = {1,364})
  cat(sprintf("\n[Step 3] Selecting third lag (S = {%d,%d})...\n",
              results$selected_lags[1], results$selected_lags[2]))
  S <- c(S, results$selected_lags[2])
  Sc <- sort(setdiff(seq_len(d), S), decreasing = TRUE)

  for (z in seq_along(Sc)) {
    j <- Sc[z]
    Sj <- sort(c(S, j), decreasing = TRUE)
    dec_S <- rev(S) # S in decreasing order

    # Frequency tables
    b_Sja <- freqTab(S = dec_S, j = j, A = A, countsTab = ct)
    b_Sj <- b_Sja %>%
      dplyr::group_by(dplyr::across(dplyr::all_of(paste0("x", Sj)))) %>%
      dplyr::summarise(Nx_Sj = sum(Nxa_Sj), .groups = "drop")
    b_S <- b_Sja %>%
      dplyr::group_by(dplyr::across(dplyr::all_of(paste0("x", dec_S)))) %>%
      dplyr::summarise(Nx_Sj = sum(Nxa_Sj), .groups = "drop")

    subx <- b_S[, -ncol(b_S)]

    # Compute ν_j
    for (t in which(b_S$Nx_Sj > 0)) {
      PIs <- PI(S = dec_S, groupTab = b_Sj, x_S = subx[t, ],
                lenX = lenX, d = d)
      dTVs <- dTV_sample(S = dec_S, j = j, lenA = lenA, base = b_Sja,
                         A_pairs = A_pairs, x_S = subx[t, ])
      PI_xS <- as.numeric(b_S[t, ncol(b_S)]/(lenX - d))
      results$nuj3[z] <- results$nuj3[z] + prod(PIs) * dTVs/PI_xS
    }
  }
  results$selected_lags[3] <- Sc[which.max(results$nuj3)]
  cat(sprintf("Selected: j = %d (ν = %.4f)\n",
              results$selected_lags[3], max(results$nuj3)))

  # Final Results
  cat("\n=== Final Selection Results ===\n")
  print(data.frame(Step = 1:3, Selected_Lag = results$selected_lags,
                   nu = c(max(results$nuj1), max(results$nuj2), max(results$nuj3))))

  return(results)
}
```


```{r sequential-selection-call, eval=FALSE}
results <- run_sequential_lag_selection(Temp12_Train)
```

```{r sequential-selection-load, echo=FALSE}
results <- results_sequential
cat("\n=== Final Selection Results ===\n")
print(data.frame(
  Step = 1:3,
  Selected_Lag = results$selected_lags,
  nu = c(max(results$nuj1), max(results$nuj2), max(results$nuj3))
))
```


### Plot Figure 4: FS sequential step analysis through $\hat{\nu}_{n,j,S}$.

```{r nu-plot, fig.width=11, fig.height=6, message=FALSE, warning=FALSE}
 par(mfrow = c(1, 3), mar = c(5, 6, 4, 2), oma = c(0, 0, 4, 0))
 palette <- c("#E41A1C", "#377EB8", "#4DAF4A")
 with(results, {
   # Graph 1
   Sc <-  364:1
   plot(1:364, rev(nuj1), type = "p", pch = 19, cex = 0.8, col = "gray70",
        ylab = "", xlab = "Lag (-j)", cex.lab = 1.8, cex.axis = 1.3,
        ylim = c(0,0.13), main = "", panel.first = grid())
   title(main = expression(paste("S = ", Ø)), cex.main = 1.5, font.main = 1)
   title(ylab = expression(widehat(nu)[n*","*j*","*S]/2), line = 3, cex.lab = 1.5)
   points(Sc[which.max(nuj1)], nuj1[which.max(nuj1)], pch = 21, bg = palette[1], cex = 1.5, lwd = 1)
   text(Sc[which.max(nuj1)], nuj1[which.max(nuj1)], labels = paste0(Sc[which.max(nuj1)]),
        pos = 3, col = palette[1], font = 2, cex = 1.4)
   # Graph 2
   Sc <- 364:2
   plot(2:364, rev(nuj2), type = "p", pch = 19, cex = 0.8, col = "gray70",
        ylab = "", xlab = "Lag (-j)", cex.lab = 1.8, cex.axis = 1.3,
        ylim = c(0,0.025),  main = "", panel.first = grid())
   title(main = "With S = {-1}", cex.main = 1.5, font.main = 1)
   title(ylab = expression(widehat(nu)[n*","*j*","*S]/2), line = 3, cex.lab = 1.5)
   points(Sc[which.max(nuj2)], nuj2[which.max(nuj2)], pch = 21, bg = palette[2], cex = 1.5, lwd = 1)
   text(Sc[which.max(nuj2)]-5, nuj2[which.max(nuj2)], labels = paste0(Sc[which.max(nuj2)]),
        pos = 3, col = palette[2], font = 2, cex = 1.4)
   # Graph 3
   Sc <- 363:2
   plot(2:363, rev(nuj3), type = "p", pch = 19, cex = 0.8, col = "gray70",
        ylab = "", xlab = "Lag (-j)", cex.lab = 1.8, cex.axis = 1.3,
        ylim = c(0,0.025), main = "", panel.first = grid())
   title(main = "With S = {-364, -1}", cex.main = 1.5, font.main = 1)
   title(ylab = expression(widehat(nu)[n*","*j*","*S]/2), line = 3, cex.lab = 1.5)
   points(Sc[which.max(nuj3)], nuj3[which.max(nuj3)], pch = 21, bg = palette[3], cex = 1.5, lwd = 1)
   text(Sc[which.max(nuj3)], nuj3[which.max(nuj3)], labels = paste0(Sc[which.max(nuj3)]),
        pos = 3, col = palette[3], font = 2, cex = 1.4)

   mtext(expression(paste("Sequential lag selection based on ", widehat(nu)[n*","*j*","*S])),
         outer = TRUE, cex = 1.6, font = 2, line = 1.4)
 })
```



## Session info
```{r session-info}
sessionInfo()
```