오늘도 데이터: R 통계 고급 시계열 문제 벡터 자기회귀

네이버 지식인에서 아래와 같은 질문이 왔다.

미국 연방기금금리 1%p 인상이 미국의 실업률과 인플레이션율에 미치는 동태적 영향을 VAR 모형과 직교화충격반응함수를 이용해 분석하시오.라는 과제를 받았습니다.. rstudio의 VAR함수와 irf 함수를 사용해야하는데 데이터는 만들어놨습니다... 근데 이데이터를 VAR함수에 적용하고 irf함수에 적용하는법을 감을 못잡겠습니다.. 도와주실분 계신가요...

아래 데이터 셋에서는 dcpi 가 인플레이션 이고 ur이 실업률이다. 데이터 셋을 보니 놀랍게도 2021년 4월 까지 있었다.

var 모형

아래의 글을 작성 하는데 아래의 링크를 활용 하였다.

https://otexts.com/fppkr/VAR.html#fn25

VAR 모형

위의 글을 인용 하여 질문의 내용을 요약 하자면 아래와 같다.

VAR 모델은 이론적이지 않다고 비판 받기도 한다.

즉, 이론적인 구조를 식에 반영하는 몇몇 경제학적 이론에 기초하여 세운 것이 아니다.

모든 변수는 시스템 안의 다른 모든 변수에 영향을 준다는 가정을 하기 때문에, 추정한 계수를 해석 하는 것이 어렵다.

하지만 VAR 모델은 여러 상황에 유용함.

분명한 해석이 필요 없을 때 관련된 변수의 모임을 예측 할때
그랑거(Granger) 인과율(causality) 검정에 기초하여) 한 변수가 다른 것을 예측이 유용한지를 판단 할때
충격반응분석(impulse response analysis), 한 변가가 다른 한 변수의 갑작 스럽지만 일시적인 변화에 반응 하는 것을 분석 할때
예측 분산 분해를 예측 할때, 다른 변수의 효과의 결과가 각 변수의 예측 분산의 비율이 될때

이번 과제는 코로나 상황에서 변화된 충격 반응 분석과 관련 된 내용 같다.

library(rio)
library(forecast)

library(vars)
library(dplyr)
library(readxl)

library(tsDyn)
library(httr)
library(dplyr)

# 파일을 불러오는 path를 url 변수로 정의      
url = "https://docs.google.com/uc?export=download&id=1aagKp808dYlXxceRAOBT2mw73Iq0etel"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))
gdp <- read_excel(tf, 1L)  # 1L은 sheet 

# 시계열 데이터 셋을 만든다.     
  gdp <- ts(gdp %>% 
              select(dcpi,ur),
              start = c(1954, 07),
              end = c(2021,04),
              frequency = 12)

  head(gdp)

##                 dcpi  ur
## Jul 1954 -0.29695620 5.8
## Aug 1954 -0.03723008 6.0
## Sep 1954 -0.14897579 6.1
## Oct 1954 -0.33569564 5.7
## Nov 1954  0.22455090 5.3
## Dec 1954 -0.03734130 5.0

이번 질문자 님께서는 아마도 위와 같이 데이터셋 만 만들 줄 아셨다면, 아주 쉽게 접근 하였을 것 같다.

아래는 인플레이션, 실업률의 벡터 자기회귀모형 구현 한 것이다. 물론 여기에서 카이제곱이 작고, 자유도가 28정도 되는 것이 가장 좋은 것으로 생각되어, 아래와 같이 그래프를 그렸다.

# dcpi가 인플레이션율, ur이 실업률   

  
# 벡터 자기회귀 
# var  AIC , HQ, SE,FPE 시차수를 고른다.
  VARselect(gdp, lag.max=8,
            type="const")[["selection"]]

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      7      4      1      7

# Portmanteau Test (asymptotic)
  var1 <- VAR(gdp, p=1, type="const")
  serial.test(var1, lags.pt=10, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 95.448, df = 36, p-value = 2.788e-07

  var2 <- VAR(gdp, p=2, type="const")
  serial.test(var2, lags.pt=10, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var2
## Chi-squared = 64.787, df = 32, p-value = 0.0005301

  var3 <- VAR(gdp, p=3, type="const")
  serial.test(var3, lags.pt=10, type="PT.asymptotic")

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var3
## Chi-squared = 53.237, df = 28, p-value = 0.002754

# 잔차에(residual) 계열상관에 대한 검증을 통과함
  forecast(var3) %>% 
    autoplot()

벡터 자기회귀

직교화충격반응함수

아래와 같이 인플레이션이 실업률에 충격을 주는지 보는 직교 좌표 이다.

 For VAR에 대한 직교충격 반응 함수  
  mod_var <- lineVar(gdp, lag = 2)
  
 irf(mod_var, impulse = "dcpi", response = c("dcpi", "ur"), boot =
       FALSE)

## 
## Impulse response coefficients
## $dcpi
##              dcpi           ur
##  [1,] 0.240343656 -0.041724820
##  [2,] 0.126424607 -0.054933135
##  [3,] 0.102308775 -0.036797312
##  [4,] 0.072963843 -0.027316829
##  [5,] 0.053808235 -0.019539348
##  [6,] 0.039321533 -0.014020958
##  [7,] 0.028807391 -0.009959866
##  [8,] 0.021091419 -0.007004504
##  [9,] 0.015445930 -0.004850446
## [10,] 0.011311922 -0.003284171
## [11,] 0.008285359 -0.002147514

# For VECM에 대한 직교충격 반응 함수
  mod_VECM <- VECM(gdp, lag = 2, estim="ML", r=1)
  irf(mod_VECM, impulse = "dcpi", response = c("dcpi", "ur"), boot =
        FALSE)

## 
## Impulse response coefficients
## $dcpi
##             dcpi           ur
##  [1,] 0.23873905 -0.042981595
##  [2,] 0.12171448 -0.059666699
##  [3,] 0.08469422 -0.057022400
##  [4,] 0.07823935 -0.028034952
##  [5,] 0.06064821 -0.014235989
##  [6,] 0.04826459 -0.006201328
##  [7,] 0.03899538  0.001868021
##  [8,] 0.03111186  0.008286613
##  [9,] 0.02494715  0.013301932
## [10,] 0.02004418  0.017341856
## [11,] 0.01610318  0.020558603