R語(yǔ)言時(shí)間序列中文教程

1、40R語(yǔ)言時(shí)間序列中文教程 2012特別聲明：R語(yǔ)言是免費(fèi)語(yǔ)言，其代碼不帶任何質(zhì)量保證，使用R語(yǔ)言所產(chǎn)生的后果由使用者負(fù)全責(zé)。前言R語(yǔ)言是一種數(shù)據(jù)分析語(yǔ)言，它是科學(xué)的免費(fèi)的數(shù)據(jù)分析語(yǔ)言，是凝聚了眾多研究人員心血的成熟的使用范圍廣泛全面的語(yǔ)言，也是學(xué)習(xí)者能較快受益的語(yǔ)言。 在R語(yǔ)言出現(xiàn)之前，數(shù)據(jù)分析的編程語(yǔ)言是SAS。當(dāng)時(shí)SAS的功能比較有限。在貝爾實(shí)驗(yàn)室里，有一群科學(xué)家討論提到，他們研究過(guò)程中需要用到數(shù)據(jù)分析軟件。SAS的局限也限制了他們的研究。于是他們想，我們貝爾實(shí)驗(yàn)室的研究歷史要比SAS長(zhǎng)好幾倍，技術(shù)力量也比SAS強(qiáng)好幾倍，且貝爾實(shí)驗(yàn)室里并不缺乏訓(xùn)練有素的專業(yè)編程人員，那么，我們貝爾實(shí)驗(yàn)室

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

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

4、語(yǔ)言樣品比較應(yīng)用舉例>>，如果學(xué)習(xí)者比較熟悉R語(yǔ)言，還可以閱讀優(yōu)秀時(shí)間序列讀物Ecomometrics in R，也可以上QuickR 網(wǎng)站。目錄R語(yǔ)言時(shí)間序列分析1前言1目錄21.運(yùn)用R語(yǔ)言研究JJ數(shù)據(jù)32.運(yùn)用R語(yǔ)言研究空氣污染193.運(yùn)用R語(yǔ)言研究自動(dòng)回歸移動(dòng)平均集成模型224.運(yùn)用R語(yǔ)言研究全球變暖理論315.運(yùn)用R語(yǔ)言研究非獨(dú)立誤差與線性回歸326.運(yùn)用R語(yǔ)言研究估計(jì)波動(dòng)的頻率367.運(yùn)用R語(yǔ)言研究厄爾尼諾頻率398.運(yùn)用R語(yǔ)言研究太陽(yáng)黑子周期頻率401.運(yùn)用R語(yǔ)言研究JJ數(shù)據(jù) 學(xué)習(xí)R言語(yǔ)時(shí)間序列分析程序操作，需要從最基礎(chǔ)、最簡(jiǎn)單的學(xué)起，例如在命令窗口中，輸入并執(zhí)行2+2

