#################Prognose in einem Interventionsmodell####################### #################Erstellen Sie eine Prognose für die Zeitreie resex.txt für 12 weitere Zeitpunkte #installation of residential telephone extensions #see p. 43-44 in Diagnostic checks in time series / Li, Wai Keung, 1953 #zwei Ausreißer bei 83 und 84 (Nov und Dec 1972) re <- scan("//ad.uni-hamburg.de/basis/mit/bm01/BAB5975/Documents/Lehre/ZRA/Anwendungen/intervention/resex/resex.txt") ts.plot(re) #Ausreißer bei 83!!! #Wir verwenden eine Inteventionsmodell re.ts <- ts(re, start=c(1966,1), frequency=12) ts.plot(re.ts, xlab="Jahr", ylab="Inst. privater Tel.anschlüsse (in Tsd.)") ts.plot(re.ts, xlab="Jahr", ylab="Inst. privater Tel.anschlüsse (in Tsd.)", type = "o") #################Modell für den Abschnitt 1:82 der Zeireihe bestimmen re.log <- ts(log(re)) ts.plot(re.log) ts.plot(re.log[1:82]) par(mfrow = c(1,2)) acf(re.log[1:82], lag.max = 20, xlab=expression(tau), main = "ACF mit Bartlett-Grenzen", ci.type = "ma") pacf(re.log[1:82], lag.max = 20, xlab=expression(tau), main = "PACF mit White-Noise-Grenzen") source("//ad.uni-hamburg.de/basis/mit/bm01/BAB5975/Documents/Lehre/ZRA/Anwendungen/tsutil.r") vartable(re.log[1:82], 12) # eine einfache Dif. und eine saisonale Dif. re.df <- diff(diff(re.log[1:82], 12),1) par(mfrow=c(1,1)) ts.plot(re.df) par(mfrow = c(1,2)) acf(re.df, lag.max = 20, xlab=expression(tau), main = "ACF mit Bartlett-Grenzen", ci.type = "ma") pacf(re.df, lag.max = 20, xlab=expression(tau), main = "PACF mit White-Noise-Grenzen") out <- arima(re.log[1:82], order = c(1,1,1),seasonal = list(order = c(0,1,0),period = 12), include.mean = T) out #=> alfa1 nicht signifikant => MA[1] out <- arima(re.log[1:82], order = c(0,1,1),seasonal = list(order = c(0,1,0), period = 12), include.mean = T) out par(mfrow=c(1,1)) plot(out$residuals) #1. QQ-Diagramm qqnorm(out$residuals) abline(0,sd(out$residuals)) #2. ACF und PACF der Residuen par(mfrow = c(1,2)) acf(out$residuals, lag.max = 20, xlab=expression(tau), main = "ACF mit Bartlett-Grenzen", ci.type = "ma") pacf(out$residuals, lag.max = 20, xlab=expression(tau), main = "PACF mit White-Noise-Grenzen") #=> SARIMA[0,1,1][1,1,0](12) out <- arima(re.log[1:82], order = c(0,1,1),seasonal = list(order = c(1,1,0), period = 12), include.mean = T) out par(mfrow = c(1,2)) acf(out$residuals, lag.max = 20, xlab=expression(tau), main = "ACF mit Bartlett-Grenzen", ci.type = "ma") pacf(out$residuals, lag.max = 20, xlab=expression(tau), main = "PACF mit White-Noise-Grenzen") #=> OK ###################### Intervention festlegen it <- rep(0, 89) it[83]<- 1 ######################Parameterschätzung ts <- c(re.log [1:82], NA, re.log [84:89]) m <- mean(ts, na.rm =TRUE) ts.zr <- ts-m ts.zr[83] <- 0 out <- arima(ts.zr, order = c(0,1,1), seasonal = list(order = c(1,1,0), period = 12), include.mean = F, xreg = it) out # Wir schätzen den Wert -.3931 + 9.75 = 9.357 für re.log[83] ########################Prognose mit Arima mit Regressor pout <- predict(out, n.ahead = 12, newxreg = rep(0,12)) ts.plot(re.log, type = "o", xlim = c(0,100)) lines(pout$pred+m, col = "red", type = "o")