时间序列分析作业讲解

《时间序列分析与应用》

课程作业

地震数据（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号点等。究其原因，可能是因为所选本组数据相对平缓数据变化相对稳定，因此较难拟合出趋势较明显的模型。