R ile Enerji Analizi – Bölüm 6- doğrusal Regresyon ile Elektrik Talep Tahmini Denemesi

R ile Enerji Analizi – Bölüm 6- Doğrusal Regresyon ile Elektrik Talep Tahmini Denemesi

Özet: Bu bölümde, teorisine fazlaca girmeden, Rob J. Hyndman tarafından hazırlanan Forecast paketi ile bazı tahmin modellerini inceleyeceğiz, eğer teorik detaylarına girmek isterseniz, yazının içindeki linkleri takip edebilirsiniz.

Doğrusal regresyon modeli kurmadan önce bir çok test yapılması gerekir. Bu testler R’da ve ekonometri kitaplarında mevcut. Bu kısımda daha çok “lm” fonksiyonu ve özelliklerine değineceğiz.

Bu Bölümdeki Fonksiyonlar

read.csv

rep(değişken, tekrar sayısı)

gecik<-function(x, k){ c(rep(NA,k),x[1: ( length(x)-k)]) }

diff()

lm()

AIC()
6.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=".")

Şimdi “veri” yazdığımız her noktada alt detaylara erişebiliriz.

6.1. Forecast paketi

Zaman serileri ve modelleri, Türkiye’de üniversitelerde gözlemlediğim kadarı ile roket bilimine yakın bir zorluktaymış gibi öğretiliyor. Çok fazla gramerden pratiğe pek vakit kalmadığı gibi, dersleri bitirenler arasında modelleme yapabilen ise pek de fazla olmuyor.

Bu yüzden bu bölümde Hyndman’ın “forecast” paketini kullanacağız. Bu paket ile ilgili gerek Hyndman’ın sayfalarında, gerek internette bir sürü örnek bulabilirsiniz. Ama bu paketi kullanmamızın sebepleri:

• 
  Bir çok metodu kullanma imkanı tanıması (Exponential Smoothing, ARIMA, Neural Networks)
  
• 
  Otomatik olarak modeli belirlemesi
  
• 
  Kullanımının basit olması
  
• 
  Farklı frekansların bir arada olduğu zaman serilerini kullanma imkanı vermesi
  

Sayılabilir. Ayrıca Hyndman’ın videoları da internetten erişilebilir.

Bakınız linkler: http://robjhyndman.com/ , “Forecasting: principles and practice (Açık Kitap)” https://www.otexts.org/fpp/ ,

Önce “forecast” paketini R Studio’dan yükleyelim.

install.packages("forecast")

Eğer eksik bir paket belirtirse ki normalde tüm paketleri de yükler: install.packages(“zoo”) yazarak forecast’in kullandığı “zoo” paketini de kullanabilirsiniz.

Artık “forecast” paketini kullanmaya hazırız.

library("forecast")


6.2. Konuya Girmeden Önce – Fonksiyon Tanımlaması- Zaman Gecikmesi (Lag)

Özellikle doğrusal regresyon çalışanlar için en çok kullanılan iki konudan biri zaman gecikmesi, diğeri de farklar konusudur.

R’da zaman gecikmesi “lag”ler konusunda internet gruplarında da farklı tartışmalar var. Ben o yüzden kolayca lag oluşturan ama sonra değişik sorunlara maruz kalınan komutların bir kısmını yazayım.

• 
  lag(zaman_serisi, -gecikme_suresi, na.pad=TRUE) : burada bir zaman_serisi, gecikme_suresi kadar geciktirilirken, boş hücreler “na.pad=TRUE” komutu gereğince NA, R’da tanımlı Not Available sabiti ile doldurulur.
  
