Zaman serisi analizi

12. Zaman serisi analizi#

Zaman serisi, zaman içinde değişen bir sistemden alınan ölçümler dizisidir. Önceki bölümlerde kullandığımız regresyon gibi araçların çoğu zaman serilerinde de kullanılabilir; ayrıca bu veri türü için özellikle yararlı yöntemler vardır.

Örnek olarak iki veri kümesini inceleyeceğiz: 2001–2024 arasında ABD’deki yenilenebilir elektrik üretimi ve aynı dönemdeki hava durumu verileri. Bir zaman serisini uzun dönemli eğilim ile yinelenen mevsimsel bileşene ayırma yöntemleri geliştireceğiz. Eğilimleri uydurup öngörmek için doğrusal regresyon kullanacak ve resmî adı “otoregresif bütünleşik hareketli ortalama”, kısa adı ARIMA olan yaygın bir modeli deneyeceğiz.

Think Stats kitabının üçüncü baskısını Bookshop.org ve Amazon üzerinden edinebilirsiniz (bunlar satış ortaklığı bağlantılarıdır). Ücretsiz çevrimiçi sürümü beğeniyorsanız bana bir kahve ısmarlamayı düşünebilirsiniz.

Bu notebook’u Colab’da çalıştırmak için buraya tıklayın.

Kod hücresi içeriğini gizle

 
from os.path import basename, exists


def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve

        local, _ = urlretrieve(url, filename)
        print("İndirildi: " + local)


download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/nb/thinkstats.py")

Kod hücresi içeriğini gizle

 
try:
    import empiricaldist
except ImportError:
    %pip install empiricaldist

Kod hücresi içeriğini gizle

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

from thinkstats import decorate

plt.rcParams["figure.dpi"] = 300

12.1. Elektrik#

Zaman serisi örneği olarak ABD Enerji Bilgi İdaresi’nin 2001–2024 arasında yenilenebilir kaynaklardan aylık toplam elektrik üretimini içeren veri kümesini kullanacağız. İndirme yönergeleri bu bölümün notebook’unda yer alır.

Aşağıdaki hücre, 17 Eylül 2024’te https://www.eia.gov/electricity/data/browser/ adresinden indirdiğim verileri indirir.

Kod hücresi içeriğini gizle

 
filename = "Net_generation_for_all_sectors.csv"
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/" + filename)

Verileri yükledikten sonra çalışmayı kolaylaştıracak biçime getirmek için bazı dönüşümler yapmalıyız.

elec = (
    pd.read_csv("Net_generation_for_all_sectors.csv", skiprows=4)
    .drop(columns=["units", "source key"])
    .set_index("description")
    .replace("--", np.nan)
    .transpose()
    .astype(float)
)

Yeniden biçimlendirilmiş veri kümesinde her sütun gigavatsaat (GWh) cinsinden aylık toplamlar dizisidir. Elektrik kaynaklarını, yani “sektörleri” gösteren sütun etiketleri şöyledir.

elec.columns

Index(['Net generation for all sectors', 'United States',
       'United States : all fuels (utility-scale)', 'United States : nuclear',
       'United States : conventional hydroelectric',
       'United States : other renewables', 'United States : wind',
       'United States : all utility-scale solar', 'United States : geothermal',
       'United States : biomass',
       'United States : hydro-electric pumped storage',
       'United States : all solar',
       'United States : small-scale solar photovoltaic'],
      dtype='object', name='description')

İndisteki etiketler ay ve yılları gösteren dizelerdir; ilk 12’si şöyledir.

elec.index[:12]

Index(['Jan 2001', 'Feb 2001', 'Mar 2001', 'Apr 2001', 'May 2001', 'Jun 2001',
       'Jul 2001', 'Aug 2001', 'Sep 2001', 'Oct 2001', 'Nov 2001', 'Dec 2001'],
      dtype='object')

Bu dizeleri Pandas Timestamp nesneleriyle değiştirirsek verilerle çalışmak kolaylaşır. Ocak 2001’den başlayan Timestamp dizisini date_range ile oluşturabiliriz. "ME" frekans kodu “ay sonu” anlamına gelir ve her ayın son gününü üretir.

elec.index = pd.date_range(start="2001-01", periods=len(elec), freq="ME")
elec.index[:6]

DatetimeIndex(['2001-01-31', '2001-02-28', '2001-03-31', '2001-04-30',
               '2001-05-31', '2001-06-30'],
              dtype='datetime64[ns]', freq='ME')

İndis artık NumPy’de tanımlanan datetime64[ns] veri türünde bir DateTimeIndex’tir. 64, her etiketin 64 bit kullandığını; ns, nanosaniye hassasiyetini gösterir.

12.2. Ayrıştırma#

İlk örnekte ABD’deki nükleer reaktörlerin elektrik üretiminin Ocak 2001–Haziran 2024 arasında nasıl değiştiğine bakıp seriyi uzun dönemli eğilim ve dönemsel bileşene ayıracağız. Aylık toplam üretimler şöyledir.

actual_options = dict(color="C0", lw=1, alpha=0.6)
trend_options = dict(color="C1", alpha=0.6)
pred_options = dict(color="C2", alpha=0.6, ls=':')
model_options = dict(color="gray", alpha=0.6, ls='--')

nuclear = elec["United States : nuclear"]
nuclear.plot(label="nükleer", **actual_options)

decorate(ylabel="GWh")
_images/cea5ccdc7ea7c446533eea9a130098a4a023cc416ac43f596e4fc4b3995d4889.png

Bazı artış ve azalışlar var gibi görünüyor; ancak aydan aya büyük değişimler nedeniyle açıkça görmek zordur. Uzun dönemli eğilimi belirginleştirmek için rolling ve mean ile hareketli ortalama hesaplayabiliriz.

