The tsbugs package for R

My tsbugs package has gone up on CRAN. The functions in the tsbugs package are aimed to automate the writing of time series models to run in WinBUGS (Lunn et al., 2000) or OpenBUGS (Lunn et al., 2009). I created these functions a while back when I was doing some work on model averaging for time series models. I found it a lot easier to build R functions to write the BUGS models than the more error-inducing process of copy and pasting BUGS scripts, and then making slight alterations to create new models. It also allowed me to add arguments to specify different lag lengths, prior distributions, variance assumptions and data lengths. I decided not to write a vignette for CRAN, as it would have involved doing some estimation in BUGS via R2WinBUGS or R2OpenBUGS and running into some problems when submitted the package. Instead I thought I would post a quick guide here.

Autoregressive Models

The ar.bugs command builds a BUGS script for autoregressive (AR) models ready to use in R2OpenBUGS (Sturtz et al., 2005). For example, consider the LakeHuron data.

> LH <- LakeHuron 
> par(mfrow=c(2,1)) 
> plot(LH, ylab="Level in feet") 
> plot(diff(LH), ylab="Differenced Level")

We can construct a AR(1) model for this data (after differencing the data to obtain a stationary mean) as such:

> library("tsbugs")
> ar1 <- ar.bugs(y=diff(LH), ar.order=1)
> print(ar1)

model{
#likelihood
for(t in 2:97){
    y[t] ~ dnorm(y.mean[t], isigma2)
}
for(t in 2:97){
    y.mean[t] <- phi1*y[t-1]
}
#priors
phi1 ~ dnorm(0,1)
isigma2 ~ dgamma(0.000001,0.000001)
sigma <- pow(isigma2,-0.5)
}

The ar.bugs function allows for alternative specifications for prior distributions, forecasts and the inclusion of mean term:

> ar2 <- ar.bugs(y=diff(LH), ar.order=2, ar.prior="dunif(-1,1)", 
                 var.prior="dgamma(0.001,0.001)", k = 10, mean.centre = TRUE)
> print(ar2)

model{
#likelihood
for(t in 3:107){
    y[t] ~ dnorm(y.mean[t], isigma2)
}
for(t in 3:107){
    y.mean[t] <- phi0 + phi1*(y[t-1]-phi0) + phi2*(y[t-2]-phi0)
}
#priors
phi0 ~ dunif(-1,1)
phi1 ~ dunif(-1,1)
phi2 ~ dunif(-1,1)
sigma2 ~ dgamma(0.001,0.001)
isigma2 <- pow(sigma2,-1)
#forecast
for(t in 98:107){
    y.new[t] <- y[t]
}
}

The tsbugs objects can be used with R2OpenBUGS to easily run models from R. This is made even easier using the inits and nodes functions (also in the tsbugs package). For example:

> writeLines(ar2$bug, "ar2.txt") 
> library("R2OpenBUGS")
> ar2.bug <- bugs(data = ar2$data,
                  inits = list(inits(ar2)),
                  param = c(nodes(ar2, "prior")$name, "y.new"),
                  model = "ar2.txt", 
                  n.iter = 11000, n.burnin = 1000, n.chains = 1)

Note, 1) the model is written to a .txt file (as required by R2OpenBUGS), 2) the data used is part of the tsbugs object. The ar.bugs command cleans the data and adds missing values at the end of the series for foretasted values, 3) the initial values offered by the inits function are very crude, and with more complicated data or models, users might be better off specifying there own list of initial values.
The parameter traces and posterior distributions can be plotted using the coda package:

> library("coda")
> param.mcmc <- as.mcmc(ar2.bug$sims.matrix[,nodes(ar2, "prior")$name])
> plot(param.mcmc[,1:4])

The fanplot package can be used to plot the entire series of posterior predictive distributions. We may also plot (after deriving using the diffinv function) the posterior predictive distributions of the lake level:

> # derive future level
> ynew.mcmc <- ar2.bug$sims.list$y.new
> lhnew.mcmc <- apply(ynew.mcmc, 1, diffinv, xi = tail(LH,1))
> lhnew.mcmc <- t(lhnew.mcmc[-1,])
> # derive percentiles
> library("fanplot")
> k0 <- end(LH)[1]
> ynew.pn <- pn(ynew.mcmc, start=k0+1)
> lhnew.pn <- pn(lhnew.mcmc, start=k0+1)
> # plot
> par(mfrow=c(2,1)) 
> plot(diff(LH), type="n", xlim=k0+c(-50,10), ylab="Differenced Level")
> fan(ynew.pn)
> lines(diff(LH))
> plot(LH, type="n", xlim=k0+c(-50,10), ylab="Level")
> fan(lhnew.pn, cex=0.8)
> lines(LH)

Stochastic Volatility (SV) Models

The sv.bugs command builds a BUGS script for SV models ready to use in R2OpenBUGS. For example, consider the svpdx data.

