plot.ts( cbind(cmort, part, tempr))

pairs( cbind( cmort, part, tempr ))

trend = time(cmort)
temp = tempr - mean(tempr)
temp2 = temp^2

fit01 = lm( cmort ~ trend + temp + temp2 + part )
summary(fit01)
## 
## Call:
## lm(formula = cmort ~ trend + temp + temp2 + part)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.0760  -4.2153  -0.4878   3.7435  29.2448 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.831e+03  1.996e+02   14.19  < 2e-16 ***
## trend       -1.396e+00  1.010e-01  -13.82  < 2e-16 ***
## temp        -4.725e-01  3.162e-02  -14.94  < 2e-16 ***
## temp2        2.259e-02  2.827e-03    7.99 9.26e-15 ***
## part         2.554e-01  1.886e-02   13.54  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.385 on 503 degrees of freedom
## Multiple R-squared:  0.5954, Adjusted R-squared:  0.5922 
## F-statistic:   185 on 4 and 503 DF,  p-value: < 2.2e-16
r = resid( fit01 )
plot(r, type="l")

acf2(r)

##         ACF  PACF
##  [1,]  0.34  0.34
##  [2,]  0.44  0.36
##  [3,]  0.28  0.08
##  [4,]  0.28  0.06
##  [5,]  0.16 -0.05
##  [6,]  0.12 -0.05
##  [7,]  0.07 -0.02
##  [8,]  0.01 -0.05
##  [9,]  0.03  0.02
## [10,] -0.05 -0.06
## [11,] -0.02  0.00
## [12,]  0.00  0.06
## [13,] -0.04 -0.02
## [14,] -0.02  0.00
## [15,]  0.01  0.04
## [16,] -0.05 -0.08
## [17,] -0.01  0.00
## [18,] -0.03 -0.01
## [19,] -0.06 -0.06
## [20,] -0.03  0.03
## [21,] -0.03  0.02
## [22,] -0.05 -0.03
## [23,]  0.00  0.05
## [24,] -0.02 -0.01
## [25,] -0.01  0.00
## [26,] -0.05 -0.06
## [27,]  0.03  0.06
## [28,] -0.11 -0.11
## [29,] -0.02  0.00
## [30,] -0.10 -0.05
## [31,] -0.02  0.05
## [32,] -0.09 -0.04
## [33,] -0.03  0.04

The residuals are correlated. Based on ACF and PACF, it looks like AR(2) might be a good model. But, we’ll try a few other alternatives as well.

