はじめに

こんにちは、事業会社で働いているデータサイエンティストです。

皆さんは業務で「時系列データの分類」をお願いされることがありますか？
実は、ビジネスの現場では、時系列を分類しなければならない場面が少なくありません。

たとえば、投資家とのコミュニケーションの中で、「御社の業績と最も相関の強いマーケット指標は何ですか？」と聞かれることはよくあります。このような場面では、自社の売上などのKGI・KPIと、各マーケット指標の時系列を分類・比較することで、「自社の業績と同じ動きをする指標」を特定できます。これにより、投資家の疑問に明確に答えられ、信頼を得ることにもつながります。

また、多数の事業を展開している企業であれば、まず地域ごと・事業ごとの売上データを分類することで、事業環境の見通しが明瞭になります。たとえば、新宿オフィスと同じトレンドを持つ別の拠点があれば、「同様のビジネス環境である可能性が高い」と仮定できます。その結果、拡張時には新宿のメンバーを優先的にアサインすることで、移行がスムーズになるかもしれません。

最後の例に入る前に、まず強調しておきたいのは、本記事は特定の手法を使って有価証券への投資を促すものではなく、読者の資産配分に対して責任を負うものでもありません。しかし、株価の時系列を分類することで、以下のような応用が可能です。まず、同じグループに分類された銘柄だけに集中投資せず、異なるグループへ分散投資することでリスクを抑える。また、逆に、あるグループに成長余地があると判断した場合、そのグループの銘柄に重点投資する。つまり、「多数の時系列を分類すること」は、より良い判断や戦略に直結する技術なのです。

時系列データを「分類する」という行為は、単なる分析手法ではなく、ビジネスの現場で具体的な意思決定に直結する力を持っています。今回の記事では、ビジネスの現場で直面するような、大量の時系列を分類するベイズモデルを紹介します。

モデル説明

本記事では、ディリクレ過程（Dirichlet Process）を用いた時系列モデルによって、時系列データを分類する手法を紹介します。

時系列データの分類においては、「いくつのグループに分ければよいのか？」という問いが避けて通れません。しかし、現実のビジネスや社会現象のデータにおいて、あらかじめ適切なグループ数を知っていることは稀です。むしろ、分析の目的そのものが「見えざる構造を明らかにすること」にあることも多いでしょう。

そこで有力なのが、ディリクレ過程を基礎としたベイズ非パラメトリック手法です。ディリクレ過程を用いることで、グループ数そのものを推定の対象とすることができ、分析者が恣意的にグループ数を指定する必要がありません。この特性により、仮定を最小限に抑えたまま柔軟な分類が可能になります。

さらに、ディリクレ過程は、観測データの構造に応じて自然にグループ数を調整する性質を持つため、モデルの汎化性能や客観性が高いという利点もあります。分析者の主観に頼らず、データそのものに語らせる。

では、モデルの説明に入ります。

まず、ハイパーパラメータ$\alpha$をガンマ分布からサンプリングします：

$$
\alpha \sim Gamma(0.001, 0.001)
$$

次に、棒折り過程を用いてグループpの所属確率$p_{p}$をサンプリングします。グループは無限にある設定なので、理論上この操作は無限回繰り返します：

$$
\pi_{p} \sim Beta(1, \alpha)
$$

$$
p_{p} = \pi_{p} \prod_{l=1}^{p - 1} (1 - \pi_{l})
$$

要するに、全ての時系列の中の$p_{p}$の割合の時系列が、グループpに属するということです。

各グループはトレンド$\mu_{p}$をもちます：

$$
\phi_{begin, p} \sim Gamma(0.001, 0.001)
$$

$$
\phi_{shift, p} \sim Gamma(0.001, 0.001)
$$

$$
\mu_{p, 1} \sim Normal(0, \phi_{begin, p})
$$

$$
\mu_{p, t} \sim Normal(\mu_{p, t-1}, \phi_{shift, p}), \space \text{for t > 1}
$$

