R語言解讀自回歸模型

時(shí)間序列是金融分析中常用到的一種數(shù)據(jù)格式，自回歸模型是分析時(shí)間序列數(shù)據(jù)的一種基本的方法。通過建立自回歸模型，找到數(shù)據(jù)自身周期性的規(guī)律，從而幫助我們理解金融市場的發(fā)展變化。

在時(shí)間序列分析中，有一個(gè)常用的模型包括AR，MA，ARMA，ARIMA，ARCH，GARCH，他們的主要區(qū)別是適用條件不同，且是層層遞進(jìn)的，后面的一個(gè)模型解決了前一個(gè)模型的某個(gè)固有問題。本文以AR模型做為開始，將對時(shí)間序列分析體系，進(jìn)行完整的介紹，并用R語言進(jìn)行模型實(shí)現(xiàn)。

由于本文為非統(tǒng)計(jì)的專業(yè)文章，所以當(dāng)出現(xiàn)與教課書不符的描述，請以教課書為準(zhǔn)。本文力求用簡化的語言，來介紹自回歸模型的知識，同時(shí)配合R語言的實(shí)現(xiàn)。

目錄

自回歸模型介紹

用R語言構(gòu)建自回歸模型

模型識別ACF/PACF

模型預(yù)測1. 自回歸模型(AR)

自回歸模型(Autoregressive model)，簡稱AR模型，是統(tǒng)計(jì)上一種處理時(shí)間序列的方法，用來描述當(dāng)前值與歷史值之間的關(guān)系，用變量自身的歷史時(shí)間數(shù)據(jù)對自身進(jìn)行預(yù)測，自回歸模型必須滿足平穩(wěn)性的要求。比如，時(shí)間序列數(shù)據(jù)集X 的歷史各期數(shù)據(jù)從X1至Xt-1，假設(shè)它們?yōu)榫€性關(guān)系，可以對當(dāng)期Xt的表現(xiàn)進(jìn)行預(yù)測。X的當(dāng)期值等于一個(gè)或數(shù)個(gè)落后期的線性組合，加常數(shù)項(xiàng)，加隨機(jī)誤差。

p階自回歸過程的公式定義：

字段解釋：

Xt是當(dāng)期的X的表現(xiàn)

c是常數(shù)項(xiàng)

p是階數(shù)，i為從1到p的值

φi是自相關(guān)系數(shù)

t為時(shí)間周期

εt是均值為0，標(biāo)準(zhǔn)差為δ 的隨機(jī)誤差，同時(shí)δ是獨(dú)立于t的

對于一階自回模型，用AR(1)來表示，簡化后的公式為：

自回歸是從線性回歸分析中發(fā)展而來，只是把自變量x對因變量y的分析，變成自變量x對自身的分析。如果你需要了解線性回歸的知識，請參考文章R語言解讀一元線性回歸模型。

自回歸模型的限制

自回歸模型是用自身的數(shù)據(jù)來進(jìn)行預(yù)測，但是這種方法受到一定的限制：

必須具有平穩(wěn)性，平穩(wěn)性要求隨機(jī)過程的隨機(jī)特征不隨時(shí)間變化。

必須具有自相關(guān)性，如果自相關(guān)系數(shù)(φi)小于0.5，則不宜采用，否則預(yù)測結(jié)果極不準(zhǔn)確。

自回歸只適用于預(yù)測與自身前期相關(guān)的現(xiàn)象，即受自身歷史因素影響較大的現(xiàn)象。對于受其他因素影響的現(xiàn)象，不宜采用自回歸，可以改用向量自回歸模型。

平穩(wěn)性時(shí)間序列的特點(diǎn)

平穩(wěn)性要求產(chǎn)生時(shí)間序列Y的隨機(jī)過程的隨機(jī)特征不隨時(shí)間變化，則稱過程是平穩(wěn)的；假如該隨機(jī)過程的隨機(jī)特征隨時(shí)間變化，則稱過程是非平穩(wěn)的。

