Verinin yüklenmesi:

veri <-read.csv("https://raw.githubusercontent.com/gungorrbaris/zaman-serisi-analizi-R/main/data/Table_7.6_Electricity_End_Use.csv")
knitr::kable(head(veri,n=10), align = "c")

Verinin Düzenlenmesi:

Tarih değişkenini Ay-Yıl formatına getirelim:

tarih = seq(from = as.Date("1973-01-01"), to = as.Date("2021-11-01"), by = 'month')
veri$Month <- tarih
veri$Month <-format(veri$Month, "%m-%Y")
knitr::kable(head(veri), align = "c")

Değişken isimlerini türkçeleştirelim:

colnames(veri) <- c("Tarih","Elektrik_Tuketimi")
knitr::kable(head(veri), align = "c")

Zaman Serisinin Oluşturulması

veri_ts <- ts(data = veri$Elektrik_Tuketimi, start = c(1973, 01), end = c(2021, 11), frequency = 12 )

Verimizde gözlemler aylık olduğu için frekans değeri (frequency) 12 olarak belirlenmelidir.

Zaman Serisi Grafiğinin Çizilmesi

ts.plot(veri_ts,gpars=list(xlab="Zaman", ylab="Elektrik Tüketimi"))

Zaman Serisi Grafiğine bakıldığında:

• Pozitif yönlü artış vardır. Bu yüzden Serinin artan bir trende sahip olduğunu söyleyebiliriz.

• Seride düzenli dalgalanmaların olduğu görülmektedir. Bu seride mevsimsellik olabileceğini gösterir. Daha kesin sonuç için ACF grafiğini çizelim.

ACF Grafiğinin Çizilmesi

! BİLGİ !

Gecikmeli seri: Serilerdeki verilerin dönem kaydırılması işlemi ile ortaya çıkan serilerdir. Zt bir zaman serisi değişkeni ise; bir dönem gecikmeli zaman serisi Zt-1, iki dönem gecikmeli Zt-2 ve k dönem gecikmeli Zt-k ile gösterilir.

Otokorelasyon katsayısı: Zaman serileri ile bu serilerin gecikmeli serileri arasındaki ilişkileri verir.

ACF: Otokorelasyon fonksiyonu (Autocorrelation function) demektir. Otokorelasyon katsayısı değerlerinden oluşur.

Seride trendin olduğunu anlamak için ACF grafiğindeki ilk dört gecikmeye bakmak yeterlidir. İlk dört gecikme sınırlar dışındaysa seri için “trende sahiptir” diyebiliriz.

library(fpp)

## Zorunlu paket yükleniyor: forecast

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Zorunlu paket yükleniyor: fma

## Zorunlu paket yükleniyor: expsmooth

## Zorunlu paket yükleniyor: lmtest

## Zorunlu paket yükleniyor: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Zorunlu paket yükleniyor: tseries

Acf(veri_ts,lag.max = 42,  ylim=c(-1,1), lwd=5,main="Orijinal veri ACF Grafiği")

ACF Grafiğine bakıldığında:

• Tüm gecikmeler (4 gecikme olması yeterlidir) sınırlar dışında olduğu için zaman serisinin trende sahip olduğu söylenebilir. Bu yüzden serinin farkı alınıp durağan hale getirilmelidir.

• Trend ile birlikte düzenli dalgalanmaların olduğu gözükse de baskın bir mevsimsellikten şu an için bahsedemeyiz.

Birinci Farkın Alınması

! BİLGİ !

Fark Alınması: Zaman serisinin akışkanlı bir şekilde son değerlerinden belli bir dönem önceki değerlerinin çıkarılması işlemidir. Bu işlem sayesinde serideki trend ya da mevsimsel dalgalanmaları yok etmek mümkün olmaktadır.

veri_birinci_fark <- diff(veri_ts, differences = 1)

Farkı alınmış seri için ACF:

Acf(veri_birinci_fark,lag.max = 42,  ylim=c(-1,1), lwd=5,main="Birinci Fark sonrası ACF Grafiği")

MEVSİMSELLİK

