时间序列大作业报告＋代码（以新冠新增人数为例）

作业要求：

1：本作业允许你下载任何数据集，并进行时间序列分析。

2：要求下载的数据集至少200个数据

3：评分标准:

a:报告格式是否规范;

b:是否涵盖所学的方法(数据预处理，平稳建模，确定性分析，随机分析等)

c:数学表达式是否正确;

d:结果是否进行分析;数据集是否为最新数据等。

4:本作业要求介绍你下载的数据集的背景、目的、意义，及出处。作业内容限制在 20 页左右。

因为是帮朋友做的，自己没上过这门课，报告中的数学描述部分可能有不少错误，所以请慎重copy。 

报告如下；代码请翻到最后，如有疑问，请邮箱1902946954@qq.com。

摘要

        新冠疫情逐渐进入新阶段，通过研究新冠每日新增人数的统计，可以帮助政府和公共卫生机构更好地准备资源和应对措施，以便在下一次新冠疫情或者其他公卫紧急事件爆发时更快地做出反应并控制疫情的传播。因此本文研究2020年3月1日到2023年3月1日的美国每日新增感染人数。

        首先对数据进行预处理。首先对数据进行7日滑动平均来去掉0值，在统计数据的各特征统计量，并使用时序图、自相关图、单位根检验法进行数据的平稳性检验分析，最后使用纯随机性检验。

        然后进行确定性分析。使用线性趋势拟合法和指数平滑法。在综合分析，使用乘法模型对数据进行分解，绘制季节指数，趋势拟合，残差图，最后再对随机项进行纯随机性检验。我们会发现该数据较为复杂，并不能用简单的线性或二次函数的趋势，乘法模型来拟合。

        接着进行随机性分析。首先通过一阶-多步差分提取线性趋势和季节周期，并使用ARIMA(1,1,1)模型来拟合并做十步的预测，其中ARIMA(1,1,1)模型的rmse和mae值都较小，拟合效果较好。

目录

1.数据介绍... 3

1.1 数据来源... 3

1.2数据介绍... 3

1.3 研究意义... 3

2 数据预处理... 4

2.1数据读取与可视化... 4

2.2 数据统计量... 5

2.3 平稳性检测... 5

2.3.1 时序图检验... 5

2.3.2自相关图检验... 6

2.3.3 单位根检验... 6

2.3.4 纯随机性检验... 7

2.4 总结... 7

3 确定性分析... 7

3.1 趋势分析... 7

3.1.1 趋势拟合法... 8

3.1.2 平滑法... 9

3.1.2.1移动平均法... 9

3.1.2.2 指数平滑法... 11

3.2 综合分析... 11

3.2.1乘法模型... 12

4 随机分析... 15

4.1 差分运算... 15

4.1.1 一阶差分... 16

4.1.2 一阶-12步差分... 16

4.2 ARIMA模型... 17

5总结... 19

6参考文献... 19

1.数据介绍

1.1 数据来源

    本文数据为美国从2020年3月1日到2023年3月14日的每日新确认的新冠肺炎病例人数。数据来自网站“Our World In Data”，网址为：COVID-19 Data Explorer - Our World in Data

    Our World in Data是一个由牛津大学马丁学派全球变化研究所创建的学术数据网站。该网站致力于收集、整理和呈现有关全球问题的数据和研究成果，以便更好地了解和解决这些问题。它的数据涵盖了许多领域，包括经济、环境、健康、能源、食品和农业等。这些数据被呈现在易于理解和使用的交互式图表和可视化工具中，让人们可以更加直观地了解全球问题的现状和趋势。

1.2数据介绍

    该新冠数据集包含世界各国2020年3月1日到2023年3月14日有关新增患者、死亡人数、ICU人数、疫苗接种情况、各变异病毒传播人数等数据。我们取其中“美国每日新增新冠肺炎病例人数”作为研究对象。