trend = nuclear.rolling(window=12).mean()

window=12, 12 aylık örtüşen aralıklar seçer: ilk aralık ilk ölçümden başlayan 12, ikinci aralık ikinci ölçümden başlayan 12 ölçümü içerir. Her aralığın ortalama üretimini hesaplarız.

Özgün verilerle birlikte sonuçlar şöyledir.

nuclear.plot(label="nükleer", **actual_options)
trend.plot(label="eğilim", **trend_options)
decorate(ylabel="GWh")
_images/0eaf7f0efbc3e863d8684d5ed731b7863d43e7b218d49cd205520b3df80c5af1.png

Eğilim hâlâ oldukça değişkendir. Daha uzun pencereyle daha fazla yumuşatabilirdik; şimdilik 12 aylık pencereyi kullanacağız.

Özgün veriden eğilimi çıkarırsak uzun dönemli ortalaması yaklaşık sabit olan “eğilimden arındırılmış” zaman serisini elde ederiz. Şöyle görünür.

detrended = (nuclear - trend).dropna()
detrended.plot(label="eğilimden arındırılmış", **actual_options)
decorate(ylabel="GWh")
_images/b370e3908de2b672c3c51ad77e330d37bda6adda52f3a9195349edd5cfc5758f.png

Yıllık yinelenen bir örüntü vardır; elektrik talebinin kışın ısınma, yazın klima nedeniyle mevsimlere göre değiştiği düşünülürse bu anlamlıdır. Yıllık örüntüyü açıklamak için indisteki datetime nesnelerinin ay bölümünü seçip verileri aya göre gruplar ve ortalama üretimi hesaplarız. Aylık ortalamalar şöyledir.

monthly_averages = detrended.groupby(detrended.index.month).mean()
monthly_averages.plot(label="aylık ortalama", **actual_options)
decorate(ylabel="GWh")
_images/1bb45b7dbccebabaf88c0968c710fc77ed7a4c557d8c45d963e0a9509394c0c3.png

X ekseninde 1 Ocak, 12 Aralık ayıdır. Elektrik üretimi en soğuk ve en sıcak aylarda en yüksek, Nisan ve Ekim’de en düşüktür.

monthly_averages ile verinin mevsimsel bileşenini oluşturabiliriz. Bu, nuclear ile aynı uzunlukta ve her ayın öğesi o ayın ortalaması olan bir seridir. Şöyle görünür.

seasonal = monthly_averages[nuclear.index.month]
seasonal.index = nuclear.index
seasonal.plot(label="mevsimsel", **actual_options)
decorate(ylabel="GWh")
_images/9664f1c4a9abdb0b480ff57ff3a6bca52009d009520109b25ba91612419d4d4d.png

Her 12 aylık dönem diğerleriyle aynıdır.

Eğilim ile mevsimsel bileşenin toplamı her ayın beklenen değerini temsil eder.

expected = trend + seasonal

Özgün seriyle karşılaştırması şöyledir.

expected.plot(label="beklenen", **pred_options)
nuclear.plot(label="gerçek", **actual_options)
decorate(ylabel="GWh")
_images/62df1fa2347c1892f6cc3de4edc6e318cd40899211cddc958b2387b7e847cd55.png

Bu toplamı özgün seriden çıkarırsak her ayın beklenen değerden sapmasını temsil eden artık bileşeni elde ederiz.

resid = nuclear - expected
resid.plot(label="artık", **actual_options)
decorate(ylabel="GWh")
_images/3613d240a73e50b6be3cf765129568e8ee4d6cb61313f7e8c6cf3fadf8a4254f.png

Artığı, enerji üretimini etkileyip uzun dönemli eğilim veya mevsimsel bileşenle açıklanmayan her şeyin toplamı olarak düşünebiliriz. Buna hava durumu, bakımdaki donanım ve belirli olaylardan doğan talep değişiklikleri dâhildir. Artık çok sayıda öngörülemez, bazen bilinemez etmenin toplamı olduğundan onu çoğu zaman rastgele nicelik kabul ederiz.

Artıkların dağılımı şöyledir.

from thinkstats import plot_kde

plot_kde(resid.dropna())
decorate(xlabel="Artık (GWh)", ylabel="Yoğunluk")
_images/4ad66ef8cb963dc745761cd464f94f1624c0f89c4308b74466b1e12886f739a3.png

Normal dağılımın çan eğrisine benzer; bu, çok sayıda rastgele katkının toplamı olduğu varsayımıyla tutarlıdır.

Modelin özgün seriyi ne kadar iyi açıkladığını nicelleştirmek için artık varyansının özgün seri varyansına göre ne kadar küçük olduğunu gösteren belirleme katsayısını hesaplayabiliriz.

rsquared = 1 - resid.var() / nuclear.var()
rsquared

np.float64(0.9054559977517084)

\(R^2\) yaklaşık 0,92’dir; uzun dönemli eğilim ile mevsimsel bileşen serideki değişkenliğin %92’sini açıklar. Bu değer önceki bölümdekilerden çok yüksektir; özellikle modeli verilere benzeyecek biçimde kurduğumuz böyle bir durumda zaman serilerinde bu yaygındır.

İzlediğimiz sürece mevsimsel ayrıştırma denir. StatsModels bunu yapan seasonal_decompose işlevini sağlar.

from statsmodels.tsa.seasonal import seasonal_decompose

decomposition = seasonal_decompose(nuclear, model="additive", period=12)

model="additive" argümanı toplamsal modeli belirtir; seri eğilim, mevsimsel bileşen ve artık toplamına ayrılır. Çarpımsal modeli birazdan göreceğiz. period=12, mevsimsel bileşenin 12 ay sürdüğünü belirtir.

Sonuç üç bileşeni içeren bir nesnedir. Bu bölümün notebook’u bunları çizen bir işlev sağlar.

