廖海燕

摘 要：從比利時(shí)皇家天文臺(tái)（the Royal Observatory of Belgium）的太陽(yáng)黑子指數(shù)數(shù)據(jù)中心（the Sunspot Index Data center）的網(wǎng)站獲得了1700—2013每年的太陽(yáng)黑子數(shù)的數(shù)據(jù)。利用R軟件結(jié)合時(shí)間序列建模方法對(duì)觀測(cè)值進(jìn)行了分析和建模，并利用該模型對(duì)未來(lái)的太陽(yáng)黑子數(shù)進(jìn)行了預(yù)測(cè)。

關(guān)鍵詞：太陽(yáng)黑子數(shù)；R軟件；fBasics添加包；ARIMA模型

中圖分類號(hào)：P182.41 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：2095-6835（2014）11-0121-03

太陽(yáng)黑子是太陽(yáng)光球上出現(xiàn)的一種臨時(shí)現(xiàn)象，它們?cè)诳梢?jiàn)光下會(huì)比周圍區(qū)域黑暗。太陽(yáng)黑子很少單獨(dú)活動(dòng)，而是經(jīng)常成群出現(xiàn)，活躍時(shí)會(huì)對(duì)地球的磁場(chǎng)產(chǎn)生影響，對(duì)各類電子產(chǎn)品和電器造成損害，還會(huì)使地球南北極和赤道的大氣環(huán)流作經(jīng)向流動(dòng)，進(jìn)而造成惡劣天氣。鑒于太陽(yáng)黑子的危害性，有必要對(duì)其數(shù)量進(jìn)行預(yù)測(cè)，進(jìn)而對(duì)其進(jìn)行有效的控制。

從比利時(shí)皇家天文臺(tái)的太陽(yáng)黑子指數(shù)數(shù)據(jù)中心網(wǎng)站獲得從1700—2013期間每年的太陽(yáng)黑子數(shù)數(shù)據(jù)，時(shí)間跨度為314年，即共獲得了314個(gè)觀測(cè)值。下面將利用R軟件結(jié)合時(shí)間序列建模的方法對(duì)這些數(shù)據(jù)進(jìn)行分析，并建立ARIMA預(yù)測(cè)模型。

1 數(shù)據(jù)的預(yù)處理

一共有314個(gè)數(shù)據(jù)，第一個(gè)數(shù)據(jù)是2013年的太陽(yáng)黑子數(shù)，最后一個(gè)數(shù)據(jù)是1700年的太陽(yáng)黑子數(shù)。為了方便時(shí)間序列分析，需要把數(shù)據(jù)的順序倒過(guò)來(lái)。

由于是年度數(shù)據(jù)，所以每個(gè)觀測(cè)值中的日期屬性并不重要，只需提取觀測(cè)值即可，輸入代碼：s1=sunspot2014[，2]，得到命名為s1的去日期數(shù)據(jù)。

將按照時(shí)間順序排列的太陽(yáng)黑子數(shù)的數(shù)據(jù)對(duì)象命名為s3.

2 基本統(tǒng)計(jì)分析

輸入代碼：library（fBasics），加載R的添加包fBasics。為了對(duì)常用統(tǒng)計(jì)量進(jìn)行計(jì)算，輸入代碼：basicStats（s3），由輸出結(jié)果可以得到以下幾個(gè)結(jié)論：①在314個(gè)數(shù)據(jù)中，無(wú)缺失值，且記錄的最小太陽(yáng)黑子數(shù)數(shù)值為0.輸入代碼：which（s3==0）+1700-1，由輸出結(jié)果可知，1711年、1712 年和1810年記錄的太陽(yáng)黑子數(shù)為0.②記錄的最大太陽(yáng)黑子數(shù)數(shù)值為190.2，輸入代碼：which（s3==190.2）+1700-1，由輸出可知，記錄的太陽(yáng)黑子數(shù)的最大值是1957年。③大約1/4的數(shù)據(jù)小于16，而1/4的數(shù)據(jù)大于69，平均每年的太陽(yáng)黑子數(shù)為49.5，標(biāo)準(zhǔn)差是40.3.④輸入代碼which（s3==40）+1700-1，由輸出結(jié)果可知，1725年、1741年和1747年的太陽(yáng)黑子數(shù)是中位數(shù)40，除了這三年的數(shù)據(jù)外，將其余數(shù)據(jù)以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）

為了觀察數(shù)據(jù)的自相關(guān)函數(shù)，輸入代碼：

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）

