R/StartingValues.R

In LeopoldoCatania/GAS: Generalized Autoregressive Score Models

StartingNu <- function(vY) {
  dStart = c(unmapVec_C(12.0, LowerNu(), UpperNu()))
  dNu = try(exp(solnp(dStart, function(x, vY) {
    dNu =  c(Map_Vec(x, LowerNu(), UpperNu()))
    - sum(dt(vY, dNu, log = TRUE))
  }, vY = vY, control = list(trace = 0))$pars) + LowerNu(), silent = TRUE)

  if (is(dNu, "try-error"))
    dNu = 7.0

  if (dNu > UpperNu())
    dNu = UpperNu() - 1.0
  if (dNu <= LowerNu())
    dNu = LowerNu() + 0.02

  return(dNu)
}

StaticStarting_Uni <- function(vY, Dist, iK) {

  if (Dist == "std") {

    dMu = mean(vY)
    dNu = 8.0
    dPhi2 = var(vY) * (dNu - 2.0)/dNu

    vTheta = c(dMu, dPhi2, dNu)

  }

  if (Dist == "sstd") {

    dMu = mean(vY)
    dNu = 8.0
    dSigma = sd(vY)
    dXi = 1.0

    vTheta = c(dMu, dSigma, dXi, dNu)

  }

  if (Dist == "norm") {

    dMu = mean(vY)
    dSigma2 = var(vY)

    vTheta = c(dMu, dSigma2)

  }
  if (Dist == "snorm") {

    dMu = mean(vY)
    dSigma = sd(vY)
    dXi = 1.0

    vTheta = c(dMu, dSigma, dXi)

  }
  if (Dist == "ast" | Dist == "ast1") {

    dNu1 = StartingNu(vY)

    if (Dist == "ast") {

      dNu2 = dNu1 * 1.5

      if (dNu2 > UpperNu()) {
        dNu2 = UpperNu() - 1.0
      }

    } else {

      dNu2 = dNu1

    }

    dAlpha = 0.5
    dMu = mean(vY)
    dSigma = sd(vY)

    if (Dist == "ast") {

      vTheta = c(dMu, dSigma, dAlpha, dNu1, dNu2)

    } else {

      vTheta = c(dMu, dSigma, dAlpha, dNu1)

    }
  }

  if (Dist == "ghskt") {

    dMu = mean(vY)
    dNu = 10.0
    dSigma = sd(vY)
    dBetaBar = 1e-4

    vTheta = c(dMu, dSigma, dBetaBar, dNu)

  }

  if (Dist == "poi") {

    vTheta = mean(vY)

  }
  if (Dist == "ber") {

    dPi = sum(vY)/length(vY)

    vTheta = c(dPi)

  }
  if (Dist == "gamma") {

    dMean = mean(vY)
    dSigma2 = var(vY)

    dBeta = dMean/dSigma2
    dAlpha = dMean^2/dSigma2

    vTheta = c(dAlpha, dBeta)

  }

  if (Dist == "exp") {

    vTheta = c(1.0/mean(vY))

  }
  if (Dist == "beta") {

    dMean = mean(vY)
    dSigma2 = var(vY)

    dAlpha = dMean * (dMean * (1.0 - dMean)/dSigma2 - 1.0)
    dBeta = (1.0 - dMean) * (dMean * (1.0 - dMean)/dSigma2 - 1.0)

    vTheta = c(dAlpha, dBeta)
  }
  if (Dist == "ald") {

    dTheta = mean(vY)
    dSigma = sd(vY)
    dKappa = 1.0

    vTheta = c(dTheta, dSigma, dKappa)

  }

  if (Dist == "negbin") {

    dMu = mean(vY)
    dVar  = var(vY)

    if (dVar > dMu) {

      dPi = 1.0 / (1.0 + (dVar - dMu)/dMu)
      dNu = dMu^2/(dVar - dMu)

    } else {
      dPi = 0.5
      dNu = 1.0
    }

    vTheta = c(dPi, dNu)

  }

  if (Dist == "skellam") {

    dMu = mean(vY)
    dSigma2 = var(vY)

    vTheta = c(dMu, dSigma2)

  }

  vTheta_tilde = as.numeric(UnmapParameters_univ(vTheta, Dist, iK = iK))

  names(vTheta_tilde) = FullNamesUni(Dist)

  return(vTheta_tilde)
}

