第12章: 時系列分析とmodeltime

library(tidyverse) library(lubridate) library(timetk) library(modeltime) # 複雑な売上時系列データの作成 set.seed(123) # 2019年1月から2023年12月までの日次データ 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( # 基本トレンド（年間成長率5%） trend = 1000 + as.numeric(date - start_date) * 0.1, # 年次季節性（年末年始、夏のボーナス時期） yearly_seasonal = 200 * sin(2 * pi * yday(date) / 365) + 150 * cos(4 * pi * yday(date) / 365), # 週次季節性（平日高、週末低） weekly_seasonal = case_when( wday(date) %in% c(1, 7) ~ -100, # 週末 wday(date) == 6 ~ -50, # 金曜日 TRUE ~ 50 # 平日 ), # 特別な日（祝日効果など） special_events = case_when( month(date) == 12 && day(date) %in% c(24, 25) ~ 300, # クリスマス month(date) == 1 && day(date) == 1 ~ -200, # 元日 month(date) == 7 && day(date) %in% c(15:25) ~ 100, # 夏のセール TRUE ~ 0 ), # ランダムノイズ noise = rnorm(n(), 0, 80), # 最終的な売上 sales = trend + yearly_seasonal + weekly_seasonal + special_events + noise ) %>% mutate(sales = pmax(sales, 100)) %>% # 最小値設定 select(date, sales) print(head(sales_ts, 10)) # 基本的な時系列統計 ts_summary <- sales_ts %>% summarise( start_date = min(date), end_date = max(date), total_days = n(), mean_sales = mean(sales), median_sales = median(sales), sd_sales = sd(sales), min_sales = min(sales), max_sales = max(sales) ) print("時系列データ概要:") print(ts_summary) # 月次集計 monthly_sales <- sales_ts %>% mutate( year = year(date), month = month(date), yearmon = floor_date(date, "month") ) %>% group_by(yearmon) %>% summarise( monthly_total = sum(sales), monthly_avg = mean(sales), days_count = n(), .groups = 'drop' ) print("月次売上（最初の12ヶ月）:") print(head(monthly_sales, 12)) # 年次成長率の計算 yearly_growth <- sales_ts %>% mutate(year = year(date)) %>% group_by(year) %>% summarise(yearly_total = sum(sales), .groups = 'drop') %>% mutate( growth_rate = (lag(yearly_total, 1) - yearly_total) / lag(yearly_total, 1) * 100 ) %>% filter(!is.na(growth_rate)) print("年次成長率:") print(yearly_growth) # 週次・月次パターンの分析 pattern_analysis <- sales_ts %>% mutate( weekday = wday(date, label = TRUE, abbr = FALSE), month_name = month(date, label = TRUE, abbr = FALSE) ) weekday_pattern <- pattern_analysis %>% group_by(weekday) %>% summarise( avg_sales = mean(sales), median_sales = median(sales), .groups = 'drop' ) %>% arrange(desc(avg_sales)) print("曜日別売上パターン:") print(weekday_pattern) monthly_pattern <- pattern_analysis %>% group_by(month_name) %>% summarise( avg_sales = mean(sales), median_sales = median(sales), .groups = 'drop' ) print("月別売上パターン:") print(monthly_pattern)

# データの分割（訓練用とテスト用） splits <- sales_ts %>% time_series_split(assess = "3 months", cumulative = TRUE) train_data <- training(splits) test_data <- testing(splits) print(paste("訓練期間:", min(train_data$date), "to", max(train_data$date))) print(paste("テスト期間:", min(test_data$date), "to", max(test_data$date))) # モデル1: ARIMA model_arima <- arima_reg() %>% set_engine("auto_arima") %>% fit(sales ~ date, data = train_data) # モデル2: 指数平滑法（ETS） model_ets <- exp_smoothing() %>% set_engine("ets") %>% fit(sales ~ date, data = train_data) # モデル3: 線形回帰（トレンド + 季節性） model_lm <- linear_reg() %>% set_engine("lm") %>% fit(sales ~ as.numeric(date) + wday(date, label = TRUE) + month(date, label = TRUE), data = train_data) # モデル4: Prophet model_prophet <- prophet_reg() %>% set_engine("prophet") %>% fit(sales ~ date, data = train_data) # モデル5: 季節性ナイーブ model_snaive <- naive_reg() %>% set_engine("snaive") %>% fit(sales ~ date, data = train_data) # modeltimeテーブルの作成 models_tbl <- modeltime_table( model_arima, model_ets, model_lm, model_prophet, model_snaive ) print("モデル一覧:") print(models_tbl) # 予測の実行 forecast_tbl <- models_tbl %>% modeltime_forecast( new_data = test_data, actual_data = train_data, keep_data = TRUE ) print("予測結果（サンプル）:") print(head(forecast_tbl, 10)) # 精度評価 accuracy_tbl <- models_tbl %>% modeltime_accuracy(test_data) print("モデル精度比較:") print(accuracy_tbl) # 最良モデルの選択と将来予測 best_model_id <- accuracy_tbl %>% filter(mae == min(mae)) %>% pull(.model_id) print(paste("最良モデル ID:", best_model_id)) # 将来予測（3ヶ月先まで） future_forecast <- models_tbl %>% filter(.model_id == best_model_id) %>% modeltime_refit(sales_ts) %>% modeltime_forecast(h = "3 months", actual_data = sales_ts) print("将来予測結果:") future_summary <- future_forecast %>% filter(.key == "prediction") %>% summarise( start_date = min(.index), end_date = max(.index), avg_forecast = mean(.value), total_forecast = sum(.value) ) print(future_summary) # 季節性分解の実行 decomposition <- sales_ts %>% plot_seasonal_diagnostics(.date_var = date, .value = sales, .feature_set = c("week", "month.lbl")) # 季節性統計の計算 seasonal_stats <- sales_ts %>% mutate( month = month(date, label = TRUE), wday = wday(date, label = TRUE) ) %>% group_by(month) %>% summarise( seasonal_index = mean(sales) / mean(sales_ts$sales), avg_sales = mean(sales), .groups = 'drop' ) %>% arrange(desc(seasonal_index)) print("月別季節指数:") print(seasonal_stats)

