時(shí)間序列分析方法及其進(jìn)展*

趙 志 周 倩 張晉昕△

·綜述·

時(shí)間序列分析方法及其進(jìn)展*

趙 志1周 倩2張晉昕1△

在醫(yī)學(xué)科研工作中，按某種(相等或不相等的)時(shí)間間隔對(duì)客觀事物進(jìn)行動(dòng)態(tài)觀察，由于隨機(jī)因素的影響，各次觀察的指標(biāo)X1,X2,X3,…,Xi，…都是隨機(jī)變量，這種按時(shí)間順序排列的一系列隨機(jī)變量(或其觀測(cè)值)稱為時(shí)間序列。

例如流行病學(xué)家會(huì)關(guān)注在某地區(qū)觀察到的流感樣病例數(shù)、登革熱病例數(shù)隨時(shí)間的變化，時(shí)間序列模型可以考慮時(shí)間對(duì)病例數(shù)的影響，亦可納入不同的流行病學(xué)影響因素來預(yù)測(cè)病例數(shù)，以及探討疫情季節(jié)性特征。臨床中連續(xù)時(shí)間內(nèi)多次血壓的測(cè)量能夠評(píng)估藥物治療高血壓的效果；人腦的功能性核磁共振造影(fMRI)中能夠采集一系列時(shí)間序列，據(jù)以研究不同實(shí)驗(yàn)條件下大腦對(duì)于刺激的反應(yīng)[1]。生物學(xué)家也會(huì)對(duì)基因表達(dá)譜中隱含的一些重要模式規(guī)律感興趣，以期與相應(yīng)的表觀或疾病關(guān)聯(lián)[2]。時(shí)間序列分析方法是處理這些問題的有力工具。本文將對(duì)時(shí)間序列分析的經(jīng)典模型及最新研究進(jìn)展進(jìn)行綜述。

時(shí)間序列的基本特征

1.平穩(wěn)性

時(shí)間序列是一種特殊的隨機(jī)過程，根據(jù)觀測(cè)記錄對(duì)隨機(jī)過程的結(jié)構(gòu)進(jìn)行統(tǒng)計(jì)推斷時(shí)，通常需要作出某些假設(shè)，其中最重要的假設(shè)是平穩(wěn)性，即決定過程特性的統(tǒng)計(jì)規(guī)律不隨著時(shí)間的變化而變化。從一定意義上講，過程位于統(tǒng)計(jì)的平衡點(diǎn)上。特別地，對(duì)一切時(shí)滯k和時(shí)點(diǎn)t1,t2,…,tn，都有Xt1,Xt2,…,Xtn與Xt1-k,Xt2-k,…,Xtn-k的聯(lián)合分布相同，則稱過程{Xt}為嚴(yán)平穩(wěn)過程[3](strictly stationary process)。

判斷一個(gè)隨機(jī)過程是否是嚴(yán)平穩(wěn)的，這在實(shí)際應(yīng)用中十分困難。因?yàn)樗筮^程的各階矩都是時(shí)不變的(time-invariant)，即只與時(shí)間間隔有關(guān)，與初始時(shí)間無關(guān)，并且各階矩可以無窮大。而實(shí)際問題中，如果數(shù)據(jù)的方差無窮大，意味著變異無限大，那么就無法對(duì)數(shù)據(jù)進(jìn)行建模來預(yù)測(cè)或控制。因此，當(dāng)把條件適當(dāng)放寬，只限定存在常均值、有限且時(shí)不變的二階矩，便得到了寬平穩(wěn)過程(weakly stationary process)。若隨機(jī)過程的均值(一階距)和協(xié)方差(二階矩)存在，且滿足：

(1)E[X(t)]=μ,?t∈T;

(2)E[Xt+τ-μ][Xt-μ]=γ(τ),?t,t+τ∈T;

則稱{X(t),t∈T}為寬平穩(wěn)過程，均值μ為常數(shù)，γ(τ)是X(t)的協(xié)方差函數(shù)，它與時(shí)間t無關(guān)。

判斷時(shí)間序列是否平穩(wěn)是建模的前提，一般可通過繪制時(shí)序圖來判斷。若時(shí)序圖中的曲線圍繞一定均值呈現(xiàn)波動(dòng)變化相似的情景，則可判定序列是平穩(wěn)的。更準(zhǔn)確地，可以通過單位根檢驗(yàn)判斷平穩(wěn)性。

