test_that("AR1 time-varying works", {
  skip_on_cran()
  skip_on_ci()
  set.seed(1)
  predictor_dat <- data.frame(
    X = runif(4000), Y = runif(4000),
    year = rep(seq_len(800), each = 5)
  )
  mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.2)

  sigma_V_hats <- numeric(10L)
  rho_hats <- numeric(10L)

  set.seed(1234)
  for (j in seq_along(rho_hats)) {
    print(j)
    # https://kaskr.github.io/adcomp/classdensity_1_1AR1__t.html
    ar1_scale <- 0.8
    rho <- 0.7
    sigma <- sqrt(1 - rho^2)
    devs <- rnorm(length(unique(predictor_dat$year))) * ar1_scale
    x0 <- rnorm(1) * ar1_scale
    B <- numeric(length(unique(predictor_dat$year)))
    B[1] <- rho * x0 + sigma * devs[1]
    for (i in seq(2, length(B))) {
      B[i] <- rho * B[i-1] + sigma * devs[i]
    }

    sim_dat <- sdmTMB_simulate(
      formula = ~ 0 + as.factor(year),
      data = predictor_dat,
      time = NULL,
      mesh = mesh,
      family = gaussian(),
      range = 0.5,
      sigma_E = NULL,
      phi = 0.1,
      sigma_O = 0,
      seed = j,
      B = B
    )

    sim_dat$year <- predictor_dat$year
    m <- sdmTMB(observed ~ 1, data = sim_dat, mesh = mesh,
      time = "year", spatial = 'off', time_varying_type = "ar1",
      spatiotemporal = 'off', time_varying = ~ 1)

    s <- as.list(m$sd_report, "Estimate")

    sigma_V_hats[j] <- exp(s$ln_tau_V)[1,1]
    m121 <- function(x) 2 * plogis(x) - 1
    rho_hats[j] <- m121(s$rho_time_unscaled)[1,1]
  }

  plot(s$b_rw_t[,,1], B)
  abline(0, 1)
  expect_gt(cor(B, s$b_rw_t[,,1]), 0.99)

  expect_equal(round(mean(sigma_V_hats), 2L), ar1_scale)
  expect_equal(round(mean(rho_hats), 2L), rho)

  ss <- simulate(m)
  sim_dat$obs2 <- ss[,1]
  m <- sdmTMB(obs2 ~ 0, data = sim_dat, mesh = mesh,
    time = "year", spatial = 'off', time_varying_type = "ar1",
    spatiotemporal = 'off', time_varying = ~ 1)

  s <- as.list(m$sd_report, "Estimate")
  expect_equal(dim(ss), c(4000L, 1L))
  expect_gt(exp(s$ln_tau_V)[1,1], 0.75)
  expect_gt(m121(s$rho_time_unscaled)[1,1], 0.65)
})


test_that("RW with mean-zero (rw0) time-varying works", {
  skip_on_cran()
  skip_on_ci()
  set.seed(1)
  predictor_dat <- data.frame(
    X = runif(4000), Y = runif(4000),
    year = rep(seq_len(800), each = 5)
  )
  mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.2)

  sigma_V_hats <- numeric(12L)
  sigma_true <- 0.3

  set.seed(1234)
  for (j in seq_along(sigma_V_hats)) {
    print(j)
    devs <- rnorm(length(unique(predictor_dat$year)))
    B <- numeric(length(unique(predictor_dat$year)))
    B[1] <- sigma_true * devs[1]
    for (i in seq(2, length(B))) {
      B[i] <- B[i-1] + sigma_true * devs[i]
    }

    sim_dat <- sdmTMB_simulate(
      formula = ~ 0 + as.factor(year),
      data = predictor_dat,
      time = NULL,
      mesh = mesh,
      family = gaussian(),
      range = 0.5,
      sigma_E = NULL,
      phi = 0.1,
      sigma_O = 0,
      seed = j,
      B = B
    )

    sim_dat$year <- predictor_dat$year
    m <- sdmTMB(observed ~ 1, data = sim_dat, mesh = mesh,
      time = "year", spatial = 'off', time_varying_type = "rw0",
      spatiotemporal = 'off', time_varying = ~ 1)

    s <- as.list(m$sd_report, "Estimate")

    sigma_V_hats[j] <- exp(s$ln_tau_V)[1,1]
  }

  plot(s$b_rw_t[,,1], B)
  abline(0, 1)
  expect_gt(cor(B, s$b_rw_t[,,1]), 0.99)
  expect_equal(round(mean(sigma_V_hats), 2L), sigma_true)
})