# R/erfeVecR.R In erfe: Fits Expectile Regression for Panel Fixed Effect Model

#### Documented in erfeVecR

#' Dexpectilize a vector according the a single asymmetric point
#' @description This function is part of the erfe package. It estimates
#'  the ERFE model for a panel dataset and for a
#'  single asymmetric point \eqn{\tau \in (0, 1)}.
#'  When \eqn{\tau=0.5} the function estimate the classical
#'  within-transformation estimator and its sandwich
#'  covariance matrix.
#' @return Return a list of objects related to the erfe model
#' such as the asymmetric point, the coefficient-estimate,
#' the standard deviation, the estimated covariance.
#' @references Barry, Amadou, Oualkacha, Karim, and Charpentier
#'  Arthur. (2022). \emph{Weighted asymmetric least squares
#'   regression with fixed-effects}.
#'  arXiv preprint arXiv:2108.04737
#' @param xmat Numeric vector to de-expectilize.
#' @param yvec Numeric vector of individual asymmetric weight.
#' @param panSizeVec Numeric vector to individual panel size.
#' @param asym Numeric vector to individual panel size.
#' @param id Numeric vector to individual panel size.
#' @examples
#' set.seed(13)
#' temps_obs <- 5
#' n_subj <- 50
#' sig <- diag(rep(1,temps_obs))
#' id <- rep(1:n_subj, each=temps_obs)
#' rvec <- c(mvtnorm::rmvnorm(n_subj, sigma = sig))
#' fvec <- (1 + rep(rnorm(n_subj) , each=temps_obs))
#' xmat <- cbind(rt(n_subj*temps_obs, df=2, ncp=1.3),
#'  1.2 * fvec + rnorm(n_subj * temps_obs, mean = 0.85, sd = 1.5) )
#' yvec <- 0.6*xmat[, 1] + xmat[, 2] + fvec + rvec
#' asym <- 0.5
#' panSizeVec <- unname(unlist(lapply(split(id, id), function(x) length(x))))
#' erfeVecR(xmat, yvec, panSizeVec, asym, id)
#' @export
#' @importFrom mvtnorm rmvnorm
#' @importFrom stats lsfit
#' @importFrom Matrix bdiag
erfeVecR <- function(xmat, yvec, panSizeVec, asym, id) {
xncol <- ncol(xmat)
betaEstimate <- rep(NA, xncol)
it <- 1
dw1 <- 1
w1 <- 0 * yvec + 0.5
while ( dw1 != 0 && it < 50 ) {
ydexpect <- dexpectilizeVecR(yvec, w1, panSizeVec)
xdexpect <- dexpectilizeMatR(xmat, w1, panSizeVec)
lsfit_coef <- stats::lsfit(xdexpect, ydexpect,
wt = w1, intercept = FALSE)\$coef
betaEstimate <- lsfit_coef
z1 <- c(xdexpect %*% betaEstimate)
w01 <- w1
w1[] <- asym
w1[ !(ydexpect > z1) ] <- 1 - asym
dw1 <- sum(w1 != w01, na.rm = TRUE)
it <- it + 1
if (it == 49)
warning("IWLS weights did not converge
after 50 iterations.")
}
res <- c(yvec - z1)
weight_res <- c(w1) * res
split_id <- split(id, id, drop = TRUE)
bdiag_sigma <- Matrix::bdiag(lapply(split_id, function(x)
{ tcrossprod(weight_res[which(id==x)]) } ) )
mat_bread <- solve(t(xmat) %*% diag(c(w1)) %*% xmat)
mat_meat <- t(xmat) %*% bdiag_sigma %*% xmat
denom_dof <- length(weight_res) - ncol(xmat) - length(split_id)
dof <- length(weight_res) / denom_dof
ecart_type <- sqrt(diag(as.matrix(cov_mat)))
out <- list(asymPoint = asym,
coefEst = betaEstimate,
standardError = ecart_type,
asymWeight = c(w1),
Residuals = res,
covMat = cov_mat)
out
}

## Try the erfe package in your browser

Any scripts or data that you put into this service are public.

erfe documentation built on Dec. 16, 2022, 5:15 p.m.