言葉で説明すると、各グループが持つトレンドは、一次の自己回帰モデルの形で表現できるということです。

最後に、時系列$s$のグループ番号$eta_{s}$を$p$からサンプリングします：

$$
\eta_{s} \sim Categorical(p)
$$

続いて、時系列は$\eta_{s}$グループのトレンドを平均とする独立した正規分布からサンプリングされます：

$$
Series_{s} \sim Normal(\mu_{\eta_{s}}, \sigma)
$$

モデルに含まれるグループ番号パラメータ $\eta$ は離散的な変数であるため、Stanのような一般的なベイズモデリング言語では直接扱うことができません。そのため、モデリングの際には $\eta$ を積分によって除去し、事後分布の推定後に generated quantities ブロックで再サンプリングするというアプローチを取る必要があります。また、理論上はグループ数は無限に存在し得ますが、現実のコンピュータは無限長の配列を扱うことができません。そのため、実装上は「無限大」をあらかじめ十分に大きな有限の数値で近似し、計算可能な範囲に収める必要があります。

モデルの実装コードはこちらです：

classify_ts.stan

functions {
  vector stick_breaking(vector breaks){
    int length = size(breaks) + 1;
    vector[length] result;
    
    result[1] = breaks[1];
    real summed = result[1];
    for (d in 2:(length - 1)) {
      result[d] = (1 - summed) * breaks[d];
      summed += result[d];
    }
    result[length] = 1 - summed;
    
    return result;
  }
}
data {
  int date_type;
  int variable_type;
  int group_type;
  
  array[variable_type] int variable_begin;
  array[variable_type] int variable_end;
  
  int N;
  array[N] int date;
  vector[N] value;
}
parameters {
  real<lower=0> group_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[group_type - 1] group_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  vector<lower=0>[group_type] phi_begin;
  vector<lower=0>[group_type] phi_shift;
  real<lower=0> sigma;
  
  array[group_type] vector[date_type] mu;
}
transformed parameters {
  simplex[group_type] group;
  
  group = stick_breaking(group_breaks);
}
model {
  group_alpha ~ gamma(0.001, 0.001);
  
  group_breaks ~ beta(1, group_alpha);
  
  phi_begin ~ gamma(0.001, 0.001);
  phi_shift ~ gamma(0.001, 0.001);
  sigma ~ gamma(0.001, 0.001);
  
  for (i in 1:group_type){
    mu[i, 1] ~ normal(0, phi_begin[i]);
    for (j in 2:date_type){
      mu[i, j] ~ normal(mu[i, j - 1], phi_shift[i]);
    }
  }
  
  for (i in 1:variable_type){
    vector[group_type] case_when = log(group);
    for (g in 1:group_type){
      case_when[g] += normal_lupdf(value[variable_begin[i]:variable_end[i]] | mu[g, date[variable_begin[i]:variable_end[i]]], sigma);
    }
    target += log_sum_exp(case_when);
  }
}
generated quantities {
  array[variable_type] vector[group_type] eta;
  for (i in 1:variable_type){
    eta[i] = rep_vector(0.0, group_type);
    vector[group_type] case_when = log(group);
    for (g in 1:group_type){
      case_when[g] += normal_lpdf(value[variable_begin[i]:variable_end[i]] | mu[g, date[variable_begin[i]:variable_end[i]]], sigma);
    }
    eta[i, categorical_rng(softmax(case_when))] = 1.0;
  }
}

データ推定

ここからは、実際のデータ分析に進みます。使用するのは、アメリカのマクロ経済に関連する時系列データです。このデータには一部欠損値が含まれていますが、本モデルの設計上、問題なく取り扱うことが可能です。

というのも、グループ番号$\eta$が与えられると、各時系列の値は条件付きで独立に分布しているため、欠損している部分をそのまま無視しても、観測されたデータのみを用いて尤度関数を構築することができます。これにより、欠損値の補完を事前に行わなくても、自然な形でモデリングと推定が進められます。