平穩(wěn)性是由樣本時(shí)間序列所得到的擬合曲線，在未來的一段期間內(nèi)能順著現(xiàn)有的形態(tài)能一直地延續(xù)下去；如果數(shù)據(jù)非平穩(wěn)，則說明樣本擬合曲線的形態(tài)不具有延續(xù)的特點(diǎn)，也就是說擬合出來的曲線將不符合當(dāng)前曲線的形態(tài)。

隨機(jī)變量Yt的均值和方差均與時(shí)間t無關(guān)

隨機(jī)變量Yt和Ys的協(xié)方差只與時(shí)間差（步長）t-s有關(guān)

對于平穩(wěn)時(shí)間序列在數(shù)學(xué)上有比較豐富的處理手段，非平穩(wěn)的時(shí)間序列通過差分等手段轉(zhuǎn)化為平穩(wěn)時(shí)間序列處理

2. 用R語言構(gòu)建自回歸模型

了解了自回歸模型的定義，我們就可以用R語言來模擬一下自回歸模型的構(gòu)建和計(jì)算過程。

生成一個(gè)隨機(jī)游走的數(shù)據(jù)集，滿足平穩(wěn)性的要求。

# 隨機(jī)游走的數(shù)據(jù)集> set.seed(0)> x<-w<-rnorm(1000) # 生成符合正態(tài)分布N(0,1)的數(shù)據(jù)> for(t in 2:1000) x[t]<-x[t-1]+w[t]> tsx<-ts(x) # 生成ts時(shí)間序列的數(shù)據(jù)集# 查看數(shù)據(jù)集> head(tsx,15) [1] 1.2629543 0.9367209 2.2665202 3.5389495 3.9535909 2.4136409 [7] 1.4850739 1.1903534 1.1845862 3.5892396 4.3528331 3.5538238[13] 2.4061668 2.1167053 1.8174901> plot(tsx) # 生成可視化圖形 > a<-ar(tsx);a # 進(jìn)行自回歸建模Call:ar(x = tsx)Coefficients: 1 0.9879 Order selected 1 sigma^2 estimated as 1.168

數(shù)據(jù)的如圖展示：

自相關(guān)系數(shù)為0.9879 ，這是一個(gè)非常強(qiáng)的自相關(guān)性，所以上述的數(shù)列符合自相關(guān)的特性。

R語言中ar()函數(shù)提供了多種自相關(guān)系數(shù)的估計(jì)，包括"yule-walker", "burg", "ols", "mle", "yw"，默認(rèn)是用yule-walker方法，常用的方法還有最小二乘法(ols)，極大似然法(mle)。

我們用最小二乘法，來進(jìn)行參數(shù)估計(jì)。

> b<-ar(tsx,method = "ols");bCall:ar(x = tsx, method = "ols")Coefficients: 1 0.9911 Intercept: -0.017 (0.03149) Order selected 1 sigma^2 estimated as 0.9906

用最小二乘法的計(jì)算結(jié)果，則自相關(guān)系統(tǒng)數(shù)為0.9911，截距為-0.017。只有使用最小二乘法進(jìn)行參數(shù)估計(jì)的時(shí)候，才會有截距。

我們用極大似然法，來進(jìn)行參數(shù)估計(jì)。

> d<-ar(tsx,method = "mle");dCall:ar(x = tsx, method = "mle")Coefficients: 1 0.9904 Order selected 1 sigma^2 estimated as 0.9902

用極大似然法計(jì)算結(jié)果，則自相關(guān)系統(tǒng)數(shù)為0.9904。對于上面3種估計(jì)方法，自相關(guān)系數(shù)的值都是很接近的。

3. 模型識別ACF/PACF

