Loading [MathJax]/jax/output/CommonHTML/jax.js
본문 바로가기

Time Series Analaysis/Time Series Analysis

Linear Gaussian State Space Model


본 글에서는 State Space Model 중 Gaussian State Space Model에 대해 살펴보고자 합니다.

 

[참고자료]


1. State Space Model(SSM)

상태 공간 모형(State Space Model)은 시계열 데이터 모형으로 크게 상태 방정식(state equation)관측 방정식(observation equation)으로 이루어집니다. 

[Figure 1] State Space Model

▶ 상태 방정식(state equation)

  • xt:  state process(hidden or latent process)로 관측되지 않는 데이터(""unobserved data"")입니다.
  • xt에 대한 가정: Markov Process 

Markov Process과거와 현재 상태가 주어졌을 때의 미래 상태의 조건부 확률 분포가 과거 상태와는 독립적으로 현재 상태에 의해서만 결정된다는 것을 의미합니다. 즉, t+1 시점은 오직 t 시점에만 의존합니다(dependent).

P(xt+1|xt,xt1,,x0)=P(xt+1|xt)

 

따라서, xt에 대한 모형 설정으로 AR(1) 모형 혹은 ARMA(1,1) 모형 등을 생각할 수 있다.

본 글에서는 ARMAX(1,1) 모형으로 설정하겠습니다.

※ 시계열 통계 모형 참고자료(AR, MA, ARMA, ARMAX)

xt=Φxt1+γut+wt

 

▶ 관측 방정식(observation equation)

yt=Atxt+Γut+vt

yt는 우리가 실제로 관측한 데이터입니다.
우리는 state vector xtRp에 대해서 직접적으로 관측할 수는 없지만,
observation vector ytRq를 통해 간접적으로 관측할 수 있습니다. 

 

▶ State Space Model의 목표

State Space Model의 목표는 observation vector y1:s={y1,y2,,ys}를 기반으로 xt를 추정하는 것입니다.

  • s<t
  • s=t
  • s>t
A primary aim of state space model, would be to produce estimators for the underlying unobserved signal xt, given the data y1:s={y1,,yn}(⇒ observed data)

 

Noise vector wtvt가 정규분포(Normal distribution)을 따르는 경우,
우리는 Gaussian State Space Model이라고 지칭합니다.


2. Gaussian State Space Model(GSSM)

▶ 상태 방정식(state equation)

xt=Φxt1+γut+wtRp

  • wtiidNp(0,Q)Rp,QRp×p
  • utRp: fixed input, exogenous variables(파생변수)
  • x0Np(μ0,Σ0)

정규분포의 성질에 의해 xt,t{1,2,T}는 모두 정규분포를 따릅니다.

 

▶ 관측 방정식(observation equation)

yt=Atxt+Γut+vtRq

  • AtRq×p: measurement or observation matrix
  • vtiidNq(0,R)Rq,RRq×q

마찬가지로 정규분포의 성질에 의해 yt,t{1,2,,T}도 정규분포를 따릅니다.

 

▶ Assumption

  • x0,{wt},{vt}는 모두 독립이다.
    • 이 가정은 필수는 아닙니다.

▶ Estimation of parameter θ

현재 우리가 추정해야 하는 파라미터는 θ=(Φ,γ,At,Γ,Q,Σ0,R)입니다.

파라미터 추정 방법으로 Expectation Maximize(EM) 알고리즘과 Markov Chain Monte Carlo(MCMC) 방법이 존재합니다.

현재 모수 θ 추정까지 완료된 상태에서
State Space Model의 목적인 y1:s={y1,y2,,ys}가 주어졌을 때 xt를 추정하는 방법에 대해 살펴보고자 합니다.

 

▶ Filtering v.s. Smoothing v.s. Forecasting

인덱스 st에 따른 다음 3가지의 문제를 정의하고자 합니다.

  1. s<t인 경우: forecasting or prediction
  2. s=t인 경우: filtering 
  3. s>t인 경우: smoothing

즉, P(xt+1|y1:t)prediction  distribution, P(xt+1|yt+1)filtered distribution, P(xt+1|y1:T)smoothed distribution이라고 지칭하겠습니다.

 

또한, 우리는 식 (1)과 식 (2)에서 정규분포를 가정했기 때문에 prediction distribution, filtered distribution, smoothed distribution 모두 정규분포를 따릅니다.

  • xt+1|y1:tN
  • xt+1|yt+1N
  • xt+1|y1:TN

▶ Kalman Filter v.s Smoothing

[1]. Kalman Filter

