R/TwoSimpleUnitRoots.R

In YueyunZhu/UnitRootsofPIAR: Eliminate unit roots for periodically integrated autoregressive models

#' @title Function to compute the coefficients of periodic integrated filter
#' @description Function to compute the coefficients of periodic integrated filter.
#'              Assume there are two simple unit roots in the process.
#' @param Xsim similarity matrix for Jordan decomposition of the multi-companion matrix F
#' @param d period
#' @param eigval eigen value of F
#'
#' @return  coefficients of $(1-\alpha_{1s}L-\alpha_{2s}L^2)$
#' @export
#'
#' @examples
#' library(mcompanion)
#' Fd <- sim_mc(dim = 4, mo = 4, mo.col = 2, eigval = c(1, 1), len.block = c(1, 1))
#' icoef <- piar2unit.simple(Xsim = Fd$eigvec, d = 4, eigval = c(1, 1))
piar2unit.simple <- function(Xsim, d, eigval) {
  vecX1 <- Xsim[, 1]
  vecX2 <- Xsim[, 2]
  alpha1 <- rep(NA, d)
  alpha2 <- rep(NA, d)
  Delta <- matrix(NA, d, d)
  for (i in 1:d) {
    {
      for (j in 1:d) {
        Delta[i, j] <- vecX1[i] * vecX2[j] - vecX1[j] * vecX2[i]
      }
    }
  }
  alpha2[1] <- (eigval[1] * vecX1[d] * vecX2[1] - eigval[2] * vecX1[1] * vecX2[d])/(Delta[2, 1])
  alpha1[1] <- (eigval[1] * vecX1[d] - alpha2[1] * vecX1[2])/vecX1[1]
  alpha2[2] <- (eigval[1] * eigval[2] * (Delta[d - 1, d]))/
    (eigval[2] * vecX1[1] * vecX2[d] - eigval[1] * vecX1[d] * vecX2[1])
  alpha1[2] <- (eigval[1] * vecX1[d - 1] - alpha2[2] * vecX1[1])/(eigval[1] * vecX1[d])
  if (d > 2) {
    for (i in 3:d) {
      alpha1[i] <- Delta[d - i + 1, d - i + 3]/Delta[d - i + 2, d - i + 3]
      alpha2[i] <- Delta[d - i + 2, d - i + 1]/Delta[d - i + 2, d - i + 3]
    }
    return(par = list(alpha1=alpha1, alpha2=alpha2))
  } else {
    return(par = list(alpha1=alpha1, alpha2=alpha2))
  }
}