Python - Zaman Serisi İnceleme ve ARIMA ile Tahmin

Bu yazı kapsamında zaman serisi analizine giriş ve ARIMA modeli ile zamana bağımlı bir değişkenin değerinin nasıl tahmin edilebileceğine dair bir örnek paylaşacağım. Kullanacağımız örnek veri setini indirmek için tıklayınız.

Kaynak : https://www.analyticsvidhya.com/blog/2016/02/time-series-forecasting-codes-python/

import pandas as pd
import numpy as np
import matplotlib.pylab as plt      
        
data = pd.read_csv('AirPassengers.csv')

Örnek verimizde 2 kolon mevcut; yıl-ay bilgisi ve yolcu sayısı. Veri setimiz zamana göre değişen yolcu sayılarını içermekte, zaman serisi analizi için .csv veriyi okurken ay değişkenini datetime formatına çevirip index olarak ayarlayabiliriz.

dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month',date_parser=dateparse)
print (data.head())







Gördüğünüz gibi zaman bilgisini index olarak ayarladık ve artık sadece bir sütunumuz var. Şimdi o sütunu seçelim ve ayrı bir değişkene aktaralım. ( data değişkeni DataFrame  tipinde, ts değişkeni pandas Series tipinde )

 ts = data['#Passengers']

Artık herhangi bir zamandaki veriye şu şekilde erişebiliriz ya da tarih aralığı verebiliriz.

ts['1949-01-01']

ts['1957-01-01':'1957-05-01']

Ya da sadece belirtilen bir yıldaki verileri kolayca çekebiliriz.

Şimdi ise ts içerisindeki tüm değerleri grafik olarak gösterelim.

 plt.plot(ts)




Grafikte gördüğünüz gibi 1949-1960 yılları arasındaki yolcu sayılarında yıl içerisinde değişim olmakla birlikte sürekli bir artış söz konusu. Bu durum zaman serilerinde genel eğilim (trend) olarak adlandırılır (Trend: zaman serilerinde uzun vadede gösterilen kararlı durum). Verileri ay ay incelediğimizde ise yıl içerisinde sezonlara göre dalgalanmalar mevcut. Zaman serilerinde durağanlık zaman içerisinde varyansın ve ortalamanın sabit olmasıdır. Verilerde bir durağanlık analizi yapılmak isteniyorsa, veri üzerinde belirli büyüklükte (mesela son 12 ay) bir kayan pencere gezdirerek son 12 değerin ortalama ve varyansını hesaplatıp inceleyebiliriz. Bunun için bir fonksiyon tanımlayalım.

 def test_stationarity(timeseries):
    
    #Determing rolling statistics
    rolmean = pd.rolling_mean(timeseries, window=12)
    rolstd = pd.rolling_std(timeseries, window=12)

    #Plot rolling statistics:
    orig = plt.plot(timeseries, color='blue',label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)

test_stationarity(ts)




Yukardaki grafiğe göre standart sapmadaki değişim çok az olmakla birlikte, ortalama da ise düzgün (doğrusala yakın) bir artış söz konusu. Grafikte görüldüğü gibi yıl içerisindeki değişimler bir hayli fazla. Log dönüşümü ile yüksek değerleri küçüklere nispeten daha fazla azaltıp, verimizi daha az dalgalı hale getirelim. Yine bu dönüştürülmüş veriye rolling_mean işlemi (pencere gezdirip ortalama alma) uygulayalım.




ts_log = np.log(ts)
moving_avg = pd.rolling_mean(ts_log,12)
plt.plot(ts_log)
plt.plot(moving_avg, color='red')




Pencere kaydırma işlemi ile son 12 değerin ortalamasını aldık, bu sebeple ilk 11 değer null oldu. Şimdi log alınmış değerlerden, moving_avg değerlerini çıkaralım ve null değerleri silelim.

ts_log_moving_avg_diff = ts_log - moving_avg
ts_log_moving_avg_diff.dropna(inplace=True)
test_stationarity(ts_log_moving_avg_diff)

