时间序列分析与生存模型

Posted 2021-05-02 智能算法研究学习1688

时间序列分析与生存模型

时间序列分析

时间序列分析方法是经济领域研究的主要工具之一,它用合适的模型描述历史数据随时间变化的规律,并预测经济变量值.而ARMA模型是适用于任何序列的发展形态的一种高级预测方法,它描述时间序列的动态性和发展变化规律.文中通过ARMA模型分析时间序列的随机性和平稳性,以一种商品月度销售额具体分析,用sas软件检验模型的可行性,并预测应用.结果表明,模拟值和真实值接近,在实际应用中预测值的准确对于指导商家的战略决策起重要作用.

基本思想: 某些事件序列是依赖于时间的一族时间变量，构成该时序的单个序列值虽然具有不确定性，但整个序列的变化具有一定的规律性，可以用相应的数学模型近似描述。通过对该数学模型的分析和调研，能够更本质地认识时间序列结构和特征，达到最小方差意义下的最优预测。

 

1 时间序列与时间序列分析

在生产和科学研究中，对某一个或者一组变量 x(t)x(t) 进行观察测量，将在一系列时刻 t1,t2,⋯,tnt1,t2,⋯,tn 所得到的离散数字组成的序列集合，称之为时间序列。 
时间序列分析是根据系统观察得到的时间序列数据，通过曲线拟合和参数估计来建立数学模型的理论和方法。时间序列分析常用于国民宏观经济控制、市场潜力预测、气象预测、农作物害虫灾害预报等各个方面。

2 时间序列建模基本步骤

1.     获取被观测系统时间序列数据；

2.     对数据绘图，观测是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，化为平稳时间序列；

3.     经过第二步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF ，通过对自相关图和偏自相关图的分析，得到最佳的阶层 p 和阶数 q

4.     由以上得到的d、q、pd、q、p ，得到ARIMA模型。然后开始对得到的模型进行模型检验。

3 ARIMA实战解剖

原理大概清楚，实践却还是会有诸多问题。相比较R语言，Python在做时间序列分析的资料相对少很多。下面就通过Python语言详细解析后三个步骤的实现过程。 
文中使用到这些基础库： pandas,numpy,scipy,matplotlib,statsmodelspandas,numpy,scipy,matplotlib,statsmodels。对其调用如下

from __future__ importprint_function

import pandas as pd

import numpy as np

from scipy import  stats

importmatplotlib.pyplot as plt

import statsmodels.apias sm

fromstatsmodels.graphics.api import qqplot

3.1 获取数据

这里我们使用一个具有周期性的测试数据，进行分析。 
数据如下：

dta=[10930,10318,10595,10972,7706,6756,9092,10551,9722,10913,11151,8186,6422, 
6337,11649,11652,10310,12043,7937,6476,9662,9570,9981,9331,9449,6773,6304,9355, 
10477,10148,10395,11261,8713,7299,10424,10795,11069,11602,11427,9095,7707,10767, 
12136,12812,12006,12528,10329,7818,11719,11683,12603,11495,13670,11337,10232, 
13261,13230,15535,16837,19598,14823,11622,19391,18177,19994,14723,15694,13248, 
9543,12872,13101,15053,12619,13749,10228,9725,14729,12518,14564,15085,14722, 
11999,9390,13481,14795,15845,15271,14686,11054,10395]

dta=pd.Series(dta)

dta.index =pd.Index(sm.tsa.datetools.dates_from_range('2001','2100'))

dta.plot(figsize=(12,8))

3.2 时间序列的差分dd

ARIMA 模型对时间序列的要求是平稳型。因此，当你得到一个非平稳的时间序列时，首先要做的即是做时间序列的差分，直到得到一个平稳时间序列。如果你对时间序列做d次差分才能得到一个平稳序列，那么可以使用ARIMA(p,d,q)模型，其中d是差分次数。

fig =plt.figure(figsize=(12,8))

ax1=fig.add_subplot(111)

diff1 = dta.diff(1)

diff1.plot(ax=ax1)

 
一阶差分的时间序列的均值和方差已经基本平稳，不过我们还是可以比较一下二阶差分的效果

