BSL (version 3.0.0)

ma2: An MA(2) model

Description

In this example we wish to estimate the parameters of a simple MA(2) time series model. We provide the data and tuning parameters required to reproduce the results in An et al. (2019).

Usage

data(ma2)
ma2_sim(theta, T)
ma2_sum(x)
ma2_prior(theta)
ma2_logPrior(theta)

Arguments

theta

A vector of proposed model parameters, \(\theta_1\) and \(\theta_2\).

T

The number of observations.

x

Observed or simulated data in the format of a vector of length \(T\).

A simulated dataset

An example ``observed'' dataset and the tuning parameters relevant to that example can be obtained using data(ma2). This ``observed'' data is a simulated dataset with \(\theta_1 = 0.6\), \(\theta_2=0.2\) and \(T=50\). Further information about this model and the specific choices of tuning parameters used in BSL and BSLasso can be found in An et al. (2019).

• data: A time series dataset, in the form of a vector of length \(T\)

• sim_options: A list containing \(T=50\)

• start: A vector of suitable initial values of the parameters for MCMC

• cov: The covariance matrix of a multivariate normal random walk proposal distribution used in the MCMC, in the form of a 2 by 2 matrix

Details

This example is based on estimating the parameters of a basic MA(2) time series model of the form

$$y_t = z_t + \theta_1 z_{t-1} + \theta_2 z_{t-2},$$

where \(t=1,\ldots,T\) and \(z_t ~ N(0,1)\) for \(t=-1,0,\ldots,T\). A uniform prior is used for this example, subject to the restrictions that \(-2<\theta_1<2\), \(\theta_1+\theta_2>-1\) and \(\theta_1-\theta_2<1\) so that invertibility of the time series is satisfied. The summary statistics are simply the full data.

References

An, Z., South, L. F., Nott, D. J. & Drovandi, C. C. (2019). Accelerating Bayesian synthetic likelihood with the graphical lasso. Journal of Computational and Graphical Statistics. https://doi.org/10.1080/10618600.2018.1537928

Examples

# NOT RUN {
# Loading the data for this example
data(ma2)
model <- BSLModel(fnSim = ma2_sim, fnSum = ma2_sum, simArgs = ma2$sim_options, theta0 = ma2$start, 
                  fnLogPrior = ma2_logPrior)
true_ma2 <- c(0.6,0.2)

# Performing BSL (reduce the number of iterations M if desired)
resultMa2BSL <- bsl(y = ma2$data, n = 500, M = 300000, model = model, covRandWalk = ma2$cov,
                     method = 'BSL', verbose = TRUE)
show(resultMa2BSL)
summary(resultMa2BSL)
plot(resultMa2BSL, which = 1, thetaTrue = true_ma2, thin = 20)

# Performing uBSL (reduce the number of iterations M if desired)
resultMa2uBSL <- bsl(y = ma2$data, n = 500, M = 300000, model = model, covRandWalk=ma2$cov,
                     method = 'uBSL', verbose = TRUE)
show(resultMa2uBSL)
summary(resultMa2uBSL)
plot(resultMa2uBSL, which = 1, thetaTrue = true_ma2, thin = 20)

# Performing tuning for BSLasso
ssy <- ma2_sum(ma2$data)
lambda_all <- list(exp(seq(-3,0.5,length.out=20)), exp(seq(-4,-0.5,length.out=20)),
                   exp(seq(-5.5,-1.5,length.out=20)), exp(seq(-7,-2,length.out=20)))
set.seed(100)
sp_ma2 <- selectPenalty(ssy = ssy, n = c(50, 150, 300, 500), lambda_all, theta = true_ma2, 
                        M = 100, sigma = 1.5, model = model, method = 'BSL', shrinkage = 'glasso')
sp_ma2
plot(sp_ma2)

# Performing BSLasso with a fixed penalty (reduce the number of iterations M if desired)
resultMa2BSLasso <- bsl(y = ma2$data, n = 300, M = 250000, model = model, covRandWalk=ma2$cov,
                        method = 'BSL', shrinkage = 'glasso', penalty = 0.027, verbose = TRUE)
show(resultMa2BSLasso)
summary(resultMa2BSLasso)
plot(resultMa2BSLasso, which = 1, thetaTrue = true_ma2, thin = 20)

# Performing semiBSL (reduce the number of iterations M if desired)
resultMa2SemiBSL <- bsl(y = ma2$data, n = 500, M = 300000, model = model, covRandWalk=ma2$cov,
                        method = 'semiBSL', verbose = TRUE)
show(resultMa2SemiBSL)
summary(resultMa2SemiBSL)
plot(resultMa2SemiBSL, which = 1, thetaTrue = true_ma2, thin = 20)

# Plotting the results together for comparison
# plot using the R default plot function
par(mar = c(5, 4, 1, 2), oma = c(0, 1, 2, 0))
combinePlotsBSL(list(resultMa2BSL, resultMa2uBSL, resultMa2BSLasso, resultMa2SemiBSL), which = 1,
                thetaTrue = true_ma2, thin = 20, label = c('bsl', 'uBSL', 'bslasso', 'semiBSL'),
                col = c('black', 'red', 'blue', 'green'), lty = 1:4, lwd = 1)
mtext('Approximate Univariate Posteriors', outer = TRUE, cex = 1.5)

# plot using the ggplot2 package
combinePlotsBSL(list(resultMa2BSL, resultMa2uBSL, resultMa2BSLasso, resultMa2SemiBSL), which = 2,
    thetaTrue = true_ma2, thin = 20, label = c('bsl', 'ubsl', 'bslasso', 'semiBSL'),
    options.color = list(values=c('black', 'red', 'blue', 'green')),
    options.linetype = list(values = 1:4), options.size = list(values = rep(1, 4)),
    options.theme = list(plot.margin = grid::unit(rep(0.03,4), 'npc'),
        axis.title = ggplot2::element_text(size=12), axis.text = ggplot2::element_text(size = 8),
        legend.text = ggplot2::element_text(size = 12)))
# }
# NOT RUN {
# }