用时间序列分析提升回归模型的稳健性 — 第一部分

利用时间序列分析提升回归模型的稳健性 — 第一篇

新加坡公屋转售价格案例研究。在这个故事中，我们展示了如何将时间序列分析与线性回归相结合，提高模型的预测能力。

介绍

距离家1.5小时的地方，新加坡总是让我着迷。被更大的邻国包围，这个小小的国家打破了规律。从独立初期的卑微开始，它现在成为了一个重要的金融中心和国际港口。尽管取得了成功，新加坡面临着一个重大挑战——土地稀缺。这个问题导致房价异常高涨。

为了应对这个问题，新加坡政府于1960年建立了新加坡房屋发展局（HDB），致力于为居民提供适当的住房。通过该机构，政府已成功为80%以上的人口提供了质量好、卫生条件良好、价格合理的公共住房，详见HDB官方网站。

认识到HDB的重要性，本文将深入研究基于公开可用数据集的新加坡HDB转售价格，利用数据驱动的方法进行分析。由于这个数据集中包含大量与转售价格和相关变量的信息，因此具有构建回归模型的潜力，令人感兴趣。

令人惊讶的是，在VoAGI和Kaggle中，我找到的出色分析中，很少有人深入讨论这个数据集的时间序列特性。我发现这个分析空白挑战着我。

• 揭示卡通动画的未来：AnimeInbet在线动画描线中的创新
• 生成型人工智能是否值得其环境足迹？’ (Shēngchéngxíng rén​gōng zhì​néng shìfǒu zhídé qí huánjìng zújì?)
• “ConvNet是否卷土重来？揭示其在大规模Web数据集和匹配视觉转换器上的性能”

因此，本文的目标是无缝地将时间序列分析和预测纳入回归模型中。希望本文的观点还能为现有模型增加新的见解，提高鲁棒性。

工具

我在分析过程中使用的编程语言是Python，还使用了以下库（但不限于）：

import pandas as pd  # 结构化数据操作import numpy as np # 科学与数值计算import scipy # 科学与数值计算import matplotlib.pyplot as plt # 数据可视化import seaborn as sns # 数据可视化import statsmodels.api as sm # 统计分析import datetime as dt # 操作日期和时间数据

数据提取、清洗和预处理

（关于这个分析的从数据提取到构建回归模型的完整代码，请随时查看我的GitHub存储库）

可以从这个网站上下载公开可用且免费使用的HDB转售价格数据集的CSV文件。一旦我们获得了CSV文件，我们可以将数据加载到Pandas DataFrame中，得到以下表格：

# 将数据集加载到DataFrame中 # 我下载数据集的位置在这里：https://data.gov.sg/dataset/resale-flat-prices# 我选择了2017年1月至2023年8月的时间范围data = pd.read_csv('Housing_Resale_Dataset.csv')data.head()

原始数据集包含160,795行（交易数量）和11列（变量数量）。由于我们的重点是预测HDB转售价格，列’resale_price’成为了因变量，或目标变量，而其他列成为独立变量，或特征。虽然数据集已经相对整洁，但我们进行了一些数据清洗和预处理，包括列重命名、提取新特征、去重和重新排列，得到了最终的干净数据集，如下所示：

cleandata.head()

从上面的图片中，我们可以看到我们可以从“resale_date”列（之前名为“month”且在图片上不可见）提取信息来创建两个新特征，“resale_year”和“resale_month”，这样我们在清理后的数据集中有13列。与此同时，行数略有减少，变为160,454行，这是由于删除了重复项的结果。

探索性数据分析

接下来，我们将创建可视化图表来揭示我们数据中最重要的一些信息。虽然我主要在Jupyter Notebook中使用Matplotlib和Seaborn创建图表进行直接分析和灵活性，但在博客的这个部分中，我还使用由Tableau生成的图像以提供精美且易读的演示。

价格分布

首先，让我们来看一下下面的新加坡组屋HDB转售价格的直方图:

cleandata['resale_price'].describe()

使用40,000作为分组间隔，我们可以观察到“resale_price”分布呈右偏态，平均值为新币486,519，中位数为新币455,000。

各镇的平均价格

使用条形图，我们可以发现不同镇的组屋平均价格。紫色柱子表示平均价格最高的五个镇。从这张图中，您还会注意到不同镇之间在平均价格上存在一些显著差异。

月平均价格

由于我们在此博客中的主要目标是展示时间序列分析如何增强回归模型，因此展示数据可视化的时间序列特征是至关重要的。因此，下面是2017年1月至2023年8月的组屋月平均价格。请注意，蓝色虚线表示实际的月平均价格，而紫色线表示6个月移动平均线。