fig =plt.figure(figsize=(12,8))

ax2=fig.add_subplot(111)

diff2 = dta.diff(2)

diff2.plot(ax=ax2)

 
可以看出二阶差分后的时间序列与一阶差分相差不大，并且二者随着时间推移，时间序列的均值和方差保持不变。因此可以将差分次数d设置为1。 
其实还有针对平稳的检验，叫“ADF单位根平稳型检验”，以后再更。

3.3 合适的p,qp,q

现在我们已经得到一个平稳的时间序列，接来下就是选择合适的ARIMA模型，即ARIMA模型中合适的p,qp,q。 
第一步我们要先检查平稳时间序列的自相关图和偏自相关图。

dta= dta.diff(1)#我们已经知道要使用一阶差分的时间序列，之前判断差分的程序可以注释掉

fig =plt.figure(figsize=(12,8))

ax1=fig.add_subplot(211)

fig =sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)

ax2 =fig.add_subplot(212)

fig =sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax2)

其中lags 表示滞后的阶数，以上分别得到acf 图和pacf 图 
 
通过两图观察得到： 
* 自相关图显示滞后有三个阶超出了置信边界； 
* 偏相关图显示在滞后1至7阶（lags 1,2,…，7）时的偏自相关系数超出了置信边界，从lag 7之后偏自相关系数值缩小至0 
则有以下模型可以供选择： 
1. ARMA(0,1)模型：即自相关图在滞后1阶之后缩小为0，且偏自相关缩小至0，则是一个阶数q=1的移动平均模型； 
2. ARMA(7,0)模型：即偏自相关图在滞后7阶之后缩小为0，且自相关缩小至0，则是一个阶层p=3的自回归模型； 
3. ARMA(7,1)模型：即使得自相关和偏自相关都缩小至零。则是一个混合模型。 
4. …还可以有其他供选择的模型 
现在有以上这么多可供选择的模型，我们通常采用ARMA模型的AIC法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。不仅仅包括AIC准则，目前选择模型常用如下准则： 
* AIC=-2 ln(L) + 2 k 中文名字：赤池信息量 akaike information criterion 
* BIC=-2 ln(L) + ln(n)*k 中文名字：贝叶斯信息量 bayesian information criterion 
* HQ=-2 ln(L) + ln(ln(n))*k hannan-quinn criterion 
构造这些统计量所遵循的统计思想是一致的，就是在考虑拟合残差的同时，依自变量个数施加“惩罚”。但要注意的是，这些准则不能说明某一个模型的精确度，也即是说，对于三个模型Ａ，Ｂ，Ｃ，我们能够判断出Ｃ模型是最好的，但不能保证Ｃ模型能够很好地刻画数据，因为有可能三个模型都是糟糕的。

arma_mod20 =sm.tsa.ARMA(dta,(7,0)).fit()

print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)

arma_mod30 =sm.tsa.ARMA(dta,(0,1)).fit()

print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)

arma_mod40 =sm.tsa.ARMA(dta,(7,1)).fit()

print(arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic)

arma_mod50 =sm.tsa.ARMA(dta,(8,0)).fit()

print(arma_mod50.aic,arma_mod50.bic,arma_mod50.hqic)

 
可以看到ARMA(7,0)的aic，bic，hqic均最小，因此是最佳模型。

3.4 模型检验

在指数平滑模型下，观察ARIMA模型的残差是否是平均值为0且方差为常数的正态分布（服从零均值、方差不变的正态分布），同时也要观察连续残差是否（自）相关。

3.4.1 我们对ARMA(7,0)模型所产生的残差做自相关图

fig =plt.figure(figsize=(12,8))

ax1 =fig.add_subplot(211)

fig =sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)

ax2 = fig.add_subplot(212)

fig =sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)

3.4.2 做D-W检验

