基于TAR模型的太阳黑子非线性时间序列预测

1、基于TAR模型的太阳黑子非线性时间序列预测摘要：太阳黑子数目的变化对地球的气候、农业、通信、导航等方面影响巨大因此对太阳黑子数目进行预测具有十分重要的意义。本文对1945-2005年的太阳黑子数据建立基于不同时间段的门限自回归模型(TAR)，分析太阳黑子时间序列的变动特征并对未来10年的太阳黑子数进行预测。从模型诊断结果可以得出：TAR(2；3，5)模型能很好地拟合该太阳黑子的非线性时间序列，相应的预测值也比较精确。关键词：太阳黑子 非线性时间序列 TAR模型 预测0 引言太阳黑子的太阳活动中最基本的现象，它是在太阳的光球层桑发生的一种太阳活动，太阳黑子是表示太阳活动强弱的一项重要指标，它是典

2、型的复杂时间序列，地磁变化、大气运动、气候异常、海洋活动、等都和太阳黑子数的变化有着不同程度的关系。对太阳黑子活动进行有效的预测以此来分析地球环境的变化有着十分重要的价值。因此，历来世界各国都十分重视对太阳黑子活动的预测工作，以便能够采取防范措施，避免意外的灾难性事故发生。任晶等（2014）建立了基于相空间重构的神经网络和神经网络的太阳黑子时间序列预测模，并在MATLAB环境下进行预测仿真，仿真结果表明，建立的模型预测精度较好。向昌盛等（2011）提出了一种相空间重构和最小二乘支持向量机（LSSVM）参数的联合优化方法，实验结果表明联合优化方法预测精度比较好，而且优化速度更快。对于太阳黑子的预

3、测文献中，运用向量自回归（TAR）模型进行预测的还比较少。本文对1945-2005年的太阳黑子数据建立基于不同时间段的门限自回归模型(TAR)，分析太阳黑子时间序列的变动特征并对未来10年的太阳黑子数进行预测。1 TAR模型门限自回归模型作为一类非线性模型，能够解释金融数据中的非线性性质。它首先是由Tong(1980)提出的。门限自回归模型设定某一特定的时点，时间序列的运动方式从一种机制跳跃到了另一种机制，同时这种跳跃是离散的。门限自回归模型在拟合实际数据时具有较好的性质，但是由于建立门限自回归模型的步骤比较复杂，直到Ruey S.Tsay (1989)提出了相对来说比较简易的建模及检验方法后

4、，这类模型才被人们广泛地应用。一般地，对于时间序列称为满足一个k阶门限自回归模型(TAR)，其门限变量为，假设初始值是已知的，如果其满足下式：其中，；表示门限；k表示段数，是正整数，表示第段，；要求门限变量在空间上可测（即：是内元素的可测函数），在这里是时刻t-1之前可用信息的域；参数d被称为延迟变量，也是正整数；是均值为0，方差为，独立的序列，实际上对于相同的，服从独立同分布。们组上述条件的模型通常记为：。在实际应用中，由Tong(1983，1990)提出了各种状态下涉及若干含有分离高阶AR(p)过程的不同状态的TAR过程。两状态的TAR模型的一般形式：两个子模型的自回归阶数不必相同，延迟参

5、数d可以大于最大的自回归阶数。这里存在由的值定义的两个可分离状态。如果门限已知，则依据的值是否在门限之上或之下，分离观测值。然后用OLS估计方法估计每段方程式，其滞后期长度d根据AR模型的方法确定，因此可以用t检验对单个系数进行检验。在门限值未知的情况下，根据Chan（1993）提出的方法来获得门限r的超一致估计，为了确保在门限两边有适当数量的观测值，从检索中排除了按大小排列的最高和最底部分各10%的值。然后，运用程序是r取每一个观测值估计TAR模型，其中相应残差平方和最小的回归方程含有门限的一致估计。2 实证分析2.1 数据来源及预处理本文数据来源于R软件TSA包中数据文件名为spots的数

