จากบทความที่แล้ว ตามลิงก์ด้านล่าง

https://algoaddict.com/blog/54520/arima

เราได้พูดถึง Basic of ARIMA เอาไว้ ซึ่งถึงแม้ว่าตัวโมเดลจะไม่ได้ถูกออกแบบมาโดยเฉพาะสำหรับการทำนายมากนัก แต่ก็เป็นโมเดล Time-series ที่สามารถนำมาประยุกต์ใช้ในการทำนายข้อมูลได้ ดังนั้น วันนี้เราจึงได้นำตัวอย่างการประยุกต์ใช้ ARIMA พื้นฐานในการทำนายอย่างง่ายมาให้ดูกัน โดยตัวอย่างนี้ใช้ จะเป็นข้อมูล GDP ของประเทศไทย ดังนี้

Data: GDP of Thailand during (1960 - 2017)

File name: GDP_Thailand.csv


1) มาดูข้อมูลดิบกันก่อน


ก่อนอื่นมาโหลดข้อมูล และ เนื่องจากข้อมูลของเราเรียงจากปีปัจจุบันไปยังอดีต แต่เพื่อให้ง่ายต่อความเข้าใจ จากนั้น เราจะได้พล็อตดูแนวโน้มข้อมูลกันต่อไป

1.1) Import libraries ที่จำเป็น


import pandas as pd
import matplotlib.pyplot as plt
import numpy as np


1.2) อ่านข้อมูล และ ตั้งชื่อ คอลัมน์ว่า GDP


data = pd.read_csv('GDP_Thailand.csv', parse_dates = True, index_col = 0)
data.columns=['GDP']


1.3) ใส่ข้อมูลไว้ในตัวแปรชื่อ thGDP เพื่อใช้งานในอนาคต


thGDP = data['GDP']


1.4) ทดลองพล็อตข้อมูลดู


## ***** Testing for stationarity ***** ##
thGDP.plot()
plt.title('Thailand GDP 1960 - 2015')
plt.grid(True)



ig 1: พล็อตข้อมูลดิบ


2) เช็คความเป็น Stationary ของข้อมูล

ตามที่ได้กล่าวไว้ในบทความก่อนหน้า อันดับแรกเราจะต้องทำการตรวจสอบข้อมูลก่อนว่าข้อมูลนั้นเป็น Stationary หรือไม่ ซึ่งในการตรวจสอบเราอาจจะดูจากสายตาได้ (ซึ่งจากกราฟในขั้นตอนที่ 1 เห็นได้ชัดว่าข้อมูลมี Trend ชัดเจน ซึ่งมีความเป็นไปได้สูงมากที่ข้อมูลจะไม่เป็น stationary) อย่างไรก็ตาม เราไม่สามารถสังเกตุได้ด้วยตาเปล่าเสมอไป ดังนั้น เราจะใช้ค่าสถิติ “Dickey-Fuller Test” ในการตรวจสอบ

ในที่นี้เราจะสร้างฟังก์ชันง่ายๆ เพื่อใช้ตรวจสอบความเป็น Stationary ดังนี้


from statsmodels.tsa.stattools import adfuller

def stationarity_test(timeseries):

print('Augment Dicky Fuller Test:')
adftest = adfuller(timeseries, autolag='AIC')
adfoutput = pd.Series(adftest[0:4], index=['Test Statistic','P-Value','#Lags','Number of Observations'])

print('Test_Statistic :' + str(adfoutput[0]))
print('p-values :' + str(adfoutput[1]))
print('Number of observation :' + str(adfoutput[3]))

จากนั้น เรียกใช้ฟังก์ชั่นด้านบนเพื่อทดสอบความเป็น stationary ของข้อมูลได้ โดยการส่ง Timeseries ที่ต้องการตรวจสอบไปคำนวณ ดังนี้

stationarity_test(thGDP)

ผลลัพธ์ที่สำคัญที่เราจะได้จากฟังก์ชัน มีดังนี้

Fig 2. ผลที่ได้จากการเทส Stationarity


