20. R code for 06.~08.

저번 포스팅에서는 AR, MA, ARMA 모델을 R code로 구현하는 방법에 대해 알아보았다.

 

이번 시간에는 주어진 시계열 데이터가 정상성을 갖는지에 대한 검정과,

 

ARIMA 모델을 R code로 구현하는 방법에 대해 알아본다.

 

---

 

우선 다음과 같이 NVIDIA의 월간 주가 데이터를 가져온다.

 

> data <- read.csv('data19.csv', header = T)
> data <- data[,1:2]
> head(data)
        date     close
1 2014-08-31 0.4649273
2 2014-09-30 0.4410236
3 2014-10-31 0.4670786
4 2014-11-30 0.5033823
5 2014-12-31 0.4812979
6 2015-01-31 0.4608935

2014년 8월부터 2022년 12월까지의 데이터를 시각화한다.

> nvda <- ts(data$close, start = c(2014,8), freq = 12)
> nvda1 <- window(nvda, end = c(2022,12))
> plot(nvda1, xlab='', ylab='')

> pmax <- floor(12*(length(nvda1)/100)^(1/4))
> fit <- ur.df(nvda1, type = 'trend', lags = (pmax - 1))
> summary(fit)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

(...)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -0.053998   0.414129  -0.130   0.8967  
z.lag.1       0.142413   0.087490   1.628   0.1084  
tt           -0.004952   0.013973  -0.354   0.7242  
z.diff.lag1  -0.379520   0.179659  -2.112   0.0385 *
z.diff.lag2  -0.148561   0.210095  -0.707   0.4820  
z.diff.lag3  -0.035753   0.213071  -0.168   0.8673  
z.diff.lag4  -0.153174   0.231090  -0.663   0.5098  
z.diff.lag5   0.287419   0.233201   1.232   0.2222  
z.diff.lag6  -0.159862   0.237911  -0.672   0.5040  
z.diff.lag7   0.161472   0.264729   0.610   0.5440  
z.diff.lag8  -0.508627   0.263694  -1.929   0.0581 .
z.diff.lag9   0.220560   0.267851   0.823   0.4133  
z.diff.lag10 -0.167435   0.269119  -0.622   0.5360  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

위의 결과에서 t-ratio의 절댓값이 1.6보다 큰 lag8을 pmax로 설정하고 Dickey-Fuller 검정을 실시한다.

 

> adfTest(nvda1, lags = 8, type = 'ct')

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 8
  STATISTIC:
    Dickey-Fuller: 1.8448
  P VALUE:
    0.99

5%의 유의수준 하에서 단위근 검정의 귀무가설을 기각하지 못하기 때문에 $x_t$가 ARIMA(p,1,q) 모델로 표현될 수 있음을 의미한다.

 

> aic <- function(p,q){
+   fit <- Arima(nvda1, order = c(p,1,q), include.constant = T)
+   fit$aic
+ }
> aic(0,0)
[1] 290.696
> aic(0,1)
[1] 292.6605
> aic(0,2)
[1] 294.6407
> aic(0,3)
[1] 291.5097
> aic(1,0)
[1] 292.6615
> aic(1,1)
[1] 294.66
> aic(1,2)
[1] 294.4015
> aic(1,3)
[1] 289.5037
> aic(2,0)
[1] 294.6182
> aic(2,1)
[1] 294.2128
> aic(2,2)
[1] 292.4198
> aic(2,3)
[1] 289.5294
> aic(3,0)
[1] 291.4181
> aic(3,1)
[1] 289.7972
> aic(3,2)
[1] 289.9142
> aic(3,3)
[1] 291.6297

AIC(1,3)의 값이 제일 작으므로 ARIMA(1,1,3) 모델을 생성한다.

 

> fit <- Arima(nvda1, order = c(1,1,3), include.constant = T)
> fit
Series: nvda1 
ARIMA(1,1,3) with drift 

Coefficients:
         ar1      ma1      ma2     ma3   drift
      0.7337  -0.8160  -0.0593  0.3567  0.3855
s.e.  0.1502   0.1805   0.1794  0.1419  0.2193

sigma^2 = 1.442:  log likelihood = -138.75
AIC=289.5   AICc=290.54   BIC=304.37

결정된 모델을 식으로 나타내면 다음과 같다.

 

$$ (\Delta x_t - 0.3855) = 0.7337(\Delta x_{t-1} - 0.3855) + \varepsilon_t - 0.8160 \varepsilon_{t-1} - 0.0593 \varepsilon_{t-2} + 0.3567 \varepsilon_{t-3} .$$

 

결정된 ARIMA 모델로 예측을 실시하고 실젯값과 비교해본다.

 

> fcst <- forecast(fit, 6)$mean
> real <- window(nvda, start = c(2022,1))
> plot(fcst, xlab = '', ylab = '', ylim = c(15,40))
> lines(real, lty=2)
> legend('topleft', c('Observations', 'Forecasts'), lty = c(1,2), inset = 0.01)

이전 AR, MA, ARMA 모델과 마찬가지로 충분한 설명력을 가지고 있지 않다.

 

---

 

오늘은 정상성을 가지지 않는 데이터를 다루는 모델인 ARIMA 모델을 R code로 구현하는 방법에 대해 알아보았다.

 

저번 포스팅에서도 말했지만, 주가는 기본적으로 random walk를 가정하기 때문에 가장 기본적인 선형모델인 ARIMA 모델로 예측을 실시하기에는 무리가 있다.

 

하지만 R code로 어떻게 구현하는지 흐름을 보기에는 충분하다고 생각한다.

 

오늘은 여기까지.