随时间变化python画图 随时间变化的数据

目录

• 时间序列分析

• 基本问题

• 时间序列的组成部分
• 时间序列的平稳性

• 平稳性
• 白噪声
• 随机行走
• 时间序列的零均值化和平稳化
• 时间序列的平稳性检验

• 预测准确度的衡量

• 统计量
• 验证步骤

• ARIMA模型

• 模型简介
• 识别估计与预测

• 识别
• 估计检测
• 预测

时间序列分析

历史数据往往以时间序列的形式呈现出来，这些数据时随着时间的变化而变化的，反映了事物、现象在时间上的发展变动情况，是相同事物或现象在不同时刻或时期所形成的数据，称之为时间序列数据，简称时间序列或时间数列。

前面颜色就的大部分数据都是反映若干事物或现象在同一时刻或时间上所处的状态或特征，或者反映其与时间无关的特征，这些数据反映了事物或现象之间存在的内在数值联系，称之为横截面数据。

基本问题

时间序列（time series）是将某一个变量或指标在不同时间上的不同数值按照时间的先后顺序排列而成的数列，也被称为时间数列。通常可以用 \(X_1,X_2,\cdots,X_t\)

时间序列的组成部分

由于受到各种偶然因素的影响，时间数列往往表现出某种随机性，彼此之间存在着统计上的依赖关系。

一个时间序列可以由4部分组成：长期趋势、季节变动、循环波动和不规则变动。

长期趋势：事物或现象在较长时间内持续发展变化的一种趋势或状态。

季节变动：事物或现象在一年内随季节更换形成的有规律变动。

循环波动：事物或现象周而复始的变动。循环变动不同于长期趋势或季节变动，它是无固定规律的交替波动。

不规则变动：无法一个上述组成部分解释或不可控的随机变动。

为了更加深入研究时间序列的规律，往往可以将一个时间序列用乘法模式或加法模式分解为上述的4各组成部分。

示例

import pandas as pd
import statsmodels.api as sm
import scipy.stats as stats
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# 收集了商场销售额度月度数据，进行时间序列分析

# 指定为黑体中文字体，防止中文乱码
plt.rcParams["font.sans-serif"] = ["Heiti TC"]
# 解决保存图像是负号'-'显示为方块的问题
plt.rcParams['axes.unicode_minus'] = False

sales = pd.read_sas('./data/sales_monthly.sas7bdat')
print(sales.head(5))
#        DATE       Month  Sales
# 0 2001-01-01  b'2001Jan'  814.0
# 1 2001-02-01  b'2001Feb'  774.8
# 2 2001-03-01  b'2001Mar'  782.8
# 3 2001-04-01  b'2001Apr'  772.0
# 4 2001-05-01  b'2001May'  817.6

# 数据趋势图
sales.index = pd.Index(pd.date_range('1/2001', '9/2008', freq='1M'))  # 设置日前索引
sales['Sales'].plot()
plt.show()
# 从图中可以看出，销售额度月度数据总体上呈直线上升趋势，但是在上升过程中还有上下波动情况。

时间序列的平稳性

按照不同的性质和特征，可以对时间序列进行分类。从统计特性上来看，试卷序列可以分为平稳时间序列和非平稳时间序列。

平稳性

如果一个时间序列的概率分布与时间 \(t\) 无关，则称该序列为严格的平稳时间序列（stationary time series）。如果时间序列的一、二阶矩存在，而且对任意时刻 \(t\)

严格平稳时间序列要求比较严格，通常情况下，如果不明确提出严格平稳，所谓的平稳即宽平稳，其特征即均值和协方差不随时间变化而变化。

宽平稳也叫弱平稳性（广义平稳时间序列）

$\[E(y_t) = \mu \\ Cov(y_s,y_t) = Cov(y_{s+z},y_{t+z})=\sigma^2 \]$

在平稳的保证下，对历史时序数据进行分析的参数估计结果也比较稳定，可以直接用于对未来时许数据的预测。此外，非平稳时间序列在分析时，还可能出现本来没有什么关系的变量之间出现伪回归的情况。因此，平稳性时合理进行时间序列分析和预测的重要保证。

白噪声

平稳时间序列有一种特殊情况，即分布不随时间变化而变化，其具有零均值和同方差性，且协方差为零，即白噪声。白噪声序列可以用于对时序模型拟合进行检验。