2.可逆性

可以假設(shè)隨機(jī)序列是隨機(jī)沖擊的線性組合，這種隨機(jī)序列就是用一般線性隨機(jī)模型來描述的[4]。例如對(duì)于隨機(jī)過程{Xt}有

Xt=εt-θεt-1

(1)

若將εt用現(xiàn)在和過去的Xt表示，則(1)式變?yōu)?/p>

εt=(1-θB)-1Xt

(2)

其中，B表示滯后算子，定義為BmXt=Xt-m。進(jìn)一步，若|θ|<1，(2)式可寫成無窮級(jí)數(shù)

(3)

這時(shí)就稱序列{Xt}是可逆的。

然而，如果|θ|>1，式(2)可寫成

(4)

即只能用過程的現(xiàn)在和將來值來刻畫當(dāng)前指標(biāo)值，這樣的模型不符合常理，被稱作不滿足可逆性條件。

時(shí)域分析方法的主流模型

常用的時(shí)間序列模型有很多，如移動(dòng)平均模型、條件異方差模型、非線性時(shí)間序列模型等。這些模型是從時(shí)域(time domain)角度對(duì)序列的將來值用過去值建模，要求相鄰時(shí)間點(diǎn)序列的相關(guān)性能夠被過去值很好地刻畫。

1.ARIMA模型

ARIMA模型(autoregressive integrated moving average model)又稱博克斯-詹金斯(Box-Jenkins)模型[5]，是由美國(guó)統(tǒng)計(jì)學(xué)家G.E.P.Box和G.M.Jenkins于1970年首次系統(tǒng)提出的，該模型有三種基本模式：自回歸模型、移動(dòng)平均模型與自回歸移動(dòng)平均模型。當(dāng)序列為平穩(wěn)序列，其形式為ARMA模型：

Xt=φ1Xt-1+φ2Xt-2+…+φpXt-p+εt-θ1εt-1-

θ2εt-2-…-θqεt-q

(5)

當(dāng)序列為非平穩(wěn)序列，但通過d次差分可使序列平穩(wěn)時(shí)，采用的模型稱作自回歸求和移動(dòng)平均模型，即ARIMA(p,d,q)模型。鑒于三種基本模式可視作ARIMA模型的特例，故ARIMA模型又常被作為這一族模型的總稱。對(duì)于判斷此模型是否平穩(wěn)，可通過DF(Dickey-Fuller法)或ADF(augmented Dickey-Fuller法)單位根檢驗(yàn)，并經(jīng)過考察θ(B)的根符合“|θ|<1”判定序列滿足可逆性。

ARIMA模型能夠?yàn)閭魅静』蛳嚓P(guān)醫(yī)療衛(wèi)生設(shè)施方面的預(yù)測(cè)工作提供有效的指導(dǎo)。例如，2003年SARS爆發(fā)期間，Earnest和Chen等[6]利用ARIMA模型對(duì)新加坡的一家專科醫(yī)院每日病床占用數(shù)進(jìn)行建模，所建立的ARMA(1,3)模型能夠較好地預(yù)測(cè)未來3天的病床需求數(shù)。

2.GARCH模型

經(jīng)典的ARIMA模型假定過程的當(dāng)前值依賴于過去值、當(dāng)前擾動(dòng)和過去擾動(dòng)，而擾動(dòng)項(xiàng)εt為白噪聲。但現(xiàn)實(shí)中擾動(dòng)不一定是理想的白噪聲，當(dāng)前的擾動(dòng)可能會(huì)受到前一期甚至前幾期序列蘊(yùn)涵的信息影響，ARCH模型[7](autoregressive conditional heteroskedasticity model)由此而來。ARCH(q) 定義為：

Xt=φ1Xt-1+εt

這是對(duì)ARCH模型的推廣，稱之為GARCH(p,q)模型(generalizedautoregressiveconditionalheteroskedasticitymodel)。

條件異方差模型可看作是利用現(xiàn)在和過去的數(shù)據(jù)，基于AR模型對(duì)條件方差過程進(jìn)行建模，通過該模型可以預(yù)測(cè)序列未來的波動(dòng)大小。

