在真實的資料中，許多現象都具有時間上的先後順序，而非彼此獨立的觀測值；例如股票的每日收盤價、氣象站每小時的溫度紀錄、工廠感測器每秒回傳的振動數值，以及電商平台每月的銷售量等。這類依照時間順序排列的資料，我們稱之為時間序列(time series)。時間序列資料的核心挑戰在於，相鄰的觀測值之間往往不是獨立的，而是存在某種程度的相關性，例如昨天的股價會影響今天的股價，上個月的銷售量也可能反映出這個月的需求趨勢。因此，無論是預測未來的數值、偵測異常的時間點，還是理解資料背後的週期規律，都需要專門針對這種「時間相依性」設計的分析方法。本篇將從時間序列的基本組成與特性出發，逐步介紹常用的分析與預測模型，以及實務上的注意事項。

在分析或預測時間序列之前，通常可以先嘗試將其分解為趨勢、季節性與殘差三個成分的和，以便分別理解各部分的行為。以下範例，會以太陽黑子的年活動量資料為例，示範如何進行此分解，並比較 statsmodels 函式庫與自行以移動平均計算的結果：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.datasets import sunspots
from statsmodels.tsa.seasonal import seasonal_decompose

data = sunspots.load_pandas().data
y = data["SUNACTIVITY"].to_numpy()
years = data["YEAR"].to_numpy()
period = 11

result = seasonal_decompose(y, model="additive", period=period)

pad = period // 2
trend_manual = np.array([
	np.mean(y[i - pad: i + pad + 1]) for i in range(pad, len(y) - pad)
])
years_valid = years[pad: len(y) - pad]
detrended = y[pad: len(y) - pad] - trend_manual
seasonal_manual = np.zeros(len(detrended))
for phase in range(period):
	idx = np.arange(phase, len(detrended), period)
	seasonal_manual[idx] = np.mean(detrended[idx])
residual_manual = detrended - seasonal_manual

fig, axes = plt.subplots(4, 2)
titles_left = ["Original (statsmodels)", "Trend", "Seasonal", "Residual"]
titles_right = ["Original (manual)", "Trend", "Seasonal", "Residual"]

data_left = [y, result.trend, result.seasonal, result.resid]
data_right = [y, trend_manual, seasonal_manual, residual_manual]
x_left = [years, years, years, years]
x_right = [years, years_valid, years_valid, years_valid]

for i in range(4):
	axes[i, 0].plot(x_left[i], data_left[i])
	axes[i, 1].plot(x_right[i], data_right[i])
	axes[i, 0].set_title(titles_left[i])
	axes[i, 1].set_title(titles_right[i])

plt.show()

在上述範例中：

• 將週期設為 11，是因為我們已經知道太陽黑子的變化週期大約是十一年。如果你對資料的週期大小還沒有足夠的把握，有可能需要多加嘗試，找出適當的值。
• 「model="additive"」是指，趨勢、季節性與殘差三個成分之間的關係，是相加後等於原訊號。
• 自行計算的步驟，是先用一個長度等於週期的 window，做移動平均的計算來當作趨勢；將訊號減去趨勢，再對每個週期內的同一位置取平均，就可以計算週期變化；而剩餘的部分就是殘差。其中，以移動平均計算趨勢時，每個時間點的趨勢值，是以該點為中心，前後各取若干點來平均。以太陽黑子為例，每個點的趨勢值，是用它前面 5 個點、它自己、以及後面 5 個點，共 11 個點來平均。

前述的加法模型，假設季節性振幅不會隨趨勢改變。然而在某些資料中，振幅會隨著整體數值的增大而放大，例如航空旅客人數在 1949 年至 1960 年間的資料，涵蓋了航空業快速成長的年代，讓旅客量不僅有明顯的上升趨勢，且季節性的波動幅度也隨之擴大。此時，改用乘法模型會更為合適，即將資料分解為趨勢、季節性與殘差三個成分的積。以下會以航空旅客人數資料，示範如何進行乘法分解，並同樣比較 statsmodels 函式庫與自行計算的結果：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.seasonal import seasonal_decompose

data = get_rdataset("AirPassengers").data
y = data["value"].to_numpy().astype(float)
period = 12

result = seasonal_decompose(y, model="multiplicative", period=period)

