모델의 차수를 찾은 다음(즉, p,d,q 값), 다음과 같은 매개변수 c,ϕ1,...ϕp,θ1,…,θq 을 추정할 필요가 있다.
R에서 ARIMA 모델을 계산시, MLE를 사용한다. 이 방법은 관찰한 데이터를 얻는 확률을 최대화하는 매개변수의 값을 찾는다.
다음은 수리통계학에서 나오는 정의이다.
가능도함수최대가능도 추정량
ARIMA 모델에서 MLE는 다음과 같은 양을 최소화하는 최소제곱(least squares)추정과 비슷하다.
ARIMA 모델이 회귀 모델을 추정하는 것보다 훨씬 더 복잡하고, 서로 다른 소프트웨어가 서로 다른 추정 기법과 최적화 알고리즘을 사용하기 때문에 살짝 다른 결과를 낼 수 있다는 것에 주목하자.
실제로, R은 데이터의 로그 가능도(log likelihood) 값을 알려줄 것이다. 즉, 추정한 모델에서 나온ㄱ ㅘㄴ측 데이터의 확률의 로그를 말한다. 다음은 수리통계학에서 나오는 log likelihood 이다.
로그 가능도함수
정보 기준
회귀에서 예측변수(predict)를 선택할 때 유용했던 아카이케(Akaike)의 정보 기준이 ARIMA 모델에서 차수를 결정할 때도 유용하다.
AIC는 다음과 같이 쓸 수 있다.
위에서 L은 데이터의 가능도, c is not 0 이면 k=1이고, c is 0 이면 k=0이다. 괄호 안의 마지막 항이 (σ^2와 잔차(residual)의 분산을 포함하는) 모델의 매개변수 개수라는 것에 주목하자.
ARIMA 모델에 대해, 수정된 AIC는 다음과 같다.
또한 베이지안 정보 기준은 다음과 같이 쓸 수 있다.
AIC,AICc 혹은 BIC를 최소화하여 좋은 모델을 얻을 수 있다. 여기서는 AICc를 사용한다.
이러한 정보 기준이 모델의 적절한 차분 차수(d)를 고를 때 별로 도움이 되지 않는 경향이 있고, p와 q값을 고를 때만 도움이 된다는 것은 중요한 점이다. 차분을 구하는 것을 통해 likelihood를 계산하는 데이터가 바뀌기 때문에, 서로 다른 차수로 차분을 구한 모델의 AIC 값을 비교할 수 없게 된다. 그렇기에 d를 고르기 위해 다른 방법을 사용해야 하고, 그 후 p와 q를 고르기 위해 AICc를 사용할 수 있다.
위 그래프는 미국 소비 지출 분기별 백분율 변화를 나타낸다. 분기별 시계열이지만, 계절성 패턴이 나타나지 않는 것 같다.
따라서 비-계절성 ARIMA 모델로 맞춰보자.
다음 코드는 모델을 자동으로 선택할 때 사용한다.
fit <- auto.arima(uschange[,"Consumption"], seasonal=FALSE)
#> Series: uschange[, "Consumption"]
#> ARIMA(1,0,3) with non-zero mean
#>
#> Coefficients:
#> ar1 ma1 ma2 ma3 mean
#> 0.589 -0.353 0.085 0.174 0.745
#> s.e. 0.154 0.166 0.082 0.084 0.093
#>
#> sigma^2 estimated as 0.35: log likelihood=-164.8
#> AIC=341.6 AICc=342.1 BIC=361
아래는 ARIMA(1,0,3) 모델이다.
여기에서 c=0.745*(1-0.589)=0.307 이고 Et는 다음과 같이 0.592=Route(0.350)이다.
이러한 표준 편차를 갖는 백색잡음(white noise)이다.
fit %>% forecast(h=10) %>% autoplot(include=80) +
ggtitle("0이 아닌 평균을 가지는 ARIMA(2,0,2)로부터 얻은 예측값") + ylab("소비")
미국 소비 지출의 분기별 백분율 변동의 예측값
다음과 같은 예측값을 나타내었다.
ARIMA 모델 이해하기
auto.arima() 함수는 유용하지만, 모든 입력을 자동으로 결정하게 두면 약간 위험할 수 있다.
그리고 자동으로 모델을 고르게 두더라도 모델이 대략적으로 작동하는 방식은 공부해볼만 하다.
상수 c는 이러한 모델에서 얻은 장기 예측값에 중요한 영향을 준다.
c=0이고d=0이면, 장기 예측값이 0에 가까워질 것이다.
c=0이고d=1이면, 장기 예측값이 0이 아닌 상수에 가까워질 것이다.
c=0이고d=2이면, 장기 예측값이 직선 형태로 나타나게 될 것이다.
c≠0이고d=0이면, 장기 예측값이 데이터의 평균에 가까워질 것이다.
c≠0이고d=1이면, 장기 예측값이 직선 형태로 나타나게 될 것이다.
c≠0이고d=2이면, 장기 예측값이 2차 곡선 추세로 나타나게 될 것이다.
d 값은 예측 구간(prediction interval)에도 영향을 준다. d값이 클 수록, 예측 구간의 크기가 더욱 급격하게 늘어난다.
d=0에서, 장기 예측 표준 편차가 과거 데이터의 표준 편차에 가까워질 것이고, 따라서 모든 예측 구간은 실제적으로 같게 될 것이다.
위 그래프는 이러한 행동을 나타낸다. 여기에서 d=0이고 c is not equal 0이다. 위 그래프에서, 예측 구간은 마지막 몇 개의 예측 수평선(forecast horizon)에 대한 경우가 거의 같고, 점 예측값(point forecast)은 데이터의 평균과 같다.
p 값은 데이터에서 주기(cycles)가 나타날 때 중요하다. 주기적 예측값을 얻기 위해서는, 매개변수에 대한 몇 가지 추가적인 조건과 함께 p>=2 이어야 한다. AR(2) 모델의 경우, ϕ21+4ϕ2<0이면 주기적인 행동이 나타난다.
이 경우에, 주기(cycle)의 평균 기간은 다음과 같다.
ACF와 PACF 그래프
보통은 단순하게 시간 그래프(time plot)만 보고나서 어떤 p와 q값이 데이터에 맞는지 이야기할 수 없다. 하지만, 적절한 pp와 q값을 결정하기 위해서 때때로 ACF 그래프와 PACF 그래프를 이용하면 가능하다.
서로 다른 k값에 대해, yt와 yt-k의 관계를 측정하는 자기상관값(autocorrelation)을 나타내는 ACF그래프를 다시 떠올려보자. yt와 yt-1이 상관관계가 있다면, yt-1과 yt-2에도 상관관계가 있어야 한다. 하지만 yt와 yt-2는(yt를 예측하는데 사용될 수 있는) yt-2에 담긴 어떤 새로운 정보 때문이 아니라 단순히 두 값 모두 yt-1과 관련이 있기 때문에 상관관계를 가질 수도 있다.
이러한 문제를 극복하기 위해, 부분 자기상관값들(partial autocorrelations)을 이용할 수 있다.
이 값은 시차 1,2,3,.....,k-1 의 효과를 제거한 후의 yt와 yt-k 사이의 관계를 측정한다. 그래서 첫 번째 부분 자기상관은 제거할 부분이 없어서 첫 번째 자기상관과 같다. 각 부분 자기상관은 자기회귀 모델의 마지막 계수처럼 측정할 수 있다. 구체적으로, k번째 부분 자기상관 계수 αk를 계산하는 알고리즘들이 존재한다.
다음 그림들은 미국 소비 데이터의 ACF와 PACF 그래프를 나타낸다.
부분 자기상관값은 보통의 자기상관값처럼 +-1.96/Route(T)로 같은 임계 값을 갖는다.
ggACF(uschange[,"Consumption"])
미국 소비 데이터에서 분기별 백분율 변동의 ACF
ggPacf(uschange[,"Consumption"])
미국 소비 데이터에서 분기별 백분율 변동의 PACF
데이터가 ARIMA(p,d,0)나 ARIMA(0,d,q) 모델에서 왔다면, p나 q 값을 결정할 때 ACF와 PACF 그래프가 유용할 수 있다.
p와 q가 모두 양수라면, 적절한 p와 q값을 찾을 때 이러한 그래프는 도움이 되지 않는다.
차분을 구한 데이터의 ACF와 PACF 그래프가 다음과 같은 패턴을 나타내면, 데이터는 ARIMA(p,d,0) 모델을 따를 수도 있다.
ACF가 지수적으로 감소하거나 사인 함수 모양인 경우
PACF 그래프에서 시차pp에 뾰족한 막대가 유의미하게 있지만, 시차pp이후에는 없을 때.
차분을 구한 데이터의 ACF와 PACF 그래프가 다음과 같은 패턴을 나타내면 데이터는 ARIMA(0,d,q)모델을 따를 수도있다.
PACF가 지수적으로 감소하거나 사인 함수 모양인 경우;
ACF 그래프에서 시차qq에 뾰족한 막대가 유의미하게 있지만, 시차qq이후에는 없을 때.
위 ACF그래프에서 시차 4에서 나타나는 거의 유의미한 뾰족한 막대 이전에 3개의 뾰족한 막대가 있는 것을 알 수 있다. PACF에서 3개의 유의미한 뾰족한 막대가 있고, 그 이후에는 (22 뒤처짐에서 경계를 한 번 벗어나는 것을 제외하고) 유의미하게 경계를 벗어나는 것을 무시할 수 있다.
결국, 뾰족한 막대가 유의미하게 될 확률은 1/20이고, 각 그래프에서 시차 22에서 나타나는 뾰족한 막대를 나타냈다.
첫 3개의 뾰족한 막대에서 나타나는 패턴은 PACF가 감소하는 경향을 통해 ARIMA(3,0,0)에서 예측했던 것이다.
그래서 ACF와 PACF을 보면 이 경우에는 ARIMA(3,0,0) 모델이 적절할 수 있다고 생각할 수 있다.
회귀에서 목표 에상 변수(forecast variable)의 과거 값을 이용하는 대신에, 이동 평균 모델은 회귀처럼 보이는 모델에서 과거 예측 오차(forecast error)을 이용한다.
여기서 E(t)는 백색잡음(White noise)이다. 이것을 q차 이동 평균 모델인 MA(q) 모델이라고 부르자.
물론, E(t)의 값을 관찰하지 않기 때문에, 실제로는 일반적인 회귀가 아니다.
Y(t)의 각 값을 과거 몇 개의 예측 오차(forecast error)의 가중 이동 평균으로 생각할 수 있다는 것에 주목하자. 하지만 이동 평균 평활과 헷갈리면 안된다. 이동 평균 모델은 미래 값을 예측할 때 사용한다. 이동 평균 평활은 과거 값의 추세-주기를 측정할 때 사용한다.
매개변수를 다르게 설정한 이동 평균 모델로부터 얻은 데이터의 두가지 예
MA(1): Y(t)=20+E(t)+0.8E(t-1)
MA(2): Y(t)=E(t)-E(t-1)+0.8E(t-2)
*두 가지 경우 모두, E(t)는 평균이 0이고 분산이 1인 정규 분포를 따르는 백색 잡음(White noise)이다.
위 그림은 MA(1) 모델과 MA(2) 모델로 얻은 몇몇 데이터를 나타낸다. 매개변수 θ1,…,θqθ1,…,θq 을 바꾸면 다른 시계열 패턴이 나타난다. 자기회귀 모델을 이용하는 경우처럼, 오차항 E(t)의 분산은 시계열의 패턴이 아니라 눈금만 맞출 것이다.
정상성을 나타내는 어떤 AR(p) 모델을 MA(∞) 모델로 쓸 수 있다. 예를 들어, 반복하여 대입하면, AR(1) 모델에 대해 다음과 같이 나타낼 수 있다.
-1<ϕ1<1 에 대해, k 가 커질 수록 ϕ1^ k의 값이 작아질 것이다. 그래서 결국 다음과 같은 식을 얻게 된다.
이는 MA(∞)의 과정이다.
MA 매개변수에 대한 몇몇 제한조건을 도입하면 반대 결과도 성립한다. 그러면 MA 모델을 가역적(invertible)이라고 부른다. 즉, 어떤 가역적인 MA(q) 과정을 AR(∞) 과정으로 쓸 수 있다. 가역적 모델은 단순하게 MA 모델을 AR 모델로 바꿀 수 있도록 하는 것만은 아니다. 몇 가지 수학적 특징도 갖고 있다.
예를 들어보자.
MA(1) 과정
다음과 같은 식을 생각해보자. 이것은 AR(∞)로 표현하면, 가장 최근의 오차는 현재와 과거 관측값의 선형 함수로 쓸 수 있다.
|θ|>1 이면, 가중치의 시작(lag) 값이 증가함에 따라 증가하고, 따라서 더 멀리 떨어진 관측값일수록 현재 오차에 미치는 영향이 커진다.
|θ|=1 이면, 가중치가 크기에 대해서는 상수, 멀리 떨어진 관측값과 가까운 관측값 모두 동일하게 영향을 준다.
위 예 AR(1), AR(2) 경우 모두 그럴듯하지 않기 때문에, |θ|<1 가 필요하며, 결국 가장 최근 관측값이 멀리 떨어진 관측값 보다 더 큰 가중치를 갖게 된다. 따라서, |θ|<1일때 과정은 가역적(invertible)이다.
다른 모델에 대한 가역성(invertibility) 제한조건은 정상성(stationarity) 제한조건과 비슷하다.
MA(1) 모델의 경우:−1<θ1<1−1<θ1<1.
MA(2) 모델의 경우:−1<θ2<1,θ2+θ1>−1,θ1−θ2<1.
q>=3에 대해서는 더 복잡한 조건이 성립한다. 여기서, 파이썬, R에서 모델을 다룰 때, 이러한 제한조건을 처리해준다.
ARIMA 모형은 대표적인 통계적 시계열 예측 모형으로, 현재 값을 과거 값과 과거 예측 오차를 통해 설명한다. 이 모형을 적합하려면 시계열 Yt가 정상성 조건을 만족하는 정상 시계열(Stationary series)이어야 한다. (정확하게는 약정상성, weak stationary). 또한, 주어진 시계열 자료에 적합한 ARIMA(p,d,q)를 결정하는 절차는 Box-Jenkins method를 따르며, ACF(Autocorrelation function: 자기상관함수)로 시계열의 특성을 파악하고 적절한 차수의 ARIMA 모형을 선택한다.
Box-Jenkins method
정상성(Stationary)
우리는 시간의 순서에 따라 기록되지 않은 일반적 자료들을 분석할 때, Random sample들에 iid가정을 한다.
*iid: All samples are independent and identically distributed
시간에 종속되어있는 시게열은 상식적으로 iid 가정을 할 수 없다. 그래서 이러한 시계열 자료에 대해 예측 모형을 적합하고 통계적 검정을 하기 위해서는 분석을 단순화 시킬 수 있는 새로운 가정이 필요하다.
이중 가장 중요한 것이 시계열 모형의 확률적 성질이 시간에 따라 변하지 않는다고 가정하는 정상성 가정이다.
ARIMA모형은 해당 시계열이 약정상성(Weak Stationary)를 만족한다고 가정하며, 약정상성을 만족해야 좋은 fitting과 predict 성능을 보여줄 수 있다. 시계열 Yt가 다음의 세 조건을 만족할 때 약정상성을 가진다고 표현한다. 시계열 분서겡서 말하는 정상 시계열(Stationary series)은 약 정상성을 가지는 시계열을 말한다.
-평균이 모든 시점 t에서 동일하다 : 추세, 계절성, 순환성 등의 패턴이 보이지 않게 된다.
-분산이 모든 시점 t에서 동일하다: 자료 변화의 폭이 일정하게 된다.
-Yt와 Y(t-h) 간의 Covariation(즉, Yt의 자기 공분산 함수)이 모든 시점 t에 대해서 동일하다: 시간에 따라 상이한 자기상관적 패턴을 보이지 않게 된다.
즉, 풀어설명하자면 추세와 계절성이 있는 시계열은 정상성을 나타내는 시계열이 아니다. 일반적으로 정상성을 나타내는 시계열은 장기적으로 볼 때 예측할 수 있는 패턴을 나타내지 않을 것이다. (어떤 주기적인 행동이 있을 수 있더라도) 시간 그래프는 시계열이 일정한 분산을 갖고 대략적으로 평평하게 될 것을 나타낼 것이다.
그렇다면 다음 시계열 그래프중 어떤 것이 정상성을 나타낼까?
a) 거래일 200일 동안의 구글 주식 가격
b) 거래일 200일 동안의 구글 주식 가격의 일일 변동
c) 미국의 연간 파업 수
d) 미국에서 판매되는 새로운 단독 주택의 월별 판매액
e) 미국에서 계란 12개의 연간 가격(고정 달러)
f) 호주 빅토리아 주에서 매월 도살한 돼지의 전체 수
g) 캐나다 북서부 맥킨지 강 지역에서 연간 포획된 시라소니의 전체 수
h) 호주 월별 맥주 생산량
i) 호주 월별 전기 생산량
분명하게 계절성이 보이는 d),h),i)는 후보가 되지 못한다. 추세가 있고 수준이 변하는 a),c),e),f),i)도 후보가 되지 못한다. 분산이 증가하는 i)도 후보가 되지 못한다. 그렇다면 b)와 g)만 정상성을 나타내는 시계열 후보로 남았다.
얼핏 보면 시계열 g)에서 나타나는 뚜렷한 주기(cycle)때문에 정상성을 나타내는 시계열이 아닌 것처럼 보일 수 있다. 하지만 이러한 주기는 불규칙적(aperiodic)이다.
즉, 먹이를 구하기 힘들만큼 시라소니 개체수가 너무 많이 늘어나 번식을 멈춰서, 개체수가 작은 숫자로 줄어들고, 그 다음 먹이를 구할 수 있게 되어 개체수가 다시 늘어나는 식이다.
장기적으로 볼 때 , 이러한 주기의 시작이나 끝은 예측할 수 없다. 따라서 이 시계열은 정상성을 나타내느 시계열이다.
차분(differencing)
구글 주식 가격의 ACF(왼쪽), 일별 변동(오른쪽)
처음 그림에서 패널(a) 의 구글 주식 가격이 정상성을 나타내는 시계열이 아니었지만 (b)의 일별 변화는 정상성을 나타냈다느 것에 주목해보자. 이 그림은 정상성을 나타내지 않는 시계열을 정상성을 나타내도록 만드는 한 가지 방법을 나타낸다.
연이은 관측값들의 차이를 계산하는 것이다. 이를 차분(differencing)이라 부른다.
로그변환은 시계열의 분산 변화를 일정하게 만드는데 도움이 될 수 있다. 차분(differencing)은 시계열의 수준에서 나타내느 변화를 제거하여 시계열의 평균 변화를 일정하게 만드는데 도움이 될 수 있다. 결과적으로 추세나 계절성이 제거(또는 감소) 된다.
정상성을 나타내지 않는 시계열을 찾아낼 때 데이터의 시간 그래프를 살펴보는 것만큼, ACF 그래프도 유용하다. 정상성을 나타내지 않는 데이터에서는 ACF가 느리게 감소하지만, 정상성을 나타내는 시계열에서는, ACF가 비교적 빠르게 0으로 떨어질 것이다. 그리고 정상성을 나타내지 않는 데이터에서 r1은 종종 큰 양수 값을 갖는다.
차분을 구한 구글 주식 가격의 ACF는 단순히 white noise 시계열처럼 생겼다. 95% limitation 바깥에 자기상관(autocorrelation)값이 없고, 융-박스(Ljung-Box) Q* 통계는 h=10 에 대해 0.355라는 p-value를 갖는다. 이 결과는 구글 주식 가격의 일별
Box.test(diff(goog200), lag=10, type="LJung-Box")
#Box-Ljung test
#data: diff(goog200)
#X-Squared=11, df=10, p-value=0.4
변동이 기본적으로는 이전 거래일의 데이터와 상관이 없는 무작위적인 양이라는 것을 말해준다.
확률보행 모델
차분(difference)을 구한 시게열은 원래의 시계열에서 연이은 관측값의 차이이고, 다음과 같이 쓸 수 있다.
첫 번째 관측값에 대한 차분 y'1을 계산할 수 없기 때문에 차분을 구한 시계열은 T-1개의 값만 가질 것이다.
차분을 구한 시계열이 white noise이면, 원래 시계열에 대한 모델은 다음과 같이 쓸 수 있다.
여기서 Error(t)는 white noise를 의미한다. 이것을 정리하면 '확률보행(random walk)' 모델을 얻는다.
확률보행(random walk) 모델은 정상성을 나타내지 않은 데이터, 특별히 금융이나 경제 데이터를 다룰 때 널리 사용되고 있다. 확률보행에는 보통 다음과 같은 특징을 가진다.
-누가 봐도 알 수 있는 긴 주기를 갖는 상향 또는 하향 추세가 있다.
-갑작스럽고 예측할 수 없는 방향 변화가 있다.
미래 이동을 예측할 수 없고 위로 갈 확률이나 아래로 갈 확률이 정확하게 같기 때문에 확률보행 모델에서 얻어낸 예측값은 마지막 관측값과 같다. 따라서, 확률보행 모델은 단순(naive) 예측값을 뒷받침한다.
밀접하게 연관된 모델은 차분값이 0이 아닌 평균값을 갖게 한다. 그러면,
c값은 연이은 관측값의 차이의 평균이다. c가 양수이면, 평균 변화는 y(t)값에 따라 증가한다. 따라서 y(t)는 위쪽 방향으로 이동하는 경향을 나타낼 것이다. 하지만 c가 음수이면, y(t)는 아래쪽 방향으로 이동하는 경향을 나타낼 것이다.
2차 차분
가끔 차분(difference)을 구한 데이터가 정상성(stationarity)이 없다고 보일 수도 있다. 정상성을 나타내는 시계열을 얻기 위해 데이터에서 다음과 같이 한 번 더 차분을 구하는 작업이 필요할 수 있다.
이 경우에는, y(t)''는 T-2개의 값을 가질 것이다. 그러면, 원본 데이터의 '변화에서 나타나는 변화'를 모델링하게 되는 셈이다. 실제 상황에서는, 2차 차분 이상으로 구해야 하는 경우는 거의 일어나지 않는다.
계절성 차분
계절성 차분(seasonal differencing)은 관측치와, 같은 계절의 이전 관측값과의 차이를 말한다. 따라서
여기에서 m은 계절의 개수이다. m 주기 시차 뒤의 관측을 빼기 때문에 시차 m 차분이라고 부르기도 한다.
계절성으로 차분을 구한 데이터가 white nois로 보이면, 원본 데이터에 대해 적절한 모델은 다음과 같다.
이 모델에서 낸 예측값은 관련 있느 계절의 마지막 관측값과 같다. 즉, 이 모델은 계절성 단순(seasonal naive) 예측값을 나타낸다.
위 그림의 아래쪽 패널은 호주에서 팔린 A10 약물(당뇨병 약)의 월별 처방전의 수에 로그를 취하여 계절성 차분을 구한 결과를 나타낸다. 변환과 차분을 통해 시계열이 정상성을 나타내는 것처럼 보인다.
cbind("판매량 (백만 달러)"= a10,
"월별 로그 판매량"= log(a10),
"로그 눈금에서 연간 변동"=diff(log(a10),12))%>%
autoplot(facets=TRUE)+
xlab("연도")+ylab("")+
ggtitle("당뇨병 약 판매량")
보통의 차분과 계절성 차분을 구분하기 위해, 때때로 보통의 차분을 시차 1에서 차분을 구한다는 의미로 "1차 차분(first difference)"라고 부른다.
위 그림에서 나타낸 것처럼, 계절성을 나타내는 데이터를 얻기 위해 계절성 차분과 1차 차분 둘 다 구하는 것이 필요하기도 한다. 여기서는 데이터를 먼저 로그로 변환하고(두번째 패널), 그 후 계절성 차분을 계산했다(세번째 패널). 데이터에서 여전히 정상성이 보이지 않는 것 같아서, 1차 차분을 더 많이 계산했다.
어떤 차분(difference)을 구할지 정할 때는 주관적인 요소가 어느정도 들어간다. 위 약물 그림의 계절성 차분(seasonal difference)을 구한 데이터는 월별 순 전기 생산량의 계절성 차분을 구한 데이터와는 큰 차이를 나타내지 않는다. 후자의 경우, 계절성 차분을 구한 데이터로 결정했어야 했고, 차분을 더 구하지 않아도 되었다. 이에 반면해 전자의 경우에는 데이터의 정상성(stationarity)이 충분히 나타나지 않아서 차분을 더 구해야 했다. 차분을 구하는 것에 대한 몇가지 형식적인 검정을 아래에서 다루지만, 모델링 과정에서 항상 몇 가지 선택이 존재하고, 분석하는 사람마다 다른 선택을 할 수 있다.
y'(t)=y(t)-y(t-m)가 계절성 차분(seasonal difference)을 구한 시계열을 나타낸다면, 두 번 차분한 시계열은 다음과 같다.
계절성 차분과 1차 차분을 둘 다 적용할 때, 어떤 것을 먼저 적용하더라도 차이는 없다. 결과는 같을 것이다. 하지만, 데이터에 계절성 패턴이 강하게 나타나면, 계절성 차분을 먼저 계산하는 것을 추천한다. 왜냐면, 때때로 결과 시계열에서 정상성이 나타나기도 해서 이런 경우 1차 차분을 구할 필요가 없게 되기 때문이다. 1차 차분을 먼저 계산했다면, 여전히 남아있는 계절성이 나타날 것이다.
차분을 구했다면, 차분 값이 해석 가능할 것이라는 것은 중요하다. 첫 번째 차분값은 한 관측값과 그 다음 관측값 사이의 변화이다. 계절성 차분값은 한 해와 그 다음 해 사이의 변화이다. 다른 시차값(lagged value)은 직관적으로 해석하기가 쉽지 않기 때문에 사용하지 않는 것이 좋다.
단위근검정
단위근검정(unit root tests)은 더 객관적으로 차분을 구하는 것이 필요할 지 결정하기 위해 사용하는 한 가지 방법이다. 차분을 구하는 것이 필요한지 결정하는 상황을 위해 설계된 통계적 가설 검정들이 존재한다. 사용할 수 있는 단위근검정은 다양하고, 서로 다른 가정에 기초하고 있으며, 상반되는 답을 낼 수도 있다. 분석 과정에서 KPSS(Kwiatkowski-Phillips-Schmidt-Shin)검정을 사용한다. 이 검정은 데이터에 정상성이 나타난다는 것이 귀무가설이고, 귀무 가설이 거짓이라는 증거를 찾으려고 한다. 결과적으로 ,작은 p-value는 차이를 구하는 것이 필요하다는 것을 나타낸다. 검정은 다음과 같다.
예로는 구글 주식 가격 데이터에 적용해보자.
library(urca)
goog %>% ur.kpss() %>% summary()
#>
#> #######################
#> # KPSS Unit Root Test #
#> #######################
#>
#> Test is of type: mu with 7 lags.
#>
#> Value of test-statistic is: 10.72
#>
#> Critical value for a significance level of:
#> 10pct 5pct 2.5pct 1pct
#> critical values 0.347 0.463 0.574 0.739
검정 통계량은 1% 임계값보다 훨씬 크다. 이것은 귀무가설이 기각 된다는 것을 의미한다. 즉, 데이터가 정상성을 가지고 있지 않다. 데이터에 차분을 수행할 수 있고 검정을 다시 적용할 수 있다.
goog %>% diff() %>% ur.kpss() %>% summary()
#>
#> #######################
#> # KPSS Unit Root Test #
#> #######################
#>
#> Test is of type: mu with 7 lags.
#>
#> Value of test-statistic is: 0.0324
#>
#> Critical value for a significance level of:
#> 10pct 5pct 2.5pct 1pct
#> critical values 0.347 0.463 0.574 0.739
이번에는 검정 통계가 작고(0.0324), 정상성이 나타나는 데이터에서 볼 수 있는 것처럼 범위 안에 잘 들어간다. 따라서 차분을 구한 데이터가 정상성을 나타낸다고 결론내릴 수 있다.
1차 차분의 적당한 횟수를 결정하기 위해 여러번의 KPSS 검정을 사용하는 이 과정을 ndiffs()로 수행할 수 있다.
ndiffs(goog)
#> [1] 1
위 검정들에서 본 것처럼, google데이터가 정상성을 나타내도록 하려면 한번의 차분(difference)이 필요하다.
계절성 차분(seasonal difference)이 필요한지 결정하기 위한 비슷한 함수는 nsdiffs()이다. 이 함수는 필요한 계절성 차분의 적당한 횟수를 결정하기 위해 계절성 강도 측정량을 사용한다. Fs<0.64이면, 계절성 차분이 필요 없다고 알려주고, 이외의 경우에는 하나의 계절성 차분이 필요하다 알려준다.
시계열(Time-Series) 데이터란 시간의 흐름에 따라 순차적으로(sequentially) 기록된 데이터를 가리킨다.
관측된 시계열 데이터를 분석하여 미래를 예측하는 문제가 바로 시계열 예측 문제이다.
시계열 예측 문제는 흔하게 접하는 문제로써 주로 경제 지표를 예측하거나, 어떤 상품의 수요를 예측하는 문제에 이르기까지 다양한 어플리케이션을 가지고 있다.
특히, 예측된 결과를 바탕으로 여러 정책이나 비즈니스 전략을 결정하는 과정에 활용되기에, 실제 비즈니스 영역에서는 시계열 예측 문제가 매우 중요하게 여겨지고 있다.
일례로 McKinsey Global Institute의 연구에 따르면, 시계열 데이터가 텍스트나 이미지 데이터 보다 더 큰 잠재적 가치를 가지고 있다고 보고있다.
데이터 유형에 따른 잠재적 인공지능 활용 가치(McKinsey Global Institute analysis)
시계열 예측 문제의 주요 챌린지
무엇보다 시계열 문제가 4가지 어려운 점이 있다.
1.예측값은 완벽할 수 없으며, 항상 변동의 가능성을 내포.
따라서 시계열 예측에 있어 불확실성(uncertainty)이 고려되어야 한다. 다시 말해 예측값에 대해 적합한 확률 분포 모델(probability distribution)이 고려되어야 한다. 단순히 어떤 에측 값을 제공한다기 보다는 그 예측값의 불확실성(또는 신뢰도)에 대해 정량적인 평가가 함께 고려되어야 예측 모델로써 더욱 가치가 있다.
2. 시계열 데이터에 숨겨진 여려 형태의 패턴을 찾아 학습하는 문제.
보통 trend,seasonality,cycle 등의 유형으로 패턴을 구분하곤 한다. 대게 시계열 데이터에는 이러한 패턴이 복잡하게 섞여있고 데이터의 길이가 불충분 하거나, 노이즈, 아웃라이어 데이터로 인해 손쉽게 패턴을 구분해 내어 찾기가 어렵다.
3. 다변량 시계열 데이터(multiple time-series)를 다뤄야 하는 경우가 많아지고 있다.
과거에는 주로 다변량의 시계열 데이터(unvariate time-series)의 분석과 예측 문제가 주로 다뤄졌지만 근래에는 다변량 시계열 데이터를 다뤄야 하는 경우가 많아지고 있다.
순차적으로 관찰된 수많은 변수들에 대해, 이 변수들간의 상관 관계를 학습할 수 있어야 한다. 가령, 어떤 상품의 수요가 다른 상품들의 수요 변화에 영향을 받거나, 어떤 지역의 택시 수요가 다른 지역의 택시 수요와 일정한 시간 차이를 두고 상관 관계를 보일 수 있다. 이 변수가 적게는 수백에서 많게는 수백만에 이를 수 있으므로, 이를 효율적으로 처리하고 학습할 수 있는 알고리즘이 매우 중요해 지고 있다.
4. 시계열 데이터에는 노이즈 데이터 또는 관찰되지 못한 기간이 종종 존재한다.
이는 측정하고 데이터를 기록하는 과정에서의 오류나 예측치 못한 상황으로 인해 발생할 수 있다. 예를 들어 상품의 품절로 인하여 장기간 판매량이 없는 경우, 해당 상품에 대한 실제 수요를 판매량으로 유추할 수 없는 경우이다.
시계열 예측 문제에 있어서는 이에 대한 적절한 전처리 과정이 매우 중요하다.
시계열 예측 모델의 평가
일반적으로 supervised machine learning 모델은 데이터 셋을 학습 데이터와 테스트 데이터로 랜덤하게 나누어 평가를 한다. 하지만 시계열 예측 문제는 평가와 학습 대상의 데이터를 특정 시간 기준으로 엄밀하게 분리시켜, data leak이 없이 평가가 이루어져야 한다.
가령 시계열 데이터의 중간 구간을 잘라내 테스트용으로, 그 이전과 이후를 학습용으로 사용하면 정확한 평가가 이루어질 수 없다. 학습 데이터는 반드시 테스트 데이터 보다 이전에 관찰된 것이어야 한다.
모델링 과정에서는 최적의 성능을 가지되, robust한 모델을 찾기 위해 cross-validation을 통해 성능을 평가한다. 제한된 데이터 셋으로부터 여러 쌍의 학습 및 테스트 데이터를 샘플링해서 모델의 성능이 기대치에 부합하며, 안정적인지를 확인하게 된다. 시계열 예측 모델에서는 앞서 언급한 chronological testing 기준을 지키며, cross-validation을 할 수 있는 방법으로 크게 두가지가 고려된다.
1. Sliding window
슬라이딩 윈도우(Sliding Window) 알고리즘은 배열이나 리스트의 요소의 일정 범위의 값을 비교할 때 사용하면 유용한 알고리즘이다.
예를 들어 정수로 이루어진 배열 [2,4,7,10,8,4,5,6,7,1] 에서 길이가 3인 서브배열의 합계가 가장 큰 서브배열은 무엇일까?
서브 배열 합계 체크
특정 길이의 하위 배열의 최대 합계값을 구하는 단순한 방식을 먼저 살펴보자.
단순하게 for 문으로 모든 배열 요소를 특정 길이만큼 순회하며 합계를 구해서 최대 값을 구하는 단순하고 비효율적인 방법이다.
import sys
def max_sub_array(arr, k):
maxsum = -sys.maxsize - 1
arraysize = len(arr)
for i in range(arraysize - k + 1):
current_sum = 0
for j in range(i, i + k):
current_sum += arr[j]
maxsum = max(maxsum, current_sum)
return maxsum
if __name__ == '__main__':
print(max_sub_array([2, 4, 7, 10, 8, 4, 5, 6, 7, 1], 3))
25 (7 + 10 + 8)
설명이 크게 필요 없는 단순한 코드이다.
물론 위와 같이 작성해도 데이터가 방대하지 않다면 큰 차이는 없을 것이다.
그렇다면 좀 더 효율적인 코드를 작성할 방법을 찾아보자.
우선 기존 단순한 방식에서 서브 배열의 요소를 순회하다 보면 중복되는 요소들이 존재한다.
서브배열의 공통된 요소 예시
위 이미지 처럼 범위가 5인 서브배열을 탐색하는 경우 0~4 범위의 서브배열과 1~5 범위의 서브배열은 공통적으로 1~4 범위가 중복된다.
중복되는 요소(공통 요소)들을 재상용하는 방법이 슬라이딩 윈도우 알고리즘의 핵심이다.
def max_sub_array(arr, k):
window_sum = 0
max_sum = 0
window_start = 0
for window_end in range(len(arr)):
window_sum += arr[window_end] # 슬라이딩 인덱스 범위 요소 합산
# 슬라이딩 윈도우의 범위가 k 보다 커진 경우
if window_end >= (k - 1):
max_sum = max(max_sum, window_sum)
window_sum -= arr[window_start] # 슬라이드 윈도우 범위를 벗어난 요소를 합계에서 제거
window_start += 1 # 슬라이드 윈도우 시작 인덱스 증가
return max_sum
if __name__ == '__main__':
print(max_sub_array([2, 4, 7, 10, 8, 4], 3))
25 (7 + 10 + 8)
이해를 쉽게 하기 위해 기존 배열보다 길이가 짧은 배열로 변경했다. 각 루프마다 다음과 같은 식으로 진행된다.
슬라이딩 알고리즘1슬라이딩 알고리즘2
우선 슬라이딩 윈도우를 관리하기 위한 변수들로 window_start, window_end, window_sum 등을 사용한다.
-window_start: 슬라이딩 윈도우 시작 인덱스
-window_end: 슬라이딩 윈도우 끝 인덱스
-window_sum: 슬라이딩 윈도우 합계
슬라이딩 윈도우는 하위 윈도우의 범위 k 와 동일하게 유지하며 범위 내의 요소는 모두 합산하고 범위 밖으로 벗어난 요소들은 빼준다. 다만 슬라이딩 윈도우의 범위가 k 보다 커질때까지 (window_end>=k-1)는 요소들을 합산하기만 하고 이후부터는 범위 밖의 요소들을 차감시켜준다.
이렇게 하면 기존처럼 매번 서브배열의 합계를 구할 필요 없이 모든 요소를 1번만 순회하면서 최대값을 구하게 되므로 빠르고 효율적으로 처리된다.
위 sliding window 알고리즘을 통해 고정된 사이즈 2개의 window를 움직여가며 학습,테스트 데이터를 추출하여 cross-validation을 하는 방법이다.
2. Expanding window
위 sliding window를 이해하였다면 이는 간단하게 이해할 수 있다. 데이터 길이가 충분치 않을때 데이터를 확장해 가며 학습,테스트 데이터를 추출하는 알고리즘이다.
Cross-validation for time-series forecasting model(Source: https://eng.uber.com/forecasting-introduction/)
시계열 예측 모델의 평가는 문제에 따라 다양한 지표가 활용될 수 있다. Root Mean Squared Error(RMSE) 또는 Mean Absolute Percentage Error(MAPE)가 많이 사용되는 지표이다. 일반적으로 baseline 예측 모델(e.g.최근 평균값으로 예측)과의 비교가 매우 유용하다. 머신러닝을 도입하기 이전의 다소 단순한 방식과의 비교를 통해 개발된 예측 모델의 도입 효과를 가늠해 보곤한다. 예측 모델의 평가는 특정 성능 지표 뿐만 아니라 residual error의 분포, 시간축에 따른 Error의 propagation 패턴을 분석하여 모델의 bias가 있는지 혹은 overfitting 여부 등에 대한 검토가 반드시 필요하다.
주요한 시계열 패턴
전통적으로 시계열 데이터를 분석할 때 보통 trend, seasonality,cycles 그리고 random noise의 패턴으로 분류하여 분석한다.
Trend 는 전체적으로 긴구간에 걸쳐 일정하게 증가 또는 감소하는 패턴을 가리킨다.
Seasonal component는 규칙적인 주기에 따라 반복되는 패턴을 나타내며, 주기가 규칙적이지 않지만 wave형태로 반복되는 패턴을 cyclical component로 분류하여 분석한다.
Trend,seasonality,cycles를 제외한 그외 불규칙한 fluctuation은 random component로 분류한다.
Decomposition of time-series data(Source: Forecasting: Principles and Practices)
일반적을 전통적인 시계열 분석 모델에서는 이 4가지 유형이 패턴이 linear한 조합으로 이루어진 모델(additive model)을 가정한다. 간혹 seasonal 또는 cyclic 패턴의 변동폭이 시계열 값의 크기에 따라서 함께 변하는 경우 multiplicative model을 가정하기도 한다. 분석 모델 안에서 패턴이 어떻게 구성되느냐에 따라 , 개별 component에 대한 수학적 모델 또한 매우 다양하다. 주어진 시계열 데이터에 맞는 모델을 찾아가는 과정이 예측 모델을 만들어 가는 과정이라고 볼 수 있다.
대표적인 전통 시계열 예측 모델: ARIMA
전통적으로 시계열 데이터에서 패턴을 추출하고 그 패턴을 이용하여 미래를 예측하는 모델로 Autoregressive integrated moving average(ARIMA)가 있다.
ARIMA모델은 크게 세가지 component로 조합되어 만들어진다. 먼저 Autoregressive는 예측 시점(t-0)이 과거 p 개의 데이터 포인트(t-1,t-2,...,t-p)와 연관이 있다는 의미이다. 일련의 과거 관측값을 이용하여 미래의 값을 예측하는 regression 모델인 것이다.
Moving average는 과거의 예측 에러를 예측 대상과 연관 시키는 component이다. 가령 autoregressive 모델로 미래를 예측함에 있어, 과거의 에러들을 고려하여 그 예측값을 보정하는 역할을 한다고 볼 수 있다. 마지막으로 integrated의 의미는 non-stationary time-series의 데이터를 differencing을 통해 시간에 따라 데이터의 통계적 특성(e.g. mean, std)이 변하지 않는 stationary time-series 데이터로 변환하는 모델을 의미한다. 적합한 ARIMA 예측 모델을 만들기 위해서는 differencing order을 조절하여 데이터로부터 non-stationary를 반드시 제거해야 한다.
ARIMA 모델은 추가적인 component를 고려함으로써 다양한 형태로 변형, 발전될 수 있다. 예를 들어 추가적인 seasonal component를 고려하는 SARIMA가 있으며, univariate time-series가 아닌 추가적인 covariates를 함께 고려하여 예측할 수 있는 ARIMAX 모델이 존재한다. 물론 이 두가지를 모두 혼합한 SARIMAX모델도 사용되곤 한다.
