r语言：时间序列ARMA基础学习

1、r语言：时间序列ARMA基础学习26 二月, 2013,oldlee11,R语言与数据挖掘,时间序列分析,?010203040506070809101112131415161718192021222324252627282930313233343536# 术语 #白噪音：其均值=0，并且独立分布的（同时间无关）#稳定时间序列：任意j-i时间段的序列：其均值相等#自相关系数acf图：研究yt同yt-l序列之间的相关性# 在纯的ma(q)序列下，acf图形表现为q+1以后的自相关系数约为0（虚线内）#偏相关系数pacf图：在yt同yt-l之间的序列固定的情况下，研究研究yt同yt-l序列之间的相关

2、性# 在纯的ar(p)序列下，pacf图形表现为p+1以后的偏相关系数约为0（虚线内）#扩展相关系数图eacf：如果yt不是纯的ar或ma，而是arma（混合体），无法通过acf确定q，也不能通过pacf确认p，需要通过eacf确认p和q# 模拟产生ma ar arma 序列 # MA 时间序列的模拟试验：产生一个ma时间序列y.ma-function(a1,a2,a3=0,a4=0,num=200,pic=TRUE)#MA滑动平均时间序列的模拟（也可以使用filter函数）e-rnorm(num,0,1)#模拟白噪声,均值=0result-0result1-e1result2-e2-a1*e

3、1result3-e3-a1*e2-a2*e1result4-e4-a1*e3-a2*e2-a3*e1for(t in 5:num) resultt-et-a1*et-1-a2*et-2-a3*et-3-a4*et-4 #构造一个ma型时间序列if(pic=TRUE)#画图形dev.new()ts.plot(result,main=paste(y.mat=et-,a1,*et-1-,a2,*et-2-,a3,*et-3-,a4,*et-4的时间序列散点图)dev.new()lag.plot(result, 9, do.lines=FALSE)dev.new()par(mfrow=c(2,1)a

