太陽黑子數(shù)的ARIMA模型

廖海燕

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

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

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

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

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

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

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

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

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

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

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

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

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

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)性，需要對季節(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相對應(yīng)，矩陣的列與MA的階q相對應(yīng)。一般來說，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，各系數(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相對應(yīng)。

為了對模型的殘差進(jì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混成檢驗，得到p值為0.882 8，由此可知，模型的殘差的各階相關(guān)函數(shù)為0.

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

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

。

6 結(jié)束語

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

參考文獻(xiàn)

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

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

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

〔編輯：王霞〕endprint

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

矩陣的行與AR的階p相對應(yīng)，矩陣的列與MA的階q相對應(yīng)。一般來說，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，各系數(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相對應(yīng)。

為了對模型的殘差進(jì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混成檢驗，得到p值為0.882 8，由此可知，模型的殘差的各階相關(guān)函數(shù)為0.

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

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

。

6 結(jié)束語

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

參考文獻(xiàn)

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

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

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

〔編輯：王霞〕endprint

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

矩陣的行與AR的階p相對應(yīng)，矩陣的列與MA的階q相對應(yīng)。一般來說，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，各系數(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相對應(yīng)。

為了對模型的殘差進(jì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混成檢驗，得到p值為0.882 8，由此可知，模型的殘差的各階相關(guān)函數(shù)為0.

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

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

。

6 結(jié)束語

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

參考文獻(xiàn)

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

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

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

〔編輯：王霞〕endprint