Estimate composite stochastic processes

This vignette shows how to estimate composite stochastic models using gmwm2(). We generate data from known models, fit composite candidates, and visualize the results.

We consider a zero-mean stochastic process \(\{Y_t\}_{t=1}^n\) generated from a composite model parameterized by \(\boldsymbol{\gamma} \in \boldsymbol{\Gamma}\).

Given a parametric model with parameters \(\boldsymbol{\gamma}\), the GMWM estimator computes \(\hat{\boldsymbol{\gamma}}\) by solving: \[\begin{equation} \hat{\boldsymbol{\gamma}} = \underset{\boldsymbol{\gamma} \in \boldsymbol{\Gamma}}{\arg\min} \left\{\hat{\boldsymbol{\nu}} - \boldsymbol{\nu}(\boldsymbol{\gamma})\right\}^{\top} \boldsymbol{\Omega} \left\{\hat{\boldsymbol{\nu}} - \boldsymbol{\nu}(\boldsymbol{\gamma})\right\}, \end{equation}\]

where \(\hat{\boldsymbol{\nu}}\) is the empirical wavelet variance of the observed series and \(\boldsymbol{\nu}(\boldsymbol{\gamma})\) is the theoretical wavelet variance implied by the model.

library(gmwmx2)
boxplot_mean_dot <- function(x, ...) {
  boxplot(x, ...)
  x_mat <- as.matrix(x)
  mean_vals <- colMeans(x_mat, na.rm = TRUE)
  points(seq_along(mean_vals), mean_vals, pch = 16, col = "black", cex = 1.5)
}

Example 1 (White noise + AR(1))

n <- 10000
sigma2_wn = 25
phi_ar1 = 0.99
sigma2_ar1 = 1
mod1 <- wn(sigma2 = sigma2_wn) + ar1(phi = phi_ar1, sigma2 = sigma2_ar1)
y1 <- generate(mod1, n = n)
plot(y1)
fit1 <- gmwm2(y1, model = wn() + ar1())
fit1
plot(fit1)

Monte Carlo simulation

B = 100
mat_res = matrix(NA, nrow = B, ncol = 3)
for(b in seq(B)){
  y <- generate(mod1, n = n)
  fit = gmwm2(y, model = wn() + ar1())
  mat_res[b,] = c(fit$theta_domain$`White Noise_1`, fit$theta_domain$`AR(1)_2`)
  # cat("Done with b =", b, "\n")
}
par(mfrow = c(1, 3))
boxplot_mean_dot(mat_res[, 1], main = expression(sigma[WN]^2), ylab = "Estimated Value")
abline(h = sigma2_wn)
boxplot_mean_dot(mat_res[, 2], main = expression(phi[AR1]), ylab = "Estimated Value")
abline(h = phi_ar1)
boxplot_mean_dot(mat_res[, 3], main = expression(sigma[AR1]^2), ylab = "Estimated Value")
abline(h = sigma2_ar1)
par(mfrow = c(1, 1))

Example 2 (White noise + stationary powerlaw)

sigma2_wn = 5
kappa_pl = -0.9
sigma2_pl = 1
mod2 <- wn(sigma2_wn) + pl(kappa = kappa_pl, sigma2 = sigma2_pl)
y2 <- generate(mod2, n = n)
plot(y2)
fit2 <- gmwm2(y2, model = wn() + pl())
fit2
plot(fit2)

Monte Carlo simulation

B = 100
mat_res = matrix(NA, nrow = B, ncol = 3)
for(b in seq(B)){
  y <- generate(mod2, n = n)
  fit2 = gmwm2(y, model = wn() + pl())
  mat_res[b,] = c(fit2$theta_domain$`White Noise_1`, fit2$theta_domain$`Stationary PowerLaw_2`)
  # cat("Done with b =", b, "\n")
}
par(mfrow = c(1, 3))
boxplot_mean_dot(mat_res[, 1], main = expression(sigma[WN]^2), ylab = "Estimated Value")
abline(h = sigma2_wn)
boxplot_mean_dot(mat_res[, 2], main = expression(kappa[PL]), ylab = "Estimated Value")
abline(h = kappa_pl)
boxplot_mean_dot(mat_res[, 3], main = expression(sigma[PL]^2), ylab = "Estimated Value")
abline(h = sigma2_pl)
par(mfrow = c(1, 1))

Example 3 (White noise + AR(1) + random walk)

sigma2_wn = 5
phi_ar1 = 0.98
sigma2_ar1 = 1
sigma2_rw = 0.1
mod3 <- wn(sigma2_wn) + ar1(phi = phi_ar1, sigma2 = sigma2_ar1) + rw(sigma2_rw)
y3 <- generate(mod3, n = n)
plot(y3)
fit3 <- gmwm2(y3, model = wn() + ar1() + rw())
fit3        
plot(fit3)

Monte Carlo simulation

B = 100
mat_res = matrix(NA, nrow = B, ncol = 4)
for(b in seq(B)){
  y <- generate(mod3, n = n)
  fit3 = gmwm2(y, model = wn() + ar1() + rw())
  mat_res[b,] = c(fit3$theta_domain$`White Noise_1`, fit3$theta_domain$`AR(1)_2`, fit3$theta_domain$`Random Walk_3`)
  # cat("Done with b =", b, "\n")
}
par(mfrow = c(1, 4))
boxplot_mean_dot(mat_res[, 1], main = expression(sigma[WN]^2), ylab = "Estimated Value")
abline(h = sigma2_wn)
boxplot_mean_dot(mat_res[, 2], main = expression(phi[AR1]), ylab = "Estimated Value")
abline(h = phi_ar1)
boxplot_mean_dot(mat_res[, 3], main = expression(sigma[AR1]^2), ylab = "Estimated Value")
abline(h = sigma2_ar1)
boxplot_mean_dot(mat_res[, 4], main = expression(sigma[RW]^2), ylab = "Estimated Value")
abline(h = sigma2_rw)
par(mfrow = c(1, 1))
sigma2_wn = 20
sigma2_rw = 0.1
mod4 <- wn(sigma2_wn) + rw(sigma2_rw)
y4 <- generate(mod4, n = n)
plot(y4)
fit4 <- gmwm2(y4, model = wn() + rw())
fit4
plot(fit4)

Monte Carlo simulation

B = 100
mat_res = matrix(NA, nrow = B, ncol = 2)
for(b in seq(B)){
  y <- generate(mod4, n = n)
  fit4 = gmwm2(y, model = wn()  + rw())
  mat_res[b,] = c(fit4$theta_domain$`White Noise_1`, fit4$theta_domain$`Random Walk_2`)
  # cat("Done with b =", b, "\n")
}
par(mfrow = c(1, 2))
boxplot_mean_dot(mat_res[, 1], main = expression(sigma[WN]^2), ylab = "Estimated Value")
abline(h = sigma2_wn)
boxplot_mean_dot(mat_res[, 2], main = expression(sigma[RW]^2), ylab = "Estimated Value")
abline(h = sigma2_rw)
par(mfrow = c(1, 1))