Rのパッケージ:KFAS

1. サンプルデータの作成

まず、シミュレーションのためのサンプルデータを作成します。 underlying（潜在的な）状態 alpha がランダムウォーク（正規分布に従う誤差項が加わっていくモデル）で変動し、観測値 y がその状態 alpha を指数変換したものをパラメータ（期待値）に持つポアソン分布から生成されるものとします。

• 状態方程式: α_{t+1} = α_t + η_t, ここで η_t ~ N(0, Q)
• 観測方程式: y_t | α_t ~ Poisson(exp(α_t))

# 必要なパッケージを読み込みます
library(KFAS)
library(ggplot2)
library(dplyr)

# 乱数のシードを固定して再現性を確保します
seed <- 20251007
set.seed(seed)

# 時系列の長さを設定します
n <- 150

# 状態の分散（真の値）を設定します
# この値が大きいほど、状態（潜在トレンド）は大きく変動します
true_q <- 0.01

# 状態ベクトルαを格納する変数を用意します
alpha <- numeric(n)

# 状態の初期値を設定します
alpha[1] <- 1.5

# 状態方程式に従って、潜在的な状態alphaを生成します
# α_t = α_{t-1} + η_t, η_t ~ N(0, true_q)
for (t in 2:n) {
  alpha[t] <- alpha[t - 1] + rnorm(1, mean = 0, sd = sqrt(true_q))
}

# 観測方程式に従って、観測値yを生成します
# y_t ~ Poisson(lambda = exp(alpha_t))
y <- rpois(n, lambda = exp(alpha))

# 生成したデータをデータフレームにまとめます
sim_data <- data.frame(
  time = 1:n,
  y = y,
  true_alpha = alpha
)

# サンプルデータの可視化
sim_data %>%
  tidyr::gather(, , colnames(.)[-1]) %>%
  ggplot(mapping = aes(x = time, y = value, colour = key)) +
  geom_line() +
  geom_point() +
  theme_bw() +
  facet_wrap(. ~ key, scales = "free_y") +
  theme(legend.position = "none") +
  labs(title = "生成したサンプルデータ")

2. KFASによる状態空間モデルの構築と推定

次に、作成したサンプルデータ y を用いて、KFASで状態空間モデルを構築し、パラメータと潜在状態を推定します。

ステップ1: モデルの定義

SSModel関数を使ってモデルを定義します。観測値 y がポアソン分布に従うローカルレベルモデル（ランダムウォーク・トレンドモデル）を仮定します。

• y ~ SSMtrend(1, Q = NA): トレンド成分を定義します。degree = 1はローカルレベルモデルを意味します。状態の分散Qは未知ですので、NAとして推定対象にします。
• distribution = "poisson": 観測値がポアソン分布に従うことを指定します。

ステップ2: パラメータの最尤推定 fitSSM関数を用いて、モデルの未知パラメータ（この場合はQ）を最尤法で推定します。initsでパラメータの初期値を設定します。

ステップ3: フィルタリングと平滑化 KFS関数を使い、推定されたパラメータを持つモデルに対してカルマンフィルタと状態平滑化を実行します。これにより、各時点における状態 alpha の推定値が得られます。

# ポアソン分布を仮定したローカルレベルモデルを定義します
# 状態の分散Qは未知(NA)とし、データから推定させます
model <- SSModel(y ~ SSMtrend(1, Q = list(matrix(NA))),
  distribution = "poisson"
)

cat("--- モデルの構造 ---\n")
print(model)

--- モデルの構造 ---
Call:
SSModel(formula = y ~ SSMtrend(1, Q = list(matrix(NA))), distribution = "poisson")

State space model object of class SSModel

Dimensions:
[1] Number of time points: 150
[1] Number of time series: 1
[1] Number of disturbances: 1
[1] Number of states: 1
Names of the states:
[1]  level
Distributions of the time series:
[1]  poisson

Object is a valid object of class SSModel.

# fitSSM関数でモデルの未知パラメータ(log(Q))を最尤推定します
# initsにはlog(Q)の初期値を設定します
fit <- fitSSM(model, inits = log(0.01), method = "BFGS")

cat("--- パラメータの最尤推定 ---\n")
cat("- 推定された状態の分散Q:", exp(fit$optim.out$par), "\n")
cat("- 真の値Q:", true_q)

