随着COVID-19大流行的持续,许多国家正面临疫苗短缺的问题。通过对全球COVID-19疫苗接种进度的调查,了解到全球已接种超过9.5亿剂疫苗,平均每100人中有12人接种。因此,对各国疫苗接种率进行预测分析,以确定特定国家何时能实现全民接种,显得尤为重要和有用。
本文将介绍如何使用时间序列分析来预测不同国家的COVID-19疫苗接种率,并分享用于创建预测模型的Python代码。在阅读本文的同时,可以打开Jupyter Notebook,使用Auto-Regressive Integrated Moving Average(ARIMA)进行简单的时间序列分析。
ARIMA是“自回归积分滑动平均”的缩写。由于本文面向时间序列分析的初学者,因此从简单的分析开始是更好的选择。尽管称之为“简单”,但它是一个非常有用的模型类别,可以帮助预测大量时间序列,其中描述了给定时间序列从过去值(称为滞后和滞后预测误差)的时间序列,并用于预测值(预测未来值)。
时间序列中的过去值称为AR项(自回归项),预测误差称为MA项(移动平均项)。本次分析的数据来自牛津大学“Our World in Data”项目。时间序列的第一天是该国家开始接种疫苗的那一天。本次分析使用了大约前75天的数据。用于ARIMA时间序列分析的样本数据可以从以下GitHub仓库下载。
让开始吧...
import pandas as pd
# 加载数据集
df = pd.read_csv('share-people-vaccinated-covid.csv', index_col='Day', parse_dates=True)
print('数据形状', df.shape)
df.head()
df
# 由于数据集中没有唯一国家标识符,可以添加如下:
df = df.sort_values('Entity')
df['eid'] = (df.groupby(['Entity']).cumcount() == 0).astype(int)
df['eid'] = df['eid'].cumsum()
ndf = df.loc[df['eid'] == 137]
# 检查威尔士的数据框架
sndf = ndf.sort_values('people_vaccinated_per_hundred')
sndf
从存储威尔士(或选择的国家)疫苗接种率的列中获取值。希望使时间序列平稳,使其不依赖于分析时间序列的时间。非平稳时间序列可能具有趋势和/或季节性。这些因素会影响时间序列的值。因此,在开始时间序列分析之前,有必要去除它们。要使非平稳时间序列变为平稳,可以进行差分操作以稳定序列的均值,计算连续的差分(即,连续观测之间的差值或从当前观测值中减去过去的一个值)。这有助于去除(或最小化)趋势和季节性。
有时,时间序列的一阶差分可能无法去除所有趋势和季节性,需要再次对数据进行差分(时间序列的二阶差分)以获得平稳时间序列。在分析自相关图后,进行了二阶差分以使时间序列平稳。
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
plt.rcParams.update({'figure.figsize':(9,7), 'figure.dpi':120})
# 原始序列
fig, axes = plt.subplots(3, 2, sharex=True)
axes[0, 0].plot(val); axes[0, 0].set_title('原始序列')
plot_acf(val, ax=axes[0, 1])
# 一阶差分
axes[1, 0].plot(np.diff(val)); axes[1, 0].set_title('一阶差分')
plot_acf(np.diff(val), ax=axes[1, 1])
diffval = np.diff(val)
# 二阶差分
axes[2, 0].plot(np.diff(diffval)); axes[2, 0].set_title('二阶差分')
plot_acf(np.diff(diffval), ax=axes[2, 1])
plt.show()
没有季节性的ARIMA模型表示为ARIMA(p,d,q),其中p是AR项的数量,d是需要的差分次数(在这个例子中,假设需要d=2)以使时间序列变为平稳,q是MA项的数量。然后找到最佳的ARIMA模型来拟合。在这里,观察Akaike信息准则(AIC)以确定最佳ARIMA模型的阶数。
要找到差分的最佳值(之前假设d=2),可以使用增强的Dickey-Fuller(ADF)测试。
from statsmodels.tsa.arima_model import ARIMA
import pmdarima as pm
model = pm.auto_arima(val, start_p=1, start_q=1,
test='adf', # 使用adftest找到最佳'd'
max_p=3, max_q=3, # 最大p和q
m=1, # 序列的频率
d=None, # 让模型确定'd'
seasonal=False, # 无季节性
start_P=0,
D=0,
trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
print(model.summary())
选择的最佳模型是:ARIMA(p,d,q) = ARIMA(3,2,2)。正如所假设的,需要进行二阶差分以使时间序列在这种情况下变得平稳。
作为下一步,执行残差的诊断检查如下:
model.plot_diagnostics(figsize=(7,5))
plt.show()
已经完成了所有必要的步骤。然后现在可以预测疫苗接种率了。
# 预测
n_periods = 50
fc, confint = model.predict(n_periods=n_periods, return_conf_int=True)
index_of_fc = np.arange(len(val), len(val)+n_periods)
# 制作绘图系列
fc_series = pd.Series(fc, index=index_of_fc)
lower_series = pd.Series(confint[:, 0], index=index_of_fc)
upper_series = pd.Series(confint[:, 1], index=index_of_fc)
# 绘图
plt.plot(val)
plt.plot(fc_series, color='darkgreen')
plt.fill_between(lower_series.index,
lower_series,
upper_series,
color='k', alpha=.15)
plt.title("疫苗接种率预测")
plt.show()