이전 관측치 y1:t를 사용하지 않고,
새로운 데이터 yt+1가 주어졌을 때, filter distribution P(xt|yt)P(xt+1|yt+1)로 업데이트 할 때 Kalman Filter 방법론을 사용합니다.

xt+1|yt+1가 정규분포를 따르는 것은 알고 있으며 정규분포의 모수인 평균분산을 Kalman Filter를 통해 구하고자 하는 것

 

[2]. (Rausch - Tung - Streinbel, RTS) Smoothing

P(xt+1|y1:T)가 주어졌을 때, P(xt|y1:T)를 예측할 때 Smoothing 방법론을 사용합니다.

즉, 이전 hidden state에 대해 예측하는 것입니다.

xt+1|y1:T의 평균과 분산을 구하고자 하는 것

 

먼저, Kalman Filter에 대해 알아보겠습니다.


2.1 Kalman Filter

▶Define problem

Calculate distribution of xt+1|yt+1

xt+1|yt+1N(?,?)

 

▶Kalman Filter의 구성요소

Kalman Filter는 크게 두 단계로 이루어져 있습니다.

  • "Time Update"
  • "Measurement Upate"

2.1.1 Time Update

▶ Purpose

마지막으로 업데이트 된 filtered distribution P(xt|y1:t)가 주어졌을 때,

prediction distribution P(xt+1|y1:t)을 계산하는 것이 Time Update의 목적입니다.

 

▶ Notation

xt|y1:s에 대한 기댓값과 공분산을 다음과 같이 명시하겠습니다.

ˆxt|s=E(xt|y1:s)

Pt1,t2|s=E[(xt1|sˆxt1|s)(xt2|sˆxt2|s)|y1:s]=E[(xt1|sˆxt1|s)(xt2ˆxt2|s)|y1:s]

 

▶ Compute mean of xt+1|y1:t

ˆxt+1|t=E[xt+1|y1:t]=E[Φxt+γut+wt]((1))=E[Φxt|y1:t]+E[γut|y1:t]+E[wt|y1:t]=Φˆxt|t+γut

  • ˆxt|t=E[xt|y1:t]
  • ut는 고정된 값이므로, E[γut|y1:t]=γut
  • wty1:tE[wt|y1:t]=E[wt]=0

▶ Compute covariance of xt+1|y1:t

Pt+1|t=Var(xt+1|y1:t)=Var(Φxt+γut+wt|y1:t)=Var[Φxt|y1:t]+Var[wt|y1:t]=ΦVar[xt|y1:t]ΦT+Q=ΦPt|tΦT+Q

  • Pt|t=Var[xt|y1:t]
  • xtwtCov(xt,wt|y1:t)=0
  • wty1:tVar(wt|y1:t)=Var(wt)=Q

▶ Distribution of xt+1|y1:t

xt+1|y1:tN(ˆxt+1|t,Pt+1|t)

  • ˆxt+1=Φˆxt|t+γut
  • Pt+1|t=ΦPt|tΦT+Q

2.1.2 Measurement Update

▶ Purpose

마지막으로 업데이트 된 prediction distribution P(xt+1|y1:t)가 주어졌을 때,

filtered distribution P(xt+1|yt+1)을 계산하는 것이 Measurement Update의 목적입니다.

 

▶ Step

Measurement Update는 크게 2단계로 이루어집니다.

  • calculate joint distribution of yt+1,xt+1|y1:t
  • calculate conditional distribution of xt+1|yt+1

First step: calculate the joint distribution of yt+1,xt+1|y1:t

우리가 계산해야 하는 값은 다음과 같습니다.

  • E[yt+1|y1:t], Cov(yt+1,yt+1|y1:t) unknown value
  • E[xt+1|y1:t]=ˆxt+1|t, Cov(xt+1,xt+1|y1:t)=Pt+1|t → known value(Time Update에서 구한 값)
  • Cov(yt+1,xt+1|y1:t)  → unknown value

[1]. calculate the E[yt+1|y1:t], Cov(yt+1,yt+1|y1:t)

ˆyt+1|t=E[yt+1|y1:t]=E[At+1xt+1+Γut+1+vt+1|y1:t]=At+1E[xt+1|y1:t]+Γut+1=At+1ˆxt+1|t+Γut+1

 

Cov(yt+1,yt+1|y1:t)=E[(yt+1ˆyt+1|t)(yt+1ˆyt+1|t)T|y1:t]=E[(At+1xt+1+vt+1At+1ˆxt+1|t)(At+1xt+1+vt+1At+1ˆxt+1|t)T|y1:t]=E[(At+1xt+1At+1ˆxt+1|t)(At+1xt+1At+1ˆxt+1|t)T|y1:t]+E[vt+1vt+1|y1:t]=ATt+1Pt+1|tATt+1+R

