Simulation of MA(1)-process with \( \theta = 0.8 \) (In ggplot2)

Setup

I am working my way through Bernhard Pfaff's Analysis of Integrated and Cointegrated Time Series with R. The book contains some excellent examples. In this RPub I'll show you how to change the code of Example 1.1 to produce Figure 1.6. I'll also show you how to do it in ggplot2.

You can see RPubs similar to this one for the following examples:
1.1
1.2

Quickly Plotting Time Series

I really like the general layout that Pfaff uses of displaying a time series and its auto- and partial auto-correlation plots together in one figure. To make this faster I have written the following custom function.

quickTSPlots<-function(TS, ylab = '', ylim=c(-1,1), ...)
{
  op <- par(no.readonly=TRUE) # We need to remember the original
  #       parameters so we can return the settings back as they were

  layout(matrix(c(1,1,2,3), 2, 2, byrow=TRUE)) # This gives the nice
  #       layout of the time series on top with the correlation plots
  #       below it

  plot(TS,ylab=ylab, ...)
  acf(TS, main='Autocorrelations', ylab=ylab,
      ylim=ylim, ci.col = "black")
  pacf(TS, main='Partial Autocorrelations', ylab=ylab,
       ylim=ylim, ci.col = "black")

  par(op) # Return the settings as they were before


}

Simulating and Plotting the MA(1)-process

We can very easily adapt the code from Example 1.1 to accomadate an MA(1)-process instead of an AR(1)-process.

set.seed(123456)
y.sim <- arima.sim(n = 100, list(ma=0.8), innov=rnorm(100))

quickTSPlots(y.sim)

plot of chunk MAprocess

Interestingly I find that there is always more than one significant partial autocorrelation lag when I run this simulation. In Pfaff's book the only significant partial autocorrelation lag is the first. At this point in time my knowledge of time series analysis is not deep enough to explain this.

Putting it in ggplot2

Start by making sure you have ggplot2 and gridExtra loaded.

require(ggplot2)
## Loading required package: ggplot2
require(gridExtra)
## Loading required package: gridExtra
## Loading required package: grid

Now let's make a new version of our quickTSPlots function that will accomodate ggplot2.

quickTSPlots<-function(TS, ylab = '', ylim=c(-1,1),ggplot2=FALSE, ...)
{
  if (ggplot2) {
    TS.df <- data.frame(time=c(1:length(TS)),value=TS)
    timeSeriesPlot <- ggplot(TS.df,aes(x=time,y=value)) + geom_line()
    TS.acf<-acf(TS, plot=FALSE)
    TS.pacf<-pacf(TS, plot=FALSE)
    ci <- 0.95 # Indicates 95% confidence level
    clim0 <- qnorm((1 + ci)/2)/sqrt(TS.acf$n.used)
    clim <- c(-clim0,clim0)

    hline.data <- data.frame(z=c(0,clim),
                       type=c("base","ci","ci"))

    acfPlot <- ggplot(data.frame(lag=TS.acf$lag,acf=TS.acf$acf)) +
    geom_hline(aes(yintercept=z,colour=type,linetype=type),hline.data) +
    geom_linerange(aes(x=lag,ymin=0,ymax=acf)) +
    scale_colour_manual(values = c("black","black")) +
    scale_linetype_manual(values =c("solid","dashed")) +
    opts(title="Autocorrelations")

    pacfPlot <- ggplot(data.frame(lag=TS.pacf$lag,pacf=TS.pacf$acf)) +
    geom_hline(aes(yintercept=z,colour=type,linetype=type),hline.data) +
    geom_linerange(aes(x=lag,ymin=0,ymax=pacf)) +
    scale_colour_manual(values = c("black","black")) +
    scale_linetype_manual(values =c("solid","dashed")) +
    opts(title="Partial Autocorrelations")


    grid.arrange(timeSeriesPlot, arrangeGrob(acfPlot, pacfPlot, ncol=2),
                 ncol=1) 




  } else {
    op <- par(no.readonly=TRUE) # We need to remember the original
    #       parameters so we can return the settings back as they were

    layout(matrix(c(1,1,2,3), 2, 2, byrow=TRUE)) # This gives the nice
    #       layout of the time series on top with the correlation plots
    #       below it

    plot(TS,ylab=ylab, ...)
    acf(TS, main='Autocorrelations', ylab=ylab,
        ylim=ylim, ci.col = "black")
    pacf(TS, main='Partial Autocorrelations', ylab=ylab,
         ylim=ylim, ci.col = "black")

    par(op) # Return the settings as they were before
  }

}

If the ggplot2 portion of the code above is a bit difficult then take a peak at my RPub for Example 1.1. I came across the arrangeGrob function above in this helpful question at stackoverflow.com

Now let's test the new function out on our MA(1)-process.

  quickTSPlots(y.sim, ggplot2=TRUE)

plot of chunk ggMA


So that's that! If you enjoyed this then please check out my other RPubs. Additionally take a peak over at my blog at http://www.quantifyingQuality.com