まずは標準化とマスター表の作成などの前処理を実施します：

data("fred_md", package = "BVAR")

ts_df <- fred_md |>
  tibble::rownames_to_column(var = "date") |>
  tibble::tibble() |>
  tidyr::pivot_longer(!date, names_to = "variable", values_to = "value") |>
  tidyr::drop_na() |>
  dplyr::group_by(variable) |>
  dplyr::mutate(
    value = (value - mean(value))/sd(value),
    date = as.Date(date)
  ) |>
  dplyr::ungroup() |>
  dplyr::arrange(variable, date)

date_master <- ts_df |>
  dplyr::select(date) |>
  dplyr::distinct() |>
  dplyr::arrange(date) |>
  dplyr::mutate(
    date_id = dplyr::row_number()
  )

variable_master <- ts_df |>
  dplyr::select(variable) |>
  dplyr::distinct() |>
  dplyr::arrange(variable) |>
  dplyr::mutate(
    variable_id = dplyr::row_number()
  )

# このbegin_end_dfの技法はまた別記事で紹介します！
begin_end_df <- ts_df |>
  dplyr::mutate(
    rid = dplyr::row_number()
  ) |>
  dplyr::summarise(
    begin = min(rid),
    end = max(rid),
    .by = variable
  )

ts_df_with_id <- ts_df |>
  dplyr::left_join(
    date_master, by = "date"
  ) |>
  dplyr::left_join(
    variable_master, by = "variable"
  )

data_list <- list(
  date_type = nrow(date_master),
  variable_type = nrow(variable_master),

  # 無限大を10で近似
  group_type = 10,
  
  variable_begin = begin_end_df$begin,
  variable_end = begin_end_df$end,
  
  N = nrow(ts_df_with_id),
  date = ts_df_with_id$date_id,
  value = ts_df_with_id$value
)

続いて、モデルのコンパイルを実施し：

m_cts_init <- cmdstanr::cmdstan_model("classify_ts.stan")

変分推論を用いたモデル推定を実施します：

> m_cts_estimate <- m_cts_init$variational(
     seed = 1234567,
     data = data_list
 )
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.009975 seconds 
1000 transitions using 10 leapfrog steps per transition would take 99.75 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Success! Found best value [eta = 1] earlier than expected. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100       -96721.959             1.000            1.000 
   200       -93816.248             0.515            1.000 
   300       -92887.571             0.347            0.031 
   400       -92643.757             0.261            0.031 
   500       -92594.586             0.209            0.010   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  15.8 seconds.

データ量は90010とかなり多いですが、推定自体は15秒で終わりました。最後に、結果を整形して保存します：

m_cts_summary <- m_cts_estimate$summary()

推定結果可視化

まず、各グループの所属確率を確認しましょう：

> m_cts_summary |>
     dplyr::filter(stringr::str_detect(variable, "^group\\["))
# A tibble: 10 × 7
   variable       mean     median       sd        mad           q5      q95
   <chr>         <dbl>      <dbl>    <dbl>      <dbl>        <dbl>    <dbl>
 1 group[1]  0.756     0.760      0.0426   0.0414     0.682        0.821   
 2 group[2]  0.00687   0.00430    0.00869  0.00387    0.000798     0.0214  
 3 group[3]  0.00694   0.00398    0.00917  0.00390    0.000672     0.0221  
 4 group[4]  0.224     0.220      0.0415   0.0412     0.161        0.296   
 5 group[5]  0.00314   0.00108    0.00651  0.00135    0.0000496    0.0130  
 6 group[6]  0.00139   0.000327   0.00345  0.000459   0.00000508   0.00609 
 7 group[7]  0.000589  0.0000783  0.00168  0.000113   0.000000822  0.00292 
 8 group[8]  0.000235  0.0000217  0.00100  0.0000315  0.000000165  0.000938
 9 group[9]  0.0000921 0.00000690 0.000333 0.0000101  0.0000000303 0.000426
