python資料分析之時間序列分析詳情

2022-08-11 14:01:36

前言

時間序列分析是基於隨機過程理論和數理統計學方法：

每日的平均氣溫
每天的銷售額
每月的降水量

時間序列分析主要通過statsmodel庫的tsa模組完成：

根據時間序列的散點圖，自相關函數和偏自相關函數圖識別序列是否平穩的非隨機序列，如果是非隨機序列，觀察其平穩性
對非平穩的時間序列資料採用差分進行平滑處理
根據識別出來的特徵建立相應的時間序列模型
引數估計，檢驗是否具有統計意義
假設檢驗，判斷模型的殘差序列是否為白噪聲序列
利用已通過檢驗的模型進行預測

時間序列的相關檢驗

白噪聲檢驗

如果為白噪聲資料(即獨立分佈的亂資料)，說明其沒有任何有用的資訊

## 輸出高清影象
%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)

## 匯入會使用到的相關庫
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")

## 讀取時間序列資料,該資料包含:X1為飛機乘客資料,X2為一組亂資料
df = pd.read_csv("data/chap6/timeserise.csv")
## 檢視資料的變化趨勢
df.plot(kind = "line",figsize = (10,6))
plt.grid()
plt.title("時序資料")
plt.show()

## 白噪聲檢驗Ljung-Box檢驗
## 該檢驗用來檢查序列是否為隨機序列，如果是隨機序列，那它們的值之間沒有任何關係
## 使用LB檢驗來檢驗序列是否為白噪聲，原假設為在延遲期數內序列之間相互獨立。
lags = [4,8,16,32]
LB = sm.stats.diagnostic.acorr_ljungbox(df["X1"],lags = lags,return_df = True)
print("序列X1的檢驗結果:n",LB)
LB = sm.stats.diagnostic.acorr_ljungbox(df["X2"],lags = lags,return_df = True)
print("序列X2的檢驗結果:n",LB)

## 如果P值小於0.05，說明序列之間不獨立，不是白噪聲

'''
序列X1的檢驗結果:
         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的檢驗結果:
       lb_stat  lb_pvalue
4    1.822771   0.768314
8    8.452830   0.390531
16  15.508599   0.487750
32  28.717743   0.633459
'''

在延遲階數為[4,6,16,32]的情況下，序列X1的LB檢驗P值均小於0.05，即該資料不是隨機的。有規律可循，有分析價值，而序列X2的LB檢驗P值均大於0.05，該資料為白噪聲，沒有分析價值

平穩性檢驗

時間序列是否平穩，對選擇預測的數學模型非常關鍵

如果資料是平穩的，可以使用自迴歸平均移動模型（ARMA）

如果資料是不平穩的，可以使用差分移動自迴歸平均移動模型（ARIMA）

## 序列的單位根檢驗，即檢驗序列的平穩性
dftest = adfuller(df["X2"],autolag='BIC')
dfoutput = pd.Series(dftest[0:4], index=['adf','p-value','usedlag','Number of Observations Used'])
print("X2單位根檢驗結果: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單位根檢驗結果:n",dfoutput)

## 對X1進行一階差分後的序列進行檢驗
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一階差分單位根檢驗結果:n",dfoutput)

## 一階差分後 P值大於0.05, 小於0.1，可以認為其是平穩的
'''
X2單位根檢驗結果:
 adf                           -1.124298e+01
p-value                        1.788000e-20
usedlag                        0.000000e+00
Number of Observations Used    1.430000e+02
dtype: float64
X1單位根檢驗結果:
 adf                              0.815369
p-value                          0.991880
usedlag                         13.000000
Number of Observations Used    130.000000
dtype: float64
X1一階差分單位根檢驗結果:
 adf                             -2.829267
p-value                          0.054213
usedlag                         12.000000
Number of Observations Used    130.000000
dtype: float64
'''

序列X2的檢驗P值小於0.05，說明X2是一個平穩時間序列（該序列是白噪聲，白噪聲序列是平穩序列）

序列X1的檢驗P值遠大於0.05，說明不平穩，而其一階差分後的結果，P值大於0.05，但小於0.1，可以認為平穩

針對資料的平穩性檢驗，還可以使用KPSS檢驗，其原假設該序列是平穩的，該檢驗可以用kpss()函數完成

## KPSS檢驗的原假設為：序列x是平穩的。

