Rのパッケージ:dlm

order

• モデルの次数を指定する整数です。これは状態ベクトル \(\boldsymbol{\theta}_t\) の次元数を決定し、モデルの複雑さを定義します。
• order = 1: ランダムウォークモデル（またはローカルレベルモデル）を意味します。
  • 状態はレベル（水準）\(\mu_t\) のみです。
  • 状態方程式: \(\mu_t = \mu_{t-1} + w_t\)
  • 時系列の基本的な水準が時間とともにランダムに変動する様子をモデル化します。
• order = 2: ローカル線形トレンドモデルを意味します。
  • 状態はレベル \(\mu_t\) とトレンド（傾き）\(\delta_t\) の2つです。
  • 状態方程式:
    • \(\mu_t = \mu_{t-1} + \delta_{t-1} + w_{\mu, t}\) （レベルは前回のレベル＋前回のトレンド分だけ変化）
    • \(\delta_t = \delta_{t-1} + w_{\delta, t}\) （トレンド自体もランダムに変化）
  • 時系列のレベルと、その変化率であるトレンドの両方が時間とともに変動する、より柔軟なモデルです。
• order = 3: レベル、トレンドに加え、トレンドの変化率（加速度）までを状態に含むモデルになります。
• デフォルト値: 2

dV

• 観測ノイズの分散 \(V_t\) を指定する数値です。
  • 観測方程式 \(y_t = \boldsymbol{F}_t \boldsymbol{\theta}_t + v_t\) における、ノイズ項 \(v_t\) の分散を意味します。
  • これは、モデルで説明しきれないランダムな変動の大きさを表します。この値が大きいほど、観測値はモデルが示す基本的な動き（レベルやトレンド）から大きくばらつきます。
  • この引数で設定された値は、作成されるdlmオブジェクトのVコンポーネントに格納されます。
  • 実際には観測ノイズの分散は未知であることがほとんどですので、dlmMLE関数その他によりデータから最尤推定する必要があります。その場合、このdVの初期値は推定の出発点として使われますが、最終的には推定された値で上書きされます。
• デフォルト値: 1

dW

• システムノイズ（状態遷移ノイズ）の分散を指定するベクトルです。
  • 状態方程式 \(\boldsymbol{\theta}_t = \boldsymbol{G}_t \boldsymbol{\theta}_{t-1} + \boldsymbol{w}_t\) における、ノイズ項 \(\boldsymbol{w}_t\) の共分散行列 \(W_t\) の対角成分を設定します。
  • dWはorderと同じ長さのベクトルで指定します。dW[i]は、状態ベクトルのi番目の要素が時間変化する際のランダムな変動の大きさを表します。
  • 例えば order = 2 の場合、dW = c(dW_level, dW_trend) となり、dW_levelがレベルのシステムノイズ分散、dW_trendがトレンドのシステムノイズ分散に対応します。
• デフォルト値: c(rep(0, order - 1), 1)
  • order = 1 の場合: dW = 1。レベルのシステムノイズ分散が1であることを意味します。
  • order = 2 の場合: dW = c(0, 1)。これは、レベルのシステムノイズ分散が0、トレンドのシステムノイズ分散が1であることを意味します。レベルの変動は完全にトレンド \(\delta_{t-1}\) によって決まり、レベル自体にはランダムな変動が加わらない（\(w_{\mu, t}=0\)）という、やや強い仮定を置いています。

m0

• 初期状態ベクトルの事前分布の平均 \(\boldsymbol{m}_0 = E[\boldsymbol{\theta}_0]\) を指定するベクトルです。
  • これは、分析を開始する時点（\(t=0\)）での各状態（レベルやトレンド）の期待値を意味します。
  • orderと同じ長さのベクトルで指定します。
  • 例えば、時系列データの最初の値が100あたりであれば、レベルの初期値として m0 = c(100, 0) のように設定することも考えられます（レベルの初期値100、トレンドの初期値0）。
  • この引数で設定された値は、作成されるdlmオブジェクトのm0コンポーネントに格納されます。