1.3 研究意义

    通过对时间序列数据的分析，我们可以发现疫情的高峰期和低谷期，以及疫情发展的趋势和速度。这些信息可以帮助政府和公共卫生机构更好地制定疫情防控策略，包括隔离措施、疫苗接种计划等。

        此外，时间序列分析还可以帮助我们预测未来的疫情形势。通过对历史数据的分析，我们可以建立预测模型，预测未来的感染人数和疫情发展趋势。这可以帮助政府和公共卫生机构更好地准备资源和应对措施，以便在疫情爆发时更快地做出反应并控制疫情的传播。[3]

2 数据预处理

2.1数据读取与可视化

    用python绘制时序图，横坐标为时间（因日期文本太长，因此报告中所有图横坐标的日期都用自然数来代替）纵坐标为每日新增人数，如图1：

 

图1：美国新冠每日新增人数数据时序图

    结果分析：我们会发现横坐标10000到12000的数据中有较多值为0（因为新冠后期医疗机构统计数据的频率变慢，因此数据不能覆盖每一天）因此我们用python对数据进行七天滚动平均来进行线性插值。滑动平均后的人数存在小数，我们用round函数取整。

对处理后的数据可视化如图2：

 

图2：经滑动平均处理的时序图

2.2 数据统计量

    用python的numpy库中的函数进行统计量分析，计算该数据集的最大值、最小值、均值、标准差和方差，如表1：

变量名	最小值	最大值	均值	标准差	方差
人数	0	809735	84304.18	110221.57	12148795984

   

2.3 平稳性检测

    平稳性是时间序列分析中的一个重要概念，它指的是时间序列的均值、方差和自协方差不随时间的推移而发生显著变化的性质。如果一个时间序列是平稳的，那么它的统计特征（如均值、方差等）在不同时间段内应该是相似的，这种性质对于时间序列的建模和预测非常重要。

    下面通过不同方法进行平稳性检测。

2.3.1 时序图检验

    分析：我们在2.1中已得到时序图（如图2）。从时序图上直接看，数据总体上为先上升后下降的趋势。有较明显以年和季节为周期的周期性，各统计值随着时间推移发生变化，因此没有平稳性特征。

2.3.2自相关图检验

 图3：自相关图

    该自相关图在不同的滞后阶数下波动较大，且在较多滞后阶数上的自相关系数显著不为 0，那么时间序列可能存在趋势或季节性，不是平稳的。

2.3.3 单位根检验

    我们使用 statsmodels 库的 adfuller() 函数进行单位根检验，将结果存储在 result变量中。最后，打印出检验结果，包括检验统计量的值和关键值：

ADF Statistic: -2.144931

p-value: 0.45566

Critical Values:

1%: -3.436

5%: -2.864

10%: -2.568

    分析：ADF Statistic 的值为 -2.144931，大于 10% 的临界值-2.568，不能拒绝原假设，即该序列存在单位根非平稳性。

        同时，p 值为 0.45566，大于 0.05 的显著性水平，不能拒绝原假设，即该序列存在单位根非平稳性。

2.3.4 纯随机性检验

    使用 statsmodels 库的 runs() 函数进行纯随机性检验，将结果存储在 result 变量中。最后，打印出检验结果，包括检验统计量（Runs）和 p 值，结果如下：

Runs: 4

P-value: 0.485727

分析：在显著性水平为 0.05 的情况下，p值大于 0.05 的显著性水平。即我们不能拒绝原假设，认为该序列是具有纯随机性的。

   

2.4 总结

通过时序图，自相关系数图，单位根检验，纯随机性检验分析，可以得出美国每日新增新冠人数数据序列为非平稳且为非纯随机性序列。需要进行后续的确定性分析和随机性分析建立拟合模型。

3 确定性分析

3.1 趋势分析

     3.1.1 趋势拟合法

    通过时序图的分析可以得到，该序列的长期趋势有大致的先增后减的趋势特征，可以初步使用二阶函数拟合，具体公式如s下：

Y=ax2+bx+c

    我么使用 Python 中的 NumPy 和 SciPy 库来进行二次函数拟合，结果如下：

拟合结果：a=-0.393, b=507.685, c=-29926.160

    拟合效果如下图：

 

