Posts

Pythonで株価予測 - 線形回帰とランダムフォレストで未来を予測

はじめに

株価は予測できるのか？

この問いには明確な答えはありませんが、機械学習を使って統計的な予測は可能です。

この記事では、Pythonのscikit-learnを使って、基本的な株価予測モデルを構築します。

注意：予測の限界

⚠️ 重要な免責事項

  • 過去のパターンが未来を保証しない
  • ブラックスワン（予測不可能なイベント）がある
  • このモデルは学習目的であり、実際の取引には使用しないこと

準備

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler

# データ取得
ticker = "AAPL"
df = yf.download(ticker, period="5y")

print(f"データ期間: {df.index[0]} ～ {df.index[-1]}")
print(f"データ数: {len(df)}")

特徴量エンジニアリング

テクニカル指標を特徴量として作成

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

def create_features(df):
    """
    特徴量を作成する関数
    """
    data = df.copy()
    
    # 終値の変化率（リターン）
    data['Return'] = data['Close'].pct_change()
    
    # 移動平均
    data['MA5'] = data['Close'].rolling(window=5).mean()
    data['MA20'] = data['Close'].rolling(window=20).mean()
    data['MA50'] = data['Close'].rolling(window=50).mean()
    
    # 移動平均との差
    data['Close_MA5'] = data['Close'] - data['MA5']
    data['Close_MA20'] = data['Close'] - data['MA20']
    
    # RSI
    delta = data['Close'].diff()
    gain = (delta.where(delta > 0, 0)).rolling(window=14).mean()
    loss = (-delta.where(delta < 0, 0)).rolling(window=14).mean()
    rs = gain / loss
    data['RSI'] = 100 - (100 / (1 + rs))
    
    # MACD
    ema12 = data['Close'].ewm(span=12, adjust=False).mean()
    ema26 = data['Close'].ewm(span=26, adjust=False).mean()
    data['MACD'] = ema12 - ema26
    data['MACD_Signal'] = data['MACD'].ewm(span=9, adjust=False).mean()
    
    # ボラティリティ（過去20日の標準偏差）
    data['Volatility'] = data['Return'].rolling(window=20).std()
    
    # 出来高の変化
    data['Volume_Change'] = data['Volume'].pct_change()
    data['Volume_MA20'] = data['Volume'].rolling(window=20).mean()
    
    # 高値・安値との位置
    data['High_Low_Range'] = (data['Close'] - data['Low']) / (data['High'] - data['Low'])
    
    # 翌日の終値（予測対象）
    data['Target'] = data['Close'].shift(-1)
    
    return data

# 特徴量作成
df_features = create_features(df)

# 最初の50日と最後の1日は削除（NaNがあるため）
df_features = df_features.iloc[50:-1]

print("作成された特徴量:")
print(df_features.columns.tolist())
print(f"\n特徴量の統計:")
print(df_features.describe())

線形回帰モデル

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

# 特徴量とターゲットを分離
feature_cols = ['Return', 'Close_MA5', 'Close_MA20', 'RSI', 'MACD', 
                'Volatility', 'Volume_Change', 'High_Low_Range']

X = df_features[feature_cols]
y = df_features['Target']

# 訓練データとテストデータに分割（時系列なのでシャッフルしない）
split_point = int(len(X) * 0.8)
X_train, X_test = X.iloc[:split_point], X.iloc[split_point:]
y_train, y_test = y.iloc[:split_point], y.iloc[split_point:]

# スケーリング
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 線形回帰モデル
lr_model = LinearRegression()
lr_model.fit(X_train_scaled, y_train)

# 予測
lr_pred = lr_model.predict(X_test_scaled)

# 評価
lr_rmse = np.sqrt(mean_squared_error(y_test, lr_pred))
lr_r2 = r2_score(y_test, lr_pred)

print("=== 線形回帰モデル ===")
print(f"RMSE: ${lr_rmse:.2f}")
print(f"R²: {lr_r2:.4f}")

# 特徴量の重要度
feature_importance = pd.DataFrame({
    'feature': feature_cols,
    'coefficient': lr_model.coef_
}).sort_values('coefficient', key=abs, ascending=False)

print("\n特徴量の係数:")
print(feature_importance)

ランダムフォレストモデル

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# ランダムフォレストモデル
rf_model = RandomForestRegressor(
    n_estimators=100,
    max_depth=10,
    min_samples_split=5,
    random_state=42
)

rf_model.fit(X_train_scaled, y_train)

# 予測
rf_pred = rf_model.predict(X_test_scaled)

# 評価
rf_rmse = np.sqrt(mean_squared_error(y_test, rf_pred))
rf_r2 = r2_score(y_test, rf_pred)

print("\n=== ランダムフォレストモデル ===")
print(f"RMSE: ${rf_rmse:.2f}")
print(f"R²: {rf_r2:.4f}")

