# NOT RUN {
## Note: The following are toy examples, for a proper analysis we recommend to run
## at least 1000 iterations and check the convergence of the Markov chain
## Example 1 (daily data, d = 3, local level + seasonal + covariates)
y <- cbind(seq(0.5,100,by=0.5)*0.1 + rnorm(200),
seq(100.25,150,by=0.25)*0.05 + rnorm(200),
seq(1,200,by=1)*(-0.01) + rnorm(200, 0, 0.5))
dates <- seq.Date(from = as.Date('2019-01-10'),by = "days", length.out = 200)
# Adding a fictional intervention and four covariates. To illustrate the
# functioning of Bayesian model selection, one covariate is assumed to be
# unrelated to y.
int.date <- as.Date('2019-04-01')
y.new <- y; y.new[dates >= int.date, ] <- y.new[dates >= int.date, ]*1.3
x1 <- y[,1]*0.5 + y[,2]*0.3 + y[,3]*0.1
x2 <- y[,2]*0.1 + rnorm(dim(y)[1],0,0.5)
x3 <- y[,3]*1.2 + rnorm(dim(y)[1],0,0.5)
x4 <- rnorm(dim(y)[1], 5, 10)
X <- cbind(x1, x2, x3, x4)
# Model definition
causal.1 <- CausalMBSTS(y.new, components = c("trend", "seasonal"), seas.period = 7,
X = X, dates = dates, int.date = int.date,
s0.r = 0.1*diag(3), s0.eps = 0.1*diag(3), niter = 50,
burn = 5, horizon = as.Date(c('2019-04-08','2019-07-28')))
## Plotting
plot(causal.1, int.date = int.date, type = 'inclusion.probs', prob = 0.1)
# as expected, x4 is rarely included in the model
oldpar <- par(no.readonly = TRUE)
par(mar = c(2,2,2,2))
par(mfrow=c(2,3))
plot(causal.1, int.date = int.date, type = c('impact', 'forecast'))
par(mfrow=c(3,4))
plot(causal.1, type = 'ppchecks', int.date = int.date)
par(oldpar)
## Example 2
set.seed(1)
t <- seq(from = 0,to = 4*pi, length.out=300)
y <- cbind(3*sin(2*t)+rnorm(300), 2*cos(2*t) + rnorm(300))
dates <- seq.Date(from = as.Date("2015-01-01"), by = "week", length.out=300)
int.date <- as.Date("2020-02-27")
y[dates >= int.date,] <- y[dates >= int.date,]+2
# Model definition
causal.2 <- CausalMBSTS(y, components = c("trend", "cycle"), cycle.period = 75,
dates = dates, int.date = int.date,
s0.r = 0.01*diag(2), s0.eps = 0.1*diag(2),
niter = 100, burn = 10)
# Plotting
oldpar <- par(no.readonly = TRUE)
par(mfrow=c(2,4))
plot(causal.2, type = 'ppchecks', int.date = int.date)
par(mfrow=c(2,2))
plot(causal.2, type = c('impact','forecast'), int.date = int.date)
par(oldpar)
# }
Run the code above in your browser using DataLab