StaticStarting_Multi <- function(mY, Dist, iN) {

  vEmpRho = build_vR(cor(t(mY)), iN)
  vEmpPhi = UnMapR_C(vEmpRho, iN)

  vEmpMu = apply(mY, 1L, mean)
  vEmpSigma = apply(mY, 1L, sd)

  if (Dist == "mvnorm") {
    vMu_tilde = vEmpMu
    vSigma_tilde = log(vEmpSigma)
    vPw = c(vMu_tilde, vSigma_tilde, vEmpPhi)
  }
  if (Dist == "mvt") {

    dNu = StartingNu(c(mY))

    vMu_tilde = vEmpMu
    vSigma_tilde = log(vEmpSigma * (dNu - 2.0)/dNu)
    vPw = c(vMu_tilde, vSigma_tilde, vEmpPhi, log(dNu - LowerNu()))
  }

  return(as.numeric(vPw))

}

UniGAS_Starting <- function(vY, iT, iK, Dist, ScalingType, GASPar, fn.optimizer) {

  StaticFit = StaticMLFIT(vY, Dist, fn.optimizer)
  vUncValues = StaticFit$optimiser$pars

  names(vUncValues) = paste("kappa", 1:iK, sep = "")

  if (iK > 1L) {

    vA = starting_vA_Uni(vY, vUncValues, mB = diag(rep(0.9, iK)), dA_foo = 1e-06, iT, iK, Dist,
                         ScalingType = ScalingType, GASPar)
    vB = starting_vB_Uni(vY, vUncValues, dB_foo = 0.9, mA = diag(vA), iT, iK, Dist, ScalingType = ScalingType,
                         GASPar)
    vKappa = (diag(iK) - diag(vB)) %*% vUncValues
    names(vKappa) = paste("kappa", 1:iK, sep = "")

  } else {

    vA = starting_vA_Uni(vY, vUncValues, mB = matrix(0.9, iK, iK), dA_foo = 1e-06, iT, iK, Dist,
                         ScalingType = ScalingType, GASPar)
    vB = starting_vB_Uni(vY, vUncValues, dB_foo = 0.9, mA = matrix(vA, iK, iK), iT, iK, Dist, ScalingType = ScalingType,
                         GASPar)
    vKappa = (1.0 - vB) * vUncValues
    names(vKappa) = paste("kappa", 1:iK, sep = "")

  }
  vA = unmapVec_C(vA, LowerA(), UpperA())
  names(vA) = paste("a", 1:iK, sep = "")
  vB = unmapVec_C(vB, LowerB(), UpperB())
  names(vB) = paste("b", 1:iK, sep = "")

  return(list(vPw = c(vKappa, vA, vB), StaticFit = StaticFit))

}

starting_vA_Uni <- function(vY, vUncValues, mB, dA_foo, iT, iK, Dist, ScalingType, GASPar) {

  seq_alpha = c(seq(1e-04, 5.5, length.out = 30L))

  vKappa = (diag(iK) - mB) %*% vUncValues

  mA = matrix(0, iK, iK)

  diag(mA)[unlist(GASPar)] = dA_foo

  dAlpha_best = dA_foo

  for (i in 1:iK) {
    if (GASPar[[i]]) {
      for (l in 1:length(seq_alpha)) {
        dLLK_foo = try(GASFilter_univ(vY, vKappa, mA, mB, iT, iK, Dist, ScalingType)$dLLK, silent = TRUE)
        if (is.numeric(dLLK_foo) & !is.nan(dLLK_foo)) {
          mA[i, i] = seq_alpha[l]
          dLLK_post = try(GASFilter_univ(vY, vKappa, mA, mB, iT, iK, Dist, ScalingType)$dLLK,
                          silent = TRUE)
          if (is.numeric(dLLK_post) & !is.nan(dLLK_post)) {
            if (dLLK_post > dLLK_foo)
              dAlpha_best = seq_alpha[l]
          }
          mA[i, i] = dAlpha_best
        }
      }
    }
  }

  return(diag(mA))
}

