inst/doc/estimate_linear_models_with_dependent_errors.R

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

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  fig.width = 10,
  fig.height = 6,
  fig.align = "center"
)

# so as not to execute this vignette during package development
knitr::opts_chunk$set(eval = FALSE)

## ----message=FALSE,warning=FALSE----------------------------------------------
# library(gmwmx2)
# library(wv)
# library(dplyr)
# 
# 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")
# }

## -----------------------------------------------------------------------------
# n = 15*365
# X = matrix(NA, nrow = n, ncol = 4)
# # intercept
# X[, 1] = 1
# # trend
# X[, 2] = 1:n
# # annual sinusoid
# omega_1 <- (1 / 365.25) * 2 * pi
# X[, 3] <- sin((1:n) * omega_1)
# X[, 4] <- cos((1:n) * omega_1)
# beta = c(0.5 , 0.00004, 0.005,  0.0008)
# 
# 
# # visualize the deterministic signal
# plot(x = X[, 2], y = X %*% beta, type = "l",
#      main = "Signal", xlab = "Time", ylab = "Signal", las = 1)

## -----------------------------------------------------------------------------
# phi_ar1 = 0.975
# sigma2_ar1 =  5e-05
# sigma2_wn =  7e-04
# eps = gmwmx2::generate(ar1(phi = phi_ar1, sigma2 = sigma2_ar1) + wn(sigma2_wn), n = n, seed = 123)$series
# plot(wv::wvar(eps))

## -----------------------------------------------------------------------------
# y = X %*% beta + eps
# plot(X[, 2], y, type = "l", main = "Simulated y (WN + AR1)")

## -----------------------------------------------------------------------------
# fit = gmwmx2(X = X, y = y, model = wn() + ar1() )
# print(fit)
# 

## -----------------------------------------------------------------------------
# B = 100
# mat_res = matrix(NA, nrow=B, ncol=19)
# for(b in seq(B)){
#   eps = generate(ar1(phi=phi_ar1, sigma2=sigma2_ar1) + wn(sigma2_wn), n=n, seed = (123 + b))$series
#   y = X %*% beta + eps
#   fit = gmwmx2(X = X, y = y, model = wn() + ar1() )
#   # misspecified model assuming white noise as the stochastic model
#   fit2 = lm(y~X[,2] + X[,3] + X[,4])
# 
#   mat_res[b, ] = c(fit$beta_hat, fit$std_beta_hat,
#                    summary(fit2)$coefficients[,1],
#                    summary(fit2)$coefficients[,2],
#                    fit$theta_domain$`AR(1)_2`,
#                    fit$theta_domain$`White Noise_1`)
#   # cat("Iteration ", b, " completed.\n")
# }

## -----------------------------------------------------------------------------
# # compute empirical coverage
# mat_res_df = as.data.frame(mat_res)
# colnames(mat_res_df) = c("gmwmx_beta0_hat", "gmwmx_beta1_hat",
#                          "gmwmx_beta2_hat", "gmwmx_beta3_hat",
#                          "gmwmx_std_beta0_hat", "gmwmx_std_beta1_hat",
#                          "gmwmx_std_beta2_hat", "gmwmx_std_beta3_hat",
#                          "lm_beta0_hat", "lm_beta1_hat", "lm_beta2_hat", "lm_beta3_hat",
#                          "lm_std_beta0_hat", "lm_std_beta1_hat", "lm_std_beta2_hat", "lm_std_beta3_hat",
#                          "phi_ar1","sigma_2_ar1" ,"sigma_2_wn")
# 