Birinci Farkı alınmış seri için ACF Grafiğine bakıldığında:

• Yine tüm gecikmeler sınırlar dışında olduğu için birinci farkı alınmış zaman serisinin de trende sahip olduğu söylenebilir. Fakat bir öncekine göre bu ACF grafiğinde baskın bir mevsimsellik görülüyor. Dalgalanmalar ve gecikmelerdeki düzenli sıçramalar bunu destekliyor.

Trend Bileşeni ve Periyodun Bulunması

Trend Bileşeni:

Serinin Trend Bileşeni Regresyon analizi yardımı ile oluşturulur.

veri_trend2 <- tslm(veri_ts~trend)

• Burada dikkat edilmesi gereken, denkleme orijinal serinin eklenmiş olmasıdır.

• Oluşturulan Regresyon denklemindeki fitted.values değişkeni dikkate alınır.

• fitted.values değişkeni orijinal serinin trend bileşeni olur.

head(veri_trend2[["fitted.values"]],n=36)

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1973 149788.4 150157.2 150526.1 150895.0 151263.9 151632.8 152001.6 152370.5
## 1974 154214.9 154583.8 154952.7 155321.5 155690.4 156059.3 156428.2 156797.0
## 1975 158641.4 159010.3 159379.2 159748.1 160117.0 160485.8 160854.7 161223.6
##           Sep      Oct      Nov      Dec
## 1973 152739.4 153108.3 153477.1 153846.0
## 1974 157165.9 157534.8 157903.7 158272.6
## 1975 161592.5 161961.3 162330.2 162699.1

Periyodun Bulunması:

Orijinal zaman Serisi ve trend bileşeninin farkını alarak periyodu bulabiliriz.

periyot_trend2 <- veri_ts - veri_trend2[["fitted.values"]]

Serinin periyoduna sahip mevsimsel bileşen serisinin ACF değerleri ile periyodu bulalım:

Acf(periyot_trend2,lag.max = 42,  ylim=c(-1,1), lwd=5,plot=FALSE)

## 
## Autocorrelations of series 'periyot_trend2', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.720  0.265 -0.021  0.043  0.320  0.455  0.318  0.040 -0.025  0.246 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.673  0.916  0.653  0.218 -0.057  0.012  0.287  0.416  0.281  0.009 -0.057 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.204  0.622  0.853  0.602  0.185 -0.078 -0.006  0.260  0.387  0.257 -0.015 
##     33     34     35     36     37     38     39     40     41     42 
## -0.079  0.177  0.587  0.814  0.567  0.158 -0.097 -0.029  0.230  0.354

ACF değerlerine bakıldığında:

• En büyük değerin 12. gecikme (0.916), en büyük ikinci gecikmenin ise 24. gecikme (0.853) olduğu gözükmektedir.

• 12 ve 24 gecikme arasında toplam 12 değer olduğu için; birinci farkı alınmış serinin periyodunun 12 olduğunu söyleyebiliriz.

Buradaki başlıkta periyot bulma işlemini yapmış ve aynı sonuca ulaşmıştık.

Mevsimsel, döngüsel veya düzensiz dalgalanmaları yok etme ya da belli bir miktar düzleştirme amacıyla basit hareketli ortalama ya da merkezsel hareketli ortalama işlemlerine ihtiyaç duyulmaktadır.

Ayrıştırma Yönteminde merkezsel hareketli ortalamayı mevsimsel bileşeni bulmak için kullanacağız

Merkezsel Hareketli Ortalama Hesabı

veri_trend_2 <- ma(veri_ts, order = 12, centre = TRUE)

Burada order = germe sayısı, bir önceki adımda bulduğumuz periyot ile aynı olmalıdır.

Mevsimsel Bileşenin Bulunması

Mevsimsel bileşene ulaşmak için yukarıda hesapladığımız merkezsel hareketli ortalama serisini, orijinal seriden bölmemiz gerekiyor.

