時序分析與預(yù)測完全指南

本文作者： skura

2019-08-15 18:45

導(dǎo)語：了解移動平均、指數(shù)平滑、平穩(wěn)性、自相關(guān)、時間序列，將這些技術(shù)應(yīng)用于兩個項目中。

時序分析與預(yù)測完全指南

無論我們是想預(yù)測金融市場的趨勢還是用電量，時間都是我們模型中必須考慮的一個重要因素。例如，預(yù)測一天中什么時候會出現(xiàn)用電高峰是很有趣的，可以以此為依據(jù)調(diào)整電價或發(fā)電量。

輸入時間序列。時間序列只是按時間順序排列的一系列數(shù)據(jù)點(diǎn)。在時間序列中，時間往往是獨(dú)立變量，其目標(biāo)通常是預(yù)測未來。

然而，在處理時間序列時，還有一些其他因素會發(fā)揮作用。

它是靜止的嗎？

有季節(jié)性嗎？

目標(biāo)變量是否自相關(guān)？

在這篇文章中，我將介紹時間序列的不同特征，以及我們?nèi)绾螌λ鼈冞M(jìn)行建模才能獲得準(zhǔn)確的預(yù)測。

時序分析與預(yù)測完全指南

預(yù)測未來是困難的

自相關(guān)

通俗地說，自相關(guān)是觀測值之間的相似度，它是觀測值之間時間滯后的函數(shù)。

時序分析與預(yù)測完全指南

自相關(guān)示例

上面是一個自相關(guān)的例子。仔細(xì)觀察，你會發(fā)現(xiàn)第一個值和第 24 個值具有很高的自相關(guān)性。同樣，第 12 個值和第 36 個觀測值也高度相關(guān)。這意味著我們將在每 24 個時間單位中找到一個非常相似的值。

注意，這個圖看起來像正弦函數(shù)。這是季節(jié)性的征兆，你可以通過在上面的圖中找到 24 小時的周期來找到它的價值。

季節(jié)性

季節(jié)性是指周期性波動。例如，白天的用電量高，晚上的用電量低，或者圣誕節(jié)期間的在線銷售額增加，節(jié)后銷售再次放緩。

時序分析與預(yù)測完全指南

季節(jié)性示例

如你所見，每天都有明顯的季節(jié)性。每天晚上，你都會看到一個高峰，最低點(diǎn)出現(xiàn)在每天的開始和結(jié)束。
記住，如果季節(jié)性是滿足正弦函數(shù)的，它也可以從自相關(guān)圖中推導(dǎo)出來。簡單地看一下周期，它給出了季節(jié)的長度。

平穩(wěn)性

平穩(wěn)性是時間序列的一個重要特征。如果時間序列的統(tǒng)計性質(zhì)不隨時間變化，則稱其為平穩(wěn)的。換句話說，它有不變的均值和方差，協(xié)方差不隨時間變化。

時序分析與預(yù)測完全指南

平穩(wěn)過程示例

再看上面的圖，我們看到上面的過程是平穩(wěn)的，平均值和方差不會隨時間變化。

通常，股票價格不是一個平穩(wěn)的過程，因為我們可能會看到一個增長的趨勢，或者，其波動性可能會隨著時間的推移而增加（這意味著方差正在變化）。

理想情況下，我們需要一個用于建模的固定時間序列。當(dāng)然，不是所有的都是平穩(wěn)的，但是我們可以通過做不同的變換，使它們保持平穩(wěn)。

如何測試過程是否平穩(wěn)

你可能已經(jīng)注意到在上圖的標(biāo)題「Dickey-Fuller」。這是我們用來確定時間序列是否穩(wěn)定的統(tǒng)計測試。

在不討論 Dickey-Fuller 測試的技術(shù)特性的情況下，它檢測了單位根是否存在空假設(shè)。

如果是，則 P>0，并且過程不是平穩(wěn)的。

