TimeSeries with R

发表于 | 阅读次数： 
1
2
3
4
5
6
7
8
9
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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#一、数据读取
ts(data,start,end,frequency)
##文件是矩阵时按列读入，可以通过先转置再线性化为向量，代码如下
D=read.table('zz.txt')
D=t(D)#转置
D=as.vector(D)#将矩阵线性化为向量
#数据读取结束

#二、描述时间序列
ts.plot(D)#ts(D);plot(D)

#三、自相关和偏相关
acf(D,lag,type,plot=True)
pacf()

#四、平稳性检验
#（1）时间序列图检验
#（2）自相关图检验
#1.确定性时间序列
x=seq(0.1,5,0.02)
y=0.8*x+2
D=jitter(y,amount=0.2)#加入随机扰动
D=ts(D)
par(mfrow=c(1,3))
plot(D)
#以上是时间序列图检验，下面进行自相关图检验                                                                           
acf(D)
pacf(D)

#2.周期性时间序列
x=seq(0.1,15,0.1)
y=sin(x)
par(mfrow=c(3,1))
D=jitter(y,amount=0.2)
ts.plot(D)
acf(D)
pacf(D)

#3.伪周期时间序列
x=seq(0.1,20,0.1)
y=100*sin(x)/x
D=jitter(y,amount=3)
par(mfrow=c(3,1))
plot(D)
acf(D)
pacf(D)

#五、纯随机性检验
#H0:是白噪声序列(acf,pacf基本为0)
D=rnorm(300,0)
D=ts(D)
par(mfrow=c(3,1))
plot(D)
acf(D)
pacf(D)
#接下来利用Box.test()进行检验
#result=0;lag=0;lb=0;p=0;
for(i in 1:12)
{	btest=Box.test(D,type="Ljung-Box",lag=i)
	lag[i]=i
	lb[i]=btest$statistic
	p[i]=btest$p.value
	result=cbind(lag,lb,p)
}
print(result)

#六、ARMA模型拟合与识别 
#arima.sim(model,n),model是模型 ，n是序列的长度
 
#1.AR(1):X(t)=0.8*X(t-1)+E(t)
D=arima.sim(list(order=c(1,0,0),ar=0.8),n=200)#order的三个参数分别对应AR,I,MA的阶数
plot(D)
par(mfrow=c(1,2))
acf(D)	
pacf(D)

#2.AR(2):X(t)=0.8*X(t-1)+0.1X(t-2)+E(t)
D=arima.sim(list(order=c(2,0,0),ar=c(0.8，0.1)),n=200)#order的三个参数分别对应AR,I,MA的阶数
plot(D)
par(mfrow=c(1,2))
acf(D)	
pacf(D)

#3.MA(2):x(t)=E(t)-0.7*E(t-1)-0.2*E(t-2)
D=arima.sim(list(order=c(0,0,2),ma=c(0.7,0.2)),n=200)
plot(D)
par(mfrow=c(1,2))
acf(D)
pacf(D)

#ARMA(2,2):X(t)=0.3*X(t-1)+0.4*X(t-2)+E(t)-0.5*E(t-1)-0.4*E(t-2)
D=arima.sim(list(order=c(2,0,2),ar=c(0.3,0.4)ma=c(0.5,0.4)),n=200)
plot(D)
par(mfrow=c(1,2))
acf(D)
pacf(D)

#七、ARMA模型的估计
D=arima.sim(list(order=c(2,0,2),ar=c(-0.5,0.4),ma=c(0.5,0.4)),n=1000)
ARIMA=arima(D,order=c(2,0,2),method="ML")
summary(ARIMA)
ARIMA$coef

#八、ARMA模型的预测

D=arima.sim(list(order=c(2,0,2),ar=c(0.5,-0.2),ma=c(0.5,0.2)),n=100)
plot(D)
par(mfrow=c(1,2))
acf(D)
pacf(D)
ARIMA=arima(D,order=c(2,0,2),method="ML"，include.mean=FALSE)
summary(ARIMA)
resid=ARIMA$residual
Box.test(resid,lag=6)
Box.test(resid,lag=12)
pre=predict(ARIMA,n.ahead=10)
U=pre$pred+1.96*pre$se#置信上限
L=pre$pred-1.96*pre$se#置信下限
ts.plot(pre$pred,col=1:2)
lines(U,col="blue",lty="dashed")
lines(L.col="blue",lty="dashed")

凡希 wechat
喜欢所以热爱，坚持干货分享，欢迎订阅我的微信公众号
呐，请我吃辣条