5、 等于4的R語(yǔ)言命令。2+2 1 4 執(zhí)行完2+2 等于4的R語(yǔ)言命令后，我們可以開(kāi)始時(shí)間序列，即學(xué)著把玩johnson & Johnson 數(shù)據(jù)。下載jj.dat或執(zhí)行下面語(yǔ)句。這個(gè)數(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") # read the data讀取數(shù)據(jù)jj <- scan("/stoffer

6、/tsa2/data/jj.dat") # read the data another way第二種方法讀取數(shù)據(jù)scan("/stoffer/tsa2/data/jj.dat") -> jj # and another第三種方法讀取數(shù)據(jù) 使用R語(yǔ)言的人，有的喜歡使用<-，有的喜歡使用->，大多數(shù)醫(yī)療系統(tǒng)的工作者喜歡用=，正因?yàn)槿绱瞬庞昧松厦娣N不同讀取數(shù)據(jù)的方法。 讀取數(shù)據(jù)后，鍵入并執(zhí)行jj，數(shù)據(jù)在窗口便會(huì)有如下顯示：jj 1 0.71 0.63 0.85 0.44 5 0.61 0.69 0.92 0

7、.55 . . . . . . . . . . 77 14.04 12.96 14.85 9.99 81 16.20 14.67 16.02 11.61 jj中有84個(gè)數(shù)據(jù)被稱作對(duì)象。下面命令可以顯示所有對(duì)象。objects()如果使用matlab，你會(huì)認(rèn)為jj是一個(gè)84行1列的向量，但實(shí)際上不是這樣。jj有次序，有長(zhǎng)度，但沒(méi)維度，R稱這些對(duì)象為向量，要小心區(qū)別。在R里，矩陣有維度，但向量沒(méi)維度。這都是程序語(yǔ)言的一些概念。jj1 # the first element列中第一個(gè)數(shù)據(jù) 1 0.71jj84 # the last element列中最后一個(gè)數(shù)據(jù) 1 11.61jj1:4 # the

8、 first 4 elements列中第一至第四個(gè)數(shù)據(jù) 1 0.71 0.63 0.85 0.44jj-(1:80) # everything EXCEPT the first 80 elements列中除第80個(gè)以外的所有數(shù)據(jù) 1 16.20 14.67 16.02 11.61length(jj) # the number of elements 有多少個(gè)數(shù)據(jù) 1 84dim(jj) # but no dimensions .但沒(méi)維度 NULLnrow(jj) # . no rows 沒(méi)行 NULLncol(jj) # . and no columns沒(méi)列 NULL#如果你要把jj轉(zhuǎn)變?yōu)橐粋€(gè)

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

10、令里，scan可以被read.table替代。用read.table讀取數(shù)據(jù)可生成matrix對(duì)象，還可以給每列起名字。 下面學(xué)習(xí)一下read.table, data frames, 和時(shí)間序列對(duì)象。輸入命令后，窗口會(huì)有如下顯示：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 生成的是有維

11、度的向量，read.table生成的則是帶有維度的數(shù)據(jù)架構(gòu)。 讀取jj數(shù)據(jù)的最后要領(lǐng)。如果這個(gè)數(shù)據(jù)是從1960年第三個(gè)季度開(kāi)始的，所需輸入命令則為ts(x,start=c(1960,3),frequency=4)；如果是一個(gè)每月每月的數(shù)據(jù)，例如數(shù)據(jù)是從1984年6月開(kāi)始的，需要輸入的命令則為ts(x,start=c(1984,6),frequency=12)。 輸入命令后，轉(zhuǎn)變后的時(shí)間序列對(duì)象為：jj Qtr1 Qtr2 Qtr3 Qtr4 1960 0.71 0.63 0.85 0.44 1961 0.61 0.69 0.92 0.55 . . . . . . . . . . 1979 14

12、.04 12.96 14.85 9.99 1980 16.20 14.67 16.02 11.61 注意到區(qū)別了嗎？時(shí)間信息，也就是4個(gè)不同的季度的數(shù)據(jù)被加載到里面了。 進(jìn)行時(shí)間數(shù)據(jù)分析后，窗口會(huì)有如下顯示：time(jj) Qtr1 Qtr2 Qtr3 Qtr4 1960 1960.00 1960.25 1960.50 1960.75 1961 1961.00 1961.25 1961.50 1961.75 . . . . . . . . . . . . 1979 1979.00 1979.25 1979.50 1979.75 1980 1980.00 1980.25 1980.50 198

13、0.75 接下來(lái)輸入如下組合命令。(jj = ts(scan("/stoffer/tsa2/data/jj.dat"), start=1960, frequency=4) 然后進(jìn)行對(duì)數(shù)據(jù)繪圖：plot(jj, ylab="Earnings per Share", main="J & J") 輸入以上命令后，可以看到如下結(jié)果： 再輸入下面的命令，看看區(qū)別。plot(jj, type="o", col="blue", lty="dash

14、ed")plot(diff(log(jj), main="logged and diffed") 下面利用操作plot.ts和ts.plot兩個(gè)相關(guān)命令，顯示區(qū)別。x = -5:5 # sequence of integers from -5 to 5y = 5*cos(x) # guesspar(mfrow=c(3,2) # multifigure setup: 3 rows, 2 cols#- plot:plot(x, main="plot(x)")plot(x, y, main="plot(x,y)")#- plot.

15、ts:plot.ts(x, main="plot.ts(x)")plot.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- x and y are ts objects #- the help files ? and help() are the same:?plot.tshelp(ts.plot)?par #

16、 might as well skim the graphical parameters help file while you're here從窗口中的顯示可以看出，如果數(shù)據(jù)是時(shí)間序列對(duì)象，使用plot()命令就足夠了；如果數(shù)據(jù)是平常序列，使用plot.ts()也可以做時(shí)間繪圖。 不過(guò)，把jj數(shù)據(jù)放在一張圖上，數(shù)據(jù)會(huì)隨著時(shí)間的變化上上下下跳動(dòng)，能從整體上反應(yīng)上升或者下降的趨勢(shì)。上文中用紅色光滑的曲線代表上升的趨勢(shì)，簡(jiǎn)單明了。這需要將過(guò)濾和光滑的技巧使用在jj數(shù)據(jù)上。在這里，我們用對(duì)稱的移動(dòng)平均值來(lái)達(dá)到過(guò)濾和光滑的目的。下面使用公式：fjj(t) = jj(t-2) + ¼

17、jj(t-1) + ¼ jj(t) + ¼ jj(t+1) + jj(t+2) 除此之外，lowess的過(guò)濾平滑技巧（藍(lán)色曲線）也要使用在jj數(shù)據(jù)中。具體操作如下圖：k = c(.5,1,1,1,.5) # k is the vector of weights用于對(duì)稱移動(dòng)平均的系數(shù)(k = k/sum(k) 1 0.125 0.250 0.250 0.250 0.125fjj = filter(jj, sides=2, k) # ?filter for help but you knew that already使用對(duì)稱移動(dòng)平均plot(jj)lines(fjj, col=

18、"red") # adds a line to the existing plot稱移動(dòng)平均的繪圖lines(lowess(jj), col="blue", lty="dashed")#lowess 的繪圖 操作后，窗口會(huì)顯示下面結(jié)果：看完jj數(shù)據(jù)，我們就需要開(kāi)始具體分析。第一步，我們把所有jj數(shù)據(jù)都取log值。第二步，我們把log值做差，即使用log值數(shù)列中第二值減去第一值，第三值減去第二值，第四值減去第三值等等。如果做差處理前數(shù)列里有n個(gè)數(shù)值，處理后的結(jié)果中將有n-1個(gè)數(shù)值。dljj = diff(log(jj) # differ

19、ence the logged data做log和差的處理plot(dljj) # plot it if you haven't already對(duì)結(jié)果制圖shapiro.test(dljj) # test for normality 測(cè)試結(jié)果的正態(tài)分布的性質(zhì) Shapiro-Wilk normality test data: dljj W = 0.9725, p-value = 0.07211處理完畢以上兩步，我們接下來(lái)就要將柱形分布圖和QQ圖放在一起。這兩個(gè)圖的本質(zhì)仍舊在于測(cè)試數(shù)據(jù)正態(tài)分布的性質(zhì)。數(shù)據(jù)正態(tài)分布的性質(zhì)是整個(gè)統(tǒng)計(jì)學(xué)構(gòu)架堅(jiān)實(shí)的基礎(chǔ)，如果這個(gè)性質(zhì)的存在比較可信、通過(guò)了所有測(cè)試

20、，統(tǒng)計(jì)分析中得出的結(jié)論就比較可信、就通得過(guò)所有測(cè)試。當(dāng)然如果這個(gè)性質(zhì)在數(shù)據(jù)中不存在，我們需要用其它的技巧來(lái)處理。詳細(xì)的，參看R語(yǔ)言樣品比較應(yīng)用的實(shí)例。以上操作，在窗口中有如下顯示：par(mfrow=c(2,1) # set up the graphics 設(shè)置為兩圖的輸出hist(dljj, prob=TRUE, 12) # histogram柱形分布圖 lines(density(dljj) # smooth it - ?density for details柱形分布圖的曲線 qqnorm(dljj) # normal Q-Q plot QQ圖 qqline(dljj) # add a l

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

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

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

24、關(guān)系。其他幾個(gè)延遲存在著負(fù)面關(guān)系。 我們可以利用這4和8的正面關(guān)系來(lái)進(jìn)行預(yù)測(cè)，即用現(xiàn)有數(shù)據(jù)計(jì)算接下來(lái)的第4個(gè)或者是第8個(gè)數(shù)據(jù)。 下面我們來(lái)看ACF和PACF圖表。 確定我們已經(jīng)觀察到的正面和負(fù)面關(guān)系。par(mfrow=c(2,1) # The power of accurate observation is commonly called cynicism # by those who have not got it. - George Bernard Shawacf(dljj, 20) # ACF to lag 20 - no graph shown. keep readingpacf(d

25、ljj, 20) # PACF to lag 20 - no graph shown. keep reading# !NOTE! acf2 on the line below is NOT available in R. details follow the graph belowacf2(dljj) # this is what you'll see below在上圖中，ACF和PACF橫坐標(biāo)的標(biāo)記是1,2,3,4,5. 但因?yàn)閿?shù)據(jù)是季度性的，每年有4個(gè)季度所以1,2,3,4,5的標(biāo)記代表的4,8,12,16,20的延遲。當(dāng)然，如果我們不喜歡上面橫坐標(biāo)的標(biāo)記，也可以將dljj更改為t

26、s(dljj, freq=1); 也就是說(shuō) acf(ts(dljj, freq=1), 20)。因?yàn)樵谏厦鍭CF圖表中橫坐標(biāo)1代表的是延遲4，橫坐標(biāo)2代表的是延遲6，其相應(yīng)的豎線代表的就是延遲4和8的正面關(guān)系。 接下來(lái)，下面我們介紹結(jié)構(gòu)拆析。在前面R代碼中，我們?cè)鴮⑺衘j數(shù)據(jù)進(jìn)行了lag變型。在變型后的數(shù)據(jù)中，存在著上升趨勢(shì)，季節(jié)的影響和每一時(shí)間點(diǎn)產(chǎn)生的隨機(jī)的誤差。根據(jù)這一數(shù)據(jù)圖，我們能夠把趨勢(shì)、季節(jié)和誤差從變型后的jj數(shù)據(jù)中拆析出來(lái)，分別研究，或者分別進(jìn)行繪圖，以便于單獨(dú)檢查。下面等式代表將要使用的數(shù)學(xué)模型：Log(jj)=趨勢(shì)+季節(jié)+誤差 log(jj) = trend + season

27、 + error結(jié)構(gòu)拆析的R命令是stl(), 下面語(yǔ)句中stl命令中輸入的是lag變型后的jj數(shù)據(jù)。其中的“per”輸入指的是使用季節(jié)循環(huán)來(lái)進(jìn)行拆析。stl語(yǔ)句在這里生成了一個(gè)叫dog的R物件，然后Plot語(yǔ)句將dog物件進(jìn)行繪圖。具體操作如下圖” plot(dog <- stl(log(jj), "per") 窗口會(huì)出現(xiàn)下面所示：上圖中有四行R的繪圖，其中第一行代表原來(lái)的log（jj）的數(shù)據(jù)。此數(shù)據(jù)可以看到總體的上升趨勢(shì)還存在著一定季節(jié)循環(huán)性的變化。第二行繪圖代表拆析后季節(jié)循環(huán)的作用。第三行繪圖代表拆析后將季節(jié)循環(huán)作用清除剩余的上升趨勢(shì)，此數(shù)據(jù)清楚地看到那種循環(huán)性

28、變化已經(jīng)不存在，剩余的只是趨勢(shì)。第四行繪圖代表將季節(jié)循環(huán)作用和總體的趨勢(shì)從數(shù)據(jù)中清除后所剩余的隨機(jī)產(chǎn)生的誤差。 如果我們需要對(duì)數(shù)據(jù)的誤差進(jìn)行一些常規(guī)檢測(cè)，例如進(jìn)行正態(tài)分布檢測(cè)，繪制QQ圖，還有繪制柱形圖。我們所需要的具體誤差數(shù)據(jù)被存在叫做dog$time.series,3的數(shù)列里。即叫dog的物件中有個(gè)叫time.series的數(shù)據(jù)矩陣，誤差就被存儲(chǔ)在這個(gè)數(shù)據(jù)矩陣的第三列里。$指調(diào)取dog物件中的time.series數(shù)據(jù)矩陣。,3指數(shù)據(jù)矩陣中第三列。如果要對(duì)這一數(shù)列的誤差值進(jìn)行ACF的分析，只需要執(zhí)行命令acf(dog$time.series,3)。 再接下來(lái)，我們對(duì)log變型后的jj數(shù)據(jù)進(jìn)

29、行線性回歸模型分析。與上面結(jié)構(gòu)拆析不同的是，我們?cè)谶@里使用四個(gè)季度來(lái)量化季節(jié)循環(huán)對(duì)數(shù)據(jù)的影響。一年中有四個(gè)季度，也是我們所使用數(shù)據(jù)所代表的。這個(gè)jj數(shù)據(jù)是某一家公司的季度收入數(shù)據(jù)，從上面繪圖中我們就可以看到，每一年第三季度就會(huì)出現(xiàn)一個(gè)收入高峰，隨之而來(lái)第四季度收入就會(huì)跌入低谷。然后在一季度和二季度收入又會(huì)逐漸上升。這也就是說(shuō)，每一季度對(duì)這家公司收入的影響都是不一樣的。 具體考慮到這種季度之間的不同，我們可使用如下數(shù)學(xué)模型： log(jj)= *time + 1*Q1 + 2*Q2 + 3*Q3 + 4*Q4 + 這個(gè)數(shù)學(xué)模型的意思是： log(jj)=趨勢(shì)*時(shí)間+一季度的影響*一季度+二季度的

30、影響*二季度+三季度的影響*三季度+四季度的影響*四季度+誤差 上面的模型代表的就是總體上升趨勢(shì)，1 2 3 4代表的是四個(gè)季度的影響。有一個(gè)非常有趣的問(wèn)題，上面模型是把所有四個(gè)季度的趨勢(shì)都加在了一起，其結(jié)果卻是某單一季度的收入。四個(gè)季度的和如何能夠與一個(gè)季度相等問(wèn)題就出在Q1 Q2 Q3 Q4 上。因?yàn)镼被我們稱作指示性函數(shù)。函數(shù)的意思就是數(shù)據(jù)進(jìn)，數(shù)據(jù)出，也就是說(shuō)把一個(gè)數(shù)據(jù)輸入到一個(gè)函數(shù)中，那個(gè)函數(shù)就會(huì)輸出一個(gè)結(jié)果。 以上面的Q1函數(shù)為例，Q1只能輸出兩種結(jié)果，1 和0. Q1所需要的輸入是四種1，2，3，4， 代表四個(gè)季度。把1輸入到Q1函數(shù)中時(shí)，Q1函數(shù)輸出的結(jié)果為1，當(dāng)把2,3,4輸入

