An Example using ARIMA in Prediction using Thailand GDP (1)

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

https://qtmlresearch.com/2017/06/07/arima-part-1-the-basic/

 

Data : GDP of Thailand during (1960 = 2015)

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

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

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

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

# reverse data from past to present
data = data.reindex(index=data.index[::-1])

thGDP = data['GDP']

## ***** Testing for stationarity ***** ##
thGDP.plot()
plt.title('Thailand GDP 1960 - 2015')
plt.grid(True)
Screen Shot 2017-06-16 at 17.48.36
Fig. 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)

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

Screenshot from 2017-06-17 15-55-17

Fig. 2

วิธีการที่พื้นฐานที่สุดที่จะอ่านผล เราสามารถดูได้จากค่า 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 สิ่ง คือ Tred (แนวโน้มการเปลี่ยนแปลของข้อมูล) และ 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 differencing)

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

thGDP_log = np.log(thGDP)

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

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

figure_1-1

Fig. 3

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

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

sta_result = stationarity_test(thGDP_log)

Screenshot from 2017-06-17 16-05-33

Fig. 4

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

3.2 – Removing those trends from the data

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

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

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

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

figure_1.png

Fig. 5

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

sta_result = stationarity_test(thGDP_log_diff)

Screenshot from 2017-06-17 16-08-50

Fig. 6

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

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

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

https://qtmlresearch.com/2017/06/07/arima-part-1-the-basic/

 

13062103_233130583716199_6224283738921716075_n
Enjoy STATISTIC
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s