在上面的例子中，我們默認(rèn)是用一階的自回歸模型AR(1)，進(jìn)行程序?qū)崿F(xiàn)的。在實(shí)際應(yīng)用中，自回歸模型AR時(shí)間序列的階數(shù)P是未知的，必須根據(jù)實(shí)際數(shù)據(jù)來決定，就要對AR模型定階數(shù)。常的方法就是利用自相關(guān)函數(shù)(ACF)和偏自相關(guān)函數(shù)(PACF)來確定自回歸模型的階數(shù)。在ACF/PACF不能確定的情況下，還需要用AIC(Aikaike info Criterion)、BIC(Bayesian information criterion)的信息準(zhǔn)則函數(shù)來確定階數(shù)。

自回歸模型的確立過程，是通過確定階數(shù)，參數(shù)估計(jì)，再次確定階數(shù)的方法進(jìn)行判斷。自相關(guān)函數(shù)ACF，用來確定采用自回歸模型是否合適。如果自相關(guān)函數(shù)具有拖尾性，則AR模型為合適模型。偏自相關(guān)函數(shù)PACF用來確定模型的階數(shù)，如果從某個(gè)階數(shù)之后，偏自相關(guān)函數(shù)的值都很接近0，則取相應(yīng)的階數(shù)作為模型階數(shù)，偏自相關(guān)函數(shù)通過截尾性確定階數(shù)。

1. 自相關(guān)函數(shù)ACF(autocorrelation function)

將一個(gè)有序的隨機(jī)變量序列與其自身相比較，這就是自相關(guān)函數(shù)在統(tǒng)計(jì)學(xué)中的定義。每個(gè)不存在相位差的序列，都與其自身相似，即在此情況下，自相關(guān)函數(shù)值最大。如果序列中的組成部分相互之間存在相關(guān)性（不再是隨機(jī)的），則由以下相關(guān)值方程所計(jì)算的值不再為零，這樣的組成部分為自相關(guān)。

自相關(guān)函數(shù)反映了同一序列在不同時(shí)序的取值之間的相關(guān)程序。

ACF的公式為：

字段解釋

Pk,為ACF的標(biāo)準(zhǔn)誤差

t,為數(shù)據(jù)集的長度

k,為滯后，取值從1到t-1，表示相距 k個(gè)時(shí)間間隔的序列值之間的相關(guān)性

Yt,為樣本在t時(shí)期的值

Yt-k,為樣本在t-k時(shí)期的值

μ,為樣本的均值

所得的自相關(guān)值Pk的取值范圍為[-1,1],1為最大正相關(guān)值，-1則為最大負(fù)相關(guān)值，0為不相關(guān)。

根據(jù)上面公式，我們可以手動計(jì)算出tsx數(shù)據(jù)集的ACF值

> u<-mean(tsx) #均值> v<-var(tsx) #方差> # 1階滯后> p1<-sum((x[1:length(tsx)-1]-u)*(x[2:length(tsx)]-u))/((length(tsx)-1)*v);p1[1] 0.9878619> # 2階滯后> p2<-sum((x[1:(length(tsx)-2)]-u)*(x[3:length(tsx)]-u))/((length(tsx)-1)*v);p2[1] 0.9760271> # 3階滯后> p3<-sum((x[1:(length(tsx)-3)]-u)*(x[4:length(tsx)]-u))/((length(tsx)-1)*v);p3[1] 0.9635961

同時(shí)，我們可以用R語言中的acf()函數(shù)來計(jì)算，會打印前30個(gè)滯后的ACF值。

> acf(tsx)$acf, , 1 [,1] [1,] 1.0000000 [2,] 0.9878619 [3,] 0.9760271 [4,] 0.9635961 [5,] 0.9503371 [6,] 0.9384022 [7,] 0.9263075 [8,] 0.9142540 [9,] 0.9024862[10,] 0.8914740[11,] 0.8809663[12,] 0.8711005[13,] 0.8628609[14,] 0.8544984[15,] 0.8462270[16,] 0.8384758[17,] 0.8301834[18,] 0.8229206[19,] 0.8161523[20,] 0.8081941[21,] 0.8009467[22,] 0.7942255[23,] 0.7886249[24,] 0.7838154[25,] 0.7789733[26,] 0.7749697[27,] 0.7709313[28,] 0.7662547[29,] 0.7623381[30,] 0.7604101[31,] 0.7577333

