时间序列预测

发布日期: 2023-09-02

更新日期: 2023-11-15

文章字数: 4.7k

阅读时长: 21 分

阅读次数:

时间序列

时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。对时间序列进行观察研究，找寻它的发展规律，预测它将来的走势就是时间序列分析。时间序列经过合理函数变换均可认为是三个部分叠加而成：趋势项部分、周期项部分、随机噪声项部分。

时间序列根据所研究的依据不同，可有不同的分类：

按所研究的对象数量分为一元时间序列和多元时间序列；
按时间的连续性分为离散时间序列和连续时间序列；
按序列的统计特性分为平稳时间序列和非平稳时间序列。如果一个时间序列的概率分布与时间 $t$ 无关，则称该序列为严格的（狭义的）平稳时间序列。如果序列的一、二阶矩存在，而且对任意时刻满足：均值为常数；协方差为时间间隔的函数，称该序列为宽平稳时间序列，也叫广义平稳时间序列；
按时间序列的分布规律来分为高斯型时间序列和非高斯型时间序列。

1. 移动平均法（Moving Average）

一次移动平均是一种简单的平滑方法，通过计算过去 $N$ 个时间点的观测值的平均值来估计当前时刻的值。但由于越远时刻的预测，误差越大，因此，一次移动平均仅应用于后一个时刻的预测。

性质： 主要用于平滑数据，减小随机波动。
趋势： 适用于较弱的趋势。
周期性： 不擅长处理明显的周期性。

设观测序列为 $y_1,y_2,\cdots,y_T$ ，取移动平均的项数 $N<T$ 。一次移动平均的计算公式为

$M^{(1)}_t =\frac{1}{N}\left(y_t+y_{t-1}+\cdots+y_{t-N+1}\right)= \frac{1}{N}\sum^{N-1}_{i=0}y_{t-i},\quad t=N,N+1,\cdots$

则有

$M^{(1)}_t = M^{(1)}_{t-1}+\frac{1}{N}(y_t-y_{t-N})$

$t+1$ 时刻的预测值 $\hat{y}_{t+1}=M^{(1)}_t$ ，其预测标准误差为

$S=\sqrt{\frac{\sum^{T}_{t=N+1}(\hat{y}_t-y_t)^2}{T-N}}$

def simple_moving_average(data, window_size=5):
    moving_average = np.convolve(data, np.ones(window_size)/window_size, mode='valid')
    return moving_average

当预测变量的基本趋势发生变化时，一次移动平均法不能迅速适应这种变化。当时间序列的变化为线性趋势时，一次移动平均法的滞后偏差使预测值偏低，不能进行合理的趋势外推。二次移动平均值的计算公式为

$M^{(2)}_t = \frac{1}{N}\left(M^{(1)}_t+M^{(1)}_{t-1}+\cdots+M^{(1)}_{t-N+1}\right)=M^{(2)}_{t-1}+\frac{1}{N}(M^{(1)}_t-M^{(1)}_{t-N})$

当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均方法建立预测模型。当预测目标的基本趋势与某一线性模型相吻合时，常用二次移动平均法。但序列同时存在线性趋势与周期波动时，可用趋势移动平均法建立预测模型

$\hat{y}_{T+m}=a_T+b_Tm,\quad m=1,2,\cdots$

其中 $a_T=2M^{(1)}_T-M^{(2)}_T,b_T=\frac{2}{N}(M^{(1)}_T-M^{(2)}_T)$ 。

2. 指数平滑法（Exponential Smoothing）

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

性质： 能够捕捉数据的趋势和季节性。
趋势： 适用于有一定趋势的数据。
周期性： 对季节性有一定的适应性。

设观测序列为 $y_1,y_2,\cdots,y_T$ ， $\alpha$ 为加权系数， $0<\alpha<1$ 。一次指数平滑的预测公式为

$\hat{y}_{t+1}=S^{(1)}_t =\alpha y_t+(1-\alpha)S^{(1)}_{t-1}=S^{(1)}_{t-1}+\alpha(y_t-S^{(1)}_{t-1})$

可以进一步展开

$S^{(1)}_{t}=\cdots=\alpha\sum^{\infty}_{j=0}(1-\alpha)^jy_{t-j}$

其中每一项的加权系数符合指数递减规律。

在进行指数平滑时，加权系数的选择是很重要的。 $\alpha$ 的大小规定了在新预测值中新数据和原预测值所占的比重。 $\alpha$ 值越大，新数据所占的比重就愈大，原预测值所占的比重就愈小，反之亦然。可以这样理解，将预测公式改写为

$\hat{y}_{t+1}=\hat{y}_{t}+\alpha(y_t-\hat{y}_{t})$