白噪声是联系横截面数据和纵向数据的纽带。严格来讲，白噪声是具有独立同分布的数据序列，即没有特定随时间变化特征的满足平稳性条件的数据。

白噪声数据在时间序列研究中之所以重要，是因为所有时间序列的技术都是要将一组数据通过一系列过程尽量变为一个白噪声数据，这一系列过程就被称之为滤子。

白噪声数据的特点是对其的点预测 将其方差不依赖于我们想要预测到多远，而只与样本数据的均值和方差有关 。如果我们有个白噪声过程 \(y_t, t=1,...,T\)，而我们要预测 \(T+s\) 期末数据的大小，则其最优期望值为样本均值 \(\bar{y}\)，而预测的 \(\alpha\)

$\[\bar{y}\pm t_{T-1,1-\alpha/2}\sqrt{(1+1/T)s_y} \]$

其中\(s_y\) 为样本方差根，而 \(t_{T-1,1-\alpha/2}\)则是自由度为\(T-1\) 的T-分布统计量 \(\alpha\)

随机行走

白噪声时间序列的累加和就构成一个随机行走时间序列。随机行走有如下特点：首先这种时间序列数据是非平稳的，其均值和方差都随着时间而变化。对随机行走时间序列取一阶差分作为滤子，过滤后的时间序列则为白噪声序列。

随机行走模型是一类非常重要的时间序列模型。因为其为对应的白噪声的时间序列累加和，所以每个时间点上该变量的期望和方差分别为

$\[E(y_t) = y_0 + t_\mu \\ Var(y_t) = t\sigma^2 \]$

其中，\(\mu,\sigma^2\) 是对应的白噪声序列的期望均值和方差，而 \(y_0\)

$\[y_{T+s} = y_T + s\hat{\mu} \pm 2\hat{\sigma}\sqrt{s} \]$

其中，\(y_T\) 是已知随机行走时间序列的末尾值，\(s\) 是要预测的未来时间间隔，\(\hat{\mu}, \hat{\sigma}\)

import numpy as np
import  matplotlib.pyplot  as plt

# 生成白噪声和随机行走
np.random.seed(1291)
z = np.random.normal(0.1, 2, 100)
y = np.cumsum(z)

fig, ax1 = plt.subplots()
plt.plot(z, label="White Noise")
plt.plot(y, label="Random Walk")
plt.legend()
plt.savefig('./data/TimeSeries/WhiteNoise.png')

plt.show()

mean1 = np.round(np.mean(y[:20]), 4);
mean2 = np.round(np.mean(y[-20:]), 4)
std1 = np.round(np.std(y[:20]), 4);
std2 = np.round(np.std(y[-20:]), 4)

print("前20个数据点的均值为%.4f, 标准差为%.4f" % (mean1, std1))
print("后20个数据点的均值为%.4f, 标准差为%.4f" % (mean2, std2))

plt.close()
plt.show()

时间序列的零均值化和平稳化

在日常生活中，社会经济现象的特征随着时间推移，大部分都会表现为上升或下降趋势的非平稳时间序列。因此可以考虑对时间序列进行变化，使非平稳序列转化为平稳序列。

零均值化：对均值不为零的时间序列进行转化，使其均值为零的数据转换过程。通常可用时间序列中的每一个数值 \(X_t\)

平稳化：对非平稳的时间序列进行转化，使之成为平稳时间序列的数据转换过程。通常可用每一个数值减去其前面的一个数值，即 \(X_t-X_{t-1}\) 差分的方法。差分方法还可以是每一个数值减去其前面的、任何间隔为 \(s\) 的一个数值，即\(X_t-X_{t-s}\)。

对原始数据进行一次差分的过程被称为一阶差分。如果数列经过一阶差分之后还是非平稳数列，则可进行二阶差分或高阶差分。在一般情况下，非平稳的时间序列在经过一阶差分或二阶差分之后都可以平稳化。

在有些情况下，还可以通过函数的形式进行零均值化和平稳化，如对时间序列的数值取对数后再进行差分。具体函数具体分析。

# 进行一阶差分，绘制时序图
sales['Sales'].diff(1).plot()
plt.show()

时间序列的平稳性检验

除了利用类似时间序列图和差分图及性能时间序列平稳性的粗略判断外，还可以领用样本自相关系数及其图形或统计检验的方法进行进一步判断

• 自相关系数和自相关图检验(ACF/PACF)