## -----------------------------------------------------------------------------
# par(mfrow = c(1, 4))
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta0_hat")],las=1,
#         names = c("beta0"),
#         main = expression(beta[0]), ylab = "Estimate")
# abline(h = beta[1], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta1_hat")],las=1,
#         names = c("beta1"),
#         main = expression(beta[1]), ylab = "Estimate")
# abline(h = beta[2], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta2_hat")],las=1,
#         names = c("beta2"),
#         main = expression(beta[2]), ylab = "Estimate")
# abline(h = beta[3], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta3_hat")],las=1,
#         names = c("beta3"),
#         main = expression(beta[3]), ylab = "Estimate")
# abline(h = beta[4], col = "black", lwd = 2)
# 
# 
# par(mfrow = c(1, 3))
# 
# boxplot_mean_dot(mat_res_df$phi_ar1,las=1,
#         names = c("phi_ar1"),
#         main = expression(phi["AR1"]), ylab = "Estimate")
# abline(h = phi_ar1, col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df$sigma_2_ar1,las=1,
#         names = c("phi_ar1"),
#         main = expression(sigma["AR1"]^2), ylab = "Estimate")
# abline(h = sigma2_ar1, col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df$sigma_2_wn,las=1,
#         names = c("phi_ar1"),
#         main = expression(sigma["WN"]^2), ylab = "Estimate")
# abline(h = sigma2_wn, col = "black", lwd = 2)
# par(mfrow = c(1, 1))

## -----------------------------------------------------------------------------
# zval = qnorm(0.975)
# mat_res_df$upper_ci_gmwmx_beta1 = mat_res_df$gmwmx_beta1_hat + zval * mat_res_df$gmwmx_std_beta1_hat
# mat_res_df$lower_ci_gmwmx_beta1 = mat_res_df$gmwmx_beta1_hat - zval * mat_res_df$gmwmx_std_beta1_hat
# # empirical coverage of gmwmx beta
# dplyr::between(rep(beta[2], B), mat_res_df$lower_ci_gmwmx_beta1, mat_res_df$upper_ci_gmwmx_beta1) %>% mean()

## -----------------------------------------------------------------------------
# # do the same for lm beta
# mat_res_df$upper_ci_lm_beta1 = mat_res_df$lm_beta1_hat + zval * mat_res_df$lm_std_beta1_hat
# mat_res_df$lower_ci_lm_beta1 = mat_res_df$lm_beta1_hat - zval * mat_res_df$lm_std_beta1_hat
# dplyr::between(rep(beta[2], B), mat_res_df$lower_ci_lm_beta1, mat_res_df$upper_ci_lm_beta1) %>% mean()

## -----------------------------------------------------------------------------
# kappa_pl <- -0.75
# sigma2_wn <- 2e-07
# sigma2_pl <- 2.25e-06
# eps_pl = generate(wn(sigma2_wn) + pl(kappa = kappa_pl, sigma2 = sigma2_pl),
#                   n = n, seed = 123)
# plot(eps_pl$series, type = "l")
# plot(wv::wvar(eps_pl$series))
# y_pl = X %*% beta + eps_pl$series
# plot(X[, 2], y_pl, type = "l", main = "Simulated y (WN + PL)")

## -----------------------------------------------------------------------------
# fit_pl = gmwmx2(X = X, y = y_pl, model = wn() + pl())
# fit_pl

## -----------------------------------------------------------------------------
# B_pl = 100
# mat_res_pl = matrix(NA, nrow = B_pl, ncol = 19)
# for(b in seq(B_pl)){
#   eps = generate(wn(sigma2_wn) + pl(kappa = kappa_pl, sigma2 = sigma2_pl),
#                  n = n, seed = (123 + b))$series
#   y = X %*% beta + eps
#   fit = gmwmx2(X = X, y = y, model = wn() + pl())
#   fit2 = lm(y~X[,2] + X[,3] + X[,4])
# 
#   mat_res_pl[b, ] = c(fit$beta_hat, fit$std_beta_hat,
#                       summary(fit2)$coefficients[,1],
#                       summary(fit2)$coefficients[,2],
#                       fit$theta_domain$`Stationary PowerLaw_2`,
#                       fit$theta_domain$`White Noise_1`)
#   # cat("Iteration ", b, " \n")
# }
# 
# mat_res_pl_df = as.data.frame(mat_res_pl)
# colnames(mat_res_pl_df) = c("gmwmx_beta0_hat", "gmwmx_beta1_hat",
#                             "gmwmx_beta2_hat", "gmwmx_beta3_hat",
#                             "gmwmx_std_beta0_hat", "gmwmx_std_beta1_hat",
#                             "gmwmx_std_beta2_hat", "gmwmx_std_beta3_hat",
#                             "lm_beta0_hat", "lm_beta1_hat", "lm_beta2_hat", "lm_beta3_hat",
#                             "lm_std_beta0_hat", "lm_std_beta1_hat",
#                             "lm_std_beta2_hat", "lm_std_beta3_hat",
#                             "kappa_pl","sigma2_pl" ,"sigma2_wn")

