【数学建模】时间序列分析_时间序列自相关

时间序列分析

时间序列分析大致可分成 三大部分：描述过去、分析规律和预测未来。而根据时间和数值性质不同，可分为两类：

1. 时期序列：数值要素反映现象在一定时期内发展的结果
2. 时点序列：数值要素反映现象在一定时点上的瞬时水平

时间序列分解：

1. 长期趋势T：表现在相当长的一段时间内，收到长期趋势影响因素的影响
2. 季节趋势S：季节是广义的，可以按月，季，周为单位，指由于季节转变使指标发生周期性变动
3. 循环变动C：通常以若干年为周期，变现为波浪式的周期变动
4. 不规则变动I：由于某些随机因素导致的数值变化
   

数据：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA

# 1. 生成时间序列数据
np.random.seed(42)
n_samples = 100
time_steps = np.arange(n_samples)
ts_data = np.random.randn(n_samples).cumsum()

plt.figure()
plt.plot(ts_data)
plt.xlabel('时间')
plt.ylabel('值')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

1. 自相关系数（Autocorrelation Coefficient）： 自相关系数是衡量时间序列在不同滞后阶数下的相关性的统计指标。在时间序列中，可能存在当前时刻和之前某些时刻之间的相关性，自相关系数就是用来度量这种相关性的。自相关系数的范围在 -1 到 1 之间，其中 1 表示完全正相关，-1 表示完全负相关，0 表示无相关性。

   自相关系数可以用来检测时间序列数据中的季节性、周期性和其他模式。常用的自相关系数包括一阶自相关系数（ACF(1)）、二阶自相关系数（ACF(2)）等，表示时间序列与自身在不同滞后阶数下的相关性。

2. 滞后阶数（Lag Order）： 滞后阶数是指在计算自相关系数时，选择的时间间隔或滞后步数。它表示当前时刻的观测值与之前第几个时刻的观测值进行比较，用于计算相关系数。选择适当的滞后阶数对于准确描述时间序列的相关性非常重要。

ADF平稳性检验

ADF( Augmented Dickey-Fuller test )检验是一种常用的平稳性检验方法。它基于单位根检验，通过检验时间序列的单位根是否存在来判断序列的平稳性。如果拒绝原假设（存在单位根），则序列不平稳。

# 2. ADF检验
def adf_test(timeseries):
    result = adfuller(timeseries)
    print("ADF 统计量:", result[0])
    print("p-value:", result[1])
    if result[1] <= 0.05:
        print("平稳")
    else:
        print("不平稳")

print("ADF检验结果:")
adf_test(ts_data)
1
2
3
4
5
6
7
8
9
10
11
12

Ljung-Box纯随机性检验

Ljung-Box检验的原假设是序列中的自相关系数为零，即序列呈现纯随机性。检验的统计量计算了一系列滞后阶数下的样本自相关系数的平方和，然后与自由度为滞后阶数的卡方分布进行比较。

# 3. 纯随机数检验
def pure_random_test(timeseries):
    n = len(timeseries)
    acf_vals = np.correlate(timeseries, timeseries, mode='full') / np.sum(timeseries**2)
    lag = int(np.floor(n/2))
    if acf_vals[lag] <= 0.05:
        print("不纯随机")
    else:
        print("纯随机")

print("\n纯随机数检验:")
pure_random_test(ts_data)
1
2
3
4
5
6
7
8
9
10
11
12

指数平滑模型

分为季节性模型和非季节性模型

一元时间序列分析模型

时间序列的平稳性

若时间序列{x_t}满足以下三个条件：

• $$E(x_t)=E(x_{t-s})=u(均值固定为常数)$$

• $$Var(x_i)=Var(x_{t-s})={\sigma }^2(方差存在且为常数)$$

• $$Cov(x_t,x_{t-s})={\gamma }_s(协方差只和间隔s有关，与t无关)$$

则称{x_t}为协方差平稳，又称弱平稳

如果对任意的t和h，多为随机变量$x_{t_1},x_{t_2},...,x_{t_k}$和$x_{t_1+h},x_{t_2+h},...x_{t_k+h}$的联合分布相同，则称{x_t}为严格平稳

若时间序列{x_t}满足以下三个条件：

• $$E(x_t)=E(t_{t-s})=0$$

• $$Var(x_i)=Var(x_{t-s})={\sigma }^2(方差存在且为常数)$$

• $$Cov(x_t,x_{t-s})={\gamma }_s$$

则称{x_t}为白噪声序列

差分方程

将某个时间序列变量表示为该变量的滞后项、时间和其他变量的函数，这个函数方程被称为差分方程

• $$y_t={\alpha }_0+{\alpha }_1{y}_{t-1}+{\alpha }_2{y}_{t_2}+...+{\alpha }_p{y}_{y-p}+{\epsilon }_t(自回归AR(p)模型)$$

• $$y_t={\epsilon }_t+{\beta }_1{\epsilon }_{t-1}+{\beta }_2{\epsilon }_{t-2}+...+{\beta }_q{\epsilon }_{t-q}(移动平均MA(q)模型)$$

