铁路客货运量时间序列预测：基于ARIMA/SARIMA建模实战-平芜编程栈

简介：针对铁路货运量与客运量进行时序建模预测的Python项目，是一份面向毕业设计、课程设计及期末大作业场景的高分参考。代码中附有详细注释，建模流程完整，初学者也能按步骤理解思路，并快速完成本地部署。整个资源包共包含50个文件，核心由17个Python脚本和4个R脚本构成，同时提供运输量、时序建模等Excel/CSV数据集，以及19张结果图表和说明文档，并附有依赖文件与环境说明，压缩包仅2.05MB，轻量易用。目前已有268人学习使用，可作为铁路运输量时序分析、预测算法选型与结果评估的典型范例。读者不仅能获得可运行的源代码与模型，还能借助评估表格及配套笔记快速复现预测流程，适合迁移应用于其他货运、客运数据分析场景，帮助提升项目完整度与答辩说服力。

1. 铁路货运量/客运量预测：为什么不能直接拿原始序列进模型

铁路货运量和客运量不是能被“一拍脑袋”的回归喂饱的问题。客运量受节假日、调图、高铁新线开通这类事件冲击，货运量又跟宏观经济周期和运输结构调整绑定，两者都带有明显的年度重复模式，同时趋势项会随着经济总量缓慢抬升。这种数据形态放进普通的线性回归或者直接扔给LSTM裸跑，经常出现训练集误差很好看、测试集却系统性偏移的情况，原因是模型把非平稳序列里的“水平漂移”当成了规律。时序建模预测在铁路运输场景里之所以是必选项，是因为它先把趋势、季节、残差拆开，再对自相关结构显式建模，预测结果可解释，也能做出可靠的置信区间。这套资源内置了完整的Python和R代码、运输量.xlsx和时序建模.xlsx等多个格式数据源，代码注释覆盖到每一步，毕业设计、课程设计和期末大作业都能直接复用。熟悉统计建模的读者可以直接跳到第3章看定阶逻辑，想先跑通流程的从第2章开始，半小时内能出第一版预测结果。

2. 数据读取、缺失值修补与平稳性检验：开始建模前的必修课

2.1 多格式数据源的统一载入

项目里同时出现了运输量.xlsx、时序建模.xlsx、ads.csv、currency.csv，这是典型的课程设计数据组织方式：Excel用于人眼检查，CSV用于程序快速读取。实际建模前必须先把它们统一成同一个DataFrame结构，否则后面所有的绘图和建模函数都会因为索引或列名不一致报错。

import pandas as pd import numpy as np df = pd.read_excel("运输量.xlsx", parse_dates=[0], index_col=0) # 如果日期读进来是文本，例如“202301”，需要显式指定格式 df.index = pd.to_datetime(df.index, format="%Y%m") # 统一列名，避免原始文件中英文混杂 df.columns = ["freight", "passenger"] # 先去重，再排序，确保时间索引严格递增 df = df[~df.index.duplicated(keep="first")].sort_index() # 检查时间索引是否连续，连续才能做后续差分和季节分解 print(df.index.is_monotonic_increasing) print(df.isnull().sum())

这里parse_dates=[0]表示把第一列解析成时间，index_col=0让时间列直接成为索引。format="%Y%m"只适用于形如202301的月份文本，如果原始数据是2023-01-01，则不需要这个参数。is_monotonic_increasing返回True才能保证后续diff()和shift()的语义正确。

得到连续时间索引后，还要看缺失值分布。铁路运输量数据偶尔会因为统计口径调整出现一两个空值，这是最常见的坑。处理方式按缺失形态选择，不要一律删行。

缺失场景	推荐处理方式	适用条件
单个缺失点	`df['freight'].interpolate(method='linear')`	序列局部平稳，无剧烈跳变
连续缺失超过2个点	用上一年同月值填充	月度数据且季节性强
突然出现异常低值	用前后三周中位数替代	节假日或疫情等异常扰动

连续缺失时“线性插值”会把断点拉成一条斜线，破坏季节形状，这时候用上年同月值更合理。插值完成后建议再看一眼时间序列图，确认没有人为制造的不自然跳变。

2.2 ADF单位根检验：为什么非平稳序列不能直接建模

平稳性是ARIMA类模型的隐含前提。一个平稳序列的均值、方差和自协方差不随时间改变，模型才能用历史自相关去推断未来。铁路货运量通常带有上升趋势和年度周期，直接从图形上看就不是平稳序列，但不能只靠肉眼判断，需要做ADF单位根检验。