6、据，是1945年至2005年美国一年一度的网络观测的太阳活动的加权平均值。为降低数据的异方差性，对原始数据进行开放变换，使用开方变换后的数据进行建模。本文所有计算都借助R语言完成。2.2 非线性检验非线性时间序列分析的研究热情大约源于20世纪70年代后期对实际数据表现出的非线性动态进行建模的要求。非线性时间序列模型通常展示出非常丰富的动态结构。实际上，很简单的一个确定性的非线性差分方程可以有以下意义上的混沌解：时间序列的解对初始值敏感，且看起来基于相关分析是难以将其与白噪声序列加以区分的。因此，非线性时间序列分析也许可以提供更精确的预测，在状态空间的某些部分可能是非常重要的，从而对所研究数据的

7、动态给出新颖的启示。门限模型是处理非线性时间序列的方法，所以建立门限模型之前要检验时间序列的是否具有非线性。2.2.1 图解法非线性检验在ARIMA建模中，新息（误差）过程通常设定为服从独立同分布的正态分布。这一正态误差假设意味着平稳时间序列同时也是正态过程，即时间序列观测的任意有限集合服从联合正态分布。在作统计推断时，采用正态假设主要是为了方便。实际当中，可能要考虑具有非正态新息的ARIMA模型。若维持正态误差假设，则非线性时间序列通常为非正态分布。因此，对于非线性的考察，可通过检验时间序列观测的任意有限集合是否服从联合正态分布来进行。二元正态分布的散点图应该与距离中心越远密度越低的椭圆形数

8、据云相似。若异于这个模式（例如，数据云中存在大洞），则有可能表明数据是非正态的，而相关的过程是非线性的。图1太阳黑子序列滞后回归图图1是太阳黑子序列的滞后回归图。从图中可以看出，1阶至4阶滞后的散点图的中心明显存在大洞，表明数据必为非正态的。另外，对2阶至3阶滞后的回归函数的估计显出了很强的非线性，表明非线性的数据机制。下面用更精确的方式检验时间序列的非线性。2.2.2 基于Keenan和Tsay的非线性检验在时间序列分析中，提出了若干用于建立非线性模型必要性评估的检验方法。其中某些检验，诸如已在Keenan（1985）、Tsay（1986）及Luukkonen等（1988）中有所研究的，可归

9、为用于特定非线性可选方案的拉格朗日乘子检验。这里主要用Keenan检验和Tsay检验。表1 非线性检验检验方法统计量P值检验结果Keenan检验16.280.0002拒绝Tsay检验3.60.0009拒绝表1是基于Keenan和Tsay的非线性检验统计表，由表1可以看出，Keenan和Tsay检验的p值分别是0.0002和0.0009，拒绝线性假设，说明太阳黑子时间序列是非线性的，与上面图解法非线性检验一致。2.3 TAR模型的估计对于线性ARIMA模型，AR的阶数可以通过最小化AIC来估计，对于固定的和，TAR模型本质上分别拟合阶数为和的两个AR模型，因此AIC成为其中除之外，参数的个数等于

10、。现在，用最小化AIC(MAIC)估计法估计参数时所最小化的AIC满足如下约束条件：在某些可保证任何区域都有足够数据用于估计的区间内搜索门限参数。基于把门限参数作为额外参数这一直观的想法，如上所得最小的AIC加2后的结果称为估计门限模型的名义AIC。因为门限参数一般会给模型增加很多的灵活性，所以也可能给模型增加不止一个自由度。从渐进角度看，这可能等价于给模型增加了三个自由度。在对TAR模型进行估计之前，要首先确定最大阶数p，这里用偏自相关图结合AIC确定的AR模型的阶数。开方后太阳黑子时间序列如图2所示。从图2可以看出，AR模型的阶数为4或者5，再结合AR模型的AIC值确定最佳阶数。图2 AR