否則，p=0，無效假設(shè)被拒絕，過程被認(rèn)為是平穩(wěn)的。

例如，下面的過程不是平穩(wěn)的。請注意為什么平均值不隨時間變化。

時序分析與預(yù)測完全指南

非平穩(wěn)過程示例

時間序列建模

有很多方法可以模擬時間序列來進(jìn)行預(yù)測。在此，我將介紹：

移動平均
指數(shù)平滑
ARIMA

移動平均

移動平均模型可能是最簡單的時間序列建模方法。這個模型簡單來說就是，下一個值是所有過去值的平均值。

雖然很簡單，但是這個模型的效果可能好到出乎意料，它代表了一個好的起點(diǎn)。

否則，移動平均值可用于識別數(shù)據(jù)中有趣的趨勢。我們可以定義一個窗口來應(yīng)用移動平均模型來平滑時間序列，并突出不同的趨勢。

時序分析與預(yù)測完全指南

24 小時窗口上的移動平均值示例

在上面的圖中，我們將移動平均模型應(yīng)用于一個 24 小時窗口。綠線平滑了時間序列，我們可以看到 24 小時內(nèi)有 2 個峰值。

當(dāng)然，窗口越長，趨勢就越平滑。下面是一個較小窗口上移動平均值的示例。

時序分析與預(yù)測完全指南

12 小時窗口上的移動平均值示例

指數(shù)平滑

指數(shù)平滑使用與移動平均相似的邏輯，但這次，對每個觀測值分配了不同的遞減權(quán)重。換言之，離現(xiàn)在的時間距離越遠(yuǎn)，觀察結(jié)果的重要性就越低。

在數(shù)學(xué)上，指數(shù)平滑表示為：

時序分析與預(yù)測完全指南

指數(shù)平滑表達(dá)式

這里，alpha 是一個平滑因子，它的值介于 0 和 1 之間。它決定了之前觀測值的權(quán)重下降的速度。

時序分析與預(yù)測完全指南

指數(shù)平滑示例

在上面的圖中，深藍(lán)色線表示時間序列的指數(shù)平滑，平滑系數(shù)為 0.3，而橙色線表示平滑系數(shù)為 0.05。

如你所見，平滑因子越小，時間序列就越平滑。這是有意義的，因為當(dāng)平滑因子接近 0 時，我們接近移動平均模型。

雙指數(shù)平滑

當(dāng)時間序列中存在趨勢時，使用雙指數(shù)平滑。在這種情況下，我們使用這種技術(shù)，它只是指數(shù)平滑的兩次遞歸使用。

數(shù)學(xué)公式為：

時序分析與預(yù)測完全指南

雙指數(shù)平滑表達(dá)式

這里，beta 是趨勢平滑因子，它的值介于 0 和 1 之間。

下面，你可以看到 alpha 和 beta 的不同值如何影響時間序列的形狀。

時序分析與預(yù)測完全指南

雙指數(shù)平滑示例

三指數(shù)平滑

該方法通過添加季節(jié)平滑因子來擴(kuò)展雙指數(shù)平滑。當(dāng)然，如果你注意到時間序列中的季節(jié)性，這很有用。

在數(shù)學(xué)上，三指數(shù)平滑表示為：

時序分析與預(yù)測完全指南

三指數(shù)平滑表達(dá)式

其中 gamma 是季節(jié)平滑因子，L 是季節(jié)長度。

季節(jié)性差分自回歸滑動平均模型（SARIMA）

SARIMA 實際上是簡單模型的組合，可以生成一個復(fù)雜的模型，該模型可以模擬具有非平穩(wěn)特性和季節(jié)性的時間序列。

首先，我們得到了自回歸模型 AR(p)。這基本上是時間序列對自身的回歸。在這里，我們假設(shè)當(dāng)前值依賴于它以前的值，并且有一定的滯后。它采用一個表示最大滯后的參數(shù) p。為了找到它，我們查看了部分自相關(guān)圖，在此之后大部分滯后并不顯著。

