时间序列分析与预测：从基础到实战-CSDN博客

本文链接：https://blog.csdn.net/m0_62153576/article/details/147185096

在当今的数据驱动时代，时间序列分析和预测发挥着至关重要的作用。它有着广泛的应用场景，比如分析零售连锁店的销售情况，能帮助商家了解销售趋势，合理安排库存；检测服务器流量中的异常，及时发现网络攻击或系统故障；预测股票市场，为投资者提供决策参考。接下来，我们将深入探讨时间序列分析的相关内容。

准备数据集

我们选用德国耶拿的天气数据进行实验。可以使用命令 wget https://storage.googleapis.com/tensorflow/tf-keras-datasets/jena_climate_2009_2016.csv.zip 从命令行界面下载数据，解压后会得到 CSV 格式的数据，用 Pandas 读取即可。本教程主要使用其中的温度数据（摄氏度），数据记录规律，每 10 分钟记录一次，跨度超过 24 小时。以下是读取数据并绘制温度序列图的代码：

import pandas as pd 
import matplotlib.pyplot as plt 
df = pd.read_csv('jena_climate_2009_2016.csv') 
time = pd.to_datetime(df.pop('Date Time'), format='%d.%m.%Y %H:%M:%S') 
series = df['T (degC)'] 
series.index = time 
print(df)
series.plot() 
plt.show()

时间序列的属性

时间序列是按时间顺序依次获取的观测值序列，某个变量随时间变化，我们分析其中的模式，并根据序列过去的变化进行预测。时间序列有几个重要属性：

趋势：变量随时间的整体上升或下降趋势。例如，随着全球气候变暖，多年来的平均气温可能呈现上升趋势。
季节性：时间序列中的周期性成分，某种模式每隔几个时间单位重复出现。就像一年中四季更替，温度会有周期性的变化；一天中，夜晚温度通常低于白天。
残差：时间序列中的噪声成分。

将时间序列分解为这些成分，可以让我们深入了解其变化规律。

滚动统计

可视化滚动统计是理解时间序列的有效方法。我们以 2600（约一个月，即 30 天）为窗口绘制滚动统计，并与数据集的均值和所有数据的最佳拟合线进行比较。以下是绘制滚动均值和滚动标准差的代码：

import numpy as np 
plt.plot(series.index, np.array([series.mean()] * len(series))) 
x = np.arange(len(series)) 
y = series.values 
m, c = np.polyfit(x, y, 1) 
plt.plot(series.index, m*x + c) 
series.rolling(3600).mean().plot() 
plt.legend(['mean', 'regression line', 'rolling mean']) 
plt.ylabel('Temp') 
plt.show()

roll_std = series.rolling(3600).std() 
roll_std.dropna(inplace=True) 
plt.plot(roll_std.index, np.array([roll_std.mean()] * len(roll_std))) 
x = np.arange(len(roll_std)) 
y = roll_std.values 
m, c = np.polyfit(x, y, 1) 
plt.plot(roll_std.index, m*x + c) 
roll_std.plot() 
plt.legend(['rolling std mean', 'rolling std regression', 'rolling std']) 
plt.ylabel('Temp') 
plt.show()

从绘制的图表可以明显看出，温度值有上升趋势，且数据具有很强的季节性。滚动标准差有微弱的下降趋势，但滚动标准差的值本身变化很大。

自相关

相关分析可以通过假设高斯分布，使用皮尔逊系数来比较两个变量之间的关系。自相关则是信号与其延迟副本之间的相关性，是延迟的函数。我们尝试找出时间序列与其时间滞后版本之间的相关性，并观察这些值随时间滞后增加的变化情况。以下是计算自相关的代码：

def correlation(x, y): 
    x_norm = x - x.mean() 
    y_norm = y - y.mean() 
    return np.sum(x_norm * y_norm) / np.sqrt(np.sum(x_norm ** 2) * np.sum(y_norm ** 2))