与相关系数类似，自相关系数实际上是构成时序的各个组成元素的相关系数，即通过历史数据和未来数据的相关性，可以得知不同时期数据之间的相关程度。其取值范围在-1到1之间，其绝对值越接近1，说明时间序列的自相关程度越高。

对于一个时序数据总体而言，在给定正整数 \(p\) 的情况下，可以考察 \(X_t,X_{t+p}\) 之间的相关系数 \(\rho_{p}\) 来度量时间间隔为 \(p\) 的两部分数据之间的相关性。因此，依据样本数据，可以定义时间序列的 \(p\)

$\[r_p = \frac{\sum_{t=1}^{n=p}{(X_t-\bar{X})(X_{t+p}-\bar{X})}}{\sum_{t=1}^{n}{(x_t-\bar{X})^2}},p=1,2,3,\cdots \]$

运用自相关图判定时间序列平稳性的一般准则是：若时间序列的自相关系数基本上（一般 \(p>3\)

除了绘制自相关系数图进行主管的平稳性检验之外，还可以进行自相关系数的白噪声检验。如果白噪声检验结果显著，则表明时间序列总体自相关显著，即表现为非平稳。当所有的白噪声检验结果均不显著时，则时间序列是平稳的。

示例

# 自相关系数和自相关图检验
from statsmodels.tsa.stattools import acf, pacf

ts_d1_ACF = pd.DataFrame(acf(sales['Sales'].diff(1).iloc[1:92]), columns=['ACF'])
ts_d1_ACF['PACF'] = pd.DataFrame(pacf(sales["Sales"].diff(1).iloc[1:92]))
res = ts_d1_ACF.head(10).T
print(res)
#         0         1         2         3  ...         6         7         8         9
# ACF   1.0 -0.213246 -0.267816  0.045357  ... -0.099197  0.131948 -0.020434 -0.047949
# PACF  1.0 -0.215615 -0.335942 -0.123429  ... -0.262224 -0.048721 -0.154343 -0.126597

# 自相关系数图
fig = plt.figure(figsize=(10, 6))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(sales['Sales'].diff(1).iloc[1:92], lags=24, ax=ax1)
# 注意要去掉一阶差分数据sales['Sales'].diff(1)的第一个空值再进行计算
plt.show()
# 时间序列的自相关系数基本都落入置信区间，其逐渐趋于0，故该数据经过一阶差分之后可以认为是平稳的。

# 白噪声检验
r, q, p = sm.tsa.acf(sales['Sales'].diff(1).iloc[1:92].values.squeeze(), qstat=True)
mat = np.c_[range(1, 41), r[1:], q, p]
table = pd.DataFrame(mat, columns=['lag', 'AC', 'Q', 'Prob(>Q)'])
LB_result = table.loc[[5, 11, 17]]
LB_result.set_index('lag', inplace=True)
print(LB_result)
#             AC          Q  Prob(>Q)
# lag                                
# 6.0  -0.099197  12.490298  0.051883
# 12.0  0.198715  18.940955  0.089963
# 18.0 -0.037085  21.678168  0.246587
# 在\alpha=0.05的条件下，白噪声检验的p值均大于0.05，表明白噪声检验不显著，即经过一阶差分的时间序列是平稳的。

• 单位根检验（ADF）

仅从图形描述来对时间序列平稳性进行判断的准确性毕竟有限，一般还可以考察了使用单位根检验的方法对时序数据的平稳性进行检验

一个时间序列如果能通过差分的方式平稳化，则可称其具有单位根，即当一个时间序列具有单位根时时非平稳的。其原假设和备择假设如下

\(H_0\)：时间序列具有单位根；\(H_1\)：时间序列是平稳序列（即没有单位根）

示例

# 单位根检验
from statsmodels.tsa.stattools import adfuller


def DFTest(sales, regression, maxlag, autolag='AIC'):
    # 便于月度DF检验输出结果，对输出内容做标注
    print('ADF-Test Result')
    dftest = adfuller(sales, regression=regression, maxlag=maxlag, autolag=autolag)
    # regression根据模型形式指定为：c(仅常数，默认)；ct(常数和长期趋势)；ctt(常数，线性和二次曲线趋势)；nc(无常数无趋势)
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', 'Lags Used', 'nobs'])
    for k, v in dftest[4].items():
        dfoutput['Critical Value at {}'.format(k)] = v
    print(dfoutput)