31、Q1函數(shù)時(shí)，Q1函數(shù)輸出的結(jié)果為0. 與Q1函數(shù)類似，Q2函數(shù)的輸入也是1,2,3,4， 但只有輸入為2時(shí)，Q2函數(shù)的輸出才為1，當(dāng)輸入為1,3,4 時(shí)，Q2函數(shù)的輸出為0. Q3函數(shù)輸入為1,2,3,4，只有當(dāng)輸入為3時(shí)，輸出為1，輸入其他數(shù)據(jù)時(shí)，輸出為0.Q4函數(shù)的輸入為1,2,3,4，只有當(dāng)輸入為4時(shí)，輸出為1，其他數(shù)據(jù)時(shí)，輸出為0.我們?cè)倩氐缴厦娴哪Ｐ?，?dāng)一個(gè)數(shù)據(jù)是從第一季度中記錄下來(lái)的，Q1給出數(shù)值1，Q2給出數(shù)值0，Q3給出數(shù)值0，Q4給出的數(shù)值0。因?yàn)檫@時(shí)Q2，Q3，Q4都是0，二季度，三季度，四季度的影響被0相乘后也變成了0. 所以在第一季度Q1為1，而其他的為0.我們就只考慮

32、了一季度的影響，其他季度的影響不存在。同理，當(dāng)季度為二、三、四時(shí)也有類似結(jié)果。下面是建立這個(gè)線性模型的R語(yǔ)句，只有頭三行是用來(lái)生成線性模型的，第四條語(yǔ)句summary()用來(lái)輸出模型參數(shù)數(shù)值。具體操作以及顯示如下： Q = factor(rep(1:4,21) # make (Q)uarter factors that's repeat 1,2,3,4, 21 timestrend = time(jj)-1970 # not necessary to "center" time, but the results look nicerreg = lm(log(jj)0

33、+trend+Q, na.action=NULL) # run the regression without an intercept#- the na.action statement is to retain time series attributessummary(reg) Call:lm(formula = log(jj) 0 + trend + Q, na.action = NULL)Residuals: Min 1Q Median 3Q Max -0.29318 -0.09062 -0.01180 0.08460 0.27644 Coefficients: Estimate St

