【日次GDP】カルマンフィルタの実装（R）

おはこんばんにちは。かなり久しぶりの投稿となってしまいました。決して研究をさぼっていたのではなく、BVARのコーディングに手こずっていました。あと少しで完成します。さて、今回はBVARやこの前のGiannnone et a (2008)のような分析でも大活躍のカルマンフィルタを実装してしまいたいと思います。このブログではパッケージソフトに頼らず、基本的に自分で一から実装を行い、研究することをポリシーとしていますので、これから頻繁に使用するであろうカルマンフィルタを関数として実装してしまうことは非常に有益であると考えます。今回はRで実装をしましたので、そのご報告をします。

カルマンフィルタとは？

まず、カルマンフィルタに関する簡単な説明を行います。非常にわかりやすい記事があるので、そちらを読んでいただいたほうがより分かりやすいかと思います。

qiita.com

カルマンフィルタとは、状態空間モデルを解くアルゴリズムの一種です。状態空間モデルとは、手元の観測可能な変数から観測できない変数を推定するモデルであり、以下のような形をしています。

ここで、はg×1ベクトルの観測可能な変数(観測変数)、はg×k係数行列、はk×1ベクトルの観測不可能な変数(状態変数)、はk×k係数行列です。また、は観測変数の誤差項、は状態変数の誤差項です。これらの誤差項はそれぞれに従います（は分散共分散行列）。は定数項です。1本目の式は観測方程式、2本目の式は遷移方程式と呼ばれます。
状態空間モデルを使用する例として、しばしば池の魚の数を推定する問題が使用されます。今、池の中の魚の全数が知りたいとして、その推定を考えます。観測時点毎に池の中の魚をすべて捕まえてその数を調べるのは現実的に困難なので、一定期間釣りをして釣れた魚をサンプルに全数を推定することを考えます。ここで、釣れた魚は観測変数、池にいる魚の全数は状態変数と考えることができます。今、経験的に釣れた魚の数と全数の間に以下のような関係があるとします。

これが観測方程式になります。また、魚の全数は過去の値からそれほど急速には変化しないと考えられるため、以下のようなランダムウォークに従うと仮定します。

これが遷移方程式になります。あとは、これをカルマンフィルタアルゴリズムを用いて計算すれば、観測できない魚の全数を推定することができます。
このように状態空間モデルは非常に便利なモデルであり、また応用範囲も広いです。例えば、販売額から潜在顧客数を推定したり、クレジットスプレッドやトービンのQ等経済モデル上の概念として存在する変数を推定する、BVARのようにVARや回帰式の時変パラメータ推定などにも使用されます。

カルマンフィルタアルゴリズムの導出

さて、非常に便利な状態空間モデルの説明はこれくらいにして、カルマンフィルタの説明に移りたいと思います。カルマンフィルタは状態空間モデルを解くアルゴリズムの一種であると先述しました。つまり、他にも状態空間モデルを解くアルゴリズムは存在します。カルマンフィルタアルゴリズムは一般に誤差項の正規性の仮定を必要としないフィルタリングアルゴリズムであり、観測方程式と遷移方程式の線形性の仮定さえあれば、線形最小分散推定量となります。カルマンフィルタアルゴリズムの導出にはいくつかの方法がありますが、今回はこの線形最小分散推定量としての導出を行います。まず、以下の３つの仮定を置きます。

①初期値は正規分布に従う確率ベクトルである(はの推定値)。
②誤差項は全てので互いに独立で、初期値ベクトルと無相関である（）。
③②より、、。

まず、期の情報集合が既知の状態でのとの期待値（予測値）を求めてみましょう。上述した状態空間モデルと誤差項の期待値がどちらもゼロである事実を用いると、以下のように計算することができます。

ここで、次に、これらの分散を求めます。

ここで、期の情報集合が得られたとします（つまり、観測値を入手）。カルマンフィルタでは、期の情報である観測値を用いてを以下の方程式で更新します。

つまり、観測値との期待値（予測値）の差をあるウェイト（k×g行列）でかけたもので補正をかけるわけです。よって、観測値と予測値が完全に一致していた場合は補正は行われないことになります。ここで重要なのは、ウエイトをどのように決めるのかです。は更新後の状態変数の分散を最小化するよう決定します。これが、カルマンフィルタが線形最小分散推定量である根拠です。最小化にあたっては以下のベクトル微分が必要になりますので、おさらいをしておきましょう。今回使用するのは以下の事実です。