# 对原始数据进行单位根检验
DFTest(sales['Sales'], regression='nc', maxlag=6, autolag='AIC')
print(37 * '-')
# 对原始数据的一阶差分数据进行单位根检验
DFTest(sales['Sales'].diff(1).iloc[1:92], regression='nc', maxlag=5, autolag='AIC')
# ADF-Test Result
# Test Statistic            0.902719
# p-value                   0.901743
# Lags Used                 6.000000
# nobs                     85.000000
# Critical Value at 1%     -2.592546
# Critical Value at 5%     -1.944575
# Critical Value at 10%    -1.614030
# dtype: float64
# -------------------------------------
# ADF-Test Result
# Test Statistic          -6.460678e+00
# p-value                  9.494200e-10
# Lags Used                5.000000e+00
# nobs                     8.500000e+01
# Critical Value at 1%    -2.592546e+00
# Critical Value at 5%    -1.944575e+00
# Critical Value at 10%   -1.614030e+00
# dtype: float64
# 原始数据的p值比较大，无法拒绝原假设，即原始序列具有单位根，即非平稳序列
# 一阶差分后的序列p值非常小，可以拒绝原假设，即一阶差分后的序列不具有单位根，即平稳序列

• KPSS

KPSS检验是一种较新的检验方式，与ADF检验相反，其默认假设待验证的时间序列是平稳的，如果得到的统计量p值较大 ，则无法拒绝原假设，说明这个时间序列是平稳的；反之则是不平稳的。

# 在statsmodels的0.8以上版本才有
# 检查版本
print(statsmodels.version.full_version)

from statsmodels.tsa.stattools import kpss

预测准确度的衡量

统计量

时间序列模型通常用来对未来的值进行预测，衡量预测的值的准确性的统计量如下所示

• 平均误差

平均误差能较好地衡量现有模型是否有很好描述的线性趋势

$\[ME = \frac{1}{T_2}\sum_{t=T_1+1}^{T_1+T_2}{e_t} \]$

• 平均百分比误差

平均百分比误差也用于衡量是否有短期趋势没有被模型很好地描述，不过它是以相对误差的形式来表达的

$\[MPE =\frac{1}{T_2}\sum_{t=T_1+1}^{T_1+T_2}{\frac{e_t}{y_t}} \]$

• 均方差

均方差相对于平均误差来说，能侦测出线性趋势之外更多的没有被模型描述的数据模式，比如周期性等，因此更常用

$\[MSE = \frac{1}{T_2}\sum_{t=T_1+1}^{T_1+T_2}{e_t^2} \]$

• 平均绝对误差

平均绝对误差在衡量模型的准确度方面和均方差有类似的效果，只是对异常值相对来说稳健性更高

$\[MAE = \frac{1}{T_2}\sum_{t=T_1+1}^{T_1+T_2}{||e_t||} \]$

• 平稳绝对百分比误差

MAPE结合了MAE和MPE的优点，能较好地侦测线性趋势之外的更多的数据模式，并以相对误差的形式表达

验证步骤

• 将长度为 \(T=T_1+T_2\) 的样本时间序列分为两个字序列，其中前面一个 \((t=1,\cdots,T_1)\) 子序列用于模型训练，后面一个子序列 \((t=T_1+1,\cdots,T)\)
• 用第一个子序列训练一个待验证模型
• 使用上一步训练的模型，使用时间范围为 \(t=1,\cdots, T_1\) 的因变量来预测未来 \(T_1+1, \cdots, T\) 时间段的因变量值\(\hat{y_t}\)，即对用于模型验证部分的子序列因变量使用待验证模型进行拟合
• 使用上一步拟合的因变量值和对应的实际因变量值，计算单步预测误差 \(e_t = u_t-\hat{y_t}\)，然后采用一种或者多种衡量模型准确度的统计量来计算综合预测能力

可以对每一个待验证模型都执行以上2～4步骤，选取综合预测能力最好，即统计量值最小的那个待选模型。

ARIMA模型

时间序列分析的ARIMA建模过程也叫做Box-Jenkins方法，它主要是在对时间序列分析的基础上，通过选择适当的模型进行预测

模型简介

ARIMA模型也叫做整合自回归移动平均模型（autoregressive integrated moving-average model），其模型可分为自回归模型（AR模型）、移动平均模型（MA模型）和自回归移动平均模型（ARMA模型）。

