Rによる時系列データのBox-Cox変換

本ポストはこちらの続きです。よって、seed値は同じく 20250517 とします。

RによるBox-Cox変換における最適ラムダの算出
対数尤度 を利用して Box-Cox変換 における 最適ラムダ(λ) の算出を R で試みます。始めにサンプルデータを作成します。なお、サンプルデータは意図的に「分散が非定常」かつ「季節性(周期性)の無い」時系列データとします。librar...
www.saecanet.com
2025.05.17

Rによる時系列データの Box-Cox変換 を試みます。

改めてサンプルデータを作成します。

seed <- 20250517
set.seed(seed)
n <- 200
t_idx <- 1:n
trend <- 50 + 0.5 * t_idx
noise_sd_ratio <- 0.15
noise <- rnorm(n, mean = 0, sd = trend * noise_sd_ratio)
x_original_values <- trend + noise
# Box-Cox変換のため、値が確実に正であることを保証させる
x_original_values[x_original_values <= 0] <- 1e-6
x_ts <- ts(x_original_values, frequency = 1) # 季節性なしのデータとして扱う

関数 forecast::BoxCox.lambda を利用して最適なラムダ(λ)を求めます。

optimal_lambda_loglik <- forecast::BoxCox.lambda(x_ts, method = "loglik")
cat(paste("Optimal lambda (loglik method):", round(optimal_lambda_loglik, 4), "\n"))
Optimal lambda (loglik method): 0.45

Box-Cox変換 のための関数を作成します。

boxcox_transform <- function(y, lambda) {
  original_attributes <- attributes(y) # tsオブジェクトの属性を保持
  y_numeric <- as.numeric(y) # 計算のために数値ベクトルに変換

  # Box-Cox変換の計算
  if (abs(lambda) < 1e-9) { # lambdaがほぼ0の場合、対数変換
    transformed_y <- log(y_numeric)
  } else { # lambdaが0でない場合
    transformed_y <- (y_numeric^lambda - 1) / lambda
  }

  # 元のオブジェクトがtsなら、tsオブジェクトとして返す
  if (is.ts(y)) {
    # tsオブジェクトの属性（開始時期、周期など）を復元
    transformed_y <- ts(transformed_y,
      start = original_attributes$tsp[1],
      frequency = original_attributes$tsp[3]
    )
  }
  return(transformed_y)
}

作成した関数を利用してサンプル時系列データを Box-Cox変換 します。

x_ts_transformed <- boxcox_transform(x_ts, optimal_lambda_loglik)
x_ts_transformed
Time Series:
Start = 1 
End = 200 
Frequency = 1 
  [1] 12.303560 10.926115  9.895295 11.212919 11.181870 12.216017 11.197002
  [8] 11.522548  9.791767  8.765698 10.754899 10.582257 12.581710 11.013034
 [15] 11.306413 12.897183 11.662039 10.142130 13.192650 13.823372 13.262558
 [22] 13.627527 10.941004 13.988421 12.158101 12.917781 11.529698 13.099224
 [29] 11.761910 13.067547 11.758352 12.437319 12.191684  9.995317 12.566340
 [36] 12.973881 13.131552 12.522335 12.217904 12.682857 11.244514 12.346608
 [43] 11.410475 12.336399 12.397760 14.005780 14.667806 14.221040 13.581120
 [50] 12.323514 12.880274 11.658165 14.131319 13.198142 14.695949 12.607398
 [57] 14.387814 14.643078 12.350418 14.021409 14.288380 13.001766 13.529733
 [64] 14.977056 12.662699 14.296725 14.844133 13.416750 13.966334 14.419613
 [71] 11.526020 12.713276 15.015883 15.381004 14.076720 15.059598 13.078936
 [78] 13.742992 14.425897 14.894895 15.390410 13.566456 13.822698 15.662622
 [85] 15.610601 13.247973 15.741638 15.145323 14.714575 13.559551 17.387243
 [92] 11.967950 15.912261 15.268174 12.837342 15.705460 15.179007 11.597971
 [99] 16.316404 15.614148 14.474377 15.414680 14.488281 16.835206 13.973556
[106] 15.465861 14.414259 17.050159 16.314383 15.462131 17.304628 15.402270
[113] 16.805953 16.112487 15.655637 17.283437 17.370022 16.310301 13.999779
[120] 14.778950 17.823199 15.346389 16.536468 15.947433 16.327611 15.625809
[127] 15.420474 15.170425 14.084622 15.534046 14.515386 17.176438 16.272183
[134] 16.005614 16.895622 14.541354 17.782362 17.948177 17.694553 19.000831
[141] 17.864323 16.295407 18.889447 14.169382 13.178991 15.869686 16.833604
[148] 17.835879 15.245122 15.529272 16.752181 17.959172 17.244867 19.088067
[155] 16.209960 16.481083 18.736243 20.842727 15.416968 16.806621 15.372214
[162] 18.262131 17.449263 17.989684 17.418979 17.370260 18.943023 14.190743
[169] 15.506217 18.062476 17.928363 17.942794 18.644096 18.570261 17.939075
[176] 16.247411 17.495822 18.345460 19.854506 19.164663 17.903038 19.684803
[183] 20.035654 18.004475 17.612414 18.847308 19.072240 18.648185 16.207550
[190] 19.855366 16.770800 19.244883 18.992126 18.892618 18.911704 18.619563
[197] 19.970070 19.078027 18.625075 19.767949

