CUHK- Shenzhen STAT 401 - Seasonal and Related Models

What is this all about?Data/Theory goals of the topicSeasonal modelsSeasonal unit roots and related discussionPeriodic modelsSome final notesReferencesSeasonal and Related ModelsVladas Pipiras, STOR @ UNC-CHFebruary, 2022What is this all about?Consider the following time series and its analysis.ts0 <- scan("airpass.txt")ts <- log(ts0)par(mfrow = c(1, 2))plot.ts(ts0)plot.ts(ts)Timets00 20 40 60 80 100 120 140100 300 500Timets0 20 40 60 80 100 120 1405.0 5.5 6.0 6.5n <- length(ts)tt <- seq(1,n)tt2 <- ttˆ2cos12 <- cos(2*pi*tt/12)sin12 <- sin(2*pi*tt/12)cos6 <- cos(2*pi*tt/6)sin6 <- sin(2*pi*tt/6)cos4 <- cos(2*pi*tt/4)sin4 <- sin(2*pi*tt/4)rg <- lm(ts ~ tt + tt2 + cos12 + sin12 + cos6 + sin6 + cos4 + sin4)summary(rg)#### Call:## lm(formula = ts ~ tt + tt2 + cos12 + sin12 + cos6 + sin6 + cos4 +## sin4)#### Residuals:1## Min 1Q Median 3Q Max## -0.173540 -0.035102 0.002494 0.038911 0.128827#### Coefficients:## Estimate Std. Error t value Pr(>|t|)## (Intercept) 4.736e+00 1.417e-02 334.223 < 2e-16 ***## tt 1.319e-02 4.508e-04 29.247 < 2e-16 ***## tt2 -2.148e-05 3.011e-06 -7.132 5.45e-11 ***## cos12 -1.415e-01 6.582e-03 -21.499 < 2e-16 ***## sin12 -4.927e-02 6.594e-03 -7.471 8.93e-12 ***## cos6 -2.275e-02 6.582e-03 -3.457 0.000732 ***## sin6 7.872e-02 6.584e-03 11.957 < 2e-16 ***## cos4 2.731e-02 6.582e-03 4.149 5.87e-05 ***## sin4 -8.706e-03 6.582e-03 -1.323 0.188144## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### Residual standard error: 0.05584 on 135 degrees of freedom## Multiple R-squared: 0.9849, Adjusted R-squared: 0.984## F-statistic: 1100 on 8 and 135 DF, p-value: < 2.2e-16ts2 <- residuals(rg)par(mfrow = c(1, 2))plot.ts(ts)lines(fitted.values(rg),col="red")plot.ts(ts2)Timets0 20 40 60 80 100 120 1405.0 5.5 6.0 6.5Timets20 20 40 60 80 100 120 140−0.15 −0.05 0.05par(mfrow = c(1, 2))acf(ts2,lag.max = 30)pacf(ts2,lag.max = 30)20 5 10 15 20 25 30−0.4 0.0 0.4 0.8LagACFSeries ts20 5 10 15 20 25 30−0.3 −0.1 0.1 0.3LagPartial ACFSeries ts2ts20 <- diff(ts,lag = 12)ts2a <- diff(ts20)par(mfrow = c(1, 2))plot.ts(ts20)plot.ts(ts2a)Timets200 20 40 60 80 100 1200.0 0.1 0.2 0.3Timets2a0 20 40 60 80 100 120−0.15 −0.05 0.05 0.15par(mfrow = c(1, 2))acf(ts2a,lag.max = 30)pacf(ts2a,lag.max = 30)0 5 10 15 20 25 30−0.4 0.0 0.4 0.8LagACFSeries ts2a0 5 10 15 20 25 30−0.3 −0.1 0.1LagPartial ACFSeries ts2a3Note that the resulting series exhibit strong (partial) correlations at some lags which are multiples of theperiod s. That is, bρ(s), bρ(2s), . . ., bα(s), bα(2s), . . .,are large.Data/Theory goals of the topicIntroduce (seasonal) models that can capture the above behavior. Discuss a few related issues, and alsoanother model with seasonal (periodic) features.Seasonal modelsMotivating seasonal models: Consider time series data for r years and 12 months:Year\Month 1 2 . . . 121 Y1Y2. . . Y122 Y13Y14. . . Y24...............r Y12(r−1)+1Y12(r−1)+2. . . Y12rThere arej= 1, . . . ,12 column seriesYj+12t,t= 0, . . . , r −1. Assume the same AR model for each:j = 1, . . . , 12,Yj+12t− Φ1Yj+12(t−1)− . . . − ΦPYj+12(t−P )= Uj+12twith {Uj+12t} ∼ WN(0, σ2). This is the same asΦ(B12)Yt= Utwith Φ(z) = 1 − Φ1z − . . . − ΦPzP. Terminology: between-year model.B Example: P = 1, Φ1= 0.7.To incorporate dependence between months, take another AR model for {Ut}:φ(B)Ut= Ztwith φ(z) = 1 − φ1z − . . . − φpzpand {Zt} ∼ WN(0, σ2Z). Combining the two models:φ(B)Φ(B12)Yt= Zt.What one obtains is known as a seasonal ARIMA model.Definition d, Dnon-negative integers:{Xt}is a seasonal ARIMA(p, d,0)×(P, D,0)stime series with periods if the differenced time series Yt= (I − B)d(I − Bs)DXtis a causal ARMA time series defined byφ(B)Φ(Bs)Yt= Ztwith φ(z) = 1 − φ1z − . . . − φpzp, Φ(z) = 1 − Φ1z − . . . − ΦPzPand {Zt} ∼ WN(0, σ2Z).MA parts can be added similarly.B Note the multiplicative nature of SARIMA models.Example: SARIMA(1, 0, 0) × (1, 0, 0)12with φ1= 0.7, Φ1= 0.6.library(astsa)set.seed(123)ts.e1 <- sarima.sim(ar=0.7, sar=0.6, S=12,n=400)ts.e1 <- as.numeric(ts.e1)4par(mfrow = c(1, 3))plot.ts(ts.e1)acf(ts.e1,lag.max = 30)pacf(ts.e1,lag.max = 30)Timets.e10 100 200 300 400−4 −2 0 2 40 5 10 15 20 25 300.0 0.2 0.4 0.6 0.8 1.0LagACFSeries ts.e10 5 10 15 20 25 30−0.2 0.0 0.2 0.4 0.6LagPartial ACFSeries ts.e1set.seed(1)ts.e2 <- arima.sim(list(order=c(13,0,0),ar=c(.7,rep(0,10),.6,-.7*.6)),n=400)par(mfrow = c(1, 1))plot.ts(ts.e2)Timets.e20 100 200 300 400−4 −2 0 2par(mfrow = c(1, 2))acf(ts.e2,lag.max = 30)pacf(ts.e2,lag.max = 30)library(forecast)50 5 10 15 20 25 300.0 0.4 0.8LagACFSeries ts.e20 5 10 15 20 25 30−0.2 0.2 0.4 0.6LagPartial ACFSeries ts.e2ar_acf<-ARMAacf(ar=c(.7,rep(0,10),.6,-.7*.6),lag.max=30,pacf=F)ar_pacf<-ARMAacf(ar=c(.7,rep(0,10),.6,-.7*.6),lag.max=30,pacf=T)par(mfrow = c(1, 2))plot(ar_acf, type='h')plot(ar_pacf, type='h')0 5 10 15 20 25 300.2 0.4 0.6 0.8 1.0Indexar_acf0 5 10 15 20 25 30−0.4 0.0 0.4Indexar_pacfH Homework 5 will ask you to explore similarly another SARIMA model.A number of R functions used earlier (model selection, fitting, forecasting, etc.) allow for seasonal models.library(forecast)set.seed(123)ts.e1 <- sarima.sim(ar=0.7, sar=0.6, S=12,n=400)auto.mod <- auto.arima(ts.e1, max.p = 2, max.q = 2, max.d = 0,start.p = 0, start.q = 0,max.P = 2, max.Q = 2, max.D = 0,start.P = 0, start.Q = 0,allowdrift = FALSE, ic = "aic")auto.mod## Series: ts.e1## ARIMA(1,0,0)(1,0,0)[12] with zero mean#### Coefficients:## ar1 sar16## 0.6061 0.5298## s.e. 0.0396 0.0426#### sigma^2 estimated as 0.9612: log likelihood=-560.87## AIC=1127.73 AICc=1127.79 BIC=1139.71sarima.model <- Arima(ts.e1,order=c(1,0,0),seasonal=c(1,0,0),method = "ML",include.mean = TRUE)sarima.model## Series: ts.e1## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean#### Coefficients:## ar1 sar1 mean## 0.6057 0.5295 0.0851## s.e. 0.0396 0.0426 0.2541#### sigma^2 estimated as 0.9634: log likelihood=-560.81## AIC=1129.62 AICc=1129.72 BIC=1145.59h = 12ts.e1.forecast <- forecast(auto.mod, h)plot(ts.e1.forecast)Forecasts from ARIMA(1,0,0)(1,0,0)[12] with zero mean0 5 10 15 20 25 30 35−4 −2 0 2 4Back to the time series motivating this topic.ts0 <- scan("airpass.txt")ts1 <- log(ts0)ts <- ts(data=ts1,frequency=12)auto.ts <- auto.arima(ts, max.p = 2, max.q = 2, max.d = 1,start.p = 0, start.q = 0,max.P = 2, max.Q = 2, max.D = 1,start.P = 0, start.Q = 0,allowdrift = FALSE, ic = "aic")auto.ts## Series: ts7## ARIMA(0,1,1)(0,1,1)[12]####