时间序列

时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。分析时间序列的方法构成数据分析的一个重要领域，即时间序列分析。时间序列根据所研究的依据不同，可有不同的分类。

1．按所研究的对象的多少分，有一元时间序列和多元时间序列。

2．按时间的连续性可将时间序列分为离散时间序列和连续时间序列两种。

3．按序列的统计特性分，有平稳时间序列和非平稳时间序列。如果一个时间序列的概率分布与时间t无关，则称该序列为严格的（狭义的）平稳时间序列。如果序列的 一、二阶矩存在，而且对任意时刻t满足：

（1）均值为常数

（2）协方差为时间间隔 $\small \tau$ 的函数。则称该序列为宽平稳时间序列，也叫广义平稳时间序列。我们以后所研究的时间序列主要是宽平稳时间序列。

4．按时间序列的分布规律来分，有高斯型时间序列和非高斯型时间序列。（1）均值为常数

（2）协方差为时间间隔$\tau$的函数。则称该序列为宽平稳时间序列，也叫广义平稳时间序列。我们以后所研究的时间序列主要是宽平稳时间序列。

时间序列的组成部分

时间序列预测技术就是通过对预测目标自身时间序列的处理，来研究其变化趋势的。一个时间序列往往是以下几类变化形式的叠加或耦合。我们常认为一个时间序列可以分解为以下四大部分：

（1）长期趋势变动。它是指时间序列朝着一定的方向持续上升或下降，或停留在某一水平上的倾向，它反映了客观事物的主要变化趋势。

（2）季节变动。

（3）循环变动。通常是指周期为一年以上，由非季节因素引起的涨落起伏波形相似的波动。

（4）不规则变动。通常它分为突然变动和随机变动。

　三种时间序列模型

设

$设:\\ \begin{cases} T_t &长期趋势项\\ S_t &季节变动趋势项\\ C_t &循环变动趋势项\\ R_t &随机干扰项 \end{cases}$

加法模型

$y_t=T_t+S_t+R_t+C_t$

乘法模型

$Y_t=T_t*S_t*C_t*R_t$

混合模型

$\begin{cases} Y_t=T_t*S_t+R_t\\ Y_t=S_t+T_t*C_t*R_t\\ \end{cases} \\ Y_y是观测目标的观测记录,\ E(R_t)=0 \quad E(R_t^2)=\sigma^2$

可以说这三个模型是适用于所有时间序列的

移动平均法

移动平均法是根据时间序列资料逐渐推移，依次计算包含一定项数的时序平均数，以反映长期趋势的方法。当时间序列的数值由于受周期变动和不规则变动的影响，起伏较大，不易显示出发展趋势时，可用移动平均法，消除这些因素的影响，分析、预测序列的长期趋势。移动平均法有简单移动平均法，加权移动平均法，趋势移动平均法等。

简单移动平均法

$设观测序列 y_1..y_T,取移动平均项N<T,一次简单移动平均值计算公式为：\\ M_t^{(1)}=\frac{1}{N}(y_t+...+y_{t-N+1})=\frac{1}{N}(y_{t-1}+...+y_{t-N})+\frac{1}{N}(y_t-y_{t-N})=M_{t-1}^{(1)}+\frac{1}{N}(y_t-y_{t-N})$

当预测目标的基本趋势是在某一水平上下波动时，可用一次简单移动平均方法建立预测模型

$ \hat y{t+1}=M_t^{(1)}=\frac{1}{N}(y_t+…+y{t-N+1}),t\in[N,T] $

其预测标准误差为$S=\sqrt{\frac{\sum_{t=N+1}^{T}(\hat y_t-y_t)^2}{T-N}}$

最近 N 期序列值的平均值作为未来各期的预测结果。一般 N 取值范围：5 ≤ N ≤ 200 。当历史序列的基本趋势变化不大且序列中随机变动成分较多时， N 的取值应较大一些。否则 N 的取值应小一些。在有确定的季节变动周期的资料中，移动平均的项数应取周期长度。选择最佳 N 值的一个有效方法是，比较若干模型的预测误差。预测标准误差最小者为好。

加权移动平均法

在简单移动平均公式中，每期数据在求平均时的作用是等同的。但是，每期数据所包含的信息量不一样，近期数据包含着更多关于未来情况的信息。因此，把各期数据等同看待是不尽合理的，应考虑各期数据的重要性，对近期数据给予较大的权重，这就是加权移动平均法的基本思想。

