tests/testthat/testNumerical.R

In helske/KFAS: Kalman Filter and Smoother for Exponential Family State Space Models

context("Numerical issues")

test_that("Test numerical issues",{
  tol <- 1e-6
  
  y <-  1 / exp(1:250)
  x <- c(1, y[-250])
  model <- SSModel(y ~ -1 + SSMcustom(Z = array(x, c(1, 1, 250)), T = 0, P1 = NA, Q = NA), H = 0)
  
  updatefn <- function(pars, model){
    model["Q",etas = 1] <- model["P1", states = 1] <- exp(pars)
    model
  }
  expect_warning(fit <- fitSSM(model = model, updatefn = updatefn,
                               inits = 0, method = "BFGS"), NA)
  expect_warning(out <- KFS(fit$model), NA)
  expect_true(all(out$F > 0))
  expect_equal(rep(exp(-1), 250), as.numeric(coef(out)))
  expect_equal(0, min(out$V))
  expect_equal(0, max(out$V))
  
  y <-  exp(1:250)
  x <- c(1, y[-250])
  model <- SSModel(y ~ -1 + SSMcustom(Z = array(x, c(1, 1, 250)), T = 0, P1 = NA, Q = NA), H = 0)
  
  updatefn <- function(pars, model){
    model["Q",etas = 1] <- model["P1", states = 1] <- exp(pars)
    model
  }
  expect_warning(fit <- fitSSM(model = model, updatefn = updatefn,
                               inits = 0, method = "BFGS"), NA)
  expect_warning(out <- KFS(fit$model), NA)
  expect_true(all(out$F > 0))
  expect_equal(rep(exp(1), 250), as.numeric(coef(out)))
  expect_equal(0, min(out$V))
  expect_equal(0, max(out$V))
  
  y <-  c(exp(1:250 + rnorm(250)), 1 / exp(1:250 + rnorm(250)))
  x <- c(1, y[-500])
  model <- SSModel(y ~ -1 + SSMcustom(Z = array(x, c(1, 1, 500)), T = 0, P1 = NA, Q = NA), H = 0)
  
  updatefn <- function(pars, model){
    model["Q",etas = 1] <- model["P1", states = 1] <- exp(pars)
    model
  }
  expect_warning(fit <- fitSSM(model = model, updatefn = updatefn,
                               inits = 0, method = "BFGS"), NA)
  expect_warning(out <- KFS(fit$model), NA)
  expect_true(all(out$F > 0))
  expect_equal(y/x, as.numeric(coef(out)))
  expect_equal(0, min(out$V))
  expect_equal(0, max(out$V))
  
  y <- c(2.37375526896038, 0.484382763212699, 0.103390991781719, 0.0318635379310975,
         0.00828792079520903, 0.00153946155449415, 0.000301154179764616,
         4.69712084908583e-05, 3.76402566253327e-06, 7.45679822731237e-09,
         -4.79127859471886e-11, 6.98648705955961e-12, -1.00417499195876e-12,
         1.26051794711504e-13, -3.20702012984292e-14, 6.18033553762622e-15,
         -1.84481158083898e-15, 3.26863125217215e-16, -3.5284201289666e-17,
         2.36395697506618e-18, -6.92742107749554e-20, -8.77878686939385e-22,
         -1.19500786691541e-22, -1.50630197243851e-23, -3.34636969232657e-24,
         -7.96696616958478e-25, -2.20460059174297e-25, -5.80636254492355e-26,
         -1.62283713164705e-26, -3.89630187347358e-27, -7.65190352174139e-28,
         -8.16966714420973e-29, -5.33388991322569e-30, -1.21520471338705e-31,
         1.2649544736029e-33, -1.55703175170465e-34, 2.38264944133605e-35,
         -4.10986347249571e-36, 7.94913694346775e-37, -1.92614292526717e-37,
         3.51989768624077e-38, -5.05154871659138e-39, 9.34466397920635e-40,
         -7.56430326238475e-41, -3.05693230494392e-42, -2.67174949012045e-44,
         -4.13949781450607e-45, -6.11119054914636e-46, -1.22801610501944e-46,
         -2.76803280116606e-47, -7.31144664065532e-48, -1.93726128368889e-48,
         -5.54464145295472e-49, -1.76601788089866e-49, -4.01627339184745e-50,
         -7.59022386085151e-51, -8.36583756696038e-52, -3.17794222502326e-53,
         -3.91839200240157e-55, 7.4825810032649e-57, -7.85738322493243e-58,
         1.46101607031008e-58, -2.57935926446874e-59, 4.19176143504855e-60,
         -7.08214683603497e-61, 1.47527215452851e-61, -1.91709438431942e-62,
         1.87544440454051e-63, 2.04613490688312e-65, 2.09853743248338e-66,
         1.74212391380961e-67, 3.18497678876087e-68, 7.11526882827013e-69,
         2.1213409671844e-69, 5.8718717895392e-70, 2.13453010217269e-70,
         7.42558169209693e-71, 2.06452658841298e-71, 6.1769680972564e-72,
         1.31550698985567e-72, 2.62563134827229e-73, 9.08272465548599e-75,
         1.73124107862413e-76, -9.44208179051201e-78, 4.00647571435298e-79,
         -3.88035674132018e-80, 4.54048352200356e-81, -5.82602591442553e-82,
         8.35180098240723e-83, -1.41883135784571e-83, 1.08530922456612e-84,
         -8.11354927154577e-86, -2.27995301915367e-87, -2.43711722034897e-88,
         -4.44653147227123e-89, -9.1784555786123e-90, -2.43523671946446e-90,
         -7.19606453989095e-91, -3.05540428315533e-91, -1.21223662471548e-91,
         -3.62080301487568e-92, -1.17815112994887e-92, -3.54185759441524e-93,
         -9.39722009745252e-94, -1.73442103553917e-94, -2.29186899966492e-95,
         -1.80564717147145e-96, -2.59730633059331e-98, 1.66278932172853e-99,
         -1.54419579745688e-100, 1.99845930282088e-101, -2.57018963744413e-102,
         3.72902594267587e-103, -4.60600152822731e-104, 2.73666214498223e-105,
         2.35310108382603e-107, 1.04812015244585e-108, 1.31721036975048e-109
  )
  Z <- array(0, dim = c(1, 2, 118))
  Z[1, 1, ] <- c(9.469247, y[-118])
  model <- SSModel(y ~ -1 + SSMcustom(Z = Z, T = diag(c(0, 1)), Q = diag(c(NA, 0)),
                                      P1 = diag(c(NA, 0)), a1 = c(0, 1)), H = NA)
  updatefn <- function(pars, model){
    model["H"] <- exp(pars[1])
    model["Q",etas = 1] <- model["P1", states = 1] <- exp(pars[2])
    model
  }
  expect_warning(fit <- fitSSM(model = model, updatefn = updatefn,
                            inits = c(-4, -4), method = "BFGS"), NA)
  expect_warning(out <- KFS(fit$model), NA)
  expect_equal(out$logLik, 14530.62, tolerance = tol, check.attributes = FALSE)
  expect_equal(fit$model["H",1], 0, tolerance = tol, check.attributes = FALSE)
  expect_equal(fit$model["Q", eta = 1], 0.03556689, tolerance = tol,
               check.attributes = FALSE)

})