R/estimation.R

In yymmhaha/PackPaper1: Package related to my first paper

#' @importFrom spams spams.fistaTree spams.normalize
#' @importFrom purrr map_dbl
#' @importFrom magrittr %>%
#' @importFrom glmnet cv.glmnet

NULL
#'
#' 2 step procedure: Variable selection and estimation and prediction improvement
#'
#' @description Implements hierarchical variable selection and plain lasso (glmnet) selection. The hierarchy comes from SPAMS paper
#' (to be described). Gives beta for different lags. The second step is also implemented here.
#'
#' @param X.shift a matrix. Shifted X matrix
#' @param Y.single a vector. Response vector
#' @param Y.shift a matrix. Shifted response vector (for second step)
#' @param lagX a number. Lags of X to be estimated. In our example it's 2,3 or 5
#'
#' @return a list of components
#'
#' \item{beta.hat}{Estimated regression coefficients for specified lags, using hierarchical penalty}
#' \item{beta.lasso}{Estimated regression coefficients for specified lags using lasso}
#' \item{phi.hat}{Estimated AR coefficients for residual fitting}
#' \item{beta.refit}{Combination of phi and beta}
#'
#'
#' @export
#'
#' @examples
#' modelsimple <- estimation(X.shift = DATsimple$Xtr,
#' Y.single = DATsimple$Ytr.single,
#' Y.shift = DATsimple$Ytr,
#' lagX = parameter.of.estimation$lag.to.estX)



estimation <- function(X.shift, Y.single, Y.shift, lagX){

  W0 <- matrix(c(0), nrow = ncol(X.shift), ncol = 1)
  tree <- autogrouper(X.shift, lx = lagX)$tree  # X is the subset of training data with correct lags


  # ---------------- choose the best lambda
  # here it's the choose.lambda with intercept

  optimal.lambda <- choose.lambda(X.shift, Y.single, W0 = W0, tree = tree) %>%  .$optimal.lambda

  meanMatrix <- matrix(rep(colMeans(X.shift), nrow(X.shift)),  nrow(X.shift), ncol(X.shift), byrow = T) # every row same
  Xtr.hier <- X.shift - meanMatrix   # centered
  Xtr.hier <- spams.normalize(Xtr.hier)  # rescale such that they have unit L2 norm. but this is different from unit variance

  # ---------------- use the chosen lambda
  # normalised X and scaled Y
  res <- spams.fistaTree(as.matrix(scale(Y.single)), Xtr.hier, W0, tree, TRUE, numThreads = -1,
                         verbose = F, lambda1 = optimal.lambda, it0 = 10,
                         max_it = 200, L0 = 0.1, tol = 1e-5,
                         intercept = F, pos = F, compute_gram = T,
                         loss = 'square', regul = 'tree-l2')  #extra attention to the regul!!!
  beta.hat <- res[[1]]  # dense p*n matrix (here p*1 as we only have one response)



  # ------------ alternatively, check with glmnet
  res.lasso <- cv.glmnet(X.shift, Y.single, alpha = 1)
  beta.lasso <- res.lasso %>% coef(., s = .$lambda.1se) %>% .[-1]



  # ============== RESIDUAL FITTING ============= #

  residual <- Y.single - X.shift %*% beta.hat   # order: true - fitted
  refit <- cv.glmnet(Y.shift, residual, alpha = 1)  # give all 5 lags in y, let lasso choose
  phi.hat <- refit %>% coef(., s = .$lambda.1se) %>% .[-1]
  beta.refit <- c(phi.hat, beta.hat)  # 5 + 21

  return(list(beta.hat = beta.hat,
              beta.lasso = beta.lasso,
              beta.refit = beta.refit,
              phi.hat = phi.hat))

}