R/lagspec.R

In mpiktas/midasr: Mixed Data Sampling Regression

##' Normalized Exponential Almon lag MIDAS coefficients
##'
##' Calculate normalized exponential Almon lag coefficients given the parameters and required number of coefficients.
##' @param p parameters for Almon lag
##' @param d number of the coefficients
##' @param m the frequency, currently ignored.
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @examples
##'
##' ##Load data
##' data("USunempr")
##' data("USrealgdp")
##'
##' y <- diff(log(USrealgdp))
##' x <- window(diff(USunempr),start=1949)
##' t <- 1:length(y)
##'
##' midas_r(y~t+fmls(x,11,12,nealmon),start=list(x=c(0,0,0)))
##'
##' @details Given unrestricted MIDAS regression
##'
##' \deqn{y_t=\sum_{h=0}^d\theta_{h}x_{tm-h}+\mathbf{z_t}\beta+u_t}
##'
##' normalized exponential Almon lag restricts the coefficients \eqn{theta_h} in the following way:
##'
##' \deqn{\theta_{h}=\delta\frac{\exp(\lambda_1(h+1)+\dots+
##' \lambda_r(h+1)^r)}{\sum_{s=0}^d\exp(\lambda_1(s+1)+\dots+\lambda_r(h+1)^r)}}
##'
##' The parameter \eqn{\delta} should be the first element in vector \code{p}. The degree of
##' the polynomial is then decided by the number of the remaining parameters.
##' @importFrom stats poly
##' @export
nealmon <- function(p, d, m) {
  i <- 1:d
  plc <- poly(i, degree = length(p) - 1, raw = TRUE) %*% p[-1]
  as.vector(p[1] * exp(plc) / sum(exp(plc)))
}


##' Produces weights for aggregates based MIDAS regression
##'
##' Suppose a weight function \eqn{w(\beta,\theta)} satisfies the following equation:
##' \deqn{w(\beta,\theta)=\beta g(\theta)}
##'
##' The following combinations are defined, corresponding to structure types
##' \code{A}, \code{B} and \code{C} respectively:
##' \deqn{(w(\beta_1,\theta_1),...,w(\beta_k,\theta_k))}
##' \deqn{(w(\beta_1,\theta),...,w(\beta_k,\theta))}
##' \deqn{\beta(w(1,\theta),...,w(1,\theta)),}
##'
##' where \eqn{k} is the number of low frequency lags, i.e. \eqn{d/m}. If the latter
##' value is not whole number, the error is produced.
##'
##'
##' The starting values \eqn{p} should be supplied then as follows:
##' \deqn{(\beta_1,\theta_1,...,\beta_k,\theta_k)}
##' \deqn{(\beta_1,...,\beta_k,\theta)}
##' \deqn{(\beta,\theta)}
##'
##'
##' @title Weights for aggregates based MIDAS regressions
##' @param p parameters for weight functions, see details.
##' @param d number of high frequency lags
##' @param m the frequency
##' @param weight the weight function
##' @param type type of structure, a string, one of A, B or C.
##' @return a vector of weights
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
amweights <- function(p, d, m, weight = nealmon, type = c("A", "B", "C")) {
  type <- match.arg(type)
  hf <- d %/% m

  if (d %% m != 0) stop("Number of high frequency lags should be a multiple of frequency")
  if (type == "A") {
    return(as.vector(apply(matrix(seq_len(length(p)), ncol = hf), 2, function(x) weight(p[x], m, m))))
  }
  if (type == "B") {
    theta <- p[-1:-hf]
    return(as.vector(sapply(p[1:hf], function(beta) weight(c(beta, theta), m))))
  }
  if (type == "C") {
    return(p[1] * as.vector(apply(matrix(2:length(p), ncol = hf), 2, function(x) weight(c(1, p[x]), m, m))))
  }
}


##' Gradient function for normalized exponential Almon lag weights
##'
##' Gradient function for normalized exponential Almon lag weights
##' @param p hyperparameters for Almon lag
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return the gradient matrix
##' @author Vaidotas Zemlys
##' @export
nealmon_gradient <- function(p, d, m) {
  i <- 1:d
  pl <- poly(i, degree = length(p) - 1, raw = TRUE)
  eplc <- exp(pl %*% p[-1])[, , drop = TRUE]
  ds <- apply(pl * eplc, 2, sum)
  s <- sum(eplc)
  cbind(eplc / s, p[1] * (pl * eplc / s - eplc %*% t(ds) / s^2))
}

