平穩時間序列以及MATLAB相關工具箱學習筆記

概念

（1）平穩序列

即序列的均值是個常數，與序列長度、起始位置無關。
直觀看上去，該序列類似于圍繞某一個值上下波動（該值為平穩序列均值）。

（2）平穩白噪聲序列

所謂噪聲，可以理解為一種非周期性的擾動。由于其自協方差函數值為0，于是其各序列值沒有任何相關關系，所以說平穩白噪聲序列是一段沒有記憶的平穩序列。
在MATLAB中可用如下代碼生成高斯分布下的平穩白噪聲序列：

elps=randn(1,1000); plot(elps)

基本的平穩時間序列

（1）$AR(p)$序列

需要說明的是：
1.$p$為模型階數，需要進行定階。
比說：隨機產生時間序列$X_t+0.6X_{t?1}+0.2X_{t?2}={\xi }_t$中
$p=2$,則該序列為$AR(2)$序列。

2.$\psi =({\psi }_1,{\psi }_2,?,{\psi }_p{)}^T$為未知參數，需要對其進行估計。

（2）$MA(q)$序列

同樣：
1.$q$為模型階數，需要進行定階。
2.$\theta =({\theta }_1,{\theta }_2,?,{\theta }_q{)}^T$為未知參數，需要對其進行估計。

（3）$ARMA(p,q)$序列

1.$p,q$為模型階數，需要進行定階。
2.$\theta =({\theta }_1,{\theta }_2,?,{\theta }_q{)}^T,\psi =({\psi }_1,{\psi }_2,?,{\psi }_p{)}^T$為未知參數，需要對其進行估計。

平穩時間序列建模流程

下面逐步分析流程并演示代碼：

（1）平穩性檢驗

先導入一組數據，生成圖像：

由圖可知，其變動范圍較大，均值可能為一常數（肉眼判斷）在定義上有可能滿足平穩序列的條件。
于是，引出了平穩性檢驗的第一種方法：$1.時序圖檢驗$。
除了這種方法，更一般的方法為$2.自相關檢驗$
平穩序列通常具有短期相關性，該性質用自相關系數來描述就是：$隨著延遲期數的增加，平穩序列的自相關系數會很快衰減為0$。

MATLAB檢驗自相關系數的函數為：

figure(1) r11=autocorr(a) %計算自相關函數 stem(r11)%繪制經線圖 title('自相關系數') figure(2) r12=parcorr(a) %計算偏相關函數 stem(r12)%繪制經線圖 title('偏相關系數')

對源數據計算可得圖像：

平穩系序列的自相關系數具有：截尾性、拖尾性、指數衰減性。根據以上兩圖可以進行直觀理解。從以上圖示，$該序列都滿足平穩性的判別標準$。

但如上圖所示，該序列的自相關系數分布圖明顯不展現截尾性、拖尾性、指數衰減性的性質，該序列為不平穩序列。
$3.Daniel檢驗法$。

MATLAB程序：

x0=[80.8 94.0 88.4 101.5 110.3 121.5 134.7 142.7]; n=length(x0); [xsort,ind]=sort(x0) Rt(ind)=1:n %計算秩 t=1:n; Qs=1-6/(n*(n^2-1))*sum((t-Rt).^2) %計算 Qs 的值 T=Qs*sqrt(n-2)/sqrt(1-Qs^2) %計算 T 統計量的值 t_0=tinv(0.975,n-2) %計算上 alpha/2 分位點