Moran和Solomon[8]使用GARCH模型對(duì)澳大利亞和新西蘭1995-2009年間重癥監(jiān)護(hù)室(ICU)成人月死亡數(shù)據(jù)進(jìn)行建模，從而對(duì)其實(shí)現(xiàn)有效的監(jiān)控。Modarres和Ouarda等[9]對(duì)加拿大蒙特利爾地區(qū)的氣候條件和髖部骨折率建立了多元GARCH模型，文章指出雪的深度、氣溫、氣壓以及白晝長(zhǎng)度與髖部骨折率存在較大的關(guān)聯(lián)。

3.非線性時(shí)間序列

條件異方差模型在序列與擾動(dòng)項(xiàng)之間建立了一種非線性的關(guān)系，同樣，也可以在序列之間建立非線性的函數(shù)關(guān)系，如

Xt=f1(Xt-1)+…+fp(Xt-p)+εt

(6)

其中，fi(Xt-i)可以是已知的非線性函數(shù)，i=1,2,…p，可取fi(Xt-i)=1/Xt-i。當(dāng)fi(Xt-i)為未知的函數(shù)時(shí)，模型(6)稱為可加自回歸模型(additiveautoregressivemodel)[10]，表示為Xt～AAP(p)，這類非線性關(guān)系未知的模型統(tǒng)稱為非參數(shù)時(shí)間序列模型。

另一種使用較為廣泛的非線性時(shí)間序列模型是門限自回歸模型(thresholdautoregressivemodel)。例如流感所引起的病死人數(shù)的增加速度往往比減少速度慢，并且當(dāng)病死人數(shù)達(dá)某一值時(shí)，過程會(huì)從一個(gè)方程(模式)突然切換到另一個(gè)方程。此時(shí)模型可設(shè)定為

其中ε1t和ε2t是白噪聲，d表示模式改變的時(shí)刻，稱之為延遲參數(shù)，c為臨界值。

頻域分析

頻域分析的思想認(rèn)為所有時(shí)間序列都可以看作是不同周期成分疊加的結(jié)果。周期成分廣泛存在于生物醫(yī)學(xué)的時(shí)間序列數(shù)據(jù)中，大到人群疾病流行強(qiáng)度的周期波動(dòng)，小到細(xì)胞新陳代謝的周期生長(zhǎng)，如生物醫(yī)學(xué)信號(hào)處理中的心電圖、腦電圖、醫(yī)院月度門診量等[11-12]。具有明顯周期成分的時(shí)間序列一般會(huì)在時(shí)域圖中顯示出周期性質(zhì)來，但更豐富的周期信息常常蘊(yùn)含于序列內(nèi)部，通過肉眼對(duì)時(shí)域圖的判讀難以發(fā)現(xiàn)并量化它，需要通過特定的方法將這種周期信息提取出來。早在 1929 年統(tǒng)計(jì)學(xué)家 R.A.Fisher就對(duì)時(shí)間序列周期性檢驗(yàn)進(jìn)行過研究，他運(yùn)用傅立葉變換獲得時(shí)間序列周期圖并提出基于周期圖法的 Fisherg統(tǒng)計(jì)量用于檢測(cè)峰值[13]。傅立葉變換是頻域分析的基本工具。之后統(tǒng)計(jì)學(xué)家對(duì)長(zhǎng)短不同、背景噪聲大小不同的周期性檢驗(yàn)方法不斷予以改進(jìn)[14]。

定性時(shí)間序列

定性時(shí)間序列(categorical time series) 又稱分類時(shí)間序列，是指每個(gè)時(shí)間點(diǎn)觀測(cè)值的取值范圍為有限狀態(tài)空間C={C1，C2，…，Cn}的時(shí)間序列，其取值只能表示狀態(tài)或者類別。定性時(shí)間序列廣泛存在于各個(gè)領(lǐng)域，如生物醫(yī)學(xué)、行為學(xué)、流行病學(xué)、遺傳學(xué)等[15]。DNA序列、睡眠狀態(tài)監(jiān)測(cè)就是由有限狀態(tài)空間形成的定性時(shí)間序列[16-17]。