在下面的例子中，p 的值是 4。

時序分析與預(yù)測完全指南

部分自相關(guān)圖示例

然后，我們添加移動平均模型 MA(q)。這需要一個參數(shù) q，它代表自相關(guān)圖上那些滯后不顯著的最大滯后。

下圖中，q 為 4。

時序分析與預(yù)測完全指南

自相關(guān)圖示例

之后，我們添加整合順序 I(d)。參數(shù) d 表示使序列平穩(wěn)所需的差異數(shù)。

最后，我們添加最后一部分：季節(jié)性 S(P, D, Q, s)，其中 S 只是季節(jié)的長度。此外，這里要求參數(shù) P 和 Q 與 p 和 q 相同，但用于季節(jié)部分。最后，D 是季節(jié)整合的順序，表示從系列中刪除季節(jié)性所需的差異數(shù)量。

綜合起來，我們得到了 SARIMA(p, d, q)(P, D, Q, s) 模型。

要注意的是：在用 SARIMA 建模之前，我們必須對時間序列進(jìn)行轉(zhuǎn)換，以消除季節(jié)性和任何非平穩(wěn)行為。

這是一個很好的理論！讓我們在第一個項目中應(yīng)用上面討論的技術(shù)。

我們將想辦法預(yù)測一家公司的股票價格。現(xiàn)在，預(yù)測股票價格幾乎是不可能的。然而，這仍然是一個有趣的練習(xí)，它將是一個很好的來實踐我們所學(xué)到知識的方法。

項目 1：股票價格預(yù)測

我們將利用 New Germany Fund（GF）的歷史股價來預(yù)測未來五個交易日的收盤價。

你可以在這里獲取數(shù)據(jù)集和資料。

像往常一樣，我強(qiáng)烈推薦你動手編碼！啟動你的筆記本，我們開始吧！

時序分析與預(yù)測完全指南

你絕不會因為這個項目而發(fā)財

導(dǎo)入數(shù)據(jù)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error

from scipy.optimize import minimize
import statsmodels.tsa.api as smt
import statsmodels.api as sm

from tqdm import tqdm_notebook

from itertools import product

