ma2 {BSL} R Documentation

## 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). The journal article An et al. (2022) provides a full description of how to use this package for the toad example.

### Usage

data(ma2)

ma2_sim(theta, TT)

ma2_sim_vec(n, theta, TT)

ma2_sum(x, epsilon = 0, delta = 1)

ma2_prior(theta)


### Arguments

 theta A vector of proposed model parameters, \theta_1 and \theta_2. TT The number of observations. n The number of simulations to run with the vectorised simulation function. x Observed or simulated data in the format of a vector of length TT. epsilon The skewness parameter in the sinh-arcsinh transformation. delta The kurtosis parameter in the sinh-arcsinh transformation.

### 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,TT and z_t \sim N(0,1) for t=-1,0,\ldots,TT. 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.

### Functions

• ma2_sim: Simulates an MA(2) time series.

• ma2_sim_vec: Simulates n MA(2) time series with a vectorised simulation function.

• ma2_sum: Returns the summary statistics for a given data set. The skewness and kurtosis of the summary statistics can be controlled via the \epsilon and \delta parameters. This is the sinh-arcsinnh transformation of Jones and Pewsey (2009). By default, the summary statistics function simply returns the raw data. Otherwise, the transformation is introduced to motivate the “semiBSL” method.

• ma2_prior: Evaluates the (unnormalised) log prior, which is uniform subject to several restrictions related to invertibility of the time series.

### 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 TT=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 TT

• sim_args: A list containing TT=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 \times 2 matrix

### Author(s)

Ziwen An, Leah F. South and Christopher Drovandi

### References

An Z, South LF, Drovandi CC (2022). “BSL: An R Package for Efficient Parameter Estimation for Simulation-Based Models via Bayesian Synthetic Likelihood.” Journal of Statistical Software, 101(11), 1–33. doi: 10.18637/jss.v101.i11.

An Z, South LF, Nott DJ, Drovandi CC (2019). “Accelerating Bayesian Synthetic Likelihood With the Graphical Lasso.” Journal of Computational and Graphical Statistics, 28(2), 471–475. doi: 10.1080/10618600.2018.1537928.

Jones MC, Pewsey A (2009). “Sinh-arcsinh distributions.” Biometrika, 96(4), 761–780. ISSN 0006-3444.

### Examples

## Not run:
# Load the data for this example and set up the model object
data(ma2)
model <- newModel(fnSimVec = ma2_sim_vec, fnSum = ma2_sum, simArgs = ma2$sim_args, theta0 = ma2$start, fnLogPrior = ma2_prior)
thetaExact <- c(0.6, 0.2)

# reduce the number of iterations M if desired for all methods below
# Method 1: standard BSL
resultMa2BSL <- bsl(y = ma2$data, n = 500, M = 300000, model = model, covRandWalk = ma2$cov,
method = "BSL", verbose = 1L)
show(resultMa2BSL)
summary(resultMa2BSL)
plot(resultMa2BSL, thetaTrue = thetaExact, thin = 20)

# Method 2: unbiased BSL
resultMa2uBSL <- bsl(y = ma2$data, n = 500, M = 300000, model = model, covRandWalk=ma2$cov,
method = "uBSL", verbose = 1L)
show(resultMa2uBSL)
summary(resultMa2uBSL)
plot(resultMa2uBSL, thetaTrue = thetaExact, thin = 20)

# Method 3: BSLasso (BSL with glasso shrinkage estimation)
# tune the penalty parameter fisrt
ssy <- ma2_sum(ma2$data) lambdaAll <- list(exp(seq(-5.5,-1.5,length.out=20))) set.seed(100) penaltyGlasso <- selectPenalty(ssy = ssy, n = 300, lambdaAll, theta = thetaExact, M = 100, sigma = 1.5, model = model, method = "BSL", shrinkage = "glasso") penaltyGlasso plot(penaltyGlasso) resultMa2BSLasso <- bsl(y = ma2$data, n = 300, M = 250000, model = model, covRandWalk=ma2$cov, method = "BSL", shrinkage = "glasso", penalty = 0.027, verbose = 1L) show(resultMa2BSLasso) summary(resultMa2BSLasso) plot(resultMa2BSLasso, thetaTrue = thetaExact, thin = 20) # Method 4: BSL with Warton's shrinkage and Whitening # estimate the Whtieing matrix and tune the penalty parameter first W <- estimateWhiteningMatrix(20000, model, method = "PCA", thetaPoint = ma2$start)
gammaAll <- list(seq(0.3, 0.8, 0.02))
set.seed(100)
penaltyWarton <- selectPenalty(ssy = ssy, n = 300, gammaAll, theta = thetaExact,
M = 100, sigma = 1.2, model = model, method = "BSL", shrinkage = "Warton",
whitening = W)
penaltyWarton
plot(penaltyWarton, logscale = FALSE)

resultMa2Whitening <- bsl(y = ma2$data, n = 300, M = 250000, model = model, covRandWalk=ma2$cov,
method = "BSL", shrinkage = "Warton", whitening = W,
penalty = 0.52, verbose = 1L)
show(resultMa2Whitening)
summary(resultMa2Whitening)
plot(resultMa2Whitening, thetaTrue = thetaExact, thin = 20)

# Method 5: semiBSL, the summary statistics function is different from previous methods
model2 <- newModel(fnSimVec = ma2_sim_vec, fnSum = ma2_sum, simArgs = ma2$sim_args, sumArgs = list(epsilon = 2), theta0 = ma2$start, fnLogPrior = ma2_prior)
sim <- simulation(model, n = 1e4, theta = ma2$start, seed = 1) # run a short simulation plot(density(sim$ssx[, 1])) # the first marginal summary statistic is right-skewed
resultMa2SemiBSL <- bsl(y = ma2$data, n = 500, M = 200000, model = model2, covRandWalk=ma2$cov,
method = "semiBSL", verbose = 1L)
show(resultMa2SemiBSL)
summary(resultMa2SemiBSL)
plot(resultMa2SemiBSL, thetaTrue = thetaExact, thin = 20)

# Method 6: BSL with consideration of model misspecification (mean adjustment)
resultMa2Mean <- bsl(y = ma2$data, n = 500, M = 200000, model = model, covRandWalk=ma2$cov,
method = "BSLmisspec", misspecType = "mean", verbose = 1L)
show(resultMa2Mean)
summary(resultMa2Mean)
plot(resultMa2Mean, thetaTrue = thetaExact, thin = 20)

# Method 7: BSL with consideration of model misspecification (variance inflation)
resultMa2Variance <- bsl(y = ma2$data, n = 500, M = 200000, model = model, covRandWalk=ma2$cov,
method = "BSLmisspec", misspecType = "variance", verbose = 1L)
show(resultMa2Variance)
summary(resultMa2Variance)
plot(resultMa2Variance, thetaTrue = thetaExact, thin = 20)

# Plotting the results together for comparison
# plot using the R default plot function
oldpar <- par()
par(mar = c(5, 4, 1, 2), oma = c(0, 1, 2, 0))
combinePlotsBSL(list(resultMa2BSL, resultMa2uBSL, resultMa2BSLasso, resultMa2SemiBSL), which = 1,
thetaTrue = thetaExact, 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 = thetaExact, 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)))
par(mar = oldpar$mar, oma = oldpar$oma)

## End(Not run)



[Package BSL version 3.2.4 Index]