Box-Jenkins方法的基本思想是用时间序列的过去值和现在值的线性组合来预测其未来的值。即将随时间推移而形成的系列数据视为一个随机序列，把时间序列作为一组仅依赖于时间 \(t\)

• AR模型

AR模型即自回归模型（autoregressive model），其具体表现为某个观测值 \(X_t\) 与其滞后 \(p\)

$\[X_t= \varphi_1X_{t-1}+\varphi_2X_{t-2}+\cdots+\varphi_pX_{t-p}+a_t \]$

其中，\(X_t\) 为零均值平稳序列，\(a_t\) 为随机误差项。常把模型简记为 \(AR(p)\)。

对于\(AR(p)\)模型而言，有其基本假设：

1. 假设 \(X_t\) 仅与 \(X_{t-1},X_{t-2},\cdots,X_{t-p}\)
2. 在 \(X_{t-1},X_{t-2},\cdots,X_{t-p}\) 已知的条件下，\(X_t\) 与 \(X_{t-p-1},X_{t-p-2},\cdots\)
3. \(a_t\)

• MA模型

MA模型即移动平均模型（moving-average model），其具体表现为某个观测值 \(X_t\) 与先前 \(t-1, t-2, t-q\) 个时刻进入系统的 \(q\) 个随机误差项即 \(a_t, a_{t-1},\cdots,a_{t-q}\)

$\[X_t= a_t-\theta_1a_{t-1}-\theta_2a_{t-2}\cdots-\theta_qa_{t-q} \]$

常把模型简记为 \(MA(q)\)。

对于\(MA(q)\)模型，\(X_t\) 仅与 \(a_t,a_{t-1},\cdots,a_{t-q}\) 有关，而与 \(a_{t-q-1},a_{t-q-2}，\cdots\) 无关，且 \(a_t\)

• ARMA模型

ARMA模型即自回归移动平均模型（auto-regressive moving average model），即观测值 \(X_t\) 不仅与其以前 \(p\) 个时刻的自身观测值有关，而且还与其以前时刻进入系统的 \(q\)

$\[X_t = \varphi_1X_{t-1}+\varphi_2X_{t-2}+\cdots+\varphi_pX_{t-p}+a_t + a_t-\theta_1a_{t-1}-\theta_2a_{t-2}\cdots-\theta_qa_{t-q} \]$

显然 \(ARMA(p,q)\) 模型便是 \(AR(p)\) 与 \(MA(q)\) 的组合，\(ARMA(p,0)\) 便是 \(AR(p)\) 模型，\(ARMA(0, q)\) 便是 \(MA(q)\)

在进行ARMA建模之前，分析的时间序列必须满足平稳性条件。非平稳的时间序列可通过差分法使之平稳化并进行平稳性检验。时间序列通过平稳化之后，便可建立ARMA模型进行分析，待模型进行参数估计之后，再通过数据交换的可逆性，使得模型参数估计结果适应平稳化之前的数据。通过这个过程建立的模型称之为整合的ARMA模型，即ARIMA模型。如果对原始数据进行了 \(d\) 次差分，用差分数据所建立的 \(ARMA(p,q)\) 可以记为 \(ARMA(p,d,q)\)。

识别估计与预测

建立ARMA模型的基本前提是保证时间序列的平稳性，ARIMA建模过程则是把非平稳时间序列平稳化，再建立ARMA模型。

• 建模基本流程

1.可视化待建模的序列数据
2.使用ADF/KPSS测试确定将数据平稳化所需的差分次数(d)
3.使用ACF/PACF确定移动平均对应的预测误差项(q)和自回归项个数(p)
4.ARIMA(p,d,q)拟合
5.对ARMA(p,q)模型进行参数估计即显著性检验
6.利用显著的模型对时间序列进行预测

• 优化流程

1.对于使用ACF/PACF确定移动平均对应的预测误差项和自回归项个数，一般从一项开始
2.对于拟合好的ARIMA模型，将预测误差项和自回归项分别减少一个再拟合
3.根据AIC或者BIC判断模型相对简单的AR或者MA模型是否有改进
4.对自回归和移动平均项个数递增一个，逐次检验

• 建模遵循一般原则

