时间序列模型(ARIMA) — — 预测未来(含 python 和 Matlab 的完整代码)

之前的两篇笔记用 Matlab 和 python 实现了对股票的预测，但是实现效果并不理想，近两天又翻阅了其他的一些资料，对预测的代码进行了改进。

望月12138

3262人浏览 · 2024-08-29 21:08:52

望月12138 · 2024-08-29 21:08:52 发布

文章目录

前言
一、python 代码实现
二、Matlab 代码实现
三、python 和 Matlab 的完整代码
- python
- Matlab

前言

【Matlab】时间序列模型(ARIMA)
【python】时间序列模型(ARIMA)

之前的两篇笔记用 Matlab 和 python 实现了对股票的预测，但是实现效果并不理想，近两天又翻阅了其他的一些资料，对预测的代码进行了改进。

一、python 代码实现

python 的 statsmodels 库中预测函数大致分为两个，一个是 predict 函数，另一个是 forecast 函数。两个函数的使用上有所不同。

`statsmodels` 中的 `predict` 和 `forecast`

`predict`

用途: predict 函数通常用于生成模型拟合期间的数据预测。它可以用于生成拟合数据的未来预测，也可以用于生成训练数据期间的拟合值。
用法: predict(start, end)，其中 start 和 end 参数定义了预测的时间范围。

`forecast`

用途: forecast 函数专门用于生成未来的预测值。它通常在模型已经拟合好的情况下使用，用于生成未来多个时间点的预测。
用法: forecast(steps)，其中 steps 参数定义了要预测的未来时间点数量。

选择使用哪个函数取决于你的具体需求。如果你需要预测未来的时间点数据，使用 forecast 更合适；如果你需要对现有数据或拟合期间的数据进行预测，使用 predict 更合适。

下面的链接有更加详细的解释，供读者参考：

https://www.statsmodels.org/dev/examples/notebooks/generated/statespace_forecasting.html

实际运用

因为我们需要对未来的时间点进行预测，所以选用 forecast 函数，但是如果直接调用 forecast_values = fit.forecast(steps=3) 实现预测未来3个时间点的值，会发现预测出的三个值是一样的，显然结果不对。

后来是参考了国外一个机器学习的教学网站：

https://machinelearningmastery.com/make-manual-predictions-arima-models-python/

了解到，为了提高预测的准确度，可以使用滚动预测。

滚动预测

滚动预测方法的使用方式如下：

测试数据集中的每个时间点都会迭代。
在每次迭代中，都会根据所有可用的历史数据训练新的 ARIMA 模型。
该模型用于预测第二天的情况。
将存储预测，并从测试集中检索“真实”观察结果，并将其添加到历史记录中以供下次迭代使用。

简单来说，就是先通过训练集拟合出来的模型往后预测一天，然后将测试集中这天的真实值加入到训练集里面，再进行拟合预测，反复迭代。直到全部预测出训练集中的数据。

向后预测的代码

p = 1
d = 1
q = 0

history = list(train.values)
forecast = list()
for t in range(len(test.values)):
	model = sm.tsa.ARIMA(history, order=(p,d,q))
	model_fit = model.fit()
	output = model_fit.forecast()
	yhat = output[0]
	forecast.append(yhat)
	obs = test[t]
	history.append(obs) # 这里实现了滚动预测，应该是ARIMA只能预测下一步而不能多步

predict = pd.Series(forecast, index=test.index)
# 绘制结果
plt.figure(figsize=(12, 6))
plt.plot(train.index, train, label='Train')
plt.plot(test.index, test, label='Test')
plt.plot(predict.index, predict, label='Predict', linestyle='--')
# plt.plot(predict_sunspots, label='Fit')
plt.legend()
plt.xticks(rotation=45)
plt.show()

运行结果：

在这里插入图片描述

可以发现，预测的结果比之前好了很多。

二、Matlab 代码实现

在之前的博客中有存疑的地方，现已解决。

ADF检验

首先是 ADF 检验部分，在 python 中 ADF 检验的结果为：

timeseries_adf :  (np.float64(-2.902677813259885), np.float64(0.04503748919120268), 0, 61, {'1%': np.float64(-3.542412746661615), '5%': np.float64(-2.910236235808284), '10%': np.float64(-2.5927445767266866)}, np.float64(-227.51423828600673))

可以看到，第一个数据只是小于了 '10%'后的数据，并没有很好的拒绝原假设，所以可以理解为该数据不平稳。

在 Matlab 中，检验的结果为：

在这里插入图片描述

同样没有拒绝原假设。