from statsmodels.tsa.stattools import adfuller result = adfuller(df["freight"].dropna()) print("ADF统计量:", result[0]) print("p值:", result[1]) print("1%临界值:", result[4]["1%"]) print("5%临界值:", result[4]["5%"])

ADF检验的原假设是“序列存在单位根，即非平稳”。当p值大于0.05时，不能拒绝原假设，说明序列非平稳，需要差分处理。铁运量数据的p值通常远大于0.05，和用眼睛看到的结果一致。对客运量也要分别执行一次检验，因为两者受不同因素驱动，差分阶数可能不一样。

如果ADF统计量已经小于5%临界值，但p值还不够小，可以尝试对数据取对数后再检验。对数变换能把指数型增长压成近似线性趋势，也能缓解异方差。我一般会先看序列标准差是否随均值增大，如果明显相关，就优先取对数。

2.3 季节性分解：把趋势、季节和残差拆开

平稳性检验只回答“能不能建模”，却不能告诉你是哪种成分在制造非平稳。月度铁路运输量数据几乎肯定有12个月的周期，所以下一步做季节分解。

from statsmodels.tsa.seasonal import seasonal_decompose decomp = seasonal_decompose(df["freight"].dropna(), model="additive", period=12) decomp.trend.plot() decomp.seasonal.plot() decomp.resid.plot()

model="additive"表示把序列拆成“趋势+季节+残差”的加和关系，适合季节波动幅度不随水平变化的数据。如果观察折线图发现每年的波动高度随运量总量放大，则改用model="multiplicative"，它拆的是乘积关系。period=12对月度数据基本是定死的，季报数据才用period=4。

分解之后可以看到残差里如果还留有明显的周期脉冲，说明除了常规季节外还有节假日等事件效应。这时候继续用纯ARIMA会欠拟合，可以把节假日哑变量作为外生变量放进SARIMAX，这也是项目里ads.csv这类文件存在的意义。不要急着把残差扔掉，残差是后续模型诊断的原材料。

3. 从ACF/PACF到ARIMA(SARIMA)定阶：用数据而不是感觉选参数

3.1 ACF和PACF怎么看：截尾、拖尾与差分阶数

差分完成后的序列需要确定AR项阶数p和MA项阶数q。自相关函数ACF度量当前值与滞后值的线性相关性，偏自相关函数PACF则是在排除中间滞后项影响后的净相关。两者的截尾位置直接指向p和q。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 1, figsize=(10, 6)) plot_acf(df["freight_diff12"].dropna(), ax=axes[0], lags=24) plot_pacf(df["freight_diff12"].dropna(), ax=axes[1], lags=24, method="ywm") plt.tight_layout() plt.show()

这里的freight_diff12是在普通一阶差分基础上再做12步季节差分，用于消除年度周期。method="ywm"是Yule-Walker的改进版，对月度数据和小样本更稳定，新版statsmodels里如果不指定，有时会提示使用ywm替代旧的yw。看ACF时重点看滞后1到滞后2处是否突然跌进置信带，看PACF时重点看滞后1到滞后3是否截尾：

实际观察到的现象	推荐的阶数	理由
PACF在滞后1处截尾，ACF拖尾	p=1, q=0	AR(1)过程的自相关呈指数衰减
ACF在滞后1处截尾，PACF拖尾	p=0, q=1	MA(1)过程只有一阶记忆
ACF在滞后12处仍有显著尖峰	增加季节项P或Q	年度季节残余未被完全消除
两者都缓慢衰减	差分阶数不足	先增加d或D，不要急着调p/q

仅仅靠看图定阶容易犹豫，特别是样本量不足时置信带很宽。我通常把ACF/PACF图当成第一步筛选，最终阶数交给信息准则去比。

3.2 用AIC/BIC自动定阶：pmdarima与网格搜索

pmdarima把auto_arima封装成了类似R语言forecast::auto.arima的接口，会自动尝试多组(p,d,q)和(P,D,Q,m)，按AIC或BIC选最优。它比手动点图快，且不容易漏掉局部最优组合。

from pmdarima import auto_arima auto = auto_arima( df["freight"].dropna(), seasonal=True, m=12, start_p=0, max_p=5, start_q=0, max_q=5, d=1, D=1, stepwise=True, information_criterion="aic", trace=True ) print(auto.order) print(auto.seasonal_order)

d=1表示允许一阶差分，D=1表示考虑季节差分，m=12对应月度周期。stepwise=True走逐步搜索路径，速度比全网格快很多，但对比较接近的候选模型可能漏掉最优组合。如果样本量不大（少于300条），可以设stepwise=False做全搜索，时间也可接受。trace=True会把每次尝试的AIC打印出来，方便答辩时展示你是如何选阶的。

