I am using the bsts package to analyze several time series, to find out whether the values in the series are increasing, decreasing or remaining stable along the time period.
I see I can extract a trend component from the model, which looks like it is increasing (code below). But how can I formally test the direction of the trend? Should I remove the local linear trend and instead do a regression against time? Or is there a way to answer that question using the trend component already extracted?
Here is one of the time series:
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2013 5.454450 5.351702 5.450414 5.490382 5.367216 5.404927 5.267974 5.211403 5.439326 5.394489 5.425333 5.413367
2014 5.365040 5.359957 5.388980 5.272279 5.337346 5.383114 5.211803 5.567671 5.446918 5.427461 5.490386 5.429216
2015 5.522443 5.432766 5.581413 5.307201 5.452850 5.426128 5.372678 5.524919 5.400167 5.453443 5.596288 5.424220
2016 5.633506 5.553917 5.553649 5.444297 5.551022 5.461216 5.503183 5.426033 5.454331 5.643385 5.575780 5.486073
2017 5.656180 5.411885 5.477615 5.437005 5.434144 5.404344 5.481005 5.437255 5.352800 5.485022 5.441534 5.472936
2018 5.366347 5.381583 5.401431 5.479976 5.439315 5.319484 5.421672 5.319448 5.574673 5.472161 5.479900 5.539380
2019 5.390139 5.429426 5.613977 5.343529 5.487940 5.624361 5.381285 5.366611 5.565688 5.503865 5.491821 5.486168
2020 5.475343 5.493866 6.010556 5.690947 5.656557 5.420500 5.453484 5.566972 5.369799 5.435967 5.613358 5.409345
And the code I used to extract the trend (from https://multithreaded.stitchfix.com/blog/2016/04/21/forget-arima/), y is the time series above:
ss <- AddLocalLinearTrend(list(), y)
ss <- AddSeasonal(ss, y, nseasons = 12)
bsts.model <- bsts(y, state.specification = ss, niter = 500, ping=0, seed=2016)
burn <- SuggestBurn(0.1, bsts.model)
components <- cbind.data.frame(
colMeans(bsts.model$state.contributions[-(1:burn),"trend",]),
colMeans(bsts.model$state.contributions[-(1:burn),"seasonal.12.1",]),
as.Date(time(Y)))
names(components) <- c("Trend", "Seasonality", "Date")
components <- melt(components, id="Date")
names(components) <- c("Date", "Component", "Value")
ggplot(data=components, aes(x=Date, y=Value)) + geom_line() +
theme_bw() + theme(legend.title = element_blank()) + ylab("") + xlab("") +
facet_grid(Component ~ ., scales="free") + guides(colour=FALSE) +
theme(axis.text.x=element_text(angle = -90, hjust = 0))
This can be done by inspecting the draws more closely.
Let's take this apart:
colMeans(bsts.model$state.contributions[-(1:burn), "trend",])
bsts.model$state.contributions is an array of shape (iter, metric, obs).
bsts.model$state.contributions[-(1:burn), "trend", ] Turns this into a (non_burnin_iter, obs) matrix.
The call to colMeans gets the average trend but you can extract any quantile or shares of properties as well:
m <- bsts.model$state.contributions[-(1:burn), "trend", ]
# Average trend, what you are getting now.
colMeans(m)
# Posterior probability that the effect of the local trend is positive
# Probably what you want though note that it's conceptually not the
# same as a p-value. It's more what people think 1-pvalue is.
colMeans(m > 0)
# Lower bound of the 95% credible interval of the trend.
# Useful if you have a region or practical equivalence (ROPE).
apply(m, 2, quantile, p = 0.025)