ここで、はベクトル（それぞれn×1ベクトル、1×nベクトル）、はn×nの対称行列です。まず、１つ目から証明していきます。とします。
このとき、なので、

次に２つ目です。とします。このとき、

よって、

さて、準備ができたので、更新後の状態変数の分散を求めてみましょう。

１回目の式変形で、に上述した更新式を代入し、２回目の式変形で観測方程式と上で計算したを代入しています。さて、更新後の状態変数の分散をの関数として書き表すことができたので、これをで微分し、0と置き、を最小化するを求めます。先述した公式で、、、とすると（は分散共分散行列の和なので対称行列）、

ここから、を解きなおすと、

突然、が出てきました。これは観測変数の予測値の分散でした。一方、は状態変数の予測値の分散を観測変数のスケールに調整したものです（観測空間に写像したもの）。つまり、カルマンゲインは状態変数と観測変数の予測値の分散比となっているのです。観測変数にはノイズがあり、観測方程式はいつも誤差０で満たされるわけではありません。また、状態方程式にも誤差項が存在します。状態の遷移も100%モデル通りにはいかないということです。期の観測変数が得られたとして、それをどれほど信頼して状態変数を更新するかは観測変数のノイズが状態変数のノイズに比べてどれほど小さいかによります。つまり、相対的に観測方程式が遷移方程式よりも信頼できる場合には状態変数を大きく更新するのです。このように、カルマンフィルタでは、観測方程式と遷移方程式の相対的な信頼度によって、更新の度合いを毎期調整しています。その度合いが分散比であり、カルマンゲインだというわけです。ちなみに欠損値が発生した場合には、観測変数の分散を無限大にし、状態変数の更新を全く行わないという対処を行います。観測変数に信頼がないので当たり前の処置です。この場合は遷移方程式を100%信頼します。これがカルマンフィルタのコアの考え方になります。
更新後の分散を計算しておきます。

では、最終的に導出されたアルゴリズムをまとめたいと思います。

初期値が所与の元で、まず状態変数の予測値を計算します。その結果を用いて、次は観測変数の予測値を計算し、カルマンゲインを得ます。期の観測可能なデータを入手したら、更新方程式を用いてを更新します。これをサンプル期間繰り返していくことになります。ちなみに、遷移方程式の誤差項と定数項がなく、遷移方程式のパラメータが単位行列のカルマンフィルタは逐次最小自乗法と一致します。つまり、新しいサンプルを入手するたびにOLSをやり直す推計方法ということです（今回はその証明は勘弁してください）。

Rで実装する。

以下がRでの実装コードです。

#Kalman filter