def plot_decomposition(original, decomposition):
    plt.figure(figsize=(6, 5))

    ax1 = plt.subplot(4, 1, 1)
    plt.plot(original, label="Özgün", color="C0", lw=1)
    plt.ylabel("Özgün")

    plt.subplot(4, 1, 2, sharex=ax1)
    plt.plot(decomposition.trend, label="Eğilim", color="C1", lw=1)
    plt.ylabel("Eğilim")

    plt.subplot(4, 1, 3, sharex=ax1)
    plt.plot(decomposition.seasonal, label="Mevsimsel", color="C2", lw=1)
    plt.ylabel("Mevsimsel")

    plt.subplot(4, 1, 4, sharex=ax1)
    plt.plot(decomposition.resid, label="Artık", color="C3", lw=1)
    plt.ylabel("Artık")

    plt.tight_layout()

plot_decomposition(nuclear, decomposition)
_images/b2e026d8cd6c0cb2fe2d3bcb23e484f568c5f8c2d3f4678d9c3f0972803e9586.png

Uygulama ayrıntılarından doğan küçük farklarla birlikte sonuçlar kendi hesapladıklarımıza benzer.

Bu ayrıştırma zaman serisinin yapısını anlamayı sağlar; sonraki bölümde göreceğimiz gibi öngörü üretmekte de yararlıdır.

12.3. Öngörü#

Geleceği öngörmek için mevsimsel ayrıştırma sonuçlarını kullanabiliriz. Göstermek üzere seriyi, öngörü üretmekte kullanacağımız eğitim serisine ve doğruluğu sınayacağımız test serisine ayıran aşağıdaki işlevi kullanacağız.

def split_series(series, n=60):
    training = series.iloc[:-n]
    test = series.iloc[-n:]
    return training, test

n=60 olduğunda test serisi Temmuz 2019’dan başlayan beş yıllık dönemi kapsar.

training, test = split_series(nuclear)
test.index[0]

Timestamp('2019-07-31 00:00:00')

Haziran 2019’da olduğunuzu ve nükleer elektrik üretimi için beş yıllık öngörü istendiğini varsayın. Eğitim verileriyle bir model kurup öngörüler üreteceğiz. Eğitim verisinin mevsimsel ayrıştırmasıyla başlayalım.

decomposition = seasonal_decompose(training, model="additive", period=12)
trend = decomposition.trend

Şimdi eğilime doğrusal model uyduracağız. Açıklayıcı months değişkeni serinin başlangıcından itibaren geçen ay sayısıdır.

import statsmodels.formula.api as smf

months = np.arange(len(trend))
data = pd.DataFrame({"trend": trend, "months": months}).dropna()
results = smf.ols("trend ~ months", data=data).fit()

Sonuçların özeti şöyledir.

from thinkstats import display_summary

display_summary(results)
coef std err t P>|t| [0.025 0.975]
Intercept 6.482e+04 131.524 492.869 0.000 6.46e+04 6.51e+04
months 10.9886 1.044 10.530 0.000 8.931 13.046
R-kare: 0.3477

\(R^2\) yaklaşık 0,35’tir; modelin veriye pek iyi uymadığını düşündürür. Uyum doğrusunu çizerek bunu daha iyi görebiliriz. Eğitim ve test verilerinin beklenen değerlerini predict ile hesaplayacağız.

months = np.arange(len(training) + len(test))
df = pd.DataFrame({"months": months})
pred_trend = results.predict(df)
pred_trend.index = nuclear.index

Eğilim bileşeni ve doğrusal model şöyledir.

trend.plot(**trend_options)
pred_trend.plot(label="doğrusal model", **model_options)
decorate(ylabel="GWh")
_images/507c96c4ed96b5ce44d63a0179f27822bde64c12cc38ba024d6c87e3ff4c7691.png

Doğrusal modelin yakalayamadığı çok şey vardır; yine de genel olarak artan bir eğilim görülüyor.

Şimdi ayrıştırmanın mevsimsel bileşenini kullanarak aylık ortalamaların Series nesnesini hesaplayacağız.

seasonal = decomposition.seasonal
monthly_averages = seasonal.groupby(seasonal.index.month).mean()

Uyum doğrusundaki tarihleri monthly_averages içinde arayarak mevsimsel bileşeni öngörebiliriz.

pred_seasonal = monthly_averages[pred_trend.index.month]
pred_seasonal.index = pred_trend.index

Son olarak öngörüleri üretmek için mevsimsel bileşeni eğilime ekleyeceğiz.

pred = pred_trend + pred_seasonal

Eğitim verileri ve öngörüler şöyledir.

pred.plot(label="öngörü", **pred_options)
training.plot(label="eğitim", **actual_options)
decorate(ylabel="GWh")
_images/1cf5432dc7f5f54ad149485933665c4be33125c300a1c23b7972ad7ebfb9e42a.png

Öngörüler eğitim verisine makul ölçüde uyar; uzun dönemli eğilimin süreceği varsayımıyla tahmin de makul görünüyor.

Gelecekten geriye bakarak bu öngörünün ne kadar doğru çıktığına bakalım. Temmuz 2019’dan başlayan beş yıllık dönemin öngörülen ve gerçek değerleri şöyledir.

forecast = pred[test.index]
forecast.plot(label="öngörülen", **pred_options)
test.plot(label="gerçek", **actual_options)
decorate(ylabel="GWh")
_images/72be2a7cc7b4666143a06a7df2332f570598ba8a75e43161ca2856b92fa6f723.png

Öngörünün ilk yılı oldukça iyidir; ancak nükleer üretim 2020’de —muhtemelen COVID-19 salgını nedeniyle— beklenenden düşük kalmış ve uzun dönemli eğilime dönmemiştir.

