[시계열] R에서 ARIMA 모델을 이용한 시계열 분석

시계열 데이터는 일정 기간 동안 일련의 시간 간격으로서 수집된 데이터 지점이다. 시계열 데이터 분석이란 이용 가능한 데이터를 분석하여 데이터의 패턴이나 추세를 파악하여 미래 가치를 예측함으로써 보다 효과적이고 비즈니스 의사결정을 최적화하는 데 도움이 되는 것을 의미한다.

 

시계열 분석 방법

 

 

 

게다가 , 타임 시리즈 분석은 다음과 같이 분류될 수 있다 :

• 1. 모수와 비모수
• 2. 선형과 비선형
• 3. 단변량과 다변량

타임시리즈 분석에 사용되는 기법들은 다음과 같다 : 

• 1. ARIMA 모델
• 2. Box-Jenkins 다변량 모델
• 3. Holt winters 지수 스무딩 (single, double and triple)

 

ARIMA 모델링

ARIMA는 AutoRegulation Integrated Moving Average의 약칭이다. 자기회귀(AR) 용어는 미분된 시리즈의 시차를 가리키고, 이동 평균(MA) 용어는 오류의 시차를 가리키며, I 는 시계열을 정상화(stationary) 하는데 사용되는 미분의 차수입니다.

ARIMA 모델의 가정 

• 1. Data 는 정상성(stationary)을 띄어야 합니다.  – 정상상태는 시리즈의 속성이 캡처된 시간에 따라 달라지지 않는다는 것을 의미한다. 백색 잡음(white noise) 시리즈와 주기적 동작이 있는 시리즈도 정상 시리즈로 간주할 수 있다. 
• 2. Data 는 단변량(univariate)이어야 합니다. – ARIMA 하나의 변수에 대해서만 작동합니다. 자기-회귀(Auto-regression)란 과거 값에 대한 모든 회귀식을 의미합니다. 

 

ARIMA 모델링을 하기 위한 단계는 다음과 같습니다: 

• 1. 탐색적 분석(Exploratory analysis)
• 2. 모델 적합 (Fit the model)
• 3. 진단적 조치 (Diagnostic measures)

 

R을 이용한 시계열 데이터 모델링의 첫번째 단계는 이용가능한 데이터를 시계열 데이터 형식으로 변환하는 것입니다. 그렇게 하기 위해서, 우리는 R에서 다음의 명령어를 수행할 필요가 있습니다. :

tsData = ts(RawData, start = c(2011,1), frequency = 12)

 

여기서, RawData는 우리가 시계열 시리즈로 바꿔야할 단변량(univariate) 데이터 입니다.  start에는 데이터의 시작점을 주도록 하자, 이 데이터의 경우에는 2011월 1월 이다. 그리고, 월별 데이터 이기 때문에 'frequency=12' 이다.  

 

실제 데이터 셋은 다음과 같이 생겼다. 

우리는 그래프 자체에서 데이터 포인트가 일부 특이치들과 함께 전반적인 상승 추세를 따른다는 것을 유추할 수 있다. 이제 우리는 데이터의 정확한 비-정상성, 계절성을 알아내기 위해 분석을 해야 한다.

탐색적 분석

• 1. 순차적인 의존성을 조사하기 위함 자기 상관 분석 : 과거의 어떤 값이 현재 값과 상관관계가 있는지 추정하는데 사용된다. ARIMA 모델에 대한 p,d,q 추정치를 제공한다. 
• 2. 주기적 동작을 조사하기 위한 스펙트럼 분석 : 시계열의 변동을 주기적 구성 요소로 설명할 수 있는 방법을 설명할 수 있는 방법을 묘사하기 위해 수행된다. 빈도 영역 분석(Frequency Domain Analysis)라고도 한다. 이를 통해 잡음이 심한 환경에서 주기적인 구성요소를 분리할 수 있다. 
• 3. 트렌드 추정과 분해 : 계절적 조정에 사용된다. 그것은 관찰된 시계열로부터 이들 각각이 특정한 특성을 갖는 다수의 구성 요소 시리즈(원래 시리즈를 재구성하는 데 사용될 수 있음)를 구성하려고 한다. 

