R语言：时间序列ARIMA模型(三)
题记：本文是个人的读书笔记，仅用于学习交流使用，本文将深入研究时间序列技术。
01
解决什么问题？
前面介绍了时间序列的经典分析方法，本节介绍另一个应用广泛的随机序列分析方法，常用的就是ARIMA模型自回归移动平均模型。
02
ARIMA模型
若是一个平稳时间序列，则该模型是一个具有p阶自回归和q阶移动平均的混合模型，用ARMA(p,q)表示。如果是非平稳时间序列，要进行差分，多一个参数d，表示为ARIMA(p,d,q)。
关于阶数的识别，可以采用第一节中的一般规则识别。本文采用自动识别模型阶数auto.arima()，之后采用arima()建模。时间序列模型的拟合一般是一下几个步骤：①模型阶数识别；②建模估计参数；③残差检验；④模型预测。
①模型阶数识别
数据集为数据集为某市2006-2008年大气SO2日平均浓度。
death <- read.csv('death_gz.csv',as.is = T)dth.ts <- ts(death$so2,frequency = 365,start = c(2006,1))plot.ts(dth.ts)auto.arima(dth.ts,stepwise=F,approximation=FALSE)
auto.arima()返回最优的拟合模型参数，这该函数默认设置p,d,q的值，根据aic，bic 遍历所有的模型，选择aic最小的为最优模型。计算得出阶数为(0,1,4)最合适。
②建模估计参数
dth.ts.arma <- arima(dth.ts,c(0,1,4))dth.ts.arma
返回的结果coef是模型参数估计值，sigma^2是拟合模型的残差方差值(537.1)，aic为10002.62.这个模型没有回报参数是否具有显著意义，可以手动计算，例如第一个t=（-0.259)/0.03 = -8.63，p <0.05，其他以此类推，所有的系数都有统计学意义。
③残差检验
par(mfrow=c(3,1),mar=c(2,3,0.5,1),mgp=c(1.8,0.5,0))plot(resid(dth.ts.arma))acf(resid(dth.ts.arma))pacf(resid(dth.ts.arma))Box.test(resid(dth.ts.arma)) ##p-value = 0.7125
结果解释
1.从残差序列图可看出(上图1)，在时间轴的1/3和2/3处比较大。
2.残差的自相关图(上图2)，已无显著自相关函数。从偏自相关函数图(上图3)，前2个偏自相关函数已无统计学意义，从滞后第3阶起的偏自相关函数中有一些具有统计学意义，不过他们的值都不大(小于0.1)。用Box.test()函数也可以做残差检验，P>0.05，残差没有相关性。
④模型预测
dth.ts.arma.for <- forecast(dth.ts.arma,365) ##向前365天预测值par(mar=c(4,4,1.4,0.5));plot(dth.ts.arma.for,col='green') #将原观测值和预测值画出来lines(dth.ts.arma.for$fitted,col='blue') ##将原时间点的回代值画出来U <- dth.ts.arma.for$lower[,2] ##95% ##预测值的上线L <- dth.ts.arma.for$upper[,2] ##95% 预测值的下线dth.ts.arma.for$mean[1:20] ###前20天的播报值
365天的预测值比较平坦，在S02的平均值上。在预测值上下两条向右延伸的平行线为95%和80%的可信区间的上下界值。
05
总结
本文作为时间序列的第三部分，主要介绍ARIMA模型的分析方法，详细介绍aima()函数的使用和模型的预测，模型的检验等。
06
参考文献
彭晓武等.R软件及其环境流行病学应用.中国环境出版社
Jonathan D.Cryer .时间序列分析及应用.机械工业出版社
Cory Leismester.精通机器学习：基于R.人民邮电出版社
作者介绍：医疗大数据统计分析师,擅长R语言。
欢迎各位在后台留言，恳请斧正！
更多阅读：