def mean_absolute_percentage_error(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

import warnings
warnings.filterwarnings('ignore')

%matplotlib inline

DATAPATH = 'data/stock_prices_sample.csv'

data = pd.read_csv(DATAPATH, index_col=['DATE'], parse_dates=['DATE'])
data.head(10)

首先，我們導(dǎo)入一些庫，這些庫將在整個分析過程中都會用到。此外，我們用平均百分比誤差（MAPE）作為我們的誤差度量。

然后，我們導(dǎo)入數(shù)據(jù)集，排在前十的是：

時序分析與預(yù)測完全指南

數(shù)據(jù)集的前 10 個條目

正如你所看到的，我們有一些關(guān)于 New Germany Fund (GF) 不同股票的數(shù)據(jù)。此外，我們還有一個關(guān)于當(dāng)天信息的數(shù)據(jù)，但我們只需要當(dāng)天結(jié)束（EOD）時的股票信息。

數(shù)據(jù)清洗

data = data[data.TICKER != 'GEF']
data = data[data.TYPE != 'Intraday']

drop_cols = ['SPLIT_RATIO', 'EX_DIVIDEND', 'ADJ_FACTOR', 'ADJ_VOLUME', 'ADJ_CLOSE', 'ADJ_LOW', 'ADJ_HIGH', 'ADJ_OPEN', 'VOLUME', 'FREQUENCY', 'TYPE', 'FIGI']

data.drop(drop_cols, axis=1, inplace=True)

data.head()

首先，我們刪除不需要的條目。

然后，我們刪除不需要的列，因為我們只想關(guān)注股票的收盤價。

如果預(yù)覽數(shù)據(jù)集，則應(yīng)該看到的是：

時序分析與預(yù)測完全指南

清洗后的數(shù)據(jù)集

令人驚嘆！我們準(zhǔn)備好進(jìn)行探索性數(shù)據(jù)分析了！

探索性數(shù)據(jù)分析（EDA）

# Plot closing price

plt.figure(figsize=(17, 8))
plt.plot(data.CLOSE)
plt.title('Closing price of New Germany Fund Inc (GF)')
plt.ylabel('Closing price ($)')
plt.xlabel('Trading day')
plt.grid(False)
plt.show()

我們繪制數(shù)據(jù)集整個時間段的收盤價。

你將會得到：

時序分析與預(yù)測完全指南

New Germany Fund （GF）收盤價

很明顯，你看到的不是一個平穩(wěn)的過程，很難判斷是否有某種季節(jié)性。

移動平均

讓我們使用移動平均模型來平滑我們的時間序列。為此，我們將使用一個輔助函數(shù)，該函數(shù)將在指定的時間窗口上運(yùn)行移動平均模型，并繪制結(jié)果平滑曲線：

def plot_moving_average(series, window, plot_intervals=False, scale=1.96):

    rolling_mean = series.rolling(window=window).mean()

    plt.figure(figsize=(17,8))
    plt.title('Moving average\n window size = {}'.format(window))
    plt.plot(rolling_mean, 'g', label='Rolling mean trend')

    #Plot confidence intervals for smoothed values
    if plot_intervals:
        mae = mean_absolute_error(series[window:], rolling_mean[window:])
        deviation = np.std(series[window:] - rolling_mean[window:])
       lower_bound = rolling_mean - (mae + scale * deviation)
       upper_bound = rolling_mean + (mae + scale * deviation)
       plt.plot(upper_bound, 'r--', label='Upper bound / Lower bound')
       plt.plot(lower_bound, 'r--')

    plt.plot(series[window:], label='Actual values')
    plt.legend(loc='best')
    plt.grid(True)

#Smooth by the previous 5 days (by week)
plot_moving_average(data.CLOSE, 5)

#Smooth by the previous month (30 days)
plot_moving_average(data.CLOSE, 30)

#Smooth by previous quarter (90 days)
plot_moving_average(data.CLOSE, 90, plot_intervals=True)

使用5天的時間窗口，我們得到：

時序分析與預(yù)測完全指南

上一個交易周的平滑曲線

如你所見，我們幾乎看不到趨勢，因為它太接近實際曲線。讓我們看看上個月和上個季度的平滑處理結(jié)果。

時序分析與預(yù)測完全指南

上個月（30 天前）的平滑曲線

時序分析與預(yù)測完全指南

按上一季度（90 天）平滑

現(xiàn)在更容易發(fā)現(xiàn)趨勢。請注意，30 天和 90 天的趨勢圖在末尾顯示一條向下的曲線。這意味著股票可能在接下來的幾天內(nèi)會下跌。

指數(shù)平滑

現(xiàn)在，讓我們用指數(shù)平滑來看看它是否能獲得更好的趨勢。

def exponential_smoothing(series, alpha):

result = [series[0]] # first value is same as series
for n in range(1, len(series)):
result.append(alpha * series[n] + (1 - alpha) * result[n-1])
return result

def plot_exponential_smoothing(series, alphas):

plt.figure(figsize=(17, 8))
for alpha in alphas:
plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))
plt.plot(series.values, "c", label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Exponential Smoothing")
plt.grid(True);

plot_exponential_smoothing(data.CLOSE, [0.05, 0.3])

這里，我們使用 0.05 和 0.3 作為平滑因子的值。當(dāng)然你也可以嘗試其他值，看看結(jié)果如何。

時序分析與預(yù)測完全指南

指數(shù)平滑

如您所見，alpha 值 0.05 平滑了曲線，同時剔除了大部分向上和向下的趨勢。

現(xiàn)在，讓我們使用雙指數(shù)平滑。

雙指數(shù)平滑

def double_exponential_smoothing(series, alpha, beta):

result = [series[0]]
for n in range(1, len(series)+1):
if n == 1:
level, trend = series[0], series[1] - series[0]
if n >= len(series): # forecasting
value = result[-1]
else:
value = series[n]
last_level, level = level, alpha * value + (1 - alpha) * (level + trend)
trend = beta * (level - last_level) + (1 - beta) * trend
result.append(level + trend)
return result

def plot_double_exponential_smoothing(series, alphas, betas):

plt.figure(figsize=(17, 8))
for alpha in alphas:
for beta in betas:
plt.plot(double_exponential_smoothing(series, alpha, beta), label="Alpha {}, beta {}".format(alpha, beta))
plt.plot(series.values, label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Double Exponential Smoothing")
plt.grid(True)

plot_double_exponential_smoothing(data.CLOSE, alphas=[0.9, 0.02], betas=[0.9, 0.02])

你將得到：

時序分析與預(yù)測完全指南

雙指數(shù)平滑

同樣，用不同的 α 和 β 組合進(jìn)行實驗，以獲得更好的曲線。

建模

如前所述，我們必須將序列轉(zhuǎn)換為一個平穩(wěn)的過程，以便對其進(jìn)行建模。因此，讓我們應(yīng)用 Dickey-Fuller 測試來看看它是否是一個平穩(wěn)的過程：

def tsplot(y, lags=None, figsize=(12, 7), syle='bmh'):

    if not isinstance(y, pd.Series):
        y = pd.Series(y)

   with plt.style.context(style='bmh'):
       fig = plt.figure(figsize=figsize)
       layout = (2,2)
       ts_ax = plt.subplot2grid(layout, (0,0), colspan=2)
       acf_ax = plt.subplot2grid(layout, (1,0))
       pacf_ax = plt.subplot2grid(layout, (1,1))

       y.plot(ax=ts_ax)
       p_value = sm.tsa.stattools.adfuller(y)[1]
       ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))
       smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
       smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
       plt.tight_layout()

