The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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))

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.