Rによる季節成分とトレンド成分への分解

# サンプルデータのパラメータ設定
seed <- 20250508
n_years <- 5
frequency <- 12 # 月次データ
n_months <- n_years * frequency

# 基本となる季節パターン
base_seasonal_pattern <- numeric(12)
base_seasonal_pattern[12] <- 3.0 # December
base_seasonal_pattern[6] <- 2.0 # June
base_seasonal_pattern[3] <- 1.0 # March
base_seasonal_pattern[9] <- -2.5 # September
base_seasonal_pattern[1] <- -0.5
base_seasonal_pattern[2] <- -0.6
base_seasonal_pattern[4] <- -0.3
base_seasonal_pattern[5] <- -0.7
base_seasonal_pattern[7] <- -0.4
base_seasonal_pattern[8] <- -0.5
base_seasonal_pattern[10] <- -0.2
base_seasonal_pattern[11] <- -0.3

# 季節成分の作成 (振幅が年々変化)
seasonal_component_varying <- numeric(n_months)
# 年ごとの振幅スケールファクター (ここでは 0.8から5.2へ線形に増加)
amplitude_scale_factors <- seq(from = 0.8, to = 5.2, length.out = n_years)

for (i in 0:(n_years - 1)) {
  year_start_index <- i * frequency + 1
  year_end_index <- (i + 1) * frequency

  # 現在の年の振幅スケール
  current_scale <- amplitude_scale_factors[i + 1]

  # スケールを適用した季節パターン
  scaled_pattern_for_year <- base_seasonal_pattern * current_scale

  # 各年のパターンが合計0になるように中心化
  centered_pattern_for_year <- scaled_pattern_for_year - mean(scaled_pattern_for_year)

  seasonal_component_varying[year_start_index:year_end_index] <- centered_pattern_for_year
}

# トレンド成分の作成
trend_component_varying_season <- numeric(n_months)
current_trend_value <- 5
for (i in 0:(n_years - 1)) {
  year_start_index <- i * frequency + 1
  year_end_index <- (i + 1) * frequency
  start_val_for_year <- current_trend_value
  if (i %% 2 == 0) {
    end_val_for_year <- start_val_for_year + 2
    trend_component_varying_season[year_start_index:year_end_index] <- seq(from = start_val_for_year, to = end_val_for_year, length.out = frequency)
  } else {
    end_val_for_year <- start_val_for_year - 2
    trend_component_varying_season[year_start_index:year_end_index] <- seq(from = start_val_for_year, to = end_val_for_year, length.out = frequency)
  }
  current_trend_value <- end_val_for_year
}

# ノイズ成分の作成
set.seed(seed)
noise_component_varying_season <- rnorm(n_months, mean = 0, sd = 0.5)

# 全成分を合成
observed_data_varying_season <- seasonal_component_varying + trend_component_varying_season + noise_component_varying_season

# 外れ値の追加
data_ts_varying_season_outliers_val <- observed_data_varying_season
set.seed(seed)
outlier_indices <- sample(1:n_months, 3)
outlier_magnitudes <- c(10, -12, 15)
data_ts_varying_season_outliers_val[outlier_indices[1]] <- data_ts_varying_season_outliers_val[outlier_indices[1]] + outlier_magnitudes[1]
data_ts_varying_season_outliers_val[outlier_indices[2]] <- data_ts_varying_season_outliers_val[outlier_indices[2]] + outlier_magnitudes[2]
data_ts_varying_season_outliers_val[outlier_indices[3]] <- data_ts_varying_season_outliers_val[outlier_indices[3]] + outlier_magnitudes[3]

# tsオブジェクトに変換
data_ts_varying_season_outliers <- ts(data_ts_varying_season_outliers_val, frequency = frequency, start = c(2020, 1))
true_trend_ts_varying_season <- ts(trend_component_varying_season, frequency = frequency, start = c(2020, 1))

# サンプルデータのプロット
plot(data_ts_varying_season_outliers,
  main = "サンプルデータ", ylab = NULL,
  ylim = range(c(data_ts_varying_season_outliers, true_trend_ts_varying_season))
)
lines(true_trend_ts_varying_season, col = "blue", lwd = 2)
points(time(data_ts_varying_season_outliers)[outlier_indices],
  data_ts_varying_season_outliers[outlier_indices],
  col = "red", pch = 19, cex = 1.5
)
legend("topright",
  legend = c("観測データ", "真のトレンド", "外れ値"),
  col = c("black", "blue", "red"),
  lty = c(1, 1, NA), pch = c(NA, NA, 19), lwd = c(1, 2, NA), cex = 0.8
)
abline(v = seq(2020, 2020 + n_years, by = 1), col = "grey", lty = 2)

# 生成された季節成分の確認
plot(ts(seasonal_component_varying, frequency = 12, start = c(2020, 1)), main = "生成された季節成分", ylab = NULL)