E8 · 동적 모형 I: 충격의 전파, 누적 분산, 순차 관측#
1. 결국 사라지는 충격이 처음에는 커질 수 있는가#
두 부문의 생산 편차를 기준 생산량 단위로 \(s_t=(s_{1t},s_{2t})^T\)에 담겠습니다. 외부 충격이 없는 한 기간의 변화는
둘째 부문은 매 기간 절반으로 줄지만 그 값이 첫째 부문의 다음 생산에 더해집니다. \(s_0=e_2\)라면 \(s_1=(1,1/2)^T\), \(s_2=(1,1/4)^T\), \(s_3=(3/4,1/8)^T\)입니다. 첫 반응의 길이는 \(\sqrt5/2>1\)입니다. 고윳값은 둘 다 \(1/2\)인데도 처음에는 커집니다.
\(A=aI+N\), \(a=1/2\), \(N^2=0\)이므로 이항 전개에서 두 항만 남아
따라서 \(A^ke_2=(k2^{1-k},2^{-k})^T\)이고 결국 0입니다. \(k\)라는 다항 인자가 초기 증폭을, \(2^{-k}\)가 장기 감쇠를 설명합니다. 첫째 부문의 누적 반응은 \(\sum_{k\ge1}k2^{1-k}=4\), 둘째는 \(\sum_{k\ge0}2^{-k}=2\)입니다. 직접 \((I-A)^{-1}=\begin{pmatrix}2&4\\0&2\end{pmatrix}\)의 둘째 열과 같습니다.
그림 152 가로축은 이산 기간이다. 점을 잇는 선은 시점 사이 동학을 뜻하지 않는다. 아래 패널은 같은 고윳값을 가진 \(A\)와 \((1/2)I\)의 거듭제곱 작용소노름이다.#
2. 시차를 상태로 만들기#
AR(2) 모형 \(y_t=(6/5)y_{t-1}-(8/25)y_{t-2}+\varepsilon_t\)는 \(s_t=(y_t,y_{t-1})^T\)로 쌓으면
동반행렬의 특성식은 \(\lambda^2-(6/5)\lambda+8/25=(\lambda-4/5)(\lambda-2/5)\)입니다. 시차 다항식 \(1-(6/5)z+(8/25)z^2=(1-4z/5)(1-2z/5)\)의 근은 \(5/4,5/2\), 즉 비영 고윳값의 역수입니다.
VAR(\(p\))도 \(p\)개 시차를 쌓으면 상태 크기가 원래 변수 수의 \(p\)배인 VAR(1)이 됩니다. 상태행렬의 고윳값은 단위원 안, 지연다항식의 유한 근은 단위원 밖이어야 안정적입니다. 최고 시차 행렬이 특이하면 영 고윳값과 다항식 차수 감소가 있으므로 모든 근이 유한한 일대일 역수라고 말하지 않습니다. 또한 단위원 위의 복소근 모두가 고윳값 1은 아닙니다.
일반 관측·충격 지도가 \(y_t=Cs_t\), \(s_{t+1}=As_t+B\varepsilon_{t+1}\)이면 \(k\)기간 충격반응은 \(CA^kB\)입니다. \(A\)가 대각화될 때만 단순 지수 모드들의 합으로 쓸 수 있습니다. Jordan 블록이 있으면 앞 예처럼 다항 인자가 붙습니다. 여러 모드가 부호를 바꾸며 합쳐지면 하나의 고윳값에서 계산한 반감기가 전체 반응의 반감기와 같지 않습니다.
3. 안정성을 네 가지 방법으로 읽기#
유한차원에서는 \(\rho(A)<1\), \(A^k\to0\), \(\sum_{k\ge0}\|A^k\|<\infty\), 그리고 어떤 \(P\succ0\)에 \(A^TPA-P\prec0\)가 동치입니다. 마지막 식은 상태의 가중 제곱길이가 매 단계 감소한다는 뜻입니다. 유클리드 길이는 잠시 증가해도 이 가중 길이는 감소할 수 있습니다.
보조정리와 정리 3. 노름의 선택과 Gelfand 공식의 \(\lim_k\|A^k\|^{1/k}=\rho(A)\)는 장기 지수율을 말합니다. \(\|A\|\)나 처음 몇 번의 \(\|A^k\|\)를 \(\rho(A)^k\)로 대체하지 않습니다. \(\rho(A)<q<1\)을 정하면 충분히 큰 \(k\)에서 \(\|A^k\|\le q^k\)이고, 유한한 앞부분을 상수 \(C\)에 흡수하면 모든 \(k\)에 \(Cq^k\) 상한을 얻습니다.
안정성은 잡음이 계속 들어오는 모형의 정상 분산을 만드는 조건이기도 합니다. 다만 초기 상태가 임의이면 처음부터 정상분포가 되는 것은 아닙니다. 다음 절의 공분산으로 초기화하거나 초기 영향이 사라지는 극한을 보아야 합니다.
4. 매 기간 충격의 공분산을 더하기#
\(s_{t+1}=As_t+w_{t+1}\), \(Ew_t=0\), \(\operatorname{Cov}w_t=Q\), 서로 다른 시점 및 초기 상태와 무상관이라고 합시다. 정상 공분산은
1절의 \(A\)와 \(Q=I\)에 \(P=\begin{pmatrix}p&r\\r&s\end{pmatrix}\)를 대입합니다. 곱의 성분은
따라서 \(s=s/4+1\)에서 \(s=4/3\), \(r=r/4+s/2\)에서 \(r=8/9\), \(p=p/4+r+s+1\)에서 \(p=116/27\)입니다. 즉
둘째 부문의 충격이 첫째로 전달되므로 첫째 분산이 더 크고 양의 교차공분산이 생깁니다. 이것은 모형의 전달 구조가 만든 결과입니다.
그림 153 초기 공분산을 0으로 두고 \(P_{k+1}=AP_kA^T+I\)를 반복한다. 점선은 정확한 정상 분산이다. 세로축 단위는 생산 편차 단위의 제곱이다.#
열 순서 vec를 쓰면 \((I-A\otimes A)\operatorname{vec}P=\operatorname{vec}Q\)입니다. 유일해 조건은 모든 고윳값 곱 \(\lambda_i(A)\lambda_j(A)\ne1\)입니다. 안정성은 이를 보장하지만 필요조건은 아닙니다. 예를 들어 스칼라 \(A=2,Q=1\)에는 유일해 \(P=-1/3\)이 있습니다. 대수적 유일해가 공분산으로 허용되는 해라는 뜻은 아닙니다.
직접 Kronecker 선형계는 미지수가 \(d^2\)개여서 조밀 저장 \(O(d^4)\), 일반 소거 \(O(d^6)\)입니다. Schur 기반 Stein 해법은 원래 \(d\times d\) 행렬에서 \(O(d^3)\) 연산과 \(O(d^2)\) 저장으로 처리합니다. 작은 크기의 실행시간으로 이 복잡도를 확정하거나 특정 크기에서 반드시 메모리 실패한다고 예언하지 않습니다. 구현 경로와 컴퓨터에 따라 교차점은 달라집니다.
5. 불완전한 관측을 받으면 무엇을 고치는가#
상태는 직접 보이지 않고 \(y_t=Hs_t+v_t\)를 본다고 하자. 예측 평균 \(m^-\)와 오차 공분산 \(P^-\)가 주어졌고 관측잡음은 평균 0, 과거와 현재 예측오차에 무상관, 공분산 \(R\succ0\)입니다. 새 정보는 관측 자체보다 예측에서 벗어난 혁신
입니다. 상태 오차와 혁신의 교차공분산은 \(P^-H^T\). 최소 선형예측의 계수는 \(K=P^-H^TS^{-1}\)이므로
독립 Gaussian 모형에서는 이것이 조건부 평균·공분산입니다. 비정규에서는 일반적으로 최적 선형 예측이지 전체 조건부 평균은 아닙니다.
수치 예로 \(P^-=\begin{pmatrix}2&1\\1&2\end{pmatrix}\), \(H=(1,0)\), \(R=1\), \(m^-=0\), \(y=3\)을 택하겠습니다. \(S=3\), \(K=(2/3,1/3)^T\)이고 새 평균은 \((2,1)^T\)입니다. 첫 상태만 보았지만 둘째도 교차공분산 때문에 변합니다. 새 공분산은
감소량은 PSD 계수 1 행렬입니다. 첫 상태의 분산은 줄지만 잡음이 있으므로 0은 아닙니다. 이후 예측은 \(m^-_{t+1}=Am_t^+\), \(P^-_{t+1}=AP_t^+A^T+Q\)입니다.
그림 154 타원은 확률질량의 고정 95% 영역이 아니라 \((s-m)^TP^{-1}(s-m)=1\)의 경계다. 점선 수직선은 관측값 3이며 잡음 때문에 사후 평균이 그 선에 놓일 필요는 없다.#
6. 혁신 우도와 순차 직교화#
매 시점 \(\nu_t=y_t-\widehat y_{t|t-1}\)는 과거 관측의 선형공간에 직교합니다. 따라서 이전 혁신과도 무상관입니다. 관측 전체에서 혁신 전체로 가는 변환은 대각이 항등인 블록 하삼각입니다. 행렬식은 1입니다. Gaussian이면 무상관 혁신은 독립이고
여기서는 매 시점 관측 차원을 \(m\)으로 고정했습니다. 결측 시점이나 특이한 혁신 공분산은 별도 처리해야 합니다. 로그행렬식 항은 상태추정의 정확도뿐 아니라 불확실성의 크기까지 우도에 반영합니다.
스칼라 확인으로 \(A=1/2\), \(Q=1\), \(R=1\), 첫 예측분산 \(P_1^-=1\), 관측 \((1,0)\)을 택하겠습니다. 첫 \(S_1=2\), \(K_1=1/2\), \(m_1^+=1/2\), \(P_1^+=1/2\). 둘째 예측은 \(m_2^-=1/4\), \(P_2^-=9/8\)라 \(S_2=17/8\), \(\nu_2=-1/4\)입니다. 결합 관측 공분산은
직접 역을 곱한 이차형식은 \(9/17\). 혁신으로 계산해도 \(1/2+(1/16)/(17/8)=1/2+1/34=9/17\)입니다. 순차 필터와 전체 Gaussian 우도가 같은 계산임을 작은 자료에서 확인했습니다.
7. 관측가능성, 제어가능성, 최소실현#
상태 차원 \(d\)에 \(\mathcal C=[B,AB,\ldots,A^{d-1}B]\), \(\mathcal O=[C^T,(CA)^T,\ldots,(CA^{d-1})^T]^T\)라 합시다. 입력으로 만들 수 있는 상태는 \(\operatorname{col}\mathcal C\), 영입력 관측에서 구별하지 못하는 상태는 \(\ker\mathcal O\)입니다. 유한 Gram은 \(W_c=\mathcal C\mathcal C^T\), \(W_o=\mathcal O^T\mathcal O\)라 각 계수는 해당 지도와 같습니다.
1절의 \(A\), \(B=e_2\), \(C=e_1^T\)에서는
로 둘 다 가역입니다. 반대로 \(B=e_1\)이면 둘째 상태에 충격을 넣을 수 없고, \(C=e_2^T\)이면 첫째 상태를 볼 수 없습니다. 모든 상태가 보이지 않는다고 모든 모수가 식별 불가능하다고 바로 결론내리지 않습니다. 모수 변화가 어느 전달함수를 바꾸는지까지 확인해야 E4의 정보행렬과 연결됩니다.
최소실현은 영 초기상태에서 같은 입력-출력 전달을 갖는 실현 중 상태 차원이 최소라는 뜻입니다. 임의 초기상태를 별도 출력으로 취급하면 제거 가능한 상태의 의미도 달라집니다.
8. Joseph 형태와 제곱근 계산#
짧은 갱신식 \(P^-=KSK^T+P^+\)에서 큰 두 항을 빼면 반올림에 민감합니다. 임의 이득 \(K\)에 대한 실제 갱신오차 공분산은
두 항이 PSD인 Joseph 형태입니다. 최적 \(K\)에서만 이것이 앞의 짧은 식과 같습니다. 임의 \(K\)에 \(P^--KSK^T\)를 쓰면 틀립니다. Joseph 식도 부동소수점에서 엄밀한 PSD 증명서를 자동 주지는 않으며, 인수 형태로 저장하면 더 직접적으로 구조를 보존합니다.
\(P^-=LL^T\), \(R=L_RL_R^T\)이면 \(D=[(I-KH)L,KL_R]\)에 \(P^+=DD^T\)입니다. \(D^T=QR\)로 QR을 하면 새 인수는 \(R^T\)입니다. 예측도 \([AL,L_Q]\)의 전치 QR로 합니다. 여기서 행렬 \(R\)의 기호는 QR 삼각인수일 때만 사용하며 관측잡음 공분산과 문맥을 구별합니다. 분산을 극단적으로 줄인 실험에서 어떤 구현이 최초로 음의 고윳값을 내는지는 입력·정밀도에 달렸으므로 실패 시점을 강제로 만들지 않습니다.
9. 계산으로 확인하기#
import numpy as np
import scipy.linalg as la
A=np.array([[.5,1.],[0.,.5]]); Q=np.eye(2)
P=la.solve_discrete_lyapunov(A,Q)
Pk=la.solve(np.eye(4)-np.kron(A,A),Q.reshape(-1,order='F')).reshape((2,2),order='F')
assert np.allclose(P,[[116/27,8/9],[8/9,4/3]])
assert np.allclose(P,Pk)
assert la.norm(P-A@P@A.T-Q)<1e-12
print('stationary covariance:',P)
stationary covariance: [[4.2962963 0.88888889]
[0.88888889 1.33333333]]
def measurement_update(m,L,y,H,R):
P=L@L.T; S=H@P@H.T+R
K=la.solve(S,H@P,assume_a='pos').T
D=np.column_stack([(np.eye(len(m))-K@H)@L,K@la.cholesky(R,lower=True)])
_,Rt=la.qr(D.T,mode='economic')
return m+K@(y-H@m),Rt.T,K,S
P0=np.array([[2.,1.],[1.,2.]])
m,L,K,S=measurement_update(np.zeros(2),la.cholesky(P0,lower=True),np.array([3.]),np.array([[1.,0.]]),np.ones((1,1)))
assert np.allclose(m,[2.,1.])
assert np.allclose(L@L.T,[[2/3,1/3],[1/3,5/3]])
Cov=np.array([[2.,.5],[.5,2.25]]); yy=np.array([1.,0.])
assert np.isclose(np.linalg.det(Cov),2*17/8)
assert np.isclose(yy@la.solve(Cov,yy),.5+1/34)
print('posterior mean and covariance:',m,L@L.T)
posterior mean and covariance: [2. 1.] [[0.66666667 0.33333333]
[0.33333333 1.66666667]]
10. 연습과 전체 풀이#
문제 1 — 안정성. AR(2)의 두 유한 시차근과 동반 고윳값을 구하라.
풀이. 특성식은 \((\lambda-4/5)(\lambda-2/5)\), 시차 다항식은 \((1-4z/5)(1-2z/5)\)입니다. 고윳값 \(.8,.4\)는 단위원 안이고 시차근 \(1.25,2.5\)는 밖입니다. 따라서 초기 상태 영향이 0으로 갑니다.
문제 2 — 대수와 공분산. \(A=2,Q=1\)의 Stein 해를 해석하라.
풀이. \(P=4P+1\)이므로 \(P=-1/3\)입니다. 선형방정식은 유일해지만 음수이므로 공분산은 아닙니다. 급수 \(\sum4^k\)가 발산하여 정상 분산 해석도 불가능합니다.
문제 3 — 관측. \(H=(0,1)\)로 바꾸면 5절 예제의 결과는?
풀이. 대칭으로 \(S=3\), \(K=(1/3,2/3)^T\), \(m^+=(1,2)^T\)입니다. 공분산은 \(\begin{pmatrix}5/3&1/3\\1/3&2/3\end{pmatrix}\). 관측한 둘째 상태의 분산이 더 많이 줄어듭니다.
문제 4 — Joseph. \(K=0\)이면 어떤 갱신식이 맞는가?
풀이. 평균은 관측을 사용하지 않고 공분산은 \(P^-\). 일반 \(K\)에서는 Joseph 식을 전개해 \(P^--KHP^--P^-H^TK^T+KSK^T\)를 얻습니다. \(KS=P^-H^T\)인 최적 \(K\)에서만 \(P^--KSK^T\)로 줄어듭니다.
문제 5 — 최소실현. 예제의 \(B=e_1,C=e_2^T\)에서 전달함수는?
풀이. 모든 \(k\)에 \(A^ke_1=2^{-k}e_1\)이므로 \(CA^kB=0\). 영 초기상태에서는 입력이 출력에 전혀 나타나지 않습니다. 전달의 최소 상태 차원은 0이며 두 차원 실현은 최소가 아닙니다.
문제 6 — 우도. 6절의 두 관측 우도를 쓰라.
풀이. \(-2\ell=2\log(2\pi)+\log(17/4)+9/17\). 직접 공분산 방식과 혁신 방식 모두 이 값입니다. 정규 상수를 생략하면 같은 표본의 모수 비교에는 무관하지만 실제 밀도값은 달라집니다.
11. 지금까지의 내용을 수학의 언어로 정리해 봅시다#
충격이 사라지는 조건, 충격을 누적한 분산, 관측으로 줄어드는 오차를 각각 거듭제곱·행렬방정식·직교사영의 명제로 옮기자. 앞에서 계산한 두 부문의 모형은 일반 증명의 각 조건이 왜 필요한지 보여 주는 예입니다.
정리 1. 안정성의 네 동치#
유한차원 실수 행렬에 3절의 네 조건이 동치입니다.
증명. \(\rho(A)<1\)이면 보조정리와 정리 3. 노름의 선택과 Gelfand 공식에 의해 임의 \(\rho(A)<q<1\)에 \(\|A^k\|\le Cq^k\)이고 노름급수가 수렴합니다. 급수 수렴이면 항이 0으로 갑니다. \(A^k\to0\)일 때 복소 고유벡터 \(v\ne0\)에 \(A^kv=\lambda^kv\to0\)이므로 모든 \(|\lambda|<1\)입니다. 실수 행렬은 복소화해 같은 결론을 얻습니다.
안정하면 \(P=\sum_{k\ge0}(A^T)^kA^k\)가 절대수렴합니다. 첫 항이 \(I\)라 \(P\succ0\)이고 지표 이동으로 \(A^TPA-P=-I\). 역으로 \(R=P-A^TPA\succ0\)이면 \(R\succeq\delta P\)인 \(\delta>0\)가 있습니다. 이는 \(P^{-1/2}RP^{-1/2}\)의 최소 고윳값입니다. \(A^TPA\succeq0\)이므로 \(\delta\le1\). \(\delta<1\)이면 \(\|A^kx\|_P^2\le(1-\delta)^k\|x\|_P^2\)이고, \(\delta=1\)이면 \(A=0\). 노름 동치로 \(A^k\to0\). ∎
정리 2. Stein 방정식#
모든 행렬 \(Q\)에 \(P=APA^T+Q\)의 유일해가 존재할 필요충분조건은 \(\lambda_i(A)\lambda_j(A)\ne1\)입니다. 안정한 \(A\), \(Q\succeq0\)이면 \(P=\sum_{k\ge0}A^kQ(A^T)^k\succeq0\)입니다.
증명. vec 항등식으로 선형사상은 \(I-A\otimes A\)입니다. 복소 삼각화 \(A=UTU^*\)를 사용하면 \(A\otimes A\)는 \(T\otimes T\)와 닮고 그 대각은 모든 \(\lambda_i\lambda_j\). 따라서 가역성 조건이 정확히 제시한 조건입니다. 안정하면 급수 각 항의 노름은 \(C^2q^{2k}\|Q\|\) 이하라 수렴합니다. 부분합은 PSD이고 PSD 원뿔은 닫혀 있어 극한도 PSD입니다. 지표 이동으로 방정식을 만족합니다. 두 해의 차 \(D\)는 \(D=A^kD(A^T)^k\)이므로 \(k\to\infty\)에서 0입니다. ∎
정리 3. 혁신 사영과 제곱근 갱신#
5절의 무상관 가정과 \(R\succ0\) 아래 제시한 \(K,m^+,P^+\)는 최적 선형예측입니다. 임의 이득의 공분산은 Joseph 식이며, 최적 공분산과의 차이는 \((K-K_*)S(K-K_*)^T\)입니다.
증명. 예측오차 \(e\)의 갱신오차는 \(e-K(He+v)\)입니다. 외적의 기댓값을 전개하고 \(E[ev^T]=0\)을 쓰면 Joseph 식을 얻습니다. \(K_*=P^-H^TS^{-1}\)에 \(K=K_*+D\)를 대입하면 \(K_*S=P^-H^T\) 때문에 교차항이 사라지고 \(P^+(K)=P^+(K_*)+DSD^T\). 따라서 모든 선형결합의 오차분산이 최소입니다. 갱신오차와 혁신의 공분산은 \(P^-H^T-K_*S=0\)입니다. Gaussian이면 무상관 두 블록은 독립이므로 선형예측이 조건부 평균입니다. 인수 \(D_0=[(I-KH)L,KL_R]\)의 곱은 Joseph 식과 정확히 같습니다. \(D_0^T=QR\)에서 \(Q^TQ=I\)이므로 \(D_0D_0^T=R^TR\)입니다. ∎
정리 4. 혁신 우도와 최소실현#
Gaussian 관측, 양정치 혁신 공분산 아래 6절의 우도 분해가 성립합니다. 유한차원 영 초기상태 실현은 제어가능·관측가능할 때 그리고 그때에만 최소입니다.
증명. 과거 관측공간에 대한 사영 잔차는 그 공간에 직교합니다. 이전 혁신은 그 공간에 속하므로 서로 무상관입니다. Gaussian이라 독립이며 각 혁신의 밀도를 곱합니다. 관측에서 혁신으로의 블록 하삼각 변환은 대각 항등이므로 Jacobian은 1입니다. 각 Gaussian 밀도의 로그를 더하면 식을 얻습니다.
최소실현의 경우 Cayley–Hamilton으로 \(A^d\) 이상의 거듭제곱은 앞의 \(d\)개 선형결합입니다. 따라서 \(\operatorname{col}\mathcal C\)는 \(A\)-불변이고 입력이 만드는 모든 상태를 포함합니다. 제어불가능하면 이 공간으로 제한하여 차원을 줄입니다. 관측불가능한 \(\ker\mathcal O\)도 \(A\)-불변이며 \(C\)가 그 공간을 0으로 보내므로 그 몫공간에서 같은 전달을 얻습니다. 따라서 최소이면 두 조건이 필요합니다.
반대로 둘 다 완전계수이면 블록 Hankel \(\mathcal H=(CA^{i+j}B)_{i,j=0}^{d-1}=\mathcal O\mathcal C\)는 계수 \(d\)입니다. \(\mathcal O\)는 단사, \(\mathcal C\)는 전사이므로 그 곱의 계수가 \(d\)라는 결론입니다. 같은 전달을 갖는 \(d'\)차원 실현도 같은 Hankel을 두 지도의 곱으로 표현하므로 \(d'\ge\operatorname{rank}\mathcal H=d\)입니다. 따라서 원래 실현이 최소입니다. ∎
정리 5. 동반행렬과 지연다항식#
\(d\)변수 VAR(\(p\))의 동반행렬 \(\mathcal A\)에 \(\det(\lambda I_{dp}-\mathcal A)=\lambda^{dp}\det(I_d-A_1/\lambda-\cdots-A_p/\lambda^p)\)입니다. 우변은 전개 후 다항식으로 이해합니다.
증명. \(\lambda\ne0\)에서 동반행렬의 아래 블록 방정식은 \(v_j=\lambda^{-(j-1)}v_1\)입니다. 행렬식에는 아래 \((p-1)d\) 블록의 소거로 \(\lambda^{(p-1)d}\)가 생기고, 첫 Schur 보원은 \(\lambda I_d-A_1-A_2/\lambda-\cdots-A_p/\lambda^{p-1}\)입니다. 이 블록에서 \(\lambda^d\)를 묶으면 제시한 식입니다. 양쪽은 전개 후 다항식이며 모든 비영 \(\lambda\)에서 같아 항등적으로 같습니다. 비영 고윳값과 유한 지연근의 역수 대응 및 중복도 일치가 따릅니다. ∎
안정한 전이에서는 충격의 누적 공분산과 필터 갱신을 계산할 수 있습니다. 다음 장에서는 단위근과 제약이 남는 경우를 살펴보고, 공통 추세와 안정한 관계를 분리합니다. E9로 이어 읽기.