Durağan

library(WDI)
df = WDI(indicator=c(un='SL.UEM.TOTL.ZS', enf='FP.CPI.TOTL.ZG' ), country=c('TR'),  start=1992, end=2020)
library(dynlm)
## Zorunlu paket yükleniyor: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Attaching package: 'zoo'
df.ts <- ts(df, start=c(1992), end=c(2020),frequency=1)
df.ts<-df.ts[,c("un","enf")]
plot(df.ts[,"un"], ylab="işsizlik")

plot(df.ts[,"enf"], ylab="enflasyon")

gecikmeliun <- data.frame(cbind(df.ts[,"un"], lag(df.ts[,"un"],-1)))
names(gecikmeliun) <- c("un","un1gec")
head(gecikmeliun)
##     un un1gec
## 1 8.51     NA
## 2 8.96   8.51
## 3 8.58   8.96
## 4 7.64   8.58
## 5 6.63   7.64
## 6 6.84   6.63
acf(df.ts[,"un"])

acf(df.ts[,"enf"])

Enflasyon ise 5’ci gecikmeye kadar pozitif otokorelasyon gösterir. işsizliğin ilk farkı alınmış serisini yaratalım ve Phillips eğrisi regresyonunu yapalım. \[inf=β0+β1Δun+u\]

enflasyon <- df.ts[,"enf"]
deltaissizlik <- diff(df.ts[,"un"])
plot(enflasyon)

plot(deltaissizlik)

phillips.reg <- dynlm(enflasyon~deltaissizlik)
summary(phillips.reg)
## 
## Time series regression with "ts" data:
## Start = 1993, End = 2020
## 
## Call:
## dynlm(formula = enflasyon ~ deltaissizlik)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -30.44 -24.05 -20.79  27.40  71.41 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     33.002      6.331   5.213 1.92e-05 ***
## deltaissizlik   -2.103      5.163  -0.407    0.687    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.2 on 26 degrees of freedom
## Multiple R-squared:  0.00634,    Adjusted R-squared:  -0.03188 
## F-statistic: 0.1659 on 1 and 26 DF,  p-value: 0.6871
uhat <- resid(phillips.reg)
plot(uhat)

Hata payının zamana göre değişiklik gösterdiğini görüyoruz. Bu durum hata payı sabit ortalama değere sahip olmalı varsayımımıza aykırı. Bu yüzden anlamlılık düzeyini ölçmek için bulunan hata payları da yanlış.

acf(uhat)

Hata payı beşinci gecikmesine kadar otokorelasyona sahip. Atokorelasyonu ölçmek için kullandığımız testlerden biri Breusch-Godfrey testidir. Örneğimizi ele alalım. \[inf=β0+β1Δun+u\] Hata payı u’nun otokorelasyona sahip olduğunu düşünüyoruz. Yani hata payının gecikmeli haliyle korelasyonu varsa, gecikmeli halinin önündeki katsayılar anlamlı ve sıfırdan farklı olmalı. \[utα0+α1ut-1+ε\] bu durumda eğer sıfıra eşit değilse. ve arasında korelasyon mevcuttur diyebiliriz. Bu yüzde hipotez \[H0:α1=0\] Eğer daha çoklu geciklemeri test etmek istersek F testi veya Chi kare istatitiğini kullanabiliriz. R’da Breusch-Godfrey testi şu şekilde yapılabilir.

library(lmtest)
bgtest(phillips.reg, order=1, type="F")
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  phillips.reg
## LM test = 131.8, df1 = 1, df2 = 25, p-value = 1.844e-11
library(lmtest)
bgtest(phillips.reg, order=1, type="Chisq")
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  phillips.reg
## LM test = 23.536, df = 1, p-value = 1.226e-06

Bu iki testte de olaılık değerlerinin (p-value) çok küçük olduğunu görüyoruz. Bu hipotezimizin doğru olma ihtimalinin çok düşük olduğunu gösteriyor. Dolayısıyla otokorelasyonu yoktur hipotezini tamin ettiğimiz reddediyoruz.