定性時(shí)間序列的分析方法早期主要集中在時(shí)域分析部分，常見的有馬爾科夫鏈分析法、連接函數(shù)法等[17]。隨著定性時(shí)間序列的周期性研究的應(yīng)用開展，頻域分析方法也開始出現(xiàn)并得到改進(jìn)。包括基于變換的功率譜分析法、基于序列啞元化的譜封分析法(spectral envelope analysis)[18]、基于非平穩(wěn)序列的小波分析法[19]等。其中，譜封分析法適用于各種類型時(shí)間序列的周期性分析，它是定性時(shí)間序列周期性分析中檢驗(yàn)效果較好的一種方法。

多元時(shí)間序列

研究某一地區(qū)由于上呼吸道感染、肺炎、哮喘所引起的就診人數(shù)時(shí)間序列，如果出現(xiàn)異常的天氣條件，各類疾病的門診量時(shí)間序列就會(huì)同時(shí)發(fā)生改變，建模中不能忽視可能存在的關(guān)聯(lián)[20]。

關(guān)于一元時(shí)間序列的大部分基本理論和方法都可以推廣到多元時(shí)間序列，最常用的多元時(shí)間序列模型是向量自回歸模型(vector autoregressive model,VAR)。考慮p階的多元時(shí)間序列Xt滿足VAR(p) 模型[21-22]，則

其中，Xt可以是上述例子中門診量的3維列向量，φ0是3維常數(shù)列向量，而φi是3×3維的系數(shù)矩陣，εt是3維的白噪聲過程。

多元情況下也存在向量ARMA模型，但其參數(shù)過多，往往存在模型的可識(shí)別性問題。針對(duì)非平穩(wěn)多元時(shí)間序列，一般需要檢驗(yàn)序列的線性協(xié)整(cointegration)并建立誤差修正模型(error correction model)。另外，多元時(shí)間序列中還存在一種重要的模型結(jié)構(gòu)關(guān)系——格蘭杰因果關(guān)系(Granger causality)，它對(duì)于研究不同序列之間的繼起關(guān)系具有重要價(jià)值。

研究熱點(diǎn)與展望

時(shí)間序列模型的研究是為了不斷適應(yīng)實(shí)踐中遇到的新問題，建立更合理有效的模型。近年來在多元時(shí)間序列方面研究較多，該領(lǐng)域主要的發(fā)展動(dòng)向是：(1)由一元非線性模型轉(zhuǎn)向多元非線性模型；(2)非平穩(wěn)的多元時(shí)間序列與小波分析等工具的結(jié)合；(3)高維多元時(shí)間序列的降維問題；(4)時(shí)間序列與統(tǒng)計(jì)過程控制相結(jié)合，可用于癥狀監(jiān)測(cè)，對(duì)于檢出疾病暴發(fā)等異常變動(dòng)具有重要價(jià)值。

一般的研究對(duì)象往往處在紛繁復(fù)雜的系統(tǒng)中。例如傳染病發(fā)病人數(shù)的變化很可能與氣候因素、醫(yī)療衛(wèi)生條件等密切相關(guān)，這些序列各自都可能呈現(xiàn)非平穩(wěn)且非線性的特性，普通的VAR模型或非參數(shù)時(shí)間序列模型也不再適用。此時(shí)，建立起合理有效的多元非線性模型引起研究者的廣泛關(guān)注，主要研究角度是從非參數(shù)協(xié)整和半?yún)?shù)協(xié)整[23-26]對(duì)多元時(shí)間序列的非平穩(wěn)性進(jìn)行修正。它們對(duì)于非平穩(wěn)多元時(shí)間序列的推廣應(yīng)用具有重要意義。

針對(duì)非平穩(wěn)的多元時(shí)間序列，另一研究方向是從頻域角度出發(fā)與小波分析相結(jié)合[27]。小波分析不僅對(duì)非平穩(wěn)性具有很好的容忍性，還能夠同時(shí)結(jié)合時(shí)域分析與頻域分析的優(yōu)勢(shì)，使過程的信息得到更全面的利用。另外，還可以將小波變換與信息論相結(jié)合，利用小波熵來提取序列蘊(yùn)涵的信息[28-29]。