• デフォルト値: rep(0, order)
  • すべての状態の初期値の期待値が0であることを意味します。後述のC0を大きく設定することで、この初期値の影響はデータによって更新されます。

C0

• 初期状態ベクトルの事前分布の共分散行列 \(\boldsymbol{C}_0 = \text{Var}(\boldsymbol{\theta}_0)\) を指定する行列です。
  • これは、分析を開始する時点（\(t=0\)）での各状態の不確実性の大きさを表します。
  • order行 order列の正方行列で指定します。
  • この引数で設定された値は、作成されるdlmオブジェクトのC0コンポーネントに格納されます。
• デフォルト値: 1e+07 * diag(nrow = order)
  • これは、対角成分が非常に大きな値（\(10^7\)）で、非対角成分が0の対角行列です。
  • 対角成分（分散）が非常に大きいということは、初期状態に関する事前知識がほとんどなく、不確実性が非常に高いことを意味します（散漫事前分布）。
  • このように設定することで、カルマンフィルタは初期状態の事前情報にほとんど依存せず、観測データに基づいて状態を推定するようになります。

frequency

• 季節性の周期を指定する整数です。この関数において唯一、指定が必須の引数です。
• 具体例:
  • 月次データの場合: frequency = 12
  • 四半期データの場合: frequency = 4
  • 週次データ（曜日効果）の場合: frequency = 7
• 状態ベクトルの次元数（frequency - 1）や、状態遷移行列\(\boldsymbol{G}\)の形が本引数の値に基づいて自動的に設定されます。

dV

• 観測ノイズの分散 \(V_t\) を指定する数値です。
  • これは、モデルが観測値\(y_t\)を説明する際に生じる誤差（\(y_t = \text{季節効果}_t + v_t\) における \(v_t\)）の分散を意味します。
• dlmModSeasは通常、dlmModPoly（トレンドモデル）など他のモデルと足し合わせて使用されます。dlmでは、複数のモデルを + で結合した場合、観測ノイズの分散Vは、原則として最初に指定されたモデルの値が使われるか、結合後のモデル全体に対して明示的に設定されます。
• よって、dlmModSeasでdVを指定しても、最終的なモデルではその値が使われないことがあります。この引数は主に、季節モデルを単独で使用する場合や、モデルの仕様を形式的に整えるために設けられていると考えられます。
• デフォルト値: 1

dW

• システムノイズの分散を指定するベクトルです。これは、季節パターン自体が時間と共にどれだけ変化するかを制御する引数です。
  • 状態方程式 \(\boldsymbol{\theta}_t = \boldsymbol{G}_t \boldsymbol{\theta}_{t-1} + \boldsymbol{w}_t\) における、ノイズ項\(\boldsymbol{w}_t\)の共分散行列\(W_t\)の対角成分を設定します。ベクトルの長さは frequency - 1 である必要があります。
• 設定による挙動の違い:
  • 動的季節モデル (Dynamic Seasonal Model): dWにゼロでない値を設定すると、季節パターンが時間と共に変化することを許容します。
    • デフォルト値 c(1, rep(0, frequency - 2)): この設定では、システムノイズの最初の要素にのみ分散を与え、残りはゼロにしています。これにより、季節パターン全体が滑らかに変化していくようなモデルになります。例えば、ある年の「12月の売上のピーク」が、翌年には少し大きくなったり小さくなったりするような変動を捉えることができます。
  • 静的季節モデル (Static / Fixed Seasonal Model): dWをすべてゼロに設定した場合（例: dW = rep(0, frequency - 1)）、季節パターンは時間を通じて変化しないと仮定されます。
    • この設定は、データ全体で一つの共通した季節パターンが存在すると仮定するものです。最初に推定された季節パターンが、将来にわたって全く同じ形で繰り返されることになります。
• パラメータ推定: dWの値を事前に決めるのが難しい場合、dlmMLEを用いてデータから推定することが可能です。

m0

• 初期状態ベクトルの事前分布の平均 \(\boldsymbol{m}_0\) を指定するベクトルです。
  • 分析を開始する時点（\(t=0\)）での、frequency - 1個の季節効果の期待値を意味します。
  • この引数で設定された値は、作成されるdlmオブジェクトのm0コンポーネントに格納されます。