mevsim4 <- veri_ts / veri_trend_2
head(mevsim4,n=36)

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1973        NA        NA        NA        NA        NA        NA 1.0689556
## 1974 1.0087181 0.9580513 0.9534005 0.9350386 0.9374657 0.9994565 1.0640906
## 1975 1.0374040 0.9912432 0.9738812 0.9386246 0.9219252 0.9840898 1.0557923
##            Aug       Sep       Oct       Nov       Dec
## 1973 1.1044922 1.1045693 1.0253052 0.9688891 0.9613348
## 1974 1.1066813 1.0717819 0.9840378 0.9443436 0.9791481
## 1975 1.0948663 1.0681237 0.9567380 0.9356111 0.9880502

Mevsim serisinin ortalamaları

Mevsimsel bileşendeki hata teriminin yok edilebilmesi için her bir periyottaki dönemlerin ortalama değerleri hesaplanır.

donemort4<-t(matrix(data=mevsim4, nrow = 12))

## Warning in matrix(data = mevsim4, nrow = 12): veri uzunluğu [587] dize sayısının
## [12] böleni veya tam böleni değil

head(donemort4)

##          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     [,7]
## [1,]       NA        NA        NA        NA        NA        NA 1.068956
## [2,] 1.008718 0.9580513 0.9534005 0.9350386 0.9374657 0.9994565 1.064091
## [3,] 1.037404 0.9912432 0.9738812 0.9386246 0.9219252 0.9840898 1.055792
## [4,] 1.069853 1.0185930 0.9743572 0.9362177 0.9276360 0.9679462 1.051004
## [5,] 1.079322 1.0292605 0.9681632 0.9130324 0.9150998 0.9809471 1.076628
## [6,] 1.071610 1.0282985 0.9848617 0.9023607 0.9056363 0.9858898 1.062838
##          [,8]     [,9]     [,10]     [,11]     [,12]
## [1,] 1.104492 1.104569 1.0253052 0.9688891 0.9613348
## [2,] 1.106681 1.071782 0.9840378 0.9443436 0.9791481
## [3,] 1.094866 1.068124 0.9567380 0.9356111 0.9880502
## [4,] 1.080428 1.047667 0.9584520 0.9537730 1.0252505
## [5,] 1.096220 1.054184 0.9727635 0.9224478 0.9863869
## [6,] 1.087688 1.072675 0.9741863 0.9254711 0.9789783

Periyottaki dönemlerin ortalama değerleri:

colMeans(donemort4, na.rm = T)

##  [1] 1.0497120 0.9724287 0.9561884 0.8991415 0.9234883 1.0151612 1.1214735
##  [8] 1.1383244 1.0533474 0.9596636 0.9170578 0.9916020

Periyottaki dönemlerin ortalama değerlerinin toplamı:

sum(colMeans(donemort4, na.rm = T))

## [1] 11.99759

Bu değerin 0 olması gerekmektedir. 0 olmadığı için ortalamaların ortalaması alınmalıdır. Bu ortalama değer diğer tüm ortalama değerlere bölünür ve mevsimsel endeks serisi bulunur.

Periyottaki dönemlerin ortalama değerlerinin ortalaması:

mean(colMeans(donemort4, na.rm = T))

## [1] 0.9997991

Mevsimsel Endeksin Bulunması

endeks2 <- colMeans(donemort4, na.rm = T) / mean(colMeans(donemort4, na.rm = T))

endeks2

##  [1] 1.0499230 0.9726242 0.9563805 0.8993222 0.9236739 1.0153652 1.1216988
##  [8] 1.1385531 1.0535591 0.9598564 0.9172421 0.9918013

Bu 12 mevsimsel endeks değeri periyottaki dönemlere dikkat edilerek (ilk döneme karşılık M1; ikinci döneme karşılık M2; üçüncü döneme karşılık M3, … , on ikinci döneme karşılık M12 gelecek şekilde) seri boyunca tekrar tekrar yazılabilmesi için aşağıdaki adım gerçekleştirilir.

Mevsimsel Endeksin ortalaması: >> Mevsimsel endekin ortalaması yani Dönem ortalamalarının ortalaması 1’e eşit olmalıdır.