Öngörü doğruluğunu nicelleştirmek için aşağıdaki işlevin hesapladığı ortalama mutlak yüzde hatayı (MAPE) kullanacağız.

def MAPE(predicted, actual):
    ape = np.abs(predicted - actual) / actual
    return np.mean(ape) * 100

Bu örnekte öngörüler ortalama %3,81 sapar.

MAPE(forecast, test)

np.float64(3.811940747879257)

Bölümün ilerleyen kısmında bu örneğe dönüp farklı modelle daha iyi sonuç alıp alamayacağımıza bakacağız.

12.4. Çarpımsal model#

Önceki bölümdeki toplamsal model, serinin uzun dönemli eğilim, mevsimsel bileşen ve artığın toplamı olduğunu varsayar; bu da mevsimsel bileşen ile artıkların büyüklüğünün zamanla değişmediği anlamına gelir.

Bu varsayımı ihlal eden örnek olarak 2014’ten bu yana küçük ölçekli güneş elektriği üretimine bakalım.

solar = elec["United States : small-scale solar photovoltaic"].dropna()
solar.plot(label="güneş", **actual_options)
decorate(ylabel="GWh")
_images/da02a2292738618b2a62bf85311cfda105cb3daa872fa54b4c2e3b94e4526633.png

Bu dönemde toplam üretim birkaç kat artmış, mevsimsel değişimin büyüklüğü de belirgin biçimde yükselmiştir.

Mevsimsel ve rastgele değişimin büyüklüklerinin eğilimle orantılı olduğunu varsayarsak serinin üç bileşenin çarpımı olduğu başka bir model elde ederiz.

Bu çarpımsal modeli denemek için seriyi eğitim ve test kümelerine ayıracağız.

training, test = split_series(solar)

Ardından model="multiplicative" argümanıyla seasonal_decompose çağıracağız.

decomposition = seasonal_decompose(training, model="multiplicative", period=12)

Sonuçlar şöyledir.

plot_decomposition(training, decomposition)
_images/43425f8eca9d148d98d9ec7c436e27620ac546ca486ed0b61c607236b3d8283f.png

Mevsimsel ve artık bileşenler artık çarpanlardır. Mevsimsel bileşen eğilimin yaklaşık %25 altı ile %25 üstü arasında değişir. İlk dönemdeki birkaç büyük değer dışında artık bileşen çoğunlukla iki yönde de %5’ten küçüktür. Bileşenleri şöyle çıkarabiliriz.

trend = decomposition.trend
seasonal = decomposition.seasonal
resid = decomposition.resid

Bu modelin \(R^2\) değeri çok yüksektir.

rsquared = 1 - resid.var() / training.var()
rsquared

np.float64(0.9999999992978134)

Güneş paneli üretimi büyük ölçüde maruz kaldığı güneş ışığına bağlıdır; üretimin yıllık döngüyü bu kadar yakından izlemesi anlamlıdır.

Uzun dönemli eğilimi öngörmek için karesel model kullanacağız.

months = range(len(training))
data = pd.DataFrame({"trend": trend, "months": months}).dropna()
results = smf.ols("trend ~ months + I(months**2)", data=data).fit()

Patsy formülündeki I(months**2) alt dizisi modele karesel terim ekler; ayrıca hesaplamamız gerekmez. Sonuçlar şöyledir.

display_summary(results)
coef std err t P>|t| [0.025 0.975]
Intercept 766.1962 13.494 56.782 0.000 739.106 793.286
months 22.2153 0.938 23.673 0.000 20.331 24.099
I(months ** 2) 0.1762 0.014 12.480 0.000 0.148 0.205
R-kare: 0.9983

Doğrusal ve karesel terimlerin p-değerleri çok küçüktür; karesel modelin eğilim hakkında doğrusal modelden fazla bilgi yakaladığını düşündürür. \(R^2\) de çok yüksektir.

Modelle geçmiş ve gelecekteki eğilimin beklenen değerini hesaplayabiliriz.

months = range(len(solar))
df = pd.DataFrame({"months": months})
pred_trend = results.predict(df)
pred_trend.index = solar.index

Şöyle görünür.

pred_trend.plot(label="karesel model", **model_options)
trend.plot(**trend_options)
decorate(ylabel="GWh")
_images/2f518204211146e2bef6687fdccc527bb26f7f5faea7be397ad6356ddd934d76.png

Karesel model geçmiş eğilime iyi uyar. Şimdi gelecekteki mevsimsel değişimi mevsimsel bileşenle öngörebiliriz.

monthly_averages = seasonal.groupby(seasonal.index.month).mean()
pred_seasonal = monthly_averages[pred_trend.index.month]
pred_seasonal.index = pred_trend.index

Son olarak geçmiş değerler için geriye dönük öngörüleri, gelecek için öngörüleri hesaplamak üzere eğilim ile mevsimsel bileşeni çarparız.

pred = pred_trend * pred_seasonal

Eğitim verileriyle birlikte sonuç şöyledir.

training.plot(label="eğitim", **actual_options)
pred.plot(label="öngörü", **pred_options)
decorate(ylabel="GWh")
_images/3f08f5fa2d4d9ebb6c1a03141ab950b70f5ee9a8c49ed0d576fab619c7209e6b.png

Geriye dönük öngörüler eğitim verisine iyi uyar ve gelecek öngörüleri makul görünüyor. Gerçekte ne kadar doğru çıktıklarına bakalım. Test verileriyle birlikte öngörüler şöyledir.

future = pred[test.index]
future.plot(label="öngörü", **pred_options)
test.plot(label="gerçek", **actual_options)
decorate(ylabel="GWh")
_images/05697e92e222cbdbf0baf2d421a31815453b6e20d93b7d7781a95e3e70d7ced0.png

İlk üç yıl öngörüler çok iyidir. Sonrasında gerçek büyüme beklentiyi aşmış gibi görünüyor.

