Rで解析：時系列のトレンドを線形フィルタで取り出す「rjd3filters」パッケージ

「rjd3filters」パッケージは、時系列データから線形フィルタでトレンドと循環成分を取り出すパッケージです。ヘンダーソン移動平均などのフィルタを作り、その係数やゲインを見比べてから、手元の時系列に当てることができます。系列の端、つまり先の観測値がまだ無い直近の時点でも、遅れを抑えたフィルタでトレンドを推計できます。渡すのは ts 形式の時系列で、季節調整済みの月次・四半期データを想定しています。JDemetra+ 3.x のJavaライブラリを使うため、Java 21 以上が必要です。景気や需要の転換点を早めに捉えたいときに役立てられるのではないかと考えます。

パッケージバージョンは2.4.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

RとPythonで学ぶ[実践的]データサイエンス&機械学習【増補改訂版】 | 有賀友紀, 大橋俊介

RとPythonを併走させながらデータサイエンスと機械学習を学ぶ。両方の言語を行き来する人に向く。

パッケージのインストール

下記コマンドを実行してください。事前に Java 21 以上をインストールしておく必要があります。

# パッケージのインストール
install.packages("rjd3filters")


# パッケージの読み込み
library("rjd3filters")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

線形フィルタ：前後の観測値に係数を掛けて足し合わせ、その時点の値を推計する方法。移動平均がその代表で、係数の並びがフィルタの性質を決めます。前後の観測値を同じ数ずつ使うものを対称フィルタ、先の観測値が足りない系列の端で使うものを非対称フィルタと呼びます。

ゲイン：周期ごとに、フィルタを通した後に変動がどれだけ残るかを表す値。トレンドを取り出すフィルタは、長い周期では1に近く、短い周期（ノイズ）では0に近いほど良いとされます。

以下では、京丹後市の月別宿泊者数指数（季節調整済みを想定した架空データ、2016年1月〜2025年12月）を使います。

# 乱数の種を固定する
set.seed(2026)


# 月数を指定する
n <- 120


# 経過月を作る
month_no <- 1:n


# 緩やかな増加と約4年周期の波に、月ごとのゆらぎを足す
tango <- 100 + 0.15 * month_no + 6 * sin(2 * pi * month_no / 48) + rnorm(n, sd = 2.5)


# 月次の時系列にする
tango <- ts(tango, start = c(2016, 1), frequency = 12)


# 先頭を確認する
head(tango)

単純移動平均の作成：simple_maコマンド

係数がすべて同じ移動平均を作ります。作った移動平均どうしは「+」「*」「/」「^」で組み合わせられるので、2×12項の中心化移動平均や3×3項の移動平均もそのまま作れます。

# 3項の単純移動平均を作る
ma3 <- simple_ma(3, -1)
ma3

# 12項の移動平均を2つずらして平均し、2×12項の中心化移動平均を作る
m2x12 <- (simple_ma(12, -6) + simple_ma(12, -5)) / 2
m2x12

# 3項の移動平均を2回重ねて、3×3項の移動平均を作る
ma3^2

局所多項式フィルタの作成：lp_filterコマンド

局所多項式の当てはめに基づくフィルタを作ります。対称フィルタに加えて、先の観測値が0〜horizon個しか無いときの非対称フィルタもまとめて返します。kernel = “Henderson” で、季節調整で使われるヘンダーソン移動平均になります。

# ヘンダーソン型の13項（horizon = 6）の局所多項式フィルタを作る
henderson <- lp_filter(horizon = 6, kernel = "Henderson", endpoints = "LC")
henderson

＜おすすめのRに関する書籍です＞

Rで学ぶプログラミングの基礎の基礎 | 舟尾暢男

統計以前のプログラミングの考え方から始めるRの入門書。つまずきやすい部分を丁寧に扱う。

フィルタ係数のプロット：plot_coefコマンド

フィルタの係数をプロットします。q で、先の観測値をいくつ使えるときのフィルタを並べるかを選びます。