如果我们仔细观察蓝色曲线，我们可以看到它呈现出锯齿状的模式，这可能表明我们的数据集呈现出月度季节性。另一方面，紫色线显示了数据的趋势。从紫色线上，我们可以看出我们的数据集具有非平稳的正趋势，即月平均价格在2018年中期下降并在2019年初趋于稳定，但从2020年中期开始，平均价格开始上涨。

每月交易量

下面的图表显示了2017年1月至2023年8月的每月交易量。该图还显示了交易数据呈现季节性，在每年的12月和1月左右，每月交易量通常会下降。我们还可以看到在2020年第一季度，交易量骤降，这很可能是由于COVID-19封锁导致的。

时间序列分析

如引言所述，本文的重点是了解HDB公寓的定价。因此，作为时间序列分析的先决条件，首先，我们将通过创建Pandas系列以月均价格为指标，构建我们的时间序列数据集。

price_time = cleandata.groupby('resale_date')['resale_price'].mean()price_time.head()

分解

通过对月均价格可视化的观察，我们可以观察到数据的趋势和可能的季节性。因此，接下来，我们想要做的是将数据分解为其组成部分，即趋势、季节性和残差，使用时间序列分解。

import statsmodels.api as sm# 我们使用'加法'模型，因为季节性似乎不会随着时间的推移而放大# 我们添加了6个外推期间，使趋势延伸到前6个月和后6个月price_decomp = sm.tsa.seasonal_decompose(price_time, model = 'additive', extrapolate_trend=6) price_decomp.plot();

通过简单移动平均法，我们已经看到了时间序列数据的趋势。同时，通过进行分解，我们还可以得知数据在12个月的周期内存在季节性。从上图可见，季节性相对于平均价格来说较小，范围在-4,000至2,500之间。为了更容易观察，我们将季节性图表分离出来，并将线型图改为条形图：

price_season = price_decomp.seasonal # 季节性 price_season_month = price_season.reset_index()price_season_month['resale_month'] = price_season_month['resale_date'].dt.strftime('%b') # 提取月份信息months_order = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']price_season_month['resale_month'] = pd.Categorical(price_season_month['resale_month'], categories= months_order, ordered = True) # 创建月份数据的顺序pm = price_season_month.groupby('resale_month')['seasonal'].mean()pm.plot(kind = 'bar')

数据分割

在建立可评估可靠性的模型之前，我们首先采用一种称为数据分割的方法。这个过程涉及将数据分成两部分 — 训练数据集和测试数据集。我们使用训练数据集构建时间序列模型，然后使用测试数据集评估模型的性能。因此，在时间序列建模的背景下，训练数据集代表“过去”，而测试数据集代表“未来”。

# 我们使用对数转换是因为这种方式能够使模型表现更好
transform = np.log(price_time) # 对价格时间进行对数转换

# 训练集-测试集拆分
price_test = price_time['Sep 2022':] # 测试数据集（最后12个月）
price_train = price_time[:'Aug 2022'] # 训练数据集（前面的数据）

对于我们的数据集，时间跨度从2017年1月到2023年8月。在这种情况下，我们选择了最后12个月的数据作为测试数据集，选择了前面的数据作为训练数据集。有人可能会想选择12个月的原因是什么。虽然这种选择可能有些随意，但它有实际意义。12个月的时间窗口足够长，让我们观察到我们预测的行为。另一方面，这段时间足够短，使我们的预测不会与实际数据相差太多。

时间序列建模
根据我们之前获得的信息，我们可以构建一个时间序列模型。最终，我们的目标是利用对模型的理解来预测数据集在未来的行为。这个过程被称为预测。

有许多不同的时间序列建模方法。经典方法包括移动平均、指数平滑和ARIMA（自回归综合移动平均），而更高级的技术包括Facebook Prophet和LSTM网络。在这个特定的背景下，经典方法效果相当不错，而在这些方法中，季节性自回归综合移动平均（SARIMA）是表现最佳的方法。

要构建一个SARIMA模型，我们需要确定ARIMA模型的参数p、d和q，以及季节性的参数P、D、Q和s。有关这些参数的解释可以在这里和这里阅读。虽然有更严谨的方法来确定这些参数，但在这个特定的实例中，我们使用了网格搜索的实用方法，并根据最低的AIC（阿卡伊信息准则）得分选择模型。

