# Simulation study

This vignette reproduce a design considered in Cucci et al. (2022) and implement a Monte Carlo simulation comparing the performance of the GMWMX-1, the GMWMX-2 and the MLE implemented in Hector. Note that this file is not intended to be run on a personal computer but should rather be executed on a high performance computing cluster.

    # libraries
    library(simts)
    library(gmwmx)
    library(dplyr)

    check_if_theta_in_ci <- function(theta_vec, ci_mat) {
      vec_emp_coverage <- vector(mode = "logical", length = length(theta_vec))
      for (i in seq(length(theta_vec))) {
        if (dplyr::between(theta_vec[i], ci_mat[i, 1], ci_mat[i, 2])) {
          vec_emp_coverage[i] <- T
        } else {
          vec_emp_coverage[i] <- F
        }
      }
      return(as.numeric(vec_emp_coverage))
    }

    # we consider example the model considered in model 3
    phase <- 0.45
    amplitude <- 2.5
    sigma2_wn <- 15
    sigma2_powerlaw <- 10
    d <- 0.4
    bias <- 0
    trend <- 5 / 365.25
    cosU <- amplitude * cos(phase)
    sinU <- amplitude * sin(phase)

Let us consider 5 years of daily data

    # consider n years of daily observations
    year <- 20
    n <- year * 365

We consider a gaussian White noise and a Power Law process as the stochastic model of the residuals \(\varepsilon\).

Note that the functions that enable to generate stochastic models that include Power Law process, Matèrn process or Fractional Gaussian noise are (for now) only available from the development version of the package `simts` that can be easily installed with:

    install.packages("devtools")
    devtools::install_github("SMAC-Group/simts")

    # define model for generating gaussian white noise + PLP
    model_gaussian_wn_plp <- WN(sigma2 = sigma2_wn) + PLP(sigma2 = sigma2_powerlaw, d = d)

We consider three location shifts (jumps) at three different heights.

    # generate data
    # define time at which there are jumps
    jump_vec <- c(200, 300, 500)
    jump_height <- c(10, 15, 20)

We define a seed for reproducibility

    # define seed
    myseed <- 123

    # add trend, gaps and sin
    nbr_sin <- 1 # define number of sinusoidal process

and create the matrix \(\boldsymbol{A}\) with function `create_A_matrix()`.

    # define A matrix
    A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)

We set the vector of functional parameters \(\boldsymbol{x}_0\).

    # define beta
    x_0 <- c(bias, trend, cosU, sinU, jump_height)

We define the number of simulations to perform

    # define number of Monte Carlo simulation
    n_simu <- 1000

    # define number of parameter estimated (depends on the model)
    # bias + trend + height * nbr of jump vec + 2* nbr of sin process + wn + powerlaw parameters
    nbr_param_check_coverage <- 4
    nbr_param <- 2 + length(jump_vec) + 2 * nbr_sin + 3
    dim_mat_results <- (nbr_param + nbr_param_check_coverage) * 3 + 1 * 3

    # define matrix of results
    mat_results <- matrix(NA, ncol = dim_mat_results, nrow = n_simu)

