Time Series Analysis

Overview

時系列分析は、時間を通じて収集されたデータポイントを調査し、予測と時間的動力学を理解するためのパターン、トレンド、季節性を特定します。

When to Use

• 過去のトレンドに基づいて将来の値を予測する
• データ内の季節性と周期的パターンを検出する
• 売上、株価、またはウェブサイトトラフィックの時間経過に伴うトレンドを分析する
• 自己相関と時間的依存関係を理解する
• 信頼区間付きで時間ベースの予測を行う
• データをトレンド、季節、残差成分に分解する

Core Components

• Trend: 長期的な方向性の動き
• Seasonality: 固定間隔での繰り返しパターン
• Cyclicity: 長期的な振動（非固定周期）
• Stationarity: 時間を通じて一定の平均と分散
• Autocorrelation: 過去の値との相関

Key Techniques

• Decomposition: トレンド、季節、残差成分の分離
• Differencing: データを定常化する
• ARIMA: 自己回帰和分移動平均モデル
• Exponential Smoothing: 過去の値の加重平均
• SARIMA: 季節的ARIMAモデル

Implementation with Python

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# Create sample time series data
dates = pd.date_range('2020-01-01', periods=365, freq='D')
values = 100 + np.sin(np.arange(365) * 2*np.pi / 365) * 20 + np.random.normal(0, 5, 365)
ts = pd.Series(values, index=dates)

# Visualize time series
fig, axes = plt.subplots(2, 2, figsize=(14, 8))

axes[0, 0].plot(ts)
axes[0, 0].set_title('Original Time Series')
axes[0, 0].set_ylabel('Value')

# Decomposition
decomposition = seasonal_decompose(ts, model='additive', period=30)
axes[0, 1].plot(decomposition.trend)
axes[0, 1].set_title('Trend Component')

axes[1, 0].plot(decomposition.seasonal)
axes[1, 0].set_title('Seasonal Component')

axes[1, 1].plot(decomposition.resid)
axes[1, 1].set_title('Residual Component')

plt.tight_layout()
plt.show()

# Test for stationarity (Augmented Dickey-Fuller)
result = adfuller(ts)
print(f"ADF Test Statistic: {result[0]:.6f}")
print(f"P-value: {result[1]:.6f}")
print(f"Critical Values: {result[4]}")

if result[1] <= 0.05:
    print("Time series is stationary")
else:
    print("Time series is non-stationary - differencing needed")

# First differencing for stationarity
ts_diff = ts.diff().dropna()
result_diff = adfuller(ts_diff)
print(f"\nAfter differencing - ADF p-value: {result_diff[1]:.6f}")

# Autocorrelation and Partial Autocorrelation
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

plot_acf(ts_diff, lags=40, ax=axes[0])
axes[0].set_title('ACF')

plot_pacf(ts_diff, lags=40, ax=axes[1])
axes[1].set_title('PACF')

plt.tight_layout()
plt.show()

# ARIMA Model
arima_model = ARIMA(ts, order=(1, 1, 1))
arima_result = arima_model.fit()
print(arima_result.summary())

# Forecast
forecast_steps = 30
forecast = arima_result.get_forecast(steps=forecast_steps)
forecast_df = forecast.conf_int()
forecast_mean = forecast.predicted_mean

# Plot forecast
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index[-90:], ts[-90:], label='Historical')
ax.plot(forecast_df.index, forecast_mean, label='Forecast', color='red')
ax.fill_between(
    forecast_df.index,
    forecast_df.iloc[:, 0],
    forecast_df.iloc[:, 1],
    color='red', alpha=0.2
)
ax.set_title('ARIMA Forecast with Confidence Interval')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

# Exponential Smoothing
exp_smooth = ExponentialSmoothing(
    ts, seasonal_periods=30, trend='add', seasonal='add', initialization_method='estimated'
)
exp_result = exp_smooth.fit()

# Model diagnostics
fig = exp_result.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

# Custom moving average analysis
window_sizes = [7, 30, 90]
fig, ax = plt.subplots(figsize=(12, 5))

ax.plot(ts.index, ts.values, label='Original', alpha=0.7)

for window in window_sizes:
    ma = ts.rolling(window=window).mean()
    ax.plot(ma.index, ma.values, label=f'MA({window})')

