1、1,慣性儀器測試與數(shù)據(jù)分析,西北工業(yè)大學(xué) 自動化學(xué)院 嚴恭敏 2015-09,2,第六章 時間序列分析,主要內(nèi)容： 一、隨機過程的基本概念 二、ARMA模型及其特點 三、ARMA建模分析,3,第六章 時間序列分析,4,一、隨機過程的基本概念,1、隨機向量,在概率論中，隨機變量 用來描述隨機事件，可分為續(xù)型隨機變量，離散型隨機變量。,（1）隨機變量,使用概率密度函數(shù)、概率分布函數(shù)、特征函數(shù)及數(shù)字特征（均值、方差和矩等）等數(shù)學(xué)語言描述。,（2）隨機向量,由 N 個隨機變量組成一組向量：,均值向量,方差矩陣,5,一、隨機過程的基本概念,2、隨機過程與時間序列,電阻熱噪聲電壓,A）每只 電阻電壓隨時間

2、是一條隨機波動的曲線,B）在同一特定 時刻各個電阻的電壓值各不相同,樣本曲線（軌跡、現(xiàn)實）,隨機變量 取值,A）所有樣本函數(shù)的集合 構(gòu)成了一個隨機過程；,（3）隨機過程 定義：,6,一、隨機過程的基本概念,2、隨機過程與時間序列,隨機變量取值,（4）隨機過程分類與時間序列,時間參數(shù)取值,連續(xù),離散,連續(xù),離散,時間序列,使用A/D轉(zhuǎn)換器對熱電阻電壓進行等間隔 采樣，假設(shè)A/D分辨率足夠高忽略量化誤差影響,將采樣周期歸一化處理，常將時間序列記為,時間序列就是按照時間的先后順序記錄的一列有序數(shù)據(jù)，這些數(shù)據(jù)由于受到各種偶然因素的影響，往往表現(xiàn)出某種隨機性，但彼此之間又存在一定的相關(guān)性； 時間序列分析

3、就是對時間序列進行觀察、研究，揭示其蘊含的內(nèi)在規(guī)律，進而根據(jù)變化規(guī)律預(yù)測走勢或?qū)嵤┛刂啤?時間序列分析方法大體可分為時域和頻域兩種分析方法。時域分析方法主要從序列自相關(guān)的角度揭示時間序列的發(fā)展規(guī)律；頻域分析方法也稱為頻譜分析，從頻率角度揭示時間序列的規(guī)律。,7,一、隨機過程的基本概念,2、隨機過程與時間序列,（5）時間序列的數(shù)字特征,均值序列,自協(xié)方差函數(shù),自相關(guān)系數(shù)函數(shù),方差函數(shù),自相關(guān)函數(shù),恒等式,8,一、隨機過程的基本概念,3、平穩(wěn)性與各態(tài)遍歷性,（1）嚴平穩(wěn)過程（狹義平穩(wěn)過程） 隨機過程的所有概率統(tǒng)計性質(zhì)不隨時間原點推移而變化，要求過于苛刻，不利于理論分析和實際應(yīng)用。,（2）寬平穩(wěn)過程

4、（廣義平穩(wěn)過程或二階矩平穩(wěn)過程 ） 在所有時刻上，均值序列和方差序列都是常值，且方差有限，即 和 ； 自相關(guān)函數(shù)與時間起點無關(guān)，而只與時間間隔 有關(guān)，即,例如，白噪聲序列 是理論分析中一種理想化的最基本的平穩(wěn)序列,正態(tài)分布，高斯白噪聲常記作,偶對稱,9,一、隨機過程的基本概念,3、平穩(wěn)性與各態(tài)遍歷性,（3）各態(tài)遍歷平穩(wěn)過程 對于平穩(wěn)過程，實際工作中通常很難取得足夠多的樣本用來分析隨機過程的總體特性，有時也是沒有必要的，所以常常只用少量甚至一個樣本函數(shù)進行分析，這就涉及到一個樣本函數(shù)的特性能否代表和估計隨機過程總體特性的問題。,若滿足所有樣本函數(shù)在某一固定時刻的一階和二階統(tǒng)計特性與單一樣本函數(shù)在

5、長時間的統(tǒng)計特性一致則稱為各態(tài)遍歷平穩(wěn)隨機過程，即,集總平均 = 時間平均,提出各態(tài)遍歷性目的：計算方便，但針對實際問題證明困難， 經(jīng)驗上：主要物理條件隨時間基本不變，各樣本隨機影響因素基本相同。,簡記作 ，可表示 1.時間序列總體 2.某一樣本 3. 時刻樣本值,10,一、隨機過程的基本概念,3、平穩(wěn)性與各態(tài)遍歷性,舉例,11,二、ARMA模型及其特點,建立模型意義： A）獲得一些重要的模型參數(shù)，有助于深入了解研究對象，為進一步改進研究對象提供依據(jù)； B）通過建立研究對象的數(shù)學(xué)表達式，是更好地發(fā)揮研究對象的使用性能的基礎(chǔ)，特別是在使用現(xiàn)代最優(yōu)控制和最優(yōu)估計理論解決實際問題時，對傳感器進行隨機

