library(TMB)
library(ggplot2)
set.seed(123)
TMBで、状態空間モデルのあてはめをやってみました。ただ、モデル式をそのままコード化したものでは、最適化計算がうまく収束しませんでした。
準備
いつものように、TMBとggplot2を読み込みます。擬似乱数も固定します。
データ
ナイル川のデータを使用します。今回は、数値計算が安定するように流量の値を1000で割っておきました。元の値の単位が 108 m3 なので、変換後は 1011 m3 ということになります。また、tsクラスのデータをデータフレームに変換しておきました。
<- data.frame(Year = 1871:1970,
Nile_df Y = c(Nile) / 1000)
<- ggplot(Nile_df, aes(x = Year, y = Y)) +
p geom_line()
plot(p)
モデル
状態空間モデルにあてはめます。ここでは1回差分で線形正規のローカルレベルモデルとします。
システムモデル
\[ \alpha_{t} = \alpha_{t-1} + \eta_{t}, \quad \eta_{t} \sim N(0, \sigma_{\eta}^2) \]
観測モデル
\[ Y_{t} = \alpha_{t} + \epsilon_{t}, \quad \epsilon_{t} \sim N(0, \sigma_{\epsilon}^2) \]
ローカルレベルモデルをTMBで記述します。まずは上のモデル式をそのままコードにしてみました。
ssm1.cpp
// State space model
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_VECTOR(Y);
DATA_SCALAR(m0);
DATA_SCALAR(C0);
PARAMETER_VECTOR(alpha);
PARAMETER(log_sigma_eta);
PARAMETER(log_sigma_eps);
int n = Y.size();
Type ll = dnorm(alpha(0), m0, sqrt(C0), true);
for (int t = 1; t < n; t++) {
ll += dnorm(alpha(t), alpha(t - 1), exp(log_sigma_eta), true);
}
for (int t = 0; t < n; t++) {
ll += dnorm(Y(t), alpha(t), exp(log_sigma_eps), true);
}
return -ll;
}
コンパイルと最適化
モデルをコンパイルして、できたライブラリをロードします。
<- "ssm1"
model_name file.path("models", paste(model_name, "cpp", sep = ".")) |>
compile()
## Note: Using Makevars in /Users/hiroki/.R/Makevars
## using C++ compiler: 'Apple clang version 15.0.0 (clang-1500.3.9.4)'
## using SDK: ''
## [1] 0
file.path("models", dynlib(model_name)) |>
dyn.load()
データと、パラメータの初期値を用意して最適化を実行します。
しかし、このモデルでは収束しませんでした。nlminb
関数のeval.max
とiter.max
を増やしてもダメでした。
<- list(Y = Nile_df$Y, m0 = 0, C0 = 1000)
data <- list(alpha = rep(1, nrow(Nile_df)),
parameters log_sigma_eta = 0,
log_sigma_eps = 0)
<- MakeADFun(data, parameters, DLL = model_name)
obj <- nlminb(obj$par, obj$fn, obj$gr,
opt control = list(eval.max = 100,
iter.max = 100))
結果です。convergence
が1なので、収束していません。message
にも出ていますが。
print(opt)
## $par
## alpha alpha alpha alpha alpha
## 1.0677169 1.0673157 1.0666821 1.0662118 1.0651758
## alpha alpha alpha alpha alpha
## 1.0635592 1.0619342 1.0602309 1.0585951 1.0557058
## alpha alpha alpha alpha alpha
## 1.0523299 1.0491020 1.0462711 1.0427704 1.0395511
## alpha alpha alpha alpha alpha
## 1.0361627 1.0327649 1.0289340 1.0256462 1.0222775
## alpha alpha alpha alpha alpha
## 1.0184435 1.0143659 1.0098995 1.0044455 0.9988440
## alpha alpha alpha alpha alpha
## 0.9923429 0.9850908 0.9779326 0.9701240 0.9629786
## alpha alpha alpha alpha alpha
## 0.9563019 0.9495046 0.9433231 0.9374785 0.9318846
## alpha alpha alpha alpha alpha
## 0.9267077 0.9218265 0.9175848 0.9134062 0.9087132
## alpha alpha alpha alpha alpha
## 0.9041769 0.8996056 0.8957796 0.8932439 0.8908095
## alpha alpha alpha alpha alpha
## 0.8887567 0.8867447 0.8836058 0.8814222 0.8791589
## alpha alpha alpha alpha alpha
## 0.8773187 0.8757948 0.8745801 0.8732096 0.8720267
## alpha alpha alpha alpha alpha
## 0.8714135 0.8707897 0.8706811 0.8708152 0.8705659
## alpha alpha alpha alpha alpha
## 0.8702835 0.8707520 0.8709854 0.8712781 0.8715063
## alpha alpha alpha alpha alpha
## 0.8713424 0.8711903 0.8713660 0.8707772 0.8708092
## alpha alpha alpha alpha alpha
## 0.8710026 0.8720633 0.8727300 0.8738704 0.8752079
## alpha alpha alpha alpha alpha
## 0.8766028 0.8777318 0.8786572 0.8800693 0.8808317
## alpha alpha alpha alpha alpha
## 0.8818143 0.8831746 0.8846187 0.8864180 0.8871512
## alpha alpha alpha alpha alpha
## 0.8881416 0.8889625 0.8896293 0.8900972 0.8904886
## alpha alpha alpha alpha alpha
## 0.8907659 0.8907845 0.8907253 0.8899792 0.8894567
## alpha alpha alpha alpha alpha
## 0.8879432 0.8872327 0.8861510 0.8854612 0.8852800
## log_sigma_eta log_sigma_eps
## -5.1123217 -2.1532164
##
## $objective
## [1] -450.4888
##
## $convergence
## [1] 1
##
## $iterations
## [1] 88
##
## $evaluations
## function gradient
## 100 88
##
## $message
## [1] "function evaluation limit reached without convergence (9)"
カルマンフィルタの尤度を使ったモデル
単純なコード化ではダメだったので、カルマンフィルタの尤度を使用することにしました。尤度の式は野村(2016)を参考にしました。最適化のためなら定数部分は省略してもよさそうですが、いちおうちゃんと尤度を計算するようにしています。
ssm2.cpp
// State space model using Kalman filter
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_VECTOR(Y);
DATA_SCALAR(a0);
DATA_SCALAR(P0);
PARAMETER(log_sigma_eta);
PARAMETER(log_sigma_eps);
int n = Y.size();
vector<Type> P(n);
vector<Type> F(n);
vector<Type> K(n);
vector<Type> L(n);
vector<Type> a(n);
vector<Type> v(n);
Type ll = -n / 2.0 * log(2.0 * M_PI);
P(0) = P0;
a(0) = a0;
for (int t = 0; t < n; t++) {
v(t) = Y(t) - a(t);
F(t) = P(t) + pow(exp(log_sigma_eps), 2.0);
K(t) = P(t) / F(t);
L(t) = 1.0 - K(t);
if (t < n - 1) {
P(t + 1) = P(t) * L(t) + pow(exp(log_sigma_eta), 2.0);
a(t + 1) = a(t) + K(t) * v(t);
}
ll += - 0.5 * (log(F(t)) + pow(v[t], 2.0) / F(t));
}
return -ll;
}
コンパイルと最適化
モデルをコンパイルして、できたライブラリをロードします。
<- "ssm2"
model_name file.path("models", paste(model_name, "cpp", sep = ".")) |>
compile()
file.path("models", dynlib(model_name)) |>
dyn.load()
データと、パラメータ(システムモデルと観測モデルの誤差の標準偏差の対数)の初期値を用意して最適化を実行します。
<- list(Y = Nile_df$Y, a0 = 0, P0 = 1000)
data <- list(log_sigma_eta = 1,
parameters log_sigma_eps = 1)
<- MakeADFun(data, parameters, DLL = model_name)
obj <- nlminb(obj$par, obj$fn, obj$gr)
opt ## outer mgc: 54.811
## outer mgc: 47.67067
## outer mgc: 30.13291
## outer mgc: 9.976196
## outer mgc: 4.905411
## outer mgc: 5.577532
## outer mgc: 2.403372
## outer mgc: 0.5632631
## outer mgc: 0.5531079
## outer mgc: 0.1161105
## outer mgc: 0.0002061353
## outer mgc: 3.300339e-06
結果です。こちらはうまく収束しました。
print(opt)
## $par
## log_sigma_eta log_sigma_eps
## -3.261528 -2.096579
##
## $objective
## [1] -46.94871
##
## $convergence
## [1] 0
##
## $iterations
## [1] 11
##
## $evaluations
## function gradient
## 18 12
##
## $message
## [1] "relative convergence (4)"
DLM
比較のため、DLMパッケージを使用して、システムモデルと観測モデルの誤差分散の最尤推定をやってみました。
library(dlm)
##
## 次のパッケージを付け加えます: 'dlm'
## 以下のオブジェクトは 'package:ggplot2' からマスクされています:
##
## %+%
<- function(x) {
buildFun dlmModPoly(order = 1, dV = exp(x[1]), dW = exp(x[2]),
m0 = 0, C0 = 1000)
}<- dlmMLE(Nile_df$Y, parm = c(0, 0), buildFun)
fit <- buildFun(fit$par) dlmNile
結果です。V(dlmNile)
が観測モデルの誤差分散、W(dlmNile)
がシステムモデルの誤差分散になります。
V(dlmNile)
## [,1]
## [1,] 0.01509853
W(dlmNile)
## [,1]
## [1,] 0.001469168
先ほどのTMBの結果を比較してみます。dlmの結果とほぼ同じでした。
exp(opt$par["log_sigma_eps"])^2
## log_sigma_eps
## 0.01509853
exp(opt$par["log_sigma_eta"])^2
## log_sigma_eta
## 0.001469171
平滑化曲線
TMBで推定された誤差標準偏差の値を使ってdlmで平滑化した値を求め、グラフにしました。
<- buildFun(c(2 * opt$par["log_sigma_eps"],
dlm2 2 * opt$par["log_sigma_eta"]))
<- dlmSmooth(Nile_df$Y, dlm2)
smoothed data.frame(Year = 1871:1970,
Y = Nile_df$Y,
theta = dropFirst(smoothed$s)) |>
ggplot() +
geom_line(aes(x = Year, y = Y)) +
geom_line(aes(x = Year, y = theta), colour = "red")
参考文献
- 野村俊一 (2016) カルマンフィルタ—Rを使った時系列予測と状態空間モデル—. 共立出版, 154pp.