时间序列-SARIMAX模型代码和结果检测

一、SARIMAX模型代码示例

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.seasonal import seasonal_decompose
from sklearn.metrics import mean_squared_error

# 加载数据
# sales_data = pd.read_csv('your_sales_data.csv')
# 为了演示，我们创建一个简单的示例 DataFrame
sales_data = pd.DataFrame({
    'date': pd.date_range(start='2024-01-01', end='2024-06-30'),
    'sales': [np.random.normal(loc=1000, scale=300, size=1)[0] for _ in range(181)]
})
sales_data.set_index('date', inplace=True)

# 检查数据平稳性
def test_stationarity(timeseries):
    result = adfuller(timeseries)
    print(f'ADF Statistic: {result[0]}')
    print(f'p-value: {result[1]}')

test_stationarity(sales_data['sales'])

# 季节性分解
decomposition = seasonal_decompose(sales_data['sales'], model='additive', period=30)
decomposition.plot()
plt.show()

# 绘制ACF和PACF图
plot_acf(sales_data['sales'], lags=40)
plt.show()
plot_pacf(sales_data['sales'], lags=40)
plt.show()

# 确定 `d` 和 `D`
# 假设 `d` 为1，`D` 为1

# 确定 `p`, `q`, `P`, `Q`
# 基于ACF和PACF图，我们初步设定 `p`=1, `q`=1, `P`=1, `Q`=1

# 网格搜索最优参数
def optimize_sarima(parameters_list, d, D, s, endog):
    results = []
    best_aic = float('inf')
    for param in parameters_list:
        try:
            model = SARIMAX(endog=endog, order=(param[0], d, param[1]), seasonal_order=(param[2], D, param[3], s))
            results.append(model.fit(disp=-1))
        except:
            continue
        aic = results[-1].aic
        # Save best model, AIC and parameters
        if aic < best_aic:
            best_model = results[-1]
            best_aic = aic
            best_param = param
    return (best_param, best_model)

# 定义参数范围
p = d = q = range(0, 2)
P = D = Q = range(0, 2)
s = 30

parameters = [(x1, x2, x3, x4) for x1 in p for x2 in d for x3 in P for x4 in Q]

# 优化SARIMA模型
best_param, best_model = optimize_sarima(parameters, d=1, D=1, s=s, endog=sales_data['sales'])
print("Best SARIMA Parameters:", best_param)
print("Best Model AIC:", best_model.aic)

# 预测
forecast = best_model.get_forecast(steps=30)
forecast_ci = forecast.conf_int()
ax = sales_data['sales'].plot(label='observed', figsize=(14, 7))
forecast.predicted_mean.plot(ax=ax, label='Forecast', alpha=.7)
ax.fill_between(forecast_ci.index,
                forecast_ci.iloc[:, 0],
                forecast_ci.iloc[:, 1], color='k', alpha=.2)
ax.set_xlabel('Date')
ax.set_ylabel('Sales')
plt.legend()
plt.show()

二、评估模型的有效性

1. Ljung-Box 检验结果

   • 如果只有第5阶及之后的滞后期的 p 值小于 0.05，而之前的滞后期的 p 值均大于 0.05，这可能意味着模型已经捕捉到了短期的相关性，但未能完全处理长期的相关性。
   • 这种情况可能表明模型需要进一步改进，比如增加季节性部分的阶数，或者考虑其他模型结构。

2. 残差直方图

   • 如果残差直方图大致符合正态分布，这表明模型的预测误差是随机的，这对于构建置信区间是有利的。
   • 正态分布的残差通常意味着模型的假设得到了满足，这对于许多统计推断方法是重要的。

3. 图形诊断

   • 使用 ACF 图来辅助判断。如果 ACF 图显示除了第5阶及之后的滞后期之外，大多数滞后期的条形都在置信区间内，那么模型可能是合理的。
   • 仔细查看 ACF 图中第5阶及之后的滞后期，看看是否存在明显的模式或趋势。

4. 综合评价

   • 结合 AIC、BIC 以及 MAE、MSE 等其他统计准则来综合评价模型的性能。
   • 如果其他指标显示模型表现良好，那么即使存在一些小问题，模型仍然可能是有用的。

