平稳时间序列分析：ARMA模型
ARMA模型拟合指令：
UDDI
arima(x,order=c(p,d,q),method=c("ML","CSS"),an)
其中:
x—带估计序列；
Order—指定模型阶数；其中，P为⾃回归阶数；d为差分阶数，若为平稳序列，则不需要差分，d=0；q为移动平均阶数。
method—估计⽅法。method =CSS-ML，系统默认的是条件最⼩⼆乘估计和极⼤似然估计的混合⽅法；method =ML，极⼤似然估计；method =CSS，条件最⼩⼆乘估计；
案例：选择合适的模型拟合1950-2008年我国油路及农村投递线路每年新增⾥程数序列。
（1）读取数据
file.choo
摄氏零度#弹出对话框，可以让我们选择⽂件的位置。R的运⾏结果：
#"C:\\Urs\\w\\Documents\\ "
d=read.table("C:\\Urs\\w\\Documents\\")#读⼊数据。
# 注意：读取txt格式⽂件时⽤“read.table”指令，header默认为FALSE，也即认为第⼀⾏就是数据；所以如果第⼀⾏不是数据⽽是标签，要改写为header=T/TU RE；读取csv格式⽂件时⽤“read.csv”，header默认为Ture，也即认为第⼀⾏是标题⾏；
licheng=ts(d,start=1950,end=2008,freq=1)#转化为时间序列数据。年度数据则freq=1; 半年度数据则freq=2; 季度数据则freq=4; ⽉度数据则freq=12; 周度数据则freq=52; ⽇度数据则freq=365;
（2）序列的平稳性检验：
时序图，⾃相关图，ADF检验
plot(licheng)#作时序图。⼤致围绕0值上下波动，且波动有界，该序列⼤致为平稳时间序列。
图1 时序图
acf(licheng,14)#作⾃相关图，⾃相关系数具有短期相关特性，因此进⼀步认为该序列为平稳时间序列。且具有拖尾特性。
pacf(licheng,14)#作偏⾃相关图，2阶截尾
图2 ⾃相关图
图3 偏⾃相关图
单位根检验的做法：
library(fBasics)
library(fUnitRoots)
#下载并安装“fUnitRoots”包，并⽤library函数进⾏调⽤。该程序包中的adfTest函数可以很便捷的进⾏单位根检验。
for(i in1:3)print(adfTest(licheng,lag=i,type="c"))
对“licheng”序列进⾏检验类型为类型2（⽆趋势，有常数项），延迟阶数分别取1,2,3的单位根检验。R的运⾏结果如下：
Title:（延迟阶数为1时的ADF检验）
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order:1
STATISTIC:
Dickey-Fuller:-6.9674
P VALUE:
0.01
Title:（延迟阶数为2阶时的ADF检验）
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order:2
STATISTIC:
小炒猪皮
Dickey-Fuller:-5.4821
P VALUE:
0.01
Title:（延迟阶数为3时的ADF检验）
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order:3
STATISTIC:
Dickey-Fuller:-4.2721
P VALUE:
0.01
adfTest检验，表明⽆论延迟阶数取1,2或3， 统计量对应的P值都为0.01，都拒绝原假设：序列⾮平稳。
表明序列为平稳序列。结合前⽂序列平稳性的时序图检验和⾃相关图检验，认为该序列为平稳序列。
（3）序列的⽩噪声检验：LB统计量
纯随机性检验,p值都⼩于5%，序列为⾮⽩噪声序列。
（4）平稳⾮⽩噪声序列的拟合：ARMA（p,q）模型
针对平稳⾮⽩噪声序列，可以运⽤ARMA类模型进⾏拟合。
根据ACF拖尾，PACF2阶截尾，建议采⽤AR（2）模型。
arima(licheng, order = c(2,0,0),method="ML")
⽤AR(2)模型拟合，这⾥选择参数估计⽅法method=“ML”。
若估计⽅法为条件最⼩⼆乘法CSS时则不显⽰AIC。
R运⾏结果：
Call:苦思冥想的意思
arima(x = licheng, order = c(2,0,0), method ="ML")
Coefficients:
ar1 ar2 intercept
0.7185-0.529411.0223
< 0.10830.1067 3.0906
sigma^2 estimated as 365.2: log likelihood =-258.23, aic =524.46
#注意：若要⽤中⼼化的AR(2)模型拟合，即假定不含截距项，则要补充语句an = F，即指令为：arima(licheng, order = c(2,0,0),method="ML", an = F)
运⾏该指令，得到不含截距项的AR（2）模型的AIC=532.03。
所以选择⾮中⼼化的模型。
根据⾮中⼼化模型的参数估计结果，得到均值为
=11.0223，⽽
=0.7185，=-0.5294，计算截距项
=8.9380。即此时的AR（2）模型为：
（5）平稳⾮⽩噪声序列拟合模型的检验
模型的检验1：模型的显著性检验。
a= arima(licheng, order = c(2,0,0),method="ML")
#把所有的拟合信息都存在⼀个叫“a”的⽂件（仓库）⾥。
r=a$residuals
#⽤r来保存残差
对残差序列进⾏纯随机性检验。
注意：fitdf的值等于（p+q），number of degrees of freedom to be subtracted if x is a ries of residuals。当检验的序列是残差序列时，需要加上命令fitdf，表⽰减去的⾃由度。
满足作文R运⾏结果：
data: r
X-squared =2.3416, df =4, p-value =0.6732
表明当检验阶数为6阶时，残差序列是⽩噪声序列。
R运⾏结果：
Box-Ljung test
data: r
X-squared =3.2655, df =10, p-value =0.9745
表明当检验阶数为12阶时，残差序列依然是⽩噪声序列。
综合检验阶数为6阶的检验结果，模型显著有效。
模型的检验2：参数的显著性检验。
注意：在时间序列中对参数显著性的要求与回归模型不同，我们更多的是考察模型整体的好坏，⽽不是参数。所以，R中的arima拟合结果中没有给出参数的t统计量值和p值。
如需计算参数的t统计量值和p值，利⽤下⾯的公式。
（1）参数 的t统计量值=
如在本例中，ar1即参数的估计值为0.7185，其估计值的标准误差为0.1083，则其t统计量=
。
其他的参数的t统计量的求法也是⼀样的。
（2）P值的求法。
卡文迪什以参数的t检验统计量的P值的计算为例。
p=pt(6.5420,df=56,lower.tail = F)*2
pt()为求t分布的p值的函数；
6.5420为t统计量的绝对值；
df为⾃由度=数据个数-参数个数，此处，数据有59个，参数3个，所以⾃由度为56；
lower.tail = F表⽰所求p值为P[T > t]即右侧概率，如不加⼊这个参数表⽰所求p值为P[T <=t]；
乘2表⽰p值是双侧的。
R运⾏的结果为：
p=pt(6.5420,df=56,lower.tail = F)*2
p
[1]6.94276e-09
由于该P值⼩于0.05，表明参数显著。
类似的，对ar(2)参数和常数项的显著性进⾏检验。表明模型通过了参数的显著性检验。
（6）基于拟合模型的平稳⾮⽩噪声序列的预测
接下来，基于该AR（2）模型进⾏未来5期的预测：
预测⽅法1：直接⽤predict函数进⾏预测。
predict(arima(licheng, order = c(2,0,0)), n.ahead =5)
#预测未来5期
R运⾏结果如下：
$pred #pred为未来5期的预测值
Time Series:
Start =2009
End =2013
Frequency =1
[1]9.4655156.2151948.39273611.67810812.885903
$ #表⽰未来5期的预测标准差
Time Series:
Start =2009
End =2013
Frequency =1
[1]19.1095323.5309223.5322624.6832625.22913
置信区间的求法：以95%的置信区间的求解为例
licheng.forecast= predict(arima(licheng, order = c(2,0,0)), n.ahead =5)
#将未来5期预测过程中的产⽣的结果保存在licheng.forecast变量中
U = licheng.forecast$pred +1.96* licheng.forecast$
L = licheng.forecast$pred -1.96* licheng.forecast$
#算出未来5期预测值的95%置信区间
R的运⾏结果：
U
Time Series:
Start =2009
End =2013
Frequency =1
[1]46.9202052.3358054.5159660.0572962.33500
L
Time Series:
Start =2009
名胜景点
蝴蝶b图片End =2013
Frequency =1
[1]-27.98917-39.90541-37.73049-36.70108-36.56319
所以，未来5期即2009-2013年的新增⾥程的预测结果为：
预测年份预测值95％置信区间
20099.465515（-27.98917,46.92020）
20106.215194（-39.90541,52.33580）
20118.392736（-37.73049,54.51596）
201211.678108（-36.70108,60.05729）
201312.885903（-36.56319,62.33500）
ts.plot(licheng,licheng.forecast$pred,col=1:2)
#作时序图，含预测。
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
#在时序图中作出95%置信区间，并指定线型为短划线。lty：连线的线型（1.实线，2.虚线，3.点线，4.点虚线，5.长虚线，6.双虚线）
