기초에서 알 수있는 시계열 분석 3. 칼만 필터 (로컬 레벨 모델 구현)

전회：기초로부터 아는 시계열 분석 2.칼만 필터(수리적인 도출)

3.1 모델 정의

로컬 레벨 모델이란, 다음과 같은 랜덤 워크를 나타내는 상태 공간 모델을 말한다. 천이 이행 열이 단위 행렬.

\begin{cases}
{\bf x}_t&={\bf x}_{t-1}+{\bf x}_t\ \ \ \ {\bf w}_t\sim\mathcal{N}({\bf 0},W)\\
{\bf y}_t&={\bf x}_t+{\bf v}_t\ \ \ \ \ {\bf v}_t\sim\mathcal{N}({\bf 0},V)
\end{cases}

#---
#　モデルの定義
#---

mod = dlmModPoly(order = 1)

dlm 라이브러리의 dlmModPoly() 함수로 모델을 정의했다. 이 함수는 다항식 모델을 설정하는 것으로, order=1로 하는 것으로 차수 1의 다항식 모델(로컬 레벨 모델)이 된다. 그 밖에도 시스템 노이즈나 관측 노이즈, 사전 분포도 설정할 수 있다. 시스템 노이즈와 관측 노이즈는 모델의 파라미터이므로,

3.2 파라미터 식별

시스템 노이즈와 관측 노이즈를 최대 우도 추정한다.

#---
#　パラメータの特定
#---

build_dlm = function(par){
  # モデルのパラメータを設定する関数

  mod$W[1,1] = exp(par[1]) # 負の領域を探索しないよう指数変換を行う
  mod$V[1,1] = exp(par[2])

  return(mod)
}

fit_dlm = dlmMLE(y = Nile, parm = c(0,0), build = build_dlm)
mod     = build_dlm(fit_dlm$par)

모델을 정의, 구축하는 함수 build_dlm을 만들어, dlmMLE로 파라미터의 최대 우도 추정치를 탐색시켜, 추정 결과를 기초로 모델을 구축하고 있다. fit_dlm$par이 돌려주는 값은 지수 변환 전의 값이므로 그대로 대입해도 좋다.

3.3 필터링

필터링을 수행한다.

#---
#　フィルタリングの実施
#---

dlmFiltered_obj = dlmFilter(y = Nile, mod = mod)

# フィルタ分布の平均、標準偏差の導出
m               = dropFirst(dlmFiltered_obj$m) # 先頭の事前分布の平均をドロップ
m_sdev          = sqrt(
  dropFirst( # 先頭を事前分布の分散をドロップ
    as.numeric(
      dlmSvd2var(dlmFiltered_obj$U.C, dlmFiltered_obj$D.C) # 特異値を分散に戻す
    ))
  )
m_quant         = list(m + qnorm(0.025, sd = m_sdev), m + qnorm(0.975, sd = m_sdev)) # 95%信頼区間

# 結果のプロット
oldpar = par(no.readonly = T) 
par(family = "HiraginoSans-W4")
ts.plot(cbind(Nile, m, do.call('cbind', m_quant)),
        col = c('lightgray', 'black', 'black', 'black'),
        lty = c('solid', 'solid', 'dashed', 'dashed'))
legend(legend = c('観測値', '平均(フィルタリング分布)', '95%区間'),
       lty    = c('solid', 'solid', 'dashed', 'dashed'),
       col    = c('lightgray', 'black', 'black', 'black'),
       x      = 'topright', text.width = 32, cex = 0.6)
par(oldpar)

dlmFilter() 함수로 필터링을 실시한다. 이 함수는 관측치 y와 모델 mod 외에

m: 필터 분포의 평균

U.C.,D.C.: 필터 분포의 공분산 행렬의 특이값 분해. dlmSvd2var에서 공분산 행렬로 돌아갑니다.

등을 반환합니다.

3.4 예측

10기 앞까지의 예측을 실시한다.

#---
#　予測
#---

dlmForecasted_obj = dlmForecast(mod = dlmFiltered_obj, nAhead = 10)

# 予測分布の平均、標準偏差の導出
a                 = ts(data = dlmForecasted_obj$a, start = c(1971,1))
a_sdev            = sqrt( as.numeric( dlmForecasted_obj$R ) )
a_quant = list(a + qnorm(0.025, sd = a_sdev), a + qnorm(0.975, sd = a_sdev))

# 結果のプロット
oldpar = par(no.readonly = T) 
par(family = "HiraginoSans-W4")
ts.plot(cbind(Nile, a, do.call('cbind', a_quant)),
        col = c('lightgray', 'black', 'black', 'black'),
        lty = c('solid', 'solid', 'dashed', 'dashed'))
legend(legend = c('観測値', '平均(予測分布)', '95%区間'),
       lty    = c('solid', 'solid', 'dashed', 'dashed'),
       col    = c('lightgray', 'black', 'black', 'black'),
       x      = 'topright', text.width = 26, cex = 0.6)
par(oldpar)

dlmForecast() 함수로 예측을 실시한다. 이 함수는

a: 예측 분포의 평균 벡터

R : 예측 분포의 공분산 행렬

등을 반환합니다.

다음 번 : 기초로부터 알 수 있는 시계열 분석 4.칼만 필터(로컬 트렌드 모델의 실장)