11、模型偏自相关图表2 最佳滞后阶数选择滞后阶数AIC值3198.264193.865187.816190.03由表2可以看出，AR(5)模型的AIC值最小，所以最佳滞后阶数选为5比较合适。使用MAIC方法，大致在第10百分位数和第90百分位数之间搜索门限，表3给出了时估计TAR模型的名义AIC值。当时，名义AIC值最小，所以对延迟的估计为2。表3 开方后（太阳黑子）的序列之拟合TAR模型的名义AIC AIC1149.95.88552110.56.06353124.66.60254126.28.04355150.58.1545MAIC法在下区域所选择的阶数是3，在上区域所选阶数是5，使用落在每个区

12、域中的数据，应用普通最小二乘（OLS)估计各区域的子模型。因此，对于噪声方差更为无偏的估计可由有效样本容量归一化的区域内误差平方的残差和估计。这里的有效样本容量等于对应区域中对的数据数量减去相关子模型的自回归参数的个数（包括截距项）。第个区域中的“无偏”噪声方差及其相应的极大似然估计量间的关系。用下式表示：其中是第个子模型的自回归阶数。另外，近似服从个自由度的分布。表4是对模型参数估计的结果。表4 对太阳黑子数据拟合的TAR(2；3，5)模型：MAIC法估计值标准误差t统计量P值26.058下区域9.521.287.440.001.050.176.170.00-1.200.29-4.130.0

13、0-0.560.24-2.280.040.8211上区域5.690.866.600.000.370.113.470.000.380.123.250.00-0.070.11-0.590.56-0.360.11-3.340.00-0.110.07-1.500.140.2141由参数估计表5可以得出拟合的TAR(2；3，5)模型方程：此估计结果说明被估计序列存在门限特征。2.4 模型诊断2.4.1 骨架分析研究骨架的长期（渐进）行为是评价非线性模型的一种方法，省略模型的噪声项得到其骨架，即用0代替噪声项。骨架可能发散到无穷，也可能收敛到极限点、极限循环或者奇异吸引子。平稳ARIMA模型的骨架总是收敛

14、到某个极限点，平稳非线性模型的骨架有可能展示出非常复杂的动态。对太阳黑子序列拟合TAR(2；3，5)模型，期骨架如图3所示。由图3可以看出，拟合的模型是平稳的，并且骨架收敛到极限点。骨架明显的长期稳定性意味着的拟合TAR(2；3，5)模型平稳。图3 太阳黑子序列TAR(2；3，5)模型的骨架2.4.2 残差分析上面对拟合的模型进行了骨架分析，这里讨论用残差分析进行模型诊断的统计方法。原始的残差定义为自真实数据减去拟合值，其中第个拟合值是利用给定的过去的值对条件均值所做的估计，即残差由下式给出：这些残差与得自拟合子模型的原始残差相同。用适当的标准差来规范化原始残差，即得到标准化的残差：即用下（上

15、）区域的噪声标准差估计规范下（上）区域的原始残差。正如线性情况里的一样，标准残差的时间序列图看上去应为随机的，原因在于如果TAR模型是真实的数据机制，即若TAR模型被正确地识别，那么标准残差应该是渐进独立同分布的。图4 太阳黑子序列TAR模型的残差分析图4是太阳黑子序列拟合TAR(2；3，5)模型的混合检验。最上面的图是标准残差的时间序列图，可以看出残差都在±2倍标准差内。中间的图是标准残差的ACF图，从ACF图可以看出标准残差无自相关。最下面的图是混合检验，由混合检验图可以看出所有的P值都大于0.05，检验结果是不显著的，说明拟合TAR(2；3，5)模型是稳健的。图5给出了标准残差

16、的正态得分QQ图，该图中显示为直线，所以误差服从正态分布。结合以上分析，建立TAR(2；3，5)模型很好地拟合了太阳黑子序列。图5标准残差的正态QQ图2.5 预测时间序列未来取值的不确定性（比如）完全取决于给定当前和过去数据，的条件概率分布，称为步向前预测分布。对于具有正态误差的ARIMA模型，其所有预测分布都是正态的，这就大大简化了对预测置信区间的计算，因为只要找到预测分布的均值和方差即可。但通常非线性模型的预测分布不是正态的，也不可解。因而预测区间可能要通过模拟来计算。这一模拟方法可以通过一阶非线性自回归模型来说明：给定，有。因此，可以从误差分布抽样得到，并计算，从而由所得一步向前预测分布