tsplot(data.CLOSE, lags=30)

# Take the first difference to remove to make the process stationary
data_diff = data.CLOSE - data.CLOSE.shift(1)

tsplot(data_diff[1:], lags=30)

你將看到：

時序分析與預(yù)測完全指南

通過 DickeyFuller 測試，時間序列是非平穩(wěn)的。另外，從自相關(guān)圖來看，我們發(fā)現(xiàn)它似乎沒有明顯的季節(jié)性。

因此，為了消除高度自相關(guān)并使過程穩(wěn)定，讓我們?nèi)〉谝粋€差異（代碼塊中的第 23 行）。我們簡單地用一天的滯后時間減去時間序列，得到：

時序分析與預(yù)測完全指南

令人驚嘆的！我們的序列現(xiàn)在是平穩(wěn)的，可以開始建模了！

SARIMA

#Set initial values and some bounds
ps = range(0, 5)
d = 1
qs = range(0, 5)
Ps = range(0, 5)
D = 1
Qs = range(0, 5)
s = 5

#Create a list with all possible combinations of parameters
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)

# Train many SARIMA models to find the best set of parameters
def optimize_SARIMA(parameters_list, d, D, s):
"""
Return dataframe with parameters and corresponding AIC

parameters_list - list with (p, q, P, Q) tuples
d - integration order
D - seasonal integration order
s - length of season
"""

results = []
best_aic = float('inf')

for param in tqdm_notebook(parameters_list):
try: model = sm.tsa.statespace.SARIMAX(data.CLOSE, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
except:
continue

aic = model.aic

#Save best model, AIC and parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])

result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
#Sort in ascending order, lower AIC is better
result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)

return result_table

result_table = optimize_SARIMA(parameters_list, d, D, s)