• $$y_t={\alpha }_0+\sum _{i=1}^{p=1}{\alpha }_i{y}_{t-i}+{\epsilon }_t+\sum _{i=1}^q{\beta }_i{\epsilon }_{t-i}(自回归移动平均ARMA(p,q)模型)$$

• $$y_t=a+b{y}_{t-1}+c{z}_t+d{z}_{t-1}+{\epsilon }_t$$

• $$y_t=a+b{y}_{t-1}+ct+{\epsilon }_t$$

${\epsilon }_t$不做说明一般指白噪声序列

差分方程的特征方程

差分方程的齐次部分：只包含该变量自身和它的滞后项的式子
$$y_t={\alpha }_0+\sum _{i=1}^p{\alpha }_i{y}_{t-i}+{\epsilon }_t+\sum _{i=1}^q{\beta }_i{\epsilon }_{t-i}(自回归移动平均ARMA(p,q)模型)$$
齐次部分：$y_t={\sum }_{i=1}^p{\alpha }_i{y}_{t-i}$

将其次部分转化为特征方程：令y_t=x^t带入：
$$x^t=\sum _{i=1}^p{\alpha }_i{x}^{t-i}\Rightarrow x^p={\alpha }_1{x}^{p-1}+{\alpha }_2{x}^{p-2}+..+{\alpha }_p$$
特征方程是一个p阶多项式，对应可求出p个解(虚根和实根都有)，p个解的模长大小决定了ARMA(p,q)模型的{y_i}是否平稳

滞后算子

用符号L表示滞后算子lag operator：$L^i{y}_t=y_{t-i}$，满足如下性质：
$$\begin{gathered} (1)LC=C(C为常数) \\ (2)(L^i+L^j)y_t=y_{t-i}+y_{t-j} \\ (3)L^i{L}^j{y}_t=y_{t-i-j} \end{gathered}$$

• 一阶差分
  $$\Delta y_t=y_t-y_{t-1}=(1-L)y_t$$

• 二阶差分
  $${\Delta }^2{y}_t=\Delta (\Delta y_t)=\Delta (y_t-y_{t-1})=(1-L{)}^2{y}_t$$

• n阶差分
  $${\Delta }^n{y}_t=(1-L{)}^n{y}_t$$

p,d,q

1. 自回归阶数 §： 自回归阶数表示模型中考虑的历史观测值的数量，也就是当前值与过去多少个时刻的值有关。例如，ARIMA(p, d, q)中的p就表示自回归阶数。当p为0时，模型不考虑历史观测值，只根据当前时刻的值进行预测；当p为1时，模型根据当前时刻和前一个时刻的值进行预测，以此类推。

自回归模型基于时间序列的过去观测值，将其作为预测当前时刻值的依据。增加p的值可以考虑更多过去时刻的观测值，但过多的自回归项可能会导致模型过拟合。

2. 差分阶数 (d)： 差分阶数表示在进行建模之前，需要对时间序列进行多少次差分以使其平稳。平稳性是ARIMA模型的前提，它指的是序列的均值、方差和自协方差都不随时间变化。如果时间序列不是平稳的，就需要通过差分来减少或消除趋势、季节性等。

当d为0时，说明原始时间序列已经平稳，无需差分；当d为1时，说明进行一次一阶差分；以此类推。通常，需要适当的差分次数才能将时间序列转换为平稳序列，以便构建ARIMA模型。

3. 滑动平均阶数 (q)： 滑动平均阶数表示模型中考虑的过去滑动平均项的数量，也就是当前值与过去多少个时刻的平均值有关。例如，ARIMA(p, d, q)中的q就表示滑动平均阶数。和自回归类似，滑动平均模型也用于考虑过去观测值对当前时刻值的影响。

增加q的值可以考虑更多过去时刻的平均值，但同样要注意过多的滑动平均项可能导致模型复杂性过高。

AR§

AR§是p阶的自回归模型，将自己的1到p阶滞后项视为自变量进行回归，只适用于平稳的时间序列。
$$y_t={\alpha }_0+{\alpha }_1{y}_{t-1}+{\alpha }_2{y}_{t_2}+...+{\alpha }_p{y}_{y-p}+{\epsilon }_t$$
自回归只能适用于预测与自身前期相关的经济现象，即受自身历史因素影响较大的经济现象，如矿的开采量，各种自然资源产量等。对于受社会因素影响较大的经济现象， 不宜采用自回归，而应使用可纳入其他变量的向量自回归模型（多元时间序列）。

将齐次部分转化为特征方程后，有：
$$x^p={\alpha }_1{x}^{p-1}+{\alpha }_2{x}^{p-2}+...+{\alpha }_p$$
这样可以求出p个解：

1. 若p个解模长均小于1，则{y_t}序列平稳，AR§模型平稳
2. 若p个解有k个解模长等于1，则{y_t}为k阶单位根过程，可经过k阶差分变为平稳序列
3. 若任意一个解模长大于1，则{y_t}为爆炸过程

MA(q)

q阶移动平均过程
$$y_t={\epsilon }_t+{\beta }_1{\epsilon }_{t-1}+{\beta }_2{\epsilon }_{t-2}+...+{\beta }_q{\epsilon }_{t-q}$$