17、可以获取的一个实现。独立地重复B次，该过程，比如1000次，可自一步向前预测分布得到B个随机样本。可以用这B个样本的均值作为一步向前预测均值的估计。通过递推非线性自回归模型，对任意整数，易推广模拟法来求出步向前预测分布：其中是从误差分布抽样得到的个值的随机样本。该过程可以重复B次，得到步向前预测分布的随机样本，使用该样本可以计算的预测区间或者其他任意预测的描述统计量。图6 太阳黑子序列预测对开方后的太阳黑子序列基于用的拟合TAR(2；3，5)模型来计算其预测区间，如图6所示，其中中间的虚线是预测分布的中位数，下面和上面的虚线是预测分布的第2.5百分位数和第97.5百分位数。这里模拟的规模为10

18、000次，预测分布的中位数可以作为点预测量。预测区间将趋于拟合TAR模型平稳分布的第2.5百分位数和第97.5百分位数。3 结论TAR模型是针对非线性时间序列数据的一种建模方法，本文对1945-2005年的太阳黑子数据建立基于不同时间段的门限自回归模型(TAR)，分析太阳黑子时间序列的变动特征并对未来10年的太阳黑子数进行预测。模型的骨架分析表明运用TAR(2；3，5)模型拟合太阳黑子时间序列数据比较稳定，而从模型拟合的残差分析中可以看出，残差不存在自相关、没有异常值并且服从正态分布，更进一步说明了所选取的模型很好地拟合了太阳和黑子序列数据。总之，太阳黑子时间序列数据表现出来的非线性使得门限自

19、回归模型(VAR)可以对其进行分析、预测。而从模型诊断结果来看，VAR模型很好地拟合太阳黑子序列，并可以进行预测，这表明VAR模型在太阳黑子预测的应用方面是合适的。4 参考文献1 靳晓婷等. 汇改后人民币汇率波动的非线性特征研究基于门限自回归TARJ.财经研究，2008(09)2 俞毅. GDP 增长与能源消耗的非线性门限对中国传统产业省际转移的实证分析J.产业经济，2010(12)3 段忠东.房地产价格与通货膨胀、产出的非线性关系基于门限模型的实证研究J.金融研究，2012(08)4 张五六.两部门生产函数门限模型及应用以能源消费与经济增长关系为例J.数理统计与管理，2010(06)5 任晶

20、等.基于相空间重构的神经网络太阳黑子预测模型J.计算机仿真，2014(01)6 向昌盛等.混沌时间序列预测模型参数的联合优化J.信息与控制，2011(05)7 (美)克莱尔(Cryer，J.D)等著，潘红宇等译.时间序列分析及应用：R语言M，北京：机械工业出版社，2011.1 5 附录R语言代码：library(TSA)data(spots)lagplot(sqrt(spots)#太阳黑子序列滞后回归图Keenan.test(sqrt(spots)#Keenan非线性检验Tsay.test(sqrt(spots)#Tsay非线性检验#AR模型滞后阶数选择pacf(sqrt(spots)m6=a

21、rma(sqrt(spots),order=c(6,0,0);summary(m6)m5=arma(sqrt(spots),order=c(5,0,0);summary(m5)m4=arma(sqrt(spots),order=c(4,0,0);summary(m4)m3=arma(sqrt(spots),order=c(3,0,0);summary(m3)#选择TAR模型最大的自回归阶数AICM=NULLfor(d in 1:5)spots.tar=tar(y=sqrt(spots),p1=5,p2=5,d=d,a=0.1,b=0.9)AICM=rbind(AICM,c(d,spots.tar$AIC,signif(spots.tar$thd,5),spots.tar$p1,spots.tar$p2)colnames(AICM)=c('d','nominal AIC','r