时间序列分析方法及其进展*

赵 志 周 倩 张晋昕△

·综述·

时间序列分析方法及其进展*

赵 志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