데이터에 대해서 어떤 EDA를 수행하기 전에, 우리는 시계열 데이터의 세가지 구성요소를 이해할 필요가 있다: 

• 추세 (Trend): 데이터의 장기적인 증가 또는 감소를 추세라고 한다. 그것은 반드시 선형적인 것은 아니다. 그것은 시간의 경과에 따른 데이터의 기본 패턴이다. 
• 계절성 (Seasonal): 시리즈가 계절적 요인에 의해 영향을 받는 경우(즉, 분기, 월 또는 요일) 시리즈에는 계절성이 존재한다. 그것은 항상 고정되고 알려진 기간이다. 예: -크리스마스 기간 동안 매출이 급 상승하는 경우 등.
• 주기 (Cyclic): 일정 기간이 아닌 데이터가 상승하거나 하락할 때 우리는 이를 주기적 패턴이라고 부른다. 예를들어 - 이러한 변동 지속기간은 일반적으로 최소 2년이다. 

우리는 이 시계열 시리즈의 구성요소를 찾기 위해 다음의 R 코드를 사용할 수 있다 : 

 

components.ts = decompose(tsData)
plot(components.ts)

결과는 다음과 같다 :

우리는 4가지 구성 요소를 얻었다 :

• 관측값 ( Observed ) – 실제 데이터 플롯 
• 추세 ( Trend ) – 데이터 점들의 전반적인 상향 또는 하향 움직임
• 계절성 ( Seasonal )  – 데이터 점들의 월별/년별 패턴 
• 임의값 ( Random ) – 데이터의 설명할 수 없는 부분 

이 4개의 그래프를 자세히 관찰하면 데이터가 주로 정상성 및 계절성 등 ARIMA 모델링의 모든 가정을 만족하는지 확인할 수 있다. 

다음으로 ARIMA를 위한 비-정상적인 부분을 제거해야 한다. 여기서 논의를 위해 데이터의 계절적 부분도 제거하겠다. 계절 부분은 분석에서 제거하고 나중에 추가할 수도 있고, ARIMA 모델 자체에서 관리할 수도 있다.

 

정상성을 달성하기 위해서는 : 

• 데이터를 미분한다  – 연속 관측치 간의 차이를 계산한다. 
• 시계열 데이터를 로그 또는 제곱근을 취하여, 비-일정한 분산을 안정화 한다. 
• 데이터에 추세가 포함된 경우 데이터에 곡선 유형을 적합 시킨 다음 해당 적합치에서 잔차를 모형화 한다. 
• 단위 루트 테스트 - 이 테스트는 추세 데이터에 사용되어야 하는 일차 미분 또는 회귀 분석을 통해 데이터를 정상화 할 수 있다는 것을 발견하고자 사용된다. Kwiatkowski-Phillips-Schmidt-Shin (KPSS) 테스트에서 작은 p-값은 미분이 필요하다는 것을 암시한다. 

다음은 단위 루트 테스트를 위한 R 코드 이다 : 

 

library("fUnitRoots")
urkpssTest(tsData, type = c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)
tsstationary = diff(tsData, differences=1)
plot(tsstationary)

결과는 다음과 같다 : 

 

비- 정상성을 제거한 후에는 다음과 같다 : 

다양한 플롯들과 함수들이 계절성을 탐지하는데 도움을 준다 : 

• 계절 하위 시계열 플롯 
• 다중 상자 플롯 
• 자기 상관 플롯 
• ndiffs() 는 시계열을 비-계절적으로 만드는데 필요한 첫번째 미분의 수를 결정하는데 사용된다. 