可知，移动平均模型具有可逆性：一阶MA可以转换为无穷阶自回归AR。
类似的，可逆性条件下，MA(q)可以转化为无穷阶自回归过程

平稳性：只要q是常数，MA(q)一定平稳。

ARMA(p,q)

自回归移动平均模型(Autoregressive Moving Average)，将MA和AR结合起来

$$y_t={\alpha }_0+\sum _{i=1}^{p=1}{\alpha }_i{y}_{t-i}+{\epsilon }_t+\sum _{i=1}^q{\beta }_i{\epsilon }_{t-i}$$
平稳性只和自回归AR§部分有关

ACF自相关系数/PACF偏自相关系数

ACF表示时间序列与其自身在不同滞后阶数下的相关性。它衡量了时间序列当前时刻与之前各个时刻之间的相关程度。ACF函数的图像可以帮助我们观察时间序列的周期性和季节性，以及确定最合适的滞后阶数。

1. 相关系数：X和Y是两个随机变量，$\rho =\frac{cov(X,Y)}{\sqrt{Var(x)}\sqrt{Var(y)}}$表示X和Y的相关系数(线性关系)

2. 自相关系数(autocorrelation)：假定时间序列平稳，定义自相关系数为：
   $${\rho }_s=\frac{cov(x,x_{t-s})}{\sqrt{Var(x_t)}\sqrt{Var(x_{t-s})}}=\frac{{\gamma }_s}{{\gamma }_0},{\gamma }_0=Cov(x_t,x_t)=Var(x_t)$$
   将$\rho$视为s的函数，则称该函数为自相关函数，函数图像为自相关图
   由于真实的自相关系数未知，只能依赖样本数据进行估计
   对于白噪声序列：
   ρ s = { 1 , s = 0 0 , s ≠ 0 \large\rho_s=

   $$\begin{cases} 1, s=0\\ 0, s\ne0 \end{cases}$$ ρs​=⎩ ⎨ ⎧​1,s=00,s=0​

PACF表示两个时刻之间的相关性，消除了中间时刻的影响，即在考虑其他滞后阶数的影响时，直接比较两个时刻之间的关系。PACF可以帮助我们确定时间序列的阶数（order），即在AR（自回归）或MA（移动平均）模型中需要考虑的滞后阶数。

# 4. 绘制ACF图和PACF图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
plot_acf(ts_data, ax=ax1)
plot_pacf(ts_data, ax=ax2)
plt.show()
1
2
3
4
5

AIC/BIC准则

Akaike information criterion:AIC=2(模型中参数的个数)-2ln(模型的极大似然函数值)

Bayesian information criterion:BIC=ln(T)(模型参数的个数)-2ln(模型的极大似然函数值)

• 样本个数：T
• 模型中参数的个数n：反应模型复杂程度
• 模型的极大似然函数值：反应模型对于数据解释程度

使AIC或BIC最小，BIC对模型复杂程度的惩罚系数更大，因此BIC往往比AIC选择的模型更简洁

模型检验

估计完成时间序列模型后，我们需要对残差进行白噪声检验，如果残差是白噪声，则说明我们选取的模型能完全识别出时间序列数据的规律，即模型可接受；如果残差不是白噪声，则说明还有部分信息没有被模型所识别，我们需要修正模型来识别这一部分的信息

Q检验：
$$\begin{gathered} H_0:{\rho }_1={\rho }_2=...={\rho }_s=0,H_1:{\rho }_i(i=1,2,...s)至少有一个不为0 \\ 若H_0成立，统计量Q=T(T+2)\sum _{k=1}^s\frac{r_k^2}{T-k}\sim {\chi }_{s-n}^2 \end{gathered}$$
p<0/05：拒绝原假设，模型需要修正

• T:样本个数
• n：未知参数个数，例如ARMA(p,q):n=p+q+1
• s一般取8,16,24等，SPSS里取的是18

ARIMA(p,d,q)

$$\begin{gathered} y_t^{\prime }t={\alpha }_0+\sum _{i=1}^{p=1}{\alpha }_i{y}_{t-i}^{\prime }+{\epsilon }_t+\sum _{i=1}^q{\beta }_i{\epsilon }_{t-i} \\ 且y_t^{\prime }={\Delta }^d{y}_t=(1-L{)}^d{y}_t \end{gathered}$$

SARIMA(p,d,q)

# 5. 构建ARIMA模型并进行预测

# p 是自回归阶数，d 是差分阶数（在这里设置为0，即不进行差分），q 是滑动平均阶数。
p, d, q = 1, 0, 1  
model = ARIMA(ts_data, order=(p, d, q))
model_fit = model.fit()
# 预测的步数
forecast_steps = 5
forecast =  model_fit.forecast(steps=forecast_steps) 

# 6. 预测结果可视化
plt.figure()
plt.plot(time_steps, ts_data, label='原始数据')
plt.plot(np.arange(n_samples, n_samples+forecast_steps), forecast, label='预测数据', color='red')
plt.xlabel('时间')
plt.ylabel('值')
plt.legend()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18