第四章 平稳序列的拟合与预测
第四章 平穩序列的擬合與預測
4.1 建模步驟
當我們拿到一組數據,想要通過這組數據來探究過去發生的規律,研究現在發生的強度并預測未來事情的發展走向,就需要從這些數據中挖掘出事物發生的特征,比如頻率、強度、集中值等等,來提取“經驗”,對過去的作出反應,彌補或者提取收益,而后建立一個合適的模型,來預測未來事物的發展,掌握未來動態,并及時調整,做好準備。
對于時間序列數據,當我們拿到數據時,我們首先要對數據進行預處理,事實上,所有的數據都是如此,經過預處理使之成為合適的,可以進行預測的,有用的數據。
時間序列數據的預處理是要進行平穩性檢驗,只有平穩的時間序列數據才是可以利用的,有效的數據,才能進行預測。還要進行純隨機檢驗,確保其是非白噪聲序列。
4.2 模型識別
若數據是平穩的非白噪聲序列,則可以考慮建模,但時間序列模型不止一種,具體要用哪種模型要根據數據的特征來判斷,即先計算樣本的ACF和PACF。根據樣本的自相關系數和偏自相關系數的截尾和拖尾性來選擇合適的模型擬合觀察值序列。
定階的基本原則如下:
| 拖尾 | p階截尾 | AR(p)模型 |
| q階截尾 | 拖尾 | MA(q)模型 |
| 拖尾 | 拖尾 | ARMA(p,q)模型 |
但是由于樣本的隨機性,樣本的相關系數不會呈現出理論上的截尾的完美情況,應該截尾的樣本自相關系數或偏自相關系數就會呈現出小值震蕩。同時,又因為平穩時間序列通常具有的短期相關性,隨著延遲階數趨近于無窮,ρk^\hat{\rho_k}ρk?^?與?kk^\hat{\phi_{kk}}?kk?^?都會衰減至零值附近作小值波動。
這就給我們留下了一個問題,判斷拖尾的標準情況是什么?其實并沒有,主要還是依靠分析人員的主觀經驗。根據諸多統計學家的證明,樣本自相關系數是總體自相關系數的有偏估計,當k足夠大時,根據平穩序列自相關系數呈負指數衰減,有ρk→0\rho_k\rightarrow0ρk?→0,而當樣本容量n充分大時,樣本的自相關系數和偏自相關系數都是近似于正態分布,根據正態分布的性質,可以利用2倍標準差范圍輔助判斷。
如果樣本自相關系數或偏自相關系數在最初的d階明顯超過2倍標準差范圍,而后幾乎95%的自相關系數都落在2倍標準差的范圍以內,而且由非零自相關系數衰減為小值波動的過程非常突然,這時,通常視為自相關系數截尾。截尾階數為d。
如果有超過5%的樣本自相關系數落入2倍標準差范圍之外,或者由顯著非零的自相關系數衰減為小值波動的過程非常連續或比較緩慢,這時,通常視為自相關系數不截尾。
4.3 參數估計
對于非中心化的ARMA(p,q)模型,有xt=μ+Θq(B)Φq(B)?tx_t=\mu+\frac{\Theta_q(B)}{\Phi_q(B)}\epsilon_txt?=μ+Φq?(B)Θq?(B)??t?,式中?t~WN(0,σ?2)\epsilon_t\thicksim WN(0,\sigma^2_{\epsilon})?t?~WN(0,σ?2?)【WN白噪聲】。模型中共含有p+q+2個未知參數:θ1,θ2,?,θq,?1,?2,?,?p,μ,σ?2\theta_1,\theta_2,\cdots,\theta_q,\phi_1,\phi_2,\cdots,\phi_p,\mu,\sigma_{\epsilon}^2θ1?,θ2?,?,θq?,?1?,?2?,?,?p?,μ,σ?2?。
參數μ\muμ是序列均值,通常采用一階中心距的估計方法,用樣本均值估計總體均值即可得到它的估計值。
再將原序列中心化,則待估參數減少至p+q+1個。對這些參數有三種估計方法:距估計、極大似然估計、最小二乘估計。
4.3.1 距估計
用p+q個樣本自相關系數估計總體自相關系數,構造p+q個方程組成的Yule-Walker方程組
{ρ1(?1,?,?p,θ1,?,θq)=ρ1^?ρp+q(?1,?,?p,θ1,?,θq)=ρp+q^\begin{cases} \rho_1(\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q)=\hat{\rho_1}\\ \vdots\\ \rho_{p+q}(\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q)=\hat{\rho_{p+q}}\\ \end{cases} ????????ρ1?(?1?,?,?p?,θ1?,?,θq?)=ρ1?^??ρp+q?(?1?,?,?p?,θ1?,?,θq?)=ρp+q?^??
從中解出的參數值?1^,?,?p^,θ1^,?,θq^\hat{\phi_1},\cdots,\hat{\phi_p},\hat{\theta_1},\cdots,\hat{\theta_q}?1?^?,?,?p?^?,θ1?^?,?,θq?^?就是?1,?,?p,θ1,?,θq\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q?1?,?,?p?,θ1?,?,θq?的距估計。
尤爾—沃克方程(Yule-Walker equation)是描述自回歸序列參數與其協方差函數之間關系的方程。
所以,距估計其實就是建立多個自相關系數與參數的方程,而后解出參數(因為ρ\rhoρ可估或者說是已知)。
4.3.2 極大似然估計
在極大似然準則下(事件既然發生,則其概率就一定盡可能大),認為樣本來自使該樣本出現大概率最大的總體。因此未知參數的極大似然估計就是使得似然函數(即聯合密度函數)達到最大的參數值。
L(β1^,?,βk^;x1,?,xn)=max{p(x1,?,xn);β1,?,βk}L(\hat{\beta_1},\cdots,\hat{\beta_k};x_1,\cdots,x_n)=max\{p(x_1,\cdots,x_n);\beta_1,\cdots,\beta_k\} L(β1?^?,?,βk?^?;x1?,?,xn?)=max{p(x1?,?,xn?);β1?,?,βk?}
但是,使用極大似然估計必須已知總體的分布函數,而在時間序列分析中,序列總體分布通常是未知的。為了便于分析計算,通常假設序列服從多元正態分布。
獨立的多元正態分布:
以上在理論上是可行的,但實際上上式存在p+q個超越方程,通常需要利用迭代法才能求出未知參數的極大似然估計值。實際操作中可以借助計算機。
4.3.3 最小二乘估計
在ARMA(p,q)模型場合,記
β~=(?1,?,?p,θ1,?,θq)′Ft(β~)=?1xt?1+?+?pxt?p?θ1?t?1???θq?t?q\tilde{\beta}=(\phi_1,\cdots,\phi_p,\theta_1,\cdots,\theta_q)'\\ F_t(\tilde{\beta})=\phi_1x_{t-1}+\cdots+\phi_px_{t-p}-\theta_1\epsilon_{t-1}-\cdots-\theta_q\epsilon_{t-q} β~?=(?1?,?,?p?,θ1?,?,θq?)′Ft?(β~?)=?1?xt?1?+?+?p?xt?p??θ1??t?1????θq??t?q?
殘差項為:
?t=xt?Ft(β~)\epsilon_t=x_t-F_t(\tilde{\beta}) ?t?=xt??Ft?(β~?)
殘差平方和為:
Q(β~)=∑t=1n?t2=∑t=1n(xt??1xt?1+?+?pxt?p?θ1?t?1???θq?t?q)2Q(\tilde{\beta})=\sum_{t=1}^{n}\epsilon_t^2\\ =\sum_{t=1}^{n}(x_t-\phi_1x_{t-1}+\cdots+\phi_px_{t-p}-\theta_1\epsilon_{t-1}-\cdots-\theta_q\epsilon_{t-q})^2 Q(β~?)=t=1∑n??t2?=t=1∑n?(xt???1?xt?1?+?+?p?xt?p??θ1??t?1????θq??t?q?)2
使殘差平方和達到最小的那組參數值即β~\tilde{\beta}β~?的最小二乘估計。
在實際中,最常用的是條件最小二乘估計方法。它假定過去未觀測到的序列值等于零,于是有
?t=Φ(B)Θ(B)xt=xt?∑i=1tπixt?i\epsilon_t=\frac{\Phi(B)}{\Theta(B)}x_t=x_t-\sum_{i=1}^{t}\pi_ix_{t-i} ?t?=Θ(B)Φ(B)?xt?=xt??i=1∑t?πi?xt?i?
于是殘差平方和為:
Q(β~)=∑t=1n?t2=∑t=1n[xt?∑i=1tπixt?i]2Q(\tilde{\beta})=\sum_{t=1}^{n}\epsilon_t^2\\ =\sum_{t=1}^{n}[x_t-\sum_{i=1}^{t}\pi_ix_{t-i}]^2 Q(β~?)=t=1∑n??t2?=t=1∑n?[xt??i=1∑t?πi?xt?i?]2
4.4 模型檢驗
4.4.1 模型的顯著性檢驗
模型的顯著性檢驗主要是檢驗模型的有效性。且一個模型是否顯著有效主要看它提取的信息是否充分。也即是殘差項應該是白噪聲序列。
4.4.2 參數的顯著性檢驗
參數的顯著性檢驗就是檢驗每一個未知參數是否顯著非零。
4.5 模型優化
有時同一個序列可以有不同的擬合模型,且都顯著有效,那么我們怎么選擇呢?
這時就要用到AIC準則和SBC準則,即自小信息量準則,一般用軟件擬合時會直接給出。
AIC=?2ln(模型的極大似然函數值)+2(模型中未知參數個數)AIC=-2ln(模型的極大似然函數值)+2(模型中未知參數個數) AIC=?2ln(模型的極大似然函數值)+2(模型中未知參數個數)
AIC準則的局限性在于受到樣本容量的限制,當樣本容量趨于無限大的時候,由AIC準則選擇的模型不收斂于真實模型,它通常比真實模型所含的未知參數個數還要多。
BIC準則也稱SBC準則彌補了這個不足。
SBC=?2ln(模型的極大似然函數值)+(lnn)(模型中未知參數個數)SBC=-2ln(模型的極大似然函數值)+(ln n)(模型中未知參數個數) SBC=?2ln(模型的極大似然函數值)+(lnn)(模型中未知參數個數)
4.6 序列預測
根據ARMA(p,q)模型的平穩性和可逆性,可以用傳遞形式和逆轉形式等價描述該序列:
xt=∑i=0∞Gi?t?i?t=∑j=0∞Ijxt?jx_t=\sum_{i=0}^{\infty}{G_i\epsilon_{t-i}}\\ \epsilon_t=\sum_{j=0}^{\infty}{I_jx_{t-j}} xt?=i=0∑∞?Gi??t?i??t?=j=0∑∞?Ij?xt?j?
其中{Gi}\{G_i\}{Gi?}是Green函數值;{Ij}\{I_j\}{Ij?}為逆轉函數值。
將兩式合并,有
xt=∑i=0∞∑j=0∞GiIjxt?i?jx_t=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}G_iI_jx_{t-i-j} xt?=i=0∑∞?j=0∑∞?Gi?Ij?xt?i?j?
簡記為
xt=∑i=0∞Cixt?1?ix_t=\sum_{i=0}^{\infty}C_ix_{t-1-i} xt?=i=0∑∞?Ci?xt?1?i?
?l?1\forall l\geqslant1?l?1時,對任意未來t+l時刻,該時刻的序列值可以表示為它的歷史數據xt+l?1,?,xt+1,xt,?x_{t+l-1},\cdots,x_{t+1},x_t,\cdotsxt+l?1?,?,xt+1?,xt?,?的線性函數
xt+l=∑i=0∞Cixt+l?1?ix_{t+l}=\sum_{i=0}^{\infty}{C_i}x_{t+l-1-i} xt+l?=i=0∑∞?Ci?xt+l?1?i?
但有趣的是xtx_txt?前的歷史信息已知,xt+1,?,xt+l?1x_{t+1},\cdots,x_{t+l-1}xt+1?,?,xt+l?1?期的信息是未知的,但根據線性函數的可加性,所有的未知歷史信息都可以用已知歷史信息的線性函數表示出來。
如xt+2=∑i=0∞Cixt+1?i=C0xt+1+∑i=0∞Ci+1xt?ix_{t+2}=\sum_{i=0}^{\infty}C_ix_{t+1-i}=C_0x_{t+1}+\sum_{i=0}^{\infty}C_{i+1}x_{t-i}xt+2?=∑i=0∞?Ci?xt+1?i?=C0?xt+1?+∑i=0∞?Ci+1?xt?i?,但是因為xt+1x_{t+1}xt+1?是未知的,而其又可以表示為xt+1=∑i=0∞Cixt?ix_{t+1}=\sum_{i=0}^{\infty}C_ix_{t-i}xt+1?=∑i=0∞?Ci?xt?i?,故有
xt+2=∑i=0∞(C0Ci+Ci+1)xt?ix_{t+2}=\sum_{i=0}^{\infty}(C_0C_i+C_{i+1})x_{t-i} xt+2?=i=0∑∞?(C0?Ci?+Ci+1?)xt?i?
由此xt+2x_{t+2}xt+2?最終表示為已知歷史信息xt,xt?1,?x_t,x_{t-1},\cdotsxt?,xt?1?,?的線性函數。同理,對于未來任意l時刻的序列值xt+lx_{t+l}xt+l?最終都可以用已知歷史信息xt,xt?1,?x_t,x_{t-1},\cdotsxt?,xt?1?,?表示出來,并將該線性函數形式記為:
xt^(l)=∑i=0∞D^ixt?i\hat{x_t}(l)=\sum_{i=0}^{\infty}{\hat{D}_ix_{t-i}} xt?^?(l)=i=0∑∞?D^i?xt?i?
x^t(l)\hat{x}_t(l)x^t?(l)也稱為序列{xt}\{x_t\}{xt?}的第l步預測值。
代碼:
A7<-read.csv(“A1_7.CSV”)
A8<-read.csv(“A1_8.csv”)
A9<-read.csv(“A1_9.csv”)
a7<-ts(A7number,start=1900,frequency=1)a8<?ts(A8number,start = 1900,frequency = 1) a8<-ts(A8number,start=1900,frequency=1)a8<?ts(A8overshort,start = 1,frequency = 365)
a9<-ts(A9$change,start = 1880,frequency = 1)
a<-diff(a9,lag = 1)#對數據進行一階差分
#install.packages(“aTSA”)
library(aTSA)
#class(A7)
adf.test(a7)#平穩性檢驗
acf(a7,lag=25)#繪制自相關圖
acf(a7,lag=25)KaTeX parse error: Expected 'EOF', got '#' at position 4: acf#?顯示自相關系數 拖尾 pac…pacf#顯示偏自相關系數 一階截尾
plot(a8)
adf.test(a8)#ADF檢驗
Box.test(a8,lag = 6)
Box.test(a8,lag = 12)
for (k in 1:2) {print(Box.test(a8,lag = 6k,type = “Ljung-Box”))
}
acf(a8,lag=20)
pacf(a8)
plot(a9)
adf.test(a9)
plot(a)
acf(a)
pacf(a)
for (k in 1:2) {print(Box.test(a9,lag = 6k,type = “Ljung-Box”))
}
par(mfrow=c(1,2))
acf(a9)
pacf(a9)
###參數估計,使用函數arima()擬合
fit1<-arima(a7,order = c(1,0,0),include.mean=T,method = “CSS-ML”)#極大似然估計與條件最小二乘估計混合
fit1
fit1<-arima(a7,order = c(1,0,0),method = “CSS”)
fit2<-arima(a8,order = c(0,0,1),method = “CSS”)#條件最小二乘估計
fit2
fit3<-arima(a,order = c(1,0,1))
fit3
#####模型檢驗,使用aTSA包中的ts.diag函數
library(aTSA)
ts.diag(fit1)#模型的顯著性檢驗
t=abs(fit1coef)/sqrt(diag(fit1coef)/sqrt(diag(fit1coef)/sqrt(diag(fit1var.coef))
pt(t,length(a7)-length(fit1$coef),lower.tail = F)
###模型預測,使用forecast包中的forecast函數
install.packages(“forecast”)
#library(forecast)
forecast::forecast##因前面打開aTSA包
fore1<-forecast::forecast(fit1,h=10)
fore1
plot(fore1,lty=2)
lines(fore1$fitted,col=4)
未完…
總結
以上是生活随笔為你收集整理的第四章 平稳序列的拟合与预测的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: CS免杀
- 下一篇: 超级SIM卡 SEID号读取 手机NFC