《时间序列分析与应用》
课程作业
地震数据（COP.BHZ-24）时间序列分析
一．前言
本次作业选取了第24号文件，共1440个数据。
截取前1200个数据进行理分析，然后建立模型。
之后再对数据进行预测，然后对1200之后的30个数据进行更新，将更新结果与原观测值进行比对分析，最后得出结论。
二．数据处理
1. 数据读取与画图
首先将文件“COP.BHZ.txt”保存到E盘根目录下，以便于读取。
用scan()函数将数据读入，并保存到sugar2文件中。
如图1所示。
图1 数据读取
然后，画出该时间序列图。
横轴表示时间，单位是*10ms，纵轴表示高程，单位是um。
代码及图示如图2、图3所示。
图2 时序图代码
图3 前1200个数据散点图
2. 平稳性检验
从图中看出，该组数据随时间变化基本平稳，仅有小幅波动。
最高点与最低点相差也仅在250um之内。
通过adf.test()函数可以验证该假设，可以看出该序列是平稳的（stationary）。
如图4所示。
然后用求平均函数mean()求出这1200个数据的平均值a，可以从图5看到结果。
图4 平稳性检验结果
图5 求平均值
然后，将原始数据减去平均值，得到一组零均值的新数据，命名为sugar3。
3. 数据建模分析
接下来绘制震前数据的自相关函数和偏自相关函数图像，初步判断其大概符合什么模型。
图6为画出图像的代码，新序列sugar3的ACF、PACF图像如下所示。
图6 ACF、PACF、EACF图像代码
图7 ACF图
图8 PACF图
从ACF、PACF图可以看出，序列一阶之后相关性较强，虽然在第19阶滞后处有超限的情况，但从总体来看，两个图都是拖尾的情况。
因此要借助于EACF 图来做进一步判断。
扩展自相关函数EACF图如下。
图9 EACF图
3 模型识别
由EACF图可以看出此时间序列符合ARMA(0,1)或ARMA(2,2)，根据以上信息尚不能明确判断出具体的模型，要建立确定的模型，就需要排除上述模型中的一种，用模型诊断的方法可以实现。
模型诊断，或模型评价，涉及检验模型的拟合优度，并且如果拟合程度很差，要给出适当的调整建议。
模型诊断的方法有两种：分析拟合模型的残差和分析过度参数化的模型。
下面先使用残差法。
3.1 ARMA(0,1)模型诊断
图10 ARMA(0,1)模型3.1.1 残差法
图11 画残差图的代码
图12 ARMA(0,1)模型残差图3.1.2 分位数图法
图13 分位数法代码
图14 ARMA(0,1)模型的分位数-分位数图3.2 ARMA(2,2)模型诊断
图15 ARMA(2,2)模型
3.2.1 残差法
图16 画残差图代码
图17 ARMA(2,2)模型残差图
3.2.2 分位数图法
图18 分位数法代码
图19 ARMA(0,1)模型的分位数-分位数图
从两种待定模型的残差图和分位数-分位数图看出：两种模型都较好符合，难以比较优劣，因此较难取舍。
下面使用第二种模型诊断的方法，过度拟合法。
3.3 利用过度拟合法进行模型诊断
图20 ARMA(0,1)情况过度拟合
根据过度拟合时应遵守的原则：第一，在拟合时不能同时增加AR和MA部分的阶数；第二应按残差分析建议的方向来扩展模型。
本文中，拟合了MA(1)模型，残差在2阶滞后处仍存在明显相关性，则尝试MA(2)、MA(3)模型，而不是ARMA(1,1)模型。
由此得到以上三种模型的参数估计。
比较前两种模型，可以看出，arima(0,0,2)中的新增的参数为0.0437，不显著地不为零。
而且两种模型的共同的参数相差很小。
再结合第三种模型，同样可以看出新增参数为-0.0096，不显著地不为零，且此三种模型的共同参数相差无几，因此可以判断出ARMA(0,1)模型是较合适的。
另外，结合以上三种模型的对数似然值和AIC值，可以看出，ARMA(0,2)模型AIC值较小，且对数似然值较大，因此选取ARMA(0,2)模型。
下面来诊断ARMA(2,2)模型是否合适。
图21 ARMA(2,2)情况过度拟合
将ARIMA(2,0,2)模型同下面两种模型比较可以看出：虽然额外的参数不显著地不为零，但是它们公同的参数发生了显著的变化，因此可以得出结论：通过过度拟合的诊断，ARMA(2,2)模型不合适。
通过上面的叙述，得出了二阶滑动平均MA(2)模型，即ARIMA(0,0,2)模型，接下来，对该模型的参数进行估计，使用的估计方法是条件平方和估计法、极大似然估计法，通过对不同的方法得出的参数估值进行比较，得到最优估值。
图22 参数估计图
可以看出，两种方法得出的参数估值、方差、对数似然值等相差无几，因此选用哪种方法所得的参数估值都可以，在这里选用对数似然值较大的CSS方法。
可得模型表达式：Y t + 0.0004= e t – 0.0982e t-1 – 0.0436e t-2 – 1.8783
4 模型预测
时间序列建模的主要目标之一，是预测该序列未来的取值。
图23 时间序列预测图
从图中可知，预测结果并不好，好像哪里出了问题。
（可能哪里代码有误，暂时还没有发现。
）下面利用新建模型更新原观测数值30步，即1201至1230。
图24 预测更新代码
图25 预测更新结果
从上图更新的30步观测值可以看出，虽然拟合模型的趋势（黄色）大致与实际观测结果的走向（黑色）大致相同，但是有些点上还是有不小的拟合误
差，例如2号、10号、27号点等。
究其原因，可能是因为所选本组数据相对平缓数据变化相对稳定，因此较难拟合出趋势较明显的模型。
5 心得体会
通过一学期的课程学习，我感觉收货颇丰。
首先，在学习上我不仅学会了如何对一组时间序列进行基本的模型识别、参数估计、模型诊断预测更新等理论方法，还学会了如何用R语言去具体实现。
同时也锻炼了自己的自学及钻研能力。
其次，在课堂上，张老师、李老师除了交给我们课本上的理论知识外，还教导我们为人处世的道理，特别是张老师，从老师的身上学到了如何去尊重别人等基本的却是容易被人忽视的基本礼节。
学术固然重要，礼节却更可贵。
参考书目：
[1] Jonathan D.Cryer, Kung-Sik Chan.时间序列分析与应用[M].机械工业出版社.
[2] 王福林,王吉权,吴昌友,吴秋峰.实数遗传算法的改进研究[J]. 生物数学学报. 2006(01)
[3] 孙晓云,高鑫,王鹏.新型并行遗传算法及其在参数估计中的应用[J]. 计算机工程与应用. 2005(19)
[4] 田小梅,龚静.实数编码遗传算法的评述[J]. 湖南环境生物职业技术学院学报. 2005(01)
[5] 陈晓梅,杨成祥.遗传进化算法在时间序列建模中的应用[J]. 计算机工程与应用. 2005(05)
[6] 来鹏,陈平.基于遗传算法的时序混合模型的参数估计[J]. 山西师范大学学报(自然科学版). 2004(03)
[7] 孙靖,程大章.基于季节性时间序列模型的空调负荷预测[J]. 电工技术学报. 2004(03)
[8] 安德洪,柳湘月,刘嘉焜,许树荆.基于季节ARIMA模型的电力负荷建模与预报[J]. 天津大学学报. 2004(02)
[9] 孙晓云,蔡远利.利用改进遗传算法的参数估计[J]. 自动化技术与应用. 2004(01)
[10] 何大阔,王福利.一种提高遗传算法全局收敛性的方法[J]. 东北大学学报. 2003(06)。