Mevsimsel ayrıştırma güneş üretimini modelleyip öngörmede iyi, önceki nükleer üretim örneğinde ise etkisizdi. Sonraki bölümde başka bir yaklaşımı, otoregresyonu deneyeceğiz.

12.5. Otoregresyon#

Otoregresyonun ilk düşüncesi geleceğin geçmişe benzeyeceğidir. İncelediğimiz serilerde belirgin yıllık döngü vardır; gelecek haziran için öngörüde iyi başlangıç noktası geçen hazirandır.

Bunun ne kadar işe yaradığını görmek için nükleer üretim serisinde ardışık yılların aynı ayları arasındaki “yıldan yıla” farkları hesaplayalım.

diff = (nuclear - nuclear.shift(12)).dropna()
diff.plot(label="yıldan yıla farklar", **actual_options)
decorate(ylabel="GWh")
_images/43c4a7042b0939a7614ffb0035ff63bac75aa092eeb3b8a8ff63a11e6a37cfed.png

Bu farkların büyüklüğü özgün seriden çok küçüktür. Bu, otoregresyonun ikinci düşüncesini önerir: özgün değerler yerine farkları öngörmek daha kolay olabilir.

Fark dizisinin ardışık öğeleri arasında korelasyon olup olmadığına bakalım. Varsa geçmiş değerlerden geleceği öngörmekte kullanabiliriz.

Farkları ilk sütuna, 1, 2 ve 3 ay kaydırılmış sürümlerini sonraki sütunlara koyan bir DataFrame oluşturacağım. Diziler geciktirildiği için sütunların adları lag1, lag2, lag3’tür.

df_ar = pd.DataFrame({"diff": diff})
for lag in [1, 2, 3]:
    df_ar[f"lag{lag}"] = diff.shift(lag)

df_ar = df_ar.dropna()

Sütunlar arasındaki korelasyonlar şöyledir.

df_ar.corr()[["diff"]]
diff
diff 1.000000
lag1 0.562212
lag2 0.292454
lag3 0.222228

Bunlara gecikmeli korelasyon ya da otokorelasyon denir; “oto” öneki seriyi kendisiyle korelasyonladığımızı gösterir. diff ile lag1 korelasyonuna, dizideki ardışık öğelerin korelasyonu olduğu için seri korelasyon denir.

Korelasyonlar öngörüye yardımcı olabilecek kadar güçlüdür; çoklu regresyona koyalım. Aşağıdaki işlev ilk sütunu yanıt, diğerlerini açıklayıcı değişken kabul eden Patsy formülü oluşturur.

def make_formula(df):
    """Sütun adlarından Patsy formülü oluşturur."""
    y = df.columns[0]
    xs = " + ".join(df.columns[1:])
    return f"{y} ~ {xs}"

Önceki üç değerden sıradaki değeri öngören doğrusal modelin sonuçları şöyledir.

formula = make_formula(df_ar)
results_ar = smf.ols(formula=formula, data=df_ar).fit()
display_summary(results_ar)
coef std err t P>|t| [0.025 0.975]
Intercept 24.2674 114.674 0.212 0.833 -201.528 250.063
lag1 0.5847 0.061 9.528 0.000 0.464 0.706
lag2 -0.0908 0.071 -1.277 0.203 -0.231 0.049
lag3 0.1026 0.062 1.666 0.097 -0.019 0.224
R-kare: 0.3239

Şimdi predict ile geçmiş değerler için geriye dönük öngörüler üretebiliriz. Verilerle karşılaştırması şöyledir.

pred_ar = results_ar.predict(df_ar)
pred_ar.plot(label="öngörüler", **pred_options)
diff.plot(label="farklar", **actual_options)
decorate(ylabel="GWh")
_images/63852a07067f2b2bbe06765dfd94b189148e5f46be16def5c29b4b62218ceaa9.png

Öngörüler bazı yerlerde iyidir; ancak \(R^2\) yalnızca yaklaşık 0,319’dur ve iyileştirme payı vardır.

resid_ar = (diff - pred_ar).dropna()
R2 = 1 - resid_ar.var() / diff.var()
R2

np.float64(0.3190252265690783)

Öngörüleri iyileştirmenin bir yolu bu modelin artıklarını hesaplayıp artıkları öngörmek için başka bir model kullanmaktır; bu, otoregresyonun üçüncü düşüncesidir.

12.6. Hareketli ortalama#

Haziran 2019’da olduğunuzu ve Haziran 2020 için öngörü istendiğini varsayın. İlk tahmininiz bu yılki değerin gelecek yıl yineleneceği olabilir.

Mayıs 2020’de tahmini güncellemeniz istenirse son üç ayın sonuçlarını ve önceki bölümdeki otokorelasyon modelini kullanarak yıldan yıla farkı öngörebilirsiniz.

Son birkaç ayın öngörülerinin sürekli düşük kaldığını görürseniz gelecek ayınkini yukarı yönlü düzeltebilirsiniz. Temel varsayım, yakın geçmişteki öngörü hatalarının gelecektekileri öngörmesidir.

Bunu sınamak için ilk sütunda otoregresyon artıklarını, diğerlerinde gecikmeli sürümlerini içeren bir DataFrame oluşturabiliriz. Burada 1 ve 6 aylık gecikmeler kullanacağım.

df_ma = pd.DataFrame({"resid": resid_ar})

for lag in [1, 6]:
    df_ma[f"lag{lag}"] = resid_ar.shift(lag)

df_ma = df_ma.dropna()

Artıkların otoregresyon modelini ols ile kurabiliriz. Bu bölüme, öngörü değişkenliğini hareketli ortalamaya benzer biçimde azalttığı için “hareketli ortalama” denir. Terimi çok açıklayıcı bulmasam da gelenekseldir.

