Topic 3 Part 1 – Model Selection STAT 505 – Time Series Analysis Lecture: TR 13:30-14:45 [1] Order Selection for an ARMA Models Using ACF and PACF to identify the order of AR(p), MA(q) models. Characteristics of ACF and PACF Process ACF PACF AR(p) Tails off as exponential decay or damped sine wave Cuts off after lag p MA(q) Cuts off after lag q Tails off as exponential decay of damped sine wave ARMA(p,q) Tails off Tails off [2] > par(mfrow=c(3,2)) > ACF = ARMAacf(ar=-0.7,lag.max=24) > PACF = ARMAacf(ar=-0.7, lag.max=24,pacf=TRUE) > plot(ACF, type=”h”, xlab=”lag”, ylim=c(-1,1),main=”True ACF for ARMA(1,0)”); > abline(h=0) > plot(PACF, type=”h”, xlab=”lag”, ylim=c(-1,1),main=”True PACF for ARMA(1,0)”); > abline(h=0) > ACF = ARMAacf(ma=-0.7,lag.max=24) > PACF = ARMAacf(ma=-0.7, lag.max=24,pacf=TRUE) > plot(ACF, type=”h”, xlab=”lag”, ylim=c(-1,1),main=”True ACF for ARMA(0,1)”); > abline(h=0) > plot(PACF, type=”h”, xlab=”lag”, ylim=c(-1,1),main=”True PACF for ARMA(0,1)”); > abline(h=0) > ACF = ARMAacf(ar=-0.7,ma=-0.5,lag.max=24) > PACF = ARMAacf(ar=-0.7,ma=-0.7, lag.max=24,pacf=TRUE) > plot(ACF, type=”h”, xlab=”lag”, ylim=c(-1,1),main=”True ACF for ARMA(1,1)”); > abline(h=0) > plot(PACF, type=”h”, xlab=”lag”, ylim=c(-1,1),main=”True PACF for ARMA(1,1)”); > abline(h=0) [3] 5 10 15 20 25 1. 0 0. 0 1. 0 True ACF for ARMA(1,0) lag AC F 5 10 15 20 1. 0 0. 0 1. 0 True PACF for ARMA(1,0) lag PA CF 5 10 15 20 25 1. 0 0. 0 1. 0 True ACF for ARMA(0,1) lag AC F 5 10 15 20 1. 0 0. 0 1. 0 True PACF for ARMA(0,1) lag PA CF 5 10 15 20 25 1. 0 0. 0 1. 0 True ACF for ARMA(1,1) lag AC F 5 10 15 20 1. 0 0. 0 1. 0 True PACF for ARMA(1,1) lag PA CF [4] [5] There are commonly two approaches: Try-and-error and information criteria 1. Try-and-error to select ARMA models For the order selection of ARMA models, we may use try-and-error to find the order of p and q. For any pair of (p, q), check the residuals to see whether the ACF and PACF display the patterns for a AR or MA model. Then continue the fitting process until the residual is a white noise. According to the large-sample result, if the residuals Zt ( obtained from the modeling of a time series {Xt} ) is a iid sequence, then the components of ACF ρk and PACF φ mm of {Zt} follows approximately normal(0, 1/n) distribution. That is for k > 0, φ kk will fall between the bounds ±1.96n 1/2 with probability close to 0.95. This suggests that φ kk < 1.96n 1/2 for all k > 0, as a preliminary estimator of p and q. R-code example for a data set, say W1, we fit ARMA(2,1) model: library(stats) w1.fit <- arima(data, order=c(2, 0, 1)) #The output of this code provides the log-likelihood and AIC value. [6] # You can also compute BIC values. data_res <- w1.fit$residual plot(data_res) acf(data_res) pacf(data_res) tsdiag(data_res) # We may use this to analyze model residual Box.test(data_res, lag = 20, type=`Ljung') 2. Using information criteria: A more systematic approach to order selection is to find the values of p and q that minimize the certain information criteria such as AIC or BIC. (a) AIC (or AICc) and BIC AIC = 2k 2(log likelihood) where k = p+ q + 1 is the number of parameters. AICc(l) = 2(p+ q + 1)n/(n p q 2) 2(log likelihood) BIC(l) = k ln(n) 2(log likelihood) where n is the number of observations. (b) To apply the AIC (or AICc) and BIC criteria, we choose the value of p and q to minimize AIC or BIC. [7] (c) AIC and BIC are usually given in the model output. > Dirdata=”~/OneDrive – University of Calgary//MyCoursesThierry/STAT505/Winter2019/DataandITSM/ITSM2000-2/” > sunspot=scan(paste(Dirdata,”sunspots.tsm”,sep=””)) > sunspot.ts=ts(data=sunspot, frequency = 1,start = 1770,end = 1869) > par(mfrow=c(1,3)) > plot.ts(sunspot.ts,lwd=1.5) > acf(sunspot.ts,lag.max = 50,lwd=1.5) > pacf(sunspot.ts,lag.max = 50,lwd=1.5) > AIC1=BIC1=matrix(0,4,4) > for (i in 1:4){ for (j in 1:4){ + fit=arima(sunspot.ts, order = c(i-1, 0,j-1)) + AIC1[i,j]=AIC(fit); + BIC1[i,j]=BIC(fit); }} > AIC1 [,1] [,2] [,3] [,4] [1,] 1010.9298 906.9696 852.8261 843.4655 [2,] 903.8076 848.8172 838.9291 839.7595 [3,] 837.2348 833.0532 834.6754 836.6731 [4,] 834.9896 834.7797 836.6897 838.4875 > BIC1 [,1] [,2] [,3] [,4] [1,] 1016.1401 914.7851 863.2468 856.4913 [8] [2,] 911.6231 859.2379 851.9549 855.3905 [3,] 847.6555 846.0791 850.3064 854.9093 [4,] 848.0155 850.4107 854.9259 859.3288 > idxminA=which(AIC1 == min(AIC1), arr.ind = TRUE) > idxminA row col [1,] 3 2 > idxminB=which(BIC1 == min(BIC1), arr.ind = TRUE) > idxminB row col [1,] 3 2 Time su n sp ot .ts 1780 1800 1820 1840 1860 0 50 10 0 15 0 0 10 20 30 40 50 0. 2 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Lag AC F Series sunspot.ts 0 10 20 30 40 50 0. 5 0. 0 0. 5 Lag Pa rti al A CF Series sunspot.ts > fit=arima(sunspot.ts, order = c(idxminA[1]-1, 0,idxminA[2]-1)) > library(lmtest) # install.packages(“lmtest”) [9] > coeftest(fit) z test of coefficients: Estimate Std. Error z value Pr(>|z|) ar1 1.22482 0.11310 10.8291 < 2.2e-16 *** ar2 -0.56008 0.10830 -5.1716 2.321e-07 *** ma1 0.38467 0.13342 2.8832 0.003937 ** intercept 48.46201 6.01870 8.0519 8.152e-16 *** --- Signif. codes: 0 > tsdiag(fit) [10] Standardized Residuals Time 1780 1800 1820 1840 1860 2 0 2 0 5 10 15 20 0. 2 0. 4 1. 0 Lag AC F ACF of Residuals 2
4 6 8 10 0. 0 0. 4 0. 8 p values for Ljung Box statistic lag p va lu e [11] You may also use R code: auto.arima (from the r package forecast) function for modeling. For seasonal models, we will discuss later on. You can still use the same function for models with seasonality. > library(forecast) > auto.arima(sunspot.ts) Series: sunspot.ts ARIMA(2,0,1) with non-zero mean Coefficients: ar1 ar2 ma1 mean 1.2248 -0.5601 0.3847 48.4620 s.e. 0.1131 0.1083 0.1334 6.0187 sigma^2 estimated as 222.7: log likelihood=-411.53 AIC=833.05 AICc=833.69 BIC=846.08 [12] Case Study: Analysis of the monthly accidental deaths (Modified from Example 6.5.4) > Dirdata=”~/OneDrive – University of Calgary/MyCoursesThierry/STAT505/Winter2019/DataandITSM/ITSM2000-2/” > death=scan(paste(Dirdata,”deaths.tsm”,sep=””)) > death.ts=ts(data=death, frequency = 12,start = c(1973,1),end = c(1978,12)) > par(mfrow=c(3,1)) > plot.ts(death.ts,main=”The monthly accidental deaths data, 1973-1978″ + ,ylab=”(thousands)”,lwd=1.5) > acf(death.ts,lag.max = 48,xlab=”Year”,lwd=1.5) > pacf(death.ts,lag.max = 48,xlab=”Year”,lwd=1.5) [13] The monthly accidental deaths data, 1973 1978 Time (th ou sa nd s) 1973 1974 1975 1976 1977 1978 1979 70 00 10 00 0 0 1 2 3 4 0. 4 0. 2 0. 8 Year AC F Series death.ts 0 1 2 3 4 0. 4 0. 2 Year Pa rti al A CF Series death.ts [14] > NoSeaso.death=diff(death.ts,lag=12) > # We use the operator 1-B^{12} to remove the seasonality > par(mfrow=c(3,1)) > plot.ts(NoSeaso.death,main=”The differenced series (1-B^{12})X_t” + ,ylab=”(thousands)”,lwd=1.5) > acf(NoSeaso.death,lag.max = 60,lwd=1.5) > pacf(NoSeaso.death,lag.max = 60,lwd=1.5) > #The seasonal component evident previously is absent from the graph, > # However there appears to be a nondecreasing trend. It can be removed > # by applying the operator 1-B [15] The differenced series (1 B^{12})X_t Time (th ou sa nd s) 1974 1975 1976 1977 1978 1979 10 00 0 0 1 2 3 4 5 0. 2 0. 4 1. 0 Lag AC F Series NoSeaso.death 0 1 2 3 4 5 0. 2 0. 2 0. 6 Lag Pa rti al A CF Series NoSeaso.death [16] > NoTrend.Seaso.death=diff(NoSeaso.death,lag=1) > par(mfrow=c(3,1)) > plot.ts(NoTrend.Seaso.death,main=”The differenced series (1-B)(1-B^{12})X_t” + ,ylab=”(thousands)”,lwd=1.5) > acf(NoTrend.Seaso.death,lag.max = 48,xlab=”Year”,lwd=1.5) > pacf(NoTrend.Seaso.death,lag.max = 48,xlab=”Year”,lwd=1.5) [17] The differenced series (1 B)(1 B^{12})X_t Time (th ou sa nd s) 1974 1975 1976 1977 1978 1979 50 0 50 0 0 1 2 3 4 0. 4 0. 2 0. 8 Year AC F Series NoTrend.Seaso.death 0 1 2 3 4 0. 3 0. 0 Year Pa rti al A CF Series NoTrend.Seaso.death [18] > fit1=arima(NoTrend.Seaso.death,order=c(0,0,1), + seasonal = list(order=c(0,0,1),period=12)) > library(lmtest) > coeftest(fit1) z test of coefficients: Estimate Std. Error z value Pr(>|z|) ma1 -0.49626 0.13512 -3.6726 0.0002401 *** sma1 -0.61460 0.19324 -3.1805 0.0014700 ** intercept 21.02204 11.99207 1.7530 0.0796029 . — Signif. codes: 0 > tsdiag(fit1,lwd=1.5) > print(AIC(fit1)) [1] 856.5324 > print(BIC(fit1)) [1] 864.8426 [19] Standardized Residuals Time 1974 1975 1976 1977 1978 1979 2 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0. 2 0. 4 1. 0 Lag AC F ACF of Residuals 2 4 6 8 10 0. 0 0. 4 0. 8 p values for Ljung Box statistic lag p va lu e [20] > fit2=arima(NoTrend.Seaso.death,order=c(1,0,1), + seasonal = list(order=c(1,0,1),period=12)) > coeftest(fit2) z test of coefficients: Estimate Std. Error z value Pr(>|z|) ar1 0.447582 0.120267 3.7216 0.000198 *** ma1 -0.999911 0.074372 -13.4448 < 2.2e-16 *** sar1 0.232117 0.161734 1.4352 0.151235 sma1 -0.998135 0.695532 -1.4351 0.151268 intercept 20.750229 3.430254 6.0492 1.456e-09 *** --- Signif. codes: 0 > tsdiag(fit2,lwd=1.5) > print(AIC(fit2)) [1] 855.1901 > print(BIC(fit2)) [1] 867.6553 [21] Standardized Residuals Time 1974 1975 1976 1977 1978 1979 1 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0. 2 0. 4 1. 0 Lag AC F ACF of Residuals 2 4 6 8 10 0. 0 0. 4 0. 8 p values for Ljung Box statistic lag p va lu e [22] > fit3=arima(NoTrend.Seaso.death,order=c(1,0,1), + seasonal = list(order=c(1,0,0),period=12)) > coeftest(fit3) z test of coefficients: Estimate Std. Error z value Pr(>|z|) ar1 0.431752 0.122144 3.5348 0.0004081 *** ma1 -0.999997 0.055902 -17.8883 < 2.2e-16 *** sar1 -0.399054 0.121840 -3.2752 0.0010557 ** intercept 20.611825 3.306490 6.2337 4.554e-10 *** --- Signif. codes: 0 > tsdiag(fit3,lwd=1.5) > print(AIC(fit3)) [1] 857.9437 > print(BIC(fit3)) [1] 868.3314 [23] Standardized Residuals Time 1974 1975 1976 1977 1978 1979 1 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0. 2 0. 4 1. 0 Lag AC F ACF of Residuals 2 4 6 8 10 0. 0 0. 4 0. 8 p values for Ljung Box statistic lag p va lu e [24] > fit4=arima(NoTrend.Seaso.death,order=c(0,0,1), + seasonal = list(order=c(1,0,0),period=12)) > coeftest(fit4) z test of coefficients: Estimate Std. Error z value Pr(>|z|) ma1 -0.51869 0.13490 -3.8450 0.0001205 *** sar1 -0.35165 0.12671 -2.7753 0.0055155 ** intercept 22.46605 16.43886 1.3666 0.1717374 — Signif. codes: 0 > tsdiag(fit4,lwd=1.5) > print(AIC(fit4)) [1] 860.6085 > print(BIC(fit4)) [1] 868.9187 [25] Standardized Residuals Time 1974 1975 1976 1977 1978 1979 2 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0. 2 0. 4 1. 0 Lag AC F ACF of Residuals 2 4 6 8 10 0. 0 0. 4 0. 8 p values for Ljung Box statistic lag p va lu e