其中 $\alpha$ 的大小体现了修正的幅度， $\alpha$ 值越大，修正幅度越大。如果时间序列的波动不大，比较平稳，则 $\alpha$ 应该取小一点以减少修正幅度，使预测模型能够包含较长的时间序列信息。如果时间序列具有迅速且明显的变动倾向，则 $\alpha$ 应取大一点，使得预测模型的灵敏度高一点，以便迅速更上数据的变化。在实际中，多取几个 $\alpha$ 值测试，预测误差小就用哪个。

此外，还要确定初始值 $S^{(1)}_0$ 。当时间序列的数据较多，比如在 20 个以上时，初始值对以后的预测值影响很少，可选用第一期数据为初始值。如果时间序列的数据较少，在 20 个以下时，初始值对以后的预测值影响很大这时，就必须认真研究如何正确确定初始值。一般以最初几期实际值的平均值作为初始值。

def exponential_smoothing(data, alpha=0.2):

    smoothed_data = [data[0]]  # 初始值为第一个观测值
    #smoothed_data = [data[:2]/len(data[:2])]

    for t in range(1, len(data)):
        smoothed_value = alpha * data[t] + (1 - alpha) * smoothed_data[-1]
        smoothed_data.append(smoothed_value)

    return smoothed_data

一次指数平滑法虽然克服了移动平均法的缺点。但当时间序列的变动出现直线趋势时，用一次指数平滑法进行预测，仍存在明显的滞后偏差。因此，也必须加以修正。再作二次指数平滑，利用滞后偏差的规律建立直线趋势模型。二次指数平滑法的计算公式为

$\begin{cases} & S^{(1)}_t =\alpha y_t+(1-\alpha)S^{(1)}_{t-1}\\ & S^{(2)}_t =\alpha S^{(1)}_t+(1-\alpha)S^{(2)}_{t-1} \end{cases}$

当时间序列从某时刻 $T$ 开始具有直线趋势时

$\hat{y}_{T+m}=a_T+b_Tm,\quad m=1,2,\cdots$

其中 $a_T=2S^{(1)}_T-S^{(2)}_T,b_T=\frac{\alpha}{1-\alpha}(S^{(1)}_T-S^{(2)}_T)$ 。

import numpy as np
import matplotlib.pyplot as plt

def double_exponential_smoothing(data, alpha, beta):

    n = len(data)
    smoothed_data = [data[0]]
    trend = [data[1] - data[0]]

    for t in range(1, n):
        smoothed_value = alpha * data[t] + (1 - alpha) * (smoothed_data[-1])
        trend_value = beta * smoothed_value + (1 - beta) * trend[-1]

        smoothed_data.append(smoothed_value)
        trend.append(trend_value)

    return smoothed_data, trend

# 示例
data = np.array([3, 5, 7, 2, 8, 10, 11, 65, 72, 81, 99, 100, 150])
alpha = 0.3
beta = 0.3
smoothed_data,trend = double_exponential_smoothing(data, alpha, beta)

at = 2 * smoothed_data[-1] - trend[-1]
bt = alpha / (1 - alpha) * (smoothed_data[-1] - trend[-1])
m = np.array([1, 2])
yh = at + bt * m
print("预测值为：", yh)

# 绘制原始数据和平滑后的数据
plt.plot(data, marker='o', label='Original Data')
plt.plot(np.arange(len(data),len(data)+2),yh, marker='s', linestyle='dashed', label='Smoothed Data')
plt.title('Double Exponential Smoothing')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.show()

3. 季节系数法

这里提到的季节，可以是自然季节，也可以是某种产品的销售季节等。比如，空调、季节性服装的生产与销售所产生的数据等。对于季节性时间序列的预测，要从数学上完全拟合其变化曲线是非常困难的。但预测的目的是为了找到时间序列的变化趋势，尽可能地做到精确。季节系数法的具体计算步骤如下：

收集 $m$ 年的每年各季度（ $n$ 个季度）或者各月份的时间序列样本数据
计算所有季度的算数平均值 $\bar{\alpha}=\frac{1}{mn}\sum^{m}_{i=1}\sum^{n}_{j=1}a_{ij}$
计算同季度的算数平均值 $\bar{\alpha}_{·j}=\frac{1}{m}\sum^{m}_{i=1}a_{ij}$
计算季度系数 $b_j=\bar{\alpha}_{·j}/\bar{\alpha}$
预测。当时间序列是按照季度列出的时候，先求出下一年的年加权平均（这里可以用上面的任意一种预测方法，也可以让 $\{w_i\}$ 为自然数列）

$y_{m+1}=\frac{\sum^{m}_{i=1}w_iy_i}{\sum^{m}_{i=1}w_i}$

接着预测下一年的季度预测值 $y_{m+1,j}=b_j\bar{y}_{m+1}=b_j\frac{y_{m+1}}{n}$ 。