在加权移动平均法中， w*t 的选择，同样具有一定的经验性。一般的原则是：近期数据的权数大，远期数据的权数小。至于大到什么程度和小到什么程度，则需要按照预测者对序列的了解和分析来确定。

指数平滑法

一次移动平均实际上认为最近 N 期数据对未来值影响相同，都加权$\frac{1}{N}$ ；而 N 期以前的数据对未来值没有影响，加权为 0。但是，二次及更高次移动平均数的权数却不是 $\frac{1}{N}$ ，且次数越高，权数的结构越复杂，但永远保持对称的权数，即两端项权数小，中间项权数大，不符合一般系统的动态性。一般说来历史数据对未来值的影响是随时间间隔的增长而递减的。所以，更切合实际的方法应是对各期观测值依时间顺序进行加权平均作为预测值。指数平滑法可满足这一要求，而且具有简单的递推形式。指数平滑法根据平滑次数的不同，又分为一次指数平滑法、二次指数平滑法和三次指数平滑法等。

一次指数平滑法（SES）

$设时间序列y_1,y_2...y_t,\alpha 为加权系数，\alpha \in[0,1] 一次指数平滑公式为\\ S_t^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}=S_{t-1}^{(1)}+\alpha(y_t-S_{t-1}^{(1)}) \quad (13)$

上式（13）是由移动平均公式改进而来的。由式（1）知，移动平均数的递推公式为$M{t-1}^{(1)}+\frac{1}{N}(y_t-y{t-N})$

以$M{t-1}^{(1)}$作为$y{t-N}$的最佳估计，则有$M{t-1}^{(1)}=\frac{y_t}{N}+(1-\frac{1}{N})M{t-1}^{(1)} 令\alpha=\frac{1}{N},S{t}^{(1)}$代替$M{t}^{(1)}$则可得(13)

$St^{(1)}=\alpha\sum{j=0}^{\infty} (1-\alpha)^j y{j-1}$ 此式表明 $S_t^{(1)}$是全部历史数据的加权平均，加权系数分别为$α,α(1−α),α(1−α)^2$ ；显然有$\sum{j=0}^{\infty} \alpha(1-\alpha)^j=1$ 由于加权系数符合指数规律，又具有平滑数据的功能，故称为指数平滑。以这种平滑值进行预测，就是一次指数平滑法。预测模型为:

$\hat y{t+1}=S{t}^{(1)} \ \hat y_1=S_0^{(1)}$ 即

$\hat y_{t+1}=\alpha y_t+(1-\alpha)\hat y_t$ ！！！最终式子

预测模型$\hat y_{n+1}=\alpha y_n+(1-\alpha)\hat y_n$

加权系数选择

在进行指数平滑时，加权系数的选择是很重要的。α 值应根据时间序列的具体性质在 0～1 之间选择。具体如何选择一般可遵循下列原则：①如果时间序列波动不大，比较平稳，则α 应取小一点，如（0.1～0.5）。以减少修正幅度，使预测模型能包含较长时间序列的信息；②如果时间序列具有迅速且明显的变动倾向，则α 应取大一点，如（0.6～0.8）。使预测模型灵敏度高一些，以便迅速跟上数据的变化。在实用上，类似移动平均法，多取几个α 值进行试算，看哪个预测误差小，就采用哪个。

初始值设定

用一次指数平滑法进行预测，除了选择合适的α 外，还要确定初始值 $S_0^{(1)}$。初始值是由预测者估计或指定的。当时间序列的数据较多，比如在 20 个以上时，初始值对以后的预测值影响很少，可选用第一期数据为初始值。如果时间序列的数据较少，在 20个以下时，初始值对以后的预测值影响很大，这时，就必须认真研究如何正确确定初始值，一般以最初几期实际值的平均值作为初始值。

二次指数平滑法

一次指数平滑法虽然克服了移动平均法的缺点。但当时间序列的变动出现直线趋势时，用一次指数平滑法进行预测，仍存在明显的滞后偏差。因此，也必须加以修正。修正的方法与趋势移动平均法相同，即再作二次指数平滑，利用滞后偏差的规律建立直线趋势模型。这就是二次指数平滑法。二次指数平滑考虑了序列的baseline和趋势，相当于加强版的一阶指数平滑法

$S_{t}^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}\\ S_{t}^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_t^{(2)}\\$