• デフォルト値: rep(0, frequency - 1)
  • 初期の季節効果の期待値はすべて0であると仮定します。後述のC0を大きく設定することで、この初期値の影響はデータによってすぐに更新されます。

C0

• 初期状態ベクトルの事前分布の共分散行列 \(\boldsymbol{C}_0\) を指定する行列です。
  • 分析を開始する時点（\(t=0\)）での季節効果の不確実性の大きさを表します。
  • frequency - 1行 frequency - 1列の正方行列で指定します。
• デフォルト値: 1e+07 * diag(nrow = frequency - 1)
  • 対角成分が非常に大きな値（\(10^7\)）である対角行列です。これは、初期の季節パターンについては何も情報がない（不確実性が非常に高い）ということを意味します（散漫事前分布）。
  • このように設定することで、カルマンフィルタは初期状態の事前情報にほとんど依存せず、観測データに基づいて季節パターンを学習します。

y

• 分析対象の時系列データを指定します。dlmMLEはこのデータ y の尤度を最大化するパラメータを探します。
• データ形式:
  • tsオブジェクト: Rで時系列データを扱う標準的な形式です。開始時点や周期といった情報が含まれているため、最も扱いやすい形式です。
  • 数値ベクトル (vector): 単変量の時系列データの場合に利用できます。
  • 行列 (matrix): 各列が異なる変数を表す、多変量の時系列データの場合に利用します。行が時間の進行に対応します。
• この引数は必ず指定する必要があります。

parm

• 推定したいパラメータの初期値を数値ベクトルで指定します。
• dlmMLEは内部で数値最適化アルゴリズム（デフォルトではoptim関数）を用いて、最適なパラメータを探します。parmは、その探索を開始する出発点となります。適切な初期値を与えることで、アルゴリズムがより速く、またより良い解（尤度が最も高い点）に収束する可能性が高まります。
• DLMで推定するパラメータが分散である場合、分散は常に正の値を取る必要がありますが、最適化アルゴリズムはパラメータ空間を制約なしに探索しようとすることがあります。
• この問題を回避するための一つのテクニックとして、対数スケールで最適化を行います。つまり、build関数の中で exp(par) のようにパラメータを指数変換して分散として用いるのです。こうすることで、par自体が負の値を取っても、exp(par)は常に正となり、分散の制約を満たします。
• このテクニックを用いる場合、parmには分散の対数を取った値を初期値として設定します。例えば、分散の初期値を1としたい場合は log(1) = 0 ですので、parm = 0 と設定します。

build

• dlmモデルを構築するための関数を指定します。このbuild関数は、dlmMLEが最適化の各ステップで呼び出す、ユーザー定義の「設計図」のようなものです。
• 関数の仕様:
  1. 引数: parmで指定したのと同じ長さのパラメータベクトルを、唯一の引数として受け取る必要があります。
  2. 返り値: dlmオブジェクトを返す必要があります。
• 内部的な動作: dlmMLEは以下の処理を繰り返します。
  1. 最適化アルゴリズムが、新しいパラメータの候補値 p を生成します。
  2. dlmMLEは、ユーザーが定義したbuild関数を build(p) のように呼び出します。
  3. build関数は、パラメータ p を使ってdlmオブジェクト（例えば、VやdWにexp(p)を設定したモデル）を構築し、返します。
  4. dlmMLEは、その返されたdlmオブジェクトとデータyを使って尤度を計算します。
  5. 最適化アルゴリズムは、計算された尤度を基に、さらに尤度が高くなりそうな次のパラメータ候補値を探索します。
• この引数により、ユーザーはbuild関数を自由に記述することで、カスタムモデルのパラメータ推定も可能になります。

method

