A_CORRELATE 함수를 이용한 자기상관(Autocorrelation) 기법의 활용 [2]

지난 회 게시물에 이어서 오늘은 갯수가 좀 더 많은 데이터에 대하여 A_CORRELATE 함수를 적용하는 예제를 다뤄보겠습니다. 다음과 같이 51개의 값들로 구성된 가상의 1차원 배열을 생성합니다. 여기서는 SIN, RANDOMN 함수를 적절히 이용하여 어느 정도 주기성을 보이는 가상 데이터를 만들어보았습니다.

x = FINDGEN(51)

y = SIN(x*20*!DTOR)*(1+RANDOMN(1, N_ELEMENTS(x))*0.2)+5

nn = N_ELEMENTS(y)

이제 이 데이터에 대하여 A_CORRELATE 함수를 적용하여 자기상관계수를 산출해보고, 그 결과를 그림으로 표출해봅시다. 지난 회 게시물에서도 언급했듯이, A_CORRELATE 함수를 사용하려면 대상 1차원 배열 및 Lag 값을 인자로 투입해야 합니다. 일단 먼저 Lag의 값을 2로 가정하고 자기상관계수를 산출하여 출력하는 과정을 다음과 같이 작성해 봅시다.

lag = 2

result = A_CORRELATE(y, lag)

PRINT, 'AC Coefficient =', result

이 내용을 실행하면 결과값은 약 0.71 정도로 나오게 됩니다. 물론 lag의 값을 바꾸면 그에 따른 자기상관계수 값을 산출할 수 있습니다.그러면 여기서 끝내지말고 그림을 좀 그려보기로 합시다. 이를 위하여 다음과 같은 예제 코드를 작성해보았습니다.

str = 'AC Coefficient = '+STRING(result, FORMAT='(F0.2)')

y_shifted = SHIFT(y, -lag)

win = WINDOW(DIMENSIONS=[600, 500], /NO_TOOLBAR)

pl = PLOT(x, y, YRANGE=[2, 10], /NODATA, MARGIN=0.1, $

  XTITLE='X', YTITLE='Y', FONT_SIZE=12, TITLE=str, /CURRENT)

plo1 = PLOT(x, y, SYMBOL='circle', /SYM_FILLED, COLOR='black', $

  NAME='Data', /OVERPLOT)

plo2 = PLOT(x, y_shifted, SYMBOL='circle', /SYM_FILLED, COLOR='red', $

  NAME='Data shifted', /OVERPLOT)

lgd = LEGEND(POSITION=[50, 10], /DATA)

arr = ARROW([0, nn-1-lag], [8, 8], ARROW_STYLE=3, THICK=2, COLOR='blue', /DATA)

여기서는 그래픽창을 띄우고 PLOT 함수를 사용하여 원본 데이터인 y와 -lag만큼 이동된 y_shifted를 서로 비교할 수 있도록 중첩하여 표출하였습니다. 그리고 앞서 산출된 자기상관계수의 값을 문자화하여 타이틀에 표시되도록 하였고, LEGEND 함수를 사용하여 범례도 표시되도록 하였습니다. 그리고 ARROW 함수를 사용하여 수평 방향 화살표를 표시하였는데, 이것은 지난 회 게시물에서 언급했던 "복제 데이터를 한쪽 방향으로 이동시킨 상태에서 원본 데이터와 위치가 겹치는" 부분을 표시하기 위함입니다. 이 내용을 실행하면 다음과 같은 그림을 얻게 됩니다.

이 그림은 lag가 2인 경우에 대한 결과이며, lag의 값을 다양하게 바꿔보면 여러가지 경우들에 대한 그림을 얻을 수 있습니다. 만약 lag의값이 5, 10인 경우에는 다음과 같은 결과들을 얻게 됩니다.

이 그림들을 보면, lag의 값을 다양하게 시도해보면서 그때마다 얻어지는 자기상관계수 값들을 모아보면 뭔가 또 다른 플롯을 그릴 수 있겠다는 생각이 듭니다. 그래서 다음과 같이  다양한 lag 값들로 구성된 lags라는 배열을 생성하여 A_CORRELATE 함수에 투입하면 그 결과 역시 배열의 형태로 얻어집니다. 여기서는 lag의 값을 0, 1, 2, ..., 29, 30과 같은 값들로 구성된 배열을 생성하여 이를 적용한 결과로 산출되는 자기상관계수 값들로 구성된 results라는 배열을 얻고 이 결과를 플롯의 형태로 표출해보았습니다.

lags = [0:30]

results = A_CORRELATE(y, lags)

win_c = WINDOW(DIMENSIONS=[600, 500], /NO_TOOLBAR)

pl_c = PLOT(lags, results, YRANGE=[-1, 1], MARGIN=0.12, $

  SYMBOL='circle', /SYM_FILLED, COLOR='green', $

  XTITLE='Lag', YTITLE='AC Coefficient', FONT_SIZE=12, $

  TITLE='Autocorrelation Plot', /CURRENT)

이와 같이 다양한 lag 값들을 적용하여 얻은 자기상관계수들을 그래프로 표출한 것을 Autocorrelation Plot 또는 다른 용어로 Correlogram이라고 합니다. 이 그림을 보면 lag가 0일 때에는 당연히 자기상관계수가 1.0이지만 lag가 커지면서 다양한 값들이 산출되는데 잘 보면 패턴이 존재합니다. 그런데 1.0 이후로 자기상관계수가 가장 높게 나타나는 경우는 lag가 17 또는 18 정도일 때로 나타나는데요. 이 값이 결국 원래 데이터인 y의 변동주기에 해당됩니다. 실제로 lag의 값이 17, 18인 경우를 앞선 예제 코드를 사용하여 그려보면 다음과 같습니다.

IDL에서 자기상관 기법을 적용한 분석의 방법은 대략 이와 같습니다. 이 내용을 참조하여 여러분이 갖고 계신 데이터에도 적용해보시면 좋을 것 같습니다. 실제로 저도 이 방법을 태양의 장기간 흑점수 데이터에 적용하여 태양활동주기를 계산해보았는데요. 1749년 1월부터 2017년 10월까지 약 268년 동안의 월별 태양흑점수 데이터를 위와 같은 방식으로 분석한 최종 결과를 그림만 보여드리면 다음과 같습니다.

이 그림에 표시되어 있듯이 변동주기가 약 128개월 정도로 산출되는데, 이는 10.7년 즉 약 11년에 해당됩니다. 흔히 잘 알려져있는 태양활동주기의 값과 비슷하게 나옵니다. 물론 변동성을 가진 데이터로부터 변동주기를 얻는 방법은 다양하게 존재합니다. 제가 예전에 비슷한 작업을 Power Spectrum 기법을 이용하여 해본 내용을 올린 적이 있는데, 여기서 나온 결과와 유사합니다. 참고로 이 내용은 여기를 누르시면 보실 수 있습니다.