การพยากรณ์ข้อมูลอนุกรมเวลาด้วยเทคนิค ARIMA ด้วย Python

Published in

QUANT I LOVE U

5 min readApr 21, 2020

ในปี 1970 มีการเสนอเทคนิคการพยากรณ์ข้อมูลที่ชื่อว่า Autoregressive integrated moving average หรือ ARIMA ขึ้นมา

โดยเทคนิค ARIMA เป็นหนึ่งในการวิเคราะห์ข้อมูลแบบอนุกรมเวลา คือการอาศัยข้อมูลในอดีตเพื่อกำหนดรูปแบบของข้อมูลและพยากรณ์ข้อมูลในอนาคต

องค์ประกอบของแบบจำลอง ARIMA(p,d,q)

องค์ประกอบของแบบจำลอง ARIMA ประกอบด้วย 3 ส่วน Autoregressive process หรือ AR(p), Integrated(d) และ Moving average process: MA(q)

Autoregressive process

Autoregressive process หรือ AR(p) เป็นกระบวนการอธิบายตัวแปล Y ด้วยค่าของตัวแปล Y(Yt-1, Yt-2, Yt-3,…, Yt-p) เอง

Y ที่เราใช้สามารถมีได้หลาย order ตัวอย่างดังนี้

AR(1): First order autoregressive

AR(2): Second order autoregressive

AR(p): (p) order autoregressive

order ที่มากขึ้นคือการใช้จำนวน Y มากขึ้น ซึ่งมีข้อดีและข้อเสียการที่มี Y เยอะ ถ้าจำนวน Y น้อยมีข้อดีคือสามารถพยากรณ์เมื่อ Y มีการเปลี่ยนเทรนได้ดี แต่มีข้อเสียหาก Y หลังๆ มีค่าผิดปกติจะทำให้การพยากรณ์ผิดพลาด

หากจำนวน Y สูงมีข้อดีคือไม่อ่อนไหวต่อค่าปกติ และข้อเสียคือไม่สามารถพยากรณ์เมื่อเปลี่ยนเทรนได้

Integrated(d)

d คือจำนวนความแตกต่างที่จำเป็นในการทำให้อนุกรมเวลาคงที่

Moving average process: MA(q)

คือการดูโมเดลว่าพึ่งพาระหว่างการสังเกตและข้อผิดพลาดที่เหลือจากแบบจำลองค่าเฉลี่ยเคลื่อนที่ที่ใช้กับการสังเกตที่ยังไง

MA(1): First Order Moving Average

MA(2): Second Order Moving Average

MA(q): (q) Order Moving Average

ตัวอย่างสมการพยากรณ์

เริ่มต้นด้วยการ Import liberly

from dateutil.parser import parseimport matplotlib as mplimport matplotlib.pyplot as pltimport seaborn as snsimport numpy as npimport pandas as pdimport statsmodels.formula.api as smfimport statsmodels.api as sm

Import data

df df =  pd.read_csv('https://raw.githubusercontent.com/AileenNielsen/TimeSeriesAnalysisWithPython/master/data/AirPassengers.csv', parse_dates=['Month'])df.rename(columns={'#Passengers':'Passengers'},inplace=True)df.head()

เราลองสร้างกราฟดูข้อมูลกันก่อน

df['Passengers'].plot(figsize = (20, 5), legend = True);

# Plot Two Side View
fig, ax = plt.subplots(1, 1, figsize=(16,5), dpi= 120)
x = df['Month'].valuesy1 = df['Passengers'].valuesplt.fill_between(x, y1=y1, y2=-y1, alpha=0.9, linewidth=2, color='skyblue')plt.ylim(-700, 700)plt.title('Air Passengers (Two Side View)', fontsize=20)plt.hlines(y=0, xmin=np.min(df.Month), xmax=np.max(df.Month), linewidth=0.5)plt.show();

การทดสอบ Stationary

เราจำเป็นต้องทำให้ข้อมูลเป็น Stationary หรือมีความนิ่ง เนื่องจากคำว่า ‘Auto Regressive’ ใน ARIMA หมายถึงเป็นรูปแบบการถดถอยเชิงเส้นที่ใช้ความล่าช้าของตัวเองเป็นตัวทำนาย ARIMA แบบจำลองการถดถอยเชิงเส้นเหมือนแบบจำลองเชิงเส้นทั่วไปคือทำงานได้ดีเมื่อตัวพยากรณ์ไม่สัมพันธ์และเป็นอิสระจากกัน

rolling_mean = df.rolling(window = 12).mean()rolling_std = df.rolling(window = 12).std()fig, ax = plt.subplots(1, 1, figsize=(16,5), dpi= 120)plt.plot(df, color = 'blue', label = 'Original')plt.plot(rolling_mean, color = 'red', label = 'Rolling Mean')plt.plot(rolling_std, color = 'black', label = 'Rolling Std')plt.legend(loc = 'best')plt.title('Rolling Mean & Rolling Standard Deviation')plt.show()

result = adfuller(df['Passengers'])print('ADF Statistic: {}'.format(result[0]))print('p-value: {}'.format(result[1]))print('Critical Values:')for key, value in result[4].items():print('\t{}: {}'.format(key, value))

ผลที่ได้ออกมาจะเห็นได้ว่าค่า P value ค่อนข้างสูง ดังนั้นอาจกล่าวได้ว่าข้อมูลเป็น non stationary และเพื่อความสะดวกเราจะเขียนฟังชั่นเพื่อตรวจสอบความเป็นstationary