所以后来对代码进行改进时，我也在 python 代码里对数据进行了一阶差分，确保了数据的平稳性。

确定 p、q 的值

一开始使用 AIC 准则和 BIC 准则获取的结果如下：

在这里插入图片描述

显然是不太对劲的，怀疑是代码部分的问题。

之后我又研究了一下 Matlab 中 aicbic(logL, numParam, n) 函数。

`aicbic(logL, numParam, n)` 函数

在 ARIMA 模型中，numParams 代表模型的参数数量。具体来说，ARIMA 模型由三个主要部分组成：AR（AutoRegressive，自回归部分）、I（Integrated，差分部分）和 MA（Moving Average，移动平均部分）。模型的形式通常表示为 ARIMA(p, d, q)，其中：

p 是自回归项的阶数。
d 是差分次数。
q 是移动平均项的阶数。

详细解释

定义 aicbic 函数：
- 该函数接收对数似然值、参数数量和样本数量来计算 AIC 和 BIC。
示例数据和参数：
- LogL 是对数似然值。
- p, d, q 是 ARIMA 模型的参数。
- include_constant 是一个布尔值，指示是否包括常数项。
- n 是样本数量。
计算参数数量：
- 使用公式 numParams = p + q + (1 if include_constant else 0)。
计算 AIC 和 BIC 并打印结果：
- 使用 aicbic 函数计算并打印 AIC 和 BIC 值。

所以使用网格搜索法确定 p、q 的值，调用aicbic(logL, numParam, n) 函数时，numParams = p + q

修改代码之后的结果：

在这里插入图片描述

BIC 准则的结果相对更加准确些，模型复杂性也较低，所以确定出 p = 1，q = 0。

向后预测的代码

在 Matlab 中我也同样使用了滚动预测的思想，代码如下：

history = train;
predict = zeros(length(test),1);

for i = 1:length(test)
    % 简化模型为 ARIMA(1,1,0)
    model = arima(1,1,0);
    % 估计模型参数
    md1 = estimate(model, history, 'Display', 'off');
    % 进行预测
    [Y, YMSE] = forecast(md1, 1, 'Y0', history);
    % 存储预测结果
    predict(i) = Y;
    % 更新历史数据
    history = [history; test(i)];
end
origin_close = close_data(1:127);
origin_date = date_data(1:127);
% 绘制预测结果与真实值的比较
figure('Position', [100, 100, 1200, 700]); 
plot(origin_date,origin_close, test_date, predict, train_date, train);
legend('真实值','预测值','训练值');
title('ARIMA 模型预测结果');
xlabel('时间');
ylabel('值');

运行结果：

在这里插入图片描述
可以看出，预测结果还不错。

三、python 和 Matlab 的完整代码

python

# %%
import pandas as pd
import matplotlib.pyplot as plt

# 导入数据

ChinaBank = pd.read_csv('ChinaBank.csv',index_col = 'Date',parse_dates=['Date'])

ChinaBank.head()

# %%
# 提取Close列

ChinaBank.index = pd.to_datetime(ChinaBank.index)
sub = ChinaBank.loc['2014-01':'2014-06','Close']

sub.head()

# %%
# 划分训练测试集

train = sub.loc['2014-01':'2014-04']
test = sub.loc['2014-04':'2014-06']

#查看训练集的时间序列与数据(只包含训练集)
plt.figure(figsize=(12,6))
plt.plot(train)
plt.xticks(rotation=45) #旋转45度
plt.show()

# %%
from statsmodels.tsa.stattools import adfuller as ADF

timeseries_adf = ADF(train.tolist())

# 打印单位根检验结果
print('timeseries_adf : ', timeseries_adf)

# %%
# 一阶差分
train_dif1 =  train.diff(1).fillna(0)

# %%
# 打印单位根检验结果
print('一阶差分 : ',ADF(train_dif1.tolist()))

# %%
train_dif1.plot(subplots=True,figsize=(12,6))
plt.show()

# %%
# 参数确定
import statsmodels.api as sm

# 绘制
fig = plt.figure(figsize=(12,7))

ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(train_dif1, lags=20,ax=ax1)
ax1.xaxis.set_ticks_position('bottom') # 设置坐标轴上的数字显示的位置，top:显示在顶部  bottom:显示在底部
#fig.tight_layout()

ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(train_dif1, lags=20, ax=ax2)
ax2.xaxis.set_ticks_position('bottom')
#fig.tight_layout()
plt.show()

# %%
train_results = sm.tsa.arma_order_select_ic(train, ic=['aic', 'bic'], trend='n', max_ar=6, max_ma=6)

