基于python的數據分解-趨勢-季節性-波動變化

系列文章目錄

前言

時間序列數據的分解，一般分為趨勢項，季節變化項和隨機波動項。可以基于加法或者乘法模型。季節變化呈現出周期變化，因此也叫季節效應(周期）。

一、數據分解步驟

（1）估計時間序列的長期趨勢，一種是通過數據平滑方式進行估計；一種是通過模擬回歸方程進行估計；
（2）去掉時間序列數據的長期趨勢。加法模型則減去，乘法模型即除去；
（3）去掉長期趨勢的時間序列數據，估計時間序列的季節變化；
（4）剩下的即為隨機波動項；

二、使用步驟

1.數據分解方法

import os
#移動平均法
os.chdir("D:/Pythonmatlab學習資料/") 
#改變工作目錄
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
import numpy as npfrom statsmodels.tsa.seasonal import seasonal_decomposetraveller_df = pd.read_csv("NZTravellersDestination.csv", usecols=['Date','China'], parse_dates=['Date'], index_col='Date')deco_muti = seasonal_decompose(traveller_df, model='mutiplicative', extrapolate_trend='freq')new,(ax1,ax2,ax3,ax4) = plt.subplots(4, 1, sharex=True, figsize=(12,8), dpi=150)
ax1.plot(deco_muti.observed, color='r')
ax1.set_ylabel(ylabel="Observed", fontsize=15)
ax2.plot(deco_muti.trend, color='b')
ax2.set_ylabel(ylabel="Trend", fontsize=15)
ax3.plot(deco_muti.seasonal, color='g')
ax3.set_ylabel(ylabel="Seasonal", fontsize=15)
ax4.plot(deco_muti.resid, color='b')
ax4.set_ylabel(ylabel="Resid", fontsize=15)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
plt.tight_layout(); plt.savefig(fname='fig/4_1.png')deco_value = pd.concat([deco_muti.trend, deco_muti.seasonal, deco_muti.resid, deco_muti.observed], axis=1)
deco_value.columns = ['trend', 'season', 'resid', 'actual_values']
deco_value.head()
##%

2.分解項計算方法


import statsmodels.formula.api as smf
#表示線性擬合計算趨勢
df = np.loadtxt("elec_prod.txt")
t = np.arange(1,397)df_t = np.vstack((df,t)).swapaxes(0,1).astype(int)
model_data = pd.DataFrame(df_t,columns=['df','t'])results_f = smf.ols('df~t',data=model_data).fit()
print(results_f.summary().tables[1])
print('std = ',np.std(results_f.resid))fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(model_data, linestyle="-", color='red')
ax.plot(t,1.423e+05 + 499.2576*t, color='blue')
ax.set_ylim((130000, 410000))
ax.set_ylabel(ylabel="Electricity", fontsize=17)
ax.set_xlabel(xlabel="Time", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_2.png')##%表示曲線擬合計算趨勢
df = pd.read_excel('ningxiaGDP.xlsx').rename(columns={'t':'t1'})
df = pd.DataFrame(df['t1'].values**2,columns=['t2']).join(df)results_f = smf.ols('gdp~ 0 + t1+ t2', data=df).fit()
print(results_f.summary().tables[1])
print('std = ', np.std(results_f.resid))from scipy.optimize import curve_fitdf = pd.read_excel('ningxiaGDP.xlsx')
t = df['t'].values
gdp = df['gdp'].valuesdef func(x, b,c):return b*x + c*x**2popt, pcov = curve_fit(func,t,gdp,p0=(1.0,1.0))
print(popt)b = popt[0]
c = popt[1]
residuals = gdp - func(t, b, c)
print(np.std(residuals))t = np.arange(1996, 2016)fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
ax.scatter(y=gdp, x=t, color='blue')
ax.plot(t, results_f.predict())
ax.xaxis.set_major_locator(ticker.MultipleLocator(3))
ax.set_ylabel(ylabel="寧夏地區生產總值",fontsize=17)
ax.set_xlabel(xlabel="時間", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_3.png')##%表示5期移動平均擬合，移動平均法計算趨勢
from statsmodels.tsa.seasonal import seasonal_decompose
nile_ar = np.loadtxt("Nile.txt"); Date = np.arange(1871, 1971)
nile_df = pd.DataFrame({"Date":Date, "Nile":nile_ar})
nile_df.index = nile_df["Date"]nile_df['5-period Moving Avg'] = nile_df['Nile'].rolling(5).mean()fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
nile_df['Nile'].plot(ax=ax, color='b', marker="o", linestyle='--')
nile_df['5-period Moving Avg'].plot(ax=ax, color='r')
ax.legend(loc=1,labels=['Index','Moving average'], fontsize=13)
ax.set_ylabel(ylabel="Index", fontsize=17)
ax.set_xlabel(xlabel="Time", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_4.png')#%%二次移動平均法計算趨勢
gdp_df = pd.read_csv('JDGDP.csv')gdp_df['Moving_Avg_1'] = gdp_df['JDGDP'].rolling(4).mean()
gdp_df['Moving_Avg_2'] = gdp_df['Moving_Avg_1'].rolling(4).mean()gdp_df['at'] = 2*gdp_df['Moving_Avg_1'] - gdp_df['Moving_Avg_2']fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
gdp_df['JDGDP'].plot(ax=ax, color='b',marker="o",linestyle='--')
gdp_df['at'].plot(ax=ax, color='r')
ax.xaxis.set_major_locator(ticker.MultipleLocator(3))
ax.legend(loc=2,labels=['季度GDP','兩次移動平均'], fontsize=13)
ax.set_ylabel(ylabel="中國季度國內生產總值", fontsize=17)
ax.set_xlabel(xlabel="時間", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_5.png')##%指數平滑法，考慮近期變化對現在影響較大，而遠期的變化對現在影響要小一些
from statsmodels.tsa.api import SimpleExpSmoothing
df = np.loadtxt("retail_price_index.txt")
index = pd.date_range(start="1990", end="2021", freq="A")
retail_df = pd.Series(df, index)fit1 = SimpleExpSmoothing(retail_df, initialization_method="heuristic").fit(smoothing_level=0.2, optimized=False)
fcast1 = fit1.forecast(3).rename(r"α=0.2")fit2 = SimpleExpSmoothing(retail_df, initialization_method="heuristic").fit(smoothing_level=0.6, optimized=False)
fcast2 = fit2.forecast(3).rename(r"α=0.6")fit3 = SimpleExpSmoothing(retail_df, initialization_method="estimated").fit()
fcast3 = fit3.forecast(3).rename(r"α=
"% fit3.model.params["smoothing_level"] )
fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(retail_df, marker="o", color="black")
ax.plot(fit1.fittedvalues, marker="8", color="green",linestyle="-.")
(line1,) = ax.plot(fcast1, marker="8", color="green",linestyle="-.")
ax.plot(fit2.fittedvalues, marker="s", color="red",linestyle=":")
(line2,) = ax.plot(fcast2, marker="s", color="red",linestyle=":")
ax.plot(fit3.fittedvalues, marker="p", color="blue",linestyle="--")
(line3,) = ax.plot(fcast3, marker="p", color="blue",linestyle="--")
plt.legend([line1, line2, line3], [fcast1.name, fcast2.name, fcast3.name], fontsize=15)
ax.set_ylabel(ylabel="商品零售價格指數", fontsize=17)
ax.set_xlabel(xlabel="時間", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_6.png')##%對于含有季節變動部分的時間序列為Holt-Winters指數平滑法
from statsmodels.tsa.api import ExponentialSmoothing
df = np.loadtxt("QGDP.txt")
index = pd.date_range(start="2000", end="2021", freq="Q")
QGDP_df = pd.Series(df, index)fit = ExponentialSmoothing(QGDP_df, seasonal_periods=4, trend="add", seasonal="mul",initialization_method="estimated").fit()
simulations = fit.simulate(8, repetitions=1000, error="mul")fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(QGDP_df, marker="o", color="black")
ax.plot(fit.fittedvalues, marker="o", color="blue", linestyle=":")
ax.plot(simulations, marker="o", color="blue", linestyle=":")
ax.set_ylabel(ylabel="國內季度生產總值累計值", fontsize=17)
ax.set_xlabel(xlabel="時間",fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_8.png')##%季節效應，季節指數，用簡單平均法計算的周期內各時期季節性影響的相對數。
df = np.loadtxt("tempdub.txt")
index = pd.date_range(start="1964", end="1976", freq="M")
tempdub_df = pd.Series(df, index)fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(tempdub_df, marker="o", color="blue")
ax.set_ylabel(ylabel="杜比克市月平均氣溫",fontsize=17)
ax.set_xlabel(xlabel="時間",fontsize=17)
plt.xticks(fontsize=15);plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_9.png')t = np.arange(1,13)
SI = np.array([0.36,0.45,0.7,1.00,1.26,1.46,1.55,1.50,1.32,1.10,0.79,0.51])
Season_index = pd.Series(SI, t)fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(Season_index, marker="o", color="blue")
ax.set_ylabel(ylabel="季節指數",fontsize=17)
ax.set_xlabel(xlabel="時間",fontsize=17)
plt.xticks(fontsize=15);plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/4_10.png')