plot(jj)

plot(log(jj))

logjj = log(jj)
Q = factor( cycle( jj ))
trend = time(jj) - 1970
fit01 = lm( logjj ~ 0 + trend + Q )
summary(fit01)
## 
## Call:
## lm(formula = logjj ~ 0 + trend + Q)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29318 -0.09062 -0.01180  0.08460  0.27644 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## trend 0.167172   0.002259   74.00   <2e-16 ***
## Q1    1.052793   0.027359   38.48   <2e-16 ***
## Q2    1.080916   0.027365   39.50   <2e-16 ***
## Q3    1.151024   0.027383   42.03   <2e-16 ***
## Q4    0.882266   0.027412   32.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1254 on 79 degrees of freedom
## Multiple R-squared:  0.9935, Adjusted R-squared:  0.9931 
## F-statistic:  2407 on 5 and 79 DF,  p-value: < 2.2e-16
ts.plot( cbind(logjj, fitted(fit01)), col=1:2, type="o", main="JJ, quarterly earnings per share", ylab="Dollars")

plot.ts(residuals(fit01), type="o", main="JJ earnings, residuals", ylab="Dollars")

acf(residuals(fit01), main="ACF, JJ residuals")

pacf(residuals(fit01), main="PACF, JJ residuals")

The residuals do not look like white noise. There is too much autocorrelation.

Let’s try seasonal ARIMA, with differenced data at the annual frequency.

fit02 = sarima( logjj, 0,0,0, 0,1,0, 4, details=FALSE )

AIC( fit01, fit02$fit )
## Warning in AIC.default(fit01, fit02$fit): models are not all fitted to the
## same number of observations
##           df       AIC
## fit01      6 -103.6153
## fit02$fit  2 -146.4593
BIC( fit01, fit02$fit )
## Warning in BIC.default(fit01, fit02$fit): models are not all fitted to the
## same number of observations
##           df        BIC
## fit01      6  -89.03043
## fit02$fit  2 -141.69529
plot( residuals( fit02$fit ))

acf( residuals( fit02$fit ))

pacf( residuals( fit02$fit ))

Definitely still some autocorrelation in the residuals…

Let’s try a few more models…

fit04 = sarima( logjj, 1,0,0, 0,1,0, 4, details=FALSE )

fit05 = sarima( logjj, 0,0,1, 0,1,0, 4, details=FALSE )

fit06 = sarima( logjj, 1,0,1, 0,1,0, 4, details=FALSE )

fit07 = sarima( logjj, 2,0,0, 0,1,0, 4, details=FALSE )

fit08 = sarima( logjj, 2,0,1, 0,1,0, 4, details=FALSE )

fit09 = sarima( logjj, 1,0,2, 0,1,0, 4, details=FALSE )

AIC( fit01, fit02$fit, fit04$fit, fit05$fit, fit06$fit, fit07$fit, fit08$fit, fit09$fit )
## Warning in AIC.default(fit01, fit02$fit, fit04$fit, fit05$fit, fit06$fit, :
## models are not all fitted to the same number of observations
##           df       AIC
## fit01      6 -103.6153
## fit02$fit  2 -146.4593
## fit04$fit  3 -151.7843
## fit05$fit  3 -149.3911
## fit06$fit  4 -150.9915
## fit07$fit  4 -152.3325
## fit08$fit  5 -150.8571
## fit09$fit  5 -153.7561
BIC( fit01, fit02$fit, fit04$fit, fit05$fit, fit06$fit, fit07$fit, fit08$fit, fit09$fit )
## Warning in BIC.default(fit01, fit02$fit, fit04$fit, fit05$fit, fit06$fit, :
## models are not all fitted to the same number of observations
##           df        BIC
## fit01      6  -89.03043
## fit02$fit  2 -141.69529
## fit04$fit  3 -144.63819
## fit05$fit  3 -142.24504
## fit06$fit  4 -141.46339
## fit07$fit  4 -142.80434
## fit08$fit  5 -138.94695
## fit09$fit  5 -141.84600
fit04
## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
##     reltol = tol))
## 
## Coefficients:
##          ar1  constant
##       0.3155    0.0386
## s.e.  0.1137    0.0037
## 
## sigma^2 estimated as 0.008136:  log likelihood = 78.89,  aic = -151.78
## 
## $degrees_of_freedom
## [1] 82
## 
## $ttable
##          Estimate     SE t.value p.value
## ar1        0.3155 0.1137  2.7758  0.0068
## constant   0.0386 0.0037 10.5049  0.0000
## 
## $AIC
## [1] -3.763872
## 
## $AICc
## [1] -3.736491
## 
## $BIC
## [1] -4.705996
fit07
## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
##     reltol = tol))
## 
## Coefficients:
##          ar1     ar2  constant
##       0.2609  0.1887    0.0381
## s.e.  0.1165  0.1171    0.0045
## 
## sigma^2 estimated as 0.007873:  log likelihood = 80.17,  aic = -152.33
## 
## $degrees_of_freedom
## [1] 81
## 
## $ttable
##          Estimate     SE t.value p.value
## ar1        0.2609 0.1165  2.2398  0.0278
## ar2        0.1887 0.1171  1.6116  0.1109
## constant   0.0381 0.0045  8.5138  0.0000
## 
## $AIC
## [1] -3.772874
## 
## $AICc
## [1] -3.743037
## 
## $BIC
## [1] -4.686059
fit09
## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1, 
##     reltol = tol))
## 
## Coefficients:
##          ar1     ma1     ma2  constant
##       0.1830  0.0890  0.3578    0.0382
## s.e.  0.2416  0.2156  0.1200    0.0043
## 
## sigma^2 estimated as 0.007526:  log likelihood = 81.88,  aic = -153.76
## 
## $degrees_of_freedom
## [1] 80
## 
## $ttable
##          Estimate     SE t.value p.value
## ar1        0.1830 0.2416  0.7573  0.4511
## ma1        0.0890 0.2156  0.4128  0.6809
## ma2        0.3578 0.1200  2.9813  0.0038
## constant   0.0382 0.0043  8.9644  0.0000
## 
## $AIC
## [1] -3.794191
## 
## $AICc
## [1] -3.761224
## 
## $BIC
## [1] -4.678438

BIC prefers \(\text{ARIMA}(1,0,0)\times(0,1,0)_{12}\). AIC prefers a couple of larger models, but notice that the additional parameters are not significant at the 0.05 level. I would tend to go with the BIC. The diagnostics look good.

```