R ile Enerji Analizi – Bölüm 7- Üssel Düzgünleştirme, Mevsimsellik ve Trend (ETS, STL) ile Elektrik Talep Tahmini Denemesi

 

Özet: Bu bölümde, 6. Bölümün devamı şeklinde teorisine fazlaca girmeden, Rob J. Hyndman tarafından hazırlanan Forecast paketi ile bazı tahmin modellerini inceleyeceğiz. Bu bölümde özellikle üssel düzgünleştirme (exponential smoothing), Loess ile mevsimsel ve trend ayrıştırma yöntemi (STL, Seasonal and Trend decomposition using Loess) ile otomatik model tahminlerinden sonra “forecast(“ komutu ile modellerimizi ileri taşıyarak, kestirim yapmaya çalışacağız.

Bu Bölümdeki Fonksiyonlar

read.csv

ts(veri, start=başladığı tarih, frequency= periyodu)

decompose(veri seti)

stl(veri seti, parametreler)

ses(veri seti)

forecast(model, h=tahmin miktarı)

hw(veri seti)

stlf(veri seti)

 

 

7.0. Veriyi yüklemek

Her zamanki gibi veri dosyamızı tekrar yükleyelim.

veri=read.csv("http://www.barissanli.com/calismalar/dersler/r/elektrik-talep.csv", header=TRUE,sep=";",dec=".")

Ortalama tüketim verisini de ot değişkenine atayalım.

ot=veri$ortalama_tuketim

Þimdi “veri” yazdığımız her noktada alt detaylara erişebiliriz. Ders 6’da olduğu gibi forecast paketini bilgisayarınıza kurduğunuza emin olun. Daha önce kurmadıysanız:

install.packages(“forecast”)

Ä°le kurabilirsiniz.

7.1. Zaman Serileri

Eğer zaman serileri kullanacak isek, önce yukarıda yüklediğimiz ot değişkenini bir zaman serisi olarak belirtmekte fayda var. Bu ise gayet basit bir şekilde yapılabiliyor. Bunun için “ts” komutunu kullanacağız.

library("forecast")

ot=veri$ortalama_tuketim

ott=ts(ot)

 

Veri her ne kadar zaman serisi olarak basitçe tanımlanıyor olsa da, zaman serisini biraz daha detaylı tanımlayabiliriz. Hatırlanacağı üzere “veri” değişkenimiz 1 Ocak 2007’den başlayarak, günlük olarak devam eden bir veri setiydi. Bu yüzden zaman serimizi biraz daha detaylı tanımlayalım. Bunun için, 2007 başından başladığı için “start=2007”, 365 günde bir devrettiği için de “frequency=365” parametrelerini ekleyelim.

ott=ts(ot,start=2007,frequency=365)

plot(ott)

Bir seriyi zaman serisi olarak tanımladıktan sonra, çizdirdiğimizde otomatik olarak x ekseninin yıllara dönüştüğünü görebilirsiniz.

Artık zaman serimiz hazır olduğuna göre çalışmaya başlayabiliriz.

7.2. Seriyi Bileşenlerine Ayırmak – Decompose

(Konuyla ilgili açıklamaları https://www.otexts.org/fpp/6/1 adresinden bulabilirsiniz.)

Daha önceki bölümlerde decompose ile bir veriyi bileşenlerine ayırmıştık.

plot(decompose(ott))

 

 

Þimdi ise bir diğer yöntem olan stl’dir. “Seasonal and Trend decomposition using Loess” baş harflerinden oluşan STL yönteminde Loess ise doğrusal olmayan ilişkileri hesaplama yöntemidir.

Bir de stl ile deneyelim.

plot(stl(ott,s.window = "periodic"))

 

Tabii, daha başka yöntemler de var. Fakat şimdilik bu iki komut işimizi görecektir. Fakat dikkat edin, Türkiye elektrik talep serisi hem 365 günlük bir periyoda hem de 354 günlük “hicri” bir periyorda daha sahiptir.

7.3. Üssel Düzgünleştirme – Exponential Smoothing

(Konuyla ilgili açıklamaları https://www.otexts.org/fpp/7/1 adresinden bulabilirsiniz.)

Exponential SmoothingÂ’de basit bir exponential smoothing deneyecek isek,

fit1<-ses(ott)

 

yeterli olacaktır. Hesaplanan modeli ileri götürmek/tahmin yapmak için ise “forecast” fonksiyonunu kullanırız. Burada forecast fonksiyonunda “h=10” değeri ile de 10 adım ileriye kestirim yaparız..

 

Çizdirmek için ise, “plot” alanını sınırlandırmamız gerekecek.

Bu basit exponential smoothing yöntemi çok başarılı gözükmüyor olabilir. Çünkü mevsimselliğin olduğu yerlerde çok da başarılı değil, fakat daha kısa zaman aralıkları var ise o zaman daha başarılı olabilir. Bu yüzden bir de “ottw” yani haftalık bir değişken tanımlayalım. (Burada 24’ten büyük periyodlarda sorun olmaktadır)

ottw=ts(ott,frequency=7)

 

Þimdi kestirim yapmayı deneyelim.

ottw=ts(ott,frequency=7)

fit2<-ses(ottw)

plot(forecast(fit2,h=300),xlim=c(255,265))

Dikkat ederseniz fonksiyonumuzu yıllıktan haftalığa çevirdiğimizde x ekseni de değişmektedir.

Þimdi bir de Holt Winters yöntemini deneyelim. Yine haftalık frekanslı verimizi kullanıyoruz. Bunun için “hw” komutunu kullanacağız.

fit3<-hw(ottw)

plot(forecast(fit3,h=300),xlim=c(255,265))

Başlıktan da anlayabileceğiniz üzere Holt Winters’in “additive” yöntemi ile bir grafik oluştu. Fakat Holt Winters’i “multiplicative” denemek istersek, 

seasonal=”multiplicative”

parametresini ekleriz.

fit3<-hw(ottw,seasonal = "multiplicative")

plot(forecast(fit3,h=300),xlim=c(255,265))

 

7.4. Üssel Düzgünleştirme Modelinin Otomatik Belirlenmesi

Modeller, testler derken kendimizi bir sürü karmaşanın içinde bulabiliriz. “forecast” paketindeki bir komut tüm bunları otomatik hale getirmekte ve bizim yerimize en başarılı modeli seçmektedir.

Önce modeli hesaplatalım

fit4<-ets(ottw)

plot(forecast(fit4,h=300))

Biraz daha yakından bakalım:

plot(forecast(fit4,h=300),xlim=c(258,285))

“ets” komutunun otomatik olarak seçtiği parametreler, grafiğin başlığında gözükmektedir. Uzun sezonsal eğilimleri olan veriler için uygun gözükmese de, kısa periyotlu verilerde istenene yakın sonuçlar vermektedir.

7.4. STL – Seasonal and Trend decomposition using Loess

Yukarıda exponential smoothing’de görüldüğü üzere uzun periyotlu mevsimsel trendlere pek de uygun olmadığı durumlar var. Bunun yerine yine “forecast” paketindeki stl ve stlf komutlarını denemekte fayda var.

Veri olarak 365 periyotlu ott verisini kullanacağız. Burada önceki bölümde kullandığımız stl yerine stl forecast olarak stlf’i kullanacağız. “stlf”’e sadece veri setini vererek sonucu gözlemleyeceğiz.

fit5<-stlf(ott)

plot(forecast(fit5,h=300))

 

Görüldüğü üzere stlf komutu ile çok daha başarılı bir forecast elde etmişe benziyoruz.

Peki bu kestirim sonuçlarını nasıl göreceğiz. Bunun için modeli kestirdikten sonra forecast komutu ile başlayan kısmı plot’a aktarmak yerine bir başka değişkene atalım. Sonra bu verinin içeriğine bakalım.

verif=forecast(fit5,h=300)

str(verif)

Ekrana tüm alt değişkenler sığmadığından yer veremediğimiz, daha altlarda bir “mean” değişkeni var.

Eğer bu değişkeni çağırırsak da sadece hesaplanan tahmini görebiliriz.

plot(verif$mean)

Verilerimizin:

-          Kullandığımız deÄŸiÅŸkenleri xÂ’de yani verif$xÂ’de

-          Hesaplanan kestirim deÄŸerleri de meanÂ’de yani verif$meanÂ’de

İki seriyi birleştirirmek için de sadece yeni bir vektör tanımlayarak iki veriyi arka arkaya virgülle sıralıyoruz.

tumveri=c(verif$x,verif$mean)

plot (tumveri,type="l")

Görüldüğü üzere tüm bir seriyi bir araya getirdiğimizde kestirimi yapılan kısmın, tarihsel veriden daha az varyansı olduğu düşünülebilir.

Þimdi tümveriyi bir zaman serisi olarak 2007’den başlayan ve günlük değişen bir veri olarak değerlendirerek bir kez daha çizelim.

plot (ts(tumveri,start=2007,frequency=365),type="l")

 

Elde ettiğimiz bu veriyi şimdi bir kez daha bileşenlerine ayıralım, bakalım matematiksel hesaplamalar nasıl bir trend ve mevsimsellik hesabı yapmış.

tumverit=ts(tumveri,start=2007,frequency=365)

plot(decompose(tumverit))

 

Beklediğimiz gibi modelin 2012’den önceki gerçekleşmelerindeki mevsimsellik aynen alınarak, trend aynı şekilde devam ettirilmiş. Fakat random yani gerçek verideki hesaplanamayan kısım yerine hemen hemen hiçbirşey eklenememiş ve dolayısıyla model kestirimi daha “sade” olmuştur.

 

Sonuç

Bu bölümde ise exponential smoothing ve mevsimse ve trend yöntemleri incelenmiştir. Modelimiz görüldüğü kadarı ile daha stabil gözükmekle birlikte, “forecast” paketinin bize sağladığı en önemli kolaylık modelin otomatik seçilmesini sağlayan ets() ve stlf() komutları olmuştur.

Eğer fonksiyonların sonuçlarını değiştirmek isterseniz, çıktıları bir değişkene alarak alt değişkenlerden değişiklikler yapabilirsiniz. Mesela trend’ini değiştirme imkanınız var.

Yine içeride geçen bir çok kısaltmayı uzun uzudayı anlatmadık. Bunların hepsi Hyndman’ın kitabında https://www.otexts.org/fpp/7/6 yer almaktadır.

Özellikle “7.6 A taxonomy of exponential smoothing methods”’da gördüğümüz A,Ad,N ne demek açıklamaları bulabilirsiniz.

Tahmin yaparken en önemli konulardan biri de “matematiksel modellere körü körüne saplanmamaktadır”. Sonunda en iyi tahminler uzman görüşü + iyi matematiksel modellerle oluşturulmaktadır.

 

Komutlar

veri=read.csv("http://www.barissanli.com/calismalar/dersler/r/elektrik-talep.csv", header=TRUE,sep=";",dec=".")

 

library("forecast")

ot=veri$ortalama_tuketim

ott=ts(ot)

ot_t=ts(ot,start=2007,frequency=365)

 

plot(ott)

 

plot(decompose(ott))

plot(stl(ott,s.window = "periodic"))

 

fit1<-ses(ott)

plot(forecast(fit1,h=300))

 

fit1<-ses(ott)

plot(forecast(fit1,h=300))

forecast(fit1,h = 10)

 

plot(forecast(fit1,h=300),xlim=c(2011.9,2012.2))

 

ottw=ts(ott,frequency=7)

fit2<-ses(ottw)

plot(forecast(fit2,h=300),xlim=c(255,265))

 

fit3<-hw(ottw)

plot(forecast(fit3,h=300),xlim=c(255,265))

 

fit3<-hw(ottw,seasonal = "multiplicative")

plot(forecast(fit3,h=300),xlim=c(255,265))

 

fit4<-ets(ottw)

plot(forecast(fit4,h=300))

 

fit4<-ets(ottw)

plot(forecast(fit4,h=300),xlim=c(258,285))

 

fit5<-stlf(ott)

plot(forecast(fit5,h=300))

 

verif=forecast(fit5,h=300)

str(verif)

plot(verif$mean)

 

tumveri=c(verif$x,verif$mean)

plot (tumveri,type="l")

 

plot (ts(tumveri,start=2007,frequency=365),type="l")

 

tumverit=ts(tumveri,start=2007,frequency=365)

plot(decompose(tumverit))