醫(yī)學(xué)研究中還存在一種具有很多變量或特征的時(shí)間序列資料，如fMRI影像數(shù)據(jù)中，每個(gè)像素就是一個(gè)變量或特征。對(duì)于高維的多元時(shí)間序列，為了便于統(tǒng)計(jì)分析以及實(shí)際意義的解釋，往往需要進(jìn)行降維處理，主成分分析[30]、因子模型[31-32]、LASSO[33]等降維方法在高維多元時(shí)間序列中使用較為廣泛。例如對(duì)于fMRI數(shù)據(jù)使用主成分分析和格蘭杰因果分析[34]，可以探討不同腦區(qū)之間的功能性聯(lián)系，為臨床診斷和治療提供參考。公共衛(wèi)生領(lǐng)域?qū)W者可以利用時(shí)間序列模型方法進(jìn)行疾病的長(zhǎng)期趨勢(shì)預(yù)測(cè)與疫情監(jiān)測(cè)，歸納出疾病演變的規(guī)律；借助統(tǒng)計(jì)過程控制方法，提高監(jiān)測(cè)的準(zhǔn)確性[35]。目前的熱點(diǎn)是變點(diǎn)分析(change-point detection)，找出樣本序列的分布或特征突然發(fā)生變化的某個(gè)或某些時(shí)間點(diǎn)，來提示可能存在的疫情爆發(fā)并加以預(yù)警[36-38]；或評(píng)價(jià)公共衛(wèi)生干預(yù)手段是否可以起效，估計(jì)起效的遲滯時(shí)間長(zhǎng)度。時(shí)間序列分析技術(shù)為前瞻性地了解疫情、主動(dòng)地控制疾病流行提供了有效支持。

[1]Shumway RH,Stoffer DS.Time Series Analysis and Its Application With R Example.2nd Edition.北京:世界圖書出版社,2011.

[2]Vivian T,Liew WC,Yan H.Periodicity analysis of DNA microarray gene expression time series profiles in mouse segmentation clock data.Statistics and Its Interface,2010,3(3):413-418.

[3]Cryer JD,Chan KS.時(shí)間序列分析及應(yīng)用：R語言.北京:機(jī)械工業(yè)出版社,2011.

[4]Box GE,Jenkins GM,Reinsel GC.時(shí)間序列分析:預(yù)測(cè)與控制.北京:機(jī)械工業(yè)出版社,2011.

[5]方積乾,陸盈.現(xiàn)代醫(yī)學(xué)統(tǒng)計(jì)學(xué).北京:人民衛(wèi)生出版社,2002.

[6]Earnest A,Chen MI,Ng D,et al.Using autoregressive integrated moving average(ARIMA) models to predict and monitor the number of beds occupied during a SARS outbreak in a tertiary hospital in Singapore.BMC Health Serv Res,2005,5(36):1-8.

[7]Enders W.Applied Econometric Time Series.New Jersey:John Wiley & Sons,2004.

[8]Moran JL,Solomon PJ.Statistical process control of mortality series in the Australian and New Zealand Intensive Care Society(ANZICS) adult patient database:implications of the data generating process.BMC Medical Res Methodol,2013,13(66):1-12.

[9]Modarres R,Ouarda TBMJ,Vanasse A,et al.Modeling climate effects on hip fracture rate by the multivariate GARCH model in Montreal region,Canada.Int J Biometeorol,2014,58(5):921-930.

[10]Fan JQ,Yao QW.Nonlinear Time Series:Nonparametric and Parametric Methods.北京:科學(xué)出版社,2006.

[11]Stevenson NJ,O’Toole JM,Rankine LJ,et al.A nonparametric feature for neonatal EEG seizure detection based on a representation of pseudo-periodicity.Medical Eng & Phys,2011,34(4):437-446.

[12]周倩,張晉昕.時(shí)間序列周期性檢驗(yàn)方法研究進(jìn)展.中國(guó)衛(wèi)生統(tǒng)計(jì),2013,30(3):445-447.

[13]Fisher RA.Tests of significance in harmonic analysis.Proc Math Phys & Eng Sci,1929,125(796):54-59.

[14]Krafty RT,Xiong S,Stoffer DS,et al.Enveloping spectral surfaces:covariate dependent spectral analysis of categorical time series.J Time Series Anal,2011,33(5):797-806.

