太阳黑子数的ARIMA模型

廖海燕

摘 要：从比利时皇家天文台（the Royal Observatory of Belgium）的太阳黑子指数数据中心（the Sunspot Index Data center）的网站获得了1700—2013每年的太阳黑子数的数据。利用R软件结合时间序列建模方法对观测值进行了分析和建模，并利用该模型对未来的太阳黑子数进行了预测。

关键词：太阳黑子数；R软件；fBasics添加包；ARIMA模型

中图分类号：P182.41 文献标识码：A 文章编号：2095-6835（2014）11-0121-03

太阳黑子是太阳光球上出现的一种临时现象，它们在可见光下会比周围区域黑暗。太阳黑子很少单独活动，而是经常成群出现，活跃时会对地球的磁场产生影响，对各类电子产品和电器造成损害，还会使地球南北极和赤道的大气环流作经向流动，进而造成恶劣天气。鉴于太阳黑子的危害性，有必要对其数量进行预测，进而对其进行有效的控制。

从比利时皇家天文台的太阳黑子指数数据中心网站获得从1700—2013期间每年的太阳黑子数数据，时间跨度为314年，即共获得了314个观测值。下面将利用R软件结合时间序列建模的方法对这些数据进行分析，并建立ARIMA预测模型。

1 数据的预处理

一共有314个数据，第一个数据是2013年的太阳黑子数，最后一个数据是1700年的太阳黑子数。为了方便时间序列分析，需要把数据的顺序倒过来。

由于是年度数据，所以每个观测值中的日期属性并不重要，只需提取观测值即可，输入代码：s1=sunspot2014[，2]，得到命名为s1的去日期数据。

将按照时间顺序排列的太阳黑子数的数据对象命名为s3.

2 基本统计分析

输入代码：library（fBasics），加载R的添加包fBasics。为了对常用统计量进行计算，输入代码：basicStats（s3），由输出结果可以得到以下几个结论：①在314个数据中，无缺失值，且记录的最小太阳黑子数数值为0.输入代码：which（s3==0）+1700-1，由输出结果可知，1711年、1712 年和1810年记录的太阳黑子数为0.②记录的最大太阳黑子数数值为190.2，输入代码：which（s3==190.2）+1700-1，由输出可知，记录的太阳黑子数的最大值是1957年。③大约1/4的数据小于16，而1/4的数据大于69，平均每年的太阳黑子数为49.5，标准差是40.3.④输入代码which（s3==40）+1700-1，由输出结果可知，1725年、1741年和1747年的太阳黑子数是中位数40，除了这三年的数据外，将其余数据以40为界分为两部分，大约一半小于40，一半大于40.

3 建模分析

输入代码：

plot（s3，type=h，xlab=""，ylab=""，axes=FALSE）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Time"，cex=2，side=1，line=3）

mtext（"s3"，cex=2，side=2，line=2.5）

（mm=mean（s3））

（m1=rep（mm，nn））

lines（1700：2013，m1，lty=2）

为了观察数据的自相关函数，输入代码：

acf（s3，lag=length（s3）-1，xlab=""，ylab=""，axes=FALSE，main=""）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Lag"，cex=2，side=1，line=3）

mtext（"ACF"，cex=2，side=2，line=2.5）

太阳黑子数具有非常明显的周期性，且周期在11～12之间，为此，对时间序列数据进行周期为11的季节差分，输入代码：s4=diff（s3，11）。为了检测季节差分后数据的平稳性，输入代码：library（fUnitRoots），以加载R的fUnitRoots添加包，再输入代码：adfTest（s4，lags=11，type='ct'），并用扩展的Dickey-Fuller单位根进行检验，结果显示p值为0.030 47，由此可得出季节差分后的数据具有平稳性的结论。

为了观察季节差分后数据的变化趋势，输入代码：

plot（s4，xlab=""，ylab=""，axes=FALSE）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Time"，cex=2，side=1，line=3）

mtext（"s4"，cex=2，side=2，line=2.5）

为了观察季节差分后数据的自相关函数，输入代码：

acf（s4，lag=length（s4）-1，xlab=""，ylab=""，axes=FALSE，main=""）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Lag"，cex=2，side=1，line=3）

mtext（"ACF"，cex=2，side=2，line=2.5）

为了进一步消除数据序列的前后相关性，需要对季节差分后的数据进行常规差分，输入代码：s5=diff（s4），得到被命名为s5的新数据序列。endprint

为了确定ARIMA模型的阶，使用扩展的自相关函数EACF，为此，输入代码：library（TSA），以加载R的TSA添加包，然后输入代码：eacf（s5，11，11），输出的矩阵图如图5所示。

矩阵的行与AR的阶p相对应，矩阵的列与MA的阶q相对应。一般来说，ARIMA模型的阶（p，q）位于“o”组成的三角形的左上角顶点处。

4 模型的建立和诊断

输入代码：library（forecast），以加载R的forecast添加包。为了建立AIC最佳的ARIMA模型，输入代码：（m1=auto.arima（s3）），建立的模型如下（对象名为m1）：

（1-1.802 7*+1.238 4*）（1-）=（1--1.625 7*+0.745*）

模型的估计值为237.6，各系数的标准差依次为：0.084 2，0.130 9，0.077 2，0.059 6，0.059 5，模型的AIC为2 615.51，模型的阶的选定要与图5相对应。

为了对模型的残差进行检验，输入代码：

plot（m1$residuals，type=p，xlab=""，ylab=""，axes=FALSE）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Year"，cex=2，side=1，line=3）

mtext（"m1$residuals"，cex=2，side=2，line=2.5）

（level=rep（0，nn））

lines（1700：2013，level，lty=2）

