plot(gnp)

plot(log(gnp))

z=diff(log(gnp))
plot(z)

acf(z)

pacf(z)

Looks like maybe MA(2) or AR(1) or something. Let’s try fitting a few.

(fit01 = arima( z, order=c(0,0,1) )) 
## 
## Call:
## arima(x = z, order = c(0, 0, 1))
## 
## Coefficients:
##          ma1  intercept
##       0.2719     0.0083
## s.e.  0.0549     0.0008
## 
## sigma^2 estimated as 9.305e-05:  log likelihood = 715.3,  aic = -1424.6
(fit02 = arima( z, order=c(0,0,2) )) 
## 
## Call:
## arima(x = z, order = c(0, 0, 2))
## 
## Coefficients:
##          ma1     ma2  intercept
##       0.3028  0.2035     0.0083
## s.e.  0.0654  0.0644     0.0010
## 
## sigma^2 estimated as 8.919e-05:  log likelihood = 719.96,  aic = -1431.93
(fit03 = arima( z, order=c(0,0,3) )) 
## 
## Call:
## arima(x = z, order = c(0, 0, 3))
## 
## Coefficients:
##          ma1     ma2     ma3  intercept
##       0.3208  0.2478  0.0909     0.0083
## s.e.  0.0662  0.0718  0.0701     0.0010
## 
## sigma^2 estimated as 8.853e-05:  log likelihood = 720.78,  aic = -1431.55
(fit10 = arima( z, order=c(1,0,0) )) 
## 
## Call:
## arima(x = z, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.3467     0.0083
## s.e.  0.0627     0.0010
## 
## sigma^2 estimated as 9.03e-05:  log likelihood = 718.61,  aic = -1431.22
(fit11 = arima( z, order=c(1,0,1) )) 
## 
## Call:
## arima(x = z, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.4632  -0.1306     0.0083
## s.e.  0.1272   0.1350     0.0010
## 
## sigma^2 estimated as 8.994e-05:  log likelihood = 719.04,  aic = -1430.08
(fit12 = arima( z, order=c(1,0,2) )) 
## 
## Call:
## arima(x = z, order = c(1, 0, 2))
## 
## Coefficients:
##          ar1     ma1     ma2  intercept
##       0.2407  0.0761  0.1623     0.0083
## s.e.  0.2066  0.2026  0.0851     0.0010
## 
## sigma^2 estimated as 8.877e-05:  log likelihood = 720.47,  aic = -1430.95
(fit13 = arima( z, order=c(1,0,3) )) 
## 
## Call:
## arima(x = z, order = c(1, 0, 3))
## 
## Coefficients:
##           ar1     ma1     ma2     ma3  intercept
##       -0.0441  0.3638  0.2607  0.1006     0.0083
## s.e.   0.4126  0.4067  0.1390  0.1111     0.0010
## 
## sigma^2 estimated as 8.852e-05:  log likelihood = 720.78,  aic = -1429.56
(fit20 = arima( z, order=c(2,0,0) ))
## 
## Call:
## arima(x = z, order = c(2, 0, 0))
## 
## Coefficients:
##          ar1     ar2  intercept
##       0.3180  0.0820     0.0083
## s.e.  0.0667  0.0668     0.0011
## 
## sigma^2 estimated as 8.968e-05:  log likelihood = 719.36,  aic = -1430.72
(fit21 = arima( z, order=c(2,0,1) )) 
## 
## Call:
## arima(x = z, order = c(2, 0, 1))
## 
## Coefficients:
##           ar1     ar2     ma1  intercept
##       -0.1420  0.2512  0.4613     0.0083
## s.e.   0.4469  0.1465  0.4560     0.0010
## 
## sigma^2 estimated as 8.929e-05:  log likelihood = 719.84,  aic = -1429.67
(fit22 = arima( z, order=c(2,0,2) )) 
## 
## Call:
## arima(x = z, order = c(2, 0, 2))
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  intercept
##       1.3459  -0.7378  -1.0634  0.5621     0.0083
## s.e.  0.1374   0.1540   0.1873  0.1971     0.0008
## 
## sigma^2 estimated as 8.649e-05:  log likelihood = 723.29,  aic = -1434.57
(fit23 = arima( z, order=c(2,0,3) )) 
## 
## Call:
## arima(x = z, order = c(2, 0, 3))
## 
## Coefficients:
##          ar1      ar2      ma1     ma2      ma3  intercept
##       1.3239  -0.4418  -1.0291  0.2645  -0.1393     0.0084
## s.e.  0.1777   0.1798   0.1789  0.1606   0.0877     0.0005
## 
## sigma^2 estimated as 8.653e-05:  log likelihood = 723.23,  aic = -1432.46
(fit30 = arima( z, order=c(3,0,0) ))
## 
## Call:
## arima(x = z, order = c(3, 0, 0))
## 
## Coefficients:
##          ar1     ar2      ar3  intercept
##       0.3271  0.1178  -0.1104     0.0083
## s.e.  0.0665  0.0697   0.0666     0.0010
## 
## sigma^2 estimated as 8.857e-05:  log likelihood = 720.73,  aic = -1431.45
(fit31 = arima( z, order=c(3,0,1) )) 
## 
## Call:
## arima(x = z, order = c(3, 0, 1))
## 
## Coefficients:
##          ar1      ar2      ar3      ma1  intercept
##       1.0905  -0.1251  -0.1739  -0.7881     0.0084
## s.e.  0.3028   0.1529   0.0777   0.3113     0.0006
## 
## sigma^2 estimated as 8.682e-05:  log likelihood = 722.89,  aic = -1433.78
(fit32 = arima( z, order=c(3,0,2) )) 
## 
## Call:
## arima(x = z, order = c(3, 0, 2))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2  intercept
##       1.6742  -1.3104  0.2351  -1.4045  1.0000     0.0083
## s.e.  0.0666   0.0974  0.0656   0.0196  0.0234     0.0009
## 
## sigma^2 estimated as 8.22e-05:  log likelihood = 726.62,  aic = -1439.25
(fit33 = arima( z, order=c(3,0,3) )) 
## 
## Call:
## arima(x = z, order = c(3, 0, 3))
## 
## Coefficients:
##          ar1     ar2      ar3      ma1      ma2     ma3  intercept
##       0.3610  0.5871  -0.7044  -0.0750  -0.4759  0.5141     0.0083
## s.e.  0.1454  0.1036   0.1541   0.1885   0.1251  0.1973     0.0008
## 
## sigma^2 estimated as 8.628e-05:  log likelihood = 723.55,  aic = -1431.09
AIC( fit01, fit02, fit03, fit10, fit11, fit12, fit13, fit20, fit21, fit22, fit23, fit30, fit31, fit32, fit33)
##       df       AIC
## fit01  3 -1424.600
## fit02  4 -1431.929
## fit03  5 -1431.552
## fit10  3 -1431.221
## fit11  4 -1430.080
## fit12  5 -1430.948
## fit13  6 -1429.564
## fit20  4 -1430.724
## fit21  5 -1429.673
## fit22  6 -1434.571
## fit23  7 -1432.461
## fit30  5 -1431.453
## fit31  6 -1433.781
## fit32  7 -1439.248
## fit33  8 -1431.090
BIC( fit01, fit02, fit03, fit10, fit11, fit12, fit13, fit20, fit21, fit22, fit23, fit30, fit31, fit32, fit33)
##       df       BIC
## fit01  3 -1414.392
## fit02  4 -1418.319
## fit03  5 -1414.539
## fit10  3 -1421.013
## fit11  4 -1416.469
## fit12  5 -1413.935
## fit13  6 -1409.148
## fit20  4 -1417.113
## fit21  5 -1412.660
## fit22  6 -1414.155
## fit23  7 -1408.642
## fit30  5 -1414.440
## fit31  6 -1413.365
## fit32  7 -1415.429
## fit33  8 -1403.869

AIC prefers ARMA(3,2), BIC prefers AR((1). This is not unusual. AIC tends to overfit.

Let’s look at some diagnostics for the AR(1).

r = resid(fit10)
lags = 10
acf(r)

Box.test( r, 10, fitdf=2 )
## 
##  Box-Pierce test
## 
## data:  r
## X-squared = 10.445, df = 8, p-value = 0.2352
qqnorm(r)
qqline(r)

jarque.bera.test(r)
## 
##  Jarque Bera Test
## 
## data:  r
## X-squared = 23.573, df = 2, p-value = 7.607e-06

There is no evidence of seasonality. The AR(2) does a good job of eliminating autocorrelation, but the residuals don’t look very normal. We fail to reject the null hypothesis that the residuals are white noise based on the Box-Pierce test, but we reject that they are normall distributed based on the Jarque-Bera test.