def autocorrelation(x, k): 
    val_0 = np.sum((x - x.mean()) ** 2) / len(x) 
    val_k = np.sum((x[:-k] - x.mean()) * (x[k:] - x.mean())) / len(x) 
    return val_k / val_0

偏自相关

偏自相关函数可以计算时间序列与其具有时间滞后的同一序列之间的偏相关性。与普通相关性不同的是，偏相关性可以控制其他滞后值的影响。我们可以通过建立线性方程来计算不同滞后值之间的依赖关系。在 Python 中，可以使用 statsmodels 包中的函数来计算偏自相关并绘制相关图：

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf 
def get_acf_pacf_plots(df): 
    fig, ax = plt.subplots(2, figsize=(12,6)) 
    ax[0] = plot_acf(df, ax=ax[0]) 
    ax[1] = plot_pacf(df, ax=ax[1]) 

get_acf_pacf_plots(series[5::6]) 
plt.show()

季节性

从多个图表中可以看出，数据具有季节性元素。直觉上，温度数据也应该有季节性成分，每天温度会有波动，夜晚低于白天，一年中也会有周期性变化。为了验证这一点，可以查看序列的傅里叶变换。傅里叶变换可以将基于振幅的序列转换为基于频率的序列，它是复值函数，将每个序列表示为复平面中正弦波的叠加。以下是绘制傅里叶变换幅度的代码：

from scipy.fft import fft 
fft_vals = fft(series[5::6].values) 
f_per_dataset = np.arange(0, len(fft_vals)) 
n_samples_h = len(series[5::6]) 
hours_per_year = 24*365.2524 
years_per_dataset = n_samples_h/(hours_per_year) 
f_per_year = f_per_dataset/years_per_dataset 
plt.step(f_per_year, np.abs(fft_vals)) 
plt.xscale('log') 
plt.xticks([1, 365.2524], labels=['1/Year', '1/day']) 
_ = plt.xlabel('Frequency (log scale)') 
plt.show()

从绘制的图表中可以看到，在 1/年和 1/天附近的值出现异常峰值，验证了我们之前的直觉。这种方法是一种基本的频谱分析，还有许多复杂的频谱密度估计方法，如周期图、巴特利特方法、韦尔奇方法等，并且常常会使用窗函数（如汉明窗和布莱克曼 - 哈里斯窗）来平滑时间序列。

时间序列的平稳性

当时间序列的统计特性（如均值、方差和自相关性）不随时间变化时，该时间序列是平稳的。像 ARIMA 及其变体等方法在建模时假设时间序列是平稳的，如果时间序列不平稳，这些方法的效果会大打折扣。不过，我们可以通过一些测试来确定时间序列是否平稳，还可以使用一些变换将非平稳时间序列转换为平稳时间序列。

自相关函数（ACF）和偏自相关函数（PACF）图

观察时间序列的均值和方差随时间的变化，可以初步判断其非平稳性，并推测可以应用的平稳化机制。从之前的图表可以猜测，天气时间序列远非平稳。还可以查看 ACF 和 PACF 图，如果图上的值迅速下降，那么时间序列可能接近平稳。以下是按天和按 31 天取数据绘制 ACF 和 PACF 图的代码：

series2 = series[119::120] 
get_acf_pacf_plots(series2) 
plt.show()

series2 = series2[30::31] 
get_acf_pacf_plots(series2) 
plt.show()

从图中可以看到，有几个点落在蓝色区域之外，这表明时间序列不平稳。

迪基 - 富勒检验（Dickey - Fuller Test）

单位根检验用于检验自回归模型中是否存在单位根。原检验将时间序列视为滞后 1 的自回归模型，单位根的存在证明时间序列不平稳。主要有三种版本的检验：

单位根检验： $\Delta y_{t} = \delta_{y_{t - 1}} + u_{t}$
带漂移的单位根检验： $\Delta y_{t} = a_{0} + \delta_{y_{t - 1}} + u_{t}$
带漂移和确定性时间趋势的单位根检验： $\Delta y_{t} = a_{0} + a_{1}t + \delta_{y_{t - 1}} + u_{t}$