# 異常値検出 anomaly_detection <- sales_ts %>% plot_anomaly_diagnostics(.date_var = date, .value = sales, .alpha = 0.05, .max_anomalies = 0.1) # 異常値の特定 anomalies <- sales_ts %>% tk_anomaly_diagnostics(.date_var = date, .value = sales) %>% filter(anomaly == "Yes") print(paste("検出された異常値数:", nrow(anomalies))) print("異常値の例:") print(head(anomalies, 5)) # 外部変数（例：広告費、気温など）を含むデータの作成 enhanced_sales <- sales_ts %>% mutate( # 広告費（週末に多く投入） ad_spend = case_when( wday(date) %in% c(1, 7) ~ rnorm(n(), 50000, 10000), TRUE ~ rnorm(n(), 30000, 8000) ), # 気温（季節性あり） temperature = 20 + 15 * sin(2 * pi * yday(date) / 365) + rnorm(n(), 0, 3), # 競合他社のキャンペーン（ランダム） competitor_campaign = rbinom(n(), 1, 0.1), # 祝日フラグ is_holiday = case_when( month(date) == 1 && day(date) %in% c(1:3) ~ 1, month(date) == 5 && day(date) %in% c(3:5) ~ 1, month(date) == 12 && day(date) %in% c(23:25) ~ 1, TRUE ~ 0 ) ) %>% mutate(ad_spend = pmax(ad_spend, 10000)) # 外部変数を含む高度なモデル enhanced_splits <- enhanced_sales %>% time_series_split(assess = "3 months", cumulative = TRUE) enhanced_train <- training(enhanced_splits) enhanced_test <- testing(enhanced_splits) # 外部変数を含む線形回帰モデル model_enhanced_lm <- linear_reg() %>% set_engine("lm") %>% fit(sales ~ as.numeric(date) + wday(date, label = TRUE) + month(date, label = TRUE) + ad_spend + temperature + competitor_campaign + is_holiday, data = enhanced_train) # 外部変数を含むランダムフォレストモデル model_enhanced_rf <- rand_forest(trees = 500) %>% set_engine("ranger") %>% set_mode("regression") %>% fit(sales ~ as.numeric(date) + wday(date, label = TRUE) + month(date, label = TRUE) + ad_spend + temperature + competitor_campaign + is_holiday, data = enhanced_train) # 拡張モデルの評価 enhanced_models <- modeltime_table( model_enhanced_lm, model_enhanced_rf ) enhanced_accuracy <- enhanced_models %>% modeltime_accuracy(enhanced_test) print("外部変数を含むモデルの精度:") print(enhanced_accuracy) # 特徴量重要度（ランダムフォレストモデル） if (require(vip, quietly = TRUE)) { feature_importance <- model_enhanced_rf %>% extract_fit_parsnip() %>% vip::vi() %>% arrange(desc(Importance)) print("特徴量重要度:") print(feature_importance) } # 残差分析 residuals_analysis <- enhanced_models %>% modeltime_forecast(new_data = enhanced_test, actual_data = enhanced_train) %>% filter(.model_desc == "RANGER", .key == "prediction") %>% left_join(enhanced_test %>% select(date, actual_sales = sales), by = c(".index" = "date")) %>% mutate(residuals = actual_sales - .value) residual_stats <- residuals_analysis %>% summarise( mean_residual = mean(residuals, na.rm = TRUE), sd_residual = sd(residuals, na.rm = TRUE), rmse = sqrt(mean(residuals^2, na.rm = TRUE)), mae = mean(abs(residuals), na.rm = TRUE) ) print("残差統計:") print(residual_stats) # 予測区間の計算 prediction_intervals <- enhanced_models %>% filter(.model_id == 2) %>% # Random Forestモデル modeltime_forecast( new_data = enhanced_test, actual_data = enhanced_train, conf_interval = 0.95 ) interval_summary <- prediction_intervals %>% filter(.key == "prediction") %>% summarise( avg_prediction = mean(.value), avg_lower = mean(.conf_lo), avg_upper = mean(.conf_hi), interval_width = mean(.conf_hi - .conf_lo) ) print("予測区間の統計:") print(interval_summary)