Welcome to ClientVPS Mirrors

Estimate composite stochastic processes

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

Need a high-speed mirror for your open-source project?
Contact our mirror admin team at info@clientvps.com.

This archive is provided as a free public service to the community.
Proudly supported by infrastructure from VPSPulse , RxServers , BuyNumber , UnitVPS , OffshoreName and secure payment technology by ArionPay.