##' Normalized beta probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized beta probability density function specification
##' @param p parameters for normalized beta probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
nbeta <- function(p, d, m) {
  nbetaMT(c(p[1:3], 0), d, m)
}

##' Gradient function for normalized beta probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized beta probability density function specification of MIDAS weights.
##' @param p parameters for normalized beta probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
nbeta_gradient <- function(p, d, m) {
  nbetaMT_gradient(c(p[1:3], 0), d, m)[, 1:3]
}


##' Normalized beta probability density function MIDAS weights specification (MATLAB toolbox compatible)
##'
##' Calculate MIDAS weights according to normalized beta probability density function specification.
##' Compatible with the specification in MATLAB toolbox.
##' @param p parameters for normalized beta probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
nbetaMT <- function(p, d, m) {
  eps <- .Machine$double.eps
  xi <- (1:d - 1) / (d - 1)
  xi[1] <- xi[1] + eps
  xi[d] <- xi[d] - eps
  nb <- xi^(p[2] - 1) * (1 - xi)^(p[3] - 1) #nolint
  if (sum(nb) < eps) {
    if (abs(p[4]) < eps) {
      rep(0, d)
    } else {
      p[1] * rep(1 / d, length(nb))
    }
  } else {
    w <- (nb / sum(nb) + p[4])
    p[1] * w / sum(w)
  }
}

##' Gradient function for normalized beta probability density function MIDAS weights specification
##' (MATLAB toolbox compatible)
##'
##' Calculate gradient function for normalized beta probability density function specification of MIDAS weights.
##' @param p parameters for normalized beta probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
nbetaMT_gradient <- function(p, d, m) {
  eps <- .Machine$double.eps
  xi <- (1:d - 1) / (d - 1)
  xi[1] <- xi[1] + eps
  xi[d] <- xi[d] - eps
  nb <- xi^(p[2] - 1) * (1 - xi)^(p[3] - 1) #nolint
  snb <- sum(nb)
  wnb <- 1 + d * p[4]
  if (snb > eps) {
    nba <- nb * log(xi)
    nbb <- nb * log(1 - xi)
    a <- nba / snb - nb * sum(nba) / snb^2
    b <- nbb / snb - nb * sum(nbb) / snb^2
    cbind((nb / snb + p[4]) / wnb, p[1] * a / wnb, p[1] * b / wnb, p[1] * (1 - nb / snb * d) / wnb^2)
  }
  else {
    gres <- matrix(0, nrow = d, ncol = 4)
    if (abs(p[4]) > eps) gres[, 1] <- 1 / d
    gres
  }
}

##' Almon polynomial MIDAS weights specification
##'
##' Calculate Almon polynomial MIDAS weights
##' @param p parameters for Almon polynomial weights
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
almonp <- function(p, d, m) {
  i <- 1:d
  plc <- poly(i, degree = length(p) - 1, raw = TRUE) %*% p[-1] + p[1]
  as.vector(plc)
}

##' Gradient function for Almon polynomial MIDAS weights
##'
##' Calculate gradient for Almon polynomial MIDAS weights specification
##' @param p vector of parameters for Almon polynomial specification
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Vaidotas Zemlys
##' @export
almonp_gradient <- function(p, d, m) {
  i <- 1:d
  plc <- poly(i, degree = length(p) - 1, raw = TRUE)
  cbind(1, plc)
}

##' Step function specification for MIDAS weights
##'
##' Step function specification for MIDAS weights
##' @param p vector of parameters
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @param a vector of increasing positive integers indicating the steps
##' @return vector of coefficients
##' @author Vaidotas Zemlys
##' @export
polystep <- function(p, d, m, a) {
  if (length(a) != length(p) - 1) stop("The number of steps should be number of parameters minus one")
  if (min(a) <= 1 | max(a) >= d) stop("The steps are out of bounds")
  a <- c(0, a, d)
  rep(p, times = diff(a))
}
##' Gradient of step function specification for MIDAS weights
##'
##' Gradient of step function specification for MIDAS weights
##' @param p vector of parameters
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @param a vector of increasing positive integers indicating the steps
##' @return vector of coefficients
##' @author Vaidotas Zemlys
##' @export
polystep_gradient <- function(p, d, m, a) {
  if (length(a) != length(p) - 1) stop("The number of steps should be number of parameters minus one")
  if (min(a) <= 1 | max(a) >= d) stop("The steps are out of bounds")
  a <- c(0, a, d)
  da <- diff(a)
  r <- cbind(a[-length(a)], da)
  apply(r, 1, function(x) {
    res <- rep(0, d)
    res[(x[1] + 1):(x[1] + x[2])] <- 1
    res
  })
}