> # plot
> plot(svpdx$pdx, type = "l", xaxt = "n", xlab = "Time", ylab = "Return")
> # x-axis
> svpdx$rdate <- format(svpdx$date, format = "%b %Y")
> mth <- unique(svpdx$rdate)
> qtr <- mth[seq(1,length(mth),3)]
> axis(1, at = match(qtr, svpdx$rdate), labels = qtr, cex.axis = 0.55)
> axis(1, at = match(mth, svpdx$rdate), labels = FALSE, tcl = -0.2)

We can construct a AR(0)-SV model for this data, and also obtain posterior simulations using the sv.bugs command:

> y <- svpdx$pdx
> sv0 <- sv.bugs(y, sim=TRUE)
> print(sv0)

model{
#likelihood
for(t in 1:945){
    y[t] ~ dnorm(y.mean[t], isigma2[t])
    isigma2[t] <- exp(-h[t])
    h[t] ~ dnorm(h.mean[t], itau2)
}
for(t in 1:945){
    y.mean[t] <- 0
}
for(t in 1:1){
    h.mean[t] <- psi0
}
for(t in 2:945){
    h.mean[t] <- psi0 + psi1*(h[t-1]-psi0)
}
#priors
psi0 ~ dnorm(0,0.001)
psi1.star ~ dunif(0,1)
psi1 <- 2*psi1.star-1
itau2 ~ dgamma(0.01,0.01)
tau <- pow(itau2,-0.5)
#simulation
for(t in 1:945){
    y.mean.c[t] <- cut(y.mean[t])
    isigma2.c[t] <- cut(isigma2[t])
    y.sim[t] ~ dnorm(y.mean.c[t],isigma2.c[t])
}
}

This model closely matches those presented in Meyer and Yu (2002). There are further options in the tsbugs package to incorporate different priors that do not involve transformations such as those for psi1 above. Using R2OpenBUGS we can fit the model,

> # decent initial value for variance in first period
> init <- inits(sv0, warn=FALSE)
> init$psi0 <- log(var(y))
> # write bug
> writeLines(sv0$bug, "sv0.txt")
> # might take a while to compile
> sv0.bug <- bugs(data = sv0$data,
                  inits = list(init),
                  param = c(nodes(sv0, "prior")$name,"y.sim","h"),
                  model = "sv0.txt", 
                  n.iter = 11000, n.burnin = 1000, n.chains = 1)

The volatility and time-dependent standard deviations estimates can be easily derived,

> h.mcmc <- sv0.bug$sims.list$h
> sigma.mcmc <- sqrt(exp(h.mcmc))

Which allows us to directly view the estimated volatility process or the time-dependent standard deviation using the fanplot package ,

> # derive percentile objects
> h.pn <- pn(h.mcmc)
> sigma.pn <- pn(sigma.mcmc)
> # plot
> par(mfrow=c(2,1), mar=c(2,4,2,1)) 
> plot(NULL, type = "l", xlim = c(1, 945), xlab = "Time", xaxt = "n", 
             ylim = range(h.pn[,5:95]), ylab="Volatility") 
> # add axis
> axis(1, at = match(qtr, svpdx$rdate), labels = qtr, cex.axis = 0.55) 
> axis(1, at = match(mth, svpdx$rdate), labels = FALSE, tcl = -0.2) 
> # fan
> fan(h.pn, ln=c(1,10,50,90,99)) 
> fan.txt(h.pn, pn.l = c(1,10,50,90,99)) 
> plot(NULL, type = "l", xlim = c(1, 945), xlab = "Time", xaxt = "n", 
             ylim = range(sigma.pn[,1:95]), ylab="Standard Deviation") 
> axis(1, at = match(qtr, svpdx$rdate), labels = qtr, cex.axis = 0.55) 
> axis(1, at = match(mth, svpdx$rdate), labels = FALSE, tcl = -0.2) 
> fan(sigma.pn, ln=c(1,10,50,90,99)) 
> fan.txt(sigma.pn, pn.l = c(1,10,50,90,99))

We can also plot the posterior simulations from the model:

> # derive percentiles
> y.mcmc <- sv0.bug$sims.list$y.sim
> y.pn <- pn(y.mcmc)
> # plot
> plot(NULL, type = "l", xlim = c(1, 945), xlab = "Time", xaxt = "n", 
             ylim = range(y.pn[,5:95]), ylab="Return")
> axis(1, at = match(qtr, svpdx$rdate), labels = qtr, cex.axis = 0.55)
> axis(1, at = match(mth, svpdx$rdate), labels = FALSE, tcl = -0.2)
> fan(y.pn)
> lines(y)

Random Variance (RV) Shift Models