def seasonal_forecast(data):
    m, n = data.shape

    # 计算季节系数
    seasonal_factors = data.mean(axis=0) / data.mean()

    # 创建权重数组，假设权重与年份有关
    weights = np.arange(1, m+1)

    # 计算下一年的预测值
    yearly_forecast = weights.dot(data.sum(axis=1)) / weights.sum()

    # 计算预测年份的季度平均值
    quarterly_forecast = yearly_forecast / n

    # 计算季度预测值
    seasonal_forecast = quarterly_forecast * seasonal_factors

    return seasonal_forecast

4. ARIMA模型（Autoregressive Integrated Moving Average）

ARIMA模型和 ARMA模型是时间序列分析中常用的模型，它们包括了 AR (Auto Regressive) 模型和 MA (Moving Average) 模型。

性质： 适用于平稳或近似平稳的时间序列。
趋势： 可以处理不同阶数的趋势。
周期性： 对季节性有一定的适应性。

下面对这些模型进行简要介绍：

AR§ 模型 (Auto Regressive Model of Order p):
- 定义： AR§ 模型是指在时间序列中，每个观测值是其过去 $p$ 个观测值的线性组合，加上一个噪声项。
- 表达式： 一般表示为 $X_t = c + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \ldots + \phi_p X_{t-p} + \varepsilon_t$ ，其中 $\phi_i$ 是参数， $c$ 是常数， $\varepsilon_t$ 是零均值的平稳白噪声误差项。
MA(q) 模型 (Moving Average Model of Order q):
- 定义： MA(q) 模型是指在时间序列中，每个观测值是其过去 $q$ 个白噪声误差项的线性组合。
- 表达式： 一般表示为 $X_t = c + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \ldots + \theta_q \varepsilon_{t-q} + \varepsilon_t$ ，其中 $\theta_i$ 是参数， $c$ 是常数， $\varepsilon_t$ 是零均值的白噪声误差项。
ARMA(p,q) 模型 (Auto Regressive Moving Average Model of Order p,q):
- 定义： ARMA(p,q) 模型是 AR§ 模型和 MA(q) 模型的结合，即同时考虑过去 $p$ 个观测值和过去 $q$ 个白噪声误差项的影响。
- 表达式： 一般表示为 $X_t = c + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \ldots + \phi_p X_{t-p} + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \ldots + \theta_q \varepsilon_{t-q} + \varepsilon_t$ ，其中 $\phi_i$ 和 $\theta_i$ 是参数， $c$ 是常数， $\varepsilon_t$ 是零均值的白噪声误差项。
ARIMA 模型 (Auto Regressive Integrated Moving Average Model):
- 定义： ARIMA 模型是 ARMA 模型的一种扩展，引入了差分的概念，用来处理非平稳时间序列。
- 表达式： ARIMA(p,d,q)，其中 $d$ 表示差分的次数， $p$ 是 AR 模型的阶数， $q$ 是 MA 模型的阶数。ARIMA 模型通过差分将非平稳时间序列转化为平稳时间序列，然后应用 ARMA 模型。

在总体上，AR§ 模型关注时间序列的自回归结构，MA(q) 模型关注白噪声误差项的影响，而 ARMA(p,q) 模型综合了两者。ARIMA 模型在此基础上引入了差分，使其更适用于非平稳时间序列的建模。

平稳时间序列：序列的统计特性不随时间的平移而变化，即均值和协方差不随时间的平移而变化。

一般来说有两种实现方法：手动建模(statsmodels)和自动建模(pmdarima)

4.1 基本流程

平稳性检测：
- 自相关图检验：观察时间序列的自相关图，检查是否存在明显的趋势和周期性。
- 单位根检验（DF检验、ADF检验）：判断时间序列是否具有单位根，即是否是非平稳的。
白噪声检验：
- Ljung-Box检验：用于检验时间序列的残差序列是否是白噪声，即是否存在序列之间的相关性。
平稳性处理方法：
- 对数变换：通过取对数，可以减小数据的振动幅度，使其更加平滑。
- 平滑技术：移动平均（一段时间内的均值均为估计值）、指数平均（通过变权来提高最近值的权重）
- 差分技术：对等周期间隔的数据进行求减，可以用于减小数据的趋势。
- 分解技术：将时间序列分离成不同的成分，例如长期趋势、季节趋势和随机成分，以便更好地理解和建模数据。
识别ARIMA模型的参数：

通过观察ACF（自相关函数）和PACF（偏自相关函数）来选择ARIMA模型的参数。ACF帮助确定移动平均部分的阶数（q），PACF帮助确定自回归部分的阶数（p）。
拟合ARIMA模型、评估、预测

