search

steve-the-bayesian / BOOM

A C++ library for Bayesian modeling, mainly through Markov chain Monte Carlo, but with a few other methods supported. BOOM = "Bayesian Object Oriented Modeling". It is also the sound your computer makes when it crashes.
GNU Lesser General Public License v2.1
35 stars 14 forks source link

Semilocal Linear Trend Possible Issue #52

Closed alexiosg closed 4 years ago

alexiosg commented 4 years ago

Hi Steve,

I'm trying to understand how the semilocal linear trend is implemented. Based on the documentation, It seems that delta[t+1] = phi x delta[t] + D x (1 - phi)

So I have tried to replicate the states in the example below:

` data("AirPassengers") library(bsts) library(dlm)

spec <- AddSemilocalLinearTrend(list(), y = as.numeric(AirPassengers)) spec <- AddSeasonal(spec, y = as.numeric(AirPassengers), nseasons = 12) opt <- BstsOptions() opt$save.full.state <- TRUE

mod <- bsts(as.numeric(AirPassengers), spec, niter = 2000, model.options = opt)

Now to create the state matrices phi <- mean(mod$trend.slope.ar.coefficient) D <- mean(mod$trend.slope.mean) Z <- matrix(c(1,0,0,1,rep(0,10)), nrow = 1) G <- matrix(0, ncol = 14, nrow = 14) G[1,1:2] <- 1 G[2,2:3] <- c(phi, 1 - phi) G[3, 3] <- 1 G[4,4:14] <- -1 diag(G[5:14, 4:13]) <- 1

W <- diag(0, 14, 14) W[1,1] <- mean(mod$trend.level.sd^2) W[2,2] <- mean(mod$trend.slope.sd^2) W[4,4] <- mean(mod$sigma.seasonal.12^2)

V <- matrix(mean(mod$sigma.obs^2),1,1)

set initial state to mean of initial state in object m0 <- colMeans(mod$full.state[,,1]) C0 <- diag(1e12, 14,14) C0[3,3] <- 0 use dlm to filter the data and states dmod <- dlm(list(FF = Z, GG = G, V = V, W = W, m0 = m0, C0 = C0)) f <- dlmFilter(as.numeric(AirPassengers), dmod) compare DLM 1-step ahead state estimates, with bsts full state and simple approach for first prediction cbind(coredata(f$a)[1,],colMeans(mod$full.state[,,2]),matrix(G %*% colMeans(mod$full.state[,,1]), ncol = 1)) ` [1,] 144.85 144.96 144.85 [2,] -1.03 -6.47 -1.03 [3,] 1.75 1.75 1.75 [4,] -27.18 -27.24 -27.18 [5,] -39.33 -39.33 -39.33 [6,] -34.42 -34.42 -34.42 [7,] -43.73 -43.73 -43.73 [8,] -13.83 -13.83 -13.83 [9,] 23.92 23.92 23.92 [10,] 50.87 50.87 50.87 [11,] 55.34 55.34 55.34 [12,] 33.11 33.11 33.11 [13,] 3.64 3.64 3.64 [14,] -1.90 -1.90 -1.90

So my question relates to the slope (row 2) (and row 2 of the G matrix)...am I doing something wrong in the way I am representing the slope component equation?

Thanks!

Alexios

alexiosg commented 4 years ago

So, after a little more digging it seems that the "full.state" are the smoothed not the filtered states, so this would explain the disparity. Closing this as it seems resolved.