式中$S_{t}^{(1)}$ 为一次指数的平滑值； $S_t^{(2)}$ 为二次指数的平滑值。当时间序列$y_t$ ，从某时期开始具有直线趋势时，类似趋势移动平均法，可用直线趋势模型进行预测。还有一个趋势模型

$\hat y_{t+m}=a_t+b_tm\\ \begin{cases} a_t=2S_t^{(1)}-S_t^{(2)}\\ b_t=\frac{\alpha}{1-\alpha}(S_t^{(1)}-S_t^{(2)}) \end{cases}$

同理可得三次指数平滑法（Holt’s Winters方法）

$\begin{cases} S_{t}^{(1)}=\alpha y_t+(1-\alpha)S_{t-1}^{(1)}\\ S_{t}^{(2)}=\alpha S_t^{(1)}+(1-\alpha)S_t^{(2)}\\ S_{t}^{(3)}=\alpha S_t^{(2)}+(1-\alpha)S_t^{(3)}\\ \end{cases}$

预测模型为

$\hat y_{t+m}=a_t+b_tm+C_tm^2\\ \begin{cases} a_t=3S_t^{(1)}-3S_t^{(2)}+S_t^{(3)}\\ b_t=\frac{\alpha}{2(1-\alpha)^2}((6-5\alpha )S_t^{(1)}-2(5-4\alpha)S_t^{(2)}+(4-3\alpha )S_t^{(3)})\\ C_t=\frac{\alpha}{2(1-\alpha)^2}(S_t^{(1)}-2S_t^{(2)}+S_t^{(3)}) \end{cases}$

但是三阶指数平滑法可以针对有趋势也有季节性的序列。当一个序列在每个固定的时间间隔中都出现某种重复的模式，就称之具有季节性特征，而这样的一个时间间隔称为一个季节(理解：比如说在一个周内，销量呈现出重复的模式)。一个季节的长度k为它所包含的序列点个数。

Holt特性曲线趋势法

$\begin{cases} Forecast\ eq:\hat y_{t+h|t}=\zeta_t+hb_t\\ Level\ eq:\zeta_t=ay_t+(1-a)(\zeta_{t-1}+b_{t-1})\\ Trend\ eq:b_t=\beta(\zeta_t-\zeta_{t-1}+(1-\beta)b_{t-1}) \end{cases}$

当曲线趋势成指数级增加或下降则可选择相乘模型。

Holt-winters 季节性预测模型

Holt-winter线性趋势方法的问题在于，趋势在未来是恒定的，无限增加或减少。对于长期预测视野，这可能会有问题。因此，阻尼趋势法是一种增加阻尼参数的方法，以使趋势在未来收敛到恒定值（它使趋势变平）。参数由φ代替。

$\begin{cases} Seasonal\ eq: S_t=\gamma (y_t-L_t)+(1-\gamma)S_{t-s} \\ Level\ eq:\zeta_t=ay_t+(1-a)(\zeta_{t-1}+b_{t-1})\\ Trend\ eq:b_t=\beta(\zeta_t-\zeta_{t-1}+(1-\beta)b_{t-1})\\ Forecast\ eq:\hat y_{t+h|t}=\zeta_t+hb_t\\ \end{cases}$

指数平滑法功能解析

具体可参考model~~是书上讲的不全吗~~

	trend	Seasonlity
SES	×	×
Hlot	√	×
Hlot-Winter	√	√

平滑方法中不同分量的平滑参数

level	trend	damping	seasonality
$\alpha$	$\beta$	$\varphi$	$\gamma$

例题

通过某市1976~1987年电器销售额来预测1988年电器销售额。

code

#-*- coding:utf-8 -*-
# @Author : Dummerfu
# @Contact : https://github.com/dummerchen 
# @Time : 2020/12/16 21:59

import matplotlib as mpl
import numpy as np
from matplotlib import pyplot as plt
from statsmodels.tsa import holtwinters as ht
import statsmodels as sm
mpl.rcParams['font.sans-serif'] = 'SimHei'
mpl.rcParams['axes.unicode_minus'] = False

def exponential_smoothing(data,alpha):
    s=[]
    s.append((data[0]+data[1]+data[2]+data[3])/4)
    for i in range(1,len(data)):
        s.append((alpha*data[i-1]+(1-alpha)*s[i-1]))
    return s

