Bunnies and Kitties
Jeff Pratt
4/26/2019
Hare -> Lynx
The data
This is the classic data set which contains the number of lynx and hare pelts (in thousands) that were traded with the Hudson Bay Trading Company between 1845 and 1935. It is available within the ecostudy package and I received it from Andy Bunn. The two species are related in that the lynx preys on the hare. An increase in hare would lead to an increase in lynx as there is more food available which would in turn reduce the hare population and in turn the lynx population. I would predict that there will be a stong cyclical nature (see the previous sentence and graph below) which will be evident in the ACF plot. There should be a strong MA model available for both the hare and the lynx populations due to their strong dependence on each other.
bk <- readRDS("LynxHare.rds")
str(bk)## Time-Series [1:91, 1:2] from 1845 to 1935: 19.6 19.6 19.6 12 28 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:2] "Hare" "Lynx"head(bk)## Hare Lynx
## [1,] 19.58 30.09
## [2,] 19.60 45.15
## [3,] 19.61 49.15
## [4,] 11.99 39.52
## [5,] 28.04 21.23
## [6,] 58.00 8.42plot(bk,main="Annual numbers of pelts in tousands", xlab="Year")
acf(bk[,1])
pacf(bk[,1])
acf(bk[,2])
pacf(bk[,2])
The models
I am going to use arima to fit an ARIMA(1,1) model to each univariate series seperately. This is just a jumping off point; if it works that is great but I am expecting to have to adjust it.
bunnies.arima11<-arima(bk[,1], order = c(1,0,1))
bunnies.arima11##
## Call:
## arima(x = bk[, 1], order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.4772 0.4079 44.3594
## s.e. 0.1166 0.1128 7.1154
##
## sigma^2 estimated as 651: log likelihood = -424.3, aic = 856.59kitties.arima11<-arima(bk[,2], order = c(1,0,1))
kitties.arima11##
## Call:
## arima(x = bk[, 2], order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.6539 0.5754 28.4511
## s.e. 0.0834 0.0762 4.5270
##
## sigma^2 estimated as 94.44: log likelihood = -336.85, aic = 681.71tsdiag(bunnies.arima11)
tsdiag(kitties.arima11)
We can see that this explains a good amount of the data, but there is still a pattern in the residuals. I will borrow a bit of EDA from Andy again. Below I run a multitude models and will look at the Akaike information criterion and the bayesian information criterion to see which model fits best. A lower score denotes a better fitting model.
y<-bk[,1]
arma10 <- arima(y,order = c(1,0,0))
arma20 <- arima(y,order = c(2,0,0))
arma01 <- arima(y,order = c(0,0,1))
arma02 <- arima(y,order = c(0,0,2))
arma11 <- arima(y,order = c(1,0,1))
arma21 <- arima(y,order = c(2,0,1))
arma12 <- arima(y,order = c(1,0,2))
arma22 <- arima(y,order = c(2,0,2))
arma00 <- arima(y,order = c(0,0,0))
armaBIC <- BIC(arma10,arma20,arma01,arma02,arma11,arma21,arma12,arma22,arma00)
armaBIC## df BIC
## arma10 3 872.8444
## arma20 4 860.2439
## arma01 3 874.1646
## arma02 4 866.4231
## arma11 4 866.6334
## arma21 5 857.4906
## arma12 5 869.2140
## arma22 6 861.9490
## arma00 2 920.3515armaAIC <- AIC(arma10,arma20,arma01,arma02,arma11,arma21,arma12,arma22,arma00)
armaAIC## df AIC
## arma10 3 865.3118
## arma20 4 850.2005
## arma01 3 866.6320
## arma02 4 856.3796
## arma11 4 856.5900
## arma21 5 844.9363
## arma12 5 856.6597
## arma22 6 846.8838
## arma00 2 915.3298We see that for the hare a ARMA(2,1) model fits best. We should look at what works best for the lynx. I would hope it is a similar model structure.
y<-bk[,2]
arma10 <- arima(y,order = c(1,0,0))
arma20 <- arima(y,order = c(2,0,0))
arma01 <- arima(y,order = c(0,0,1))
arma02 <- arima(y,order = c(0,0,2))
arma11 <- arima(y,order = c(1,0,1))
arma21 <- arima(y,order = c(2,0,1))
arma12 <- arima(y,order = c(1,0,2))
arma22 <- arima(y,order = c(2,0,2))
arma00 <- arima(y,order = c(0,0,0))
armaBIC <- BIC(arma10,arma20,arma01,arma02,arma11,arma21,arma12,arma22,arma00)
armaBIC## df BIC
## arma10 3 723.0004
## arma20 4 654.8141
## arma01 3 725.1711
## arma02 4 689.6670
## arma11 4 691.7516
## arma21 5 653.3290
## arma12 5 681.9763
## arma22 6 657.8332
## arma00 2 800.8479armaAIC <- AIC(arma10,arma20,arma01,arma02,arma11,arma21,arma12,arma22,arma00)
armaAIC## df AIC
## arma10 3 715.4679
## arma20 4 644.7707
## arma01 3 717.6385
## arma02 4 679.6235
## arma11 4 681.7082
## arma21 5 640.7747
## arma12 5 669.4220
## arma22 6 642.7681
## arma00 2 795.8262Great, the same structure fits well for this as well. Let’s look at the models and their diagnostics.
bunnies.arima21<-arima(bk[,1], order = c(2,0,1))
bunnies.arima21##
## Call:
## arima(x = bk[, 1], order = c(2, 0, 1))
##
## Coefficients:
## ar1 ar2 ma1 intercept
## 1.3811 -0.7120 -0.5747 45.1433
## s.e. 0.1016 0.0784 0.1240 3.2422
##
## sigma^2 estimated as 557.6: log likelihood = -417.47, aic = 844.94kitties.arima21<-arima(bk[,2], order = c(2,0,1))
kitties.arima21##
## Call:
## arima(x = bk[, 2], order = c(2, 0, 1))
##
## Coefficients:
## ar1 ar2 ma1 intercept
## 1.4851 -0.8468 -0.3392 28.1029
## s.e. 0.0652 0.0571 0.1185 1.4796
##
## sigma^2 estimated as 58.24: log likelihood = -315.39, aic = 640.77tsdiag(bunnies.arima21)
tsdiag(kitties.arima21)
The diagnostics look much better than (1,1) and the coefficients can offer some interesting conclusions. I don’t quite feel confident explaining the coefficient values so I am going to run some simulations. The wizards on stackexchange explain a bit but a plot helps me understand.
n <- 100
epsilon <- rnorm(n=n)
y <- numeric(length = n)
phi <- 1.4
phi2<-0.0
theta <- 0.0
for(i in 3:n){
y[i] <- phi * y[i-1] +phi2 * y[i-2]+ theta * epsilon[i-1] + epsilon[i]
}
plot(y,type="l")
acf(y)
pacf(y)
So without a negative phi2, phi1 > 1 leads to a non-stationary rocket/drill. Adding in a negative AR(2) process will bring balance to the system and center it around a mean? Sum(phi1+phi2)!>=1
n <- 100
epsilon <- rnorm(n=n)
y <- numeric(length = n)
phi <- 1.4
phi2<- -0.7
theta <- 0.0
for(i in 3:n){
y[i] <- phi * y[i-1] +phi2 * y[i-2]+ theta * epsilon[i-1] + epsilon[i]
}
plot(y,type="l")
acf(y)
pacf(y)
Wow I must have a great teacher. I expect adding in a negative moving average will result in coefficients with lower amplitudes.
n <- 100
epsilon <- rnorm(n=n)
y <- numeric(length = n)
phi <- 1.4
phi2<- -0.7
theta <- -0.6
for(i in 3:n){
y[i] <- phi * y[i-1] +phi2 * y[i-2]+ theta * epsilon[i-1] + epsilon[i]
}
plot(y,type="l")
acf(y)
pacf(y)
Dead station
Describing Random
The way that I best understand the lack of stationarity in a random walk ts is through the use of images.
This image is from https://www.softcover.io/read/bf34ea25/math_for_finance/random_walks and while it does not show the distribtutions, it is still useful for me.
It shows the multiple iterations of a single model with conditional probabilites at each time set. That is why the model is not stationary. A stationary model draws from a single distribution for the entire model while a random walk draws from a specific distribution for each point in time based on the previous value. It has a probability to express as stationary like some of the iterations in the image above. They stay centered around a mean of 0 but after only 100 iterations many have already marched away from 0.
Math?
Thank God I read the book. 2.2.6 and 3.1.2
My attempt at making a stationary process out of a random walk would be to model ϵt as a function of time. yt - yt−1=ϵt allows you describe the ‘white noise’. This sounded wrong to me but my friend investopedia backed me up. I didn’t know the word to describe the process that I had created by doing this but searching for difference stationary opened up Box-Pandora’s box of differencing and other things I don’t comprehend. My best attempt is that the ‘difference’ is the result of subtracting yt−1 from yt and could continue to be ‘differenced’ until stationarity is achieved.
One could plot the values for ϵt as a function of time and this is where my biggest problem of uderstanding is. Won’t ϵt have a trend? If it was part of a normal distrubtion then the model would appear stationary anyways?
Why?
Stationarity is desirable because a described model with a uniform mean, variance, and ARMA characteristics can be used to forecast and predict while a random model would require information about prior conditions to predict. That would increase the uncertainty quite a lot and reduce the utility of the model.