[15]Freeman WJ,Viana DPG.Relation of olfactory EEG to behavior:time series analysis.Behav Neurosci,1986,100(5):753-763.

[16]Stoffer DS,Tyler DE,Mcdougall AJ.Spectral analysis for categorical time series:scaling and the spectral envelope.Biometrika,1993,83(3):611-622.

[17]Mcgee M,Ensor K.Tests for harmonic components in the spectra of categorical time series.J Time Series Anal,1998,19(3):309-323.

[18]Stoffer DS,Tyler DE.Matching sequences:Cross-spectral analysis of categorical time series.Biometrika,1998,85(1):201-213.

[19]Rasheed F,Alshalalfa M,Alhajj R.Adaptive machine learning technique for periodicity detection in biological sequences.Int J Neural Syst,2009,19(1):11-24.

[20]Diggle P.Time Series:An Biostatistical Introduction.New York:Oxford University Press,1990.

[21]Tsay RS.Multivariate Time Series Analysis:With R and Financial Applications.New Jersey:John Wiley & Sons,2014.

[22]倪延延,張晉昕.向量自回歸模型擬合與預(yù)測(cè)效果評(píng)價(jià).中國(guó)衛(wèi)生統(tǒng)計(jì),2014,31(1):53-56.

[23]Gu JP,Liang ZW.Testing cointegration relationship in a semiparametric varying coefficient model.J Econometrics,2014,178(1):57-70.

[24]Boswijk HP,Lucas A.Semi-nonparametric cointegration testing.J Econometrics,2002,108(2):253-280.

[25]Wang QY,Phillips CB.Structural Nonparametric Cointegrating Regression.Econometrics ,2009,77(6):1901-1948.

[26]Kasparis I,Phillips CB.Dynamic misspecification in nonparametric cointegrating regression.J Econometrics,2012,168(2):270-284.

[27]Maharaj EA,Alonso AM.Discriminant analysis of multivariate time series:Application to diagnosis based on ECG signals.Comput Stat & Data Anal,2014,70:67-87.

[28]Chen JK,Li GQ.Tsallis Wavelet Entropy and Its Application in Power Signal Analysis.Entropy,2014,16(6):3009-3025.

[29]Fraiwan L,Lweesy K,Khasawneh N,et al.Classification of Sleep Stages Using Multi-wavelet Time Frequency Entropy and LDA.Methods Inf Med,2010,49(3):230-237.

[30]Li H.Asynchronism-based principal component analysis for time series data mining.Expert Syst Appl,2014,41(6):2842-2850.

[31]Lam C,Yao QW.Factor modelling for high-dimensional time series:inference for the number of factors.Annals Stat,2012,40(2):694-726.

[32]Pan JZ,Yao QW.Modelling multiple time series via common factors.Biometrika,2008,95(2):365-379.

[33]Hsu NJ,Hung HL,Chang YM.Subset selection for vector autoregressive processes using Lasso.Comput Stat & Data Anal,2008,52(7):3645-3657.

[34]Zhou ZY,Chen YH,Ding MZ,et al.Analyzing Brain Networks With PCA and Conditional Granger Causality.Hum Brain Mapp ,2009,30(7):2197-2206.

[35]詹思延,李立明,吳系科.流行病學(xué)進(jìn)展.第12卷.北京:人民衛(wèi)生出版社,2010.

[36]Chen MP,Shang N,Winston CA,et al.A Bayesian analysis of the 2009 decline in tuberculosis morbidity in the United States.Stat Med,2012,31(27):3278-3284.

[37]Kass-Hout TA,Xu ZH,McMurray P,et al.Application of change point analysis to daily influenza-like illness emergency department visits.J Am Med Inform Assoc,2012,19(6):1075-1081.

[38]Riffenburgh RH,Cummins KM.A simple and general change-point identifier.Stat Med,2006,25(6):1067-1077.

(責(zé)任編輯：郭海強(qiáng))

國(guó)家自然科學(xué)基金(30872182)

1.中山大學(xué)公共衛(wèi)生學(xué)院醫(yī)學(xué)統(tǒng)計(jì)與流行病學(xué)系(510080)

2.中山大學(xué)附屬第一醫(yī)院流行病學(xué)研究室

△通信作者：張晉昕，Email：zhjinx@mail.sysu.edu.cn