def second_exponential_smoothing(data,alpha):
    s1=exponential_smoothing(data,alpha)
    s2=exponential_smoothing(s1,alpha)

    a=[]
    b=[]
    for i in range(1,len(s1)+1):
        a.append(2*s1[i]-s2[i])
        b.append((s1[i]-s2[i])*alpha/(1-alpha))
    return a,b

def cubic_index_exponential_smoothing(data,alpha):
    s1=exponential_smoothing(data,alpha)
    s2=exponential_smoothing(s1,alpha)
    s3=exponential_smoothing(s2,alpha)

    a=[]
    b=[]
    c=[]
    pre=[]
    for i in range(0,len(data)):
        a.append(3*s1[i]-3*s2[i]+s3[i])
        b.append(((6-5*alpha)*s1[i]-(10-8*alpha)*s2[i]+(4-3*alpha)*s3[i]*alpha)/(2*(1-alpha)**2))
        c.append(s1[i]-2*s2[i]+s3[i])

    pre.append(a[0])
    for i in range(len(a)):
        pre.append(a[i]+b[i]+c[i])

    plt.plot(t,pre,marker='o',label='cubic_index_exponential_smoothing')
    return a,b,c



if __name__ == "__main__":

    data2 = [253993,275396.2,315229.5,356949.6,400158.2,442431.7,495102.9,570164.8,640993.1,704250.4,767455.4,781807.8,776332.3,794161.7,834177.7,931651.5,1028390,1114914]

    data=[20.4,20.6,25.72,34.61,51.77,55.92,80.65,131.11,148.58,162.67,232.26]
    t_ = []
    for i in range(len(data)):
        t_.append(i)

    f=plt.figure('exponential_smoothing')


    t=t_+[i for i in range(len(data)+1,len(data)+5)]
    print(t)
    ses=ht.SimpleExpSmoothing(data).fit()
    plt.plot(t,list(ses.fittedvalues)+list(ses.forecast(4)),color='y',marker='o',label='SES')

    #有加法模型和指数模型，这里用指数模型
    holt=ht.Holt(data,exponential=True).fit(smoothing_level=0.8,smoothing_trend=0.2,optimized=False)
    plt.plot(t,list(holt.fittedvalues)+list(holt.forecast(4)),color='r',marker='o',label='holt')

    holt_winter=ht.ExponentialSmoothing(data,seasonal_periods=4,trend='add',seasonal='add').fit(use_boxcox=True)
    plt.plot(t,list(holt_winter.fittedvalues)+list(holt_winter.forecast(4)),color='g',marker='o',label='holt_winter')

    plt.plot(t_, data, color='b', label='raw', marker='o')
    plt.legend()
    plt.show()

指数平滑法总结

指数平滑是当今行业中应用最广泛、最成功的预测方法之一。

但是由各预测函数和图像可知：

一次只能预测序列后一个值，二次只能预测直线，三次只能预测二次曲线，局限性极强。

并且$\hat y{t+m}对应的是y{t+m-1} 而不是y_{t+m}$

比如m=1时 $\hat y{1+1}对应的是y{1}而不是y_2$

表示没看到差分指数平滑法的解析，是不常用吗❓

平稳时间序列模型

这里的平稳是指序列的统计特性不随时间平移而变化，即均值和协方差不随时间的平移变化。~~简单点理解为图像看的有一种‘美感’~~

这里对ARMA进行介绍，而ARIMA进行主要说明

在此之前还需自行了解 ~~不了解也无妨~~

ACF:自相关系数
PACF:偏自相关系数
白噪声检验
单位根检验

AR(p)序列

MA(p)序列

自回归移动平均序列（ARMA）

ARMA包括

AR序列自回归序列(Auto Regressive Model)
MA序列移动平均序列(Moving Average Model)

ARMA(p)序列

设{X~t~ ,t = 0,±1,±2…} 是零均值平稳序列，满足下列模型：$Xt-\varphi X{t-1}-…-\varphi X{t-p}=\epsilon_t-\theta_1\epsilon{t-1}-…-\thetaq\epsilon{t-q}$ （70）其中 t ε 是零均值、方差是 σ~ε~^2^ 的平稳白噪声，则称X~t~是阶数为p, q的自回归滑动平均序列，简记为ARMA( p, q) 序列。当 q = 0时，它是 AR( p) 序列；当 p = 0 时，它为MA(q) 序列。

差分自回归移动平均模型（ARIMA）

