精华内容
下载资源
问答
  • 我的序列有432个点,一阶差分后有431个点,已经平稳了,p和q也选好了都是1,我想要的是原序列432个点的拟合值 ARIMA输出的拟合结果是一阶差分的拟合值,为什么只有430个点呢?要怎么才能还原成432个原序列的拟合值...
  • 时间序列分析的主要目的是根据已有的历史数据对未来进行预测。如餐饮销售预测可以看做是基于时间序列的短期数据预测, 预测的对象时具体菜品的销售量。1.时间序列算法:常见的时间序列模型;​ 2.时序模型的预处理1. ...

    时间序列分析的主要目的是根据已有的历史数据对未来进行预测。如餐饮销售预测可以看做是基于时间序列的短期数据预测, 预测的对象时具体菜品的销售量。

    1.时间序列算法:

    常见的时间序列模型;

    3a323e95be0459ede8083277252fa576805.png3900054

    2.时序模型的预处理

    1. 对于纯随机序列,也称为白噪声序列,序列的各项之间没有任何的关系, 序列在进行完全无序的随机波动, 可以终止对该序列的分析。

    2. 对于平稳非白噪声序列, 它的均值和方差是常数。ARMA 模型是最常用的平稳序列拟合模型。

    3. 对于非平稳序列, 由于它的方差和均值不稳定, 处理方法一般是将其转化成平稳序列。 可以使用ARIMA 模型进行分析。

    对平稳性的检验:

    1.时序图检验:根据平稳时间序列的均值和方差都是常数的特性,平稳序列的时序图显示该序列值时钟在一个参数附近随机波动,而且波动的范围是有界的。如果有明显的趋势或者周期性, 那它通常不是平稳序列。

    2.自相关图检验:平稳序列具有短期相关性, 这个性质表明对平稳序列而言, 通常 只有近期的序列值得影响比较明显, 间隔越远的过去值对现在的值得影响越小。 而非平稳序列的自相关系数衰减的速度比较慢。

    3.单位根检验:单位根检验是指检验序列中是否存在单位根, 如果存在单位根, 那就是非平稳时间序列。 目前最常用的方法就是单位根检验。

    原假设是 非平稳序列过程, 备择假设是 平稳序列, 趋势平稳过程

    60cace347a2d398eb88d1029a705ac53ba9.jpg3900054

    2353a4222c8e9afbd744ba67069414b06a8.jpg3900054

    2b595ff7a9aa3ebabeb1e0cf2c5fd81fcbd.jpg3900054

    cfe8bd136e2a5ce7666c637230c7954cedc.jpg3900054

    edc2fddcac9c00cf19febac9e12d6b627c2.jpg3900054

    3.时间序列分析:

    •平稳性:

    •平稳性要求经由样本时间序列所得到的拟合曲线,在未来一段时间内仍然沿着现有的形态‘惯性’地延续下去。

    •平稳性要求序列的均值和方差不发生明显的变化。

    •弱平稳:期望和相关系数(依赖性)不变,未来某个时刻t 的值,Xt 要依赖于它过去的信息。

    •差分法:时间序列在 T 与 T-1 时刻的差值(使用差分使其满足平稳性),一般差分1,2 阶就可以了。

    •AR(自回归模型):

    •描述当前值与历史值之间的关系, 用变量自身的历史时间数据对自身进行预测。自回归模型必须满足平稳性的要求。

    公式定义:2dfaa46f24f085d1b6a2b9670d327a9b23d.png3900054

    edcc4d03f15352494070a9df2ebfc112cc2.png3900054

    d177b440a92e5fdf926a5415956f8b2e3bd.png3900054

    自回归模型的限制:

    1.自回归模型是使用自身的数据进行预测的

    2.必须具有平稳性

    3.必须具有相关性,如果相关性小于 0.5 , 则不宜使用

    4.自回归模型只适用于预测与自身前期相关的预测。

    •MA(移动平均模型):

    •移动平均模型关注的是自回归模型中的误差项的累加

    •移动平均法能有效地消除预测中的随机波动。

    b7f08ccb401ae11212f2a5b2bf024e567b4.png3900054

    •ARMA(自回归平均模型):

    •自回归和移动平均的结合。

    a897f4fa4ff82fad75702f6e6ceb64efe5a.png3900054

    •ARIMA(p,d,q)差分自回归移动平均模型(Autoregressive Integrated Moving Average Model ,简称ARIMA)

    •AR 是自回归, p 是自回归项, MA 是移动平均, q 为移动平均项, d 为时间序列称为平稳时 所做的差分次数。

    •原理: 将非平稳时间序列转换成平稳时间序列, 然后将因变量仅对它的滞后值(p阶)以及随机误差项的现值和滞后值进行回顾所建立的模型。

    •ARIMA 建模流程:

    •1.将序列平稳化(差分法确定 d)

    •2.p 和 q 阶数的确定(ACF 和 PACF 确定)

    •3.建立模型 ARIMA (p , d , q )

    f6304dafaf7f3998948504de7369813d21f.png3900054

    使用ARIMA 模型对某餐厅的销售数据进行预测

    #使用ARIMA 模型对非平稳时间序列进行建模操作

    #差分运算具有强大的确定性的信息提取能力, 许多非平稳的序列差分后显示出平稳序列的性质, 这是称这个非平稳序列为差分平稳序列。

    #对差分平稳序列可以还是要ARMA 模型进行拟合, ARIMA 模型的实质就是差分预算与 ARMA 模型的结合。

    #coding=gbk

    #使用ARIMA 模型对非平稳时间序列记性建模操作

    #差分运算具有强大的确定性的信息提取能力, 许多非平稳的序列差分后显示出平稳序列的性质, 这是称这个非平稳序列为差分平稳序列。

    #对差分平稳序列可以还是要ARMA 模型进行拟合, ARIMA 模型的实质就是差分预算与 ARMA 模型的结合。

    #导入数据 import pandas as pd filename = r'D:\datasets\arima_data.xls' data = pd.read_excel(filename, index_col = u'日期') #画出时序图 import matplotlib.pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] #定义使其正常显示中文字体黑体 plt.rcParams['axes.unicode_minus'] = False #用来正常显示表示负号 # data.plot() # plt.show()

    3900054

    c1ee862f6855926a853da4e08cea27684b6.png3900054

    #画出自相关性图

    from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

    # plot_acf(data)

    # plt.show() #平稳性检测 from statsmodels.tsa.stattools import adfuller print('原始序列的检验结果为:',adfuller(data[u'销量'])) #原始序列的检验结果为: (1.8137710150945268, 0.9983759421514264, 10, 26, {'1%': -3.7112123008648155, # '10%': -2.6300945562130176, '5%': -2.981246804733728}, 299.46989866024177) #返回值依次为:adf, pvalue p值, usedlag, nobs, critical values临界值 , icbest, regresults, resstore #adf 分别大于3中不同检验水平的3个临界值,单位检测统计量对应的p 值显著大于 0.05 , 说明序列可以判定为 非平稳序列

    3900054

    aaf2a3db27e23e9e7c7e3f89e84ea3f9834.png3900054

    #对数据进行差分后得到 自相关图和 偏相关图

    D_data = data.diff().dropna()

    D_data.columns = [u'销量差分']

    D_data.plot() #画出差分后的时序图

    # plt.show()

    plot_acf(D_data) #画出自相关图 # plt.show() plot_pacf(D_data) #画出偏相关图 # plt.show() print(u'差分序列的ADF 检验结果为: ', adfuller(D_data[u'销量差分'])) #平稳性检验 #差分序列的ADF 检验结果为: (-3.1560562366723537, 0.022673435440048798, 0, 35, {'1%': -3.6327426647230316, # '10%': -2.6130173469387756, '5%': -2.9485102040816327}, 287.5909090780334) #一阶差分后的序列的时序图在均值附近比较平稳的波动, 自相关性有很强的短期相关性, 单位根检验 p值小于 0.05 ,所以说一阶差分后的序列是平稳序列

    3900054

    eb76299b3bd2b76b2c01ed1b15280848b35.png3900054

    746e9b83c00f5bbc8e47b09b5233bb48e16.png3900054

    cf1a96328506f13d26fd0822fea6f573c69.png3900054

    #对一阶差分后的序列做白噪声检验

    from statsmodels.stats.diagnostic import acorr_ljungbox

    print(u'差分序列的白噪声检验结果:',acorr_ljungbox(D_data, lags= 1)) #返回统计量和 p 值 # 差分序列的白噪声检验结果: (array([11.30402222]), array([0.00077339])) p值为第二项, 远小于 0.05 #对模型进行定阶 from statsmodels.tsa.arima_model import ARIMA pmax = int(len(D_data) / 10) #一般阶数不超过 length /10 qmax = int(len(D_data) / 10) bic_matrix = [] for p in range(pmax +1): temp= [] for q in range(qmax+1): try: temp.append(ARIMA(data, (p, 1, q)).fit().bic) except: temp.append(None) bic_matrix.append(temp) bic_matrix = pd.DataFrame(bic_matrix) #将其转换成Dataframe 数据结构 p,q = bic_matrix.stack().idxmin() #先使用stack 展平, 然后使用 idxmin 找出最小值的位置 print(u'BIC 最小的p值 和 q 值:%s,%s' %(p,q)) # BIC 最小的p值 和 q 值:0,1 #所以可以建立ARIMA 模型,ARIMA(0,1,1) model = ARIMA(data, (p,1,q)).fit() model.summary2() #生成一份模型报告 model.forecast(5) #为未来5天进行预测, 返回预测结果, 标准误差, 和置信区间

    3900054

    利用模型向前预测的时期越长, 预测的误差就会越大,这是时间预测的典型特点。

    9aa5eeae09b9f164dcc407a630f080b2777.png3900054

    展开全文
  • 一阶差分后的序列:nan,a2−a1,a3−a2,a4−a3,a5−a4,a6−a5,a7−a6nan,a_2-a_1,a_3-a_2,a_4-a_3,a_5-a_4,a_6-a_5,a_7-a_6nan,a2​−a1​,a3​−a2​,a4​−a3​,a5​−a4​,a6​−a5​,a7​−a6​ 一阶差分后的...

    一、思路:

    1、差分

    原始序列: a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 a_1,a_2,a_3,a_4,a_5,a_6,a_7 a1,a2,a3,a4,a5,a6,a7
    一阶差分后的序列: n a n , a 2 − a 1 , a 3 − a 2 , a 4 − a 3 , a 5 − a 4 , a 6 − a 5 , a 7 − a 6 nan,a_2-a_1,a_3-a_2,a_4-a_3,a_5-a_4,a_6-a_5,a_7-a_6 nana2a1,a3a2,a4a3,a5a4,a6a5,a7a6

    一阶差分后的序列: n a n , n a n , a 3 − a 2 − ( a 2 − a 1 ) , a 4 − a 3 − ( a 3 − a 2 ) , a 5 − a 4 − ( a 4 − a 3 ) , a 6 − a 5 − ( a 5 − a 4 ) , a 7 − a 6 − ( a 6 − a 5 ) nan,nan,a_3-a_2-(a_2-a_1),a_4-a_3-(a_3-a_2),a_5-a_4-(a_4-a_3),a_6-a_5-(a_5-a_4),a_7-a_6-(a_6-a_5) nan,nan,a3a2(a2a1),a4a3(a3a2),a5a4(a4a3),a6a5(a5a4),a7a6(a6a5)

    2、还原

    2.1 一阶差分还原

    一阶差分还原我们必须有一个原始值 a 1 a_1 a1,可以看为测试集的最后一个值,接下来还原就很简单

    • 一阶差分结果求累加: n a n , a 2 − a 1 , a 3 − a 1 , a 2 − a 1 , a 3 − a 1 , a 4 − a 1 , a 5 − a 1 , a 6 − a 1 , a 7 − a 1 nan,a_2-a_1,a_3-a_1,a_2-a_1,a_3-a_1,a_4-a_1,a_5-a_1,a_6-a_1,a_7-a_1 nan,a2a1,a3a1,a2a1,a3a1,a4a1,a5a1,a6a1,a7a1
    • 对累加后的序列加上原始值 a 1 a_1 a1 a 2 , a 3 , a 4 , a 5 , a 6 , a 7 a_2,a_3,a_4,a_5,a_6,a_7 a2,a3,a4,a5,a6,a7

    2.2 二阶差分还原

    二阶差分还原我们必须有两个原始值 a 1 , a 2 a_1,a_2 a1a2,可以看为测试集的最后两个值,接下来还原就很简单

    • 先还原到一阶差分后的结果
      • 二阶差分结果求累加: n a n , n a n , a 3 − a 2 − ( a 2 − a 1 ) , a 4 − a 3 − ( a 2 − a 1 ) , a 5 − a 4 − ( a 2 − a 1 ) , a 6 − a 5 − ( a 2 − a 1 ) , a 7 − a 6 − ( a 2 − a 1 ) nan,nan,a_3-a_2-(a_2-a_1),a_4-a_3-(a_2-a_1),a_5-a_4-(a_2-a_1),a_6-a_5-(a_2-a_1),a_7-a_6-(a_2-a_1) nan,nan,a3a2(a2a1),a4a3(a2a1),a5a4(a2a1),a6a5(a2a1),a7a6(a2a1)
      • 对累加后的序列加上原始值 a 2 − a 1 a_2-a_1 a2a1 n a n , n a n , a 3 − a 2 , a 4 − a 3 , a 5 − a 4 , a 6 − a 5 , a 7 − a 6 nan,nan,a_3-a_2,a_4-a_3,a_5-a_4,a_6-a_5,a_7-a_6 nannan,a3a2,a4a3,a5a4,a6a5,a7a6
    • 重复2.1的步骤

    二、python实现

    import pandas as pd
    
    
    def diff1_reduction(y_pred_test_1, y_train_real):
        '''
        :param y_pred_test_1: 一阶差分后测试集的预测结果
        :param y_train_real: 训练集的输出
        :return: 测试集的真实预测结果(未经差分)
        '''
        y_test_real = y_pred_test_1.cumsum() + y_train_real.values[-1]
        return y_test_real.dropna()
    
    
    def diff2_reduction(y_pred_test_2, y_train_real):
        '''
        :param y_pred_test_1: 二阶差分后测试集的预测结果
        :param y_train_real: 训练集的输出
        :return: 测试集的真实预测结果(未经差分)
        '''
    
        data_diff_1_reduction = y_pred_test_2.cumsum() + (
                    y_train_real.values[-1] - y_train_real.values[-2])  # 先还原到y_pred_test_1
        y_test_real = diff1_reduction(data_diff_1_reduction, y_train_real)
        return y_test_real
    
    
    if __name__ == '__main__':
        data = pd.DataFrame([1, 3, 4, 6, 8, 9, 11, 13, 12])
        data_diff_1 = data.diff(1)  # [nan,  2.,  1.,  2.,  2.,  1.,  2.,  2., -1.]
        data_diff_2 = data_diff_1.diff()  # [nan, nan, -1.,  1.,  0., -1.,  1.,  0., -3.]
    
        # 为了更好地代入ARIMA模型的预测过程
        X_train = data[:2]  # [1, 3]
        X_test = data[2:]  # [4, 6, 8, 9, 11, 13, 12]
    
        y_train_real = data[:2]  # 训练集的真实值
        y_test_real = data[2:]  # 测试集的真实值 [4, 6, 8, 9, 11, 13, 12]
        y_pred_test_1 = data_diff_1[2:]  # 测试集一阶差分后的预测值
        y_pred_test_2 = data_diff_2[2:]  # 测试集二阶差分后的预测值
    
        """一阶差分还原"""
        y_test_real_1 = diff1_reduction(y_pred_test_1, y_train_real)  # [4, 6, 8, 9, 11, 13, 12]
    
        """二阶差分还原"""
        y_test_real_2 = diff2_reduction(data_diff_2, y_train_real)  # [4, 6, 8, 9, 11, 13, 12] #
    
    
    展开全文
  • 什么是时间序列时间序列简单的说就是各时间点上形成的数值序列时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里不...

    什么是时间序列

    时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的(这里不考虑含外生变量的时间序列)。

    为什么用python

    用两个字总结“情怀”,爱屋及乌,个人比较喜欢python,就用python撸了。能做时间序列的软件很多,SAS、R、SPSS、Eviews甚至matlab等等,实际工作中应用得比较多的应该还是SAS和R,前者推荐王燕写的《应用时间序列分析》,后者推荐“

    环境配置

    python推荐直接装Anaconda,它集成了许多科学计算包,有一些包自己手动去装还是挺费劲的。statsmodels需要自己去安装,这里我推荐使用0.6的稳定版,0.7及其以上的版本能在github上找到,该版本在安装时会用C编译好,所以修改底层的一些代码将不会起作用。

    时间序列分析

    1.基本模型

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

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

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

    2.pandas时间序列操作

    大熊猫真的很可爱,这里简单介绍一下它在时间序列上的可爱之处。和许多时间序列分析一样,本文同样使用航空乘客数据(AirPassengers.csv)作为样例。

    数据读取:

    # -*- coding:utf-8 -*- import numpy as np import pandas as pd

    from datetime import datetime

    import matplotlib.pylab as plt

    # 读取数据,pd.read_csv默认生成DataFrame对象,需将其转换成Series对象

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

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

    ts = df['x'] # 生成pd.Series对象

    # 查看数据格式

    ts.head()

    ts.head().index

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

    ts['1949-01-01'] ts[datetime(1949,1,1)]

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

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

    ts['1949']

    切片操作:

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

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

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

    3. 平稳性检验

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

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

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

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

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

    常数方差

    常数自协方差

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

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

    # -*- coding:utf-8 -*- from statsmodels.tsa.stattools import adfuller import pandas as pd import matplotlib.pyplot as plt import numpy as np from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

    # 移动平均图 def draw_trend(timeSeries, size): f = plt.figure(facecolor='white') # 对size个数据进行移动平均 rol_mean = timeSeries.rolling(window=size).mean() # 对size个数据进行加权移动平均 rol_weighted_mean = pd.ewma(timeSeries, span=size) 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() def draw_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 ''' def testStationarity(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] = value return dfoutput # 自相关和偏相关图,默认阶数为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()

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

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

    3. 平稳性处理

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

    a. 对数变换

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

    ts_log = np.log(ts) test_stationarity.draw_ts(ts_log)

    b. 平滑法

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

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

    test_stationarity.draw_trend(ts_log, 12)

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

    c. 差分

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

    作者大概的意思是说预测方法中并没有解决高于2阶的差分,有没有感觉很牵强,不过没关系,我们有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)

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

    d. 分解

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

    from statsmodels.tsa.seasonal import seasonal_decompose

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

    trend = decomposition.trend

    seasonal = decomposition.seasonal

    residual = decomposition.resid

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

    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)

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

    ts_diff_2 = ts_diff_1.diff(1) ts_diff_2.dropna(inplace=True)

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

    from statsmodels.tsa.arima_model import ARMA

    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()

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

    6. 完善ARIMA模型

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

    # 差分操作

    def diff_ts(ts, d):

    global shift_ts_list # 动态预测第二日的值时所需要的差分序列

    global last_data_shift_list

    shift_ts_list = []

    last_data_shift_list = []

    tmp_ts = ts

    for i in d:

    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)

    return tmp_ts # 还原操作

    def predict_diff_recover(predict_value, d):

    if isinstance(predict_value, float):

    tmp_data = predict_value

    for i in range(len(d)):

    tmp_data = tmp_data + last_data_shift_list[-i-1]

    elif isinstance(predict_value, np.ndarray):

    tmp_data = predict_value[0]

    for i in range(len(d)):

    tmp_data = tmp_data + last_data_shift_list[-i-1]

    else:

    tmp_data = predict_value

    for i in range(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])

    model = arima_model(diffed_ts)

    model.certain_model(1, 1)

    predict_ts = model.properModel.predict()

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

    log_recover = np.exp(diff_recover_ts)

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

    对于个数不多的时序数据,我们可以通过观察自相关图和偏相关图来进行模型识别,倘若我们要分析的时序数据量较多,例如要预测每只股票的走势,我们就不可能逐个去调参了。这时我们可以依据BIC准则识别模型的p, q值,通常认为BIC值越小的模型相对更优。这里我简单介绍一下BIC准则,它综合考虑了残差大小和自变量的个数,残差越小BIC值越小,自变量个数越多BIC值越大。个人觉得BIC准则就是对模型过拟合设定了一个标准(过拟合这东西应该以辩证的眼光看待)。

    def proper_model(data_ts, maxLag): init_bic = sys.maxint init_p = 0 init_q = 0 init_properModel = None for p in np.arange(maxLag): for q in np.arange(maxLag): model = ARMA(data_ts, order=(p, q)) try: results_ARMA = model.fit(disp=-1, method='css') except: continue bic = results_ARMA.bic if bic < init_bic: init_p = p init_q = q init_properModel = results_ARMA init_bic = bic return init_bic, init_p, init_q, init_properModel

    相对最优参数识别结果:BIC: -1090.44209358 p: 0 q: 1 , RMSE:11.8817198331。我们发现模型自动识别的参数要比我手动选取的参数更优。

    7.滚动预测

    所谓滚动预测是指通过添加最新的数据预测第二天的值。对于一个稳定的预测模型,不需要每天都去拟合,我们可以给他设定一个阀值,例如每周拟合一次,该期间只需通过添加最新的数据实现滚动预测即可。基于此我编写了一个名为arima_model的类,主要包含模型自动识别方法,滚动预测的功能,详细代码可以查看附录。数据的动态添加:

    from dateutil.relativedelta import relativedelta

    def _add_new_data(ts, dat, type='day'):

    if type == 'day': new_index = ts.index[-1] + relativedelta(days=1) elif type == 'month': new_index = ts.index[-1] + relativedelta(months=1) ts[new_index] = dat def add_today_data(model, ts, data, d, type='day'): _add_new_data(ts, data, type) # 为原始序列添加数据 # 为滞后序列添加新值 d_ts = diff_ts(ts, d) model.add_today_data(d_ts[-1], type) def forecast_next_day_data(model, type='day'): if model == None: raise ValueError('No model fit before') fc = model.forecast_next_day_value(type) return predict_diff_recover(fc, [12, 1])

    现在我们就可以使用滚动预测的方法向外预测了,取1957年之前的数据作为训练数据,其后的数据作为测试,并设定模型每第七天就会重新拟合一次。这里的diffed_ts对象会随着add_today_data方法自动添加数据,这是由于它与add_today_data方法中的d_ts指向的同一对象,该对象会动态的添加数据。

    ts_train = ts_log[:'1956-12']

    ts_test = ts_log['1957-1':]

    diffed_ts = diff_ts(ts_train, [12, 1])

    forecast_list = []

    for i, dta in enumerate(ts_test):

    if i%7 == 0:

    model = arima_model(diffed_ts)

    model.certain_model(1, 1)

    forecast_data = forecast_next_day_data(model, type='month')

    forecast_list.append(forecast_data)

    add_today_data(model, ts_train, dta, [12, 1], type='month')

    predict_ts = pd.Series(data=forecast_list, index=ts['1957-1':].index)

    log_recover = np.exp(predict_ts)

    original_ts = ts['1957-1':]

    动态预测的均方根误差为:14.6479,与前面样本内拟合的均方根误差相差不大,说明模型并没有过拟合,并且整体预测效果都较好。

    8. 模型序列化

    在进行动态预测时,我们不希望将整个模型一直在内存中运行,而是希望有新的数据到来时才启动该模型。这时我们就应该把整个模型从内存导出到硬盘中,而序列化正好能满足该要求。序列化最常用的就是使用json模块了,但是它是时间对象支持得不是很好,有人对json模块进行了拓展以使得支持时间对象,这里我们不采用该方法,我们使用pickle模块,它和json的接口基本相同,有兴趣的可以去查看一下。我在实际应用中是采用的面向对象的编程,它的序列化主要是将类的属性序列化即可,而在面向过程的编程中,模型序列化需要将需要序列化的对象公有化,这样会使得对前面函数的参数改动较大,我不在详细阐述该过程。

    总结

    与其它统计语言相比,python在统计分析这块还显得不那么“专业”。随着numpy、pandas、scipy、sklearn、gensim、statsmodels等包的推动,我相信也祝愿python在数据分析这块越来越好。与SAS和R相比,python的时间序列模块还不是很成熟,我这里仅起到抛砖引玉的作用,希望各位能人志士能贡献自己的力量,使其更加完善。实际应用中我全是面向过程来编写的,为了阐述方便,我用面向过程重新罗列了一遍,实在感觉很不方便。原本打算分三篇来写的,还有一部分实际应用的部分,不打算再写了,还请大家原谅。实际应用主要是具体问题具体分析,这当中第一步就是要查询问题,这步花的时间往往会比较多,然后再是解决问题。以我前面项目遇到的问题为例,当时遇到了以下几个典型的问题:1.周期长度不恒定的周期成分,例如每月的1号具有周期性,但每月1号与1号之间的时间间隔是不相等的;2.含有缺失值以及含有记录为0的情况无法进行对数变换;3.节假日的影响等等。

    展开全文
  • 时间序列的模型主要分四种:自回归模型AR,移动回归模型MA,两者的结合移动自回归模型ARMA,以及差分过的差分移动自回归模型ARIMA。1、AR模型:Xt时刻的值等于自回归系数乘上对应时刻的数值,ut为时间序列的随机游走...

    今天给大家分享python实现时间序列的案例。时间序列的模型主要分四种:自回归模型AR,移动回归模型MA,两者的结合移动自回归模型ARMA,以及差分过的差分移动自回归模型ARIMA。

    1、AR模型:

    Xt时刻的值等于自回归系数乘上对应时刻的数值,ut为时间序列的随机游走。

    2、MA模型:

    Xt时刻的数值为每个时刻的白噪声的系数的加权和。当自回归和移动回归结合就是ARMA。

    3、ARMA模型:

    自回归移动平均模型由两部分组成:自回归部分和移动平均部分,因此包含两个阶数,可以表示为ARMA(p,q),p是自回归阶数,q为移动平均阶数,回归方程表示为:

    ARMA既体现了每个时刻数据对于预测数据的影响也提现了随机游走因素对于预测数据的影响,模型应用比较常见。

    4、ARIMA模型:

    差分移动自回归模型,由于前三个模型都有时间序列平稳的假设在,如果时间序列存在明显的上升或者下降趋势,模型预测的效果大大折扣,对于这些有明显下降或者上升趋势的数据集,可以使用差分的方式使它们平稳,之后使用ARMA拟合。假设模型经过d次差分通过了时间序列平稳的检验,ARMA的系数为p,q,ARIMA模型为ARIMA(p,d,q)。

    案例:

    from __future__ import print_function

    import pandas as pd

    import numpy as np

    from scipy import stats

    import matplotlib.pyplot as plt

    import statsmodels.api as sm

    from statsmodels.graphics.api import qqplot

    #时间序列

    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=np.array(dta,dtype=np.float)

    #生成时间序列并画图

    dta=pd.Series(dta)

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

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

    有上升趋势需要差分

    #一阶差分

    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)

    可以观察一阶差分平稳,做单位根检验

    #一阶单位根检验

    sm.tsa.stattools.adfuller(diff1[1:])

    #二阶单位根检验

    sm.tsa.stattools.adfuller(diff2[2:])

    #选择使用一阶差分的时间序列

    diff1= 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)

    自相关图和偏相关图

    根据上图,猜测有以下模型可以供选择:

    1)ARMA(0,1)模型:即自相关图在滞后1阶之后缩小为0,且偏自相关缩小至0,则是一个阶数q=1的移动平均模型;

    2.)ARMA(7,0)模型:即偏自相关图在滞后7阶之后缩小为0,且自相关缩小至0,则是一个阶层p=7的自回归模型

    3.)ARMA(7,1)模型:即使得自相关和偏自相关都缩小至零。则是一个混合模型。

    4) 还可以有其他供选择的模型 ARMA(8,0)

    #不同的值的p,q拟合

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

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

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

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

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

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

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

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

    看不同的信息量:

    1619.191769832876 1641.6900568658484 1628.2644016400693

    1657.2172636424214 1664.716692653412 1660.241474244819

    1605.686728785649 1630.6848254889517 1615.7674307936416

    1597.9359854086617 1622.9340821119645 1608.0166874166543

    选择 ARMA(8,0)

    resid = arma_mod80.resid

    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)

    plt.show()

    残差的ACF和PACF图,可以看到序列残差基本为白噪声。

    德宾-沃森检验:

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

    QQplot:

    print(stats.normaltest(resid))

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

    ax = fig.add_subplot(111)

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

    #plt.show()

    预测:

    predict_dta = arma_mod80.predict('2090', '2150', dynamic=True)

    print(predict_dta)

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

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

    fig = arma_mod80.plot_predict('2090', '2150', dynamic=True, ax=ax, plot_insample=False)

    plt.show()

    时间序列的初步建模并预测了之后的数据。

    如果对数据分析和挖掘可以关注我的公众号:论数据如何炼成

    展开全文
  • 时间序列模型:严格来说包含4个要素,Trend/趋势、Circle/循环、Seasonal /季节性和不规则要素。...在平稳化时间序列数据中,差分/differencing是种用得广&受欢迎的方法。 笔记的目的是为了理解: 平稳的时间...
  • 在应用ARIMA模型时,要保证以下几点:时间序列数据是相对稳定的,总体基本不存在一定的上升或者下降趋势,如果不稳定可以通过差分的方式来使其变稳定。非线性关系处理不好,只能处理线性关系判断时序数据稳定基本...
  • 一、什么是时间序列时间序列简单的说就是各时间点上形成的数值序列时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的...
  • 在应用ARIMA模型时,要保证以下几点:时间序列数据是相对稳定的,总体基本不存在一定的上升或者下降趋势,如果不稳定可以通过差分的方式来使其变稳定。非线性关系处理不好,只能处理线性关系判断时序数据稳定基本...
  • 一、什么是时间序列时间序列简单的说就是各时间点上形成的数值序列时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的...
  • python时间序列分析

    万次阅读 多人点赞 2018-03-05 22:45:49
    本文转载自博客园大神“大熊猫淘沙”的一篇文章——python时间序列分析。 文章写的生动有趣干货满满,特此收藏转载一下。原文地址:https://www.cnblogs.com/foley/p/5582358.html 1. 什么是时间序列 1.1 环境...
  • 1 时间序列时间序列分析在生产和科学研究中,对某一个或者一组变量 进行观察测量,将在一系列时刻 所得到的离散数字组成的序列集合,称之为时间序列时间序列分析是根据系统观察得到的时间序列数据,通过曲线拟合...
  • 在应用ARIMA模型时,要保证以下几点:时间序列数据是相对稳定的,总体基本不存在一定的上升或者下降趋势,如果不稳定可以通过差分的方式来使其变稳定。非线性关系处理不好,只能处理线性关系判断时序数据...
  • 01 引言上一篇推文《Python量化基础:时间序列的自相关性与平稳性》着重介绍了时间序列的一些基础...时间序列经典模型主要有自回归模型AR,移动回归模型MA,移动自回归模型ARMA,以及差分移动自回归模型ARIMA,今天...
  • 1, pandas生成时间一般采用date_range操作,这个之前的博客已经详细的讲解过,这里就不在阐述2, pandas的数据重采样什么是数据重采样?就好比原来一堆统计数据是按照天来进行统计的,持续一年;那我们能不能看月...
  • airpassenger.csv:https://pan.baidu.com/s/1mi7dp6w...时间序列分析时间序列分析是根据系统观察得到的时间序列数据,通过曲线拟合和参数估计来建立数学模型的理论和方法。时间序列分析常用于国民宏观经济控制、市...
  • python时间序列数据预测教程之 arima

    千次阅读 2020-04-14 23:27:01
    最近接触时间序列较多,在借鉴很多人的知识之后,特此总结一下。目前关于时间序列数据分析预测大致有三种主流方法: 1、ARIMA系列方法 2、facebook开源的Prophet模型 3、LSTM时间序列预测 本系列将用python进行实现...
  • Python时间序列案例分析实战--奶牛产奶量预测

    千次阅读 多人点赞 2017-10-23 17:53:55
    然后使用test_stationarity函数进行原始数据以及一阶差分数据的平稳性检验,用于满足后边ARIMA模型对于平稳性的要求;之后对数据进行差分,使数据满足平稳性;之后使用网格搜索的方法对最优参数进行初步确定,之后...
  • Python时间序列分析

    万次阅读 2017-05-08 14:15:04
    什么是时间序列 时间序列简单的说就是各时间点上形成的数值序列时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是,时间序列分析并不是关于时间的回归,它主要是研究自身的变化规律的...
  • Python时间序列分析--ARIMA模型实战案例

    万次阅读 多人点赞 2020-12-22 10:09:46
    本文将介绍使用Python来完成时间序列分析ARIMA模型的完整步骤与流程 时间序列分析概念 **《时间序列分析》**是统计学中的一个非常重要的分支,是以概率论与数理统计为基础、计算机应用为技术支撑,迅速发展起来的...
  • python 时间序列分析

    千次阅读 2017-03-24 18:29:20
    1 时间序列时间序列分析 在生产和科学研究中,对某一个或者一组变量 x(t) 进行观察测量,将在一系列时刻 t1,t2,⋯,tn 所得到的离散数字组成的序列集合,称之为时间序列时间序列分析是根据系统观察得到的...
  • Python时间序列处理之ARIMA模型

    千次阅读 2018-08-22 00:39:35
    ARIMA模型 ARIMA模型的全称是自回归移动平均模型...时间序列数据是相对稳定的,总体基本不存在一定的上升或者下降趋势,如果不稳定可以通过差分的方式来使其变稳定。 非线性关系处理不好,只能处理线性关系 判...
  • 更新:2020/09/02:新增一篇写的很好的关于视频监控... 介绍异常检测(Anomaly detection)是目前时序数据分析最成熟的应用之一,定义是从正常的时间序列中识别不正常的事件或行为的过程。有效的异常检测被广泛用于...
  • python时间序列处理1-预处理时间序列平稳性判断时序图检验自相关图检验纯随机性检验LB检验Q检验 时间序列平稳性判断 时序图检验 import pandas as pd import numpy as np import seaborn as sns import ...

空空如也

空空如也

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

python时间序列一阶差分

python 订阅