python數(shù)據(jù)分析之時(shí)間序列分析詳情
更新時(shí)間：2022年08月11日 10:20:33 作者：???KaneBurger???????
這篇文章主要介紹了python數(shù)據(jù)分析之時(shí)間序列分析詳情，時(shí)間序列分析是基于隨機(jī)過程理論和數(shù)理統(tǒng)計(jì)學(xué)方法，具體詳細(xì)內(nèi)容介紹，需要的小伙伴可以參考一下
前言
時(shí)間序列分析是基于隨機(jī)過程理論和數(shù)理統(tǒng)計(jì)學(xué)方法：
- 每日的平均氣溫
- 每天的銷售額
- 每月的降水量
時(shí)間序列分析主要通過statsmodel庫的tsa模塊完成：
- 根據(jù)時(shí)間序列的散點(diǎn)圖，自相關(guān)函數(shù)和偏自相關(guān)函數(shù)圖識(shí)別序列是否平穩(wěn)的非隨機(jī)序列，如果是非隨機(jī)序列，觀察其平穩(wěn)性
- 對非平穩(wěn)的時(shí)間序列數(shù)據(jù)采用差分進(jìn)行平滑處理
- 根據(jù)識(shí)別出來的特征建立相應(yīng)的時(shí)間序列模型
- 參數(shù)估計(jì)，檢驗(yàn)是否具有統(tǒng)計(jì)意義
- 假設(shè)檢驗(yàn)，判斷模型的殘差序列是否為白噪聲序列
- 利用已通過檢驗(yàn)的模型進(jìn)行預(yù)測
時(shí)間序列的相關(guān)檢驗(yàn)
白噪聲檢驗(yàn)
如果為白噪聲數(shù)據(jù)(即獨(dú)立分布的隨機(jī)數(shù)據(jù))，說明其沒有任何有用的信息
## 輸出高清圖像 %config InlineBackend.figure_format = 'retina' %matplotlib inline ## 圖像顯示中文的問題 import matplotlib matplotlib.rcParams['axes.unicode_minus']=False import seaborn as sns sns.set(font= "Kaiti",style="ticks",font_scale=1.4)
## 導(dǎo)入會(huì)使用到的相關(guān)庫
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import *
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.tsa.api import SimpleExpSmoothing,Holt,ExponentialSmoothing,AR,ARIMA,ARMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pmdarima as pm
from sklearn.metrics import mean_absolute_error
import pyflux as pf
from fbprophet import Prophet
## 忽略提醒
import warnings
warnings.filterwarnings("ignore")## 讀取時(shí)間序列數(shù)據(jù),該數(shù)據(jù)包含:X1為飛機(jī)乘客數(shù)據(jù),X2為一組隨機(jī)數(shù)據(jù)
df = pd.read_csv("data/chap6/timeserise.csv")
## 查看數(shù)據(jù)的變化趨勢
df.plot(kind = "line",figsize = (10,6))
plt.grid()
plt.title("時(shí)序數(shù)據(jù)")
plt.show()
## 白噪聲檢驗(yàn)Ljung-Box檢驗(yàn)
## 該檢驗(yàn)用來檢查序列是否為隨機(jī)序列，如果是隨機(jī)序列，那它們的值之間沒有任何關(guān)系
## 使用LB檢驗(yàn)來檢驗(yàn)序列是否為白噪聲，原假設(shè)為在延遲期數(shù)內(nèi)序列之間相互獨(dú)立。
lags = [4,8,16,32]
LB = sm.stats.diagnostic.acorr_ljungbox(df["X1"],lags = lags,return_df = True)
print("序列X1的檢驗(yàn)結(jié)果:\n",LB)
LB = sm.stats.diagnostic.acorr_ljungbox(df["X2"],lags = lags,return_df = True)
print("序列X2的檢驗(yàn)結(jié)果:\n",LB)
## 如果P值小于0.05，說明序列之間不獨(dú)立，不是白噪聲
'''
序列X1的檢驗(yàn)結(jié)果:
lb_stat lb_pvalue
4 427.738684 2.817731e-91
8 709.484498 6.496271e-148
16 1289.037076 1.137910e-264
32 1792.523003 0.000000e+00
序列X2的檢驗(yàn)結(jié)果:
lb_stat lb_pvalue
4 1.822771 0.768314
8 8.452830 0.390531
16 15.508599 0.487750
32 28.717743 0.633459
'''在延遲階數(shù)為[4,6,16,32]的情況下，序列X1的LB檢驗(yàn)P值均小于0.05，即該數(shù)據(jù)不是隨機(jī)的。有規(guī)律可循，有分析價(jià)值，而序列X2的LB檢驗(yàn)P值均大于0.05，該數(shù)據(jù)為白噪聲，沒有分析價(jià)值
平穩(wěn)性檢驗(yàn)
時(shí)間序列是否平穩(wěn)，對選擇預(yù)測的數(shù)學(xué)模型非常關(guān)鍵
如果數(shù)據(jù)是平穩(wěn)的，可以使用自回歸平均移動(dòng)模型（ARMA）
如果數(shù)據(jù)是不平穩(wěn)的，可以使用差分移動(dòng)自回歸平均移動(dòng)模型（ARIMA）
## 序列的單位根檢驗(yàn)，即檢驗(yàn)序列的平穩(wěn)性
dftest = adfuller(df["X2"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X2單位根檢驗(yàn)結(jié)果:\n",dfoutput)
dftest = adfuller(df["X1"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X1單位根檢驗(yàn)結(jié)果:\n",dfoutput)
## 對X1進(jìn)行一階差分后的序列進(jìn)行檢驗(yàn)
X1diff = df["X1"].diff().dropna()
dftest = adfuller(X1diff,autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X1一階差分單位根檢驗(yàn)結(jié)果:\n",dfoutput)
## 一階差分后 P值大于0.05, 小于0.1，可以認(rèn)為其是平穩(wěn)的
'''
X2單位根檢驗(yàn)結(jié)果:
adf -1.124298e+01
p-value 1.788000e-20
usedlag 0.000000e+00
Number of Observations Used 1.430000e+02
dtype: float64
X1單位根檢驗(yàn)結(jié)果:
adf 0.815369
p-value 0.991880
usedlag 13.000000
Number of Observations Used 130.000000
dtype: float64
X1一階差分單位根檢驗(yàn)結(jié)果:
adf -2.829267
p-value 0.054213
usedlag 12.000000
Number of Observations Used 130.000000
dtype: float64
'''序列X2的檢驗(yàn)P值小于0.05，說明X2是一個(gè)平穩(wěn)時(shí)間序列（該序列是白噪聲，白噪聲序列是平穩(wěn)序列）
序列X1的檢驗(yàn)P值遠(yuǎn)大于0.05，說明不平穩(wěn)，而其一階差分后的結(jié)果，P值大于0.05，但小于0.1，可以認(rèn)為平穩(wěn)
針對數(shù)據(jù)的平穩(wěn)性檢驗(yàn)，還可以使用KPSS檢驗(yàn)，其原假設(shè)該序列是平穩(wěn)的，該檢驗(yàn)可以用kpss()函數(shù)完成
## KPSS檢驗(yàn)的原假設(shè)為：序列x是平穩(wěn)的。
## 對序列X2使用KPSS檢驗(yàn)平穩(wěn)性
dfkpss = kpss(df["X2"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X2 KPSS檢驗(yàn)結(jié)果:\n",dfoutput)
## 接受序列平穩(wěn)的原假設(shè)
## 對序列X1使用KPSS檢驗(yàn)平穩(wěn)性
dfkpss = kpss(df["X1"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1 KPSS檢驗(yàn)結(jié)果:\n",dfoutput)
## 拒絕序列平穩(wěn)的原假設(shè)
## 對序列X1使用KPSS檢驗(yàn)平穩(wěn)性
dfkpss = kpss(X1diff)
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1一階差分KPSS檢驗(yàn)結(jié)果:\n",dfoutput)
## 接受序列平穩(wěn)的原假設(shè)
'''
X2 KPSS檢驗(yàn)結(jié)果:
kpss_stat 0.087559
p-value 0.100000
usedlag 14.000000
dtype: float64
X1 KPSS檢驗(yàn)結(jié)果:
kpss_stat 1.052175
p-value 0.010000
usedlag 14.000000
dtype: float64
X1一階差分KPSS檢驗(yàn)結(jié)果:
kpss_stat 0.05301
p-value 0.10000
usedlag 14.00000
dtype: float64
'''ARIMA(p,d,q)模型
## 檢驗(yàn)ARIMA模型的參數(shù)d
X1d = pm.arima.ndiffs(df["X1"], alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法對序列X1的參數(shù)d取值進(jìn)行預(yù)測,d = ",X1d)
X1diffd = pm.arima.ndiffs(X1diff, alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法對序列X1一階差分后的參數(shù)d取值進(jìn)行預(yù)測,d = ",X1diffd)
X2d = pm.arima.ndiffs(df["X2"], alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法對序列X2的參數(shù)d取值進(jìn)行預(yù)測,d = ",X2d)
'''
使用KPSS方法對序列X1的參數(shù)d取值進(jìn)行預(yù)測,d = 1
使用KPSS方法對序列X1一階差分后的參數(shù)d取值進(jìn)行預(yù)測,d = 0
使用KPSS方法對序列X1的參數(shù)d取值進(jìn)行預(yù)測,d = 0
'''針對SARIMA模型，還有一個(gè)季節(jié)性平穩(wěn)性參數(shù)D
## 檢驗(yàn)SARIMA模型的參數(shù)季節(jié)階數(shù)D
X1d = pm.arima.nsdiffs(df["X1"], 12, max_D=2)
print("對序列X1的季節(jié)階數(shù)D取值進(jìn)行預(yù)測,D = ",X1d)
X1diffd = pm.arima.nsdiffs(X1diff, 12, max_D=2)
print("序列X1一階差分后的季節(jié)階數(shù)D取值進(jìn)行預(yù)測,D = ",X1diffd)
'''
對序列X1的季節(jié)階數(shù)D取值進(jìn)行預(yù)測,D = 1
序列X1一階差分后的季節(jié)階數(shù)D取值進(jìn)行預(yù)測,D = 1
'''自相關(guān)和偏相關(guān)分析
## 對隨機(jī)序列X2進(jìn)行自相關(guān)和偏相關(guān)分析可視化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X2"],"r-")
plt.grid()
plt.title("X2序列波動(dòng)")
ax = fig.add_subplot(1,3,2)
plot_acf(df["X2"], lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(df["X2"], lags=60,ax = ax)
plt.grid()
plt.tight_layout()
plt.show()
在圖像中滯后0表示自己和自己的相關(guān)性，恒等于1。不用于確定p和q。
## 對非隨機(jī)序列X1進(jìn)行自相關(guān)和偏相關(guān)分析可視化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X1"],"r-")
plt.grid()
plt.title("X1序列波動(dòng)")
ax = fig.add_subplot(1,3,2)
plot_acf(df["X1"], lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(df["X1"], lags=60,ax = ax)
plt.ylim([-1,1])
plt.grid()
plt.tight_layout()
plt.show()
## 對非隨機(jī)序列X1一階差分后的序列進(jìn)行自相關(guān)和偏相關(guān)分析可視化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(X1diff,"r-")
plt.grid()
plt.title("X1序列一階差分后波動(dòng)")
ax = fig.add_subplot(1,3,2)
plot_acf(X1diff, lags=60,ax = ax)
plt.grid()
ax = fig.add_subplot(1,3,3)
plot_pacf(X1diff, lags=60,ax = ax)
plt.grid()
plt.tight_layout()
plt.show()
ARMA(p,q)中，自相關(guān)系數(shù)的滯后，對應(yīng)著參數(shù)q；偏相關(guān)系數(shù)的滯后對應(yīng)著參數(shù)p。
## 時(shí)間序列的分解
## 通過觀察序列X1，可以發(fā)現(xiàn)其既有上升的趨勢，也有周期性的趨勢，所以可以將該序列進(jìn)行分解
## 使用乘法模型分解結(jié)果(通常適用于有增長趨勢的序列)
X1decomp = pm.arima.decompose(df["X1"].values,"multiplicative", m=12)
## 可視化出分解的結(jié)果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
show=False)
ax[0].set_title("乘法模型分解結(jié)果")
plt.show()
## 使用加法模型分解結(jié)果(通常適用于平穩(wěn)趨勢的序列)
X1decomp = pm.arima.decompose(X1diff.values,"additive", m=12)
## 可視化出分解的結(jié)果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
show=False)
ax[0].set_title("加法模型分解結(jié)果")
plt.show()
移動(dòng)平均算法
## 數(shù)據(jù)準(zhǔn)備 ## 對序列X1進(jìn)行切分，后面的24個(gè)數(shù)據(jù)用于測試集 train = pd.DataFrame(df["X1"][0:120]) test = pd.DataFrame(df["X1"][120:]) ## 可視化切分后的數(shù)據(jù) train["X1"].plot(figsize=(14,7), title= "乘客數(shù)量數(shù)據(jù)",label = "X1 train") test["X1"].plot(label = "X1 test") plt.legend() plt.grid() plt.show()
print(train.shape) print(test.shape) df["X1"].shape ''' (120, 1) (24, 1) (144,) '''
簡單移動(dòng)平均法
## 簡單移動(dòng)平均進(jìn)行預(yù)測
y_hat_avg = test.copy(deep = False)
y_hat_avg["moving_avg_forecast"] = train["X1"].rolling(24).mean().iloc[-1]
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["moving_avg_forecast"].plot(style="g--o", lw=2,
label="移動(dòng)平均預(yù)測")
plt.legend()
plt.grid()
plt.title("簡單移動(dòng)平均預(yù)測")
plt.show()
## 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("預(yù)測絕對值誤差:",mean_absolute_error(test["X1"],y_hat_avg["moving_avg_forecast"]))
'''
預(yù)測絕對值誤差: 82.55208333333336
'''簡單指數(shù)平滑法
## 數(shù)據(jù)準(zhǔn)備
y_hat_avg = test.copy(deep = False)
## 模型構(gòu)建
model1 = SimpleExpSmoothing(train["X1"].values).fit(smoothing_level=0.15)
y_hat_avg["exp_smooth_forecast1"] = model1.forecast(len(test))
model2 = SimpleExpSmoothing(train["X1"].values).fit(smoothing_level=0.5)
y_hat_avg["exp_smooth_forecast2"] = model2.forecast(len(test))
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["exp_smooth_forecast1"].plot(style="g--o", lw=2,
label="smoothing_level=0.15")
y_hat_avg["exp_smooth_forecast2"].plot(style="g--s", lw=2,
label="smoothing_level=0.5")
plt.legend()
plt.grid()
plt.title("簡單指數(shù)平滑預(yù)測")
plt.show()
## 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("smoothing_level=0.15,預(yù)測絕對值誤差:",
mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast1"]))
print("smoothing_level=0.5,預(yù)測絕對值誤差:",
mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast2"]))
smoothing_level=0.15,預(yù)測絕對值誤差: 81.10115706423566
smoothing_level=0.5,預(yù)測絕對值誤差: 106.813228720506
霍爾特(Holt)線性趨勢法
## 數(shù)據(jù)準(zhǔn)備
y_hat_avg = test.copy(deep = False)
## 模型構(gòu)建
model1 = Holt(train["X1"].values).fit(smoothing_level=0.1, smoothing_slope = 0.05)
y_hat_avg["holt_forecast1"] = model1.forecast(len(test))
model2 = Holt(train["X1"].values).fit(smoothing_level=0.1, smoothing_slope = 0.25)
y_hat_avg["holt_forecast2"] = model2.forecast(len(test))
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["holt_forecast1"].plot(style="g--o", lw=2,
label="Holt線性趨勢法(1)")
y_hat_avg["holt_forecast2"].plot(style="g--s", lw=2,
label="Holt線性趨勢法(2)")
plt.legend()
plt.grid()
plt.title("Holt線性趨勢法預(yù)測")
plt.show()
## 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("smoothing_slope = 0.05,預(yù)測絕對值誤差:",
mean_absolute_error(test["X1"],y_hat_avg["holt_forecast1"]))
print("smoothing_slope = 0.25,預(yù)測絕對值誤差:",
mean_absolute_error(test["X1"],y_hat_avg["holt_forecast2"]))
smoothing_slope = 0.05,預(yù)測絕對值誤差: 54.727467142360275
smoothing_slope = 0.25,預(yù)測絕對值誤差: 69.79052992788556
Holt-Winters季節(jié)性預(yù)測模型
## 數(shù)據(jù)準(zhǔn)備
y_hat_avg = test.copy(deep = False)
## 模型構(gòu)建
model1 = ExponentialSmoothing(train["X1"].values,
seasonal_periods=12, # 周期性為12
trend="add", seasonal="add").fit()
y_hat_avg["holt_winter_forecast1"] = model1.forecast(len(test))
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat_avg["holt_winter_forecast1"].plot(style="g--o", lw=2,
label="Holt-Winters")
plt.legend()
plt.grid()
plt.title("Holt-Winters季節(jié)性預(yù)測模型")
plt.show()
## 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("Holt-Winters季節(jié)性預(yù)測模型,預(yù)測絕對值誤差:",
mean_absolute_error(test["X1"],y_hat_avg["holt_winter_forecast1"]))
Holt-Winters季節(jié)性預(yù)測模型,預(yù)測絕對值誤差: 30.06821059070873
ARIMA模型
- 注意針對乘客數(shù)據(jù)X1，使用AR模型或者ARMA模型進(jìn)行預(yù)測，并不是非常的合適，
- 這里是使用AR和ARMA模型進(jìn)行預(yù)測的目的主要是為了和更好的模型預(yù)測結(jié)果進(jìn)行對比
## 使用AR模型對乘客數(shù)據(jù)進(jìn)行預(yù)測 ## 經(jīng)過前面序列的偏相關(guān)系數(shù)的可視化結(jié)果，使用AR(2)模型可對序列進(jìn)行建模 ## 數(shù)據(jù)準(zhǔn)備 y_hat = test.copy(deep = False) ## 模型構(gòu)建 ar_model = ARMA(train["X1"].values,order = (2,0)).fit() ## 輸出擬合模型的結(jié)果 print(ar_model.summary()) ## AIC=1141.989；BIC= 1153.138；兩個(gè)系數(shù)是顯著的
## 查看模型的擬合殘差分布
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1)
plt.plot(ar_model.resid)
plt.title("AR(2)殘差曲線")
## 檢查殘差是否符合正太分布
ax = fig.add_subplot(1,2,2)
sm.qqplot(ar_model.resid, line='q', ax=ax)
plt.title("AR(2)殘差Q-Q圖")
plt.tight_layout()
plt.show()
## 預(yù)測未來24個(gè)數(shù)據(jù)，并輸出95%置信區(qū)間
pre, se, conf = ar_model.forecast(24, alpha=0.05)
## 整理數(shù)據(jù)
y_hat["ar2_pre"] = pre
y_hat["ar2_pre_lower"] = conf[:,0]
y_hat["ar2_pre_upper"] = conf[:,1]
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat["ar2_pre"].plot(style="g--o", lw=2,label="AR(2)")
## 可視化出置信區(qū)間
plt.fill_between(y_hat.index, y_hat["ar2_pre_lower"],
y_hat["ar2_pre_upper"],color='k',alpha=.15,
label = "95%置信區(qū)間")
plt.legend()
plt.grid()
plt.title("AR(2)模型")
plt.show()
# 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("AR(2)模型預(yù)測的絕對值誤差:",
mean_absolute_error(test["X1"],y_hat["ar2_pre"]))
AR(2)模型預(yù)測的絕對值誤差: 165.79608244918572
可以發(fā)現(xiàn)使用AR(2)的預(yù)測效果并不好
ARMA模型
## 嘗試使用ARMA模型進(jìn)行預(yù)測 ## 根據(jù)前面的自相關(guān)系數(shù)和偏相關(guān)系數(shù),為了降低模型的復(fù)雜讀,可以使用ARMA(2,1) ## 數(shù)據(jù)準(zhǔn)備 y_hat = test.copy(deep = False) ## 模型構(gòu)建 arma_model = ARMA(train["X1"].values,order = (2,1)).fit() ## 輸出擬合模型的結(jié)果 print(arma_model.summary()) ## AIC=1141.989；BIC= 1153.138；兩個(gè)系數(shù)是顯著的
## 查看模型的擬合殘差分布
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1)
plt.plot(arma_model.resid)
plt.title("ARMA(2,1)殘差曲線")
## 檢查殘差是否符合正太分布
ax = fig.add_subplot(1,2,2)
sm.qqplot(arma_model.resid, line='q', ax=ax)
plt.title("ARMA(2,1)殘差Q-Q圖")
plt.tight_layout()
plt.show()
## 預(yù)測未來24個(gè)數(shù)據(jù)，并輸出95%置信區(qū)間
pre, se, conf = arma_model.forecast(24, alpha=0.05)
## 整理數(shù)據(jù)
y_hat["arma_pre"] = pre
y_hat["arma_pre_lower"] = conf[:,0]
y_hat["arma_pre_upper"] = conf[:,1]
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat["arma_pre"].plot(style="g--o", lw=2,label="ARMA(2,1)")
## 可視化出置信區(qū)間
plt.fill_between(y_hat.index, y_hat["arma_pre_lower"],
y_hat["arma_pre_upper"],color='k',alpha=.15,
label = "95%置信區(qū)間")
plt.legend()
plt.grid()
plt.title("ARMA(2,1)模型")
plt.show()
# 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("ARMA模型預(yù)測的絕對值誤差:",
mean_absolute_error(test["X1"],y_hat["arma_pre"]))
ARMA模型預(yù)測的絕對值誤差: 147.26531763335154
針對ARMA模型自動(dòng)選擇合適的參數(shù)
## 自動(dòng)搜索合適的參數(shù)
model = pm.auto_arima(train["X1"].values,
start_p=1, start_q=1, # p,q的開始值
max_p=12, max_q=12, # 最大的p和q
d = 0, # 尋找ARMA模型參數(shù)
m=1, # 序列的周期
seasonal=False, # 沒有季節(jié)性趨勢
trace=True,error_action='ignore',
suppress_warnings=True, stepwise=True)
print(model.summary())## 使用ARMA(3,2)對測試集進(jìn)行預(yù)測
pre, conf = model.predict(n_periods=24, alpha=0.05,
return_conf_int=True)
## 可視化ARMA(3,2)的預(yù)測結(jié)果，整理數(shù)據(jù)
y_hat["arma_pre"] = pre
y_hat["arma_pre_lower"] = conf[:,0]
y_hat["arma_pre_upper"] = conf[:,1]
## 可視化出預(yù)測結(jié)果
plt.figure(figsize=(14,7))
train["X1"].plot(figsize=(14,7),label = "X1 train")
test["X1"].plot(label = "X1 test")
y_hat["arma_pre"].plot(style="g--o", lw=2,label="ARMA(3,2)")
## 可視化出置信區(qū)間
plt.fill_between(y_hat.index, y_hat["arma_pre_lower"],
y_hat["arma_pre_upper"],color='k',alpha=.15,
label = "95%置信區(qū)間")
plt.legend()
plt.grid()
plt.title("ARMA(3,2)模型")
plt.show()
# 計(jì)算預(yù)測結(jié)果和真實(shí)值的誤差
print("ARMA模型預(yù)測的絕對值誤差:",
mean_absolute_error(test["X1"],y_hat["arma_pre"]))
ARMA模型預(yù)測的絕對值誤差: 158.11464180972925
可以發(fā)現(xiàn)使用自動(dòng)ARMA(3,2)模型的效果并沒有ARMA(2,1)的預(yù)測效果好
時(shí)序數(shù)據(jù)的異常值檢測
可以將突然增大或突然減小的數(shù)據(jù)無規(guī)律看作異常值
## 使用prophet檢測時(shí)間序列是否有異常值
## 從1991年2月到2005年5月，每周提供美國成品汽車汽油產(chǎn)品的時(shí)間序列（每天數(shù)千桶）
## 數(shù)據(jù)準(zhǔn)備
data = pm.datasets.load_gasoline()
datadf = pd.DataFrame({"y":data})
datadf["ds"] = pd.date_range(start="1991-2",periods=len(data),freq="W")
## 可視化時(shí)間序列的變化情況
datadf.plot(x = "ds",y = "y",style = "b-o",figsize=(14,7))
plt.grid()
plt.title("時(shí)間序列數(shù)據(jù)的波動(dòng)情況")
plt.show()
## 對該數(shù)據(jù)建立一個(gè)時(shí)間序列模型
np.random.seed(1234) ## 設(shè)置隨機(jī)數(shù)種子
model = Prophet(growth="linear",daily_seasonality = False,
weekly_seasonality=False,
seasonality_mode = 'multiplicative',
interval_width = 0.95, ## 獲取95%的置信區(qū)間
)
model = model.fit(datadf) # 使用數(shù)據(jù)擬合模型
forecast = model.predict(datadf) # 使用模型對數(shù)據(jù)進(jìn)行預(yù)測
forecast["y"] = datadf["y"].reset_index(drop = True)
forecast[["ds","y","yhat","yhat_lower","yhat_upper"]].head()| ds | y | yhat | yhat_lower | yhat_upper | |
|---|---|---|---|---|---|
| 0 | 1991-02-03 | 6621.0 | 6767.051491 | 6294.125979 | 7303.352309 |
| 1 | 1991-02-10 | 6433.0 | 6794.736479 | 6299.430616 | 7305.414252 |
| 2 | 1991-02-17 | 6582.0 | 6855.096282 | 6352.579489 | 7379.717614 |
| 3 | 1991-02-24 | 7224.0 | 6936.976642 | 6415.157617 | 7445.523000 |
| 4 | 1991-03-03 | 6875.0 | 6990.511503 | 6489.781400 | 7488.240435 |
## 根據(jù)模型預(yù)測值的置信區(qū)間"yhat_lower"和"yhat_upper"判斷樣本是否為異常值
def outlier_detection(forecast):
index = np.where((forecast["y"] <= forecast["yhat_lower"])|
(forecast["y"] >= forecast["yhat_upper"]),True,False)
return index
outlier_index = outlier_detection(forecast)
outlier_df = datadf[outlier_index]
print("異常值的數(shù)量為:",np.sum(outlier_index))
'''
異常值的數(shù)量為: 38
'''## 可視化異常值的結(jié)果
fig, ax = plt.subplots()
## 可視化預(yù)測值
forecast.plot(x = "ds",y = "yhat",style = "b-",figsize=(14,7),
label = "預(yù)測值",ax=ax)
## 可視化出置信區(qū)間
ax.fill_between(forecast["ds"].values, forecast["yhat_lower"],
forecast["yhat_upper"],color='b',alpha=.2,
label = "95%置信區(qū)間")
forecast.plot(kind = "scatter",x = "ds",y = "y",c = "k",
s = 20,label = "原始數(shù)據(jù)",ax = ax)
## 可視化出異常值的點(diǎn)
outlier_df.plot(x = "ds",y = "y",style = "rs",ax = ax,
label = "異常值")
plt.legend(loc = 2)
plt.grid()
plt.title("時(shí)間序列異常值檢測結(jié)果")
plt.show()
異常值大部分都在置信區(qū)間外
到此這篇關(guān)于python數(shù)據(jù)分析之時(shí)間序列分析詳情的文章就介紹到這了,更多相關(guān)python時(shí)間序列分析內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！
相關(guān)文章
解決使用Pandas 讀取超過65536行的Excel文件問題
這篇文章主要介紹了解決使用Pandas 讀取超過65536行的Excel文件問題，具有很好的參考價(jià)值，希望對大家有所幫助。一起跟隨小編過來看看吧2020-11-11
這篇文章主要為大家介紹了Python列表list的使用全面解析，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進(jìn)步，早日升職加薪2023-06-06
python學(xué)習(xí)--使用QQ郵箱發(fā)送郵件代碼實(shí)例
這篇文章主要介紹了python使用QQ郵箱發(fā)送郵件，文中通過示例代碼介紹的非常詳細(xì)，對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2019-04-04
今天小編就為大家分享一篇Tensorflow累加的實(shí)現(xiàn)案例，具有很好的參考價(jià)值，希望對大家有所幫助。一起跟隨小編過來看看吧2020-02-02
從零開始學(xué)習(xí)Python與BeautifulSoup網(wǎng)頁數(shù)據(jù)抓取
想要從零開始學(xué)習(xí)Python和BeautifulSoup網(wǎng)頁數(shù)據(jù)抓?。勘局改蠈槟闾峁┖唵我锥闹笇?dǎo),讓你掌握這兩個(gè)強(qiáng)大的工具,不管你是初學(xué)者還是有經(jīng)驗(yàn)的開發(fā)者,本指南都能幫助你快速入門并提升技能,不要錯(cuò)過這個(gè)機(jī)會(huì),開始你的編程之旅吧！2024-01-01
這篇文章主要介紹了Python爬取豆瓣視頻信息代碼實(shí)例,文中通過示例代碼介紹的非常詳細(xì)，對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2019-11-11
今天小編就為大家分享一篇python中多層嵌套列表的拆分方法，具有很好參考價(jià)值，希望對大家有所幫助。一起跟隨小編過來看看吧2018-07-07
PyQt5設(shè)置登錄界面及界面美化的實(shí)現(xiàn)
這篇文章主要介紹了PyQt5設(shè)置登錄界面及界面美化的實(shí)現(xiàn)，文中通過示例代碼介紹的非常詳細(xì)，對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2021-03-03
這篇文章主要介紹了安裝2019Pycharm最新版本的教程詳解,非常不錯(cuò)，具有一定的參考借鑒價(jià)值，需要的朋友可以參考下2019-10-10
這篇文章主要為大家詳細(xì)介紹了python如何把嵌套列表轉(zhuǎn)變成普通列表，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2018-03-03