• 尤度を最大化するために内部で使用する、数値最適化アルゴリズムを指定します。
• dlmMLEは、Rの標準的な汎用最適化関数である stats::optim を利用しています。よって、method引数にはoptimで利用可能なアルゴリズム名を文字列で指定します。
• デフォルト値: "L-BFGS-B"
  • これは、ボックス制約付き準ニュートン法です。パラメータの取りうる値に上限や下限（ボックス制約）を設定できる特徴があります。
• その他の選択肢: "Nelder-Mead", "BFGS", "CG", "SANN" など、optimがサポートする他の手法も指定できます。

...

• 説明: optim関数に渡す、その他の追加引数を指定するために使用します。
• 最適化プロセスの挙動を細かく制御するためのcontrolリストを渡すことが出来ます。
• 具体例:
  • control = list(maxit = 1000): 最大反復回数をデフォルトの100回から1000回に増やす。
  • control = list(trace = 1): 最適化の過程（反復ごとのパラメータ値や目的関数値）を表示する。
  • control = list(reltol = 1e-12): 収束判定の相対的な許容誤差をより厳しくする。

  # 使用例
  dlmMLE(y, parm, build, control = list(trace = 1, maxit = 1000))

debug

• デバッグ情報を出力するかどうかを指定する論理値（TRUEまたはFALSE）です。
• debug = TRUEの場合:
  • dlmMLEがbuild関数を呼び出すたびに、その際に渡されたパラメータベクトルと、計算された負の対数尤度（optimは最小化を行うため）の値を画面に出力します。
  • これは、control = list(trace = ...)よりも詳細な情報を提供します。
• 利用シーン:
  • build関数が特定のパラメータ値でエラーを返す場合、どの値で問題が起きているのかを特定するのに役立ちます。
  • 最適化が予期せぬ挙動を示す場合に、内部で何が起きているかを調査する手がかりとなります。
• デフォルト値: FALSE（デバッグ情報は表示されません）。

# --- モデル構築関数の定義 ---
# パラメータ(ノイズの分散)を引数として、dlmオブジェクトを返す関数を定義します
# これにより、最尤法で最適なパラメータを探すことができます
build_sales_model <- function(par) {
  # ローカル線形トレンドモデル (レベル + トレンド)
  trend_model <- dlmModPoly(order = 2, dV = 0, dW = exp(par[1:2]))

  # 12か月周期の季節モデル
  # 季節性の変動はないと仮定し、dWは0に固定します
  seasonal_model <- dlmModSeas(frequency = 12, dV = 0, dW = rep(0, 11))

  # モデルを結合します
  model <- trend_model + seasonal_model

  # 観測ノイズの分散を設定します
  model$V <- exp(par[3])

  return(model)
}

# --- パラメータの最尤推定 ---
# dlmMLEを用いて、モデルのパラメータを推定します
# parmにはパラメータの初期値を設定します
# buildには上で定義したモデル構築関数を指定します
mle_fit <- dlmMLE(
  y = sales_ts,
  parm = c(log(1), log(1), log(10)),
  build = build_sales_model
)

# 収束の確認 (0なら成功)
cat("最尤推定の収束コード:", mle_fit$convergence, "(0は成功を示します)\n")

# 推定されたパラメータを使って、最終的なモデルを構築します
sales_model_est <- build_sales_model(mle_fit$par)

# --- 推定された分散パラメータの確認 ---
# dlmでは分散が負にならないように対数スケールで最適化を行うため、exp()で元に戻します
estimated_variances <- exp(mle_fit$par)
cat("\n推定された分散パラメータ:\n")
cat("レベルのシステムノイズ分散 (W[1,1]):", estimated_variances[1], "\n")
cat("トレンドのシステムノイズ分散 (W[2,2]):", estimated_variances[2], "\n")
cat("観測ノイズ分散 (V):", estimated_variances[3], "\n")

最尤推定の収束コード: 0 (0は成功を示します)

推定された分散パラメータ:
レベルのシステムノイズ分散 (W[1,1]): 284.6552 
トレンドのシステムノイズ分散 (W[2,2]): 3.30265 
観測ノイズ分散 (V): 17.69443 

最尤推定によって、データに最もフィットするノイズの分散が計算されました。これらの値が設定されたモデルを用いて、次のステップに進みます。