From c8d0895e7f3465d50d236b2064b36d08c10e5404 Mon Sep 17 00:00:00 2001
From: Lionel Voirol <lionelvoirol@hotmail.com>
Date: Thu, 9 Nov 2023 20:54:11 +0100
Subject: [PATCH] add scripts of demo generating signal

---
 scripts/generate_signal.R | 93 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 93 insertions(+)
 create mode 100644 scripts/generate_signal.R

diff --git a/scripts/generate_signal.R b/scripts/generate_signal.R
new file mode 100644
index 0000000..d605c31
--- /dev/null
+++ b/scripts/generate_signal.R
@@ -0,0 +1,93 @@
+# clean ws
+rm(list=ls())
+
+# Load pkg
+library(simts)
+library(gmwmx)
+library(progress)
+
+# we consider example the model considered in model 3
+phase <- 0.45
+amplitude <- 2.5
+sigma2_wn <- 10
+gamma2_rw = .5
+bias <- 0
+trend <- 5 / 365.25
+cosU <- amplitude * cos(phase)
+sinU <- amplitude * sin(phase)
+
+# consider n years of daily observations
+year <- 5
+n <- year * 365
+
+
+# define model for generating gaussian white noise + RW
+model_gaussian_wn_rw <- WN(sigma2 = sigma2_wn) + RW(gamma2 = gamma2_rw)
+
+# generate data
+# define time at which there are jumps
+jump_vec <- c(200, 300, 500)
+jump_height <- c(10, 15, 20)
+
+# define seed
+myseed <- 123
+
+# add trend, gaps and sin
+nbr_sin <- 1 # define number of sinusoidal process
+
+# define A matrix
+A <- create_A_matrix(1:n, jump_vec, n_seasonal = 2)
+
+# define beta
+x_0 <- c(bias, trend, cosU, sinU, jump_height)
+
+# define number of Monte Carlo simulation
+n_simu <- 100
+
+# mat_res = matrix(ncol=18, nrow=n_simu)
+
+mat_res = matrix(ncol=9, nrow=n_simu)
+
+pb <- progress_bar$new(total = n_simu)
+
+for(b in seq(n_simu)){
+  
+  # b = 4
+  # fix seed for reproducibility
+  set.seed(myseed + b)
+  
+  # generate residuals from a Gaussian White noise + Power law process
+  eps <- simts::gen_gts(model = model_gaussian_wn_rw, n = n)
+  # plot(wv::wvar(eps))
+  
+  # create time series
+  yy <- A %*% x_0 + eps
+  
+  # create gnssts
+  gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
+  
+  # gmwmx 1 step
+  fit_gmwmx <- estimate_gmwmx(
+    x = gnssts_obj,
+    model = "wn+rw",
+    theta_0 = c(1, 1),
+    n_seasonal = 1,maxit = 5000,
+    k_iter = 1
+  )
+  
+
+  
+  mat_res[b, ] = c(  fit_gmwmx$beta_hat,
+                     fit_gmwmx$theta_hat)
+  
+  pb$tick()
+  
+}
+
+
+
+# plot empirical distribution of estimated stochastic parameters
+boxplot(mat_res[, c(8)])
+abline(h = sigma2_wn)
+boxplot(mat_res[, c(9)])
+abline(h =gamma2_rw)