Time Series–Prediction model

Reference book: R in Action
For further study, the following books are recommended:

Time Series: A tutorial book
A Little Book of R for Time Series: Another tutorial book with R code
Forecast: Principles and Practice: Highly recommended

指数预测模型: ets function from forecast package

pros: 短期预测性能良好

library(forecast)

Single exponential model: 拟合常数水平项和随机项

观测值: Y_t = level + irregular_t

fitANN=ets(nhtemp, model = "ANN")
fitANN

## 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

预测值(1-step ahead forecast): Y_t+1 = \(\sum_0^t\){c_i*Y_t-i}
其中 c_i = \(\alpha * (1 - \alpha)^i\); \(0 \leq \alpha \leq 1\); \(\sum_0^t{c_i} = 1\)

forecast(fitANN, 1)

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

时序值, 预测值, 预测区间

accuracy(fitANN)

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

准确性度量

Double exponential model (Holt exponential smoothing): 拟合水平项和趋势项

观测值: Y_t = level + slope * t + irregular_t

fitAAN=ets(nhtemp, model = "AAN")

Trible exponential model (Holt-Winters exponential smoothing): 拟合水平项, 趋势项以及季节效应

观测值: Y_t = level + slope * t + s_t + irregular_t

fitAAA=ets(log(AirPassengers), model = "AAA")
fitAAA

## 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

用原始尺度预测

pred=forecast(fitAAA, 5) ##预测未来5个月

p=cbind(exp(pred$mean), exp(pred$lower), exp(pred$upper))
colnames(p)=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

ARIMA预测模型: arima function from forecast package

!!要求序列具有平稳性

建立ARIMA模型的步骤

1.变换序列的值, 确保序列方差为常数

常用变换: 对数变换, Box-Cox变换

2.通过d次差分移除趋势项

library(forecast)
ndiffs(ts) # Find the optimal d
dts=diff(ts, differences = d) # 对序列进行差分

3.通过ADF(Augmented Dickey-Fuller)统计检验来验证平稳性检验, 由tseries package中的adf.test()实现

4.确定p, q的值, 使得序列中的每个观测值可由过去p个观测值和q个残差的线性组合表示

Acf(dts) # 绘制ACF图
Pacf(dts) # 绘制PACF图

选择ARIMA模型的方法
模型	ACF自相关	PACF偏自相关
ARIMA(p, d, 0)	逐渐减小到零	p阶后减小到零
ARIMA(0, d, q)	q阶后减小到零	逐渐减小到零
ARIMA(p, d, q)	逐渐减小到零	逐渐减小到零

5.使用arima(ts, order=c(p, d, q))拟合模型

6.根据AIC或者准确性度量来选择最佳模型

7.模型的独立正态分布检验, 由Q-Q图及Box.test()实现

An example of Nile

library(forecast)
library(tseries)

plot(Nile) # 方差稳定但存在趋势项

ndiffs(Nile)

## [1] 1

dNile=diff(Nile) # 一次差分
plot(dNile)

adf.test(dNile)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dNile
## Dickey-Fuller = -6.5924, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Acf(dNile)

Pacf(dNile)

由上两图可知, q取1, p取0

fit=arima(Nile, order=c(0, 1, 1))
fit

## 
## Call:
## arima(x = Nile, order = c(0, 1, 1))
## 
## Coefficients:
##           ma1
##       -0.7329
## s.e.   0.1143
## 
## sigma^2 estimated as 20600:  log likelihood = -632.55,  aic = 1269.09

accuracy(fit) # 准确性度量

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

qqnorm(fit$residuals)
qqline(fit$residuals) # Q-Q Plot

Box.test(fit$residuals, type = "Ljung-Box")

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

p-value不显著, 故认为残差的自相关系数为0

forecast(fit, 3) # 预测后三个月

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1971       798.3673 614.4307 982.3040 517.0605 1079.674
## 1972       798.3673 607.9845 988.7502 507.2019 1089.533
## 1973       798.3673 601.7495 994.9851 497.6663 1099.068

plot(forecast(fit, 3), xlab="Year", ylab="Annual Flow")