R語言時間序列中文教程教學(xué)文稿

1、Good is good, but better carries it.精益求精，善益求善。R語言時間序列中文教程R語言時間序列中文教程2012特別聲明：R語言是免費(fèi)語言，其代碼不帶任何質(zhì)量保證，使用R語言所產(chǎn)生的后果由使用者負(fù)全責(zé)。前言R語言是一種數(shù)據(jù)分析語言，它是科學(xué)的免費(fèi)的數(shù)據(jù)分析語言，是凝聚了眾多研究人員心血的成熟的使用范圍廣泛全面的語言，也是學(xué)習(xí)者能較快受益的語言。在R語言出現(xiàn)之前，數(shù)據(jù)分析的編程語言是SAS。當(dāng)時SAS的功能比較有限。在貝爾實(shí)驗室里，有一群科學(xué)家討論提到，他們研究過程中需要用到數(shù)據(jù)分析軟件。SAS的局限也限制了他們的研究。于是他們想，我們貝爾實(shí)驗室的研究歷史要比SA

2、S長好幾倍，技術(shù)力量也比SAS強(qiáng)好幾倍，且貝爾實(shí)驗室里并不缺乏訓(xùn)練有素的專業(yè)編程人員，那么，我們貝爾實(shí)驗室為什么不自己編寫數(shù)據(jù)分析語言，來滿足我們應(yīng)用中所需要的特殊要求呢？于是，貝爾實(shí)驗室研究出了S-PLUS語言。后來，新西蘭奧克蘭大學(xué)的兩位教授非常青睞S-PLUS的廣泛性能。他們決定重新編寫與S-PLUS相似的語言，并且使之免費(fèi)，提供給全世界所有相關(guān)研究人員使用。于是，在這兩位教授努力下，一種叫做R的語言在奧克蘭大學(xué)誕生了。R基本上是S-PLUS的翻版，但R是免費(fèi)的語言，所有編程研究人員都可以對R語言做出貢獻(xiàn)，且他們已經(jīng)將大量研究成果寫成了R命令或腳本，因而R語言的功能比較強(qiáng)大，比較全面。研

3、究人員可免費(fèi)使用R語言，可通過閱讀R語言腳本源代碼，學(xué)習(xí)其他人的研究成果。筆者曾有幸在奧克蘭大學(xué)受過幾年熏陶，曾經(jīng)向一位統(tǒng)計系的老師提請教過一個數(shù)據(jù)模擬方面的問題。那位老師只用一行R語句就解答了。R語言的強(qiáng)大功能非常令人驚訝。為了進(jìn)一步推廣R語言，為了方便更多研究人員學(xué)習(xí)使用R語言，我們收集了R語言時間序列分析實(shí)例，以供大家了解和學(xué)習(xí)使用。當(dāng)然，這是非常簡單的模仿練習(xí)，具體操作是，用復(fù)制粘貼把本材料中R代碼放入R的編程環(huán)境；材料中藍(lán)色背景的內(nèi)容是相關(guān)代碼和相應(yīng)輸出結(jié)果。經(jīng)過反復(fù)模仿，學(xué)習(xí)者便能熟悉和學(xué)會。需要提醒學(xué)習(xí)者的是：建議學(xué)習(xí)者安裝了R語言編程，再繼續(xù)閱讀本材料；執(zhí)行R命令時，請刪除命令

4、的中文注解，沒使用過在命令中加入中文；如果學(xué)習(xí)者是初次接觸R或者Splus，建議先閱讀，如果學(xué)習(xí)者比較熟悉R語言，還可以閱讀優(yōu)秀時間序列讀物EcomometricsinR，也可以上QuickR網(wǎng)站。目錄TOCo1-3huHYPERLINKl_Toc15658R語言時間序列分析PAGEREF_Toc156581HYPERLINKl_Toc8339前言PAGEREF_Toc83391HYPERLINKl_Toc18551目錄PAGEREF_Toc185512HYPERLINKl_Toc304941.運(yùn)用R語言研究JJ數(shù)據(jù)PAGEREF_Toc304943HYPERLINKl_Toc251242.運(yùn)

5、用R語言研究空氣污染PAGEREF_Toc2512419HYPERLINKl_Toc159543.運(yùn)用R語言研究自動回歸移動平均集成模型PAGEREF_Toc1595422HYPERLINKl_Toc29364.運(yùn)用R語言研究全球變暖理論P(yáng)AGEREF_Toc293631HYPERLINKl_Toc47715.運(yùn)用R語言研究非獨(dú)立誤差與線性回歸PAGEREF_Toc477132HYPERLINKl_Toc138836.運(yùn)用R語言研究估計波動的頻率PAGEREF_Toc1388336HYPERLINKl_Toc208207.運(yùn)用R語言研究厄爾尼諾頻率PAGEREF_Toc2082039HYPER

6、LINKl_Toc68738.運(yùn)用R語言研究太陽黑子周期頻率PAGEREF_Toc6873401.運(yùn)用R語言研究JJ數(shù)據(jù)學(xué)習(xí)R言語時間序列分析程序操作，需要從最基礎(chǔ)、最簡單的學(xué)起，例如在命令窗口中，輸入并執(zhí)行2+2等于4的R語言命令。2+214執(zhí)行完2+2等于4的R語言命令后，我們可以開始時間序列，即學(xué)著把玩johnson&Johnson數(shù)據(jù)。下載jj.dat或執(zhí)行下面語句。這個數(shù)據(jù)已被人上傳到因特網(wǎng)中。R所需要做的只是將網(wǎng)址進(jìn)行掃描就可以將數(shù)據(jù)讀取進(jìn)入R的編程環(huán)境中。下面有3種不同讀取數(shù)據(jù)的方法：jj=scan(/stoffer/tsa2/data/jj.dat)#readthedata讀取