34、d. Error t value Pr(>|t|) trend 0.167172 0.002259 74.00 <2e-16 *Q1 1.052793 0.027359 38.48 <2e-16 *Q2 1.080916 0.027365 39.50 <2e-16 *Q3 1.151024 0.027383 42.03 <2e-16 *Q4 0.882266 0.027412 32.19 <2e-16 *-Signif. codes: 0 '*' 0.001 '*' 0.01 '*' 0.05 '.&#

35、39; 0.1 ' ' 1 Residual standard error: 0.1254 on 79 degrees of freedomMultiple R-squared: 0.9935, Adjusted R-squared: 0.9931 F-statistic: 2407 on 5 and 79 DF, p-value: < 2.2e-16 上面語(yǔ)句第一行用來(lái)生成所需要的Q函數(shù)。第二行用來(lái)生成線性回歸模型中的時(shí)間輸入，然后存儲(chǔ)于叫做trend的數(shù)列中。第三行語(yǔ)句是用lm()命令來(lái)建立線性模型。在R語(yǔ)言中，用來(lái)分隔進(jìn)行預(yù)測(cè)的變量和被預(yù)測(cè)的變量。在上面語(yǔ)句第三行中

36、，log（jj）是被預(yù)測(cè)的變量，它被放置在的左邊。而且0，trend, 和Q函數(shù)都是用來(lái)進(jìn)行預(yù)測(cè)的變量。它們被放在右邊。其目的就是要用0TREND和Q function來(lái)預(yù)測(cè)Log（jj）。現(xiàn)在，我們復(fù)習(xí)一下什么叫做線性。線性是在坐標(biāo)平面上的一條直線。坐標(biāo)平面上有兩個(gè)數(shù)軸，橫軸叫做x軸，豎軸叫做y軸。在平面上畫(huà)一條直，可使用截距和斜率定義它。如果這條直線通過(guò)原點(diǎn)，其截距為0. 在右邊的0是我們模型中的截距，也就是說(shuō)我們模型中的直線通過(guò)原點(diǎn)。如果我們決定使用有截距的線性模型，即模型中的直線不通過(guò)原點(diǎn)，就要把第三條R語(yǔ)句中的0改為1.這樣，R語(yǔ)言就可以將截距計(jì)算出來(lái)。 執(zhí)行summary（）命令后

