時間序列分析方法及其進展*

趙 志 周 倩 張晉昕△

·綜述·

時間序列分析方法及其進展*

趙 志1周 倩2張晉昕1△

在醫學科研工作中，按某種(相等或不相等的)時間間隔對客觀事物進行動態觀察，由于隨機因素的影響，各次觀察的指標X1,X2,X3,…,Xi，…都是隨機變量，這種按時間順序排列的一系列隨機變量(或其觀測值)稱為時間序列。

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

時間序列的基本特征

1.平穩性

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

判斷一個隨機過程是否是嚴平穩的，這在實際應用中十分困難。因為它要求過程的各階矩都是時不變的(time-invariant)，即只與時間間隔有關，與初始時間無關，并且各階矩可以無窮大。而實際問題中，如果數據的方差無窮大，意味著變異無限大，那么就無法對數據進行建模來預測或控制。因此，當把條件適當放寬，只限定存在常均值、有限且時不變的二階矩，便得到了寬平穩過程(weakly stationary process)。若隨機過程的均值(一階距)和協方差(二階矩)存在，且滿足：

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

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

則稱{X(t),t∈T}為寬平穩過程，均值μ為常數，γ(τ)是X(t)的協方差函數，它與時間t無關。

判斷時間序列是否平穩是建模的前提，一般可通過繪制時序圖來判斷。若時序圖中的曲線圍繞一定均值呈現波動變化相似的情景，則可判定序列是平穩的。更準確地，可以通過單位根檢驗判斷平穩性。

2.可逆性

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

Xt=εt-θεt-1

(1)

若將εt用現在和過去的Xt表示，則(1)式變為

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

(2)

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

(3)

這時就稱序列{Xt}是可逆的。

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

(4)

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

時域分析方法的主流模型

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

1.ARIMA模型

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

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

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

(5)

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

ARIMA模型能夠為傳染病或相關醫療衛生設施方面的預測工作提供有效的指導。例如，2003年SARS爆發期間，Earnest和Chen等[6]利用ARIMA模型對新加坡的一家專科醫院每日病床占用數進行建模，所建立的ARMA(1,3)模型能夠較好地預測未來3天的病床需求數。

2.GARCH模型

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

Xt=φ1Xt-1+εt

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

條件異方差模型可看作是利用現在和過去的數據，基于AR模型對條件方差過程進行建模，通過該模型可以預測序列未來的波動大小。

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

3.非線性時間序列

條件異方差模型在序列與擾動項之間建立了一種非線性的關系，同樣，也可以在序列之間建立非線性的函數關系，如

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

(6)

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

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

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

頻域分析

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

定性時間序列

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

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

多元時間序列

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

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

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

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

研究熱點與展望

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

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

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

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

[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.時間序列分析及應用：R語言.北京:機械工業出版社,2011.

[4]Box GE,Jenkins GM,Reinsel GC.時間序列分析:預測與控制.北京:機械工業出版社,2011.

[5]方積乾,陸盈.現代醫學統計學.北京:人民衛生出版社,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.北京:科學出版社,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]周倩,張晉昕.時間序列周期性檢驗方法研究進展.中國衛生統計,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]倪延延,張晉昕.向量自回歸模型擬合與預測效果評價.中國衛生統計,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]詹思延,李立明,吳系科.流行病學進展.第12卷.北京:人民衛生出版社,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.

(責任編輯：郭海強)

國家自然科學基金(30872182)

1.中山大學公共衛生學院醫學統計與流行病學系(510080)

2.中山大學附屬第一醫院流行病學研究室

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