starting_vB_Uni <- function(vY, vUncValues, dB_foo, mA, iT, iK, Dist, ScalingType, GASPar) {

  seq_beta = c(seq(0.5, 0.98, length.out = 30L))

  mB = matrix(0, iK, iK)

  # diag(mB)[unlist(GASPar)] = dB_foo
  diag(mB) = dB_foo

  vKappa = (diag(iK) - mB) %*% vUncValues

  dB_best = dB_foo

  for (i in 1:iK) {
    if (GASPar[[i]]) {
      for (l in 1:length(seq_beta)) {
        dLLK_foo = try(GASFilter_univ(vY, vKappa, mA, mB, iT, iK, Dist, ScalingType)$dLLK, silent = TRUE)
        if (is.numeric(dLLK_foo) & !is.nan(dLLK_foo)) {
          mB[i, i] = seq_beta[l]
          vKappa = (diag(iK) - mB) %*% vUncValues
          dLLK_post = try(GASFilter_univ(vY, vKappa, mA, mB, iT, iK, Dist, ScalingType)$dLLK,
                          silent = TRUE)
          if (is.numeric(dLLK_post) & !is.nan(dLLK_post)) {
            if (dLLK_post > dLLK_foo)
              dB_best = seq_beta[l]
          }
          mB[i, i] = dB_best
          vKappa = (diag(iK) - mB) %*% vUncValues
        }
      }
    }
  }

  return(diag(mB))
}

StartingValues_mvnorm <- function(mY, iT, iN, iK, GASPar, ScalingType, ScalarParameters) {

  vEmpRho = build_vR(cor(t(mY)), iN)
  vEmpPhi = UnMapR_C(vEmpRho, iN)

  vEmpMu = apply(mY, 1L, mean)
  vEmpSigma = apply(mY, 1L, sd)

  vUncValues = c(vEmpMu, log(vEmpSigma), vEmpPhi)
  names(vUncValues) = paste("kappa.", mvnormParNames(iN), sep = "")
  #

  if (ScalarParameters) {

    mA = starting_mA_Multi_Scalars(mY, vUncValues, mB = diag(rep(0.9, iK)), dA_foo = 0.001, iT,
                                   iK, iN, "mvnorm", ScalingType, GASPar)
    mB = starting_mB_Multi_Scalars(mY, vUncValues, dB_foo = 0.9, mA, iT, iK, iN, "mvnorm", ScalingType,
                                   GASPar)

    vKappa = (diag(iK) - mB) %*% vUncValues
    names(vKappa) = paste("kappa.", mvnormParNames(iN), sep = "")

    vA = diag(mA)[c(1L, iN + 1L, 2L * iN + 1L, 2L * iN + iN * (iN - 1L)/2L + 1L)]
    vB = diag(mB)[c(1L, iN + 1L, 2L * iN + 1L, 2L * iN + iN * (iN - 1L)/2L + 1L)]

    vA = vA[!is.na(vA)]  # nas arise when the distribution do not have shape pars
    vB = vB[!is.na(vB)]  # nas arise when the distribution do not have shape pars

  } else {

    vA = starting_vA_Multi(mY, vUncValues, mB = diag(rep(0.9, iK)), dA_foo = 0.01, iT, iK, iN, "mvnorm",
                           ScalingType, GASPar)
    vB = starting_vB_Multi(mY, vUncValues, dB_foo = 0.9, mA = diag(vA), iT, iK, iN, "mvnorm", ScalingType,
                           GASPar)

    vKappa = (diag(iK) - diag(vB)) %*% vUncValues
    names(vKappa) = paste("kappa.", mvnormParNames(iN), sep = "")
  }

  vA = unmapVec_C(vA, LowerA(), UpperA())
  names(vA) = paste("a.", mvnormParNames(iN, ScalarParameters), sep = "")
  vB = unmapVec_C(vB, LowerB(), UpperB())
  names(vB) = paste("b.", mvnormParNames(iN, ScalarParameters), sep = "")

  pw = c(vKappa, vA, vB)
  return(pw)
}

