diff --git a/NAMESPACE b/NAMESPACE
index 7bcc23c..4916296 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -84,6 +84,7 @@ export(rlassologitEffects)
 export(tsls)
 importFrom(Formula,model.part)
 importFrom(methods,is)
+importFrom(stats,.lm.fit)
 importFrom(stats,as.formula)
 importFrom(stats,binomial)
 importFrom(stats,coef)
diff --git a/R/pkg-package.R b/R/pkg-package.R
index 2d998dd..48804ec 100644
--- a/R/pkg-package.R
+++ b/R/pkg-package.R
@@ -33,6 +33,7 @@
 #' @importFrom stats cor
 #' @importFrom stats glm
 #' @importFrom stats lm
+#' @importFrom stats .lm.fit
 #' @importFrom stats model.frame
 #' @importFrom stats model.matrix
 #' @importFrom stats model.response
diff --git a/R/rlasso.R b/R/rlasso.R
index 8e18fcd..32ee478 100644
--- a/R/rlasso.R
+++ b/R/rlasso.R
@@ -221,7 +221,9 @@ rlasso.default <- function(x, y, post = TRUE, intercept = TRUE, model = TRUE,
   XX <- crossprod(x)
   Xy <- crossprod(x, y)
   
-  startingval <- init_values(x,y)$residuals
+  init_vals <- init_values(x,y) 
+  
+  startingval <- init_vals$residuals
   pen <- lambdaCalculation(penalty = penalty, y = startingval, x = x)
   lambda <- pen$lambda
   Ups0 <- Ups1 <- pen$Ups0
@@ -235,13 +237,13 @@ rlasso.default <- function(x, y, post = TRUE, intercept = TRUE, model = TRUE,
     #coefTemp <- LassoShooting.fit(x, y, lambda, XX = XX, Xy = Xy)$coefficients
     #xn <- t(t(x)/as.vector(Ups1))
     if (mm==1 && post) {
-      coefTemp <- LassoShooting.fit(x, y, lambda/2, XX = XX, Xy = Xy)$coefficients
+      coefTemp <- LassoShooting.fit(x, y, lambda/2, XX = XX, Xy = Xy, beta.start = init_vals$coefficients)$coefficients
       #lasso.reg <- glmnet::glmnet(xn, y, family = c("gaussian"), alpha = 1,
       #                            lambda = lambda0/(2*n)/2, standardize = FALSE, intercept = FALSE)
       #lasso.reg <- glmnet::glmnet(x, y, family = c("gaussian"), alpha = 1,
       #                           lambda = lambda0/(2*n)/2, standardize = FALSE, intercept = FALSE,  penalty.factor = Ups1)
     } else {
-      coefTemp <- LassoShooting.fit(x, y, lambda, XX = XX, Xy = Xy)$coefficients
+      coefTemp <- LassoShooting.fit(x, y, lambda, XX = XX, Xy = Xy, beta.start = init_vals$coefficients)$coefficients
       #lasso.reg <- glmnet::glmnet(xn, y, family = c("gaussian"), alpha = 1,
       #                           lambda = lambda0/(2*n), standardize = FALSE, intercept = FALSE)
       #lasso.reg <- glmnet::glmnet(x, y, family = c("gaussian"), alpha = 1,
@@ -280,8 +282,10 @@ rlasso.default <- function(x, y, post = TRUE, intercept = TRUE, model = TRUE,
     
     # refinement variance estimation
     if (post) {
-      reg <- lm(y ~ -1 + x1)
-      coefT <- coef(reg)
+      #reg <- lm(y ~ -1 + x1)
+      #coefT <- coef(reg)
+      reg <- .lm.fit(x1, y)
+      coefT <- reg$coefficients
       coefT[is.na(coefT)] <- 0
       e1 <- y - x1 %*% coefT
       coefTemp[ind1] <- coefT
@@ -492,14 +496,22 @@ lambdaCalculation <- function(penalty = list(homoscedastic = FALSE, X.dependent.
     # }
     xehat <- x*ehat
     psi <- apply(xehat, 2, function(x) mean(x^2))
-    tXehattpsi <- t(t(xehat)/sqrt(psi))
+    #tXehattpsi <- t(t(xehat)/sqrt(psi))
+    #
+    #  for (l in 1:R) {
+    #    g <- matrix(rep(rnorm(n), each = p), ncol = p, byrow = TRUE)
+    #    #sim[l] <- n * max(2 * colMeans(x * ehat* g))
+    #    sim[l] <- n * max(2 * abs(colMeans(tXehattpsi * g)))
+    #  }
     
-      for (l in 1:R) {
-        g <- matrix(rep(rnorm(n), each = p), ncol = p, byrow = TRUE)
-        #sim[l] <- n * max(2 * colMeans(x * ehat* g))
-        sim[l] <- n * max(2 * abs(colMeans(tXehattpsi * g)))
-      }
+    Xehattpsi <- t(xehat)/sqrt(psi)
+    for (l in 1:R) {
+      g <- rnorm(n)
+      sim[l] <- max(2* abs(Xehattpsi %*% g))
+    }
     
+    #g <- matrix(rnorm(n*R), ncol = R, byrow=FALSE)
+    #sim <- apply(2 * abs(Xehattpsi %*% g), 2, max) 
     
     lambda0 <- penalty$c * quantile(sim, probs = 1 - penalty$gamma)
     Ups0 <- 1/sqrt(n) * sqrt(t(t(y^2) %*% (x^2)))