--- パラメータの最尤推定 ---
- 推定された状態の分散Q: 0.007581761 
- 真の値Q: 0.01

# KFS関数でフィルタリングと平滑化を実行します
# smoothing = "state"で状態の平滑化、"mean"で観測値の期待値の平滑化を行います
kfs_results <- KFS(fit$model, smoothing = c("state", "mean"))
cat("--- 平滑化の結果 ---\n")
print(kfs_results)

--- 平滑化の結果 ---
Smoothed values of states and standard errors at time n = 150:
       Estimate  Std. Error
level  1.1996    0.2101    

# 平滑化された状態（α_tの推定値）を取得します
smoothed_alpha <- kfs_results$alphahat

# 平滑化された期待値（exp(α_t)の推定値）を取得します
smoothed_lambda <- fitted(kfs_results)

# 結果をデータフレームに追加します
sim_data$smoothed_alpha <- as.numeric(smoothed_alpha)
sim_data$smoothed_lambda <- as.numeric(smoothed_lambda)

cat("--- シミュレーション結果の一部を確認 ---\n")
str(sim_data)

--- シミュレーション結果の一部を確認 ---
'data.frame':   150 obs. of  5 variables:
 $ time           : int  1 2 3 4 5 6 7 8 9 10 ...
 $ y              : int  4 5 9 3 1 3 5 2 5 6 ...
 $ true_alpha     : num  1.5 1.44 1.42 1.49 1.67 ...
 $ smoothed_alpha : num  1.47 1.47 1.47 1.43 1.4 ...
 $ smoothed_lambda: num  4.34 4.35 4.34 4.18 4.06 ...

3. 結果の可視化と解釈

最後に、シミュレーションの結果を可視化して、モデルが潜在的な状態を推定できているか確認します。

# ggplot2を使用して結果をプロットします
g <- ggplot(sim_data, aes(x = time)) +
  geom_line(aes(y = y, color = "観測値 (y)"), linewidth = 0.1) +
  geom_line(aes(y = smoothed_lambda, color = "平滑化された期待値 (KFAS)"), linewidth = 1.2) +
  geom_line(aes(y = exp(true_alpha), color = "真の期待値"), linewidth = 1.2) +
  labs(
    title = "KFASによる非ガウス状態空間モデルのシミュレーション結果",
    subtitle = "観測値(ポアソン分布)と潜在トレンドの推定",
    x = "時間",
    y = "値",
    color = "系列"
  ) +
  scale_color_manual(values = c(
    "観測値 (y)" = "black",
    "平滑化された期待値 (KFAS)" = "blue",
    "真の期待値" = "red"
  )) +
  theme_bw() +
  theme(legend.position = "bottom")

print(g)

Figure 2 には、以下の3つの要素がプロットされています。

• 観測値 (y): 黒い実線で示された、ポアソン分布から生成したサンプルのカウントデータです。
• 平滑化された期待値 (Smoothed Lambda): 青い実線で示された、KFASによって推定された各時点でのポアソン分布の期待値（λ）です。観測値のばらつきの裏にあるトレンドになります。
• 真の期待値 (True Lambda): 赤い実線で示された、シミュレーションで最初に設定した真の期待値 exp(alpha) です。

4. SSModelの引数

args(SSModel)

function (formula, data, H, u, distribution, tol = .Machine$double.eps^0.5) 
NULL

SSModel関数は、状態空間モデルの構造を定義するための関数です。以下にその引数を確認します。

• formula モデルの構造を定義する数式オブジェクトです。左辺に観測値ベクトル（時系列データ）を、右辺にモデルの成分を指定します。右辺には、SSMtrend（トレンド成分）、SSMseasonal（季節成分）、SSMcycle（周期成分）といったKFASパッケージが提供する関数や、通常の回帰モデルのような説明変数を記述できます。今回のシミュレーションではy ~ SSMtrend(1, Q = list(matrix(NA)))と指定しました。

• data formulaで指定した変数が含まれるデータフレームやリスト、環境を指定します。

• H 観測誤差の分散共分散行列を指定します。これは観測方程式 y_t = Z_t * α_t + ε_t における誤差項 ε_t の分散にあたります。distributionが "gaussian"（正規分布）の場合に主に用いられます。今回のシミュレーションのような非ガウスモデルでは、分散は状態に依存して決まるため、通常この引数は使用しません。NAを指定すると、この分散パラメータが最尤法による推定対象となります。