StartingValues_mvt <- function(mY, iT, iN, iK, GASPar, ScalingType, ScalarParameters, fn.optimizer) {

  StaticFit = StaticMLFIT_Multiv(mY, "mvt", fn.optimizer)

  vUncValues = StaticFit$optimiser$pars
  names(vUncValues) = paste("kappa.", mvtParNames(iN), sep = "")

  if (ScalarParameters) {

    mA = starting_mA_Multi_Scalars(mY, vUncValues, mB = diag(rep(0.9, iK)), dA_foo = 0.001, iT,
                                   iK, iN, "mvt", ScalingType, GASPar)
    mB = starting_mB_Multi_Scalars(mY, vUncValues, dB_foo = 0.9, mA, iT, iK, iN, "mvt", ScalingType,
                                   GASPar)

    vKappa = (diag(iK) - mB) %*% vUncValues
    names(vKappa) = paste("kappa.", mvtParNames(iN), sep = "")

    vA = diag(mA)[c(1L, iN + 1L, 2L * iN + 1L, 2L * iN + iN * (iN - 1L)/2L + 1L)]
    vB = diag(mB)[c(1L, iN + 1L, 2L * iN + 1L, 2L * iN + iN * (iN - 1L)/2L + 1L)]

    vA = vA[!is.na(vA)]  # nas arise when the distribution do not have shape pars
    vB = vB[!is.na(vB)]  # nas arise when the distribution do not have shape pars

  } else {

    vA = starting_vA_Multi(mY, vUncValues, mB = diag(rep(0.9, iK)), dA_foo = 0.01, iT, iK, iN, "mvt",
                           ScalingType, GASPar)
    vB = starting_vB_Multi(mY, vUncValues, dB_foo = 0.9, mA = diag(vA), iT, iK, iN, "mvt", ScalingType,
                           GASPar)

    vKappa = (diag(iK) - diag(vB)) %*% vUncValues
    names(vKappa) = paste("kappa.", mvtParNames(iN), sep = "")
  }

  vA = unmapVec_C(vA, LowerA(), UpperA())
  names(vA) = paste("a.", mvtParNames(iN, ScalarParameters), sep = "")
  vB = unmapVec_C(vB, LowerB(), UpperB())
  names(vB) = paste("b.", mvtParNames(iN, ScalarParameters), sep = "")

  pw = c(vKappa, vA, vB)

  return(pw)
}

MultiGAS_Starting <- function(mY, iT, iN, iK, Dist, GASPar, ScalingType, ScalarParameters, fn.optimizer) {

  if (Dist == "mvnorm")
    vPw = StartingValues_mvnorm(mY, iT, iN, iK, GASPar, ScalingType, ScalarParameters)
  if (Dist == "mvt")
    vPw = StartingValues_mvt(mY, iT, iN, iK, GASPar, ScalingType, ScalarParameters, fn.optimizer)

  return(vPw)
}

starting_vA_Multi <- function(mY, vUncValues, mB, dA_foo, iT, iK, iN, Dist, ScalingType, GASPar) {

  seq_alpha = c(seq(1e-04, 0.5, length.out = 30L))

  vKappa = (diag(iK) - mB) %*% vUncValues

  vBool = FixedDynamicPar_Multi(Dist, iN, GASPar)

  mA = matrix(0, iK, iK)
  diag(mA)[vBool] = dA_foo

  dAlpha_best = dA_foo

  for (i in 1:iK) {

    if (vBool[i]) {

      for (l in 1:length(seq_alpha)) {

        dLLK_foo = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                       silent = TRUE)

        if (is.numeric(dLLK_foo) & !is.nan(dLLK_foo)) {

          mA[i, i] = seq_alpha[l]

          dLLK_post = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                          silent = TRUE)

          if (is.numeric(dLLK_post) & !is.nan(dLLK_post)) {

            if (dLLK_post > dLLK_foo)
              dAlpha_best = seq_alpha[l]

          }

          mA[i, i] = dAlpha_best

        }
      }
    }
  }

  return(diag(mA))
}

