時系列予測の魔法

# 時系列分析エコシステムの読み込み library(tidyverse) library(lubridate) # 時系列分析パッケージ library(timetk) library(modeltime) # サンプルデータの作成 set.seed(123) # 期間の設定 start_date <- as.Date("2019-01-01") end_date <- as.Date("2023-12-31") dates <- seq(start_date, end_date, by = "day") # 時系列データの作成 sales_ts <- tibble(date = dates) %>% mutate(sales = 1000 + rnorm(n(), 0, 100)) # データの確認 glimpse(sales_ts)

# 基本統計 ts_summary <- sales_ts %>% summarise(mean_sales = mean(sales)) print(ts_summary) # 月次集計 monthly_sales <- sales_ts %>% mutate(yearmon = floor_date(date, "month")) %>% group_by(yearmon) %>% summarise(monthly_total = sum(sales), .groups = 'drop') print(head(monthly_sales))

# データ分割 splits <- sales_ts %>% time_series_split(assess = "3 months", cumulative = TRUE) train_data <- training(splits) test_data <- testing(splits) # ARIMAモデル model_arima <- arima_reg() %>% set_engine("auto_arima") %>% fit(sales ~ date, data = train_data) # 指数平滑化モデル model_ets <- exp_smoothing() %>% set_engine("ets") %>% fit(sales ~ date, data = train_data)

library(anomalize) library(timetk) # 季節性を含むデータの作成 seasonal_sales <- tibble( date = seq.Date(as.Date("2020-01-01"), as.Date("2023-12-31"), by = "day") ) %>% mutate( day_of_year = yday(date), # 基本トレンド trend = 1000 + row_number() * 0.1, # 季節性（年次） seasonal = 200 * sin(2 * pi * day_of_year / 365), # 週次パターン weekly = 50 * sin(2 * pi * wday(date) / 7), # ノイズ noise = rnorm(n(), 0, 30), # 売上総計 sales = trend + seasonal + weekly + noise ) # 季節性分解 decomposition <- seasonal_sales %>% tk_ts(start = 2020, freq = 365.25) %>% decompose() %>% sw_tidy() glimpse(decomposition) # 異常検知 anomalies <- seasonal_sales %>% time_decompose(sales, method = "stl") %>% anomalize(remainder, method = "iqr") %>% time_recompose() # 異常値の確認 anomaly_count <- anomalies %>% filter(anomaly == "Yes") %>% nrow() cat("検出された異常値の数:", anomaly_count)

library(prophet) library(randomForest) # データ分割（訓練・検証・テスト） train_end <- as.Date("2022-12-31") valid_end <- as.Date("2023-06-30") train_data <- seasonal_sales %>% filter(date <= train_end) valid_data <- seasonal_sales %>% filter(date > train_end, date <= valid_end) test_data <- seasonal_sales %>% filter(date > valid_end) # 1. ARIMAモデル arima_model <- arima_reg() %>% set_engine("auto_arima") %>% fit(sales ~ date, data = train_data) # 2. Prophetモデル prophet_model <- prophet_reg( seasonality_yearly = TRUE, seasonality_weekly = TRUE ) %>% set_engine("prophet") %>% fit(sales ~ date, data = train_data) # 3. 指数平滑化モデル ets_model <- exp_smoothing( seasonal_period = 365 ) %>% set_engine("ets") %>% fit(sales ~ date, data = train_data) # 4. 機械学習アプローチ（特徴量エンジニアリング） ml_features <- train_data %>% tk_augment_timeseries_signature() %>% select(sales, date, index.num, year, quarter, month, wday, mday, yday, week, hour, minute, second) rf_model <- rand_forest(trees = 500) %>% set_engine("randomForest") %>% set_mode("regression") %>% fit(sales ~ ., data = select(ml_features, -date))

# modeltimeテーブル作成 model_table <- modeltime_table( arima_model, prophet_model, ets_model, rf_model ) # 精度評価（検証データ） accuracy_results <- model_table %>% modeltime_accuracy(valid_data) print(accuracy_results) # 最良モデルの選択（MAEが最小） best_model_id <- accuracy_results %>% arrange(mae) %>% slice(1) %>% pull(.model_id) cat("最良モデル ID:", best_model_id) # 未来予測（3ヶ月先） forecast_results <- model_table %>% modeltime_forecast( h = "3 months", actual_data = seasonal_sales ) # 予測結果の可視化 forecast_plot <- forecast_results %>% plot_modeltime_forecast(.interactive = FALSE) print(forecast_plot)

library(plotly) library(DT) # KPI計算 kpi_summary <- seasonal_sales %>% mutate( year_month = floor_date(date, "month") ) %>% group_by(year_month) %>% summarise( 月間売上 = sum(sales), 平均日売上 = mean(sales), 売上標準偏差 = sd(sales), 変動係数 = 売上標準偏差 / 平均日売上, .groups = 'drop' ) %>% mutate( 前月比成長率 = (月間売上 / lag(月間売上) - 1) * 100, 売上トレンド = case_when( 前月比成長率 > 5 ~ "📈 上昇", 前月比成長率 < -5 ~ "📉 下降", TRUE ~ "➡️ 横ばい" ) ) # 最新3ヶ月のサマリー recent_kpi <- kpi_summary %>% tail(3) print(recent_kpi) # 予測精度のビジネス価値計算 business_impact <- accuracy_results %>% mutate( 予測精度ランク = rank(mae), コスト削減潜在力 = (max(mae) - mae) / max(mae) * 100, 推奨度 = case_when( 予測精度ランク == 1 ~ "🏆 最優秀", 予測精度ランク == 2 ~ "🥈 優秀", 予測精度ランク == 3 ~ "🥉 良好", TRUE ~ "⚠️ 要改善" ) ) %>% select(.model_desc, mae, mape, 推奨度, コスト削減潜在力) print(business_impact)