10 group[10] 0.0000785 0.00000398 0.000346 0.00000584 0.0000000152 0.000308

基本的にグループ1（75.6%）とグループ4（22.4%）以外のグループの所属確立はゼロに近いので、これからの可視化では無視します。

ここで、ggplot2でグループ1とグループ4のトレンドを描画します。

g_trend <- m_cts_summary |>
  dplyr::filter(stringr::str_detect(variable, "^mu\\[")) |>
  dplyr::mutate(
    id = variable |>
      purrr::map(
        \(x){
          as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
        }
      )
  ) |>
  tidyr::unnest_wider(id, names_sep = "_") |>
  dplyr::filter(id_1 %in% c(1, 4)) |>
  dplyr::left_join(
    date_master, by = c("id_2" = "date_id")
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean, group = as.factor(id_1), color = as.factor(id_1))) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95, group = as.factor(id_1), fill = as.factor(id_1)), alpha = 0.3) +
  ggplot2::scale_x_date(
    breaks = seq(min(date_master$date), max(date_master$date), by = "5 year"),
    date_labels = "%Y"
  ) + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(
    title = "trend",
    color = "group", 
    fill = "group"
  )

続いて、グループ1とグループ4に分類された時系列の図をそれぞれ描画します：

eta <- m_cts_summary |>
  dplyr::filter(stringr::str_detect(variable, "^eta\\[")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10)

g_group_1 <- ts_df_with_id |>
  dplyr::left_join(
    variable_master |>
      dplyr::mutate(
        id = .env$eta |>
          nrow() |>
          seq_len() |>
          purrr::map_int(
            \(i){
              which.max(eta[i,])
            }
          )
      ) |>
      dplyr::select(variable_id, id),
    by = "variable_id"
  ) |>
  dplyr::filter(id == 1) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = value, color = variable), alpha = 0.2) + 
  ggplot2::scale_x_date(
    breaks = seq(min(date_master$date), max(date_master$date), by = "5 year"),
    date_labels = "%Y"
  ) + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(
    title = "group 1"
  ) +
  ggplot2::theme(legend.position = "none")


g_group_4 <- ts_df_with_id |>
  dplyr::left_join(
    variable_master |>
      dplyr::mutate(
        id = .env$eta |>
          nrow() |>
          seq_len() |>
          purrr::map_int(
            \(i){
              which.max(eta[i,])
            }
          )
      ) |>
      dplyr::select(variable_id, id),
    by = "variable_id"
  ) |>
  dplyr::filter(id == 4) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = value, color = variable), alpha = 0.2) + 
  ggplot2::scale_x_date(
    breaks = seq(min(date_master$date), max(date_master$date), by = "5 year"),
    date_labels = "%Y"
  ) + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(
    title = "group 4"
  ) +
  ggplot2::theme(legend.position = "none")

最後に、patchworkを使って、3枚の画像をきれいに整形して出します：

library(patchwork)

(g_group_1/g_group_4) | g_trend

まず、右側の図に示されているのは、グループ1とグループ4の時系列トレンドです。注目すべき点として、グループ1のトレンドは一貫して上昇傾向にあり、非常に安定しています。一方で、グループ4のトレンドは変動が激しく、とくに直近ではリーマンショック後に大きく下落している様子が見て取れます。

この結果から推察できるのは、グループ1に属する時系列は、アメリカ経済の長期的・構造的な潜在力を反映している可能性が高いということです。それに対して、グループ4の時系列は、短期的なショックや外部要因によって大きく影響を受けやすい指標で構成されていると考えられます。

左側の2枚の図は、それぞれグループ1とグループ4に分類された時系列の推移を示しています。図を比較すると明らかなように、グループ1に属する時系列は全体として安定した動きを見せている一方で、グループ4に属する時系列は変動が非常に大きく、不安定な傾向が見られます。

