R/hmoment.R

#' @include hspec.R
NULL

setGeneric("get_lambda0", function(object, ...) standardGeneric("get_lambda0"))
# Get the long-run expectation of intensities - lambda matrix
#
# This function is not crucial but useful when you want to set lambda0 by default.
#
setMethod(
  f = "get_lambda0",
  signature(object = "hspec"),
  definition = function(object, inter_arrival = inter_arrival, type = type, mark = mark,
                        N = N, Nc = Nc,
                        mu = mu, alpha = alpha, beta = beta,
                        lambda = NULL, lambda_component = NULL, lambda_component_n = NULL,
                        ...){

    dimens <- object@dimens

    # parameter setting
    if (is.function(object@mu)){
      if(length(formals(object@mu)) == 1){
        mu0 <- evalf(object@mu)
      } else {
        mu0 <- mu(n = 1, inter_arrival = inter_arrival, type = type, mark = mark,
                  N = N, Nc = Nc, lambda = lambda, lambda_component = lambda_component,
                  lambda_component_n = lambda_component_n,
                  alpha = alpha, beta = beta)
      }
    } else{
      mu0 <- object@mu
    }
    if (is.function(object@alpha)){
      alpha <- evalf(object@alpha)
    } else{
      alpha <- object@alpha
    }
    if (is.function(object@beta)){
      beta <- evalf(object@beta)
    } else{
      beta <- object@beta
    }


    if (dimens == 1){
      lamb0 <- (mu[1]*beta[1]/(beta[1]- alpha[1]) - mu[1])/2
      LAMBDA0 <- matrix(lamb0, nrow=1)

    } else {
      # need to handle when the matrix is singular
      LAMBDA0 <- matrix(rep(0, dimens^2), nrow=dimens)
      LAMBDA0 <- tryCatch({
          LAMBDA_st <- solve(diag(dimens) - alpha / beta) %*% mu0
          matrix(rep(LAMBDA_st, dimens), nrow=dimens, byrow=T) * alpha / beta
        },
        error = function(e){
          warning("Due to the singualr martrix in caculcation, set initial value of lambda to mu.")
          matrix(rep(0, dimens^2), nrow=dimens)
        }
      )
    }
    LAMBDA0
  }
)


#' Compute Hawkes volatility
#'
#' This function computes Hawkes volatility. Only works for bi-variate Hawkes process.
#'
#' @param object \code{\link{hspec-class}}
#' @param horizon Time horizon for volatility.
#' @param inter_arrival Inter-arrival times of events which includes inter-arrival for events that occur in all dimensions. Start with zero.
#' @param type  A vector of dimensions. Distinguished by numbers, 1, 2, 3, and so on. Start with zero.
#' @param mark A vector of mark (jump) sizes. Start with zero.
#' @param dependence Dependence between mark and previous sigma-algebra.
#' @param lambda_component0 A matrix of the starting values of lambda component.
#' @param ... Further arguments passed to or from other methods.
#'
#' @rdname hvol
#' @export
setGeneric("hvol", function(object, horizon = 1,
                           inter_arrival = NULL,
                           type = NULL,
                           mark = NULL,
                           dependence = FALSE,
                           lambda_component0 = NULL, ...) standardGeneric("hvol"))