計算結果：
$T=11.0235>t_0=2.4469$該序列為$非平穩性序列$。
其次，若已知序列的類型($AR(p)$序列、$MA(q)$序列、$ARMA(p,q)$序列），可采用$4.單位根法$進行平穩性檢驗。
$\theta (B)及\psi （B）$的根均在單位圓內，則該序列平穩。
例子：已知隨機隨機時間序列：
$X_t+0.6X_{t?1}+0.2X_{t?2}={\xi }_t$
則$\theta (B)=1+0.6B+0.2B^2$

roots([0.2 0.6 1])

（2）構造差分序列

需要指出，并非是所有的非平穩序列都能進行差分使其平穩化，能進行差分的序列須為$帶有一定趨勢$的非平穩差分序列。

da=diff(a); %計算1階差分 r21=autocorr(da) %計算自相關函數 r22=parcorr(da) %計算偏相關函數 figure(1) stem(r21) title('1階差分自相關') figure(2) stem(r22) title('1階差分自相關')

以上兩幅圖為原平穩序列差分后的結果，由圖可知平穩序列的差分序列不會喪失其平穩特性，反倒越加平穩。若1階差分滿足不了，可采用2階甚至多階。

（3）模型定階與擬合（參數估計）

一般地，若已知該序列的類型，可采用直接指定的方式定階。還是這個例子：
$X_t+0.6X_{t?1}+0.2X_{t?2}={\xi }_t$

p=3;%指定大概率會的模型階數 elps=randn(10000,1); x(1:2)=0; for i=3:10000x(i)=-0.6*x(i-1)-0.2*x(i-2)+elps(i); %產生模擬數據 end x=x'; %x必須為列向量 for i=1:pmodel{i}=ar(x,i);%采用工具箱里面的AR函數進行擬合 AIC(i)=aic(model{i});%根據AIC準則，確定最終 end N=min(AIC); P_great=find(AIC==N);

更一般地，對于一組未知類型的時間序列（已知某化工生產濃度數據、已知1997-2012年我國人均國內生產總值、1974-1981年布料銷售量$?$，可采用試探模型（即通過構造{$p=1:m,q=1:n$}的方法)進行定階。這里，我們可以假設模型為計量經濟學領域的CARCH模型：

這里我們需要用到MATLAB的兩個常用模型指定函數：$$model=arima{(}^{\prime }MALag{s}^{\prime },[MA模型滯后項]{,}^{\prime }ARLag{s}^{\prime },[AR模型滯后項]{,}^{\prime }M{A}^{\prime },MA模型滯后項系數{,}^{\prime }A{R}^{\prime },MA模型滯后項系數{,}^{\prime }Constan{t}^{\prime },常數項）$$
$X_t+0.6X_{t?1}+0.2X_{t?2}={\xi }_t$
可用如下語句進行指定：

model=arima('ARLags',[1 2 3],'MAlags',1','Constant',0)

另外，指定模型之后需要進行擬合，具體函數如下：

輸入：Md1為上面用arima指定的模型；y為需要擬合的時間序列；
輸出：EstMd1為擬合得到的模型方程；EstParamCov為返回與估計參數；logL為目標函數值；info為匯總信息。
一段完整的定階與擬合MATLAB程序如下：

for i=0:pfor j=0:qif i==0 && j==0continueelseif i==0ToEstMd=arima('MALags',1:j);elseif j==0ToEstMd=arima('ARLags',1:i);elseToEstMd=arima('ARLags',1:i,'MALags',1:j); end[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,d);%模型擬合numParams=sum(any(EstParamCov));%計算擬合參數個數[aic(i,j),bic(i,j)]=aicbic(logL,numParams,n);end end N=min(min(bic)); [p_great,q_great]=find(N==bic);

（4）模型預報

這里也將有大篇幅的理論推導，這里將不再敘述。
AR預報：
ARMA預報：

forecast_=forecast(EstMd,5,'Y0',d)

這里，EstMd為擬合出來的模型，5為需要預報的值，’Y0‘為預測參量，d為序列。

（5）模型檢驗

模型檢驗常用的方法有：
時間序列檢驗
殘差的 Portmanteau 檢驗
Lagrange 乘數檢驗
MATLAB殘差檢驗的程序如下：

res=infer(EstMd,d);%計算殘差h=lbqtest(res);%模型檢驗

總結

以上是生活随笔為你收集整理的平稳时间序列以及MATLAB相关工具箱学习笔记的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。