# Regression model with time-varying coefficients
set.seed(1)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
b1 <- 1 + cumsum(rnorm(n, sd = 0.5))
b2 <- 2 + cumsum(rnorm(n, sd = 0.1))
y <- 1 + b1 * x1 + b2 * x2 + rnorm(n, sd = 0.1)
Z <- rbind(1, x1, x2)
H <- 0.1
T <- diag(3)
R <- diag(c(0, 1, 0.1))
a1 <- rep(0, 3)
P1 <- diag(10, 3)
# updates the model given the current values of the parameters
update_fn <- function(theta) {
R <- diag(c(0, theta[1], theta[2]))
dim(R) <- c(3, 3, 1)
list(R = R, H = theta[3])
}
# prior for standard deviations as half-normal(1)
prior_fn <- function(theta) {
if(any(theta < 0)) {
log_p <- -Inf
} else {
log_p <- sum(dnorm(theta, 0, 1, log = TRUE))
}
log_p
}
model <- ssm_ulg(y, Z, H, T, R, a1, P1,
init_theta = c(1, 0.1, 0.1),
update_fn = update_fn, prior_fn = prior_fn,
state_names = c("level", "b1", "b2"),
# using default values, but being explicit for testing purposes
C = matrix(0, 3, 1), D = numeric(1))
out <- run_mcmc(model, iter = 5000)
out
sumr <- summary(out, variable = "state", times = 1:n)
sumr$true <- c(b1, b2, rep(1, n))
library(ggplot2)
ggplot(sumr, aes(x = time, y = Mean)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), alpha = 0.5) +
geom_line() +
geom_line(aes(y = true), colour = "red") +
facet_wrap(~ variable, scales = "free") +
theme_bw()
# Perhaps easiest way to construct a general SSM for bssm is to use the
# model building functionality of KFAS:
library("KFAS")
model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = 5e-4)+
SSMseasonal(period = 12, sea.type = "trigonometric", Q = 0) +
log(PetrolPrice) + law, data = Seatbelts, H = 0.005)
# use as_bssm function for conversion, kappa defines the
# prior variance for diffuse states
model_bssm <- as_bssm(model_kfas, kappa = 100)
# define updating function for parameter estimation
# we can use SSModel and as_bssm functions here as well
# (for large model it is more efficient to do this
# "manually" by constructing only necessary matrices,
# i.e., in this case a list with H and Q)
prior_fn <- function(theta) {
if(any(theta < 0)) -Inf else sum(dnorm(theta, 0, 0.1, log = TRUE))
}
update_fn <- function(theta) {
model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
SSMseasonal(period = 12,
sea.type = "trigonometric", Q = theta[2]^2) +
log(PetrolPrice) + law, data = Seatbelts, H = theta[3]^2)
# the bssm_model object is essentially list so this is fine
as_bssm(model_kfas, kappa = 100, init_theta = init_theta,
update_fn = update_fn, prior_fn = prior_fn)
}
init_theta <- rep(1e-2, 3)
names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y")
model_bssm <- update_fn(init_theta)
# \donttest{
out <- run_mcmc(model_bssm, iter = 10000, burnin = 5000)
out
# }
# Above the regression coefficients are modelled as
# time-invariant latent states.
# Here is an alternative way where we use variable D so that the
# coefficients are part of parameter vector theta. Note however that the
# first option often preferable in order to keep the dimension of theta low.
updatefn2 <- function(theta) {
# note no PetrolPrice or law variables here
model_kfas2 <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
SSMseasonal(period = 12, sea.type = "trigonometric", Q = theta[2]^2),
data = Seatbelts, H = theta[3]^2)
X <- model.matrix(~ -1 + law + log(PetrolPrice), data = Seatbelts)
D <- t(X %*% theta[4:5])
as_bssm(model_kfas2, D = D, kappa = 100)
}
prior2 <- function(theta) {
if(any(theta[1:3] < 0)) {
-Inf
} else {
sum(dnorm(theta[1:3], 0, 0.1, log = TRUE)) +
sum(dnorm(theta[4:5], 0, 10, log = TRUE))
}
}
init_theta <- c(rep(1e-2, 3), 0, 0)
names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y", "law", "Petrol")
model_bssm2 <- updatefn2(init_theta)
model_bssm2$theta <- init_theta
model_bssm2$prior_fn <- prior2
model_bssm2$update_fn <- updatefn2
# \donttest{
out2 <- run_mcmc(model_bssm2, iter = 10000, burnin = 5000)
out2
# }
Run the code above in your browser using DataLab