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 ， matplotlib和statsmodels库。 warnings和itertools库包含在标准Python库集中，因此您不需要安装它们。

像其他Python包一样，我们可以用pip安装这些要求。 

我们现在可以安装pandas ， statsmodels和数据绘图包matplotlib 。 它们的依赖也将被安装：

pip install pandas numpy statsmodels matplotlib

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

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

要开始使用我们的数据，我们将启动Jupyter Notebook：

jupyter notebook

要创建新的笔记本文件，请从右上角的下拉菜单中选择新建 > 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()



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

要了解有关时间序列预处理的更多信息，请参阅“ 使用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()



我们的主要关切是确保我们的模型的残差是不相关的，并且平均分布为零。 如果季节性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 = " 可以省略，直接写条件表达式 



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

量化预测的准确性

这是很有必要的。 我们将使用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()



我们再次通过计算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()



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

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

结论

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

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

更改动态预测的开始日期，以了解其如何影响预测的整体质量。
	
尝试更多的参数组合，看看是否可以提高模型的适合度。
	
选择不同的指标以选择最佳模型。 例如，我们使用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()


 

SARIMAX (seasonal autoregressive integrated moving average with exogenous regressor)是一种常见的时间序列预测方法，可以分为趋势部分和周期性部分；每个部分又可以分为自回归、差分和平滑部分。
趋势稳定性检测：Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test
null-hypothesis: 时间序列趋势稳定。significance: 0.05. 选择KPSS和不是Dickey Fuller由于KPSS的null-hypothesis是趋势稳定，所以接受条件相对于Dickey Fuller更宽松，差分阶数更少。

$y_t = \beta^{\prime} D_t + \mu_t + u_t$

$\mu_t = \mu_{t-1} + \epsilon_t,$

$\epsilon_t \sim WN(0, \sigma^2_\epsilon)$
其中D为确定时间序列趋势/常数。u为随机漫步。epsilon为残差。
LM statistic:

$KPSS = (T^-2 \sum_{t=1}^{T} \hat S_t^2) / \hat\lambda^2$

其中$\hat S_t = \sum_{j=1}^{t} \hat u_j$

$\hat{u}$ 是yt 拟合Dt的残差，$ \hat \lambda^2$ 是var(ut)的预估值。问题归结为拉格朗日乘数法证明 $\sigma^2_\epsilon = 0$
季节稳定性检测：Canova-Hansen方法

null-hypothesis: 时间序列季节稳定性。significance: 0.05
测试频率：给定待测最大频率m以下所有 2pi/m整数倍。
时间序列yi的拟合：

$y_i = \mu + x^{\prime}_i \beta + S_i + e_i$

$S_i = \sum_{j=1}^{m/2} f^{\prime}_{ji} \gamma_j$

其中

$f^{\prime}_{ji} = (cos((j/q) \pi i),  sin((j/q) \pi i))$
LM statistic:

$LM = \frac{1}{n^2} trace( (\hat \Omega ^f)^-1 \sum_{i=1}^{n} \hat F_i \hat F^{\prime}_i)$
其中

$\hat F_i = \sum_{t=1}^{i}f_t \hat e_t$

$\hat \Omega^f = \sum_{k=-m}^{m} w (\frac{k}{m}) \frac{1}{n} \sum_{i} f_{i+k} \hat e_{i+k} f^{\prime}_i \hat e_i$
用根据待检测频率m计算f’ji向量：用sample data示例代码如下：
import numpy as np
import pandas as pd

def seasonalDummy(tsArray, frequency):
    n = len(tsArray)
    m = frequency
    #if m == 1: tsArray.reshape([n, m])
    tt = np.arange(1, n + 1, 1)
    mat = np.zeros([n, 2 * m], dtype = float)
    for i in np.arange(0, m):
        mat[:, 2 * i] = np.cos(2.0 * np.pi * (i + 1) * tt / m)
        mat[:, 2 * i + 1] = np.sin(2.0 * np.pi * (i + 1) *tt / m)
    return mat[:, 0:(m-1)]

