471 lines
20 KiB
Python
471 lines
20 KiB
Python
"""
|
||
================================================================================
|
||
국제 유가 충격의 국내 소비자물가 품목별 전이 시차 분석
|
||
- VAR 모형과 충격반응함수(IRF)를 이용한 다변량 시계열 분석
|
||
- 분석 기간: 2015.01 ~ 2024.12 (월별, 120개 관측치)
|
||
================================================================================
|
||
|
||
[분석 구조]
|
||
독립변수 : WTI 국제유가 (로그)
|
||
종속변수 : 에너지 / 식품 / 서비스 품목별 CPI (로그)
|
||
통제변수 : 원/달러 환율 (로그)
|
||
|
||
품목별로 [유가, 환율, 해당 품목 CPI] 3변수 VAR을 각각 추정하고,
|
||
유가 충격에 대한 각 품목의 반응(IRF)을 비교하여 전이 시차를 분석한다.
|
||
"""
|
||
|
||
import io
|
||
import urllib.request
|
||
import warnings
|
||
from pathlib import Path
|
||
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
from matplotlib import rc
|
||
|
||
from statsmodels.tsa.stattools import adfuller, grangercausalitytests
|
||
from statsmodels.tsa.api import VAR
|
||
|
||
warnings.filterwarnings("ignore") # 통계 경고 메시지 숨김 (결과 해석에 불필요)
|
||
|
||
# ------------------------------------------------------------------------------
|
||
# 0. 한글 폰트 설정 (그래프 제목·라벨을 한글로 표시하기 위함)
|
||
# ------------------------------------------------------------------------------
|
||
|
||
output_path = Path("./outputs")
|
||
output_path.mkdir(parents=True, exist_ok=True)
|
||
|
||
rc('font', family='Nanum Gothic')
|
||
plt.rcParams["font.family"] = "Nanum Gothic"
|
||
plt.rcParams["axes.unicode_minus"] = False
|
||
plt.rcParams["figure.dpi"] = 300
|
||
plt.rcParams["savefig.dpi"] = 300
|
||
|
||
# 분석 기간 상수
|
||
START, END = "2015-01", "2024-12"
|
||
|
||
|
||
# ==============================================================================
|
||
# 1. 데이터 수집
|
||
# ==============================================================================
|
||
def fetch_csv(url):
|
||
"""URL에서 CSV를 받아 문자열로 반환하는 헬퍼 함수"""
|
||
req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"})
|
||
return urllib.request.urlopen(req, timeout=30).read().decode()
|
||
|
||
|
||
def load_wti():
|
||
"""
|
||
WTI 국제유가 (일별 → 월평균)
|
||
출처: GitHub datasets/oil-prices (원출처: 미국 EIA)
|
||
"""
|
||
url = "https://raw.githubusercontent.com/datasets/oil-prices/main/data/wti-daily.csv"
|
||
df = pd.read_csv(io.StringIO(fetch_csv(url)))
|
||
df["Date"] = pd.to_datetime(df["Date"])
|
||
# 일별 데이터를 '월초(MS) 기준 월평균'으로 집계
|
||
monthly = df.set_index("Date").resample("MS").mean()
|
||
return monthly.loc[START:END, "Price"].rename("wti")
|
||
|
||
|
||
def load_fx():
|
||
"""
|
||
원/달러 환율 (월별)
|
||
출처: GitHub datasets/exchange-rates (원출처: 미국 연준 계열)
|
||
"""
|
||
url = "https://raw.githubusercontent.com/datasets/exchange-rates/main/data/monthly.csv"
|
||
df = pd.read_csv(io.StringIO(fetch_csv(url)))
|
||
df = df[df["Country"] == "South Korea"].copy()
|
||
df["Date"] = pd.to_datetime(df["Date"])
|
||
fx = df.set_index("Date")["Exchange rate"]
|
||
return fx.loc[START:END].rename("fx")
|
||
|
||
|
||
def load_cpi():
|
||
"""
|
||
품목별 소비자물가지수 (월별)
|
||
출처: 통계청 KOSIS, 지출목적별 소비자물가지수(2020=100)
|
||
로컬에 업로드된 CSV 파일을 읽어서 처리한다.
|
||
|
||
[품목 매핑]
|
||
에너지 <- '04 주택, 수도, 전기 및 연료'
|
||
식품 <- '01 식료품 및 비주류음료'
|
||
서비스 <- '11 음식 및 숙박'
|
||
"""
|
||
path = "./지출목적별_소비자물가지수_품목포함__2020100__20260713213203.csv"
|
||
df = pd.read_csv(path, encoding="utf-8")
|
||
|
||
# 분석에 사용할 지출목적별 항목 지정
|
||
mapping = {
|
||
"04 주택, 수도, 전기 및 연료": "energy",
|
||
"01 식료품 및 비주류음료": "food",
|
||
"11 음식 및 숙박": "service",
|
||
}
|
||
|
||
# YYYY.MM 형식의 월별 컬럼만 추출 (연간 집계 컬럼 '2023' 등은 제외)
|
||
month_cols = [c for c in df.columns if "." in str(c) and len(str(c)) == 7]
|
||
|
||
records = {}
|
||
for kor_name, eng_name in mapping.items():
|
||
row = df[df["지출목적별"] == kor_name][month_cols] # 해당 품목 행 추출
|
||
series = row.iloc[0].astype(float) # 값을 실수로 변환
|
||
series.index = pd.to_datetime(series.index, format="%Y.%m") # 인덱스를 날짜로
|
||
records[eng_name] = series
|
||
|
||
cpi = pd.DataFrame(records)
|
||
# ★ KOSIS CSV는 월 컬럼 순서가 뒤섞여 있으므로 반드시 날짜순 정렬
|
||
cpi = cpi.sort_index()
|
||
return cpi.loc[START:END] # 분석 기간으로 자름
|
||
|
||
|
||
# ==============================================================================
|
||
# 2. 데이터 병합 및 전처리
|
||
# ==============================================================================
|
||
def build_dataset():
|
||
"""세 데이터를 하나의 월별 데이터프레임으로 병합하고 로그 변환한다."""
|
||
wti = load_wti()
|
||
fx = load_fx()
|
||
cpi = load_cpi()
|
||
|
||
# 인덱스(월)를 기준으로 병합
|
||
data = pd.concat([wti, fx, cpi], axis=1)
|
||
data = data.dropna() # 결측 있는 월 제거
|
||
|
||
print("=" * 60)
|
||
print("[1] 데이터 병합 완료")
|
||
print(f" 기간: {data.index[0].date()} ~ {data.index[-1].date()}")
|
||
print(f" 관측치: {len(data)}개월, 변수: {list(data.columns)}")
|
||
print(data.head(3).round(2))
|
||
print()
|
||
|
||
# 로그 변환: 탄력성(% 반응) 해석 + 분산 안정화
|
||
# IRF 해석이 '유가 1% 충격 -> CPI x% 반응' 형태가 됨
|
||
log_data = np.log(data)
|
||
return data, log_data
|
||
|
||
|
||
# ==============================================================================
|
||
# 3. 단위근 검정 (ADF Test)
|
||
# ==============================================================================
|
||
def adf_test(series, name):
|
||
"""
|
||
ADF(Augmented Dickey-Fuller) 검정
|
||
귀무가설 H0: 단위근이 있다 (= 비정상 시계열)
|
||
p < 0.05 이면 H0 기각 -> 정상(stationary) 시계열
|
||
"""
|
||
result = adfuller(series.dropna(), autolag="AIC")
|
||
pval = result[1]
|
||
verdict = "정상 (I(0))" if pval < 0.05 else "비정상 (단위근 존재)"
|
||
print(f" {name:12s} | ADF통계량={result[0]:7.3f} | p값={pval:6.4f} | {verdict}")
|
||
return pval
|
||
|
||
|
||
def run_unit_root_tests(log_data):
|
||
"""수준(level)과 1차 차분에 대해 각각 ADF 검정을 수행한다."""
|
||
print("=" * 60)
|
||
print("[2] 단위근 검정 (ADF Test)")
|
||
print("-" * 60)
|
||
print(" (a) 수준 변수 (log level)")
|
||
for col in log_data.columns:
|
||
adf_test(log_data[col], col)
|
||
|
||
print("\n (b) 1차 차분 변수 (log difference)")
|
||
diff = log_data.diff().dropna()
|
||
for col in diff.columns:
|
||
adf_test(diff[col], f"d.{col}")
|
||
print()
|
||
return diff
|
||
|
||
|
||
# ==============================================================================
|
||
# 4. 품목별 VAR 모형 추정 + IRF + FEVD + Granger 인과성
|
||
# ==============================================================================
|
||
def analyze_item(log_data, item, item_kor, maxlags=12, irf_periods=12):
|
||
"""
|
||
특정 품목에 대해 [유가, 환율, 품목CPI] 3변수 VAR을 추정하고
|
||
유가 충격에 대한 해당 품목 CPI의 반응(IRF)을 분석한다.
|
||
|
||
[모형 설정 — 수준(level) VAR]
|
||
물가·유가·환율은 대부분 I(1) 비정상 시계열이지만, 차분하면 변수 간
|
||
장기 균형관계(공적분) 정보가 사라진다. 특히 서비스 물가는 추세가 강해
|
||
차분해도 정상성이 확보되지 않는다. 따라서 Sims, Stock & Watson(1990)이
|
||
제시한 대로 '로그 수준' 변수로 VAR을 추정한다. 변수 간 공적분이 존재하면
|
||
수준 VAR의 충격반응함수(IRF)는 일치추정량이 되므로, 유가-물가 전이 연구에서
|
||
널리 사용되는 표준적 접근이다. 계절성은 월별 계절더미로 통제한다.
|
||
|
||
Parameters
|
||
----------
|
||
item : 품목 영문명 ('energy'/'food'/'service')
|
||
item_kor : 그래프에 쓸 한글명 ('에너지'/'식품'/'서비스')
|
||
"""
|
||
print("=" * 60)
|
||
print(f"[3] VAR 분석 — {item_kor} 품목")
|
||
print("-" * 60)
|
||
|
||
# (1) 분석 대상 변수 구성: 유가 -> 환율 -> 품목CPI 순서
|
||
# Cholesky 분해 순서 = 외생성이 강한 변수부터 (유가가 가장 외생적)
|
||
cols = ["wti", "fx", item]
|
||
sub = log_data[cols].copy() # 로그 '수준' 변수 (차분하지 않음)
|
||
|
||
# (2) 월별 계절더미 생성 (1~11월, 12월은 기준월로 제외)
|
||
# 식품·에너지 CPI의 계절성을 통제하기 위함
|
||
seasonal = pd.get_dummies(sub.index.month, prefix="m", drop_first=True)
|
||
seasonal.index = sub.index
|
||
seasonal = seasonal.astype(float)
|
||
|
||
# (3) VAR 모형 생성 및 최적 시차 선택 (AIC 기준)
|
||
model = VAR(sub, exog=seasonal) # 계절더미를 외생변수로 포함
|
||
sel = model.select_order(maxlags=maxlags)
|
||
lag = sel.aic # AIC가 최소가 되는 시차
|
||
if lag == 0:
|
||
lag = 1 # 최소 1시차는 확보
|
||
print(f" 최적 시차(AIC 기준): {lag}")
|
||
|
||
results = model.fit(lag) # VAR 추정 (계절더미는 모형에 이미 포함)
|
||
|
||
# (4) 그랜저 인과성 검정: 유가 -> 품목CPI 방향
|
||
print(f" 그랜저 인과성 검정 (유가 → {item_kor}):")
|
||
gc = grangercausalitytests(sub[[item, "wti"]], maxlag=lag, verbose=False)
|
||
pvals = [round(gc[i + 1][0]["ssr_ftest"][1], 4) for i in range(lag)]
|
||
print(f" 시차별 p값: {pvals} (0.05 미만이면 인과성 유의)")
|
||
|
||
# (5) 충격반응함수(IRF) 계산
|
||
irf = results.irf(irf_periods) # irf_periods개월까지 반응 추적
|
||
|
||
# 유가('wti') 충격에 대한 품목CPI(item)의 반응 추출
|
||
# irf.irfs shape = (기간, 반응변수, 충격변수)
|
||
wti_idx = cols.index("wti")
|
||
item_idx = cols.index(item)
|
||
response = irf.irfs[:, item_idx, wti_idx] # 유가충격 -> 품목 반응
|
||
cum_response = irf.cum_effects[:, item_idx, wti_idx] # 누적 반응
|
||
|
||
# (6) 전이 지표 3종 계산
|
||
# [중요] 우리 변수는 로그 '수준'이라 IRF는 누적되어 우상향한다.
|
||
# 단일 지표(peak 시점)는 오해를 부르므로 3가지를 함께 본다.
|
||
# ① onset : 충격 직후 1개월차 전이량 → '얼마나 빨리 반응이 시작되나'
|
||
# ② mean_lag : 전이량 무게중심(가중평균 시점) → '전이의 평균 시차'
|
||
# ③ cumulative : 12개월 누적 전이량 → '전이가 실제로 유의미한가'
|
||
increments = np.diff(response) # 월별 증가분(=수준 IRF 기울기)
|
||
months = np.arange(1, len(increments) + 1) # 1~12개월
|
||
|
||
onset = increments[0] # ① 1개월차 즉각 반응
|
||
pos = np.where(increments > 0, increments, 0) # 양(+) 전이량만 (음수 왜곡 방지)
|
||
mean_lag = (np.sum(months * pos) / np.sum(pos)) if np.sum(pos) > 0 else float("nan") # ②
|
||
cumulative = float(np.sum(increments)) # ③ 12개월 누적 전이량
|
||
|
||
print(f" ▶ onset(1개월 즉각반응) = {onset:+.5f}")
|
||
print(f" ▶ 평균 전이 시차(mean lag) = {mean_lag:.2f}개월")
|
||
print(f" ▶ 12개월 누적 전이량 = {cumulative:+.5f}")
|
||
print()
|
||
|
||
return {
|
||
"item": item,
|
||
"item_kor": item_kor,
|
||
"results": results,
|
||
"irf": irf,
|
||
"response": response,
|
||
"cum_response": cum_response,
|
||
"onset": onset,
|
||
"mean_lag": mean_lag,
|
||
"cumulative": cumulative,
|
||
"lag": lag,
|
||
"cols": cols,
|
||
}
|
||
|
||
|
||
# ==============================================================================
|
||
# 5. 시각화
|
||
# ==============================================================================
|
||
def plot_raw_series(data):
|
||
"""원본 시계열 4개(유가·환율·품목별 CPI) 추이 그래프"""
|
||
fig, axes = plt.subplots(2, 2, figsize=(12, 7))
|
||
# fig.suptitle("원자료 시계열 추이 (2015–2024)", fontsize=15, fontweight="bold")
|
||
|
||
plots = [
|
||
("wti", "WTI 국제유가 (달러/배럴)", axes[0, 0]),
|
||
("fx", "원/달러 환율", axes[0, 1]),
|
||
("energy", "에너지 CPI", axes[1, 0]),
|
||
("food", "식품 CPI", axes[1, 1]),
|
||
]
|
||
for col, title, ax in plots:
|
||
ax.plot(data.index, data[col], color="#333333", linewidth=1.2)
|
||
ax.set_title(title, fontsize=12)
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.tight_layout()
|
||
plt.savefig("./outputs/01_원자료_시계열.png", bbox_inches="tight")
|
||
plt.close()
|
||
print(" 저장: 01_원자료_시계열.png")
|
||
|
||
|
||
def plot_irf_with_ci(analyses):
|
||
"""
|
||
품목별 유가충격 IRF를 개별 패널로 그리고 95% 신뢰구간을 함께 표시한다.
|
||
신뢰구간이 0을 포함하지 않는 구간 = 통계적으로 유의한 반응.
|
||
"""
|
||
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5), sharey=True)
|
||
# fig.suptitle("유가 충격에 대한 품목별 소비자물가 반응 (95% 신뢰구간 포함)",
|
||
# fontsize=14, fontweight="bold")
|
||
|
||
colors = {"energy": "#c0392b", "food": "#e67e22", "service": "#16a085"}
|
||
for ax, a in zip(axes, analyses):
|
||
irf = a["irf"]
|
||
wti_idx = a["cols"].index("wti")
|
||
item_idx = a["cols"].index(a["item"])
|
||
|
||
resp = a["response"]
|
||
# 신뢰구간(표준오차 기반, 약 95% = ±1.96 SE)
|
||
stderr = irf.stderr()[:, item_idx, wti_idx]
|
||
months = np.arange(len(resp))
|
||
lower = resp - 1.96 * stderr
|
||
upper = resp + 1.96 * stderr
|
||
|
||
ax.plot(months, resp, color=colors[a["item"]], linewidth=1.8, marker="o", markersize=3)
|
||
ax.fill_between(months, lower, upper, color=colors[a["item"]], alpha=0.15) # 신뢰구간 음영
|
||
ax.axhline(0, color="gray", linewidth=0.8, linestyle="--")
|
||
ax.axvline(a["mean_lag"], color="gray", linewidth=0.8, linestyle=":") # 평균 전이 시차 표시
|
||
ax.set_title(f"{a['item_kor']} (평균 전이 시차: {a['mean_lag']:.1f}개월)", fontsize=12)
|
||
ax.set_xlabel("경과 개월")
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
axes[0].set_ylabel("소비자물가 반응")
|
||
plt.tight_layout()
|
||
plt.savefig("./outputs/04_품목별_IRF_신뢰구간.png", bbox_inches="tight")
|
||
plt.close()
|
||
print(" 저장: 04_품목별_IRF_신뢰구간.png")
|
||
|
||
|
||
def plot_irf_comparison(analyses):
|
||
"""품목별 유가충격 IRF를 한 그래프에 겹쳐서 전이 시차를 비교"""
|
||
fig, ax = plt.subplots(figsize=(11, 6))
|
||
|
||
colors = {"energy": "#c0392b", "food": "#e67e22", "service": "#16a085"}
|
||
for a in analyses:
|
||
months = range(len(a["response"]))
|
||
ax.plot(months, a["response"],
|
||
marker="o", markersize=4, linewidth=1.5,
|
||
color=colors[a["item"]],
|
||
label=f"{a['item_kor']} (누적 {a['cumulative']:+.3f})")
|
||
|
||
ax.axhline(0, color="gray", linewidth=0.8, linestyle="--") # 0 기준선
|
||
# ax.set_title("유가 충격에 대한 품목별 소비자물가 반응 (충격반응함수)",
|
||
# fontsize=14, fontweight="bold")
|
||
ax.set_xlabel("유가 충격 이후 경과 개월")
|
||
ax.set_ylabel("소비자물가 반응 (로그 차분)")
|
||
ax.legend(title="품목", fontsize=11)
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.tight_layout()
|
||
plt.savefig("./outputs/02_품목별_IRF비교.png", bbox_inches="tight")
|
||
plt.close()
|
||
print(" 저장: 02_품목별_IRF비교.png")
|
||
|
||
|
||
def plot_monthly_increments(analyses):
|
||
"""
|
||
품목별 '월별 증가분'(수준 IRF의 기울기 = 그 달의 순수 전이량)을 그린다.
|
||
수준 IRF는 우상향해서 전이 시차가 안 보이지만, 증가분으로 보면
|
||
'언제 전이가 가장 활발한가'(=전이 시차)가 봉우리로 명확히 드러난다.
|
||
"""
|
||
fig, ax = plt.subplots(figsize=(11, 6))
|
||
|
||
colors = {"energy": "#c0392b", "food": "#e67e22", "service": "#16a085"}
|
||
for a in analyses:
|
||
inc = np.diff(a["response"]) # 월별 증가분
|
||
months = range(1, len(inc) + 1) # 1개월차부터
|
||
ax.plot(months, inc,
|
||
marker="o", markersize=4, linewidth=1.5,
|
||
color=colors[a["item"]],
|
||
label=f"{a['item_kor']} (평균 전이 시차: {a['mean_lag']:.1f}개월)")
|
||
|
||
ax.axhline(0, color="gray", linewidth=0.8, linestyle="--")
|
||
# ax.set_title("유가 충격의 월별 전이량 (수준 IRF의 기울기)",
|
||
# fontsize=14, fontweight="bold")
|
||
ax.set_xlabel("유가 충격 이후 경과 개월")
|
||
ax.set_ylabel("월별 전이량 (그 달의 물가 상승 기여분)")
|
||
ax.legend(title="품목", fontsize=11)
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.tight_layout()
|
||
plt.savefig("./outputs/05_월별_전이량.png", bbox_inches="tight")
|
||
plt.close()
|
||
print(" 저장: 05_월별_전이량.png")
|
||
|
||
|
||
def plot_cumulative_irf(analyses):
|
||
"""품목별 누적 반응(cumulative IRF) 비교"""
|
||
fig, ax = plt.subplots(figsize=(11, 6))
|
||
|
||
colors = {"energy": "#c0392b", "food": "#e67e22", "service": "#16a085"}
|
||
for a in analyses:
|
||
months = range(len(a["cum_response"]))
|
||
ax.plot(months, a["cum_response"],
|
||
marker="s", markersize=4, linewidth=1.5,
|
||
color=colors[a["item"]], label=a["item_kor"])
|
||
|
||
ax.axhline(0, color="gray", linewidth=0.8, linestyle="--")
|
||
# ax.set_title("유가 충격에 대한 품목별 소비자물가 누적 반응",
|
||
# fontsize=14, fontweight="bold")
|
||
ax.set_xlabel("유가 충격 이후 경과 개월")
|
||
ax.set_ylabel("누적 반응 (로그 차분 누적)")
|
||
ax.legend(title="품목", fontsize=11)
|
||
ax.grid(True, alpha=0.3)
|
||
|
||
plt.tight_layout()
|
||
plt.savefig("./outputs/03_품목별_누적반응.png", bbox_inches="tight")
|
||
plt.close()
|
||
print(" 저장: 03_품목별_누적반응.png")
|
||
|
||
|
||
# ==============================================================================
|
||
# 6. 메인 실행 흐름
|
||
# ==============================================================================
|
||
def main():
|
||
# --- 1) 데이터 준비 ---
|
||
data, log_data = build_dataset()
|
||
|
||
# 병합 데이터를 CSV로 저장 (보고서 부록/재현용)
|
||
data.to_csv("./outputs/00_병합데이터.csv", encoding="utf-8-sig")
|
||
|
||
# --- 2) 단위근 검정 ---
|
||
run_unit_root_tests(log_data)
|
||
|
||
# --- 3) 품목별 VAR 분석 ---
|
||
items = [
|
||
("energy", "에너지"),
|
||
("food", "식품"),
|
||
("service", "서비스"),
|
||
]
|
||
analyses = [analyze_item(log_data, eng, kor) for eng, kor in items]
|
||
|
||
# --- 4) 시각화 ---
|
||
print("=" * 60)
|
||
print("[4] 그래프 생성")
|
||
print("-" * 60)
|
||
plot_raw_series(data)
|
||
plot_irf_comparison(analyses)
|
||
plot_irf_with_ci(analyses)
|
||
plot_monthly_increments(analyses)
|
||
plot_cumulative_irf(analyses)
|
||
print()
|
||
|
||
# --- 5) 최종 요약 ---
|
||
print("=" * 70)
|
||
print("[5] 전이 지표 종합 (세 지표를 함께 해석)")
|
||
print("-" * 70)
|
||
print(f" {'품목':8s} | {'onset(즉각반응)':14s} | {'평균전이시차':12s} | {'12개월 누적':12s} | VAR시차")
|
||
print(" " + "-" * 64)
|
||
for a in analyses:
|
||
print(f" {a['item_kor']:8s} | {a['onset']:+.5f} | "
|
||
f"{a['mean_lag']:>6.2f}개월 | {a['cumulative']:+.4f} | lag={a['lag']}")
|
||
print()
|
||
print(" [해석 지침] 단일 지표는 오해를 부르므로 세 지표를 함께 본다:")
|
||
print(" · onset → 반응이 얼마나 빨리 시작되나 (전이 속도)")
|
||
print(" · mean_lag → 전이의 무게중심 (평균 시차)")
|
||
print(" · 누적 → 전이가 실제로 유의미한가 (양수여야 실질 전이)")
|
||
print(" ※ 식품은 mean_lag가 짧아도 누적≈0 → 전이 없음(반짝 후 반전)")
|
||
print("=" * 70)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main() |