# 特徴量の重要度
rf_importance = pd.DataFrame({
    'feature': feature_cols,
    'importance': rf_model.feature_importances_
}).sort_values('importance', ascending=False)

print("\n特徴量の重要度:")
print(rf_importance)

予測結果の可視化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))

# 実際の価格と予測価格
dates = y_test.index

ax1.plot(dates, y_test.values, label='Actual', color='black', linewidth=2)
ax1.plot(dates, lr_pred, label='Linear Regression', color='blue', alpha=0.7)
ax1.plot(dates, rf_pred, label='Random Forest', color='green', alpha=0.7)
ax1.set_title('Stock Price Prediction')
ax1.set_ylabel('Price (USD)')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 予測誤差
lr_error = y_test.values - lr_pred
rf_error = y_test.values - rf_pred

ax2.plot(dates, lr_error, label='Linear Regression Error', color='blue', alpha=0.5)
ax2.plot(dates, rf_error, label='Random Forest Error', color='green', alpha=0.5)
ax2.axhline(y=0, color='red', linestyle='--', alpha=0.5)
ax2.fill_between(dates, lr_error, alpha=0.2, color='blue')
ax2.set_title('Prediction Error')
ax2.set_ylabel('Error (USD)')
ax2.set_xlabel('Date')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('stock_prediction.png', dpi=150)
plt.show()

翌日の価格予測

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 最新のデータで翌日を予測
latest_features = X.iloc[-1:].values
latest_features_scaled = scaler.transform(latest_features)

lr_next_day = lr_model.predict(latest_features_scaled)[0]
rf_next_day = rf_model.predict(latest_features_scaled)[0]
actual_latest = df['Close'].iloc[-1]

print("=== 翌日の価格予測 ===")
print(f"現在の終値: ${actual_latest:.2f}")
print(f"線形回帰予測: ${lr_next_day:.2f} ({(lr_next_day/actual_latest-1)*100:+.2f}%)")
print(f"ランダムフォレスト予測: ${rf_next_day:.2f} ({(rf_next_day/actual_latest-1)*100:+.2f}%)")

モデルの評価と注意点

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 方向性の正解率（上がるか下がるか）
actual_direction = np.diff(y_test.values) > 0
lr_direction = np.diff(lr_pred) > 0
rf_direction = np.diff(rf_pred) > 0

lr_accuracy = np.mean(actual_direction == lr_direction)
rf_accuracy = np.mean(actual_direction == rf_direction)

print("\n=== 方向性予測の精度 ===")
print(f"ランダム（ベースライン）: 50.00%")
print(f"線形回帰: {lr_accuracy*100:.2f}%")
print(f"ランダムフォレスト: {rf_accuracy*100:.2f}%")

まとめ

モデルのポイント:

read more
Posts

Pythonで株価のボラティリティを予測する - GARCHモデル入門

はじめに

株価の**ボラティリティ（変動性）**は、リスク管理やオプション価格の決定に重要です。

「ボラティリティは常に変化する」という特徴を捉えるために、GARCHモデルが広く使われています。

この記事では、GARCHモデルの基礎とPython実装を解説します。

GARCHモデルとは？

ボラティリティの特徴

株価のボラティリティには以下の特徴があります：

  1. クラスタリング: 高いボラティリティの後に高いボラティリティが続きやすい
  2. 平均回帰: 長期的には平均レベルに戻る
  3. レバレッジ効果: 下落時の方がボラティリティが上がりやすい

GARCH(1,1)モデル

GARCH（Generalized Autoregressive Conditional Heteroskedasticity）モデルは、 条件付き分散（ボラティリティ）を以下の式でモデル化します：

$$\sigma_t^2 = \omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2$$

ここで：

  • $\sigma_t^2$ = 時点$t$での条件付き分散
  • $\omega$ = 定数項
  • $\alpha$ = ARCH項の係数（前日のリターンの影響）
  • $\beta$ = GARCH項の係数（前日のボラティリティの影響）
  • $r_{t-1}$ = 前日のリターン

制約条件:

  • $\omega > 0$
  • $\alpha \geq 0, \beta \geq 0$
  • $\alpha + \beta < 1$（安定性条件）

Pythonでの実装

ライブラリのインストール

1

pip install arch yfinance pandas numpy matplotlib

Step 1: データ取得と前処理

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from arch import arch_model

# S&P500のデータを取得
ticker = "^GSPC"
data = yf.download(ticker, start="2020-01-01", end="2024-12-31")['Close']

# 日次リターンを計算
returns = 100 * data.pct_change().dropna()  # パーセンテージ表記

print(f"データ期間: {returns.index[0]} ～ {returns.index[-1]}")
print(f"観測数: {len(returns)}")
print(f"平均リターン: {returns.mean():.4f}%")
print(f"リターンの標準偏差: {returns.std():.4f}%")