print('AIC', train_results.aic_min_order)
print('BIC', train_results.bic_min_order)

# %%
# 模型检验
#根据以上求得
p = 1
d = 1
q = 0

model = sm.tsa.ARIMA(train, order=(p,d,q))
results = model.fit()
resid = results.resid #获取残差

#绘制
#查看测试集的时间序列与数据(只包含测试集)
fig, ax = plt.subplots(figsize=(12, 5))

ax = sm.graphics.tsa.plot_acf(resid, lags=40, ax=ax)

plt.show()

# %%
start_index = train.index[2]
end_index = train.index[-1]
predict_sunspots = results.predict(start_index, end=end_index)
print(predict_sunspots)

#查看训练集的时间序列与数据(只包含训练集)
plt.figure(figsize=(12,6))
plt.plot(train)
plt.xticks(rotation=45) #旋转45度
plt.plot(predict_sunspots)
plt.show()

# %%
p = 1
d = 1
q = 0

history = list(train.values)
forecast = list()
for t in range(len(test.values)):
	model = sm.tsa.ARIMA(history, order=(p,d,q))
	model_fit = model.fit()
	output = model_fit.forecast()
	yhat = output[0]
	forecast.append(yhat)
	obs = test[t]
	history.append(obs) # 这里实现了滚动预测，应该是ARIMA只能预测下一步而不能多步

predict = pd.Series(forecast, index=test.index)
# 绘制结果
plt.figure(figsize=(12, 6))
plt.plot(train.index, train, label='Train')
plt.plot(test.index, test, label='Test')
plt.plot(predict.index, predict, label='Predict', linestyle='--')
# plt.plot(predict_sunspots, label='Fit')
plt.legend()
plt.xticks(rotation=45)
plt.show()

Matlab

clc;clear
%% 数据读取
% 读取 CSV 文件
filename = 'ChinaBank.csv';
data = readtable(filename);

% 读取文件中的两列
close_data = data.Close;
date_data = data.Date;

% 划分训练测试集
train = close_data(1:62);
test = close_data(63:127);
train_date = date_data(1:62);
test_date = date_data(63:127);

% 定义候选模型阶数范围
maxP = 8;
maxQ = 8;
n = length(train);

% 初始化结果存储
aicValues = NaN(maxP, maxQ);
bicValues = NaN(maxP, maxQ);

% 迭代计算所有候选模型的AIC和BIC值
for p = 0:maxP
    for q = 0:maxQ
        try
            Mdl = arima(p,1,q);
            [~,~,logL] = estimate(Mdl, train, 'Display', 'off');
            numParam = p + q; % p个AR参数, q个MA参数
            [aicValues(p+1, q+1),bicValues(p+1, q+1)] = aicbic(logL, numParam, n);
        catch
            % 忽略无法估计的模型
            continue;
        end
    end
end

% 找到AIC最小值对应的(p, q)
[minAIC, idxAIC] = min(aicValues(:));
[pAIC, qAIC] = ind2sub(size(aicValues), idxAIC);
pAIC = pAIC - 1;
qAIC = qAIC - 1;

% 找到BIC最小值对应的(p, q)
[minBIC, idxBIC] = min(bicValues(:));
[pBIC, qBIC] = ind2sub(size(bicValues), idxBIC);
pBIC = pBIC - 1;
qBIC = qBIC - 1;

fprintf('AIC选择的模型阶数: p = %d, q = %d\n', pAIC, qAIC);
fprintf('BIC选择的模型阶数: p = %d, q = %d\n', pBIC, qBIC);

model = arima(1,1,0);
md1 = estimate(model, train, 'Display', 'off');

% 检查残差的自相关性
residuals = infer(md1, train);
figure;
autocorr(residuals);
title('Residuals Autocorrelation');

history = train;
predict = zeros(length(test),1);

for i = 1:length(test)
    % 简化模型为 ARIMA(1,1,0)
    model = arima(1,1,0);
    % 估计模型参数
    md1 = estimate(model, history, 'Display', 'off');
    % 进行预测
    [Y, YMSE] = forecast(md1, 1, 'Y0', history);
    % 存储预测结果
    predict(i) = Y;
    % 更新历史数据
    history = [history; test(i)];
end
origin_close = close_data(1:127);
origin_date = date_data(1:127);
% 绘制预测结果与真实值的比较
figure('Position', [100, 100, 1200, 700]); 
plot(origin_date,origin_close, test_date, predict, train_date, train);
legend('真实值','预测值','训练值');
title('ARIMA 模型预测结果');
xlabel('时间');
ylabel('值');