kalmanfiter <- function(y,I,t,z,c=0,R=NA,Q=NA,d=0,S=NA,h=NA,a_int=NA,sig_int=NA){
  #-------------------------------------------------------------------
  # Implemention of Kalman filter
  #   y - observed variable
  #   I - the number of unobserved variable
  #   t - parameter of endogenous variable in state equation
  #   z - parameter of endogenous variable in observable equation
  #   c - constant in state equaion
  #   R - parameter of exogenous variable in state equation
  #   Q - var-cov matrix of exogenous variable in state equation
  #   d - constant in observable equaion
  #   S - parameter of exogenous variable in observable equation
  #   h - var-cov matrix of exogenous variable in observable equation
  #   a_int - initial value of endogenous variable
  #   sig_int - initial value of variance of endogenous variable
  #-------------------------------------------------------------------
  
  library(MASS)
    
  # 1.Define Variable
  if (class(y)!="matrix"){
    y <- as.matrix(y)
  }
  N <- NROW(y) # sample size
  L <- NCOL(y) # the number of observable variable 
  a_pre <- array(0,dim = c(I,1,N)) # prediction of unobserved variable
  a_fil <- array(0,dim = c(I,1,N+1)) # filtered of unobserved variable
  sig_pre <- array(0,dim = c(I,I,N)) # prediction of var-cov mat. of unobserved variable
  sig_fil <- array(0,dim = c(I,I,N+1)) # filtered of var-cov mat. of unobserved variable
  y_pre <- array(0,dim = c(L,1,N)) # prediction of observed variable
  F_pre <- array(0,dim = c(L,L,N)) # auxiliary variable 
  k <- array(0,dim = c(I,L,N)) # kalman gain
  
  if (is.na(a_int)==TRUE){
    a_int <- matrix(0,nrow = I,ncol = 1)
  }
  if (is.na(sig_int)==TRUE){
    sig_int <- diag(1,nrow = I,ncol = I)
  }
  if (is.na(R)==TRUE){
    R <- diag(1,nrow = I,ncol = I)
  }
  if (is.na(Q)==TRUE){
    Q <- diag(1,nrow = I,ncol = I)
  }
  if (is.na(S)==TRUE){
    S <- matrix(1,nrow = L,ncol = L)
  }
  if (is.na(h)==TRUE){
    H <- array(0,dim = c(L,L,N))
    for(i in 1:N){
      diag(H[,,i]) = 1
    }
  }else if (class(h)!="array"){
    H <- array(h,dim = c(NROW(h),NCOL(h),N))
  }
  
  # fill infinite if observed data is NA
  for(i in 1:N){
    miss <- is.na(y[i,])
    diag(H[,,i])[miss] <- 1e+32
  }
  y[is.na(y)] <- 0  

  # 2.Set Initial Value
  a_fil[,,1] <- a_int
  sig_fil[,,1] <- sig_int
  
  # 3.Implement Kalman filter
  for (i in 1:N){
    if(class(z)=="array"){
      Z <- z[,,i]
    }else{
      Z <- z
    }
    a_pre[,,i] <- t%*%a_fil[,,i] + c
    sig_pre[,,i] <- t%*%sig_fil[,,i]%*%t(t) + R%*%Q%*%t(R)
    y_pre[,,i] <- Z%*%a_pre[,,i] + d
    F_pre[,,i] <- Z%*%sig_pre[,,i]%*%t(Z) + S%*%H[,,i]%*%t(S)
    k[,,i] <- sig_pre[,,i]%*%t(Z)%*%ginv(F_pre[,,i])
    a_fil[,,i+1] <- a_pre[,,i] + k[,,i]%*%(y[i,]-y_pre[,,i])
    sig_fil[,,i+1] <- sig_pre[,,i] - k[,,i]%*%F_pre[,,i]%*%t(k[,,i])
  }
  
  # 4.Aggregate results
  result <- list(a_pre,a_fil,sig_pre,sig_fil,y_pre,k)
  names(result) <- c("state prediction", "state filtered", "state var prediction", 
                     "state var filtered", "observable prediction", "kalman gain")
  return(result)
}

案外簡単に書けるもんですね。これを使って、Giannone et al (2008)をやり直してみます。過去記事はこちら。

osashimix.hatenablog.com

データセットは前回記事と変わりません。データの収集方法はこちら。

osashimix.hatenablog.com

以下、分析用のRコードです。

#------------------------
# Giannone et. al. 2008 
#------------------------

# ファクターを計算
f <- 3　# the number of factors
z <- scale(dataset1) # データセットを標準化
for (i in 1:nrow(z)){
  eval(parse(text = paste("S_i <- z[i,]%*%t(z[i,])",sep = "")))
  if (i==1){
    S <- S_i
  }else{
    S <- S + S_i
  }
}
S <- (1/nrow(z))*S
gamma <- eigen(S)
D <- diag(gamma$values[1:f])
V <- gamma$vectors[,1:f] #固有ベクトルを抽出
F_t <- matrix(0,nrow(z),f)
for (i in 1:nrow(z)){
  eval(parse(text = paste("F_t[",i,",]<- z[",i,",]%*%V",sep = ""))) #主成分分析にてファクターを計算（カルマンフィルタの初期値に使用）
}
lambda_hat <- V #観測方程式の係数行列
R <- diag(diag(cov(z-z%*%V%*%t(V)))) #観測方程式の分散共分散行列 

a <- matrix(0,f,f)
b <- matrix(0,f,f)
for(t in 2:nrow(z)){
  a <- a + F_t[t,]%*%t(F_t[t-1,])
  b <- b + F_t[t-1,]%*%t(F_t[t-1,])
}
b_inv <- solve(b)
A_hat <- a%*%b_inv #VAR推計

e <- numeric(f)
for (t in 2:nrow(F_t)){
  e <- e + F_t[t,]-F_t[t-1,]%*%A_hat
}
H <- t(e)%*%e
Q <- diag(1,f,f)
Q[1:f,1:f] <- H #遷移方程式の分散共分散行列

a <- which(dataset$publication == "1980-01-01")
dataset2 <- dataset[a:nrow(dataset),]
rownames(dataset2) <- dataset2$publication
dataset2 <- dataset2[,-2]
zz <- scale(dataset2)
p <- ginv(diag(nrow(kronecker(A_hat,A_hat)))-kronecker(A_hat,A_hat))