def get_stationarity(timeseries):# rolling statisticsrolling_mean = timeseries.rolling(window=12).mean()rolling_std = timeseries.rolling(window=12).std()# rolling statistics plotoriginal = plt.plot(timeseries, color='blue', label='Original')mean = plt.plot(rolling_mean, color='red', label='Rolling Mean')std = plt.plot(rolling_std, color='black', label='Rolling Std')plt.legend(loc='best')plt.title('Rolling Mean & Standard Deviation')plt.show(block=False)# Dickey–Fuller test:result = adfuller(timeseries['Passengers'])print('ADF Statistic: {}'.format(result[0]))print('p-value: {}'.format(result[1]))print('Critical Values:')for key, value in result[4].items():print('\t{}: {}'.format(key, value))

การปรับข้อมูลให้เรียบมีหลักๆ 3 วิธีคือ ลบด้วยค่าเฉลี่ย 12 เดือน ลบด้วยค่าก่อนหน้า และ take log

วิธีการแรกในการปรับข้อมูลให้เรียบคือการลบด้วยค่าเฉลี่ยเคลื่อนที่ 12เดือนย้อนหลังออกไป เพื่อขจัด ฤดูกาลกับแนวโน้มนั้นเอง

rolling_mean = df_log.rolling(window=12).mean()df_log_minus_mean = df_log - rolling_meandf_log_minus_mean.dropna(inplace=True)get_stationarity(df_log_minus_mean)

เห็นได้ว่า P value ที่ออกมาค่อข้างตํ่าแล้ว ที่นี้เราลองวิธีที่สองtake log กันบ้าง

rolling_mean_exp_decay = df_log.ewm(halflife=12, min_periods=0, adjust=True).mean()df_log_exp_decay = df_log - rolling_mean_exp_decaydf_log_exp_decay.dropna(inplace=True)get_stationarity(df_log_exp_decay)

เห็นได้ว่า P value น้อยลงกว่าเดิม

สุดท้ายลบด้วยค่าก่อนหน้า

df_log_shift = df_log - df_log.shift()df_log_shift.dropna(inplace=True)get_stationarity(df_log_shift)

สร้าง ARIMA Model

ได้เข้าหัวข้อหลักของเราสักที วิธีเเรกของเราใช้ statsmodels

สิ่งที่เราต้องใส่คือ datafram และ order โดย order คือ number of time lags(p,d,q) เราจะใส่ order อย่างไรจริงๆสามารภหาได้จากการดู Auto correlation แต่สำหรับผมใช้อีกวิธีคือลองใส่ๆไปใน modeldก่อน ตัวอย่างแรกให้ order=(1,1,2)

# 1,1,2 ARIMA Modelmodel = ARIMA(df.Passengers, order=(1,1,2))model_fit = model.fit(disp=0)print(model_fit.summary())

ตารางมีข้อมูลค่อนข้างเยอะ แต่ตัวที่เราสนใจจริงๆแล้วคือช่อง coef และ ‘P>|z|’

‘coef’ คือน้ำหนักของตัวแปรที่เกี่ยวข้อง

‘P>|z|’ คือระดับความเชื่อมั่น

# 1,1,1 ARIMA Modelmodel = ARIMA(df.Passengers, order=(1,1,1))model_fit = model.fit(disp=0)print(model_fit.summary())

จะเห็นได้ว่า ‘P>|z|’ ดีกว่าผมเลือกอันนี้ (อันนี้เป็นวิธีส่วนตัวอาจจะไม่ค่อยถูกตามกลักการสักท่าใหร่) ลองพล็อต residual errors ว่าข้อมูลใช้ได้ใหม

residuals = pd.DataFrame(model_fit.resid)fig, ax = plt.subplots(1,2)residuals.plot(title="Residuals", ax=ax[0])residuals.plot(kind='kde', title='Density', ax=ax[1])plt.show()

ตรวจสอบอีกครั้งค่าที่ได้จากแบบจำลองกับค่าจริง

# Actual vs Fittedmodel_fit.plot_predict(dynamic=False)plt.show()

ถือว่าไกล้เคียงปิดงานได้

pmdarima

อีกวิธีสำหรับวันนี้คือใช้ pmdarima ข้อดีของวิธีนี้คือผมไม่ต้องมาเดา order ให้ใส่ order ตํ่าสุดสูงสุดไปแทน model จะเลือกให้เรา (วิธีการข้อผมคือถ้าข้อมูลไม่มากกว่าใส่กว้างๆ )เริ่มจากการเทรน model

# Fit auto_arima function to AirPassengers datasetstepwise_fit = auto_arima(df['Passengers'], start_p = 1, start_q = 1,max_p = 3, max_q = 3, m = 12,start_P = 0, seasonal = True,d = None, D = 1, trace = True,error_action ='ignore', # we don't want to know if an order does not worksuppress_warnings = True, # we don't want convergence warnings# To print the summarystepwise_fit.summary()

ตารางของวิธีนี้จะเบี่ยวๆนิดหนึ่งแต่ใช้การอ่านเหมือนกับวิธีเเรก

สุดท้าย

ในการใช้งานจริงของผมก็ดูเพียงแนวโน้มเท่านั้นว่าประมาณใหนควรทำอย่างไรสำหรับผมไม่ค่อยอ่านค่าตัวเลขแบบเฉพาะเจาะจง เพราะอ่านไปก็ไม่ถูกอยู่ดี จึงจะสร้างกราฟมาดูแนวโน้ม

# Forecast for the next 3 yearsforecast = result.predict(start = len(df),end = (len(df)-1) + 3 * 12,typ = 'levels').rename('Forecast')# Plot the forecast valuesdf['Passengers'].plot(figsize = (12, 5), legend = True)forecast.plot(legend = True)