太陽(yáng)黑子數(shù)具有非常明顯的周期性，且周期在11～12之間，為此，對(duì)時(shí)間序列數(shù)據(jù)進(jìn)行周期為11的季節(jié)差分，輸入代碼：s4=diff（s3，11）。為了檢測(cè)季節(jié)差分后數(shù)據(jù)的平穩(wěn)性，輸入代碼：library（fUnitRoots），以加載R的fUnitRoots添加包，再輸入代碼：adfTest（s4，lags=11，type='ct'），并用擴(kuò)展的Dickey-Fuller單位根進(jìn)行檢驗(yàn)，結(jié)果顯示p值為0.030 47，由此可得出季節(jié)差分后的數(shù)據(jù)具有平穩(wěn)性的結(jié)論。

為了觀察季節(jié)差分后數(shù)據(jù)的變化趨勢(shì)，輸入代碼：

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）

為了觀察季節(jié)差分后數(shù)據(jù)的自相關(guān)函數(shù)，輸入代碼：

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）

為了進(jìn)一步消除數(shù)據(jù)序列的前后相關(guān)性，需要對(duì)季節(jié)差分后的數(shù)據(jù)進(jìn)行常規(guī)差分，輸入代碼：s5=diff（s4），得到被命名為s5的新數(shù)據(jù)序列。endprint

為了確定ARIMA模型的階，使用擴(kuò)展的自相關(guān)函數(shù)EACF，為此，輸入代碼：library（TSA），以加載R的TSA添加包，然后輸入代碼：eacf（s5，11，11），輸出的矩陣圖如圖5所示。

矩陣的行與AR的階p相對(duì)應(yīng)，矩陣的列與MA的階q相對(duì)應(yīng)。一般來(lái)說(shuō)，ARIMA模型的階（p，q）位于“o”組成的三角形的左上角頂點(diǎn)處。

4 模型的建立和診斷

輸入代碼：library（forecast），以加載R的forecast添加包。為了建立AIC最佳的ARIMA模型，輸入代碼：（m1=auto.arima（s3）），建立的模型如下（對(duì)象名為m1）：

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

模型的估計(jì)值為237.6，各系數(shù)的標(biāo)準(zhǔn)差依次為：0.084 2，0.130 9，0.077 2，0.059 6，0.059 5，模型的AIC為2 615.51，模型的階的選定要與圖5相對(duì)應(yīng)。

為了對(duì)模型的殘差進(jìn)行檢驗(yàn)，輸入代碼：

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）

該模型的殘差分布狀態(tài)比較理想。

輸入代碼：Box.test（m1$residuals，type=Ljung），進(jìn)行Box-Ljung混成檢驗(yàn)，得到p值為0.882 8，由此可知，模型的殘差的各階相關(guān)函數(shù)為0.

5 樣本外的四步預(yù)測(cè)

對(duì)建立的模型m1進(jìn)行樣本外的四步預(yù)測(cè)，即預(yù)測(cè)2014—2017期間的太陽(yáng)黑子數(shù)。輸入代碼：predict（m1，4），得到2014—2017太陽(yáng)黑子數(shù)的預(yù)測(cè)值及其標(biāo)準(zhǔn)差，

。

6 結(jié)束語(yǔ)

可以對(duì)太陽(yáng)黑子數(shù)的ARIMA模型進(jìn)行進(jìn)一步的改進(jìn)，使對(duì)太陽(yáng)黑子歷史數(shù)據(jù)分析和對(duì)未來(lái)太陽(yáng)黑子數(shù)的預(yù)測(cè)更加科學(xué)和準(zhǔn)確，為相關(guān)部門制訂安全的控制措施提供參考依據(jù)。

參考文獻(xiàn)

[1]李素蘭.數(shù)據(jù)分析與R軟件[M].北京：科學(xué)出版社，2013.

[2]Paul Teetor.R語(yǔ)言經(jīng)典實(shí)例[M].北京：機(jī)械工業(yè)出版社，2013.

[3]Jonathan D.Cryer，Kung-Sik Chan.時(shí)間序列分析及應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2011.

〔編輯：王霞〕endprint

為了確定ARIMA模型的階，使用擴(kuò)展的自相關(guān)函數(shù)EACF，為此，輸入代碼：library（TSA），以加載R的TSA添加包，然后輸入代碼：eacf（s5，11，11），輸出的矩陣圖如圖5所示。

矩陣的行與AR的階p相對(duì)應(yīng)，矩陣的列與MA的階q相對(duì)應(yīng)。一般來(lái)說(shuō)，ARIMA模型的階（p，q）位于“o”組成的三角形的左上角頂點(diǎn)處。

4 模型的建立和診斷

輸入代碼：library（forecast），以加載R的forecast添加包。為了建立AIC最佳的ARIMA模型，輸入代碼：（m1=auto.arima（s3）），建立的模型如下（對(duì)象名為m1）：

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