德宾-沃森（Durbin-Watson）检验。德宾-沃森检验,简称D-W检验，是目前检验自相关性最常用的方法，但它只使用于检验一阶自相关性。因为自相关系数ρ的值介于-1和1之间，所以 0≤DW≤４。并且DW＝O＝＞ρ＝１　　即存在正自相关性 
DW＝４＜＝＞ρ＝－１　即存在负自相关性 
DW＝２＜＝＞ρ＝０　　即不存在（一阶）自相关性 
因此，当DW值显著的接近于O或４时，则存在自相关性，而接近于２时，则不存在（一阶）自相关性。这样只要知道ＤＷ统计量的概率分布，在给定的显著水平下，根据临界值的位置就可以对原假设Ｈ０进行检验。

print(sm.stats.durbin_watson(arma_mod20.resid.values))

检验结果是2.02424743723，说明不存在自相关性。

3.4.3 观察是否符合正态分布

这里使用QQ图，它用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布。QQ图细节，下次再更。

resid =arma_mod20.resid#残差

fig =plt.figure(figsize=(12,8))

ax =fig.add_subplot(111)

fig = qqplot(resid,line='q', ax=ax, fit=True)

3.4.4 Ljung-Box检验

Ljung-Box test是对randomness的检验,或者说是对时间序列是否存在滞后相关的一种统计检验。对于滞后相关的检验，我们常常采用的方法还包括计算ACF和PCAF并观察其图像，但是无论是ACF还是PACF都仅仅考虑是否存在某一特定滞后阶数的相关。LB检验则是基于一系列滞后阶数，判断序列总体的相关性或者说随机性是否存在。 
时间序列中一个最基本的模型就是高斯白噪声序列。而对于ARIMA模型，其残差被假定为高斯白噪声序列，所以当我们用ARIMA模型去拟合数据时，拟合后我们要对残差的估计序列进行LB检验，判断其是否是高斯白噪声，如果不是，那么就说明ARIMA模型也许并不是一个适合样本的模型。

r,q,p = sm.tsa.acf(resid.values.squeeze(),qstat=True)

data =np.c_[range(1,41), r[1:], q, p]

table =pd.DataFrame(data, columns=['lag', "AC", "Q","Prob(>Q)"])

print(table.set_index('lag'))

 
检验的结果就是看最后一列前十二行的检验概率（一般观察滞后1~12阶），如果检验概率小于给定的显著性水平，比如0.05、0.10等就拒绝原假设，其原假设是相关系数为零。就结果来看，如果取显著性水平为0.05，那么相关系数与零没有显著差异，即为白噪声序列。

3.5 模型预测

模型确定之后，就可以开始进行预测了，我们对未来十年的数据进行预测。

predict_sunspots =arma_mod20.predict('2090', '2100', dynamic=True)

print(predict_sunspots)

fig, ax =plt.subplots(figsize=(12, 8))

ax =dta.ix['2001':].plot(ax=ax)

predict_sunspots.plot(ax=ax)

 
前面90个数据为测试数据，最后10个为预测数据；从图形来，预测结果较为合理。至此，本案例的时间序列分析也就结束了。

参考文献与推荐阅读

1.     statsmodels–statistics in python

2.     时间序列分析—(ARIMA模型)

3.     Arima预测模型（R语言）

4.     介绍QQplot

5.     LBQ检验

6.     经管之家

生存模型

Cox-regression 模型

logistic回归，与线性回归并成为两大回归，应用范围一点不亚于线性回归，甚至有青出于蓝之势。因为logistic回归太好用了，而且太有实际意义了。解释起来直接就可以说，如果具有某个危险因素，发病风险增加2.3倍，听起来多么地让人通俗易懂。线性回归相比之下其实际意义就弱了。logistic回归与线性回归恰好相反，因变量一定要是分类变量，不可能是连续变量。分类变量既可以是二分类，也可以是多分类，多分类中既可以是有序，也可以是无序。二分类logistic回归有时候根据研究目的又分为条件logistic回归和非条件logistic回归。条件logistic回归用于配对资料的分析，非条件logistic回归用于非配对资料的分析，也就是直接随机抽样的资料。无序多分类logistic回归有时候也成为多项logit模型，有序logistic回归有时也称为累积比数logit模型。

