Random_forest(dft.4)_TAIEX_0628

加入 beta,beta_120

beta_120的離散大

from google.colab import drive
drive.mount('/content/drive')

!pip install yfinance scikit-learn matplotlib statsmodels
import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, roc_curve, roc_auc_score
from sklearn.decomposition import PCA
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint
import statsmodels.api as sm

# 定義 MSCI30 指數成分股代碼和大盤指數
tickers = ["2330.TW", "2454.TW", "2317.TW", "2412.TW", "1303.TW", "2882.TW", "3008.TW", "2308.TW", "1402.TW",
           "1216.TW", "2881.TW", "2891.TW", "2382.TW", "2409.TW", "1802.TW", "1101.TW", "3045.TW", "2324.TW",
           "2105.TW", "2880.TW", "2887.TW", "2885.TW", "4904.TW", "2603.TW", "2884.TW", "2886.TW", "2357.TW",
           "2344.TW", "4938.TW", "2888.TW", "^TWII"]

# 下載股票數據
data = yf.download(tickers, start="2021-01-01", end="2024-06-28")

# 使用前向填充處理缺失值
data = data.ffill()

# 提取調整後收盤價
adj_close = data['Adj Close']

# 計算日變動率
daily_change = adj_close.pct_change()

# 計算 Beta 值
def calculate_beta(stock_returns, market_returns, window):
    cov_matrix = stock_returns.rolling(window).cov(market_returns)
    var_market = market_returns.rolling(window).var()
    beta = cov_matrix.div(var_market, axis=0)
    return beta

# 市場回報率
market_returns = daily_change["^TWII"]

# 計算 Beta 和 Beta_120
beta = daily_change.apply(lambda x: calculate_beta(x, market_returns, 60))
beta_120 = daily_change.apply(lambda x: calculate_beta(x, market_returns, 120))

# 整合所有特徵變數
features = pd.DataFrame()
for ticker in tickers[:-1]:  # 除去 "^TWII"
    features[f'{ticker}_beta'] = beta[ticker]
    features[f'{ticker}_beta_120'] = beta_120[ticker]

# 使用前向填充處理缺失值
features = features.ffill()

# 增加技術指標特徵，只保留 MA7，刪除 MA21、MACD、BB_upper 和 BB_lower 變數
for ticker in tickers[:-1]:
    features[f'{ticker}_MA7'] = adj_close[ticker].rolling(window=7).mean()
    features[f'{ticker}_RSI14'] = (100 - (100 / (1 + adj_close[ticker].pct_change().rolling(window=14).mean())))

# 增加前一天的 TAIEX 值
features['Previous_TAIEX'] = adj_close['^TWII'].shift(1)

# 使用前向填充處理缺失值
features = features.ffill()

# 確保没有 NaN 值
features = features.fillna(0)

# 設置目標變數
target = adj_close["^TWII"]

# 構建線性回歸模型
X = features  # 使用之前處理好的特徵
y = target    # 使用調整後的收盤價作為目標變數

# 添加常數項
X = sm.add_constant(X)

# 構建模型
model = sm.OLS(y, X).fit()

# 顯示模型總結
print(model.summary())

# 數據標準化
scaler = MinMaxScaler(feature_range=(0, 1))
scaled_features = scaler.fit_transform(features)
scaled_target = scaler.fit_transform(target.values.reshape(-1, 1))

# 確認没有出現nan值
print("Checking for NaN values in scaled_features and scaled_target:")
print(np.isnan(scaled_features).sum())
print(np.isnan(scaled_target).sum())

# PCA 降維
pca = PCA(n_components=20)
scaled_features_pca = pca.fit_transform(scaled_features)

# 創建時間序列數據
def create_dataset(X, y, time_step=1):
    Xs, ys = [], []
    for i in range(len(X)-time_step):
        Xs.append(X[i:(i+time_step), :])
        ys.append(y[i + time_step])
    return np.array(Xs), np.array(ys)

# 定義時間步長
time_step = 60

# 創建訓練和測試數據集
train_size = int(len(scaled_features_pca) * 0.8)
X_train, y_train = create_dataset(scaled_features_pca[:train_size], scaled_target[:train_size], time_step)
X_test, y_test = create_dataset(scaled_features_pca[train_size:], scaled_target[train_size:], time_step)

# 將數據調整為2D以適應隨機森林
X_train_2d = X_train.reshape(X_train.shape[0], -1)
X_test_2d = X_test.reshape(X_test.shape[0], -1)

# 建立隨機森林模型
model = RandomForestRegressor(random_state=42)

# 使用 RandomizedSearchCV 調整隨機森林的參數
param_dist = {
    'n_estimators': randint(50, 100),  # 減少樹木數量
    'max_depth': [10, 20, None],
    'min_samples_split': randint(2, 5),  # 減少樣本拆分數量
    'min_samples_leaf': randint(1, 3),  # 減少葉子數量
    'bootstrap': [True, False]
}