We perform the simulation, saving for each Monte Carlo run, the estimated parameters of the functional and stochastic model and the empirical coverage on the functional parameters for the GMWMX-1, the GMWMX-2 as well as the MLE implemented in Hector.

    for (simu_b in seq(n_simu)) {

      # fix seed for reproducibility
      set.seed(myseed + simu_b)

      # generate residuals from a Gaussian White noise + Power law process
      eps <- simts::gen_gts(model = model_gaussian_wn_plp, n = n)

      # create time series
      yy <- A %*% x_0 + eps

      # create gnssts
      gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)

      # MLE
      fit_simu_b_mle <- estimate_hector(
        x = gnssts_obj,
        model = "wn+powerlaw",
        n_seasonal = 1
      )

      # gmwmx 1 step
      fit_simu_b_gmwm_1_step <- estimate_gmwmx(
        x = gnssts_obj,
        model = "wn+powerlaw",
        theta_0 = c(0.1, 0.1, 0.1),
        n_seasonal = 1,
        ci = T,
        k_iter = 1
      )

      # gmwmx 2 steps
      fit_simu_b_gmwm_2_step <- estimate_gmwmx(
        x = gnssts_obj,
        model = "wn+powerlaw",
        theta_0 = c(0.1, 0.1, 0.1),
        n_seasonal = 1,
        ci = T,
        k_iter = 2
      )

      # define vector of estimators
      fit_mle <- c(fit_simu_b_mle$beta_hat, fit_simu_b_mle$theta_hat)
      fit_gmwm_1_step <- c(fit_simu_b_gmwm_1_step$beta_hat, fit_simu_b_gmwm_1_step$theta_hat)
      fit_gmwm_2_step <- c(fit_simu_b_gmwm_2_step$beta_hat, fit_simu_b_gmwm_2_step$theta_hat)

      # compute coverage for each method
      alpha <- .05
      z_val <- qnorm(1 - alpha / 2)

      mat_ci_mle_beta <- matrix(c(
        fit_simu_b_mle$beta_hat - z_val * fit_simu_b_mle$beta_std,
        fit_simu_b_mle$beta_hat + z_val * fit_simu_b_mle$beta_std
      ),
      byrow = F, ncol = 2
      )

      mat_ci_gmwm_1_step_beta <- matrix(c(
        fit_simu_b_gmwm_1_step$beta_hat - z_val * fit_simu_b_gmwm_1_step$beta_std,
        fit_simu_b_gmwm_1_step$beta_hat + z_val * fit_simu_b_gmwm_1_step$beta_std
      ),
      byrow = F, ncol = 2
      )

      mat_ci_gmwm_2_step_beta <- matrix(c(
        fit_simu_b_gmwm_2_step$beta_hat - z_val * fit_simu_b_gmwm_2_step$beta_std,
        fit_simu_b_gmwm_2_step$beta_hat + z_val * fit_simu_b_gmwm_2_step$beta_std
      ),
      byrow = F, ncol = 2
      )

      # save empirical coverage
      inside_ci_mle <- check_if_theta_in_ci(x_0, mat_ci_mle_beta)[1:4]
      inside_ci_gmwm_1 <- check_if_theta_in_ci(x_0, mat_ci_gmwm_1_step_beta)[1:4]
      inside_ci_gmwm_2 <- check_if_theta_in_ci(x_0, mat_ci_gmwm_2_step_beta)[1:4]

      # define vector of estimated parameters as object
      res <- c(
        fit_mle, fit_gmwm_1_step,
        fit_gmwm_2_step,
        inside_ci_mle,
        inside_ci_gmwm_1, inside_ci_gmwm_2,
        fit_simu_b_mle$estimation_time[3], fit_simu_b_gmwm_1_step$estimation_time[3], fit_simu_b_gmwm_2_step$estimation_time[3]
      )

      # save in matrix of results
      mat_results[simu_b, ] <- res

      # print status
      cat(paste("Completed simulation", simu_b, "\n",sep = " "))
    }

    # define names of mat results
    name_param_functionnal <- c("bias", "trend", "cosU", "sinU")
    colnames(mat_results) <- c(
        paste("mle", names(fit_mle), sep = "_"),
        paste("gmwm_1_step", names(fit_gmwm_1_step), sep = "_"),
        paste("gmwm_2_step", names(fit_gmwm_2_step), sep = "_"),
        paste0("mle_inside_ci_", name_param_functionnal),
        paste0("gmwm_1_inside_ci_", name_param_functionnal),
        paste0("gmwm_2_inside_ci_", name_param_functionnal),
        c("time_mle", "time_gmwm_1", "time_gmwm_2")
      )