时间序列分析:ARIMA 建模示例
1. 一般流程
对时间序列的分析和建模一般有一下几个步骤:
- 查看时间序列的图像:
1.1. 如果具有趋势,可能需要对序列进行差分;
1.2. 如果方差随时间变化,可能需要对序列进行转换,如 $\log$ 变换。先进行 $\log$ 变换再做差分的变换称为 Log-Return。 - 查看 ACF 图像,判断滑动平均的阶数。
- 查看 PACF 图像,判断自回归阶数。
- 选择若干组参数进行建模,可以使用如下指标来选择模型:
- SSE:评价拟合质量;
- AIC:综合评价拟合质量和模型复杂度;
- Ljung-Box Q 统计量:判断残差是否有自相关。
- 进行估计。
2. 一个例子
数据可以在这里下载到。
2.1. 绘制图像
首先读取数据并绘制图像如图 1 所示:
birth.data <- read.csv("./data/daily-total-female-births-CA.csv") birth.num <- birth.data$births birth.date <- as.Date(birth.data$date, "%Y-%m-%d") plot(birth.num ~ birth.date, type='l')
可以看到数据似乎有一些趋势。对序列进行 Ljung-Box 检验
Box.test(birth.num, lag=log(length(birth.num)))
结果为
Box-Pierce test data: birth.num X-squared = 36.391, df = 5.8999, p-value = 2.088e-06
可见 $p = 2.088 \times 10^{-6} < 0.05$,拒绝零假设,认为序列存在自相关。
2.2. 进行差分
为了移除趋势,对序列进行差分,绘制图像如图 2 所示:
birth.num.diff <- diff(birth.num) plot(birth.num.diff ~ birth.date[1:364], type='l')
可见趋势基本消失了。再进行 Ljung-Box 检验
Box.test(birth.num.diff, lag=log(length(birth.num.diff)))
结果为
Box-Pierce test data: birth.num.diff X-squared = 78.094, df = 5.8972, p-value = 7.661e-15
可见 $p = 7.661 \times 10^{-15} < 0.05$,拒绝零假设,序列仍存在自相关。
2.3. 查看 ACF 和 PACF
绘制 ACF 和 PACF 图像如图 3、图 4 所示:
acf(birth.num.diff, 50) pacf(birth.num.diff, 50)
可见 ACF 在滞后为 $1$ 处截断,PACF 像是在衰减,也像是在延迟为 $7$ 处截断。
2.4. 模型比较
结合分析,使用以下 $4$ 组模型进行拟合:
params <- list(c(0, 1, 1), c(0, 1, 2), c(7, 1, 1), c(7, 1, 2)) results <- data.frame() for (i in 1:length(params)) { param <- params[[i]] model <- arima(birth.num, order=param) test <- Box.test(model$residuals, lag=log(length(model$residuals))) result = data.frame(model=paste("ARIMA(", param[1], ",", param[2], ",", param[3], ")", sep=""), AIC=model$aic, SSE=sum(resid(model)^2), P.value=test$p.value) results <- rbind(results, result) }
结果为
model AIC SSE P.value ARIMA(0,1,1) 2462.221 18148.46 0.5333604 ARIMA(0,1,2) 2459.571 17914.65 0.9859227 ARIMA(7,1,1) 2464.883 17584.39 0.9999899 ARIMA(7,1,2) 2466.666 17574.06 0.9999929
选择 AIC 最小的模型,即 $\mathrm{ARIMA}(0,1,2)$。
2.5. 最终模型
最后使用 astsa
包的 sarima
函数进行拟合:
library(astsa) sarima(birth.num, 0, 1, 2, 0, 0, 0)
得到系数为
Coefficients: ma1 ma2 constant -0.8511 -0.1113 0.015 s.e. 0.0496 0.0502 0.015
故最终的模型为
\begin{equation}
(1 – B)X_t = 0.015 + e_t – 0.8511 e_{t-1} – 0.1113 e_{t-2}
\end{equation}
即
\begin{equation}
X_t = X_{t-1} + 0.015 + e_t – 0.8511 e_{t-1} – 0.1113 e_{t-2}
\end{equation}