자기 상관을 계산하기 위한 R 코드는 다음과 같다 : 

 

acf(tsData,lag.max=34) 

자기 상관 함수(acf())는 가능한 모든 시차에서 자기 상관을 제공한다. 지연 0에서의 자기 상관은 기본적으로 포함되며, 이 값은 항상 데이터와 자신 사이의 상관 관계를 나타내는 값 1을 차지한다. 위의 그래프에서 유추할 수 있듯이, 자기 상관성은 지연이 증가함에 따라 계속 감소하여 더 큰 시차로 구분된 관측치 사이에 선형적인 연관성이 없음을 확인시켜 준다. 

 

timeseriesseasonallyadjusted <- tsData- timeseriescomponents$seasonal
tsstationary <- diff(timeseriesseasonallyadjusted, differences=1)

데이터에서 계절성을 제거하기 위해 원래 시계열에서 계절 성분을 뺀 다음 이를 변화시켜 정지 상태로 만든다.

 

계절성을 제거하고 데이터를 정지시킨 후 다음과 같이 보이게 된다 : 

스무딩은 보통 우리가 시계열의 패턴과 추세를 더 잘 볼 수 있도록 돕기 위해 행해진다. 일반적으로 불규칙한 거칠기를 부드럽게 하여 보다 선명한 신호를 본다. 계절 데이터의 경우 추세를 파악할 수 있도록 계절성을 완화시킬 수 있다. 스무딩은 우리에게 모델을 제공하지 않지만, 시리즈의 다양한 성분을 설명하는 좋은 첫 번째 단계가 될 수 있다.

시계열을 평활하려면:

• 보통 이동 평균(단일 이동, 중심) - 각 시점에서 특정 시간 이전의 관측 값의 평균을 결정한다. 시계열에서 계절성을 제거하여 추세를 더 잘 볼 수 있또록 길이 = 계절적 범위 인 이동 평균을 사용한다. 계절적 범위는 계절성이 반복되는 기간이다. 예를들어 계절성이 매년 12월 이면 - 12개월이다. 따라서 평활 시리즈에서 각 평활 값은 전체 계절에 걸쳐 평균을 구했다. 
• 지수 가중 평균 – 각 시점에서 기하급수적으로 감소하는 가중 요소를 적용한다. 각 오래된 기준점에 대한 가중치는 기하급수적으로 감소하고 영(0)에 도달하지 않는다.

 

모델에 적합하기

데이터가 준비되고 모델링의 모든 가정을 만족하면, 데이터에 적합할 모델의 순서를 결정하기 위해, 우리는 모델의 자기 회귀, 통합 및 이동 평균 부분의 순서를 각각 참조하는 음이 아닌 정수인 p, d, q의 세 변수가 필요하다.

어떤 p와 q 값이 적절할지를 조사하기 위해서는 acf()와 pacf() 함수를 실행해야 한다.

지연 k에서의 pacff는 자기 상관 함수로, 정확히 k step 떨어져 있는 모든 데이터 지점 간의 상관 관계를 설명한 후, 해당 k 단계 간의 데이터와의 상관 관계를 설명한다. ARIMA 모델에서 자기 회귀(AR) 계수(p-값)의 수를 식별하는 데 도움이 된다.

다음은 acf() 및 pacf() 명령을 실행하는 R 코드이다.

 

acf(tsstationary, lag.max=34)
pacf(tsstationary, lag.max=34)

플롯은 다음과 같다. 

 

p 및 q 값을 정의하는 acf()의 모양:
그래프를 보고 표를 통해 어떤 모델을 선택할지, 그리고 p, d, q의 값을 결정할 수 있다.

fitARIMA <- arima(tsData, order=c(1,1,1),seasonal = list(order = c(1,0,0), period = 12),method="ML")
library(lmtest)
coeftest(fitARIMA)