• 
  “dynlm” paketini yüklerseniz, dynlm komutu içinde L(zaman_serisi, -gecikme_suresi) komutunu da kullanabilirsiniz. (tabii dynlm( komutu içinde)
  

Mantığını oluşturmak açısında ben size daha direkt bir yöntemi göstermek isterim. Bir veri serisine gecikme eklediğimizde pratikte ne yapıyoruz? Veriyi gecikme suresi kadar aşağı kaydırıyoruz, verinin en başındaki yerlere varsa geçmiş veriyi yoksa “NA” doldurup, verinin en altındaki(en yeni tarihli veri setinde de) gecikme suresi kadar veri atıyoruz.

Şimdi bunu bir seri ile deneyelim.

a=1:10

a

a_1=c(NA,a[1:9])

a_1

Burada yapılan, a_1 yani 1 zaman gecikmeli bir a vektörü oluşturarak, bu vektörü “NA” ile başlattık ve a vektörünün eski serisindeki ilk 9 veriyi de ekledik.

Peki ya iki zaman dilimi kaydırmak istersek:

a_2=c(NA,NA,a[1:8])

a_2

Algoritmik olarak:

• 
  Gecikme zamanı kadar NA ekle
  
  • 
    rep(NA,gecikme_zamani)
    
• 
  Ardından veri serisini başlat, fakat sondan gecikme zamanı kadarına kadar olan kısmı ekle
  
  • 
    veri_serisi[1: (length(veri_serisi)-gecikme_zamani)]
    

şeklinde tamamlayabiliriz.

Dolayısıyla geciktirme fonksiyonumuz:

c(rep(NA,gecikme),veri_serisi[1: (length(veri_serisi)-gecikme) ]

şeklinde olacaktır.

Kodumuzu 4 gecikme için deneyelim

a_4=c(rep(NA,4),a[1:(length(a)-4)])

a_4

Bunu istersek bir fonksiyon olarak tanımlayabiliriz. Veri setimiz “x”, gecikme süremiz de “k” olacaktır.

gecik<-function(x, k){ c(rep(NA,k),x[1: ( length(x)-k)]) }

Artık ne zaman gecik komutunu çağırırsak bize gecikmeli bir zaman serisi üretir

6.3. Farklar

Fark almak konusunda ise “diff” komutu kullanılabilir. Özellikle durağan olmayan serilerde model yapmak için ilk serilerin alınması çok önemli olabilir. Bu konuda dikkat etmeniz gereken

• 
  lag=3 veya direkt olarak veri_serisinden sonra 3 girerseniz, 3 önceki veri ile farkını verir
  
• 
  differences=3 derseniz de 3. Farkı verir
  

?diff

a=a^3

diff(a,1)

a

diff(a,lag=1)

diff(a,differences=1)

diff(a,differences=2)

diff(a,lag=2)

diff(a,differences=3)

diff(a,3)

Mesela aşağıda, differences=3 parametresi ile 3 defa verinin yinelemeli 1. Farkları birbirine eşittir (6 6 6 6 6 6 6)

Peki fark neden önemlidir, talep serilerinde genelde veri seti durağan olmayabilir (enerji , elektrik) gibi, yani veri seti sürekli yükselen artan bir eğilim gösteriyordur. Fakat birinci farkı veya ikinci farkı, normal dağılıma daha yakın olabilir.

Mesela, veri setimizdeki ortalama tüketim verisine iki şekilde bakalım.

ot=veri$ortalama_tuketim

plot(ot,col="blue",type="l")

Yukarıdaki yükselen talep eğrisinin ilk farkını aldığımızda da aşağıdaki grafik ile karşılaşırız.

plot(diff(ot,1),col="blue",type="l")

İlk farkın histogramına bakarsak da , aslında ilk farkın yeterli olmadığını düşünebiliriz.

hist(diff(ot,1),col="blue",breaks = 300)

Yedinci farka bakarsak ise daha normal bir dağılımla karşılaşırız.

hist(diff(ot,7),col="blue",breaks = 300)

Görüldüğü üzere, ihtiyacımız olan araçlarımızı da gördük. Artık doğrusal regresyon ile konumuza başlayabiliriz.

6.4. Doğrusal Regresyon – Yeniden

Bu noktaya geldiğimizde:

• 
  elimizde ot değişkenine yüklenmiş veri$ortalama_tuketim
  
• 
  library(forecast) ile yüklenmiş “forecast” kütüphanesi
  
• 
  Ayrıca “gecik” isimli yukarıda tanımladığımız bir fonksiyon var.
  

Daha önceki bölümlerde lm(y~x) şeklinde fonksiyonumuzu görmüştük.

NOT: bundan sonraki kısımlarda, ekonometrik testler yapılmadan model yapılması yanlış sonuç verebilir. Bir bilene danışın, ekonometri anlatımı bu öğreticinin konusu değildir.

Şimdi ise lm(y~x1+x2+x3) şeklinde fonksiyonumuzu kullanacağız. Fonksiyonu bu şekilde yazarak aslında y bağımlı değişkenini x1,x2,x3’ler cinsinden ifade edecek bir modeli arıyoruz ki

Y= c+ c1*x1 + c2*x2 + c3 * x3

Olacak şekilde c (yani intercept), c1,c2,c3’leri program hesaplasın. Bu katsayılar öyle bir hesaplansın ki, x1,x2,x3’ün Y’ye uzaklıklarının karelerinin toplamı minimum olsun.

Biz de elektrik talebinin bir doğrusal regresyonuna bakalım.

lm(ot~gecik(ot,1))

Görüldüğü üzere modelimizin katsayıları hemen hesaplandı.

Peki modelimiz ne kadar iyi?

Bunun için önce orijinal verilerimizi gri ile çizdirdikten sonra üzerine modelin hesapladığı “y” değerlerini çizdirelim

plot(ot,col="grey",type="l")

m1=lm(ot~gecik(ot,1))

lines(fitted(m1),col="red",type="l")

Görüldüğü üzere modelimiz gayet güzel oturmuş gözüküyor.

Şimdi doğrusal regresyonumuza bir parametre daha ekleyelim. O da bir hafta öncesinin değeri olsun.

plot(ot,col="blue",type="l")

m2=lm(ot~gecik(ot,1)+gecik(ot,7))

lines(fitted(m2),col="red",type="l")

Peki hangi model daha iyi bir model oldu?

6.5. Model Sonuçları

Model sonuçlarının hangisi daha iyi oldu konusunda, summary komutunu kullanabilirsiniz.

summary(m1)

Görüldüğü üzere yanına *** konan değişkenler sonuca en çok etki edenler. Genelde burada “t value” en büyük, “Pr(>|t|)” ise 0.05 den küçük olan değişkenler alınır.

Aynı şekilde “Adjusted R-Squared” değeri ise ne kadar 1’e yakın ise o kadar iyidir, denebilir. Bu yüzden bu değeri yüksek istiyoruz.

İlgilenenler için  Akaike's ‘An Information Criterion AIC komutu da işe yarar. Bu değişkeni de

AIC(m1)
Şeklinde çağırabilirsiniz.

Model sonuçlarını değerlendirirken bir diğer önemli konu ise tortu/kalıntı yani residual’lara bakmaktır. Kalıntıların mümkün olduğunca random/rastgele dağılması çok önemlidir.

Biz de şimdi bu iki modelin kalıntılarına bakalım. Önce mavi olarak birinci modelin kalıntılarını çizelim sonra da ikinci modelinkileri. Bunu da “resid” fonksiyonu ile yapacağız. (m1 ve m2 yukarıda tanımlanmış lm modelleri)

plot(resid(m1),col="blue",type="l")

lines(resid(m2),col="red",type="l")

Sanki ikinci modelin sonucunda daha az kalıntı oluşmuş gibi gözüküyor.

6.6. Modeli İleri Doğru Hareket Ettirirsek

Bu kısım bilmeyenler için biraz “advanced” gelebilir. Ama mantığımız çok basit. Hem m1, m2 modelimiz var. “m1” modeli bir gün öncesini alarak geliyor, “m2” modeli ise hem gün öncesini hem de bir hafta öncesini dikkate alıyor.

Yapacağımız şey sırası ile:

• 
  Önce m1f ve m2f (m1 forecast ve m2 forecast kısaltması) iki veri serisi tanımlayıp bunları 0 ile dolduralım.
  
  • 
    m1f=rep(0,100)
    
  • 
    m2f=rep(0,100)
    
• 
  ot’deki (ortalama_tuketim) son 8 günü m1f ve m2f’in ilk sekiz değişkenine atayalım.
  
  • 
    lng=length(ot)
    
  • 
    tmp=ot[(lng-7):lng]
    
  • 
    m1f[1:8]=tmp[1:8]
    
  • 
    m2f[1:8]=tmp[1:8]
    
• 
  sonra 100 defa(100 adım ötesine kadar görmek için), her adımda geçtiğimiz 8 veriyi kullanarak bir adım ilerisini hesaplayalım (predict ile)
  
  • 
    for(syc in 1:100){
    
  • 
    tmp2=predict(m1,newdata=data.frame(ot=m1f[syc:(syc+7)]))
    
• 
  Her adımda elde ettiğimiz serinin son değerini m1f ve m2f’in sonuna ekleyelim.
  
  • 
    m1f[syc+7]=tmp2[8]
    
  • 
    
    
  • 
    tmp2=predict(m2,newdata=data.frame(ot=m2f[syc:(syc+7)]))
    
  • 
    m2f[syc+7]=tmp2[8]
    
  • 
    
    
  • 
    }
    

Ardından da m1f ve m2f’i çizdirelim.

Görüldüğü üzere sadece gün öncesini içeren modeli ileri taşıdığımızda model hızla sönümleniyor (mavi), bir hafta öncesini de dahil ettiğimiz model ise daha geç sönümleniyor(kırmızı).

Sonuç

Bu bölümde ekonometrik testlere girmeden model oluşumu ve modelin sonuçlarını tahmin ederek ileri doğru tahminler yapmaya çalıştık. Model analizimiz çok da doğru sayılamayacağından, modeli alıp, son haftadan ileri götürdükçe modelimizin saptığını gördük.

Bu bölümde doğrusal regresyonu gördük, bundan sonraki bölümlerde ise diğer yöntemleri de inceleyeceğiz.

Not: Aşağıdaki kodların en altında “Bonus Material olarak” 365 gün öncesini de baz alan bir model için kodlar var, bir önceki bölümden farklı değil, sadece 100 adım yerine 700 adım ötesine gidiyor.

Komutlar

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

## Program

ot=veri$ortalama_tuketim
install.packages("forecast")

library("forecast")

a=1:10

a

a_1=c(NA,a[1:9])

a_1

a_1=c(NA,NA,a[1:8])

a_2=c(NA,NA,a[1:8])

a_2

rep(NA,3)

a_4=c(rep(NA,4),a[1:(length(a)-4)])

a_4
#------------------------------------------------------------------
#gecik fonksiyonu

gecik<-function(x, k){ c(rep(NA,k),x[1: ( length(x)-k)]) }
#------------------------------------------------------------------
gecik(a,2)

gecik(a,5)

?diff

a=a^3

diff(a,1)

a

diff(a,lag=1)

diff(a,differences=1)

diff(a,differences=2)

diff(a,lag=2)

diff(a,differences=3)

diff(a,3)

diff(a,3)

diff(a,differences=3)

diff(a,1)

diff(diff(a,1),1)

diff(diff(diff(a,1),1),1)

ot=veri$ortalama_tuketim

plot(ot,col="blue",type="l")

plot(diff(ot,1),col="blue",type="l")

hist(diff(ot,1),col="blue",breaks = 300)

hist(diff(ot,7),col="blue",breaks = 300)
plot(ot,col="blue",type="l")

m1=lm(ot~gecik(ot,1))

lines(fitted(m1),col="red",type="l")
summary(m1)

AIC(m1)

AIC(m2)
plot(resid(m1),col="blue",type="l")

lines(resid(m2),col="red",type="l")
#------------------------------------------------------------------

# programımız

m1f=rep(0,100)

m2f=rep(0,100)

lng=length(ot)

tmp=ot[(lng-7):lng]

m1f[1:8]=tmp[1:8]

m2f[1:8]=tmp[1:8]
for(syc in 1:100){

tmp2=predict(m1,newdata=data.frame(ot=m1f[syc:(syc+7)]))

m1f[syc+7]=tmp2[8]

tmp2=predict(m2,newdata=data.frame(ot=m2f[syc:(syc+7)]))

m2f[syc+7]=tmp2[8]

}

plot(m1f,type="l",col="blue")

lines(m2f,type="l",col="red")

#------------------------------------------------------------------


# --------------------------------------------------------------------------

# Bonus material

# lm modeli ile üretilen fonksiyonun 2 yıl ileri uzatılmlış hali

m1=lm(ot~gecik(ot,1)+gecik(ot,7)+gecik(ot,14))

m2=lm(ot~gecik(ot,1)+gecik(ot,7)+gecik(ot,14)+gecik(ot,21)+gecik(ot,28)+gecik(ot,364)+gecik(ot,365))
m1f=rep(0,1000)

m2f=rep(0,1000)

lng=length(ot)

tmp=ot[(lng-365):lng]

m1f[1:366]=tmp[1:366]

m2f[1:366]=tmp[1:366]
for(syc in 1:725){

tmp2=predict(m1,newdata=data.frame(ot=m1f[syc:(syc+365)]))

m1f[syc+366]=tmp2[366]

tmp2=predict(m2,newdata=data.frame(ot=m2f[syc:(syc+365)]))

m2f[syc+366]=tmp2[366]

}

plot(m1f,type="l",col="blue",ylim=c(15000,40000))

lines(m2f,type="l",col="red")