对于自相关性，一般可以通过增加自回归系项或者移动平均部分里面的预测误差项个数来消除。一般的原则是如果未消除的自相关是正自相关关系，即ACF图里面第一项是正值，则使用增加自回归项的方法较好；而如果未消除的自相关是负自相关关系，则增加预测误差项的方法更为合适。这是因为一般而言，差分方法对于消除正相关关系非常有效，但是同时也会额外引入反向的相关关系，这时候会出现过度差分的情况，需要额外引入一个预测误差项来消除负相关关系，这也是为什么在上面的建模步骤里预先引入预测误差项建模，而不是先引入自回归项开始建模，也就是先拟合一个ARIMA(0, 1, 1)模型再看看ARIMA(1,1,0)模型，通常ARIMA(0,1,1)模型会比ARIMA(1,1,0)模型拟合效果好一些。

识别差分项的原则

1.如果建模的序列的自相关系数一直衍生到很长的滞后项（比如10或者更多滞后项），则获得平稳序列所需的差分次数较多
2.如果滞后一项的自相关系数为0或者为负，或者所有的 自相关系数都很小，则该序列不需要更多的差分来获取平稳性。通常而言，如果滞后一项的自相关性为-0.5或更小，则很可能该序列被过度差分了，这是需要注意的。
3.最优的差分项个数通常对应于差分后拥有最小标准差的时间序列
4.如果元序列不需要进行差分，则假定原序列是平稳的。一阶差分则意味着原序列有一个为常数的平均趋势，二阶差分则意味着原序列有一个依时间变化的趋势。
5.对不需要进行差分的时间序列建模时通常包含一个常数项。如果对一个需要一阶差分的时间序列进行建模，则只有在该时间序列包含非0的平均趋势的时候才需要包含常数项。而对一个需要进行二阶差分的时间序列进行建模时则通常不用包含常数项。

识别自回归或者预测误差项的原则

1.如果差分后的序列的PACF显示为sharp Cutoff或者滞后一项的自相关为正，则说明该序列差分不足，这时候可以对模型增加一个或者多个自相关项，增加个数通常为PACF截断的地方
2.如果差分后的序列的ACF显示为急剧截断或者滞后一项的自相关为负自相关，则说明该序列差分过度，这时候可以对模型增加一个或者多个预测误差项，增加个数通常为ACF截断的地方
3.自回归项和预测误差项有可能会相互抵消，因此如果一个两种要素都包含的ARIMA模型对数据拟合的很好，则通常可以试一试一个少一个自回归项或者少一个预测误差项的模型。一般来说，同时包含多个自回归和多个预测误差项的ARIMA模型都会过度拟合
4.如果自回归项的系数和接近1，即自回归部分有单位根现象，那么这时候应该将自回归项减少一个，同时增加一次差分操作
5.如果预测误差项的系数和接近1，即移动平均部分有单位根现象，那么这时候应该将预测误差项减少一个，同时减少一个差分操作
6.自回归或者移动平均部分有单位根通常也表现为长期预测不稳定。

识别模型的季节性

1.如果一个时间序列有很强的季节性，则必须使用一次季节周期作为差分，否则模型会认为季节性会随着时间逐渐消除。但是使用季节周期做差分不能超过一次，如果使用了季节周期做差分，则非季节周期的差分最多也只能再进行一次。
2.如果一个适当差分之后的序列的自相关系数在第s个滞后上 仍然表现为正 ，而s为季节性周期包含的时间段数，则在模型里添加一个季节性自回归项。如果这个自相关系数为负，则添加一个季节性预测误差项。通常情况下，如果已经使用了季节周期做差分，则第二种情况更为常见，而第一种情况通常是还没使用季节性周期做差分。如果季节性周期很规律，则使用差分是比引入一个季节性自回归项更好的方法，分析师应该尽量避免在模型里同时引入季节性自回归和季节性预测误差项，否则模型会过度拟合，甚至在拟合过程本身会出现不收敛的情况。

识别

ARMA模型的识别主要是针对确定 \(p,q\)

• 利用自相关系数图、偏自相关系数图进行模型识别

自相关系数描述的是时间序列观测值与过去值之间的相关性；偏自相关系数（PACF）则为在给定中间观测值的条件下，观测值与前面某个间隔的观测值的相关系数。

截尾 是指在自相关关系图或偏自相关系图中，西相关系数或偏自相关系数在滞后的前几期处于置信区间之外，而之后的系数基本上都落入置信区间内，且逐渐趋于零的情况。通常把相关系数图在滞后第 \(p\) 期后截尾的情况叫做 \(p\)