pad = period // 2
trend_manual = np.array([
	np.mean(y[i - pad: i + pad + 1]) for i in range(pad, len(y) - pad)
])
x_valid = np.arange(pad, len(y) - pad)
detrended = y[pad: len(y) - pad] / trend_manual
seasonal_manual = np.zeros(len(detrended))
for phase in range(period):
	idx = np.arange(phase, len(detrended), period)
	seasonal_manual[idx] = np.mean(detrended[idx])
residual_manual = detrended / seasonal_manual

fig, axes = plt.subplots(4, 2)
titles_left = ["Original (statsmodels)", "Trend", "Seasonal", "Residual"]
titles_right = ["Original (manual)", "Trend", "Seasonal", "Residual"]

x = np.arange(len(y))
data_left = [y, result.trend, result.seasonal, result.resid]
data_right = [y, trend_manual, seasonal_manual, residual_manual]
x_left = [x, x, x, x]
x_right = [x, x_valid, x_valid, x_valid]

for i in range(4):
	axes[i, 0].plot(x_left[i],  data_left[i])
	axes[i, 1].plot(x_right[i], data_right[i])
	axes[i, 0].set_title(titles_left[i])
	axes[i, 1].set_title(titles_right[i])

plt.show()

在上述範例中：

• 因為這筆資料的週期，和我們採用的分析的方法都不同，所以跟第一個範例相比，需要將週期從 11 改為 12，並將分解模型從加法改為乘法。
• 承上，若需要自行計算時，可以簡單地把相減改成相除。但若你的資料中有特殊狀況（例如某天完全沒客人上門），則有可能遇到除以 0 的問題，此時可以在分母加上一個極小值，或者在分母為 0 時填 NaN 或其他特定值。
• X 軸的意義改為月份的整數索引。如果你希望顯示年月，則可能需要搭配 pandas 做額外處理，會比較方便。

如果你要進行時間序列資料的預測，則有很多種模型可以選擇。以下先介紹傳統統計模型中的 ARIMA (AutoregRessive Integrated Moving Average)。ARIMA 是一個經典且廣泛使用的預測模型，適合單變數、具有趨勢但不含明顯季節性的時間序列；若資料同時含有季節性，則可進一步擴展為 SARIMA。以下範例會以航空旅客人數資料為例，示範資料切割、滾動預測，以及常見評估指標的計算：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.datasets import get_rdataset
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, mean_squared_error

data = get_rdataset("AirPassengers").data
y = np.log(data["value"].to_numpy().astype(float))
n = len(y)

test_size = 24
train = y[:n-test_size]
test = y[n-test_size :]

predictions = []
history = list(train)
for i in range(test_size):
	model = ARIMA(history, order=(2, 1, 2))
	result = model.fit()
	pred = result.forecast(steps=1)[0]
	predictions.append(pred)
	history.append(test[i])

predictions = np.array(predictions)

y_true = np.exp(test)
y_pred = np.exp(predictions)

mae = mean_absolute_error(y_true, y_pred)
rmse = mean_squared_error(y_true, y_pred) ** 0.5
mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
print(f"MAE:  {mae:.2f}")
print(f"RMSE: {rmse:.2f}")
print(f"MAPE: {mape:.2f}%")

x_all = np.arange(n)
x_train = np.arange(n - test_size)
x_test = np.arange(n - test_size, n)

plt.plot(x_train, np.exp(train), label="Train")
plt.plot(x_test, y_true, label="Test (ground truth)")
plt.plot(x_test, y_pred, label="Prediction")
plt.legend()
plt.show()

在上述範例中：

• 取對數計算是因為 AirPassengers 資料序列的結構比較適合乘法結構，因此透過取對數轉換為加法，讓 ARIMA 較容易處理；唯評估與繪圖時，仍需透過指數函數還原。
• 由於我們不應該在模型訓練的階段，讓模型「看見未來」，因此資料切分時，必須按照時間來切。範例中只切出了訓練與測試兩個子集，你若有需求，也可以切成訓練、驗證、測試三個子集。
• ARIMA 模型的三個參數，通常以 p、d、q 表示，說明如下：
  • p：AR (AutoregRessive)的階數，代表用過去的數值來預測當前值。
  • d：差分次數，代表對資料去除趨勢，使資料趨於平穩。
  • q：MA (Moving Average)的階數，代表用過去的預測誤差，來修正當前預測。
  範例中的三個參數為依照經驗選擇。而在實務上，可以用 ACF (AutoCorrelation Function) 以及 AIC/BIC 等統計指標來輔助選擇。