Step 2: リターンの可視化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 株価
axes[0].plot(data.index, data, label='S&P 500')
axes[0].set_title('S&P 500 Index')
axes[0].set_ylabel('Price')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# リターン
axes[1].plot(returns.index, returns, label='Daily Returns', alpha=0.7)
axes[1].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[1].set_title('Daily Returns (%)')
axes[1].set_ylabel('Return (%)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('returns_analysis.png', dpi=150)
plt.show()

Step 3: GARCH(1,1)モデルの推定

1
2
3
4
5
6
7

# GARCH(1,1)モデルを構築
model = arch_model(returns, vol='Garch', p=1, q=1, dist='normal')

# モデルをフィット
result = model.fit(update_freq=5)

print(result.summary())

Step 4: 結果の解釈

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# 推定されたパラメータ
params = result.params
omega = params['omega']
alpha = params['alpha[1]']
beta = params['beta[1]']

print(f"\n=== 推定パラメータ ===")
print(f"ω (定数項): {omega:.6f}")
print(f"α (ARCH項): {alpha:.4f}")
print(f"β (GARCH項): {beta:.4f}")
print(f"α + β (持続性): {alpha + beta:.4f}")
print(f"長期分散: {omega / (1 - alpha - beta):.4f}")

# 持続性の解釈
persistence = alpha + beta
if persistence > 0.95:
    print("\n高い持続性: ショックは長期間影響を残します")
elif persistence > 0.9:
    print("\n中程度の持続性: ショックは中期的に影響を残します")
else:
    print("\n低い持続性: ショックは短期間で消散します")

Step 5: 条件付きボラティリティの可視化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# 条件付きボラティリティを取得
conditional_volatility = result.conditional_volatility

fig, ax = plt.subplots(figsize=(12, 6))

ax.plot(returns.index, conditional_volatility, label='Conditional Volatility', color='red')
ax.set_ylabel('Volatility (%)')
ax.set_xlabel('Date')
ax.set_title('GARCH(1,1) Estimated Volatility')
ax.legend()
ax.grid(True, alpha=0.3)

# ボラティリティが高かった時期を強調
high_vol_periods = conditional_volatility[conditional_volatility > conditional_volatility.quantile(0.95)]
ax.scatter(high_vol_periods.index, high_vol_periods.values, color='orange', s=20, zorder=5)

plt.tight_layout()
plt.savefig('garch_volatility.png', dpi=150)
plt.show()

Step 6: ボラティリティの予測

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

# 将来のボラティリティを予測
forecast = result.forecast(horizon=30)

# 予測された分散
forecast_variance = forecast.variance.values[-1, :]
forecast_volatility = np.sqrt(forecast_variance)

# 予測期間の日付を生成
last_date = returns.index[-1]
forecast_dates = pd.date_range(start=last_date, periods=31, freq='B')[1:]

# プロット
fig, ax = plt.subplots(figsize=(12, 6))

# 過去のボラティリティ
ax.plot(returns.index[-60:], conditional_volatility[-60:], 
       label='Historical Volatility', color='blue')

# 予測ボラティリティ
ax.plot(forecast_dates, forecast_volatility, 
       label='Forecasted Volatility', color='red', linestyle='--')

# 信頼区間（仮に±2標準偏差）
ax.fill_between(forecast_dates, 
               forecast_volatility * 0.8, 
               forecast_volatility * 1.2,
               alpha=0.2, color='red', label='Confidence Interval')

ax.set_ylabel('Volatility (%)')
ax.set_xlabel('Date')
ax.set_title('Volatility Forecast (Next 30 Days)')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('volatility_forecast.png', dpi=150)
plt.show()

print(f"\n=== 30日先のボラティリティ予測 ===")
print(f"現在のボラティリティ: {conditional_volatility.iloc[-1]:.2f}%")
print(f"30日後の予測: {forecast_volatility[-1]:.2f}%")

モデルの診断

残差の検定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# 標準化残差
std_residuals = result.resid / conditional_volatility

# 残差の自己相関を検定
from statsmodels.stats.diagnostic import acorr_ljungbox

ljung_box = acorr_ljungbox(std_residuals**2, lags=10, return_df=True)
print("\n=== 残差のLjung-Box検定（ボラティリティの残差）===")
print(ljung_box)

# p値が0.05以上なら、残差に自己相関なし（モデルは適切）

他のモデルとの比較

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# EGARCHモデル（非対称性を考慮）
egarch_model = arch_model(returns, vol='EGarch', p=1, q=1)
egarch_result = egarch_model.fit(disp='off')

# モデル比較
print("\n=== モデル比較 ===")
print(f"GARCH(1,1) AIC: {result.aic:.2f}")
print(f"EGARCH(1,1) AIC: {egarch_result.aic:.2f}")

if result.aic < egarch_result.aic:
    print("GARCH(1,1)がより良いフィット")
else:
    print("EGARCH(1,1)がより良いフィット（レバレッジ効果あり）")

よくあるエラーと解決法

エラー1: “ConvergenceWarning”

原因: モデルが収束しない

read more