Simulated AR(2) time series

Below I will simulate an AR(2) time series using random normally-distributed data for the error term.

n <- 500
epsilon <- rnorm(n = n,mean = 0,sd = 1) # white noise
x <- numeric(length = n)
phi <- 0.5
phi2 <- 0.2
for(i in 3:n){
  x[i] <- phi * x[i-1] + phi2 * x[i-2] + epsilon[i] }
plot(x,type="l")

The above figure is the AR(2) time series. Below I show the acf and pacf plots.

acf(x)

pacf(x)

The above figures look what what I expected them to. The acf plot shows significant autocorrelation out to 6 lags, but the partial acf plot shows that it is truly an AR(2) process, and that the significant autocorrelation to further lags shown in the acf plot is a result of the the propagation of autocorrelation from the first two lags. The plot shows correlations that are close to my phi values. Below I will try using different values for phi and phi2, and compare the acf figures.

phi <- 0.9
phi2 <- -0.4
for(i in 3:n){
  x[i] <- phi * x[i-1] + phi2 * x[i-2] + epsilon[i] }
acf(x)

pacf(x)

Again, the above plots are what I expected. The partial acf plot shows significant AR(2) process with correlations roughly what I set as phi and phi2.

Tree-ring data

Below I show code that I used for my tree-ring data. More information about these processes can be found in last week’s rpub note. I call dplR, read in my ring width data, attempt to detrend the data by fitting a negative exponential curve, and create a mean value chronology. I then create acf and pacf plots for my data.

library(dplR)
## Warning: package 'dplR' was built under R version 3.1.2
mwa <- read.rwl('MWA_308series.rwl',long=TRUE)
mwa.rwi <- detrend(rwl = mwa, method = "ModNegExp")
mwa.crn <- chron(mwa.rwi, prefix = "MWA")
acf(mwa.crn$MWAstd)

pacf(mwa.crn$MWAstd)

The partial acf plot shows significant autocorrelation out to lag 6 for the bristlecone pine data. This makes sense to me biologically. The trees can hold their needles for years and photosynthates can become stored up. Because of this, growth in one year is dependent on growth in previous years. I then fit an AR model to the data to see which order of model would fit the data best. The AR function uses AIC, which is a criterion to create an accurate but hopefully simplistic model.

mwa.ar <- ar(mwa.crn$MWAstd)
mwa.ar$order
## [1] 8
plot(mwa.ar$aic)

mwa.ar <- ar(mwa.crn$MWAstd,order.max = 5)

The AR function fit an AR(8) model, but after looking at a plot of the AIC values, it seems like an AR(5) model would be sufficient, and a simpler model. Below, I show the values for phi1-phi5.

mwa.ar$ar
## [1] 0.29363605 0.19502136 0.11006081 0.09470540 0.05501646