#はじめに

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装
2018/3/17、amazonで購入。
第3部のVARモデルとARCH・GARCHモデルはスキップ

#単位根検定の種類

ADF検定(原系列に単位根がある) {urca::ur.df}
KPSS検定(原系列に単位根がない) {urca::ur.kpss}
PO検定(共和分がない) {urca::ca.po}
DW検定(残差に自己相関がない) {lmtest::dwtest}

#時系列データに回帰分析を適用する際のフローチャート

参考　時系列データへの回帰分析

単位根検定
├単位根あり
│└共和分検定：PO検定
│　├共和分あり→普通の回帰(OLS)
│　└共和分無し→差分系列へ回帰(OLS)
└単位根無し
　└残差の自己相関の検定：DW検定
　　├自己相関あり→一般化最小二乗法(GLS)
　　└自己相関無し→普通の回帰(OLS)

#ARMAの特性多項式の解

(1)AR過程の特性多項式の全ての解の絶対値が1より大
・MA(∞)過程で表現可能
・AR過程は定常となる
・MA過程は定義から元々定常

(2)MA過程の特性多項式の全ての解の絶対値が1より大
・AR(∞)過程で表現可能
・MA過程は反転可能となる (現在および過去の$x_t$から現在の$\epsilon_t$が計算可能であること)
・AR過程は定義から元々反転可能

任意のMA過程に対して、それと同一の平均、同一の自己相関関数を持つMA過程が複数存在する。
但し反転可能なものは1つしかない。

#散漫カルマンフィルター

Diffuse Kalman Filter
状態の分散の初期値を無限大とし、対数尤度の計算からは除外する。
対数尤度を最大化する過程誤差と観測誤差の分散を求める。

#KFAS

ローカルレベルモデル

build_kfas = SSModel(Nile~SSMtrend(degree=1,Q=NA),H=NA)
fit_kfas = fitSSM(build_kfas,inits=c(1,1))

ローカル線形トレンドモデル

build_kfas = SSModel(Nile~SSMtrend(degree=2,Q=list(NA,NA)),H=NA)
fit_kfas = fitSSM(build_kfas,inits=c(1,1,1))

ローカル線形トレンド季節性モデル

build_kfas = SSModel(AirPassengers~SSMtrend(degree=2,Q=list(NA,NA))+SSMseasonal(period=12,Q=NA),H=NA)
fit_kfas = fitSSM(build_kfas,inits=c(1,1,1,1))

#Stan

ARモデル

AR.stan

data {
  int<lower=1> T;
  real y[T];
}

parameters {
  real mu;
  real alpha;
//  real beta;
  real<lower=0> sigma;
}

model {
  real eps;
  mu ~ normal(0, 10);
  alpha ~ normal(0, 2);
//  beta ~ normal(0, 2);
  sigma ~ cauchy(0, 5);
  eps = 0;
  for(t in 2:T){
    eps = y[t] - (mu + alpha * y[t-1]);
    eps ~ normal(0, sigma);
  }
}

ARMAモデル

ARMA.stan

data {
  int<lower=1> T;
  real y[T];
}

parameters {
  real mu;
  real alpha;
  real beta;
  real<lower=0> sigma;
}

model {
  real eps;
  mu ~ normal(0, 10);
  alpha ~ normal(0, 2);
  beta ~ normal(0, 2);
  sigma ~ cauchy(0, 5);
  eps = 0;
  for(t in 2:T){
    eps = y[t] - (mu + alpha * y[t-1] + beta * eps);
    eps ~ normal(0, sigma);
  }
}

ちなみにサンプルのARMAモデルは少し違っている。

ARMA.stan

data {
  int<lower=1> T;
  real y[T];
}

parameters {
  real mu;
  real alpha;
  real beta;
  real<lower=0> sigma;
}

model {
  real eps;
  mu ~ normal(0, 10);
  alpha ~ normal(0, 2);
  beta ~ normal(0, 2);
  sigma ~ cauchy(0, 5);
  eps = y[1] - (mu + alpha * mu);
  eps ~ normal(0, sigma);
  for(t in 2:T){
    eps = y[t] - (mu + alpha * y[t-1] + beta * eps);
    eps ~ normal(0, sigma);
  }
}

ARMA(P,Q)モデル
最終的にはこんな感じにしてみました。

ARMA(P,Q).stan

// ARMA(P,Q): y[t] = mu + alpha[1]*y[t-1] + beta[1]*e[t-1] + e[t]

data {
  int<lower=1> T;
  real y[T];
  int P;
  int Q;
}

parameters {
  real mu;
  real alpha[P];
  real beta[Q];
  real<lower=0> sigma;
  real state[P];
}

transformed parameters {
  real prd[T];
  for(t in 1:P){
    prd[t] = state[t];
  }
  for(t in (P+1):T){
    prd[t] = mu;
    for(p in 1:P){
      prd[t] += alpha[p] * y[t-p];
    }
    for(q in 1:min(Q,t-1)){
      prd[t] += beta[q] * (y[t-q] - prd[t-q]);
    }
  }
}

model {
  mu ~ normal(0, 10);
  alpha ~ normal(0, 2);
  beta ~ normal(0, 2);
  sigma ~ cauchy(0, 5);
  for(t in 1:T){
    y[t] ~ normal(prd[t], sigma);
  }
}