The rv.bugs command builds a BUGS script for random variance shift models, similar to that of McCulloch and Tsay (1993 ready to use in R2OpenBUGS. Consider the ew data.

> r <- ts(ew[2:167]/ew[1:166]-1, start=1841)
> y <- diff(r) 
> plot(y, ylab="Difference in Population Growth Rate")

We can create a BUGS script to fit a RV model to this data, including posterior simulations, using the rv.bugs command:

> rv0 <- rv.bugs(y, sim=TRUE)
> print(rv0)

model{
#likelihood
for(t in 1:165){
    y[t] ~ dnorm(y.mean[t], isigma2[t])
    isigma2[t] <- exp(-h[t])
    h[t] <- 2*lsig[t]
}
for(t in 1:165){
    y.mean[t] <- 0
}
lsig[1] <- -0.5*log(isig02)
for(t in 2:165){
    lsig[t] <- lsig[t-1]+(delta[t]*beta[t])
    delta[t] ~ dbern(epsilon)
    beta[t] ~ dnorm(0,ilambda2)
}
#priors
isig02 ~ dgamma(0.000001,0.000001)
sig0 <- pow(isig02,-0.5)
epsilon ~ dbeta(1, 100)
ilambda2 ~ dgamma(0.01,0.01)
lambda <- pow(ilambda2,-0.5)
#simulation
for(t in 1:165){
    y.mean.c[t] <- cut(y.mean[t])
    isigma2.c[t] <- cut(isigma2[t])
    y.sim[t] ~ dnorm(y.mean.c[t],isigma2.c[t])
}
}

and then run the script in R2OpenBUGS:

> # decent inital value for variance in first period
> init <- inits(rv0, warn=FALSE)
> init$isig02<-sd(y)^-2
> # write bug
> writeLines(rv0$bug,"rv0.txt")
> # might take a while to compile
> rv0.bug <- bugs(data = rv0$data,
                  inits = list(init),
                  param = c(nodes(rv0, "prior")$name, "y.sim", "h", "delta", "beta"),
                  model = "rv0.txt", 
                  n.iter = 11000, n.burnin = 1000, n.chains = 1)

We can plot the posterior simulations from the model using the fanplot package:

> # derive percentiles
> y0 <- tsp(y)[1]
> y.mcmc <- rv0.bug$sims.list$y.sim
> y.pn <- pn(y.mcmc, start = y0)
> # plot
> plot(NULL, type = "n", xlim=tsp(y.pn)[1:2], ylim = range(y.pn[,5:95]), ylab="Difference in Population Growth Rate", xlab="Time")
> fan(y.pn, txt=NA)
> fan.txt(y.pn, pn.r = seq(10,90,40))
> fan.txt(y.pn, pn.l = seq(10,90,10))
> lines(y)

Alongside the posterior distributions of the volatility and standard deviations,

> # derive percentiles
> h.mcmc <- rv0.bug$sims.list$h
> h.pn <- pn(h.mcmc, start = y0)
> sigma.pn <- pn(sims = sqrt(exp(h.mcmc)), start = y0)
> # plots
> par(mfrow=c(2,1))
> plot(NULL, type = "n&quot;, xlim =tsp(h.pn)[1:2], ylim = range(h.pn[,5:95]), ylab=&quot;Volatility&quot;, xlab=&quot;Time&quot;)
> fan(h.pn, ln=c(1,10,50,90,99))
> fan.txt(h.pn, pn.l = c(1,10,50,90,99))
> plot(NULL, type = "n", xlim =tsp(sigma.pn)[1:2], ylim = range(sigma.pn[,1:95]), ylab="Standard Deviation", xlab="Time")
> fan(sigma.pn, ln=c(1,10,50,90,99))
> fan.txt(sigma.pn, pn.l = c(1,50,99))

The posterior distributions of the multiplier effect of the shift in variance (beta[t] in the BUGS model) can also be plotted. Note, when there is no variance shift, the posterior of the beta[t] is similar to the prior distribution.

> # derive percentiles
> beta.mcmc <- rv0.bug$sims.list$beta
> beta.pn <- pn(beta.mcmc, start = y0)
> # plots
> plot(NULL, type = "n", xlim =tsp(beta.pn)[1:2], ylim = range(beta.pn), ylab="Variance Multiplier", xlab="Time")
> fan(beta.pn)
> fan.txt(beta.pn, pn.l = seq(10,90,10))

Further Work

I hope to develop the package in the future with a few more functions to produce some MA and GARCH models. I am also thinking about ways to incorporate non-normal time series data.

References

Lunn, D. J., A. Thomas, N. Best, and D. Spiegelhalter (2000). WinBUGS – A Bayesian modelling framework: Concepts, structure, and extensibility. Statistics and Computing 10 (4), 325–337.
Lunn, D. et al. (2009). The BUGS project: Evolution, critique and future directions. Statistics in Medicine 28.25, 3049-3067
Meyer, R. and J. Yu (2002). BUGS for a Bayesian analysis of stochastic volatility models. Econometrics Journal, 3, 2, pp. 198–215.
McCulloch, R. E. and R. S. Tsay (1993). Bayesian Inference and Prediction for Mean and Variance Shifts in Autoregressive Time Series. Journal of the American Statistical Association 88 (423), 968-978.
Sturtz, S., U. Ligges, and A. Gelman (2005). R2WinBUGS: a package for running WinBUGS from R. Journal of Statistical Software 12 (3).

About these ads