Otokorelasyon mevcut olduğu için, regresyonun gerçek anlamlılık düzeylerini bulmak için başka bir hata payı değeri kullanmalıyız. Bunun için heteroskadastisite ve otokorelasyon düzelten (heteroskedasticity and autocorrelation consistent, HAC) standart hatalar kullanacağız. Bu aynı zamanda Newey-West hata payı olarak bilinir.

library(sandwich)
coeftest(phillips.reg, vcov.=vcovHAC(phillips.reg))
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    33.0018    10.0679  3.2779 0.002968 **
## deltaissizlik  -2.1029     4.9319 -0.4264 0.673334   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(phillips.reg, vcov.=NeweyWest(phillips.reg))
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    33.0018    28.0141  1.1780   0.2495
## deltaissizlik  -2.1029     9.0586 -0.2321   0.8182
coeftest(phillips.reg, vcov.=kernHAC(phillips.reg))
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    33.0018    30.5887  1.0789   0.2905
## deltaissizlik  -2.1029     8.8534 -0.2375   0.8141

Farklı varyans kovaryans matriksleri kullanıldığında, katsayıların değişmediğini, sadece standart hataların değiştiğini gözlemleyebilirsiniz.

Durağan Olmayan

library(WDI)
gsyh <- WDI(country=c("US", "TR","JP"), indicator=c("NY.GDP.PCAP.CD"), start=1960, end=2020)
names(gsyh) <- c("iso2c", "Ülke", "KisiBasiGSYH", "Sene")
head(gsyh)
##   iso2c  Ülke KisiBasiGSYH Sene
## 1    JP Japan     40193.25 2020
## 2    JP Japan     40777.61 2019
## 3    JP Japan     39808.17 2018
## 4    JP Japan     38891.09 2017
## 5    JP Japan     39400.74 2016
## 6    JP Japan     34960.64 2015
library(ggplot2)
ggplot(gsyh, aes(Sene, KisiBasiGSYH, color=Ülke, linetype=Ülke)) + geom_line()

TR <- cbind(gsyh$KisiBasiGSYH[gsyh$Ülke == "Turkey"], gsyh$Sene[gsyh$Ülke == "Turkey"])
TR <- TR[order(TR[,2]),]
TR
##             [,1] [,2]
##  [1,]   509.4240 1960
##  [2,]   283.8283 1961
##  [3,]   309.4466 1962
##  [4,]   350.6630 1963
##  [5,]   369.5835 1964
##  [6,]   386.3581 1965
##  [7,]   444.5495 1966
##  [8,]   481.6937 1967
##  [9,]   526.2135 1968
## [10,]   571.6178 1969
## [11,]   489.9304 1970
## [12,]   455.1049 1971
## [13,]   558.4209 1972
## [14,]   686.4901 1973
## [15,]   927.7992 1974
## [16,]  1136.3756 1975
## [17,]  1275.9566 1976
## [18,]  1427.3718 1977
## [19,]  1549.6444 1978
## [20,]  2079.2203 1979
## [21,]  1564.2472 1980
## [22,]  1579.0738 1981
## [23,]  1402.4064 1982
## [24,]  1310.2557 1983
## [25,]  1246.8245 1984
## [26,]  1368.4017 1985
## [27,]  1510.6763 1986
## [28,]  1705.8944 1987
## [29,]  1745.3649 1988
## [30,]  2021.8595 1989
## [31,]  2794.3505 1990
## [32,]  2735.7076 1991
## [33,]  2842.3700 1992
## [34,]  3180.1876 1993
## [35,]  2270.3373 1994
## [36,]  2897.8666 1995
## [37,]  3053.9472 1996
## [38,]  3144.3857 1997
## [39,]  4499.7375 1998
## [40,]  4116.1706 1999
## [41,]  4337.4780 2000
## [42,]  3142.9210 2001
## [43,]  3687.9561 2002
## [44,]  4760.1040 2003
## [45,]  6101.6321 2004
## [46,]  7456.2961 2005
## [47,]  8101.8569 2006
## [48,]  9791.8825 2007
## [49,] 10941.1721 2008
## [50,]  9103.4741 2009
## [51,] 10742.7750 2010
## [52,] 11420.5555 2011
## [53,] 11795.6335 2012
## [54,] 12614.7816 2013
## [55,] 12157.9904 2014
## [56,] 11006.2795 2015
## [57,] 10894.6034 2016
## [58,] 10589.6677 2017
## [59,]  9454.3484 2018
## [60,]  9121.5152 2019
## [61,]  8536.4333 2020
TR <- ts(TR[,1], start=min(gsyh$Sene), end=max(gsyh$Sene))
plot(TR, ylab="Kişi başı GSYH", xlab="Sene")