##' Normalized Gompertz probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized Gompertz probability density function specification
##' @param p parameters for normalized Gompertz probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
gompertzp <- function(p, d, m) {
  i <- 1:d / d
  gm <- exp(p[3] * i - p[2] * exp(p[3] * i))
  p[1] * gm / sum(gm)
}

##' Gradient function for normalized Gompertz probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized Gompertz probability density function specification of MIDAS weights.
##' @param p parameters for normalized Gompertz probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
gompertzp_gradient <- function(p, d, m) {
  i <- 1:d / d
  gm <- exp(p[3] * i - p[2] * exp(p[3] * i))
  dp2 <- -gm * exp(i * p[3])
  dp3 <- -gm * i * (p[2] * exp(i * p[3]) - 1)
  cbind(
    gm, p[1] * (dp2 - gm * sum(dp2) / sum(gm)),
    p[1] * (dp3 - gm * sum(dp3) / sum(gm))
  ) / sum(gm)
}

##' Normalized Nakagami probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized Nakagami probability density function specification
##' @param p parameters for normalized Nakagami probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
nakagamip <- function(p, d, m) {
  i <- 1:d / d
  ng <- i^(2 * p[2] - 1) * exp(-p[2] / p[3] * i^2) #nolint
  p[1] * ng / sum(ng)
}

##' Gradient function for normalized Nakagami probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized Nakagami probability density function specification of MIDAS weights.
##' @param p parameters for normalized Nakagami probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
nakagamip_gradient <- function(p, d, m) {
  i <- 1:d / d
  ng <- i^(2 * p[2] - 1) * exp(-p[2] / p[3] * i^2) #nolint
  dp2 <- ((2 * log(i) * p[3] - i^2) / p[3]) * ng
  dp3 <- (p[2] * i^2 / p[3]^2) * ng
  cbind(
    ng, (dp2 - ng * sum(dp2) / sum(ng)) * p[1],
    (dp3 - ng * sum(dp3) / sum(ng)) * p[1]
  ) / sum(ng)
}

##' Normalized log-Cauchy probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized log-Cauchy probability density function specification
##' @param p parameters for normalized log-Cauchy probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
lcauchyp <- function(p, d, m) {
  i <- 1:d / d
  lc <- 1 / (i * ((log(i) - p[2])^2 + p[3]^2))
  p[1] * lc / sum(lc)
}

##' Gradient function for normalized log-Cauchy probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized log-Cauchy probability density function specification of MIDAS weights.
##' @param p parameters for normalized log-Cauchy probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
lcauchyp_gradient <- function(p, d, m) {
  i <- 1:d / d
  lc <- 1 / (i * ((log(i) - p[2])^2 + p[3]^2))
  dp2 <- 2 * lc^2 * i * (log(i) - p[2])
  dp3 <- -2 * lc^2 * p[3] * i
  cbind(
    lc, p[1] * (dp2 - lc * sum(dp2) / sum(lc)),
    p[1] * (dp3 - lc * sum(dp3) / sum(lc))
  ) / sum(lc)
}

