Simulation of AR(1)-process with \( \phi \) = 0.9 (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 use ggplot2.

The Existing Example

As it is the example demonstrates how to simulate an AR(1)-process. Additionally the code runs an autocorrelation test and partial autocorrelation test.

  set.seed(123456)
  y <- arima.sim(n = 100, list(ar = 0.9), innov=rnorm(100))
  op <- par(no.readonly=TRUE)
  layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE))
  plot.ts(y, ylab='')
  acf(y, main='Autocorrelations', ylab='',
      ylim=c(-1, 1), ci.col = "black")
  pacf(y, main='Partial Autocorrelations', ylab='',
       ylim=c(-1, 1), ci.col = "black")

plot of chunk basicCode

  par(op)

Converting To `ggplot2`

The Time Series Plot

ggplot2 prefers data in a data.frame. So let's change our ts object to a data.frame

  y.df <- data.frame(time=c(1:length(y)),value=y)

Now let's do the time series plot

  require(ggplot2)

## Loading required package: ggplot2

  timeSeriesPlot <- ggplot(y.df,aes(x=time,y=value)) + geom_line()

The Autocorrelation Plot

This will require a bit more work. Notice there are two parts to the autocorrelation plot. First, the autocorrelations are plotted as vertical lines for each lag. Second, the dashed lines indicate a 95% confidence interval. To get started we will store the results of the autocorrelation function.

  y.acf<-acf(y, plot=FALSE)
  y.pacf<-pacf(y, plot=FALSE)

If you look here you will see that the confidence intervals are calculated using the following code.

  ci <- 0.95 # Indicates 95% confidence level
  clim0 <- qnorm((1 + ci)/2)/sqrt(y.acf$n.used)
  clim <- c(-clim0,clim0)

Again we will need to use data frames to make the plot. Note that we will plot three horizontal lines:

Lower Confidence Level (dashed and blue)
Upper Confidence Level (dashed and blue)
Horizontal Axis (solid and black)

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

  acfPlot <- ggplot(data.frame(lag=y.acf$lag,acf=y.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","blue")) +
    scale_linetype_manual(values =c("solid","dashed")) +
    opts(title="Autocorrelations")

The partial autocorrelation plot is a very similar process.




  pacfPlot <- ggplot(data.frame(lag=y.pacf$lag,pacf=y.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","blue")) +
    scale_linetype_manual(values =c("solid","dashed")) +
    opts(title="Partial Autocorrelations")

To put it together we have

  require(gridExtra)

## Loading required package: gridExtra

## Loading required package: grid

  timeSeriesPlot

plot of chunk allPlots


  grid.arrange(acfPlot,pacfPlot,ncol=2)