bsm_ng  R Documentation 
Constructs a nonGaussian basic structural model with local level or local trend component, a seasonal component, and regression component (or subset of these components).
bsm_ng(
y,
sd_level,
sd_slope,
sd_seasonal,
sd_noise,
distribution,
phi,
u,
beta,
xreg = NULL,
period,
a1 = NULL,
P1 = NULL,
C = NULL
)
y 
A vector or a 
sd_level 
Standard deviation of the noise of level equation.
Should be an object of class 
sd_slope 
Standard deviation of the noise of slope equation.
Should be an object of class 
sd_seasonal 
Standard deviation of the noise of seasonal equation.
Should be an object of class 
sd_noise 
A prior for the standard deviation of the additional noise
term to be added to linear predictor, defined as an object of class

distribution 
Distribution of the observed time series. Possible
choices are 
phi 
Additional parameter relating to the nonGaussian distribution.
For negative binomial distribution this is the dispersion term, for gamma
distribution this is the shape parameter, and for other distributions this
is ignored. Should an object of class 
u 
A vector of positive constants for nonGaussian models. For Poisson, gamma, and negative binomial distribution, this corresponds to the offset term. For binomial, this is the number of trials. 
beta 
A prior for the regression coefficients.
Should be an object of class 
xreg 
A matrix containing covariates with number of rows matching the
length of 
period 
Length of the seasonal pattern.
Must be a positive value greater than 2 and less than the length of the
input time series. Default is 
a1 
Prior means for the initial states (level, slope, seasonals). Defaults to vector of zeros. 
P1 
Prior covariance matrix for the initial states (level, slope, seasonals).Default is diagonal matrix with 100 on the diagonal. 
C 
Intercept terms for state equation, given as a m x n or m x 1 matrix. 
An object of class bsm_ng
.
# Same data as in Vihola, Helske, Franks (2020)
data(poisson_series)
s < sd(log(pmax(0.1, poisson_series)))
model < bsm_ng(poisson_series, sd_level = uniform(0.115, 0, 2 * s),
sd_slope = uniform(0.004, 0, 2 * s), P1 = diag(0.1, 2),
distribution = "poisson")
out < run_mcmc(model, iter = 1e5, particles = 10)
summary(out, variable = "theta", return_se = TRUE)
# should be about 0.093 and 0.016
summary(out, variable = "states", return_se = TRUE,
states = 1, times = c(1, 100))
# should be about 0.075, 2.618
model < bsm_ng(Seatbelts[, "VanKilled"], distribution = "poisson",
sd_level = halfnormal(0.01, 1),
sd_seasonal = halfnormal(0.01, 1),
beta = normal(0, 0, 10),
xreg = Seatbelts[, "law"],
# default values, just for illustration
period = 12L,
a1 = rep(0, 1 + 11), # level + period  1 seasonal states
P1 = diag(1, 12),
C = matrix(0, 12, 1),
u = rep(1, nrow(Seatbelts)))
set.seed(123)
mcmc_out < run_mcmc(model, iter = 5000, particles = 10, mcmc_type = "da")
mcmc_out$acceptance_rate
theta < expand_sample(mcmc_out, "theta")
plot(theta)
summary(theta)
library("ggplot2")
ggplot(as.data.frame(theta[,1:2]), aes(x = sd_level, y = sd_seasonal)) +
geom_point() + stat_density2d(aes(fill = ..level.., alpha = ..level..),
geom = "polygon") + scale_fill_continuous(low = "green", high = "blue") +
guides(alpha = "none")
# Traceplot using as.data.frame method for MCMC output
library("dplyr")
as.data.frame(mcmc_out) >
filter(variable == "sd_level") >
ggplot(aes(y = value, x = iter)) + geom_line()
# Model with slope term and additional noise to linear predictor to capture
# excess variation
model2 < bsm_ng(Seatbelts[, "VanKilled"], distribution = "poisson",
sd_level = halfnormal(0.01, 1),
sd_seasonal = halfnormal(0.01, 1),
beta = normal(0, 0, 10),
xreg = Seatbelts[, "law"],
sd_slope = halfnormal(0.01, 0.1),
sd_noise = halfnormal(0.01, 1))
# instead of extra noise term, model using negative binomial distribution:
model3 < bsm_ng(Seatbelts[, "VanKilled"],
distribution = "negative binomial",
sd_level = halfnormal(0.01, 1),
sd_seasonal = halfnormal(0.01, 1),
beta = normal(0, 0, 10),
xreg = Seatbelts[, "law"],
sd_slope = halfnormal(0.01, 0.1),
phi = gamma_prior(1, 5, 5))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.