##' HAR(3)-RV model MIDAS weights specification
##'
##' MIDAS weights for Heterogeneous Autoregressive model of Realized Volatilty (HAR-RV).
##' It is assumed that month has 20 days.
##' @title HAR(3)-RV model MIDAS weights specification
##' @param p parameters for Almon lag
##' @param d number of the coefficients
##' @param m the frequency, currently ignored.
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @references Corsi, F., \emph{A Simple Approximate Long-Memory Model of Realized Volatility},
##' Journal of Financial Econometrics Vol. 7 No. 2 (2009) 174-196
##' @export
harstep <- function(p, d, m) {
  if (d != 20) stop("HAR(3)-RV process requires 20 lags")
  out <- rep(0, 20)
  out[1] <- p[1] + p[2] / 5 + p[3] / 20
  out[2:5] <- p[2] / 5 + p[3] / 20
  out[6:20] <- p[3] / 20
  out
}
##' Gradient function for HAR(3)-RV model MIDAS weights specification
##'
##' MIDAS weights for Heterogeneous Autoregressive model of Realized Volatilty (HAR-RV).
##' It is assumed that month has 20 days.
##' @title Gradient function for HAR(3)-RV model MIDAS weights specification
##' @param p parameters for Almon lag
##' @param d number of the coefficients
##' @param m the frequency, currently ignored.
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @references Corsi, F., \emph{A Simple Approximate Long-Memory Model of Realized Volatility},
##' Journal of Financial Econometrics Vol. 7 No. 2 (2009) 174-196
##' @export
harstep_gradient <- function(p, d, m) {
  if (d != 20) stop("HAR(3)-RV process requires 20 lags")
  out <- matrix(0, ncol = 3, nrow = d)
  out[1, 1] <- 1
  out[1:5, 2] <- 1 / 5
  out[1:20, 3] <- 1 / 20
  out
}


#' Calculates the MIDAS coefficients for generalized exponential MIDAS lag specification
#'
#' Generalized exponential MIDAS lag specification is a generalization of exponential
#' Almon lag. It is defined as a product of first order polynomial with exponent
#' of the second order polynomial. This spefication was used by V. Kvedaras and
#' V. Zemlys (2012).
#'
#' @title Generalized exponential MIDAS coefficients
#' @param p a vector of parameters
#' @param d number of coefficients
#' @param m the frequency, currently ignored
#'
#' @return a vector of coefficients
#' @export
#' @author Virmantas Kvedaras, Vaidotas Zemlys
#' @references Kvedaras V., Zemlys, V. \emph{Testing the functional constraints on
#' parameters in regressions with variables of different frequency} Economics Letters 116 (2012) 250-254
genexp <- function(p, d, m) {
  i <- (1:d - 1) / 100
  pol <- p[3] * i + p[4] * i^2
  (p[1] + p[2] * i) * exp(pol)
}

#' Calculates the gradient of generalized exponential MIDAS lag specification
#'
#' Generalized exponential MIDAS lag specification is a generalization of exponential
#' Almon lag. It is defined as a product of first order polynomial with exponent
#' of the second order polynomial. This spefication was used by V. Kvedaras and
#' V. Zemlys (2012).
#'
#' @title Gradient of generalized exponential MIDAS coefficient generating function
#' @param p a vector of parameters
#' @param d number of coefficients
#' @param m the frequency, currently ignored
#'
#' @return a vector of coefficients
#' @export
#' @author Virmantas Kvedaras, Vaidotas Zemlys
#' @references Kvedaras V., Zemlys, V. \emph{Testing the functional constraints on parameters in
#' regressions with variables of different frequency} Economics Letters 116 (2012) 250-254
genexp_gradient <- function(p, d, m) {
  i <- (1:d - 1) / 100
  pol <- p[3] * i + p[4] * i^2
  pol0 <- p[1] + p[2] * i
  epl <- exp(pol)
  cbind(epl, epl * i, pol0 * epl * i, pol0 * epl * i^2)
}

is_weight_normalized <- function(wf, init, chk = rnorm(5), tol = 1e-8, nc = 1) {
  tst <- lapply(c(0, rnorm(chk)), function(x) {
    ss <- init
    ss[nc] <- ss[nc] + x
    try(abs(sum(wf(ss)) - ss[nc]))
  })

  fail <- sapply(tst, inherits, "try-error")
  got_na <- fail
  got_na[!fail] <- sapply(tst[!fail], is.na)

  if (all(fail) | all(got_na)) stop("Error in checking normalization, all candidate starting values produced errors")

  eq <- unlist(tst[!got_na])
  if (all(eq < tol)) {
    return(TRUE)
  } else {
    return(FALSE)
  }
}

find_normalisation_coefficient <- function(wf, init) {
  tst <- lapply(seq_len(length(init)), function(i) try(is_weight_normalized(wf, init, nc = i)))
  tst <- sapply(tst, function(l) if (inherits(l, "try-error")) {
    return(FALSE)
  } else {
    return(l)
  })
  which(tst)
}