7、數(shù)據(jù)jjjj#andanother第三種方法讀取數(shù)據(jù)使用R語言的人，有的喜歡使用，大多數(shù)醫(yī)療系統(tǒng)的工作者喜歡用=，正因為如此才用了上面種不同讀取數(shù)據(jù)的方法。讀取數(shù)據(jù)后，鍵入并執(zhí)行jj，數(shù)據(jù)在窗口便會有如下顯示：jj10.710.630.850.4450.610.690.920.557714.0412.9614.859.998116.2014.6716.0211.61jj中有84個數(shù)據(jù)被稱作對象。下面命令可以顯示所有對象。objects()如果使用matlab，你會認(rèn)為jj是一個84行1列的向量，但實(shí)際上不是這樣。jj有次序，有長度，但沒維度，R稱這些對象為向量，要小心區(qū)別。在R里，矩陣有維度，

8、但向量沒維度。這都是程序語言的一些概念。jj1#thefirstelement列中第一個數(shù)據(jù)10.71jj84#thelastelement列中最后一個數(shù)據(jù)111.61jj1:4#thefirst4elements列中第一至第四個數(shù)據(jù)10.710.630.850.44jj-(1:80)#everythingEXCEPTthefirst80elements列中除第80個以外的所有數(shù)據(jù)116.2014.6716.0211.61length(jj)#thenumberofelements有多少個數(shù)據(jù)184dim(jj)#butnodimensions.但沒維度NULLnrow(jj)#.norows

9、沒行NULLncol(jj)#.andnocolumns沒列NULL#如果你要把jj轉(zhuǎn)變?yōu)橐粋€向量，執(zhí)行如下操作后，維度為84行1列。jj=as.matrix(jj)dim(jj)1841然后把jj轉(zhuǎn)變?yōu)橐粋€時間序列對象。jj=ts(jj,start=1960,frequency=4)#ts()命令這個數(shù)據(jù)是從1960年開始，個個季度的收入，frequency=4指四個季度。R語言的優(yōu)勢在于可用一條命令做很多事，即可以把前面的命令放在一起打包執(zhí)行。其操作如下：jj=ts(scan(/stoffer/tsa2/data/jj.dat),start=1960,frequency=4)在上面命令里，

10、scan可以被read.table替代。用read.table讀取數(shù)據(jù)可生成matrix對象，還可以給每列起名字。下面學(xué)習(xí)一下read.table,dataframes,和時間序列對象。輸入命令后，窗口會有如下顯示：jj=ts(read.table(/stoffer/tsa2/data/jj.dat),start=1960,frequency=4)help(read.table)help(ts)help(data.frame)需要注意的是，Scan和read.table不一樣。Scan生成的是有維度的向量，read.table生成的則是帶有維度的數(shù)據(jù)架構(gòu)。讀取jj數(shù)據(jù)的最后要領(lǐng)。如果這個數(shù)據(jù)是

11、從1960年第三個季度開始的，所需輸入命令則為ts(x,start=c(1960,3),frequency=4)；如果是一個每月每月的數(shù)據(jù)，例如數(shù)據(jù)是從1984年6月開始的，需要輸入的命令則為ts(x,start=c(1984,6),frequency=12)。輸入命令后，轉(zhuǎn)變后的時間序列對象為：jjQtr1Qtr2Qtr3Qtr419600.710.630.850.4419610.610.690.920.55197914.0412.9614.859.99198016.2014.6716.0211.61注意到區(qū)別了嗎？時間信息，也就是4個不同的季度的數(shù)據(jù)被加載到里面了。進(jìn)行時間數(shù)據(jù)分析后，窗口

