R/wrapKFAS.R

In sarima: Simulation and Prediction with Seasonal ARIMA Models

## 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)
}