tsArray = np.array([  4.00195672,   4.99944175,   5.99861146,   7.00000213,
         8.00124207,   9.00059539,  10.00029542,  10.99871969,
        11.99728933,  12.99965231,  14.00148869,  15.0006378 ,
        16.00117112,  17.00159081,  17.99848509,  18.99957668,
        20.00066721,  20.99932292,  21.99992471,  23.00099164,
        24.00127222,  25.00014385,  26.00014191,  27.00037435,
        27.9985619 ,  28.99949718,  29.99844772,  30.99911627,
        31.99965086,  33.00211019,  34.00240288,  34.99889972,
        36.00240406,  37.0002379 ,  37.99958145,  38.99825111,
        39.99932529,  40.9998374 ,  42.00034236,  43.00206289,
        43.9994545 ,  45.00141283,  46.00041818,  47.00132581,
        48.00216031,  48.99812424,  50.00060522,  51.00049892,
        51.99817633,  52.9997362 ])
frequency = 6
pd.DataFrame(seasonalDummy(tsArray, frequency)).head(2)

Canova-Hansen用sample data示例代码如下：
from numpy.linalg import lstsq as lsq

n = len(tsArray)
frec = np.ones(int((frequency + 1) / 2), dtype = int)
ltrunc = int(np.round(frequency * np.power(n / 100.0, 0.25)))

$y_i = \mu + x^{\prime}_i \beta + S_i + e_i$

$S_i = \sum_{j=1}^{m/2} f^{\prime}_{ji} \gamma_j$
其中$f^{\prime}_{ji} = (cos((j/q) \pi i),  sin((j/q) \pi i))$
# create dummy column f'ji
r1 = seasonalDummy(tsArray, frequency)
#create intercept column for regression
r1wInterceptCol = np.column_stack([np.ones(r1.shape[0], dtype = float), r1])

# residual ei:
result = lsq(a = r1wInterceptCol, b = tsArray)
residual = tsArray - np.matmul(r1wInterceptCol, result[0])

long-run covariance matrix：$\Omega =  lim_{n \to \infty}\frac{1}{n}E(F_n F_n^{\prime})$

在ei可能有serial correlation的情况下可以用estimate

$\hat{\Omega} = \sum_{k=-m}^{m}w(\frac{k}{m})\frac{1}{n}\sum_{i}F_{i+k}F_i^{\prime}$
fhat = np.zeros([n, frequency - 1], dtype = float)
fhataux = np.zeros([n, frequency - 1], dtype = float)

for i in np.arange(0, frequency - 1):
    fhataux[:, i] = r1[:, i] * residual

for i in np.arange(0, n):
    for j in np.arange(0, frequency - 1):
        mySum = sum(fhataux[0:(i + 1), j])
        fhat[i, j] = mySum

$w(\cdot)$ 是kernel function，来自rob j. Hyndman的forecast包:
wnw = np.ones(ltrunc, dtype = float) - np.arange(1, ltrunc + 1, 1) / (ltrunc + 1)

$\hat{\Omega}^f$的计算：
Ne = fhataux.shape[0]
omnw = np.zeros([fhataux.shape[1], fhataux.shape[1]], dtype = float)
for k in np.arange(0, ltrunc):
    omnw = omnw + np.matmul(fhataux.T[:, (k+1):Ne], fhataux[0:(Ne-(k+1)), :]) * float(wnw[k])
cross = np.matmul(fhataux.T, fhataux)
omnwplusTranspose = omnw + omnw.T
omfhat = (cross + omnwplusTranspose) / float(Ne)

Generalized Hannan’s model：通过设置矩阵A选择需要测试的频率的子集。在程序中用的A = eye:
sq = np.arange(0, frequency - 1, 2)
frecob = np.zeros(frequency - 1, dtype = int)    
for i in np.arange(0, len(frec)):
    if (frec[i] == 1) & (i + 1 == int(frequency / 2.0)):
        frecob[sq[i]] = 1
    if (frec[i] == 1) & (i + 1 < int(frequency / 2.0)):
        frecob[sq[i]] = 1
        frecob[sq[i] + 1] = 1

