SARIMA 模型 Python 代码 带数据库 sarima模型的预测

SARIMA时间序列模型预测城市房价数据

数据清洗

文件中含有大量城市的房价数据，考虑到此次为学习性质的练习，为了节省数据处理的繁琐步骤。我截取了北京的2010-2021房价数据作为样例，并将价格的数据格式改为数值，去除多余的逗号

数据导入

#导入数据
df = pd.read_csv('Beijingpricedata.csv', encoding='utf-8', index_col='Timing')
#将日期转化为标准时间格式
i=0;
for str in df.index:
    df.index.values[i]=datetime.datetime.strptime(str, '%y-%b').strftime('%Y-%m')
    i+=1;
df.index = pd.to_datetime(df.index)  # 将字符串索引转换成时间索引TravelDate

ts = df['北京'] # 生成pd.Series对象

数据平稳性检测

因为采用ARIMA模型训练数据，要求数据是稳定的，所以需要对数据进行平稳性检测

画出原始数据图表

可以看到价格数据具有逐步上升的趋势

#原始数据图
def draw_ts(timeseries):
    f = plt.figure(facecolor='white')
    timeseries.plot(color='blue')
    plt.title('housing price changes')
    plt.show()
draw_ts(ts)

移动平均图

我们可以发现数据的移动平均值有越来越大的趋势，是不稳定的。

