Baslangic
R Studio ve ilk kodlar
Elektrik Talep Analizi - (Histogram)
Elektrik Talep Analizi - (Filtreleme)
Elektrik Talep Analizi - (Verinin Bilesenleri)
Elektrik Talep Analizi - (Programlamaya Giris)
Dogrusal Regresyon ile Elektrik Talebi
Ussel Duzgunlestirme ve STL ile Elektrik Talebi
Elektrik Talebinde Farkli Periyodlar (ARIMA ve TBATS)
2015-2024 - Aylik Dogalgaz Talebi
NOAA'dan gunluk sicaklik verilerine erisim
TCMB(Merkez Bankasi) verilerine erisim
EIA (ABD Eneji Bilgi Dairesi) verilerine erisim
Dunya Bankasi verileri, otomatik modelleme, ggplot2 ve tema
Merkez Bankasi ile EIA verileri birarada: Akaryakit fiyatlari
Enflasyon ve Akaryakit Iliskisini inceleme
R Dersleri - Elektrik Talebi - Dorusal Regresyon

R ile Enerji Analizi B繹l羹m 6- Dorusal 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 inceleyeceiz, eer teorik detaylar覺na girmek isterseniz, yaz覺n覺n i癟indeki linkleri takip edebilirsiniz.

Dorusal regresyon modeli kurmadan 繹nce bir 癟ok test yap覺lmas覺 gerekir. Bu testler Rda ve ekonometri kitaplar覺nda mevcut. Bu k覺s覺mda daha 癟ok lm fonksiyonu ve 繹zelliklerine deineceiz.

 

Bu B繹l羹mdeki Fonksiyonlar

read.csv

rep(deiken, 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 eriebiliriz.

6.1. Forecast paketi

Zaman serileri ve modelleri, T羹rkiyede 羹niversitelerde g繹zlemlediim kadar覺 ile roket bilimine yak覺n bir zorluktaym覺 gibi 繹retiliyor. ok fazla gramerden pratie 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 olduu zaman serilerini kullanma imkan覺 vermesi

Say覺labilir. Ayr覺ca Hyndman覺n videolar覺 da internetten eriilebilir.

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

nce forecast paketini R Studiodan y羹kleyelim.

install.packages("forecast")

 

Eer eksik bir paket belirtirse ki normalde t羹m paketleri de y羹kler: install.packages(zoo) yazarak forecastin 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 dorusal regresyon 癟al覺anlar i癟in en 癟ok kullan覺lan iki konudan biri zaman gecikmesi, dieri de farklar konusudur.

Rda zaman gecikmesi lagler konusunda internet gruplar覺nda da farkl覺 tart覺malar var. Ben o y羹zden kolayca lag oluturan ama sonra deiik 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 gereince NA, Rda 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覺 oluturmak a癟覺s覺nda ben size daha direkt bir y繹ntemi g繹stermek isterim. Bir veri serisine gecikme eklediimizde pratikte ne yap覺yoruz? Veriyi gecikme suresi kadar aa覺 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羹 oluturarak, bu vekt繹r羹 NA ile balatt覺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

o   rep(NA,gecikme_zamani)

-          Ard覺ndan veri serisini balat, fakat sondan gecikme zaman覺 kadar覺na kadar olan k覺sm覺 ekle

o   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 duraan 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 aa覺da, differences=3 parametresi ile 3 defa verinin yinelemeli 1. Farklar覺 birbirine eittir (6 6 6 6 6 6 6)

Peki fark neden 繹nemlidir, talep serilerinde genelde veri seti duraan olmayabilir (enerji , elektrik) gibi, yani veri seti s羹rekli y羹kselen artan bir eilim 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 erisinin ilk fark覺n覺 ald覺覺m覺zda da aa覺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 dorusal regresyon ile konumuza balayabiliriz.

6.4. Dorusal Regresyon Yeniden

Bu noktaya geldiimizde:

-          elimizde ot deikenine 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 deildir.

imdi ise lm(y~x1+x2+x3) eklinde fonksiyonumuzu kullanaca覺z. Fonksiyonu bu ekilde yazarak asl覺nda y ba覺ml覺 deikenini x1,x2,x3ler cinsinden ifade edecek bir modeli ar覺yoruz ki

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

Olacak ekilde c (yani intercept), c1,c2,c3leri program hesaplas覺n. Bu katsay覺lar 繹yle bir hesaplans覺n ki, x1,x2,x3羹n Yye uzakl覺klar覺n覺n karelerinin toplam覺 minimum olsun.

File:Linear least squares example2.png

Biz de elektrik talebinin bir dorusal 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 deerlerini 癟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 dorusal regresyonumuza bir parametre daha ekleyelim. O da bir hafta 繹ncesinin deeri 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 deikenler sonuca en 癟ok etki edenler. Genelde burada t value en b羹y羹k, Pr(>|t|) ise 0.05 den k羹癟羹k olan deikenler al覺n覺r.

Ayn覺 ekilde Adjusted R-Squared deeri ise ne kadar 1e yak覺n ise o kadar iyidir, denebilir. Bu y羹zden bu deeri y羹ksek istiyoruz.

襤lgilenenler i癟in  Akaike's An Information Criterion AIC komutu da ie yarar. Bu deikeni de

AIC(m1)

 

eklinde 癟a覺rabilirsiniz.

Model sonu癟lar覺n覺 deerlendirirken bir dier 繹nemli konu ise tortu/kal覺nt覺 yani residuallara bakmakt覺r. Kal覺nt覺lar覺n m羹mk羹n olduunca 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覺 olumu gibi g繹z羹k羹yor.

6.6. Modeli 襤leri Doru 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.

o   m1f=rep(0,100)

o   m2f=rep(0,100)

-          otdeki (ortalama_tuketim) son 8 g羹n羹 m1f ve m2fin ilk sekiz deikenine atayal覺m.

o   lng=length(ot)

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

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

o   m2f[1:8]=tmp[1:8]

-          sonra 100 defa(100 ad覺m 繹tesine kadar g繹rmek i癟in), her ad覺mda ge癟tiimiz 8 veriyi kullanarak bir ad覺m ilerisini hesaplayal覺m (predict ile)

o   for(syc in 1:100){

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

-          Her ad覺mda elde ettiimiz serinin son deerini m1f ve m2fin sonuna ekleyelim.

o   m1f[syc+7]=tmp2[8]

o  

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

o   m2f[syc+7]=tmp2[8]

o  

o   }

Ard覺ndan da m1f ve m2fi 癟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 ettiimiz model ise daha ge癟 s繹n羹mleniyor(k覺rm覺z覺).

Sonu癟

Bu b繹l羹mde ekonometrik testlere girmeden model oluumu ve modelin sonu癟lar覺n覺 tahmin ederek ileri doru tahminler yapmaya 癟al覺t覺k. Model analizimiz 癟ok da doru say覺lamayaca覺ndan, modeli al覺p, son haftadan ileri g繹t羹rd羹k癟e modelimizin sapt覺覺n覺 g繹rd羹k.

Bu b繹l羹mde dorusal regresyonu g繹rd羹k, bundan sonraki b繹l羹mlerde ise dier y繹ntemleri de inceleyeceiz.

Not: Aa覺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覺 deil, 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")