plot.CausalMBSTS {CausalMBSTS}

Plotting function for object of class CausalMBSTS

Description

Given an object of class 'CausalMBSTS', the function draws: i) the plot of the estimated (pointwise) causal effect; ii) the original time series plotted against the predicted counterfactual; iii) posterior predictive checks; iv) regressor inclusion probabilities (only for models with a regression component).

Usage

## S3 method for class 'CausalMBSTS'
plot(x, int.date, type = c("impact", "forecast", "ppchecks"), prob = NULL, ...)

Arguments

x

Object of class 'CausalMBSTS'

int.date

Date of the intervention.

type

A character string indicating the type of plot to be produced. Possible values in 'c('impact', 'forecast', 'ppchecks', 'inclusion.probs')'. See Details for further explanation.

prob

Regressors inclusion probabilities above 'prob' are plotted. Optional, only required for type = 'inclusion.prob'.

...

Arguments passed to other methods (currently unused).

Details

Option 'impact' for parameter type plots the general causal effect at every time points in the post period. Multiple plots will be generated, corresponding to each combination of time series and horizon (if specified). Option 'forecast' plots the observed time series against the predicted counterfactual, one plot per each combination of time series and horizon (if specified). 'ppchecks' draws posterior predictive checks for the model estimated in the pre-period. They include four plots generated for each time series (and horizon). The plots are (1) density of posterior mean vs. density of the data before intervention, (2) Histogram of maximum in-sample forecasts and Bayes p-value, (3) QQ-plot of residuals, and (4) ACF of residuals. Option 'inclusion.probs' plots the regressors' inclusion probabilities above 'prob'.

Value

NULL, invisibly.

Examples

## 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)