Implementing the \(\text{ARIMA}(1, 1, 0)\) model using rgeneric
Jeff Eaton (jeffrey.eaton@imperial.ac.uk)
Abstract
To write!
library(INLA)## Loading required package: Matrix## Loading required package: foreach## Loading required package: parallel## Loading required package: sp## This is INLA_22.05.07 built 2022-05-07 09:58:26 UTC.
## - See www.r-inla.org/contact-us for how to get help.
## - To enable PARDISO sparse library; see inla.pardiso()library(Matrix)
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forcats)1 R-INLA \(\text{AR}(1)\) example using rgeneric
From https://inla.r-inla-download.org/r-inla.org/doc/vignettes/rgeneric.pdf.
inla.rgeneric.ar1.model <- function(
cmd = c("graph", "Q", "mu", "initial", "log.norm.const", "log.prior", "quit"),
theta = NULL) {
## for reference and potential storage for objects to
## cache, this is the environment of this function
## which holds arguments passed as `...` in
## `inla.rgeneric.define()`.
interpret.theta = function() {
return(list(prec = exp(theta[1L]),
rho = 2 * exp(theta[2L])/(1 + exp(theta[2L])) - 1))
}
envir = parent.env(environment())
graph = function(){ return (Q()) }
Q = function() {
p = interpret.theta()
i = c(1L, n, 2L:(n - 1L), 1L:(n - 1L))
j = c(1L, n, 2L:(n - 1L), 2L:n)
x = p$prec/(1 - p$rho^2) *
c(1L, 1L, rep(1 + p$rho^2, n - 2L),
rep(-p$rho, n - 1L))
return (sparseMatrix(i = i, j = j, x = x, giveCsparse = FALSE))
}
mu = function() { return(numeric(0)) }
log.norm.const = function() {
p = interpret.theta()
prec.i = p$prec / (1.0 - p$rho^2)
val = n * (- 0.5 * log(2*pi) + 0.5 * log(prec.i)) +
0.5 * log(1.0 - p$rho^2)
return (val)
}
log.prior = log.prior = function() {
p = interpret.theta()
val = dgamma(p$prec, shape = 1, rate = 1, log=TRUE) + theta[1L] +
dnorm(theta[2L], mean = 0, sd = 1, log=TRUE)
return (val)
}
initial = function() {
return (rep(1, 2))
}
quit = function() { return (invisible()) }
## sometimes this is useful, as argument 'graph' and 'quit'
## will pass theta=numeric(0) (or NULL in R-3.6...) as
## the values of theta are NOT
## required for defining the graph. however, this statement
## will ensure that theta is always defined.
if (!length(theta)) theta = initial()
val = do.call(match.arg(cmd), args = list())
return (val)
}
n <- 100
rho <- 0.9
x <- arima.sim(n, model = list(ar = rho)) * sqrt(1-rho^2)
y <- x + rnorm(n, sd = 0.1)
model <- inla.rgeneric.define(inla.rgeneric.ar1.model, n = n)## Warning in sparseMatrix(i = i, j = j, x = x, giveCsparse = FALSE): 'giveCsparse' has been deprecated; setting 'repr = "T"' for youformula <- y ~ -1 + f(idx, model = model)
r <- inla(formula, data = data.frame(y, idx = 1:n))## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'rgeneric' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.summary(r)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ",
## " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call,
## inla.arg = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug =
## debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 2.99, Running = 1.72, Post = 0.0214, Total = 4.73
## Random effects:
## Name Model
## idx RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for the Gaussian observations 1.83e+04 1.77e+04 1232.297 1.29e+04 65308.63 NA
## Theta1 for idx 6.26e-01 2.26e-01 0.151 6.38e-01 1.04 NA
## Theta2 for idx 1.92e+00 2.70e-01 1.431 1.91e+00 2.49 NA
##
## Marginal log-Likelihood: -74.64
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')fformula = y ~ -1 +
f(idx, model = "ar1",
hyper = list(prec = list(prior = "loggamma", param = c(1,1)),
rho = list(prior = "normal", param = c(0,1))))
rr = inla(fformula, data = data.frame(y, idx = 1:n))
plot(inla.smarginal(rr$marginals.hyperpar[[2]]),
type="l", lwd=5, col="red", xlab="stdev", ylab="density")
lines(inla.tmarginal(exp, r$internal.marginals.hyperpar[[2]]),
col="yellow")
plot(inla.smarginal(rr$marginals.hyperpar[[3]]),
type="l", lwd=5, col="red", xlab="rho", ylab="density")
lines(inla.tmarginal(function(x) 2*exp(x)/(1+exp(x))-1,
r$internal.marginals.hyperpar[[3]]),
col="yellow")
round(rbind(native = rr$cpu.used,
rgeneric = r$cpu.used), digits = 3)## Pre Running Post Total
## native 2.989 0.240 0.012 3.241
## rgeneric 2.992 1.716 0.021 4.7292 R-INLA \(\text{ARIMA}(1, 1, 0)\) using rgeneric
\(\alpha \sim \text{ARIMA}(1, 1, 0)\) is equivalent to \(\gamma ~ AR(1)\) where \(\gamma = D\alpha\) and \(D\) is the first order difference matrix. Let \(\tilde Q\) be the \(\text{AR}(1)\) precision matrix. Then \(\gamma' \tilde Q \gamma = \alpha' D' \tilde Q D \alpha\). Thus the precision matrix for \(\text{ARIMA}(1, 1, 0)\) is simply \(Q = D' \tilde Q D\). \(\det(Q) = n \det(\tilde Q)\). I don’t know why yet… The normalising constant has coefficient \(n - 1\).
inla.rgeneric.arima110.model <- function(cmd = c("graph", "Q", "mu", "initial",
"log.norm.const", "log.prior", "quit"),
theta = NULL) {
## for reference and potential storage for objects to
## cache, this is the environment of this function
## which holds arguments passed as `...` in
## `inla.rgeneric.define()`.
interpret.theta = function() {
return(list(prec = exp(theta[1L]),
rho = 2 * exp(theta[2L])/(1 + exp(theta[2L])) - 1))
}
envir = parent.env(environment())
graph = function(){ return (Q()) }
Q = function() {
n.ar1 <- n-1 # dimension of AR1 on first order differences
D1 = Matrix::diff(Matrix::Diagonal(n)) # first order difference operator
p = interpret.theta()
## Define both uperr and lower tri because D'*Q*D calculation
i = c(1L, n.ar1, 2L:(n.ar1 - 1L), 1L:(n.ar1 - 1L), 2L:n.ar1)
j = c(1L, n.ar1, 2L:(n.ar1 - 1L), 2L:n.ar1, 1L:(n.ar1 - 1L))
x = p$prec/(1 - p$rho^2) *
c(1L, 1L, rep(1 + p$rho^2, n.ar1 - 2L),
rep(-p$rho, n.ar1 - 1L),
rep(-p$rho, n.ar1 - 1L))
Q.ar1 <- sparseMatrix(i = i, j = j, x = x, giveCsparse = FALSE)
return (t(D1) %*% Q.ar1 %*% D1)
}
mu = function() { return(numeric(0)) }
log.norm.const = function() {
p = interpret.theta()
prec.i = p$prec / (1.0 - p$rho^2)
val = (n - 1) * (- 0.5 * log(2*pi) + 0.5 * log(prec.i)) +
0.5 * log(1.0 - p$rho^2) +
log(n) # + log(n) because det(Q) = n * det(Q.ar1)
return (val)
}
log.prior = log.prior = function() {
p = interpret.theta()
val = dgamma(p$prec, shape = 1, rate = 1, log=TRUE) + theta[1L] +
dnorm(theta[2L], mean = 0, sd = 1, log=TRUE)
return (val)
}
initial = function() {
return (rep(1, 2))
}
quit = function() { return (invisible()) }
## sometimes this is useful, as argument 'graph' and 'quit'
## will pass theta=numeric(0) (or NULL in R-3.6...) as
## the values of theta are NOT
## required for defining the graph. however, this statement
## will ensure that theta is always defined.
if (!length(theta)) theta = initial()
val = do.call(match.arg(cmd), args = list())
return (val)
}Simulate data
set.seed(1)
n <- 100
rho <- 0.8
c <- 0.3The parameter mean = <> to arima.sim() passes a mean value to rnorm() used to simulate the innovations.
To obtain a process with marginal slope c, pass the value mean = c * (1 - rho) / sqrt(1 - rho^2).
x <- arima.sim(n, model = list(order = c(1, 1, 0), ar = rho), mean = c * (1 - rho) / sqrt(1 - rho^2))
x <- x[-1] * sqrt(1 - rho^2)
y <- x + rnorm(n, sd = 0.1)Define R-INLA \(\text{ARIMA}(1,1,0)\) model object.
arima110_model <- inla.rgeneric.define(inla.rgeneric.arima110.model, n = n)## Warning in sparseMatrix(i = i, j = j, x = x, giveCsparse = FALSE): 'giveCsparse' has been deprecated; setting 'repr = "T"' for youFit R-INLA model to full data
formula <- y ~ idx_lin + f(idx, model = arima110_model, constr = TRUE)
fit <- inla(formula, data = data.frame(y, idx_lin = 1:n, idx = 1:n),
control.predictor = list(compute = TRUE))
summary(fit)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ",
## " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call,
## inla.arg = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug =
## debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 2.89, Running = 1.96, Post = 0.0153, Total = 4.87
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -1.320 7.479 -16.137 -1.322 13.509 NA 0
## idx_lin 0.609 0.148 0.316 0.609 0.903 NA 0
##
## Random effects:
## Name Model
## idx RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for the Gaussian observations 1.61e+04 1.67e+04 1081.609 1.09e+04 6.07e+04 NA
## Theta1 for idx 5.13e-01 2.08e-01 0.079 5.22e-01 8.95e-01 NA
## Theta2 for idx 1.57e+00 2.44e-01 1.120 1.56e+00 2.08e+00 NA
##
## Marginal log-Likelihood: -93.67
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')fit$summary.fitted.values %>%
mutate(idx = 1:n) %>%
ggplot(aes(idx, mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
geom_ribbon(color = NA, alpha = 0.3) +
geom_line() +
geom_point(aes(y = y), color = "red3")
Truncate last 25% projection data
y_trunc <- y
y_trunc[round(n*0.75):n] <- NA
fit_trunc <- inla(y_trunc ~ idx_lin + f(idx, model = arima110_model, constr = TRUE),
data = data.frame(y_trunc, idx_lin = 1:n, idx = 1:n),
control.predictor = list(compute = TRUE))
summary(fit_trunc)##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, ",
## " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, control.compute = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla, control.fixed = control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call,
## inla.arg = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, safe = FALSE, debug =
## debug, ", " .parent.frame = .parent.frame)")
## Time used:
## Pre = 2.72, Running = 1.9, Post = 0.0153, Total = 4.63
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -1.104 8.412 -17.729 -1.139 15.722 NA 0
## idx_lin 0.600 0.177 0.244 0.601 0.948 NA 0
##
## Random effects:
## Name Model
## idx RGeneric2
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Precision for the Gaussian observations 2.07e+04 2.01e+04 1391.694 1.47e+04 73703.55 NA
## Theta1 for idx 4.76e-01 2.38e-01 -0.025 4.89e-01 0.91 NA
## Theta2 for idx 1.62e+00 2.84e-01 1.104 1.60e+00 2.22 NA
##
## Marginal log-Likelihood: -70.94
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')fit_trunc$summary.fitted.values %>%
mutate(idx = 1:n) %>%
ggplot(aes(idx, mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
geom_ribbon(color = NA, alpha = 0.3) +
geom_line() +
geom_point(aes(y = y), color = "red3")
2.1 Compare ARIMA(1,1,0) + c with other models
- \(\text{ARIMA}(1, 1, 0)\) (no constant)
- \(\text{RW}(2)\)
- \(\text{AR}(2)\)
- \(\text{AR}(2) + c\)
fit_arima110 <- inla(y_trunc ~ f(idx, model = arima110_model, constr = TRUE),
data = data.frame(y_trunc, idx_lin = 1:n, idx = 1:n),
control.predictor = list(compute = TRUE),
control.compute = list(config = TRUE))
fit_rw2 <- inla(y_trunc ~ f(idx, model = "rw2"),
data = data.frame(y_trunc, idx_lin = 1:n, idx = 1:n),
control.predictor = list(compute = TRUE))
fit_ar2 <- inla(y_trunc ~ f(idx, model = "ar", order = 2),
data = data.frame(y_trunc, idx_lin = 1:n, idx = 1:n),
control.predictor = list(compute = TRUE),
control.compute = list(config = TRUE))
fit_ar2_c <- inla(y ~ idx_lin + f(idx, model = "ar", order = 2),
data = data.frame(y_trunc, idx_lin = 1:n, idx = 1:n),
control.predictor = list(compute = TRUE))
fitted_all <- bind_rows(
mutate(fit_trunc$summary.fitted.values, idx = row_number(), model = "ARIMA(1,1,0) + c"),
mutate(fit_arima110$summary.fitted.values, idx = row_number(), model = "ARIMA(1,1,0)"),
mutate(fit_rw2$summary.fitted.values, idx = row_number(), model = "RW2"),
mutate(fit_ar2$summary.fitted.values, idx = row_number(), model = "AR(2)"),
mutate(fit_ar2_c$summary.fitted.values, idx = row_number(), model = "AR(2) + c")
) %>%
mutate(model = fct_inorder(model))
fitted_all %>%
ggplot(aes(idx, mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
geom_ribbon(aes(fill = model, color = model), alpha = 0.2, linetype = "dotted", size = 0.25) +
geom_line(aes(color = model)) +
geom_point(aes(x = idx, y = y), color = "red3", inherit.aes = FALSE,
data = data.frame(idx = seq_along(y), y)) +
scale_color_brewer(palette = "Set1") +
scale_fill_brewer(palette = "Set1")
3 Original computing environment
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Monterey 12.6
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.5.1 dplyr_1.0.9 ggplot2_3.3.6 INLA_22.05.07 sp_1.4-7 foreach_1.5.2 Matrix_1.4-1 knitr_1.39
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.2 xfun_0.31 TMB_1.9.1 bslib_0.3.1 purrr_0.3.4 splines_4.2.0 lattice_0.20-45 colorspace_2.0-3 vctrs_0.4.1 generics_0.1.2 htmltools_0.5.2 yaml_2.3.5
## [13] mgcv_1.8-40 utf8_1.2.2 rlang_1.0.2 jquerylib_0.1.4 pillar_1.7.0 glue_1.6.2 withr_2.5.0 DBI_1.1.2 RColorBrewer_1.1-3 lifecycle_1.0.1 stringr_1.4.0 MatrixModels_0.5-0
## [25] munsell_0.5.0 gtable_0.3.0 codetools_0.2-18 evaluate_0.15 labeling_0.4.2 fastmap_1.1.0 fansi_1.0.3 highr_0.9 scales_1.2.0 jsonlite_1.8.0 farver_2.1.0 Deriv_4.1.3
## [37] digest_0.6.29 stringi_1.7.6 bookdown_0.26 grid_4.2.0 cli_3.3.0 tools_4.2.0 magrittr_2.0.3 sass_0.4.1 tibble_3.1.7 crayon_1.5.1 pkgconfig_2.0.3 ellipsis_0.3.2
## [49] assertthat_0.2.1 rmarkdown_2.14 iterators_1.0.14 R6_2.5.1 nlme_3.1-157 compiler_4.2.0