37、，我們就可得到具體模型參數(shù)的數(shù)值。模型中的趨勢(shì)，即估計(jì)值為0.167172，一季度的影響記為a1，估計(jì)值為1.052793，第二季度的影響記為a2，估計(jì)值為1.080916，第三季度的影響記為a3，估計(jì)值為 1.151024，第四季度的影響記為a4，估計(jì)值為 0.882266。 隨后，我們將所有估計(jì)值插入上面模型中，可得到下面等式： Log(jj)=0.167172*時(shí)間+1.052793*Q1+ 1.080916*Q2+ 1.151024 *Q3+0.882266*Q4+誤差 然后，我們就可以進(jìn)行預(yù)測(cè)。比如說(shuō)，在第200個(gè)時(shí)間所處的季度為4，其預(yù)測(cè)值就可被計(jì)算為： Log(jj)= 0.16

38、7172*200+1.052793*0+ 1.080916*0+ 1.151024 *0 +0.882266*1=Now check out what happened. Look at a plot of the observations and their fitted values: 在上面，我們已經(jīng)知道如何使用計(jì)算出來(lái)的模型對(duì)某一時(shí)間點(diǎn)進(jìn)行預(yù)測(cè)。但是，我們的預(yù)測(cè)是否符合在實(shí)際中觀察到的時(shí)間序列數(shù)據(jù)，我們需要將實(shí)際當(dāng)中觀察到的數(shù)據(jù)與模型計(jì)算出來(lái)的數(shù)據(jù)放在一起，繪成圖表進(jìn)行比較。 接下來(lái)要進(jìn)行的工作是，將預(yù)測(cè)的數(shù)據(jù)和實(shí)觀察的數(shù)據(jù)進(jìn)行比較。下面第一條語(yǔ)句是對(duì)實(shí)際觀察到的log(jj)數(shù)據(jù)進(jìn)行

