SIMULASI AR(1)

Simulating an AR(p) process

set.seed(456)
## list description for AR(1) model with small coef
AR.sm <- list(order = c(1, 0, 0), ar = 0.1, sd = 0.1)
## list description for AR(1) model with large coef
AR.lg <- list(order = c(1, 0, 0), ar = 0.9, sd = 0.1)
## simulate AR(1)
AR1.sm <- arima.sim(n = 50, model = AR.sm)
AR1.lg <- arima.sim(n = 50, model = AR.lg)

Now let’s plot the 2 simulated series.

## setup plot region
par(mfrow = c(1, 2))
## get y-limits for common plots
ylm <- c(min(AR1.sm, AR1.lg), max(AR1.sm, AR1.lg))
## plot the ts
plot.ts(AR1.sm, ylim = ylm, ylab = expression(italic(x)[italic(t)]), 
    main = expression(paste(phi, " = 0.1")))
plot.ts(AR1.lg, ylim = ylm, ylab = expression(italic(x)[italic(t)]), 
    main = expression(paste(phi, " = 0.9")))

set.seed(123)
## list description for AR(1) model with small coef
AR.pos <- list(order = c(1, 0, 0), ar = 0.5, sd = 0.1)
## list description for AR(1) model with large coef
AR.neg <- list(order = c(1, 0, 0), ar = -0.5, sd = 0.1)
## simulate AR(1)
AR1.pos <- arima.sim(n = 50, model = AR.pos)
AR1.neg <- arima.sim(n = 50, model = AR.neg)

OK, let’s plot the 2 simulated series.

## setup plot region
par(mfrow = c(1, 2))
## get y-limits for common plots
ylm <- c(min(AR1.pos, AR1.neg), max(AR1.pos, AR1.neg))
## plot the ts
plot.ts(AR1.pos, ylim = ylm, ylab = expression(italic(x)[italic(t)]), 
    main = expression(paste(phi[1], " = 0.5")))
plot.ts(AR1.neg, ylab = expression(italic(x)[italic(t)]), main = expression(paste(phi[1], 
    " = -0.5")))

Correlation structure of AR(p) processes

set.seed(123)
## the 4 AR coefficients
ARp <- c(0.7, 0.2, -0.1, -0.3)
## empty list for storing models
AR.mods <- list()
## loop over orders of p
for (p in 1:4) {
    ## assume SD=1, so not specified
    AR.mods[[p]] <- arima.sim(n = 10000, list(ar = ARp[1:p]))
}

Now that we have our four AR( p ) models, lets look at plots of the time series, ACF’s, and PACF’s.

## setup plot region
par(mfrow = c(4, 3))
## loop over orders of p
for (p in 1:4) {
    plot.ts(AR.mods[[p]][1:50], ylab = paste("AR(", p, ")", sep = ""))
    acf(AR.mods[[p]], lag.max = 12)
    pacf(AR.mods[[p]], lag.max = 12, ylab = "PACF")
}

SIMULASI MA(1)

Simulating an MA(q) process

set.seed(123)
## list description for MA(1) model with small coef
MA.sm <- list(order = c(0, 0, 1), ma = 0.2, sd = 0.1)
## list description for MA(1) model with large coef
MA.lg <- list(order = c(0, 0, 1), ma = 0.8, sd = 0.1)
## list description for MA(1) model with large coef
MA.neg <- list(order = c(0, 0, 1), ma = -0.5, sd = 0.1)
## simulate MA(1)
MA1.sm <- arima.sim(n = 50, model = MA.sm)
MA1.lg <- arima.sim(n = 50, model = MA.lg)
MA1.neg <- arima.sim(n = 50, model = MA.neg)

with their associated plots.

## setup plot region
par(mfrow = c(1, 3))
## plot the ts
plot.ts(MA1.sm, ylab = expression(italic(x)[italic(t)]), main = expression(paste(theta, 
    " = 0.2")))
plot.ts(MA1.lg, ylab = expression(italic(x)[italic(t)]), main = expression(paste(theta, 
    " = 0.8")))
plot.ts(MA1.neg, ylab = expression(italic(x)[italic(t)]), main = expression(paste(theta, 
    " = -0.5")))

#Correlation structure of MA(q) processes

set.seed(123)
## the 4 MA coefficients
MAq <- c(0.7, 0.2, -0.1, -0.3)
## empty list for storing models
MA.mods <- list()
## loop over orders of q
for (q in 1:4) {
    ## assume SD=1, so not specified
    MA.mods[[q]] <- arima.sim(n = 1000, list(ma = MAq[1:q]))
}

Now that we have our four MA(q) models, lets look at plots of the time series, ACF’s, and PACF’s.

## set up plot region
par(mfrow = c(4, 3))
## loop over orders of q
for (q in 1:4) {
    plot.ts(MA.mods[[q]][1:50], ylab = paste("MA(", q, ")", sep = ""))
    acf(MA.mods[[q]], lag.max = 12)
    pacf(MA.mods[[q]], lag.max = 12, ylab = "PACF")
}

SIMULASI ARMA(1,1)

Fitting ARMA(p,q) models with arima()

set.seed(123)
## ARMA(2,2) description for arim.sim()
ARMA22 <- list(order = c(2, 0, 2), ar = c(-0.7, 0.2), ma = c(0.7, 
    0.2))
## mean of process
mu <- 5
## simulated process (+ mean)
ARMA.sim <- arima.sim(n = 10000, model = ARMA22) + mu
## estimate parameters
arima(x = ARMA.sim, order = c(2, 0, 2))
## 
## Call:
## arima(x = ARMA.sim, order = c(2, 0, 2))
## 
## Coefficients:
##           ar1     ar2     ma1     ma2  intercept
##       -0.7079  0.1924  0.6912  0.2001     4.9975
## s.e.   0.0291  0.0284  0.0289  0.0236     0.0125
## 
## sigma^2 estimated as 0.9972:  log likelihood = -14175.92,  aic = 28363.84

Searching over model orders

## empty list to store model fits
ARMA.res <- list()
## set counter
cc <- 1
## loop over AR
for (p in 0:3) {
    ## loop over MA
    for (q in 0:3) {
        ARMA.res[[cc]] <- arima(x = ARMA.sim, order = c(p, 0, 
            q))
        cc <- cc + 1
    }
}
## Warning in arima(x = ARMA.sim, order = c(p, 0, q)): possible convergence
## problem: optim gave code = 1
## get AIC values for model evaluation
ARMA.AIC <- sapply(ARMA.res, function(x) x$aic)
## model with lowest AIC is the best
ARMA.res[[which(ARMA.AIC == min(ARMA.AIC))]]
## 
## Call:
## arima(x = ARMA.sim, order = c(p, 0, q))
## 
## Coefficients:
##           ar1     ar2     ma1     ma2  intercept
##       -0.7079  0.1924  0.6912  0.2001     4.9975
## s.e.   0.0291  0.0284  0.0289  0.0236     0.0125
## 
## sigma^2 estimated as 0.9972:  log likelihood = -14175.92,  aic = 28363.84