太阳黑子数量的预测

太阳黑子数量的预测

这篇博客是学习zonination/sunspots的笔记。展示了每日太阳表面黑子的数量。最后我还尝试对每年平均的太阳黑子数进行了建模,验证了太阳黑子的周期确实大致为11年。

绘图

该图表数据来自:http://www.sidc.be/silso/INFO/sndtotcsv.php,数据实时更新。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
download.file("http://www.sidc.be/silso/INFO/sndtotcsv.php", "每日黑子的数量.csv")

library(readr)
library(ggplot2)
spots <- read_delim("每日黑子的数量.csv", ";", escape_double = F, col_names = F, trim_ws = T)

# 列名分别是:
# 年
# 月
# 日
# 小数日期
# 每日太阳黑子数
# 标准偏差
# 观察次数
# 确定/临时指标
names(spots) <- c("y", "m", "d", "y.frac", "ssn", "sd", "n.obs", "def")

# 把-1换成缺失值
spots$ssn <- ifelse(spots$ssn == -1, NA, spots$ssn)
library(cowplot)

ggplot(spots, aes(y.frac, ssn)) +
geom_point(alpha = 0.1, size = 0.01, colour = 'orangered') +
scale_x_continuous(breaks = seq(1800, 2100, 50),
minor_breaks = seq(1800, 2100, 10)) +
scale_y_continuous(breaks = seq(0, 600, 100),
minor_breaks = seq(0, 600, 20)) +
labs(title = "图:每日太阳黑子数量",
x = "日期", y = "太阳黑子数量",
subtitle = "太阳表面黑子数量",
caption = "数据来源:http://www.sidc.be") +
theme_bw(base_size = 20, base_family = "STSong") +
theme(plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm")) +
theme(plot.title = element_text(hjust = 0.5)) + theme(plot.subtitle = element_text(hjust = 0.5))

从长期来看,太阳黑子的数量确实有很强的周期性,记得高中地理书上说的是11年为周期。另外一个感慨的是,西方人做记录好认真,这个数据从1818年开始的。。。

建模

由于有些年份的缺失值很多,可以把这些年份找出来删除掉。下面的代码是为了统计每年的观测值的数目:

1
2
3
library(dplyr)
year_count <- summarise(group_by(spots[which(!is.na(spots$ssn)),], y), count = length(ssn))
plot(year_count)

可以看到1849年之后的观测值基本稳定在365和366上,不再有缺失值,所以仅取用1849-2017年的观测数据进行建模:

1
2
3
4
5
d1 <- subset(spots, spots$y >= 1849)
d1 <- subset(d1, d1$y != 2018)
d1 <- subset(d1, !is.na(d1$ssn))
# 计算每年的平均值
d1 <- summarise(group_by(d1, y), mean = mean(ssn))

平均后的序列是这样的:

1
2
ggplot(data = d1, aes(x = y, y = mean)) +
geom_line()

显然是平稳的,下面再对该序列建立一个季节ARMA模型,以往的研究表明太阳黑子的周期是11年,所以取11阶差分进行建模,第一步是定阶:

1
2
3
4
ssn <- d1$mean
ssn = ts(ssn, start = c(1849))
sdssn <- diff(ssn, 11)
acf(sdssn)

ACF图表明MA部分取2阶比较合适,再看AR部分:

1
pacf(sdssn)


同样我们也试试2阶的AR部分。

11阶差分后的序列是这样的:

1
plot(sdssn)

再进行ADF平稳性检验:

1
2
library(fUnitRoots)
adfTest(sdssn)

结果:

1
2
3
4
5
6
7
8
9
Augmented Dickey-Fuller Test

Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -6.0149
P VALUE:
0.01

p值表明这个序列是平稳的。

forecast包的auto.arima()函数运行的结果也是ARMA(2, 2):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
library(forecast)
m1 <- auto.arima(sdssn)
summary(m1)
# 结果
Series: sdssn
ARIMA(2,0,2) with zero mean

Coefficients:
ar1 ar2 ma1 ma2
1.1747 -0.5640 -0.3691 0.2683
s.e. 0.1621 0.1259 0.1803 0.1001

sigma^2 estimated as 896.5: log likelihood=-759.76
AIC=1529.52 AICc=1529.91 BIC=1544.83

Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -1.262856 29.5596 22.03454 5.280506 199.0357 0.8639828 0.01424979

这个模型是这样的:
$$
(1 - 1.17B + 1 + 0.56B^2)sdssn_t = (1 - 0.37B + 0.27B^2)\varepsilon_n
$$

再对模型检验:

1
tsdiag(m1)

检验图由三个子图构成,第一幅是模型拟合后的标准化残差,一般要求不能出现过多的极端值(即要满足正态性假设);第二幅图是残差的自相关图,如果残差是正态的,那么自相关函数应该一直为0(或者不显著);第三幅图是对残差的Ljung-Box检验,散点位于虚线之上表明无法拒绝残差无自相关的原假设。因此这个模型是充分的。

下面为了减小估计误差,直接对原序列(均值序列)拟合季节ARMA模型:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
m2 <- arima(ssn, order = c(2, 0, 2), seasonal = list(order = c(1, 0, 1), period = 11))
summary(m2)

# 结果
Call:
arima(x = ssn, order = c(2, 0, 2), seasonal = list(order = c(1, 0, 1), period = 11))

Coefficients:
ar1 ar2 ma1 ma2 sar1 sma1 intercept
1.4315 -0.7525 -0.4151 0.1652 0.5983 -0.2952 84.5579
s.e. 0.0955 0.0788 0.1219 0.0825 0.1788 0.1963 7.3059

sigma^2 estimated as 607.7: log likelihood = -783.26, aic = 1582.53

Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.2218944 24.65169 18.93943 -41.19816 59.43233 0.6167939
ACF1
Training set -0.01870418

拟合的模型结果为:
$$
(1 - 1.43B + 0.75B^2)(1 - 0.60B^{12})(ssd_t - 84.56) = \\ (1 - 0.42B + 0.16B^2)(1 - 0.30B^{12})ε_n
$$

同样进行模型检验:

1
tsdiag(m2)

同样充分。

预测

例如,做个100年的预测:

1
2
library(forecast)
plot(forecast(m2, h = 100))

如果觉得这个图太丑了,可以直接自己画:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
forecast <- as.data.frame(forecast(m2, h = 100))
forecast$y <- as.numeric(rownames(forecast))
colnames(forecast) <- c('pe', 'lo80', 'hi80', 'lo95', 'hi95', 'year')
ggplot() +
geom_line(data = d1, aes(y, mean)) +
geom_ribbon(data = forecast, aes(x = year, ymin = lo95, ymax = hi95), alpha = 0.8, fill = "#92b4f4") +
geom_ribbon(data = forecast, aes(x = year, ymin = lo80, ymax = hi80), alpha = 0.8, fill = "#5e7ce2") +
geom_line(data = forecast, aes(x = year, y = pe),
colour = '#fc4a1a', linetype = 6) +
scale_x_continuous(limits = c(1849, 2117),
breaks = seq(1850, 2100, 50)) +
scale_y_continuous(limits = c(-70, 300),
breaks = c(-100, 0, 100, 200, 300)) +
theme_bw(base_family = 'STSong', base_size = 20) +
labs(title = '图:未来100年每年太阳黑子平均数量的预测',
caption = '数据来源:http://www.sidc.be',
subtitle = '周期为11年',
x = "年份",
y = "每年的平均数量") +
theme(plot.title = element_text(hjust = 0.5)) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
theme(plot.margin = grid::unit(c(0.5, 0.5, 0.5, 0.5), "cm")) +
theme(panel.grid.minor.y = element_blank()) +
theme(panel.grid.minor.y = element_blank()) +
draw_line(x = c(2017, 2017), y = c(-70, 300), colour = '#fc4a1a') +
draw_label(x = 2022, y = 250, label = '2017年', colour = '#fc4a1a', fontfamily = 'STSong', hjust = 0)

典型的时间序列的预测结果啦,长期预测会趋于平均。但是上面的分析表明黑子数量的周期大致为11年。

# R

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×