拖尾 是指在自相关系图或偏自关系数图中的系数有指数型、正弦型或震荡型衰减的波动，并不会都落在置信区间内。

利用自相关系数图、偏自相关系数图进行模型识别，主要依据如下原则

模型识别	ACF(自相关)	PACF(偏自相关)
AR(p)	衰减趋于零(几何型或振荡型)	P阶后截尾
MA(q)	q阶后截尾	拖尾
ARMA(p,q)	q阶后拖尾	p阶后拖尾

注意：利用图形进行定阶只是一种模型识别的辅助手段。实际上需要根据数据情况反复对模型进行调整。

示例

# 绘制ACF图和PACF图
fig = plt.figure(figsize=(10, 6))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(sales['Sales'].diff(1).iloc[1:92].dropna(), lags=24, ax=ax1)  # ACF图
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(sales['Sales'].diff(1).iloc[1:92].dropna(), lags=24, ax=ax2)  # PACF图
plt.show()
# 由图可知，ACF图在滞后期p=2之后截尾，PACF图随着滞后期扩大，拖尾趋势较明显。
# 可把模型初步定为MA(2)，由于由一阶差分数据而得，故ARIMA(0,1,2)

• 利用最小信息准则进行模型识别

$可以计算 \(ARMA(p,q)\)$

模型选择AIC和BIC：选择更简单的模型

AIC：赤池信息准则(Akaike Information Criterion)

$\[AIC= 2k-2\ln(L) \]$

BIC：贝叶斯信息准则(Bayesian Information Criterion)

$\[BIC = k\ln(n) - 2\ln(L) \]$

其中，\(k\) 为模型参数个数，\(n\) 为样本数量，\(L\)

示例

# 最小信息准则
order_p, order_q, bic = [], [], []
model_order = pd.DataFrame()
for p in range(4):
    for q in range(4):
        arma_model = sm.tsa.ARMA(sales['Sales'].diff(1).iloc[1:92].dropna(), (p, q)).fit()
        order_p.append(p)
        order_q.append(q)
        bic.append(arma_model.bic)
        print('The BIC of ARMA({0},{1}) is {2}'.format(p, q, arma_model.bic))

model_order['p'] = order_p
model_order['q'] = order_q
model_order['BIC'] = bic
P = list(model_order['p'][model_order['BIC'] == model_order['BIC'].min()])
Q = list(model_order['q'][model_order['BIC'] == model_order['BIC'].min()])
print('The best model is ARMA({0}, {1})'.format(P[0], Q[0]))
# The best model is ARMA(0, 2)
# 其中BIC指数最小的是1028.24,其对应模型为ARMA(0, 2),故本例的模型为MA(2)

估计检测

在对时间序列模型进行识别并确定模型的具体形式之后，便可以利用样本数据进行模型参数的估计并对估计结果进行检验。对于ARMA模型，可以对其拟合程度和参数估计显著性等方面进行检验。此外，对于一个适当的ARMA模型，还应当抱枕其残差项无自相关性，即对残差进行白噪声检验。如果模型残差项非白噪声，则需要重新对模型进行调整或识别。

示例

# 模型估计检验
# 1.估计值显著性
model = sm.tsa.ARMA(sales['Sales'].diff(1).iloc[1:92].dropna(), (0, 2)).fit(method='css')
# method有css-mle(默认),mle(极大似然),css(最小二乘)
params = model.params  # 参数估计结果
tvalues = model.tvalues  # t统计量
pvalues = model.pvalues  # p值
result_mat = pd.DataFrame({'Estiamte': params, 't-values': tvalues, 'p-values': pvalues})
print(result_mat)
#              Estiamte  t-values  p-values
# const        1.725820  1.368120  0.171275
# ma.L1.Sales -0.436521 -4.320340  0.000016
# ma.L2.Sales -0.382629 -3.875835  0.000106
# 此处ma.L1.Sales和ma.L2.Sales的Estiamte为模型X_t = a_t-\theta_1a_{t-1}-\theta_2a_{t-2}中-\theta的结果，
# 不是通常理解的\theta的估计结果
# 输出模型的AIC信息指数和方差
print('AIC', model.aic)
print('Variance Estimates:', model.sigma2)
# AIC 1018.4162954572903
# Variance Estimates: 3887.7999374786523
# 针对MA(2)模型X_t = a_t-\theta_1a_{t-1}-\theta_2a_{t-2},其估计值项用const表示。当ARMA
# 模型参数为空时(p,q均为0),const的值为所分析序列的样本均值。模型的两个参数-\theta_1,-\theta_2分别为-4.320340和-3.875835
# 由于ARMA估计过程中的标准误差是基于大数定律的，因此t值和对应的p值在小样本条件下不一定可靠。
# 本例经过一阶差分之后样本量91为大样本，从MA(2)模型两个参数估计值非常显著

