library(ggplot2)
library(dplyr)
library(dygraphs)
library(xts)
library(plotly)
library(TSA)
library(hrbrthemes)

Random Walk

# 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 Manual

#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

Simulasi Model AR dan MA

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='')

Mesimulasikan Model dengan Parameter

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]

Mensimulasikan PACF

#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]