如果不想引入额外依赖，也可以自己写两层循环遍历p和q，用SARIMAX拟合后取AIC最小值。注意statsmodels的AIC默认基于极大似然估计，和pmdarima计算口径一致，可以互相验证。

3.3 SARIMA建模：参数含义与外生变量扩展

定阶完成后，直接进入SARIMAX建模。SARIMAX是ARIMA在状态空间框架下的实现，支持季节项、趋势项和外生回归变量，比老的ARIMA类灵活得多。

from statsmodels.tsa.statespace.sarimax import SARIMAX model = SARIMAX( df["freight"].dropna(), order=(1, 1, 1), seasonal_order=(1, 1, 1, 12), trend="c" ) res = model.fit(disp=False) print(res.summary())

order=(p,d,q)里的p是自回归阶数，q是移动平均阶数，d是差分次数。seasonal_order=(P,D,Q,m)对应对季节分量的自回归、差分、移动平均以及一个完整的季节周期长度。trend="c"表示允许一个常数项，等价于在模型里包含漂移项。铁道运输量序列在有下降或上升年份时，trend="c"能让长期预测不偏到一个固定水平上。

res.summary()会给出每个系数的z统计量和p值，p值大于0.05的项可以考虑去掉再重新拟合。这里有一个新手常犯错的地方：不要看到p-value大于0.05就把差分阶数d调成2，差分过度反而会让模型丢失水平信息，得出来的预测区间会宽到失去意义。差分阶数的确认要回到ADF检验，而不是看系数显著性。

4. 训练/测试集划分、滚动回测与误差诊断：让预测结果经得起答辩

4.1 时间序列划分：绝对不能用随机打乱的训练集

交叉验证里的train_test_split默认随机抽样，这在时间序列里是致命的：模型会看见未来数据，训练损失极低，测试损失却一泡污。时间序列预测必须用“前一段训练，后一段验证”的切分方式。

train_size = int(len(df) * 0.8) train, test = df.iloc[:train_size].copy(), df.iloc[train_size:].copy() model = SARIMAX( train["freight"], order=(1, 1, 1), seasonal_order=(1, 1, 1, 12) ) res = model.fit(disp=False)

用80%做训练、20%做测试是比较常见的比例。若数据本身较短，比如只有10年120条月度记录，我会把测试集压缩到最后12条，也就是一整年的长度，保证测试集里包含完整的季节周期。如果测试集只有两三个月，测出来的误差会被某个偶然月份主导，不具参考价值。

4.2 多步预测与误差度量：MAE、RMSE、MAPE怎么读

拟合完成后，用get_forecast做多步预测，一次性预测整个测试集长度。这里与滚动预测不同，所有预测都只基于训练期内信息，最能反映真实上线后的行为。

forecast = res.get_forecast(steps=len(test)) mean = forecast.predicted_mean ci = forecast.conf_int() from sklearn.metrics import mean_absolute_error, mean_squared_error mae = mean_absolute_error(test["freight"], mean) rmse = np.sqrt(mean_squared_error(test["freight"], mean)) mape = np.mean(np.abs( (test["freight"] - mean) / np.clip(test["freight"], 1, None) )) * 100 print(f"MAE={mae:.2f}, RMSE={rmse:.2f}, MAPE={mape:.2f}%")

np.clip(test["freight"], 1, None)是为了防止某个月份客运量为0时除零。MAPE对真实值接近0的场景极度敏感，铁运量一般不会为零，但在其他业务里要注意这个坑。RMSE对离群值惩罚更大，如果后期有突发性大客流，RMSE会比MAE敏感得多。三个指标同时看，不要只挑好看的那一个写进报告。

指标	计算逻辑	业务上的意义
MAE	误差绝对值的平均	平均偏离多少万吨/万人
RMSE	误差平方后开根号	大误差被放大，反映最差情况
MAPE	相对误差百分比	跨量级比较，答辩展示最直观
MASE	与朴素预测对比	小于1说明比用上一年同月值强

4.3 残差白噪声检验：Ljung-Box才是模型是否合格的关键

预测误差再小，如果残差里还有明显的自相关，说明模型没有把时间依赖结构提取干净，后续预测区间也会失真。Ljung-Box检验是残差诊断的标准做法。

from statsmodels.stats.diagnostic import acorr_ljungbox resid = res.resid.dropna() lb = acorr_ljungbox(resid, lags=[6, 12, 24], return_df=True) print(lb[["lb_pvalue"]].round(3))