## -----------------------------------------------------------------------------
# par(mfrow = c(1, 4))
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta0_hat")],las=1,
#         names = c("beta0"),
#         main = expression(beta[0]), ylab = "Estimate")
# abline(h = beta[1], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta1_hat")],las=1,
#         names = c("beta1"),
#         main = expression(beta[1]), ylab = "Estimate")
# abline(h = beta[2], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta2_hat")],las=1,
#         names = c("beta2"),
#         main = expression(beta[2]), ylab = "Estimate")
# abline(h = beta[3], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta3_hat")],las=1,
#         names = c("beta3"),
#         main = expression(beta[3]), ylab = "Estimate")
# abline(h = beta[4], col = "black", lwd = 2)
# 
# par(mfrow = c(1, 1))

## -----------------------------------------------------------------------------
# par(mfrow = c(1, 3))
# 
# boxplot_mean_dot(mat_res_pl_df$kappa_pl,las=1,
#         main = expression(kappa["PL"]), ylab = "Estimate")
# abline(h = kappa_pl, col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_pl_df$sigma2_pl,las=1,
#         main = expression(sigma["PL"]^2), ylab = "Estimate")
# abline(h = sigma2_pl, col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_pl_df$sigma2_wn,las=1,
#         names = c("phi_ar1"),
#         main = expression(sigma["WN"]^2), ylab = "Estimate")
# abline(h = sigma2_wn, col = "black", lwd = 2)
# par(mfrow = c(1, 1))

## -----------------------------------------------------------------------------
# zval = qnorm(0.975)
# mat_res_pl_df$upper_ci_gmwmx_beta1 = mat_res_pl_df$gmwmx_beta1_hat + zval * mat_res_pl_df$gmwmx_std_beta1_hat
# mat_res_pl_df$lower_ci_gmwmx_beta1 = mat_res_pl_df$gmwmx_beta1_hat - zval * mat_res_pl_df$gmwmx_std_beta1_hat
# # empirical coverage of gmwmx beta
# dplyr::between(rep(beta[2], B_pl), mat_res_pl_df$lower_ci_gmwmx_beta1, mat_res_pl_df$upper_ci_gmwmx_beta1) %>% mean()
# 
# # do the same for lm beta
# mat_res_pl_df$upper_ci_lm_beta1 = mat_res_pl_df$lm_beta1_hat + zval * mat_res_pl_df$lm_std_beta1_hat
# mat_res_pl_df$lower_ci_lm_beta1 = mat_res_pl_df$lm_beta1_hat - zval * mat_res_pl_df$lm_std_beta1_hat
# dplyr::between(rep(beta[2], B_pl), mat_res_pl_df$lower_ci_lm_beta1, mat_res_pl_df$upper_ci_lm_beta1) %>% mean()

## -----------------------------------------------------------------------------
# # fix stochastic parameters
# sigma2_wn_fl <- 8e-07
# sigma2_fl <- 1e-06
# eps_fl = generate(wn(sigma2_wn_fl) + flicker(sigma2 = sigma2_fl),
#                   n = n, seed = 123)
# plot(eps_fl$series, type = "l")
# plot(wv::wvar(eps_fl$series))
# y_fl = X %*% beta + eps_fl$series
# plot(X[, 2], y_fl, type = "l", main = "Simulated y (WN + Flicker)")

## -----------------------------------------------------------------------------
# fit_fl = gmwmx2(X = X, y = y_fl, model = wn() + flicker())
# fit_fl