mean(endeks2)

## [1] 1

• Dönem ortalamalarının ortalaması 1’e eşittir. Bu hata yapmadığımızı gösterir.

İndeks Değerlerini Seri Boyunca Yazdırma İşlemi:

indeks2<-  matrix(data = endeks2, nrow = 587)

## Warning in matrix(data = endeks2, nrow = 587): veri uzunluğu [12] dize sayısının
## [587] böleni veya tam böleni değil

head(indeks2)

##           [,1]
## [1,] 1.0499230
## [2,] 0.9726242
## [3,] 0.9563805
## [4,] 0.8993222
## [5,] 0.9236739
## [6,] 1.0153652

Hata Terimli Trend Bileşeninin Bulunması

trendhata2 <- veri_ts / indeks2

Hata Terimli Trend Bileşeninine ait zaman serisi grafiği:

plot.ts(trendhata2)

>>Elde edilen bu trent serisinde hata terimi de olduğundan seriye ‘’trendhata2’’ adı verilir. Bu seriyi hata teriminden arındırabilmek amacıyla trendhata serisine doğrusal regresyon uygulanır.

trend2<-tslm(ts(trendhata2)~trend)

• Burada fitted.values değişkeni dikkate alınır.

• fitted.values değişkeni orijinal serinin trend bileşeni olur.

head(trend2[["fitted.values"]],n=24)

## Time Series:
## Start = 1 
## End = 24 
## Frequency = 1 
##  [1] 150201.9 150569.1 150936.3 151303.5 151670.7 152037.9 152405.1 152772.3
##  [9] 153139.5 153506.7 153873.9 154241.1 154608.2 154975.4 155342.6 155709.8
## [17] 156077.0 156444.2 156811.4 157178.6 157545.8 157913.0 158280.2 158647.3

Tahmin ve Hata Serilerinin Bulunması

Tahmin Bileşenlerini Bulalım:

tahmin = mevsimsel endeks + saf trend serisi

tahmin2<- indeks2*trend2[["fitted.values"]]
tahmin2 <- ts(data = tahmin2, start = c(1973, 01), end = c(2021, 11), frequency = 12 )
plot(tahmin2,main = "Tahmin Serisi Grafiği")

Hata Bileşenlerini Bulalım:

hata = saf trend serisi - tahmin

hata2 <- ts(veri_ts) - ts(tahmin2)
plot(hata2,main = "Tahmin serisi için Hata Serisi Grafiği")

Modelin Güvenirliği

Çarpımsal modelin orijinal seri üzerinde geçerli bir model olup olmadığını kontrol edelim.

Orijinal seri ile tahmin serisinin uyumu:

plot( window(veri_ts), 
      xlab="Zaman", ylab="",lty=1, col=4, lwd=2,main="Orijinal Seri ve Tahmin Serisinin Birlikte Grafiği")
lines( window(tahmin2) ,lty=3,col=2,lwd=3)
legend("topleft",c(expression(paste(OrijinalSeri )),
                   expression(paste(Tahmin ))),
       lwd=c(2,2),lty=c(1,3), cex=0.6, col=c(4,2))

Zaman Serisi Grafiğine bakıldığında:

• Tahmin ile orijinal seri arasında sapmanın yüksek olduğunu dolayısıyla uyum göstermediğini söyleyebiliriz.

• Fakat emin olmak için hataların akgürültü olup olmadığına bakabiliriz.

Acf(hata2, lag.max = 42,  ylim=c(-1,1), lwd=5,main="Hatalar için ACF Grafiği")

Box.test(hata2,type = "Ljung",lag = 42)

## 
##  Box-Ljung test
## 
## data:  hata2
## X-squared = 12215, df = 42, p-value < 2.2e-16

SONUÇ

Bir önceki adımda görüldüğü üzere;

• Hatalar Arasında ilişki olduğunu yani hataların akgürültü olmadığını gördük.

• Sonuç olarak ABD’deki aylık elektrik tüketimi serisinin analizi için Çarpımsal Ayrıştırma yöntemi uygun değildir.