11.1 복잡한 계절성
지금까지, 분기별이나 월별 데이터 같이 비교적 단순한 계절성 패턴을 다루었습니다. 하지만, 더 높은 빈도의 시계열에서는 종종 더욱 복잡한 계절성 패턴이 나타납니다. 예를 들면, 일별 데이터에는 일별 패턴이 있을 수 있고 연간 패턴도 있을 수 있습니다. 시간별 데이터에는 보통 3가지 종류의 계절성이 나타납니다: 일별 패턴, 주별 패턴, 그리고 연간 패턴. 주별 데이터 평균적으로 \(365.25/7\approx 52.179\)의 계절성 주기를 갖는 연간 패턴이 보통 나타날 수 있기 때문에 주별 데이터를 다루는 것은 어려울 수 있습니다.
이러한 다중 계절성 패턴은 높은 빈도로 기록된 데이터에서 흔하게 나타납니다. 다중 계절성 패턴이 나타날 수 있는 다른 예제에는 콜센터의 전화량 기록, 일별 병원 방문자 수, ATM에서 현금 인출 요청 수, 전력과 수도 사용량, 웹 사이트 방문자 등이 있습니다.
지금까지 다룬 대부분의 기법은 이러한 복잡한 계절성을 다룰 수 없습니다. R의 ts 클래스는 보통 정수 값을 가진다고 가정하는 1가지 계절성만 다룰 수 있습니다.
이러한 시계열을 다루기 위해, 다중 계절성 시계열을 다루는 msts 클래스를 사용할 것입니다. 이것을 이용하여 관련 있을 수 있는 주기 값 전체를 정할 수 있습니다. 이 클래스는 정수가 아닌 주기도 다룰 수 있을 정도로 유연합니다.
유연성이 있더라도, 이러한 주기 전부를 넣을 필요가 없을 수 있습니다. 데이터에서 가장 잘 나타나는 하나만 필요할 수 있습니다. 예를 들면, 180일 데이터만 가지고 있을 때, 연도 계절성을 무시할 수 있습니다. 자연 현상을 측정한 데이터라면(예: 온도), 아마도 안전하게 주별 계절성도 무시할 수 있을 것입니다.
그림 11.1 의 윗 부분은 33주 동안 주중 오전 7시부터 오후 9시 5분까지 5분 간격으로 소매 은행에 걸려오는 통화의 수를 나타냅니다. 아래 부분은 같은 시계열의 처음 3주를 나타냅니다. 주기 169인 (하루에 5분 시간 간격이 169개 있습니다) 일별 계절성 패턴이 강하게 나타나고, 주기 \(169 \times 5=845\)인 주별 계절성 패턴이 약하게 나타납니다. (월요일에는 통화의 수가 주중 나머지 날보다 높게 나타나는 경향이 있습니다.) 더 긴 시계열 데이터를 이용할 수 있다면, 연간 계절성 패턴도 관측할 수 있겠습니다.
p1 <- autoplot(calls) +
ylab("통화량") + xlab("주차") +
scale_x_continuous(breaks=seq(1,33,by=2))
p2 <- autoplot(window(calls, end=4)) +
ylab("통화량 ") + xlab("주차") +
scale_x_continuous(minor_breaks = seq(1,4,by=0.2))
gridExtra::grid.arrange(p1,p2)
Figure 11.1: 북미 어떤 상업 은행에서 주중 오전 7시부터 오후 9시 5분 사이에 처리된 5분 통화량. 위쪽 패널은 2003년 3월 3일부터 시작하여 164일 동안의 데이터를 나타냅니다.
다중 계절성 주기를 고려하는 STL
mstl() 함수는 여러 계절성을 고려하도록 stl()을 변형한 것입니다. 결과값으로 다중 계절성 성분, 추세, 나머지 성분을 돌려줍니다.
calls %>% mstl() %>%
autoplot() + xlab("통화량")
Figure 11.2: 통화량 데이터에 다양한 STL 기법을 적용한 것.
두 가지 계절성 패턴이 나타납니다. 하루에 대한 것(3번째 패널)과 주에 대한 것(4번째 패널)입니다. 이 그래프를 적절하게 해석하려면, 수직 눈금을 살펴보는 것이 중요합니다. 이 경우에는, 데이터에서 작은 추세가 나타나고 주별 계절성이 약해서 다른 성분에 비해 추세와 주별 계절성이 상대적으로 좁은 범위로 나타납니다.
이러한 분해를 계절성 단순(naïve) 기법을 이용하여 계절 성분 각각을 예측할 때와 ETS (또는 몇몇 사용자가 정한 다른 기법)를 이용하여 계절성으로 조정된 데이터를 예측할 때 사용할 수 있습니다. stlf() 함수가 이 작업을 자동으로 할 것입니다.
calls %>% stlf() %>%
autoplot() + xlab("주차")
Figure 11.3: 통화량 데이터에 다양한 STL 기법을 적용한 것.
다중 계절성 주기를 고려하는 동적 조화 회귀
앞의 장에서 한 것처럼 푸리에 항을 이용하여 다중 계절성을 다룰 수 있습니다(5.4 절과 9.5 절을 참조하시길 바랍니다). 다중 계절성이 있기 때문에, 각 계절성 주기에 푸리에 항을 추가할 필요가 있습니다. 이 경우에, 계절성 주기는 169와 845이고, 따라서 \(k=1,2,\dots\)에 대한 푸리에 항의 형태는 다음과 같습니다.
\[
\sin\left(\frac{2\pi kt}{169}\right), \quad
\cos\left(\frac{2\pi kt}{169}\right), \quad
\sin\left(\frac{2\pi kt}{845}\right), \quad \text{and} \quad
\cos\left(\frac{2\pi kt}{845}\right),
\]
fourier() 함수로 이러한 항을 생성할 수 있습니다.
ARMA 오차 구조를 고려하는 동적 조화 회귀 모델을 맞출 것입니다. AICc를 최소화하기 위해 각 계절성 주기에 대한 푸리에 항의 전체 숫자가 선택될 것입니다. 예측값과 예측 구간이 양수로 유지되는지 확인하기 위해 로그 변환(lambda=0)을 사용할 것입니다.
fit <- auto.arima(calls, seasonal=FALSE, lambda=0,
xreg=fourier(calls, K=c(10,10)))
fit %>%
forecast(xreg=fourier(calls, K=c(10,10), h=2*169)) %>%
autoplot(include=5*169) +
ylab("통화량") + xlab("주차")
Figure 11.4: 동적 조화 회귀를 통화량 데이터에 적용하여 얻은 예측값.
이것은 40개의 매개변수를 포함하는 큰 모델입니다: 4개의 ARMA 계수, 주기 169에 대한 20개의 푸리에 계수, 주기 845에 대한 16 개의 푸리에 계수. 주기 169의 항과 겹치기 때문에 (\(845=5\times169\)) 주기 845에 대한 모든 푸리에 항은 사용하지 않습니다.
TBATS 모델
De Livera, Hyndman, & Snyder (2011) 가 개발한 또 다른 접근 방식에서는 지수 평활 상태 공간 모델(exponential smoothing state space model)과 박스-칵스(Box-Cox) 변환을 고려하는 푸리에 항의 조합을 완벽히 자동화된 방식으로 사용합니다. 다른 자동화된 모델링 체계가 그렇듯이, 안 좋은 결과를 내는 경우가 있을 수 있습니다만, 몇 가지 상황에서 유용한 접근 방식이 될 수 있습니다.
조화 회귀 항은 계절성 패턴이 변하지 않고 주기적으로 반복해서 나타나도록 합니다만, TBATS 모델은 계절성이 TBATS 모델에서 시간에 따라 느리게 변할 수 있다는 점에서 동적 조화 회귀와는 다릅니다. 하지만 TBATS 모델의 한 가지 단점은 특별히 긴 시계열을 추정하는데 느릴 수 있다는 것입니다. 따라서 시간을 아끼기 위해 calls 데이터의 일부분을 다룰 것입니다.
calls %>%
subset(start=length(calls)-2000) %>%
tbats() -> fit2
fc2 <- forecast(fit2, h=2*169)
autoplot(fc2, include=5*169) +
ylab("통화량") + xlab("주차")
Figure 11.5: TBATS 모델을 통화량 데이터에 적용하여 얻은 예측값.
여기에서 나타나는 예측구간은 너무 넓습니다. 안타깝게도 TBATS 모델에서는 종종 이런 상황이 벌어지기도 합니다.
공변량이 있는 복잡한 계절성
동적 조화 회귀(dynamic harmoic regression) 모델에는 공변량(covariate)이 들어갈 수 있긴해도, TBATS 모델은 공변량을 허용하지 않습니다. 이러한 모델을 보통 전력 수요 모델링에 적용합니다.
그림 11.6은 2014년에 호주 빅토리아 주에서 매 30분마다 전력 수요를 빅토리아 주의 가장 큰 도시인 멜버른에 대해 같은 기간의 온도에 따라 나타낸 것입니다.
colnames(elecdemand) <- c("수요","WorkDay","기온")
autoplot(elecdemand[,c("수요","기온")],
facet=TRUE) +
scale_x_continuous(minor_breaks=NULL,
breaks=2014+
cumsum(c(0,31,28,31,30,31,30,31,31,30,31,30))/365,
labels=month.abb) +
xlab("시간") + ylab("")
Figure 11.6: 2014년 호주 빅토리아의 30분 간격 전력 수요와 해당 기온.
colnames(elecdemand) <- c("Demand","WorkDay","Temperature") 전력 수요를 온도에 따라 그리면 (그림 11.7), 기온이 낮을 때는 (난방 때문에) 수요가 증가하고, 기온이 높을 때는 (냉방 때문에) 수요가 증가하는 둘 사이의 비선형 관계가 나타납니다.
elecdemand %>%
as.data.frame() %>%
ggplot(aes(x=Temperature, y=Demand)) + geom_point() +
xlab("기온 (섭씨)") +
ylab("수요 (GW)")
Figure 11.7: 빅토리아주의 30분 간격 전력 수요를 빅토리아주에서 가장 큰 도시인 멜버른의 기온과 함께 나타낸 것.
(섭씨 18도에 있는 매듭을 포함하는) 온도의 단계별 선형 함수와 일별 계절성 패턴을 허용하기 위한 회귀 항을 고려하는 회귀 모델을 맞출 것입니다.
cooling <- pmax(elecdemand[,"Temperature"], 18)
fit <- auto.arima(elecdemand[,"Demand"],
xreg = cbind(fourier(elecdemand, c(10,10,0)),
heating=elecdemand[,"Temperature"],
cooling=cooling))예측변수(predictor variable)의 미래 값이 필요하기 때문에 이러한 모델로 예측하는 것은 어려운 일입니다. 푸리에 항의 미래값은 계산하기 쉽습니다만, 미래 온도값은 물론 알 수 없습니다. 한 주 앞까지만 예측하는데 관심이 있다면, 기상학적인 모델로부터 얻은 기온 예측값을 사용할 수 있습니다. 다른 방법으로는, 시나리오 예측을 사용하여 (4.5 절), 가능한 온도 패턴을 넣는 방법도 있겠습니다. 대신에, 시나리오 예측을 사용하여 가능한 기온 패턴을 넣을 수도 있습니다. 다음 예제에서는, 가능한 미래 수요값을 생성하기 위해 마지막 2주의 온도 값을 반복하여 사용합니다.
temps <- subset(elecdemand[,"Temperature"],
start=NROW(elecdemand)-2*48+1)
fc <- forecast(fit,
xreg=cbind(fourier(temps, c(10,10,0)),
heating=temps, cooling=pmax(temps,18)))
autoplot(fc, include=14*48)
Figure 11.8: 동적 조화 회귀 모델을 30분 간격 전력 수요 데이터에 적용하여 얻은 예측값.
단기 예측값이 그럴듯해보여도, 이건 복잡한 과정을 설명하기 위해 대강 만든 모델입니다. 잔차(residual)는 이 모델이 잡아내지 못한 정보가 여전히 많다는 것을 나타냅니다.
checkresiduals(fc)
Figure 11.9: 동적 조화 회귀 모델에 대한 잔차 진단.
#>
#> Ljung-Box test
#>
#> data: Residuals from Regression with ARIMA(5,1,4) errors
#> Q* = 738412, df = 3455, p-value <2e-16
#>
#> Model df: 49. Total lags used: 3504이 모델의 좀 더 세련된 형태는 훨씬 더 나은 예측값을 제공합니다. 이 내용은 Hyndman & Fan (2010) 와 Fan & Hyndman (2012) 에 있습니다.