# 先の観測値が0・3・6個のときの係数を並べてプロットする
plot_coef(henderson, q = c(0, 3, 6), legend = TRUE, lwd = 2)

ゲインのプロット：plot_gainコマンド

フィルタのゲインを周波数ごとにプロットします。横軸は0〜πの周波数で、左ほど長い周期です。

# 直近の時点（q = 0）と対称フィルタ（q = 6）のゲインを比べる
plot_gain(henderson, q = c(0, 6), legend = TRUE, lwd = 2)

時系列データのフィルタリング：filterコマンド

作ったフィルタを時系列に当てます。非対称フィルタもまとめて渡すと、系列の端まで値が埋まります。読み込むと stats パッケージの filter コマンドを上書きする点に注意してください。

# ヘンダーソン移動平均で宿泊者数指数のトレンドを取り出す
tango_trend <- filter(tango, henderson)

# 元の指数とトレンドを重ねてプロットする
plot(ts.union(tango, tango_trend),
  plot.type = "single",
  col = c("gray50", "red"), lwd = 2, ylab = "宿泊者数指数"
)

＜おすすめのRに関する書籍です＞

R／RStudioでやさしく学ぶプログラミングとデータ分析 【ダウンロード特典付き！】

プログラミング未経験からでも読める入門書。RとRStudioの操作を一つずつ確認しながら分析の初歩まで進む。

www.amazon.co.jp

RKHSフィルタの作成：rkhs_filterコマンド

再生核ヒルベルト空間（RKHS）の考え方でフィルタを作ります。系列の端で使う非対称フィルタの幅を、asymmetricCriterion で選んだ基準が最も良くなるように決めます。

# 遅れ（Timeliness）を基準にRKHSフィルタを作る
rkhs <- rkhs_filter(horizon = 6, asymmetricCriterion = "Timeliness")

# ヘンダーソン移動平均と同じ q で係数をプロットして見比べる
plot_coef(rkhs, q = c(0, 3, 6), legend = TRUE, lwd = 2)

FSTフィルタの作成：fst_filterコマンド

忠実度（Fidelity）・滑らかさ（Smoothness）・遅れ（Timeliness）の3つの基準の重みを決めて、フィルタを1本作ります。lags と leads で、使う過去と先の観測値の数を決めます。

# 先の観測値を2つ使い、滑らかさと遅れを0.3ずつ重視したFSTフィルタを作る
fst2 <- fst_filter(lags = 6, leads = 2, smoothness.weight = 0.3, timeliness.weight = 0.3)
fst2

# 同じく先の観測値を2つ使うヘンダーソン移動平均と、トレンドの推計を比べる
fst_ma <- filter(tango, fst2)
lpp_ma <- filter(tango, henderson[, "q=2"])

# 元の指数と2つの推計を重ねてプロットする
plot(ts.union(tango, fst_ma, lpp_ma),
  plot.type = "single",
  col = c("gray50", "red", "blue"), lwd = 2, ylab = "宿泊者数指数"
)

トレンドの信頼区間の算出：confint_filterコマンド

フィルタで取り出したトレンドに信頼区間を付けます。結果は推定値・下限・上限の3列の時系列です。

# ヘンダーソン移動平均で取り出したトレンドに90%信頼区間を付ける
tango_ci <- confint_filter(tango, henderson, level = 0.9)

# 直近の6か月を確認する
tail(tango_ci, 6)

# トレンドを実線、信頼区間を破線でプロットする
plot(tango_ci,
  plot.type = "single",
  col = c("red", "black", "black"), lty = c(1, 2, 2), lwd = 2,
  ylab = "宿泊者数指数"
)

＜おすすめのRに関する書籍です＞

Rによるスクレイピング入門 | 石田基広, 市川太祐, 瓜生真也, 湯谷啓明

Web上の情報をRで収集する方法を、HTMLやAPIの基礎から解説。例題を動かしながら覚えられる。

この記事が誰かの役に立ちますように。