Artık modelinin özeti şöyledir.

formula = make_formula(df_ma)
results_ma = smf.ols(formula=formula, data=df_ma).fit()
display_summary(results_ma)
coef std err t P>|t| [0.025 0.975]
Intercept -14.0016 114.697 -0.122 0.903 -239.863 211.860
lag1 0.0014 0.062 0.023 0.982 -0.120 0.123
lag6 -0.1592 0.063 -2.547 0.011 -0.282 -0.036
R-kare: 0.0247

\(R^2\) oldukça küçüktür; bu terim pek yardımcı olmayacak gibi görünüyor. Ancak 6 aylık gecikmenin p-değeri küçüktür; şanstan beklenenden fazla bilgi sağladığını düşündürür.

Modelle artıklar için geriye dönük öngörüler üretebiliriz.

pred_ma = results_ma.predict(df_ma)

Yıldan yıla farkların geriye dönük öngörülerini üretmek için ikinci modelin düzeltmesini ilk modelin sonuçlarına ekleriz.

pred_diff = pred_ar + pred_ma

İki model toplamının \(R^2\) değeri yaklaşık 0,332’dir; hareketli ortalama düzeltmesi olmadan elde edilen 0,319’dan biraz iyidir.

resid_ma = (diff - pred_diff).dropna()
R2 = 1 - resid_ma.var() / diff.var()
R2

np.float64(0.3315101001391231)

Şimdi bu farkları özgün değerlerin geriye dönük öngörülerini üretmekte kullanacağız.

12.7. Otoregresyonla geriye dönük öngörü#

Yıldan yıla farkları özgün serinin indisiyle hizalı bir Series içine koyarak başlayacağız.

pred_diff = pd.Series(pred_diff, index=nuclear.index)

NaN değerlerini isna ile denetlediğimizde yeni Series’in ilk 21 öğesinin eksik olduğunu görürüz.

n_missing = pred_diff.isna().sum()
n_missing

np.int64(21)

Çünkü yıldan yıla fark için seriyi 12 ay, ilk otoregresyon için farkları 3 ay, ikinci model için ilk model artıklarını 6 ay kaydırdık. Her kaydırmada baştan birkaç değer kaybederiz; toplam 21’dir.

Geriye dönük öngörü üretmeden önce özgün serinin ilk 21 öğesini yeni bir Series içine kopyalayarak başlangıç koşullarını sağlamalıyız.

pred_series = pd.Series(index=nuclear.index, dtype=float)
pred_series.iloc[:n_missing] = nuclear.iloc[:n_missing]

İndis 21’den sonuna kadar öğeleri dolduran aşağıdaki döngüyü çalıştırabiliriz. Her öğe önceki yılın değeri ile öngörülen yıldan yıla farkın toplamıdır.

for i in range(n_missing, len(pred_series)):
    pred_series.iloc[i] = pred_series.iloc[i - 12] + pred_diff.iloc[i]

İlk 21 değeri kusursuz “öngörmüş” sayılmamak için kopyaladığımız öğeleri yeniden NaN ile değiştireceğiz.

pred_series[:n_missing] = np.nan

Özgün seriyle karşılaştırıldığında geriye dönük öngörüler şöyledir.

pred_series.plot(label="öngörülen", **pred_options)
nuclear.plot(label="gerçek", **actual_options)
decorate(ylabel="GWh")
_images/cf019402afdf59c25db7111eefa04705bf020d436b262424b6f7e1919c67f127.png

Oldukça iyi görünüyorlar; \(R^2\) yaklaşık 0,86’dır.

resid = (nuclear - pred_series).dropna()
R2 = 1 - resid.var() / nuclear.var()
R2

np.float64(0.8586566911201015)

Bu öngörüleri hesaplamakta kullandığımız modele SARIMA denir; ARIMA ailesinin bir üyesidir. Kısaltmadaki parçalar model öğelerini gösterir.

  • S, ilk adımda bir mevsimsel dönem aralıklı değerlerin farkını aldığımız için seasonal (mevsimsel).

  • AR, farklardaki gecikmeli korelasyonları modelleyen autoregression (otoregresyon).

  • I, pred_series hesaplamasındaki yinelemeli süreç kalkülüsteki integrasyona benzediği için integrated (bütünleşik).

  • MA, ilk modelin artıklarıyla kurduğumuz ikinci otoregresyonun geleneksel adı moving average (hareketli ortalama).

ARIMA modelleri zaman serileri için güçlü ve esnek araçlardır.

12.8. ARIMA#

StatsModels’ın “time series analysis” (zaman serisi analizi) kısaltması olan tsa kitaplığı, ARIMA modellerini uydurup öngörü üreten ARIMA işlevini içerir.

Önceki bölümlerde geliştirdiğimiz SARIMA modelini uydurmak için işlevi order ve seasonal_order adlı iki demetle çağıracağız. Önceki modele karşılık gelen order şöyledir.

order = ([1, 2, 3], 0, [1, 6])

order değerleri şunları belirtir:

  • AR modeline eklenecek gecikmeler; burada ilk üçü.

  • Ardışık öğeler arasında kaç kez fark alınacağı; mevsimsel fark kullandığımız için burada 0.

  • MA modeline eklenecek gecikmeler; burada birinci ve altıncı.

seasonal_order değerleri şöyledir.

seasonal_order = (0, 1, 0, 12)

İlk ve üçüncü öğelerin 0 olması mevsimsel AR ya da MA bulunmadığını gösterir. İkinci öğe 1’dir; mevsimsel fark hesaplandığını belirtir. Son öğe mevsimsel dönemdir.

Modeli ARIMA ile şöyle oluşturup uydururuz.

import statsmodels.tsa.api as tsa