6、測量誤差建模分析具有重要意義。,（1）ARMA(p,q)、MA(q)與AR(p)模型定義,以零均值高斯白噪聲序列 作為時間序列分析的最基本組成單元，一般各態(tài)遍歷平穩(wěn)時間序列 使用白噪的線性組合聲來表示。,特征根,12,二、ARMA模型及其特點,A) ARMA(p,q) 模型（自回歸-滑動平均模型 ）,觀測值 與既往p個觀測 存在相關(guān)性，并且除 外與既往q個噪聲 也存在相關(guān)性,B) MA(q) 模型（滑動平均模型 ）,C) AR (p) 模型（自回歸模型 ）,滑動平均系數(shù),自回歸系數(shù),13,二、ARMA模型及其特點,（2）MA(q) 模型特點,自協(xié)方差函數(shù),自相關(guān)系數(shù)函數(shù),q步截尾,矩陣形式,相

7、關(guān)分析：已知 求,模型辨識：已知 求,14,二、ARMA模型及其特點,（3）AR(p) 模型特點,自協(xié)方差函數(shù),自相關(guān)系數(shù)函數(shù),矩陣形式,簡記為 ，即尤爾沃克方程（Y-W方程）,相關(guān)分析：已知 求,模型辨識：已知 求,15,二、ARMA模型及其特點,AR(p)自相關(guān)系數(shù)函數(shù)拖尾性,如 已知，當 時,研究表明， 按負指數(shù)函數(shù)衰減，理論上是無限延伸趨于0的，這種性質(zhì)稱為拖尾性。,為了判斷AR（p）過程的階數(shù)，引入偏自相關(guān)系數(shù)函數(shù)定義,其中,含義：扣除中間量 的影響后， 與 之間的相關(guān)性。,為最佳線性估計系數(shù)；,16,二、ARMA模型及其特點,按定義不好計算，研究發(fā)現(xiàn) 恰好與k階Y-W方程的解系數(shù)

8、完全相同,顯然 是 的函數(shù)。,容易驗證這兩個特例 和 ，并且 。,因此，AR(p)過程的偏自相關(guān)系數(shù)函數(shù)是p步截尾的，這是用它作為過程階數(shù)判斷的重要標志。,偏自相關(guān)系數(shù)函數(shù)的定義和計算方法也適用于MA(q)過程，但它按負指數(shù)函數(shù)衰減，也就是說MA(q)偏自相關(guān)系數(shù)函數(shù)具有拖尾性質(zhì),AR(p),MA(q),對偶性,自相關(guān)系數(shù)q步截尾,偏自相關(guān)系數(shù)拖尾,自相關(guān)系數(shù)拖尾,偏自相關(guān)系數(shù)p步截尾,17,二、ARMA模型及其特點,例6.2-1 假設(shè)AR(1)模型 ，白噪聲 ，試求該模型的自相關(guān)系數(shù)函數(shù)和偏自相關(guān)系數(shù)函數(shù)。,解：,首先,其次大于等于1階時由遞推方程,自協(xié)方差,綜合,18,二、ARMA模型及其

9、特點,直接由偏自相關(guān)系數(shù)函數(shù)性質(zhì) 和 得,當k1時,可以驗證,AR(1)過程也稱為一階馬爾科夫過程（Markov過程），主要特點是當前時刻觀測值僅與相鄰的前一時刻觀測值存在相關(guān)性。,序列采樣周期， 相關(guān)時間常數(shù)， 反相關(guān)時間常數(shù)（單位1/s）， 序列相關(guān)長度，含義,當 時 ，此時AR(1)過程與白噪聲非常接近。,19,二、ARMA模型及其特點,在AR(1)中令 ，考慮隨機過程模型,由于對應(yīng)傳遞函數(shù)的分母特征根 在單位圓上，因此該模型不屬于平穩(wěn)過程。 特點：當前觀測值完全由上一時刻觀測值加上現(xiàn)時噪聲決定隨機游走。,假設(shè),則,（一階馬爾科夫）,（隨機游走）,（白噪聲）,非平穩(wěn)過程,三者關(guān)系：,20

10、,二、ARMA模型及其特點,% MATLAB AR(1)仿真 a1 = 0.95; % AR(1)參數(shù) sigma = 1; % 噪聲均方差 N = 100; % 仿真長度 wn = sigma*randn(N,1); % 白噪聲序列 x = wn; % 初始值 for n = 2:N x(n) = a1*x(n-1) + wn(n); end plot(x,b), grid on, hold on % 數(shù)據(jù)圖 xlabel(n); ylabel(x);,AR(1)過程的MATLAB仿真程序,21,(1) a1=0.37 (2) a1=0.95 (3) a1=0.0 (4) a1=1.0,22

