精华内容
下载资源
问答
  • 时间序列分析 sas实例

    2008-09-25 19:52:40
    里面有三个时间序列最常用的三个实例,在建模中有用..
  • SAS时间序列分析

    2014-04-03 10:06:11
    SAS时间序列分析中的运用,很详细,武汉大学出版
  • sas时间序列分析

    2011-07-21 22:58:04
    sas时间序列分析,时间序列是把反映现象发展水平的统计指标数值,按照时间先后顺序排列起来所形成的一组统计数字序列。时间序列又称动态数列或时间数列。时间序列分析就是利用这组数列,应用数理统计方法加以处理,...
  • SAS学习系列38.时间序列分析Ⅱ非平稳时间序列的确定性分析.doc
  • SAS系统讲义-平稳时间序列分析
  • 本篇以美国1980年-2015年月度失业率为对象,进行一个更为完善的、有季节效应的非平稳时间序列分析流程。 理论支持: 拿到非平稳时间序列,首先进行的就是差分法消除趋势性,然后根据情况判断拟合季节加法模型或乘法...

    前言:前一篇介绍了对平稳时间序列的分析方法和流程,在没有考虑季节效应的情况下,模型建立的并不成功。本篇以美国1980年-2015年月度失业率为对象,进行一个更为完善的、有季节效应的非平稳时间序列分析流程。
    理论支持:
    拿到非平稳时间序列,首先进行的就是差分法消除趋势性,然后根据情况判断拟合季节加法模型或乘法模型,最后进行模型检验。常用的模型有两种:ARIMA和因素分解模型。

    1. ARIMA(加、乘法)模型,本篇分析采用。
    2. 因素分解模型:序列收三个因素影响:长期趋势,季节效应,随机波动。剔除前两者后留下随机波动。使用方法为简单中心移动平均法提取趋势效应,加法乘法提取季节效应。数学运算涉及较多,步骤繁琐,此处不多赘述。

    步骤:
    一. 数据录入,做出时序图

    data unem;
    input rate@@;
    time=intnx('month','01jan1980'd, _n_ -1);
    format time monyy.;
    cards;
    …………
    ;
    run;
    proc gplot data=unem;
    plot rate*time/ vaxis=2 to 13 by 0.5;
    symbol v=none c=red i=join;
    run;
    

    在这里插入图片描述
    初步观察发现,序列有明显的趋势性和周期性。
    二. 对原序列进行1阶12步差分,做出差分后时序图并进行ADF、白噪声检验

    data unem;
    input rate@@;
    difrate=dif12(dif(rate));
    time=intnx('month','01jan1980'd, _n_ -1);
    format time monyy.;
    cards;
    ………
    proc gplot;
    plot rate*time difrate*time;
    symbol v=none c=red i=join;
    proc arima data=unem;
    identify var=difrate stationarity=(adf);
    run;
    

    在这里插入图片描述
    在这里插入图片描述
    由ADF和白噪声检验得知,差分后的序列为平稳非白噪声序列。
    三. 模型定阶,选择加法还是乘法模型
    通常来说,加法模型适用于序列的季节效应,趋势效应和随机波动彼此之间很容易分开。但实践中更常见的情况是,序列的季节效应、长期趋势效应和随机波动之间存在复杂的交互影响关系,简单的加法模型不足以充分提取相关关系,此时应当使用乘法模型。
    在这里插入图片描述
    此例中,自相关图显示延迟12阶自相关系数显著大于二倍标准差,说明差分后的序列具有显著的季节效应,此外延迟2、3、5阶也大于二倍标准差,意味着差分后序列还具有短期相关性。在通过加法模型无法充分提取的情况下,我们尝试使用乘法模型,构造原理如下:

    1. 用低阶ARMA(p,q)模型提取序列短期相关性。
    2. 当序列有季节效应,季节效应本身又具有相关性的时候,用以周期步长为单位的ARMA(P,Q)S 提取季节相关性。
    3. 拟合模型实际上为ARMA(p,q)和ARMA(P,Q)S的乘积,记为ARMA(p,d,q) * ARMA(P,D,Q)S

    定阶分析过程:
    首先观察1阶12步差分后系列12阶以内的自相关系数和偏自相关系数的特征,以确定短期相关模型。自相关图和偏自相关图显示12阶以内均不截尾,考虑使用ARMA(1,1)模型提取短期自相关信息。
    再考虑季节相关性,检查延迟12阶、24阶等以周期长度(12)为单位的自相关系数和偏自相关系数的特征。自相关图12阶自相关系数显著,但是24阶系数在二倍标准差范围之内。而偏自相关图显示12阶和24阶都显著,且12阶到24阶之间没有超出二倍标准差范围。因而可以认为季节自相关特征是自相关系数截尾、偏自相关系数拖尾。此时用以12步为周期的ARMA(0,1)12模型提取季节自相关信息。
    综上,最终要拟合的乘法模型为ARIMA(1,1,1)*(0,1,1)12。

    四. 参数估计
    模型定阶之后,参数估计就简单很多了。

    /*注释:乘法模型常见格式为ARIMA(p,1,q)*(m,1,n)s, sas拟合命令为:
    identify var=x(1,s); 
    estimate p=(p)(ms)q=(q)(ns)
    第一句命令要求系统进行一阶s步差分,这里我们前面差分过了因此不用写。
    第二局命令要求系统对差分后序列拟合非季节效应模型ARMA(p,q)与季节效应模型ARMA(m,n)s。合并起来就完成了乘法模型的拟合。*/
    estimate p=(1) q=(1)(12) noint;
    run;
    

    五. 模型检验

    1. 可以发现参数显著性检验通过
      在这里插入图片描述
    2. 残差自相关检验中,P值在小于等于延迟24阶时都不大于5%显著性水平,因而不能认为残差序列已经是白噪声序列,即该模型拟合不显著,结合下面残差相关性检验结果进行分析原因。
      在这里插入图片描述

    六、问题原因及处理:
    传统的纯随机性检验都是借助LB检验统计量进行的,而LB检验统计量是在序列满足方差齐性的假定下构造的。当序列存在异方差属性时,LB统计量不在近似服从卡方分布。所以在条件异方差存在的场合,白噪声检验结果不再准确。通常现象就是残差序列的相关系数很小,近似白噪声序列,但是LB检验结果P值很小。因此,在异方差可能存在的场合,如果自相关系数很小(<0.2),则可以认为残差序列近似为白噪声序列。
    残差自回归性检验和异方差性检验:
    使用model语句,让系统建立序列关于时间的线性回归模型,检验残差序列5阶延迟的自相关性并输出DW检验的p值,同时对残差序列进行异方差检验:

    proc autoreg data=unem;
    model difrate=time/nlag=5 dwprob archtest; 
    run;
    

    DW检验结果,R方很小,模型不显著,证明残差序列不具有显著的自相关性,这与我们之前分析的结果相同(相关系数<0.2,残差序列可以认为是白噪声)。
    在这里插入图片描述
    参数估计也表明回归模型不显著,因而总结为残差序列不具有自相关性
    在这里插入图片描述

    异方差性(ARCH)检验:
    在这里插入图片描述
    图中Q统计量和LM统计量的P值均小于0.05显著性水平,因而可以认定该序列方差非齐。
    剩下的步骤为构建GARCH或ARCH模型以及检验,最后预测,内容较多,放在下一篇讨论。

    展开全文
  • 时间序列分析步骤及sas代码

    万次阅读 多人点赞 2015-07-13 20:59:50
    • 1. 时间序列数据的预处理:平稳性... 平稳时间序列数据分析 • 3. 非平稳时间序列数据分析    一. 时间序列数据的预处理 (1) 平稳性检验   程序: data a; input  x @@; time=_n_; cards; 17.4 20 17.9

    •         1. 时间序列数据的预处理:平稳性检验、纯随机性检验

    •         2. 平稳时间序列数据分析

    •         3. 非平稳时间序列数据分析 

     

    一.           时间序列数据的预处理

    (1)  平稳性检验

     

    程序

    data a;input x @@;

    time=_n_;

    cards;

    17.4   20  17.9   14.1   12.9   13.6   14.9

    18.2   16.8   16.9   17.6   18.9   19.3   17.7

    15.6   15  16.8   18.5   19.5   19  17.5

    14.5   14.3   16.2   16.5   16  17.5   19.6

    20  19.3   17  15.4   15.4   16.9   18.2

    17.7   17  16.8   15.2   14.5   16  17.1

    ;

    procgplot;

    plot x*time;

    symbolv=diamondi=joinc=blue;

    procarimadata=a;

    identifyvar=x;

    run;

     

    程序运行结果显示为:

     

    1.       时序图  2.自相关系数图

    结合时序图(在一条平行于x轴的线上起伏变化)和自相关系数图(在7阶以内收敛于0,两倍标准误差)

     

    稳定时间序列:

    (2)纯随机性检验(白噪声检验)

    原假设为白噪声序列,需拒绝原假设,则所有p值<0.05,则满足随机性。

     

    (3)模型识别

    将对应语句修改为:identify var=factory nlag=18 minic p=(0:5) q=(0:5);

    运行结果显示最小信息值及其对应模型为RA()或MA()。

     

    (4)模型拟合

    estimate p=a q=b method=ml;(a,b为对应p,q值。必要时需采用疏系数,如:p=(2,4,5))

     

    (5)预测并作出拟合图

    forecast id=time lead=5 out=results;/*lead预测期数,id指定身份变量,out预测结果存入某数据集*/

     

    proc gplot data=results;

    plot factory*time=1  forecast*time=2  l95*time=3u95*time=3/overlay;

    symbol1 v=star i=join c=black;

    symbol2 v=none i=join c=red;

    symbol3 v=none i=join c=green;

    run;

     

     

    非稳定时间序列

    一.       差分运算

    data a;input x@@;

    dif1=dif(x);

    time=_n_;

    cards;

    1.05   -0.84  -1.42  0.2 2.81   6.72   5.4 4.38

    5.52   4.46   2.89   -0.43  -4.86  -8.54  -11.54 -16.22

    -19.41 -21.61 -22.51 -23.51 -24.49 -25.54 -24.06 -23.44

    -23.41 -24.17 -21.58 -19 -14.14 -12.69 -9.48  -10.29

    -9.88  -8.33  -4.67  -2.97  -2.91  -1.86  -1.91  -0.8

    ;

    procgplot;plot x*time dif1*time;symbolc=blacki=joinv=star;

    procarima;

    identifyvar=x(1)nlag=22minicp=(0:5)q=(0:5);

    estimatep=1 noint method=ml;

    forecastlead=3id=timeout=out;/*forecastΪԤ²â¹Ø¼ü´Ê£¬leadΪԤ²â½×Êý*/

    procgplotdata=out;  plot x*time=1 forecast*time=2 l95*time=3 u95*time=3/overlay;

    symbol1c=blacki=nonev=star; symbol2c=redi=joinv=none;symbol3c=bluei=joinv=none;

    run;

     

    二.简单季节模型和乘积季节模型

     

    一阶差分12步+乘积季节模型

     

    data a;input x@@;

    dif1_12=dif12(dif(x));

    time=intnx('quarter','1jan1948'd,_n_-1);

    format time year4.;

    cards;

    /*数据省略*/

    ;

    proc gplot; plot x*time dif1_12*time; symbol c=blacki=join v=star;

    proc arima; identify var=x(1,12);

    estimate p=1 q=(1)(12) noint;

    forecast lead=0 id=time out=out;

    proc gplot data=out;

    plot x*time=1 forecast*time=2 /overlay;

    symbol1 c=blacki=none v=star; symbol2 c=redi=join v=none;

    run;

    展开全文
  • 为分析火控精度时间序列,采用SAS 语言模块,分别从时间序列数据的预处理、时间序列建模、时间序列的预报方面,讨论了时间序列分析问题,并结合火控精度实验中的数据实例探讨了SAS 的优越性和实用性。结果表明,利用SAS ...
  • 介绍SAS软件中相关系数,ARMA模型,Garch模型等常用代码的用法
  • SAS系统讲义-时间序列分析1-3实例
  • 一、什么是时间序列时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的...

    一、什么是时间序列

    时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。

    在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里不考虑含外生变量的时间序列)。

    环境配置

    python作为科学计算的利器,当然也有相关分析的包:statsmodels中tsa模块,当然这个包和SAS、R是比不了,但是python有另一个神器:pandas!pandas在时间序列上的应用,能简化我们很多的工作。这两个包pip就能安装。

    数据准备

    许多时间序列分析一样,本文同样使用航空乘客数据(AirPassengers.csv)作为样例。下载链接。

    用pandas操作时间序列

    #-*- coding:utf-8 -*-

    importnumpy as npimportpandas as pdfrom datetime importdatetimeimportmatplotlib.pylab as plt#读取数据,pd.read_csv默认生成DataFrame对象,需将其转换成Series对象

    df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='Month')

    df.index= pd.to_datetime(df.index) #将字符串索引转换成时间索引

    ts = df['Passengers'] #生成pd.Series对象#查看数据格式

    print(ts.head())print(ts.head().index)

    1365470-20200130213601204-441660966.png

    不知道为啥,时间隔1000,不过没什么影响。

    查看某日的值既可以使用字符串作为索引,又可以直接使用时间对象作为索引

    ts['2049-01-01']

    ts[datetime(2049,1,1)]

    两者的返回值都是第一个序列值:112

    如果要查看某一年的数据,pandas也能非常方便的实现

    ts['2049']

    切片操作:

    ts['2049-1' : '2049-6']

    注意时间索引的切片操作起点和尾部都是包含的,这点与数值索引有所不同

    pandas还有很多方便的时间序列函数,在后面的实际应用中在进行说明。

    二、时间序列分析

    1. 基本模型

    自回归移动平均模型(ARMA(p,q))是时间序列中最为重要的模型之一,它主要由两部分组成: AR代表p阶自回归过程,MA代表q阶移动平均过程,其公式如下:

    1365470-20200130213931237-1472383726.png

    1365470-20200130213936859-1941353511.png

    依据模型的形式、特性及自相关和偏自相关函数的特征,总结如下:

    1365470-20200130213959360-1759676234.png

    在时间序列中,ARIMA模型是在ARMA模型的基础上多了差分的操作。

    2. 平稳性检验

    我们知道序列平稳性是进行时间序列分析的前提条件,很多人都会有疑问,为什么要满足平稳性的要求呢?在大数定理和中心定理中要求样本同分布(这里同分布等价于时间序列中的平稳性),而我们的建模过程中有很多都是建立在大数定理和中心极限定理的前提条件下的,如果它不满足,得到的许多结论都是不可靠的。以虚假回归为例,当响应变量和输入变量都平稳时,我们用t统计量检验标准化系数的显著性。而当响应变量和输入变量不平稳时,其标准化系数不在满足t分布,这时再用t检验来进行显著性分析,导致拒绝原假设的概率增加,即容易犯第一类错误,从而得出错误的结论。

    平稳时间序列有两种定义:严平稳和宽平稳

    严平稳顾名思义,是一种条件非常苛刻的平稳性,它要求序列随着时间的推移,其统计性质保持不变。对于任意的τ,其联合概率密度函数满足:

    1365470-20200130214100419-558255082.png

    严平稳的条件只是理论上的存在,现实中用得比较多的是宽平稳的条件。

    宽平稳也叫弱平稳或者二阶平稳(均值和方差平稳),它应满足:

    常数均值

    常数方差

    常数自协方差

    平稳性检验:观察法和单位根检验法

    基于此,我写了一个名为test_stationarity的统计性检验模块,以便将某些统计检验结果更加直观的展现出来。

    #-*- coding:utf-8 -*-

    from statsmodels.tsa.stattools importadfullerimportpandas as pdimportmatplotlib.pyplot as pltimportnumpy as npfrom statsmodels.graphics.tsaplots importplot_acf, plot_pacf#移动平均图

    defdraw_trend(timeSeries, size):

    f= plt.figure(facecolor='white')#对size个数据进行移动平均

    rol_mean = timeSeries.rolling(window=size).mean()#对size个数据进行加权移动平均

    rol_weighted_mean = pd.DataFrame.ewm(timeSeries, span=size).mean()

    timeSeries.plot(color='blue', label='Original')

    rolmean.plot(color='red', label='Rolling Mean')

    rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')

    plt.legend(loc='best')

    plt.title('Rolling Mean')

    plt.show()defdraw_ts(timeSeries):

    f= plt.figure(facecolor='white')

    timeSeries.plot(color='blue')

    plt.show()'''Unit Root Test

    The null hypothesis of the Augmented Dickey-Fuller is that there is a unit

    root, with the alternative that there is no unit root. That is to say the

    bigger the p-value the more reason we assert that there is a unit root'''

    deftestStationarity(ts):

    dftest=adfuller(ts)#对上述函数求得的值进行语义描述

    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])for key,value in dftest[4].items():

    dfoutput['Critical Value (%s)'%key] =valuereturndfoutput#自相关和偏相关图,默认阶数为31阶

    def draw_acf_pacf(ts, lags=31):

    f= plt.figure(facecolor='white')

    ax1= f.add_subplot(211)

    plot_acf(ts, lags=31, ax=ax1)

    ax2= f.add_subplot(212)

    plot_pacf(ts, lags=31, ax=ax2)

    plt.show()

    观察法,通俗的说就是通过观察序列的趋势图与相关图是否随着时间的变化呈现出某种规律。所谓的规律就是时间序列经常提到的周期性因素,现实中遇到得比较多的是线性周期成分,这类周期成分可以采用差分或者移动平均来解决,而对于非线性周期成分的处理相对比较复杂,需要采用某些分解的方法。下图为航空数据的线性图,可以明显的看出它具有年周期成分和长期趋势成分。平稳序列的自相关系数会快速衰减,下面的自相关图并不能体现出该特征,所以我们有理由相信该序列是不平稳的。

    1365470-20200130214258727-1235545443.png

    1365470-20200130214303800-696996578.png

    单位根检验:ADF是一种常用的单位根检验方法,他的原假设为序列具有单位根,即非平稳,对于一个平稳的时序数据,就需要在给定的置信水平上显著,拒绝原假设。ADF只是单位根检验的方法之一,如果想采用其他检验方法,可以安装第三方包arch,里面提供了更加全面的单位根检验方法,个人还是比较钟情ADF检验。以下为检验结果,其p值大于0.99,说明并不能拒绝原假设。

    1365470-20200130214436436-1646626229.png

    3. 平稳性处理

    由前面的分析可知,该序列是不平稳的,然而平稳性是时间序列分析的前提条件,故我们需要对不平稳的序列进行处理将其转换成平稳的序列。

    a. 对数变换

    对数变换主要是为了减小数据的振动幅度,使其线性规律更加明显(我是这么理解的时间序列模型大部分都是线性的,为了尽量降低非线性的因素,需要对其进行预处理,也许我理解的不对)。对数变换相当于增加了一个惩罚机制,数据越大其惩罚越大,数据越小惩罚越小。这里强调一下,变换的序列需要满足大于0,小于0的数据不存在对数变换。

    ts_log =np.log(ts)

    test_stationarity.draw_ts(ts_log)

    1365470-20200131112947728-324036261.png

    b. 平滑法

    根据平滑技术的不同,平滑法具体分为移动平均法和指数平均法。

    移动平均即利用一定时间间隔内的平均值作为某一期的估计值,而指数平均则是用变权的方法来计算均值

    test_stationarity.draw_trend(ts_log, 12)

    1365470-20200131113022780-494670320.png

    从上图可以发现窗口为12的移动平均能较好的剔除年周期性因素,而指数平均法是对周期内的数据进行了加权,能在一定程度上减小年周期因素,但并不能完全剔除,如要完全剔除可以进一步进行差分操作。

    c. 差分

    时间序列最常用来剔除周期性因素的方法当属差分了,它主要是对等周期间隔的数据进行线性求减。前面我们说过,ARIMA模型相对ARMA模型,仅多了差分操作,ARIMA模型几乎是所有时间序列软件都支持的,差分的实现与还原都非常方便。而statsmodel中,对差分的支持不是很好,它不支持高阶和多阶差分。

    不过没关系,我们有pandas。我们可以先用pandas将序列差分好,然后在对差分好的序列进行ARIMA拟合,只不过这样后面会多了一步人工还原的工作。

    diff_12 = ts_log.diff(12)

    diff_12.dropna(inplace=True)

    diff_12_1= diff_12.diff(1)

    diff_12_1.dropna(inplace=True)

    test_stationarity.testStationarity(diff_12_1)

    1365470-20200131113244379-491062215.png

    从上面的统计检验结果可以看出,经过12阶差分和1阶差分后,该序列满足平稳性的要求了。

    d. 分解

    所谓分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程,它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件一样,statsmodels也支持两类分解模型,加法模型和乘法模型,这里我只实现加法,乘法只需将model的参数设置为”multiplicative”即可。

    from statsmodels.tsa.seasonal importseasonal_decompose

    decomposition= seasonal_decompose(ts_log, model="additive")

    trend=decomposition.trend

    seasonal=decomposition.seasonal

    residual=decomposition.resid

    trend.plot()#seasonal.plot()#residual.plot()

    plt.show()

    1365470-20200131113538638-335942066.png

    得到不同的分解成分后,就可以使用时间序列模型对各个成分进行拟合,当然也可以选择其他预测方法。我曾经用过小波对时序数据进行过分解,然后分别采用时间序列拟合,效果还不错。由于我对小波的理解不是很好,只能简单的调用接口,如果有谁对小波、傅里叶、卡尔曼理解得比较透,可以将时序数据进行更加准确的分解,由于分解后的时序数据避免了他们在建模时的交叉影响,所以我相信它将有助于预测准确性的提高。

    4. 模型识别

    在前面的分析可知,该序列具有明显的年周期与长期成分。对于年周期成分我们使用窗口为12的移动平进行处理,对于长期趋势成分我们采用1阶差分来进行处理。

    rol_mean = ts_log.rolling(window=12).mean()

    rol_mean.dropna(inplace=True)

    ts_diff_1= rol_mean.diff(1)

    ts_diff_1.dropna(inplace=True)

    test_stationarity.testStationarity(ts_diff_1)

    1365470-20200131115155296-1404135668.png

    观察其统计量发现该序列在置信水平为95%的区间下并不显著,我们对其进行再次一阶差分。再次差分后的序列其自相关具有快速衰减的特点,t统计量在99%的置信水平下是显著的,这里我不再做详细说明。

    ts_diff_2 = ts_diff_1.diff(1)

    ts_diff_2.dropna(inplace=True)

    1365470-20200131115446144-1883712969.png

    数据平稳后,需要对模型定阶,即确定p、q的阶数。观察上图,发现自相关和偏相系数都存在拖尾的特点,并且他们都具有明显的一阶相关性,所以我们设定p=1, q=1。下面就可以使用ARMA模型进行数据拟合了。这里我不使用ARIMA(ts_diff_1, order=(1, 1, 1))进行拟合,是因为含有差分操作时,预测结果还原老出问题,至今还没弄明白。

    from statsmodels.tsa.arima_model importARMA

    model= ARMA(ts_diff_2, order=(1, 1))

    result_arma= model.fit( disp=-1, method='css')

    5. 样本拟合

    模型拟合完后,我们就可以对其进行预测了。由于ARMA拟合的是经过相关预处理后的数据,故其预测值需要通过相关逆变换进行还原。

    predict_ts =result_arma.predict()#一阶差分还原

    diff_shift_ts = ts_diff_1.shift(1)

    diff_recover_1=predict_ts.add(diff_shift_ts)#再次一阶差分还原

    rol_shift_ts = rol_mean.shift(1)

    diff_recover=diff_recover_1.add(rol_shift_ts)#移动平均还原

    rol_sum = ts_log.rolling(window=11).sum()

    rol_recover= diff_recover*12 - rol_sum.shift(1)#对数还原

    log_recover =np.exp(rol_recover)

    log_recover.dropna(inplace=True)

    我们使用均方根误差(RMSE)来评估模型样本内拟合的好坏。利用该准则进行判别时,需要剔除“非预测”数据的影响。

    ts = ts[log_recover.index] #过滤没有预测的记录

    plt.figure(facecolor='white')

    log_recover.plot(color='blue', label='Predict')

    ts.plot(color='red', label='Original')

    plt.legend(loc='best')

    plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))

    plt.show()

    1365470-20200131115629840-2081403479.png

    观察上图的拟合效果,均方根误差为11.8828,感觉还过得去。

    6. 完善ARIMA模型

    前面提到statsmodels里面的ARIMA模块不支持高阶差分,我们的做法是将差分分离出来,但是这样会多了一步人工还原的操作。基于上述问题,我将差分过程进行了封装,使序列能按照指定的差分列表依次进行差分,并相应的构造了一个还原的方法,实现差分序列的自动还原。

    #差分操作

    defdiff_ts(ts, d):globalshift_ts_list#动态预测第二日的值时所需要的差分序列

    globallast_data_shift_list

    shift_ts_list=[]

    last_data_shift_list=[]

    tmp_ts=tsfor i ind:

    last_data_shift_list.append(tmp_ts[-i])print(last_data_shift_list)

    shift_ts=tmp_ts.shift(i)

    shift_ts_list.append(shift_ts)

    tmp_ts= tmp_ts -shift_ts

    tmp_ts.dropna(inplace=True)returntmp_ts#还原操作

    defpredict_diff_recover(predict_value, d):ifisinstance(predict_value, float):

    tmp_data=predict_valuefor i inrange(len(d)):

    tmp_data= tmp_data + last_data_shift_list[-i-1]elifisinstance(predict_value, np.ndarray):

    tmp_data=predict_value[0]for i inrange(len(d)):

    tmp_data= tmp_data + last_data_shift_list[-i-1]else:

    tmp_data=predict_valuefor i inrange(len(d)):try:

    tmp_data= tmp_data.add(shift_ts_list[-i-1])except:raise ValueError('What you input is not pd.Series type!')

    tmp_data.dropna(inplace=True)return tmp_data

    现在我们直接使用差分的方法进行数据处理,并以同样的过程进行数据预测与还原。

    #差分

    diffed_ts = diff_ts(ts_log, d=[12, 1])#预测

    from statsmodels.tsa.arima_model importARMA

    model= ARMA(diffed_ts, order=(1, 1))

    result_arma= model.fit( disp=-1, method='css')

    predict_ts=result_arma.predict()#还原

    predict_ts =model.properModel.predict()

    diff_recover_ts= predict_diff_recover(predict_ts, d=[12, 1])

    log_recover=np.exp(diff_recover_ts)#绘制结果

    ts = ts[log_recover.index] #过滤没有预测的记录

    plt.figure(facecolor='white')

    log_recover.plot(color='blue', label='Predict')

    ts.plot(color='red', label='Original')

    plt.legend(loc='best')

    plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))

    plt.show()

    1365470-20200131121655038-2134940672.png

    是不是发现这里的预测结果和上一篇的使用12阶移动平均的预测结果一模一样。这是因为12阶移动平均加上一阶差分与直接12阶差分是等价的关系,后者是前者数值的12倍,这个应该不难推导。

    7. BIC准则

    8. 滚动预测

    9. 模型序列化

    三、总结

    与其它统计语言(SAS和R)相比,python在统计分析这块还显得不那么“专业”。随着numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推动,我相信也祝愿python在数据分析这块越来越好。

    展开全文
  • 一、什么是时间序列时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的...

    一、什么是时间序列

    时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。

    在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里不考虑含外生变量的时间序列)。

    环境配置

    python作为科学计算的利器,当然也有相关分析的包:statsmodels中tsa模块,当然这个包和SAS、R是比不了,但是python有另一个神器:pandas!pandas在时间序列上的应用,能简化我们很多的工作。这两个包pip就能安装。

    数据准备

    许多时间序列分析一样,本文同样使用航空乘客数据(AirPassengers.csv)作为样例。下载链接。

    用pandas操作时间序列

    #-*- coding:utf-8 -*-

    importnumpy as npimportpandas as pdfrom datetime importdatetimeimportmatplotlib.pylab as plt#读取数据,pd.read_csv默认生成DataFrame对象,需将其转换成Series对象

    df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='Month')

    df.index= pd.to_datetime(df.index) #将字符串索引转换成时间索引

    ts = df['Passengers'] #生成pd.Series对象#查看数据格式

    print(ts.head())print(ts.head().index)

    1365470-20200130213601204-441660966.png

    不知道为啥,时间隔1000,不过没什么影响。

    查看某日的值既可以使用字符串作为索引,又可以直接使用时间对象作为索引

    ts['2049-01-01']

    ts[datetime(2049,1,1)]

    两者的返回值都是第一个序列值:112

    如果要查看某一年的数据,pandas也能非常方便的实现

    ts['2049']

    切片操作:

    ts['2049-1' : '2049-6']

    注意时间索引的切片操作起点和尾部都是包含的,这点与数值索引有所不同

    pandas还有很多方便的时间序列函数,在后面的实际应用中在进行说明。

    二、时间序列分析

    1. 基本模型

    自回归移动平均模型(ARMA(p,q))是时间序列中最为重要的模型之一,它主要由两部分组成: AR代表p阶自回归过程,MA代表q阶移动平均过程,其公式如下:

    1365470-20200130213931237-1472383726.png

    1365470-20200130213936859-1941353511.png

    依据模型的形式、特性及自相关和偏自相关函数的特征,总结如下:

    1365470-20200130213959360-1759676234.png

    在时间序列中,ARIMA模型是在ARMA模型的基础上多了差分的操作。

    2. 平稳性检验

    我们知道序列平稳性是进行时间序列分析的前提条件,很多人都会有疑问,为什么要满足平稳性的要求呢?在大数定理和中心定理中要求样本同分布(这里同分布等价于时间序列中的平稳性),而我们的建模过程中有很多都是建立在大数定理和中心极限定理的前提条件下的,如果它不满足,得到的许多结论都是不可靠的。以虚假回归为例,当响应变量和输入变量都平稳时,我们用t统计量检验标准化系数的显著性。而当响应变量和输入变量不平稳时,其标准化系数不在满足t分布,这时再用t检验来进行显著性分析,导致拒绝原假设的概率增加,即容易犯第一类错误,从而得出错误的结论。

    平稳时间序列有两种定义:严平稳和宽平稳

    严平稳顾名思义,是一种条件非常苛刻的平稳性,它要求序列随着时间的推移,其统计性质保持不变。对于任意的τ,其联合概率密度函数满足:

    1365470-20200130214100419-558255082.png

    严平稳的条件只是理论上的存在,现实中用得比较多的是宽平稳的条件。

    宽平稳也叫弱平稳或者二阶平稳(均值和方差平稳),它应满足:

    常数均值

    常数方差

    常数自协方差

    平稳性检验:观察法和单位根检验法

    基于此,我写了一个名为test_stationarity的统计性检验模块,以便将某些统计检验结果更加直观的展现出来。

    #-*- coding:utf-8 -*-

    from statsmodels.tsa.stattools importadfullerimportpandas as pdimportmatplotlib.pyplot as pltimportnumpy as npfrom statsmodels.graphics.tsaplots importplot_acf, plot_pacf#移动平均图

    defdraw_trend(timeSeries, size):

    f= plt.figure(facecolor='white')#对size个数据进行移动平均

    rol_mean = timeSeries.rolling(window=size).mean()#对size个数据进行加权移动平均

    rol_weighted_mean = pd.DataFrame.ewm(timeSeries, span=size).mean()

    timeSeries.plot(color='blue', label='Original')

    rolmean.plot(color='red', label='Rolling Mean')

    rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')

    plt.legend(loc='best')

    plt.title('Rolling Mean')

    plt.show()defdraw_ts(timeSeries):

    f= plt.figure(facecolor='white')

    timeSeries.plot(color='blue')

    plt.show()'''Unit Root Test

    The null hypothesis of the Augmented Dickey-Fuller is that there is a unit

    root, with the alternative that there is no unit root. That is to say the

    bigger the p-value the more reason we assert that there is a unit root'''

    deftestStationarity(ts):

    dftest=adfuller(ts)#对上述函数求得的值进行语义描述

    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])for key,value in dftest[4].items():

    dfoutput['Critical Value (%s)'%key] =valuereturndfoutput#自相关和偏相关图,默认阶数为31阶

    def draw_acf_pacf(ts, lags=31):

    f= plt.figure(facecolor='white')

    ax1= f.add_subplot(211)

    plot_acf(ts, lags=31, ax=ax1)

    ax2= f.add_subplot(212)

    plot_pacf(ts, lags=31, ax=ax2)

    plt.show()

    观察法,通俗的说就是通过观察序列的趋势图与相关图是否随着时间的变化呈现出某种规律。所谓的规律就是时间序列经常提到的周期性因素,现实中遇到得比较多的是线性周期成分,这类周期成分可以采用差分或者移动平均来解决,而对于非线性周期成分的处理相对比较复杂,需要采用某些分解的方法。下图为航空数据的线性图,可以明显的看出它具有年周期成分和长期趋势成分。平稳序列的自相关系数会快速衰减,下面的自相关图并不能体现出该特征,所以我们有理由相信该序列是不平稳的。

    1365470-20200130214258727-1235545443.png

    1365470-20200130214303800-696996578.png

    单位根检验:ADF是一种常用的单位根检验方法,他的原假设为序列具有单位根,即非平稳,对于一个平稳的时序数据,就需要在给定的置信水平上显著,拒绝原假设。ADF只是单位根检验的方法之一,如果想采用其他检验方法,可以安装第三方包arch,里面提供了更加全面的单位根检验方法,个人还是比较钟情ADF检验。以下为检验结果,其p值大于0.99,说明并不能拒绝原假设。

    1365470-20200130214436436-1646626229.png

    3. 平稳性处理

    由前面的分析可知,该序列是不平稳的,然而平稳性是时间序列分析的前提条件,故我们需要对不平稳的序列进行处理将其转换成平稳的序列。

    a. 对数变换

    对数变换主要是为了减小数据的振动幅度,使其线性规律更加明显(我是这么理解的时间序列模型大部分都是线性的,为了尽量降低非线性的因素,需要对其进行预处理,也许我理解的不对)。对数变换相当于增加了一个惩罚机制,数据越大其惩罚越大,数据越小惩罚越小。这里强调一下,变换的序列需要满足大于0,小于0的数据不存在对数变换。

    ts_log =np.log(ts)

    test_stationarity.draw_ts(ts_log)

    1365470-20200131112947728-324036261.png

    b. 平滑法

    根据平滑技术的不同,平滑法具体分为移动平均法和指数平均法。

    移动平均即利用一定时间间隔内的平均值作为某一期的估计值,而指数平均则是用变权的方法来计算均值

    test_stationarity.draw_trend(ts_log, 12)

    1365470-20200131113022780-494670320.png

    从上图可以发现窗口为12的移动平均能较好的剔除年周期性因素,而指数平均法是对周期内的数据进行了加权,能在一定程度上减小年周期因素,但并不能完全剔除,如要完全剔除可以进一步进行差分操作。

    c. 差分

    时间序列最常用来剔除周期性因素的方法当属差分了,它主要是对等周期间隔的数据进行线性求减。前面我们说过,ARIMA模型相对ARMA模型,仅多了差分操作,ARIMA模型几乎是所有时间序列软件都支持的,差分的实现与还原都非常方便。而statsmodel中,对差分的支持不是很好,它不支持高阶和多阶差分。

    不过没关系,我们有pandas。我们可以先用pandas将序列差分好,然后在对差分好的序列进行ARIMA拟合,只不过这样后面会多了一步人工还原的工作。

    diff_12 = ts_log.diff(12)

    diff_12.dropna(inplace=True)

    diff_12_1= diff_12.diff(1)

    diff_12_1.dropna(inplace=True)

    test_stationarity.testStationarity(diff_12_1)

    1365470-20200131113244379-491062215.png

    从上面的统计检验结果可以看出,经过12阶差分和1阶差分后,该序列满足平稳性的要求了。

    d. 分解

    所谓分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程,它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件一样,statsmodels也支持两类分解模型,加法模型和乘法模型,这里我只实现加法,乘法只需将model的参数设置为”multiplicative”即可。

    from statsmodels.tsa.seasonal importseasonal_decompose

    decomposition= seasonal_decompose(ts_log, model="additive")

    trend=decomposition.trend

    seasonal=decomposition.seasonal

    residual=decomposition.resid

    trend.plot()#seasonal.plot()#residual.plot()

    plt.show()

    1365470-20200131113538638-335942066.png

    得到不同的分解成分后,就可以使用时间序列模型对各个成分进行拟合,当然也可以选择其他预测方法。我曾经用过小波对时序数据进行过分解,然后分别采用时间序列拟合,效果还不错。由于我对小波的理解不是很好,只能简单的调用接口,如果有谁对小波、傅里叶、卡尔曼理解得比较透,可以将时序数据进行更加准确的分解,由于分解后的时序数据避免了他们在建模时的交叉影响,所以我相信它将有助于预测准确性的提高。

    4. 模型识别

    在前面的分析可知,该序列具有明显的年周期与长期成分。对于年周期成分我们使用窗口为12的移动平进行处理,对于长期趋势成分我们采用1阶差分来进行处理。

    rol_mean = ts_log.rolling(window=12).mean()

    rol_mean.dropna(inplace=True)

    ts_diff_1= rol_mean.diff(1)

    ts_diff_1.dropna(inplace=True)

    test_stationarity.testStationarity(ts_diff_1)

    1365470-20200131115155296-1404135668.png

    观察其统计量发现该序列在置信水平为95%的区间下并不显著,我们对其进行再次一阶差分。再次差分后的序列其自相关具有快速衰减的特点,t统计量在99%的置信水平下是显著的,这里我不再做详细说明。

    ts_diff_2 = ts_diff_1.diff(1)

    ts_diff_2.dropna(inplace=True)

    1365470-20200131115446144-1883712969.png

    数据平稳后,需要对模型定阶,即确定p、q的阶数。观察上图,发现自相关和偏相系数都存在拖尾的特点,并且他们都具有明显的一阶相关性,所以我们设定p=1, q=1。下面就可以使用ARMA模型进行数据拟合了。这里我不使用ARIMA(ts_diff_1, order=(1, 1, 1))进行拟合,是因为含有差分操作时,预测结果还原老出问题,至今还没弄明白。

    from statsmodels.tsa.arima_model importARMA

    model= ARMA(ts_diff_2, order=(1, 1))

    result_arma= model.fit( disp=-1, method='css')

    5. 样本拟合

    模型拟合完后,我们就可以对其进行预测了。由于ARMA拟合的是经过相关预处理后的数据,故其预测值需要通过相关逆变换进行还原。

    predict_ts =result_arma.predict()#一阶差分还原

    diff_shift_ts = ts_diff_1.shift(1)

    diff_recover_1=predict_ts.add(diff_shift_ts)#再次一阶差分还原

    rol_shift_ts = rol_mean.shift(1)

    diff_recover=diff_recover_1.add(rol_shift_ts)#移动平均还原

    rol_sum = ts_log.rolling(window=11).sum()

    rol_recover= diff_recover*12 - rol_sum.shift(1)#对数还原

    log_recover =np.exp(rol_recover)

    log_recover.dropna(inplace=True)

    我们使用均方根误差(RMSE)来评估模型样本内拟合的好坏。利用该准则进行判别时,需要剔除“非预测”数据的影响。

    ts = ts[log_recover.index] #过滤没有预测的记录

    plt.figure(facecolor='white')

    log_recover.plot(color='blue', label='Predict')

    ts.plot(color='red', label='Original')

    plt.legend(loc='best')

    plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))

    plt.show()

    1365470-20200131115629840-2081403479.png

    观察上图的拟合效果,均方根误差为11.8828,感觉还过得去。

    6. 完善ARIMA模型

    前面提到statsmodels里面的ARIMA模块不支持高阶差分,我们的做法是将差分分离出来,但是这样会多了一步人工还原的操作。基于上述问题,我将差分过程进行了封装,使序列能按照指定的差分列表依次进行差分,并相应的构造了一个还原的方法,实现差分序列的自动还原。

    #差分操作

    defdiff_ts(ts, d):globalshift_ts_list#动态预测第二日的值时所需要的差分序列

    globallast_data_shift_list

    shift_ts_list=[]

    last_data_shift_list=[]

    tmp_ts=tsfor i ind:

    last_data_shift_list.append(tmp_ts[-i])print(last_data_shift_list)

    shift_ts=tmp_ts.shift(i)

    shift_ts_list.append(shift_ts)

    tmp_ts= tmp_ts -shift_ts

    tmp_ts.dropna(inplace=True)returntmp_ts#还原操作

    defpredict_diff_recover(predict_value, d):ifisinstance(predict_value, float):

    tmp_data=predict_valuefor i inrange(len(d)):

    tmp_data= tmp_data + last_data_shift_list[-i-1]elifisinstance(predict_value, np.ndarray):

    tmp_data=predict_value[0]for i inrange(len(d)):

    tmp_data= tmp_data + last_data_shift_list[-i-1]else:

    tmp_data=predict_valuefor i inrange(len(d)):try:

    tmp_data= tmp_data.add(shift_ts_list[-i-1])except:raise ValueError('What you input is not pd.Series type!')

    tmp_data.dropna(inplace=True)return tmp_data

    现在我们直接使用差分的方法进行数据处理,并以同样的过程进行数据预测与还原。

    #差分

    diffed_ts = diff_ts(ts_log, d=[12, 1])#预测

    from statsmodels.tsa.arima_model importARMA

    model= ARMA(diffed_ts, order=(1, 1))

    result_arma= model.fit( disp=-1, method='css')

    predict_ts=result_arma.predict()#还原

    predict_ts =model.properModel.predict()

    diff_recover_ts= predict_diff_recover(predict_ts, d=[12, 1])

    log_recover=np.exp(diff_recover_ts)#绘制结果

    ts = ts[log_recover.index] #过滤没有预测的记录

    plt.figure(facecolor='white')

    log_recover.plot(color='blue', label='Predict')

    ts.plot(color='red', label='Original')

    plt.legend(loc='best')

    plt.title('RMSE: %.4f'% np.sqrt(sum((log_recover-ts)**2)/ts.size))

    plt.show()

    1365470-20200131121655038-2134940672.png

    是不是发现这里的预测结果和上一篇的使用12阶移动平均的预测结果一模一样。这是因为12阶移动平均加上一阶差分与直接12阶差分是等价的关系,后者是前者数值的12倍,这个应该不难推导。

    7. BIC准则

    8. 滚动预测

    9. 模型序列化

    三、总结

    与其它统计语言(SAS和R)相比,python在统计分析这块还显得不那么“专业”。随着numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推动,我相信也祝愿python在数据分析这块越来越好。

    展开全文
  • 什么是时间序列时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里不...
  • python时间序列分析

    千次阅读 2018-04-04 09:51:49
    什么是时间序列 时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里...
  • •1. 时间序列数据的预处理:平稳性... 平稳时间序列数据分析•3. 非平稳时间序列数据分析一.时间序列数据的预处理(1)平稳性检验程序:dataa;inputx @@;time=_n_;cards;17.420 17.9 14.1 12.9 13.6 14.918.216.8...
  • 笔者通过自学,在借鉴了CSDN社区多位大佬的教程和说明,以及参考由中国人民大学出版社出版的《应用时间序列分析》之后,借助SAS统计软件,选取北京市1980-2009年年降水量为对象,进行建模分析。第一次实战操作,笔法...
  • 第五章 非平稳时间序列分析 程序 《应用时间序列分析》第五章所有SAS程序
  • 本课程实验目的和基本要求:培养学生运用所学数学知识,并利用计算机等现代化手段来解决实际问题的综合能力。...掌握SAS/ETS模块进行时间序列分析的一些基本方法和技巧;并逐步了解科学研究的基本思维过程及方法。
  • 时间序列分析 王燕 ppt 代码 全部打包, 注意针对SAS的源代码,不是SPSS等其他类型
  • 平稳序列建模 零、基本概念 1.两种方法工具 →→{\to}差分运算 →→{\to}延迟算子 2. 三种模型 →→{\to}AR模型 p—偏自相关系数 →→{\to}MA模型 q—自相关系数 →→{\to}ARMA模型 p, q ...
  • Python时间序列分析

    万次阅读 2017-05-08 14:15:04
    时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里不考虑含外生...
  • 数学建模之时间序列分析

    千次阅读 2020-07-12 14:48:54
    时序图检验:根据平稳时间序列均值、方差为常数的性质,平稳序列的时序图应该显示出该序列始终在一个常数值附近随机波动,而且波动的范围有界、无明显趋势及周期特征。 自相关图检验:平稳序列通常具有短期相关性。...
  • 本书的独特之处在于它与统计软件包SAS(统计分析系统)计算环境的集成。 通过多元回归假设基本的应用统计数据。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 3,193
精华内容 1,277
关键字:

时间序列分析sas