#Set parameters that give the lowest AIC (Akaike Information Criteria)
p, q, P, Q = result_table.parameters[0]

best_model = sm.tsa.statespace.SARIMAX(data.CLOSE, order=(p, d, q),
seasonal_order=(P, D, Q, s)).fit(disp=-1)

print(best_model.summary())

現(xiàn)在，對于 SARIMA，我們首先定義一些參數(shù)值的范圍，以生成 p, q, d, P, Q, D, s 的所有可能組合的列表。

現(xiàn)在，在上面的代碼單元中，我們有 625 種不同的組合！我們將嘗試每種組合，并訓(xùn)練 SARIMA，以便找到性能最佳的模型。這可能需要一些時間，具體多長時間取決于計算機(jī)的處理能力。

完成后，我們將輸出最佳模型的摘要，你將看到：

時序分析與預(yù)測完全指南

令人驚嘆！最后，我們預(yù)測未來五個交易日的收盤價，并評估模型的 MAPE。

在這種情況下，有一個 0.79% 的 MAPE，這是非常好的！

將預(yù)測價格與實際數(shù)據(jù)進(jìn)行比較

# Make a dataframe containing actual and predicted prices
comparison = pd.DataFrame({'actual': [18.93, 19.23, 19.08, 19.17, 19.11, 19.12],
'predicted': [18.96, 18.97, 18.96, 18.92, 18.94, 18.92]},
index = pd.date_range(start='2018-06-05', periods=6,))

#Plot predicted vs actual price

plt.figure(figsize=(17, 8))
plt.plot(comparison.actual)
plt.plot(comparison.predicted)
plt.title('Predicted closing price of New Germany Fund Inc (GF)')
plt.ylabel('Closing price ($)')
plt.xlabel('Trading day')
plt.legend(loc='best')
plt.grid(False)
plt.show()

現(xiàn)在，為了將我們的預(yù)測與實際數(shù)據(jù)進(jìn)行比較，我們從雅虎財務(wù)（YahooFinance）獲取財務(wù)數(shù)據(jù)并創(chuàng)建一個數(shù)據(jù)框架。

然后，我們繪出曲線，看看我們與實際收盤價的差距有多大：

時序分析與預(yù)測完全指南

預(yù)計值和實際收盤價比較

我們的預(yù)測似乎有點(diǎn)偏離。事實上，預(yù)測價格很平穩(wěn)，這意味著我們的模型可能表現(xiàn)不佳。

當(dāng)然，這不是因為我們的程序，而是因為預(yù)測股票價格基本上是不可能的。

從第一個項目開始，我們學(xué)習(xí)了在使用 SARIMA 建模之前平滑時間序列的整個過程。

現(xiàn)在，讓我們介紹一下 Facebook 的 Prophet。它是一個在 python 和 r 中都可用的預(yù)測工具。該工具幫助生成高質(zhì)量的預(yù)測。

讓我們看看如何在第二個項目中使用它！

項目2-使用 Prophet 預(yù)測空氣質(zhì)量

標(biāo)題說明了一切：我們將使用 Prophet 來幫助我們預(yù)測空氣質(zhì)量！

時序分析與預(yù)測完全指南

導(dǎo)入數(shù)據(jù)

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt

%matplotlib inline

DATAPATH = 'data/AirQualityUCI.csv'

data = pd.read_csv(DATAPATH, sep=';')
data.head()

打印出前五行：

時序分析與預(yù)測完全指南

如你所見，數(shù)據(jù)集包含有關(guān)不同氣體濃度的信息。每天隔一個小時記錄一次。

如果更深入地研究數(shù)據(jù)集，會發(fā)現(xiàn)有許多值 -200 的實例。當(dāng)然，負(fù)濃度是沒有意義的，所以我們需要在建模前清洗數(shù)據(jù)。

數(shù)據(jù)清洗與特征工程