模型的估計(jì)值為237.6，各系數(shù)的標(biāo)準(zhǔn)差依次為：0.084 2，0.130 9，0.077 2，0.059 6，0.059 5，模型的AIC為2 615.51，模型的階的選定要與圖5相對(duì)應(yīng)。

為了對(duì)模型的殘差進(jìn)行檢驗(yàn)，輸入代碼：

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）

該模型的殘差分布狀態(tài)比較理想。

輸入代碼：Box.test（m1$residuals，type=Ljung），進(jìn)行Box-Ljung混成檢驗(yàn)，得到p值為0.882 8，由此可知，模型的殘差的各階相關(guān)函數(shù)為0.

5 樣本外的四步預(yù)測(cè)

對(duì)建立的模型m1進(jìn)行樣本外的四步預(yù)測(cè)，即預(yù)測(cè)2014—2017期間的太陽(yáng)黑子數(shù)。輸入代碼：predict（m1，4），得到2014—2017太陽(yáng)黑子數(shù)的預(yù)測(cè)值及其標(biāo)準(zhǔn)差，

。

6 結(jié)束語(yǔ)

可以對(duì)太陽(yáng)黑子數(shù)的ARIMA模型進(jìn)行進(jìn)一步的改進(jìn)，使對(duì)太陽(yáng)黑子歷史數(shù)據(jù)分析和對(duì)未來(lái)太陽(yáng)黑子數(shù)的預(yù)測(cè)更加科學(xué)和準(zhǔn)確，為相關(guān)部門制訂安全的控制措施提供參考依據(jù)。

參考文獻(xiàn)

[1]李素蘭.數(shù)據(jù)分析與R軟件[M].北京：科學(xué)出版社，2013.

[2]Paul Teetor.R語(yǔ)言經(jīng)典實(shí)例[M].北京：機(jī)械工業(yè)出版社，2013.

[3]Jonathan D.Cryer，Kung-Sik Chan.時(shí)間序列分析及應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2011.

〔編輯：王霞〕endprint

為了確定ARIMA模型的階，使用擴(kuò)展的自相關(guān)函數(shù)EACF，為此，輸入代碼：library（TSA），以加載R的TSA添加包，然后輸入代碼：eacf（s5，11，11），輸出的矩陣圖如圖5所示。

矩陣的行與AR的階p相對(duì)應(yīng)，矩陣的列與MA的階q相對(duì)應(yīng)。一般來(lái)說(shuō)，ARIMA模型的階（p，q）位于“o”組成的三角形的左上角頂點(diǎn)處。

4 模型的建立和診斷

輸入代碼：library（forecast），以加載R的forecast添加包。為了建立AIC最佳的ARIMA模型，輸入代碼：（m1=auto.arima（s3）），建立的模型如下（對(duì)象名為m1）：

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

模型的估計(jì)值為237.6，各系數(shù)的標(biāo)準(zhǔn)差依次為：0.084 2，0.130 9，0.077 2，0.059 6，0.059 5，模型的AIC為2 615.51，模型的階的選定要與圖5相對(duì)應(yīng)。

為了對(duì)模型的殘差進(jìn)行檢驗(yàn)，輸入代碼：

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）

該模型的殘差分布狀態(tài)比較理想。

輸入代碼：Box.test（m1$residuals，type=Ljung），進(jìn)行Box-Ljung混成檢驗(yàn)，得到p值為0.882 8，由此可知，模型的殘差的各階相關(guān)函數(shù)為0.

5 樣本外的四步預(yù)測(cè)

對(duì)建立的模型m1進(jìn)行樣本外的四步預(yù)測(cè)，即預(yù)測(cè)2014—2017期間的太陽(yáng)黑子數(shù)。輸入代碼：predict（m1，4），得到2014—2017太陽(yáng)黑子數(shù)的預(yù)測(cè)值及其標(biāo)準(zhǔn)差，

。

6 結(jié)束語(yǔ)

可以對(duì)太陽(yáng)黑子數(shù)的ARIMA模型進(jìn)行進(jìn)一步的改進(jìn)，使對(duì)太陽(yáng)黑子歷史數(shù)據(jù)分析和對(duì)未來(lái)太陽(yáng)黑子數(shù)的預(yù)測(cè)更加科學(xué)和準(zhǔn)確，為相關(guān)部門制訂安全的控制措施提供參考依據(jù)。

參考文獻(xiàn)

[1]李素蘭.數(shù)據(jù)分析與R軟件[M].北京：科學(xué)出版社，2013.

[2]Paul Teetor.R語(yǔ)言經(jīng)典實(shí)例[M].北京：機(jī)械工業(yè)出版社，2013.

[3]Jonathan D.Cryer，Kung-Sik Chan.時(shí)間序列分析及應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2011.

〔編輯：王霞〕endprint