ax.set_title('Moving Averages')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

# Seasonal subseries plot
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, month in enumerate(range(1, 5)):
    month_data = ts[ts.index.month == month]
    axes[i // 2, i % 2].plot(month_data.values)
    axes[i // 2, i % 2].set_title(f'Month {month} Pattern')

plt.tight_layout()
plt.show()

# Forecast accuracy metrics
def calculate_forecast_metrics(actual, predicted):
    mae = np.mean(np.abs(actual - predicted))
    rmse = np.sqrt(np.mean((actual - predicted) ** 2))
    mape = np.mean(np.abs((actual - predicted) / actual)) * 100
    return {'MAE': mae, 'RMSE': rmse, 'MAPE': mape}

metrics = calculate_forecast_metrics(ts[-30:], forecast_mean[:30])
print(f"\nForecast Metrics:\n{metrics}")

# Additional analysis techniques

# Step 10: Seasonal subseries plots
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, season in enumerate([1, 2, 3, 4]):
    seasonal_ts = ts[ts.index.month % 4 == season % 4]
    axes[i // 2, i % 2].plot(seasonal_ts.values)
    axes[i // 2, i % 2].set_title(f'Season {season}')
plt.tight_layout()
plt.show()

# Step 11: Granger causality (for multiple series)
from statsmodels.tsa.stattools import grangercausalitytests

# Create another series for testing
ts2 = ts.shift(1).fillna(method='bfill')

try:
    print("\nGranger Causality Test:")
    print(f"Test whether ts2 Granger-causes ts:")
    gc_result = grangercausalitytests(np.column_stack([ts.values, ts2.values]), maxlag=3)
except Exception as e:
    print(f"Granger causality not performed: {str(e)[:50]}")

# Step 12: Autocorrelation and partial autocorrelation analysis
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

acf_values = acf(ts.dropna(), nlags=20)
pacf_values = pacf(ts.dropna(), nlags=20)

# Step 13: Seasonal strength
def seasonal_strength(series, seasonal_period=30):
    seasonal = seasonal_decompose(series, model='additive', period=seasonal_period)
    var_residual = np.var(seasonal.resid.dropna())
    var_seasonal = np.var(seasonal.seasonal)
    return 1 - (var_residual / (var_residual + var_seasonal)) if (var_residual + var_seasonal) > 0 else 0

ss = seasonal_strength(ts)
print(f"\nSeasonal Strength: {ss:.3f}")

# Step 14: Forecasting with uncertainty
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(ts.index[-60:], ts.values[-60:], label='Historical', linewidth=2)

# Multiple horizon forecasts
for steps_ahead in [10, 20, 30]:
    try:
        fc = arima_result.get_forecast(steps=steps_ahead)
        fc_mean = fc.predicted_mean
        ax.plot(pd.date_range(ts.index[-1], periods=steps_ahead+1)[1:],
               fc_mean.values, marker='o', label=f'Forecast (+{steps_ahead})')
    except:
        pass

ax.set_title('Multi-step Ahead Forecasts')
ax.set_xlabel('Date')
ax.set_ylabel('Value')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Step 15: Model comparison summary
print("\nTime Series Analysis Complete!")
print(f"Original series length: {len(ts)}")
print(f"Trend strength: {1 - np.var(decomposition.resid.dropna()) / np.var((ts - ts.mean()).dropna()):.3f}")
print(f"Seasonal strength: {ss:.3f}")

Stationarity

• Stationary: 時間を通じて平均、分散、自己相関が一定
• Non-stationary: トレンドまたは季節パターンが存在
• Solution: 差分、対数変換、またはトレンド除去

Model Selection

• ARIMA: 単変量予測に最適
• SARIMA: 季節成分を含む
• Exponential Smoothing: よりシンプル、トレンドに最適
• Prophet: 祝日と変化点を処理

Evaluation Metrics

• MAE: 平均絶対誤差
• RMSE: 二乗平均平方根誤差
• MAPE: 平均絶対パーセント誤差

Deliverables

• 分解分析チャート
• 定常性テスト結果
• ACF/PACFプロット
• 診断を含むフィットされたモデル
• 信頼区間付き予測
• 精度メトリクスの比較