forked from kabacoff/RiA2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Ch15 Time series.R
128 lines (95 loc) · 3.29 KB
/
Ch15 Time series.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#-----------------------------------------#
# R in Action (2nd ed): Chapter 15 #
# Time series #
# requires forecast, tseries packages #
# install.packages("forecast", "tseries") #
#-----------------------------------------#
par(ask=TRUE)
# Listing 15.1 - Creating a time series object in R
sales <- c(18, 33, 41, 7, 34, 35, 24, 25, 24, 21, 25, 20,
22, 31, 40, 29, 25, 21, 22, 54, 31, 25, 26, 35)
tsales <- ts(sales, start=c(2003, 1), frequency=12)
tsales
plot(tsales)
start(tsales)
end(tsales)
frequency(tsales)
tsales.subset <- window(tsales, start=c(2003, 5), end=c(2004, 6))
tsales.subset
# Listing 15.2 - Simple moving averages
library(forecast)
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2))
ylim <- c(min(Nile), max(Nile))
plot(Nile, main="Raw time series")
plot(ma(Nile, 3), main="Simple Moving Averages (k=3)", ylim=ylim)
plot(ma(Nile, 7), main="Simple Moving Averages (k=7)", ylim=ylim)
plot(ma(Nile, 15), main="Simple Moving Averages (k=15)", ylim=ylim)
par(opar)
# Listing 15.3 - Seasonal decomposition using slt()
plot(AirPassengers)
lAirPassengers <- log(AirPassengers)
plot(lAirPassengers, ylab="log(AirPassengers)")
fit <- stl(lAirPassengers, s.window="period")
plot(fit)
fit$time.series
exp(fit$time.series)
par(mfrow=c(2,1))
library(forecast)
monthplot(AirPassengers, xlab="", ylab="")
seasonplot(AirPassengers, year.labels="TRUE", main="")
par(opar)
# Listing 15.4 - Simple exponential smoothing
library(forecast)
fit <- HoltWinters(nhtemp, beta=FALSE, gamma=FALSE)
fit
forecast(fit, 1)
plot(forecast(fit, 1), xlab="Year",
ylab=expression(paste("Temperature (", degree*F,")",)),
main="New Haven Annual Mean Temperature")
accuracy(fit)
# Listing 15.5 - Exponential smoothing with level, slope, and seasonal components
fit <- HoltWinters(log(AirPassengers))
fit
accuracy(fit)
pred <- forecast(fit, 5)
pred
plot(pred, main="Forecast for Air Travel",
ylab="Log(AirPassengers)", xlab="Time")
pred$mean <- exp(pred$mean)
pred$lower <- exp(pred$lower)
pred$upper <- exp(pred$upper)
p <- cbind(pred$mean, pred$lower, pred$upper)
dimnames(p)[[2]] <- c("mean", "Lo 80", "Lo 95", "Hi 80", "Hi 95")
p
# Listing 15.6 - Automatic exponential forecasting with ets()
library(forecast)
fit <- ets(JohnsonJohnson)
fit
plot(forecast(fit), main="Johnson and Johnson Forecasts",
ylab="Quarterly Earnings (Dollars)", xlab="Time")
# Listing 15.7 - Transforming the time series and assessing stationarity
library(forecast)
library(tseries)
plot(Nile)
ndiffs(Nile)
dNile <- diff(Nile)
plot(dNile)
adf.test(dNile)
# Listing 15.8 - Fit an ARIMA model
fit <- arima(Nile, order=c(0,1,1))
fit
accuracy(fit)
# Listing 15.9 - Evaluating the model fit
qqnorm(fit$residuals)
qqline(fit$residuals)
Box.test(fit$residuals, type="Ljung-Box")
# Listing 15.10 - Forecasting with an ARIMA model
forecast(fit, 3)
plot(forecast(fit, 3), xlab="Year", ylab="Annual Flow")
# Listing 15.11 - Automated ARIMA forecasting
library(forecast)
fit <- auto.arima(sunspots)
fit
forecast(fit, 3)
accuracy(fit)