Regresi Deret Waktu

package

library(dplyr)
library(readxl)
library(forecast)
library(TTR)
library(tseries)
library(lmtest) #uji-Durbin Watson
library(orcutt) #Cochrane-Orcutt
library(HoRM)#Hildreth Lu

Data

#working directory
setwd("D:/MATERI KULIAH S2 IPB/ASPRAK 1/RESPONSI 3")

#membuka file data
dataregresi <- read.csv("dataregresi.csv")
head(dataregresi)

##    y   x
## 1 32  38
## 2 49  40
## 3 50  44
## 4 39  62
## 5 38  50
## 6 55 106

Eksplorasi Data

x <- dataregresi$x
y <- dataregresi$y

#diagram pencar identifikasi model
plot(x,y,pch = 20, col = "blue", main = "Scatter Plot X vs Y",
     ylab = "Nilai Peubah Y", xlab = "Nilai Peubah X")

#korelasi x dan y
cor(x,y)

## [1] 0.6090332

Model Regresi Deret Waktu

#model regresi
model <- lm(y~x, data = dataregresi)
summary(model)

## 
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.9921  -6.0457  -0.9104   5.4266  21.1712 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.04340    4.08925  10.281 7.27e-10 ***
## x            0.20918    0.05808   3.602  0.00159 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.6 on 22 degrees of freedom
## Multiple R-squared:  0.3709, Adjusted R-squared:  0.3423 
## F-statistic: 12.97 on 1 and 22 DF,  p-value: 0.001585

Deteksi Autokorelasi

Grafik

1. Residual Plot

#sisaan dan fitted value
resi1 <- residuals(model)
fit <- predict(model)

#Diagnostik dengan eksploratif
par(mfrow = c(2,2))

qqnorm(resi1)
qqline(resi1, col = "steelblue", lwd = 2)

plot(fit, resi1, col = "steelblue", pch = 20, xlab = "Sisaan", 
     ylab = "Fitted Values", main = "Sisaan vs Fitted Values")
abline(a = 0, b = 0, lwd = 2)

hist(resi1, col = "steelblue")

plot(seq(1,24,1), resi1, col = "steelblue", pch = 20, 
     xlab = "Sisaan", ylab = "Order", main = "Sisaan vs Order")
lines(seq(1,24,1), resi1, col = "red")
abline(a = 0, b = 0, lwd = 2)

Bila plot tidak membentuk pola tertentu, maka asumsi kebebasan terpenuhi

2. ACF dan PACF Plot

#ACF dan PACF identifikasi autokorelasi
par(mfrow = c(2,1))
acf(resi1)
pacf(resi1)

Bila ACF dan PACF tidak ada yang signifikan, maka sisaan saling bebas

Uji Statistik

1. Durbin Watson Test

H0: tidak ada autokorelasi

H1: ada autokorelasi

lmtest::dwtest(model,alternative = 'two.sided') #ada autokorelasi

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.2088, p-value = 0.02727
## alternative hypothesis: true autocorrelation is not 0

p-value < 5% : Tolak H0 sehingga cukup bukti untuk menyatakan bahwa terdapat autokorelasi positif pada taraf 5%

2. Breusch-Godfrey Test

H0: tidak ada autokorelasi

H1: ada autokorelasi

lmtest::bgtest(y ~ x, data=dataregresi, order=1) #Perform Breusch-Godfrey test for first-order serial correlation

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  y ~ x
## LM test = 2.5821, df = 1, p-value = 0.1081

p-value > 5% : Tidak Tolak H0 sehingga dapat disimpulkan bahwa tidak terdapat autokorelasi pada taraf 5%

Run’s Test

#install.packages("lawstat")
library(lawstat)
lawstat::runs.test(resi1,alternative="two.sided")

## 
##  Runs Test - Two sided
## 
## data:  resi1
## Standardized Runs Statistic = -0.83485, p-value = 0.4038

Penanganan Aurokorelasi

1. Cochrane-Orcutt

#Interactive method using to solve first order autocorrelation problems. 

modelco <- orcutt::cochrane.orcutt(model)
modelco

## Cochrane-orcutt estimation for first order autocorrelation 
##  
## Call:
## lm(formula = y ~ x, data = dataregresi)
## 
##  number of interaction: 13
##  rho 0.441367
## 
## Durbin-Watson statistic 
## (original):    1.20877 , p-value: 1.364e-02
## (transformed): 1.66348 , p-value: 1.992e-01
##  
##  coefficients: 
## (Intercept)           x 
##   47.908320    0.132056

#rho optimum
rho <- modelco$rho
y

##  [1] 32 49 50 39 38 55 57 50 58 81 81 67 69 64 60 51 47 46 40 49 72 60 54 40

y[-1]

##  [1] 49 50 39 38 55 57 50 58 81 81 67 69 64 60 51 47 46 40 49 72 60 54 40

y[-24]

##  [1] 32 49 50 39 38 55 57 50 58 81 81 67 69 64 60 51 47 46 40 49 72 60 54

#transformasi terhadap y dan x
(y.trans <- y[-1]-y[-24]*rho)

##  [1] 34.87624 28.37300 16.93163 20.78667 38.22804 32.72479 24.84206 35.93163
##  [9] 55.40069 45.24924 31.24924 39.42838 33.54565 31.75249 24.51795 24.49026
## [17] 25.25573 19.69710 31.34530 50.37300 28.22155 27.51795 16.16616

