Rの関数:prophet {prophet}

# 必要なパッケージをロード
library(prophet)
library(dplyr)
library(ggplot2)

options(scipen = 999)

# 再現性のための乱数シード設定
seed <- 20250915
set.seed(seed)

# 1. データ期間の設定 (2020年1月1日 から 2023年12月31日)
ds <- seq(as.Date("2020-01-01"), as.Date("2023-12-31"), by = "day")
n <- length(ds)

# 2. 各コンポーネントの生成
# 2-1. トレンド成分 (2022年7月1日に傾きが変わる)
t1 <- 1:n
t2 <- pmax(0, t1 - which(ds == as.Date("2022-07-01")))
trend <- 100 + 0.2 * t1 + 0.6 * t2

# 2-2. 季節性成分 (年周期と週周期)
yearly_seasonality <- 20 * sin(2 * pi * (as.numeric(format(ds, "%j")) / 365.25 - 0.5))
weekly_seasonality <- 10 * case_when(
  weekdays(ds) == "土曜日" ~ -1,
  weekdays(ds) == "日曜日" ~ -1.5,
  TRUE ~ 0.2
)

# 2-3. 祝日効果データフレームの作成 (正月とゴールデンウィーク)
holidays_df <- bind_rows(
  data.frame(holiday = "正月", ds = as.Date(paste0(2020:2024, "-01-01"))),
  data.frame(holiday = "正月", ds = as.Date(paste0(2020:2024, "-01-02"))),
  data.frame(holiday = "正月", ds = as.Date(paste0(2020:2024, "-01-03"))),
  data.frame(holiday = "GW", ds = as.Date(paste0(2020:2024, "-05-03"))),
  data.frame(holiday = "GW", ds = as.Date(paste0(2020:2024, "-05-04"))),
  data.frame(holiday = "GW", ds = as.Date(paste0(2020:2024, "-05-05")))
)
holiday_effect <- ifelse(ds %in% holidays_df$ds[holidays_df$holiday == "正月"], -30, 0) +
  ifelse(ds %in% holidays_df$ds[holidays_df$holiday == "GW"], 40, 0)


# 2-4. ノイズ成分
noise <- rnorm(n, mean = 0, sd = 10)

# 3. yの値を合成
y <- trend + yearly_seasonality + weekly_seasonality + holiday_effect + noise

# 4. prophet用のデータフレームを作成
df <- data.frame(ds, y)

# 5. データの可視化
cat("作成したサンプルデータの一部を確認:\n")
print(head(df))

p <- ggplot(df, aes(x = ds, y = y)) +
  geom_line() +
  labs(title = "生成したサンプル時系列データ", subtitle = "縦鎖線はトレンドの傾きを変化させた時点", x = "日付", y = "値") +
  theme_minimal() +
  geom_vline(xintercept = as.Date("2022-07-01"), linetype = "dashed")
print(p)

set.seed(seed)

# 1. 日周期を含むサンプルデータ（時間単位）の作成
# 2023年1月1日から1月31日までの1時間ごとのデータ
ds_hourly <- seq(
  from = as.POSIXct("2023-01-01 00:00:00", tz = "Asia/Tokyo"),
  to = as.POSIXct("2023-01-31 23:00:00", tz = "Asia/Tokyo"),
  by = "hour"
)

n_hourly <- length(ds_hourly)

# トレンド、週周期、日周期、ノイズを生成
trend_hourly <- 100 + 0.02 * (1:n_hourly)
# 週末は下がる週周期
weekly_seasonality_hourly <- 10 * ifelse(weekdays(ds_hourly) %in% c("土曜日", "日曜日"), -1.5, 0.2)
# 午後3時(15時)頃にピークを迎え、午前3時頃に底になる日周期
daily_seasonality_hourly <- 20 * sin(2 * pi * (as.numeric(format(ds_hourly, "%H")) - 9) / 24)
noise_hourly <- rnorm(n_hourly, mean = 0, sd = 5)

# 全てを合成
y_hourly <- trend_hourly + weekly_seasonality_hourly + daily_seasonality_hourly + noise_hourly

df_hourly <- data.frame(ds = ds_hourly, y = y_hourly)

cat("作成した時間単位データの一部を確認:\n")
print(head(df_hourly))

# 2. prophetモデルの学習
# 時間単位データを渡すと、日周期が自動で有効になります
m5_daily <- prophet(df_hourly)

# 3. 未来のデータフレームを作成 (7日分)
# freq='hour' を指定して時間単位で未来を生成
future5_daily <- make_future_dataframe(m5_daily, periods = 7 * 24, freq = "hour")

# 4. 予測の実行
forecast5_daily <- predict(m5_daily, future5_daily)

# 5. 予測結果のプロット（最後の10日間を拡大）
p_forecast5 <- plot(m5_daily, forecast5_daily) +
  labs(title = "日周期データに対する予測結果", x = "日時", y = "予測値") +
  coord_cartesian(xlim = c(max(df_hourly$ds) - 10, max(future5_daily$ds))) + # 最後の10日間にズーム
  theme_minimal()
print(p_forecast5)

# 6. 成分のプロット
prophet_plot_components(m5_daily, forecast5_daily)