a = frecob.tolist().count(1)  # find nr of 1's in frecob
A = np.zeros([frequency - 1, a], dtype = float)
j = 0
for i in np.arange(0, frequency - 1):
    if frecob[i] == 1:
        A[i, j] = 1
        j = j + 1

LM statistic的计算 （Nyblom(1989), Hansen(1990)）：当LM statistic值超过对应自由度的critical value时，拒绝 （null hypothesis = 没有单位根）：

$LM statistic = \frac{1}{n^2} trace((A&#x27; \hat{\Omega}^f A)^{-1} A&#x27; \sum_{i=1}^{n}\hat{F_i} \hat{F_i}^{\prime}A)$
其中$\sum_{i=1}^{n}\hat{F_i}$ 和$\sum_{i=1}^{n} \hat{F_i}^{\prime}$由前面步骤中的 $\hat{f}$ 给出:
from numpy.linalg import svd

aTomfhat = np.matmul(A.T, omfhat)
tmp = np.matmul(aTomfhat, A)

machineDoubleEps = 2.220446e-16

problems = min(svd(tmp)[1]) < machineDoubleEps # svd[1] are the singular values
if problems:
    stL = 0.0
else:
    solved = np.linalg.solve(tmp, np.eye(tmp.shape[1], dtype = float))
    step1 = np.matmul(solved, A.T)
    step2 = np.matmul(step1, fhat.T)
    step3 = np.matmul(step2, fhat)
    step4 = np.matmul(step3, A)
    stL = (1.0 / np.power(n, 2.0)) * sum(np.diag(step4))

在Canova-Hansen test接受alternative hypothesis的情况下对时间序列进行lag = 待测频率的差分（季节差分）。
季节性检测：

在进行季节性时间序列稳定性检测之前，首先判断a.时间序列是否有季节性，和b.时间序列在什么频率上有季节性。结果会作为时间序列稳定性检测的参数输入。
季节性检测根据离散傅里叶变换和自相关函数的“与”关系得出结论（只有两个method都返回真值，才会判定时间序列有季节性）。
离散傅里叶变换吧时间序列从时域变为频域。变换后频域的新序列为：
$X_k = \sum_{n=0}^{N-1}x_n \cdot [cos(2 \pi kn/N) - i \cdot sin(2 \pi kn/N)]$
在待检测频率上如果能量为最大值，则返回真值。
自相关函数检测最大lag = 待检频率各阶上的correlation coefficient。时间点t和s之间的自相关性R的定义为：
$R(s, t) = \frac{E[(X_t - \mu_t)(X_s - \mu_s)}{\sigma_t \sigma_s}$
如果在待检频率上的相关系数超过双边confidence interval在0.05的临界值
clim = norm.ppf((1 + ci) / 2) / np.sqrt(n)

则method返回真值。
DFT调用了numpy.fft.fft方法。ACF调用了statsmodels.tsa.stattools.acf方法。
测量函数：根据不同情况采用不同模型测量方法。算法使用了Rob J. Hyndman的 MASE (mean absolute scaled error)。与其他测量方法的优劣对比。
定义：

R-squared:

$SS_{tot} = \sum_{i}(y_i - \bar{y})^2$

$SS_{res} = \sum_{i}(y_i - \hat{y}_i)^2$

$R^2 = 1 - \frac{SS_{res}}{SS_{tot}}$
RMSE:

$RMSE = \sqrt{\frac{\sum_{t=1}^{n} (\hat{y}_t - y_t)^2}{n}}$
MAE:

$MAE = \frac{\sum_{t=1}^{n} |\hat{y}_t - y_t|}{n}$
MAPE:

$MAPE = \frac{100}{n} \sum_{t=1}^{n} | \frac{\hat{y}_t - y_t}{y_t} |$
sMAPE:

$SMAPE = \frac{100}{n} \sum_{t=1}^{n} \frac{|\hat{y}_t - y_t|}{(|\hat{y}_t| + |y_t|)/2}$
MASE without seasonality:

$MASE = \frac{\sum_{t = 1}^{n} |\hat{y}_t - y_t|}{\frac{n}{n-1} \sum_{t=2}^{n} |y_t - y_{t-1}|}$
MASE with seasonality:

$MASE = \frac{\sum_{t = 1}^{n} |\hat{y}_t - y_t|}{\frac{n}{n-m} \sum_{t=m+1}^{n} |y_t - y_{t-m}|}$
趋势平滑：在SARIMA模型中引入时间序列的趋势作为exogenous regressor(X)，有几种算法可以选择：
Lowess (locally weighted scatterplot smoothing): 基于KNN的非参数拟合方法。代码调用了 statsmodels.api.nonparametric.lowess
RANSAC

Random sample consensus，一种robust regression方法，可以探测异常值并使拟合对于异常值的敏感度降低。代码调用sklearn.linear_model.RANSACRegressor
Weighted moving average:
指数平滑递归表达：

$W_n = (1-\alpha) * W_{n-1} + \alpha * y_n$

$W_0 = y_0$

$\alpha = 2 / (span + 1)$
调用了 pandas.ewma
关于脉冲响应：
如果有另外的（多维）exogenous regressor Xi 影响预测模型，比如类似离散脉冲波形的机会点数据：假设各个脉冲regressor之间是独立的，并不受时间序列本身的影响：可以用多元线性回归首先发现Xi 和时间序列趋势yhat的关系。
考虑到历史上时间序列对于脉冲输入的响应不同，在算法中会测试三种不同脉冲响应与时间序列的相关性，并挑选相关性最强的脉冲响应作为Xi向量。这三种分别是a. 原始脉冲，b. 经指数平滑的脉冲，c. 经指数平滑并累加的脉冲（cumulative)。
示例：
%matplotlib inline
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings('ignore')

def myEwma(x, histPeriod = 6, fcstPeriod = 18):    
    from pandas import ewma
    xfit = ewma(x, span = histPeriod)
    xpred = np.zeros(fcstPeriod)
    tmp = np.zeros(histPeriod + 1)
    tmp[:histPeriod] = x[-histPeriod:].copy()
    tmp[histPeriod] = xfit[-1]
    for i in np.arange(0, fcstPeriod):
        xpred[i] = ewma(tmp, span = histPeriod)[-1]
        tmp = shiftLeft(tmp)
        tmp[-1] = xpred[i]
    return xfit, xpred

def shiftLeft(_ar, step = 1):
    if step == 0:
        return _ar
    ar = _ar.copy()
    ar[0:-step] = ar[step:]
    ar[-step:] = 0
    return ar

original = np.array([3000,0,0,0,0,0,0,0,0,0,1000,0,0,0,0,0,0,0,0,1000,0,0,0,0,0,0,0,0,0,0])
pd.DataFrame(original).plot(title = 'original')
movingAverage = myEwma(original)[0]
pd.DataFrame(movingAverage).plot(title = 'moving average')
cumulative = movingAverage.cumsum()
pd.DataFrame(cumulative).plot(title = 'cumulative')

检测相加性：时间序列拆分成趋势，季节性和残差的方式有相加和相乘两种。
y = Trend + Seasonality + Residual, 或者
y = Trend * Seasonality * Residual
如果是相乘的情况，残差的分布是不稳定的。所以如果时间序列没有通过相加性检测，会对时间序列做对数处理，变为相加：
$ynew = log(y) = log(Trend * Seasonality * Residual) = log(Trend) + log(Seasonality) + log(Residual)$
自动搜索SARIMA参数空间 Auto ARIMA：搜索差分阶数，检测季节性的存在，并搜索能给出Akaike Information Criterion最小值的ARMA模型维度 ARMA(p, q, P, Q)
ARMA模型范式：
$X_t - \alpha_1 X_{t-1} - \dots - \alpha_{p^{\prime}} X_{t-p^{\prime}} = \epsilon_t + \theta_1 \epsilon_{t-1} + \dots + \theta_q \epsilon_{t-q}$
自动搜索参数空间返回的结果是p和q的个数。作为输入给SARIMAX核心算法
搜索参数空间的初始值来自Rob J. Hyndman的paper section3.2.
SARIMAX 带外部变量的季节性自回归差分平滑算法：预测是基于SARIMAX的state-space representation。 核心算法调用了statsmodels.api.tsa.statespace.SARIMAX.
总结
以上是SARIMAX各个功能模块的详细数学方法和编程实现。由于一些语言（包括python）的时间序列预测开源算法包里缺少比如稳定性检测、季节性差分和自动搜索参数空间的功能，需要自己依照数学公式实现。

序列分解



1、非季节性时间序列分解

移动平均MA(Moving Average)

①SAM（Simple Moving Average） 

简单移动平均，将时间序列上前n个数值做简单的算术平均。 

SMAn=(x1+x2+…xn)/n

②WMA(Weighted Moving Average) 

加权移动平均。基本思想，提升近期的数据、减弱远期数据对当前预测值的影响，使平滑值更贴近最近的变化趋势。 

用Wi来表示每一期的权重，加权移动平均的计算： 

WMAn=w1x1+w2x2+…+wnxn

R中用于移动平均的API 

install.packages(“TTR”) 

SAM(ts,n=10)

ts   时间序列数据
n    平移的时间间隔，默认值为10
WMA(ts,n=10,wts=1:n)

wts  权重的数组，默认为1:n
#install.packages('TTR')
library(TTR)

data <- read.csv("data1.csv", fileEncoding="UTF8")

plot(data$公司A, type='l')

data$SMA <- SMA(data$公司A, n=3)

lines(data$SMA)



plot(data$公司A, type='l')

data$WMA <- WMA(data$公司A, n=3, wts=1:3)

lines(data$WMA)





2、季节性时间序列分解

在一个时间序列中，若经过n个时间间隔后呈现出相似性，就说该序列具有以n为周期的周期性特征。 

分解为三个部分： 

①趋势部分 

②季节性部分 

③不规则部分 
R中用于季节性时间序列分解的API 

序列数据周期确定

freg<-spec.pgram(ts,taper=0, log=’no’, plot=FALSE)
start<-which(freqspec==max(freqspec))周期开始位置
frequency<-1/freqfreq[which(freqspec==max(freq$spec))]周期长度
序列数据分解 

decompose(ts)

data <- read.csv("data2.csv", fileEncoding = "UTF8")

freq <- spec.pgram(data$总销量, taper=0, log='no', plot=FALSE);

start <- which(freq$spec==max(freq$spec))
frequency <- 1/freq$freq[which(freq$spec==max(freq$spec))]


data$均值 <- data$总销量/data$分店数

freq <- spec.pgram(data$均值, taper=0, log='no', plot=FALSE);

start <- which(freq$spec==max(freq$spec))
frequency <- 1/freq$freq[which(freq$spec==max(freq$spec))]

plot(data$均值, type='l')

meanTS <- ts(
  data$均值[start:length(data$均值)], 
  frequency=frequency
) 
ts.plot(meanTS)

meanTSdecompose <- decompose(meanTS)
plot(meanTSdecompose)

#趋势分解
meanTSdecompose$trend
#季节性分解数据
meanTSdecompose$seasonal
#随机部分
meanTSdecompose$random

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import ExponentialSmoothing

data = pd.read_excel('时间序列预测数据集.xlsx')
# data.columns=[时间,投递人数,投递次数,工程师投递人数,工程师投递次数,招聘发布公司量,发布职位量,工程师岗位发布公司,工程师岗位发布量]

for i in data.drop('时间',axis=1).columns:
    y = data[i][:37]
    ets = ExponentialSmoothing(y, trend='add', seasonal='add', seasonal_periods=12)
    r = ets.fit()
    pred = r.predict(start=len(y), end=len(y) + 12)
    data[i][37:] = pred
    pd.DataFrame({
    'origin': data[i],
    'fitted': r.fittedvalues,
    'pred': pred
}).plot(legend=True)
    plt.title(i)
data


https://www.cnblogs.com/lfri/p/12243268.html