增强迪基 - 富勒检验是一种单尾检验，应用于滞后 $p$ 的序列。可以使用信息准则（如赤池、贝叶斯或汉南 - 奎因）来确定滞后值。以下是进行迪基 - 富勒检验的代码：

from statsmodels.tsa.stattools import adfuller 
def test_dickey_fuller_stationarity(df): 
    dftest = adfuller(df, autolag='AIC') 
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p - value', 'Number of Lags Used', 'Number of Observations Used']) 
    for key, value in dftest[4].items(): 
        dfoutput['Critical Value (%s)' % key] = value 
    if dfoutput['Critical Value (1%)'] < dfoutput['Test Statistic']: 
        print('Series is not stationary with 99% confidence. ') 
    elif dfoutput['Critical Value (5%)'] < dfoutput['Test Statistic']: 
        print('Series is not stationary with 95% confidence. ') 
    elif dfoutput['Critical Value (10%)'] < dfoutput['Test Statistic']: 
        print('Series is not stationary with 90% confidence. ') 
    else: 
        print('Series is possibly stationary. ') 
    return dfoutput 

out = test_dickey_fuller_stationarity(series) 
print(out)

KPSS 检验

KPSS 检验也是一种单位根检验，与迪基 - 富勒检验不同的是，它将存在单位根作为备择假设。以下是进行 KPSS 检验的代码：

from statsmodels.tsa.stattools import kpss 
def test_kpss(df): 
    dftest = kpss(df) 
    dfoutput = pd.Series(dftest[0:3], index=['Test Statistic', 'p - value', 'Number of Lags Used']) 
    for key, value in dftest[3].items(): 
        dfoutput['Critical Value (%s)' % key] = value 
    if abs(dfoutput['Critical Value (1%)']) < abs(dfoutput['Test Statistic']): 
        print('Series is not stationary with 99% confidence. ') 
    elif abs(dfoutput['Critical Value (5%)']) < abs(dfoutput['Test Statistic']): 
        print('Series is not stationary with 95% confidence. ') 
    elif abs(dfoutput['Critical Value (10%)']) < abs(dfoutput['Test Statistic']): 
        print('Series is not stationary with 90% confidence. ') 
    else: 
        print('Series is possibly stationary. ') 
    return dfoutput 

out = test_kpss(series) 
print(out)

我们的时间序列在不同检验中得到了相互矛盾的结果。但如果序列不平稳，将其平稳化很重要，因为自回归模型的设计基于平稳性假设。

平稳化变换

有几种方法可以使时间序列平稳：

去趋势

可以通过从序列中减去滚动均值并进行归一化来消除时间序列中的趋势。以下是消除趋势的代码：

def eliminate_trends(series): 
    roll = series.rolling(4).mean() 
    avg_diff = (series - roll) / roll 
    avg_diff.dropna(inplace=True) 
    return avg_diff 

diff = eliminate_trends(series[5::6]) 
test_dickey_fuller(diff) 
test_kpss(diff)

如果序列中有线性趋势，也可以找到回归直线并相应地去除线性趋势：

def eliminate_linear_trend(series): 
    x = np.arange(len(series)) 
    m, c = np.polyfit(x, series, 1) 
    return (series - c) / m 

diff = eliminate_linear_trend(series[5::6]) 
test_dickey_fuller(diff) 
test_kpss(diff)

差分

可以选择一个滞后值 $p$ ，将时间 $t$ 的值与时间 $t - p$ 的值相减。以下是进行差分的代码：

def difference(series, lag = 1): 
    differenced = [] 
    for x in range(lag, len(series)): 
        differenced.append(series[x] - series[x - lag]) 
    return pd.Series(differenced) 

diff = difference(series[6::5]) 
test_dickey_fuller(diff) 
test_kpss(diff)