result1 <- kalmanfiter(zz,f,A_hat,V,c=0,R=NA,Q=Q,d=0,S=NA,h=psi,a_int=F_t[1,],sig_int=matrix(p,f,f))
F_tk <- xts(t(result1$`state filtered`[,1,305:469]),order.by = as.Date(rownames(dataset2[305:469,])))
plot.zoo(F_tk,col = c("red","blue","green"),plot.type = "single",ylab = "factor",xlab = "time")

出力されたグラフが以下です。

giannoneの記事を書いた際は、元論文のMATLABコードを参考にRで書いたのですが、通常のカルマンフィルタとは観測変数の分散共分散行列の逆数の計算方法が違うらしくグラフの形が異なっています。まあでも、概形はほとんど同じですが（なので、ちゃんと動いているはず）。

カルマンスムージング

カルマンフィルタの実装は以上で終了なのですが、誤差項の正規性を仮定すれば期までの情報集合を用いて、をへ更新することができます。これをカルマンスムージングと呼びます。これを導出してみましょう。その準備として、以下のようなとの混合分布を計算しておきます。

ここで、は

ここで、条件付き多変量正規分布は以下のような分布をしていることを思い出しましょう（
多変量正規分布における条件付き確率の式と意味 - 具体例で学ぶ数学
）。

これを用いて、を計算してみましょう。

ただし、の値は観測不可能なので、上式を用いて状態変数を更新することはできません。今、わかるのは期におけるの分布のみです。ということで、をで代用し、の分布を求めてみます。では、計算していきます。ですが、

というように、も確率変数となるので、繰り返し期待値の法則と繰り返し分散の法則を使用します（こちらを参照）。

繰り返し期待値の法則

繰り返し分散の法則

よって、

となります。カルマンスムージングのアルゴリズムをまとめておきます。

カルマンスムージングの特徴的な点は後ろ向きに計算をしていく点です。つまり、期から1期に向けて計算を行っていきます。に関してはそもそもカルマンフィルタを回した時点で計算可能ですが、は期までのデータが手元にないと計算できません。今、期まで観測可能なデータが入手できたとしましょう。すると、２番目の方程式を用いて、を計算します。ちなみにはカルマンフィルタを回した時点ですでに手に入っているので、計算する必要はありません。同時に、３番目の式を用いてを計算します。そして、とを用いてを計算、というように1期に向けて後ろ向きに計算をしていくのです。さきほど、遷移方程式の誤差項と定数項がなく、遷移方程式のパラメータが単位行列のカルマンフィルタは逐次最小自乗法と一致すると書きましたが、カルマンスムージングの場合は期までのサンプルでOLSを行った結果と一致します。
Rで実装してみます。

kalmansmoothing <- function(filter){
  #-------------------------------------------------------------------
  # Implemention of Kalman smoothing
  #   t - parameter of endogenous variable in state equation
  #   z - parameter of endogenous variable in observable equation
  #   a_pre - prediction of state
  #   a_fil - filtered value of state
  #   sig_pre - prediction of var of state
  #   sig_fil - filtered value of state
  #-------------------------------------------------------------------
  
  library(MASS)
  
  # 1.Define variable
  a_pre <- filter$`state prediction`
  a_fil <- filter$`state filtered`
  sig_pre <- filter$`state var prediction`
  sig_fil <- filter$`state var filtered`
  t <- filter$`parameter of state eq`
  C <- array(0,dim = dim(sig_fil))
  a_sm <- array(0,dim = dim(a_fil))
  sig_sm <- array(0,dim = dim(sig_fil))
  N <- dim(C)[3]
  a_sm[,,N] <- a_fil[,,N]
  sig_sm[,,N] <- sig_fil[,,N]
  
  for (i in dim(C)[3]:2){
    C[,,i-1] <- sig_fil[,,i-1]%*%t(t)%*%ginv(sig_pre[,,i])
    a_sm[,,i-1] <- a_fil[,,i-1] + C[,,i-1]%*%(a_sm[,,i]-a_pre[,,i])
    sig_sm[,,i-1] <- sig_fil[,,i-1] + C[,,i-1]%*%(sig_sm[,,i]-sig_pre[,,i])%*%t(C[,,i-1])
  }
  
  
  result <- list(a_sm,sig_sm,C)
  names(result) <- c("state smoothed", "state var smoothed", "c")
}

先ほどのコードの続きでRコードを書いてみます。

result2 <- kalmansmoothing(result1)

かなりシンプルですね。ちなみにグラフにしましたが、１個目とほぼ変わりませんでした。とりあえず、今日はここまで。