本 notebook 介绍如何使用 statsmodels 的时间序列模型进行预测。
本章只适用于以下状态空间模型类：
sm.tsa.SARIMAXsm.tsa.UnobservedComponentssm.tsa.VARMAXsm.tsa.DynamicFactor
下文的数字、表格与两幅图均为官方 notebook 随附的历史示例输出，保留其原始日期与数值。本稿没有重新拟合模型，也没有重新运行或测量性能。
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
macrodata = sm.datasets.macrodata.load_pandas().data
macrodata.index = pd.period_range("1959Q1", "2009Q3", freq="Q")基本示例
先用 AR(1) 模型预测通货膨胀。在预测之前，先查看这个时间序列：
endog = macrodata["infl"]
endog.plot(figsize=(15, 5))官方示例原输出：
<Axes: >官方示例原输出：
<Figure size 1500x500 with 1 Axes>构建并估计模型
下一步是明确要用来预测的计量经济模型。这里通过 statsmodels 的 SARIMAX 类构建 AR(1) 模型。
构建模型后，用 fit 估计参数。summary 会生成几个便于阅读的结果表格。
# Construct the model
mod = sm.tsa.SARIMAX(endog, order=(1, 0, 0), trend="c")
# Estimate the parameters
res = mod.fit()
print(res.summary())官方示例原输出：
SARIMAX Results
==============================================================================
Dep. Variable: infl No. Observations: 203
Model: SARIMAX(1, 0, 0) Log Likelihood -472.714
Date: Thu, 27 Aug 2026 AIC 951.427
Time: 06:59:42 BIC 961.367
Sample: 03-31-1959 HQIC 955.449
- 09-30-2009
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
intercept 1.3962 0.254 5.488 0.000 0.898 1.895
ar.L1 0.6441 0.039 16.482 0.000 0.568 0.721
sigma2 6.1519 0.397 15.487 0.000 5.373 6.930
===================================================================================
Ljung-Box (L1) (Q): 8.43 Jarque-Bera (JB): 68.45
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 1.47 Skew: -0.22
Prob(H) (two-sided): 0.12 Kurtosis: 5.81
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
预测
结果对象的 forecast 或 get_forecast 方法可以产生样本外预测。
forecast 只提供点预测。
# The default is to get a one-step-ahead forecast:
print(res.forecast())官方示例原输出：
2009Q4 3.68921
Freq: Q-DEC, dtype: float64
get_forecast 提供更完整的结果，还支持构建置信区间。
# Here we construct a more complete results object.
fcast_res1 = res.get_forecast()
# Most results are collected in the `summary_frame` attribute.
# Here we specify that we want a confidence level of 90%
print(fcast_res1.summary_frame(alpha=0.10))官方示例原输出：
infl mean mean_se mean_ci_lower mean_ci_upper
2009Q4 3.68921 2.480302 -0.390523 7.768943
默认置信水平为 95%，可通过 alpha 参数调整。置信水平等于 (1 - alpha) × 100%。上述例子设置 alpha=0.10，对应 90% 置信水平。
指定预测步数
forecast 和 get_forecast 都接受一个表示预测步数的参数。无论是否使用日期索引，都可以传入整数，指定向前预测多少步。
print(res.forecast(steps=2))官方示例原输出：
2009Q4 3.689210
2010Q1 3.772434
Freq: Q-DEC, Name: predicted_mean, dtype: float64
fcast_res2 = res.get_forecast(steps=2)
# Note: since we did not specify the alpha parameter, the
# confidence level is at the default, 95%
print(fcast_res2.summary_frame())官方示例原输出：
infl mean mean_se mean_ci_lower mean_ci_upper
2009Q4 3.689210 2.480302 -1.172092 8.550512
2010Q1 3.772434 2.950274 -2.009996 9.554865
如果数据包含具有明确频率的 Pandas 索引，也可以直接指定预测终点日期。有关索引的更多说明，见本章末尾“索引”一节。
print(res.forecast("2010Q2"))官方示例原输出：
2009Q4 3.689210
2010Q1 3.772434
2010Q2 3.826039
Freq: Q-DEC, Name: predicted_mean, dtype: float64
fcast_res3 = res.get_forecast("2010Q2")
print(fcast_res3.summary_frame())官方示例原输出：
infl mean mean_se mean_ci_lower mean_ci_upper
2009Q4 3.689210 2.480302 -1.172092 8.550512
2010Q1 3.772434 2.950274 -2.009996 9.554865
2010Q2 3.826039 3.124571 -2.298008 9.950087
绘制数据、预测与置信区间
把实际数据、预测值与置信区间放在一张图上，通常很有帮助。绘制方法有很多，下面是一种示例：
fig, ax = plt.subplots(figsize=(15, 5))
# Plot the data (here we are subsetting it to get a better look at the forecasts)
endog.loc["1999":].plot(ax=ax)
# Construct the forecasts
fcast = res.get_forecast("2011Q4").summary_frame()
fcast["mean"].plot(ax=ax, style="k--")
ax.fill_between(
fcast.index, fcast["mean_ci_lower"], fcast["mean_ci_upper"], color="k", alpha=0.1
);官方示例原输出：
<Figure size 1500x500 with 1 Axes>如何理解预测结果
上面的预测几乎是一条直线，看起来可能并不惊艳。这是因为使用了非常简单的单变量预测模型。不过，简单模型也可能具有很强的竞争力。图中区间反映的是该模型假设下的预测不确定性，不是未来实际值必定落入其中的保证。
样本内预测与样本外预测
结果对象还提供 predict、get_prediction，二者既能返回样本内拟合值，也能做样本外预测。predict 只返回点预测，类似 forecast；get_prediction 返回更多结果，类似 get_forecast。
如果只关心样本外预测，直接使用 forecast、get_forecast 通常更方便。
交叉验证
本节所用的部分函数最早在 statsmodels 0.11.0 中引入。
一种常见做法是递归执行向前 h 步的预测，按以下流程交叉验证预测方法：
- 在训练样本上拟合模型参数。
- 从样本末尾生成向前 h 步的预测。
- 将预测与测试集比较，计算误差。
- 扩展样本，加入下一条观测，然后重复。
经济学中有时把这个过程称为伪样本外预测评估，或时间序列交叉验证。
示例
这里使用前面的通胀数据，做一个简单的递归评估。完整数据集包含 203 条观测；为了便于说明，前 80% 用作训练样本，首先只考虑向前一步的预测。
上述过程的一次迭代如下：
# Step 1: fit model parameters w/ training sample
training_obs = int(len(endog) * 0.8)
training_endog = endog[:training_obs]
training_mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend="c")
training_res = training_mod.fit()
# Print the estimated parameters
print(training_res.params)官方示例原输出：
intercept 1.162076
ar.L1 0.724242
sigma2 5.051600
dtype: float64
# Step 2: produce one-step-ahead forecasts
fcast = training_res.forecast()
# Step 3: compute root mean square forecasting error
true = endog.reindex(fcast.index)
error = true - fcast
# Print out the results
print(
pd.concat(
[true.rename("true"), fcast.rename("forecast"), error.rename("error")], axis=1
)
)官方示例原输出：
true forecast error
1999Q3 3.35 2.55262 0.79738
加入下一条观测时，可以使用结果对象的 append 或 extend。二者可以得到相同的预测，但能取得的其他结果不同：
append更完整：始终保存所有训练观测的结果，并允许基于新增观测重新拟合模型参数。默认不会重新拟合。extend更快，适合训练样本很大的情况：只保存新增观测的结果，也不允许重新拟合，只能继续使用之前估计的参数。
训练样本较小，例如不到几千条，或希望尽可能得到更好的预测时，可以使用 append。如果样本过大使该方法不可行，或能接受参数稍旧、预测可能略逊，则可以考虑 extend。这是一种取舍，不是对任何数据集准确率的保证。
第二次迭代使用 append 并重新拟合参数。注意 append 默认不重估参数，这里通过 refit=True 显式覆盖默认行为：
# Step 1: append a new observation to the sample and refit the parameters
append_res = training_res.append(endog[training_obs : training_obs + 1], refit=True)
# Print the re-estimated parameters
print(append_res.params)官方示例原输出：
intercept 1.171544
ar.L1 0.723152
sigma2 5.024580
dtype: float64
估计出的参数与第一次略有不同。有了新的 append_res 结果对象，就可以从比上次多一条观测的位置开始预测：
# Step 2: produce one-step-ahead forecasts
fcast = append_res.forecast()
# Step 3: compute root mean square forecasting error
true = endog.reindex(fcast.index)
error = true - fcast
# Print out the results
print(
pd.concat(
[true.rename("true"), fcast.rename("forecast"), error.rename("error")], axis=1
)
)官方示例原输出：
true forecast error
1999Q4 2.85 3.594102 -0.744102
把这些步骤组合起来，就得到以下递归预测评估：
# Setup forecasts
nforecasts = 3
forecasts = {}
# Get the number of initial training observations
nobs = len(endog)
n_init_training = int(nobs * 0.8)
# Create model for initial training sample, fit parameters
training_endog = endog.iloc[:n_init_training]
mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend="c")
res = mod.fit()
# Save initial forecast
forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)
# Step through the rest of the sample
for t in range(n_init_training, nobs):
# Update the results by appending the next observation
updated_endog = endog.iloc[t : t + 1]
res = res.append(updated_endog, refit=False)
# Save the new set of forecasts
forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)
# Combine all forecasts into a dataframe
forecasts = pd.concat(forecasts, axis=1)
print(forecasts.iloc[:5, :5])官方示例原输出：
1999Q2 1999Q3 1999Q4 2000Q1 2000Q2
1999Q3 2.552620 NaN NaN NaN NaN
1999Q4 3.010790 3.588286 NaN NaN NaN
2000Q1 3.342616 3.760863 3.226165 NaN NaN
2000Q2 NaN 3.885850 3.498599 3.885225 NaN
2000Q3 NaN NaN 3.695908 3.975918 4.196649
现在得到从 1999Q2 到 2009Q3 每个时间点发出的三个预测。把每个预测从对应时间点的实际 endog 值中减去，就得到预测误差。
# Construct the forecast errors
forecast_errors = forecasts.apply(lambda column: endog - column).reindex(
forecasts.index
)
print(forecast_errors.iloc[:5, :5])官方示例原输出：
1999Q2 1999Q3 1999Q4 2000Q1 2000Q2
1999Q3 0.797380 NaN NaN NaN NaN
1999Q4 -0.160790 -0.738286 NaN NaN NaN
2000Q1 0.417384 -0.000863 0.533835 NaN NaN
2000Q2 NaN 0.304150 0.691401 0.304775 NaN
2000Q3 NaN NaN -0.925908 -1.205918 -1.426649
评估预测时，通常需要均方根误差等汇总指标。先把预测误差按预测步长整理，而不是按日期排列，然后分别计算各步长的均方根误差。
# Reindex the forecasts by horizon rather than by date
def flatten(column):
return column.dropna().reset_index(drop=True)
flattened = forecast_errors.apply(flatten)
flattened.index = (flattened.index + 1).rename("horizon")
print(flattened.iloc[:3, :5])官方示例原输出：
1999Q2 1999Q3 1999Q4 2000Q1 2000Q2
horizon
1 0.797380 -0.738286 0.533835 0.304775 -1.426649
2 -0.160790 -0.000863 0.691401 -1.205918 -0.311464
3 0.417384 0.304150 -0.925908 -0.151602 -2.384952
# Compute the root mean square error
rmse = (flattened**2).mean(axis=1) ** 0.5
print(rmse)官方示例原输出：
horizon
1 3.292700
2 3.421808
3 3.280012
dtype: float64
使用 extend
extend 不会根据新增观测重新估计参数。与会重新拟合的 append(refit=True) 比较时，参数和预测可能不同。不过，下方两段完整主循环分别使用 append(..., refit=False) 和 extend，并不是重新拟合与固定参数之间的对比，应按实际代码和输出理解。
# Setup forecasts
nforecasts = 3
forecasts = {}
# Get the number of initial training observations
nobs = len(endog)
n_init_training = int(nobs * 0.8)
# Create model for initial training sample, fit parameters
training_endog = endog.iloc[:n_init_training]
mod = sm.tsa.SARIMAX(training_endog, order=(1, 0, 0), trend="c")
res = mod.fit()
# Save initial forecast
forecasts[training_endog.index[-1]] = res.forecast(steps=nforecasts)
# Step through the rest of the sample
for t in range(n_init_training, nobs):
# Update the results by appending the next observation
updated_endog = endog.iloc[t : t + 1]
res = res.extend(updated_endog)
# Save the new set of forecasts
forecasts[updated_endog.index[0]] = res.forecast(steps=nforecasts)
# Combine all forecasts into a dataframe
forecasts = pd.concat(forecasts, axis=1)
print(forecasts.iloc[:5, :5])官方示例原输出：
1999Q2 1999Q3 1999Q4 2000Q1 2000Q2
1999Q3 2.552620 NaN NaN NaN NaN
1999Q4 3.010790 3.588286 NaN NaN NaN
2000Q1 3.342616 3.760863 3.226165 NaN NaN
2000Q2 NaN 3.885850 3.498599 3.885225 NaN
2000Q3 NaN NaN 3.695908 3.975918 4.196649
# Construct the forecast errors
forecast_errors = forecasts.apply(lambda column: endog - column).reindex(
forecasts.index
)
print(forecast_errors.iloc[:5, :5])官方示例原输出：
1999Q2 1999Q3 1999Q4 2000Q1 2000Q2
1999Q3 0.797380 NaN NaN NaN NaN
1999Q4 -0.160790 -0.738286 NaN NaN NaN
2000Q1 0.417384 -0.000863 0.533835 NaN NaN
2000Q2 NaN 0.304150 0.691401 0.304775 NaN
2000Q3 NaN NaN -0.925908 -1.205918 -1.426649
# Reindex the forecasts by horizon rather than by date
def flatten(column):
return column.dropna().reset_index(drop=True)
flattened = forecast_errors.apply(flatten)
flattened.index = (flattened.index + 1).rename("horizon")
print(flattened.iloc[:3, :5])官方示例原输出：
1999Q2 1999Q3 1999Q4 2000Q1 2000Q2
horizon
1 0.797380 -0.738286 0.533835 0.304775 -1.426649
2 -0.160790 -0.000863 0.691401 -1.205918 -0.311464
3 0.417384 0.304150 -0.925908 -0.151602 -2.384952
# Compute the root mean square error
rmse = (flattened**2).mean(axis=1) ** 0.5
print(rmse)官方示例原输出：
horizon
1 3.292700
2 3.421808
3 3.280012
dtype: float64
当前 notebook 的两段完整主循环都保持参数不变：前一段使用 append(..., refit=False)，后一段使用 extend。它们随附的三个步长均方根误差完全相同，依次为 3.292700、3.421808、3.280012，因此这些代码与输出不能支持“extend 的预测稍差、每个步长误差更高”的结论。原说明仍在讨论与 append(refit=True) 的比较，与当前两段主循环的条件不一致。前面单次加入观测的示例确实使用了 append(refit=True)，应与这两个固定参数循环区分。
原作者另以 %%timeit 报告历史计时：extend 约 570ms，append(refit=True) 约 1.7s，并称 extend 比 append(refit=False) 也更快。该计时描述中的重新拟合条件不同于当前两段固定参数主循环；这里保留为原作者的历史报告，不把它当作这些现存单元的重新验证结果。本稿没有重新执行。运行时间取决于数据、模型、依赖版本和硬件，不能当成本机结果或普遍性能保证。
索引
本 notebook 一直使用带有频率信息的 Pandas 日期索引。如下所示，数据为季度频率，范围从 1959Q1 到 2009Q3。
print(endog.index)官方示例原输出：
PeriodIndex(['1959Q1', '1959Q2', '1959Q3', '1959Q4', '1960Q1', '1960Q2',
'1960Q3', '1960Q4', '1961Q1', '1961Q2',
...
'2007Q2', '2007Q3', '2007Q4', '2008Q1', '2008Q2', '2008Q3',
'2008Q4', '2009Q1', '2009Q2', '2009Q3'],
dtype='period[Q-DEC]', length=203)
多数情况下，如果数据对应明确频率的日期或时间索引，例如季度、月度等，最好使用带有正确索引的 Pandas Series。下面给出三种例子：
# Annual frequency, using a PeriodIndex
index = pd.period_range(start="2000", periods=4, freq="Y")
endog1 = pd.Series([1, 2, 3, 4], index=index)
print(endog1.index)官方示例原输出：
PeriodIndex(['2000', '2001', '2002', '2003'], dtype='period[Y-DEC]')
# Quarterly frequency, using a DatetimeIndex
index = pd.date_range(start="2000", periods=4, freq="QS")
endog2 = pd.Series([1, 2, 3, 4], index=index)
print(endog2.index)官方示例原输出：
DatetimeIndex(['2000-01-01', '2000-04-01', '2000-07-01', '2000-10-01'], dtype='datetime64[us]', freq='QS-JAN')
# Monthly frequency, using a DatetimeIndex
index = pd.date_range(start="2000", periods=4, freq="ME")
endog3 = pd.Series([1, 2, 3, 4], index=index)
print(endog3.index)官方示例原输出：
DatetimeIndex(['2000-01-31', '2000-02-29', '2000-03-31', '2000-04-30'], dtype='datetime64[us]', freq='ME')
如果数据频率不规则，应把索引替换为可以明确外推的形式。替换索引只明确了预测位置，并没有解决不规则采样间隔本身的建模问题。
original_index = pd.DatetimeIndex(
[
"2000-01-01 10:08am",
"2000-01-01 11:32am",
"2000-01-01 5:32pm",
"2000-01-02 6:15am",
]
)
new_index = np.arange(4)
endog4 = pd.Series([0.2, 0.5, -0.1, 0.1], index=new_index)
print(endog4.index)官方示例原输出：
Index([0, 1, 2, 3], dtype='int64')
然后把数据传给 statsmodels 模型类，就能得到可预期的预测索引。
mod = sm.tsa.SARIMAX(endog4)
res = mod.fit()例如，做向前一步预测：
res.forecast(1)官方示例原输出：
4 0.011866
dtype: float64新预测的索引是 4，因为输入数据使用整数索引，4 是其下一位置。
最终，数据需要有明确频率，或使用整数序列索引；也可以直接使用没有索引的 NumPy 数组。不过，如果能使用带频率的 Pandas Series，就可以采用更多预测终点指定方式，返回的结果索引也更有意义。
来源与许可
原文：Forecasting in statsmodels。正文、代码、输出和图来自该站点提供的 notebook 源，对应 statsmodels 官方仓库的 notebook。
Copyright (C) 2006, Jonathan E. Taylor；Copyright (c) 2006-2008 Scipy Developers；Copyright (c) 2009-2018 statsmodels Developers。源 notebook 属于官方仓库的 BSD 3-Clause 许可范围，许可条件和免责声明完整保留于下方及 assets/license.txt。
本稿将全部说明文字汉化，保留 29 个代码单元、28 个原始输出及两幅官方教学图，并补充了历史输出、预测区间及不规则索引的边界说明。数值、运行时间和图像来自原 notebook；本次只做静态核对，没有重运行。
许可全文（英文原文）
Copyright (C) 2006, Jonathan E. Taylor
All rights reserved.
Copyright (c) 2006-2008 Scipy Developers.
All rights reserved.
Copyright (c) 2009-2018 statsmodels Developers.
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
a. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
b. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
c. Neither the name of statsmodels nor the names of its contributors
may be used to endorse or promote products derived from this software
without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL STATSMODELS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
暂无评论内容