11、,(1) AR(1): a1=-0.95 (2) AR(1): a1=0.95 (3) MA(1): b1=-0.95 (4) MA(1): b1=0.95,23,二、ARMA模型及其特點,例6.2-2 假設(shè)AR(2)模型 ，白噪聲 ，試求該模型的自相關(guān)系數(shù)函數(shù)和偏自相關(guān)系數(shù)函數(shù)。,解：,首先,其次根據(jù)2階Y-W方程,最后高于1階時由遞推方程,直接由偏自相關(guān)系數(shù)函數(shù)性質(zhì)得 和,AR(2)的級聯(lián)表示：,狀態(tài)方程,24,二、ARMA模型及其特點,（4）ARMA(p,q) 模型特點,有限項級數(shù)之和 （自相關(guān)系數(shù)函數(shù)截尾）,無窮級數(shù)之和 （自相關(guān)系數(shù)函數(shù)拖尾）,ARMA(p,q)模型特點,無窮級數(shù)之和

12、,白噪聲 截尾 截尾,25,三、ARMA建模分析,傳統(tǒng)數(shù)字濾波器,處理確定性信號 仔細設(shè)計濾波器的結(jié)構(gòu)和參數(shù) 近乎完美的帶通和快速起伏性能 偏好于IIR濾波器,ARMA有色噪聲 成形濾波器,處理隨機信號 觀測時間序列樣本還是有限，隨機性和帶誤差 成形濾波器設(shè)計 “粗糙” 偏好于低階AR模型,（1）傳統(tǒng)數(shù)字濾波器與ARMA建模比較,26,三、ARMA建模分析,（2）ARMA建模流程,27,三、ARMA建模分析,（3）樣本統(tǒng)計特性,樣本均值,樣本自協(xié)方差函數(shù),設(shè) 是含N個數(shù)據(jù)的一個時間序列樣本，則,樣本自相關(guān)系數(shù)函數(shù),樣本偏自相關(guān)系數(shù)函數(shù) （Durbin遞推公式）,遞推計算順序,28,三、ARMA

13、建模分析,（4）樣本數(shù)據(jù)預(yù)處理（平穩(wěn)化處理）,測試數(shù)據(jù) 典型分解式,趨勢項,周期項,平穩(wěn)序列,趨勢項提取,最小二乘法,k階差分法（消除k次多項式趨勢）,周期項提取,功率譜密度（PSD）分析（檢驗是否存在周期項）,d步差分法（去除周期d）,平穩(wěn)性檢驗：特點不隨時間變化的均值、方差和自相關(guān)系數(shù)函數(shù)等統(tǒng)計量。有時序圖檢驗方法和相關(guān)圖檢驗方法。,隨機序列結(jié)果,29,三、ARMA建模分析,（5）ARMA模型識別（確定模型類別 ）,模型識別依據(jù)ARMA模型（偏）自相關(guān)系數(shù)截尾/拖尾特點。 當樣本數(shù)據(jù)量N充分大時，以有限序列樣本代替總體統(tǒng)計特性，結(jié)論：,根據(jù)學(xué)者Bartlett（1946）的研究,如果當kq

14、時 成立，則表示 是q步截尾的，判斷為MA(q)模型。,MA(2)序列的自相關(guān)系數(shù)函數(shù),30,三、ARMA建模分析,（5）ARMA模型識別（確定模型類別 ）,根據(jù)學(xué)者Quenouille（1949）的研究,如果當kp時 ，則表示 是p步截尾的，判斷為AR(p)模型。,AR(2)序列的偏自相關(guān)系數(shù)函數(shù),在模型識別時，一般原則是先進行AR(p)模型分析，力求近似為AR(p)模型，如果實在不行或者階數(shù)太高，再考慮使用MA(q)和ARMA(p,q)模型。在Matlab/financial工具箱的函數(shù)parcorr()、autocorr()、aryule()和tsmovavg()，可方便用于AR和MA時

15、間序列建模分析。,31,三、ARMA建模分析,（6）AR(p)模型參數(shù)估計,模型參數(shù),白噪聲方差,（7）后續(xù)工作,適用性檢驗； 模型優(yōu)化； 適應(yīng)性評估； ARMA模型轉(zhuǎn)化為狀態(tài)空間模型。,32,三、ARMA建模分析,（8）MA(q)模型參數(shù)估計,遞推公式：,33,三、ARMA建模分析,（9）幾種常見的低階ARMA模型的參數(shù)估計結(jié)果,34,三、ARMA建模分析,（10）舉例1,35,三、ARMA建模分析,（10）舉例2,ts = 0.001; t = (ts:ts:10); len = length(t); w = randn(length(t),1); w(len/2:end) = w(len

16、/2:end) + 3.0; b0,a0 = butter(6, 0.016); w0 = filter(b0,a0,w); b1,a1 = butter(6, 0.033); w1 = filter(b1,a1,w); b2,a2 = butter(6, 0.066); w2 = filter(b2,a2,w); b3,a3 = butter(6, 0.131); w3 = filter(b3,a3,w); b4,a4 = butter(6, 0.262); w4 = filter(b4,a4,w); figure, plot(t, w,w4,w3,w2,w1,w0,linewidth,1.5); xygo(w) legend(LP-NO,LP262,LP131,LP066,LP033,LP016),36,三、ARMA建模分析,（10）舉例2,37,三、ARMA建模分析,（10）舉例3,imu = load(stim300.txt)