• u 観測値が従う確率分布の追加パラメータを指定します。例えば、distributionが "binomial"（二項分布）の場合の試行回数や、"negative binomial"（負の二項分布）の場合の形状パラメータなどが該当します。ポアソン分布や正規分布のように追加パラメータが不要な場合は指定する必要はありません。

• distribution 観測値 y が従う確率分布を指定する文字列です。デフォルトは "gaussian" ですので、正規分布を仮定する線形ガウス状態空間モデルを扱います。ポアソン分布 ("poisson")、二項分布 ("binomial")、ガンマ分布 ("gamma") などを指定することで、非ガウス状態空間モデルを定義できます。

• tol 数値計算上の許容誤差を指定する微小な値です。モデルの行列計算などにおいて、ゼロと見なす閾値として使用されます。

5. SSMtrendの引数

args(SSMtrend)

function (degree = 1, Q, type, index, a1, P1, P1inf, n = 1, state_names = NULL, 
    ynames) 
NULL

SSMtrend関数は、SSModelのformula内で使用され、モデルにトレンド成分を追加するための関数です。ローカルレベルモデルやローカル線形トレンドモデルなどのトレンドを定義できます。

• degree トレンドの階数を整数で指定します。

  • degree = 1（デフォルト）: ローカルレベルモデル（stochastic level model）を定義します。状態はレベル（水準）のみで、レベルがランダムウォークに従って変動するモデルです。状態方程式は level_{t+1} = level_t + η_t となります。
  • degree = 2: ローカル線形トレンドモデル（local linear trend model）を定義します。状態はレベル（水準）とスロープ（傾き）の2つで、レベルとスロープの両方が時間と共に確率的に変動するモデルです。

• Q 状態方程式における誤差項（η_t）の分散共分散行列を指定します。時系列が1つの場合は、list(matrix(Q_value)) のようにリストで渡します。degree=1なら1×1、degree=2なら最大2×2の行列となります。分散の値を未知として最尤法で推定したい場合は、NAを指定します（例: Q = list(matrix(NA))）。

• type degreeの代わりに使用できる、トレンドタイプを文字列で指定する引数です。

• index 多変量モデルの場合に、このトレンド成分がどの観測値系列に対応するかを整数ベクトルで指定します。単変量モデルでは不要です。

• a1 状態の初期値 α_1（時刻t=1における状態の期待値）をベクトルで指定します。指定しない場合は、モデルから自動的に（通常は0として）設定されます。

• P1 初期状態の分散共分散行列 P_1 を指定します。これは初期状態 a1 の不確実性を表します。

• P1inf 初期状態分散のうち、分散が無限大（diffuse）であるかどうかを示す行列を0と1で指定します。1を指定した要素は分散が無限大であることを意味し、初期状態に関する情報が全くない場合（非定常な状態など）に用いられます。

• n 時系列の長さを指定します。通常はSSModelに渡されるデータから自動的に決定されるため、指定は不要です。

• state_names 状態の名前を文字列ベクトルで指定します。例えばdegree = 2の場合、c("level", "slope")のように指定できます。指定しない場合はデフォルトの名前が与えられます。

• ynames 観測値系列の名前を指定します。通常はデータから自動で取得されます。

6. KFSの引数

args(KFS)

function (model, filtering, smoothing, simplify = TRUE, transform = c("ldl", 
    "augment"), nsim = 0, theta, maxiter = 50, convtol = 1e-08, 
    return_model = TRUE, expected = FALSE, H_tol = 1e+15, transform_tol) 
NULL

KFS (Kalman Filter and Smoother) 関数は、SSModelで定義されたモデルオブジェクトに対して、カルマンフィルタリングと平滑化を実行するための関数です。

• model SSModelオブジェクト、またはfitSSMによって返されたモデルオブジェクトを指定します。分析対象となる状態空間モデルそのものです。

• filtering カルマンフィルタリングによって計算する要素を文字列ベクトルで指定します。指定できる主な要素は以下の通りです。

  • "state": フィルタ化状態 a_t（t時点までの観測値に基づくt時点の状態の推定値）とその分散 P_t。
  • "mean": フィルタ化された観測値の期待値 E(y_t | y_{1:t}) とその分散。
  • "signal": フィルタ化されたシグナル（Z_t * a_t）の推定値。
  • "disturbance": 状態攪乱項のフィルタ化推定値。 デフォルトでは何も計算しませんので、必要な要素を明示的に指定する必要があります。