from statsmodels.tsa.statespace.sarimax import SARIMAX
s = 12 # 季节性由12个月的循环组成
for p in range(1, 4):
    for d in range(0, 3):
        for q in range(0, 4):
            for P in range(1, 4):
                for D in range(0, 3):
                    for Q in range(0, 4):
                        model = SARIMAX(price_train, order=(p, d, q), seasonal_order=(P, D, Q, s))
                        results = model.fit()
                        print(f"ARIMA{p},{d},{q}-SARIMA{P},{D},{Q}-{s} - AIC: {results.aic}")
# 注意，运行网格搜索可能会占用大量计算资源并花费一些时间
# 所以得到所需的参数后记得关闭这些代码

对于这组数据，我们得到的最佳参数是ARIMA参数为3、1和3，季节性参数为1、1、3和12。在以下代码块中，我们可以看到如何使用Statsmodels库将训练数据集拟合到模型中。

p, d, q = 3, 1, 3
P, D, Q, s = 1, 1, 3, 12
sarima_model = SARIMAX(price_train, order=(p, d, q), seasonal_order=(P, D, Q, s))
sarima_result = sarima_model.fit()

预测
得到模型后，我们需要检查残差图。如下图所示，由于某些原因，我仍然无法解释，第一个和第13个索引处的残差值始终远离零，无论如何调整模型的参数。尽管我多次尝试，这个问题仍然没有解决。

sarima_result.resid.plot()

    
    


然而，剩下的数据行为正常。因此，为了解决这个问题，我决定使用第13个索引之后的数据来分析模型。这种实用的方法使我能够使用数据集更稳定的部分，并避免极端异常值的影响。

sarima_result.resid[13:].plot()

所以，在构建SARIMA模型之后，下一步是对训练和测试数据集进行预测，并将其与实际数据集一起绘制，以进行可视化评估。

模型评估

从上面的可视化结果可以看出，预测值（来自训练和测试数据集）与实际值之间有明显的相似性。预测值也在模式和趋势方面紧密跟随实际值。因此，通过视觉检查，我们可以得出结论，我们的SARIMA模型可以很好地预测训练数据集中的数据。此外，这种准确性还延伸到测试数据集。这意味着该模型在短期内（在我们的例子中是12个月）可以很好地推广到未见过的数据。

除了视觉检查外，我们还可以使用某些指标来评估我们的时序模型的准确性。其中之一是平均绝对误差（MAE）。我们通过计算预测值与数据集中的真实值之间的平均差异的绝对值来计算MAE。从本质上讲，MAE衡量了我们的模型与真实值的平均接近程度。较低的MAE意味着我们的模型的预测与实际数据接近，而较高的MAE意味着我们的模型存在较大误差。为了计算MAE，我们可以使用Scikit-learn库中提供的函数。

from sklearn.metrics import mean_absolute_error as MAE # 不要忘记反转我们的对数转换 MAE_train = MAE(price_train[13:].apply(reverse_transform), train_predict[13:].apply(reverse_transform)) MAE_test = MAE(price_test.apply(reverse_transform), test_predict.apply(reverse_transform)) print(f'训练数据集的MAE = {MAE_train:.2f}') print(f'测试数据集的MAE = {MAE_test:.2f}')

从以上计算结果可以看出，训练数据集和测试数据集的MAE分别为5748新加坡元和7549新加坡元。为了对比，我们可以将这些值与时间序列数据的平均值相比较，平均值为481,070新加坡元。这意味着MAE值与平均每月价格相比较小。这些值表明我们的模型与真实值之间只存在较小的偏差。

此外，测试数据集的MAE仅略高于训练数据集的MAE。这种微小的差异进一步证实了我们的时序模型在短期内可以很好地适应未见过的数据。

待续……

最后，我们已经完成了我们的时序模型，我们准备进入这个故事的第二部分。在下一部分中，我将向您展示如何利用我们从时间序列分析中获得的见解，并将其用于改善回归模型的稳健性。

进一步阅读：

[1] 房屋发展局，HDB历史与城镇规划，2023。[在线]. 可获取: https://www.hdb.gov.sg/cs/infoweb/about-us/history [访问日期：2023年10月30日]

[2] J. Brownlee，如何使用Python创建ARIMA时间序列预测模型，2020。[在线]. 可获取: https://machinelearningmastery.com/arima-for-time-series-forecasting-with-python/ [访问日期：2023年10月30日]

[3] J. Brownlee，《Python时间序列预测的SARIMA简介》，2019年。[在线]. 可获得: https://machinelearningmastery.com/sarima-for-time-series-forecasting-in-python/ [访问日期：2023年10月30日]