library(ggplot2)
library(dplyr)
library(dygraphs)
library(xts)
library(plotly)
library(TSA)
library(hrbrthemes)
# Deklarasi Konstanta
NSim <- 200 #NSim banyaknya simulasi
rataan <- 0 #besar rataan untuk error
sig2 <- 0.5 #besar variansi untuk error
#inisialisasi peubah acak Yt
data <- data.frame(
index = seq(1:NSim),
Yt_p = rep(0, NSim),
et_p = rep(0, NSim)
)
# Mengisi nilai Random Walk
for (i in 1:(NSim-1)) {
et_p <- rnorm(1,rataan,sig2)
data[i+1,2] <- data[i,1]+ et_p
data[i+1,3] <- et_p
}
p <- data %>% ggplot( aes(x=index, y=Yt_p)) +
geom_line(color="#69b3a2") +
ylab("Nilai Y_t") +
theme_ipsum()
p <- ggplotly(p)
p
#Menghitung Kovariansi
maxLag <- 16
time <- seq(from = as.POSIXlt("2002-01-01"), to = as.POSIXct("2018-08-30"), by = "month")
df <- data.frame(time, data) #membuat data frame baru, time+data awal
#Data Frame Kovariansi
cov <- data.frame(k = c(1:maxLag),
acf = NA,
stringsAsFactors = FALSE)
for (k in 1: maxLag) {
N <- nrow(df) - k
matkov <- matrix(c(data[1:N,2], data[(1+k):(N+k),2]), nrow = N, ncol = 2,F)
ac = cov(matkov[,1], matkov[,2], use = "everything", method = "pearson")
cov[k,1] <- k
cov[k,2] <- ac
}
cov
## k acf
## 1 1 3317.373
## 2 2 3283.307
## 3 3 3250.237
## 4 4 3217.642
## 5 5 3185.131
## 6 6 3152.848
## 7 7 3121.146
## 8 8 3088.870
## 9 9 3056.902
## 10 10 3024.615
## 11 11 2992.554
## 12 12 2960.778
## 13 13 2929.560
## 14 14 2898.557
## 15 15 2866.976
## 16 16 2836.329
#Menghitung korelasi
kor<-data.frame(k=c(1:maxLag),
korelasi=NA)
for (k in 1:maxLag) {
kor[k,2]=cov[k,2]/var(data$Yt_p)
}
kor
## k korelasi
## 1 1 0.9899074
## 2 2 0.9797419
## 3 3 0.9698740
## 4 4 0.9601476
## 5 5 0.9504460
## 6 6 0.9408128
## 7 7 0.9313529
## 8 8 0.9217219
## 9 9 0.9121825
## 10 10 0.9025480
## 11 11 0.8929809
## 12 12 0.8834991
## 13 13 0.8741836
## 14 14 0.8649322
## 15 15 0.8555085
## 16 16 0.8463634
ts.sim.ar<-arima.sim(list(order=c(1,0,0),ar
=0.3),n=100)
ts.plot(ts.sim.ar, main = 'Simulasi AR(1)')
acf(ts.sim.ar,main='ACF Simulasi AR 1')
pacf(ts.sim.ar,main='PACF Simulasi AR 1')
ts.sim.ma<-arima.sim(list(order=c(0,0,1),ma
=0.3),n=100)
ts.plot(ts.sim.ma, main = 'Simulasi Plot MA(1)')
acf(ts.sim.ma,main='')
pacf(ts.sim.ma,main='')
ts.sim.ar1<-arima.sim(list(order=c(2,0,3),
ar=c(0.3,0.6),
ma=c(0.2,0.3,0.4)),n=100)
ts.plot(ts.sim.ar1)
acf(ts.sim.ar1,main='')
pacf(ts.sim.ar1,main='')
Tn <- 45 #banyak datum
phi <- 0.5 #parameter AR
theta <-0.5 #parameter MA
sig <- 7.2 #variansi error
Y <-c(rep(0,Tn)) #bangun vector Y
e <- rnorm(Tn,0,sig) #bangkitkan data error
#AR(1)
Y[2:Tn] <- (phi*Y[1:(Tn-1)])+e[2:Tn]
#MA(1)
Y[2:Tn] <- e[2:Tn]-(theta*e[1:(Tn-1)])
#ARMA(1,1)
Y[2:Tn] <- (phi*Y[1:(Tn-1)])+ e[2:Tn]-(theta*e[1:(Tn-1)])
plot(Y)
# Mensimulasikan ACF
#ACF Manual
maxLag<-12
df<-data.frame(1:Tn,Y)
acfResults<- data.frame(k=c(rep(1:maxLag)),
acf=NA,
stringsAsFactors = FALSE)
for (k in 1: maxLag) {
N <-nrow(df)-k
acf<- sum((Y[1:N]-mean(Y))%*%(Y[(1+k):(N+k)]-mean(Y)))/sum((Y-mean(Y))^2)
acfResults[k,1] <-k
acfResults[k,2] <-acf
}
plot(acfResults[,2],
type="h",
ylim=c(-.4,.3),
xlab="lag",ylab="ACF"
,main='Plot ACF Manual')
abline(h=0)
abline(h=c(1,-1)*1.96/sqrt(Tn),
lty=2,col="blue")
#ACF Automatic
acf(Y)
acf(Y, plot=F)
##
## Autocorrelations of series 'Y', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## -0.098 -0.121 0.210 0.053 0.187 -0.054 -0.095 0.112 0.000 -0.108 -0.065
## 12 13 14 15 16
## -0.048 -0.096 -0.124 0.028 -0.090
#Bandingkan Hasilnya
Hasil<-data.frame(k=c(rep(1:maxLag)),acf_manual=NA, acf_R=NA, stringsAsFactors = FALSE)
Hasil[,1] <-c(rep(1:maxLag))
Hasil[,2] <-acf(Y,plot=F)
Hasil[,3] <-acfResults[,2]
#Perhitungan PACF Manual
maxLag <- 12
df <- data.frame(1:Tn,Y)
PacfResults <- data.frame(k=c(1:maxLag),pacf=NA,stringsAsFactors = FALSE)
rho = acfResults[,2]
phi = matrix(0, nrow=maxLag,ncol=maxLag)
phi[1,1]=rho[1]
for (i in 2:maxLag){
phi[i,i]=(rho[i]-(sum(phi[i-1,1:i-1]*rho[(i-1):1])))/(1-sum(phi[i-1,1:i-1]*rho[1:i-1]))
for (j in 1:maxLag){
if (j<i){
phi[i,j]=phi[i-1,j]-(phi[i,i]*phi[i-1,i-j])
}
else {
0
}
}
PacfResults[,2] <-diag(phi)
}
plot(PacfResults[,2],
type="h",
ylim=c(-.3,.3),
xlab="lag",
ylab="PACF",
main='Plot PACF Manual')
abline(h=0)
abline(h=c(1,-1)*1.96/sqrt(length(phi)),
lty=2,col="blue")
#Perhitungan PACF dengan R
pacf(Y)
pacf(df, plot=FALSE)
##
## Partial autocorrelations of series 'df', by lag
##
## , , X1.Tn
##
## X1.Tn Y
## 0.935 ( 1) 0.052 ( -1)
## -0.050 ( 2) 0.335 ( -2)
## -0.015 ( 3) 0.171 ( -3)
## -0.048 ( 4) 0.241 ( -4)
## -0.028 ( 5) 0.159 ( -5)
## -0.023 ( 6) 0.050 ( -6)
## -0.020 ( 7) -0.097 ( -7)
## -0.050 ( 8) 0.147 ( -8)
## -0.005 ( 9) 0.242 ( -9)
## -0.108 ( 10) -0.137 (-10)
## -0.043 ( 11) -0.081 (-11)
## 0.032 ( 12) 0.177 (-12)
## -0.005 ( 13) -0.078 (-13)
##
## , , Y
##
## X1.Tn Y
## -0.033 ( 1) -0.104 ( 1)
## 0.107 ( 2) -0.152 ( 2)
## -0.059 ( 3) 0.204 ( 3)
## -0.083 ( 4) 0.056 ( 4)
## -0.015 ( 5) 0.259 ( 5)
## 0.096 ( 6) -0.053 ( 6)
## -0.033 ( 7) -0.111 ( 7)
## 0.122 ( 8) -0.077 ( 8)
## -0.121 ( 9) 0.031 ( 9)
## -0.031 ( 10) -0.114 ( 10)
## -0.098 ( 11) -0.155 ( 11)
## -0.019 ( 12) -0.200 ( 12)
## 0.040 ( 13) -0.161 ( 13)
#Bandingkan Hasilnya
Hasil<-data.frame(k=c(rep(1:maxLag)),
pacf_manual=NA,
pacf_R=NA,
stringsAsFactors = FALSE)
Hasil[,1] <-c(rep(1:maxLag))
Hasil[,2] <-pacf(Y,plot=F)
Hasil[,3] <-PacfResults[,2]