图4：二次函数拟合效果图

    结果分析：提取了大致先增后减的趋势，但是该数据不平稳且较为复杂，简单的二次函数显然不能较好地拟合数据，需要用更加强大的模型进行进一步拟合。

3.1.2 平滑法

    平滑法是一种常用的趋势分析方法，它的基本思想是通过对时间序列数据进行加权平均，去除随机扰动，从而得到趋势分量。其中最常用的平滑法包括移动平均法（Moving Average，MA）和指数平滑法（Exponential Smoothing，ES）。

移动平均法是一种简单的平滑法，它的思想是将一定时期内的数据取平均值作为平滑后的值，以消除随机波动，突出趋势变化。移动平均法的优点是计算简单，易于理解和应用。但是，它对于非常规数据（如季节性波动等）的处理能力较弱，且滞后效应明显。

指数平滑法是一种加权平均法，它的思想是对历史数据进行加权平均，使得近期数据的权重较大，远期数据的权重较小，以适应数据的变化。指数平滑法的优点是对近期数据的反应较快，而且适用于非常规数据的处理。但是，指数平滑法对于噪声数据的处理能力较弱，且需要选择合适的平滑系数和初始值。

3.1.2.1移动平均法

假定在一个比较短的时间间隔里，序列值之间的差异主要是由随机波动造成的。根据这种假定，我们可以用一定时间间隔内的平均值作为某一期的估计值。我们在移动平均模型中将间隔长度设置为30和7，具体公式如下：

xt=1k(xt+xt-1+⋯+xt-n+1)

    得到的结果可视化后如下：

 

图5：移动平均法（窗口为30天）

 

图6：移动平均法（窗口为7天）

分析：移动平均拟合效果很好，特别是窗口选取的较小时，但是不便于用于预测，而且在窗口选取较大时有明显的滞后效应。

3.1.2.2 指数平滑法

    我们用python pandas库中的series.ewd函数进行指数平滑法的拟合，窗口span选为30。拟合结果为下图：（蓝色为原数据，橙色为滑动平均、绿色为指数平滑）

 

图7：指数平滑法法（窗口为7天）

    分析：指数平滑法的优点是对近期数据的反应较快，但是对于噪声的去除没有移动平均好，因此一定程度上移动平均总体上对该数据拟合的更好。

   3.2 综合分析

    对于一个同时受到趋势和季节因素影响的序列，首先要判断趋势和季节之间的相互影响关系。

1. 如果季节波动的振幅不受到趋势变动的影响，那么季节与趋势之间通常没有交互影响关系，可以认为是可加关系。
2. 如果季节波动的振幅随着趋势的变化而变化，那么季节与趋势之间意味着有交互影响关系，可以认为是乘法关系。

3. 

 

    分析：通过观察该时序图可以发现：周期的振幅随着货运量趋势的改变而改变，也就是说季节与趋势之间有交互影响关系，本文所以尝试使用乘法模型。

     3.2.1乘法模型

    我们使用python的seasonal_decompose()函数返回的result对象是一个具有四个属性的对象：

result.observed：原始时间序列数据。

result.trend：时间序列数据的趋势分量。

result.seasonal：时间序列数据的季节性分量。

result.resid：时间序列数据的残差分量。

    趋势分量、季节性分量、残差分量可视化后分别为下面三幅图：

 

图8：趋势分量

图9：季节性分量

图10：残差序列

    将四个属性信息整合后如下图：

 

图11：乘法模型分析整合

    我们对乘法模型进行纯随机性检验。检验的一种常用方法是检查残差序列是否具有自相关性。如果残差序列是纯随机的，那么它们应该是独立的，不应该存在自相关性。如果残差序列存在自相关性，则表示模型未能捕捉到数据中的某些模式。用python的plot_acf()函数画出自相关图，如下图：

 

图12：乘法模型残差序列纯随机检验

    可以看到，置信区间由蓝色的上下两条水平线表示。有较多点在置信区间外，因此说明残差序列并不是纯随机的，因此乘法模型并未很好地拟合该数据。

4 随机分析

4.1 差分运算

    差分方法是一种非常简便、有效的确定性信息提取方法，当序列蕴含着显著的线性趋势，一阶差分就可以实现趋势平稳。