12、會有如下顯示：time(jj)Qtr1Qtr2Qtr3Qtr419601960.001960.251960.501960.7519611961.001961.251961.501961.7519791979.001979.251979.501979.7519801980.001980.251980.501980.75接下來輸入如下組合命令。(jj=ts(scan(/stoffer/tsa2/data/jj.dat),start=1960,frequency=4)然后進(jìn)行對數(shù)據(jù)繪圖：plot(jj,ylab=EarningsperShare,main=J&J)輸入以上命令后，可以看到如下結(jié)果：再

13、輸入下面的命令，看看區(qū)別。plot(jj,type=o,col=blue,lty=dashed)plot(diff(log(jj),main=loggedanddiffed)下面利用操作plot.ts和ts.plot兩個相關(guān)命令，顯示區(qū)別。x=-5:5#sequenceofintegersfrom-5to5y=5*cos(x)#guesspar(mfrow=c(3,2)#multifiguresetup:3rows,2cols#plot:plot(x,main=plot(x)plot(x,y,main=plot(x,y)#plot.ts:plot.ts(x,main=plot.ts(x)plo

14、t.ts(x,y,main=plot.ts(x,y)#ts.plot:ts.plot(x,main=ts.plot(x)ts.plot(ts(x),ts(y),col=1:2,main=ts.plot(x,y)#note-xandyaretsobjects#thehelpfiles?andhelp()arethesame:?plot.tshelp(ts.plot)?par#mightaswellskimthegraphicalparametershelpfilewhileyourehere從窗口中的顯示可以看出，如果數(shù)據(jù)是時間序列對象，使用plot()命令就足夠了；如果數(shù)據(jù)是平常序列，使用pl

15、ot.ts()也可以做時間繪圖。不過，把jj數(shù)據(jù)放在一張圖上，數(shù)據(jù)會隨著時間的變化上上下下跳動，能從整體上反應(yīng)上升或者下降的趨勢。上文中用紅色光滑的曲線代表上升的趨勢，簡單明了。這需要將過濾和光滑的技巧使用在jj數(shù)據(jù)上。在這里，我們用對稱的移動平均值來達(dá)到過濾和光滑的目的。下面使用公式：fjj(t)=jj(t-2)+jj(t-1)+jj(t)+jj(t+1)+jj(t+2)除此之外，lowess的過濾平滑技巧（藍(lán)色曲線）也要使用在jj數(shù)據(jù)中。具體操作如下圖：k=c(.5,1,1,1,.5)#kisthevectorofweights用于對稱移動平均的系數(shù)(k=k/sum(k)10.1250.2

16、500.2500.2500.125fjj=filter(jj,sides=2,k)#?filterforhelpbutyouknewthatalready使用對稱移動平均plot(jj)lines(fjj,col=red)#addsalinetotheexistingplot稱移動平均的繪圖lines(lowess(jj),col=blue,lty=dashed)#lowess的繪圖操作后，窗口會顯示下面結(jié)果：看完jj數(shù)據(jù)，我們就需要開始具體分析。第一步，我們把所有jj數(shù)據(jù)都取log值。第二步，我們把log值做差，即使用log值數(shù)列中第二值減去第一值，第三值減去第二值，第四值減去第三值等等。如

17、果做差處理前數(shù)列里有n個數(shù)值，處理后的結(jié)果中將有n-1個數(shù)值。dljj=diff(log(jj)#differencetheloggeddata做log和差的處理plot(dljj)#plotitifyouhaventalready對結(jié)果制圖shapiro.test(dljj)#testfornormality測試結(jié)果的正態(tài)分布的性質(zhì)Shapiro-Wilknormalitytestdata:dljjW=0.9725,p-value=0.07211處理完畢以上兩步，我們接下來就要將柱形分布圖和QQ圖放在一起。這兩個圖的本質(zhì)仍舊在于測試數(shù)據(jù)正態(tài)分布的性質(zhì)。數(shù)據(jù)正態(tài)分布的性質(zhì)是整個統(tǒng)計學(xué)構(gòu)架堅實(shí)的

18、基礎(chǔ)，如果這個性質(zhì)的存在比較可信、通過了所有測試，統(tǒng)計分析中得出的結(jié)論就比較可信、就通得過所有測試。當(dāng)然如果這個性質(zhì)在數(shù)據(jù)中不存在，我們需要用其它的技巧來處理。詳細(xì)的，參看R語言樣品比較應(yīng)用的實(shí)例。以上操作，在窗口中有如下顯示：par(mfrow=c(2,1)#setupthegraphics設(shè)置為兩圖的輸出hist(dljj,prob=TRUE,12)#histogram柱形分布圖lines(density(dljj)#smoothit-?densityfordetails柱形分布圖的曲線qqnorm(dljj)#normalQ-QplotQQ圖qqline(dljj)#addaline在Q

19、Q圖上加直線經(jīng)過測試數(shù)據(jù)后，窗口會有如下顯示：在實(shí)踐操作中，時間序列數(shù)據(jù)存在著前后關(guān)系。例如，今天股票的價格很有可能決定明天股票的價格。明天的溫度取決于今天的氣溫。做天氣預(yù)報的具體操作方法，是使用已經(jīng)存在的天氣歷史記錄，比如說今天的氣溫，昨天的氣溫，前天的氣溫等等，來預(yù)測明天的氣溫。當(dāng)然，在進(jìn)行預(yù)測之前，我們一定要看清時間序列數(shù)據(jù)中的前后關(guān)系結(jié)構(gòu)，清楚哪一個特定的歷史數(shù)據(jù)可以精確預(yù)測未來的數(shù)據(jù)。在這里，我們使用被log和求差后的dljj數(shù)據(jù)，來介紹分析時間序列數(shù)據(jù)前后關(guān)系結(jié)構(gòu)的具體技巧。在預(yù)測的實(shí)際應(yīng)用中，我們總希望用歷史數(shù)據(jù)來預(yù)測即將要產(chǎn)生的數(shù)據(jù)。事實(shí)上，已產(chǎn)生的數(shù)據(jù)相對于即將產(chǎn)生的數(shù)據(jù)，中

20、間存在著一定的延遲，也就是lag.比方說在某天凌晨12點(diǎn)開始記錄室內(nèi)溫度，每小時記一次，一共連續(xù)記錄了10個小時。凌晨12點(diǎn)的數(shù)據(jù)和凌晨3點(diǎn)的數(shù)據(jù)之間就存在著延遲。12點(diǎn)的數(shù)據(jù)比3點(diǎn)的早了3個小時，可記作lag3.3點(diǎn)的數(shù)據(jù)比12點(diǎn)的晚了3個小時，可記作lead3.我們下面來介紹關(guān)聯(lián)性。例如，冷飲的銷量與天氣溫度存在關(guān)聯(lián)性。溫度越高冷飲銷量就越高，溫度越低冷飲銷量也越低。這種關(guān)聯(lián)性稱為正面關(guān)聯(lián)性。又如，人的體重和跑步速度也存在關(guān)聯(lián)性。不過，人的體重越重，跑步速度就會越慢，體重越低，相對來講，速度就會越快。這種關(guān)聯(lián)性稱為負(fù)面關(guān)聯(lián)性。下面我們回到預(yù)測應(yīng)用上。如果現(xiàn)在收集的數(shù)據(jù)與將來的數(shù)據(jù)之間存在著

21、正面或者是負(fù)面的關(guān)聯(lián)性，我們就可以用現(xiàn)在收集的數(shù)據(jù)來預(yù)測未來的數(shù)據(jù)。因此找到現(xiàn)在收集到的數(shù)據(jù)與未來數(shù)據(jù)之間的關(guān)聯(lián)性是最關(guān)鍵的。找到這種關(guān)聯(lián)性的具體技巧被稱作延遲圖表，也就是lag.plot.我們可以在電腦里輸入如下命令：lag.plot(dljj,9,do.lines=FALSE)#whythedo.lines=FALSE?.tryleavingitout上面語句顯示了lag.plot用法，輸入的數(shù)據(jù)是被log和作差后的jj.dat數(shù)據(jù)。其中9表示我們要考慮從1到9這9個不同的延遲。值得注意的是，在下面圖表中顯示出延遲4和8顯示出了正面關(guān)系。其他幾個延遲存在著負(fù)面關(guān)系。我們可以利用這4和8的正

22、面關(guān)系來進(jìn)行預(yù)測，即用現(xiàn)有數(shù)據(jù)計算接下來的第4個或者是第8個數(shù)據(jù)。下面我們來看ACF和PACF圖表。確定我們已經(jīng)觀察到的正面和負(fù)面關(guān)系。par(mfrow=c(2,1)#Thepowerofaccurateobservationiscommonlycalledcynicism#bythosewhohavenotgotit.-GeorgeBernardShawacf(dljj,20)#ACFtolag20-nographshown.keepreadingpacf(dljj,20)#PACFtolag20-nographshown.keepreading#!NOTE!acf2onthelinebe

23、lowisNOTavailableinR.detailsfollowthegraphbelowacf2(dljj)#thisiswhatyoullseebelow在上圖中，ACF和PACF橫坐標(biāo)的標(biāo)記是1,2,3,4,5.但因為數(shù)據(jù)是季度性的，每年有4個季度所以1,2,3,4,5的標(biāo)記代表的4,8,12,16,20的延遲。當(dāng)然，如果我們不喜歡上面橫坐標(biāo)的標(biāo)記，也可以將dljj更改為ts(dljj,freq=1);也就是說acf(ts(dljj,freq=1),20)。因為在上面ACF圖表中橫坐標(biāo)1代表的是延遲4，橫坐標(biāo)2代表的是延遲6，其相應(yīng)的豎線代表的就是延遲4和8的正面關(guān)系。接下來，下面我

24、們介紹結(jié)構(gòu)拆析。在前面R代碼中，我們曾將所有jj數(shù)據(jù)進(jìn)行了lag變型。在變型后的數(shù)據(jù)中，存在著上升趨勢，季節(jié)的影響和每一時間點(diǎn)產(chǎn)生的隨機(jī)的誤差。根據(jù)這一數(shù)據(jù)圖，我們能夠把趨勢、季節(jié)和誤差從變型后的jj數(shù)據(jù)中拆析出來，分別研究，或者分別進(jìn)行繪圖，以便于單獨(dú)檢查。下面等式代表將要使用的數(shù)學(xué)模型：Log(jj)=趨勢+季節(jié)+誤差log(jj)=trend+season+error結(jié)構(gòu)拆析的R命令是stl(),下面語句中stl命令中輸入的是lag變型后的jj數(shù)據(jù)。其中的“per”輸入指的是使用季節(jié)循環(huán)來進(jìn)行拆析。stl語句在這里生成了一個叫dog的R物件，然后Plot語句將dog物件進(jìn)行繪圖。具體操作

25、如下圖”plot(dog|t|)trend0.1671720.00225974.002e-16*Q11.0527930.02735938.482e-16*Q21.0809160.02736539.502e-16*Q31.1510240.02738342.032e-16*Q40.8822660.02741232.192e-16*Signif.codes:0*0.001*0.01*0.05.0.11Residualstandarderror:0.1254on79degreesoffreedomMultipleR-squared:0.9935,AdjustedR-squared:0.9931F-st

26、atistic:2407on5and79DF,p-value:|t|)(Intercept)69.010201.3749850.1902e-16*part0.151400.028985.2252.56e-07*part40.262970.028999.0712e-16*Signif.codes:0*0.001*0.01*0.05.0.11Residualstandarderror:8.323on501degreesoffreedomMultipleR-Squared:0.3091,AdjustedR-squared:0.3063F-statistic:112.1on2and501DF,p-va

27、lue:2.2e-16如上圖所示，我們便獲得了所研究需要的數(shù)據(jù)。3.運(yùn)用R語言研究自動回歸移動平均集成模型現(xiàn)在來介紹自動回歸移動平均集成模型，即ARIMA模型。AR指自動回歸，MA指移動平均，I指集成。我們先來生成兩個隨機(jī)的自動回歸數(shù)列，記為x1和x2。它們都是ar1數(shù)列，因為生成過程中使用order=c(1,0,0)的參數(shù)輸入。X1的AR參數(shù)是+0.9，x2的AR參數(shù)是-0.9。兩個語句生成完畢后，我們以將這兩個數(shù)列繪圖出來。兩個時間序列非常不同，x1的參數(shù)為+0.9，它的圖像類似于股票價格的變動圖形；x2的參數(shù)為-0.9，它的圖像看上去像是一列聲波。窗口具體顯示如下：x1=arima.si

28、m(list(order=c(1,0,0),ar=0.9),n=100)x2=arima.sim(list(order=c(1,0,0),ar=-0.9),n=100)par(mfrow=c(2,1)plot(x1)plot(x2)在R語言中，ACF指自動關(guān)系性，ACF即Auto-CorrelationFunction的簡稱.比方說，股票價格今天的價格跟昨天的價格有關(guān)系，明天的價格會跟今天的或者昨天的價格有關(guān)系。它們之間的關(guān)系性便用ACF來衡量。PACF被稱作不完全自動關(guān)系性。自動關(guān)系性ACF中存在著線性關(guān)系性和非線性關(guān)系性。不完全自動關(guān)系性就是把線性關(guān)系性從自動關(guān)系性中消除。如果在線性關(guān)系性

29、被去除以后，兩個時間點(diǎn)之間的關(guān)系性也就是不完全關(guān)系性。當(dāng)PACF近似于0，這表明兩個時間點(diǎn)之間的關(guān)系性是完全由線性關(guān)系性所造成的。如果不完全關(guān)系性在兩個時間點(diǎn)之間不近似于0，這表明線性模型是不能夠表達(dá)這兩個時間點(diǎn)之間的關(guān)系。我們用下面的語句來檢驗x1和x2數(shù)列中的ACF和PACF。par(mfcol=c(2,2)acf(x1,20)acf(x2,20)pacf(x1,20)pacf(x2,20)對數(shù)列x1來講，ACF從lag1到lag12都高于藍(lán)色虛線，也就是說兩個時間點(diǎn)的距離在1到12之間它們的自動關(guān)系都是正面的。所有線性關(guān)系在x1數(shù)列中被清除了結(jié)果是PartialACF的x1圖表。在PAC

30、Fx1圖表上可以看到只有l(wèi)ag0值為1，其他的lag上的關(guān)系值都低于藍(lán)色虛線，近似于0，也就是說在x1數(shù)列中存在的自動關(guān)系基本上是線性關(guān)系。使用線性模型可以符合x1數(shù)列的要求。對數(shù)列x2來講，兩個時間點(diǎn)之間的ACF關(guān)系有負(fù)面的也有正面的，但經(jīng)過去除線性關(guān)系以后所有不完全自動關(guān)系都接近于0，也就是說x2數(shù)列中的自動關(guān)系也基本上是線性的。使用線性模型符合x2數(shù)列要求。下面，我們來介紹移動平均MA模型。我們先用ARIMA.SIM()命令生成一個移動平均的數(shù)據(jù)列。不同的是，order=c(0,0,1)，ma1代表生成的數(shù)列，ma是一個參數(shù)，設(shè)置為+0.8，數(shù)列中有100個數(shù)據(jù)。其余3行命令是用來設(shè)置輸

31、出窗口繪制圖表和計算ACF，PACF數(shù)值的。經(jīng)過一系列操作后，代碼及其結(jié)果如下所示：#anMA1x=arima.sim(list(order=c(0,0,1),ma=.8),n=100)par(mfcol=c(3,1)plot(x,main=(expression(MA(1)theta=.8)acf(x,20)pacf(x,20)需要注意的是，因數(shù)據(jù)都是隨機(jī)生成的，我們在重新運(yùn)行這些命令時，生成的數(shù)據(jù)會與以前所生成的數(shù)據(jù)完全不同，其結(jié)果也不一樣。除非在以上代碼的前面使用隨機(jī)鎖定的命令set.seed(1)，每次生成的就是相同的數(shù)據(jù)了。當(dāng)然也可以用set.seed(100)，和其它數(shù)字輸入都是可

32、以的。在上圖中，MA1代表了數(shù)列輸出的結(jié)果。最上面顯示的是數(shù)列的數(shù)目，橫向坐標(biāo)從1到100代表有100個數(shù)據(jù)。中間的圖顯示的是x數(shù)列的自動關(guān)系性，在lag0上有很強(qiáng)的自動關(guān)系性，因為Lag0表示沒時間延遲，當(dāng)前數(shù)據(jù)與它有完全的關(guān)系性。所以lag0上的關(guān)系性什么都不代表。而Lag1上存在著正面的關(guān)系性，即當(dāng)前數(shù)據(jù)與前一個數(shù)據(jù)的關(guān)系性是正面的。其他延遲上的自動關(guān)系性都近似于0，可以不被考慮。order=c(2,0,0)代表AR2模型。下面，我們來看生成AR2數(shù)列及分析結(jié)果。在這里，需要有兩個AR參數(shù)1和-0.9。輸入相關(guān)命令后，小窗會有如下顯示：#anAR2x=arima.sim(list(ord

33、er=c(2,0,0),ar=c(1,-.9),n=100)par(mfcol=c(3,1)plot(x,main=(expression(AR(2)phi1=1phi2=-.9)acf(x,20)pacf(x,20)接下來，我們介紹一下如何生成ARIMA（1,1,1）模型的隨機(jī)數(shù)列。其中的設(shè)置order=c(1,1,1)說明該數(shù)列是ARIMA模型，且AR的參數(shù)為0.9，MA的參數(shù)為0.5，數(shù)列長度為200。輸入相關(guān)命令后，小窗會有如下顯示：#anARIMA(1,1,1)x=arima.sim(list(order=c(1,1,1),ar=.9,ma=-.5),n=200)par(mfcol=

34、c(3,1)plot(x,main=(“expression(ARIMA(1,1,1)phi=.9theta=-.5)”)acf(x,30)#theprocessisnotstationary,sothereisnopopulationPACF.pacf(x,30)#butlookatthesamplevaluestoseehowtheydifferfromtheexamplesabove我們再介紹一下自動回歸移動平均集成模型的估計分析，即ARIMA模型。也就是說，我們手頭上有一個時間序列數(shù)據(jù)，如果知道這個數(shù)據(jù)是被ARIMA模型生成的，我們希望能估計具體的ARIMA模型的參數(shù)，可以對整個數(shù)據(jù)生

35、成的過程進(jìn)行重演和預(yù)測。下面，我們隨機(jī)生成一些ARIMA模型的數(shù)據(jù)，并且假裝我們不知道模型的參數(shù)，然后做估計練習(xí)。生成ARIMA數(shù)據(jù)的命令為arima.sim().自動回歸的參數(shù)為0.9，移動平均的參數(shù)為-0.5，估計中要假設(shè)這兩個參數(shù)的設(shè)置為未知，并進(jìn)行估計。計算估計值的命令為arima()。其中，x指隨機(jī)生成的數(shù)據(jù)。在下面的輸出中，我們可以看到AR參數(shù)的估計值為0.8465，ma的參數(shù)估計值為-0.5021，這兩個估計值與前面生成數(shù)據(jù)時設(shè)置的參數(shù)非常相似。也就是說，估計得比較精確。具體操作和顯示如下所示：x=arima.sim(list(order=c(1,0,1),ar=.9,ma=-.

36、5),n=100)#simulatesomedata(x.fit=arima(x,order=c(1,0,1)#fitthemodelandprinttheresultsCall:arima(x=x,order=c(1,0,1)Coefficients:ar1ma1intercept-NOTtheintercept-seeHYPERLINK/stoffer/tsa2/Rissues.htmlIssue1tnewRIssue10.8465-0.50210.5006s.e.0.08370.13560.3150sigma2estimatedas1.027:loglikelihood=-143.44,

37、aic=294.89下面來介紹衡量參數(shù)估計的可信度的方法。我們對估計的模型進(jìn)行檢測和評估，需要執(zhí)行的命令為tsdiag()。其中，x.fit指代估計出來的模型。這個語句輸出的結(jié)果如以下幾個圖表所示。第一個圖表代表估計模型誤差的繪圖。英文叫做StandardizedResiduals,上面有很多豎線在橫向坐標(biāo)的上下分布。如果這個估計的模型比較可信，豎線的長度是比較相似的。如果豎線的長度互相有很大出入或者根本就不同，估計模型的可信度就非常差。下面誤差繪圖中豎線的長度比較相似，都處在穩(wěn)定范圍之內(nèi)，即估計的模型沒產(chǎn)生不符合要求的誤差分布。再介紹輸出的第二張繪圖，標(biāo)題是ACFofResiduals。AC

38、F指數(shù)據(jù)點(diǎn)相互之間的關(guān)系，當(dāng)然在生成這個數(shù)據(jù)時，數(shù)據(jù)點(diǎn)之間互相獨(dú)立，并不存在任何關(guān)系。所以在這張圖上，只有位于0刻度上的豎線最高，其ACF值為1。這個0代表數(shù)據(jù)點(diǎn)與自己相比較，即數(shù)據(jù)點(diǎn)永遠(yuǎn)和它自己有關(guān)系，這種關(guān)系數(shù)值為1。其他橫向數(shù)軸上的刻度代表一個數(shù)據(jù)點(diǎn)于其他數(shù)據(jù)點(diǎn)之間的關(guān)系，這些刻度上豎線的長度幾乎等于0，即這個數(shù)據(jù)點(diǎn)與其他數(shù)據(jù)點(diǎn)沒明顯關(guān)系。這張ACF圖代表估計的模型沒造成誤差之間的任何關(guān)系。這是符合數(shù)據(jù)生成時每個數(shù)據(jù)都是獨(dú)立的這個前提的。由此可見，這ACF圖符合檢測要求。下面來介紹第三張圖，也就是Ljung-Box指標(biāo)。這個指標(biāo)可對每一個時間序列的延遲進(jìn)行顯著性的評估。這張圖的橫坐標(biāo)代表

39、時間序列的延遲，縱坐標(biāo)代表P-value，即顯著性。如果P-value十分小，就說明在其相對應(yīng)的延遲點(diǎn)上是顯著的。我們就需要拋棄所假設(shè)的模型，并且結(jié)論在所假設(shè)的模型不可信。需要注意的是，他們使用假設(shè)的模型對一個時間序列進(jìn)行估計，如果P-value是顯著的話，我們使用的模型就不可信，需要嘗試其他新模型。具體判定技巧是，P-value點(diǎn)的高度越高，我們的模型越可信。tsdiag(x.fit,gof.lag=20)#youknowtheroutine-?tsdiagfordetailsx=arima.sim(list(order=c(1,0,1),ar=.9,ma=-.5),n=100)#simul

40、atesomedata(x.fit=arima(x,order=c(1,0,1)#fitthemodelandprinttheresultstsdiag(x.fit,gof.lag=20)#youknowtheroutine-?tsdiagfordetails在對模型的參數(shù)進(jìn)行了估計后，我們接下來要對這些估計的參數(shù)進(jìn)行預(yù)測。我們需要做的是向前預(yù)測10個數(shù)據(jù)點(diǎn)。預(yù)測所使用的語句是predict().在()中輸入x.fit和n.ahead=10，其中x.fit代表前面估計出來的模型與參數(shù)，n.ahead=10代表向前預(yù)測10個數(shù)據(jù)點(diǎn)。這條語句給我們生成的預(yù)測程序項目，被命名為x.fore.在下面兩

41、條語句中，計算了兩列的新數(shù)據(jù)，分別被稱為U和L.U指最高界限，L指最低界限。最高限和最低限告訴我們未來數(shù)據(jù)有很高可能性事發(fā)生在兩個界限之間，超出或低于這兩個界限的幾率并不高。下一步使用Min和Max條語句，把數(shù)據(jù)中的最大值和最小值挑選出來，為接下來的統(tǒng)計繪圖設(shè)置圖標(biāo)界限。再使用TS.PLOT()語句.把原有的數(shù)據(jù)和向前預(yù)測10個數(shù)據(jù)點(diǎn)的數(shù)據(jù)繪在一張圖上。最后使用lines()語句，把計算的最高界限和最低界限加到繪好的圖上。這兩條界限的顏色為藍(lán)色，形式也被設(shè)置為虛線的形式。x.fore=predict(x.fit,n.ahead=10)#plottheforecastsU=x.fore$pred

42、+2*x.fore$seL=x.fore$pred-2*x.fore$seminx=min(x,L)maxx=max(x,U)ts.plot(x,x.fore$pred,col=1:2,ylim=c(minx,maxx)lines(U,col=blue,lty=dashed)lines(L,col=blue,lty=dashed)下面的繪圖是由上面的語句所完成的。圖右側(cè)除去原有的數(shù)據(jù)外，又加上了三條曲線。最上面和最下面的曲線是虛線，是計算出的最高界限和最低界限。中間的曲線為實(shí)線，是對未來的10個數(shù)據(jù)點(diǎn)進(jìn)行的預(yù)測。在計算最高、最低限時，我們使用了一個2的系數(shù)，也就是說，這兩個界限之間的幾率為95

43、%。在未來的10個數(shù)據(jù)點(diǎn)有95%的幾率是出現(xiàn)在上限和下限之間的。而出現(xiàn)在兩個界限之外的幾率小于或等于5%。4.運(yùn)用R語言研究全球變暖理論下面，來介紹R語言在全球溫度變暖研究中的應(yīng)用。我們將全球溫度變暖研究的數(shù)據(jù)叫g(shù)lobtemp2.dat.數(shù)據(jù)中有3列，第二列指全球每年溫度從1880年開始到2004年。這個數(shù)據(jù)已經(jīng)在網(wǎng)上已公布，我們只需要用下面語句讀取就行。第一條語句讀取全球溫度數(shù)據(jù)生成一個叫做u的數(shù)列矩陣；第二條語句讀取的第二列數(shù)據(jù)變成時間序列，命名為gtemp；第三條語句是將gtemp數(shù)據(jù)繪成圖來看溫度變化的具體情況。操作以后，窗口里會有如下顯示：u=read.table(/stoffer

44、/tsa2/data/globtemp2.dat)#readthedatagtemp=ts(u,2,start=1880,freq=1)#yearlytempincol2plot(gtemp)#graphtheseries(notshownhere)經(jīng)過多年，科學(xué)家們論證這個時間序列符合ARIMA(1,1,1)模型，并有0.6oC每世紀(jì)的漂移。這也就是所謂全球變暖理論。下面我們要計算出具體模型。我們有近200年的氣溫數(shù)據(jù)，且知道具體模型是ARIMA，下面我們需要估計AR和MA的具體數(shù)值。語句如下，我們從結(jié)果中可以看到AR的數(shù)值為0.2545，MA的數(shù)值為-0.7742.電腦窗口命令具體顯示如下

45、：arima(gtemp,order=c(1,1,1)Coefficients:ar1ma10.2545-0.7742s.e.0.11410.0651下面，我們就來介紹漂移參數(shù)的估計方法。具體操作如下所示：drift=1:length(gtemp)arima(gtemp,order=c(1,1,1),xreg=drift)Coefficients:ar1ma1drift0.2695-0.81800.0061s.e.0.11220.06240.0030在上面的語句中，第一行生成了一個名為drift的空數(shù)列，其長度與全球溫度數(shù)列的長度相同。而在第二個語句中，我們只要把drift數(shù)列作為xreg輸入

46、就可以了。從輸出中，我們可以看到AR參數(shù)值為0.2695，MA參數(shù)值為-0.8180，是drift值為0.0061，即漂移產(chǎn)生的值.上面提到過漂移值為+0.6，我們的估計值比提供的數(shù)值小100倍。其原因是因為我們使用的全球溫度數(shù)據(jù)是經(jīng)過被100整除的，所以估計的漂移值也相應(yīng)地小了100倍。如果將gtemp數(shù)列輸出后看一下，我們就可以明白所有溫度都是被100整除過的，就不會有什么疑惑了。具體情況，代碼顯示如下：u=read.table(/stoffer/tsa2/data/globtemp2.dat)#readthedatagtemp=ts(u,2,start=1880,freq=1)#year

47、lytempincol2plot(gtemp)drift=1:length(gtemp)arima(gtemp,order=c(1,1,1),xreg=drift)gtemp5.運(yùn)用R語言研究非獨(dú)立誤差與線性回歸在前面我們介紹的模型當(dāng)中，誤差都是相互獨(dú)立的，即誤差之間不互相影響。但如果模型中存在著互相影響的誤差，我們運(yùn)用R語言該怎么做呢？在下面例子當(dāng)中，誤差是一個自動回歸2的數(shù)列，記作AR2，意思是現(xiàn)在的誤差依賴于前兩個已經(jīng)被觀測到的誤差,誤差之間是非獨(dú)立的。具體有兩個意思，我們要學(xué)習(xí)如何從數(shù)據(jù)中確認(rèn)所有的誤差符合自動回歸模型，而且參數(shù)為2；我們要學(xué)習(xí)如何把這個AR2模型的誤差加入到線性回歸模

48、型當(dāng)中。我們使用的數(shù)據(jù)仍舊是空氣中污染物和心血管疾病死亡率的數(shù)據(jù)。在下面模型中，我們認(rèn)為，時間可以決定心血管疾病死亡人數(shù)，隨著時間向前移動，心血管疾病死亡人數(shù)會呈現(xiàn)下降趨勢。對于空氣中的污染物來講，空氣中的污染物越高，心血管疾病的死亡人數(shù)也會越高。所以在這個模型中，心血管疾病死亡人數(shù)是被預(yù)測的變量。時間變量和空氣中污染物的變量是用來進(jìn)行預(yù)測的變量。第一句命令是調(diào)用nlme，R語言程序包。其具體顯示如下：library(nlme)#loadthepackage第二句命令是將時間數(shù)列從mort數(shù)列當(dāng)中提出，用來作為進(jìn)行預(yù)測的變量稱之為trend。命令如下：trend=time(mort)#assu

49、mesmortandpartaretherefrompreviousexamples第三句命令，先建立線性回歸模型，再從模型中提取模型的誤差來研究誤差之間的相互關(guān)系性，以便確認(rèn)誤差之間是相互獨(dú)立的，還是符合前面提到的AR2模型。因為模型建立以后，我們并不知道模型的誤差是否是相互獨(dú)立的。fit.lm=lm(morttrend+part)#ols第四行命令使用上面第三行命令中生成的被稱作fit.lm線性回歸模型結(jié)果。命令Resid()是用來提取線性回歸模型中的誤差數(shù)列的。Acf（）命令是用來研究線性模型誤差的自動關(guān)系性的，也就是確定誤差之間是否相互關(guān)系。acf(resid(fit.lm)#chec

50、kacfandpacfoftheresids第五行命令同樣是將模型中誤差數(shù)列提出后使用PACF()命令進(jìn)行不完全相關(guān)性的研究。pacf(resid(fit.lm)#oruseacf2(resid(fit.lm)ifyouhaveacf2上面的pacf繪圖是第五行命令運(yùn)行的結(jié)果。在上圖中，橫向坐標(biāo)代表的延遲，被叫做lag.在lag1和lag2上出現(xiàn)了高于藍(lán)色虛線的pacf,其他lag上超出藍(lán)色虛線的距離并不明顯可以忽略。這兩個高于藍(lán)色虛線的pacf說明這個誤差數(shù)列符合AR2模型。第六行語句，使用gls()命令來對數(shù)據(jù)建立帶有AR2誤差的線性回歸模型，生成結(jié)果記為fit.gls。命令gls()在這

51、里使用時間和空氣污染物來預(yù)測心血管疾病死亡人數(shù)，然corARMA（p=2）是告訴R語言模型的誤差符合AR2模型，如果模型的誤差符合MA2模型就要改為corARMA(q=2).其中method=”ML”,代表使用極大似然方法進(jìn)行模型估計。大似然方法是一個統(tǒng)計學(xué)的核心理念，極為重要，需要使用相關(guān)的教材著重學(xué)習(xí)。具體操作如下：#residsappeartobeAR(2).nowusegls()fromnlme:fit.gls=gls(morttrend+part,correlation=corARMA(p=2),method=ML)第七行語句使用summary命令來輸出上面fit.gls模型結(jié)果的具

52、體參數(shù)值。從下面窗口顯示中，我們可以看到這個線性模型的截距為3131.5452，模型中時間對于心血管疾病死亡人數(shù)的影響，即時間的斜率值為-1.5444，也就是說隨著時間每增長一周，心血管疾病死亡率都要降低1.5444。模型中空氣污染物對心血管疾病死亡人數(shù)的影響為0.1503，也就是說空氣中污染物每上升一個單位，心血管疾病死亡人數(shù)會增加0.1503個人。下面是電腦里顯示的內(nèi)容：summary(fit.gls)Parameterestimate(s):Phi1Phi20.39805660.4134305Coefficients:ValueStd.Errort-valuep-value(Interc

53、ept)3131.5452857.21413.6531663e-04trend-1.54440.4340-3.5580214e-04part0.15030.02107.1624080e+00#residanalysis-weassumedet=1et-1+2et-2+wtwherewtisweight.事實(shí)上，ARIMA命令也可以完成上述自動關(guān)系誤差線性模型。使用ARIMA命令，其中輸入的第一個參數(shù)是mort數(shù)列，表示Mort是被預(yù)測的變量，第二個輸入的是Order=c(2,0,0)代表AR2誤差模型，xreg=cbind(trend,part)代表用時間和空氣污染物作為進(jìn)行預(yù)測的變量。下面我

54、們用ARIMA命令從新估計這個模型，所輸出的模型參數(shù)可以看出截距為3132.7085，時間的影響是-1.5449，空氣污染物的影響是0.1503。這些數(shù)值與我們上面得到的結(jié)果之間略有差異是由計算中四舍五入造成的。(fit2.gls=arima(mort,order=c(2,0,0),xreg=cbind(trend,part)Coefficients:ar1ar2intercepttrendpart0.39800.41353132.7085-1.54490.1503s.e.0.04050.0404854.66620.43280.0211sigma2estimatedas28.99:loglik

55、elihood=-1576.56,aic=3165.13Box.test(resid(fit2.gls),12,type=Ljung)#andsoon.6.運(yùn)用R語言研究估計波動的頻率我們來介紹如何運(yùn)用R語言來估計波動的頻率。波動是一種非常常見的物理現(xiàn)象。人說話發(fā)出聲波，手機(jī)通訊通過電波。我們可以用R語言對一段被記錄下來的波動進(jìn)行頻率的估計。比方說，每一年都會有流感出現(xiàn)，得流感的人數(shù)在一年中的某個時間都是不一樣的。對流感人數(shù)的記錄，就是對流感周期的記錄。我們能夠知道流感發(fā)生的頻率，有利疾病的防治。又如，股票市場其價格也是存在周期的，如果我們能計算出股票價格變化的周期，就可以預(yù)測市場的變化。我們

56、隨機(jī)生成一個波動的時間序列數(shù)據(jù)來進(jìn)行分析，闡述R語言在估計撥動的頻率。下面的語句是用來生成一些自動回歸2，即AR2波動數(shù)據(jù)的。具體操作如下：x=arima.sim(list(order=c(2,0,0),ar=c(1,-.9),n=28)#somedata第二行R語句是對設(shè)定的自動回歸模型的復(fù)雜數(shù)根進(jìn)行計算。這兩個復(fù)數(shù)根的長度都小于1，值得注意，這是自動回歸模型特有的性質(zhì)。具體操作如下所示：(u=polyroot(c(1,-1,.9)#xisAR(2)w/complexroots因為我們已經(jīng)知道數(shù)據(jù)是從自動回歸模型中生成的，根據(jù)模型條件，對于其頻率可進(jìn)行精確的計算，并不需要估計。下面語句就是對

57、這個波動頻率的計算，其模型已知為自動回顧模型。經(jīng)過四舍五入，頻率大概是0.16.Arg(u1)/(2*pi)#dominantfrequencyaround.16:下面這個語句是對繪圖輸出窗口的一些具體設(shè)置。如果不進(jìn)行這個設(shè)置，在繪圖輸出窗口中，每次只能輸出一張圖。我們需要在窗口中一次輸出3張圖，并且羅列在一起。這3張圖羅列在一起就形成了3行1列的布局，也就是下面語句中要使用(3,1)的輸入。par(mfcol=c(3,1)在估計頻率之前，需要看一看隨機(jī)生成數(shù)據(jù)的樣子。需要把這個數(shù)據(jù)繪到圖上。生成繪圖的語句如下：plot.ts(x)接下來，我們對已生成的隨機(jī)波動數(shù)據(jù)的頻率進(jìn)行估計。注意，在這里

58、不能夠使用自動回歸模型作為已知條件來進(jìn)行頻率的計算，而是要對頻率進(jìn)行估計。這里使用的R語言是SPEC.PGRAM().在下面語句中，X代表生成的隨機(jī)波動數(shù)據(jù)，SPANS=C（3,3）是對計算的側(cè)重進(jìn)行的設(shè)置。spec.pgram(x,spans=c(3,3),log=no)#nonparametricspectralestimate;alsoseespectrum()如果使用者對上面的語句有疑問，可將其介紹文件查調(diào)出來。具體使用語句如下：?spec.pgram#somehelp在生成隨機(jī)波動數(shù)據(jù)時，我們對于自動回歸模型的輸入設(shè)置為2,即AR(2)。假設(shè)我們已知模型為自動回歸，但不知道其設(shè)置為2

59、，我們也可以使用不完全的信息進(jìn)行頻率估計。但是由于不知道其設(shè)置為2，所以頻率也只能是估計出來的。具體輸入命令為：spec.ar(x,log=no)#parametricspectralestimate下面是執(zhí)行上面語句后的輸出結(jié)果。我們來看窗口中輸出的繪成圖3個圖表。最上面圖表顯示的是隨機(jī)波動的走勢。中間的圖表是在假設(shè)，用自動回歸模型的條件未知下進(jìn)行的估計。在最頂端峰相對應(yīng)的水平數(shù)軸的值就是我們所估計出來的頻率。具體來講，水平數(shù)軸上標(biāo)有刻度，從0.00到0.5.山峰最頂端所對應(yīng)的刻度在0.16左右，與前面所計算出的理論值相對應(yīng)。最下面的繪圖與中間的繪圖結(jié)論相似。峰值所對應(yīng)的水平軸刻度也是在0.16左右。但這個繪圖是在已知自