VARMAX와 상태공간 모형으로 다변량 시계열을 다루는 법
VARMAX와 상태공간 모형의 구조, 추정 방식, 진단 항목과 결측·외생 변수·구조 변화에 대응하는 시계열 모델링 방법을 정리합니다.
2026-08-14 · 최초 발행 2024-04-29
다변량 시계열에서 VARMAX와 상태공간 모형이 맡는 역할
여러 시계열이 함께 움직일 때는 각 변수의 과거만 보는 방식으로는 관계를 충분히 설명하기 어렵다. VARMAX는 변수 간 상호작용과 외생 요인을 한 모형 안에 넣고, 상태공간 모형은 관측값 뒤의 추세·계절·레짐 같은 상태를 분리해 다룬다.
VARMAX(Vector AutoRegressive Moving Average with eXogenous)는 다변량 자기회귀-이동평균 모형에 외생 변수를 포함한다. 식은 Y*t = c + Σ A_i Y*{t−i} + Σ M*j ε*{t−j} + B X_t + ε_t 형태이며, 변수 사이의 상호 인과와 공변동, 외부 요인의 영향을 동시에 반영한다. VARMA와 달리 마케팅이나 정책 같은 외생 충격을 구조 안에 포함할 수 있다.
상태공간 모형은 관측식 y*t = Z_t α_t + ε_t, 상태식 α_t = T_t α*{t−1} + R_t η_t로 표현한다. ε_t, η_t는 보통 가우시안 백색잡음이며, 칼만 필터와 스무더를 통해 추정과 예측을 수행한다. 불완전 관측, 불규칙 간격, 시간가변 계수, 계층 구조를 다루기 쉬운 프레임이다.
관계 해석과 상태 분해의 차이
VARMAX에서는 차수 p, q와 안정성이 핵심이다. 고유값이 1보다 작은지 확인하고, ACF/PACF, AIC/BIC, 교차상관, 도메인 지식을 함께 사용해 구조를 정한다. 비정상 시계열에 공적분이 있으면 VECM으로 재표현할 수 있으며, 외생 변수는 동시 효과와 지연 효과를 구분해 설계한다.
두 모형 모두 최대우도(MLE)를 기반으로 추정한다. 상태공간 모형의 칼만 필터 복잡도는 O(T · n_state^3)이고, VARMAX는 모수 제약과 공분산 구조에 따라 수치최적화를 사용한다. 결측치나 누락 구간은 상태공간 모형에서 자연스럽게 처리할 수 있으며, EM 알고리즘은 초기값과 은닉상태 추정을 보강하는 데 쓰인다.
운영 시에는 잔차의 독립성·정규성·등분산성을 Ljung-Box, Q-Q, ARCH 테스트로 확인해야 한다. 공분산 행렬의 양정성도 보장해야 한다. 구조변화(Break) 검정과 롤링 윈도 성능 모니터링을 병행하고, 축차적 차수 제어와 벌점화(regularization)로 과적합을 억제한다.
VARMAX는 충격반응함수(IRF), 분산분해(FEVD), 그랜저 인과 검정을 통해 정책이나 프로모션 시나리오를 해석하는 데 유리하다. 상태공간 모형은 스무딩 결과에서 추세·계절·레짐을 나누어 읽을 수 있고, 시간가변 계수로 구조 변화도 추적할 수 있다.
함께 움직이는 지표를 모델에 넣는 장면
리테일과 제조에서는 SKU 수요, 가격, 프로모션, 매장 이벤트를 VARMAX로 함께 모델링하고, 시나리오별 매출과 마진 민감도를 분석할 수 있다. 이 결과는 리드타임과 재고 정책을 다루는 공급망 의사결정으로 이어진다.
성장률, 물가, 금리, 환율처럼 거시지표가 얽힌 경우에는 상호작용을 포착하고 IRF로 정책 충격의 파급 경로를 평가할 수 있다. 공분산 예측은 동적 자산배분을 위한 위험예측에도 활용된다.
멀티센서 환경에서는 상태공간 모형이 신호를 융합하고 결측을 보간하며 잡음을 분리한다. 로컬레벨과 시간가변 회귀를 결합하면 이상치에 대한 견고성을 확보할 수 있고, 실시간 칼만 필터링은 스트리밍 추정에 맞는다.
진단부터 운영까지 이어지는 흐름
먼저 ADF/PP와 Johansen 테스트로 단위근과 공적분을 확인하고, 스케일 정규화와 이상치·결측 처리를 수행한다. 그랜저 인과와 교차상관 분석은 후보 지연 차수와 외생 변수 선정에 사용한다.
VARMAX에서는 (p, q)와 외생 변수의 지연 라그를 설계한 뒤 안정성 제약 아래에서 최적화한다. 상태공간 모형에서는 Z, T, R, Q, H를 정의하고 칼만 필터·스무더와 MLE 또는 EM을 적용한다.
적합 뒤에는 잔차 검정, 고유값 안정성, 파라미터 유의성을 확인한다. VARMAX라면 IRF의 안정성과 표준오차도 검토 대상이다. 필요하면 모형을 단순화하거나 시간가변·레짐전환·계층화 구조로 확장한다.
운영 단계에서는 롤링 재학습 주기를 설정한다. 예를 들어 주별로 재학습할 수 있으며, PSI와 KS로 데이터 드리프트를 감지한다. 예측 성능 SLA는 sMAPE, RMSE, MAPE, CRPS(확률예측)로 추적한다.
| 항목 | VARMAX | 상태공간(State-Space) |
|---|---|---|
| 성능(정확도) | 상호작용·외생 효과 명시 시 강점. 정상성 가정 하 안정적 | 결측·비정상·시간가변 구조에서 강점. 스무딩으로 노이즈 견고 |
| 확장성 | 차수·변수 증가 시 모수 폭증, 규제화 필요 | 상태 차원 확장 용이하나 n_state^3 비용. 희소성/분해로 완화 |
| 일관성(해석) | IRF/FEVD/그랜저로 경제적 해석 용이 | 상태 분해(추세·계절·레짐) 해석 우수, 직접적 인과 해석은 설계 의존 |
| 안정성 | 고유값<1, 역이동평균 조건 민감 | 칼만 필터 수치안정성 중요, Q/H 양정성 요구 |
| 운영 편의 | statsmodels VARMAX API 간결 | 결측/온라인 추정, 구조적 변경 대응 유리 |
Python으로 확인하는 VARMAX와 상태공간 모형
환경과 전제조건은 다음과 같다.
- Python 3.10+, statsmodels 0.14+, numpy, pandas, matplotlib
- 예시 데이터: 합성 시뮬레이션
외생 변수를 포함한 VARMAX
# pip install statsmodels pandas numpy matplotlib
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.varmax import VARMAX
np.random.seed(42)
n = 400
# 외생 변수 X: 계절 + 잡음
t = np.arange(n)
X = 0.5*np.sin(2*np.pi*t/24) + 0.1*np.random.randn(n)
# 2변량 VAR(1) + MA(1) 구조
A1 = np.array([[0.5, 0.2],
[0.1, 0.4]])
M1 = np.array([[0.3, 0.0],
[0.0, 0.2]])
B = np.array([[0.8],
[0.3]])
eps = np.random.multivariate_normal(mean=[0,0], cov=np.array([[1.0,0.3],[0.3,0.8]]), size=n)
y = np.zeros((n,2))
u = np.zeros((n,2))
for t_idx in range(1, n):
u[t_idx] = eps[t_idx] + (M1 @ eps[t_idx-1])
y[t_idx] = (A1 @ y[t_idx-1]) + (B @ np.array([X[t_idx]])).ravel() + u[t_idx]
df = pd.DataFrame(y, columns=['y1','y2'])
df['x'] = X
# 학습/검증 분할
train, test = df.iloc[:320], df.iloc[320:]
# 모델 적합: VARMAX(p=1, q=1) + exog
model = VARMAX(endog=train[['y1','y2']], exog=train[['x']], order=(1,1), trend='c')
res = model.fit(maxiter=500, disp=False)
print(res.summary())
# 다단계 예측
pred = res.get_forecast(steps=len(test), exog=test[['x']])
f_mean = pred.predicted_mean
f_ci = pred.conf_int()
# 시각화
ax = df[['y1','y2']].plot(figsize=(10,4))
f_mean.index = test.index
f_mean[['y1','y2']].plot(ax=ax, style='--')
plt.show()초기값이나 공분산 제약에서 실패하면 enforce_stationarity=False, enforce_invertibility=False를 검토할 수 있다. 외생 변수를 표준화하면 수치 안정성 확보에 도움이 된다. 차수는 p,q를 0~2 범위 격자로 두고 AIC/BIC 최소화를 기준으로 탐색하는 방식을 권장한다.
로컬레벨과 결측 처리를 포함한 상태공간 모형
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.structural import UnobservedComponents
np.random.seed(7)
n = 300
x = np.sin(np.arange(n)/12) + 0.2*np.random.randn(n)
level_true = np.cumsum(0.05 + 0.02*np.random.randn(n)) # 랜덤워크 추세
y = level_true + 0.8*x + 0.5*np.random.randn(n)
# 결측 주입
y[50:60] = np.nan
y[150] = np.nan
df = pd.DataFrame({'y': y, 'x': x})
# 로컬레벨(State-Space) + 회귀(exog)
mod = UnobservedComponents(endog=df['y'], level='local level', exog=df[['x']])
res = mod.fit(disp=False)
print(res.summary())
# 칼만 스무딩으로 추세 복원 및 결측 자동 보간
smoothed = res.smoothed_state[0] # level state
pred = res.get_forecast(steps=24, exog=pd.DataFrame({'x': np.sin(np.arange(n,n+24)/12)}))
mean = pred.predicted_mean
ci = pred.conf_int()
ax = df['y'].plot(figsize=(10,4), label='obs', alpha=0.6)
pd.Series(smoothed, index=df.index, name='level').plot(ax=ax, label='smoothed level')
mean.plot(ax=ax, style='--', label='forecast')
plt.legend()
plt.show()상태공간 구성 요소는 level, local linear trend, seasonal, ar(1) 등을 조합할 수 있다. Q(상태잡음)와 H(관측잡음)의 추정 안정성을 점검해야 하며, 실시간 필터링이 필요하면 res.filter_results를 활용해 순차 업데이트를 구현한다.
모형 복잡도와 운영 비용의 균형
변수 선택에는 그랜저 인과와 도메인 피처 엔지니어링을 함께 적용한다. 휴일이나 프로모션 라그도 이 단계에서 검토한다. 안정성은 고유값 제약과 공분산 양정성 유지로 관리하며, Cholesky 파라메터화를 사용할 수 있다. 검증에서는 롤링 원스텝·다단계 예측 평가와 예측분포 적합도(PIT, CRPS) 점검이 필요하다.
높은 차수와 시간가변 구조는 정확도를 높일 수 있지만 해석성과 운영 복잡도를 함께 키운다. 다변량 모델에서는 모수가 빠르게 늘어나므로 벌점화(L1/L2), 차원 축소(PC-VAR), 공분산 구조 가정으로 이를 제어한다. 상태공간 모형은 n_state^3 비용을 가지며, VARMAX는 변수와 차수가 늘어날수록 최적화 비용이 급증한다. 배치와 온라인 운용을 혼합하는 설계가 필요한 이유다.
동일 변수 단변량 대비 예측 오차(RMSE, sMAPE)는 5~20% 개선 가능하며, 이는 데이터와 도메인에 의존한다. 결측·잡음 환경에서는 상태공간 스무딩으로 신호대잡음비(SNR)를 개선하고 다운스트림 오탐률을 10%p 내외 낮춘 사례가 관찰됐다.
IRF와 FEVD 기반의 시나리오 분석은 정책·마케팅 의사결정의 신뢰도를 높일 수 있다. 시간가변 계수와 상태 분해는 구조 변화를 더 일찍 감지하고 리스크 관리를 강화하는 데도 활용된다. 관계 해석과 정상성이 중심이면 VARMAX를, 결측과 구조 변화가 핵심이면 상태공간 모형을 우선 평가하되, 안정성 진단·롤링 검증·모델 단순화를 운영 원칙으로 둔다.