4.1.1 一阶差分

    我们尝试用一阶差分来分析新冠每日新增序列。具体公式为：

∇t=xt-xt-1

    下面对序列进行一阶差分，结果绘制时序图如下：（红色线为原数据，蓝色线为一阶差分后的序列）

 

图13：一阶差分结果

    结果分析：因为原数据并非线性递增趋势，因此一阶差分数据并不平稳，用一阶差分效果并不好。

     4.1.2 一阶-12步差分

对于蕴含着固定周期的序列进行步长为周期长度的差分运算，通常可以较好地提取周期信息。而该数据有以12个月为周期的趋势，所以使用一阶差分提取线性趋势，加上以年为单位的阶数（365）差分提取季节趋势。

    结果如下图：（绿色线为原数据，蓝色线为一阶差分后的序列）

 图14：一阶-365步差分结果

        分析：在一阶差分提取线性趋势基础上，使用以年为单位的步数差分提取周期性，绘制时序图后发现序列相对一阶差分的时序图更加平稳。

4.2 ARIMA模型

    ARIMA模型（Autoregressive Integrated Moving Average，自回归积分移动平均模型）是一种经典的时间序列模型，常用于预测和分析时间序列数据。

    ARIMA模型的一般形式是ARIMA(p,d,q)，其中p是自回归项数，d是差分阶数，q是移动平均项数。ARIMA模型可以根据时间序列数据的特点，动态地确定p、d、q的取值，使得模型最优化地拟合时间序列数据。[4]

ARIMA模型通常包括三个步骤：

1：模型识别：根据时间序列的自相关图和偏自相关图，确定ARIMA模型的阶数。

2：参数估计：利用极大似然估计法等方法，估计ARIMA模型的参数。

3：模型检验：对ARIMA模型进行残差分析和预测检验，检验模型的拟合优度和预测能力。

我们通过python statsmodels库来拟合ARIMA模型，确定结束，估计参数，进行预测，并且最后预测十个步长的值，并使用均方根误差（RMSE）和平均绝对误差（MAE）来衡量模型的精度。RMSE和MAE越小，模型的精度越高。

Arima(1,1,1)模型拟合结果如下图：

 图15：arima模型拟合效果图

表2：arima精度衡量

均方根误差（RMSE）	平均绝对误差（MAE）
5004.524	1996.604

    分析：我们可以发现红色线条是拟合值，蓝色线条是原数据，两条线基本重合。并且均方根误差和平均绝对误差的值相比于数据都较小，因此arima(1,1,1)拟合效果较好。

源码如下：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.optimize import curve_fit
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoLocator, MaxNLocator
import statsmodels.api as sm
from sklearn.metrics import mean_squared_error, mean_absolute_error

# 读取Excel文件数据
df = pd.read_excel("数据集路径")

# 获取第一列的数据和数据长度
y_data = df.iloc[:, 1]
y_data = round(y_data)
n = len(y_data)
# # 生成横坐标数据
x_data = range(1, n+1)
#统计量
# mean=np.mean(y_data)
# print(mean)
# std=np.std(y_data)
# print(std)
# var=np.var(y_data)
# print(var)
# max=np.max(y_data)
# print(max)

# 绘制折线图
# x_2stage=list(x_data)
# y_2stage= [0 for i in range(n)]
# for i in x_2stage:
#     y_2stage[i-1]=x_2stage[i-1]*x_2stage[i-1]*(-0.393)+x_2stage[i-1]*507.685-29926.16

#滑动平均法
# y_average = pd.Series(y_data)
# ma = y_average.rolling(window=30).mean()
# es = y_average.ewm(span=30).mean()
# 绘制折线图
# plt.plot(x_data, y_data)
# plt.plot(x_data, ma)
# plt.plot(x_data, es)

# plt.plot(x_data, y_2stage)
# 设置横轴标签和纵轴标签
# plt.xlabel('date')
# plt.ylabel('New per day')

# 显示图形
# plt.show()

