📖 시계열 분석을 하는데 있어 기본적인 요소들을 공부해봅시다.
par(mfrow=c(2,3))
options(family="NanumGothic")
# 백색잡음 시계열 -----------------------------------------------------------
y_t <- rnorm(100, mean=0, sd=1)
ts.plot(y_t, main="백색잡음(White Noise)",xlab="시간(time)", col="blue", lwd=2)
abline(h=0)
# 추세를 갖는 시계열--------------------------------------------------
e <- rnorm(100, 0, 10)
trend_yt <- 0.5*seq(1,100) + e
ts.plot(trend_yt, lwd=2, col="blue", main="추세를 갖는 시계열")
abline(a=0, b=0.5)
# 랜덤워크 시계열 -----------------------------------------------------------
e <- rnorm(100)
rw_yt <- cumsum(e)
ts.plot(rw_yt, lwd=2, col="blue", main="랜덤워크 시계열")
abline(h=0)
# 분산이 커지는 시계열 -----------------------------------------------------------
e <- sapply(1:100, function(x) rnorm(1, sd = 2*x/100))
var_yt <- 0.1*seq(1,100) + e
ts.plot(var_yt, lwd=2, col="blue", main="추세와 분산이 커지는 시계열")
abline(a=0, b=0.1)
# 평균이 증가하는 시계열 -----------------------------------------------------------
mean_yt <- c(rnorm(50), rnorm(50, 5))
ts.plot(mean_yt, lwd=2, col="blue", main="평균이 변화하는 시계열")
abline(a=0, b=0)
abline(a=5, b=0)
# 분산이 변화하는 시계열 -----------------------------------------------------------
var_change_yt <- c(rnorm(150,0,1), rnorm(200,0,10), rnorm(150,0,1))
ts.plot(var_change_yt, lwd=2, col="blue", main="분산이 변화하는 시계열")
abline(a=0, b=0)
1) 로그변환: log(Yt) 을 통해 분산이 커지는 경향(혹은 변동폭이 일정하지 않는 경향)을 갖는 시계열을 안정화 시킴.
2) 차분: diff(Yt) 을 통해 차분을 하계 되면 추세를 제거하는 효과를 거둠.
3) 계절차분: diff(Yt,s) 을 통해 계절 차분을 하계 되면 계절추세를 제거하는 효과를 거둠.
white_noise_zero <- arima.sim(model = list(order = c(0, 0, 0)), n = 250)
white_noise_mean_sd <- arima.sim(model = list(order = c(0, 0, 0)), n = 250, mean =5, sd =2)
par(mfrow=c(1,2))
ts.plot(white_noise_zero, ylim=c(-10,10))
ts.plot(white_noise_mean_sd, ylim=c(-10,10))
추세/순환/계절/불규칙
요소로 분해하는 기법입니다.yt=Tt+Ct+St+It
beer <- window(fpp2::ausbeer,start=1992,end=c(2007,4)) #data
par(mfrow = c(2,2))
decompose_beer <- decompose(beer, "additive")
plot(beer, main = '원 데이터')
plot(decompose_beer$seasonal, main = '계절성')
plot(decompose_beer$trend, main = '추세성')
plot(decompose_beer$random, main = '불규칙성')
yt=Tt×Ct×St×It
par(mfrow = c(2,2))
decompose_beer <- decompose(beer, "multiplicative")
plot(beer, main = '원 데이터')
plot(decompose_beer$seasonal, main = '계절성')
plot(decompose_beer$trend, main = '추세성')
plot(decompose_beer$random, main = '불규칙성')
STL은 다양한 상황에서 사용할 수 있는 강력한 시계열 분해 기법입니다. STL은 “Seasonal and Trend decomposition using Loess(Loess를 사용한 계절성과 추세 분해)”의 약자입니다. 여기에서 Loess는 비선형 관계를 추정하기 위한 기법입니다.
장점
1) SEAT와 X11과는 다르게, STL은 월별이나 분기별 데이터를 포함하여 어떤 종류의 계절성도 다룰 수 있습니다.
2) 계절적인 성분이 시간에 따라 변해도 괜찮습니다. 계절성분의 변화율을 사용자가 조절할 수 있습니다.
3) 추세-주기의 매끄러운 정도를 사용자가 조절할 수 있습니다.
4) 가끔 있는 이상값이 추세-주기와 계절성분에 영향을 주지 않게 만들 수 있습니다. 하지만 이상값은 나머지 성분(remainder)에 영향을 줄 것입니다.
window(fpp2::ausbeer,start=1992,end=c(2007,4)) %>%
stl(t.window=13, s.window="periodic", robust=TRUE) %>%
autoplot()
# t.window : 추세-주기를 추정할 때 사용할 연이은 관측값의 개수
# s.window : 계절성분에서 각 값을 추정할 때 연이은 관측값의 개수
# 위 두 매개변수 모두 값이 홀수이어야하며, 값이 작을수록 더 급격하게 변합니다.
gglagplot(beer) +
guides(colour=guide_legend(title="분기")) +
ggtitle("")
ggAcf(beer)
ggPacf(beer)
# 위 4가지 기법 시각화-------------------------------------
ggplot2::autoplot(beer) +
autolayer(meanf(beer, h=5),
series="평균", PI=FALSE) +
autolayer(naive(beer, h=5),
series="단순", PI=FALSE) +
autolayer(snaive(beer, h=5),
series="계절성 단순", PI=FALSE) +
autolayer(rwf(beer, h=5, drift = TRUE),
series="표류", PI=FALSE) +
labs(title = '분기별 맥주 생산량 예측값',
x = '연도', y = '단위:백만 리터') +
guides(color = guide_legend(title = '예측'))
beer <- window(fpp2::ausbeer,start=1992,end=c(2007,4))
# 원 데이터는 2010 2분기까지 존재
tseries::adf.test(beer) #단위근 검정
## 귀무가설(H0): 자료에 단위근이 존재한다.
## 대립가설(H1): 자료에 단위근이 존재하지 않아 시계열 자료가 정상성을 만족한다.
> tseries::adf.test(beer)
Augmented Dickey-Fuller Test
data: beer
Dickey-Fuller = -3.8813, Lag order = 3, p-value = 0.02061
alternative hypothesis: stationary
ndiffs(데이터) : 몇 번의 차분을 해야할지 정해준다.
nsdiffs(데이터) : 몇 번의 계절 차분을 해야할지 정해준다.
auto.arima
함수를 활용하여 여러 모델의 AICc를 확인하고 가장 낮은 AIC를 가진 모델을 채택하도록 하겠습니다.auto_arima_beer <- auto.arima(beer, trace=TRUE, test="adf", ic="aicc")
summary(auto_arima_beer)
checkresiduals(auto_arima_beer) # 잔차 진단
#Ljung-BOX test까지 해줘요!
> auto_arima_beer <- auto.arima(beer, trace=TRUE, test="adf", ic="aicc")
ARIMA(2,0,2)(1,1,1)[4] with drift : Inf
ARIMA(0,0,0)(0,1,0)[4] with drift : 511.942
ARIMA(1,0,0)(1,1,0)[4] with drift : 490.9862
ARIMA(0,0,1)(0,1,1)[4] with drift : 481.6034
ARIMA(0,0,0)(0,1,0)[4] : 510.7779
ARIMA(0,0,1)(0,1,0)[4] with drift : 508.715
ARIMA(0,0,1)(1,1,1)[4] with drift : 482.738
ARIMA(0,0,1)(0,1,2)[4] with drift : 482.4688
ARIMA(0,0,1)(1,1,0)[4] with drift : 490.5411
ARIMA(0,0,1)(1,1,2)[4] with drift : 484.9341
ARIMA(0,0,0)(0,1,1)[4] with drift : 483.4529
ARIMA(1,0,1)(0,1,1)[4] with drift : 483.9154
ARIMA(0,0,2)(0,1,1)[4] with drift : 483.822
ARIMA(1,0,0)(0,1,1)[4] with drift : 481.7344
ARIMA(1,0,2)(0,1,1)[4] with drift : 486.2001
ARIMA(0,0,1)(0,1,1)[4] : 488.6937
Best model: ARIMA(0,0,1)(0,1,1)[4] with drift
> summary(auto_arima_beer)
Series: beer
ARIMA(0,0,1)(0,1,1)[4] with drift
Coefficients:
ma1 sma1 drift
-0.2593 -0.7313 -0.3983
s.e. 0.1208 0.0938 0.0991
sigma^2 estimated as 154.9: log likelihood=-236.44
AIC=480.88 AICc=481.6 BIC=489.25
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.4749186 11.7439 8.735055 -0.1229367 2.011984 0.610843
ACF1
Training set -0.01387772
> checkresiduals(auto_arima_beer)
Ljung-Box test
data: Residuals from ARIMA(0,0,1)(0,1,1)[4] with drift
Q* = 2.787, df = 5, p-value = 0.7328
Model df: 3. Total lags used: 8
fore_model1 <- forecast(auto_arima_beer, h = 10)
fore_model1
autoplot(fore_model1, lwd = 5) + autolayer(fore_model1$fitted) +
autolayer(window(fpp2::ausbeer,start=c(2008,1),end=c(2010,2)))
> fore_model1
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2008 Q1 427.6151 411.6671 443.5631 403.2248 452.0054
2008 Q2 383.9695 367.4939 400.4450 358.7723 409.1666
2008 Q3 399.9629 383.4874 416.4385 374.7658 425.1601
2008 Q4 478.5510 462.0755 495.0266 453.3539 503.7482
2009 Q1 423.6286 406.6047 440.6524 397.5929 449.6642
2009 Q2 382.3762 365.3161 399.4362 356.2850 408.4673
2009 Q3 398.3696 381.3096 415.4297 372.2785 424.4608
2009 Q4 476.9577 459.8977 494.0178 450.8666 503.0489
2010 Q1 422.0353 404.4451 439.6254 395.1334 448.9371
2010 Q2 380.7829 363.1576 398.4081 353.8274 407.7383