【机器学习实战中阶】使用SARIMAX，ARIMA预测比特币价格，时间序列预测
数据集说明
比特币价格预测（轻量级CSV）关于数据集
致谢
这些数据来自CoinMarketCap，并且可以免费使用该数据。
https://coinmarketcap.com/ 数据集:链接: 价格预测器 源代码与数据集
算法说明
SARIMAX（Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors）和 ARIMA（AutoRegressive Integrated Moving Average）都是用于时间序列分析和预测的强大工具，但它们之间存在一些关键区别。这些差异主要体现在模型的灵活性、处理季节性和外生变量的能力上。
1. 季节性成分
-
ARIMA：适用于非季节性的数据，即没有明显周期模式的时间序列。它通过自回归（AR）、差分（I）和移动平均（MA）三个部分来捕捉时间序列中的趋势和平稳性特征。
-
SARIMAX：是ARIMA模型的一种扩展，特别设计用来处理具有季节性模式的数据。除了包含ARIMA的所有组件外，SARIMAX还增加了额外的参数来描述季节性影响，如季节性自回归（SAR）、季节性差分（SI）和季节性移动平均（SMA）。因此，SARIMAX可以更好地捕捉时间序列中的长期周期波动。
2. 外生变量的支持
-
ARIMA：传统上的ARIMA模型仅依赖于时间序列本身的历史值来进行预测，不直接支持外部变量的影响。
-
SARIMAX：允许引入外生变量（exogenous variables），即那些可能影响目标时间序列但不属于该序列内部结构的因素。例如，在预测销售量时，节假日、促销活动等外部因素可以通过外生变量的形式加入到模型中，以提高预测精度。
3. 数学表示
-
ARIMA(p, d, q)：其中(p)代表自回归项的阶数，(d)代表差分阶数，(q)代表移动平均项的阶数。
-
SARIMAX(p, d, q) × (P, D, Q)_m：这里不仅包含了非季节性的(p), (d), (q)参数，还有额外的季节性参数(P)（季节性自回归阶数）、(D)（季节性差分阶数）、(Q)（季节性移动平均阶数），以及(m)为季节周期长度。
4. 应用场景
-
ARIMA：适合用于没有明显季节性的平稳或经过适当变换后变得平稳的时间序列数据。
-
SARIMAX：更适合那些展示出显著季节性和/或受其他外部因素影响较大的时间序列数据。比如，零售业销售额通常会显示出明显的月度或年度季节性变化，而天气条件、经济指标等外部因素也可能对销售产生重要影响。
综上所述，选择使用ARIMA还是SARIMAX取决于您的具体需求和数据特性。如果您的时间序列包含明显的季节性并且您有可用的外生变量信息，那么SARIMAX可能是更好的选择；否则，对于简单的非季节性数据，ARIMA已经足够强大且易于实现。
接下来我们从详细代码方案开始
# 这个 Python 3 环境附带了许多有用的分析库
# 它是由 kaggle/python 容器镜像定义的：https://github.com/kaggle/docker-python
# 例如，这里有几个有用的数据包可以加载
import numpy as np # 线性代数库
import pandas as pd # 数据处理库，用于处理 CSV 文件的输入输出（如：pd.read_csv）
import seaborn; seaborn.set() # seaborn 是一个基于 matplotlib 的数据可视化库，这里通过 seaborn.set() 来设置默认的绘图风格
from sklearn.metrics import mean_squared_error # 用于计算平均平方误差，一种常见的评估回归预测模型的方式
# 输入数据文件位于 "../input/" 目录下。
# 例如，运行此代码（通过点击运行或按下 Shift+Enter）将列出输入目录中的文件
# 任何你写入当前目录的结果都会作为输出保存。
import matplotlib.pylab as plt # 导入 matplotlib 的子模块 pyplab（pylab），用于绘制图像
%matplotlib inline # 魔法命令，使图表内嵌在 jupyter notebook 中
from matplotlib.pylab import rcParams # 从 pyplab 中导入 rcParams，用于设置图表的参数
rcParams['figure.figsize'] = 15, 6 # 设置图表的默认尺寸为宽度15英寸，高度6英寸
from statsmodels.tsa.stattools import adfuller # 引入 adfuller 模块，用于进行时间序列平稳性检验
from statsmodels.tsa.seasonal import seasonal_decompose # 引入 seasonal_decompose 模块，用于对时间序列数据进行分解
from statsmodels.tsa.stattools import acf, pacf # 引入 acf 和 pacf 模块，用于计算时间序列的自相关和偏自相关函数
from statsmodels.tsa.arima_model import ARIMA # 引入 ARIMA 模型，用于时间序列数据的预测
加载数据并研究行和列的值
# 加载 bitcoin_price_Training - Training.csv 文件，作为一个包含时间序列数据的 DataFrame 对象
data = pd.read_csv("../input/bitcoin_price_Training - Training.csv")
# 打印数据表的前五行，以查看数据格式
print(data.head(5))
# 打印数据表的最后五行，以查看最近的数据
print(data.tail(5))
# 查看每一列的数据类型
data.dtypes
# 获取数据表的信息，包括行数、列数、列名、数据类型等
data.info()
# 生成数据表的描述统计信息，如计数、均值、最小值、最大值等
data.describe()
# 重新加载数据，设置日期列为索引
data = pd.read_csv("../input/bitcoin_price_Training - Training.csv", index_col='Date')
# 再次打印前五行，检查日期是否正确设置为索引
print(data.head(5))
# 检查数据表的信息，确认索引列已正确设置
data.info()
# 将索引从字符串格式转化为日期时间格式
data.index = pd.to_datetime(data.index)
# 打印索引，确认转换成功
print(data.index)
# 再次打印数据表的前五行，查看转换后的数据
data.head(5)
将数据集按时间顺序排序，从最旧到最近
# 将数据按索引排序，从最旧的数据到最近的数据
data = data.sort_index()
# 打印排序后的前五行数据，确认排序成功
data.head()
绘制比特币收盘价的折线图
# 绘制比特币每日收盘价的折线图
data['Close'].plot()
# 设置 y 轴标签为"每日比特币价格"
plt.ylabel("每日比特币价格")
为了建立模型，我将侧重于预测比特币的收盘价，因此创建一个新的对象，排除其他列
# 从数据中提取收盘价列，创建一个新的 Series 对象
data = data['Close']
# 将原始数据按周重采样并累计，创建周比特币价格数据
weekly = data.resample('W').sum()
# 绘制比特币每周价格的折线图
weekly.plot()
# 设置 y 轴标签为"周比特币价格"
plt.ylabel('周比特币价格')
按年份统计并绘制比特币平均收盘价格
# 通过数据中的日期索引按年份进行分组，并计算每年的均值
by_year = data.groupby(data.index.year).mean()
# 绘制按年份分组的比特币平均收盘价折线图
by_year.plot()
按星期几统计并绘制比特币平均收盘价格
# 通过数据中的日期索引按星期几进行分组，这里是累加每个星期几的所有数据
by_weekday = data.groupby(data.index.dayofweek).sum()
# 将原始 0-6 的数字索引替换为对应的星期几名称
by_weekday.index = ['Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun']
# 绘制按星期几分组的比特币价格累加图
by_weekday.plot()
# 通过数据中的日期索引按星期几进行分组，这一次是计算每个星期几的收盘价平均值
by_weekday = data.groupby(data.index.dayofweek).mean()
# 再次替换数字索引为星期几名称
by_weekday.index = ['Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun']
# 绘制按星期几分组的比特币平均收盘价折线图
by_weekday.plot()
按年内第几天统计并绘制比特币平均收盘价格
# 通过数据中的日期索引按年内第几天进行分组，并计算每一天的均值
by_day = data.groupby(data.index.dayofyear).mean()
# 绘制按日内第几天分组的比特币平均收盘价折线图
by_day.plot()
按月份统计并绘制比特币平均收盘价格
# 通过数据中的日期索引按月份进行分组，并计算每个月份的均值
by_month = data.groupby(data.index.month).mean()
# 绘制按月份分组的比特币平均收盘价折线图
by_month.plot()
按季度统计并绘制比特币平均收盘价格
# 通过数据中的日期索引按季度进行分组，并计算每个季度的均值
by_quarter = data.groupby(data.index.quarter).mean()
# 绘制按季度分组的比特币平均收盘价折线图
by_quarter.plot()
按季度再次分组绘制，这次仅为了展示分组操作
# 通过数据中的日期索引按季度进行分组
by_quarter = data.groupby(data.index.quarter)
# 这次我们只调用 plot 方法，不计算平均值，用于展示不同季度的数据分布
by_quarter.plot()
绘制按季度统计的比特币平均收盘价格，以及年度的对比图
# 计算每个季度的数据平均值
by_quarter_overall = data.groupby(data.index.quarter).mean()
# 绘制季度的比特币平均收盘价折线图
by_quarter_overall.plot()
按工作日与周末统计并绘制比特币平均收盘价格
# 使用 numpy 的 np.where 函数识别交易日为工作日（1-5）或周末（5-0）
weekend = np.where(data.index.weekday < 5, 'Weekday', 'Weekend')
# 将数据按工作日周末及年份分组，计算每组的平均收盘价
by_time = data.groupby([weekend, data.index.year]).mean()
# 创建一个包含1行2列的子图，用于同时展示工作日和周末的价格对比
fig, ax = plt.subplots(1, 2, figsize=(14, 5))
# 绘制工作日的平均收盘价格
by_time.loc['Weekday'].plot(ax=ax[0], title='工作日')
# 绘制周末的平均收盘价格
by_time.loc['Weekend'].plot(ax=ax[1], title='周末')
最后，绘制比特币的价格时间序列图
# 将数据赋值给 ts
ts = data
# 绘制比特币价格的时间序列图
plt.plot(ts)
以上代码主要用于加载比特币的历史价格数据，对数据进行预处理、转换、排序，并按不同的时间单位（如周、月、季及工作日与周末）对价格数据进行分组统计，最终绘制出这些统计结果的折线图。这有助于分析比特币价格随时间的变化趋势及其季节性特征。
继续干
这一部分内容的灵感来源于：https://www.analyticsvidhya.com/blog/2016/02/time-series-forecasting-codes-python/
from statsmodels.tsa.stattools import adfuller
def test_for_stationary(input_data):
r_mean = input_data.rolling(window = 7, center=False).mean() # 计算7日移动平均值，中心设为False
r_std = input_data.rolling(window = 7, center=False).std() # 计算7日移动标准差，中心设为False
# 绘制数据
given = plt.plot(input_data, color='blue', label='原始时间序列')
rolling_mean = plt.plot(r_mean, color='red', label='7日移动平均')
rolling_std = plt.plot(r_std, color='green', label='7日移动标准差')
plt.legend(loc='best')
plt.title('7日移动平均与标准差')
plt.show(block=False) # 显示图表
# 执行Dickey-Fuller检验
print('Dickey-Fuller检验结果：')
dftest = adfuller(input_data)
dfoutput = pd.Series(dftest[0:4], index=['检验统计量', 'p值', '使用滞后数量', '使用观察数量'])
for key, value in dftest[4].items():
dfoutput['临界值 (%s)'%key] = value # 添加临界值到输出中
print(dfoutput)
test_for_stationary(ts)
总结：
这段代码定义了一个名为test_for_stationary的函数，用于检查时间序列的平稳性。函数首先计算了7日移动平均值和7日移动标准差，并绘制了这两个指标与原始时间序列的图表，以便于直观地观察时间序列的平稳性。接着，代码执行了Dickey-Fuller检验，该检验是为了确定时间序列是否有单位根，从而判断其是否平稳。检验结果包括检验统计量、p值、使用的滞后数量、使用的观察数量以及不同显著性水平下的临界值。最后，函数调用test_for_stationary(ts)对时间序列ts进行了检验。
消除趋势的数据转换
ts_logtransformed = np.log(ts) # 对时间序列进行对数转换，消除趋势
plt.plot(ts_logtransformed) # 绘制对数转换后的时间序列图表
总结：
这段代码对时间序列ts进行了对数转换，目的是消除时间序列中的趋势成分。对数转换后的结果被绘制出来，以便于观察转换效果。
ts_logtransformed.head(10) # 查看转换后的时间序列的前10行数据
Rolling_average = ts_logtransformed.rolling(window = 7, center=False).mean() # 计算对数转换后的时间序列的7日移动平均值
plt.plot(ts_logtransformed, label='对数转换后的序列')
plt.plot(Rolling_average, color='red', label='7日移动平均')
plt.legend(loc='best')
plt.title('7日移动平均')
plt.show()
Rolling_average.head(10)
总结：
这段代码计算了对数转换后的时间序列的7日移动平均值，并绘制了原对数转换序列和7日移动平均值的图表，用以观察去除趋势的效果。
log_Rolling_difference = ts_logtransformed - Rolling_average # 对数转换后的序列减去7日移动平均值
log_Rolling_difference.head(10) # 查看差分后的前10行数据
log_Rolling_difference.tail(10) # 查看差分后的后10行数据
# 用0替换上述数据框中的NAN，以避免将来出现错误
log_Rolling_difference.dropna(inplace=True)
plt.plot(log_Rolling_difference)
test_for_stationary(log_Rolling_difference)
总结：
这段代码对对数转换后的时间序列进行了7日移动平均值的差分，通过减去7日移动平均值来进一步去除趋势。差分后的结果被绘制出，并检查了其平稳性。差分后的序列应该更加接近平稳性。
使用指数加权移动平均值来改进解决方案
expwighted_avg = ts_logtransformed.ewm(halflife=7, min_periods=0, adjust=True, ignore_na=False).mean() # 计算指数加权移动平均值，半衰期设为7天
plt.plot(ts_logtransformed, label='对数转换后的序列')
plt.plot(expwighted_avg, color='red', label='指数加权移动平均值')
plt.legend(loc='best')
plt.title('指数加权移动平均值')
plt.show()
总结：
这段代码计算了对数转换后的时间序列的指数加权移动平均值，半衰期设为7天。指数加权移动平均值可以更加灵活地捕捉时间序列的变化，因为它对较近的数据点赋予了更高的权重。绘制了原对数转换序列和指数加权移动平均值的图表，用以观察去除趋势的效果。
expwighted_avg.head(10) # 查看指数加权移动平均值的前10行数据
log_expmovwt_diff = ts_logtransformed - expwighted_avg # 对数转换后的序列减去指数加权移动平均值
test_for_stationary(log_expmovwt_diff)
总结：
这段代码通过减去指数加权移动平均值来进一步去除对数转换后的时间序列中的趋势。差分后的结果被绘制出，并检查了其平稳性。差分后的序列应该更加接近平稳性。
现在时间序列已经平稳，DS（检验统计量）值小于1%的临界值。由于权重是从起始索引分配的，因此检验结果有效。
消除趋势和季节性的其他方法
这里我将实现差分和分解方法来消除趋势和季节性。
# 绘制对数转换后的时间序列
ts_logtransformed.plot()
总结：
这段代码绘制了对数转换后的时间序列，用以观察时间序列的变化趋势和季节性成分。
通过差分消除季节性
时间序列中的季节性成分可以通过差分来消除。如果存在一个月的季节性成分，那么可以通过从当前值中减去上一个月的值来消除。例如，10月1日的值减去9月1日的值，10月2日的值减去9月2日的值，以此类推。差分后的序列将不包含第一个月的数据，因为第一个月的数据无法进行差分。
# 尝试不同的季节性差分，并测试平稳性数据
ts_diff_logtrans = ts_logtransformed - ts_logtransformed.shift(7) # 通过7天差分来消除季节性
plt.plot(ts_diff_logtrans)
ts_diff_logtrans.head(10)
ts_diff_logtrans.dropna(inplace=True) # 删除缺失值
test_for_stationary(ts_diff_logtrans)
总结：
这段代码通过对数转换后的时间序列进行了7天的差分，以消除季节性成分。差分后的序列被绘制出，并检查了其平稳性。差分后的序列应该更加接近平稳性。
通过分解消除趋势和季节性
decomposition = seasonal_decompose(ts_logtransformed) # 对时间序列进行季节性分解
trend = decomposition.trend # 提取趋势成分
seasonal = decomposition.seasonal # 提取季节性成分
residual = decomposition.resid # 提取残差成分
plt.subplot(411)
plt.plot(ts_logtransformed, label='原序列')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='趋势')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal, label='季节性')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='残差')
plt.legend(loc='best')
plt.tight_layout() # 调整子图布局
plt.show()
总结：
这段代码通过seasonal_decompose函数对对数转换后的时间序列进行了季节性分解，提取了趋势成分、季节性成分和残差成分。分解后的成分分别被绘制在四个子图中，以便于观察时间序列的不同成分。最后，对残差成分进行了平稳性检验。
decomposed_TS = residual # 将残差成分作为分解后的时间序列
decomposed_TS.dropna(inplace=True) # 删除缺失值
test_for_stationary(decomposed_TS) # 检查残差成分的平稳性
预测
# ACF和PACF图
lag_acf = acf(ts_diff_logtrans, nlags=30) # 计算自相关函数（ACF），滞后30期
lag_pacf = pacf(ts_diff_logtrans, nlags=50, method='ols') # 计算偏自相关函数（PACF），滞后50期，使用普通最小二乘法
# 绘制ACF图
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0, linestyle='--', color='gray') # 绘制0线
plt.axhline(y=-1.96/np.sqrt(len(ts_diff_logtrans)), linestyle='--', color='gray') # 绘制下置信区间线
plt.axhline(y=1.96/np.sqrt(len(ts_diff_logtrans)), linestyle='--', color='gray') # 绘制上置信区间线
plt.title('自相关函数')
# 绘制PACF图
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0, linestyle='--', color='gray') # 绘制0线
plt.axhline(y=-1.96/np.sqrt(len(ts_diff_logtrans)), linestyle='--', color='gray') # 绘制下置信区间线
plt.axhline(y=1.96/np.sqrt(len(ts_diff_logtrans)), linestyle='--', color='gray') # 绘制上置信区间线
plt.title('偏自相关函数')
plt.tight_layout() # 调整子图布局
plt.show()
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from matplotlib import pyplot
pyplot.figure()
pyplot.subplot(211)
plot_acf(ts_diff_logtrans, ax=pyplot.gca(),lags=40)
pyplot.subplot(212)
plot_pacf(ts_diff_logtrans, ax=pyplot.gca(), lags=50)
pyplot.show()
总结：
这段代码计算并绘制了时间序列ts_diff_logtrans的自相关函数（ACF）和偏自相关函数（PACF）图。ACF图用于确定移动平均（MA）参数q，而PACF图用于确定自回归（AR）参数p。图表中的虚线表示置信区间，可以用来判断哪些滞后值显著。
从图表中的一些观察结果
在这个图表中，两侧的虚线表示置信区间。这些虚线可以用来确定自回归（AR）和移动平均（MA）参数p和q：
p值是PACF图中首次超过上置信区间的滞后值。在这个例子中，p=2，PACF在第2天显示了一个显著的滞后值。q值是ACF图中首次超过上置信区间的滞后值。在这个例子中，q=18，ACF在第18天显示了一个显著的滞后值。
PACF在第7天、第8天、第12天、第15天等处可能也有一些显著的滞后值，这表明差分后的数据中仍然存在一些季节性成分。我们可以通过尝试不同的ARIMA模型来识别最佳模型，通过计算每个模型的残差平方和（RSS）来比较模型的性能。RSS值越低，模型的拟合效果越好。
# 从 statsmodels.tsa.arima.model 导入 ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX as ARIMA
ts_diff_logtrans = ts_diff_logtrans.fillna(0) # 用0填充缺失值，避免将来出现错误
# 构建和拟合ARIMA模型 (p=8, d=1, q=0)
model = ARIMA(ts_logtransformed, order=(8, 1, 0))
results_AR = model.fit(disp=-1) # 拟合模型，disp=-1表示不显示迭代过程
plt.plot(ts_diff_logtrans) # 绘制差分后的序列
plt.plot(results_AR.fittedvalues, color='red', label='模型拟合值 (p=8)') # 绘制模型拟合值
RSS = results_AR.fittedvalues - ts_diff_logtrans # 计算残差
RSS.dropna(inplace=True) # 删除缺失值
plt.title('RSS: %.4f' % sum(RSS**2)) # 设置图表标题，显示RSS值
plt.legend(loc='best') # 设置图例位置
plt.show() # 显示图表
# 构建和拟合ARIMA模型 (p=2, d=1, q=0)
model = ARIMA(ts_logtransformed, order=(2, 1, 0))
results_AR = model.fit(disp=-1) # 拟合模型，disp=-1表示不显示迭代过程
plt.plot(ts_diff_logtrans) # 绘制差分后的序列
plt.plot(results_AR.fittedvalues, color='red', label='模型拟合值 (p=2)') # 绘制模型拟合值
RSS = results_AR.fittedvalues - ts_diff_logtrans # 计算残差
RSS.dropna(inplace=True) # 删除缺失值
plt.title('RSS: %.4f' % sum(RSS**2)) # 设置图表标题，显示RSS值
plt.legend(loc='best') # 设置图例位置
plt.show() # 显示图表
# 打印模型摘要
print(results_AR.summary())
总结：
这段代码构建并拟合了两个ARIMA模型，一个是p=8, d=1, q=0，另一个是p=2, d=1, q=0。每个模型的拟合值和残差都被绘制出来，并计算了残差平方和（RSS）。RSS值越低，模型的拟合效果越好。最后，打印了模型的摘要信息，包括模型参数和统计检验结果。
# 构建和拟合ARIMA模型 (p=0, d=1, q=18)
model = ARIMA(ts_logtransformed, order=(0, 1, 18))
results_MA = model.fit(disp=-1) # 拟合模型，disp=-1表示不显示迭代过程
plt.plot(ts_diff_logtrans) # 绘制差分后的序列
plt.plot(results_MA.fittedvalues, color='red', label='模型拟合值 (q=18)') # 绘制模型拟合值
RSS = results_MA.fittedvalues - ts_diff_logtrans # 计算残差
RSS.dropna(inplace=True) # 删除缺失值
plt.title('RSS: %.4f' % sum(RSS**2)) # 设置图表标题，显示RSS值
plt.show() # 显示图表
# 打印模型摘要
print(results_MA.summary())
总结：
这段代码构建并拟合了一个ARIMA模型，参数为p=0, d=1, q=18。模型的拟合值和残差被绘制出来，并计算了残差平方和（RSS）。RSS值越低，模型的拟合效果越好。最后，打印了模型的摘要信息，包括模型参数和统计检验结果。
# 绘制对数转换后的时间序列和模型残差
plt.plot(ts_logtransformed, label='log_tranfromed_data')
plt.plot(results_MA.resid, color='green', label='Residuals')
plt.title('MA Model Residual plot')
plt.legend(loc='best')
plt.show()
# 绘制残差的密度图
results_MA.resid.plot(kind='kde')
plt.title('Density plot of the residual error values')
plt.show()
# 打印残差的描述统计信息
print(results_MA.resid.describe())
总结：
这段代码绘制了对数转换后的时间序列和MA模型的残差，用以观察模型的拟合效果。残差的密度图被绘制出来，以检查残差的分布情况。最后，打印了残差的描述统计信息，包括均值、标准差、最小值、最大值等。这些统计信息可以帮助我们了解残差的分布特征。
ARIMA 组合模型
介绍部分
ARIMA 组合模型
model = ARIMA(ts_logtransformed, order=(8, 1, 18))
results_ARIMA = model.fit(trend= 'nc', disp=-1)
plt.plot(ts_diff_logtrans)
plt.plot(results_ARIMA.fittedvalues, color='red', label = 'p =8, q =18')
RSS =results_ARIMA.fittedvalues-ts_diff_logtrans
RSS.dropna(inplace=True)
plt.title('RSS: %.4f'% sum(RSS**2))
plt.legend(loc='best')
上面的图表和模型展示了不同的残差 RSS（残差平方和）。评估所有这些模型可能会很困难，因此我们对所有可能的 ARIMA 参数进行网格搜索。我们寻找不同的 p、d、q 组合，并找到最佳组合。p 的取值为 7、10、13、16、19，d 的取值为 0 到 2，模型运行 10 次，这需要一些时间来完成。
以下是一个 Python 代码，用于评估不同模型的性能，以找到最佳的 ARIMA 超参数。
代码部分
代码总结：
这段代码的目的是通过网格搜索来评估不同 ARIMA 模型的性能，找到最佳的 p 和 d 参数组合。代码的逻辑是通过遍历不同的 p 和 d 值，计算每个组合的 RSS（残差平方和），并选择 RSS 最小的组合作为最佳模型。
import warnings
# 定义函数来评估 ARIMA 模型的性能
def evaluate_arima_model(data_set, arima_order):
model = ARIMA(data_set, order=arima_order) # 创建 ARIMA 模型
results_ARIMA = model.fit(disp=-1) # 拟合模型
RSS_diff = results_ARIMA.fittedvalues - ts_diff_logtrans # 计算残差
RSS = RSS_diff**2 # 计算残差平方和
return RSS
# 定义函数来评估不同 p 和 d 值的组合
def evaluate_models(dataset, p_values, d_values):
best_score, best_cfg = float("inf"), None # 初始化最佳得分和最佳配置
for p in p_values: # 遍历 p 值
for d in d_values: # 遍历 d 值
order = (p, d, 18) # 设置 ARIMA 模型的阶数
try:
rss = evaluate_arima_model(dataset, order) # 评估模型
if rss < best_score: # 如果当前 RSS 小于最佳 RSS
best_score, best_cfg = rss, order # 更新最佳得分和最佳配置
print('ARIMA%s RSS=%.3f' % (order, rss)) # 打印当前模型的 RSS
except:
continue
print('Best ARIMA%s RSS=%.3f' % (best_cfg, best_score)) # 打印最佳模型的 RSS
p_values = range(8, 20, 3) # 设置 p 值的范围
d_values = range(0, 3) # 设置 d 值的范围
warnings.filterwarnings('ignore') # 忽略警告
evaluate_models(ts_logtransformed, p_values, d_values) # 评估模型
代码总结：
这段代码的目的是生成 ARIMA 模型的摘要，并绘制残差图。代码的逻辑是通过拟合 ARIMA 模型，生成模型的摘要信息，并绘制残差图以评估模型的性能。
# 打印 ARIMA 模型的摘要
print(results_ARIMA.summary())
# 绘制对数变换后的数据和残差图
plt.plot(ts_logtransformed, label='log_transformed_data')
plt.plot(results_ARIMA.resid, color='green', label='Residuals')
plt.title('ARIMA Model Residual plot')
plt.legend(loc='best')
# 绘制残差误差值的密度图
results_ARIMA.resid.plot(kind='kde')
plt.title('Density plot of the residual error values')
print(results_ARIMA.resid.describe())
test = pd.read_csv("../input/bitcoin_price_1week_Test - Test.csv",index_col= 'Date')
test.index = pd.to_datetime(test.index)
test = test['Close']
test = test.sort_index()
test
代码总结：
这段代码的目的是使用 ARIMA 模型进行预测和预测结果的累积求和。代码的逻辑是通过拟合 ARIMA 模型，生成预测值，并对预测值进行累积求和，最终将结果转换回原始尺度。
# 使用 ARIMA 模型进行预测
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print(predictions_ARIMA_diff.head())
# 对预测值进行累积求和
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
print(predictions_ARIMA_diff_cumsum.head())
# 将累积求和的结果转换回对数尺度
predictions_ARIMA_log = pd.Series(ts_logtransformed.iloc[0], index=ts_logtransformed.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum, fill_value=0)
predictions_ARIMA_log.head()
# 将对数尺度的预测值转换回原始尺度
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(data)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f' % np.sqrt(sum((predictions_ARIMA - data)**2) / len(data)))
dates = [pd.Timestamp('2017-08-01'), pd.Timestamp('2017-08-02'), pd.Timestamp('2017-08-03'),pd.Timestamp('2017-08-04'), pd.Timestamp('2017-08-05'), pd.Timestamp('2017-08-06'), pd.Timestamp('2017-08-07')]
forecast = pd.Series(results_ARIMA.forecast(steps=7)[0],dates)
forecast = np.exp(forecast)
print(forecast)
error = mean_squared_error(test, forecast)
print('Test MSE: %.3f' % error)
plt.plot(forecast, color ='green', label ='Predicted rates')
plt.plot(test, color = 'red', label = 'Observed from test data')
plt.title('RMSE: %.4f'% np.sqrt(sum((forecast-test)**2)/len(data)))
plt.legend(loc = 'best')
代码总结：
这段代码的目的是使用 MA 模型进行预测，并计算预测值与实际值之间的均方误差（MSE）。代码的逻辑是通过拟合 MA 模型，生成预测值，并计算预测值与实际值之间的误差。
# 使用 MA 模型进行预测
predictions_MA_diff = pd.Series(results_MA.fittedvalues, copy=True)
print(predictions_MA_diff.head())
# 对预测值进行累积求和
predictions_MA_diff_cumsum = predictions_MA_diff.cumsum()
print(predictions_MA_diff_cumsum.head())
# 将累积求和的结果转换回对数尺度
predictions_MA_log = pd.Series(ts_logtransformed.iloc[0], index=ts_logtransformed.index)
predictions_MA_log = predictions_MA_log.add(predictions_MA_diff_cumsum, fill_value=0)
predictions_MA_log.head()
# 将对数尺度的预测值转换回原始尺度
predictions_MA = np.exp(predictions_MA_log)
plt.plot(data)
plt.plot(predictions_MA)
plt.title('RMSE: %.4f' % np.sqrt(sum((predictions_MA - data)**2) / len(data)))
# 计算预测值与实际值之间的均方误差
error = mean_squared_error(test, forecast)
print('Test MSE: %.3f' % error)
以下是翻译后的文章内容，包括对代码部分的注释和总结：
# 定义日期列表
dates = [pd.Timestamp('2017-08-01'), pd.Timestamp('2017-08-02'), pd.Timestamp('2017-08-03'),
pd.Timestamp('2017-08-04'), pd.Timestamp('2017-08-05'), pd.Timestamp('2017-08-06'), pd.Timestamp('2017-08-07')]
# 使用移动平均模型对未来的7天进行预测，并将预测结果存储在 forecast 变量中
forecast = pd.Series(results_MA.forecast(steps=7)[0], dates)
# 对预测结果进行指数变换，恢复到原始数据的尺度
forecast = np.exp(forecast)
# 打印预测结果
print(forecast)
# 计算测试数据与预测数据之间的均方误差（MSE）
error = mean_squared_error(test, forecast)
# 打印测试均方误差（MSE）
print('Test MSE: %.3f' % error)
# 以下代码的目的是绘制预测结果和测试数据的图表
plt.plot(forecast, color='green', label='Predicted rates') # 绘制预测结果，颜色为绿色，标签为 Predicted rates
plt.plot(test, color='red', label='Observed from test data') # 绘制测试数据，颜色为红色，标签为 Observed from test data
# 设置图表标题，并显示均方根误差（RMSE）
plt.title('RMSE: %.4f' % np.sqrt(sum((forecast - test) ** 2) / len(data)))
plt.legend(loc='best') # 设置图例位置为最佳
消除上述时间序列中的季节性是一个需要进行大量工作的任务。可以扩展使用更多方法，如曲线拟合和差分，来识别季节性并从数据中去除季节性成分。
改进的模型可以是从前一年的相同日历月的平均币价中减去，而不是从相同日期中减去。
我们可以通过将数据集重新采样为月平均价格来开始。重新采样方法处理了闰年的概念，同时消除了偏移（例如2月只有28天）。
# 将数据重新采样为月平均价格
monthly_mean = data.resample('M').mean()
# 打印前13个月的月平均价格
print(monthly_mean.head(13))
# 绘制月平均价格的图表
monthly_mean.plot()
以下是另一种方法来评估预测性能：
# 对测试数据进行对数变换
test_logtransformed = np.log(test)
# 初始化历史数据列表，包含对数变换后的训练数据
history = [x for x in ts_logtransformed]
# 初始化预测结果列表
predictions = list()
# 遍历测试数据的每个时间点
for t in range(len(test)):
# 使用移动平均模型进行预测
output = results_MA.forecast()
yhat = output[0] # 获取预测结果
predictions.append(yhat) # 将预测结果添加到预测列表中
# 获取实际观察值并添加到历史数据列表中
obs = test_logtransformed[t]
history.append(obs)
print('predicted=%f, expected=%f' % (yhat, obs)) # 打印预测值和实际值
# 计算对数变换后的测试数据与预测数据之间的均方误差（MSE）
error = mean_squared_error(test_logtransformed, predictions)
# 打印测试均方误差（MSE）
print('Test MSE: %.3f' % error)
消除上述时间序列中的季节性是一个需要进行大量工作的任务。可以扩展使用更多方法，如曲线拟合和差分，来识别季节性并从数据中去除季节性成分。
改进的模型可以是从前一年的相同日历月的平均币价中减去，而不是从相同日期中减去。
我们可以通过将数据集重新采样为月平均价格来开始。重新采样方法处理了闰年的概念，同时消除了偏移（例如2月只有28天）。
代码总结：
这段代码的目的是使用 Prophet 模型进行时间序列预测。代码的逻辑是通过拟合 Prophet 模型，生成未来 7 天的预测值，并绘制预测结果。
from prophet import Prophet
# 准备数据
data_prophet = data.copy()
data_prophet = pd.DataFrame(data_prophet)
data_prophet.reset_index(drop=False, inplace=True)
data_prophet.columns = ['ds', 'y']
# 创建并拟合 Prophet 模型
m = Prophet()
m.fit(data_prophet)
# 生成未来 7 天的预测数据
future = m.make_future_dataframe(periods=7, freq='D')
forecast = m.predict(future)
# 绘制预测结果
m.plot(forecast)
data.plot()
# 绘制预测结果的组件
m.plot_components(forecast)
# 提取预测值并计算均方误差
forecasted_values = forecast[['ds', 'yhat']].tail(7)
forecasted_values = forecasted_values.set_index('ds')
forecasted_values.columns = ['y']
mean_squared_error(forecasted_values['y'], test)
代码总结：
这段代码的目的是使用 LSTM 模型进行时间序列预测。代码的逻辑是通过创建和训练 LSTM 模型，对时间序列数据进行预测。
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from sklearn.preprocessing import MinMaxScaler
# 创建 LSTM 模型
model = Sequential()
model.add(LSTM(50, input_shape=(1, 1)))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
# 数据预处理
scaler = MinMaxScaler(feature_range=(0, 1))
scaled_data = scaler.fit_transform(data)
# 训练模型
model.fit(scaled_data, epochs=100, batch_size=1, verbose=2)
总结
本文介绍了如何使用 ARIMA 和 MA 模型进行时间序列预测，并通过网格搜索找到最佳的超参数组合。此外，还介绍了如何使用 Prophet 和 LSTM 模型进行时间序列预测。通过这些方法，可以有效地预测时间序列数据，并评估模型的性能。
原文地址：https://blog.csdn.net/jrckkyy/article/details/145235553
免责声明：本站文章内容转载自网络资源，如侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！