Why should standard auxiliary residuals ever be NA?

[Figuera]

Hi, I am replicating a little program of mine that was built with dlm package with KFAS (which, in my experience, is much faster, congratulations).

In this program I use auxiliary residuals to detect interventions. In dlm package there is no function (as far as I know) to calculate auxiliary residuals so I had to create my own, but with KFAS i would rather use the package rstandard function.

Comparing the two functions I get, as expected, very similar results, but for some cases rstandard returns just a bunch of NA values while my function return actual numbers. This mostly happen when I am calculating "state" residuals and the variance is very close to 0.

Why should this happen? Why should a standard auxiliary residual ever be NA?

Inspecting you code is easy to see what is happening, in state_standardized function we have:

D <- sqrt((object$model$Q[, , i * tvq + 1] - object$V_eta[, , i])[1 + 0:(k - 1) * (k + 1)])
pos <- D > sqrt(.Machine$double.eps) * max(D, 0)
eta[i, pos] <- eta[i, pos] / D[pos]
eta[i, !pos] <- NA

eta[i, ] is NAbecause D is less than sqrt(.Machine$double.eps) * max(D, 0). I just fail to see the meaning of this, what are we testing? Is this to avoid some bug?

In my specific case: etahat[1,3] = 1.44e-14 and D[3] = 1.16e-13, I was expecting to get standard eta as 1.44e-14/1.16e-13 = 0.12, but I get NA.

[helske]

To be honest I have no idea why there is that kind of check in the codes. D can be zero for some states, and even nanin case of numerical issues, but I don't really remember why I have used that extra guard against small values of D.

[helske]

I now added an argument called zerotol for rstandard function which you can use to define new tolerance parameter instead of sqrt(machine epsilon). By default, this is now 0. Let me know if you still see some issues.

[Figuera]

Very fast response, I think that will do, thank you.