# 2.残差项进行检验白噪声检验
resid = model.resid
r, q, p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
mat_res = np.c_[range(1, 41), r[1:], q, p]
table_res = pd.DataFrame(mat_res, columns=['to lag', 'AC', 'Q', 'Prob(>Q)'])
LB_result_res = table_res.loc[[5, 11, 17, 23]]
LB_result_res.set_index('to lag', inplace=True)
print('残差白噪声检验结果\n', LB_result_res)
# 残差白噪声检验结果
#                AC          Q  Prob(>Q)
# to lag
# 6.0    -0.115900   2.048824  0.915154
# 12.0    0.217214   8.508128  0.744269
# 18.0    0.012164  10.126873  0.927686
# 24.0    0.058879  12.924341  0.967351
# 残差项白噪声检验的原假设为：残差是白噪声，备择假设为非白噪声
# 由上面的结果，发现各滞后期的残差项不存在自相关，可认为无法拒绝原假设，即为白噪声
# 残差项的ACF图
fig = plt.figure(figsize=(10, 6))
ax3 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid, lags=24, ax=ax3)
plt.show()
# 观察图可知各滞后期(Lag=0不考虑)的残差项自相关系数，可以看出均处于区间之内在0附近波动，是一白噪声序列
# 正态分布检验
sm.ProbPlot(resid, stats.t, fit=True).ppplot(line='45')
sm.ProbPlot(resid, stats.t, fit=True).qqplot(line='45')
plt.show()
# pp和qq图在45度线上，认为是正态分布
# 正态分布图
plt.figure()
x = pd.Series(resid)
p1 = x.plot(kind='kde')
p2 = x.hist(normed=True)
plt.grid(True)
plt.show()

预测

模型经过识别和参数估计，并进行相应的检验之后，便可利用所建立的模型进行预测。

示例

# 模型的预测
res_p = arma_model.forecast(steps=4)
# 对后续n期数据进行预测
# 参数step指定期数
# 输出3各数组，分别为预测值、预测标准误差、预测值的置信区间(置信区间可通过forecast的参数alpha进行设置，此处默认95%)
print(res_p)
# (array([-38.24694743,  -4.44512376,   1.10646676,   6.72089693]),
# array([58.1070393 , 66.93239305, 71.70980541, 71.71613542]),
# array([[-152.13465171,   75.64075685],
#        [-135.63020354,  126.73995601],
#        [-139.44216918,  141.6551027 ],
#        [-133.8401456 ,  147.28193946]]))
# 注意：上述预测值是对建模时所用的数据进行预测的，本例使用的是一阶分差，得到的预测值也是差分数值，所以第1个数组有负值且偏小
# 观测与预测值对比
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.plot(sales['Sales'].diff(1).iloc[1:92], color='blue', label='Sales')
ax.plot(arma_model.fittedvalues, color='green', label='Predicted Sales')
plt.legend(loc='lower right')
plt.show()


# 预测值与真实值之间比较吻合，且在波动的方向和幅度基本一致

# 克服预测数据与原始数据不一致问题
def forecast(step, var, modelname):
    # 参数step表示预测期数
    diff = list(modelname.predict(len(var) - 1, len(var) - 1 + step, dynamic=True))
    prediction = []
    prediction.append(var[len(var) - 1])
    seq = []
    seq.append(var[len(var) - 1])
    seq.extend(diff)
    for i in range(step):
        v = prediction[i] + seq[i + 1]
        prediction.append(v)
    prediction = pd.DataFrame({'Predicted Sales': prediction})
    return prediction[1:]  # 第一个值是原序列的最后一个值，故从第二个开始是预测值


res_p = forecast(4, sales['Sales'], arma_model)
print(res_p)
#    Predicted Sales
# 1       964.553053
# 2       960.107929
# 3       961.214396
# 4       967.935293