[2]. calcualte the Cov(yt+1,xt+1|y1:t)

Cov(yt+1,xt+1|y1:t)=E[(yt+1ˆyt+1|t)(xt+1ˆxt+1|t)|y1:t]=E[(At+1xt+1+vt+1At+1ˆxt+1|t)(xt+1ˆxt+1|t)T|y1:t]=At+1E[(xt+1ˆxt+1|t)(xt+1ˆxt+1|t)T|y1:t]=At+1Pt+1|t

 

[3]. calculate joint distribution

[xt+1yt+1]N([ˆxt+1|tAt+1ˆxt+1|t+Γut+1],[ΦPt|tΦT+QPt+1|tATt+1At+1Pt+1|tATtPt+1|tATt+R])

 

Second step: calculate the conditional distribution of xt+1|yt+1

정규분포의 결합 분포(joint distribution)에서 조건부 분포(conditional distribution)로 변환하는 공식을 이용합니다.

※ 다변수 정규분포 참고자료 → 식 (8.6.22)와 식 (8.6.25) 참고 

 

[1]. conditional mean of xt+1|yt+1

E[xt+1|yt+1]=ˆxt+1|t+Pt+1|tATt+1(ATtPt+1|tATt+R)1(yt+1At+1ˆxt+1|t+Γut+1)

 

[2]. conditional covariance of xt+1|yt+1

Cov(xt+1,xt+1|yt+1)=ΦPt+1|tΦT+QPt+1|tATt+1(ATtPt+1|tATt+R)1At+1Pt+1|t

 

[3]. conditional distribution

xt+1|yt+1N(ˆxt+1|t+1,Pt+1|t+1)

  • ˆxt+1|t+1=ˆxt+1|t+Pt+1|tATt+1(ATtPt+1|tATt+R)1(yt+1At+1ˆxt+1|t+Γut+1)
  • Pt+1|t+1=ΦPt+1|tΦT+QPt+1|tATt+1(ATtPt+1|tATt+R)1At+1Pt+1|t

 

우리가 Kalman filter를 통해 최종적으로 구하고자 했던 xt+1|yt+1의 분포를 구하였습니다.


▶ Algorithms of Kalman Filter

Kalman Filter의 알고리즘을 요약하면 다음과 같습니다.

🔶 Algorithm of Kalmann Filter
1. Initialize with ˆx0|1=μ0,P0|1=Σ0
2. Time update
ˆxt+1|t=Φˆxt|t+γut
Pt+1|t=ΦPt|tΦT+Q
3. Measurement Update
    
    E[xt+1|yt+1]=ˆxt+1|t+Pt+1|tATt+1(ATtPt+1|tATt+R)1(yt+1At+1ˆxt+1|t+Γut+1)
    
    Cov(xt+1,xt+1|yt+1)=ΦPt+1|tΦT+QPt+1|tATt+1(ATtPt+1|tATt+R)1At+1Pt+1|t

 

[참고]

🔶 Kalman Gain Matrix

Kt+1=Pt+1|tATt+1(ATtPt+1|tATt+R)1

Kt+1을 Kalman Gain Matrix라고 불리우며 Sherman-Morrison-Woodbury 이론에 따라서 다음과 같이 Kt+1을 표현할 수 있습니다.

Kt+1=(P1t+1|t+ATt+1RAt+1)1ATt+1R1

 


2.2 (Rausch - Tung - Streibel) Smoothing

▶Define problem

Given, xt+1|y1:T,calculate distribution of xt|y1:T,t<T

xt|y1:TN(?,?)

 

▶Two property

(a). xtyt+1|T|xt+1
(b). (a)번의 성질에 의해, p(xt|xt+1,y1:T)=p(xt|xt+1,y1:t)이 성립합니다.

 

Calculate distribution of xt|y1:T

크게 3단계로 이루어집니다.

  • calculate the joint distribution of $x_{t}, x_{t+1}|y_{1:t}( = x_{t}, x_{t+1}|y_{1:T})$
  • calcualte the conditional distribution of xt|xt+1,y1:t(=xt|xt+1,y1:T)
  • calculate the xt|y1:T using below two property
    • (Tower property) E[Z|X]=E[E[Z|Y,X]|X]
    • (Law of Total Conditional Variance) Cov[Z|X]=Cov[E[Z|Y,X]|X]+E[Cov[Z|Y,X]|X]

 First step: calculate the joint distribution of $x_{t}, x_{t+1}|y_{1:t}( = x_{t}, x_{t+1}|y_{1:T})$

