Selamlar, serinin üçüncü yazısına hoşgeldiniz. İlk iki yazıda, Ocak 1966 ve Ekim 1975 yılları arasında Boston’da aylık silahlı soygunların verisini incelemiş ve bu veri üzerinde 2 farklı metod kullanarak tahminlemeler gerçekleştirmiştik. Bunlar zaman serisi regresyonu ve üstel düzeltme metodlarıydı. Bu yazının konusu ise aynı veri üzerinde ARIMA modellerini kullanarak tahminleme yapmak olacak.
Bir önceki yazıda yaptığımız gibi, burada da verinin tanıtım aşamasını kısa tutup, sadece kodları vererek geçeceğiz. İndirme linki de serinin ilk yazısında bulunuyor.
Veriyi R üzerine alalım ve hatırlamak amacıyla grafiğini tekrar çizdirelim:
Data <- read.csv("MonthlyBostonArmedRobberies.csv",
stringsAsFactors = FALSE,
sep = ";")
nrow(Data)
ncol(Data)
head(Data, n = 10)
Üstteki kod bloğunu zaten konuşmuştuk.
Sonuçları hatırlayalım:
Evet verinin ilk hali bu şekildeydi, daha sonra aşağıdaki kod bloğuyla zaman serisi verisine çevirmiştik:
install.packages("fpp2")
library(fpp2)
Data_ts <- ts(data = Data$Value,
start = c(1966, 01),
end = c(1975, 10),
frequency = 12)
Verinin son hali:
Evet veriyi R’a zaman serisi olarak tanıttıktan sonra, grafiğini de aşağıdaki kodlar ile çizdirmiştik:
autoplot(Data_ts) +
ggtitle("Boston Aylık Silahlı Soygun Verisi") +
xlab("Ay") +
ylab("Sayı")
Kod bloklarını açıklamadan devam ediyorum, önceki yazılarda zaten detaylı bir biçimde üzerinde durduk.
Grafiği inceleyelim:
Evet veri setinin grafik üzerinde görünümü bu şekildeydi. Grafiği inceleyip, yorumlar yapmış ve detaylarına girmiştik. Bu yüzden yine burada bu kısımları geçecek ve yazının asıl konusuna odaklanacağız.
Yazının ana konusunun ARIMA modelleme üzerine olacağından bahsetmiştik, fakat buna geçmeden önce, üzerinde durulması gereken bazı kavramlar bulunuyor. Örneğin “durağanlık”. Bir zaman serisinin “durağan” olması demek, o serinin özelliklerinin zamandan etkilenmemesi, ondan bağımsız olması demektir. Yani gözlemler zamandan etkileniyorsa, seri durağanlığını kaybetmiş olur. Örneğin trend ya da mevsimsellik veya her ikisini birden barındıran zaman serilerini ele alalım. Eğer bir seri bunlardan herhangi birini ya da ikisini birden barındırıyorsa, durağan değildir çünkü trend ve mevsimsellik zaman serisinin gözlemlerini etkileyecek ve zamana bağlı olarak değiştirecektir. Örneğin mevsimsellik gösteren bir seride mevsimselliğin frekansı ne ise, gözlemler o frekansa geldiklerinde aynı davranışı gösterecekler ve zamandan etkilenmiş olacaklardır. Bu da serinin durağan olmadığının göstergesidir. Durağan serilere örnek verecek olursak da, ilk yazıda bahsettiğimiz beyaz gürültü(white noise) serisini örnek verebiliriz. Bu serinin hareketleri neredeyse tamamen rastgeledir ve herhangi tahmin edilebilir bir patern göstermez. Bu da beyaz gürültü serisinin durağan olduğunun kanıtıdır. Son olarak, bazı periyodik hareketler sergileyen seriler de durağan olabilir. Çünkü ilk yazıda da bahsettiğimiz gibi, periyodik hareketlerin frekansı sabit değildir ve kolay kolay tahmin edilemez.
Durağan serileri bir de görsel olarak görelim. Bunu yaparken de bizim serimiz ile karşılaştıralım ki iyice somutlaşsın:
Evet 3 farklı zaman serisini üst kısma yerleştirdik. Alt taraftaki bizim kendi verimiz, üsttekiler ise iki farklı zaman serisi. Bizim verimize baktığımızda, zaten ilk göze çarpan net bir trend var. Aynı zamanda çok aşırı olmasa da mevsimsellik de olduğunu biliyoruz. Bu yüzden bizim verimiz durağan olarak sınıflandırılamaz. Fakat üstteki iki veriyi incelediğimizde, ikisinin de herhangi bir trend veya mevsimsellik göstermediğini görebiliyoruz. İkisinin de hareketleri rastgele gibi görünüyor. Sağ üstteki verinin gösterdiği periyodik hareketler de aslında üstte bahsettiğimiz olayı bize gösteriyor. Frekansı belirsiz olan bu hareketler, tahmin edilebilir olmadığı için, verinin durağanlığını bozmaz.
Peki, durağan olmayan bir seri nasıl durağan hale getirilebilir? Burada da karşımıza “fark alma” veya “mevsimsel fark alma” metodları çıkıyor. Fark alma derken, aslında çok basit bir şekilde zaman serisinin ardışık gözlemlerinin birbirlerinden çıkarılmasından bahsediyoruz. Mevsimsel fark alma derken de, zaman serisinin bir önceki mevsimsel periyoddaki değerlerinin birbirlerinden çıkarılmasından bahsediyoruz. Örneğin, birinci çeyrek verisinden, bir önceki yılın birinci çeyrek verisinin çıkarılması. Hatırlarsanız varyansı değişken olan zaman serilerinde varyansın sabitleştirilmesi için seriye logaritma alma gibi dönüşümler uyguladığımızı ilk yazıda anlatmıştık. Aynı şekilde fark alma metodu kullanılarak da, durağan olmayan bir zaman serisi durağanlaştırılabilir. Yani dönüşümler varyansı stabilize ederken, fark alma da verinin trend ve mevsimselliğini stabilize ederek eleyebilir. Tabii bazı durumlarda sadece bir fark almak yetmeyecek, bir kez daha fark almak gerekecektir. Ya da normal fark almanın üstüne bir de mevsimsel fark almak da gerekebilir. Buna tabii ki veriye bakarak karar vereceğiz.
Bir zaman serisi üzerinde hangi türden kaç fark alma gerçekleştirilmesi gerektiğini, R üzerinde otomatik olarak belirleyen fonksiyonlar bulunuyor. Hem bunu hem de logaritma dönüşümünü bizim verimize uygularak verinin geldiği son hali görelim:
ndiffs(Data_ts)
nsdiffs(Data_ts)
cbind("Data" = Data_ts,
"Data_Log" = log(Data_ts),
"Data_Log_Diff" = diff(log(Data_ts), 1)) %>%
autoplot(facets = TRUE) +
xlab("Ay")+
ylab("Sayı") +
ggtitle("Boston Aylık Silahlı Soygun Verisi")
Kod bloğunu açıklayalım. “ndiffs()” komutu kaç normal fark alınması gerektiğini, “nsdiffs()” komutu ise kaç mevsimsel fark alınması gerektiğini bize veriye birim kök testleri uygulayarak bildiriyor. Bu testlerin detayına yazıyı uzatmamak adına burada değinmeyeceğim. Komutların sonucunu görelim:
Evet R bize zaman serisine 1 tane normal fark alma gerçekleştirmemizin onu durağan yapmak için yeterli olduğunu söylüyor. Mevsimsel farka da ihtiyaç olmadığını görebiliyoruz.
Kodun devamında ise, “cbind” komutu ile bir matris yaratıyor ve verinin kendisini, log dönüşümü uygulanmış halini ve en son da hem log dönüşümü uygulanmış hem de bir tane normal fark alma işlemi gerçekleştirilmiş halini 3 farklı kolon olarak bu matrisin içine yerleştiriyoruz. Sonrasında ise grafiğini çizdiriyoruz. Buradaki “diff()” komutu fark alırken, içerisine verdiğimiz “1” argümanı da birinci gecikmede fark alması gerektiğini yani ardışık gözlemlerin farkını almasını söylüyor. Mevsimsel fark alıyor olsaydık, örneğin çeyrek verisi için bu argümanının değerini “4” olarak atayacaktık. Sonuç aşağıdaki gibi:
Evet karşımıza çıkan grafiği inceleyelim. İlk kısımda verinin kendisi var. Burayı geçiyoruz. İkinci kısımdaki grafiği ilk yazıda görmüştük. Burada log dönüşümü uyguladık ve verinin varyansını stabilize ettik. Görüldüğü üzere verinin seviyesi değiştikçe artık değişkenlik ilk halindeki kadar değişmiyor. Verinin içerisinde oynadığı mesafe hemen hemen her seviyede birbirine yakın. Son kısımda ise log dönüşümüne ek olarak bir fark alma uygulanmış hali var. Dikkat ederseniz burada artık iyice trend ve mevsimsellik arındırılmış ve verinin hareketleri rastgele hale gelmiş durumda. Bu da aslında veriyi büyük ölçüde durağanlaştırdığımızın bir göstergesi.
Son olarak, durağan serinin bir de ACF grafiğini inceleyelim. Bu grafiğin de veri içerisindeki otokorelasyonları gösterdiğini ilk yazıda işlemiştik. Grafiği çizdirecek olursak:
ggAcf(diff(log(Data_ts), 1),
lag = 30) +
ggtitle("Fark-Otokorelasyon Fonksiyonu Grafiği")
Evet buradaki “ggAcf” komutunu ilk yazıda da kullanmıştık. Log ve fark alınmış veriyi ilk argüman olarak verdikten sonra, 30 gecikmeye kadar grafiği çizdirmesini söyledik. Sonuç:
Grafiğe baktığımızda, beyaz gürültü(white noise) serisinin otokorelasyon grafiğine çok benzediğini görüyoruz. Yani anlamlılık sınırlarının dışında sadece bir tane otokorelasyon katsayısı bulunuyor ve katsayıların gidişatı tamamen rastgele. Birbirleriyle hiçbir ilişki bulunmuyor. Verimizin kendi otokorelasyon grafiğini de getirelim ve karşılaştıralım:
Evet bu da verinin kendisinin otokorelasyon grafiği. Baktığımızda, tamamen birbirleriyle ilişkili ve yavaş yavaş azalan bir grafik görüyoruz. Üstteki ise tamamen rastgele ilerleyen ve bir sonraki adımda nereye yöneleceği belirsiz bir grafik. Bu da aslında durağan bir serinin otokorelasyon fonksiyonunun bu şekilde dağınık ve rastgele bir grafik çizeceğini gösteriyor.
Durağanlık kavramı, ARIMA modellemeyi öğrenmeden önce bilinmesi gereken temellerden biriydi. Şimdi ARIMA’nın içerisine daha detaylı olarak bakalım. Aslında ARIMA ifadesi “AutoRegressive Integrated Moving Average” ifadesinin kısaltmasıdır. Yani Aslında AR() – Auto Regressive(Otoregresif) ve MA() – Moving Average(Hareketli Ortalama) modellerinin bir entegrasyonudur. Şimdi AR() ve MA() modelleri üzerinde durarak ARIMA’nın temelini oluşturan bu ifadeleri anlayalım.
AR() modelleriyle başlayalım. Nasıl ki çoklu regresyonda tahminci değişkenlerin doğrusal bir kombinasyonunu kullanarak bağımlı değişkeni tahminliyorsak, burada da bağımlı değişkenin önceki değerlerinin lineer kombinasyonunu kullanarak bağımlı değişkeni tahminliyoruz. Yani bağımlı değişkeninin kendisini kullanarak regresyon yapıyoruz, “otoregresif” ifadesi de buradan geliyor. Bir AR(p) modeli aşağıdaki gibidir:
Görüleceği üzere burada, bağımlı değişkeni, o değişkenin gecikmeli(önceki) değerleriyle tahminliyoruz ve hepsinin bir katsayısı bulunuyor. Bu durum çoklu regresyon ile aynı. Buradaki “p” ifadesi ise, bu modelde kaç gecikmeye kadar geri gidileceğini gösteriyor. Yani AR(2) dediğimizde, 2 gecikmeye kadar giderek bağımsız değişkenin bu gecikmelerdeki değerlerini oto-regresyon modelinde kullanıyoruz. Buradaki “c” sabit bir sayıyı ifade ediyorken(regresyon denkleminde de “intercept” yani doğrunun başlangıç noktası vardı), epsilon ifadesi de tabii ki modelin hatasını temsil ediyor.
Burada gecikme değerlerinin yani bir diğer ifadeyle tahminci değişkenlerin katsayılarının 0 olduğunu düşünelim. Bu durumda oto-regresyon denkleminde sadece c sabiti ve epsilon kalacaktır. Bu da aslında bizim için beyaz gürültü(white noise) serisini temsil eden denklem demektir. Yani katsayılar sıfır olduğunda geçmiş değerlerin denkleme hiçbir etkisi yok demek oluyor ve bu da serinin tamamen beyaz gürültü olduğunun göstergesi.
İkinci durumda da bir AR(1) modeli düşünelim. Yani sadece bir gecikmedeki değeri tahminci olarak kullanıyoruz. Aynı modelde, c sabitinin değeri 0 ve tahmincinin katsayısı da 1 olursa karşımıza aşağıdaki gibi bir denklem çıkar:
Bu denklem de aslında bir “Random Walk(Rastgele Yürüyüş)” modelini temsil ediyor. Yeri gelmişken bu modeli de kısaca açıklayalım. Denklemden de anlaşılabileceği gibi, bir sonraki değer, bir önceki değer + epsilon yani hata formülünden oluşuyor. Bu da demek oluyor ki bir sonraki adımda ne olacağına dair bir fikrimiz yok. Rastgele yürüyüş modellerinde, uzun süren aşağı veya yukarı doğru trendlerden sonra ani yön değişiklikleri olur ve bu da tahmin edilmeyi çok güçleştirir. Örneğin bu durum finansal verilerde çok görülür. Gelecek gözlemin aşağı veya yukarı doğru olma olasılığı eşit olduğundan, bu tarz verilerde genellikle tahminleme “Naive” metod yani son gözlem değerine eşit olma olarak bulunur. Yani aslına bakarsanız sağlam bir tahminimiz yok elimizde.
“Otoregresif” modellerden sonra, biraz da ARIMA’nın ikinci unsuru olan “Hareketli Ortalama” modellerinden bahsedelim. Bu tarz modellerde de, tahminlemeye çalıştığımız değişkenin geçmiş değerlerini tahminciler olarak kullanmak yerine, geçmiş tahmin hatalarını modelde tahminciler olarak kullanıyoruz. Örneğin bir MA(q) modeli aşağıdaki gibi yazılabilir:
Gördüğünüz gibi geriye dönük “q” adedince geçmiş tahmin hatasını modelde kullanıyor ve geleceği tahminlemeye çalışıyoruz. Burada “c” yine sabit bir sayı iken, E(t) değeri modelimizin normal hatası, diğer epsilon değerleri ise geçmiş tahmin hataları. Bu tarz modellere hareketli ortalama modelleri denmesinin sebebiyse, “Yt” yani tahmin edilmek istenen değerin, geçmiş hata değerlerinin hareketli ortalaması olmasına eşit olmasıdır.
Evet artık ARIMA modelleme kısmına geçebiliriz. Burada da ARIMA modellerini inceledikten sonra, uygulama kısmına geçelim. ARIMA modelleri de, aynı Üstel Düzeltme de olduğu gibi, mevsimsel ve mevsimsel olmayan modeller olarak ayrılabilir. Bizim verimizde mevsimsellik olduğunu bildiğimizden, mevsimsel bir ARIMA modelinin daha uygun olacağını düşünüyoruz. Fakat yine de göstermek amacıyla, mevsimsel olmayan modeli de deneyeceğiz. Mevsimsel ARIMA modellerini aşağıdaki gibi düşünebiliriz:
Görüldüğü üzere, ARIMA modelleri üç ana unsurdan oluşuyor. AR() unsuru, MA() unsuru ve “d” yani fark alma(differencing) unsuru. Buradaki küçük “p” mevsimsel olmayan AR() sırasını gösterirken(kaç gecikmedeki gözlemlerin kullanılacağı) büyük “P” ise mevsimsel olan AR() sırasını gösteriyor(kaç mevsimsel gecikmedeki gözlemlerin kullanılacağı). Aynı şekilde küçük “q” mevsimsel olmayan MA() sırasını gösterirken(kaç gecikmedeki tahmin hatalarının kullanılacağı), büyük “Q” ise mevsimsel olan MA() sırasını gösteriyor(kaç mevsimsel gecikmedeki tahmin hatalarının kullanılacağı). Küçük “d” değeri, veriyi durağan yapana kadar kaç normal fark alma uygulanacağı, büyük “D” değeri ise mevsimsel etkiyi veriden kaldırana kadar kaç mevsimsel fark alma işleminin yapılacağını gösteriyor. Son olarak “m” ise yıl bazında gözlem frekansı. Örneğin çeyrek verisi için 4 olacaktır.
Bu mevsimsel ARIMA modeliydi. Mevsimsel olmayan ARIMA modelini ise aşağıdaki gibi düşünebiliriz:
Görüleceği üzere hem gecikme değerleri hem de geçmiş hata değerleri aynı denklemde birleşmiş durumda. Üsttekinden tek farkı ise mevsimsel kısmın(yani büyük harfli ifadelerin) bulunmaması. Yani mevsimsel olmayan ARIMA modelini ARIMA(p, d, q) olarak düşünebiliriz. Bu denklemde Yt değerlerinin üstündeki kesme işaretleri ise fark alma işlemini ifade ediyor:
Mevsimsel ve mevsimsel olmayan ARIMA modelleri kısaca bu şekilde. Çok fazla detaya girip de yazıyı da boğmaya gerek yok. Bu noktada bahsetmem gereken ek bir kavram daha bulunuyor. O da “Kısmi Otokorelasyon” ifadesi. Daha önce ACF(otokorelasyon grafiği) üzerine çalışmıştık. Fakat ARIMA modellerinde p,d,q ve P,D,Q değerlerini bulabilmek için, sadece zaman serisi grafiğine bakmak çoğu zaman yeterli değildir. Bu yüzden ACF ve PACF(Kısmi oto-korelasyon grafiği) grafiklerine bakmak gerekecektir. Kısmi otokorelasyonu şöyle açıklayabiliriz; Yt ve Yt-k gözlemleri arasındaki otokorelasyonu, aradaki tüm gözlemlerin etkilerini kaldırarak hesaplarsak, bu iki değer arasındaki kısmi otokorelasyon katsayısını bulmuş oluruz. Yani aradakileri kaldırmasaydık, sırf aradakilerin Yt ile ilişkili olmasından dolayı, Yt-k de Yt ile ilişkili gözükebilirdi. Fakat kısmi otokorelasyon katsayısını bulup aradaki gözlemlerin etkilerini kaldırdığımızda, elimizde gerçekten arada ilişki olup olmadığını gösteren ve Yt-k’nin Yt’yi tahminlemede kullanılabilecek bir bilgi barındırıp barındırmadığına dair bir işaret veren bir değer elde etmiş oluyoruz. Burada dikkat çekilmesi gereken bir nokta da, her kısmi otokorelasyon katsayısının bir AR() modelindeki son katsayı olarak tahminlenebileceğidir. Yani yine Yt ve Yt-k arasındaki kısmi otokorelasyon katsayını arıyorsak, bir AR(k) modelinin son katsayısı bize bu değerin tahmini verebilir. Çünkü AR(k) modelinin denklemini göz önüne alırsanız, son katsayı aynı zamanda Yt-k gözleminin Yt üzerindeki etkisini göstermektedir.
Sonuç olarak bu değerleri içeren PACF grafiği de, bizim işimize ARIMA modelinin değerlerini bulmada çok yardımcı olacak.
Peki ACF ve PACF grafilerini kullanarak, p,d,q ve P,D,Q değerlerine nasıl ulaşabiliriz?
Örneğin bir ARIMA(0,0,0)(0,0,1)4 modelini ele alırsak;
-> ACF grafiğinde 4. gecikmede anlamlı bir çubuk görülür, fakat diğer gecikmelerde görülmez.
-> PACF grafiğinde mevsimsel gecikmelerde katsayı çubukları üstel olarak azalır.(4, 8, 12, …)
Bu kez bir ARIMA(0,0,0)(1,0,0)4 modeline bakarsak;
-> ACF grafiğinde mevsimsel gecikmelerde üstel olarak azalma görülür.
-> PACF grafiğinde 4. gecikmede anlamlı bir çubuk görülür fakat diğer gecikmelerde görülmez.
Mevsimsel değerleri(P,D,Q) bulurken, mevsimsel gecikmelere bakmak gerekir. Fakat mevsimsel olmayan değerlerde(p,d,q) mevsimsel olmayan gecikmelere bakacağız. Bu teorik detayları şimdi uygulamaya geçtiğimizde çok rahat anlayacağız.
Tekrar verinin grafiğini gözümüzün önüne getirelim:
Grafikten de net olarak anlaşılıyor ki veri durağan değil. Açık bir yukarı yönlü trend var. Aynı zamanda mevsimsellik göstergeleri de bulunuyor. Buna ek olarak, verinin varyansı da stabil değil ve serinin farklı seviyelerinde farklı büyüklükler gösteriyor. Bu bilgilerden anlaşılıyor ki hem dönüşüm yapacak, hem de fark alma uygulayacağız. Yukarıda hem log dönüşümü uygulanmış hem de fark alınmış halini hesaplamıştık. Tekrar getirelim:
Evet gördüğümüz gibi burada verinin son hali durağan durumda ve varyansı da kısmen stabilize edilmiş. ACF ve PACF’yi görecek olursak;
Data_ts %>% log() %>% diff(lag = 1) %>% ggtsdisplay(xlab = 'Ay', ylab = 'Sayı')
Evet veriye önce log dönüşümü uyguladık, sonra da bir fark alma yaptıktan sonra “ggtsdisplay” komutuyla tüm grafiklerini döktük:
ACF grafiğinde sadece bir tane anlamlı çubuk görülüyor ve çubuklar tamamen rastgele. Burada belki daha da durağanlaştırmak adına bir fark alma işlemi uygulamak daha düşünülebilir ama ben şu an için gerek görmüyorum. Hatırlayalım ki kaç fark almamız gerektiğine dair “ndiffs()” ve “nsdiffs()” komutlarıyla da üstte çalışma yapmıştık. Sonuçlar aşağıdaki gibiydi:
Yani bu komutlar da sadece birinci gecikmede bir adet fark alma işleminin yeterli olacağını söylüyor. Şimdi üstteki ACF ve PACF grafiklerine bakarak, uygun mevsimsel ve mevsimsel olmayan değerleri bulmaya çalışalım.
Unutmayalım ki verimiz aylık olduğu için yıl bazında gözlem frekansı 12. Yani mevsimsel değerleri bulurken 12 ve katlarındaki çubuklara bakacağız. Diğerleri ise mevsimsel olmayan değerleri gösterecek. ACF grafiğinde 1. gecikmedeki anlamlı çubuk ve sonrasında herhangi bir anlamlı çubuk olmaması, bir MA(1) modelinin göstergesi diyebiliriz. Aynı şekilde PACF grafiğinde de 2. çubuktan sonra 10. çubukta da bir anlamlılık var. Biz burada hem AR(2) hem de AR(10) modellerini deneyelim. ACF grafiğindeki 12 ve 36. gecikmelerdeki neredeyse anlamlı olan çubuklar, mevsimsel bir MA(1) modelini tavsiye ediyor olabilir. Aynı şekilde PACF grafiğinde de, ACF’deki kadar belirgin olmasa da, denemeye değer bir mevsimsel AR(1) modeli olabilir çünkü buralarda anlamlılığa yaklaşan çubuklar bulunuyor. Tüm bu verileri topladığımızda, elimizde denemek için ARIMA(2,1,1)(1,0,1)[12] ve ARIMA(10,1,1)(1,0,1)[12] modelleri var. Tabii ki biz sadece bu modelleri değil, bunların küçük değişimleriyle de oluşacak olan modelleri de uygulayacak ve sonuçları göreceğiz. Aşağıdaki kod bloğuyla modellerimizi kuralım:
#ARIMA(2,1,1)(1,0,1)[12]
Model1 <- Arima(y = Data_ts,
order = c(2, 1, 1),
seasonal = list(order = c(1, 0, 1), period = 12),
lambda = 0)
#ARIMA(10,1,1)(1,0,1)[12]
Model2 <- Arima(y = Data_ts,
order = c(10, 1, 1),
seasonal = list(order = c(1, 0, 1), period = 12),
lambda = 0)
#ARIMA(2,1,1)(1,0,0)[12]
Model3 <- Arima(y = Data_ts,
order = c(2, 1, 1),
seasonal = list(order = c(1, 0, 0), period = 12),
lambda = 0)
#ARIMA(2,1,1)(0,0,1)[12]
Model4 <- Arima(y = Data_ts,
order = c(2, 1, 1),
seasonal = list(order = c(0, 0, 1), period = 12),
lambda = 0)
#ARIMA(2,1,1)(0,0,0)[12]
Model5 <- Arima(y = Data_ts,
order = c(2, 1, 1),
seasonal = list(order = c(0, 0, 0), period = 12),
lambda = 0)
#ARIMA(2,1,2)(1,0,1)[12]
Model6 <- Arima(y = Data_ts,
order = c(2, 1, 2),
seasonal = list(order = c(1, 0, 1), period = 12),
lambda = 0)
#ARIMA(2,1,2)(0,0,1)[12]
Model7 <- Arima(y = Data_ts,
order = c(2, 1, 2),
seasonal = list(order = c(0, 0, 1), period = 12),
lambda = 0)
list = list(Model1, Model2, Model3, Model4, Model5, Model6, Model7)
for (i in list) {
print(i$aicc)
}
Üstte bahsettiğimiz modelleri, birkaç değişik versiyonlarıyla beraber oluşturduk. Burada dikkat edelim ki “Arima” komutunun içine argüman olarak veriyi, mevsimsel olmayan değerlerin sırasını, ve mevsimsel değerlerin sırasını verdik. En sondaki “lambda=0” ifadesi de, aslında veriye bir log dönüşümü uyguladığımızı gösteriyor. En son tüm bu modelleri bir listeye atıp ufak bir for döngüsü ile de hepsini karşılaştırmak adına AICc değerlerini çektik:
En küçük AICc değerini veren model, ilk model yani ARIMA(2,1,1)(1,0,1)[12] modeli oldu. Bu da grafiklere bakarak bizim seçtiğimiz modeldi. Şimdi bu modelin kalıntılarını kontrol edip varsayımları sağlayıp sağlamadığına bakalım:
checkresiduals(Model1)
Evet “checkresiduals()” komutunun sonucu aşağıda:
Sonuçları inceleyelim. Grafiğe baktığımızda, ACF grafiğinde neredeyse tüm çubukların anlamlılık sınırlarını aşmadığını görüyoruz. Sadece bir tanesi ufak bir aşım yapmış durumda. Dolayısıyla kalıntılar beyaz gürültü serisine çok yakın bir seri gösteriyor diyebiliriz. Bu da bizim modelimizin kalıntı varsayımını sağladığını gösteriyor. Ek olarak, aşağıda “Ljung-Box” testinin sonucunda elde edilen p değeri de 0.7522. Yani, kalıntılarda geriye kalan herhangi bir otokorelasyon olmadığını ve kalıntıların bir beyaz gürültü serisinden çok da ayrılamadığını gösteriyor. Bu bilgiler de seçtiğimiz modelin kalıntı varsayımlarını sağladığını bize anlatıyor.
Şimdi veriyi train ve test setler olarak ayırıp, üstteki modellerin test set üzerindeki performanslarını inceleyelim. AICc bazında en iyi model ARIMA(2,1,1)(1,0,1)[12] idi, fakat bakalım test set üzerindeki performansta işler değişecek mi.
Veriyi ayırıp aynı modelleri bu kez train set üzerinde kuralım:
train_ts <- window(Data_ts, end = c(1973, 12)) test_ts <- window(Data_ts, start = c(1974, 01)) #ARIMA(2,1,1)(1,0,1)[12] Model1_t <- Arima(y = train_ts, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 1), period = 12), lambda = 0) #ARIMA(10,1,1)(1,0,1)[12] Model2_t <- Arima(y = train_ts, order = c(10, 1, 1), seasonal = list(order = c(1, 0, 1), period = 12), lambda = 0) #ARIMA(2,1,1)(1,0,0)[12] Model3_t <- Arima(y = train_ts, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 0), period = 12), lambda = 0) #ARIMA(2,1,1)(0,0,1)[12] Model4_t <- Arima(y = train_ts, order = c(2, 1, 1), seasonal = list(order = c(0, 0, 1), period = 12), lambda = 0) #ARIMA(2,1,1)(0,0,0)[12] Model5_t <- Arima(y = train_ts, order = c(2, 1, 1), seasonal = list(order = c(0, 0, 0), period = 12), lambda = 0) #ARIMA(2,1,2)(1,0,1)[12] Model6_t <- Arima(y = train_ts, order = c(2, 1, 2), seasonal = list(order = c(1, 0, 1), period = 12), lambda = 0) #ARIMA(2,1,2)(0,0,1)[12] Model7_t <- Arima(y = train_ts, order = c(2, 1, 2), seasonal = list(order = c(0, 0, 1), period = 12), lambda = 0)
İlk iki satırda veriyi train ve test set olarak ayırdıktan sonra, aynı modelleri train set üzerinde kurduk. Şimdi bu modellerden tahminleri “forecast()” komutuyla üretip, “accuracy()” komutuyla da test set performanslarını ölçeceğiz. Unutmayalım ki “Arima()” komutu direkt olarak tahmin döndürmediği için, ekstra bir adım daha uygulayarak tahminleri manuel olarak hesaplıyoruz.
Aşağıdaki kod bloğuyla gerçekleştirelim:
#Tahminlerin elde edilmesi
Fit1_t <- forecast(Model1_t, h = 22)
Fit2_t <- forecast(Model2_t, h = 22)
Fit3_t <- forecast(Model3_t, h = 22)
Fit4_t <- forecast(Model4_t, h = 22)
Fit5_t <- forecast(Model5_t, h = 22)
Fit6_t <- forecast(Model6_t, h = 22)
Fit7_t <- forecast(Model7_t, h = 22)
list2 = list(Fit1_t, Fit2_t, Fit3_t, Fit4_t, Fit5_t, Fit6_t, Fit7_t)
for (i in list2) {
print(accuracy(i, test_ts))
}
Evet yine tüm tahminleri bir listeye attık ve basit bir for döngüsüyle test set üzerindeki performanslarını hesapladık. Sonuç:
Test set üzerinde RMSE bazında hataları inceledik. Diğer hata metrikleri bazında incelemeler de yapılabilir. Yine en düşük hata değerini, dolayısıyla en iyi performansı bizim seçmiş olduğumuz ARIMA(2,1,1)(1,0,1)[12] modeli verdi. Ona çok yakın olan ve iyi performans gösteren diğer model ise sondan birinci kurduğumuz ARIMA(2,1,2)(1,0,1)[12] modeli. Fakat biz tahminlerimizi tabii ki en iyi model üzerinden elde edeceğiz.
Test üzerindeki performansı bir de grafiğe dökelim:
autoplot(Data_ts) +
autolayer(Fit1_t, series="ARIMA(2,1,1)(1,0,1)[12]", PI = FALSE) +
ggtitle("Model Tahminleri") + xlab("Yıl") +
ylab("Silahlı Soygunlar") +
guides(colour=guide_legend(title="Method"))
Sonuca bakarsak:
Evet seçtiğimiz modelin test set üzerindeki performansını görüyoruz. Unutmamak gerekir ki bazı modeller dışında, modellerin hata metrikleri de birbirine çok yakın. Yani diğer varyasyonlardan birini de seçseydik, buna benzer bir sonuç elde edebilirdik. Fakat az farkla da olsa elimizdeki en iyi sonucu veren model bu. Bir de bu model verinin gelecek değerlerini nasıl tahmin ediyor onu görelim:
Model1 %>% forecast(h=22) %>% autoplot()
Dikkat edelim ki burada Model1’i kullandık. Bu, train set ile değil verinin kendisi üzerinde kurduğumuz modeldi. Fakat tabii ki model düzeni aynı. Sonucu görelim:
Modelin veri üzerindeki tahminleri bu şekilde. Dikkat ederseniz güven aralıkları çok geniş. Bu da aslında ileride verinin geniş bir mesafede oynayabileceğini gösteriyor. %80 ve %95 güven aralıklarını görebiliyoruz. Modelin içeriğini ve ürettiği sayısal nokta tahminlerini de görelim:
Model1 Model1 %>% forecast(h=22) -> Fit Fit$mean
Modelin ürettiği tahminleri “Fit” nesnesine atadık ve “mean” metoduyla nokta tahminlerini çektik. Sonuçları görelim:
Model bilgisine baktığımızda, AR() ve MA() modellerinin katsayılarını görebiliyoruz. Buradaki katsayılar tabii ki üstteki denklemlerde gösterdiğimiz gecikme gözlemlerinin ve gecikme tahmin hatalarının katsayıları. Bunlara ek olarak, AIC, AICc ve BIC değerleri de gözüküyor.
Tahminlere geçtiğimizdeyse, ay ay tüm tahmin değerlerini görebiliyoruz. Bunlar yukarıda güven aralıklarıyla beraber grafiğini çizdirdiğimiz değerler.
Aslına bakarsanız ARIMA modelleme ve tahminleme işlemini bitirdik ve gerekli detaylara değindik. Fakat yazıyı bitirmeden önce, bahsetmek istediğim son bir şey daha var. O da üstte adım adım gerçekleştirdiğimiz prosedürü R üzerinde otomatik olarak gerçekleştiren “auto.arima” fonksiyonu. Fakat unutmamak gerekir ki otomatik prosedürler her zaman en doğru sonucu getirmez. Bu yüzden manuel olarak üstteki adımları gerçekleştirmek çok daha güvenilir. Şimdi aşağıdaki kod bloğuyla bu fonksiyonu çalıştırıp veri için hangi modeli seçeceğini görelim:
Fit_auto <- auto.arima(y = Data_ts, seasonal = TRUE, stepwise = FALSE, approximation = FALSE) Fit_auto
Kod bloğunu inceleyelim. “auto.arima” fonksiyonunu kullanarak otomatik bir ARIMA modeli kuruyoruz ve burada veriyi argüman olarak verdikten sonra, “seasonal=TRUE” ifadesiyle algoritmanın mevsimsel modelleri de gözden geçirmesini söylüyoruz. “approximation = FALSE” ifadesi, algoritmanın araştırmayı hızlandırması için yaklaşımlar kullanmasını iptal ediyor. Daha detaylı ve yavaş bir şekilde araştırmasını söylüyoruz. Son olarak, “stepwise = FALSE” argümanı ise, çok daha fazla modelin araştırılmasını sağlıyor. Sonuca bakalım:
Aloritma bizim için ARIMA(0,1,4)(0,0,1)[12] modelini drift ekleyerek seçti. Bu bizim kendimiz seçtiğimiz modelden farklı. AICc değeri 1177.64. Şimdi bu modelin test set üzerindeki performansını inceleyelim. Bakalım sonuçlar bizim seçtiğimiz modelden iyi mi.
Fit_auto_train <- Arima(y = train_ts, order = c(0, 1, 4), seasonal = list(order = c(0, 0, 1), period = 12), lambda = 0) Fit_auto_forecasts % forecast(h = 22) accuracy(Fit_auto_forecasts, test_ts)
Algoritmanın seçtiği düzende ARIMA modelini train set üzerinde kurduk ve test seti üzerinde de yaptığı tahminlerin performansını karşılaştırdık. Sonucu görelim:
Evet modelin RMSE değeri 73.28. Bizim seçtiğimiz modelinki ise 68.74 idi. Görüleceği gibi bizim seçtiğimiz model algoritmanın otomatik olarak seçtiği modelden az da olsa daha iyi performans göstermiş durumda. Tabii ki ortada büyük bir fark yok. Son olarak, bu modelin de hem test set hem de gelecek tahminlerini görelim:
autoplot(Data_ts) +
autolayer(Fit_auto_forecasts, series="ARIMA(0,1,4)(0,0,1)[12]", PI = FALSE) +
ggtitle("Model Tahminleri") + xlab("Yıl") +
ylab("Silahlı Soygunlar") +
guides(colour=guide_legend(title="Method"))
Fit_auto %>% forecast(h = 22) %>% autoplot()
Sonuçlar:
Evet test set üzerindeki performans kötü görünüyor. Gelecek gözlemlerin tahminine bakarsak, tahminlerin bir yerden sonra yükselen bir trendle düz gittiğini görebiliriz. Bu da aslında modeldeki drift etkisinden dolayı.
Bizim modelimiz ve otomatik modelin nokta tahminlerini de karşılaştıralım:
Fit_auto_f % forecast(h = 22)
cbind('Manuel Model Tahminleri' = Fit$mean,
'Otomatik Model Tahminleri' = Fit_auto_f$mean)
“cbind” komutu ile matriks haline getirik ve sonuç:
Tabloyu incelersek, sonlara doğru bizim modelimizin tahminleri daha üstel biçimde artarken, otomatik modelin tahminleri ağır bir biçimde düz olarak artıyor. Bu da bize tekrar drift etkisini gösteriyor.
Evet, serinin üçüncü ve büyük ihtimalle son yazısında, zaman serileri için çok önemli olan bir metod olan ARIMA modelleri üzerinde durduk ve iyi tahminler elde ettik. Bu yöntem bu konu için çok kritiktir ve sıkça kullanılır. İnternette de birçok kaynaktan detaylı araştırmalar yapabilirsiniz.
NOTLAR
- ARIMA metodları tahminleme işleminde sıkça kullanılır ve başarılı tahminler üretir. Fakat üzerinde sıkı çalışmak gerekmektedir.
- Manuel seçtiğimiz modelin birkaç versiyonunu denedik ve hepsinin başarısını inceledik. Burada denenen model sayısı artırılıp başka modeller de test edilebilir. Fakat aşağı yukarı aynı sonuçlar elde edilecektir. Bizim seçtiğimiz model, yüksek ihtimalle yine en iyisi çıkacaktır.
- ARIMA ve Üstel Düzeltme metodları en yaygın iki metod olduğu için, bunların karşılaştırılmasını inceleyip veri için en uygun olanın hangisi olduğuna karar verme konusunda başka kaynaklar okunması da oldukça önemlidir.
REFERANSLAR
fpp2.pdf erişimi için tıklayın