## 對序列X2使用KPSS檢驗平穩性
dfkpss = kpss(df["X2"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X2 KPSS檢驗結果:n",dfoutput)
## 接受序列平穩的原假設
## 對序列X1使用KPSS檢驗平穩性
dfkpss = kpss(df["X1"])
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1 KPSS檢驗結果:n",dfoutput)
## 拒絕序列平穩的原假設

## 對序列X1使用KPSS檢驗平穩性
dfkpss = kpss(X1diff)
dfoutput = pd.Series(dfkpss[0:3], index=["kpss_stat"," p-value"," usedlag"])
print("X1一階差分KPSS檢驗結果:n",dfoutput)
## 接受序列平穩的原假設

'''
X2 KPSS檢驗結果:
 kpss_stat     0.087559
 p-value      0.100000
 usedlag     14.000000
dtype: float64
X1 KPSS檢驗結果:
 kpss_stat     1.052175
 p-value      0.010000
 usedlag     14.000000
dtype: float64
X1一階差分KPSS檢驗結果:
 kpss_stat     0.05301
 p-value      0.10000
 usedlag     14.00000
dtype: float64
'''

ARIMA(p,d,q)模型

## 檢驗ARIMA模型的引數d
X1d = pm.arima.ndiffs(df["X1"], alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法對序列X1的引數d取值進行預測,d = ",X1d)

X1diffd = pm.arima.ndiffs(X1diff, alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法對序列X1一階差分後的引數d取值進行預測,d = ",X1diffd)

X2d = pm.arima.ndiffs(df["X2"], alpha=0.05, test="kpss", max_d=3)
print("使用KPSS方法對序列X2的引數d取值進行預測,d = ",X2d)

'''
使用KPSS方法對序列X1的引數d取值進行預測,d =  1
使用KPSS方法對序列X1一階差分後的引數d取值進行預測,d =  0
使用KPSS方法對序列X1的引數d取值進行預測,d =  0
'''

針對SARIMA模型，還有一個季節性平穩性引數D

## 檢驗SARIMA模型的引數季節階數D
X1d = pm.arima.nsdiffs(df["X1"], 12, max_D=2)
print("對序列X1的季節階數D取值進行預測,D = ",X1d)
X1diffd = pm.arima.nsdiffs(X1diff, 12, max_D=2)
print("序列X1一階差分後的季節階數D取值進行預測,D = ",X1diffd)

'''
對序列X1的季節階數D取值進行預測,D =  1
序列X1一階差分後的季節階數D取值進行預測,D =  1
'''

自相關和偏相關分析

## 對隨機序列X2進行自相關和偏相關分析視覺化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X2"],"r-")
plt.grid()
plt.title("X2序列波動")
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表示自己和自己的相關性，恆等於1。不用於確定p和q。

## 對非隨機序列X1進行自相關和偏相關分析視覺化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(df["X1"],"r-")
plt.grid()
plt.title("X1序列波動")
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()

## 對非隨機序列X1一階差分後的序列進行自相關和偏相關分析視覺化
fig = plt.figure(figsize=(16,5))
plt.subplot(1,3,1)
plt.plot(X1diff,"r-")
plt.grid()
plt.title("X1序列一階差分後波動")
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)中，自相關係數的滯後，對應著引數q；偏相關係數的滯後對應著引數p。

## 時間序列的分解
## 通過觀察序列X1，可以發現其既有上升的趨勢，也有周期性的趨勢，所以可以將該序列進行分解
## 使用乘法模型分解結果(通常適用於有增長趨勢的序列)
X1decomp = pm.arima.decompose(df["X1"].values,"multiplicative", m=12)
## 視覺化出分解的結果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
                              show=False)
ax[0].set_title("乘法模型分解結果")
plt.show()

## 使用加法模型分解結果(通常適用於平穩趨勢的序列)
X1decomp = pm.arima.decompose(X1diff.values,"additive", m=12)
## 視覺化出分解的結果
ax = pm.utils.decomposed_plot(X1decomp,figure_kwargs = {"figsize": (10, 6)},
                              show=False)
ax[0].set_title("加法模型分解結果")
plt.show()

移動平均演演算法

## 資料準備
## 對序列X1進行切分，後面的24個資料用於測試集
train = pd.DataFrame(df["X1"][0:120])
test = pd.DataFrame(df["X1"][120:])
## 視覺化切分後的資料
train["X1"].plot(figsize=(14,7), title= "乘客數量資料",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,)
'''

簡單移動平均法

## 簡單移動平均進行預測
y_hat_avg = test.copy(deep = False)
y_hat_avg["moving_avg_forecast"] =  train["X1"].rolling(24).mean().iloc[-1]
## 視覺化出預測結果
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="移動平均預測")
plt.legend()
plt.grid()
plt.title("簡單移動平均預測")
plt.show()

## 計算預測結果和真實值的誤差
print("預測絕對值誤差:",mean_absolute_error(test["X1"],y_hat_avg["moving_avg_forecast"]))
'''
預測絕對值誤差: 82.55208333333336
'''

簡單指數平滑法

## 資料準備
y_hat_avg = test.copy(deep = False)
## 模型構建
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))

## 視覺化出預測結果
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("簡單指數平滑預測")
plt.show()

## 計算預測結果和真實值的誤差
print("smoothing_level=0.15,預測絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast1"]))
print("smoothing_level=0.5,預測絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat_avg["exp_smooth_forecast2"]))

smoothing_level=0.15,預測絕對值誤差: 81.10115706423566

smoothing_level=0.5,預測絕對值誤差: 106.813228720506

霍爾特(Holt)線性趨勢法

## 資料準備
y_hat_avg = test.copy(deep = False)
## 模型構建
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))
## 視覺化出預測結果
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線性趨勢法預測")
plt.show()

## 計算預測結果和真實值的誤差
print("smoothing_slope = 0.05,預測絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat_avg["holt_forecast1"]))
print("smoothing_slope = 0.25,預測絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat_avg["holt_forecast2"]))

smoothing_slope = 0.05,預測絕對值誤差: 54.727467142360275

smoothing_slope = 0.25,預測絕對值誤差: 69.79052992788556

Holt-Winters季節性預測模型

## 資料準備
y_hat_avg = test.copy(deep = False)
## 模型構建
model1 = ExponentialSmoothing(train["X1"].values,
                              seasonal_periods=12, # 週期性為12  
                              trend="add", seasonal="add").fit()
y_hat_avg["holt_winter_forecast1"] = model1.forecast(len(test))
## 視覺化出預測結果
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季節性預測模型")
plt.show()
## 計算預測結果和真實值的誤差
print("Holt-Winters季節性預測模型,預測絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat_avg["holt_winter_forecast1"]))

Holt-Winters季節性預測模型,預測絕對值誤差: 30.06821059070873

ARIMA模型

注意針對乘客資料X1，使用AR模型或者ARMA模型進行預測，並不是非常的合適，
這裡是使用AR和ARMA模型進行預測的目的主要是為了和更好的模型預測結果進行對比

## 使用AR模型對乘客資料進行預測 

## 經過前面序列的偏相關係數的視覺化結果，使用AR(2)模型可對序列進行建模
## 資料準備
y_hat = test.copy(deep = False)
## 模型構建
ar_model = ARMA(train["X1"].values,order = (2,0)).fit()
## 輸出擬合模型的結果
print(ar_model.summary())

## AIC=1141.989；BIC= 1153.138；兩個係數是顯著的

## 檢視模型的擬合殘差分佈
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()

## 預測未來24個資料，並輸出95%置信區間
pre, se, conf = ar_model.forecast(24, alpha=0.05)  
## 整理資料
y_hat["ar2_pre"] = pre
y_hat["ar2_pre_lower"] = conf[:,0]
y_hat["ar2_pre_upper"] = conf[:,1]
## 視覺化出預測結果
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)")
## 視覺化出置信區間
plt.fill_between(y_hat.index, y_hat["ar2_pre_lower"], 
                 y_hat["ar2_pre_upper"],color='k',alpha=.15,
                 label = "95%置信區間")
plt.legend()
plt.grid()
plt.title("AR(2)模型")
plt.show()
# 計算預測結果和真實值的誤差
print("AR(2)模型預測的絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat["ar2_pre"]))

AR(2)模型預測的絕對值誤差: 165.79608244918572

可以發現使用AR(2)的預測效果並不好

ARMA模型

## 嘗試使用ARMA模型進行預測

## 根據前面的自相關係數和偏相關係數,為了降低模型的複雜讀,可以使用ARMA(2,1)

## 資料準備
y_hat = test.copy(deep = False)
## 模型構建
arma_model = ARMA(train["X1"].values,order = (2,1)).fit()
## 輸出擬合模型的結果
print(arma_model.summary())

## AIC=1141.989；BIC= 1153.138；兩個係數是顯著的

## 檢視模型的擬合殘差分佈
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()

## 預測未來24個資料，並輸出95%置信區間
pre, se, conf = arma_model.forecast(24, alpha=0.05)
## 整理資料
y_hat["arma_pre"] = pre
y_hat["arma_pre_lower"] = conf[:,0]
y_hat["arma_pre_upper"] = conf[:,1]
## 視覺化出預測結果
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)")
## 視覺化出置信區間
plt.fill_between(y_hat.index, y_hat["arma_pre_lower"], 
                 y_hat["arma_pre_upper"],color='k',alpha=.15,
                 label = "95%置信區間")
plt.legend()
plt.grid()
plt.title("ARMA(2,1)模型")
plt.show()
# 計算預測結果和真實值的誤差
print("ARMA模型預測的絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat["arma_pre"]))

ARMA模型預測的絕對值誤差: 147.26531763335154

針對ARMA模型自動選擇合適的引數

## 自動搜尋合適的引數
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模型引數
                      m=1,              # 序列的週期
                      seasonal=False,   # 沒有季節性趨勢
                      trace=True,error_action='ignore',  
                      suppress_warnings=True, stepwise=True)

print(model.summary())

## 使用ARMA(3,2)對測試集進行預測
pre, conf = model.predict(n_periods=24, alpha=0.05,
                          return_conf_int=True)
## 視覺化ARMA(3,2)的預測結果，整理資料
y_hat["arma_pre"] = pre
y_hat["arma_pre_lower"] = conf[:,0]
y_hat["arma_pre_upper"] = conf[:,1]
## 視覺化出預測結果
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)")
## 視覺化出置信區間
plt.fill_between(y_hat.index, y_hat["arma_pre_lower"], 
                 y_hat["arma_pre_upper"],color='k',alpha=.15,
                 label = "95%置信區間")
plt.legend()
plt.grid()
plt.title("ARMA(3,2)模型")
plt.show()

# 計算預測結果和真實值的誤差
print("ARMA模型預測的絕對值誤差:",
      mean_absolute_error(test["X1"],y_hat["arma_pre"]))

ARMA模型預測的絕對值誤差: 158.11464180972925

可以發現使用自動ARMA(3,2)模型的效果並沒有ARMA(2,1)的預測效果好

時序資料的異常值檢測

可以將突然增大或突然減小的資料無規律看作異常值

## 使用prophet檢測時間序列是否有異常值

## 從1991年2月到2005年5月，每週提供美國成品汽車汽油產品的時間序列（每天數千桶）

## 資料準備
data = pm.datasets.load_gasoline()
datadf = pd.DataFrame({"y":data})
datadf["ds"] = pd.date_range(start="1991-2",periods=len(data),freq="W")
## 視覺化時間序列的變化情況
datadf.plot(x = "ds",y = "y",style = "b-o",figsize=(14,7))
plt.grid()
plt.title("時間序列資料的波動情況")
plt.show()

## 對該資料建立一個時間序列模型
np.random.seed(1234)  ## 設定亂數種子
model = Prophet(growth="linear",daily_seasonality = False,
                weekly_seasonality=False,
                seasonality_mode = 'multiplicative',
                interval_width = 0.95,   ## 獲取95%的置信區間
                )
model = model.fit(datadf)     # 使用資料擬合模型
forecast = model.predict(datadf)  # 使用模型對資料進行預測
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

## 根據模型預測值的置信區間"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("異常值的數量為:",np.sum(outlier_index))
'''
異常值的數量為: 38
'''

## 視覺化異常值的結果
fig, ax = plt.subplots()
## 視覺化預測值
forecast.plot(x = "ds",y = "yhat",style = "b-",figsize=(14,7),
              label = "預測值",ax=ax)
## 視覺化出置信區間
ax.fill_between(forecast["ds"].values, forecast["yhat_lower"], 
                forecast["yhat_upper"],color='b',alpha=.2,
                label = "95%置信區間")
forecast.plot(kind = "scatter",x = "ds",y = "y",c = "k",
              s = 20,label = "原始資料",ax = ax)
## 視覺化出異常值的點
outlier_df.plot(x = "ds",y = "y",style = "rs",ax = ax,
                label = "異常值")
plt.legend(loc = 2)
plt.grid()
plt.title("時間序列異常值檢測結果")
plt.show()