The tsbugs package for R

I have/will update this post as I expanded the tsbugs package.

My tsbugs package has gone up on CRAN. I decided not to write a vignette for the submission, 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 below. Note, the package code is on my github if you would like to contribute more models to the package…

The functions in the tsbugs package are aimed to automate the writing of time series models to run in WinBUGS or OpenBUGS. 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. Below are examples for three types of time series models; autorgressive models with constant variance, stochastic volatility and random variance shift models.

Autoregressive Models

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

LH <- LakeHuron 
par(mfrow=c(2,1)) 
plot(LH, main="Level (ft)") 
plot(diff(LH), main="Differenced Level")

tsbugs1-cvdata
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)

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)

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])

tsbugs1-cvparam
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,])

# plot differenced
par(mfrow=c(2,1)) 
plot(diff(LH) ,xlim=k0+c(-50,10), main="Differenced Level")

# add fan
library("fanplot")
k0 <- end(LH)[1]
fan(ynew.mcmc, start=k0+1, rcex=0.5)

# plot undifferenced
plot(LH, xlim=k0+c(-50,10), main="Level")
fan(lhnew.mcmc, start=k0+1, rcex=0.5)

tsbugs1-cvforc

Stochastic Volatility Models

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

# plot
plot(svpdx$pdx, type = "l", 
     main = "Return of Pound-Dollar exchange rate data 
             from 2nd October 1981 to 28th June 1985", 
     cex.main = 0.8)

tsbugs1-svdata
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)

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 estimates can be easily extracted,

h.mcmc <- sv0.bug$sims.list$h

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

# plot
plot(NULL, xlim = c(1, 945)+c(0,40), ylim = c(-4,2), 
     main="Estimated Volatility from SV Model") 
# fan
fan(h.mcmc, type = "interval")

tsbugs1-svh
We can also plot the posterior simulations from the model:

# derive percentiles
y.mcmc <- sv0.bug$sims.list$y.sim
# plot
plot(NULL, type = "l", xlim = c(1, 945)+c(0,20), 
     ylim = range(y), 
     main = "Posterior Model Simulations and Data")
fan(y.mcmc)
lines(y)

tsbugs1-svysim

Random Variance Shift Models

The rv.bugs command builds a BUGS script for random variance (RV) 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, main="Difference in England and Wales 
              Population Growth Rate")

tsbugs1-rvdata
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)

and then run the script in R2OpenBUGS (this can take a couple of hours):

# 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

# plot
plot(NULL, xlim=tsp(y)[1:2]+c(-5,5), ylim = range(y), 
     main="Posterior Simulations")
fan(y.mcmc, start = y0, rlab=c(10,50,90), llab=TRUE)
lines(y)

tsbugs1-rvysim
Alongside the posterior distributions of the standard deviations,

# derive sigma
h.mcmc <- rv0.bug$sims.list$h
sigma.mcmc <- sqrt(exp(h.mcmc))
# plots
plot(NULL, xlim =tsp(y)[1:2]+c(-5,5), ylim = c(0,0.008), 
     main="Standard Deviation")
fan(sigma.mcmc, start = y0, rlab=c(5,50,95), llab = c(5,50,95))

tsbugs1-rvhsigma
The posterior distributions of the probability of a variance shift and multiplier effect of the shift in variance (delta[t] and 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.

#extract data
delta.mcmc <- rv0.bug$sims.list$delta
beta.mcmc <- rv0.bug$sims.list$beta

# plots
par(mfrow=c(2,1))
plot(NULL, xlim = tsp(y)[1:2]+c(-5,5), ylim = c(0,1),
     main="Probability of Variance Change Point")
fan(delta.mcmc, start=y0, ln = NULL, rlab = NULL)

plot(NULL, xlim = tsp(y)[1:2]+c(-5,5), ylim = c(-2,2), 
     main="Variance Multiplier")
fan(beta.mcmc, start=y0)

tsbugs1-rvbeta

About these ads
This entry was posted in BUGS, R, Research, tsbugs and tagged , , , , , , , , , , , , , . Bookmark the permalink.

11 Responses to The tsbugs package for R

  1. rtkapend says:

    Speechless! Anyway, well done mate.

    Happy new year Guy!

    Any plan to visit Southampton soon? Maybe then I can have a chance to be talked through some of the clever stuff you don’t stop doing. I’m getting closer to conduct my sensitivity analysis so expect more bugs from this end.

    Check you later,

    Rich

    ________________________________ De : Guy Abel À : rkapend@yahoo.fr Envoyé le : Mardi 15 janvier 2013 16h26 Objet : [New post] tsbugs Package

    WordPress.com gjabel posted: “My tsbugs package has gone up on CRAN. The functions are aimed to automate the writing of time series models to run in WinBUGS or OpenBUGS. I created these functions a while back when I was doing some work on model averaging for time series models. I foun”

  2. Pablo E. Verde says:

    I expected that the code for prediction to be:
    #forecasts
    for(t in 98:107){
    y.new[t] ~ dnorm(y.mean[t], isigma2)
    }
    and not the outcome code:

    #forecasts
    for(t in 98:107){
    y.new[t] <- y[t]
    }
    Is this a bug or I missed something?
    Thanks, great tool!

    Pablo

    • I think that the likelihood for the forecast period (98-107) already has been made by this part:
      #likelihood
      for(t in 3:107){
      y[t] ~ dnorm(y.mean[t], isigma2)
      }

      Then the y values for the forecast period (98-107) are copied into a new variable y.new just to make it ieasier to plot the forecast separately using fan(ynew.pn).

    • gjabel says:

      Thanks Pablo.

      I think our code for the forecasts are broadly equivalent. Both of our stochastic y.new[t] is generated from the same time dependent mean and variance as my y[t]. However, in the BUGS model I have only set up random nodes for future values once (in the first for loop of the BUGS script, where WinBUGS or OpenBUGS will see the NA’s in the data as missing) and then relabeled them in the forecast chunk of the BUGS code. This works, as using the ar.bugs function the data is modified to add on the right amount of NA’s at the end of the series so BUGS treats them as random nodes (in the top loop). You can see the data (that is used in the data argument on the R2OpenBUGS bugs command) by entering ar2$data into the R console. If I did not modify the data, then I would need a new loop for forecasted values (like the one you propose). Hope this makes sense?

  3. Thanks a lot for making this package! This certainly lowers the threshold for using MCMC modelling for a lot of users.

  4. Mudassar chanda says:

    This is very useful package and less time consuming. I am interested in non-normal time series models like Poisson. Is there any plan to extend this package to handle such data.
    Thanks

  5. clarkh85 says:

    Awesome package! Is there an extension to the package to perform multivariate SV modeling? i.e., similar to the method described in Meyer & Yu (2006). Thank you!

  6. gjabel says:

    Not at the moment. I have the code to run the models, but I have not incorporated them into the package. Just slapped all the package files up on github if you want to contribute… https://github.com/gjabel/tsbugs. The multivariate models are going to be in fn_tsbugs4.R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s