순서는 ARIMA 모델의 비계절 부분을 지정한다: (p, d, q)는 AR 순서, 차이 정도 및 MA 순서를 가리킨다.

ARIMA 모델의 계절 부분과 주기(이 경우 기본적으로 빈도(x) 즉 12로 설정됨)를 지정한다. 이 함수는 성분 순서와 주기가 있는 리스트를 요구하지만, 길이 3의 숫자 벡터를 부여하면 사양을 '차수'로 하는 적합한 리스트로 변한다.

메소드는 '최대우도(ML)' 또는 '조건부 제곱합(CSS)'을 최소화할 수 있는 적합 방법을 가리킨다. 기본값은 조건부 제곱합이다.

 

 

이것은 재귀적 과정이며 우리는 가장 최적화되고 효율적인 모델을 찾기 위해 다른 (p,d,q) 값으로 arima() 함수를 실행해야 한다.

fitarima()의 출력은 적합 계수와 각 계수에 대한 표준 오차(s.e)를 포함한다. 계수를 관찰하면 중요하지 않은 계수는 제외할 수 있다. confint() 함수를 이 목적에 사용할 수 있다.

confint(fitARIMA)

 

최적 모델을 선택

R은 ARIMA 모형을 추정하기 위해 최대우도 추정(MLE)을 사용한다. 모수 추정치를 찾을 때 p, d, q의 주어진 값에 대한 로그 우도를 최대화하여 우리가 관찰한 데이터를 얻을 확률을 최대화하려고 한다.

모델 세트에 대한 Akaike의 정보 기준(AIC)을 알아보고 AIC 값이 가장 낮은 모델을 조사하자. 또 슈바르츠 베이지안 정보 기준(BIC)을 사용해 BIC 값이 가장 낮은 모델을 조사해보자. 최대우도 추정을 사용하여 모형 모수를 추정할 때 모수를 추가함으로써 가능성을 높일 수 있으며, 이는 과다 적합을 초래할 수 있다. BIC는 모형의 파라미터 수에 대한 페널티 항을 도입하여 이 문제를 해결한다. AIC, BIC와 함께 우리는 또한 그러한 계수값들을 면밀히 관찰할 필요가 있고 우리는 그 성분들의 유의 수준에 따라 포함시킬지 여부를 결정해야 한다.

진단적 조치 

잔차의 ACF를 표시하고 portmanteau 테스트를 수행하여 선택한 모델의 잔차에서 패턴을 찾아 보자. 플롯이 백색 소음처럼 보이지 않으면 우리는 수정된 모델을 시도해 볼 필요가 있다.

잔차가 백색 노이즈처럼 보이면 예측값을 계산해보자.

Box-Ljung test

이것은 하나의 구체화된 모든 시차들의 독립성을 검증하는 테스트이다. 각각의 다른 시차들의 임의성을 테스트 하기보다, 시차의 수에 기반한 "전체적인" 무작위성을 테스트하는 것이다. 그래서 portmanteau 테스트라고 한다. 이는 적합된 ARIMA 모델의 원래 시리즈(series)가 아닌 나머지(residuals)에 적용된다. 이러한 적용에서 실제로 테스트되고 있는 가설은 ARIMA 모델의 잔차가 자기 상관 관계가 없다는 것이다. 

 

box 테스트 결과를 얻기 위한 R code 이다. : 

 

acf(fitARIMA$residuals)
library(FitAR)
boxresult-LjungBoxTest (fitARIMA$residuals,k=2,StartLag=1)
plot(boxresult[,3],main= "Ljung-Box Q Test", ylab= "P-values", xlab= "Lag")
qqnorm(fitARIMA$residuals)
qqline(fitARIMA$residuals)

Output:

잔차의 ACF는 중요한 자기상관이 없음을 보여준다. 

Ljung-Box Q 테스트를 위한 p-values는 0.05를 훨씬 상회하며, "비-중요함"을 암시한다. 

그 값들은 한 줄로 서서 사방에 널려 있지 않기 때문에 정상이다. 

 

우리의 가설을 지지하기 위한 모든 그래프들이 잔차에 대한 패턴을 보이지 않기 때문에, 우리는 다음으로 넘어가서 예측값을 계산할 수 있다. 

 

작업 흐름도 (Work flow diagram)

auto.arima() 함수:

forecast package는 두가지 함수를 제공한다: 지수적 그리고 ARIMA 모델의 자동 선택을 위한 ets() 와 auto.arima() 이다.  
R에서의 auto.arima() 함수는 AIC의 축약 버전인 단위 루트 테스트의 조합과 ARIMA 모델을 얻기 위한 MLE를 사용한다.

KPSS 테스트는 자동 ARIMA 모델링을 위한 Hyndman-Khandakar 알고리즘에서 미분의 수 (d)를 결정하기 위해 사용된다. 

p,d, 와 q는 AICc의를 최소화하면서 선택된다. 알고리즘은 단계적 검색을 사용하여 모형 공간을 가로지르며 가장 작은 AICc로 최적의 모형을 선택한다.

만약 d = 0 이라면, 상수 c는 포함이 된다; 만약 d≥1 이라면 상수 c는 0으로 설정된다. 현재 모델의 분산은 현재 모델에서 ±1로 p 및/또는 q를 변화시키고 현재 모델에서 c를 포함/제외함으로써 고려된다.

지금까지 고려된 최고의 모델(현재 모델 또는 이러한 변형 중 하나)이 새로운 현재 모델이 된다.
이제 이 과정은 더 낮은 AIC를 찾을 수 없을 때까지 반복된다.

 

 

auto.arima(tsData, trace=TRUE)

 

ARIMA 모델을 사용하여 예측하기

해당 ARIMA 모델의 매개변수는 시계열 데이터에 가장 적합한 모델을 선택한 후 시계열의 미래 값에 대한 예측을 위한 예측 모델로 사용할 수 있다.

d-값은 예측 구간에 영향을 준다 - 'd'값이 높아질 수록 예측 구간은 커진다. 장기 예측의 표준 편차가 역사적 데이터의 표준 편차를 향해 가기 때문에, d = 0 일 때, 예측 구간은 필수적으로 같다. 

다양한 모델 피팅 함수의 결과로부터 예측하는 데 사용되는 predict()라는 함수가 있다. 이는 예측하기 위해 얼마나 많은 시간 단계를 앞지르는지 명시하는 인자 n.ahead()가 필요로 한다. 

 

predict(fitARIMA,n.ahead = 5)

forecast R 패키지의 forecast.Arima() 함수는 또한 타임 시리즈의 미래 값을 예측하기 위해 사용될 수 있다. 

우리는 또한 level 인자를 사용함으로써 예측구간을 위한 신뢰 수준을 구체화할 수 있다.  

 

futurVal <- forecast.Arima(fitARIMA,h=10, level=c(99.5))
plot.forecast(futurVal)

예측 오차는 보통 평균 0과 일정한 분산을 가지고 분포하며 상관성이 없는지 확인할 필요가 있다. 진단 측정값을 사용하여 예측값이 가장 높은 적절한 모델을 찾을 수 있다.

 

 

 

 

예측은 파란색 선으로 표시되며 80% 예측 구간은 어두운 음영 영역으로, 95% 예측 구간은 밝은 음영 영역으로 표시된다.
이것은 ARIMA를 이용하여 시계열 데이터를 분석하고 기존 시계열의 값을 예측할 수 있는 전체적인 과정이다.

 

 

Click here to get the entire code.

 

References

Stationarity and differencing
Time Series and Forecasting
stats
sdstate
General seasonal ARIMA models

 

 

 

 

출처

https://datascienceplus.com/time-series-analysis-using-arima-model-in-r/

Time Series Analysis Using ARIMA Model In R