4、cf(result, 30,main=paste(y.ma自相关图,y.mat=et-,a1,*et-1-,a2,*et-2-,a3,*et-3-,a4,*et-4)pacf(result, 30,main=paste(y.ma偏自相关图,y.mat=et-,a1,*et-1-,a2,*et-2-,a3,*et-3-,a4,*et-4)resulty.ma-y.ma(0.92,0.65)结果一：绘制散点图结果二：绘制出yt同yt-1(延迟1)、yt-2(延迟2).yt-9(延迟9)的2维散点图，用以观察yt同yt-i的相关性（i=19）结果三：绘制自相关和偏自相关图（在虚线外的表示有相关性）可

5、以看到：1）在纯的ma(q)序列下，acf图形表现为q+1以后的自相关系数约为0（虚线内）2）在纯的ma(q)序列下，pacf则不规则的会大于1/-1可以通过acf图确定q的数值?01020304050607080910111213141516171819202122# AR 时间序列的模拟试验：产生一个ar时间序列y.ar-function(b1,b2,b3=0,b4=0,num=200,pic=TRUE)#AR自回归型时间序列的模拟（也可以使用filter函数）e-rnorm(num,0,1)#模拟白噪声,均值=0result-0result1-e1result2-e2+b1*result

6、1result3-e3+b1*result2+b2*result1result4-e4+b1*result3+b2*result2+b3*result1for(t in 5:num) resultt-et+b1*resultt-1+b2*resultt-2+b3*resultt-3+b4*resultt-4 #构造一个ar型时间序列if(pic=TRUE)#画图形dev.new()ts.plot(result,main=paste(y.art=et+,b1,*y.art-1+,b2,*y.art-2+,b3,*y.art-3+,b4,*y.art-4的时间序列散点图)dev.new()lag.

7、plot(result, 9, do.lines=FALSE)dev.new()par(mfrow=c(2,1)acf(result, 30,main=paste(y.ar自相关图,y.art=et+,b1,*y.art-1+,b2,*y.art-2+,b3,*y.art-3+,b4,*y.art-4)pacf(result, 30,main=paste(y.ar偏自相关图,y.art=et+,b1,*y.art-1+,b2,*y.art-2+,b3,*y.art-3+,b4,*y.art-4)resulty.ar-y.ar(0.92,-0.65)结果一：绘制散点图结果二：绘制出yt同yt-1

8、(延迟1)、yt-2(延迟2).yt-9(延迟9)的2维散点图，用以观察yt同yt-i的相关性（i=19）结果三：绘制自相关和偏自相关图（在虚线外的表示有相关性）1）在纯的ar(q)序列下，pacf图形表现为q+1以后的自相关系数约为0（虚线内）2）在纯的ar(q)序列下，acf则不规则的会大于1/-1可以通过pacf图确定p的数值?0102030405060708091011121314151617181920# ARMA 时间序列的模拟试验：产生一个arma时间序列library(TSA)y.arma-function(a1,a2,a3=0,a4=0,b1,b2,b3=0,b4=0,num

9、=200,pic=TRUE)result-y.ma(a1=a1,a2=a2,a3=a3,a4=a4,pic=F,num=num)+y.ar(b1=b1,b2=b2,b3=b3,b4=b4,pic=F,num=num)#产生自回归滑动平均时间序列ARIMAexp.str=q+1的自相关系数全部约小于0 | 逐步将接近0# AR(p) 逐步将接近0 | lag=p+1的自相关系数全部约小于0 # ARMA(p,q) 逐步将接近0 | 逐步将接近0# step2 序列的类型确定后：求p q# 如果是MA(q) 类型的，看acf图，看前几个不为0# 如果是AR(p) 类型的，看pacf图，看前几个不为

10、0# 如果是ARMA(p,q) 类型的，看eacf图，看x三角形的顶端#代码：data是向量数据library(TSA)plot.cf-function(data)dev.new()ts.plot(data,main=时间序列散点图)dev.new()lag.plot(data, 9, do.lines=FALSE)dev.new()par(mfrow=c(2,1)acf(data, 30,main=自相关图)pacf(data, 30,main=偏自相关图)eacf(data)data-y.arma(a1=0.92,a2=0.65,b1=0.92,b2=-0.65,num=100,pic=F

11、)plot.cf(data)# 产生ariam模型 # p AR的参数# q MA的参数# d 移动取差的阶数# arima(data.ts,order=c(p,d,q)#产生模型。arima只处理时间序列，不处理向量等p=2;d=0;q=2data-y.arma(a1=0.92,a2=0.65,b1=0.92,b2=-0.65,num=100,pic=F)data.ts-ts(data,freq=1,star=1)#把向量化为时间序列sol-arima(data.ts,order=c(p,d,q)#产生模型。arima只处理时间序列，不处理向量等dev.new()plot(sol,n.ahe

12、ad=30)#预测30个点# 检验预测模型的性能 # 序列的正态分布验证函数norm.test-function(input.data,alpha=0.05,pic=TRUE)if(pic=TRUE)#画图形dev.new()par(mfrow=c(2,1)qqnorm(input.data,main=qq图)qqline(input.data)hist(input.data,freq=F,main=直方图和密度估计曲线)lines(density(input.data),col=blue)#密度估计曲线x-c(round(min(input.data):round(max(input.dat

13、a)lines(x,dnorm(x,mean(input.data),sd(input.data),col=red)#正态分布曲线solalpha)print(paste(success:服从正态分布，p.value=,sol$p.value,alpha)elseprint(paste(error:不服从正态分布，p.value=,sol$p.value,=,alpha) sol#残差分析#方案一（自己编写的，不一定好）# ariam.sol由ariam()产生的模型aya.res-function(ariam.sol)#step1 ：检验正太性，如果符合正太，表示均值相同：稳定的dev.ne

14、w()norm.test(ariam.sol$residuals)#画出QQ图,并给出检验正态分布的测试结果par(mfrow=c(3,1)#step2 ：查看时序，如果有趋势变化，则mean不为相同，即不是稳定的plot(ariam.sol$residuals,main=残差序列图,type=o)#画出时序图abline(h=0)#step3 ：查看acf图，如果全部在虚线内，则表示残差和时间无关系acf(ariam.sol$residuals, 30,main=残差的自相关图:如果所有lag的系数均处于虚线内，则原始模型较好)#step4 ：使用Ljung来进行白噪声验证：表示均值全部为0

15、，即残差很小test.p-0for(i in 1:20)test.pi-Box.test(sol$residuals,i,fitdf=2,type=Ljung)$p.value#选取残差中lag=1-i的自相关系数来建立Ljung量。print(test.p)plot(c(NA,NA,test.p-c(1:2),main=残差的Ljung检验p.value:如果lag=3-以后的所有点均大于0.05，则是白噪声,ylim=c(0,1)abline(h=0.05,col=red) aya.res(sol)#方案二(好，书上的)#step1 ：检验正太性，如果符合正太，表示均值相同：稳定的dev.

16、new()norm.test(sol$residuals)#画出QQ图,并给出检验正态分布的测试结果dev.new()#使用对arima产生的模型，使用tsdiag()函数，直接画出：残差的时序图，残差的acf图，Ljung产生的p.value（和0.05比较）tsdiag(sol,gof=15,omit.initial=F)# sol是由ariam()产生的模型# 使用ariam模型进行预测 # arima.pred函数：# ariam.sol由ariam()产生的模型# nahead预测的点数# alpha计算预测区间时的置信度：1-alpha# 输出的数据为数据框：$pred是预测的估值# $max是预测的估值的最大值（1-alpha置信度，双侧估计）# $min是预测的估值的最小值（1-alpha置信度，双侧估计）arima.pred-function(ariam.sol,nahead=5,alpha=0.05)pred-predict(ariam.sol,n.ahead=nahead)#产生nahead个预测数据，包括预测的平均值$pred.和预测标准误$sek-qnorm(1-alpha/2),0,1)pred.data-data.frame(pred=pred$pred,max