#' @rdname hvol
setMethod(
  f = "hvol",
  signature(object = "hspec"),
  definition = function(object, horizon = 1,
                        inter_arrival = NULL,
                        type = NULL,
                        mark = NULL,
                        dependence = FALSE,
                        lambda_component0 = NULL,
                        ...){

    additional_argument <- list(...)
    if ("lambda0" %in% names(additional_argument)) {

      warning("lambda0 is deprecated; instead use lambda_component0.")

      lambda_component0 <- additional_argument[["lambda0"]]

    }

    h <- object

    infer_result <- infer_lambda(h, inter_arrival = inter_arrival, type = type, mark = mark,
                                 lambda_component0 = lambda_component0)

    Mu <- matrix(h@mu, nrow=2)
    Alpha <- h@alpha
    Beta <- diag(diag(h@beta))
    # Todo : if eta does not exist
    Eta <- h@eta

    up_move_idx <- type == 1
    down_move_idx <- type == 2


    #### compute Hawkes vol

    if (dependence){

      Z_1 <- sum(mark[up_move_idx] * infer_result$lambda[up_move_idx,1]) /
        sum(infer_result$lambda[up_move_idx,1])

      Z_2 <- sum(mark[down_move_idx] * infer_result$lambda[down_move_idx,2]) /
        sum(infer_result$lambda[down_move_idx,2])

      Z_sq_1 <- sum(mark[up_move_idx]^2 * infer_result$lambda[up_move_idx,1]) /
        sum(infer_result$lambda[up_move_idx,1])

      Z_sq_2 <- sum(mark[down_move_idx]^2 * infer_result$lambda[down_move_idx,2]) /
        sum(infer_result$lambda[down_move_idx,2])

      Z_ln_11 <- sum(mark[up_move_idx] * infer_result$lambda[up_move_idx,1] * infer_result$Nc[up_move_idx,1]) /
        sum(infer_result$lambda[up_move_idx,1] * infer_result$Nc[up_move_idx,1])

      Z_ln_12 <- sum(mark[down_move_idx] * infer_result$lambda[down_move_idx,1] * infer_result$Nc[down_move_idx,2]) /
        sum(infer_result$lambda[down_move_idx,1] * infer_result$Nc[down_move_idx,2])

      Z_ln_21 <- sum(mark[up_move_idx] * infer_result$lambda[up_move_idx,2] * infer_result$Nc[up_move_idx,1]) /
        sum(infer_result$lambda[up_move_idx,2] * infer_result$Nc[up_move_idx,1])

      Z_ln_22 <- sum(mark[down_move_idx] * infer_result$lambda[down_move_idx,2] * infer_result$Nc[down_move_idx,2]) /
        sum(infer_result$lambda[down_move_idx,2] * infer_result$Nc[down_move_idx,2])


      Z_nl_11 <- sum(mark[up_move_idx] * infer_result$Nc[up_move_idx,1] * infer_result$lambda[up_move_idx,1]) /
        sum(infer_result$Nc[up_move_idx,1] * infer_result$lambda[up_move_idx,1])

      Z_nl_12 <- sum(mark[down_move_idx] * infer_result$Nc[down_move_idx,1] * infer_result$lambda[down_move_idx,2]) /
        sum(infer_result$Nc[down_move_idx,1] * infer_result$lambda[down_move_idx,2])

      Z_nl_21 <- sum(mark[up_move_idx] * infer_result$Nc[up_move_idx,2] * infer_result$lambda[up_move_idx,1]) /
        sum(infer_result$Nc[up_move_idx,2] * infer_result$lambda[up_move_idx,1])

      Z_nl_22 <- sum(mark[down_move_idx] * infer_result$Nc[down_move_idx,2] * infer_result$lambda[down_move_idx,2]) /
        sum(infer_result$Nc[down_move_idx,2] * infer_result$lambda[down_move_idx,2])


      Z_ll_11 <- sum(mark[up_move_idx] * infer_result$lambda[up_move_idx,1] * infer_result$lambda[up_move_idx,1]) /
        sum(infer_result$lambda[up_move_idx,1] * infer_result$lambda[up_move_idx,1])

      Z_ll_12 <- sum(mark[down_move_idx] * infer_result$lambda[down_move_idx,1] * infer_result$lambda[down_move_idx,2]) /
        sum(infer_result$lambda[down_move_idx,1] * infer_result$lambda[down_move_idx,2])

      Z_ll_21 <- sum(mark[up_move_idx] * infer_result$lambda[up_move_idx,2] * infer_result$lambda[up_move_idx,1]) /
        sum(infer_result$lambda[up_move_idx,2] * infer_result$lambda[up_move_idx,1])

      Z_ll_22 <- sum(mark[down_move_idx] * infer_result$lambda[down_move_idx,2] * infer_result$lambda[down_move_idx,2]) /
        sum(infer_result$lambda[down_move_idx,2] * infer_result$lambda[down_move_idx,2])


      Z <- matrix(c(Z_1, Z_2,
                    Z_1, Z_2), nrow=2, byrow=TRUE)

      Z_sq <- matrix(c(Z_sq_1, Z_sq_2,
                       Z_sq_1, Z_sq_2), nrow=2, byrow=TRUE)

      Z_ln <- matrix(c(Z_ln_11, Z_ln_12,
                       Z_ln_21, Z_ln_22), nrow=2, byrow=TRUE)

      Z_nl <- matrix(c(Z_nl_11, Z_nl_12,
                       Z_nl_21, Z_nl_22), nrow=2, byrow=TRUE)

      Z_ll <- matrix(c(Z_ll_11, Z_ll_12,
                       Z_ll_21, Z_ll_22), nrow=2, byrow=TRUE)

      Z_nl <- Z


      E_lambda <- solve(Beta - Alpha + Eta - Eta * Z) %*% Beta %*% Mu

      # To find E_ll : E_lambda_lambda

      G <- (Alpha - Eta + Eta * Z) %*% diag(c(E_lambda)) %*% t(Alpha - Eta) +
        (Alpha - Eta) %*% diag(c(E_lambda)) %*% t(Eta * Z) +
        (Eta * Z_sq^(0.5)) %*% diag(c(E_lambda)) %*% t(Eta * Z_sq^(0.5))


      H <- kronecker(diag(2), Alpha - Beta) + kronecker(Alpha - Beta, diag(2)) +
        rbind(cbind(Eta %*% diag(t(Z_ll)[,1] - 1), matrix(rep(0,4), nrow=2)),
              cbind(matrix(rep(0,4), nrow=2), Eta %*% diag(t(Z_ll)[,2] - 1))) +
        rbind(cbind(t(Eta)[,1] %*% diag(t(Z_ll)[,1] - 1), matrix(rep(0,2), nrow=1)),
              cbind(matrix(rep(0,2), nrow=1), t(Eta)[,1] %*% diag(t(Z_ll)[,2] - 1)),
              cbind(t(Eta)[,2] %*% diag(t(Z_ll)[,1] - 1), matrix(rep(0,2), nrow=1)),
              cbind(matrix(rep(0,2), nrow=1), t(Eta)[,2] %*% diag(t(Z_ll)[,2] - 1)))


      C <- Beta %*% Mu %*% t(E_lambda) + t(Beta %*% Mu  %*% t(E_lambda)) + G

      E_ll <- matrix(solve(H, -c(C)), nrow=2)

      HA <- kronecker(diag(2), Alpha - Beta) +
        rbind(cbind(Eta %*% diag(t(Z)[,1] - 1), matrix(rep(0,4), nrow=2)),
              cbind(matrix(rep(0,4), nrow=2), Eta %*% diag(t(Z)[,2] - 1)))

      CA <- (Beta %*% Mu) %*% t(E_lambda) %*% diag(diag(Z))

      A <- t(matrix(solve(HA, -c(CA)), nrow=2))
      # A <- diag(diag(Z)) %*% E_lambda %*% t(E_lambda)

      CB <- Z_ll * E_ll + ((Alpha - Eta) * Z + Eta * Z_sq) %*% diag(c(E_lambda)) - t(A)

      B <- t(matrix(solve(HA, -c(CB)), nrow=2))

      E_nl <- A * horizon + B

      P1 <- Z_nl * 0.5 * A + t(Z_nl * 0.5 * A)
      P2 <- Z_nl * B  + t(Z_nl * B)
      P3 <- Z_sq * diag(c(E_lambda))

      E_NN <- P1 * horizon^2 + P2 * horizon + P3 * horizon

      variance <- matrix(c(1, -1), nrow=1) %*% (P2 * horizon + P3 * horizon) %*% matrix(c(1, -1))

      (std <- sqrt(variance))

    } else {

      ##### independent

      Z_ind_1 <- mean(mark[up_move_idx])
      Z_ind_2 <- mean(mark[down_move_idx])

      Z_sq_ind_1 <- mean(mark[up_move_idx]^2)
      Z_sq_ind_2 <- mean(mark[down_move_idx]^2)

      Z_ind <- matrix(c(Z_ind_1, Z_ind_2,
                        Z_ind_1, Z_ind_2), nrow=2, byrow=TRUE)

      Z_sq_ind <-  matrix(c(Z_sq_ind_1, Z_sq_ind_2,
                            Z_sq_ind_1, Z_sq_ind_2), nrow=2, byrow=TRUE)

      E_lambda_ind <- solve(Beta - Alpha + Eta - Eta * Z_ind) %*% Beta %*% Mu

      G_ind <- (Alpha - Eta + Eta * Z_ind) %*% diag(c(E_lambda_ind)) %*% t(Alpha - Eta) +
        (Alpha - Eta) %*% diag(c(E_lambda_ind)) %*% t(Eta * Z_ind) +
        (Eta * Z_sq_ind^(0.5)) %*% diag(c(E_lambda_ind)) %*% t(Eta * Z_sq_ind^(0.5))

      H_ind <- kronecker(diag(2), Alpha - Beta) + kronecker(Alpha - Beta, diag(2)) +
        rbind(cbind(Eta %*% diag(t(Z_ind)[,1] - 1), matrix(rep(0,4), nrow=2)),
              cbind(matrix(rep(0,4), nrow=2), Eta %*% diag(t(Z_ind)[,2] - 1))) +
        rbind(cbind(t(Eta)[,1] %*% diag(t(Z_ind)[,1] - 1), matrix(rep(0,2), nrow=1)),
              cbind(matrix(rep(0,2), nrow=1), t(Eta)[,1] %*% diag(t(Z_ind)[,2] - 1)),
              cbind(t(Eta)[,2] %*% diag(t(Z_ind)[,1] - 1), matrix(rep(0,2), nrow=1)),
              cbind(matrix(rep(0,2), nrow=1), t(Eta)[,2] %*% diag(t(Z_ind)[,2] - 1)))

      C_ind <- Beta %*% Mu %*% t(E_lambda_ind) + t(Beta %*% Mu  %*% t(E_lambda_ind)) + G_ind

      E_ll_ind <- matrix(solve(H_ind, -c(C_ind)), nrow=2)

      HA_ind <- kronecker(diag(2), Alpha - Beta) +
        rbind(cbind(Eta %*% diag(t(Z_ind)[,1] - 1), matrix(rep(0,4), nrow=2)),
              cbind(matrix(rep(0,4), nrow=2), Eta %*% diag(t(Z_ind)[,2] - 1)))

      CA_ind <- (Beta %*% Mu) %*% t(E_lambda_ind) %*% diag(diag(Z_ind))

      A_ind <- t(matrix(solve(HA_ind, -c(CA_ind)), nrow=2))

      CB_ind <- Z_ind * E_ll_ind + ((Alpha - Eta) * Z_ind + Eta * Z_sq_ind) %*% diag(c(E_lambda_ind)) - t(A_ind)

      B_ind <- t(matrix(solve(HA_ind, -c(CB_ind)), nrow=2))

      E_nl_ind <- A_ind * horizon + B_ind

      P1_ind <- Z_ind * 0.5 * A_ind + t(Z_ind * 0.5 * A_ind)
      P2_ind <- Z_ind * B_ind  + t(Z_ind * B_ind)
      P3_ind <- diag(diag(Z_sq_ind)) %*% diag(c(E_lambda_ind))

      E_NN_ind <- P1_ind * horizon^2 + P2_ind * horizon + P3_ind * horizon

      variance_ind <- matrix(c(1, -1), nrow=1) %*% (P2_ind * horizon + P3_ind * horizon) %*% matrix(c(1, -1))

      (std_ind <- sqrt(variance_ind))

    }

  }
)
ksublee/emhawkes documentation built on Feb. 2, 2023, 12:36 p.m.