ARIMA模型(包括ARMA模型、AR模型和MA模型)是预测平稳时间序列的一类通用模型。ARIMA模型由三部分组成:

序列滞后值的加权和（自回归（AR）部分）
序列滞后预测误差的加权和（移动平均值（MA）部分）
时间序列的差分(Integrated (I)部分)

ARIMA模型通常被标注为ARIMA（p，d，q），其中p表示AR部分的顺序，d表示差分的顺序（“ I”部分），q表示MA的顺序。ARIMA模型相对ARMA模型，仅多了差分操作，ARIMA模型几乎是所有时间序列软件都支持的，差分的实现与还原都非常方便。而statsmodel中，对差分的支持不是很好，它不支持高阶和多阶差分，为什么不支持，作者大概的意思是说预测方法中并没有解决高于2阶的差分，不过没关系，我们可以先用pandas将序列差分好，然后在对差分好的序列进行ARIMA拟合，只不过这样后面会多了一步人工还原的工作。

ARIMA步骤

因为ARMA不能差分所以说经常不使用（差分就变成ARIMA了)所以仅列出ARIMA的步骤，ARMA步骤类似就不再赘述了

画出自相关图 ACF
画出偏自相关图 PACF
进行单位根检验：如果检验p值大于0.05则说明时间序列是非平稳的需要进行差分~~下图是一阶差分后的adf 很接近于0了，画图也能看出来非常平稳~~，否则是平稳的。
白噪声检：验主要是检验p值是否大于0.05，大于0.05的时间序列是平稳的白噪声时间序列，p值小于0.05的是平稳的非白噪声的时间序列，是平稳的非白噪声的时间序列才可以进行下一步的ARMA分析(p<0.05才能进行ARMA分析，否则应回到上一步继续差分)
通过枚举进行参数分析求得最佳的p，q（类似灵敏度分析）
将参数带入模型，得到预测矩阵
forecas预测已有数据集，predict预测未知数据

#-*- coding:utf-8 -*-
# @Author : Dummerfu
# @Contact : https://github.com/dummerchen 
# @Time : 2020/12/17 16:58

import matplotlib as mpl
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from statsmodels.tsa.arima_model import ARIMA
from dataprocessing import data_process
from statsmodels.tsa.stattools import  adfuller as adf
from statsmodels.stats.diagnostic import acorr_ljungbox as al
from statsmodels.graphics.tsaplots import plot_pacf,plot_acf
mpl.rcParams['font.sans-serif'] = 'SimHei'
mpl.rcParams['axes.unicode_minus'] = False
import warnings
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',
                        FutureWarning)
warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',
                        FutureWarning)
warnings.filterwarnings("ignore")

def bf_arima(data):
    pq = [(0, 2, 0), (0, 2, 1), (0, 2, 2), (1, 2, 0), (1, 2, 1), (1, 2, 2), (2, 2, 0), (2, 2, 1), (2, 2, 2)]
    best_acc = float('inf')
    for i in pq:
        try:
            model = ARIMA(data, order=i)
            arima = model.fit(disp=0)

            acc = arima.aic
            if acc < best_acc:
                best_acc = acc
                best_model = arima
                best_param = i
        except Exception as e:
            #print(e)
            continue
    print(best_acc, best_param, best_model.summary())
    return best_param

if __name__ == "__main__":
    data,raw_data=data_process()
    #print(len(data['Count']),data.head())
    # arima_param=bf_arima(data['Count'])
    arima=ARIMA(data,order=(2,2,2))

    res=arima.fit(disp=0)
    print(res.summary())
    # predict 得到的是序列内的预测值，最多向后只能预测1个（没有自动还原差分
    # forecast 是向后预测，自动还原差分值
    fig,ax=plt.subplots(figsize=(10,5))
    plt.plot(raw_data,label='raw_data')
    forecast=res.forecast(steps=len(raw_data)-len(data))

    #print(forecast[0],predict)
    #res.plot_predict(start=raw_data.index[4],end=raw_data.index[-1],plot_insample=False,dynamic=False,ax=ax)

    # forecast[0].plot(ax=ax,label='forecast')
    plt.plot(pd.to_datetime(raw_data.index[700:]),forecast[0],label='forecast',color='g')
    plt.axvline(x=data.index[699],ymin=0,ymax=20000, linestyle='--', color='r', label='Start forecast')

    plt.legend(loc='best')
    plt.show()