## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 11,
  dpi = 120,
  out.width = "100%",
  warning = FALSE,
  message = FALSE
)

## ----packages-----------------------------------------------------------------
library(VARcheck)

## ----helpers------------------------------------------------------------------
fit_ar1 <- function(x) {
  n     <- length(x)
  mod   <- lm(x[-1] ~ x[-n])
  xhat  <- c(NA, predict(mod))
  res   <- x - xhat
  list(mod = mod, xhat = xhat, res = res, res_var = sd(res, na.rm = TRUE))
}

sim_ar1 <- function(fit, n, seed = 92) {
  b    <- fit$mod$coefficients
  xsim <- numeric(n)
  set.seed(seed)
  for (i in 2:n) xsim[i] <- b[1] + b[2] * xsim[i - 1] + rnorm(1, 0, fit$res_var)
  xsim
}

## ----dgms---------------------------------------------------------------------
N <- 200

# DGM 1: correctly specified, AR > 0
set.seed(2)
x1 <- numeric(N)
for (i in 2:N) x1[i] <- 0.6 * x1[i - 1] + rnorm(1)

# DGM 2: correctly specified, AR = 0 (white noise)
set.seed(3)
x2 <- numeric(N)
for (i in 2:N) x2[i] <- rnorm(1)

# DGM 3: trend
set.seed(2)
x3 <- numeric(N); x3[1] <- -15
for (i in 2:N) x3[i] <- 0.6 * x3[i - 1] + rnorm(1) - 6 + 0.06 * i
x3 <- (x3 / sd(x3)) * 2

# DGM 4: switching / non-linear
x4 <- c(rep(-2, 30), rep(2, 25), rep(-2, 40), rep(2, 10),
        rep(-2, 60), rep(2, 35)) + rnorm(N, sd = 0.5)

# DGM 5: heteroscedasticity (growing innovation variance)
set.seed(2)
x5 <- numeric(N)
for (i in 2:N) x5[i] <- 0.6 * x5[i - 1] + rnorm(1, sd = 0.1 + 0.012 * i)

# DGM 6: seasonality (weekend effect)
set.seed(4)
weekend <- rep(c(rep(0, 5), rep(1, 2)), 29)[1:N]
x6 <- numeric(N)
for (i in 2:N) x6[i] <- -1 + 2 * weekend[i] + 0.6 * x6[i - 1] + rnorm(1)

# DGM 7: state-dependent innovations
set.seed(2)
x7 <- numeric(N)
for (i in 2:N) x7[i] <- 0.6 * x7[i - 1] + rnorm(1, sd = 1 + 0.3 * x7[i - 1])

# DGM 8: non-Gaussian innovations (exponential)
set.seed(2)
x8 <- numeric(N)
for (i in 2:N) x8[i] <- -1 + 0.6 * x8[i - 1] + rexp(1)

# Fit AR(1) to each DGM
fits <- lapply(list(x1, x2, x3, x4, x5, x6, x7, x8), fit_ar1)
sims <- lapply(fits, sim_ar1, n = N)

## ----correct, fig.height = 7--------------------------------------------------
vd_correct <- new_var_data(
  empirical = cbind(x1, x2),
  predicted = cbind(fits[[1]]$xhat, fits[[2]]$xhat),
  residuals = cbind(fits[[1]]$res,  fits[[2]]$res),
  simulated = cbind(sims[[1]],      sims[[2]]),
  var_names = c("Dependence", "Independence")
)

plot_var_check(vd_correct)

## ----misspec-main, fig.height = 9---------------------------------------------
vd_main <- new_var_data(
  empirical = cbind(x4, x3, x5),
  predicted = cbind(fits[[4]]$xhat, fits[[3]]$xhat, fits[[5]]$xhat),
  residuals = cbind(fits[[4]]$res,  fits[[3]]$res,  fits[[5]]$res),
  simulated = cbind(sims[[4]],      sims[[3]],      sims[[5]]),
  var_names = c("Switching", "Trend", "Heteroscedasticity")
)

plot_var_check(vd_main)

## ----misspec-extra, fig.height = 9--------------------------------------------
vd_extra <- new_var_data(
  empirical = cbind(x6, x7, x8),
  predicted = cbind(fits[[6]]$xhat, fits[[7]]$xhat, fits[[8]]$xhat),
  residuals = cbind(fits[[6]]$res,  fits[[7]]$res,  fits[[8]]$res),
  simulated = cbind(sims[[6]],      sims[[7]],      sims[[8]]),
  var_names = c("Seasonality", "State-Dependent Innovations", "Non-Gaussian Innovations")
)

plot_var_check(vd_extra)