比較前3個(gè)值的計(jì)算結(jié)果，與我們自己的計(jì)算結(jié)果是一樣的，同時(shí)可以用R語言進(jìn)行可視化輸出。

> acf(tsx)

從上圖中看出，數(shù)據(jù)的ACF為拖尾，存在很嚴(yán)重的自相關(guān)性。接下來，這時(shí)候我們用偏自相關(guān)函數(shù)確定一下AR的階數(shù)。

2. 偏自相關(guān)函數(shù)(PACF)(partial autocorrelation function)

偏自相關(guān)函數(shù)是有自相關(guān)函數(shù)推到而來。對于一個(gè)平穩(wěn)AR(p)模型，求出滯后k自相關(guān)系數(shù)p(k)時(shí)，實(shí)際上得到并不是x(t)與x(t-k)之間單純的相關(guān)關(guān)系。因?yàn)閤(t)同時(shí)還會受到中間k-1個(gè)隨機(jī)變量x(t-1)、x(t-2)、……、x(t-k+1)的影響，而這k-1個(gè)隨機(jī)變量又都和x(t-k)具有相關(guān)關(guān)系，所以自相關(guān)系數(shù)p(k)里實(shí)際摻雜了其他變量對x(t)與x(t-k)的影響。

為了能單純測度x(t-k)對x(t)的影響，引進(jìn)偏自相關(guān)系數(shù)的概念。對于平穩(wěn)時(shí)間序列{x(t)}，所謂滯后k偏自相關(guān)系數(shù)指在給定中間k-1個(gè)隨機(jī)變量x(t-1)、x(t-2)、……、x(t-k+1)的條件下，或者說，在剔除了中間k-1個(gè)隨機(jī)變量x(t-1)、x(t-2)、……、x(t-k+1)的干擾之后，x(t-k)對x(t)影響的相關(guān)程度。

簡單來說，就是自相關(guān)系數(shù)ACF還包含了其他變量的影響，而偏自相關(guān)系數(shù)PACF是嚴(yán)格這兩個(gè)變量之間的相關(guān)性。在ACF中存在著線性關(guān)系和非線性關(guān)系，偏自相關(guān)函數(shù)就是把線性關(guān)系從自動關(guān)系性中消除。當(dāng)PACF近似于0，表明兩個(gè)時(shí)間點(diǎn)之間的關(guān)系性是完全由線性關(guān)系所造成的。

通過R語言的pacf()函數(shù)來進(jìn)行偏自相關(guān)函數(shù)計(jì)算。

> pacf(tsx)$acf, , 1 [,1] [1,] 0.987861891 [2,] 0.006463542 [3,] -0.030541593 [4,] -0.041290415 [5,] 0.047921168 [6,] -0.009774246 [7,] -0.006267004 [8,] 0.002146693 [9,] 0.028782423[10,] 0.014785187[11,] 0.019307564[12,] 0.060879259[13,] -0.007254278[14,] -0.004139848[15,] 0.015707900[16,] -0.018615370[17,] 0.037067452[18,] 0.019322565[19,] -0.048471479[20,] 0.023388065[21,] 0.027640953[22,] 0.051177900[23,] 0.028063875[24,] -0.003957142[25,] 0.034030631[26,] 0.004270416[27,] -0.029613088[28,] 0.033715973[29,] 0.092337583[30,] -0.031264028# 可視化輸出 > pacf(tsx)

從上面的這個(gè)結(jié)果分析，當(dāng)滯后為1時(shí)AR模型顯著，滯后為其他值是PACF的值接近于0不顯著。所以，對于數(shù)據(jù)集tsx來說，數(shù)據(jù)滿足AR(1)的自回歸模型。對于上文中參數(shù)估計(jì)出的1階自相關(guān)系數(shù)值是可以用的。

4. 模型預(yù)測