moving_avg  ile fark alınca değerlerdeki  sürekli artış trendi ortadan kalktı ve artık nispeten daha durağan bir hale gelmiş oldu. Zaman serisinde tahmin işlemi gerçekleştirebilmek için veriyi durağan hale getirmemiz lazım, biz de bunun için moving_avg ile çıkarma işlemi yaptık. (bunun için yazının devamında örneklerini bulacağınız farklı yöntemler de kullanılabilir).

 pd.rolling_mean ile son 12 değerin ortalamasını almış olduk ve bu durumda x ayı için kendisinden önceki 11 ayın değerleri aynı ağırlıkta hesaba katmış olduk. Oysa ki bazı durumlarda (değişkenimiz zaman olduğu için) x ayına daha yakın olan zaman dilimlerine daha fazla ağırlık vermek mantıklı olabilir. Bunun için pandas kütüphanesinin exponentially weighted moving average algoritmasını kullanan ewma() fonksiyonunu kullanabiliriz.

expwighted_avg = pd.ewma(ts_log, halflife=12)
log = plt.plot(ts_log, label='Log')
exp = plt.plot(expwighted_avg, color='red', label='Log expwighted_avg')
avg = plt.plot(moving_avg, color='green', label='Log moving_avg')
plt.legend(loc='best')
plt.title('Log - expwighted_avg - moving_avg')
plt.show(block=False)




moving_avg ile fark aldığımız gibi expweighted_avg ile de fark alıp verimizi bir başka şekilde de durağanlaştırabiliriz. Ayrıca son 12 değerin ortalamasını alıp fark almanın dışında, her bir değeri bir önceki ile çıkarıp sadece fark değerlerini de kullanabiliriz.

ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)
ts_log_diff.dropna(inplace=True)




Zaman Serisi Tahmini (ARIMA)




ARIMA ile log dönüşümü uygulanmış ts_log verisine zaman serisi analizi uygulayıp, geçmiş değerlere göre gelecek değerleri tahmin etmeye çalışalım. ARIMA modelini oluştururken 3 parametre veriyoruz, bunlar sırasıyla; p,d ve q. Ayrıntılarına yazının başında kaynak olarak belirttiğim linkten ve ARIMA modelinin dökümanından bakabilirsiniz. p: x(t) anında kaç adım önceki değerlerin tahmin işleminde göz önüne alınacağını, q: x(t) anındaki bir tahmin hatasının kaç adım önceki değerler ile moving average işlemine tabi tutulacağını , d: veriyi durağan hale getirmek için uygulanacak fark alma işlemi derecesini (yukarıda örneklerle anlattığımız: ts_log_moving_avg_diff, ts_log_ewma_diff ya da ts_log_diff) ifade eder. Eğer d parametresini 0 verirsek verimiz durağan hale getirilmeyecektir, dolayısıyla bu parametreyi 0 verdiğimiz zaman input olarak durağan hale getirilmiş ts_log_diff ya da diğer iki veriyi ARIMA modelinde kullanmamız gerekmekte. 

Yani bir başka deyişle ;

model = ARIMA(ts_log, order=(2, 1, 2))  ile 

model = ARIMA(ts_log_diff, order=(2, 0, 2))  aynı sonucu verecektir.

from statsmodels.tsa.arima_model import ARIMA       


model = ARIMA(ts_log, order=(2, 1, 2))  
results_ARIMA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))

Log dönüşümü yapılmış ts_log verisi ile bir ARIMA modeli oluşturup, eğittik ve .fit() ile tahmin edilen değerleri almış olduk. Şimdi tahmin edilen değerlerden gerçek ölçeğe dönebilmek için, ARIMA ile tahmin edilen değerlerin kümülatif toplamını bulalım

predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)

predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()




ts_log un ilk değerini baz alıp, kümülatif toplamları ekleyelim.

predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index)

predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)




Ve son olarak log alınmış değerlerin e inci değerden üstünü alalım, orjinal değerler (ts) ile tahmin edilmiş değerleri (predictions_ARIMA) kıyaslayalım

predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)

plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))







Sonuç itibariyle çok iyi bir tahmin olmasa da orjinal değişime yakın bir grafik elde etmiş olduk. Referans olarak kullandığım kaynağa ve biraz daha detaylı bilgiye linkten ulaşabilirsiniz.