## -----------------------------------------------------------------------------
# B_fl = 100
# mat_res_fl = matrix(NA, nrow = B_fl, ncol = 18)
# for(b in seq(B_fl)){
#   eps = generate(wn(sigma2_wn_fl) + flicker(sigma2 = sigma2_fl),
#                  n = n, seed = (123 + b))$series
#   y = X %*% beta + eps
#   fit = gmwmx2(X = X, y = y, model = wn() + flicker())
#   fit2 = lm(y~X[,2] + X[,3] + X[,4])
# 
#   mat_res_fl[b, ] = c(fit$beta_hat, fit$std_beta_hat,
#                       summary(fit2)$coefficients[,1],
#                       summary(fit2)$coefficients[,2],
#                       fit$theta_domain$`Flicker_2`,
#                       fit$theta_domain$`White Noise_1`)
#   # cat("Iteration ", b, " \n")
# }
# 
# mat_res_fl_df = as.data.frame(mat_res_fl)
# colnames(mat_res_fl_df) = c("gmwmx_beta0_hat", "gmwmx_beta1_hat",
#                             "gmwmx_beta2_hat", "gmwmx_beta3_hat",
#                             "gmwmx_std_beta0_hat", "gmwmx_std_beta1_hat",
#                             "gmwmx_std_beta2_hat", "gmwmx_std_beta3_hat",
#                             "lm_beta0_hat", "lm_beta1_hat", "lm_beta2_hat", "lm_beta3_hat",
#                             "lm_std_beta0_hat", "lm_std_beta1_hat",
#                             "lm_std_beta2_hat", "lm_std_beta3_hat",
#                             "sigma2_fl" ,"sigma2_wn")

## -----------------------------------------------------------------------------
# par(mfrow = c(1, 4))
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta0_hat")],las=1,
#         names = c("beta0"),
#         main = expression(beta[0]), ylab = "Estimate")
# abline(h = beta[1], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta1_hat")],las=1,
#         names = c("beta1"),
#         main = expression(beta[1]), ylab = "Estimate")
# abline(h = beta[2], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta2_hat")],las=1,
#         names = c("beta2"),
#         main = expression(beta[2]), ylab = "Estimate")
# abline(h = beta[3], col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_df[, c("gmwmx_beta3_hat")],las=1,
#         names = c("beta3"),
#         main = expression(beta[3]), ylab = "Estimate")
# abline(h = beta[4], col = "black", lwd = 2)
# 
# par(mfrow = c(1, 1))

## -----------------------------------------------------------------------------
# par(mfrow = c(1, 2))
# 
# boxplot_mean_dot(mat_res_fl_df$sigma2_fl,las=1,
#         main = expression(sigma["FL"]^2), ylab = "Estimate")
# abline(h = sigma2_fl, col = "black", lwd = 2)
# 
# boxplot_mean_dot(mat_res_fl_df$sigma2_wn,las=1,
#         names = c("phi_ar1"),
#         main = expression(sigma["WN"]^2), ylab = "Estimate")
# abline(h = sigma2_wn_fl, col = "black", lwd = 2)
# par(mfrow = c(1, 1))

## -----------------------------------------------------------------------------
# zval = qnorm(0.975)
# mat_res_fl_df$upper_ci_gmwmx_beta1 = mat_res_fl_df$gmwmx_beta1_hat + zval * mat_res_fl_df$gmwmx_std_beta1_hat
# mat_res_fl_df$lower_ci_gmwmx_beta1 = mat_res_fl_df$gmwmx_beta1_hat - zval * mat_res_fl_df$gmwmx_std_beta1_hat
# # empirical coverage of gmwmx beta
# dplyr::between(rep(beta[2], B_pl), mat_res_fl_df$lower_ci_gmwmx_beta1, mat_res_fl_df$upper_ci_gmwmx_beta1) %>% mean()
# 
# # do the same for lm beta
# mat_res_fl_df$upper_ci_lm_beta1 = mat_res_fl_df$lm_beta1_hat + zval * mat_res_fl_df$lm_std_beta1_hat
# mat_res_fl_df$lower_ci_lm_beta1 = mat_res_fl_df$lm_beta1_hat - zval * mat_res_fl_df$lm_std_beta1_hat
# dplyr::between(rep(beta[2], B_pl), mat_res_fl_df$lower_ci_lm_beta1, mat_res_fl_df$upper_ci_lm_beta1) %>% mean()