通過模型識別，我們已經(jīng)確定了數(shù)據(jù)集tsx是符合AR(1)的建模條件的，同時(shí)我們也創(chuàng)建了AR(1)模型。接下來，就可以利用這個(gè)自回測的模型的進(jìn)行預(yù)測，通過規(guī)律發(fā)現(xiàn)價(jià)值。在R語言中，我們可以用predict()函數(shù)，實(shí)現(xiàn)預(yù)測的計(jì)算。

使用AR(1)模型進(jìn)行預(yù)測，并保留前5個(gè)預(yù)測點(diǎn)。

> predict(a,10,n.ahead=5L)$predTime Series:Start = 2 End = 6 Frequency = 1 [1] 9.839680 9.681307 9.524855 9.370303 9.217627$seTime Series:Start = 2 End = 6 Frequency = 1 [1] 1.080826 1.519271 1.849506 2.122810 2.359189

上面結(jié)果中，變量$pred表示預(yù)測值，變量$se為誤差。

我可以生成可視化的圖，更直觀的看到預(yù)測的結(jié)果。

# 生成50個(gè)預(yù)測值 > tsp<-predict(a,n.ahead=50L)# 把原數(shù)據(jù)畫圖 > plot(tsx)# 把預(yù)測值和誤差畫出來> lines(tsp$pred,col='red') > lines(tsp$pred+tsp$se,col='blue')> lines(tsp$pred-tsp$se,col='blue')

圖中，黑色線為原始數(shù)據(jù)的，紅色線為預(yù)測值，藍(lán)色線為預(yù)測值的范圍。這樣我們就利用AR(1)模型，實(shí)現(xiàn)了對規(guī)律的預(yù)測計(jì)算。

上面關(guān)于預(yù)測和可視化的過程,我們是通過原生的predict()函數(shù)和plot()函數(shù)完成的。在R語言中，可以用forecast包來簡化上面的操作過程，讓代碼更少，操作更便捷。

# 加載forecast包> library('forecast')# 生成模型AR(1) > a2 <- arima(tsx, order=c(1,0,0))> tsp2<-forecast(a2, h=50)> plot(tsp2)

查看forecast()計(jì)算后的預(yù)測結(jié)果。

> tsp2 Point Forecast Lo 80 Hi 80 Lo 95 Hi 951001 -15.71590 -16.99118 -14.440628 -17.66627 -13.76553691002 -15.60332 -17.39825 -13.808389 -18.34843 -12.85820921003 -15.49181 -17.67972 -13.303904 -18.83792 -12.14569661004 -15.38136 -17.89579 -12.866932 -19.22685 -11.53587261005 -15.27197 -18.06994 -12.474000 -19.55110 -10.99284321006 -15.16362 -18.21425 -12.112996 -19.82915 -10.49809221007 -15.05631 -18.33593 -11.776682 -20.07206 -10.04055411008 -14.95001 -18.43972 -11.460312 -20.28705 -9.61297501009 -14.84474 -18.52891 -11.160567 -20.47919 -9.21028461010 -14.74046 -18.60591 -10.875013 -20.65216 -8.82876731011 -14.63718 -18.67257 -10.601802 -20.80877 -8.46559941012 -14.53489 -18.73030 -10.339486 -20.95121 -8.11857231013 -14.43357 -18.78024 -10.086905 -21.08123 -7.78591741014 -14.33322 -18.82333 -9.843112 -21.20026 -7.46619031015 -14.23383 -18.86034 -9.607319 -21.30947 -7.15819231016 -14.13538 -18.89190 -9.378864 -21.40985 -6.8609139

通過forecast()函數(shù)，直接生成了Forecast值，80%概率的預(yù)測值范圍，和95%概率的預(yù)測值范圍。

在明白了整個(gè)自回歸模型的設(shè)計(jì)思路、建模過程、檢驗(yàn)條件、預(yù)測計(jì)算、可視化展示的完整操作后，我們就可以真正地把自回歸模型用到實(shí)際的業(yè)務(wù)中。發(fā)現(xiàn)規(guī)律，發(fā)現(xiàn)價(jià)值??！