时间序列分析作业讲解.docx
- 文档编号:6338010
- 上传时间:2023-05-09
- 格式:DOCX
- 页数:13
- 大小:419.11KB
时间序列分析作业讲解.docx
《时间序列分析作业讲解.docx》由会员分享,可在线阅读,更多相关《时间序列分析作业讲解.docx(13页珍藏版)》请在冰点文库上搜索。
时间序列分析作业讲解
《时间序列分析与应用》
课程作业
地震数据(COP.BHZ-24)时间序列分析
一.前言
本次作业选取了第24号文件,共1440个数据。
截取前1200个数据进行理分析,然后建立模型。
之后再对数据进行预测,然后对1200之后的30个数据进行更新,将更新结果与原观测值进行比对分析,最后得出结论。
二.数据处理
1.数据读取与画图
首先将文件“COP.BHZ.txt"保存到E盘根目录下,以便于读取。
用scan()ffi数将数据读入,并保存到sugar?
文件中。
如图1所示。
>sugar2=scan(1COP.BHZ-024.txt")
Read1440items
>
图1数据读取
然后,画出该时间序列图。
横轴表示时间,单位是*10ms,纵轴表示高程,单位是um。
代码及图示如图2、图3所示。
>wln.grapli(widch=-l•875rheighL=2.5,poinEsize=8)
>plot(sugax2(0:
1200]fxlab=,*10nis'.yLab=,uni,r9q9)
>
图2时序图代码
a20040060080010001200
餐10ms
图3前1200个数据散点图
2.平稳性检验
从图中看出,该组数据随时间变化基本平稳,仅有小幅波动。
最高点与最低点相差也仅在250um之内。
通过adf.test()函数可以验证该假设,可以看出该序列是平稳的(stationary)o如图4所示。
然后用求平均函数mean()求出这1200个数据的平均值a,可以从图5看到结果。
>library(tsexies)
>adf・rest(sugax2[0:
1200])
Augirtent皂dDickey-FullerTest
dara:
suqax2[0:
1200]
Dickey-Fuller=一9・3423,Lagorder-10^p-value=0.01alternativehypothesis:
stationary
图4平稳性检验结果
>a=Tnean(sugar2[0:
1200]、
>a
[1]7・878333
■
图5求平均值
然后,将原始数据减去平均值,得到一组零均值的新数据,命名为sugar3o
>5UQar3:
=sugar2[0:
1200]-a
3.数据建模分析
接下来绘制震前数据的自相关函数和偏自相关函数图像,初步判断其大概符合什么模型。
图6为画岀图像的代码,新序列sugar3的ACF、PACF图像如下所示。
>win.graph(width=4・87SFheight=3,pointsize=9)
>acf(sugar3)
>win•graph(width=4・总*75『height=3『pointsize=&)
>pacf(sugax3)
>library(ISA)
>eacf(sugar3)
图6ACF、PACF、EACF图像代码
Lag
图7ACF图
Seriessugar3
图8PACF图
从ACF、PACF图可以看出,序列一阶之后相关性较强,虽然在第19阶滞后处有超限的情况,但从总体来看,两个图都是拖尾的情况。
因此要借助于EACF图来做进一步判断。
扩展自相关函数EACF图如下。
0
1
2
34
5
6
7
89
10
11
12
13
0
X
o
o
oo
o
o
o
oo
o
o
o
o
1
X
X
o
oo
o
o
o
oo
o
o
o
o
2
X
X
o
oo
o
o
o
oo
o
o
o
o
3
X
X.
X.
oo
o
Q
Q
oo
Q
Q
o
o
4
X
X.
X
□□
□
o
o
oo
o
o
o
o
5
X
X
X
XX
□
o
o
oo
o
o
o
o
6
X
X
X
□X
X
o
o
oo
o
o
o
o
7
X
X
0
X□
X
X
o
oo
o
o
o
o
图9EACF图
3模型识别
IIIEACF图可以看出此时间序列符合ARMA(O,1咸ARMA(2,2),根据以上信息尚不能明确判断出具体的模型,要建立确定的模型,就需要排除上述模型中的一种,用模型诊断的方法可以实现。
模型诊断,或模型评价,涉及检验模型的拟合优度,并且如果拟合程度很差,要给出适当的调整建议。
模型诊断的方法有两种:
分析拟合模型的残差和分析过度参数化的模型。
下面先使用残差法。
3.1ARMA(0,1)模型诊断
>sugar3=arima(sugar3zorder=c(0z0r1))
>ml・sugar3
Call:
arima(x=sugar3,order=c(QfQf1)}
Coefficients:
malintercept
0.09310.0004
s.e.0.02751.1529
sigmai八2estimatedas1335:
loglikelihood=-6020・78zale=12045・57
图10ARMA(0,l)模型
3.1.1残差法
>win・Qiraph(width=4.875,heighc=3zpointsize=8)
>plot(rstandard(ml.sugar3)rylab=1StandardizedResiduals1#o1)
>abline(h=0)
%
图11画残差图的代码
今—
Time
sCTnp-s①比PQNPJEPUEJS
图12ARMA(OJ)模型残差图
3.1.2分位数图法
>win・graph(width=4・875fheight=2.5fpoxn七sIn已=9)
>q(inorni(residuals(ml.sugar3))
>qqline(residuals(ml・sugar3)rcol=Bred,')
图13分位数法代码
NormalQ・QPlot
TheoreticalQuantiles
图14ARMA(O,1)模型的分位数-分位数图
3.2ARMA⑵2)模型诊断
>m2・sugar3=arima(sugar3Forder=c(2Z2))
>m2・sugar3
Call:
arima(x=sugar3,order=c(2Z0,2))
Coefficients:
arl
ar2
mal
ma2
intercept
-0.6134
0.0988
0.712
0.0052
・0.0131
0・3146
0・2736
0.316
0.2735
1・1940
sigma八2estimatedas1331:
loglikelihood=-6019・06,aic=12048.13
图15ARMA(乙2)模型
3.2.1残差法
>win.graph(width=4.875rheight=3fpointsize=8)
>plot(rsrandard(m2.sugar3)rylab='ScandaMdizEdResiduals1rtype=*o1)
>abline:
(h=0)
图16画残差图代码
0200400600SOO10001200
Time
图17ARMA(2,2)模型残差图
3.2.2分位数图法
>win・graph(width=4・875fheight=2・5rpointsizE=9)
>qqnorm(residuals(m2.sugar3))
>qqline(residuals(m2・3ugar3).col=1red.1)
>
图18分位数法代码
NormalQ-QPlot
-3-2-10123
TheoreticalQuan-tiles
图19ARMA(OJ)模型的分位数■分位数图
从两种待定模型的残差图和分位数-分位数图看出:
两种模型都较好符合,
难以比较优劣,因此较难取舍。
下面使用第二种模型诊断的方法,过度拟合法。
3.3利用过度拟合法进行模型诊断
3.3.1ARMA(OZ1)的过度拟合
>arima(sugar3rorder=c(0,Of1))
Call:
arima(x=sugar3zorder=c(Q,Q,1))
Coefficients:
malintercept;
0.09310・0004
0・02751.1529
estimatedas1335:
loglikelihood=-6020.78,aic=12047.57
>
>a二ima(suga二3,ozrdEhc(0,0,2))
Call:
arima(x=sugar3,order=c(0z0,2))
Coefficients:
Call:
arima(x=sugar3zorder=c(0,0,3))
Coefficients:
malma2ma3intercept
0.09790.042-0.0026-0.0003
s・^・0・02890.0290・02871・1909
图20ARMA(OJ)情况过度拟合
根据过度拟合时应遵守的原则:
第一,在拟合时不能同时增加AR和MA部分的阶数;第二应按残差分析建议的方向来扩展模型。
本文中,拟合了IVIA(l)模型,残差在2阶滞后处仍存在明显相关性,则尝试MA
(2).MA(3)模型,而不是ARMA(1Z1)模型。
由此得到以上三种模型的参数估计。
比较前两种模型,可以看出,arima(0,0,2)中的新增的参数为0.0437,不显著地不为零。
而且两种模型的共同的参数相差很小。
再结合第三种模型,同样可以看岀新增参数为-0.0096,不显著地不为零,且此三种模型的共同参数相差无儿,因此可以判断出ARMA(0,l)模型是较合适的。
另外,结合以上三种模型的对数似然值和AIC值,可以看出,ARMA(0,2)模型AIC值较小,且对数似然值较大,因此选取ARMA(0,2)模型。
下面来诊断ARMA⑵2)模型是否合适。
3.3.2ARMA(2,2)的过度拟合
>arima(sugar3zorder=c(2r0,2))
Call:
arima(x=sugar3,order=c(2Z0f2))
Coefficients:
arl
ar2
mal
ma2
intercept
-0・6134
0・0988
0.712
0・0052
-0・0131
S・E・
0・3146
0・2736
0・316
0・2735
1・1940
sigma八2estimatedas1331:
loglikelihood=-6019・O6,aic=12050・13
>
>arima(sugar3rorder=c(3r0,2))
Call:
arima(x=sugar3zorder=c(3z0,2))
Coefficients:
arl
ar2
ar3
mal
ma2
intercept
0・1539
0・9163
-Q・0933
-0・0584
-0・8974
0・1097
0・1621
0・1534
0・0299
0・1609
0・1555
1・9874
sigmaA2estimatedas1328:
loglikelihood=-6017.49zaic=12048.97
>arima(sugar3rorder=c(2f0r3))
Call:
arima(x=sugar3forder=c(2Z0,3))
Coefficients:
arl
ar2
mal
ma2
ma3
intercept
-0.4580
0.2187
0.5567
-0.1310
-0.0182
0・00丄0
1・1694
0・7703
工・1690
0・8844
0・1283
1・1961
图21ARMA(2,2)t#况过度拟合
将ARIMA(2,0,2)模型同下面两种模型比较可以看出:
虽然额外的参数不显著地不为零,但是它们公同的参数发生了显著的变化,因此可以得出结论:
通过过
度拟合的诊断,ARMA⑵2)模型不合适。
通过上面的叙述,得出了二阶滑动平均MA
(2)模型,即ARIMA(0,0,2)模型,接下来,对该模型的参数进行佔计,使用的估计方法是条件平方和估计法、极大似然估计法,通过对不同的方法得岀的参数估值进行比较,得到最优估值。
>arima(sugar3rorder=c(0,0,2)zmerhod='CSS1)
Call:
arima(x=sugar3rorder=c(0z0,2}rmethod=f,CSSrr)
Coefficients:
malma2
Q・09820・0436
3・己・Q・02890・0285
intercept-0・00041・2030
sigma^2estimatedas
1332:
partloglikelihood=-6019・61
>arima(sugar3rorder=c(0.0,2)zmerhod='ML1)
Call:
arima(x=sugar3rorder=c(0z0,2}rmethod=r,MLr,)
Coefficients:
malma2
Q・09810・0437
intercept-0・00042030
s・^・Q・02890・0286
图22参数估计图
可以看岀,两种方法得岀的参数估值、方差、对数似然值等相差无儿,因此
选用哪种方法所得的参数估值都可以,在这里选用对数似然值较大的CSS方法。
可得模型表达式:
Yt+0.0004=et-0.0982et-i-0.0436et-2-1.8783
4模型预测
时间序列建模的主要LI标之一,是预测该序列未来的取值。
020040060080010001200
Time
图23时间序列预测图
从图中可知,预测结果并不好,好像哪里出了问题。
(可能哪里代码有误,
暂时还没有发现。
)下面利用新建模型更新原观测数值30步,即1201至1230o
>r-numeric(30)
>data1=sugar2[0;1200]
>fox(iin1:
30)
+t=predict(arima.(data!
order=c(0r0F2))Fn.ahead=30)
>r[i]=t$pred|
图24预测更新代码
1:
30
O0「s£0」s)ns
图25预测更新结果
从上图更新的30步观测值可以看出,虽然拟合模型的趋势(黄色)大致与实际观测结果的走向(黑色)大致相同,但是有些点上还是有不小的拟合误差,例如2号、10号、27号点等。
究其原因,可能是因为所选本组数据相对平缓数据变化相对稳定,因此较难拟合出趋势较明显的模型。
5心得体会
通过一学期的课程学习,我感觉收货颇丰。
首先,在学习上我不仅学会了如何对一组时间序列进行基本的模型识别、参数估讣、模型诊断预测更新等理论方法,还学会了如何用R语言去具体实现。
同时也锻炼了自己的自学及钻研能力。
其次,在课堂上,张老师、李老师除了交给我们课本上的理论知识外,还教导我们为人处世的道理,特别是张老师,从老师的身上学到了如何去尊重别人等基本的却是容易被人忽视的基本礼节。
学术固然重要,礼节却更可贵。
参考书目:
[1]JonathanD.Cryer,Kung-SikChan.时间序列分析与应用[M].机械工业出版社.
[2]王福林,王吉权,吴昌友,吴秋峰.实数遗传算法的改进研究[J].生物数学学报.2006(01)
⑶孙晓云,高鑫,王鹏.新型并行遗传算法及其在参数估计中的应用[」]•计算机工程与应用.2005(19)
[4]田小梅,龚静.实数编码遗传算法的评述[J].湖南环境生物职业技术学院学报.2005(01)
⑸陈晓梅,杨成祥.遗传进化算法在时间序列建模中的应用[J].计算机工程与应用.2005(05)
⑹来鹏,陈平.基于遗传算法的时序混合模型的参数佔计[J].山西师范大学学报(自然科学版).2004(03)
[7]孙靖,程大章.基于季节性时间序列模型的空调负荷预测[J].电工技术学报.2004(03)
[8]安德洪,柳湘月,刘嘉焜,许树荆.基于季节ARIMA模型的电力负荷建模与预报[J].天津大学学报.2004(02)
⑼孙晓云,蔡远利.利用改进遗传算法的参数佔计卩]・自动化技术与应用.
2004(01)
[10]何大阔,王福利.一种提高遗传算法全局收敛性的方法卩].东北大学学报.
2003(06)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 时间 序列 分析 作业 讲解