なお、Rでは forecast::BoxCox により Box-Cox変換 が可能です。

unique(x_ts_transformed - forecast::BoxCox(x_ts, lambda = optimal_lambda_loglik))
[1] 0

変換前と変換後のサンプル時系列データを比較します。

library(ggplot2)
library(zoo)
library(patchwork)
df_comparison <- data.frame(
  time = time(x_ts),
  original = coredata(x_ts),
  transformed = coredata(x_ts_transformed)
)

p_original_ts <- ggplot(df_comparison, aes(x = time, y = original)) +
  geom_line(color = "blue") +
  ggtitle("Original Time Series") +
  ylab("Value") +
  xlab("Time") +
  theme_bw()

p_transformed_ts <- ggplot(df_comparison, aes(x = time, y = transformed)) +
  geom_line(color = "red") +
  ggtitle(paste("Box-Cox Transformed Time Series (lambda =", round(optimal_lambda_loglik, 3), ")")) +
  ylab("Transformed Value") +
  xlab("Time") +
  theme_bw()

p_original_ts + p_transformed_ts + plot_layout(ncol = 1)
Figure 1

変換の前後で分散が変化していること(安定していること)を確認します。

library(dplyr)
calculate_rolling_sd <- function(series, window_size) {
  n_series <- length(series)
  rolling_sd <- rep(NA, n_series)
  if (n_series < window_size) {
    return(rolling_sd)
  }
  for (i in window_size:n_series) {
    rolling_sd[i] <- sd(series[(i - window_size + 1):i], na.rm = TRUE)
  }
  return(rolling_sd)
}

window_size <- 20 # ローリングウィンドウのサイズ (データ長の10%程度)
df_comparison$rolling_sd_original <- calculate_rolling_sd(df_comparison$original, window_size)
df_comparison$rolling_sd_transformed <- calculate_rolling_sd(df_comparison$transformed, window_size)

p_rolling_sd <- ggplot(
  df_comparison %>% filter(!is.na(rolling_sd_original)),
  aes(x = time)
) +
  geom_line(aes(y = rolling_sd_original, color = "Original"), linewidth = 1) +
  geom_line(aes(y = rolling_sd_transformed, color = "Transformed"), linewidth = 1) +
  scale_color_manual(values = c("Original" = "blue", "Transformed" = "red")) +
  ggtitle(paste("Rolling Standard Deviation (Window =", window_size, ")")) +
  ylab("Rolling SD") +
  xlab("Time") +
  labs(color = "Series") +
  theme_bw() +
  theme(legend.position = "top")

print(p_rolling_sd)
Figure 2

一定期間（ウィンドウ）ごとのデータの標準偏差（ばらつき具合）の推移を確認しますと、青線（Original）は経時に伴い分散が変化していますが、赤線（Transformed）は 青線（Original）と比較しますと低い水準で安定しており、変換によって分散が安定化したことを確認できます。

続いて、ヒストグラムによる分布の変化を確認します。

p_hist_original <- ggplot(df_comparison, aes(x = original)) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "blue", alpha = 0.7, color = "white") +
  geom_density(color = "darkblue", linewidth = 1) +
  ggtitle("Histogram of Original Data") +
  xlab("Value") +
  ylab("Density") +
  theme_bw()

p_hist_transformed <- ggplot(df_comparison, aes(x = transformed)) +
  geom_histogram(aes(y = ..density..), bins = 20, fill = "red", alpha = 0.7, color = "white") +
  geom_density(color = "darkred", linewidth = 1) +
  ggtitle(paste("Histogram of Transformed Data (lambda =", round(optimal_lambda_loglik, 3), ")")) +
  xlab("Transformed Value") +
  ylab("Density") +
  theme_bw()

plot_comparison_hist <- p_hist_original | p_hist_transformed
print(plot_comparison_hist)
Figure 3

左のヒストグラム（Original）は、右に裾の長い（歪んだ）分布を示していますが、右のヒストグラム（Transformed）では、変換により分布がより対称的になっていることを確認できます。

最後に両者の基本統計量を確認します。

summary_original <- summary(df_comparison$original)
summary_transformed <- summary(df_comparison$transformed)

cat("--- 要約統計量の比較 ---\n")
cat("Original Data Summary:\n")
print(summary_original)
cat("\nTransformed Data Summary:\n")
print(summary_transformed)
cat(paste("SD Original:", round(sd(df_comparison$original), 2), "  SD Transformed:", round(sd(df_comparison$transformed), 2), "\n"))
cat(paste("IQR Original:", round(IQR(df_comparison$original), 2), " IQR Transformed:", round(IQR(df_comparison$transformed), 2), "\n"))
--- 要約統計量の比較 ---
Original Data Summary:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  34.87   72.76   96.61   98.63  124.44  181.19 

Transformed Data Summary:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  8.766  13.076  15.158  15.089  17.255  20.843 
SD Original: 31.98   SD Transformed: 2.58 
IQR Original: 51.68  IQR Transformed: 4.18

特に標準偏差 (sd) や四分位範囲 (IQR = Q3 – Q1) に注目すると、変換後のデータの方が値のスケールが小さくなり、ばらつきが 抑えられていることが確認できます。

ただし、スケールが異なっているのですから(意図的に変換しているのですから)、単純なSDの比較は注意が必要 です。

時系列データの Box-Cox変換 の主目的は「時間を通じたばらつきの均一化」です。

以上です。