Stan を使うことで、状態空間モデルのパラメーターを簡単に推定することができる。このページではナイル川の流量データに対して、状態空間モデルを適用し、パラメーターを推定する方法を示す。また、推定されたモデルを使用して、さらに 10 年度の流量データを予測する例も合わせて示す。

data(Nile)
plot(Nile)

状態空間モデルの方程式

状態空間モデルは、状態とよばれる観測されない潜在的な確率変数 x を導入して、その状態を使用して観測値 y を説明しようするとモデルである。

状態空間モデルでは、状態と観測値の 2 つの変数で構成されれ、時系列データが定常性を持つとき、t 時点で観測される観測値 y_t は状態 x_t から生成されるとされる。そして、状態 x_t は、状態 x_t-1 のみに依存して定まる。これらのことを方程式で表すと次のようになる。これらのことを方程式で表すと次のようになる。

状態が定常性を持つのであれば、 t 時点の状態と t - 1 時点の状態の差は、平均ゼロ、標準偏差 σ_ε の正規分布に従うと考えられる。また、状態 x から観測値・ナイル川流量 y が観測されるときに、観測ノイズ δ が生じる。ナイル川の流量は連続値であるから、観測ノイズを平均ゼロ、標準偏差 σ_δ の正規分布に従うと仮定することができる。このとき、ε および δ は、次のようにモデル化できる。

以上、まとめると、ナイル川流量に関する状態空間モデルのモデル式は次のように書ける。

Stan コード

上で求めた状態空間モデルの方程式を Stan コードで記述すると次のようになる。ただし、パラメーター推定後に、Npred 年後の流量を予測したいので、ここで generated quantities ブロックの中で、推定されたパラメーターを用いて、N + 1 年から N + Npred 年までの流量データのサンプリングを行なっている。

data {
    int N;
    int Npred;
    vector[N] Y;
}
parameters {
    real mu0;
    vector[N] mu;
    real<lower=0> sigma_obs;
    real<lower=0> sigma_sys;
}
model {
    // state equation
    mu[1] ~ normal(mu0, sigma_sys);
    for(n in 2:N) {
        mu[n] ~ normal(mu[n - 1], sigma_sys);
    }

    // observation  equation
    for(n in 1:N) {
        Y[n] ~ normal(mu[n], sigma_obs);
    }
}
generated quantities {
    vector[N + Npred] mu_hat;
    vector[Npred] y_hat;

    mu_hat[1:N] = mu;
    for (n in 1:Npred) {
        mu_hat[N + n] = normal_rng(mu_hat[N + n - 1], sigma_sys);
        y_hat[n] =  normal_rng(mu_hat[N + n], sigma_obs);
    }
}

状態空間モデルのパラメーター推定

R で上述の Stan コードを読み込み、ナイル川流量データ（Nile）を整形して代入し、パラメーター推定を行う。

library(rstan)
d <- list(Y = as.numeric(Nile), N = length(Nile), Npred = 10)
fit <- stan(file = "ss.stan", data = d)

パラメーターの推定結果として、状態の平均および分散、観測値の平均および分散、そして、N + 10 年後のサンプリング結果など多くの情報が含まれている。ここでは、状態の平均、状態のサンプリング結果、観測値のサンプリング結果だけを extract 関数で取り出す。

ms <- rstan::extract(fit, pars = c('mu', 'mu_hat', 'y_hat'))

取り出したデータを図示する。図は次のような手順で作成する。（ただ、ggplot で最後に実行した関数が、画像のもっとも上の層として描かれるので、よりきれいな図を作成するには、関数の順序を入れ変えて作図したほうがいい。例えば、geom_lin 関数を一番最後に実行するなど。）

1. 1 年目から N 年目までの状態の推定値を図示。推定された状態の平均値を赤線で、推定された 95% 信頼区間を薄い赤色の塗る。
2. N + 1年目から N + Npred 年目までの状態の予測値を図示。予測された状態の平均値を赤線で、予測された 95% 信頼区間を薄い赤い色で塗る。
3. N + 1年目から N + Npred 年目までの観測値の予測値を図示。予測された観測値の 95% 予測区間を灰色で塗る。（予測された観測値の平均値は、状態の平均値とほぼ同一であるため、描かないようにした。）

library(ggplot2)

df.past <- data.frame(x = time(Nile),
                      y = as.numeric(Nile),
                      lower = apply(ms$mu, 2, quantile, prob = 0.025),
                      median = apply(ms$mu, 2, quantile, prob = 0.500),
                      upper = apply(ms$mu, 2, quantile, prob = 0.925))
df.pred.mu <- data.frame(x = seq(time(Nile)[d$N] + 1, time(Nile)[d$N] + d$Npred),
                         lower = apply(ms$mu_hat[, (d$N + 1):(d$N + d$Npred)], 2, quantile, prob = 0.025),
                         median = apply(ms$mu_hat[, (d$N + 1):(d$N + d$Npred)], 2, quantile, prob = 0.500),
                         upper = apply(ms$mu_hat[, (d$N + 1):(d$N + d$Npred)], 2, quantile, prob = 0.925))
df.pred.y <- data.frame(x = seq(time(Nile)[d$N] + 1, time(Nile)[d$N] + d$Npred),
                        lower = apply(ms$y_hat, 2, quantile, prob = 0.025),
                        median = apply(ms$y_hat, 2, quantile, prob = 0.500),
                        upper = apply(ms$y_hat, 2, quantile, prob = 0.925))

g <- ggplot()

# observed data
g <- g + geom_line(data = data.frame(x = time(Nile), y = as.numeric(Nile)), aes(x = x, y = y))
g <- g + geom_ribbon(aes(x = x, ymin = lower, ymax = upper), fill = '#CC0C00', alpha = 0.3, data = df.past)
g <- g + geom_line(aes(x = x, y = median), color = '#CC0C00', data = df.past)

# predicted results
g <- g + geom_line(aes(x = x, y = median), color = '#CC0C00', data = df.pred.mu)
g <- g + geom_ribbon(aes(x = x, ymin = lower, ymax = upper), fill = '#CC0C00', alpha = 0.3, data = df.pred.mu)
g <- g + geom_ribbon(aes(x = x, ymin = lower, ymax = upper), fill = '#000000', alpha = 0.3, data = df.pred.y)

g <- g + xlab('Year') + ylab('Nile')
print(g)