r语言时间序列ARMA基础学习.docx
《r语言时间序列ARMA基础学习.docx》由会员分享，可在线阅读，更多相关《r语言时间序列ARMA基础学习.docx（18页珍藏版）》请在冰豆网上搜索。
r语言时间序列ARMA基础学习
r语言：
时间序列ARMA基础学习
26二月,2013, oldlee11, R语言与数据挖掘, 时间序列分析,
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
######################################### 术语#################################
#白噪音：
其均值=0，并且独立分布的（同时间无关）
#稳定时间序列：
任意j--i时间段的序列：
其均值相等
#自相关系数acf图：
研究y[t]同y[t-l]序列之间的相关性
# 在纯的ma（q）序列下，acf图形表现为q+1以后的自相关系数约为0（虚线内）
#偏相关系数pacf图：
在y[t]同y[t-l]之间的序列固定的情况下，研究研究y[t]同y[t-l]序列之间的相关性
# 在纯的ar（p）序列下，pacf图形表现为p+1以后的偏相关系数约为0（虚线内）
#扩展相关系数图eacf：
如果y[t]不是纯的ar或ma，而是arma（混合体），无法通过acf确定q，也不能通过pacf确认p，需要通过eacf确认p和q
################################################################################
###############################################
#### 模拟产生maararma序列 ####
###############################################
#####MA时间序列的模拟试验：
产生一个ma时间序列
y.ma<-function（a1,a2,a3=0,a4=0,num=200,pic=TRUE）{#MA滑动平均时间序列的模拟（也可以使用filter函数）
e<-rnorm（num,0,1）#模拟白噪声,均值=0
result<-0
result[1]<-e[1]
result[2]<-e[2]-a1*e[1]
result[3]<-e[3]-a1*e[2]-a2*e[1]
result[4]<-e[4]-a1*e[3]-a2*e[2]-a3*e[1]
for（tin5:
num）{result[t]<-e[t]-a1*e[t-1]-a2*e[t-2]-a3*e[t-3]-a4*e[t-4]}#构造一个ma型时间序列
if（pic==TRUE）{#画图形
dev.new（）
ts.plot（result,main=paste（"y.ma[t]=e[t]-",a1,"*e[t-1]-",a2,"*e[t-2]-",a3,"*e[t-3]-",a4,"*e[t-4]的时间序列散点图"））
dev.new（）
lag.plot（result,9,do.lines=FALSE）
dev.new（）
par（mfrow=c（2,1））
acf（result,30,main=paste（"y.ma自相关图,y.ma[t]=e[t]-",a1,"*e[t-1]-",a2,"*e[t-2]-",a3,"*e[t-3]-",a4,"*e[t-4]"））
pacf（result,30,main=paste（"y.ma偏自相关图,y.ma[t]=e[t]-",a1,"*e[t-1]-",a2,"*e[t-2]-",a3,"*e[t-3]-",a4,"*e[t-4]"））
}
result
}
y.ma<-y.ma（0.92,0.65）
结果一：
绘制散点图
结果二：
绘制出y[t]同y[t-1]（延迟1）、y[t-2]（延迟2）……..y[t-9]（延迟9）的2维散点图，用以观察
y[t]同y[t-i]的相关性（i=1–9）
结果三：
绘制自相关和偏自相关图（在虚线外的表示有相关性）
可以看到：
1）在纯的ma（q）序列下，acf图形表现为q+1以后的自相关系数约为0（虚线内）
2）在纯的ma（q）序列下，pacf则不规则的会大于1/-1
可以通过acf图确定q的数值
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
#####AR时间序列的模拟试验：
产生一个ar时间序列
y.ar<-function（b1,b2,b3=0,b4=0,num=200,pic=TRUE）{#AR自回归型时间序列的模拟（也可以使用filter函数）
e<-rnorm（num,0,1）#模拟白噪声,均值=0
result<-0
result[1]<-e[1]
result[2]<-e[2]+b1*result[1]
result[3]<-e[3]+b1*result[2]+b2*result[1]
result[4]<-e[4]+b1*result[3]+b2*result[2]+b3*result[1]
for（tin5:
num）{result[t]<-e[t]+b1*result[t-1]+b2*result[t-2]+b3*result[t-3]+b4*result[t-4]}#构造一个ar型时间序列
if（pic==TRUE）{#画图形
dev.new（）
ts.plot（result,main=paste（"y.ar[t]=e[t]+",b1,"*y.ar[t-1]+",b2,"*y.ar[t-2]+",b3,"*y.ar[t-3]+",b4,"*y.ar[t-4]的时间序列散点图"））
dev.new（）
lag.plot（result,9,do.lines=FALSE）
dev.new（）
par（mfrow=c（2,1））
acf（result,30,main=paste（"y.ar自相关图,y.ar[t]=e[t]+",b1,"*y.ar[t-1]+",b2,"*y.ar[t-2]+",b3,"*y.ar[t-3]+",b4,"*y.ar[t-4]"））
pacf（result,30,main=paste（"y.ar偏自相关图,y.ar[t]=e[t]+",b1,"*y.ar[t-1]+",b2,"*y.ar[t-2]+",b3,"*y.ar[t-3]+",b4,"*y.ar[t-4]"））
}
result
}
y.ar<-y.ar（0.92,-0.65）
结果一：
绘制散点图
结果二：
绘制出y[t]同y[t-1]（延迟1）、y[t-2]（延迟2）……..y[t-9]（延迟9）的2维散点图，用以观察
y[t]同y[t-i]的相关性（i=1–9）
结果三：
绘制自相关和偏自相关图（在虚线外的表示有相关性）
1）在纯的ar（q）序列下，pacf图形表现为q+1以后的自相关系数约为0（虚线内）
2）在纯的ar（q）序列下，acf则不规则的会大于1/-1
可以通过pacf图确定p的数值
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
#####ARMA时间序列的模拟试验：
产生一个arma时间序列
library（TSA）
y.arma<-function（a1,a2,a3=0,a4=0,b1,b2,b3=0,b4=0,num=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）#产生自回归滑动平均时间序列ARIMA
exp.str<-paste（"y.arma[t]=e[t]-",a1,"*e[t-1]-",a2,"*e[t-2]-",a3,"*e[t-3]-",a4,"*e[t-4]+",b1,"*y.arma[t-1]+",b2,"*y.arma[t-2]+",b3,"*y.arma[t-3]+",b4,"*y.arma[t-4]"）
if（pic==TRUE）{#画图形
dev.new（）
ts.plot（result,main=paste（exp.str,"的时间序列散点图"））
dev.new（）
lag.plot（result,9,do.lines=FALSE）
dev.new（）
par（mfrow=c（2,1））
acf（result,30,main=paste（exp.str,"的自相关图"））
pacf（result,30,main=paste（exp.str,"的偏自相关图"））
}
print（paste（exp.str,"的扩展相关图。
用于确认p和q的数值"））
eacf（result）#观察eacf图，确定p和q,###需要加载：
library（TSA）
result
}
y.arma（a1=0.92,a2=0.65,b1=0.92,b2=-0.65,num=100）
结果一：
绘制散点图
结果二：
绘制出y[t]同y[t-1]（延迟1）、y[t-2]（延迟2）……..y[t-9]（延迟9）的2维散点图，用以观察
y[t]同y[t-i]的相关性（i=1–9）
结果三：
绘制自相关和偏自相关图（在虚线外的表示有相关性）
1）在arma（q）序列下，acf和pacf图形全市不规则的
所以在arma下acf和pacf并不能用于确定p和q，此时用eacf来大致确定
AR/MA
012345678910111213
0ooxxooooooo x o o
1xoooooooooo x o o
2xoooooooooo x o o
3xoooooooooo o o o
4xoooooooooo o o o
5xoooooooooo o o o
6ooooooooooo o o o
7xoooooooooo o o o
如何使用该改图，我会以后补充（基本就是查看0的倒三角）
以上介绍了arma的一般p和q情况，下面来阐述分析一个arma的一般过程
001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
052
053
054
055
056
057
058
059
060
061
062
063
064
065
066
067
068
069
070
071
072
073
074
075
076
077
078
079
080
081
082
083
084
085
086
087
088
089
090
091
092
093
094
095
096
097
098
099
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#############################################################
#### 使用差分log对不稳定序列进行“稳定化”并计算d ####
#############################################################
#取差分，使其稳定化：
#代码：
# diff（data）#1阶差分，d=1
# diff（data,2）#2阶差分，d=2
#取log再差点，使其稳定化：
#代码：
# diff（log（data））#1阶差分，d=1
# diff（log（data）,2）#2阶差分，d=2
#此处不做代码的实验，各位可以自己尝试
###############################################
#### 通过acf图pacf图eacf来计算pq ####
###############################################
####step1通过acf和pacf确定序列的类型：
# acf pacf
# MA（q） lag>=q+1的自相关系数全部约小于0 | 逐步将接近0
# AR（p） 逐步将接近0 | lag>=p+1的自相关系数全部约小于0
# ARMA（p,q） 逐步将接近0 | 逐步将接近0
####step2序列的类型确定后：
求pq
# 如果是MA（q）类型的，看acf图，看前几个不为0
# 如果是AR（p）类型的，看pacf图，看前几个不为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）
plot.cf（data）
###############################################
#### 产生ariam模型 ####
###############################################
#pAR的参数
#qMA的参数
#d移动取差的阶数
#arima（data.ts,order=c（p,d,q））#产生模型。
arima只处理时间序列，不处理向量等
p=2;d=0;q=2
data<-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.ahead=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.data）））
lines（x,dnorm（x,mean（input.data）,sd（input.data））,col="red"）#正态分布曲线
}
sol<-shapiro.test（input.data）
if（sol$p.value>alpha）{
print（paste（"success:
服从正态分布，p.value=",sol$p.value,">",alpha））
}else{
print（paste（"error:
不服从正态分布，p.value=",sol$p.value,"<=",alpha））
}
sol
}
#####残差分析
#方案一（自己编写的，不一定好）
#ariam.sol由ariam（）产生的模型
aya.res<-function（ariam.sol）{
#step1：
检验正太性，如果符合正太，表示均值相同：
稳定的
dev.new（）
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，即残差很小
test.p<-0
for（iin1:
20）{
test.p[i]<-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.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.和预测标准误$se
k<-qnorm（（1-alpha/2）,0,1）
pred.data<-data.frame（pred=pred$pred,max<-pred$pred+k*pred$se,min<-pred$pred-k*pred$se）#95%的pred区间值（双侧估计）
pred.data
}
data<-y.arma（a1=0.92,a2=0.65,b1=0.92,b2=-0.65,num=100,pic=F）
sol1<-arima（data.ts,o
升级会员 