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_trend <- 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_trend[["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_trend <- veri_ts - veri_trend[["fitted.values"]]

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

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

## 
## Autocorrelations of series 'periyot_trend', 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_1 <- 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 çıkartmamız gerekiyor.

mevsim3 <- veri_ts - veri_trend_1
head(mevsim3,n=36)

##               Jan          Feb          Mar          Apr          May
## 1973           NA           NA           NA           NA           NA
## 1974   1244.43125  -5987.08842  -6644.94317  -9240.61225  -8877.00792
## 1975   5386.46054  -1263.04279  -3774.80125  -8881.94913 -11311.41638
##               Jun          Jul          Aug          Sep          Oct
## 1973           NA   9841.39842  14898.72092  14892.45550   3604.40875
## 1974    -77.18058   9125.58383  15241.64129  10288.53987  -2292.94225
## 1975  -2311.87663   8148.75617  13941.57687  10060.55892  -6412.73892
##               Nov          Dec
## 1973  -4435.71604  -5518.18642
## 1974  -8002.64087  -2999.57629
## 1975  -9583.86033  -1785.75867

Mevsim serisinin ortalamaları

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

donemort3<-t(matrix(data=mevsim3, nrow = 12))

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

head(donemort3)

##           [,1]      [,2]      [,3]       [,4]       [,5]        [,6]      [,7]
## [1,]        NA        NA        NA         NA         NA          NA  9841.398
## [2,]  1244.431 -5987.088 -6644.943  -9240.612  -8877.008   -77.18058  9125.584
## [3,]  5386.461 -1263.043 -3774.801  -8881.949 -11311.416 -2311.87663  8148.756
## [4,] 10481.996  2802.820 -3880.936  -9695.276 -11061.992 -4935.45200  7910.461
## [5,] 12666.531  4699.200 -5136.574 -14090.190 -13786.849 -3093.95550 12455.279
## [6,] 11804.552  4676.338 -2511.719 -16272.308 -15786.122 -2368.93296 10589.750
##          [,8]      [,9]     [,10]      [,11]     [,12]
## [1,] 14898.72 14892.455  3604.409  -4435.716 -5518.186
## [2,] 15241.64 10288.540 -2292.942  -8002.641 -2999.576
## [3,] 13941.58 10060.559 -6412.739  -9583.860 -1785.759
## [4,] 12553.10  7480.415 -6545.047  -7306.031  4008.703
## [5,] 15675.59  8854.138 -4461.604 -12721.138 -2238.178
## [6,] 14851.83 12373.875 -4417.750 -12814.542 -3622.708

Periyottaki dönemlerin ortalama değerleri:

colMeans(donemort3, na.rm = T)

##  [1]  12293.710  -8992.596 -12270.776 -27486.215 -19556.515   5455.674
##  [7]  33711.582  37332.788  13396.718 -10784.080 -22267.319  -2109.830

Periyottaki dönemlerin ortalama değerlerinin toplamı:

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

## [1] -1276.859

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ğerlerden çıkartılır ve mevsimsel endeks serisi bulunur.

Periyottaki dönemlerin ortalama değerlerinin ortalaması:

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

## [1] -106.4049

Mevsimsel Endeksin Bulunması

endeks1 <- colMeans(donemort3, na.rm = T) - mean(colMeans(donemort3, na.rm = T))

endeks1

##  [1]  12400.115  -8886.191 -12164.371 -27379.810 -19450.110   5562.079
##  [7]  33817.987  37439.193  13503.123 -10677.675 -22160.914  -2003.426

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.

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

indeks<-  matrix(data = endeks1, nrow = 587)

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

head(indeks)

##            [,1]
## [1,]  12400.115
## [2,]  -8886.191
## [3,] -12164.371
## [4,] -27379.810
## [5,] -19450.110
## [6,]   5562.079

Hata Terimli Trend Bileşeninin Bulunması

trendhata <- veri_ts - indeks

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

plot.ts(trendhata)

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

trend1<-tslm(ts(trendhata)~trend)

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

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

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

## Time Series:
## Start = 1 
## End = 24 
## Frequency = 1 
##  [1] 149853.0 150221.7 150590.3 150959.0 151327.6 151696.3 152064.9 152433.6
##  [9] 152802.2 153170.9 153539.5 153908.1 154276.8 154645.4 155014.1 155382.7
## [17] 155751.4 156120.0 156488.7 156857.3 157226.0 157594.6 157963.2 158331.9

Tahmin ve Hata Serilerinin Bulunması

Tahmin Bileşenlerini Bulalım:

tahmin = mevsimsel endeks + saf trend serisi

tahmin<- indeks+trend1[["fitted.values"]]
tahmin <- ts(data = tahmin, start = c(1973, 01), end = c(2021, 11), frequency = 12 )
plot(tahmin,main = "Tahmin Serisi Grafiği")

Hata Bileşenlerini Bulalım:

hata = saf trend serisi - tahmin

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

Modelin Güvenirliği

Toplamsal 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(tahmin) ,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(hata, lag.max = 42,  ylim=c(-1,1), lwd=5,main="Hatalar için ACF Grafiği")

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

## 
##  Box-Ljung test
## 
## data:  hata
## X-squared = 9137.7, 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 Toplumsal Ayrıştırma yöntemi uygun değildir.