無論我們是想預測金融市場的趨勢還是用電量，時間都是我們模型中必須考慮的一個重要因素。例如，預測一天中什麼時候會出現用電高峰是很有趣的，可以以此為依據調整電價或發電量。
輸入時間序列。時間序列只是按時間順序排列的一系列數據點。在時間序列中，時間往往是獨立變量，其目標通常是預測未來。
然而，在處理時間序列時，還有一些其他因素會發揮作用。
它是靜止的嗎？
有季節性嗎？
目標變量是否自相關？
在這篇文章中，我將介紹時間序列的不同特徵，以及我們如何對它們進行建模才能獲得準確的預測。
預測未來是困難的
自相關
通俗地說，自相關是觀測值之間的相似度，它是觀測值之間時間滯後的函數。
自相關示例
上面是一個自相關的例子。仔細觀察，你會發現第一個值和第 24 個值具有很高的自相關性。同樣，第 12 個值和第 36 個觀測值也高度相關。這意味著我們將在每 24 個時間單位中找到一個非常相似的值。
注意，這個圖看起來像正弦函數。這是季節性的徵兆，你可以通過在上面的圖中找到 24 小時的周期來找到它的價值。
季節性
季節性是指周期性波動。例如，白天的用電量高，晚上的用電量低，或者聖誕節期間的在線銷售額增加，節後銷售再次放緩。
季節性示例
如你所見，每天都有明顯的季節性。每天晚上，你都會看到一個高峰，最低點出現在每天的開始和結束。
記住，如果季節性是滿足正弦函數的，它也可以從自相關圖中推導出來。簡單地看一下周期，它給出了季節的長度。
平穩性
平穩性是時間序列的一個重要特徵。如果時間序列的統計性質不隨時間變化，則稱其為平穩的。換句話說，它有不變的均值和方差，協方差不隨時間變化。
平穩過程示例
再看上面的圖，我們看到上面的過程是平穩的，平均值和方差不會隨時間變化。
通常，股票價格不是一個平穩的過程，因為我們可能會看到一個增長的趨勢，或者，其波動性可能會隨著時間的推移而增加（這意味著方差正在變化）。
理想情況下，我們需要一個用於建模的固定時間序列。當然，不是所有的都是平穩的，但是我們可以通過做不同的變換，使它們保持平穩。
如何測試過程是否平穩
你可能已經注意到在上圖的標題「Dickey-Fuller」。這是我們用來確定時間序列是否穩定的統計測試。
在不討論 Dickey-Fuller 測試的技術特性的情況下，它檢測了單位根是否存在空假設。
如果是，則 P>0，並且過程不是平穩的。
否則，p=0，無效假設被拒絕，過程被認為是平穩的。
例如，下面的過程不是平穩的。請注意為什麼平均值不隨時間變化。
非平穩過程示例
時間序列建模
有很多方法可以模擬時間序列來進行預測。在此，我將介紹：
移動平均
指數平滑
ARIMA
移動平均
移動平均模型可能是最簡單的時間序列建模方法。這個模型簡單來說就是，下一個值是所有過去值的平均值。
雖然很簡單，但是這個模型的效果可能好到出乎意料，它代表了一個好的起點。
否則，移動平均值可用於識別數據中有趣的趨勢。我們可以定義一個窗口來應用移動平均模型來平滑時間序列，並突出不同的趨勢。
24 小時窗口上的移動平均值示例
在上面的圖中，我們將移動平均模型應用於一個 24 小時窗口。綠線平滑了時間序列，我們可以看到 24 小時內有 2 個峰值。
當然，窗口越長，趨勢就越平滑。下面是一個較小窗口上移動平均值的示例。
12 小時窗口上的移動平均值示例
指數平滑
指數平滑使用與移動平均相似的邏輯，但這次，對每個觀測值分配了不同的遞減權重。換言之，離現在的時間距離越遠，觀察結果的重要性就越低。
在數學上，指數平滑表示為：
指數平滑表達式
這裡，alpha 是一個平滑因子，它的值介於 0 和 1 之間。它決定了之前觀測值的權重下降的速度。
指數平滑示例
在上面的圖中，深藍色線表示時間序列的指數平滑，平滑係數為 0.3，而橙色線表示平滑係數為 0.05。
如你所見，平滑因子越小，時間序列就越平滑。這是有意義的，因為當平滑因子接近 0 時，我們接近移動平均模型。
雙指數平滑
當時間序列中存在趨勢時，使用雙指數平滑。在這種情況下，我們使用這種技術，它只是指數平滑的兩次遞歸使用。
數學公式為：
雙指數平滑表達式
這裡，beta 是趨勢平滑因子，它的值介於 0 和 1 之間。
下面，你可以看到 alpha 和 beta 的不同值如何影響時間序列的形狀。
雙指數平滑示例
三指數平滑
該方法通過添加季節平滑因子來擴展雙指數平滑。當然，如果你注意到時間序列中的季節性，這很有用。
在數學上，三指數平滑表示為：
三指數平滑表達式
其中 gamma 是季節平滑因子，L 是季節長度。
季節性差分自回歸滑動平均模型（SARIMA）
SARIMA 實際上是簡單模型的組合，可以生成一個複雜的模型，該模型可以模擬具有非平穩特性和季節性的時間序列。
首先，我們得到了自回歸模型 AR(p)。這基本上是時間序列對自身的回歸。在這裡，我們假設當前值依賴於它以前的值，並且有一定的滯後。它採用一個表示最大滯後的參數 p。為了找到它，我們查看了部分自相關圖，在此之後大部分滯後並不顯著。
在下面的例子中，p 的值是 4。
部分自相關圖示例
然後，我們添加移動平均模型 MA(q)。這需要一個參數 q，它代表自相關圖上那些滯後不顯著的最大滯後。
下圖中，q 為 4。
自相關圖示例
之後，我們添加整合順序 I(d)。參數 d 表示使序列平穩所需的差異數。
最後，我們添加最後一部分：季節性 S(P, D, Q, s)，其中 S 只是季節的長度。此外，這裡要求參數 P 和 Q 與 p 和 q 相同，但用於季節部分。最後，D 是季節整合的順序，表示從系列中刪除季節性所需的差異數量。
綜合起來，我們得到了 SARIMA(p, d, q)(P, D, Q, s) 模型。
要注意的是：在用 SARIMA 建模之前，我們必須對時間序列進行轉換，以消除季節性和任何非平穩行為。
這是一個很好的理論！讓我們在第一個項目中應用上面討論的技術。
我們將想辦法預測一家公司的股票價格。現在，預測股票價格幾乎是不可能的。然而，這仍然是一個有趣的練習，它將是一個很好的來實踐我們所學到知識的方法。
項目 1：股票價格預測
我們將利用 New Germany Fund（GF）的歷史股價來預測未來五個交易日的收盤價。
你可以在這裡獲取數據集和資料。
像往常一樣，我強烈推薦你動手編碼！啟動你的筆記本，我們開始吧！
你絕不會因為這個項目而發財
導入數據
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.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)
首先，我們導入一些庫，這些庫將在整個分析過程中都會用到。此外，我們用平均百分比誤差（MAPE）作為我們的誤差度量。
然後，我們導入數據集，排在前十的是：
數據集的前 10 個條目
正如你所看到的，我們有一些關於 New Germany Fund (GF) 不同股票的數據。此外，我們還有一個關於當天信息的數據，但我們只需要當天結束（EOD）時的股票信息。
數據清洗
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()
首先，我們刪除不需要的條目。
然後，我們刪除不需要的列，因為我們只想關注股票的收盤價。
如果預覽數據集，則應該看到的是：
清洗後的數據集
令人驚嘆！我們準備好進行探索性數據分析了！
探索性數據分析（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()
我們繪製數據集整個時間段的收盤價。
你將會得到：
New Germany Fund （GF）收盤價
很明顯，你看到的不是一個平穩的過程，很難判斷是否有某種季節性。
移動平均
讓我們使用移動平均模型來平滑我們的時間序列。為此，我們將使用一個輔助函數，該函數將在指定的時間窗口上運行移動平均模型，並繪製結果平滑曲線：
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天的時間窗口，我們得到：
上一個交易周的平滑曲線
如你所見，我們幾乎看不到趨勢，因為它太接近實際曲線。讓我們看看上個月和上個季度的平滑處理結果。
上個月（30 天前）的平滑曲線
按上一季度（90 天）平滑
現在更容易發現趨勢。請注意，30 天和 90 天的趨勢圖在末尾顯示一條向下的曲線。這意味著股票可能在接下來的幾天內會下跌。
指數平滑
現在，讓我們用指數平滑來看看它是否能獲得更好的趨勢。
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 作為平滑因子的值。當然你也可以嘗試其他值，看看結果如何。
指數平滑
如您所見，alpha 值 0.05 平滑了曲線，同時剔除了大部分向上和向下的趨勢。
現在，讓我們使用雙指數平滑。
雙指數平滑
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])
你將得到：
雙指數平滑
同樣，用不同的 α 和 β 組合進行實驗，以獲得更好的曲線。
建模
如前所述，我們必須將序列轉換為一個平穩的過程，以便對其進行建模。因此，讓我們應用 Dickey-Fuller 測試來看看它是否是一個平穩的過程：
def tsplot(y, lags=None, figsize=(12, 7), syle='bmh'):
if not isinstance(y, pd.Series):
y = pd.Series(y)
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)
ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p='.format(p_value))
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)
你將看到：
通過 DickeyFuller 測試，時間序列是非平穩的。另外，從自相關圖來看，我們發現它似乎沒有明顯的季節性。
因此，為了消除高度自相關並使過程穩定，讓我們取第一個差異（代碼塊中的第 23 行）。我們簡單地用一天的滯後時間減去時間序列，得到：
令人驚嘆的！我們的序列現在是平穩的，可以開始建模了！
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):
seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
except:
continue
aic = model.aic
#Save best model, AIC and parameters
if aicpan>
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]
seasonal_order=(P, D, Q, s)).fit(disp=-1)
print(best_model.summary())
現在，對於 SARIMA，我們首先定義一些參數值的範圍，以生成 p, q, d, P, Q, D, s 的所有可能組合的列表。
現在，在上面的代碼單元中，我們有 625 種不同的組合！我們將嘗試每種組合，並訓練 SARIMA，以便找到性能最佳的模型。這可能需要一些時間，具體多長時間取決於計算機的處理能力。
完成後，我們將輸出最佳模型的摘要，你將看到：
令人驚嘆！最後，我們預測未來五個交易日的收盤價，並評估模型的 MAPE。
在這種情況下，有一個 0.79% 的 MAPE，這是非常好的！
將預測價格與實際數據進行比較
# 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()
現在，為了將我們的預測與實際數據進行比較，我們從雅虎財務（YahooFinance）獲取財務數據並創建一個數據框架。
然後，我們繪出曲線，看看我們與實際收盤價的差距有多大：
預計值和實際收盤價比較
我們的預測似乎有點偏離。事實上，預測價格很平穩，這意味著我們的模型可能表現不佳。
當然，這不是因為我們的程序，而是因為預測股票價格基本上是不可能的。
從第一個項目開始，我們學習了在使用 SARIMA 建模之前平滑時間序列的整個過程。
現在，讓我們介紹一下 Facebook 的 Prophet。它是一個在 python 和 r 中都可用的預測工具。該工具幫助生成高質量的預測。
讓我們看看如何在第二個項目中使用它！
項目2-使用 Prophet 預測空氣品質
標題說明了一切：我們將使用 Prophet 來幫助我們預測空氣品質！
導入數據
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()
列印出前五行：
如你所見，數據集包含有關不同氣體濃度的信息。每天隔一個小時記錄一次。
如果更深入地研究數據集，會發現有許多值 -200 的實例。當然，負濃度是沒有意義的，所以我們需要在建模前清洗數據。
數據清洗與特徵工程
# 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()pan>
# 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)
在這裡，我們首先分析日期列，將其轉換為日期類型。
然後，我們把所有的測量值轉換成浮點數。
之後，我們用每天的平均值來匯總數據。
我們還有一些需要刪除的 NAN。
最後，我們按周匯總數據，這將提供一個更平滑的分析趨勢。
我們可以畫出每種化學物質濃度的趨勢。這裡，我們展示了 NOx。
NOx 濃度
氮氧化物是非常有害的，因為它們會形成煙霧和酸雨，同時也會形成細顆粒和臭氧。這些都會對健康產生不利影響，因此氮氧化物的濃度是空氣品質的一個關鍵特徵。
建模
# 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()
我們將只關注氮氧化物濃度。因此，我們刪除所有其他不相關的列。
然後，我們導入 Prophet。
Prophet 要求日期列命名為 ds，特徵列命名為 y，因此我們進行了適當的更改。
此時，我們的數據如下：
然後，我們定義一個訓練集。為此，我們將保留最後 30 個條目進行預測和驗證。
之後，我們簡單地初始化 Prophet，將模型與數據匹配，並進行預測！
你會看到：
這裡，yhat 代表預測值，yhat_lower 和 yhat_upper 分別代表預測值的下限和上限。
Prophet 讓你可以輕鬆繪製預測圖，我們得到：
NOx 濃度預測
如你所見，Prophet 只是用一條直線來預測未來的 NOx 濃度。
然後，我們檢查時間序列是否具有某些有趣的特性，例如季節性：
在這裡，Prophet 沒有發現季節性的趨勢。
通過計算模型的平均絕對百分誤差（MAPE）和平均絕對誤差（MAE）來評估模型的性能，我們發現 MAPE 為 13.86%，MAE 為 109.32，這還不錯！記住，我們根本沒有對模型進行微調。
最後，我們只需繪製預測的上限和下限：
每周 NOx 平均濃度預測
恭喜你達到目的！這篇文章很長，但內容豐富。你學會了如何強有力地分析和建模時間序列，並將你的知識應用到兩個不同的項目中。我希望你覺得這篇文章有用。
雷鋒網雷鋒網雷鋒網