Bug 14738 - Wrong sign in KalmanSmooth computing the smoothed state variance
Wrong sign in KalmanSmooth computing the smoothed state variance
Status: CLOSED FIXED
Product: R
Classification: Unclassified
Component: Models
R 2.14.0
Other Linux
: P5 normal
Assigned To: R-core
Depends on:
Blocks:
  Show dependency treegraph
 
Reported: 2011-11-23 20:19 UTC by Javier López-de-Lacalle
Modified: 2011-11-23 22:22 UTC (History)
0 users

See Also:


Attachments
Graphic displaying confidence bands based on the current and a corrected version of KalmanSmooth (16.76 KB, application/pdf)
2011-11-23 20:19 UTC, Javier López-de-Lacalle
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Javier López-de-Lacalle 2011-11-23 20:19:10 UTC
Created attachment 1250 [details]
Graphic displaying confidence bands based on the current and a corrected version of KalmanSmooth

Dear R developers:

I am using the function KalmanSmooth to replicate the top-left chart in Figure 2.2 of the book Durbin and Koopman 2001 (cited in the documentation of the function KalmanSmooth).

The following code computes the smoothed level and the 95% confidence bands based on a local level model fitted to the Nile time series.

#begin R code
y <- Nile
vars <- StructTS(y, type = "level")$coef
mod <- list(Z = cbind(1), a = y[1], P = var(y)*10000,
  T = 1, V = vars[1], h = vars[2], Pn = var(y)*10000)

ks <- KalmanSmooth(y, mod) 

ub <- ts(ks$smooth[,1] + 1.96 * sqrt(ks$var[,1,1]), start = start(y))
lb <- ts(ks$smooth[,1] - 1.96 * sqrt(ks$var[,1,1]), start = start(y))

plot(ts(ks$smooth[,1], start = start(y)), ylim =  c(456, 1410), 
main = "Nile time series. Smoothed level and confidence interval", type = "n", yaxt = "n")
axis(side = 2, at = seq(500, 1400, 150), labels = seq(500, 1400, 150), las = 1)
lines(y, col = "lightgray")
lines(ts(ks$smooth[,1], start = start(y)))
lines(ub, lty = 2, col = "blue")
lines(lb, lty = 2, col = "blue")
#end R code

The confidence bands are noticeably wider than those shown in Figure 2.2 of the reference given above.

Inspecting the code, there is an error in the sign of one of the calculations.

The location of the code I refer to is the line 335 of the file /R-2.14.0/src/library/stats/src/arima.c.

It should be:

tmp -= mm[i + p * k] * Pt[l + n*k + n*p*j];

instead of:

tmp += mm[i + p * k] * Pt[l + n*k + n*p*j];

It can be easily checked that the sign should be a minus. According to equation (4.32) in Durbin and Koopman 2001 we have:

V_t = P_t - P_t * N_{t-1} * P_t.

In the code, 'tmp' is initialized as P_t and 'mm' is N_{t-1} * P_t. Therefore, the correct sign is a minus (tmp -= ...).

I attach a graphic showing the results based on the current version of KalmanSmooth and a modified/corrected version.

Regards

javi