## Do not edit this file manually.
## It has been automatically generated from *.org sources.
makeUpdateFun_arma0 <- function(pind, qind, s2ind){
pind
qind
s2ind
## todo: check the case when both p=0 and q=0?
## this is based on an example in package KFAS
function(pars, model, estimate = TRUE){
wrk <- try(KFAS::SSMarima(KFAS::artransform(pars[pind]), KFAS::artransform(pars[qind]),
Q = exp(pars[s2ind])), silent = TRUE)
if(inherits(wrk, "try-error")){
if(estimate) 1e100 else model
}else{
model["T", "arima"] <- wrk$T
model["R", "arima"] <- wrk$R
model["P1", "arima"] <- wrk$P1
model["Q", "arima"] <- wrk$Q
if(estimate){
-logLik(model)
}else
model
}
}
}
makeArma0SSModel <- function(ar = numeric(0), ma = numeric(0), Q = 1, H = 0, y){
## todo: more care with the environment?
## 2022-01-19 moved KFAS to Suggests.
##
## however, then SSMarima is not visible and in a model formula can't
## be qualified with 'KFAS::'. So, copying it here.
SSMarima <- KFAS::SSMarima
arma0 <- KFAS::SSModel(y ~ -1 + SSMarima(ar = ar, ma = ma, Q = 1), H = 0)
p <- length(ar)
q <- length(ma)
pind <- seq_len(p)
qind <- p + seq_len(q)
s2ind <- p + q + 1
## todo: check the case when both p=0 and q=0?
likfn0 <- makeUpdateFun_arma0(pind, qind, s2ind)
list(model = arma0, updateFun = likfn0)
}
fitArma0Model <- function(arma0ss, init = NULL, method = "BFGS", hessian = TRUE, ...){
e <- environment(arma0ss$updateFun)
p <- length(e$pind)
q <- length(e$qind)
if(is.null(init)){
init <- numeric(p + q + 1)
}
arma0_fit <- optim(par = init,
fn = arma0ss$updateFun,
method = method,
model = arma0ss$model,
hessian = hessian, ...)
## turn into SS model object
modelObj <- arma0ss$updateFun(arma0_fit$par, arma0ss$model, FALSE)
ar <- modelObj$T[seq_len(p)]
ma <- modelObj$R[seq_len(q+1)[-1]] # 2:(q+1), safely
sigma2 <- as.vector(modelObj$Q)
# nn <- dim(arma0_fit$hessian)[1]
# sqrt(diag(solve(arma0_fit$hessian)))
# ## ignoring sigma2:
# sqrt(diag(solve(arma0_fit$hessian[-nn, -nn])))
fi <- solve(arma0_fit$hessian[-(p+q+1), -(p+q+1)])
## correct for phi = phi(parcor)
Jphi <- pacf2ArWithJacobian(ar2Pacf(ar))$J
Jtheta <- pacf2ArWithJacobian(ar2Pacf(ma), TRUE)$J
J2 <- dbind(Jphi, Jtheta)
## correct for the tanh transform.
# Jtanh <- diag(1/cosh(arma0_fit$par[seq_len(p+q)])^2)
# fi <- J2 %*% Jtanh %*% fi %*% Jtanh %*% t(J2)
# J <- J2 %*% Jtanh
J <- J2 / rep(cosh(arma0_fit$par[seq_len(p+q)])^2, each = p + q)
fi <- J %*% tcrossprod(fi, J) # equivalent to J %*% fi %*% t(J)
#
# fi <- solve(arma0_fit$hessian)
# J <- dbind(J2 %*% Jtanh, exp(arma0_fit$par[p+q+1]))
# fi <- J %*% fi %*% t(J)
# fi <- fi[-(p+q+1), -(p+q+1)]
se <- sqrt(diag(fi))
se.asy <- sqrt(diag(## 2022-02-14 was: InformationMatrixARMA(ar, -ma) /
## .FisherInfo has the SP convention:
.FisherInfo(-ar, ma)) /
attr(modelObj, "n", exact = TRUE))
param <- c(ar, ma)
nams <- c(paste0("ar", seq_len(p)), paste0("ma", seq_len(q)))
names(se.asy) <- names(se) <- names(param) <- nams
list(par = param, se = se, se.asy = se.asy,
optim = arma0_fit, SSModel = modelObj, vcov = fi)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.