본 글에서는 State Space Model 중 Gaussian State Space Model에 대해 살펴보고자 합니다.
[참고자료]
- Time Series Analysis and Its Applications
- Linear Gaussian State Space Model
- Stanford Lecture Note: Linear Gaussian State Space Model
1. State Space Model(SSM)
상태 공간 모형(State Space Model)은 시계열 데이터 모형으로 크게 상태 방정식(state equation)과 관측 방정식(observation equation)으로 이루어집니다.

▶ 상태 방정식(state equation)
- xt: state process(hidden or latent process)로 관측되지 않는 데이터(""unobserved data"")입니다.
- xt에 대한 가정: Markov Process
Markov Process는 과거와 현재 상태가 주어졌을 때의 미래 상태의 조건부 확률 분포가 과거 상태와는 독립적으로 현재 상태에 의해서만 결정된다는 것을 의미합니다. 즉, t+1 시점은 오직 t 시점에만 의존합니다(dependent).
P(xt+1|xt,xt−1,⋯,x0)=P(xt+1|xt)
따라서, xt에 대한 모형 설정으로 AR(1) 모형 혹은 ARMA(1,1) 모형 등을 생각할 수 있다.
본 글에서는 ARMAX(1,1) 모형으로 설정하겠습니다.
※ 시계열 통계 모형 참고자료(AR, MA, ARMA, ARMAX)
xt=Φxt−1+γut+wt
▶ 관측 방정식(observation equation)
yt=Atxt+Γut+vt
yt는 우리가 실제로 관측한 데이터입니다.
우리는 state vector xt∈Rp에 대해서 직접적으로 관측할 수는 없지만,
observation vector yt∈Rq를 통해 간접적으로 관측할 수 있습니다.
▶ 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 wt와 vt가 정규분포(Normal distribution)을 따르는 경우,
우리는 Gaussian State Space Model이라고 지칭합니다.
2. Gaussian State Space Model(GSSM)
▶ 상태 방정식(state equation)
xt=Φxt−1+γut+wt∈Rp
- wtiid∼Np(0,Q)∈Rp,Q∈Rp×p
- ut∈Rp: fixed input, exogenous variables(파생변수)
- x0∼Np(μ0,Σ0)
정규분포의 성질에 의해 xt,t∈{1,2,⋯T}는 모두 정규분포를 따릅니다.
▶ 관측 방정식(observation equation)
yt=Atxt+Γut+vt∈Rq
- At∈Rq×p: measurement or observation matrix
- vtiid∼Nq(0,R)∈Rq,R∈Rq×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
인덱스 s와 t에 따른 다음 3가지의 문제를 정의하고자 합니다.
- s<t인 경우: forecasting or prediction
- s=t인 경우: filtering
- 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:t∼N
- xt+1|yt+1∼N
- xt+1|y1:T∼N
▶ 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+1∼N(?,?)
▶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
- wt⊥y1:t→E[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]
- xt⊥wt→Cov(xt,wt|y1:t)=0
- wt⊥y1:t→Var(wt|y1:t)=Var(wt)=Q
▶ Distribution of xt+1|y1:t
xt+1|y1:t∼N(ˆ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+1−At+1ˆxt+1|t)(At+1xt+1+vt+1−At+1ˆxt+1|t)T|y1:t]=E[(At+1xt+1−At+1ˆxt+1|t)(At+1xt+1−At+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+1−At+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+1−At+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+Q−Pt+1|tATt+1(ATtPt+1|tATt+R)−1At+1Pt+1|t
[3]. conditional distribution
xt+1|yt+1∼N(ˆ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+1−At+1ˆxt+1|t+Γut+1)
- Pt+1|t+1=ΦPt+1|tΦT+Q−Pt+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+1−At+1ˆxt+1|t+Γut+1)
Cov(xt+1,xt+1|yt+1)=ΦPt+1|tΦT+Q−Pt+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=(P−1t+1|t+ATt+1RAt+1)−1ATt+1R−1
2.2 (Rausch - Tung - Streibel) Smoothing
▶Define problem
Given, xt+1|y1:T,calculate distribution of xt|y1:T,t<T
xt|y1:T∼N(?,?)
▶Two property
(a). xt⊥yt+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|yt∼N([ˆ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ΦTP−1t+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|t−LtPt+1|tLTt
[3]. calculate conditional distribution
xt|xt+1,y1:t∼N(ˆxt|t+Lt(xt+1−ˆxt+1|t),Pt|t−LtPt+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|t−LtPt+1|tLTt|y1:T]=LtPt+1|tLTt+Pt|t−LtPt+1|tLTt=Pt|t+Lt(Pt+1|T−Pt+1|t)LTt
[3].conditional distribution
xt|y1:T∼N(ˆxt|T,Pt|T)
- ˆxt|T=ˆxt|t+Lt(ˆxt+1|T−ˆxt+1|t)
- Pt|T=Pt|t+Lt(Pt+1|T−Pt+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 |