吴裕雄--天生自然 R语言开发学习：时间序列（续二）-CSDN博客

#-----------------------------------------#
# 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)
转载于:https://www.cnblogs.com/tszr/p/11176317.html