cox回归，cox回归的因变量就有些特殊，因为他的因变量必须同时有2个，一个代表状态，必须是分类变量，一个代表时间，应该是连续变量。只有同时具有这两个变量，才能用cox回归分析。cox回归主要用于生存资料的分析，生存资料至少有两个结局变量，一是死亡状态，是活着还是死亡？二是死亡时间，如果死亡，什么时间死亡？如果活着，从开始观察到结束时有多久了？所以有了这两个变量，就可以考虑用cox回归分析

 

 

在介绍Cox回归模型之前，先介绍几个有关的概念。

1.生存函数具有变量  的观察对象的生存时间  大于某时刻  的概率，

也就是相对于当前时间t，未来生存时间T的生存率。称为生存函数。生存函数  又称为累积生存率。

2. 死亡函数具有变量  的观察对象的生存时间  不大于某时刻 的概率，

称为死亡函数。死亡函数  的实际意义是当观察随访到  时刻的累积死亡率。

3. 死亡密度函数具有变量X的观察对象在某时刻t的瞬时死亡率，称为死亡密度函数。

4. 危险率(风险)函数具有变量X，且生存时间已达到 的观察对象在时刻 的瞬时死亡率，

危险率函数  实际上是一个条件瞬间死亡率 ^[2]  。

基本原理

生存分析的主要目的在于研究变量X与观察结果即生存函数(累积生存率)  

之间的关系。当  受很多因素影响，即 为向量时，传统的方法是考虑回归方程——即诸变量  对  的影响。但由于生存分析研究中的数据包含删失数据。且时间变量t通常不满足正态分布和方差齐性的要求，这就造成了用一般的回归方法研究上述关系的困难 ^[2]  。

Cox回归模型的基本形式D.R.Cox提出了Cox比例风险回归模型，它不是直接考察

  与X的关系，而是用 作为因变量，模型的基本形式为：

式中，  为自变量的偏回归系数，它是须从样本数据作出估计的参数；  是当X向量为0时， 的基准危险率，它是有待于从样本数据作出估计的量。公式(1)简称为Cox回归模型。

由于Cox回归模型对  未作任何假定，因此Cox回归模型在处理问题时具有较大的灵活性；另一方面，在许多情况下，我们只需估计出参数 (如因素分析等)，即使在

  未知的情况下，仍可估计出参数 。这就是说，Cox回归模型由于含有  ，因此它不是完全的参数模型，但仍可根据公式(1)作出参数  的估计，故Cox回归模型属于半参数模型。

公式(1)可以转化为：

Cox回归模型的假定

1. 比例风险假定 各危险因素的作用不随时间的变化而变化，即

  

不随时间的变化而变化。因此，公式(1)又称为比例风险率模型(PH Model)。这一假定是建立Cox回归模型的前提条件。

2．对数线性假定 模型中的协变量应与对数风险比呈线性关系，如公式(2)。

Cox回归模型中偏回归系数的意义，若 是非暴露组观察对象的各因素取值，  

是暴露组观察对象的各因素取值，由公式(3)就可以求出暴露组对非暴露组的相对危险度RR。

由公式(2)可见，模型中偏回归系数 的流行病学含义是在其他协变量不变的情况下，协变量  每增加一个测定单位时所引起的相对危险度的自然对数的改变量。即式中，  分别表示在不同情况下的取值。当协变量  分别取1和0时，其对应的 为

从公式(1)和公式(4)可以看出有如下关系：

若  ，则各  取值越大时， 的值越大，即  为危险因素。若

  ，则各  的取值对  的值没有影响，即 为无关因素。若

  ，则各  取值越大时 的值越小，即  为保护因素。

假设检验

Cox回归模型中的偏回归系数可以通过建立偏似然函数，利用Newton-Raphson迭代法求得。其他自变量不变的情况下，变量 每增加一个单位，相对危险度  的

 可信区间为：式中  为 的标准误。

对于回归模型的假设检验通常采用似然比检验、Wald检验和记分检验，其检验统计量均服从  分布，其自由度为模型中待检验的自变量个数。一般说来，Cox回归系数的估计和模型的假设检验计算量较大，通常需利用计算机来完成相应的计算 ^[2]  。

 

[Spark2生存分析Survival regression]

https://www.cnblogs.com/wwxbi/p/6150352.html

https://www.cnblogs.com/webRobot/p/6754160.html