(fit00 = arima( cmort, order=c(0,0,0), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(0, 0, 0), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##       intercept  trend + temp + temp2 + part
##         37.9206                       0.0241
## s.e.     8.6567                       0.0041
## 
## sigma^2 estimated as 93.43:  log likelihood = -1873.28,  aic = 3752.56
(fit10 = arima( cmort, order=c(1,0,0), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(1, 0, 0), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ar1  intercept  trend + temp + temp2 + part
##       0.7862    53.0703                       0.0170
## s.e.  0.0273     4.6462                       0.0021
## 
## sigma^2 estimated as 35.81:  log likelihood = -1630.17,  aic = 3268.35
(fit01 = arima( cmort, order=c(0,0,1), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(0, 0, 1), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ma1  intercept  trend + temp + temp2 + part
##       0.4955    46.3444                       0.0201
## s.e.  0.0286     6.3093                       0.0030
## 
## sigma^2 estimated as 61:  log likelihood = -1765.12,  aic = 3538.24
(fit20 = arima( cmort, order=c(2,0,0), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(2, 0, 0), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ar1     ar2  intercept  trend + temp + temp2 + part
##       0.4300  0.4532    53.1733                        0.017
## s.e.  0.0394  0.0396     4.7093                        0.002
## 
## sigma^2 estimated as 28.45:  log likelihood = -1571.98,  aic = 3153.95
(fit11 = arima( cmort, order=c(1,0,1), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(1, 0, 1), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ar1      ma1  intercept  trend + temp + temp2 + part
##       0.9423  -0.4514    51.5990                       0.0177
## s.e.  0.0163   0.0380     5.0385                       0.0022
## 
## sigma^2 estimated as 29.66:  log likelihood = -1582.45,  aic = 3174.9
(fit02 = arima( cmort, order=c(0,0,2), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(0, 0, 2), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ma1     ma2  intercept  trend + temp + temp2 + part
##       0.4947  0.5168    55.5202                       0.0158
## s.e.  0.0444  0.0369     5.1356                       0.0024
## 
## sigma^2 estimated as 44.14:  log likelihood = -1683.17,  aic = 3376.33
(fit30 = arima( cmort, order=c(3,0,0), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(3, 0, 0), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ar1     ar2     ar3  intercept  trend + temp + temp2 + part
##       0.4097  0.4338  0.0447    52.6029                       0.0173
## s.e.  0.0443  0.0441  0.0450     4.8148                       0.0021
## 
## sigma^2 estimated as 28.4:  log likelihood = -1571.48,  aic = 3154.96
(fit21 = arima( cmort, order=c(2,0,1), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(2, 0, 1), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ar1     ar2      ma1  intercept  trend + temp + temp2 + part
##       0.5020  0.3965  -0.0908    52.6152                       0.0173
## s.e.  0.0865  0.0736   0.0932     4.8091                       0.0021
## 
## sigma^2 estimated as 28.4:  log likelihood = -1571.52,  aic = 3155.04
(fit12 = arima( cmort, order=c(1,0,2), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(1, 0, 2), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ar1      ma1     ma2  intercept  trend + temp + temp2 + part
##       0.9173  -0.5024  0.1983    51.8192                       0.0177
## s.e.  0.0215   0.0480  0.0427     4.7529                       0.0021
## 
## sigma^2 estimated as 28.48:  log likelihood = -1572.27,  aic = 3156.54
(fit03 = arima( cmort, order=c(0,0,3), xreg=trend+temp+temp2+part))
## 
## Call:
## arima(x = cmort, order = c(0, 0, 3), xreg = trend + temp + temp2 + part)
## 
## Coefficients:
##          ma1     ma2     ma3  intercept  trend + temp + temp2 + part
##       0.5300  0.5656  0.3266    52.0839                       0.0174
## s.e.  0.0466  0.0363  0.0401     4.8427                       0.0023
## 
## sigma^2 estimated as 38.67:  log likelihood = -1649.61,  aic = 3311.22
AIC(fit00,fit10,fit01,fit20,fit11,fit02,fit30,fit21,fit12,fit03)
##       df      AIC
## fit00  3 3752.562
## fit10  4 3268.350
## fit01  4 3538.245
## fit20  5 3153.952
## fit11  5 3174.895
## fit02  5 3376.335
## fit30  6 3154.963
## fit21  6 3155.038
## fit12  6 3156.539
## fit03  6 3311.217
BIC(fit00,fit10,fit01,fit20,fit11,fit02,fit30,fit21,fit12,fit03)
##       df      BIC
## fit00  3 3765.253
## fit10  4 3285.272
## fit01  4 3555.167
## fit20  5 3175.105
## fit11  5 3196.048
## fit02  5 3397.487
## fit30  6 3180.346
## fit21  6 3180.421
## fit12  6 3181.922
## fit03  6 3336.600

AIC and BIC both suggest AR(2) as the preferred model. Let’s look at some diagnostics… (The sarima function provides diagnostics automatically, so let’s just use that.)

sarima(cmort, 2, 0, 0, xreg=trend+temp+temp2+part, details=FALSE )

## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, optim.control = list(trace = trc, REPORT = 1, 
##     reltol = tol))
## 
## Coefficients:
##          ar1     ar2  intercept   xreg
##       0.4300  0.4532    53.1733  0.017
## s.e.  0.0394  0.0396     4.7093  0.002
## 
## sigma^2 estimated as 28.45:  log likelihood = -1571.98,  aic = 3153.95
## 
## $degrees_of_freedom
## [1] 504
## 
## $ttable
##           Estimate     SE t.value p.value
## ar1         0.4300 0.0394 10.9021       0
## ar2         0.4532 0.0396 11.4573       0
## intercept  53.1733 4.7093 11.2912       0
## xreg        0.0170 0.0020  8.3705       0
## 
## $AIC
## [1] 4.363952
## 
## $AICc
## [1] 4.368124
## 
## $BIC
## [1] 3.397263

Looks good. Residuals are approximately normal. No significant autocorrelation.

```