17. 시계열 분석-2

안녕하세요~ 오늘은 시계열 분석의 두 번째 파트, 실습을 조금 더 해보려고 합니다.

바로 시작해보도록 하겠습니다.

저 또한 여기서 쓰이는 패키지는 어떻게 구성되어있는지 몰라서 예제에 쓰이는 데이터가 무엇인지 직접찾아봅니다.

함수가 궁금할 경우 앞에 ? 를, 데이터가 궁금할 경우 앞에 ??를 치면 R이 검색해줍니다.

install.packages("fpp")

library(fpp)

data(elecequip)

# elecequip은 Manufacture of electrical equipment: computer, electronic and optical products. Data adjusted by working days; Euro area (16 countries). Industry new orders index. 2005=100. 라고 나와있네요. 유럽 16개국의 컴퓨터와 같은 전기 혹은 광학제품 제조량에 대한 데이터이며, 2005년도를 100으로 index하여 비교한 데이터 같군요?...

#stl 함수는 seasonal, trend and irregular components를 분해하는 함수이네요.

#seasadj함수는 Returns seasonally adjusted data constructed by removing the seasonal component. 라고 설명되어있습니다.

#즉 계절적 요인을 조정한 데이터를 도출해주는 함수인거 같습니다.

#따라서 stl함수로 seasonal component를 분해하고 seasadj함수로 조정하는 것인 것 같습니다.

#그래서 elecequip데이터가 시계열 데이터이므로 각종 요인들을 제거한 데이터를 eeadj라는 변수에 저장하고 plot을 찍어보겠습니다.

eeadj <- seasadj(stl(elecequip, s.window="periodic"))
plot(eeadj)

#이렇게 플랏을 찍어봤는데요~ 꾸준히 상승과 하락을 반복하다가 2000년대 들어서 급증하게되고, 급락하였다가 다시 증가했다가,

 2008~2009년도에 크~게 하락하는 모습을 보이고 있습니다. 아무래도 2000년대 들어서 IT버블이 일어났기 때문에 제조량에서도 크게 증가했을까요?

 2008~2009년도에는 금융위기가 있어서 전 세계적으로 큰 풍파를 맞았었죠. 일명 '서브프라임 모기지' 사태라고도 하죠.

 아무래도 그런 영향이 있어서 그렇지 않을까 생각해봅니다.

 

tsdisplay(diff(eeadj), main="")

# tsdisplay()함수를 통해서 1차 차분한 후 정상성을 확인하고, PACF로부터 lag확인 결과, AR(3)정도가 적합할 것으로 보입니다.

  이를 확인해보기 위해서 유사한 ARIMA Order의 AIC값을 확인해봅니다.

fit_310<-Arima(eeadj, order=c(3,1,0))

fit_410<-Arima(eeadj, order=c(4,1,0))

fit_210<-Arima(eeadj, order=c(2,1,0))

fit_311<-Arima(eeadj, order=c(3,1,1))

fit_310$aic;fit_410$aic;fit_210$aic;fit_311$aic;

#즉, 전에 알아보았듯이 ARIMA(p,d,q)모형을 알아보았었죠,

#자기회귀누적이동평균모형(Autoregressive integrated moving average model), 즉 1번 차분하고, p기를 2~4까지로 잡고,

 q를 0~1까지 잡은 것입니다. 무엇이 적절한지 비정상시계열에서 forecasting을 위해 하는 작업입니다..

 여기서는 p=2, d=1, q=0일때 AIC값이 가장 높네요.
[1] 979.3314
[1] 978.9048
[1] 996.6795
[1] 978.1664
 

다음으로 넘어가서,

Acf(residuals(fit_311))

#ARIMA(3,1,1)모델의 ACF plot에서 잔차가 whith noise임을 확인합니다.

#그리고 portmanteau test를 통해서 잔차가 white noise임을 재확인 합니다.

#portmanteau test는 귀무가설이 well speicified 되었다/ 대립가설은 loosely specified입니다.

Box.test(residuals(fit_311), lag=24, fitdf=4, type="Ljung")

        Box-Ljung test

data:  residuals(fit_311)
X-squared = 20.496, df = 20, p-value = 0.4273

#그 결과, p값이 0.4273으로 크죠? 따라서 귀무가설을 기각할수 없습니다. 따라서 well specified 된 것을 알 수 있죠. 따라서 white noise임을

 재 확인합니다.

#이제 모델을 이용한 예측을 해봅니다.

#auto.arima가 위의 과정과 같은 결과를 얻는지 확인해봅니다.

library(forecast)
fit_eeadj<-auto.arima(eeadj)
plot(forecast(fit_eeadj))

# 매우 흡사합니다~ 예측도 되었죠.

# 앞선 시간에 했던 lynx도 똑같은 방식으로 해볼 수 있습니다.

AR model 의 data인 lynx를 auto.arima에 적용해봅니다.

fit_lynx<-auto.arima(lynx)

plot(lynx_fit_lynx$residuals)

lines(lynx, col=2)

summary(fit_lynx)

Series: lynx 
ARIMA(2,0,2) with non-zero mean

Coefficients:
         ar1      ar2      ma1      ma2  intercept
      1.3421  -0.6738  -0.2027  -0.2564  1544.4039
s.e.  0.0984   0.0801   0.1261   0.1097   131.9242

sigma^2 estimated as 761965:  log likelihood=-932.08
AIC=1876.17   AICc=1876.95   BIC=1892.58

Training set error measures:
                       ME          RMSE       MAE         MPE          MAPE       MASE
Training set -1.608903  853.5488  610.1112  -63.90926  140.7693  0.7343143
                    ACF1
Training set -0.01267127

# 요렇게 비슷하게 할 수 있죠.

#오늘은 시계열분석에대해서 실습을 한번 더 해봤는데요~ 전보다 시계열 들어오고나서 많이 좀 어려워진 느낌이 들죠?..

#저도 그렇습니다. ㅋㅋㅋㅋ 익숙치 않아서 그래요. 하지만 갈길이 멀기에..

#다음 시간에는 다차원척도법을 알아보도록 하겠습니다.

#BYEBYE