우리가 계산해야 하는 값은 다음과 같습니다.

  • E[xt|y1:t]=ˆxt|t,Cov(xt,xt|y1:t)=Pt|t → known value(definition of (a),(b))
  • E[xt+1|y1:t]=ˆxt+1|t,Cov(xt+1,xt+1|y1:t)=Pt+1|t → known value
  • Cov(xt,xt+1|y1:t) → unknown value

[1]. Calculate the Cov(xt,xt+1|y1:t)

Cov(xt,xt+1|y1:t)=E[(xtˆxt|t)(xt+1ˆxt+1|t)T|y1:t]=E[(xtˆxt|t)(Φxt+wtΦˆxt|t)T|y1:t]=Pt|tΦT

 

[2]. Calculate joint distribution

xt,xt+1|ytN([ˆxt|tˆxt+1|t],[Pt|tPt|tΦTΦPt|tPt+1|t])

 

 Second step: calculate the conditional distribution of xt|xt+1,y1:t(=xt|xt+1,y1:T)

[1]. calculate the conditional mean xt|xt+1,y1:t

E(xt|xt+1,y1:T)=E(xt|xt+1,y1:t)=ˆxt|t+Lt(xt+1ˆxt+1|t)

  • Lt=Pt|tΦTP1t+1|t

[2]. calculate the conditional covariance xt|xt+1,y1:t

Cov(xt|xt+1,y1:T)=Cov(xt|xt+1,y1:t)=Pt|tLtPt+1|tLTt

 

[3]. calculate conditional distribution

xt|xt+1,y1:tN(ˆxt|t+Lt(xt+1ˆxt+1|t),Pt|tLtPt+1|tLTt)

 

Third step: calculate the conditional distribution of xt|y1:T

  • "Tower property"(식 (5))와 " Law of Total Conditional Variance "(식 (6))을 이용하여 조건부 평균과 조건부 공분산을 구합니다.

[1]. calculate the conditional mean xt|y1:T

ˆxt|T=E[xt|y1:T]=E[E[xt|xt+1,y1:T]|y1:T]=E[ˆxt|t+Lt(xt+1ˆxt+1|t)|y1:T]=ˆxt|t+Lt(ˆxt+1|Tˆxt+1|t)

 

[2]. calculate the conditional covariance xt|y1:T

Pt|T=Cov[xt|y1:T]=Cov(E[xt|xt+1,y1:T]|y1:T)+E[Cov[xt|xt+1,y1:T]|y1:T]=Cov[ˆxt|t+Lt(xt+1ˆxt+1|t)|y1:T]+E[Pt|tLtPt+1|tLTt|y1:T]=LtPt+1|tLTt+Pt|tLtPt+1|tLTt=Pt|t+Lt(Pt+1|TPt+1|t)LTt

 

[3].conditional distribution

xt|y1:TN(ˆxt|T,Pt|T)

  • ˆxt|T=ˆxt|t+Lt(ˆxt+1|Tˆxt+1|t)
  • Pt|T=Pt|t+Lt(Pt+1|TPt+1|t)LTt

Algorithm of (Rausch - Tung - Streibel) Smoothing

 (Rausch - Tung - Streibel) Smoothing의 알고리즘을 요약하면 다음과 같습니다.

🔶 Algorithm of (Rausch - Tung - Streibel) Smoothing

[1]. 먼저, Kalman Filter를 통해 filtered와 one-step prediciton quantities인 ˆzt|t,ˆzt+1|t,Pt|t,Pt+1|t을 얻습니다(0,1,,T).

  • filtered quantities : ˆzt|t,Pt|t
  • one-step prediction quantities: ˆzt+1|t,Pt+1|t

[2]. ˆxT|T,PT|T으로 smoother 초기화

[3]. p(xt+1|y1:T)가 주어진 상황에서, p(xt|y1:T)


본 글에서는 observation data y1:s={y1,,ys}가 주어졌을 때, hidden state xt를 추정하는 (Gaussian) State Space Model에 대해 알아보았습니다.


 

'Time Series Analaysis > Time Series Analysis' 카테고리의 다른 글

Convolutional LSTM network(ConvLSTM)  (0) 2024.01.29
Gated Recurrent Unit(GRU)  (2) 2024.01.23
Long Short Term Memory(LSTM)  (0) 2024.01.22
Recurrent Neural Network(RNN)  (0) 2024.01.22