如果各滞后阶数的p值都大于0.05，说明残差序列是白噪声，模型已经充分提取了数据中的自相关信息。滞后6和滞后12分别关注短期和季节性残留，滞后12尤其重要，因为月度数据的年度周期最容易在这里暴露没洗干净的季节性。如果p值小于0.05，常见处理是提高p/q阶数，或者检查是否遗漏了节假日外部变量。

除了Ljung-Box，我还会画一张残差Q-Q图，看尾部是否严重偏离正态假设。ARIMA的置信区间基于正态近似，如果残差厚尾严重，预测区间会偏窄，实际覆盖概率低于宣称的90%。

5. Python与R双实现对照：从课程设计到可交付的预测流程

5.1 用R的forecast包快速复现同一套ARIMA

项目里带了R代码的note.md，这个设计很务实：Python版适合演示数据清洗和自制模型，R的forecast包则在自动化定阶和预测可视化上更顺手。答辩时用两套代码得出相近结果，本身就是对模型稳定性的证明。

library(forecast) library(readxl) df <- read_excel("运输量.xlsx") ts_freight <- ts(df$freight, start = c(2010, 1), frequency = 12) fit <- auto.arima(ts_freight, stepwise = FALSE, approximation = FALSE) fc <- forecast(fit, h = 12) plot(fc) accuracy(fc)

ts()函数创建时间序列对象，start = c(2010, 1)指从2010年1月开始，frequency = 12是月度数据。auto.arima里的stepwise = FALSE要求执行更大范围的搜索，approximation = FALSE则让每个候选模型都用精确似然估计而不是近似值。这两项会显著增加运行时间，但能降低选到次优模型的概率。

accuracy(fc)直接输出训练集和测试集的MAE、RMSE、MAPE和MASE，参数和Python端基本一一对应。R的输出里还自带MASE指标，这是Python端需要自己算的，所以R适合做最终验证展示。

5.2 两种实现的结果一致性怎么验证

Python的pmdarima和R的forecast虽然都基于Hyndman的算法思路，但在AIC计算、参数估计终止条件上仍有细微差别。两组结果不需要完全相等，但预测曲线的形状应当高度重合。

我一般会在两个环境里分别输出未来12个月的点预测值，然后算一下两组结果的平均相对偏差：

import numpy as np py_pred = np.array([...]) # Python预测结果 r_pred = np.array([...]) # R预测结果 rel_diff = np.mean(np.abs(py_pred - r_pred) / np.clip(r_pred, 1, None)) * 100 print(f"两语言预测平均相对偏差: {rel_diff:.2f}%")

如果偏差在3%以内，说明模型结论稳定，可以直接把这份对比放进毕业设计的验证章节。偏差超过5%，先检查数据对齐：R的ts起点是否正确，Python的差分阶数是否与R一致，而不是急着调参。

5.3 模型更新策略：新数据来了该全量拟合还是增量更新

上线之后的模型不会一成不变。铁路运输量每年都有新数据进来，随着线路开通和区域经济变化，之前估计的系数会逐步失灵。常见做法是每月新增一个真实值后，把模型整体重训一遍。如果每天都要更新，才考虑增量更新。

# 每来一条新数据，用append做状态更新 new_value = [test["freight"].iloc[0]] res = res.append(new_value, refit=False)

append会用新的观测值继续卡尔曼滤波状态递推，保持已估计的系数不变，只更新隐含的状态。refit=False省去重新估计参数的时间，适合快速产出短期预测。但如果长期不重估系数，模型会逐渐偏离数据，我的经验是每12个月至少全量重训一次。evalute.xlsx这类评估文件就是用来专门记录每次重训前后误差对比的，直接往里追加即可。

生产环境里更省事的方式是设置一个误差监控：每次真实值出来后，计算单步预测误差的滚动均值，如果连续三个月超过设定阈值，就触发全量重训。阈值建议用第4章里测试集的RMSE乘1.5作为警戒线。

6. 放空预测区间与模型说服力：三个值得做进毕业设计里的验证技巧

6.1 滚动预测比一次多步预测更接近真实发布场景

一次get_forecast(steps=len(test))的输出会随步长增加越来越不确定，因为误差会累积。实际业务里往往是每个月发布未来一个月或一个季度的预测，所以滚动预测更贴近应用。

history = train["freight"].tolist() preds = [] for y_true in test["freight"]: model = SARIMAX(history, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12)) fit = model.fit(disp=False) preds.append(fit.forecast(steps=1).iloc[0]) history.append(y_true) # 把真实值放进历史，模拟真实发布环境

这段代码每走一步都用包含真实值在内的历史重新拟合，预测效果通常比一次性多步预测好，因为它不断用最新信息校正状态。缺点是时间开销大，样本多时可以把重训频率改成每三步或每五步一次。