# Make dates actual dates
data['Date'] = pd.to_datetime(data['Date'])

# Convert measurements to floats
for col in data.iloc[:,2:].columns:
if data[col].dtypes == object:
data[col] = data[col].str.replace(',', '.').astype('float')

# Compute the average considering only the positive values
def positive_average(num):
return num[num > -200].mean()

# Aggregate data
daily_data = data.drop('Time', axis=1).groupby('Date').apply(positive_average)

# Drop columns with more than 8 NaN
daily_data = daily_data.iloc[:,(daily_data.isna().sum() <= 8).values]

# Remove rows containing NaN values
daily_data = daily_data.dropna()

# Aggregate data by week
weekly_data = daily_data.resample('W').mean()

# Plot the weekly concentration of each gas
def plot_data(col):
plt.figure(figsize=(17, 8))
plt.plot(weekly_data[col])
plt.xlabel('Time')
plt.ylabel(col)
plt.grid(False)
plt.show()

for col in weekly_data.columns:
plot_data(col)

在這里，我們首先分析日期列，將其轉(zhuǎn)換為日期類型。

然后，我們把所有的測量值轉(zhuǎn)換成浮點(diǎn)數(shù)。

之后，我們用每天的平均值來匯總數(shù)據(jù)。

我們還有一些需要刪除的 NAN。

最后，我們按周匯總數(shù)據(jù)，這將提供一個更平滑的分析趨勢。

我們可以畫出每種化學(xué)物質(zhì)濃度的趨勢。這里，我們展示了 NOx。

時序分析與預(yù)測完全指南

NOx 濃度

氮氧化物是非常有害的，因為它們會形成煙霧和酸雨，同時也會形成細(xì)顆粒和臭氧。這些都會對健康產(chǎn)生不利影響，因此氮氧化物的濃度是空氣質(zhì)量的一個關(guān)鍵特征。

建模

# Drop irrelevant columns
cols_to_drop = ['PT08.S1(CO)', 'C6H6(GT)', 'PT08.S2(NMHC)', 'PT08.S4(NO2)', 'PT08.S5(O3)', 'T', 'RH', 'AH']

weekly_data = weekly_data.drop(cols_to_drop, axis=1)

# Import Prophet
from fbprophet import Prophet
import logging

logging.getLogger().setLevel(logging.ERROR)

# Change the column names according to Prophet's guidelines
df = weekly_data.reset_index()
df.columns = ['ds', 'y']
df.head()

# Split into a train/test set
prediction_size = 30
train_df = df[:-prediction_size]

# Initialize and train a model
m = Prophet()
m.fit(train_df)

# Make predictions
future = m.make_future_dataframe(periods=prediction_size)
forecast = m.predict(future)
forecast.head()

# Plot forecast
m.plot(forecast)

# Plot forecast's components
m.plot_components(forecast)

# Evaluate the model
def make_comparison_dataframe(historical, forecast):
return forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical.set_index('ds'))

cmp_df = make_comparison_dataframe(df, forecast)
cmp_df.head()

def calculate_forecast_errors(df, prediction_size):

df = df.copy()

df['e'] = df['y'] - df['yhat']
df['p'] = 100 * df['e'] / df['y']

predicted_part = df[-prediction_size:]

error_mean = lambda error_name: np.mean(np.abs(predicted_part[error_name]))

return {'MAPE': error_mean('p'), 'MAE': error_mean('e')}

for err_name, err_value in calculate_forecast_errors(cmp_df, prediction_size).items():
print(err_name, err_value)

# Plot forecast with upper and lower bounds
plt.figure(figsize=(17, 8))
plt.plot(cmp_df['yhat'])
plt.plot(cmp_df['yhat_lower'])
plt.plot(cmp_df['yhat_upper'])
plt.plot(cmp_df['y'])
plt.xlabel('Time')
plt.ylabel('Average Weekly NOx Concentration')
plt.grid(False)
plt.show()