39、繪圖。第二行語(yǔ)句是對(duì)計(jì)算出的數(shù)據(jù)進(jìn)行繪圖，其中fitted(reg)是使用我們已建立好的模型來(lái)計(jì)算預(yù)期值。 plot(log(jj), type="o") # the data in black with little dots lines(fitted(reg), col=2) # the fitted values in bloody red - or use lines(reg$fitted, col=2) 輸入上面兩個(gè)命令后，窗口就會(huì)顯示如下圖表。在圖表中，黑線代表實(shí)際中觀察到的log(jj)數(shù)據(jù)。紅線代表由模型計(jì)算出來(lái)的log(jj)數(shù)據(jù)。它們看上去非常相似，但區(qū)

40、別在于模型中計(jì)算出來(lái)的數(shù)值沒(méi)誤差。. and a plot of the residuals and the ACF of the residuals:在對(duì)log(jj)建立模型時(shí)，具體的誤差也被估計(jì)了出來(lái)，并存儲(chǔ)于一個(gè)叫做residuals的數(shù)列中，可使用命令resid(reg)進(jìn)行調(diào)用。下面，plot(resid(reg)語(yǔ)句是用誤差的繪圖來(lái)確定誤差的變化是否比較小范圍之內(nèi)。acf(resid(reg),20)語(yǔ)句是用來(lái)檢驗(yàn)誤差的自動(dòng)關(guān)聯(lián)性。因?yàn)槲覀冎勒`差之間不應(yīng)當(dāng)存在任何關(guān)聯(lián)性。如果誤差之間沒(méi)關(guān)聯(lián)性存在，只有在lag0的ACF上有數(shù)值為1的關(guān)聯(lián)性。具體操作和顯示如下： par(mfro

41、w=c(2,1) plot(resid(reg) # residuals - reg$resid is same as resid(reg) acf(resid(reg),20) # acf of the resids 下面給出了相關(guān)誤差的兩個(gè)繪圖。第一個(gè)圖可以看到誤差的變化范圍在允許范圍之內(nèi)，并且沒(méi)大變化。第二張圖只有在lag0上有數(shù)值為1的關(guān)聯(lián)性，其他lag上的關(guān)聯(lián)性都非常小，即說(shuō)明誤差和誤差之間是沒(méi)關(guān)系。這符合我們模型所需要假設(shè)的。2.運(yùn)用R語(yǔ)言研究空氣污染接下來(lái)，我們講解一個(gè)延遲線性回歸模型，同樣要使用到R語(yǔ)言中l(wèi)m()命令來(lái)生成模型。我們以空氣污染研究為實(shí)例進(jìn)行講解。在實(shí)例中，我們要

42、使用空氣中的固體污染物含量來(lái)預(yù)測(cè)心血管疾病死亡人數(shù)，并且延遲四周，即一個(gè)月。心血管疾病死亡人數(shù)的數(shù)據(jù)記為cmort.dat, 空氣中固體污染物的數(shù)據(jù)記為part.dat。先將相關(guān)數(shù)據(jù)傳到網(wǎng)上，然后用R語(yǔ)言掃描其連接，進(jìn)而讀取數(shù)據(jù)。下面是R語(yǔ)言讀取兩個(gè)數(shù)據(jù)所使用的命令語(yǔ)言。 mort = ts(scan("/stoffer/tsa2/data/cmort.dat"),start=1970, frequency=52) # make these time series objects Read 508 items part = ts

43、(scan("/stoffer/tsa2/data/part.dat"),start=1970, frequency=52) Read 508 items 讀取數(shù)據(jù)后，運(yùn)用ersect()命令把mort 數(shù)據(jù)、part數(shù)據(jù)以及part經(jīng)過(guò)延遲四周的數(shù)據(jù)捆綁在了一起，使三列數(shù)據(jù)形成行與行之間的正確對(duì)應(yīng)關(guān)系，然后進(jìn)行接線性回歸分析。在這個(gè)線性模型中，我們要使用當(dāng)前空氣污染物的數(shù)據(jù)和四周以前空氣污染物數(shù)據(jù)對(duì)當(dāng)前心血管病死亡的人數(shù)進(jìn)行預(yù)測(cè)。因?yàn)槲覀兪褂盟膫€(gè)星期以前的數(shù)據(jù)，因而所建立起來(lái)的模型被稱為延遲模型，由R語(yǔ)言中的lag

44、（）命令來(lái)實(shí)現(xiàn)，延遲四個(gè)星期則用-4表示。 ded = ersect(mort,part,part4=lag(part,-4), dframe=TRUE) # tie them together in a data frame 接下來(lái)的語(yǔ)句表示的是如何生成線性模型，使用的也是lm（）命令。在語(yǔ)句中，mort數(shù)列表示被預(yù)測(cè)的變量，放在左邊，part變量和延遲了四周的part變量屬于用來(lái)進(jìn)行預(yù)測(cè)的變量，要放在右邊。上面已被捆綁好的數(shù)據(jù)被稱作ded，要跟在data = 后面。 fit = lm(mortpart+part4, data=ded, na.action=NULL) # now

45、 the regression will work 接下來(lái)使用summary命令輸出所有估計(jì)完成的參數(shù)。summary(fit) Call: lm(formula = mort part + part4, data = ded, na.action = NULL) Residuals: Min 1Q Median 3Q Max -22.7429 -5.3677 -0.4136 5.2694 37.8539 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 69.01020 1.37498 50.190 <

46、2e-16 * part 0.15140 0.02898 5.225 2.56e-07 * part4 0.26297 0.02899 9.071 < 2e-16 * - Signif. codes: 0 '*' 0.001 '*' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 8.323 on 501 degrees of freedom Multiple R-Squared: 0.3091, Adjusted R-squared: 0.3063

47、F-statistic: 112.1 on 2 and 501 DF, p-value: < 2.2e-16 如上圖所示，我們便獲得了所研究需要的數(shù)據(jù)。3.運(yùn)用R語(yǔ)言研究自動(dòng)回歸移動(dòng)平均集成模型 現(xiàn)在來(lái)介紹自動(dòng)回歸移動(dòng)平均集成模型，即ARIMA模型。AR指自動(dòng)回歸，MA指移動(dòng)平均，I指集成。 我們先來(lái)生成兩個(gè)隨機(jī)的自動(dòng)回歸數(shù)列，記為x1和x2。它們都是ar1數(shù)列，因?yàn)樯蛇^(guò)程中使用order=c(1,0,0)的參數(shù)輸入。X1的AR參數(shù)是+0.9，x2的AR參數(shù)是-0.9。 兩個(gè)語(yǔ)句生成完畢后，我們以將這兩個(gè)數(shù)列繪圖出來(lái)。兩個(gè)時(shí)間序列非常不同，x1的參數(shù)為+0.9，它的圖像類似于股票價(jià)格的

48、變動(dòng)圖形；x2的參數(shù)為-0.9， 它的圖像看上去像是一列聲波。窗口具體顯示如下：x1 = arima.sim(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語(yǔ)言中，ACF指自動(dòng)關(guān)系性，ACF即Auto-Correlation Function的簡(jiǎn)稱. 比方說(shuō)，股票價(jià)格今天的價(jià)格跟昨天的價(jià)格有關(guān)系，明天的價(jià)格會(huì)跟今天的或者昨天的價(jià)格有關(guān)系。它們之間的關(guān)系性便用ACF來(lái)衡量。PACF被稱作不完全自動(dòng)

49、關(guān)系性。自動(dòng)關(guān)系性ACF中存在著線性關(guān)系性和非線性關(guān)系性。不完全自動(dòng)關(guān)系性就是把線性關(guān)系性從自動(dòng)關(guān)系性中消除。如果在線性關(guān)系性被去除以后，兩個(gè)時(shí)間點(diǎn)之間的關(guān)系性也就是不完全關(guān)系性。當(dāng)PACF近似于0，這表明兩個(gè)時(shí)間點(diǎn)之間的關(guān)系性是完全由線性關(guān)系性所造成的。如果不完全關(guān)系性在兩個(gè)時(shí)間點(diǎn)之間不近似于0，這表明線性模型是不能夠表達(dá)這兩個(gè)時(shí)間點(diǎn)之間的關(guān)系。我們用下面的語(yǔ)句來(lái)檢驗(yàn)x1和x2數(shù)列中的ACF和PACF。par(mfcol=c(2,2)acf(x1, 20)acf(x2, 20)pacf(x1, 20)pacf(x2, 20)對(duì)數(shù)列x1來(lái)講，ACF從lag1到lag12都高于藍(lán)色虛線，也就是說(shuō)

50、兩個(gè)時(shí)間點(diǎn)的距離在1到12之間它們的自動(dòng)關(guān)系都是正面的。所有線性關(guān)系在x1數(shù)列中被清除了結(jié)果是Partial ACF的x1圖表。在PACF x1圖表上可以看到只有l(wèi)ag0值為1，其他的lag上的關(guān)系值都低于藍(lán)色虛線，近似于0，也就是說(shuō)在x1數(shù)列中存在的自動(dòng)關(guān)系基本上是線性關(guān)系。使用線性模型可以符合x(chóng)1數(shù)列的要求。對(duì)數(shù)列x2來(lái)講，兩個(gè)時(shí)間點(diǎn)之間的ACF關(guān)系有負(fù)面的也有正面的，但經(jīng)過(guò)去除線性關(guān)系以后所有不完全自動(dòng)關(guān)系都接近于0，也就是說(shuō)x2數(shù)列中的自動(dòng)關(guān)系也基本上是線性的。使用線性模型符合x(chóng)2數(shù)列要求。下面，我們來(lái)介紹移動(dòng)平均MA模型。我們先用ARIMA.SIM()命令生成一個(gè)移動(dòng)平均的數(shù)據(jù)列。不

51、同的是，order=c(0,0,1)，ma1代表生成的數(shù)列，ma是一個(gè)參數(shù)，設(shè)置為+0.8，數(shù)列中有100個(gè)數(shù)據(jù)。其余3行命令是用來(lái)設(shè)置輸出窗口繪制圖表和計(jì)算ACF，PACF數(shù)值的。 經(jīng)過(guò)一系列操作后，代碼及其結(jié)果如下所示：# an MA1 x = 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ī)生成的，我們?cè)谥匦逻\(yùn)行這些命令時(shí)，生成的數(shù)據(jù)會(huì)與以前所生成的數(shù)據(jù)完全不同，其結(jié)果也不

52、一樣。除非在以上代碼的前面使用隨機(jī)鎖定的命令set.seed(1)，每次生成的就是相同的數(shù)據(jù)了。當(dāng)然也可以用set.seed(100)，和其它數(shù)字輸入都是可以的。在上圖中，MA1代表了數(shù)列輸出的結(jié)果。最上面顯示的是數(shù)列的數(shù)目，橫向坐標(biāo)從1到100代表有100個(gè)數(shù)據(jù)。中間的圖顯示的是x數(shù)列的自動(dòng)關(guān)系性，在lag0上有很強(qiáng)的自動(dòng)關(guān)系性，因?yàn)長(zhǎng)ag0表示沒(méi)時(shí)間延遲，當(dāng)前數(shù)據(jù)與它有完全的關(guān)系性。所以lag0上的關(guān)系性什么都不代表。而Lag1上存在著正面的關(guān)系性，即當(dāng)前數(shù)據(jù)與前一個(gè)數(shù)據(jù)的關(guān)系性是正面的。其他延遲上的自動(dòng)關(guān)系性都近似于0，可以不被考慮。 order=c(2,0,0)代表AR2模型。下面，我

53、們來(lái)看生成AR2數(shù)列及分析結(jié)果。在這里，需要有兩個(gè)AR參數(shù)1和-0.9。輸入相關(guān)命令后，小窗會(huì)有如下顯示：# an AR2 x = arima.sim(list(order=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) 接下來(lái)，我們介紹一下如何生成ARIMA（1,1,1）模型的隨機(jī)數(shù)列。其中的設(shè)置order=c(1,1,1)說(shuō)明該數(shù)列是ARIMA模型，且AR的參數(shù)為0.9，MA的參數(shù)為0.5，數(shù)列長(zhǎng)度為200。輸入相

54、關(guān)命令后，小窗會(huì)有如下顯示：# an ARIMA(1,1,1) x = arima.sim(list(order=c(1,1,1), ar=.9, ma=-.5), n=200)par(mfcol=c(3,1)plot(x, main=(“expression(ARIMA(1,1,1)phi=.9theta=-.5)”)acf(x, 30) # the process is not stationary, so there is no population PACF . pacf(x, 30) # but look at the sample values to see how they di

55、ffer from the examples above 我們?cè)俳榻B一下自動(dòng)回歸移動(dòng)平均集成模型的估計(jì)分析，即ARIMA模型。也就是說(shuō)，我們手頭上有一個(gè)時(shí)間序列數(shù)據(jù)，如果知道這個(gè)數(shù)據(jù)是被ARIMA模型生成的，我們希望能估計(jì)具體的ARIMA模型的參數(shù)，可以對(duì)整個(gè)數(shù)據(jù)生成的過(guò)程進(jìn)行重演和預(yù)測(cè)。下面，我們隨機(jī)生成一些ARIMA模型的數(shù)據(jù)，并且假裝我們不知道模型的參數(shù)，然后做估計(jì)練習(xí)。 生成ARIMA數(shù)據(jù)的命令為arima.sim(). 自動(dòng)回歸的參數(shù)為0.9，移動(dòng)平均的參數(shù)為-0.5，估計(jì)中要假設(shè)這兩個(gè)參數(shù)的設(shè)置為未知，并進(jìn)行估計(jì)。計(jì)算估計(jì)值的命令為arima()。其中，x指隨機(jī)生成的數(shù)據(jù)。在下面的輸出中，我們可以看到AR參數(shù)的估計(jì)值為0.8465， ma的參數(shù)估計(jì)值為-0.5021，這兩個(gè)估計(jì)值與前面生成數(shù)據(jù)時(shí)設(shè)置的參數(shù)非常相似。也就是說(shuō)，估計(jì)得比較精確。具體操作和顯示如下所示：x = arima.sim(lis