或者
这里的y为时间序列的值,S为季节分量,T为趋势分量,n为噪声。
为了进行分解,除了选择分解类之外,还需要设置一个季节周期(例如,p=1表示年度数据,p=4表示季度数据,p=12表示月度数据等)。
前面提到的经典分解是一种非常幼稚和简单的方法。它具有明显的局限性,如线性,无法捕捉动态季节性和难以处理时间序列中的非平稳性,但是就本文作为演示,这种方法是可以的。
为了进行经典的分解,我们将使用statmodels库中的seasonal_decomposition函数,周期等于24,因为我们处理的是每小时的数据:
vars = {'t2m': 'Air Temperature', 'tp': 'Total Precipitation', 'sp': 'Surface Pressure', 'ssr': 'Surface Net Solar Radiation'}
for var in df.columns:
result = sm.tsa.seasonal_decompose(df[var], model='additive', period = 24)
results_df = pd.DataFrame({'trend': result.trend, 'seasonal': result.seasonal, 'resid': result.resid, 'observed': result.observed})
fig, ax = plt.subplots(ncols = 2, nrows = 2,figsize=(16,9))
ax[0,0].plot(df.index, results_df.trend)
ax[0,0].set_title('Trend')
ax[0,0].set_ylabel('Value')
ax[0,1].plot(df.index, results_df.seasonal)
ax[0,1].set_title('Seasonal')
ax[1,0].plot(df.index, results_df.resid)
ax[1,0].set_title('Residual')
ax[1,0].set_ylabel('Value')
ax[1,0].set_xlabel('time')
ax[1,1].plot(df.index, results_df.observed)
ax[1,1].set_title('Observed')
ax[1,1].set_xlabel('time')
opinionated.set_title_and_suptitle(vars[var], f"Dickey-Fuller test: {round(sm.tsa.stattools.adfuller(df[var])[1],5)}", position_title=[0.45,1],
position_sub_title=[0.95, 1])
plt.tight_layout()
plt.savefig(f'Seasonal_{var}.png')
plt.show()