KFAS_ARMA.R

require(KFAS) # 1.2.9
require(dplyr)
require(tidyr)
require(ggplot2)

# --- ARIMA(2, 1) with linear trend ---
# generate a dataset
set.seed(42)
t <- 100
y <- arima.sim(n = t, model = list(ar=c(.3, -0.1),  ma=.2), sd=.1) + seq(from=1, to=10, length.out = t)
plot(y)

# create model
model <- SSModel(y~ SSMtrend(2, Q=list(0, 0)) + SSMarima(ar=c(0, 0), ma=0, Q=NA), H=0)

# estimate unknown parameters
update.fn <- function(pars, model) {
  model <- SSModel(y ~ SSMtrend(2, Q=list(0, 0)) +
                     SSMarima(ar =artransform(pars[1:2]), ma=pars[3], Q = exp(pars[4])), H = 0)
  model$T["slope", "slope", 1] <- pars[5]
  return(model)
}


fit.model <- fitSSM(model, inits=rep(.1, 5), updatefn = update.fn, method="SANN")

# check the result
print(paste("Converged: ", fit.model$optim.out$convergence == 0))
print(paste("log likelihood:", round(-fit.model$optim.out$value, 3)))
fit.model$model$Z
fit.model$model$T
fit.model$model$R
fit.model$model$Q
fit.model$model$H

filter.model <- KFS(fit.model$model, smoothing = "mean")

# plot prediced values
pred <- predict(fit.model$model, n.ahead = 20, interval = "confidence", type="response", level=.95)
df <- data.frame(t=1:length(y), y=as.numeric(y)) 
df <- bind_rows(df, data.frame(pred) %>% rename(y=fit) %>% mutate(t=time(pred))) # ignore warning
ggplot(filter(df, t>= 10)) + geom_line(aes(x=t, y=y), lwd=1) + 
  geom_ribbon(aes(x=t, ymin=lwr, ymax=upr), fill="blue", alpha=.2)