random_search = RandomizedSearchCV(estimator=model, param_distributions=param_dist, n_iter=20, cv=2, n_jobs=-1, verbose=2, random_state=42)
random_search.fit(X_train_2d, y_train.ravel())

best_model = random_search.best_estimator_

# 預測
train_predict = best_model.predict(X_train_2d)
test_predict = best_model.predict(X_test_2d)

# 反標準化
train_predict = scaler.inverse_transform(train_predict.reshape(-1, 1))
test_predict = scaler.inverse_transform(test_predict.reshape(-1, 1))
y_train = scaler.inverse_transform(y_train)
y_test = scaler.inverse_transform(y_test)

# 計算 MSE 和 MAE
train_mse = mean_squared_error(y_train, train_predict)
test_mse = mean_squared_error(y_test, test_predict)
train_mae = mean_absolute_error(y_train, train_predict)
test_mae = mean_absolute_error(y_test, test_predict)

print(f"Training MSE: {train_mse:.4f}")
print(f"Testing MSE: {test_mse:.4f}")
print(f"Training MAE: {train_mae:.4f}")
print(f"Testing MAE: {test_mae:.4f}")

# 二元分類
threshold = np.median(y_test)  # 設定閾值為測試數據的中位數
y_test_class = (y_test > threshold).astype(int)
test_predict_class = (test_predict > threshold).astype(int)

# 確保有兩個類別
if len(np.unique(y_test_class)) > 1:
    # 計算 ROC 曲線和 AUC 值
    fpr, tpr, _ = roc_curve(y_test_class, test_predict)
    roc_auc = roc_auc_score(y_test_class, test_predict)

    # 繪製 ROC 曲線
    plt.figure()
    plt.plot(fpr, tpr, color='darkorange', lw=2, label='ROC curve (area = %0.2f)' % roc_auc)
    plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Receiver Operating Characteristic (ROC)')
    plt.legend(loc="lower right")
    plt.show()
else:
    print("ROC AUC score is not defined as there is only one class present in y_test_class.")

# 調整索引和預測值的長度一致
train_predict_index = target.index[time_step:time_step + len(train_predict)]
test_predict_index = target.index[train_size + time_step:train_size + time_step + len(test_predict)]

# 限制每日最大漲跌幅為 10%
def limit_change(predictions, previous_value, max_change=0.1):
    limited_predictions = []
    for pred in predictions:
        change = (pred - previous_value) / previous_value
        if change > max_change:
            pred = previous_value * (1 + max_change)
        elif change < -max_change:
            pred = previous_value * (1 - max_change)
        limited_predictions.append(pred)
        previous_value = pred
    return np.array(limited_predictions)

# 應用漲跌幅限制
train_predict = limit_change(train_predict, y_train[0])
test_predict = limit_change(test_predict, y_test[0])

# 繪製整個數據集的走勢圖
plt.figure(figsize=(14, 7))
plt.plot(target.index, target, color='blue', label='ACTUAL')
plt.plot(train_predict_index, train_predict, color='orange', label='Training')
plt.plot(test_predict_index, test_predict, color='red', label='Prediction')
plt.legend()  # 確保圖例顯示
plt.xlabel('Date')
plt.ylabel('TAIEX')
plt.title('TAIEX Prediction')
plt.grid(True)
plt.show()

# 延伸預測至未來 100 天
future_steps = 100
future_features = scaled_features_pca[-time_step:]
future_predict = []

for i in range(future_steps):
    future_features_2d = future_features.reshape(1, -1)
    pred = best_model.predict(future_features_2d)
    future_predict.append(pred[0])
    pred = np.repeat(pred, future_features.shape[1]).reshape(1, future_features.shape[1])
    future_features = np.append(future_features[1:, :], pred, axis=0)

future_predict = np.array(future_predict).reshape(-1, 1)
future_predict = scaler.inverse_transform(future_predict)
future_dates = pd.date_range(start=target.index[-1], periods=future_steps+1, inclusive='right')

# 繪製未來 100 天的預測走勢圖
plt.figure(figsize=(14, 7))
plt.plot(target.index, target, color='blue', label='ACTUAL')
plt.plot(train_predict_index, train_predict, color='orange', label='Training')
plt.plot(test_predict_index, test_predict, color='red', label='Prediction')
plt.plot(future_dates, future_predict, color='green', label='Future Prediction')
plt.legend()  # 確保圖例顯示
plt.xlabel('Date')
plt.ylabel('TAIEX')
plt.title('TAIEX Future Prediction（100D）')
plt.grid(True)
plt.show()

Time Step 調為 30