starting_mA_Multi_Scalars <- function(mY, vUncValues, mB, dA_foo, iT, iK, iN, Dist, ScalingType, GASPar) {

  seq_alpha = c(seq(1e-04, 0.5, length.out = 30L))

  vKappa = (diag(iK) - mB) %*% vUncValues

  vBool = unlist(GASPar)
  iK2 = length(GASPar)

  mA = matrix(0, iK, iK)
  diag(mA)[FixedDynamicPar_Multi(Dist, iN, GASPar)] = dA_foo

  dAlpha_best = dA_foo

  lStructure = MatrixCoefficientStructure_Multi(Dist, iN, iK, GASPar)

  for (i in 1:iK2) {
    if (vBool[i]) {
      for (l in 1:length(seq_alpha)) {
        dLLK_foo = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                       silent = TRUE)
        if (is.numeric(dLLK_foo) & !is.nan(dLLK_foo)) {
          mA[lStructure[[i]]] = seq_alpha[l]
          dLLK_post = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                          silent = TRUE)
          if (is.numeric(dLLK_post) & !is.nan(dLLK_post)) {
            if (dLLK_post > dLLK_foo)
              dAlpha_best = seq_alpha[l]
          }
          mA[lStructure[[i]]] = dAlpha_best
        }
      }
    }
  }



  return(mA)
}

starting_vB_Multi <- function(mY, vUncValues, dB_foo, mA, iT, iK, iN, Dist, ScalingType, GASPar) {

  seq_beta = c(seq(0.5, 0.98, length.out = 30L))
  vBool = FixedDynamicPar_Multi(Dist, iN, GASPar)

  mB = matrix(0, iK, iK)
  diag(mB) = dB_foo

  vKappa = (diag(iK) - mB) %*% vUncValues

  dB_best = dB_foo

  for (i in 1:iK) {
    if (vBool[i]) {
      for (l in 1:length(seq_beta)) {
        dLLK_foo = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                       silent = TRUE)
        if (is.numeric(dLLK_foo) & !is.nan(dLLK_foo)) {
          mB[i, i] = seq_beta[l]
          vKappa = (diag(iK) - mB) %*% vUncValues
          dLLK_post = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                          silent = TRUE)
          if (is.numeric(dLLK_post) & !is.nan(dLLK_post)) {
            if (dLLK_post > dLLK_foo)
              dB_best = seq_beta[l]
          }
          mB[i, i] = dB_best
          vKappa = (diag(iK) - mB) %*% vUncValues
        }
      }
    }
  }

  return(diag(mB))
}

starting_mB_Multi_Scalars <- function(mY, vUncValues, dB_foo, mA, iT, iK, iN, Dist, ScalingType, GASPar) {

  seq_beta = c(seq(0.5, 0.98, length.out = 30L))

  vBool = unlist(GASPar)
  iK2 = length(GASPar)

  mB = matrix(0, iK, iK)
  diag(mB)[FixedDynamicPar_Multi(Dist, iN, GASPar)] = dB_foo

  vKappa = (diag(iK) - mB) %*% vUncValues

  dBeta_best = dB_foo

  lStructure = MatrixCoefficientStructure_Multi(Dist, iN, iK, GASPar)

  for (i in 1:iK2) {
    if (vBool[i]) {
      for (l in 1:length(seq_beta)) {
        dLLK_foo = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                       silent = TRUE)
        if (is.numeric(dLLK_foo) & !is.nan(dLLK_foo)) {
          mB[lStructure[[i]]] = seq_beta[l]
          vKappa = (diag(iK) - mB) %*% vUncValues
          dLLK_post = try(GASFilter_multi(mY, vKappa, mA, mB, iT, iN, iK, Dist, ScalingType)$dLLK,
                          silent = TRUE)
          if (is.numeric(dLLK_post) & !is.nan(dLLK_post)) {
            if (dLLK_post > dLLK_foo)
              dBeta_best = seq_beta[l]
          }
          mB[lStructure[[i]]] = dBeta_best
          vKappa = (diag(iK) - mB) %*% vUncValues
        }
      }
    }
  }

  return(mB)
}