| 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,
|
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
|
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\epsilonand\deltaparameters. 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 lengthTT -
sim_args: A list containingTT=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\times2 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)