该模型的残差分布状态比较理想。

输入代码：Box.test（m1$residuals，type=Ljung），进行Box-Ljung混成检验，得到p值为0.882 8，由此可知，模型的残差的各阶相关函数为0.

5 样本外的四步预测

对建立的模型m1进行样本外的四步预测，即预测2014—2017期间的太阳黑子数。输入代码：predict（m1，4），得到2014—2017太阳黑子数的预测值及其标准差，

。

6 结束语

可以对太阳黑子数的ARIMA模型进行进一步的改进，使对太阳黑子历史数据分析和对未来太阳黑子数的预测更加科学和准确，为相关部门制订安全的控制措施提供参考依据。

参考文献

[1]李素兰.数据分析与R软件[M].北京：科学出版社，2013.

[2]Paul Teetor.R语言经典实例[M].北京：机械工业出版社，2013.

[3]Jonathan D.Cryer，Kung-Sik Chan.时间序列分析及应用[M].北京：机械工业出版社，2011.

〔编辑：王霞〕endprint

为了确定ARIMA模型的阶，使用扩展的自相关函数EACF，为此，输入代码：library（TSA），以加载R的TSA添加包，然后输入代码：eacf（s5，11，11），输出的矩阵图如图5所示。

矩阵的行与AR的阶p相对应，矩阵的列与MA的阶q相对应。一般来说，ARIMA模型的阶（p，q）位于“o”组成的三角形的左上角顶点处。

4 模型的建立和诊断

输入代码：library（forecast），以加载R的forecast添加包。为了建立AIC最佳的ARIMA模型，输入代码：（m1=auto.arima（s3）），建立的模型如下（对象名为m1）：

（1-1.802 7*+1.238 4*）（1-）=（1--1.625 7*+0.745*）

模型的估计值为237.6，各系数的标准差依次为：0.084 2，0.130 9，0.077 2，0.059 6，0.059 5，模型的AIC为2 615.51，模型的阶的选定要与图5相对应。

为了对模型的残差进行检验，输入代码：

plot（m1$residuals，type=p，xlab=""，ylab=""，axes=FALSE）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Year"，cex=2，side=1，line=3）

mtext（"m1$residuals"，cex=2，side=2，line=2.5）

（level=rep（0，nn））

lines（1700：2013，level，lty=2）

该模型的残差分布状态比较理想。

输入代码：Box.test（m1$residuals，type=Ljung），进行Box-Ljung混成检验，得到p值为0.882 8，由此可知，模型的残差的各阶相关函数为0.

5 样本外的四步预测

对建立的模型m1进行样本外的四步预测，即预测2014—2017期间的太阳黑子数。输入代码：predict（m1，4），得到2014—2017太阳黑子数的预测值及其标准差，

。

6 结束语

可以对太阳黑子数的ARIMA模型进行进一步的改进，使对太阳黑子历史数据分析和对未来太阳黑子数的预测更加科学和准确，为相关部门制订安全的控制措施提供参考依据。

参考文献

[1]李素兰.数据分析与R软件[M].北京：科学出版社，2013.

[2]Paul Teetor.R语言经典实例[M].北京：机械工业出版社，2013.

[3]Jonathan D.Cryer，Kung-Sik Chan.时间序列分析及应用[M].北京：机械工业出版社，2011.

〔编辑：王霞〕endprint

为了确定ARIMA模型的阶，使用扩展的自相关函数EACF，为此，输入代码：library（TSA），以加载R的TSA添加包，然后输入代码：eacf（s5，11，11），输出的矩阵图如图5所示。

矩阵的行与AR的阶p相对应，矩阵的列与MA的阶q相对应。一般来说，ARIMA模型的阶（p，q）位于“o”组成的三角形的左上角顶点处。

4 模型的建立和诊断

输入代码：library（forecast），以加载R的forecast添加包。为了建立AIC最佳的ARIMA模型，输入代码：（m1=auto.arima（s3）），建立的模型如下（对象名为m1）：

（1-1.802 7*+1.238 4*）（1-）=（1--1.625 7*+0.745*）

模型的估计值为237.6，各系数的标准差依次为：0.084 2，0.130 9，0.077 2，0.059 6，0.059 5，模型的AIC为2 615.51，模型的阶的选定要与图5相对应。

为了对模型的残差进行检验，输入代码：

plot（m1$residuals，type=p，xlab=""，ylab=""，axes=FALSE）

axis（1，cex.axis=1.5）

axis（2，cex.axis=1.5）

mtext（"Year"，cex=2，side=1，line=3）

mtext（"m1$residuals"，cex=2，side=2，line=2.5）

（level=rep（0，nn））

lines（1700：2013，level，lty=2）

该模型的残差分布状态比较理想。

输入代码：Box.test（m1$residuals，type=Ljung），进行Box-Ljung混成检验，得到p值为0.882 8，由此可知，模型的残差的各阶相关函数为0.

5 样本外的四步预测

对建立的模型m1进行样本外的四步预测，即预测2014—2017期间的太阳黑子数。输入代码：predict（m1，4），得到2014—2017太阳黑子数的预测值及其标准差，

。

6 结束语

可以对太阳黑子数的ARIMA模型进行进一步的改进，使对太阳黑子历史数据分析和对未来太阳黑子数的预测更加科学和准确，为相关部门制订安全的控制措施提供参考依据。

参考文献

[1]李素兰.数据分析与R软件[M].北京：科学出版社，2013.

[2]Paul Teetor.R语言经典实例[M].北京：机械工业出版社，2013.

[3]Jonathan D.Cryer，Kung-Sik Chan.时间序列分析及应用[M].北京：机械工业出版社，2011.

〔编辑：王霞〕endprint