Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment