2.1.1 ARIMA 모형

ARIMA(Auto-Regressive Integrated Moving Average) 모형은 현재 시계열 값을 예측하기 위해 과거의 관측값, 오차 그리고 차분 절차를 활용한다. 이 모형은 자기회귀(AR) 모델과 이동평균(MA) 모델의 특성들을 결합하여 구성되어 있다. 시계열 $Z_{t}$가 자기회귀 모형의 차수 p, 이동평균 모형의 차수 q, 그리고 차분 차수 d를 갖는 $ARIMA(p,\: d,\: q)$ 과정을 따르며 평균 $\mu$를 가진다면, 해당 시계열은 Eq. (1)의 형태로 표현된다^{(Box et al., 1970)}.

where, $W_{t}=(1-B)^{d}(Z_{t}-\mu)$,

$\phi(B)=1-\phi_{1}(B)-\cdots -\phi_{p}(B^{p})$,

$\theta(B)=1-\theta_{1}(B)-\cdots -\theta_{q}(B^{q})$,

$\epsilon_{t}\sim iid N(0,\: \sigma^{2})$

여기서, p는 자기회귀 모형의 차수로, t 시점의 값이 (t - 1), …, (t - p) 시점의 값에 영향을 주는 것을 의미한다. d는 차분을 뜻하고, t 시점의 값에 (t - d) 시점의 값을 빼서 정상 시계열 형태로 만들어준다. q는 이동평균 모형의 차수로, t 시점의 값에 연속적인 (t - 1), …, (t - p) 시점의 오차들이 영향을 주는 것을 의미한다. 그리고, $\phi(B),\: \theta(B)$는 각각 자기회귀, 이동평균 모형에 대한 다항식을 나타낸다. B는 후진연산자를 의미하며, d는 시계열 $Z_{t}$의 차분 차수를 나타낸다. ARIMA 모형은 비정상적(Nonstationary) 시계열 자료를 분석하는 데 사용되는 방법이다.

실제 ARIMA 시계열 분석은 적용공분산 정상성(Covariance Stationary) 적합 과정을 거친 후 예측을 진행하게 된다. 단변량 ARIMA 분석기법은 두 가지의 특징을 가진다. 첫째, 하나의 시계열 자료만으로 변동상태를 확인할 수 있다. 둘째, 시계열 자료 중 시간별 자료의 변동이 클 때 변동을 반영하여 예측을 진행할 수 있다는 장점을 가지고 있다.

ARIMA 모델을 시계열 자료에 적합하는 과정은 시계열 자료의 정상성(Stationarity) 여부를 검정한다. 정상성은 시계열 데이터에서 평균, 표준편차와 같은 확률적 특성들이 시간 경과에 따라 변하지 않는 상태를 의미한다. 비정상적일 때는 차분을 이용해 정상 시계열 자료로 만들어주는 과정이 선행되어야 한다. 이 과정에서 자기상관함수(Autocorrelation Function)와 부분자기상관함수(Partial Autocorrelation Function) 등의 그림을 활용할 수 있다. 위 과정 이후 잠정 모델의 AIC(Corrected Akaike’s information criterion)값을 산출하여 AIC가 최소가 되는 모형의 하이퍼 파라미터 p, d, q를 결정한다. 모델의 적합성을 평가하기 위한 통계적 진단 방법으로는 잔차의 자기상관함수(ACF)와 부분자기상관함수(PACF)를 계산하여 잔차가 독립적인지를 확인하는 잔차 분석 방법을 주로 사용한다. 최종적으로, 적합한 ARIMA 모델을 활용하여 미래의 시계열 데이터 값을 예측한다^{(Bari et al., 2015)}.

2.1.2 SARIMA 모형

계절 ARIMA 모형은 순수 계절적 특성을 가진 계절시계열이나 계절적 특성 없이 비정상성만을 지닌 시계열과는 다른 특성을 가진 시계열 자료이다. 그러므로, 계절성과 비정상성을 함께 고려하는 것이 필수적이다. 계절적 영향을 제거한 후, 분산을 안정화하기 위해 데이터 변환과 추세를 제거하기 위한 차분을 적용함으로써, 정상 시계열 데이터를 확보할 수 있다. 해당 모형을 예측 적합의 과정을 통하여 예측 결과를 도출할 수 있다. SARIMA(Seasonal ARIMA) 모형은 ARIMA 모형에 계절적 특성과 주기성을 추가한 모형이며, 시계열 자료의 구성성분들이 확률성을 가지거나 다른 구성요소에 종속적인 경우에 사용할 수 있다. 시계열 $Z_{t}$가 계절적인 주기 S를 가진다면 Eq. (2)으로 표현할 수 있다^{(Balibey et al., 2015; Bari et al., 2015; Naz, 2015; YoosefDoost et al., 2017)}.

where, $\phi_{p}(B)=1-\phi_{1}(B)-\cdots -\phi_{p}(B^{p})$,

$\theta(B)=1-\theta_{1}(B)-\cdots -\theta_{q}(B^{q})$,

$\Phi_{P}(B^{S})=1-\Phi_{1}(B)-\cdots -\Phi_{P}(B^{P})$,

$\Theta_{Q}(B^{S})=1-\Theta_{1}(B)-\cdots -\Theta_{Q}(B^{Q})$,

$\epsilon_{t}\sim iid N(0,\: \sigma^{2})$

여기서, B는 후진연산자이고, $\phi_{p}(B)$와 $\theta_{q}(B)$는 p, q 차수의 비 계절적 자기회귀, 이동평균 모형에 대한 다항식이며, $\Phi_{P}(B^{S})$와 $\Theta_{Q}(B^{S})$는 각각 P, Q 차수를 갖는 계절적 자기회귀, 이동평균 연산자, d는 비 계절적 차분 차수를 D는 계절적 차분 차수를 의미한다.

비정상적 시계열 자료에 계절효과가 추가되면 계절적 특성이 존재하는 비정상적 시계열이기 때문에 Eq. (3)과 같은 계절 ARIMA 모형을 이용한다. 계절 ARIMA 모형은 비 계절적 요소인 p, d, q와 계절적 요소인 P, D, Q로 구성되어 있다. 여기서 p는 시계열의 현재 값이 이전의 ‘p’개의 값의 어떻게 의존하는지를 의미한다. d는 비정상적인 시계열을 정상적인 시계열로 변환하기 필요한 차분의 횟수이다. q는 시계열의 현재 오차가 이전의 ‘q’개의 오차에 어떻게 의존하는지를 나타낸다. 계절적 요소인 P는 계절적 주기 내에서 시계열의 현재 값이 이전의 ‘P’개의 계절적 값에 어떻게 의존하는지를 나타낸다. D는 계절적 비정상성을 제거하기 위해 필요한 차분의 횟수이다. Q는 계절적 주기 내에서 시계열의 현재 오차가 이전의 ‘Q’개의 계절적 오차에 어떻게 의존하는지를 나타낸다.

계절 ARIMA 모형은 순수 계절적 특성을 가진 계절시계열이나 계절적 특성 없이 비정상성만을 지닌 시계열과는 다른 특성을 가진 시계열 자료이다. 그러므로, 계절성과 비정상성을 함께 고려하는 것이 필수적이다. 계절적 영향을 제거한 후, 분산을 안정화하기 위해 데이터 변환과 추세를 제거하기 위한 차분을 적용함으로써, 정상 시계열 데이터를 확보할 수 있다. 해당 모형을 예측 적합의 과정을 통하여 예측 결과를 도출할 수 있다.

SARIMA 모형 적합 과정은 정상성검정, 추정, 적합모형 선정으로 이루어진다. 첫째, 정상성 검정에서는 시계열 데이터의 평균, 분산, 자기상관함수가 시간에 따라 일정한지를 확인한다. 이 과정에서 시계열이 비정상적으로 나타날 경우, 변수 변환(Variable Transformation)과 차분(Difference)을 이용해 시계열 패턴을 정상적으로 만든다. 정상성 검정을 위해서는 그래프 분석 방법이 사용된다. 이러한 분석을 통해 차분의 차수인 ‘d’와 계절적 차분의 차수 ‘D’가 결정된다. 다음 단계에서는 ACF와 PACF 그래프, AIC 통계량 등을 활용하여 비계절적 차수 p, q와 계절적 차수 P, Q를 선택하여 초기 모델을 설정한다. 최종적으로 적합 모형 선정 과정에서는 통계적 검증을 통해 만들어진 모형이 적합한지를 판단하는 것으로 잔차 분석과 과대 적합 진단을 이용한다. 만약 모델이 통계적으로 유의미하지 않은 경우, 다시 이전 단계를 수행하며 다른 모델을 선택한다. 이후, 선택된 모델에 대한 계수 추정과 적합성 검정을 반복적으로 수행한다. 위 과정을 거쳐 최종적으로 확정된 모델을 사용하여 시계열 데이터의 예측값을 도출한다. 위 과정은 모델이 데이터에 적합하고 통계적으로 유의미한 결과를 제공할 때까지 반복 진행한다.

2.2.1 LSTM(Long Short-Term Memory)

일반적인 RNN(Recurrent Neuron Network) 구조는 Fig. 1과 같으며, 좌측의 RNN 구조를 시간 순서대로 나열하면 과거의 정보가 임의 시점인 t에 영향을 주는 루프형태를　띄고 있으며, 데이터의 가변길이인 시퀀스(Sequence)를 시간의 흐름과 연결 지을 수 있기에 시계열 데이터의 예측에 유용하다. 다만 긴 시퀀스의 자료를 사용한 경우 깊게 쌓여진 심층 RNN이 만들어지고, 심층 RNN의 경우 학습이 진행되다가 오차 경사의 기울기가 사라지고 입력층과 인접한 층에서 데이터의 정보를 기억하지 못해 예측에 이용할 수 없다.

Fig. 1. Concept of Recurrent Neural Network Structure(Lukic, 2020)

Fig. 2. Structure of LSTM(Yang et al., 2020)

RNN의 한 종류인 LSTM(Long Short-Term Memory)은 기존의 장기 의존성 문제(Long-Term Depedencies)를 해결하기 위해서 나온 모형이다. 따라서 임의 시점의 데이터뿐만 아니라, 과거 데이터와 임의 시점간의 유사성을 고려하여 미래 데이터를 예측하기 위한 모형이다. LSTM은 RNN의 문제점인 오차 경사의 기울기 소실(Gradient Vanishing) 및 최적화 오류(Optimization Hurdle)를 해결할 뿐 아니라 시퀀스를 늘린 장기적인 시간의 종속성을 감지 하는데도 효과적이다.

LSTM 구조는 Fig. 2와 같으며, 핵심요소로는 시간에 따라 상태(State)를 유지할 수 있는 메모리이동 셀(Cell)과 셀 안팎으로 데이터의 흐름을 조절하는 세 개의 비선형 게이트 구조이다. 즉, LSTM은 특정 시점의 상태($h_{t}$)를 업데이트하기 위해 셀(Cell, $C_{t}$)이라는 개념을 도입하여 입력과 현재까지의 상태를 이용하여 내부에 가지고 있는 정보를 업데이트할 것인지 아닌지를 판단한다^{(Olah, 2015)}. 셀의 들어가는 자료를 조절하기 위한 게이트의 종류로는 입력 게이트(Input Gate, $i_{t}$), 망각 게이트(Forget Gate, $f_{t}$) 그리고 출력 게이트(Output Gate, $o_{t}$) 크게 3가지로 이루어져 있다.

망각게이트는 Eq. (4)로 나타낼 수 있으며 이전 셀의 출력인 $h_{t-1}$과 현재의 입력인 $x_{t}$를 시그모이드(Sigmoid Function) 활성레이어에 적용해 0과 1사이의 값을 얻어 이 값을 현재 상태와 요소끼리 곱하게 되며, 그 과정에서 이 정보를 유지할지 제거할지를 선택하게 된다^{(Jung et al., 2018)}.

위 식에서 $\sigma$는 활성화 함수를 의미하며 $W_{f}$는 망각게이트의 가중치를 $b_{f}$는 망각게이트의 편향을 의미한다. 입력게이트 $i_{t}$에서는 Eq. (5)와 같이 셀에 어떤 정보를 선택해 저장할지를 정하게 되며 해당 과정은 크게 두 단계로 나누어진다. 첫째, 활성화 함수를 이용해 업데이트 여부를 결정하고, 둘째, 하이퍼볼릭 탄젠트(Hyperbolic Tangent, tanh) 함수를 통해 새로운 셀 상태로 업데이트하는데 사용 될 후보 셀($\widetilde{C}_{t}$)을 Eq. (6)을 통해 생성하게 된다.

여기서, $W_{i}$는 입력게이트의 가중치를 $W_{c}$는 후보 셀의 가중치를 의미한다. $b_{i}$는 입력게이트의 편향을 $b_{c}$는 후보 셀의 편향을 의미한다. 임의시점 이전의 셀 상태($C_{t-1}$)와 임의시점의 생성된 후보 셀($\widetilde{C}_{t}$)을 Eq. (7)과 같이 조합하여 임의시점의 셀 상태($C_{t}$)를 업데이트하게 된다.

마지막으로 출력게이트 $o_{t}$는 Eq. (8)은 활성화 함수인 시그모이드 함수를 이용해서 셀 상태 중 중요도에 따라 출력 여부를 결정한다. 마지막으로 Eq. (9)와 같이 하이퍼볼릭 탄젠트 함수를 사용하여 활성화된 임의 시점의 셀 상태($C_{t}$)와의 곱을 통해 특정 시점의 상태($h_{t}$)를 업데이트하게 된다.

여기서, $W_{o}$는 출력게이트의 가중치, $b_{o}$는 편향을 나타낸다.

3.2.1 ARIMA 모형 월 강수량 예측

모형의 적합성을 높이기 위해 최소의 AIC 값을 얻을 수 있는 p, d, q를 0~2의 범위에서 시행 오차법을 통해 선정 후 예측을 진행하였다^{(Hyndman et al., 2008)}. 월 강수량이 음수로 예측되면, 강수량의 특성을 고려하여 0으로 변환하였고 시행 오차법으로 선정된 ARIMA 모형의 파라미터인 p, d, q를 Fig. 6에 박스플롯과 커널 밀도 추정의 특성을 합한 바이올린 플롯을 이용하여 도시하였으며, 바이올린 플롯을 사용하여 지역별 적합한 모형의 매개변수의 빈도를 확인하였다. Fig. 7에서는 ARIMA 모형을 활용하여 2013년부터 2022년까지 예측을 진행하였다.

평가지표인 RMSE와 MAE가 가장 작게 나온 대구의 RMSE는 64.39 mm, MAE는 47.81 mm이고, 가장 크게 나온 서울의 경우 RMSE는 106.97 mm, MAE는 72.01 mm이다. Fig. 7에서는 특정 기간을 제외하고 주기함수의 형태로 예측이 진행되었음을 알 수 있다.

Fig. 6. ARIMA Model Hyperparameter(p, d, q). (a) Gangneung, (b) Gwangju, (c) Daegu, (d) Daejeon, (e) Busan, (f) Seoul, (g) Jeju, (h) Chuncheon

Fig. 7. Predicted Monthly Precipitation Using ARIMA(2013~2022). (a) Gangneung, (b) Gwangju, (c) Daegu, (d) Daejeon, (e) Busan, (f) Seoul, (g) Jeju, (h) Chuncheon

3.2.2 SARIMA 모형 월 강수량 예측

SARIMA 모형의 하이퍼 파라미터의 경우 유사 연구^{(Dimri et al., 2020)}에서 월별 시계열 자료를 예측할 때 주기성(S)를 12로 설정하였고 ARIMA 모형의 p, d, q와 계절적 특성을 고려한 P, D, Q가 합쳐진 SARIMA 모형을 활용해 월 강수량 예측을 진행한 선행 연구들을 참고하여 월별 강수량 예측 모형에 적합한 값으로 선정되는 0~2를 대입하여 시행 오차법을 통해 최소의 AIC값을 가지는 하이퍼 파라미터를 산정하였다^{(Aliyu et al., 2021; Dabral et al., 2017)}. 매개변수 산정 결과는 Fig. 8에 박스플롯과 커널 밀도 추정의 특성을 합한 바이올린 플롯을 이용하여 도시하였다. Fig. 9에서는 SARIMA 모형을 활용해 2013년부터 2022년까지 월 강수량 예측 결과를 도시하였다.

ARIMA(Fig. 7) 서울과 SARIMA(Fig. 9) 서울을 비교해보면 주기함수의 형태가 보이는 것은 유사하나 강수 최대값을 예측하는 것이 ARIMA(Fig. 7) 서울에서는 약 250 mm 정도이지만, SARIMA (Fig. 9) 서울은 약 700 mm로 ARIMA 모형에 비해 강수 최대값을 예측하기에 적합한 모형임을 알 수 있다.

Fig. 8. SARIMA Model Hyperparameter(p, d, q, P, D, Q). (a) Gangneung, (b) Gwangju, (c) Daegu, (d) Daejeon, (e) Busan, (f) Seoul, (g) Jeju, (h) Chuncheon

Fig. 9. Predicted Monthly Precipitation Using SARIMA(2013~2022). (a) Gangneung, (b) Gwangju, (c) Daegu, (d) Daejeon, (e) Busan, (f) Seoul, (g) Jeju, (h) Chuncheon

3.2.3 LSTM 모형 월 강수량 예측

LSTM 모형은 Table 1과 같이 구성하였다. Learning rate는 0.001과 0.01를 Batch size는 32, 64, 128을 Epoch는 50, 100, 200으로 모형을 구성하였을 때, 시행오착법을 통해 RMSE를 최소로 하는 하이퍼 파라미터를 선정하였다. 연도별로 선정한 하이퍼 파라미터를 Learning rate, Batch size와 Epoch 순으로 Fig. 10에 도시하였으며 모형 구성 이후 MSE(Mean Squared Error)를 손실함수로 사용해 월 강수량 예측을 진행하였다. Fig. 10은 적합한 하이퍼 파라미터를 색상으로 구분하여 해당연도별 및 지역별로 적합 모형의 하이퍼 파라미터를 직관적으로 도시하였다. LSTM 모형을 활용해 2013년부터 2022년까지 월 강수량을 예측한 결과를 Fig. 11에 도시하였다.

통계적 모형(ARIMA, SARIMA)의 예측 결과와 LSTM 모형의 예측 결과를 비교해보면 통계적 모형의 경우 예측에서 주기함수의 형태를 보였으나 LSTM 모형의 예측 결과 중 광주, 서울, 춘천을 보면 단순 주기함수에서 벗어나 강수 최대값이 발생하는 월을 예측하여 강수 최대값에 대한 예측뿐 아니라 건기 구간에서도 적절한 예측이 이루어짐을 살펴볼 수 있다.

Fig. 10. Optimal Learning Rate, Batch Sizes, and Epochs for LSTM Models by Stations and Year

Fig. 11. Predicted Monthly Precipitation Using LSTM(2013~2022). (a) Gangneung, (b) Gwangju, (c) Daegu, (d) Daejeon, (e) Busan, (f) Seoul, (g) Jeju, (h) Chuncheon

3.2.4 GBM 모형 월 강수량 예측

GBM 모형의 경우 하이퍼 파라미터를 Table 2와 같이 Learning rate는 0.001, 0.01, 0.05을 Number of estimators 50, 100, 200을 넣어 가장 작은 RMSE를 갖는 하이퍼 파라미터를 시행오차법을 이용해 모형을 적합 후 예측을 진행하였다. Fig. 12에 연도별로 적합한 모형의 하이퍼 파라미터는 색상별로 구분하여 사용되는 하이퍼 파라미터를 직관적으로 지역 및 연도별로 구분하였다.

GBM 모형을 이용한 월 강수량 예측 결과는 Fig. 13에 도시하였다. 산정 결과를 분석해보니, GBM 모형은 최대값을 예측하는 데 어려움이 있었고, 매월 최소 80 mm 이상의 강수가 내릴 것으로 예측되는 특성을 파악할 수 있다.

Fig. 12. Optimal Learning Rate and n_estimators for GBM　Models by Stations and Year

Fig. 13. Predicted Monthly Precipitation Using GBM(2013~2022). (a) Gangneung, (b) Gwangju, (c) Daegu, (d) Daejeon, (e) Busan, (f) Seoul, (g) Jeju, (h) Chuncheon