精华内容
下载资源
问答
  • 旅游季节性是一个旅游研究的常态问题.本文对近年来旅游季节性研究成果进行了综述,有助于深入了解旅游季节性影响因素以及动力机制,进而可以从实际出发提出解决旅游季节性矛盾的主要策略,期望更有利于理论服务实践活动...
  • 季节性时间序列建模与预测 ,孟玲清,王晓雨,所谓所谓时间序列,就是各种社会、经济、自然现象的数量指标按照时间次序排列起来的统计数据,其有多种构成因素,每种因素对系统
  • 趋势性、季节性、周期性

    千次阅读 2020-07-14 10:23:44
    以上都是基本概念,重点是区分季节性和周期性。事物都是盛极必衰,但是跌落谷底之后又反弹,这种往复的运动,被称之为周期性。但是像日出日落、春夏秋冬、周六日放假这种节律如果能够固定下来,那么这种固定的往复...

    时间序列有一些比较常见的特征,包括:
    1、趋势性(Trend):即在一定时间内的单调性,一般来说斜率是固定的。
    2、季节性(Seasonal):固定长度的变化,就像春夏秋冬的温度变化一样。
    3、周期性(Cyclic):与季节性很像,但是它的波动的时间频率不是固定的。
    以上都是基本概念,重点是区分季节性和周期性。事物都是盛极必衰,但是跌落谷底之后又反弹,这种往复的运动,被称之为周期性。但是像日出日落、春夏秋冬、周六日放假这种节律如果能够固定下来,那么这种固定的往复变化就被成为季节性。

    展开全文
  • 时间序列预测,非季节性ARIMA及季节性SARIMA

    万次阅读 多人点赞 2019-03-24 21:55:00
    我们将首先介绍和讨论自相关,平稳性和季节性概念,并继续应用最常用的时间序列预测方法之一,称为ARIMA。 介绍 时间序列提供了预测未来价值的机会。基于以前的价值观,可以使用时间序列来预测经济,天气和...

    Python 3中使用ARIMA进行时间序列预测的指南

    在本教程中,我们将提供可靠的时间序列预测。我们将首先介绍和讨论自相关,平稳性和季节性的概念,并继续应用最常用的时间序列预测方法之一,称为ARIMA。

     

    介绍

    时间序列提供了预测未来价值的机会。 基于以前的价值观,可以使用时间序列来预测经济,天气和能力规划的趋势,其中仅举几例。 时间序列数据的具体属性意味着通常需要专门的统计方法。

    在本教程中,我们将针对时间序列产生可靠的预测。 我们将首先介绍和讨论自相关,平稳性和季节性的概念,并继续应用最常用的时间序列预测方法之一,称为ARIMA。

    用于建模和预测时间序列未来点的Python中的一种方法被称为SARIMAX ,其代表具有eXogenous回归的季节性自动反馈集成移动平均值 。 在这里,我们将主要关注ARIMA组件,该组件用于适应时间序列数据,以更好地了解和预测时间序列中的未来点。

    时间序列预测——ARIMA(差分自回归移动平均模型)

    ARIMA(p,d,q)中,AR是"自回归",p为自回归项数;I为差分,d为使之成为平稳序列所做的差分次数(阶数);MA为"滑动平均",q为滑动平均项数,。ACF自相关系数能决定q的取值,PACF偏自相关系数能够决定q的取值。ARIMA原理:非平稳时间序列转化为平稳时间序列然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型

    基本解释:

    自回归模型(AR)

    • 描述当前值与历史值之间的关系用变量自身的历史时间数据对自身进行预测
    • 自回归模型必须满足平稳性的要求
    • 必须具有自相关性,自相关系数小于0.5不适用
    • p阶自回归过程的公式定义:

                                                  

                                           ,y t-i 为前几天的值

    PACF,偏自相关函数(决定p值),剔除了中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的干扰之后x(t-k)对x(t)影响的相关程度。

    移动平均模型(MA)

    • 移动平均模型关注的是自回归模型中的误差项的累加,移动平均法能有效地消除预测中的随机波动
    • q阶自回归过程的公式定义:

                                                  

    ACF,自相关函数(决定q值)反映了同一序列在不同时序的取值之间的相关性。x(t)同时还会受到中间k-1个随机变量x(t-1)、x(t-2)、……、x(t-k+1)的影响而这k-1个随机变量又都和x(t-k)具有相关关系,所 以自相关系数p(k)里实际掺杂了其他变量对x(t)与x(t-k)的影响

                                               

    ARIMA(p,d,q)阶数确定:

                               

                                                                                                                           截尾:落在置信区间内(95%的点都符合该规则)

                                  

                                                                                                 acf和pacf图

    平稳性要求:平稳性就是要求经由样本时间序列所得到的拟合曲线在未来的一段期间内仍能顺着现有的形态“惯性”地延续下去

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

        具体分为严平稳与弱平稳:
        严平稳:严平稳表示的分布不随时间的改变而改变。
        如:白噪声(正态),无论怎么取,都是期望为0,方差为1
        弱平稳:期望与相关系数(依赖性)不变
        未来某时刻的t的值Xt就要依赖于它的过去信息,所以需要依赖性

    因为实际生活中我们拿到的数据基本都是弱平稳数据,为了保证ARIMA模型的要求,我们需要对数据进行差分,以求数据变的平稳。

    模型评估:


    AIC:赤池信息准则(AkaikeInformation Criterion,AIC)
                                                             ???=2?−2ln(?)
    BIC:贝叶斯信息准则(Bayesian Information Criterion,BIC)
                                                             ???=????−2ln(?)
                                            k为模型参数个数,n为样本数量,L为似然函数

     

    模型残差检验:

    • ARIMA模型的残差是否是平均值为0且方差为常数的正态分布
    • QQ图:线性即正态分布

    先决条件

    本指南将介绍如何在本地桌面或远程服务器上进行时间序列分析。 使用大型数据集可能是内存密集型的,所以在这两种情况下,计算机将至少需要2GB的内存来执行本指南中的一些计算。

    要充分利用本教程,熟悉时间序列和统计信息可能会有所帮助。

    对于本教程,我们将使用Jupyter Notebook来处理数据。 如果您还没有,您应该遵循我们的教程安装和设置Jupyter Notebook for Python 3 。

    第1步 - 安装软件包

    为了建立我们的时间序列预测环境,我们先进入本地编程环境或基于服务器的编程环境:

    cd environments
    
    . my_env/bin/activate

    从这里,我们为我们的项目创建一个新的目录。 我们称之为ARIMA ,然后进入目录。 如果您将项目称为不同名称,请务必在整个指南中将您的名称替换为ARIMA

    mkdir ARIMA
    cd ARIMA
    

    本教程将需要warnings , itertools , itertools , numpy , matplotlibstatsmodels库。 warningsitertools库包含在标准Python库集中,因此您不需要安装它们。

    像其他Python包一样,我们可以用pip安装这些要求。 
    我们现在可以安装pandas , statsmodels和数据绘图包matplotlib 。 它们的依赖也将被安装:

    pip install pandas numpy statsmodels matplotlib

    在这一点上,我们现在设置为开始使用已安装的软件包。

    第2步 - 导入包并加载数据

    要开始使用我们的数据,我们将启动Jupyter Notebook:

    jupyter notebook

    要创建新的笔记本文件,请从右上角的下拉菜单中选择新建 > Python 3 :

    创建一个新的Python 3笔记本

    这将打开一个笔记本。

    最好的做法是,从笔记本电脑的顶部导入需要的库:

    import warnings
    import itertools
    import pandas as pd
    import numpy as np
    import statsmodels.api as sm
    import matplotlib.pyplot as plt
    plt.style.use('fivethirtyeight')

    我们还为我们的地块定义了一个matplotlib风格 。

    我们将使用一个名为“来自美国夏威夷Mauna Loa天文台的连续空气样本的大气二氧化碳”的数据集,该数据集从1958年3月至2001年12月期间收集了二氧化碳样本。我们可以提供如下数据:

    data = sm.datasets.co2.load_pandas()
    y = data.data

    让我们稍后再进行一些数据处理。 每周数据可能很棘手,因为它是一个很短的时间,所以让我们使用每月平均值。 我们将使用resample函数进行转换。 为了简单起见,我们还可以使用fillna()函数来确保我们的时间序列中没有缺少值。

    # The 'MS' string groups the data in buckets by start of the month
    y = y['co2'].resample('MS').mean()
    
    # The term bfill means that we use the value before filling in missing values
    y = y.fillna(y.bfill())
    
    print(y)
    Outputco2
    1958-03-01  316.100000
    1958-04-01  317.200000
    1958-05-01  317.433333
    ...
    2001-11-01  369.375000
    2001-12-01  371.020000

    我们来探索这个时间序列e作为数据可视化:

    y.plot(figsize=(15, 6))
    plt.show()

    图1:二氧化碳浓度时间序列

    当我们绘制数据时,会出现一些可区分的模式。 时间序列具有明显的季节性格局,总体呈上升趋势。

    要了解有关时间序列预处理的更多信息,请参阅“ 使用Python 3进行时间序列可视化的指南 ”,其中上面的步骤将更详细地描述。

    现在我们已经转换和探索了我们的数据,接下来我们继续使用ARIMA进行时间序列预测。

    第3步 - ARIMA时间序列模型

    时间序列预测中最常用的方法之一就是被称为ARIMA模型,它代表了A utoreg R essive综合M oving A版本 。 ARIMA是可以适应时间序列数据的模型,以便更好地了解或预测系列中的未来点。

    有三个不同的整数( p , d , q )用于参数化ARIMA模型。 因此,ARIMA模型用符号ARIMA(p, d, q) 。 这三个参数共计数据集中的季节性,趋势和噪音:

    • p是模型的自回归部分。 它允许我们将过去价值观的影响纳入我们的模型。 直观地说,这将是类似的,表示如果过去3天已经变暖,明天可能会变暖。
    • d是模型的集成部分。 这包括模型中包含差异量(即从当前值减去的过去时间点的数量)以适用于时间序列的术语。 直观地说,这将类似于说如果过去三天的温度差异非常小,明天可能会有相同的温度。
    • q是模型的移动平均部分。 这允许我们将模型的误差设置为过去以前时间点观察到的误差值的线性组合。

    在处理季节性影响时,我们利用季节性 ARIMA,表示为ARIMA(p,d,q)(P,D,Q)s 。 这里, (p, d, q)是上述非季节性参数,而(P, D, Q)遵循相同的定义,但适用于时间序列的季节分量。 术语s是时间序列的周期(季度为4 ,年度为12 ,等等)。

    由于所涉及的多个调整参数,季节性ARIMA方法可能会令人望而生畏。 在下一节中,我们将介绍如何自动化识别季节性ARIMA时间序列模型的最优参数集的过程。

    第4步 - ARIMA时间序列模型的参数选择

    当考虑使用季节性ARIMA模型拟合时间序列数据时,我们的第一个目标是找到优化感兴趣度量的ARIMA(p,d,q)(P,D,Q)s的值。 实现这一目标有许多指导方针和最佳实践,但ARIMA模型的正确参数化可能是一个需要领域专长和时间的艰苦的手工过程。 其他统计编程语言(如R提供了自动化的方法来解决这个问题 ,但尚未被移植到Python中。 在本节中,我们将通过编写Python代码来编程选择ARIMA(p,d,q)(P,D,Q)s时间序列模型的最优参数值来解决此问题。

    非季节性参数:p,d,q

    季节参数:P、D、Q

    s:时间序列的周期,年周期s=12

    我们将使用“网格搜索”来迭代地探索参数的不同组合。 对于参数的每个组合,我们使用statsmodels模块的SARIMAX()函数拟合一个新的季节性ARIMA模型,并评估其整体质量。 一旦我们探索了参数的整个范围,我们的最佳参数集将是我们感兴趣的标准产生最佳性能的参数。 我们开始生成我们希望评估的各种参数组合:

    # Define the p, d and q parameters to take any value between 0 and 2
    p = d = q = range(0, 2)
    
    # Generate all different combinations of p, q and q triplets
    pdq = list(itertools.product(p, d, q))
    
    # Generate all different combinations of seasonal p, q and q triplets
    seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
    
    print('Examples of parameter combinations for Seasonal ARIMA...')
    print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
    print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
    print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
    print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))
    Output
    Examples of parameter combinations for Seasonal ARIMA...
    SARIMAX: (0, 0, 1) x (0, 0, 1, 12)
    SARIMAX: (0, 0, 1) x (0, 1, 0, 12)
    SARIMAX: (0, 1, 0) x (0, 1, 1, 12)
    SARIMAX: (0, 1, 0) x (1, 0, 0, 12)

    我们现在可以使用上面定义的参数三元组来自动化不同组合对ARIMA模型进行培训和评估的过程。 在统计和机器学习中,这个过程被称为模型选择的网格搜索(或超参数优化)。

    在评估和比较配备不同参数的统计模型时,可以根据数据的适合性或准确预测未来数据点的能力,对每个参数进行排序。 我们将使用AIC (Akaike信息标准)值,该值通过使用statsmodels安装的ARIMA型号方便地返回。 AIC衡量模型如何适应数据,同时考虑到模型的整体复杂性。 在使用大量功能的情况下,适合数据的模型将被赋予比使用较少特征以获得相同的适合度的模型更大的AIC得分。 因此,我们有兴趣找到产生最低AIC的模型。

    下面的代码块通过参数的组合来迭代,并使用SARIMAX函数来适应相应的季节性ARIMA模型。 这里, order参数指定(p, d, q)参数,而seasonal_order参数指定季节性ARIMA模型的(P, D, Q, S)季节分量。 在安装每个SARIMAX()模型后,代码打印出其各自的AIC得分。

    warnings.filterwarnings("ignore") # specify to ignore warning messages
    
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(y,
                                                order=param,
                                                seasonal_order=param_seasonal,
                                                enforce_stationarity=False,
                                                enforce_invertibility=False)
    
                results = mod.fit()
    
                print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            except:
                continue

    由于某些参数组合可能导致数字错误指定,因此我们明确禁用警告消息,以避免警告消息过载。 这些错误指定也可能导致错误并引发异常,因此我们确保捕获这些异常并忽略导致这些问题的参数组合。

    上面的代码应该产生以下结果,这可能需要一些时间:

    Output
    SARIMAX(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.3436240402125
    SARIMAX(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.711172764114
    SARIMAX(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320026
    SARIMAX(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878315690562
    SARIMAX(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.6578978064144
    SARIMAX(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.7647912940095
    ...
    ...
    ...
    SARIMAX(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647112294245
    SARIMAX(1, 1, 1)x(1, 0, 1, 12) - AIC:327.9049123596742
    SARIMAX(1, 1, 1)x(1, 1, 0, 12) - AIC:444.12436865161305
    SARIMAX(1, 1, 1)x(1, 1, 1, 12) - AIC:277.7801413828764

    我们的代码的输出表明, SARIMAX(1, 1, 1)x(1, 1, 1, 12)产生最低的AIC值为277.78。 因此,我们认为这是我们考虑过的所有模型中的最佳选择。

    第5步 - 安装ARIMA时间序列模型

    使用网格搜索,我们已经确定了为我们的时间序列数据生成最佳拟合模型的参数集。 我们可以更深入地分析这个特定的模型。

    我们首先将最佳参数值插入到新的SARIMAX模型中:

    mod = sm.tsa.statespace.SARIMAX(y,
                                    order=(1, 1, 1),
                                    seasonal_order=(1, 1, 1, 12),
                                    enforce_stationarity=False,
                                    enforce_invertibility=False)
    
    results = mod.fit()
    
    print(results.summary().tables[1]) #详细输出,results.summary()可以输出全部的模型计算参数表
    Output
    ==============================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
    ------------------------------------------------------------------------------
    ar.L1          0.3182      0.092      3.443      0.001       0.137       0.499
    ma.L1         -0.6255      0.077     -8.165      0.000      -0.776      -0.475
    ar.S.L12       0.0010      0.001      1.732      0.083      -0.000       0.002
    ma.S.L12      -0.8769      0.026    -33.811      0.000      -0.928      -0.826
    sigma2         0.0972      0.004     22.634      0.000       0.089       0.106
    ==============================================================================

    由SARIMAX的输出产生的SARIMAX返回大量的信息,但是我们将注意力集中在系数表上 coef显示每个特征的重量(即重要性)以及每个特征如何影响时间序列 P>|z| 列通知我们每个特征重量的意义。 这里,每个重量的p值都低于或接近0.05 ,所以在我们的模型中保留所有权重是合理的

    在 fit 季节性ARIMA模型(以及任何其他模型)的情况下,运行模型诊断是非常重要的,以确保没有违反模型的假设 plot_diagnostics对象允许我们快速生成模型诊断并调查任何异常行为。

    results.plot_diagnostics(figsize=(15, 12))
    plt.show()

    图2:模型诊断

    我们的主要关切是确保我们的模型的残差是不相关的,并且平均分布为零。 如果季节性ARIMA模型不能满足这些特性,这是一个很好的迹象,可以进一步改善

    残差在数理统计中是指实际观察值与估计值(拟合值)之间的差。“残差”蕴含了有关模型基本假设的重要信息。如果回归模型正确的话, 我们可以将残差看作误差的观测值。 它应符合模型的假设条件,且具有误差的一些性质。利用残差所提供的信息,来考察模型假设的合理性及数据的可靠性称为残差分析

    在这种情况下,我们的模型诊断表明,模型残差正常分布如下:

    • 在右上图中,我们看到红色KDE线与N(0,1)行(其中N(0,1) )是正态分布的标准符号,平均值0 ,标准偏差为1 ) 。 这是残留物正常分布的良好指示。

    • 左下角的qq图显示,残差(蓝点)的有序分布遵循采用N(0, 1)的标准正态分布采样的线性趋势。 同样,这是残留物正常分布的强烈指示。

    • 随着时间的推移(左上图)的残差不会显示任何明显的季节性,似乎是白噪声。 这通过右下角的自相关(即相关图)来证实,这表明时间序列残差与其本身的滞后值具有低相关性

    这些观察结果使我们得出结论,我们的模型选择了令人满意,可以帮助我们了解我们的时间序列数据和预测未来价值

    虽然我们有一个令人满意的结果,我们的季节性ARIMA模型的一些参数可以改变,以改善我们的模型拟合。 例如,我们的网格搜索只考虑了一组受限制的参数组合,所以如果我们拓宽网格搜索,我们可能会找到更好的模型

    第6步 - 验证预测

    我们已经获得了我们时间序列的模型,现在可以用来产生预测。 我们首先将预测值与时间序列的实际值进行比较,这将有助于我们了解我们的预测的准确性。 get_prediction()conf_int()属性允许我们获得时间序列预测的值和相关的置信区间。

    pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)#预测值
    pred_ci = pred.conf_int()#置信区间

    上述规定需要从1998年1月开始进行预测。

    dynamic=False参数确保我们每个点的预测将使用之前的所有历史观测。

    我们可以绘制二氧化碳时间序列的实际值和预测值,以评估我们做得如何。 注意我们如何在时间序列的末尾放大日期索引。

    ax = y['1990':].plot(label='observed')
    pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)
    
    ax.fill_between(pred_ci.index,
                    pred_ci.iloc[:, 0],
                    pred_ci.iloc[:, 1], color='k', alpha=.2)
    #图形填充fill fill_between,参考网址:
    #https://www.cnblogs.com/gengyi/p/9416845.html
    
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
    plt.legend()
    
    plt.show()

     pred.predicted_mean 可以得到预测均值,就是置信区间的上下加和除以2

    函数间区域填充 fill_between用法:

    plt.fill_between( x, y1, y2=0, where=None, interpolate=False, step=None, hold=None,data=None,  **kwarg)

    x - array( length N) 定义曲线的 x 坐标

    y1 - array( length N ) or scalar 定义第一条曲线的 y 坐标

    y2 - array( length N )  or scalar 定义第二条曲线的 y 坐标

    where - array of bool (length N), optional, default: None 排除一些(垂直)区域被填充。注:我理解的垂直区域,但帮助文档上写的是horizontal regions

    也可简单地描述为:

    plt.fill_between(x,y1,y2,where=条件表达式, color=颜色,alpha=透明度)" where = " 可以省略,直接写条件表达式 

    图3:二氧化碳浓度静态预测

    总体而言,我们的预测与真实价值保持一致,呈现总体增长趋势

    量化预测的准确性

    这是很有必要的。 我们将使用MSE(均方误差,也就是方差),它总结了我们的预测的平均误差。 对于每个预测值,我们计算其到真实值的距离并对结果求平方。 结果需要平方,以便当我们计算总体平均值时,正/负差异不会相互抵消。

    y_forecasted = pred.predicted_mean
    y_truth = y['1998-01-01':]
    
    # Compute the mean square error
    mse = ((y_forecasted - y_truth) ** 2).mean()
    print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
    Output
    The Mean Squared Error of our forecasts is 0.07

    我们前进一步预测的MSE值为0.07 ,这是接近0的非常低的值。MSE=0是预测情况将是完美的精度预测参数的观测值,但是在理想的场景下通常不可能。

    然而,使用动态预测(dynamic=True)可以获得更好地表达我们的真实预测能力。 在这种情况下,我们只使用时间序列中的信息到某一点,之后,使用先前预测时间点的值生成预测。

    在下面的代码块中,我们指定从1998年1月起开始计算动态预测和置信区间。

    pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
    pred_dynamic_ci = pred_dynamic.conf_int()

    绘制时间序列的观测值和预测值,我们看到即使使用动态预测,总体预测也是准确的。 所有预测值(红线)与地面真相(蓝线)相当吻合,并且在我们预测的置信区间内。

    ax = y['1990':].plot(label='observed', figsize=(20, 15))
    pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)
    
    ax.fill_between(pred_dynamic_ci.index,
                    pred_dynamic_ci.iloc[:, 0],
                    pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)
    
    ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1998-01-01'), y.index[-1],
                     alpha=.1, zorder=-1)
    
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
    
    plt.legend()
    plt.show()

    图4:二氧化碳浓度动态预测

    我们再次通过计算MSE量化我们预测的预测性能

    # Extract the predicted and true values of our time-series
    y_forecasted = pred_dynamic.predicted_mean
    y_truth = y['1998-01-01':]
    
    # Compute the mean square error
    mse = ((y_forecasted - y_truth) ** 2).mean()
    print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
    OutputThe Mean Squared Error of our forecasts is 1.01

    从动态预测获得的预测值产生1.01的MSE。 这比前进一步略高,这是意料之中的,因为我们依赖于时间序列中较少的历史数据。

    前进一步和动态预测都证实了这个时间序列模型是有效的。 然而,关于时间序列预测的大部分兴趣是能够及时预测未来价值观。

    第7步 - 生成和可视化预测

    在本教程的最后一步,我们将介绍如何利用季节性ARIMA时间序列模型来预测未来的价值。 我们的时间序列对象的

    get_forecast()属性可以计算预先指定数量的步骤的预测值。

    # Get forecast 500 steps ahead in future
    pred_uc = results.get_forecast(steps=500)
    
    # Get confidence intervals of forecasts
    pred_ci = pred_uc.conf_int()

    我们可以使用此代码的输出绘制其未来值的时间序列和预测。

    ax = y.plot(label='observed', figsize=(20, 15))
    pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
    ax.fill_between(pred_ci.index,
                    pred_ci.iloc[:, 0],
                    pred_ci.iloc[:, 1], color='k', alpha=.25)
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
    
    plt.legend()
    plt.show()

    图5:时间序列和未来价值预测

    现在可以使用我们生成的预测和相关的置信区间来进一步了解时间序列并预见预期。 我们的预测显示,时间序列预计将继续稳步增长。

    随着我们对未来的进一步预测,我们对自己的价值观信心愈来愈自然。 这反映在我们的模型产生的置信区间,随着我们进一步走向未来,这个模型越来越大。

    结论

    在本教程中,我们描述了如何在Python中实现季节性ARIMA模型。 我们广泛使用了pandasstatsmodels图书馆,并展示了如何运行模型诊断,以及如何产生二氧化碳时间序列的预测。

    这里还有一些其他可以尝试的事情:

    • 更改动态预测的开始日期,以了解其如何影响预测的整体质量。
    • 尝试更多的参数组合,看看是否可以提高模型的适合度。
    • 选择不同的指标以选择最佳模型。 例如,我们使用AIC测量来找到最佳模型,但是您可以寻求优化采样均方误差

    补充:其中一些好用的方法,代码如下:

     

    filename_ts = 'data/series1.csv'
    ts_df = pd.read_csv(filename_ts, index_col=0, parse_dates=[0])
    
    n_sample = ts_df.shape[0]
    print(ts_df.shape)
    print(ts_df.head())

    # Create a training sample and testing sample before analyzing the series
    
    n_train=int(0.95*n_sample)+1
    n_forecast=n_sample-n_train
    #ts_df
    ts_train = ts_df.iloc[:n_train]['value']
    ts_test = ts_df.iloc[n_train:]['value']
    print(ts_train.shape)
    print(ts_test.shape)
    print("Training Series:", "\n", ts_train.tail(), "\n")
    print("Testing Series:", "\n", ts_test.head())

     

    def tsplot(y, lags=None, title='', figsize=(14, 8)):
        
        fig = plt.figure(figsize=figsize)
        layout = (2, 2)
        ts_ax   = plt.subplot2grid(layout, (0, 0))
        hist_ax = plt.subplot2grid(layout, (0, 1))
        acf_ax  = plt.subplot2grid(layout, (1, 0))
        pacf_ax = plt.subplot2grid(layout, (1, 1))
        
        y.plot(ax=ts_ax) # 折线图
        ts_ax.set_title(title)
        y.plot(ax=hist_ax, kind='hist', bins=25) #直方图
        hist_ax.set_title('Histogram')
        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax) # ACF自相关系数
        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax) # 偏自相关系数
        [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]
        sns.despine()
        fig.tight_layout()
        return ts_ax, acf_ax, pacf_ax
    
    tsplot(ts_train, title='A Given Training Series', lags=20);

     

                               

    # 此处运用BIC(贝叶斯信息准则)进行模型参数选择
    # 另外还可以利用AIC(赤池信息准则),视具体情况而定
    import itertools
    
    p_min = 0
    d_min = 0
    q_min = 0
    p_max = 4
    d_max = 0
    q_max = 4
    
    # Initialize a DataFrame to store the results
    results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min,p_max+1)],
                               columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])
    
    for p,d,q in itertools.product(range(p_min,p_max+1),
                                   range(d_min,d_max+1),
                                   range(q_min,q_max+1)):
        if p==0 and d==0 and q==0:
            results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
            continue
        
        try:
            model = sm.tsa.SARIMAX(ts_train, order=(p, d, q),
                                   #enforce_stationarity=False,
                                   #enforce_invertibility=False,
                                  )
            results = model.fit() ##下面可以显示具体的参数结果表
    ## print(model_results.summary())
    ## print(model_results.summary().tables[1])
    # http://www.statsmodels.org/stable/tsa.html
            # print("results.bic",results.bic)
            # print("results.aic",results.aic)
            
            results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.bic
        except:
            continue
    results_bic = results_bic[results_bic.columns].astype(float)

    results_bic如下所示:

    为了便于观察,下面对上表进行可视化:、

    fig, ax = plt.subplots(figsize=(10, 8))
    ax = sns.heatmap(results_bic,
                     mask=results_bic.isnull(),
                     ax=ax,
                     annot=True,
                     fmt='.2f',
                     );
    ax.set_title('BIC');
    //annot
    //annotate的缩写,annot默认为False,当annot为True时,在heatmap中每个方格写入数据
    //annot_kws,当annot为True时,可设置各个参数,包括大小,颜色,加粗,斜体字等
    # annot_kws={'size':9,'weight':'bold', 'color':'blue'}
    #具体查看:https://blog.csdn.net/m0_38103546/article/details/79935671

                                              

    # Alternative model selection method, limited to only searching AR and MA parameters
    
    train_results = sm.tsa.arma_order_select_ic(ts_train, ic=['aic', 'bic'], trend='nc', max_ar=4, max_ma=4)
    
    print('AIC', train_results.aic_min_order)
    print('BIC', train_results.bic_min_order)

    注:

    SARIMA总代码如下:

    
    """
    https://www.howtoing.com/a-guide-to-time-series-forecasting-with-arima-in-python-3/
    http://www.statsmodels.org/stable/datasets/index.html
    """
     
    import warnings
    import itertools
    import pandas as pd
    import numpy as np
    import statsmodels.api as sm
    import matplotlib.pyplot as plt
    plt.style.use('fivethirtyeight')
     
    '''
    1-Load Data
    Most sm.datasets hold convenient representations of the data in the attributes endog and exog.
    If the dataset does not have a clear interpretation of what should be an endog and exog, 
    then you can always access the data or raw_data attributes.
    https://docs.scipy.org/doc/numpy/reference/generated/numpy.recarray.html
    http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.resample.html
    # resample('MS') groups the data in buckets by start of the month,
    # after that we got only one value for each month that is the mean of the month
    # fillna() fills NA/NaN values with specified method
    # 'bfill' method use Next valid observation to fill gap
    # If the value for June is NaN while that for July is not, we adopt the same value
    # as in July for that in June
    '''
     
    data = sm.datasets.co2.load_pandas()
    y = data.data  # DataFrame with attributes y.columns & y.index (DatetimeIndex)
    print(y)
    names = data.names  # tuple
    raw = data.raw_data  # float64 np.recarray
     
    y = y['co2'].resample('MS').mean()
    print(y)
     
    y = y.fillna(y.bfill()) # y = y.fillna(method='bfill')
    print(y)
     
    y.plot(figsize=(15,6))
    plt.show()
     
    '''
    2-ARIMA Parameter Seletion
    ARIMA(p,d,q)(P,D,Q)s
    non-seasonal parameters: p,d,q
    seasonal parameters: P,D,Q
    s: period of time series, s=12 for annual period
    Grid Search to find the best combination of parameters
    Use AIC value to judge models, choose the parameter combination whose
    AIC value is the smallest
    https://docs.python.org/3/library/itertools.html
    http://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html
    '''
     
    # Define the p, d and q parameters to take any value between 0 and 2
    p = d = q = range(0, 2)
     
    # Generate all different combinations of p, q and q triplets
    pdq = list(itertools.product(p, d, q))
     
    # Generate all different combinations of seasonal p, q and q triplets
    seasonal_pdq = [(x[0], x[1], x[2], 12) for x in pdq]
     
    print('Examples of parameter combinations for Seasonal ARIMA...')
    print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
    print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
    print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
    print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))
     
    warnings.filterwarnings("ignore") # specify to ignore warning messages
     
    for param in pdq:
        for param_seasonal in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(y,
                                                order=param,
                                                seasonal_order=param_seasonal,
                                                enforce_stationarity=False,
                                                enforce_invertibility=False)
     
                results = mod.fit()
     
                print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            except:
                continue
     
     
    '''
    3-Optimal Model Analysis
    Use the best parameter combination to construct ARIMA model
    Here we use ARIMA(1,1,1)(1,1,1)12 
    the output coef represents the importance of each feature
    mod.fit() returnType: MLEResults 
    http://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.mlemodel.MLEResults.html#statsmodels.tsa.statespace.mlemodel.MLEResults
    Use plot_diagnostics() to check if parameters are against the model hypothesis
    model residuals must not be correlated
    '''
     
    mod = sm.tsa.statespace.SARIMAX(y,
                                    order=(1, 1, 1),
                                    seasonal_order=(1, 1, 1, 12),
                                    enforce_stationarity=False,
                                    enforce_invertibility=False)
     
    results = mod.fit()
    print(results.summary().tables[1])
     
    results.plot_diagnostics(figsize=(15, 12))
    plt.show()
     
    '''
    4-Make Predictions
    get_prediction(..., dynamic=False) 
    Prediction of each point will use all historic observations prior to it 
    http://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.mlemodel.MLEResults.get_prediction.html#statsmodels.regression.recursive_ls.MLEResults.get_prediction
    http://pandas.pydata.org/pandas-docs/stable/generated/pandas.Series.plot.html
    https://matplotlib.org/api/_as_gen/matplotlib.axes.Axes.fill_between.html
    '''
     
    pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)
    pred_ci = pred.conf_int() # return the confidence interval of fitted parameters
     
    # plot real values and predicted values
    # pred.predicted_mean is a pandas series
    ax = y['1990':].plot(label='observed')  # ax is a matplotlib.axes.Axes
    pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)
     
    # fill_between(x,y,z) fills the area between two horizontal curves defined by (x,y)
    # and (x,z). And alpha refers to the alpha transparencies 
    ax.fill_between(pred_ci.index,
                    pred_ci.iloc[:, 0],
                    pred_ci.iloc[:, 1], color='k', alpha=.2)
     
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
    plt.legend()
     
    plt.show()
     
    # Evaluation of model
    y_forecasted = pred.predicted_mean
    y_truth = y['1998-01-01':]
     
    # Compute the mean square error
    mse = ((y_forecasted - y_truth) ** 2).mean()
    print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
     
    '''
    5-Dynamic Prediction
    get_prediction(..., dynamic=True)
    Prediction of each point will use all historic observations prior to 'start' and
    all predicted values prior to the point to predict  
    '''
     
    pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
    pred_dynamic_ci = pred_dynamic.conf_int()
     
    ax = y['1990':].plot(label='observed', figsize=(20, 15))
    pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)
     
    ax.fill_between(pred_dynamic_ci.index,
                    pred_dynamic_ci.iloc[:, 0],
                    pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)
     
    ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1998-01-01'), y.index[-1],
                     alpha=.1, zorder=-1)
     
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
     
    plt.legend()
    plt.show()
     
    # Extract the predicted and true values of our time-series
    y_forecasted = pred_dynamic.predicted_mean
    y_truth = y['1998-01-01':]
     
    # Compute the mean square error
    mse = ((y_forecasted - y_truth) ** 2).mean()
    print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
     
    '''
    6-Visualize Prediction
    In-sample forecast: forecasting for an observation that was part of the data sample;
    Out-of-sample forecast: forecasting for an observation that was not part of the data sample.
    '''
     
    # Get forecast 500 steps ahead in future
    # 'steps': If an integer, the number of steps to forecast from the end of the sample.
    pred_uc = results.get_forecast(steps=500)  # retun out-of-sample forecast 
     
    # Get confidence intervals of forecasts
    pred_ci = pred_uc.conf_int()
     
    ax = y.plot(label='observed', figsize=(20, 15))
    pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
    ax.fill_between(pred_ci.index,
                    pred_ci.iloc[:, 0],
                    pred_ci.iloc[:, 1], color='k', alpha=.25)
    ax.set_xlabel('Date')
    ax.set_ylabel('CO2 Levels')
     
    plt.legend()
    plt.show()
    

     

    展开全文
  • 通常把去除季节性因素的过程叫做季节性调整。季节调整的意义由于生产活动中一般都要提前安排来年或者下个月甚至下周的生产采购和销售计划。所以对当前整体经济活动的趋势判断就显得至关重要。对于年度序列的数据而言...

    时间序列分解的含义是,将时间序列分解成若干个因子叠加的共同影响。一般认为时间序列可以分解成三个因素:趋势-循环因素、季节性因素和白噪音。按照组合方式的不同,可以分为加法模型、乘法模型和对数加法模型。通常把去除季节性因素的过程叫做季节性调整。

    季节调整的意义

    由于生产活动中一般都要提前安排来年或者下个月甚至下周的生产采购和销售计划。所以对当前整体经济活动的趋势判断就显得至关重要。对于年度序列的数据而言,由于反映的是一个完整年份的经济运行状况,因此只要比较历史数据一般就可以得到每年的趋势变化情况。可是要得到一整年的数据情况才能判断趋势的话,却往往会出现落差很大或者超出预期的情况。比如提前观察到下个月经济开始收缩,但是企业反而根据去年制订的计划加快生产,显然就不合时宜。

    正所谓计划没有变化快,所以就有了每个季度或者每月的公布数据,有助于提前进行对计划的修正。但是使用季度或者月度数据就会遇到季节效应的问题。比如说夏季的时候空调销量很好,可是如果不能剔除掉季节性影响的话,单纯与上月相比,很容易判断成销售加快,趋势上升。所以为了克服季节性因素对数据的影响,因此就必须使用季节性调整办法来对数据进行修正。当然更理想化来说,我们希望得到趋势因素。

    国家统计局公布的数据一般是同比的。虽然说与去年同期比较,相对于未经季调的环比数据有优势。但是它不能及时判断出拐点。有研究表明,使用未经处理的同比数据对经济周期判断往往滞后于实际的转折点。而且我们国家还有个“特色”,很多数据都是1、2月合计公布的,统计局说这是为了避免1月份和2月份数据产生较大误差的处理方式。但实际上嘛,个人认为当所有人都明白1月和2月数据是受到季节性影响的时候,市场就不会对此有太大反应。反而容易让人觉得,统计局找这么个理由完全是为了省点工作。当然我没法改变统计局公布数据的方式和计划,但是实在要吐吐苦水,不公布单独的一月和二月数据对我们这些要做分析数据的研究人员实在有点为难。

    关于Matlab的季节性调整方法

    唠叨那么多,回归到正题。Matlab没有专门做季节调整的函数,这是让我觉得不解的地方。也许没有人会有这个需求,因为像Eviews这些软件都有相对应的工具,按个键就好了,非常方便。但是在使用这些季节性调整工具软件的时候,很容易就会犯错误。因为对各个参数设定不了解,同时也不能按照中国的情况来进行模型修改,所以直接用来做季调跟自己说不清楚,也没法跟别人解释清楚。中国人民银行有自己开发的一套软件,叫做中国版的X-12-ARIMA,也出了一套书(文献[a])去说明这个软件。当然这些都要花大价钱买,对于广大的屌丝研究员而言,根本没有机会去用。

    在Masthworks的网站上可以找到季节性调整相关的概念,以及季节性调整的几个例子。其中比较详尽的是利用复合移动平均法做季调的例子。[b]

    这个例子基本涵盖了其余与季调相关的过程。如果对季调有所了解的朋友,应该能看出来这个例子是不考虑对极端情况以及节假日效应处理的估计,单纯连续使用移动平均法对数据进行季调,也就是经常用到的X-11方法。

    关于X-11做季调的流程步骤,可以参考下表(参考[1] P134)

    表1 X-11计算原型(乘法模型)

    使用MATLAB对CPI环比进行季节性调整

    准备数据

    本文选用未经季调过的CPI环比数据来做展示。这个数据开始于1995年的1月份,到2013年7月结束。把1994年12月31日定义为基期,基期指数设为100。将所有环比数据还原为CPI的指数序列,得到图1。

    图1 CPI走势图

    从图上看,CPI有明显的上升趋势。但是季节性因素却不明显,仅有在1997年至2003年这段时间才呈现明显的季节性规律,主要是由于上升趋势较为明显。

    1、使用2×12移动平均估计趋势-循环成分

    由上图来看,使用乘法模型对该数据建模较合适。为了去除趋势成分,参考Matlab给出的例子,使用S2×13的中心对称移动平均对其进行平滑处理以估计趋势-循环的部分。但是中心化的移动平均容易丢失原数据序列两端的数据信息。

    这段程序清楚展示了如何利用S2×13中心对称移动平均进行平滑。由于平滑处理后的数据首尾两端的信息会被丢失,此处的简单处理办法是使首尾6个数据等于邻近的数据值。

    %% 1、使用2×12移动平均估计趋势-循环成分:

    sW13 = [1/24;repmat(1/12,11,1);1/24];

    CPIs = conv(CPI,sW13,'same');

    CPIs(1:6) = CPIs(7);

    CPIs(N-5:N) = CPIs(N-6);

    可以将平滑处理后的数据与原始数据作比较,在图2上用红线表示的是平滑处理后的数据。可以看到首尾是水平无变化的。

    图2 初步平滑后的序列与原始序列对比

    2、估计季节-不规则成分:

    用原始数据除以平滑处理后的数据就得到去除趋势成分的CPI数据,它包含了周期成分和不规则的成分(白噪音)。在这个基础上可以进一步分离出周期成分,也就是季节性的成分。这只是第一步,并非是最后得到的季节成分。后面还要反复运用移动平均来做调整。

    %% 2、估计季节-不规则成分:xt =

    CPI./CPIs;

    3、对每个月份应用3×3移动平均估计预备季节成分

    Matlab要输入S3×3阶的移动平均算子必须手动输入。

    % S3x3 季调滤子

    % 中心对称系数sW3 = [1/9;2/9;1/3;2/9;1/9];

    % 非中心对称系数(用于调整首尾6个数据)

    aW3 = [.259 .407;.37 .407;.259 .185;.111 0];

    aW3是Musgrave(1964a,1964b)给出的用于估计季节因子的复合移动平均相对应的非对称滤子,但是这在X-11方法中并未被采用。对于这些滤子的合理性,尚未看到有任何出版物对其作出解释。[a]

    但是Matlab里面却调用了这些滤子,在实际运用中需要当心。

    % 使用2×12移动平均估计预备季节成分sW13 =

    [1/24;repmat(1/12,11,1);1/24];

    sb = conv(shat,sW13,'same');

    sb(1:6) = sb(s+1:s+6);

    sb(N-5:N) = sb(N-s-5:N-s);

    % 然后标准化

    s33 = shat./sb;

     

    图3 预备季调因子

    4、估计季节调整后序列

    得到预备季调因子以后,便可得到第一阶段季调后的序列

    dt = CPI./s33;

    5、用13项Henderson移动平均估计趋势-循环成分

    在X-11中,使用Henderson移动平均从季调后的序列中提取趋势-循环成分。本例中使用的是13项Henderson移动平均算子

    % Henderson 系数

    sWH = [-0.019;-0.028;0;.066;.147;.214;

    .24;.214;.147;.066;0;-0.028;-0.019];

    对应的非对称移动平均系数。但对比[a]中所给的系数略有不同。

    % 非对称系数(针对首位6个数据)

    aWH = [-.034 -.017 .045 .148 .279 .421;

    -.005 .051 .130 .215 .292 .353;

    .061 .135 .201 .241 .254 .244;

    .144 .205 .230 .216 .174 .120;

    .211 .233 .208 .149 .080 .012;

    .238 .210 .144 .068 .002 -.058;

    .213 .146 .066 .003 -.039 -.092;

    .147 .066 .004 -.025 -.042 0 ;

    .066 .003 -.020 -.016 0 0 ;

    .001 -.022 -.008 0 0 0 ;

    -.026 -.011 0 0 0 0 ;

    -.016 0 0 0 0 0 ];

    应用这些滤子便可以得到最终的趋势成分h13

    % h13就是最终的趋势成分

    first = 1:12;

    last = N-11:N;

    h13 = conv(dt,sWH,'same');

    h13(N-5:end) = conv2(dt(last),1,aWH,'valid');

    h13(1:6) = conv2(dt(first),1,rot90(aWH,2),'valid');

    对应的季节成分

    % 估计季节-不规则成分:xt =

    CPI./h13;

    图4 13项Henderson移动平均法得到的趋势因素

    6、对每个月份应用3×5移动平均估计最终季节成分

    % S3x5 季调滤子

    % 中心对称系数sW5 =

    [1/15;2/15;repmat(1/5,3,1);2/15;1/15];

    % 非中心对称系数aW5 = [.150 .250

    .293;

    .217 .250 .283;

    .217 .250 .283;

    .217 .183 .150;

    .133 .067 0;

    .067 0 0];

    ……

    % 运用2×12移动平均sW13 =

    [1/24;repmat(1/12,11,1);1/24];

    sb = conv(shat,sW13,'same');

    sb(1:6) = sb(s+1:s+6);

    sb(N-5:N) = sb(N-s-5:N-s);

    % 然后标准化得到最终季节成分s35 =

    shat./sb;

    Season = s35;

    7、估计季节调整后序列

    用原始数据除以季节因子就可以得到最终季调后的序列

    dt = CPI./s35;

    图5 季调后CPI走势

    8、估计最终趋势

    最终趋势便是前面的h13,这里给出各个成分序列的对比图。

    图6 CPI各成分序列的对比图

    9、估计最终不规则成分

    用季调后的序列除以趋势序列便得到不规则序列

    Irr = dt./h13;

    图7 CPI的不规则成分(白噪音)

    后记

    1、本文主要参考资料有二:

    a)中国人民银行调查统计司编著的《时间序列X-12-ARIMA季节调整——原理与方法》,其中对X-12-ARIMA做季调的程序做了非常详细的说明。至今统计局仍采用这一程序进行季节性调整。

    b)MATHWORKS主页上关于季调的文档。例子中做季调的对象是航空公司乘客统计的月度数据。此处的代码都基本出自于这个例子的源程序。

    2、本文使用的X-11方法一直饱受诟病,对它的批评集中在下面三点:缺乏统计理论支撑,有时候会扭曲各种构成成分,会扭曲变量之间的关系。后来Box和Jenkins(1970)有关ARIMA模型的成果对基于模型的季节调整方法做出了重大贡献,后期形成了X-11-ARIMA方法。后来再在前者基础上出现了X-12-ARIMA方法,基本思路都是在使用X-11进行季调前必须对序列进行建模,以弥补X-11方法的缺陷。

    3、往后会介绍ARIMA模型,以及如何结合ARIMA进行季调的方法。当然未来可能还会涉及到异常值处理以及节假日调整等问题。这个系统比较庞大,所以本文并不打算就此展开。后期将逐一介绍,请关注。

    展开全文
  • 本篇内容分为两部分:滞后算子和季节性模型,其中介绍前者既回顾了之前的内容,也为介绍后者做了铺垫。1 滞后算子时间序列模型中的自回归项和移动平均项都可以使用滞后算子(Lag Operator...

    本篇内容分为两部分:滞后算子和季节性模型,其中介绍前者既回顾了之前的内容,也为介绍后者做了铺垫。

    1 滞后算子

    时间序列模型中的自回归项和移动平均项都可以使用滞后算子(Lag Operator)进行方便地表示。如可以被表示为。

    • 对于常数而言,它的任意滞后项都还是它本身,即。

    • 分配律

    滞后算子满足分配律:。

    • 结合律

    滞后算子满足结合律:

    根据以上性质,对于ARMA(, ),有

    进而有

    上式是ARMA(, )模型的一个特解,它的展开式是一个MA()过程。

    滞后算子有如下几个常用的公式:

    当时,

    当时,

    当时,

    2 季节性模型

    2.1 示例数据

    许多时间序列数据都带有一定的季节性。USAccDeaths数据集是R中系统自带的数据,它储存的是美国1973-1978年逐月的意外死亡人口数。

    data <- c(USAccDeaths)
    plot(USAccDeaths, type = "l")
    points(USAccDeaths, pch = 21, col = "red", bg = "white")
    6907e460940b7d22934b53e5afcd1a8e.png

    从图中可以看出该序列存在明显的周期性。下面再绘制出它的ACF和PACF图象:

    acf(data, lag.max = 48)
    pacf(data, lag.max = 48)
    881009d323347b2d7e387e571b5fa84e.png 208f5becddc01b47453b9c80215cd4ca.png

    从ACF图象中可以看出,相差6六个月的月份之间的意外死亡人口数呈现明显负相关(),而相差12个月则呈现明显正相关(),这也表明数据序列的周期性可能与季节因子有关。

    由于许多时间序列的周期性都与季节因子相关,因此周期性和季节性两个概念在时间序列分析中常常混用。

    2.2 普通季节模型

    考虑如下两个模型:

    model.1:

    model.2:

    对于model.1,它的自相关系数的特点是:当为周期的整数倍时,,其余为0,也就是呈现出间断指数递减的态势。假定、,使用arima.sim()函数模拟该序列数据,并绘制序列走势图和ACF图象:

    y1 <- arima.sim(model = list(order = c(4,0,0),
                                 ar = c(0,0,0,0.5)),
                    n = 200)
    plot(y1, type = "l")
    acf(y1)
    720efac9cc6dfaaeed86e2ba91809e2e.png ac518d4ccedaa0bd0906f73888dd9984.png

    对于model.2,它的自相关系数在时存在一个峰值,而在其他处均为0。同样,假定、模拟该序列数据:

    y2 <- arima.sim(model = list(order = c(0,0,4),
                                 ma = c(0,0,0,0.5)),
                    n = 200)
    plot(y2, type = "l")
    acf(y2)
    af94459e5513f626b514f1383349c729.png 04a389111057ef01d5da3699809e39a6.png

    从数据集USAccDeaths所对应的ACF图象上看,除了在12的整数倍呈现峰值外,在单个周期内自相关系数也呈现递减趋势,而非直接截尾为0,这表明该序列数据在季节性之外还存在非季节性。现实中的数据也往往不会如model.1model.2所刻画的那么理想。

    可以使用以下模型尝试对USAccDeaths数据集进行拟合:

    model.3:
    model.4:

    由于上述模型中有许多滞后值的系数为0,可以使用fixed参数进行设置:

    (model.3 <- arima(data, order = c(12,0,0),
                      fixed = c(NA, rep(0,10), NA, NA)))
    ## arima(x = data, order = c(12, 0, 0), fixed = c(NA, rep(0, 10), NA, NA))
    ## 
    ## Coefficients:
    ##          ar1  ar2  ar3  ar4  ar5  ar6  ar7  ar8  ar9  ar10  ar11    ar12
    ##       0.3442    0    0    0    0    0    0    0    0     0     0  0.6300
    ## s.e.  0.0648    0    0    0    0    0    0    0    0     0     0  0.0662
    ##       intercept
    ##       9176.4060
    ## s.e.   685.9749
    ## 
    ## sigma^2 estimated as 214051:  log likelihood = -548.73,  aic = 1105.47
    
    (model.4 <- arima(USAccDeaths, order = c(1,0,12),
                      fixed = c(NA, NA, rep(0,10), NA, NA)))
    ## arima(x = USAccDeaths, order = c(1, 0, 12), fixed = c(NA, NA, rep(0, 10), NA, 
    ##     NA))
    ## 
    ## Coefficients:
    ##          ar1     ma1  ma2  ma3  ma4  ma5  ma6  ma7  ma8  ma9  ma10  ma11
    ##       0.7241  0.0028    0    0    0    0    0    0    0    0     0     0
    ## s.e.  0.0946  0.1432    0    0    0    0    0    0    0    0     0     0
    ##         ma12  intercept
    ##       0.7203  8953.6706
    ## s.e.  0.1365   328.2761
    ## 
    ## sigma^2 estimated as 235875:  log likelihood = -552.24,  aic = 1114.47

    从AIC和模型简洁性的角度上看,model.3都优于model.4。下面再比较二者的预测值与实际值的差距:

    fit.3 <- data - model.3$residuals
    fit.4 <- data - model.4$residuals
    
    plot(data, type = "l")
    lines(1:72, fit.3, col = "red") # model.3
    lines(1:72, fit.4, col = "blue") # model.4
    c19de2e82d86b1dbdc53540fd3e2235d.png

    2.3 乘法季节模型

    假设示例数据中的非季节性由AR(1)过程控制,该过程使用滞后算子可以书写成;假设季节性由model.1所示的过程控制,其中,该过程使用滞后算子可以书写成。

    同理,若示例数据中的非季节性和季节性分别由MA过程控制,那么可以分别写作和。

    在上节所示的普通季节模型中,其思想是将控制季节性和非季节性的过程进行相加。如果要考虑二者的相互影响,可以使用乘法季节模型,如:

    model.5:

    或者,

    model.6:

    上述两式展开后可以得到ARMA模型的一般形式,并且其阶数高于周期12,但其需要估计的系数并没有增加。

    实际运用中并没有必要将上述两式展开,它们可以写成ARIMA(, , )(, , )的形式,其中、、为非季节性过程的参数,、、为非季节性过程的参数,为周期。前面已经提过,对于ARMA模型,、。具体地,model.5可以写作ARIMA(1,0,0)(1,0,0),model.6可以写作ARIMA(1,0,1)(0,0,1)。

    上述形式可以直接使用arima()函数进行估计:

    (model.5 <- arima(USAccDeaths, order = c(1,0,0),
                      seasonal = list(order = c(1,0,0), period = 12)))
    ## arima(x = USAccDeaths, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), 
    ##     period = 12))
    ## 
    ## Coefficients:
    ##          ar1    sar1  intercept
    ##       0.7579  0.8501  9214.7102
    ## s.e.  0.0768  0.0498   676.1701
    ## 
    ## sigma^2 estimated as 127203:  log likelihood = -533.45,  aic = 1074.89
    
    (model.6 <- arima(USAccDeaths, order = c(1,0,1),
                      seasonal = list(order = c(0,0,1), period = 12)))
    ## arima(x = USAccDeaths, order = c(1, 0, 1), seasonal = list(order = c(0, 0, 1), 
    ##     period = 12))
    ## 
    ## Coefficients:
    ##          ar1     ma1    sma1  intercept
    ##       0.6849  0.0907  0.7060  8947.0658
    ## s.e.  0.1065  0.1332  0.1287   311.2709
    ## 
    ## sigma^2 estimated as 235882:  log likelihood = -552,  aic = 1114.01

    从结果上看,model.5要优于model.6。下面对比model.5和普通季节模型中的优胜者model.3之间的预测值差异:

    fit.5 <- data - model.5$residuals
    
    plot(data, type = "l")
    lines(1:72, fit.3, col = "red", lty = 2) # model.3
    lines(1:72, fit.5, col = "red") # model.5
    0367ae4b028dd3d3e257490e8a65cb0d.png
    • 上图中,黑色实线为真实数据,红色虚线为使用普通季节模型(model.3)预测的结果,红色实线表示使用乘法季节模型(model.5)预测的结果。可以看出后者的拟合度更高。

    展开全文
  • 概念 时间序列(Time Series)  时间序列是均匀时间间隔上的观测值序列 ... 时间写按照季节性来分类,分为季节性时间序列和非季节性时间序列 非季节性时间序列:趋势部分、不规则部分; 季节性时间...
  • 1. 这本书对Python的知识点的描述很详细,而且排版看的很舒服. 2. 几个例题:假装自己从零开始学,将一些有代表、有意思的例题抽取出来. 3. 还有自己对一部分课后复习题,全部课后上机实践题的解题思路
  • 20210208-光大证券-财政政策专题研究系列之二:详解财政收支活动的架构和季节性特点.pdf
  • 季节性对打算交易农产品或其他有较明显季节性规律的其他资产的投资者来说是一个重要概念。在进行交易决策时,这些市场在一年中的特定时期出现的比较固定的供需规律是做出投资决策的关键考虑因素...
  • 有很多兼职亚马逊人做季节性产品,但是很多都没有什么效果,为什么? 因为你的季节性产品和别人的一样,可能是固有思维让你把路走窄了,所以你选择的商品,别人都在卖。 然而,“季节性”并不一定意味着特定的季节...
  • 季节性积水的河流叫季节河,也叫时令河。好比时令病就是指一定季节流行的病。  3、理论时区和现实时区:某一种标准时适用的地区范围叫时区,分理论时区和现实时区,理论时区是按经度把全球分成24区,每区跨经度15度...
  • 最近一直在接触时间序列,所以打算写一些有关时间序列的文章,预测部分会从规则开始、到传统模型、到机器学习、再到深度学习,此外也会介绍一些时间序列的基本概念,包括自相关、平稳、滞后季节...
  • 差异和独特 。能再同类产品中脱颖而出,吸引用户的注意力,捕获用户的心。   6.2 内功效用: 指满足用户基本需求,引导用户期望需求和创造用户兴奋需求的强大内功。即 核心价值,期望价值,附加价值的完美 ...
  • 旅游流概念的研究的探讨旅游流概念的研究的探讨摘 要:文章在比较分析各类旅游流定义研究的基础上,确定了旅游流的定义,区别了旅游流与旅游者个体综合;论述了旅游流的特征和影响因素,特别提出了旅游流的周期、...
  • 机器学习之概念偏移

    千次阅读 2017-12-17 15:32:55
    简单介绍了概念偏移,做了一些笔记
  • 【时序列】时序列数据如何一步步分解成趋势(trend)季节性(seasonality)和误差(residual)- 详细理解python sm.tsa.seasonal_decompose 在做时序列分析的时候,好多教程都告诉你要把时序列分解成趋势,季节性,...
  • 概念:对过去的观察值得加权平均值进行预测的一种方法,适用于水平历史数据 一次指数平滑法:Ft+1 =aYt+(1-a)Ft Ft表示t时预测值,Yt表示t时观察值。t取1时,F1=Y1。a为平滑系数,介于0到1之间。 最终的式子展开为 ...
  • 机器学习和深度学习概念入门

    万次阅读 多人点赞 2017-06-03 11:27:28
     对于很多初入学习人工智能的学习者来说,对机器学习、深度学习、人工智能的概念和区别还不是很了解,那么接下来就给大家从概念和特点上进行阐述。先看下三者的关系。  人工智能包括了机器学习,机器学习包括...
  • 来源:数据STUDIO 作者:云朵君大家对时间序列知多少?何为时间序列、时间序列分析、时间序列分解、时间序列预测,以及时间序列预测都有哪些方法?????点击关注|设为星标|干货速递????随着...
  • 额叶脑电图(EEG)的alpha不对称在情绪、动机和精神病理学研究中得到了广泛的研究,但它是一个使用多种程序进行量化和分析的指标,而程序的多样使交叉研究的解释变得模糊。本文的目的是提供脑电图alpha不对称的记录...
  • 再保险系统涉及的概念

    千次阅读 2018-12-20 10:51:24
    1.1再保险基本概念 1.再保险:再保险(reinsurance)也称分保,是保险人在原保险合同的基础上,通过签订分保合同,将其所承保的部分风险和责任向其他保险人进行保险的行为。 2.分入公司:就是再保公司 3.分出公司:...
  • 同样,随机梯度下降法的很多变体都有很高的可能能够找到接近严格凸函数最小值的点(但并非一定能找到)。 两个凸函数的和(例如 L2 损失函数 + L1 正则化)也是凸函数。 深度模型绝不会是凸函数。值得注意的是,...
  • 概念作为风靡全球的经典教材,看到大家对它如此的热情,我无比的感动,可见到各位学习它的方法,作为已经“背”完新概念全四册的我,又深表忧虑。 两年前,我也是新概念的忠实追随者,买齐全套后,恨不得一个月...
  • OpenERP库存管理的若干概念讲解(新增库存价值) « 于: 六月 03, 2011, 09:01:50 下午 » 一、复式库存(Double-Entry Stock Management)和库存移动(Stock Move) OpenERP的库存管理采取了独特的复式...
  • 时间序列(一)基本概念

    千次阅读 2020-07-05 10:16:52
    时间序列定义基本概念时间序列类型长期趋势 T季节趋势 S*循环变动 C*不规则变动 I时间序列的类型组合判断方法: 定义 时间序列也称动态序列,是指将某种现象的指标数值按照时间顺序排列而成的数值序列。时间序列分析...
  • AI:人工智能概念之《Google发布机器学习术语表 (中英对照)》——持续更新ML、DL相关概念2018年4月! 相关文章AI:人工智能概念之《Google发布机器学习术语表 (中英对照)》——持续更新ML、DL相关概念2018年4月...
  • 【深度强化学习】强化学习的基本概念

    千次阅读 多人点赞 2020-06-21 18:31:01
    文章目录前言第一章:强化学习的基本概念学习——监督, 无监督与强化学习强化学习的体系与联系Reward 奖励AgentEnvironmentActionsObservation马尔科夫决策过程马尔科夫链马尔科夫奖励过程马尔科夫决策过程Policy...
  • 在写这篇大数据文章之前,我发现身边很多IT人对于这些热门的新技术、新趋势往往趋之若鹜却又很难说的透彻,如果你问他大数据是什么,什么是大数据概念?估计很少能说出一二三来。究其原因,一是因为大家对大数据这类...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 7,824
精华内容 3,129
关键字:

季节性的概念