太陽黑子數的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