ma2: An MA(2) model

In BSL: Bayesian Synthetic Likelihood

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 \insertCiteAn2019;textualBSL. The journal article \insertCiteAn2022;textualBSL 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, θ_1 and θ_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 + θ_1 z_{t-1} + θ_2 z_{t-2},

where t=1,…,TT and z_t \sim N(0,1) for t=-1,0,…,TT. A uniform prior is used for this example, subject to the restrictions that -2<θ_1<2, θ_1+θ_2>-1 and θ_1-θ_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 ε and δ parameters. This is the sinh-arcsinnh transformation of \insertCiteJones2009;textualBSL. 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 θ_1 = 0.6, θ_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

\insertAllCited

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)

BSL documentation built on Nov. 3, 2022, 9:06 a.m.