{rstan}「状態空間時系列分析入門」 を再現したい。やっていること / 数式はテキストを参照。

リポジトリ: https://github.com/sinhrks/stan-statespace

必要パッケージのインストール/読み込み

library(devtools)
# devtools::install_github('hoxo-m/pforeach')
# devtools::install_github('sinhrks/ggfortify')
library(rstan)
library(pforeach)
library(ggplot2)
ggplot2::theme_set(theme_bw(base_family="HiraKakuProN-W3"))
library(ggfortify)

関数定義

# モデルが収束しているか確認
is.converged <- function(stanfit) {
  summarized <- summary(stanfit)  
  all(summarized$summary[, 'Rhat'] < 1.1)
}

# 値がだいたい近いか確認
is.almost.fitted <- function(result, expected, tolerance = 0.001) {
  if (abs(result - expected) > tolerance) {
    print(paste('Result is ', result))
    return(FALSE)
  } else {
    return(TRUE)
  }
}

データの読み込み

fatalities <- read.table('../data/Norwayfinland.txt', skip = 1)
colnames(fatalities) <- c('year', 'Norwegian_fatalities',
                          'Finnish_fatalities')
norwegian_fatalities <- fatalities[['Norwegian_fatalities']]
norwegian_fatalities <- log(ts(norwegian_fatalities, start = 1970, frequency = 1))
finnish_fatalities <- fatalities[['Finnish_fatalities']]
finnish_fatalities <- log(ts(finnish_fatalities, start = 1970, frequency = 1))

モデル定義

確定的レベルと確率的傾きをもつローカル線形トレンドモデル。

model_file <- '../models/fig03_05.stan'
cat(paste(readLines(model_file)), sep = '\n')
data {
  int<lower=1> n;
  vector[n] y;
}
parameters {
  # 確定的レベル
  real mu;
  # 確率的傾き
  vector[n-1] v;
  # 傾き撹乱項
  real<lower=0> sigma_drift;
  # 観測撹乱項
  real<lower=0> sigma_irreg;
}
transformed parameters {
  vector[n] yhat;
  yhat[1] <- mu;
  for(t in 2:n) {
    yhat[t] <- yhat[t-1] + v[t-1];
  }
}
model {
  # 式 3.1
  for(t in 2:n-1)
    v[t] ~ normal(v[t-1], sigma_drift);
  for(t in 1:n)
    y[t] ~ normal(yhat[t], sigma_irreg);
}
y <- finnish_fatalities 

standata <-
  within(list(), {
    y <- as.vector(y)
    n <- length(y)
  })
stan_fit <- stan(file = model_file, chains = 0)
## 
## TRANSLATING MODEL 'fig03_05' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'fig03_05' NOW.
fit <- pforeach(i = 1:4, .final = sflist2stanfit)({
  stan(fit = stan_fit, data = standata, chains = 1, seed = i)
})
stopifnot(is.converged(fit))

yhat <- get_posterior_mean(fit, par = 'yhat')[, 'mean-all chains']
mu <- get_posterior_mean(fit, par = 'mu')[, 'mean-all chains']
v <- get_posterior_mean(fit, par = 'v')[, 'mean-all chains']
sigma_irreg <- get_posterior_mean(fit, par = 'sigma_irreg')[, 'mean-all chains']
sigma_drift <- get_posterior_mean(fit, par = 'sigma_drift')[, 'mean-all chains']
# stopifnot(is.almost.fitted(mu, 7.0133))
is.almost.fitted(mu, 7.0133)
## [1] "Result is  7.00996309884879"
## [1] FALSE
# stopifnot(is.almost.fitted(v[[1]], 0.0068482))
is.almost.fitted(v[[1]], 0.0068482)
## [1] "Result is  0.0116291036578503"
## [1] FALSE
# 誤植あるが値は正しい。
stopifnot(is.almost.fitted(sigma_irreg^2, 0.00320083))
stopifnot(is.almost.fitted(sigma_drift^2, 0.00153314))
title <- 'Figure 3.5.1. Trend of deterministic level and stochastic slope model for Finnish fatalities'
title <- paste('図 3.5.1 フィンランド事故データの',
               '確定的レベルと確率的傾きモデルのトレンド', sep = '\n')
# 原系列
p <- autoplot(y)

# stan
yhat <- ts(yhat, start = start(y), frequency = frequency(y))
p <- autoplot(yhat, p = p, ts.colour = 'blue')
p + ggtitle(title)

title <- 'Figure 3.5.2 Stochastic slope component for Finnish fatalities.'
title <- paste('図 3.5.2 フィンランド事故データの',
               '確定的レベルと確率的傾きモデルの確率的傾き要素', sep = '\n')
slope <- ts(v, start = start(y), frequency = frequency(y))
autoplot(slope) + ggtitle(title)

title <- 'Figure 3.6. Irregular component for Finnish fatalities.'
title <- '図 3.6 フィンランド事故データに対する不規則要素'

autoplot(y - yhat, ts.linetype = 'dashed') + ggtitle(title)