廖海燕
摘 要:從比利時皇家天文臺(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