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

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2003 18 33 41 7 34 35 24 25 24 21 25 20
2004 22 31 40 29 25 21 22 54 31 25 26 35

plot(tsales)

start(tsales)
[1] 2003    1


end(tsales)
[1] 2004   12


frequency(tsales)

  [1] 12

 

tsales.subset<-window(tsales,start=c(2003,5),end=c(2004,6))
tsales.subset

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2003 34 35 24 25 24 21 25 20
2004 22 31 40 29 25 21

 

 

#简单移动平均
install.packages("forecast")
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 Average (k=3)",ylim=ylim)
plot(ma(Nile,7),main="Simple Moving Average (k=7)",ylim=ylim)
plot(ma(Nile,15),main="Simple Moving Average (k=15)",ylim=ylim)
par(opar)

 

#季节性分解
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="")

 

 

#单指数平滑
library(forecast)
fit<-ets(nhtemp,model="ANN")
fit

ETS(A,N,N)

Call:
ets(y = nhtemp, model = "ANN")

Smoothing parameters:
alpha = 0.182

Initial states:
l = 50.2759

sigma: 1.1263

AIC AICc BIC
265.9298 266.3584 272.2129

 

 

forecast(fit,1)

Point  Forecast Lo 80     Hi 80   Lo 95    Hi 95
1972   51.87045 50.42708 53.31382 49.66301 54.0779

 

plot(forecast(fit,1),xlab="Year",ylab=expression(paste("Temperature (",degreee*F,")",)),main="New Haven Annual Mean Temperature")

accuracy(fit)

             ME        RMSE     MAE       MPE       MAPE     MASE      ACF1
Training set 0.1460295 1.126268 0.8951331 0.2418693 1.748922 0.7512497 -0.00653111

 

ME: Mean Error

RMSE: Root Mean Squared Error

MAE: Mean Absolute Error

MPE: Mean Percentage Error

MAPE: Mean Absolute Percentage Error

MASE: Mean Absolute Scaled Error

ACF1: Autocorrelation of errors at lag 1.

 

#有水平项，斜率以及季节性的指数模型
library(forecast)
fit<-ets(log(AirPassengers),model="AAA")
fit

ETS(A,A,A)

Call:
ets(y = log(AirPassengers), model = "AAA")

Smoothing parameters:
alpha = 0.6534
beta = 1e-04
gamma = 1e-04

Initial states:
l = 4.8022
b = 0.01
s=-0.1047 -0.2186 -0.0761 0.0636 0.2083 0.217
0.1145 -0.011 -0.0111 0.0196 -0.1111 -0.0905

sigma: 0.0359

AIC AICc BIC
-208.3619 -203.5047 -157.8750

 

 

accuracy(fit)

pred<-forecast(fit,5)
pred

              ME            RMSE       MAE          MPE       MAPE      MASE      ACF1
Training set -0.0006710596 0.03592072 0.02773886 -0.01250262 0.508256 0.2291672 0.09431354

plot(pred,main="Forecast for Air Travel",ylab="Log(AirePassengers)",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

          mean     Lo 80    Lo 95    Hi 80    Hi 95
Jan 1961 447.4958 427.3626 417.0741 468.5774 480.1365
Feb 1961 442.7926 419.1001 407.0756 467.8245 481.6434
Mar 1961 509.6958 478.7268 463.1019 542.6682 560.9776
Apr 1961 499.2613 465.7258 448.8949 535.2116 555.2788
May 1961 504.3514 467.5503 449.1688 544.0491 566.3135

 

 

#ETS函数的自动指数预测
library(forecast)
fit<-ets(JohnsonJohnson)
fit

ETS(M,A,M)

Call:
ets(y = JohnsonJohnson)

Smoothing parameters:
alpha = 0.1481
beta = 0.0912
gamma = 0.4908

Initial states:
l = 0.6146
b = 0.005
s=0.692 1.2644 0.9666 1.077

sigma: 0.0889

AIC AICc BIC
166.6964 169.1289 188.5738

 

plot(forecast(fit),main="Johnson & Johnson Forecasts",ylab="Quarterly Earnings (Dollars)",xlab="Time",flty=2)

 

#序列的变换以及稳定性评估
library(forecast)
library(tseries)
plot(Nile)

ndiffs(Nile)

[1] 1

 

dNile<-diff(Nile)
plot(dNile)

#拟合ARIMA模型
library(forecast)
fit<-arima(Nile,order=c(0,1,1))
fit

accuracy(fit)

             ME RMSE MAE MPE MAPE MASE ACF1
Training set -11.9358 142.8071 112.1752 -3.574702 12.93594 0.841824 0.1153593

 

#模型评价
qqnorm(fit$residuals)
qqline(fit$residuals)
Box.test(fit$residuals,type="Ljung-Box")

Box-Ljung test

data: fit$residuals
X-squared = 1.3711, df = 1, p-value = 0.2416

#ARIMA 模型预测
forecast(fit,3)
plot(forecast(fit,3),xlab="Year",ylab="Annual Flow")

#ARIMA自动预测
library(forecast)
fit<-auto.arima(sunspots)
fit

Series: sunspots
ARIMA(2,1,2)

Coefficients:
ar1 ar2 ma1 ma2
1.3467 -0.3963 -1.7710 0.8103
s.e. 0.0303 0.0287 0.0205 0.0194

sigma^2 estimated as 243.8: log likelihood=-11745.5
AIC=23500.99 AICc=23501.01 BIC=23530.71

 

forecast(fit,3)

Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 1984 40.43784 20.42717 60.44850 9.834167 71.04150
Feb 1984 41.35311 18.26341 64.44281 6.040458 76.66576
Mar 1984 39.79670 15.23663 64.35677 2.235319 77.35808

accuracy(fit)

ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.02672716 15.60055 11.02575 NaN Inf 0.4775401 -0.01055012

 

 

小结

这章主要讲解了怎么用R语言来进行时间序列分析，例如：模型的建立，图表的绘制，以及未来趋势的预测。这类型的数据分析完全不在程序开发的范畴了，所有的分析都是基于数理统计，这应该就是现在的数据科学方向吧。