# y_data_mul=y_data
# y_data_mul[y_data_mul == 0] = 1
# y_data_mul = pd.Series(y_data_mul)
# index = pd.date_range(start='20200301', periods=len(y_data_mul), freq='d')
# y_data_mul.index = index
# 对数/据进行季节分解
# result = seasonal_decompose(y_data_mul, model='multiplicative')
# result.plot()
# plt.show()
# result.observed.plot()
# plt.show()
# result.trend.plot()
# plt.show()
# result.seasonal.plot()
# plt.show()
# result.resid.plot()
# plt.show()
# 获取残差序列
# residuals = result.resid.dropna()

# 绘制ACF和PACF图
# plot_acf(residuals, lags=20)
# plt.show()
# plot_pacf(residuals, lags=20)
# plt.show()



# #自相关图
# # 计算数组的自相关系数
# autocorr = np.correlate(y_data, y_data, mode='full')
# # 画出自相关图
# plt.plot(autocorr)
# # 设置图形标题和轴标签
# plt.title('Autocorrelation Plot')
# plt.xlabel('Lag')
# plt.ylabel('Autocorrelation')
# # 显示图形
# plt.show()

#求单位根检验
# result = adfuller(y_data)
# print('ADF Statistic: %f' % result[0])
# print('p-value: %f' % result[1])
# print('Critical Values:')
# for key, value in result[4].items():
#     print('\t%s: %.3f' % (key, value))

#纯随机性检验
# ljungbox_result = acorr_ljungbox(y_data, lags=20)  # 返回统计量和p值，lags为检验的延迟数
# print(ljungbox_result)

# def quadratic_function(x, a, b, c):
#     return a * x**2 + b * x + c
# # 构造数据点
# x_data = list(range(1, n+1))

# # 使用 curve_fit 函数进行拟合
# popt, pcov = curve_fit(quadratic_function, x_data, y_data)
# # 输出拟合结果
# a_fit, b_fit, c_fit = popt
# print(f"拟合结果：a={a_fit:.3f}, b={b_fit:.3f}, c={c_fit:.3f}")

#差分运算
# 对数据进行一阶差分
# diff_data = y_data.diff(periods=1)
# # 对差分后的数据进行365阶差分
# diff365 = diff_data.diff(periods=365)
#
# # 输出一阶差分后的数据
# # 创建画布和子图
# fig, ax1 = plt.subplots()
#
# # 绘制第一个数据集
# ax1.plot(diff365, color='blue')
# ax1.set_xlabel('X轴')
#
# # 创建第二个Y轴
# ax2 = ax1.twinx()
#
# #绘制第二个数据集
# ax2.plot(y_data, color='green')
#
# # 自适应选择坐标轴刻度
# ax1.yaxis.set_major_locator(AutoLocator())
# ax2.yaxis.set_major_locator(MaxNLocator(integer=True))
#
# # 显示图形
# plt.show()


#arima
#y_data_arima = pd.Series(y_data)
#index = pd.date_range(start='20200301', periods=len(y_data_arima), freq='d')
# 拟合ARIMA模型
model = sm.tsa.ARIMA(y_data, order=(1, 1, 1)).fit()

# 进行未来10个时间步长的预测
forecast = model.forecast(steps=10)

# 计算预测误差和残差
error = y_data - model.fittedvalues
residuals = model.resid

# 计算RMSE和MAE
rmse = mean_squared_error(y_data, model.fittedvalues, squared=False)
mae = mean_absolute_error(y_data, model.fittedvalues)
print(rmse)
print(mae)
# 将模型结果可视化
fig, ax = plt.subplots(figsize=(12, 6))

# 绘制原始数据
ax.plot(y_data.index, y_data.values, label='Actual')

# 绘制拟合值
ax.plot(model.fittedvalues.index, model.fittedvalues.values, color='lightcoral', label='Fitted')

# 绘制预测值
ax.plot(forecast.index, forecast.values, color='palegreen', label='Forecast')

# 设置图例和标题
ax.legend()
ax.set_title('ARIMA Model Forecast')

# 显示图形
plt.show()

 请自己对照着看要哪一部分，来进行注释和取消注释。

报告中的数学方面有可能存在问题，请慎重copy。