วิธีการที่พื้นฐานที่สุดที่จะอ่านผล เราสามารถดูได้จากค่า p-values ซึ่งในที่นี้มีค่าประมาณ 0.99 ซึ่งมากกว่า 0.05 ซึ่งสอดคล้องกับกราฟที่เราสังเกตุได้ด้วยตาเปล่าว่า ข้อมูลชุดนี้ ไม่เป็น Stationary

ค่า p-value = 0.99 หมายความว่า มีความเชื่อมั่นว่าข้อมูลชุดนี้จะเป็น Stationary ที่ ประมาณ 1% เท่านั้น (100 – 99) ซึ่งในการทดสอบสมมุติฐานทางสถิติโดยปกติแล้ว เราต้องการความเชื่อมั่นที่ ประมาณ 90% – 95% นั่นคือ ค่า p-value ที่เราอยากจะได้เพื่อยืนยันความเป็น Stationary ของข้อมูล คือ ค่าที่น้อยกว่า 0.1 – 0.05 เป็นต้น

ในบทความนี้ เราได้กำหนดให้ต้องมีความเชื่อมั่น 95% ขึ้นไปเท่านั้น เราจึงจะยอมรับว่าข้อมูลนี้เป็น Stationary ดังนั้น เราจะปรับข้อมูลไปจนกว่าจะได้ค่า p-value < 0.05 เพื่อให้ข้อมูลเป็น Stationary ซึ่งจะหมายความว่า ข้อมูลนั้นพร้อมสำหรับการทำนายด้วย ARIMA แล้ว


3) ทำให้ข้อมูล Stationary

ดังที่กล่าวไว้ในในบทความที่แล้ว ความเป็นไม่เป็น Stationary ของข้อมูล time seriesโดยหลักการแล้ว ขึ้นอยู่กับ 2 สิ่ง คือ Trend (แนวโน้มการเปลี่ยนแปลของข้อมูล) และ Seasonality (การเกิดแพทเทิร์นซ้ำๆ ณ ช่วงเวลาหนึ่งๆ) ดังนั้น การแก้ไขปัญหาที่สามารถทำได้ก็คือ การกำจัด Trend และ Seasonality ออกไปจากข้อมูลนั่นเอง

การกำจัด Trend และ Seasonality จากข้อมูลสามารถทำได้หลายวิธี ตัวอย่างเช่น การใช้ Log, Cube root, Aggregation, Smoothing, Polynomial Fitting เป็นต้น อย่างไรก็ตาม ถ้าข้อมูลของเรานอกจากจะมีปัญหาเรื่อง Trend แล้ว ยังมีปัญหาเรื่อง Seasonality

ด้วย ปัญหานี้ก็จะซับซ้อนขึ้น ทำให้ในการกำจัด Sesonality อาจจะต้องใช้วีธีอื่นๆ ร่วมด้วย ยกตัวอย่างเช่น การทำ Decomposition เป็นต้น

ในที่นี้ เดี๋ยวเราจะมาลองใช้ขั้นตอนพื้นฐาน ตามที่ได้กล่าวไว้ในบทความที่แล้วเพื่อกำจัด Trend ออกไปจากข้อมูลค่ะ


3.1 – Modelling/ estimating trend (taking log and differentiating)

อันดับแรกหา Logarithm ของข้อมูล

thGDP_log = np.log(thGDP)

ลองพล็อตข้อมูลดูก่อน

thGDP_log.plot()
plt.title('Logarithm of Thailand GDP')
plt.grid(True)



Fig 3. พล็อตดาต้าหลังจากการทำ log


มองจากกราฟ สิ่งที่สังเกตุได้ชัดเจนก็คือ ความ smooth ของข้อมูล ซึ่งจะเห็นได้ชัดว่าข้อมูล smooth ขึ้นอย่างชัดเจนเทียบกับรูปใน Fig. 1

มาลองทดสอบกับ Dickey-Fuller Test ดูกันดีกว่าค่ะ ว่าข้อมูล smooth ขึ้นแบบนี้จะเพียงพอต่อการทำให้ข้อมูลเป็น stationary หรือ ไม่ โดยการเรียกใช้ฟังก์ชัน “test_stationarity” ที่เราสร้างไว้ในขั้นตอนที่ 2 อีกครั้ง

