helske / KFAS

KFAS: R Package for Exponential Family State Space Models
64 stars 17 forks source link

No effect of missing observations on seasonal estimation error variance? #48

Closed k41m4n closed 4 years ago

k41m4n commented 4 years ago

Dear Mr Helske, I have run with KFAS a stochastic level and deterministic seasonal model with missing observations as presented in the book "An Introduction to State Space Time Series Analysis" (pages 103 - 106) of J.F. Commandeur and S.J Koopman (see my code with data below). In my results, it can be seen that the level estimation error variance increases at time points with the missing observations as presented in the book. It seems, however, that in my results, the seasonal estimation error variance does not increase for those time points as it is shown in the book. Should we also expect to get increases for the seasonal estimation error variance at time points with missing observations?

`

Libraries

library(KFAS) if(!(require(rstudioapi))){install.packages('rstudioapi')} if(!(require(dplyr))){install.packages('dplyr')} library(dplyr)

Loading data and treating some observations as missing

dataUKdriversKSI <- read.table(text = " 1687 1508 1507 1385 1632 1511 1559 1630 1579 1653 2152 2148 1752 1765 1717 1558 1575 1520 1805 1800 1719 2008 2242 2478 2030 1655 1693 1623 1805 1746 1795 1926 1619 1992 2233 2192 2080 1768 1835 1569 1976 1853 1965 1689 1778 1976 2397 2654 2097 1963 1677 1941 2003 1813 2012 1912 2084 2080 2118 2150 1608 1503 1548 1382 1731 1798 1779 1887 2004 2077 2092 2051 1577 1356 1652 1382 1519 1421 1442 1543 1656 1561 1905 2199 1473 1655 1407 1395 1530 1309 1526 1327 1627 1748 1958 2274 1648 1401 1411 1403 1394 1520 1528 1643 1515 1685 2000 2215 1956 1462 1563 1459 1446 1622 1657 1638 1643 1683 2050 2262 1813 1445 1762 1461 1556 1431 1427 1554 1645 1653 2016 2207 1665 1361 1506 1360 1453 1522 1460 1552 1548 1827 1737 1941 1474 1458 1542 1404 1522 1385 1641 1510 1681 1938 1868 1726 1456 1445 1456 1365 1487 1558 1488 1684 1594 1850 1998 2079 1494 1057 1218 1168 1236 1076 1174 1139 1427 1487 1483 1513 1357 1165 1282 1110 1297 1185 1222 1284 1444 1575 1737 1763 ")[, 1] %>% log() %>% replace(c(48:62, 120:140), NA) %>% ts(start = 1969, frequency=12)

Defining model

model <- SSModel(dataUKdriversKSI ~ SSMtrend(degree=1, Q=list(matrix(NA))) + SSMseasonal(12, sea.type='dummy', Q=matrix(0)), H=matrix(NA))

ownupdatefn <- function(pars,model,...){ model$H[,,1] <- exp(pars[1]) diag(model$Q[,,1]) <- c(exp(pars[2]), 0) model }

Fitting model

w <- 2 #Number of estimated hyperparameters (i.e. disturbance variances) fit <- fitSSM(inits = log(rep(0.005, w)), model = model, updatefn = ownupdatefn, method = "Nelder-Mead") outKFS <- KFS(fit$model, smoothing = c("state", "mean", "disturbance"))

Level and seasonal estimation error variances

levEstErVar <- outKFS$V[1, 1, ] %>% ts(start = 1969, frequency=12) seasEstErVar <- outKFS$V[2, 2, ] %>% ts(start = 1969, frequency=12)

Figure 8.17. Stochastic level estimation error variance for log drivers KSI

with observations at t=48,...,62 and t=120,...,140 treated as missing

plot(levEstErVar, xlab = "", ylab = "", lty = 1) title(main = "Figure 8.17. Stochastic level estimation error variance for log drivers KSI with observations at t=48,...,62 \n and t=120,...,140 treated as missing", cex.main = 0.8) legend("topleft",leg = "level estimation error variance", cex = 0.6, lty = 1, horiz = T)

Figure 8.19. Seasonal estimation error variance for log drivers KSI

with observations at t=48,...,62 and t=120,...,140 treated as missing

plot(seasEstErVar, xlab = "", ylab = "", lty = 1) title(main = "Figure 8.17. Seasonal estimation error variance for log drivers KSI with observations at t=48,...,62 \n and t=120,...,140 treated as missing", cex.main = 0.8) legend("topleft",leg = "seasonal estimation error variance", cex = 0.6, lty = 1, horiz = T)

`

helske commented 4 years ago

I don't have that book at hand but so I can't compare the results, but in your code, you set the variance of the seasonal component to zero, i.e. seasonal component is non-stochastic. This to me would suggest that the error variance of the seasonal component shouldn't increase with missing observations?

k41m4n commented 4 years ago

I created another model in which the variance of the seasonal component was allowed to change over time, i.e. seasonal component was stochastic. In this model, I also could not observe any effect of missing observations on the error variance of the seasonal component. These are screenshots of J.F. Commandeur and S.J Koopman's book showing the reuslts for the stochastic level and deterministic seasonal model: IMG_1953 IMG_1952

These are my results for the same model: Rplot2 Rplot1

It seems that the plot patern for the level estimation error variance in the book resembles my plot patern for the covariance between level estimation errors and seasonal estimation errors (outKFS$V[1, 2, ]). In my results are also negative values. Rplot3

helske commented 4 years ago

Hmm, this is quite strange. I tested the same model using the bssm package (using as_bssm and smoother functions) and the results coincide with KFAS. bssm has an independent (almost, same author...) implementation of Kalman filter/smoother soI lean on the possibility of an error in the figure. Might be a good idea to try the same model also with say dlm package...

k41m4n commented 4 years ago

Thank you for your prompt reply and test in the bssm package. I will also try to check it using dlm.