状態空間時系列分析入門 第5章 (1/2)

{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)
  }
}

データの読み込み

ukdrivers <- read.table('../data/UKdriversKSI.txt', skip = 1)
ukdrivers <- ts(ukdrivers[[1]], start = c(1969, 1), frequency = 12)
ukdrivers <- log(ukdrivers)

ukpetrol <- read.table('../data/logUKpetrolprice.txt', skip = 1)
ukpetrol <- ts(ukpetrol, start = start(ukdrivers), frequency = frequency(ukdrivers))

モデル定義

確定的レベルと説明変数のあるローカルレベルモデル。

model_file <- '../models/fig05_01.stan'
cat(paste(readLines(model_file)), sep = '\n')

data {
  int<lower=1> n;
  vector[n] y;
  vector[n] x;
}
parameters {
  # 確定的レベル
  real mu;
  # 確定的回帰係数
  real beta;
  # 観測撹乱項
  real<lower=0> sigma_irreg;
}
transformed parameters {
  vector[n] yhat;
  for(t in 1:n) {
    yhat[t] <- mu + beta * x[t];
  }
}
model {
  # 式 5.3
  for(t in 1:n)
    y[t] ~ normal(yhat[t], sigma_irreg);
}

y <- ukdrivers
x <- ukpetrol

standata <- within(list(), {
  y <- as.vector(y)
  x <- as.vector(x)
  n <- length(y)
})

stan_fit <- stan(file = model_file, chains = 0)

## 
## TRANSLATING MODEL 'fig05_01' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'fig05_01' NOW.

fit <- pforeach(i = 1:4, .final = sflist2stanfit)({
  stan(fit = stan_fit, data = standata,
       iter = 2000, 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']
beta <- get_posterior_mean(fit, par = 'beta')[, 'mean-all chains']
sigma_irreg <- get_posterior_mean(fit, par = 'sigma_irreg')[, 'mean-all chains']

# stopifnot(is.almost.fitted(mu, 5.8787))
is.almost.fitted(mu, 5.8787)

## [1] "Result is  5.88187006315061"

## [1] FALSE

# stopifnot(is.almost.fitted(beta, -0.67166))
is.almost.fitted(beta, -0.67166)

## [1] "Result is  -0.670262729263585"

## [1] FALSE

stopifnot(is.almost.fitted(sigma_irreg^2, 0.0230137))

title <- 'Figure 5.1. Deterministic level and explanatory variable ‘log petrol price’.'
title <- '図 5.1 確定的レベルと説明変数「対数石油価格」'

# 原系列
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 <- paste('Figure 5.2. Conventional classical regression representation of ',
               'deterministic level and explanatory variable log petrol price.', sep = '\n')
title <- paste('図 5.2 確定的レベルと説明変数「対数石油価格」の',
               '伝統的な古典的回帰表現', sep = '\n')
df <- data.frame(drivers = as.vector(ukdrivers),
                petrol = as.vector(ukpetrol),
                stan = as.vector(ukpetrol) * beta + mu)

p <- ggplot(df, aes(x = petrol)) +
  geom_point(aes(y = drivers)) +
  geom_line(aes(y = stan), colour = 'blue') + 
  stat_smooth(aes(y = drivers), method = 'lm', colour = 'red',
              linetype = 'dashed', se = FALSE)
p + ggtitle(title)

title <- 'Figure 5.3. Irregular component for deterministic level model with explanatory variable ‘log petrol price’.'
title <- paste('図 5.3 説明変数「対数石油価格」のある',
               '確定的モデルの不規則要素', sep = '\n')
autoplot(y - yhat, ts.linetype = 'dashed') + ggtitle(title)