sta_result = stationarity_test(thGDP_log)

Fig 4. ผลจากการเทส Stationary


ค่า p-value = 0.57 ซึ่งมากกว่า 0.05 ซึ่งหมายความว่าข้อมูลยัง ไม่เป็น Stationary ไปขั้นตอนต่อไปกันเลยค่ะ


3.2 – กำจัดแนวโน้ม (Trends) ออกจากข้อมูล

ขั้นตอนต่อไป เราจะมาทำ Differentiating ข้อมูลกันเพื่อ โมเดล Trend และ ทำการ Remove Trends นั้น ออกจากข้อมูล ด้วยคำสั่ง

thGDP_log_diff = thGDP_log - thGDP_log.shift()
thGDP_log_diff.dropna(inplace = True)

มาพล็อตดูกันค่ะ ว่าหลังจาก differentiating แล้ว ข้อมูลของเราเปลี่ยนไปอย่างไรบ้าง

thGDP_log_diff.plot()
plt.title('Differencing of Thailand GDP')
plt.grid(True)


Fig 5. พล็อตข้อมูลหลังจากการ remove trends


จากกราฟใน Fig. 5 จะเห็นว่าการทำ differencing ไม่ได้หมายความว่าจะทำให้ข้อมูล smooth ขึ้นเหมือนการทำ logarithm แต่ประโยชน์ของการทำ differencing ก็คือ การกำจัด Trend ออกจากข้อมูลค่ะ จากกราฟนี้ เราจะเห็นได้ว่า Trend ของข้อมูลชุดนี้ ได้ถูกกำจัดออกไป จนแทบจะสังเกตุ Trend ไม่ได้เลย มาดูค่าทางสถิติกันดีกว่าค่ะ ว่าข้อมูลมีความเป็น Stationary เพียงพอที่จะนำไปทำนายหรือยัง โดยการคำนวณจากฟังก์ชัน “test_stationarity” ที่เราสร้างไว้ในขั้นตอนที่ 2 อีกครั้ง


sta_result = stationarity_test(thGDP_log_diff)


Fig 6. ผลจากการเทส Stationariy


จะเห็นว่า ข้อมูลหลังทำ Differentiating ค่า p-value มีค่าตำ่กว่า 0.05 (0.00014171) ซึ่งหมายความว่า ข้อมูลมีความเชื่อมั่นที่จะเป็น Stationary ที่ความมั่นใจสูงมากกว่า 95% ดังนั้น เราก็จะขอหยุดการทำ differentiating ไว้เพียงแค่ระดับนี้ (ถ้าข้อมูลยังไม่เป็น stationary เราสามารถทำ differentiating ต่อไปได้) และ ใช้ข้อมูลนี้ในขั้นตอนการทำนายข้อมูลโดย ARIMA ต่อไป

ตอนแรกตั้งใจว่าจะเขียนแค่บทความสั้นๆ เขียนไปเขียนมายาวอีกแล้ว เลยขอแยกออกเป็น 2 บทความแล้วกันค่ะ บทความนี้ขอให้เข้าใจการจัดการข้อมูลให้พร้อมสำหรับการทำนายด้วย ARIMA แล้วในบทความถัดไป เราก็จะจับเอาข้อมูลนี้ที่เราจัดการไว้ให้เป็น Stationary เรียบร้อยแล้ว ซึ่งก็คือ time series (thGDP_log_diff) มาลองทำนายด้วย ARIMA กัน

ลิงก์บทความตอนที่ 2

https://algoaddict.com/blog/54679/arima_th_gdp_2

ลิงก์บทความก่อนหน้า สำหรับท่านที่อ่านบทความนี้แล้วยังไม่ค่อยเข้าใจ ก็สามารถกลับไปอ่านทฤษฏีที่เคยเขียนไว้ได้ค่ะ

https://algoaddict.com/blog/54520/arima