5.2 Stationarity
It is important to test and transform (via differencing) your data to ensure stationarity when fitting an ARMA model using standard algorithms. The standard algorithms for ARIMA models assume stationarity and we will be using those algorithms. It possible to fit ARMA models without transforming the data. We will cover that in later chapters. However, that is not commonly done in the literature on forecasting with ARMA models, certainly not in the literature on catch forecasting.
Keep in mind also that many ARMA models are stationary and you do not want to get in the situation of trying to fit an incompatible process model to your data. We will see examples of this when we start fitting models to non-stationary data and random walks.
5.2.1 Look at stationarity in simulated data
We will start by looking at white noise and a stationary AR(1) process from simulated data. White noise is simply a string of random numbers drawn from a Normal distribution. rnorm()
with return random numbers drawn from a Normal distribution. Use ?rnorm
to understand what the function requires.
<- 100
TT <- rnorm(TT, mean = 0, sd = 1) # 100 random numbers
y <- par(mfrow = c(1, 2))
op plot(y, type = "l")
acf(y)
par(op)
Here we use ggplot()
to plot 10 white noise time series.
<- data.frame(t = 1:TT, y = y)
dat <- ggplot(dat, aes(x = t, y = y)) + geom_line() + ggtitle("1 white noise time series") +
p1 xlab("") + ylab("value")
<- matrix(rnorm(TT * 10), TT, 10)
ys <- data.frame(ys)
ys $id = 1:TT
ys
<- melt(ys, id.var = "id")
ys2 <- ggplot(ys2, aes(x = id, y = value, group = variable)) +
p2 geom_line() + xlab("") + ylab("value") + ggtitle("10 white noise processes")
grid.arrange(p1, p2, ncol = 1)
These are stationary because the variance and mean (level) does not change with time.
An AR(1) process is also stationary.
<- 0.8
theta <- 10
nsim <- arima.sim(TT, model = list(ar = theta))
ar1 plot(ar1)
We can use ggplot to plot 10 AR(1) time series, but we need to change the data to a data frame.
<- data.frame(t = 1:TT, y = ar1)
dat <- ggplot(dat, aes(x = t, y = y)) + geom_line() + ggtitle("AR-1") +
p1 xlab("") + ylab("value")
<- matrix(0, TT, nsim)
ys for (i in 1:nsim) ys[, i] <- as.vector(arima.sim(TT, model = list(ar = theta)))
<- data.frame(ys)
ys $id <- 1:TT
ys
<- melt(ys, id.var = "id")
ys2 <- ggplot(ys2, aes(x = id, y = value, group = variable)) +
p2 geom_line() + xlab("") + ylab("value") + ggtitle("The variance of an AR-1 process is steady")
grid.arrange(p1, p2, ncol = 1)
Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
5.2.2 Stationary around a linear trend
Fluctuating around a linear trend is a very common type of stationarity used in ARMA modeling and forecasting. This is just a stationary process, like white noise or AR(1), around an linear trend up or down.
<- 0.5
intercept <- 0.1
trend <- 0.5
sd <- 20
TT <- rnorm(TT, sd = sd) #white noise
wn <- wn + intercept #white noise witn interept
wni <- wn + trend * (1:TT) + intercept wnti
See how the white noise with trend is just the white noise overlaid on a linear trend.
<- par(mfrow = c(1, 3))
op plot(wn, type = "l")
plot(trend * 1:TT)
plot(wnti, type = "l")
par(op)
We can make a similar plot with ggplot.
<- data.frame(t = 1:TT, wn = wn, wni = wni, wnti = wnti)
dat <- ggplot(dat, aes(x = t, y = wn)) + geom_line() + ggtitle("White noise")
p1 <- ggplot(dat, aes(x = t, y = wni)) + geom_line() + ggtitle("with non-zero mean")
p2 <- ggplot(dat, aes(x = t, y = wnti)) + geom_line() + ggtitle("with linear trend")
p3 grid.arrange(p1, p2, p3, ncol = 3)
We can make a similar plot with AR(1) data. Ignore the warnings about not knowing how to pick the scale.
<- 0.8
beta1 <- arima.sim(TT, model = list(ar = beta1), sd = sd)
ar1 <- ar1 + intercept
ar1i <- ar1 + trend * (1:TT) + intercept
ar1ti <- data.frame(t = 1:TT, ar1 = ar1, ar1i = ar1i, ar1ti = ar1ti)
dat <- ggplot(dat, aes(x = t, y = ar1)) + geom_line() + ggtitle("AR1")
p4 <- ggplot(dat, aes(x = t, y = ar1i)) + geom_line() + ggtitle("with non-zero mean")
p5 <- ggplot(dat, aes(x = t, y = ar1ti)) + geom_line() + ggtitle("with linear trend")
p6
grid.arrange(p4, p5, p6, ncol = 3)
Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
5.2.3 Greek landing data
We will look at the anchovy data. Notice the two ==
in the subset call not one =
. We will use the Greek data before 1989 for the lab.
<- subset(landings, Species == "Anchovy" & Year <= 1989)$log.metric.tons
anchovy <- ts(anchovy, start = 1964) anchovyts
Plot the data.
plot(anchovyts, ylab = "log catch")
Questions to ask.
- Does it have a trend (goes up or down)? Yes, definitely
- Does it have a non-zero mean? Yes
- Does it look like it might be stationary around a trend? Maybe