rstandard.KFS {KFAS}R Documentation

Extract Standardized Residuals from KFS output

Description

Extract Standardized Residuals from KFS output

Usage

## S3 method for class 'KFS'
rstandard(
  model,
  type = c("recursive", "pearson", "state"),
  standardization_type = c("marginal", "cholesky"),
  zerotol = 0,
  expected = FALSE,
  ...
)

Arguments

model

KFS object

type

Type of residuals. See details.

standardization_type

Type of standardization. Either "marginal" (default) for marginal standardization or "cholesky" for Cholesky-type standardization.

zerotol

Tolerance parameter for positivity checking in standardization. Default is zero. The values of D <= zerotol * max(D, 0) are deemed to zero.

expected

Logical value defining the approximation of H_t in case of Gamma and negative binomial distribution. Default is FALSE which matches the algorithm of Durbin & Koopman (1997), whereas TRUE uses the expected value of observations in the equations, leading to results which match with glm (where applicable). The latter case was the default behaviour of KFAS before version 1.3.8. Essentially this is the difference between observed and expected information.

...

Ignored.

Details

For object of class KFS with fully Gaussian observations, several types of standardized residuals can be computed. Also the standardization for multivariate residuals can be done either by Cholesky decomposition L^{-1}_t residual_t, or component-wise residual_t/sd(residual_t),.

Note that the variance estimates from KFS are of form Var(x | y), e.g., V_eps from KFS is Var(\epsilon_t | Y) and matches with equation 4.69 in Section 4.5.3 of Durbin and Koopman (2012). (in case of univariate Gaussian model). But for the standardization we need Var(E(x | y)) (e.g., Var(epshat) which we get with the law of total variance as H_t - V_eps for example.

Examples

# Replication of residual plot of Section 8.2 of Durbin and Koopman (2012)
model <- SSModel(log(drivers) ~ SSMtrend(1, Q = list(NA))+
    SSMseasonal(period = 12, sea.type = "trigonometric", Q = NA),
  data = Seatbelts, H = NA)

updatefn <- function(pars, model){
  model$H[] <- exp(pars[1])
  diag(model$Q[, , 1]) <- exp(c(pars[2], rep(pars[3], 11)))
  model
}

fit <- fitSSM(model = model,
  inits = log(c(var(log(Seatbelts[, "drivers"])), 0.001, 0.0001)),
  updatefn = updatefn)

# tiny differences due to different optimization algorithms
setNames(c(diag(fit$model$Q[,,1])[1:2], fit$model$H[1]),
  c("level", "seasonal", "irregular"))

out <- KFS(fit$model, smoothing = c("state", "mean", "disturbance"))

plot(cbind(
  recursive = rstandard(out),
  irregular = rstandard(out, "pearson"),
  state = rstandard(out, "state")[,1]),
  main = "recursive and state residuals", type = "h")

# Figure 2.8 in DK2012
model_Nile <- SSModel(Nile ~
    SSMtrend(1, Q = list(matrix(NA))), H = matrix(NA))
model_Nile <- fitSSM(model_Nile, c(log(var(Nile)), log(var(Nile))),
  method = "BFGS")$model

out_Nile <- KFS(model_Nile,  smoothing = c("state", "mean", "disturbance"))

par(mfrow = c(2, 2))
res_p <- rstandard(out_Nile, "pearson")
ts.plot(res_p)
abline(a = 0, b= 0, lty = 2)
hist(res_p, freq = FALSE)
lines(density(res_p))
res_s <- rstandard(out_Nile, "state")
ts.plot(res_s)
abline(a = 0, b= 0, lty = 2)
hist(res_s, freq = FALSE)
lines(density(res_s))


[Package KFAS version 1.5.1 Index]