model = tsa.ARIMA(nuclear, order=order, seasonal_order=seasonal_order)
results_arima = model.fit()
display_summary(results_arima)
coef std err z P>|z| [0.025 0.975]
ar.L1 0.0458 0.379 0.121 0.904 -0.697 0.788
ar.L2 -0.0035 0.116 -0.030 0.976 -0.230 0.223
ar.L3 0.0375 0.049 0.769 0.442 -0.058 0.133
ma.L1 0.2154 0.382 0.564 0.573 -0.533 0.964
ma.L6 -0.0672 0.019 -3.500 0.000 -0.105 -0.030
sigma2 3.473e+06 1.9e-07 1.83e+13 0.000 3.47e+06 3.47e+06

Sonuçlar AR modelindeki üç, MA modelindeki iki gecikmenin tahmini katsayılarını ve artıkların varyansı sigma2’yi içerir.

results_arima içinden geriye dönük öngörüleri içeren fittedvalues çıkarılabilir. Kendi sonuçlarımızın başında eksik değerler bulunmasıyla aynı nedenle bunun başında da yanlış değerler vardır; bunları çıkaracağız.

fittedvalues = results_arima.fittedvalues[n_missing:]

Uydurulan değerler kendi hesapladıklarımıza benzer, ancak aynı değildir; muhtemelen ARIMA başlangıç koşullarını farklı işler.

fittedvalues.plot(label="ARIMA modeli", **pred_options)
nuclear.plot(label="gerçek", **actual_options)
decorate(ylabel="GWh")
_images/f13c890c5278081a05e133943b2aa45ea6e255c53fcfd6706467d45b8ce07440.png

\(R^2\) de benzer, ancak tam olarak aynı değildir.

resid = fittedvalues - nuclear
R2 = 1 - resid.var() / nuclear.var()
R2

np.float64(0.8262717330814415)

ARIMA işlevi modelin farklı sürümlerini denemeyi kolaylaştırır.

Alıştırma olarak order ve seasonal_order içindeki farklı değerleri deneyip daha yüksek \(R^2\) değerli model bulmaya çalışın.

12.9. ARIMA ile öngörü#

ARIMA sonucunun get_forecast yöntemi öngörü üretir. Göstermek için seriyi eğitim ve test kümelerine ayırıp aynı modeli eğitim kümesine uyduracağız.

training, test = split_series(nuclear)
model = tsa.ARIMA(training, order=order, seasonal_order=seasonal_order)
results_training = model.fit()

Sonuçla test kümesi için öngörü üretebiliriz.

forecast = results_training.get_forecast(steps=len(test))

Sonuç, forecast_mean özniteliği ile güven aralığı döndüren bir işlev içeren nesnedir.

forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int()
forecast_ci.columns = ["lower", "upper"]

Sonuçları şöyle çizip gerçek zaman serisiyle karşılaştırabiliriz.

plt.fill_between(
    forecast_ci.index,
    forecast_ci.lower,
    forecast_ci.upper,
    lw=0,
    color="gray",
    alpha=0.2,
)
plt.plot(forecast_mean.index, forecast_mean, label="öngörü", **pred_options)
plt.plot(test.index, test, label="gerçek", **actual_options)
decorate(ylabel="GWh")
_images/d47484804621e7034e7f0c331598ee238646d70acc170a52990198754bc8fd58.png

Gerçek değerler öngörülerin güven aralığında neredeyse bütünüyle kalır. Öngörülerin MAPE değeri şöyledir.

MAPE(forecast_mean, test)

np.float64(3.3817549247338285)

Öngörüler ortalama %3,38 sapar; mevsimsel ayrıştırmayla elde edilen %3,81’den biraz iyidir.

ARIMA mevsimsel ayrıştırmadan daha esnektir ve çoğu zaman daha iyi öngörüler üretir. Bu seride otokorelasyonlar özellikle güçlü olmadığından üstünlüğü sınırlıdır.

12.10. Sözlük#

  • zaman serisi (time series): Her değerin belirli bir zamanla ilişkili olduğu, çoğu zaman düzenli aralıklarla alınmış ölçümleri temsil eden veri kümesi.

  • mevsimsel ayrıştırma (seasonal decomposition): Zaman serisini uzun dönemli eğilim, yinelenen mevsimsel bileşen ve artık bileşenine ayırma yöntemi.

  • eğitim serisi (training series): Modeli uydurmakta kullanılan zaman serisi bölümü.

  • test serisi (test series): Model öngörülerinin doğruluğunu denetlemekte kullanılan bölüm.

  • geriye dönük öngörü (retrodiction): Geçmişte gözlenmiş bir değer için, çoğunlukla modeli sınamak veya doğrulamak amacıyla üretilen öngörü.

  • pencere (window): Hareketli ortalama hesaplamak için kullanılan ardışık değerler dizisi.

  • hareketli ortalama (moving average): Dalgalanmaları yumuşatmak için örtüşen pencerelerdeki değerlerin ortalamasıyla hesaplanan seri.

  • seri korelasyon (serial correlation): Zaman serisinin ardışık öğeleri arasındaki korelasyon.

  • otokorelasyon (autocorrelation): Zaman serisi ile kendisinin kaydırılmış ya da geciktirilmiş sürümü arasındaki korelasyon.

  • gecikme (lag): Seri korelasyon ya da otokorelasyondaki kaydırma miktarı.

12.11. Alıştırmalar#

12.11.1. Alıştırma 12.1#

Mevsimsel ayrıştırma örneği olarak ABD’deki aylık ortalama yüzey sıcaklıklarını modelleyelim. Our World in Data’nın veri kümesi, 1950–2024 arasında dünyanın çoğu ülkesi için “kara, deniz ve iç su yüzeylerini kapsayan, yerden 2 metre yüksekte ölçülmüş hava sıcaklığını [Santigrat]” içerir. İndirme yönergeleri bu bölümün notebook’unda yer alır.