acf(TR)

pacf(TR)

plot(TR, ylab="Kişi başı GSYH", xlab="Sene")

Bu otokorelasyon ve partial otokorelasyon fonksiyonunu regresyon formunda gösterelim. \[yt=α0+α1yt-1+u\]

alpha1 y’nin birinci gecikmesiyle arasında olan otokorelasyonu gösteriyor. Bu regresyonu örneğimiz için yapalım ve alpha katsayılarını gözlemleyelim.

library(dynlm)
Ilkgecikme <- dynlm(TR ~ L(TR, 1))
summary(Ilkgecikme)
## 
## Time series regression with "ts" data:
## Start = 1961, End = 2020
## 
## Call:
## dynlm(formula = TR ~ L(TR, 1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1903.5  -229.7   -58.7   223.2  1596.5 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 172.73437  122.19650   1.414    0.163    
## L(TR, 1)      0.99022    0.02182  45.383   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 665.2 on 58 degrees of freedom
## Multiple R-squared:  0.9726, Adjusted R-squared:  0.9721 
## F-statistic:  2060 on 1 and 58 DF,  p-value: < 2.2e-16
Ikincigecikme <- dynlm(TR ~ L(TR, 2))
summary(Ikincigecikme)
## 
## Time series regression with "ts" data:
## Start = 1962, End = 2020
## 
## Call:
## dynlm(formula = TR ~ L(TR, 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1719.8  -492.9  -260.8   275.6  2637.7 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 363.85199  180.93530   2.011   0.0491 *  
## L(TR, 2)      0.97998    0.03277  29.905   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 984.6 on 57 degrees of freedom
## Multiple R-squared:  0.9401, Adjusted R-squared:  0.939 
## F-statistic: 894.3 on 1 and 57 DF,  p-value: < 2.2e-16
Ucuncugecikme <- dynlm(TR ~ L(TR, 3))
summary(Ucuncugecikme)
## 
## Time series regression with "ts" data:
## Start = 1963, End = 2020
## 
## Call:
## dynlm(formula = TR ~ L(TR, 3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2281.4  -587.7  -239.6   450.4  3316.5 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 571.80137  234.78995   2.435   0.0181 *  
## L(TR, 3)      0.96755    0.04325  22.371   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1277 on 56 degrees of freedom
## Multiple R-squared:  0.8994, Adjusted R-squared:  0.8976 
## F-statistic: 500.4 on 1 and 56 DF,  p-value: < 2.2e-16

Bu örnekler bu şekilde çoğaltılabilir. Katsayı hep bire yakın ve anlamlılık düzeyi çok yüksek.

Şimdi bu katsayıları birbirlerini kontrol ederek yazalım. İlk yazdığımız regresyon AR(1) olarak adlandırılıyor. Sadece bir gecikme ile yapılan regresyon. \[yt=α0+α1yt-1+ut\]

Benzer bir regresyonu diğer gecikmeleri sabit tutarak katsayının anlamlığını ölçmek için AR(p) gecikmesine kadar yazabiliriz. \[yt=α0+α1yt-1+α2yt-2+α3yt-3+α4yt-4+...+αpyt-p+ut\]