• 範例中這種每次預測下一個資料點的方式，稱為滾動預測。而範例中展示的情境是，每次僅讓模型預測下一個資料點（以 AirPassengers 資料集的粒度來說，是預測下個月的旅客數）；而當要預測下下一個資料點時，此情境已允許有下個資料點的真實結果，因此會放入真實資料。若在你的情境中，希望模型一次預測未來的數個資料點，則需要改放入預測資料，但需留意預測誤差可能會隨時間累積。

除了傳統統計模型以外，深度學習模型和樹模型等模型，也可以用來對時間序列資料進行預測。以下以太陽黑子資料集，示範如何用 XGBoost 處理時間序列資料：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.datasets import sunspots
from xgboost import XGBRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

data = sunspots.load().data
y = data["SUNACTIVITY"].to_numpy().astype(float)
n = len(y)

test_size = 44
train = y[:n-test_size]
test  = y[n-test_size:]


def make_features(series, n_lags):
	X, y_out = [], []
	for i in range(n_lags, len(series)):
		X.append(series[i - n_lags: i])
		y_out.append(series[i])
	return np.array(X), np.array(y_out)


n_lags = 11
X_train, y_train = make_features(train, n_lags)

model = XGBRegressor(n_estimators=100, random_state=0)
model.fit(X_train, y_train)

predictions = []
history = list(train)
for i in range(test_size):
	x = np.array(history[-n_lags:]).reshape(1, -1)
	pred = model.predict(x)[0]
	predictions.append(pred)
	history.append(test[i])
predictions = np.array(predictions)

mae = mean_absolute_error(test, predictions)
rmse = mean_squared_error(test, predictions) ** 0.5
mape = np.mean(np.abs((test - predictions) / (test + 1e-8))) * 100
print(f"MAE:  {mae:.2f}")
print(f"RMSE: {rmse:.2f}")
print(f"MAPE: {mape:.2f}%")

years = data["YEAR"].astype(int)
plt.plot(years[:n-test_size], train, label="Train")
plt.plot(years[n-test_size:], test, label="Test (ground truth)")
plt.plot(years[n-test_size:], predictions, label="Prediction")
plt.xlabel("Year")
plt.legend()
plt.show()

在上述範例中：

• 將資料集換成太陽黑子，是因為 XGBoost 等樹模型不擅長處理 AirPassengers 資料集中，預測未來時很可能需要外插（即測試集中的值，不在訓練資料的範圍中）的狀況。但你如果在很可能外插的資料集中，有必要使用樹模型的話，則有以下方式可以處理：
  • 先對資料做對數轉換，有機會稍微緩解此問題。
  • 先對資料做差分，讓模型學習變化量，而非絕對值；預測時，再將差分結果累加還原。
  • 先去除趨勢，再讓模型學習殘差。
• 使用 XGBoost 處理時間序列時，需要自己將序列用 sliding window 重新組合成適用的特徵格式。n_lags 變數代表特徵的維度，即 window 的長度；若設定的太小，則模型無法看到足夠的歷史規律；而若設定的太大，則可能引入相關性較低的遠期資訊。一般來說，可以先考慮設為資料週期的一至三倍，特別是整數倍。
• 模型在滾動中，改為不重新訓練。實務上，你可以視運算資源等狀況，來決定是否重新訓練。
• MAPE 的分母加上 1e-8，是為了防止原始資料點為零時，會發生除以零的錯誤。但需要注意的是，因為我們在分母加上了很小的值，因此此時計算出的 MAPE，可能會是非常大的值。

此外，由於時間序列資料的特殊性，若遇到資料中有缺失值時，由於相鄰時間點之間通常具有連續性，因此以前後相鄰的有效值來填補，往往會比用整體統計值填補更為合適；而當資料來自不同來源，或取樣頻率不一致時，則可能需要進行時間對齊與重採樣，例如將日資料累加為月資料，才容易進行後續分析。近年來，也出現了所謂的「時間序列基礎模型」(Time Series Foundation Model，TSFM)，其概念類似於自然語言處理領域的大語言模型，即先以大量來自不同領域的時間序列資料進行預訓練，之後便能對未曾見過的新資料進行零樣本(zero-shot)預測，而不需要針對特定資料集重新訓練模型。目前較具代表性的選擇包括 Google 的 TimesFM、Amazon 的 Chronos-2，以及 Salesforce 的 Moirai-2，三者在各大基準測試上的排名互有高低，且仍在持續更新競爭中。如果你的情境適合使用此類模型，可以考慮將其作為快速建立預測基準的起點，再視需要進行微調。