Standardized residuals

[k41m4n]

I try to calculate myself , without using function rstandard, standardized residuals following the formulas provided in the documentation for KFAS on pages 44-45 (see my code below) I can replicate the results of function rstandard for residuals ("recursive") but not for auxiliary residuals ("pearson", "state"). It seem that I apply the formulas correctly. Are there any aditional adjustments made in function rstandard?

data(UKDriverDeaths) dataUKdriversKSI <- log(UKDriverDeaths) model <- SSModel(dataUKdriversKSI ~ SSMtrend(degree = 1, Q = list(matrix(NA))), H = matrix(NA))

ownupdatefn <- function(pars,model){ model$H[,, 1] <- exp(pars[1]) model$Q[,, 1] <- exp(pars[2]) model } fit <- fitSSM(model, inits = log(c(0.001, 0.001)), updatefn = ownupdatefn, method = "BFGS") outKFS <- KFS(fit$model, smoothing = c("state", "mean", "disturbance"))

residuals - one-ahead predictions

predResidSt <- c(rstandard(outKFS, type = "recursive"))

predResid <- c(residuals(outKFS, type = "recursive")) predResidVar <- c(outKFS$F) predResidStCal <- predResid/sqrt(predResidVar)

(diffPredResi <- data.frame(predResidSt, predResidStCal, diff = predResidSt-predResidStCal))

auxiliary residuals - irregular disturbances

irregResidSt <- c(rstandard(outKFS, type = "pearson"))

irregResid <- c(residuals(outKFS, type = "pearson")) irregResidVar <- c(outKFS$V_eps) irregResidStCal <- irregResid/sqrt(irregResidVar)

(diffIrregResid <- data.frame(irregResidSt, irregResidStCal, diff = irregResidSt-irregResidStCal))

auxiliary residuals - state disturbances

disturResidSt <- c(rstandard(outKFS, type = "state"))

disturResid <- c(residuals(outKFS, type = "state")) disturResidVar <- c(outKFS$V_eta) disturResidStCal <- disturResid/sqrt(disturResidVar)

(diffDisturResid <- data.frame(disturResidSt, disturResidStCal, diff = disturResidSt-disturResidStCal))

[helske]

The variances for the last two residuals are not correct, they should be H - V_eps and Q - v_eta.

[k41m4n]

Thank you for your reply. Could you please be so kind to explain why H - V_eps and Q - V_eta, or could you refer to some documentation or literature which could give more information on the standardization formula used? It seems that H - V_eps and Q - V_eta are not explicitly mentioned in the KFAS documentation concerning the standardization of residuals.

[helske]

The V_eps and V_eta returned by KFAS are V(eps_t | Y) and V(eta_t |Y) as in Durbin and Koopman 2012 equation 4.69, but for standardised residuals we use hat(eps_t)/sqrt(Var(hat(eps_t)), where hat(eps_t) = E(eps_t | Y). So by law of total variance Var(E(eps_t | Y)) = Var(eps_t) - E(Var(eps_t | Y)) = H_t - V_eps and similarly for state residuals. And because V_eps has a form of H_t - X, we are in the end scaling with X which matches with the equation 12 in https://www.jstatsoft.org/article/view/v041i01 (again similarly for state residuals).

Here's an example which shows that these definitions match with DK2012:

# 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")

Thanks for this, I'll add some text regarding these to the rstandard documentation.

[k41m4n]

Does it mean that E(Var(eps_t | Y)) = Var(eps_t | Y) because Y is fixed?

[helske]

Yeah, my notation was bit sloppy, but indeed eps and Y are multivariate normal and the conditional variance Var(eps|Y=y) does not depend on y so E(Var(eps | Y=y)) = Var(eps | Y=y).

Edit: Google pointed me to this pdf which has some explicit derivations for residuals which you might find helpful: https://arxiv.org/abs/1411.0045.

[k41m4n]

Thank you for your explanations which helped me to understand how the standardized auxiliary residuals are calculated.