特に注目すべきは、リーマンショック後の動きです。グループ4の時系列はこのタイミングで一斉に大きく下落しており、外部ショックに対して敏感に反応していることがわかります。

これは抽出されたトレンドと一致した結果です。

推定結果の解釈

ここでは分類結果の詳細を確認します。

グループ1には、アメリカ経済の実体的な動向を表す指標が多く含まれています。例えば以下の指標がグループ1に分類されました：

雇用関連指標：非農業部門雇用者数（PAYEMS）、労働力人口（CLF16OV）、週あたり労働時間（AWHMAN）など
産業生産・製造業指標：鉱工業生産指数（INDPRO）、各種工業生産（IP）コード（IPMANSICS, IPFINALなど）
消費関連指標：小売売上（RETAILx）、消費者支出（PCEPIなど）
インフレ・物価指標：CPI（CPIAUCSL, CPIMEDSLなど）、PPI（PPICMM）など。
貨幣・信用：マネーサプライ（M1SL, M2SL）、消費者信用（NONREVSL）など
国際収支・為替：輸出入関連（USTPU, USGOODなど）、為替レート（EXCAUSxなど）

このクラスターの構成から、モデルは「基礎的な経済活動」や「マクロ経済の循環要因（景気循環）」を一つの構造として抽出していると考えられます。つまり、GDPや雇用、消費、インフレといった、FRBやエコノミストが注視する中核的なマクロ変数を捉える共通のパターンがあることを示唆しています。

グループ4には、景気よりも政策金利や金融市場の動向、あるいは構造的トレンドや価格の調整速度が速いセクターに敏感な指標が集まっています。特徴的なカテゴリは以下の通りです：

金利・金融市場指標：FF金利（FEDFUNDS）、国債利回り（GS1, GS5, GS10）、短期金融（TB3MS, TB6MS）など
為替：ドル円（EXJPUSx）、ドルポンド（EXUSUKx）など
雇用の特定サブグループ：製造業雇用（MANEMP, DMANEMP）
住宅市場：住宅着工（HOUST）、地域別住宅建設許可（PERMITNEなど）
センチメントと信用リスク：消費者心理（UMCSENTx）、短期資金の金利（COMPAPFFx）など

このクラスターは、より金融政策やリスクセンチメントの変化に対して反応しやすいセクターや価格指標が集約されています。景気の「ストック」的な側面よりも、流動性や信用状況、金融の「ショック伝播」の経路を捉えています。

結論

いかがでしたか？

最後に、クラスタリングや分類モデルをビジネスやアカデミアの現場で用いる際に、しばしば直面する懐疑的な態度について触れておきます。

クラスタリングを回すと、決まって「それ、知ってたわ」「直感的だね」と言う方が現れます。しかし、そのような発言は、本当に意味のある直感だったのでしょうか？推定結果を見た「後で」なら、誰でもそれらしく語ることはできます。ですが、「事前に」それを体系的に導き出すことができたかどうかは、まったく別問題です。

そもそも、「知ってたわ」と感じるということ自体が、むしろモデルが正しく、自然な構造を捉えていることの証拠なのです。よくできたクラスタリングが直感にフィットするのは、それが意味のある構造を可視化しているからに他なりません。

にもかかわらず、「モデルを使わなくても同じことが分かる」と主張するのは、結果論でしかなく、科学的な態度とは言えません。データの構造を明示的に抽出し、再現可能な形で提示するという意味で、モデルには計り知れない価値があります。

直感を補強し、曖昧だった知識を形にするのがモデルの力です。「知ってたわ」で片付けてしまう人は、その貴重なプロセスを無視し、後知恵で世界を語ろうとする危うさに気づくべきです。データに向き合い、そこから新たな洞察を得ることこそが、我々がモデルを使う理由であり、真の知的営みなのです。

私たちと一緒にデータサイエンスの力を使って社会を改善していきたい方はぜひこちらのページをご確認ください：