第21章: 時系列分析とmodeltime

# modeltimeエコシステムの読み込み library(tidymodels) library(modeltime) library(modeltime.ensemble) library(timetk) library(tidyverse) # サンプル時系列データの準備 library(tidyquant) # 株価データの取得例 stock_data <- tq_get("AAPL", get = "stock.prices", from = "2020-01-01", to = "2023-12-31") %>% select(date, adjusted) %>% rename(value = adjusted) glimpse(stock_data)

# 時系列データの分割（最後の12か月をテスト用） splits <- stock_data %>% time_series_split( assess = "12 months", cumulative = TRUE ) # 分割結果の可視化 splits %>% tk_time_series_cv_plan() %>% plot_time_series_cv_plan(date, value, .interactive = TRUE) # 訓練・テストデータの取得 train_data <- training(splits) test_data <- testing(splits) cat("訓練データ期間:", min(train_data$date), "〜", max(train_data$date), "\\n") cat("テストデータ期間:", min(test_data$date), "〜", max(test_data$date), "\\n")

# 時系列クロスバリデーション設定 resamples <- stock_data %>% time_series_cv( initial = "3 years", assess = "12 months", skip = "3 months", cumulative = TRUE, slice_limit = 6 ) # CV分割の可視化 resamples %>% tk_time_series_cv_plan() %>% plot_time_series_cv_plan( date, value, .facet_ncol = 2, .interactive = TRUE ) print(resamples)

# 時系列署名の生成 ts_features <- train_data %>% tk_augment_timeseries_signature(date) %>% select(-contains("hour|minute|second|am.pm")) %>% mutate( # ラグ特徴量 lag_1 = lag_vec(value, 1), lag_7 = lag_vec(value, 7), lag_14 = lag_vec(value, 14), lag_30 = lag_vec(value, 30), # 移動平均特徴量 ma_7 = slide_vec(value, mean, .before = 6, .complete = TRUE), ma_14 = slide_vec(value, mean, .before = 13, .complete = TRUE), ma_30 = slide_vec(value, mean, .before = 29, .complete = TRUE), # 差分特徴量 diff_1 = diff_vec(value, lag = 1), diff_7 = diff_vec(value, lag = 7), # 変化率特徴量 pct_change_1 = (value / lag_vec(value, 1) - 1) * 100, pct_change_7 = (value / lag_vec(value, 7) - 1) * 100 ) %>% filter(!is.na(lag_30)) glimpse(ts_features)

# フーリエ変換特徴量の生成 fourier_features <- ts_features %>% tk_augment_fourier(date, .periods = c(365.25, 365.25/2), .K = 2) # 滞留期間特徴量 holiday_features <- fourier_features %>% tk_augment_holiday_signature(date) %>% mutate( # カスタム特徴量 is_weekend = ifelse(wday.lbl %in% c("Saturday", "Sunday"), 1, 0), is_month_end = ifelse(mday == days_in_month(date), 1, 0), quarter_start = ifelse(month %in% c(1, 4, 7, 10) & mday == 1, 1, 0), # 技術指標風の特徴量 rsi_14 = slide_vec(pct_change_1, ~{ gains <- pmax(.x, 0) losses <- pmax(-.x, 0) avg_gain <- mean(gains, na.rm = TRUE) avg_loss <- mean(losses, na.rm = TRUE) rs <- avg_gain / avg_loss rsi <- 100 - (100 / (1 + rs)) return(rsi) }, .before = 13, .complete = TRUE) ) # 特徴量重要性の確認 feature_names <- holiday_features %>% select(-date, -value) %>% names() cat("生成された特徴量数:", length(feature_names), "\\n")

# 1. ARIMAモデル model_arima <- arima_reg() %>% set_engine("auto_arima") %>% fit(value ~ date, data = train_data) # 2. Prophetモデル model_prophet <- prophet_reg( seasonality_yearly = TRUE, seasonality_weekly = FALSE, seasonality_daily = FALSE ) %>% set_engine("prophet") %>% fit(value ~ date, data = train_data) # 3. 機械学習モデル（Random Forest） recipe_ts <- recipe(value ~ date, data = train_data) %>% step_timeseries_signature(date) %>% step_rm(contains("hour|minute|second|am.pm")) %>% step_normalize(all_numeric_predictors()) %>% step_dummy(all_nominal_predictors()) model_rf <- rand_forest(trees = 500, min_n = 10) %>% set_engine("ranger") %>% set_mode("regression") wflow_rf <- workflow() %>% add_recipe(recipe_ts) %>% add_model(model_rf) %>% fit(train_data)

# modeltimeテーブルの作成 model_table <- modeltime_table( model_arima, model_prophet, wflow_rf ) # テストデータでのキャリブレーション calibration_table <- model_table %>% modeltime_calibrate(new_data = test_data) # 精度指標の計算 accuracy_table <- calibration_table %>% modeltime_accuracy() print(accuracy_table) # 予測結果の可視化 calibration_table %>% modeltime_forecast( new_data = test_data, actual_data = train_data, keep_data = TRUE ) %>% plot_modeltime_forecast( .facet_ncol = 2, .interactive = TRUE )

# 平均アンサンブルの作成 ensemble_avg <- calibration_table %>% ensemble_average(type = "mean") # 重み付きアンサンブル（精度逆数重み） weights <- 1 / accuracy_table$mae weights <- weights / sum(weights) ensemble_weighted <- calibration_table %>% ensemble_weighted(loadings = weights) # 最終モデルテーブル final_table <- modeltime_table( ensemble_avg, ensemble_weighted ) %>% combine_modeltime_tables(calibration_table) %>% modeltime_calibrate(new_data = test_data) # アンサンブル結果の評価 final_accuracy <- final_table %>% modeltime_accuracy() %>% arrange(mae) print(final_accuracy)