• smoothing 平滑化によって計算する要素を文字列ベクトルで指定します。全ての観測値（y_{1:T}）を用いて過去の状態などを再推定します。指定できる主な要素は以下の通りです。

  • "state": 平滑化状態 alphahat_t（全観測値に基づくt時点の状態の推定値）とその分散 V_t。
  • "mean": 平滑化された観測値の期待値 E(y_t | y_{1:T}) とその分散。
  • "signal": 平滑化されたシグナル（Z_t * alphahat_t）の推定値。
  • "disturbance": 状態攪乱項の平滑化推定値。
  • "epsilon": 観測誤差の平滑化推定値。 シミュレーションの例では c("state", "mean") を指定し、平滑化状態と平滑化平均を計算しました。

• simplify TRUE（デフォルト）の場合、結果の配列が不要な次元を持つことを防ぎ、出力を単純化します。

• transform 内部的な数値計算アルゴリズムの種類を指定します。数値的安定性のためのもので、通常はデフォルトのままで問題ありません。

• nsim シミュレーション平滑化（simulation smoothing）を行う際のサンプルサイズを指定します。nsim > 0 とすると、平滑化された状態や応答の事後分布から、指定した数だけサンプルを生成します。モデルからのデータ生成や、より詳細な不確実性の評価に用いることができます。

• theta, maxiter, convtol これらは非ガウス状態空間モデルの近似計算（重要度サンプリングなど）を制御するための引数です。

  • theta: 重要度分布のパラメータ。
  • maxiter: 近似計算における最大反復回数。
  • convtol: 近似計算の収束判定に用いる許容誤差。

• return_model TRUE（デフォルト）の場合、戻り値のオブジェクトに入力として使われたmodelオブジェクトが含まれます。

• expected 非ガウスモデルの場合に、近似ガウスモデルのモード (FALSE) を使うか、期待値 (TRUE) を使うかを指定します。

• H_tol, transform_tol 数値計算上の許容誤差を指定するパラメータです。通常、ユーザーが変更する必要はありません。

7. fitSSMの引数

args(fitSSM)

function (model, inits, updatefn, checkfn, update_args = NULL, 
    ...) 
NULL

fitSSM関数は、SSModelで定義された状態空間モデルに含まれる未知のパラメータ（分散など、NAと指定されたパラメータ）を最尤法によって推定するための関数です。内部ではstats::optim関数を利用して、対数尤度が最大となるようなパラメータ値を探索します。

• model SSModelで作成されたモデルオブジェクトを指定します。このモデルオブジェクトの中に、推定対象となる未知パラメータ（通常はNAで指定）が含まれている必要があります。

• inits パラメータの初期値を指定する数値ベクトルです。最適化計算（尤度最大化）を開始する際の出発点となります。KFASでは、分散パラメータは対数変換されるなど、パラメータが制約なしの値（-∞から∞）を取るように変換されてから最適化が行われます。ですので、initsにはこの変換後のスケールでの初期値を指定する必要があります。例えば、分散Qを推定する場合、log(Q)の初期値を渡します。

• updatefn パラメータの推定プロセス中に、新しいパラメータ値でモデルを更新するための関数を指定します。KFASが提供する標準的なモデル成分（SSMtrendなど）を使用している限り、ユーザーがこの関数を自ら定義する必要は通常ありません。

• checkfn 最適化の各ステップで、提案されたパラメータ値が妥当かどうかをチェックするための関数を指定します。こちらも、通常は指定する必要はありません。

• update_args updatefnに渡す追加の引数をリスト形式で指定します。

• ... 内部で使用される最適化関数 stats::optim に渡す追加の引数を指定します。これにより、最適化のアルゴリズムや制御パラメータを調整できます。以下の引数が利用されます。

  • method: 最適化アルゴリズムを指定します。デフォルトは"BFGS"です。他の手法として"L-BFGS-B"（ボックス制約付き）や"Nelder-Mead"なども利用可能です。
  • control: 最適化の挙動を制御するためのパラメータをリストで指定します（例: control = list(maxit = 1000)で最大反復回数を設定）。

以上です。