On Sat, 03 Apr 2004 21:34:25 -0400, you wrote:>Defining the following: > > set.seed(123) > > kalmanTorture <- function(iter) { > x <- arima.sim(model = list(ar=0.9, ma=0.5),n=150 ) > x[10:20] <- NA > mod <- arima(x, order=c(1,0,1) ) > for (i in 1:iter) { > smooth <- KalmanSmooth(x, mod=mod$model)$smooth > if (any(is.na(smooth))) stop("NA on iteration ",i) } >} > >and running in R1.9.0beta (equal in R1.8.1), windows XP, precompiled >binaries: > >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 147 >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 3 >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 2 >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 3 >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 1 >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 4 >> kalmanTorture(1000) >Error in kalmanTorture(1000) : NA on iteration 1 >. >. >. >Starting a new session, using the same set.seed(123) >gives another sequence of iteration numbers where failing.I can reproduce this. The problem happens in arima.c where the KalmanSmooth function is, but I can't see why yet. Duncan Murdoch
On Sun, 4 Apr 2004 16:24:22 +0200 (CEST), Kjetil Halvorsen wrote:>On Sat, 03 Apr 2004 21:34:25 -0400, you wrote: > >>Defining the following: >> >> set.seed(123) >> >> kalmanTorture <- function(iter) { >> x <- arima.sim(model = list(ar=0.9, ma=0.5),n=150 ) >> x[10:20] <- NA >> mod <- arima(x, order=c(1,0,1) ) >> for (i in 1:iter) { >> smooth <- KalmanSmooth(x, mod=mod$model)$smooth >> if (any(is.na(smooth))) stop("NA on iteration ",i) } >>} >> >>and running in R1.9.0beta (equal in R1.8.1), windows XP, precompiled >>binaries: >> >>> kalmanTorture(1000) >>Error in kalmanTorture(1000) : NA on iteration 147 >>> kalmanTorture(1000)Found it. The entries in the Mt array in KalmanSmooth weren't being initialized in the cases where data was missing. The value was being multiplied by zero, so if it was numerically valid junk it had no effect, but if it happened to be a NaN, the multiplication gave another NaN. I'll commit a patch. Duncan Murdoch