4.2 代码

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings("ignore")  # 选择过滤警告

# Load data
data = pd.read_csv('sunspots.csv')
data['year'].values.astype(int)
print("样本量为{}个".format(len(data)))
print(data.head())

# Plot original data
data['counts'].plot(marker='o')
plt.xlabel("Year", fontsize=15)
plt.ylabel("Counts", fontsize=15)
plt.legend(loc='best')
plt.grid(linestyle='--')
plt.show()

ax1=plt.subplot(121); plot_acf(data['counts'],ax=ax1,title='self-correlation')
ax2=plt.subplot(122); plot_pacf(data['counts'],ax=ax2,title='partial correlation')


# Decompose the time series to check for seasonality
ts_decomposition = seasonal_decompose(data['counts'], period=10)
ts_decomposition.plot()
plt.show()

for i in range(1,6):
    for j in range(1,6):
        md=sm.tsa.ARIMA(data['counts'],order=(i,0,j)).fit()
        print([i,j,md.aic,md.bic])

# Fit ARIMA model
order = (4, 0, 2)  # Example order, you may need to adjust
arima_model = sm.tsa.ARIMA(data['counts'], order=order)
zmd = arima_model.fit()

# Summary of ARIMA model
print(zmd.summary())


# Plot original data, fitted values, and forecasted values
plt.figure(figsize=(10, 6))
plt.plot(data['counts'], label='Original Data')
plt.plot(zmd.predict(), color='red', label='Fitted Values')
forecast = zmd.forecast(steps=10)
plt.plot(forecast, linestyle='--', marker='o', color='black', label='Forecasted Values')
plt.title('Original Data, Fitted Values, and Forecasted Values')
plt.legend()
plt.show()

import statsmodels.api as sm
import matplotlib.pylab as plt
import pmdarima as pm
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf  # 画图定阶
from statsmodels.tsa.stattools import adfuller                 # ADF检验
from statsmodels.stats.diagnostic import acorr_ljungbox        # 白噪声检验

import warnings
warnings.filterwarnings("ignore")  # 选择过滤警告

data=pd.read_csv('sunspots.csv',index_col='year')
data.index = pd.to_datetime(data.index, format='%Y')
print("样本量为{}个".format(len(data)))
print(data.head())


data.plot(marker='o')
# 设置坐标轴标签
plt.xlabel("Year", fontsize=15)
plt.ylabel("Counts", fontsize=15)
# 图例位置
plt.legend(loc='best')
# 生成网格
plt.grid(linestyle='--')
plt.show()

# 确定平稳序列还是非平稳序列

# 分解数据查看季节性   freq为周期
ts_decomposition = seasonal_decompose(data['counts'], period=10)
ts_decomposition.plot()
plt.show()


# 使用训练集的数据来拟合模型
built_arimamodel = pm.auto_arima(data,
                                 start_p=0,   # p最小值
                                 start_q=0,   # q最小值
                                 test='adf',  # ADF检验确认差分阶数d
                                 max_p=10,     # p最大值
                                 max_q=10,     # q最大值
                                 m=1,        # 季节性周期长度，当m=1时则不考虑季节性
                                 d=None,      # 通过函数来计算d
                                 stepwise=True,  # You can set stepwise to True or False based on your preference
                                 suppress_warnings=True,  # Suppress warnings during the fitting process
                                 seasonal=True,  # Enable or disable seasonality
                                 trace=True,  # Enable or disable printing additional information during the fitting process
                                 error_action='ignore',  # How to handle errors during the fitting process
                                 information_criterion='aic',
                                 )

print(built_arimamodel.summary())

built_arimamodel.plot_diagnostics(figsize=(10, 8))
plt.show()


# Plotting original data and fitted values
plt.figure(figsize=(10, 6))
plt.plot(data, label='Original Data')
plt.plot(built_arimamodel.predict_in_sample(), color='red', label='Fitted Values')
plt.title('Original Data vs Fitted Values')
plt.legend()
plt.show()

# Plotting original data, fitted values, and forecasted values for the next 10 periods
forecast, conf_int = built_arimamodel.predict(n_periods=10, return_conf_int=True)

plt.figure(figsize=(10, 6))
plt.plot(data, label='Original Data')
plt.plot(built_arimamodel.predict_in_sample(), color='red', label='Fitted Values')
plt.plot(pd.date_range(start=data.index[-1], periods=10, freq='Y'), forecast, linestyle='--', marker='o', color='black', label='Forecasted Values')
plt.fill_between(pd.date_range(start=data.index[-1], periods=10, freq='Y'),  conf_int[:, 0], conf_int[:, 1], color='green', alpha=0.2, label='Confidence Interval')
plt.title('Original Data, Fitted Values, and Forecasted Values')
plt.legend()
plt.show()