# Aşağıdaki hücre, Our World in Data'nın hazırladığı verileri indirir.
# Verileri 18 Eylül 2024'te şu adresten indirdim:
# https://ourworldindata.org/grapher/average-monthly-surface-temperature

# Kaynak, Copernicus Climate Change Service'in (2019) değiştirilmiş verilerine dayanır
# ve Our World in Data tarafından “kapsamlı biçimde işlenmiştir”.

filename = "monthly-average-surface-temperatures-by-year.csv"
download("https://github.com/TerekliTahaBerk/thinkstatstr/raw/v3/data/" + filename)

Verileri şöyle okuyabiliriz.

temp = pd.read_csv("monthly-average-surface-temperatures-by-year.csv")

temp.head()
Entity Code Year 2024 2023 2022 2021 2020 2019 2018 ... 1959 1958 1956 1954 1952 1957 1955 1953 1951 1950
0 Afghanistan AFG 1 3.300064 -4.335608 -0.322859 -1.001608 -2.560545 0.585145 1.042471 ... -2.333814 0.576404 -3.351925 -2.276692 -2.812619 -4.239172 -2.191683 -2.915993 -3.126317 -2.655707
1 Afghanistan AFG 2 1.024550 4.187041 2.165870 5.688000 2.880046 0.068664 3.622793 ... -1.545529 0.264962 0.455350 -0.304205 0.798226 -2.747945 1.999074 1.983414 -2.642800 -3.996040
2 Afghanistan AFG 3 5.843506 10.105444 10.483686 9.777976 6.916731 5.758049 10.794412 ... 5.942937 7.716459 5.090270 4.357703 4.796146 4.434027 7.066073 4.590406 3.054388 3.491112
3 Afghanistan AFG 4 11.627398 14.277164 17.227650 15.168276 12.686832 13.838840 14.321226 ... 13.752827 14.712909 11.982360 12.155265 13.119270 8.263829 10.418768 11.087193 9.682878 8.332797
4 Afghanistan AFG 5 18.957850 19.078170 19.962734 19.885902 18.884047 18.461287 18.100782 ... 17.388723 16.352045 20.125462 18.432117 17.614851 15.505956 15.599709 17.865084 17.095737 17.329062

5 rows × 78 columns

Aşağıdaki hücre ABD’nin 2001’den serinin sonuna kadarki verilerini seçip Pandas Series içine yerleştirir.

temp_us = temp.query("Code == 'USA'")
columns = [str(year) for year in range(2000, 2025)]
temp_series = temp_us.loc[:, columns].transpose().stack()
temp_series.index = pd.date_range(start="2000-01", periods=len(temp_series), freq="ME")

Şöyle görünür.

temp_series.plot(label="aylık ortalama", **actual_options)
decorate(ylabel="Yüzey sıcaklığı (°C)")
_images/9ca56037fc93340f7885800eced4805469d6a5c5a1e70168a3e14dd552938e40.png

Şaşırtıcı olmayacak biçimde güçlü bir mevsimsel örüntü vardır. 12 aylık dönemle toplamsal mevsimsel ayrıştırma hesaplayın. Eğilime doğrusal model uydurun. Bu dönemde yüzey sıcaklığındaki ortalama yıllık artış nedir? Merak ederseniz analizi başka dönemlerle ya da ülkelerle yineleyin.

12.11.2. Alıştırma 12.2#

Bu bölümde 2014–2019 arasında küçük ölçekli güneş elektriği üretimini modellemek ve 2019–2024 üretimini öngörmek için çarpımsal ayrıştırma kullandık. Şimdi aynısını şebeke ölçekli güneş enerjisi için yapalım. Zaman serisi şöyledir.

util_solar = elec["United States : all utility-scale solar"].dropna()
util_solar = util_solar[util_solar.index.year >= 2014]
util_solar.plot(**actual_options)
decorate(ylabel="GWh")
_images/0bbb981a15cc7921a7aa00224026ba928aa26e6a0449bea2599a4fea5f0783f4.png

Veriyi eğitim ve test serilerine ayırmak için split_series kullanın. Eğitim serisini 12 aylık dönemle çarpımsal olarak ayrıştırın. Eğilime doğrusal ya da karesel model uydurup mevsimsel bileşeni de içeren beş yıllık öngörü üretin. Öngörüyü test serisiyle birlikte çizip ortalama mutlak yüzde hatayı (MAPE) hesaplayın.

12.11.3. Alıştırma 12.3#

ARIMA modelinin ABD’deki hidroelektrik üretimine ne kadar iyi uyduğuna bakalım. 2001–2024 zaman serisi şöyledir.

hydro = elec["United States : conventional hydroelectric"]
hydro.plot(**actual_options)
decorate(ylabel="GWh")
_images/55f57886cecc8ba0104561fccba921b757649aa8df2ecb336bbdc51b2788d751.png

Bu verilere 12 aylık mevsimsel dönemli SARIMA modeli uydurun. Otoregresyon ve hareketli ortalama bölümlerinde farklı gecikmeler deneyip \(R^2\) değerini en yükseğe çıkaran birleşimi bulmaya çalışın. Beş yıllık öngörü üretip güven aralığıyla birlikte çizin.

NOT: Modele eklediğiniz gecikmelere bağlı olarak uydurulan değerlerin ilk 12–24 öğesi güvenilir olmayabilir. Çizmeden ya da \(R^2\) hesaplamadan önce bunları çıkarmak isteyebilirsiniz.

Think Stats: Python ile Keşifsel Veri Analizi, 3. Baskı

Telif hakkı 2024 Allen B. Downey

Kod lisansı: MIT Lisansı

Metin lisansı: Creative Commons Atıf-Gayriticari-AynıLisanslaPaylaş 4.0 Uluslararası