5. 外样本预测

   • 使用模型对外样本数据进行预测，并评估其预测性能。如果模型在外样本上的表现良好，那么这表明模型具有较好的泛化能力。

6. 模型调整

   • 考虑增加模型的复杂度，例如增加 p, q, P, 或 Q 的值，或者尝试其他模型，如考虑更高阶的模型或者使用不同的季节性周期。
   • 如果可能的话，可以尝试使用自动模型选择工具（如 auto_arima 函数）来寻找更合适的模型。

三、模型拟合结果汇总报告

`statsmodels` 库中的 `Summary` 对象包含了关于拟合模型的大量信息。可以帮助你理解模型的性能和统计显著性。以下是如何查看和解释 `Summary` 的详细信息：

### 查看 Summary

import statsmodels.api as sm

# 假设你已经拟合了一个 SARIMAX 模型
model = sm.tsa.SARIMAX(sales_data['sales'], order=(1, 1, 1), seasonal_order=(1, 1, 1, 30))
results = model.fit()

# 显示 Summary
print(results.summary())

`Summary` 通常包含以下部分：

#### 1. **Model Information**
   - **Dep. Variable**: 依赖变量（在这里是 `sales`）
   - **Method**: 估计方法（通常是最大似然法 `MLE`）
   - **Date**: 模型拟合日期
   - **Time**: 模型拟合时间
   - **No. Observations**: 观测数量
   - **Log-Likelihood**: 对数似然值
   - **AIC**: Akaike 信息准则
   - **BIC**: Bayesian 信息准则
   - **Omnibus**: 统计检验值，用于检测残差是否符合正态分布
   - **Prob(Omnibus)**: Omnibus 检验的概率值
   - **Skew**: 残差的偏度
   - **Kurtosis**: 残差的峰度
   - **Durbin-Watson**: Durbin-Watson 统计量，用于检测残差的自相关性
   - **Cond. No.**: 条件数，表示模型矩阵的条件数

#### 2. **Coefficients Table**
   - **coef**: 参数估计值
   - **std err**: 参数的标准误差
   - **z**: z-统计量，用于检验参数是否显著
   - **P>|z|**: 参数的 p 值，用于检验参数是否显著
   - **[0.025]**: 参数置信区间的下限
   - **[0.975]**: 参数置信区间的上限

#### 3. **Covariance Matrix**
   - 参数估计值之间的协方差矩阵

#### 4. **Residuals**
   - 残差的相关统计信息

### 示例

### 解释

- **Dep. Variable**: 依赖变量为 `sales`。
- **Model**: 模型类型为 `SARIMAX(1, 1, 1)x(1, 1, 1, 30)`，表示非季节性部分的 ARMA(1,1) 和季节性部分的 SARMA(1,1)。
- **Log-Likelihood**: 对数似然值为 123.456。
- **AIC**: Akaike 信息准则为 -234.911。
- **BIC**: Bayesian 信息准则为 -201.890。
- **HQIC**: Hannan-Quinn 信息准则为 -223.222。
- **Coefficients Table**:
  - `ar.L1`: 非季节性自回归系数为 0.5000，标准误差为 0.050，z-统计量为 10.000，p 值为 0.000，置信区间为 (0.403, 0.597)。
  - `ma.L1`: 非季节性移动平均系数为 -0.4000，标准误差为 0.040，z-统计量为 -10.000，p 值为 0.000，置信区间为 (-0.480, -0.320)。
  - `ar.S.L30`: 季节性自回归系数为 0.3000，标准误差为 0.030，z-统计量为 10.000，p 值为 0.000，置信区间为 (0.242, 0.358)。
  - `ma.S.L30`: 季节性移动平均系数为 -0.2000，标准误差为 0.020，z-统计量为 -10.000，p 值为 0.000，置信区间为 (-0.240, -0.160)。
  - `sigma2`: 方差项为 1.0000，标准误差为 0.100，z-统计量为 10.000，p 值为 0.000，置信区间为 (0.805, 1.195)。

通过这些信息，可以评估模型的拟合效果、参数显著性以及其他统计特性。