def draw_trend(timeseries, size):
    plt.rcParams['font.sans-serif'] = [u'SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    f = plt.figure(facecolor='white')
    # 对size个数据进行移动平均
    rol_mean = timeseries.rolling(window=size).mean()
    # 对size个数据移动平均的方差
    rol_std = timeseries.rolling(window=size).std()

    timeseries.plot(color='blue', label='Original')
    rol_mean.plot(color='red', label='Rolling Mean')
    rol_std.plot(color='black', label='Rolling standard deviation')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show()
draw_trend(ts,12)

平稳性检测

此时p值为0.623007，不能拒绝原假设，即数据存在单位根,数据是非平稳序列。

def TestStationaryAdfuller(ts, cutoff = 0.01):
    ts_test = adfuller(ts, autolag = 'AIC')
    ts_test_output = pd.Series(ts_test[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])

    for key,value in ts_test[4].items():
        ts_test_output['Critical Value (%s)'%key] = value
    print(ts_test_output)

    if ts_test[1] <= cutoff:
        print(u"拒绝原假设，即数据没有单位根,序列是平稳的。")
    else:
        print(u"不能拒绝原假设，即数据存在单位根,数据是非平稳序列。")

# 平稳性检测:
def teststationarity(ts):
    dftest = adfuller(ts)
    # 对上述函数求得的值进行语义描述
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', '#Lags Used', 'Number of Observations Used'])
    for key, value in dftest[4].items():
        dfoutput['Critical Value (%s)' % key] = value
    return dfoutput
print(teststationarity(ts))

数据分解

将时序数据分离成不同的成分：趋势部分，季节性部分，残留部分。

可以看到，数据具有上升的趋势，且具有季节性的特点。周期=12（/月）

所以考虑采用季节ARIMA模型（SRIMA）方法模型拟合数据

def decompose(timeseries):
    # 返回包含三个部分 trend（趋势部分） ， seasonal（季节性部分） 和residual (残留部分)
    decomposition = seasonal_decompose(timeseries)

    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid

    plt.subplot(411)
    plt.plot(ts, label='Original')
    plt.legend(loc='best')
    plt.subplot(412)
    plt.plot(trend, label='Trend')
    plt.legend(loc='best')
    plt.subplot(413)
    plt.plot(seasonal, label='Seasonality')
    plt.legend(loc='best')
    plt.subplot(414)
    plt.plot(residual, label='Residuals')
    plt.legend(loc='best')
    plt.tight_layout()
    plt.show()
    return trend, seasonal, residual

trend, seasonal, residual = decompose(ts)
residual.dropna(inplace=True)
draw_trend(residual, 12)
print(teststationarity(residual))

SARIMA模型

模型定阶

用图解法求解ARIMA模型的最优参数并非易事，主观性很大，而且耗费时间。所以本文进一步考虑利用网格搜索的方法系统地选择最优的参数值

#SARIMA模型参数的选择
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
pdq_x_PDQs = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
a=[]
b=[]
c=[]
wf=pd.DataFrame()
for param in pdq:
    for seasonal_param in pdq_x_PDQs:
        try:
            mod = sm.tsa.statespace.SARIMAX(ts,order=param,seasonal_order=seasonal_param,enforce_stationarity=False,enforce_invertibility=False)
            results = mod.fit()
            print('ARIMA{}x{} - AIC:{}'.format(param, seasonal_param, results.aic))
            a.append(param)
            b.append(seasonal_param)
            c.append(results.aic)
        except:
            continue
wf['pdq']=a
wf['pdq_x_PDQs']=b
wf['aic']=c
print(wf[wf['aic']==wf['aic'].min()])


mod = sm.tsa.statespace.SARIMAX(ts,
                                order=(1,1,1),
                                seasonal_order=(1,1,1,12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)
results = mod.fit()
print(results.summary())

可以看到，模型参数的最佳组合为(1,1,1)(1,1,1,12),aic最小值为1475.5377

模型检验

# 模型检验
mod = sm.tsa.statespace.SARIMAX(ts,
                                order=(1,1,1),
                                seasonal_order=(1,1,1,12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)
results = mod.fit()
print(results.summary())
# 模型诊断
results.plot_diagnostics(figsize=(15, 12))
plt.show()
# LB检验
r, q, p = sm.tsa.acf(results.resid.values.squeeze(), qstat=True)
data = np.c_[range(1, 22), r[1:], q, p]
table = pd.DataFrame(data, columns=['lag', "AC", "Q", "Prob(>Q)"])
print(table.set_index('lag'))

coef列显示每个变量的权重(即重要性)，以及每个变量如何影响时间序列。P>|z|列是对每个变量系数的检验。每个变量的P值均小于0.05，所以在0.05的显著性水平下，拒绝加假设，模型中每个变量的系数通过显著性检验，可以认为拟合的模型中包含这些变量是合理的。

由右上角的正态分布图可知，模型的残差是正态分布的。红色的KDE线紧随着N(0，1)线。其中，N(0,1)是均值为0，标准差为1的标准正态分布。这很好地说明了残差是正态分布的。而图中左下角的正太Q—Q图也说明了残差服从正态分布。

图中右下角残差的自相关图显示残差在初始位置偏差较大，其他位置不存在自相关，说明残差序列是白噪声序列。

由此，可以得出结论，SARIMAX(1，1，1)x(0，1，1，12)模型的拟合效果还不错，可以帮助我们理解原始的时间序列数据并对未来的数值做出预测。

模型预测

1）静态预测

#静态预测
predict_ts = results.predict()
predict_ts.plot(color='blue', label='Predict')
ts.plot(color='red', label='Original')
plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((predict_ts-ts)**2)/ts.size))
plt.show()

可以看到，模型在2011左右预测效果不太理想，不过之后预测效果很不错

动态预测

#动态预测
pred_dynamic = results.get_prediction(start=pd.to_datetime('2014-01-1'), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int()
ts_forecast = pred_dynamic.predicted_mean
ts_orginal =ts['2014-01-01':]
ts_forecast.plot(color='blue', label='Predict')
ts_orginal.plot(color='red', label='Original')
plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((ts_forecast-ts_orginal)**2)/ts.size))
plt.show()

可以看到，魔性的动态预测效果很差，不适合使用该模型

预测未来5年的房价

#预测未来5年的数据
forecast = results.get_forecast(steps= 60)
# 得到预测的置信区间
forecast_ci = forecast.conf_int()
ts_forecast = forecast.predicted_mean
ts_pred_concat = pd.concat([ts_forecast,forecast_ci],axis=1)
ts_pred_concat.columns = [u'预测值',u'下限',u'上限']
print(ts_pred_concat.head(60))

#绘制时间序列图
ax = ts.plot(label='origin', figsize=(20, 15))
forecast.predicted_mean.plot(ax=ax, label='predict')
ax.fill_between(forecast_ci.index,
                forecast_ci.iloc[:, 0],
                forecast_ci.iloc[:, 1], color='g', alpha=.4)
plt.xticks(fontsize = 36)
plt.yticks(fontsize = 36)
ax.set_xlabel('时间(年)',fontsize=36)
ax.set_ylabel('房价\元',fontsize=36)
plt.legend(loc = 'upper left',fontsize=20)
plt.show()

问题

模型在初始采样时间（2011年左右）拟合效果不佳

参考

https://zhuanlan.zhihu.com/p/127032260