Estimate composite stochastic processes

In gmwmx2: Estimate Functional and Stochastic Parameters of Linear Models with Correlated Residuals and Missing Data

knitr::opts_chunk$set(
  fig.width = 10,
  fig.height = 6,
  fig.align = "center"
)


knitr::opts_chunk$set(eval = FALSE)

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

gmwmx2 documentation built on June 10, 2026, 5:06 p.m.