太阳黑子的时间序列模型
太阳黑子是在太阳的光球层上发生的一种太阳活动,是太阳活动中最基本、
最明显的。太阳黑子很少单独活动,通常是成群出现。从长期的黑子相对数记录
可见,黑子相对数的平均值明显的表现出 11 年左右的周期性,最短为 9.0 年,
最长为 13.6 年。黑子相对数的年平均值的极大和极小年份,分别称为太阳活动
的极大年(峰年)和极小年(谷年)。太阳黑子活跃时会对地球的磁场产生影响,
当太阳上有大群黑子出现的时候,会出现磁暴现象使指南针乱抖动,不能正确地
指示方向;平常很善于识别方向的信鸽会迷路;无线电通讯也会受到严重阻碍,
甚至会突然中断一段时间,这些反常现象将会对飞机、轮船和人造卫星的安全航
行、还有电视传真等等方面造成严重威胁。因此,对太阳黑子数据进行分析建立
模型,进而对其数据进行预测有利于对太阳黑子活跃事件带来危害的预防与控制。
本文从比利时皇家天文台的太阳黑子指数数据中心网站获得 1700 年至 2016
年的太阳黑子年度数据,时间跨度为 317 年,共获得 317 个数据,对数进行分析,
建立季节时间序列模型,并对模型进行检验,最后对太阳黑子数据进行预测,得
到 2017 年至 2020 年的太阳黑子数量的预测值。
一、利用 R 软件对太阳黑子数据进行处理
利用 R 软件对太阳黑子数据进行分析,得到太阳黑子数据的时序图,如图
1-1:
图 1-1 的 R 程序如下:
setwd("D:/Data")
w<-read.table(file="sunspot.txt")
names(w)=c("Year","Number")
w
s1=w[,2]
s1
s3=ts(s1,star=1700)
plot(s3,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)
从上图中可以看出太阳黑子数量有明显的周期性,对太阳黑子数据进行分析
得到其条形图和自相关和偏自相关图,如图 1-2、图 1-3、图 1-4:
图 1-1 太阳黑子数量时序图
图 1-2 的 R 程序如下:
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)
nn=length(s3)
(mm=mean(s3))
(m1=rep(mm,nn))
lines(1700:2016,m1,lty=2)
图 1-3 的 R 程序如下:
acf(s3,lag=length(s3)-1,xlab="",ylab="",axes=FALSE,mian="")
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)
图 1-2 1700 年至 2016 年太阳黑子数量条形图
图 1-3 太阳黑子数据的自相关函数图
图 1-4 太阳黑子数据的偏自相关函数图
由于,研究人员通过对太阳黑子数量的长期观测可知其周期性为 11 年左右,
故我们对太阳黑子数据进行 11 阶季节差分,得到差分后的数据记为 4s ,得到差分
后的时序图 1-4:
图 1-4R 程序如下:
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)
图 1-4 11 阶差分后太阳黑子数量的时序图
从上图可以看出,进行季节差分后的数据 4s 基本平稳,对 4s 进行单位根检验
(ADF 检验)R 程序如下:
s4=diff(s3,11)
library(fUnitRoots)
adfTest(s4,lags=11,type='ct')
得到结果为:
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 11
STATISTIC:
Dickey-Fuller: -3.7971
P VALUE:
0.01954
可 以 看 到 , 在 显 著 性 水 平 0.05
的 条 件 下 , 检 验 统 计 量 所 对 应 的
,故拒绝序列具有单位根的原假设,所以我们认为序列 4s 是
p
值
=0.01954 0.05
平稳的序列。
二、利用 R 软件对平稳时间序列建立模型
对原始的太阳黑子数据已经进行消除周期性的处理得到平稳的序列 4s ,故对
平稳时间序列建立模型:
图 1-5 季节差分后太阳黑子数据的自相关函数图
图 1-6 季节差分后太阳黑子数据的偏自相关函数图
R 程序如下:
library(forecast)
m1=auto.arima(s3)
m1
得到结果:
Series: s3
ARIMA(2,0,1) with non-zero mean
Coefficients:
ar1
1.4763
s.e. 0.0481
ar2
ma1
-0.7693 -0.1857
0.0714
0.0434
mean
79.1840
3.9904
sigma^2 estimated as 659.8:
AIC=2965.99
从结果中可以得到 4s 序列最优模型,即 AIC 值,SC 值最小的模型:
log likelihood=-1477.99
AICc=2966.18
BIC=2984.78
1
0.7693
即:
1 1.4763
B
序列
三、模型检验
11
B X
t
a
t
a
1
1
t
1
B
2
1
B
2
B
2
1
11
B X
a
t
0.1857
a
1
t
t
,其中 ta 为白噪声
利用 R 软件对模型画出残差图并对其进行 Box-Ljung 检验,R 程序如下:
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:2016,level,lty=2)
Box.test(m1$residuals,type='Ljung')
结果如下:
图 1-7 季节
2AR 模型的残差图
Box-Ljung test
data: m1$residuals
X-squared = 0.013295, df = 1, p-value = 0.9082
从残差图可以看出,大部分残差数据均匀地分布在 0 的两侧,对其进行
,其对应的 P 值 为 0.9082,
Box-Ljung 检验的结果表明卡方统计量 2=0.013295
的条件下不能拒绝原假设,即认为残差序列是独立的,故
在显著性水平 =0.05
模型是适合的。
四、模型预测
利用 R 软件进行模型预测,预测 2017 年至 2020 年的太阳黑子数量,得到
结果:
R 程序如下:
predict(m1,4)
plot(forecast(m1,h=4))
结果如下:
$pred
Time Series:
Start = 2017
End = 2020
Frequency = 1
[1] 29.68974 36.41375 54.11751 75.08042
$se
Time Series:
Start = 2017
End = 2020
Frequency = 1
[1] 25.68699 41.93792 51.09031 54.02872
从结果中,我们可以得到 2017 年至 2020 年太阳黑子数量的预测值以及其标
准差,可以得到表 1-6:
表 1-6 2017 年-2020 年太阳黑子数量预测值及标准差
年份
预测值
标准差
画出加上预测值的太阳黑子数量时序图,如图 1-8 所示:
2017 年
29.68974
25.68699% 41.93792% 51.09031% 54.02872%
2020 年
75.08042
2018 年
36.41375
2019 年
54.11751
图 1-8 模型预测图