(x.trans <- x[-1]-x[-24]*rho)

##  [1]  23.2280380  26.3453032  42.5798335  22.6352199  83.9316289   3.2150534
##  [7]  29.9316289 109.0488941  79.7395004  39.0912959  51.8632579  31.6287276
## [13]  99.3388108  -6.2604996   9.0488941  41.8762425 -10.7165756  11.8208561
## [19]  28.0553864  26.1107728  -0.5374317  34.0553864  -8.5374317

#model baru
modelcorho <- lm(y.trans~x.trans)
summary(modelcorho)

## 
## Call:
## lm(formula = y.trans ~ x.trans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.454  -6.105  -1.869   5.291  20.162 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.76315    2.74122   9.763 2.94e-09 ***
## x.trans      0.13206    0.05898   2.239   0.0361 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.15 on 21 degrees of freedom
## Multiple R-squared:  0.1927, Adjusted R-squared:  0.1543 
## F-statistic: 5.013 on 1 and 21 DF,  p-value: 0.03612

lmtest::dwtest(modelcorho,alternative = 'two.sided')

## 
##  Durbin-Watson test
## 
## data:  modelcorho
## DW = 1.6635, p-value = 0.3985
## alternative hypothesis: true autocorrelation is not 0

p-value > 5% : Gagal Tolak H0 sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%

2. Hildreth-lu

# Hildreth-Lu
hildreth.lu.func<- function(r, model){
  x <- model.matrix(model)[,-1]
  y <- model.response(model.frame(model))
  n <- length(y)
  t <- 2:n
  y <- y[t]-r*y[t-1]
  x <- x[t]-r*x[t-1]
  
  return(lm(y~x))
}

#mencari rho yang meminimumkan SSE (iteratif)
r <- c(seq(0.1,0.8, by= 0.1), seq(0.9,0.99, by= 0.01))
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,40 memiliki SSE Terkecil

##     rho      SSE
## 1  0.10 1979.103
## 2  0.20 1869.163
## 3  0.30 1796.786
## 4  0.40 1761.665
## 5  0.50 1765.243
## 6  0.60 1810.768
## 7  0.70 1902.698
## 8  0.80 2045.622
## 9  0.90 2243.228
## 10 0.91 2266.095
## 11 0.92 2289.534
## 12 0.93 2313.546
## 13 0.94 2338.132
## 14 0.95 2363.293
## 15 0.96 2389.029
## 16 0.97 2415.341
## 17 0.98 2442.231
## 18 0.99 2469.697

#grafik rho dan SSE
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

#rho optimal di sekitar 0.4
r <- seq(0.3, 0.5, by= 0.01)
tab <- data.frame("rho" = r, "SSE" = sapply(r, function(i){deviance(hildreth.lu.func(i, model))}))
round(tab, 4) #0,44 memiliki SSE Terkecil

##     rho      SSE
## 1  0.30 1796.786
## 2  0.31 1791.588
## 3  0.32 1786.762
## 4  0.33 1782.309
## 5  0.34 1778.229
## 6  0.35 1774.523
## 7  0.36 1771.193
## 8  0.37 1768.241
## 9  0.38 1765.667
## 10 0.39 1763.475
## 11 0.40 1761.665
## 12 0.41 1760.241
## 13 0.42 1759.205
## 14 0.43 1758.560
## 15 0.44 1758.307
## 16 0.45 1758.452
## 17 0.46 1758.996
## 18 0.47 1759.944
## 19 0.48 1761.298
## 20 0.49 1763.063
## 21 0.50 1765.243

#grafik SSE optimum
plot(tab$SSE ~ tab$rho , type = "l")
abline(v = tab[tab$SSE==min(tab$SSE),"rho"], lty = 3)

# Model terbaik
modelhl <- hildreth.lu.func(0.44, model)
summary(modelhl)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.460  -6.101  -1.872   5.278  20.160 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.82008    2.74456   9.772  2.9e-09 ***
## x            0.13228    0.05897   2.243   0.0358 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.15 on 21 degrees of freedom
## Multiple R-squared:  0.1933, Adjusted R-squared:  0.1549 
## F-statistic: 5.031 on 1 and 21 DF,  p-value: 0.03581

# Deteksi autokorelasi
lmtest::dwtest(modelhl)

## 
##  Durbin-Watson test
## 
## data:  modelhl
## DW = 1.6625, p-value = 0.1984
## alternative hypothesis: true autocorrelation is greater than 0

p-value > 5% : Gagal Tolak H0 sehingga belum cukup bukti untuk menyatakan bahwa terdapat autokorelasi pada taraf 5%

# Transformasi Balik
cat("y = ", coef(modelhl)[1]/(1-0.44), "+", coef(modelhl)[2],"x", sep = "") #persamaan regresi setelah di transformasi ke persamaan awal

## y = 47.893+0.1322756x

#fungsi Hildreth Lu dengan library HoRM
modelhl2 <- HoRM::hildreth.lu(y, x, 0.44) #fungsi hildreth lu dengan rho tertentu
summary(modelhl2)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.460  -6.101  -1.872   5.278  20.160 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.82008    2.74456   9.772  2.9e-09 ***
## x            0.13228    0.05897   2.243   0.0358 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.15 on 21 degrees of freedom
## Multiple R-squared:  0.1933, Adjusted R-squared:  0.1549 
## F-statistic: 5.031 on 1 and 21 DF,  p-value: 0.03581