Böyle bir regresyonu gecikme 10’a kadar yapalım, yani AR(10)

AR10 <- dynlm(TR ~ L(TR, c(1:10)))
summary(AR10)
## 
## Time series regression with "ts" data:
## Start = 1970, End = 2020
## 
## Call:
## dynlm(formula = TR ~ L(TR, c(1:10)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2124.75  -277.50   -56.73   337.91  1423.78 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      236.40345  158.63707   1.490   0.1440    
## L(TR, c(1:10))1    1.03427    0.15747   6.568 7.53e-08 ***
## L(TR, c(1:10))2    0.06830    0.22724   0.301   0.7653    
## L(TR, c(1:10))3    0.12906    0.22957   0.562   0.5771    
## L(TR, c(1:10))4   -0.39876    0.23226  -1.717   0.0937 .  
## L(TR, c(1:10))5    0.25626    0.24058   1.065   0.2932    
## L(TR, c(1:10))6    0.01197    0.24336   0.049   0.9610    
## L(TR, c(1:10))7    0.02645    0.23639   0.112   0.9115    
## L(TR, c(1:10))8   -0.06448    0.24204  -0.266   0.7913    
## L(TR, c(1:10))9   -0.02921    0.24118  -0.121   0.9042    
## L(TR, c(1:10))10  -0.11036    0.19298  -0.572   0.5706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 725 on 40 degrees of freedom
## Multiple R-squared:  0.9734, Adjusted R-squared:  0.9668 
## F-statistic: 146.7 on 10 and 40 DF,  p-value: < 2.2e-16

Farkındaysanız aynı partial otokorelasyon (pacf) grafiğinde gösterdiği gibi diğer gecikmeleri kontrol edince birtek ilk gecikme anlamlılık düzeyine sahip oluyor, ve katsayı 1’e çok yakın. Bu katsayının doğru olduğunu kabul edersek regresyon fonksiyonumuz şu şekilde yazılabilir. \[yt=α0+yt-1+ut\] α1 1’e eşit olduğu için katsayının yerine bir yazabildik.yt-1’i sol tarafa atarsak \[yt-yt-1=α0+ut\] \[Δyt=α0+ut\] olarak yazabiliriz. Yani sabit bir sayı ve rastgele hata payı. O zaman α1 1’e şit olan bir zaman serisinin farkını alırsak durağan bir seriye çevirebiliriz diyebiliriz.

Random Walk Örnegi

Serilerin 200 gözlemi olsun.

n<-200
u <- ts(rnorm(n))
v <- ts(rnorm(n))
y <- ts(rep(0,n))
for (t in 2:n){
  y[t]<- y[t-1]+u[t]
}
x <- ts(rep(0,n))
for (t in 2:n){
  x[t]<- x[t-1]+v[t]
}
plot(y,type='l', ylab="y[t-1]+u[t]")

plot(x,type='l', ylab="x[t-1]+v[t]")

Spurious <- lm(y~x)
summary(Spurious)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4607  -1.9811   0.6854   2.1456   7.7723 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.60072    0.29783 -15.447   <2e-16 ***
## x            0.01722    0.05486   0.314    0.754    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.61 on 198 degrees of freedom
## Multiple R-squared:  0.0004974,  Adjusted R-squared:  -0.004551 
## F-statistic: 0.09853 on 1 and 198 DF,  p-value: 0.7539
duragan <- dynlm(d(y) ~ d(x))
summary(duragan)
## 
## Time series regression with "ts" data:
## Start = 2, End = 200
## 
## Call:
## dynlm(formula = d(y) ~ d(x))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.57857 -0.59583 -0.06965  0.65648  2.80468 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.06761    0.06482  -1.043    0.298
## d(x)         0.02091    0.06220   0.336    0.737
## 
## Residual standard error: 0.914 on 197 degrees of freedom
## Multiple R-squared:  0.0005732,  Adjusted R-squared:  -0.0045 
## F-statistic: 0.113 on 1 and 197 DF,  p-value: 0.7371