SARIMAX时间序列分析预测代码

这个只是传统的时间序列分析代码，数据集暂时没有，代码比较详细，可以直接用自己的数据进行替换

#这些都是所需要的加载包，先运行一遍，会进行报错，哪个包没有它会显示，如果没有相关的module，如缺少pandas，就找一个空白的行，输入pip install pandas
#安装好之后，重新跑这个区域的代码
import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import pmdarima as pm
from pmdarima.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
import pdb
import statsmodels.api as sm
warnings.filterwarnings('ignore')

#这一步是文件的读取
df=pd.read_excel(r'C:\Users\12464\Desktop\PVW\DATA\PVW_time\summer_30data.xlsx')#r''里面的东西是文件的所在位置，index_col='time stamp'是索引
print(df.shape)#显示数据的形状
df.head()#显示数据的前五行

PV_face_T=df['PV face T']#设定因变量及名字，[‘’]引号中的内容是原始数据里的因变量的名城
exog=df.iloc[:,[1,2]]#设定外生变量，1，2代表原始数据中的第2列、3列
print(PV_face_T)#显示PV_Back_T
print(exog)#显示exog

#绘制原始数据图、自相关图、偏自相关图
fig = plt.figure(figsize=(16,10),dpi=200)#画板大小尺寸设置
def plot(PV_face_T):#这些都不用动，只要把你需要画的因变量的名字放到plot（）里面，这个模块里的统一替换
    results = adfuller(PV_face_T)
    results_str = 'ADF test, p-value is: {}'.format(results[1])

    grid = plt.GridSpec(2, 2)
    ax1 = plt.subplot(grid[0, :])
    ax2 = plt.subplot(grid[1, 0])
    ax3 = plt.subplot(grid[1, 1])

    ax1.plot(PV_face_T)
    ax1.set_title(results_str)
    plot_acf(PV_face_T, lags=90, ax=ax2)
    plot_pacf(PV_face_T, lags=90, ax=ax3)
    plt.show()
    
plot(PV_face_T)

#数据的自相关和偏自相关绘图，和上一步重了，这里面看自相关图比较清楚
fig = plt.figure(figsize=(16,10),dpi=200)
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(PV_face_T, lags=40, ax=ax1)#自相关图
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(PV_face_T, lags=40, ax=ax2)#偏自相关图
plt.show()
#根据相关图确定周期性，或者自己的数据周期性可以在excel看出来

#绘制进行周期后的数据图像
#看自相关图可以看到数据存在周期性，周期为29，因此这里要对数据进行季节性差分，这里diff（29）代表进行周期为29的差分
#季节性差分的数字根据你自己因变量的周期来看
fig = plt.figure(figsize=(16,10),dpi=200)
PV_face_T_seasonal_diff =PV_face_T.diff(29)
plt.plot(PV_face_T_seasonal_diff)
plt.show()

#绘制周期差分后的自相关与偏自相关图
fig = plt.figure(figsize=(16,10))
ax1 = fig.add_subplot(211)
#Plots with first order differenced Sales variable
fig = sm.graphics.tsa.plot_acf(PV_face_T_seasonal_diff.dropna(), lags=70, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(PV_face_T_seasonal_diff.dropna(), lags=70, ax=ax2)
plt.show()
#从自相关图可以看到，数据仍然不平稳，还需要进行一次差分，如果自相关图没有波动性，则说明数据平稳，到此为止，否则继续

#数据不平稳，继续进行一次差分，并做出进行差分后的图
fig = plt.figure(figsize=(16,10),dpi=200)
PV_face_T_seasonal_diff1 =PV_face_T_seasonal_diff.diff(1)
plt.plot(PV_face_T_seasonal_diff1)
plt.show()

#绘制经过差分后的自相关图与偏自相关图，可以看到，没有周期性了，数据平稳
fig = plt.figure(figsize=(16,10))
ax1 = fig.add_subplot(211)
#Plots with first order differenced Sales variable
fig = sm.graphics.tsa.plot_acf(PV_face_T_seasonal_diff1.dropna(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(PV_face_T_seasonal_diff1.dropna(), lags=40, ax=ax2)
plt.show()

#设置训练集与测试集，我们要预测最后一天的数据，因此选取最后一天的数据作为测试，其他的数据用于训练，让机器通过训练的数据来预测测试的值
#最后一天的数据有29个，开始设置训练集与测试集
PV_face_T_train=PV_face_T[0:203]#因变量的训练数集
PV_face_T_test=PV_face_T[203:]#因变量的测试数集
exog_train=exog[0:203]##外生变量的训练数集
exog_test=exog[203:]#外生变量的训练数集
print(exog_test)#显示外生变量的测试集
print(PV_face_T_test)#显示因变量的测试集
print(PV_face_T_train)#显示因变量的训练集

#开始建立模型SARIMA,不带外生变量的
mod1 = sm.tsa.statespace.SARIMAX(endog=PV_face_T_train,order=[1,1,1],seasonal_order=[1,1,2,29])#这一行出了括号里的东西，其他的都不用动
#endog=PV_face_T_train，endog后面放的是因变量的训练集，order=[1,1,1],这个是经过多次差分后最后一个自相关图和偏自相关图看的
#seasonal_order=[1，1，2，29],这里面除了29改掉，其他的不用管，144改成你自己数据的周期
results1 = mod1.fit()
print(results1.summary())

fitted=results1.fittedvalues#得到拟合值
print(fitted)
forecast1=results1.forecast(steps=29)#对未来29步的数据进行滚动预测
print(forecast1)

#绘制拟合值与预测值图
plt.figure(figsize=(30,18 ),dpi=200)#设立画布大小
PV_face_T.plot(label='PV_face_T',linewidth=5)#真实值
fitted.plot(label='fit',linewidth=5)#拟合值
forecast1.plot(label='forecast',linewidth=5)#预测值
plt.legend()
plt.grid(True)
plt.show()

#只绘制预测值
plt.figure(figsize=(30,18 ),dpi=200)
PV_face_T_test.plot(label='PV_face_T_test',linewidth=6)#用测试集的数据，字体的大小自己在csdn搜索怎么改
forecast1.plot(label='forecast',linewidth=6)
plt.legend()
plt.grid(True)
plt.show()

#开始建立模型SARIMAX,带外生变量的
mod2= sm.tsa.statespace.SARIMAX(endog=PV_face_T_train,exog=exog_train,order=[1,1,1],seasonal_order=[1,1,2,29])#这一行出了括号里的东西，其他的都不用动
#endog=PV_face_T_train，endog后面放的是因变量的训练集，exogb不用动，order=[1,1,1],这个是经过多次差分后最后一个自相关图和偏自相关图看的
#seasonal_order=[1，1，2，29],这里面除了29改掉，其他的不用管，144改成你自己数据的周期
results2 = mod2.fit()
print(results2.summary())

fitted2=results2.fittedvalues#得到拟合值
print(fitted2)
forecast2=results2.forecast(steps=29,exog=exog_test)#对未来29步的数据进行预测,这里要加exog=exog_test外生变量
print(forecast2)

#绘制拟合值与预测值图
plt.figure(figsize=(30,18 ),dpi=200)#设立画布大小
PV_face_T.plot(label='PV_face_T',linewidth=5)#真实值
fitted2.plot(label='fit2',linewidth=5)#拟合值
forecast2.plot(label='forecast2',linewidth=5)#预测值
plt.legend()
plt.grid(True)
plt.show()