diff --git a/.DS_Store b/.DS_Store
index ca7b857..daf4c11 100644
Binary files a/.DS_Store and b/.DS_Store differ
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..a7c002d
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,15 @@
+language: r
+cache: packages
+sudo: true
+
+r:
+#  - oldrel
+  - release
+  - devel
+
+
+warnings_are_errors: true
+
+
+notifications:
+    slack: newstat:ZqBL05R4NoRMqPZXnuzbF4OZ
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index 4b36d29..5801194 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -3,7 +3,7 @@ Type: Package
 Title: Nonparametric Estimation of Regression Models with
     Factor-by-Curve Interactions
 Version: 1.4.0.9000
-Date: 2016-11-18
+Date: 2017-10-27
 Author: Marta Sestelo [aut, cre],
     Nora M. Villanueva [aut],
     Javier Roca-Pardinas [aut]
@@ -35,4 +35,4 @@ Imports:
     shinyjs,
     wesanderson,
     ggplot2
-RoxygenNote: 5.0.1
+RoxygenNote: 6.0.1
diff --git a/NAMESPACE b/NAMESPACE
index 8496894..d877168 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -37,7 +37,10 @@ importFrom(sfsmisc,D1D2)
 importFrom(shinyjs,colourInput)
 importFrom(shinyjs,useShinyjs)
 importFrom(stats,as.formula)
+importFrom(stats,coef)
 importFrom(stats,lm)
+importFrom(stats,model.response)
+importFrom(stats,model.weights)
 importFrom(stats,na.omit)
 importFrom(stats,predict)
 importFrom(stats,quantile)
diff --git a/NEWS b/NEWS
index 0801c74..1f08761 100644
--- a/NEWS
+++ b/NEWS
@@ -26,4 +26,18 @@ This file documents software changes since the previous edition.
 
    * the plotting functions have been changed. Users can choose between the plot.frfast function, used for base graphics, and the autoplot.frfast function, which is appropriate for ggplot2-type plot and returns objects of the ggplot class. 
    * we have merged the plottdiff function of the previous version of the package with plot.frfast and with  autoplot.frfast by means of the new argument diffwith. This new argument lets users visualize the differences between two factor’s levels.
+ 
+ 
+ 
+  # npregfast 1.4.1 (2017-10-30)
+
+  * the localtest function has a new argument: ci.level.
+  * new implementation in R of the allotest. Now it is possible
+  to choose between two statistic test: "res" (based on residuals on the 
+  null model) or "lrt" (based on the likelihood ratio test using rss0 and rss1).
+  * fixed some bugs in the formula argument of frfast function.
+  
+  
+
+  
  
\ No newline at end of file
diff --git a/R/allotest.R b/R/allotest.R
index 88ef767..367a06b 100644
--- a/R/allotest.R
+++ b/R/allotest.R
@@ -8,9 +8,26 @@
 #' @param na.action A function which indicates what should happen when the 
 #' data contain 'NA's. The default is 'na.omit'.
 #'@param nboot Number of bootstrap repeats.
-#'@param kbin Number of binning nodes over which the function
-#' is to be estimated.
 #'@param seed Seed to be used in the bootstrap procedure.
+#' @param cluster A logical value. If  \code{TRUE} (default), the
+#'  bootstrap procedure is  parallelized (only for \code{smooth = "splines"}).
+#'   Note that there are cases 
+#'  (e.g., a low number of bootstrap repetitions) that R will gain in
+#'  performance through serial computation. R takes time to distribute tasks
+#'  across the processors also it will need time for binding them all together
+#'  later on. Therefore, if the time for distributing and gathering pieces
+#'  together is greater than the time need for single-thread computing, it does
+#'  not worth parallelize.
+#'@param ncores An integer value specifying the number of cores to be used
+#' in the parallelized procedure. If \code{NULL} (default), the number of cores 
+#' to be used is equal to the number of cores of the machine - 1.
+#'@param test Statistic test to be used, based on residuals on the null model
+#'  (\code{res}) or based on the likelihood ratio test 
+#'  using rss0 and rss1 \code{lrt}.
+#' @param \ldots Other options.
+#'  
+#'  
+#' 
 #'@details In order to facilitate the choice of a model appropriate
 #' to the data while at the same time endeavouring to minimise the 
 #' loss of information,  a bootstrap-based procedure, that test whether the 
@@ -65,7 +82,18 @@
 
 
 allotest <- function(formula, data = data, na.action = "na.omit",
-                     nboot = 500, kbin = 200, seed = NULL) {
+                     nboot = 500, seed = NULL, cluster = TRUE,
+                     ncores = NULL, test = "res", ...) {
+  
+  if (isTRUE(cluster)) {
+    if (is.null(ncores)) {
+      num_cores <- detectCores() - 1
+    }else{
+      num_cores <- ncores
+    }
+    registerDoParallel(cores = num_cores)
+  }
+  
   
   ffr <- interpret.frfastformula(formula, method = "frfast")
   varnames <- ffr$II[2, ]
@@ -101,7 +129,7 @@ allotest <- function(formula, data = data, na.action = "na.omit",
     set.seed(seed)
   }
   
-  umatrix <- matrix(runif(n*nboot), ncol = nboot, nrow = n)
+  #umatrix <- matrix(runif(n*nboot), ncol = nboot, nrow = n)
   
   if (is.null(f)) 
     f <- rep(1, n)
@@ -114,23 +142,40 @@ allotest <- function(formula, data = data, na.action = "na.omit",
   for (i in etiquetas) {
     yy <- data[, 1][f == i]
     xx <- data[, 2][f == i]
+    yy[yy == 0] <- 0.0001
+    xx[xx == 0] <- 0.0001
     n <- length(xx)
     w <- rep(1, n)
-    fit <- .Fortran("allotest_", 
-                    x = as.double(xx), 
-                    y = as.double(yy), 
-                    w = as.double(w), 
-                    n = as.integer(n), 
-                    kbin = as.integer(kbin), 
-                    nboot = as.integer(nboot), 
-                   # seed = as.integer(seed),
-                    T = as.double(-1), 
-                    pvalue = as.double(-1),
-                    umatrix = as.double(umatrix),
-                    PACKAGE = "npregfast"
-                    )
     
-    res[[i]] <- list(statistic = fit$T, pvalue = fit$pvalue)
+    m <- lm(log(yy) ~ log(xx))
+    muhatg <- exp(coef(m)[1]) * xx**coef(m)[2]
+    errg <- yy - muhatg
+    
+   if(test == "res") {t <- sta_res(x = xx, y = yy)}
+   if(test == "lrt") {t <- sta_rss(x = xx, y = yy)}
+    #print(c(t1, t2))
+    
+    yboot <- replicate(nboot, muhatg + errg * 
+                         sample(c(-sqrt(5) + 1, sqrt(5) + 1)/2, size = n, 
+                                replace = TRUE, 
+                                prob = c(sqrt(5) + 1, sqrt(5) - 1)/(2 * sqrt(5))))
+    
+    if(test == "res") {
+      tboot <- unlist(foreach(i = 1:nboot) %dopar% {
+      sta_res(x = xx, y = yboot[, i])
+    })
+    }
+    
+    if(test == "lrt") {
+    tboot <- unlist(foreach(i = 1:nboot) %dopar% {
+      sta_rss(x = xx, y = yboot[, i])
+    })
+    }
+    pvalue <- mean(tboot>t)
+
+
+    res[[i]] <- list(statistic = c(t), pvalue = c(pvalue))
+    
   }
   
   
@@ -156,4 +201,40 @@ allotest <- function(formula, data = data, na.action = "na.omit",
   
   return(result)
   
-} 
\ No newline at end of file
+} 
+
+
+
+
+
+
+sta_res <- function(x, y){
+  y[y == 0] <- 0.0001
+  model <- lm(log(y) ~ log(x))
+  muhat <- exp(coef(model)[1]) * x**coef(model)[2]
+  residuo <- y - muhat
+  pred <- as.numeric(predict(gam(residuo ~ s(x))))
+  #pred <- predict(frfast(residuo ~ x, data = data.frame(x, residuo), nboot = 0), newdata = data.frame(x=x))$Estimation[,1]
+  #pred <- pred - mean(pred)
+  #rango <- max(x) - min(x)
+  #ii <- abs(x) <= (max(x)-0.1*rango)
+  t <- sum(abs(pred))
+}
+
+sta_rss <- function(x, y){
+  y[y == 0] <- 0.0001
+  model <- lm(log(y) ~ log(x))
+  m0 <- exp(coef(model)[1]) * x**coef(model)[2]
+  rss0 <- sum((y - m0)**2)
+  m1 <- as.numeric(predict(gam(y ~ s(x))))
+  rss1 <- sum((y - m1)**2)
+  #rango <- max(x) - min(x)
+  #ii <- abs(x) <= (max(x)-0.1*rango)
+  t <- (rss0 - rss1)/rss1
+}
+
+
+
+
+
+
diff --git a/R/frfast.R b/R/frfast.R
index 3ec5b3c..3503d03 100644
--- a/R/frfast.R
+++ b/R/frfast.R
@@ -219,13 +219,14 @@
 #' @importFrom doParallel registerDoParallel
 #' @importFrom parallel detectCores
 #' @importFrom foreach foreach %dopar%
+#' @importFrom stats coef model.response model.weights
 #' 
 #' @export
 
 
 
 
-frfast <- function(formula, data = data, na.action = "na.omit",
+frfast <- function(formula, data, na.action = "na.omit",
                    model = "np", smooth = "kernel", 
                    h0 = -1.0, h = -1.0, 
                    nh = 30, weights = NULL, kernel = "epanech", p = 3, 
@@ -239,9 +240,9 @@ frfast <- function(formula, data = data, na.action = "na.omit",
   if (missing(formula)) {
     stop("Argument \"formula\" is missing, with no default")
   }
-  if (missing(data)) {
-    stop("Argument \"data\" is missing, with no default")
-  }
+  #if (missing(data)) {
+  #  stop("Argument \"data\" is missing, with no default")
+  #}
   if (!(kernel %in% 1:3)) {
     stop("Kernel not suported")
   }
@@ -274,37 +275,89 @@ frfast <- function(formula, data = data, na.action = "na.omit",
   ncmax <- 5
   c2 <- NULL
   
+
   
   
-  
+ 
   if (smooth != "splines") {
     
-    ffr <- interpret.frfastformula(formula, method = "frfast")
-    varnames <- ffr$II[2, ]
-    aux <- unlist(strsplit(varnames,split = ":"))
+    cl <- match.call()
+    mf <- match.call(expand.dots = FALSE)
+    m <- match(x = c("formula", "data", "subset", "weights", "na.action", "offset"), 
+               table = names(mf), nomatch = 0L)
+    mf <- mf[c(1L, m)]
+    mf$drop.unused.levels <- TRUE
+    mf[[1L]] <- quote(stats::model.frame)
+    
+    
+    
+    
+    
+    mf <- eval(expr = mf, envir = parent.frame())
+    
+    mt <- attr(mf, "terms")
+    y <- model.response(mf, "numeric")
+    w <- as.vector(model.weights(mf))
+    if (!is.null(w) && !is.numeric(w)) 
+      stop("'weights' must be a numeric vector")
+    
+    terms <- attr(mt, "term.labels")
+    
+    
+    
+    aux <- unlist(strsplit(terms,split = ":"))
     varnames <- aux[1]
+    namef <- aux[2]
+    response <- as.character(attr(mt, "variables")[2])
     if (unlist(strsplit(varnames,split = ""))[1] == "s") {
       stop("Argument \"formula\" is wrong specified, see details of
-model specification in 'Details' of the frfast help." )
+           model specification in 'Details' of the frfast help." )
     }
-    namef <- aux[2]
-    if (length(aux) == 1) {f <- NULL}else{f <- data[ ,namef]}
-    newdata <- data
-    data <- data[ ,c(ffr$response, varnames)]
     
+
+    #newdata <- data
+    data <- mf
     
-    if (na.action == "na.omit"){ # ver la f
+    if (na.action == "na.omit"){ # ver la f, corregido
       data <- na.omit(data)
     }else{
       stop("The actual version of the package only supports 'na.omit' (observations are removed 
            if they contain any missing values)")
     }
-    #newdata <- na.omit(newdata[ ,varnames])
+    
+    if (length(aux) == 1) {f <- NULL}else{f <- data[ ,namef]}
     n <- nrow(data)
     
+
   }else{
+    
+    
     ffr <- interpret.gam(formula)
-    varnames <- ffr$pred.names[1]
+    cl <- match.call()
+    mf <- match.call(expand.dots = FALSE)
+    mf$formula <- ffr$fake.formula
+    
+    m <- match(x = c("formula", "data", "subset", "weights", "na.action", "offset"), 
+              table = names(mf), nomatch = 0L)
+   mf <- mf[c(1L, m)]
+    
+    mf$drop.unused.levels <- TRUE
+    mf[[1L]] <- quote(stats::model.frame)
+    
+    mf <- eval(expr = mf, envir = parent.frame())
+    
+    mt <- attr(mf, "terms")
+    y <- model.response(mf, "numeric")
+    w <- as.vector(model.weights(mf))
+    if (!is.null(w) && !is.numeric(w)) 
+      stop("'weights' must be a numeric vector")
+    
+    terms <- attr(mt, "term.labels")
+    
+    response <- as.character(attr(mt, "variables")[2])
+    
+
+    varnames <- terms[1]
     if (":" %in% unlist(strsplit(ffr$fake.names,split = ""))) {
       stop("Argument \"formula\" is wrong specified, see details of
            model specification in 'Details' of the frfast help." )
@@ -314,30 +367,72 @@ model specification in 'Details' of the frfast help." )
            model specification in 'Details' of the frfast help." )
     }
     
-    namef <- ffr$pred.names[2]
-    if (length(ffr$pred.names) == 1) {f <- NULL}else{f <- data[ ,namef]}
-    newdata <- data
     
-    if (length(ffr$pred.names) == 1) {
-      data <- data[ ,c(ffr$response, varnames)]
-    }else{
-      data <- data[ ,c(ffr$response, varnames, namef)]
-    }
+    
+    # newdata <- data
+    
+    # if (length(ffr$pred.names) == 1) {
+    #   data <- data[ ,c(ffr$response, varnames)]
+    # }else{
+    #   data <- data[ ,c(ffr$response, varnames, namef)]
+    # }
+    # 
+    datam <- mf
+    
     
     if (na.action == "na.omit"){
-      data <- na.omit(data)
+      datam <- na.omit(datam)
     }else{
       stop("The actual version of the package only supports 'na.omit' (observations are removed 
            if they contain any missing values)")
     }
     
-    n <- nrow(data)
+    
+    if (length(terms) == 1) {
+      f <- NULL
+      namef <- 1
+    }else{
+      namef <- terms[2]
+      f <- mf[ ,namef]
+    }
+    n <- nrow(datam)
   }
   
   
   
   
   
+  if(missing(data)) {
+    
+    response <- strsplit(response, "\\$")[[1]][2]
+    terms2 <- strsplit(terms, "\\$")
+    
+    if(length(terms) == 1){
+      formula <- as.formula(paste0(response, "~s(",terms2[[1]][2],")"))
+      #data <- data
+      names(datam) <- c(response, terms2[[1]][2])
+      varnames <- terms2[[1]][2]
+      namef <- "F"
+    }else{
+      formula <- as.formula(paste0(response, "~s(",terms2[[1]][2],", by = ", terms2[[2]][2],")"))
+      #data2 <- data
+      names(datam) <- c(response, terms2[[1]][2], terms2[[2]][2])
+      varnames <- terms2[[1]][2]
+      namef <- terms2[[2]][2]
+    }
+    
+    
+    
+  data <- datam
+
+  
+  
+  }
+
+  
+#   strsplit(formula,split = "$")
+  
+  
   if (is.null(f)) f <- rep(1.0, n)
   etiquetas <- unique(f)
   nf <- length(etiquetas)
@@ -356,7 +451,7 @@ model specification in 'Details' of the frfast help." )
     if (length(h) == 1) h <- rep(h, nf)
   }
   
-  
+  weights <- w
   if (is.null(weights)) {
     weights <- rep(1.0, n)
   }else{
@@ -393,7 +488,7 @@ model specification in 'Details' of the frfast help." )
     frfast  <- .Fortran("frfast_",
                         f = as.integer(f),
                         x = as.double(data[ ,varnames]),
-                        y = as.double(data[ ,ffr$response]),
+                        y = as.double(data[ ,response]),
                         w = as.double(weights),
                         n = as.integer(n),
                         h0 = as.double(h0),
@@ -511,7 +606,7 @@ model specification in 'Details' of the frfast help." )
                 b = frfast$b,
                 bl = frfast$binf,
                 bu = frfast$bsup,
-                name = c(ffr$response,varnames),
+                name = c(response,varnames),
                 formula = formula,
                 nh = frfast$nh,
                 r2 = r2,
@@ -527,8 +622,8 @@ model specification in 'Details' of the frfast help." )
       # grid
       xgrid <- seq(min(data[ ,varnames]), max(data[ ,varnames]), length.out = kbin)
       newd <- expand.grid(xgrid, unique(f))
-      names(newd) <- ffr$pred.names
-      
+      names(newd) <- c(varnames, namef)
+     
       # estimations
       p <- array(NA, dim = c(kbin, 3, nf))
       m <- gam(formula, weights = weights, data = data.frame(data, weights), ...)
@@ -549,7 +644,8 @@ model specification in 'Details' of the frfast help." )
                           FUN.VALUE = numeric(kfino))
       
       newdfino <- data.frame(as.vector(xgridfino), rep(unique(f), each = kfino))
-      names(newdfino) <- ffr$pred.names
+        names(newdfino) <- c(varnames, namef)
+     
       
       # max
       muhatfino <- as.vector(predict(m, newdata = newdfino, type = "response"))
@@ -619,7 +715,7 @@ model specification in 'Details' of the frfast help." )
     # bootstrap
     m <- gam(formula, weights = weights, data = data.frame(data, weights), ...)
     muhat <- as.vector(predict(m, type = "response"))
-    err <- data[, ffr$response] - muhat
+    err <- data[, response] - muhat
     err <- err - mean(err)
     yboot <- replicate(nboot, muhat + err *
                          sample(c(-sqrt(5) + 1, sqrt(5) + 1)/2, size = n,
@@ -630,7 +726,7 @@ model specification in 'Details' of the frfast help." )
     
     allboot <- foreach(i = 1:nboot) %dopar% {
       datab <- data
-      datab[, ffr$response] <- yboot[, i]
+      datab[, response] <- yboot[, i]
       aux <- mainfun(formula, data = data.frame(datab, weights), 
                      weights = weights, ...)
       return(aux)
@@ -729,7 +825,7 @@ model specification in 'Details' of the frfast help." )
                 h0 = NA,
                 fmod = f,
                 xdata = data[, varnames],
-                ydata = data[, ffr$response],
+                ydata = data[, response],
                 w = weights,
                 #fact=fact,  # Lo tuve que comentar pq me daba error
                 kbin = kbin,
@@ -753,7 +849,7 @@ model specification in 'Details' of the frfast help." )
                 b = NA,
                 bl = NA,
                 bu = NA,
-                name = c(ffr$response,varnames),
+                name = c(response,varnames),
                 formula = formula,
                 nh = nh,
                 r2 = NA,
diff --git a/R/globaltest.R b/R/globaltest.R
index bec8df3..a2fa1b4 100644
--- a/R/globaltest.R
+++ b/R/globaltest.R
@@ -300,19 +300,25 @@ globaltest <- function(formula, data = data, na.action = "na.omit",
                            kernel = as.integer(kernel), 
                            nboot = as.integer(nboot), 
                            r = as.integer(der), 
-                           T = as.double(rep(-1, 1)), 
-                           pvalor = as.double(rep(-1, 1)), 
+                           T = as.double(rep(-1, 4)), 
+                           pvalor = as.double(rep(-1, 4)), 
                            #  umatrix = as.double(umatrix)
                            #seed = as.integer(seed),
                            umatrix = array(umatrix, c(n, nboot)),
                            PACKAGE = "npregfast"
     )
     
-    if (globaltest$pvalor < 0.05) {
-      decision <- "Rejected"
+    decision <- character(4)
+    for (j in 1:4){
+    if (globaltest$pvalor[j] < 0.05) {
+      decision[j] <- "Rejected"
     } else {
-      decision <- "Acepted"
+      decision[j] <- "Accepted"
+    }
     }
+    
+    
+    
     res <- data.frame(cbind(Statistic = globaltest$T, pvalue = globaltest$pvalor), 
                       Decision = I(decision))
     # res=cbind(Statistic=round(globaltest$T,digits=4),pvalue=round(globaltest$pvalor,digits=4),Decision=I(decision))
diff --git a/R/interpret.frfastformula.R b/R/interpret.frfastformula.R
index 1763433..4f751ea 100755
--- a/R/interpret.frfastformula.R
+++ b/R/interpret.frfastformula.R
@@ -1,54 +1,54 @@
 #' @importFrom stats as.formula terms.formula 
 
 interpret.frfastformula <-
-function(formula, method = "frfast") {
-
+  function(formula, method = "frfast") {
+    
     env <- environment(formula) 
     if(inherits(formula, "character"))          
-        formula <- as.formula(formula)
+      formula <- as.formula(formula)
     tf <- terms.formula(formula)
     terms <- attr(tf, "term.labels")
-   # if(length(grep(":",terms))!=0)  stop("Symbol '*' is not allowed")
+    # if(length(grep(":",terms))!=0)  stop("Symbol '*' is not allowed")
     
     nt <- length(terms)
-   if(attr(tf, "response") > 0) {     
-    	ns <- attr(tf, "specials")$frfast - 1 # -1 for the response
-    	response <- as.character(attr(tf, "variables")[2])
-    	vtab<-attr(tf,"factors")
+    if(attr(tf, "response") > 0) {     
+      ns <- attr(tf, "specials")$frfast - 1 # -1 for the response
+      response <- as.character(attr(tf, "variables")[2])
+      vtab<-attr(tf,"factors")
     } else {
-    	ns <- attr(tf, "specials")$frfast
-    	response <- NULL
+      ns <- attr(tf, "specials")$frfast
+      response <- NULL
     }
     
     II <- list()
     k <- 0
     if(nt) {
-        for (i in 1:nt) {
-            if (i %in% ns) {
-                k = k+1                   
-                st <- eval(parse(text = terms[i]), envir = env)
-                if(method ==  "frfast" & st$cov[1] != "ONE") {
-                	stop("The function frfast does not allow for \"by\" variables")
-                }
-                II[[k]] <- st$cov
-               # h[[k]] <- st$h
-               # partial[k] <- terms[i]                                                          
-            } else {
-                k = k+1
-                II[[k]]<- c("ONE", terms[i])
-               # h[[k]] <- 0
-               # partial[k] <- terms[i]
-            }
-        }           
+      for (i in 1:nt) {
+        if (i %in% ns) {
+          k = k+1                   
+          st <- eval(parse(text = terms[i]), envir = env)
+          if(method ==  "frfast" & st$cov[1] != "ONE") {
+            stop("The function frfast does not allow for \"by\" variables")
+          }
+          II[[k]] <- st$cov
+          # h[[k]] <- st$h
+          # partial[k] <- terms[i]                                                          
+        } else {
+          k = k+1
+          II[[k]]<- c("ONE", terms[i])
+          # h[[k]] <- 0
+          # partial[k] <- terms[i]
+        }
+      }           
     }       
     II <- if(length(II)) {
-        matrix(unlist(II), nrow = 2)
+      matrix(unlist(II), nrow = 2)
     } else {
-        matrix(0, nrow = 2)
+      matrix(0, nrow = 2)
     }       
     res <- list(response = response, II = II,vtab=vtab)
     
     class(res) <- "frfast.formula"
     return(res)
-}
+  }
 
diff --git a/R/localtest.R b/R/localtest.R
index c8ab255..3e9aaad 100644
--- a/R/localtest.R
+++ b/R/localtest.R
@@ -55,6 +55,7 @@
 #'@param ncores An integer value specifying the number of cores to be used
 #' in the parallelized procedure. If \code{NULL} (default), the number of cores 
 #' to be used is equal to the number of cores of the machine - 1.
+#'@param ci.level Level of bootstrap confidence interval. Defaults to 0.95 (corresponding to 95\%). Note that the function accepts a vector of levels.
 #' @param \ldots Other options.
 #' 
 #' 
@@ -141,7 +142,7 @@ localtest <- function(formula, data = data, na.action = "na.omit",
                       nboot = 500, h0 = -1.0, h = -1.0, nh = 30, 
                       kernel = "epanech", p = 3, kbin = 100, rankl = NULL, 
                       ranku = NULL, seed = NULL, cluster = TRUE, 
-                      ncores = NULL, ...) {
+                      ncores = NULL, ci.level = 0.95, ...) {
   
   if(kernel == "gaussian")  kernel <- 3
   if(kernel == "epanech")   kernel <- 1
@@ -187,7 +188,7 @@ localtest <- function(formula, data = data, na.action = "na.omit",
   ncmax <- 5
   c2 <- NULL
   # if(is.null(seed)) seed <- -1
-  
+  nalfas <- length(ci.level)
   
   
   
@@ -323,22 +324,30 @@ localtest <- function(formula, data = data, na.action = "na.omit",
                           pcmin = as.double(rankl), # rango de busqueda maximo
                           r = as.integer(der),
                           D = as.double(rep(-1.0,1)),
-                          Ci = as.double(rep(-1.0,1)),
-                          Cs = as.double(rep(-1.0,1)),
+                          Ci = as.double(rep(-1.0,nalfas)),
+                          Cs = as.double(rep(-1.0,nalfas)),
                          # seed = as.integer(seed),
                           umatrix = as.double(umatrix),
+                         level = as.double(ci.level),
+                         nalfas = as.integer(nalfas),
                           PACKAGE = "npregfast"
     )
     
+    decision <- character(nalfas)
+    for(i in 1:nalfas){
+      if (localtest$Ci[i] <= 0 & 0 <= localtest$Cs[i]) {
+        decision[i] <- "Accepted"
+      } else {
+        decision[i] <- "Rejected"
+      }
+    }
+    
+    
     
     
-    if (localtest$Ci <= 0 & 0 <= localtest$Cs) {
-      decision <- "Acepted"
-    } else {
-      decision <- "Rejected"
-    }
     res <- cbind(d = round(localtest$D, digits = 4), Lwr = round(localtest$Ci, digits = 4), 
-                 Upr = round(localtest$Cs, digits = 4), Decision = decision)
+                 Upr = round(localtest$Cs, digits = 4), Decision = decision,
+                 Ci.Level = round(ci.level, digits = 2))
     # class(res) <- 'localtest'
     
   }else{
@@ -434,15 +443,30 @@ localtest <- function(formula, data = data, na.action = "na.omit",
     }
     
     
-    ci <- quantile(unlist(d_allboot), probs = c(0.025, 0.975), na.rm = TRUE)
+    
+    
+    decision <- character(nalfas)
+    cilower <- numeric(nalfas)
+    ciupper <- numeric(nalfas)
+    
+    for(i in 1:nalfas){
+      alpha <- 1-ci.level[i]
+      ci <- quantile(unlist(d_allboot), 
+                     probs = c(alpha/2, 1 - (alpha/2)), na.rm = TRUE)
     
     if (ci[1] <= 0 & 0 <= ci[2]) {
-      decision <- "Acepted"
+      decision[i] <- "Accepted"
     } else {
-      decision <- "Rejected"
+      decision[i] <- "Rejected"
     }
-    res <- cbind(d = round(d, digits = 4), Lwr = round(ci[1], digits = 4), 
-                 Upr = round(ci[2], digits = 4), Decision = decision)
+    cilower[i] <- ci[1]
+    ciupper[i] <- ci[2] 
+    }
+    
+    
+    res <- cbind(d = round(d, digits = 4), Lwr = round(cilower, digits = 4), 
+                 Upr = round(ciupper, digits = 4), Decision = decision,
+                 Ci.Level = round(ci.level, digits = 2))
     
     rownames(res) <- NULL
     
diff --git a/R/print.frfast.R b/R/print.frfast.R
index 1d364d8..19738b8 100644
--- a/R/print.frfast.R
+++ b/R/print.frfast.R
@@ -56,6 +56,7 @@ print.frfast <- function(x = model, ...) {
         cat("Triangular")
       if (model$kernel == 3) 
         cat("Gaussian")
+      cat("\n")
     }
     
     
diff --git a/inst/.DS_Store b/inst/.DS_Store
new file mode 100644
index 0000000..5b6ff9c
Binary files /dev/null and b/inst/.DS_Store differ
diff --git a/man/allotest.Rd b/man/allotest.Rd
index a970733..36ff4fc 100644
--- a/man/allotest.Rd
+++ b/man/allotest.Rd
@@ -5,7 +5,7 @@
 \title{Bootstrap based test for testing an allometric model}
 \usage{
 allotest(formula, data = data, na.action = "na.omit", nboot = 500,
-  kbin = 200, seed = NULL)
+  seed = NULL, cluster = TRUE, ncores = NULL, test = "res", ...)
 }
 \arguments{
 \item{formula}{An object of class \code{formula}: a sympbolic description
@@ -19,10 +19,27 @@ data contain 'NA's. The default is 'na.omit'.}
 
 \item{nboot}{Number of bootstrap repeats.}
 
-\item{kbin}{Number of binning nodes over which the function
-is to be estimated.}
-
 \item{seed}{Seed to be used in the bootstrap procedure.}
+
+\item{cluster}{A logical value. If  \code{TRUE} (default), the
+bootstrap procedure is  parallelized (only for \code{smooth = "splines"}).
+ Note that there are cases 
+(e.g., a low number of bootstrap repetitions) that R will gain in
+performance through serial computation. R takes time to distribute tasks
+across the processors also it will need time for binding them all together
+later on. Therefore, if the time for distributing and gathering pieces
+together is greater than the time need for single-thread computing, it does
+not worth parallelize.}
+
+\item{ncores}{An integer value specifying the number of cores to be used
+in the parallelized procedure. If \code{NULL} (default), the number of cores 
+to be used is equal to the number of cores of the machine - 1.}
+
+\item{test}{Statistic test to be used, based on residuals on the null model
+(\code{res}) or based on the likelihood ratio test 
+using rss0 and rss1 \code{lrt}.}
+
+\item{\ldots}{Other options.}
 }
 \value{
 An object is returned with the following elements:
@@ -59,9 +76,6 @@ library(npregfast)
 data(barnacle)
 allotest(DW ~ RC, data = barnacle, nboot = 50, seed = 130853)
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Sestelo, M. and Roca-Pardinas, J. (2011). A new approach to estimation of 
@@ -74,4 +88,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/autoplot.frfast.Rd b/man/autoplot.frfast.Rd
index f3ac9a0..225e049 100644
--- a/man/autoplot.frfast.Rd
+++ b/man/autoplot.frfast.Rd
@@ -132,4 +132,3 @@ gridExtra::grid.arrange(grobs = facs, ncol = 2, nrow = 1)
 \author{
 Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
-
diff --git a/man/barnacle.Rd b/man/barnacle.Rd
index 8c3fa8d..a800bd7 100644
--- a/man/barnacle.Rd
+++ b/man/barnacle.Rd
@@ -31,4 +31,3 @@ length-weight relationship of \eqn{Pollicipes}  \eqn{pollicipes}
 aspects of its biology and management. Journal of Shellfish Research, 
 30(3), 939--948.
 }
-
diff --git a/man/children.Rd b/man/children.Rd
index c17ae6e..3190bc9 100644
--- a/man/children.Rd
+++ b/man/children.Rd
@@ -28,4 +28,3 @@ data(children)
 head(children)
 
 }
-
diff --git a/man/critical.Rd b/man/critical.Rd
index c0d9f3e..999075c 100644
--- a/man/critical.Rd
+++ b/man/critical.Rd
@@ -51,9 +51,6 @@ critical(fit, der = 2)
 # critical(fit2, der = 1)
 # critical(fit2, der = 2)
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Sestelo, M. (2013). Development and computational implementation of 
@@ -61,4 +58,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/criticaldiff.Rd b/man/criticaldiff.Rd
index 928a4b6..f276245 100644
--- a/man/criticaldiff.Rd
+++ b/man/criticaldiff.Rd
@@ -48,9 +48,6 @@ criticaldiff(fit2)
 criticaldiff(fit2, der = 1)
 criticaldiff(fit2, der = 1, level1 = "lens", level2 = "barca")
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Sestelo, M. (2013). Development and computational implementation of 
@@ -58,4 +55,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/frfast.Rd b/man/frfast.Rd
index 06ee8c8..285aeec 100644
--- a/man/frfast.Rd
+++ b/man/frfast.Rd
@@ -4,7 +4,7 @@
 \alias{frfast}
 \title{Fitting nonparametric models}
 \usage{
-frfast(formula, data = data, na.action = "na.omit", model = "np",
+frfast(formula, data, na.action = "na.omit", model = "np",
   smooth = "kernel", h0 = -1, h = -1, nh = 30, weights = NULL,
   kernel = "epanech", p = 3, kbin = 100, nboot = 500, rankl = NULL,
   ranku = NULL, seed = NULL, cluster = TRUE, ncores = NULL, ...)
@@ -226,9 +226,6 @@ summary(fit3)
 # fit4 <- frfast(DW ~ RC : F, data = barnacle, model = "allo", nboot = 100)
 # summary(fit4)
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Huxley, J. S. (1924). Constant differential growth-ratios and their 
@@ -239,4 +236,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/globaltest.Rd b/man/globaltest.Rd
index 0b78fab..ee294c4 100644
--- a/man/globaltest.Rd
+++ b/man/globaltest.Rd
@@ -127,9 +127,6 @@ globaltest(DW ~ RC : F, data = barnacle, der = 1, seed = 130853, nboot = 100)
 # seed = 130853, der = 0, smooth = "splines")
 
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Sestelo, M. (2013). Development and computational implementation of 
@@ -137,4 +134,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/localtest.Rd b/man/localtest.Rd
index 6a5560a..22c7df0 100644
--- a/man/localtest.Rd
+++ b/man/localtest.Rd
@@ -7,7 +7,8 @@
 localtest(formula, data = data, na.action = "na.omit", der,
   smooth = "kernel", weights = NULL, nboot = 500, h0 = -1, h = -1,
   nh = 30, kernel = "epanech", p = 3, kbin = 100, rankl = NULL,
-  ranku = NULL, seed = NULL, cluster = TRUE, ncores = NULL, ...)
+  ranku = NULL, seed = NULL, cluster = TRUE, ncores = NULL,
+  ci.level = 0.95, ...)
 }
 \arguments{
 \item{formula}{An object of class \code{formula}: a sympbolic 
@@ -82,6 +83,8 @@ not worth parallelize.}
 in the parallelized procedure. If \code{NULL} (default), the number of cores 
 to be used is equal to the number of cores of the machine - 1.}
 
+\item{ci.level}{Level of bootstrap confidence interval. Defaults to 0.95 (corresponding to 95\%). Note that the function accepts a vector of levels.}
+
 \item{\ldots}{Other options.}
 }
 \value{
@@ -148,9 +151,6 @@ localtest(DW ~ RC : F, data = barnacle, der = 1, seed = 130853, nboot = 100)
 # localtest(height ~ s(age, by = sex), data = children, seed = 130853, 
 # der = 1, smooth = "splines") 
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Sestelo, M. (2013). Development and computational implementation of 
@@ -158,4 +158,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/npregfast.Rd b/man/npregfast.Rd
index f284d3c..aa13f89 100644
--- a/man/npregfast.Rd
+++ b/man/npregfast.Rd
@@ -60,9 +60,6 @@ For a listing of all routines in the NPRegfast package type:
 View a \href{http://sestelo.shinyapps.io/npregfast}{demo Shiny app}
 or see the full \href{https://github.com/sestelo/npregfast}{README} on GitHub.
 }
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
-}
 \references{
 Efron, B. (1979). Bootstrap methods: another look at the jackknife. 
 Annals of Statistics, 7, 1--26.
@@ -86,4 +83,6 @@ aspects of its biology and management. Journal of Shellfish Research,
 
 Wand, M. P. and Jones, M. C. (1995). Kernel Smoothing. Chapman & Hall, London.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/man/plot.frfast.Rd b/man/plot.frfast.Rd
index 07b4e72..30f1e9a 100644
--- a/man/plot.frfast.Rd
+++ b/man/plot.frfast.Rd
@@ -115,4 +115,3 @@ plot(fit2, fac = "barca", diffwith = "lens", der = 1)
 \author{
 Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
-
diff --git a/man/predict.frfast.Rd b/man/predict.frfast.Rd
index 993b1a5..b14e091 100644
--- a/man/predict.frfast.Rd
+++ b/man/predict.frfast.Rd
@@ -55,4 +55,3 @@ predict(fit, newdata = nd)
 \author{
 Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
-
diff --git a/man/reexports.Rd b/man/reexports.Rd
index 8fc910e..3a7999e 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -2,9 +2,10 @@
 % Please edit documentation in R/autoplot.frfast.R
 \docType{import}
 \name{reexports}
-\alias{autoplot}
 \alias{reexports}
+\alias{autoplot}
 \title{Objects exported from other packages}
+\keyword{internal}
 \description{
 These objects are imported from other packages. Follow the links
 below to see their documentation.
@@ -12,5 +13,4 @@ below to see their documentation.
 \describe{
   \item{ggplot2}{\code{\link[ggplot2]{autoplot}}}
 }}
-\keyword{internal}
 
diff --git a/man/runExample.Rd b/man/runExample.Rd
index 2e5c6b0..9b6a602 100644
--- a/man/runExample.Rd
+++ b/man/runExample.Rd
@@ -20,4 +20,3 @@ if (interactive()) {
   runExample()
 }
 }
-
diff --git a/man/summary.frfast.Rd b/man/summary.frfast.Rd
index 0a9ab44..f309e64 100644
--- a/man/summary.frfast.Rd
+++ b/man/summary.frfast.Rd
@@ -1,8 +1,8 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/summary.frfast.R
 \name{summary.frfast}
-\alias{print.frfast}
 \alias{summary.frfast}
+\alias{print.frfast}
 \title{Summarizing fits of \code{frfast} class}
 \usage{
 \method{summary}{frfast}(object = model, ...)
@@ -51,9 +51,6 @@ fit3 <- frfast(DW ~ RC, data = barnacle, model = "allo", nboot = 100)
 fit3
 summary(fit3)
 
-}
-\author{
-Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
 }
 \references{
 Sestelo, M. (2013). Development and computational implementation of 
@@ -61,4 +58,6 @@ estimation and inference methods in flexible regression models.
 Applications in Biology, Engineering and Environment. PhD Thesis, Department
 of Statistics and O.R. University of Vigo.
 }
-
+\author{
+Marta Sestelo, Nora M. Villanueva and Javier Roca-Pardinas.
+}
diff --git a/src/.gitignore b/src/.gitignore
index ff122c5..f52894f 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -1,2 +1,3 @@
 NPRegfast.so
 program2.o
+lsq.mod
diff --git a/src/program2.f90 b/src/program2.f90
index b51c8e5..bf6f8c3 100644
--- a/src/program2.f90
+++ b/src/program2.f90
@@ -17,10 +17,10 @@ subroutine allotest_(X,Y,W,n,kbin,nboot,T,pvalor,umatrix)
 !!DEC$ ATTRIBUTES DLLEXPORT::test_allo
 !!DEC$ ATTRIBUTES C, REFERENCE, ALIAS:'test_allo_' :: test_allo
 
-integer n,kbin,p,iboot,nboot,i
+integer n,kbin,p,iboot,nboot,i,j
 double precision X(n),X2(n),Y(n),Y2(n),W(n),&
 errg(n),muhatg(n),Yboot(n),h,T,Tboot,pvalor,&
-umatrix(n,nboot), aux
+umatrix(n,nboot), aux, beta(10)
 !real u, rand
 double precision u
 real,external::rnnof
@@ -28,7 +28,7 @@ subroutine allotest_(X,Y,W,n,kbin,nboot,T,pvalor,umatrix)
 
 
 
-
+w=1
 
 
 h=-1.0
@@ -44,15 +44,27 @@ subroutine allotest_(X,Y,W,n,kbin,nboot,T,pvalor,umatrix)
 
 ! Estimación Piloto 
 
+!p=1
+!call Reglineal_pred(X,Y,W,n,p,muhatg)
+
+
 p=1
-call Reglineal_pred(X,Y,W,n,p,muhatg)
+call Reglineal (X,Y,W,n,p,Beta)
 do i=1,n
- errg(i)=Y(i)-muhatg(i)
+muhatg(i)=beta(1)
+do j=1,p
+muhatg(i)=muhatg(i)+beta(j+1)*X(i)**j
+end do
 end do
 
 
-call RfastC3(X,Y,W,n,p,kbin,h,T)
+ errg=Y-muhatg
+
+
+!print *, errg(1:n)
 
+call RfastC3(X,Y,W,n,p,kbin,h,T)
+!print *, T 
 
 pvalor=0
 do iboot=1,nboot
@@ -79,35 +91,232 @@ subroutine allotest_(X,Y,W,n,kbin,nboot,T,pvalor,umatrix)
 
 
 
+subroutine allotest_sestelo_(X,Y,W,n,kbin,nboot,T,pvalor,umatrix)
+implicit none
+
+!!DEC$ ATTRIBUTES DLLEXPORT::test_allo
+!!DEC$ ATTRIBUTES C, REFERENCE, ALIAS:'test_allo_' :: test_allo
+
+integer n,kbin,p,iboot,nboot,i,j
+double precision X(n),X2(n),Y(n),Y2(n),W(n),&
+errg(n),muhatg(n),Yboot(n),h,T,Tboot,pvalor,&
+umatrix(n,nboot), aux, beta(10)
+!real u, rand
+double precision u
+real,external::rnnof
+integer,external::which_min,which_max2
+
+
+h=-1.0
+aux = 0.001
+do i=1,n
+ X2(i)=max(X(i),aux)
+ Y2(i)=max(Y(i),aux)
+end do
+
+X2=log(X2)
+Y2=log(Y2)
+
+! Estimación Piloto escala normal
+p=1
+call Reglineal (X2,Y2,W,n,p,Beta)
+
+do i=1,n
+muhatg(i)=exp(beta(1))
+do j=1,p
+muhatg(i)=muhatg(i)*X(i)**beta(j+1)
+end do
+end do
+errg=Y-muhatg ! residuos modelo alometrico
+
+
+!print *, errg(1:n)
+
+call RfastC3_sestelo(X,Y,W,n,p,kbin,h,T)
+!print *, T 
+
+pvalor=0
+do iboot=1,nboot
+ do i=1,n
+  !u=RAND()
+  !call test_random(u)
+  u = umatrix(i,iboot)
+  if (u.le.(5.0+sqrt(5.0))/10) then
+   Yboot(i)=muhatg(i)+errg(i)*(1-sqrt(5.0))/2
+  else
+   Yboot(i)=muhatg(i)+errg(i)*(1+sqrt(5.0))/2
+  end if
+ end do
+ h=-1.0
+call RfastC3_sestelo(X,Yboot,W,n,p,kbin,h,Tboot)
+if(Tboot.gt.T) pvalor=pvalor+1
+end do
+
+pvalor=pvalor/nboot
+
+end subroutine
+
+
+subroutine RfastC3_sestelo(X,Y,W,n,p,kbin,h,T)
+implicit none
+
+integer,parameter::kernel=1,nh=20
+integer n,kbin,p,i,j
+double precision X(n),Y(n),W(n),Xb(kbin),pred1(n),h,X2(n),Y2(n),&
+Pb(kbin,3),residuo(n),predg(n),T,sumw,sum2,xmin,xmax,rango,beta(10),aux
+integer,external::which_min
+
+
+
+aux = 0.001
+do i=1,n
+ X2(i)=max(X(i),aux)
+ Y2(i)=max(Y(i),aux)
+end do
+! Ajustamos el modelo lineal primero
+X2=log(X2)
+Y2=log(Y2)
+
+p=1
+call Reglineal (X2,Y2,W,n,p,Beta)
+do i=1,n
+predg(i)=exp(beta(1))
+do j=1,p
+predg(i)=predg(i)*X(i)**beta(j+1)
+end do
+end do
+
+Residuo=Y-predg
+
+! -----------------------------------------
+
+!print *, predg
+
+!do i=1,n
+ ! print (*,*) predg(i)
+!end do
+
+p=3
+!call rfast_h_alo(X,Residuo,W,n,h,p,Xb,Pb,kbin,kernel,nh)
+call Grid1d(X,W,n,Xb,kbin)
+
+
+
+call rfast_h(X,Residuo,W,n,h,p,Xb,Pb(1,1),kbin,kernel,nh)
+
+
+!stop
+
+!call Interpola_alo(Xb,Pb,kbin,X,pred1,pred2,n)
+call Interpola(Xb,Pb(1,1),kbin,X,pred1,n)
+
+!do i=1,n
+!print *, residuo(1:n)
+! print *, pred1(i)
+!end do
+
+
+ !print *, pred1(1:n)
+
+!Centro las pred1
+sumw=0
+sum2=0
+do i=1,n
+sumw=sumw+W(i)
+sum2=sum2+pred1(i)
+end do
+
+
+
+do i=1,n
+Pred1(i)=pred1(i)-(sum2/sumw)
+end do
+
+
+
+
+xmin=9999
+xmax=-xmin
+do i=1,n
+if(x(i).le.xmin) xmin=x(i)
+if(x(i).ge.xmax) xmax=x(i)
+end do
+
+rango=xmax-xmin
+
+T=0
+do i=1,n
+!if (abs(X(i)).le.xmax-(0.10*rango)) 
+T=T+abs(pred1(i))
+end do
+
+
+
+end subroutine
+
+
 
 subroutine RfastC3(X,Y,W,n,p,kbin,h,T)
 implicit none
 
 integer,parameter::kernel=1,nh=20
-integer n,kbin,p,i
-double precision X(n),Y(n),W(n),Xb(kbin),pred1(n),pred2(n),h,&
-Pb(kbin),residuo(n),predg(n),T,sumw,sum2,xmin,xmax,rango
+integer n,kbin,p,i,j
+double precision X(n),Y(n),W(n),Xb(kbin),pred1(n),h,&
+Pb(kbin,3),residuo(n),predg(n),T,sumw,sum2,xmin,xmax,rango,beta(10)
 integer,external::which_min
 
 
+
+
 ! Ajustamos el modelo lineal primero
 
 p=1
-call Reglineal_pred(X,Y,W,n,p,predg)
+!call Reglineal_pred(X,Y,W,n,p,predg)
+
+call Reglineal (X,Y,W,n,p,Beta)
 do i=1,n
-Residuo(i)=Y(i)-predg(i)
+predg(i)=beta(1)
+do j=1,p
+predg(i)=predg(i)+beta(j+1)*X(i)**j
+end do
 end do
+
+Residuo=Y-predg
+
 ! -----------------------------------------
 
-p=2 
+!print *, predg
+
+!do i=1,n
+ ! print (*,*) predg(i)
+!end do
+
+p=2
+
+
+
+!call rfast_h_alo(X,Residuo,W,n,h,p,Xb,Pb,kbin,kernel,nh)
 
+call Grid1d(X,W,n,Xb,kbin)
 
-call rfast_h_alo(X,Residuo,W,n,h,p,Xb,Pb,kbin,kernel,nh)
 
-call Interpola_alo(Xb,Pb,kbin,X,pred1,pred2,n)
 
+call rfast_h(X,Residuo,W,n,h,p,Xb,Pb(1,1),kbin,kernel,nh)
 
 
+!stop
+
+!call Interpola_alo(Xb,Pb,kbin,X,pred1,pred2,n)
+call Interpola(Xb,Pb(1,1),kbin,X,pred1,n)
+
+do i=1,n
+!print *, residuo(1:n)
+ !print *, pred1(i)
+end do
+
+
+ !print *, pred1(1:n)
+
 !Centro las pred1
 sumw=0
 sum2=0
@@ -116,10 +325,15 @@ subroutine RfastC3(X,Y,W,n,p,kbin,h,T)
 sum2=sum2+pred1(i)
 end do
 
+
+
 do i=1,n
-Pred1(i)=pred1(i)-sum2/sumw
+Pred1(i)=pred1(i)-(sum2/sumw)
 end do
 
+
+
+
 xmin=9999
 xmax=-xmin
 do i=1,n
@@ -170,6 +384,8 @@ subroutine Reglineal_pred(X,Y,W,n,p,Pred)
 
 
 
+
+
 !*********************************************************
 !		
 !Subroutine RFAST_H_alo MODIFICADA PARA EL CONTRASTE ALOMETRICO
@@ -271,7 +487,7 @@ subroutine rfast_h_alo(X,Y,W,n,h,p,Xb,Pb,kbin,kernel,nh)
 
 
 subroutine localtest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,&
-pcmax,pcmin,r,D,Ci,Cs,umatrix)
+pcmax,pcmin,r,D,Ci,Cs,umatrix,level,nalfas)
 
 
 !!DEC$ ATTRIBUTES DLLEXPORT::localtest
@@ -280,12 +496,12 @@ subroutine localtest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,&
 implicit none
 integer,parameter::kfino=1000
 integer i,n,j,kbin,p,nf,F(n),fact(nf),iboot,ir,l,k,&
-nh,nboot,kernel,r,index,posmin,posmax
+nh,nboot,kernel,r,index,posmin,posmax,nalfas
 double precision X(n),Y(n),W(n),Waux(n),xb(kbin),pb(kbin,3,nf),&
 u,h(nf),Pb_0(kbin,3),res(n),Pb_0boot(kbin,3,nboot),meanerr,P_0(n),Err(n),&
-C(3,nf),xmin(nf),xmax(nf),pcmax(nf),pcmin(nf),Ci,Cs,&
+C(3,nf),xmin(nf),xmax(nf),pcmax(nf),pcmin(nf),Ci(nalfas),Cs(nalfas),&
 Dboot(nboot),D,pmax,pasox,pasoxfino,icont(kbin,3,nf),xminc,xmaxc,h0,&
-umatrix(n,nboot)
+umatrix(n,nboot),level(nalfas)
 !REAL(4) rand 
 double precision, allocatable:: Yboot(:),muhatg(:),errg(:),errgboot(:),&
 muhatgboot(:),Xfino(:),Pfino(:),p0(:,:),pred(:),pboot(:,:,:,:),cboot(:,:,:),&
@@ -709,16 +925,39 @@ subroutine localtest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,&
 Ci=-1
 Cs=-1
 
+do i=1,nalfas
+  call ICbootstrap_beta_per(Dboot,nboot,1-level(i),Ci(i),Cs(i)) 
+end do
+
+end subroutine
+
 
-call ICbootstrap(D,Dboot,nboot,Ci,Cs) 
 
 
+
+subroutine ICbootstrap_beta_per(X,nboot,beta,li,ls)
+implicit none
+integer nboot,nalfa
+double precision X(nboot),li,ls,alfa(3),Q(3),beta
+
+alfa(1)=beta/2
+alfa(2)=0.5
+alfa(3)=1-beta/2
+nalfa=3
+call quantile (X,nboot,alfa,nalfa,Q)
+
+li=Q(1)!-Q(2)
+ls=Q(3)!-Q(2)
+
 end subroutine
 
 
 
 
 
+
+
+
 !*************************************************
 !*************************************************
 
@@ -732,18 +971,22 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
 
 implicit none
 integer i,z,n,j,kbin,p,nf,F(n),fact(nf),iboot,k,&
-nh,nboot,kernel,r,pp
+nh,nboot,kernel,r,pp,icont,ii
 double precision X(n),Y(n),W(n),Waux(n),xb(kbin),pb(kbin,3,nf),&
 h(nf),h0,hp(nf),pred1(kbin,nf),pred0(kbin),pol(n,nf),&
-u,Tboot,T,pvalor,umatrix(n,nboot)
+u,Tboot(4),pvalor(4),umatrix(n,nboot),h0i,hi(nf),hgi(nf),meanerr,T(4),&
+RSS0,RSS1,hg(nf)
 !REAL(4) rand 
 double precision, allocatable:: Yboot(:),muhatg(:),errg(:),errgboot(:),&
-muhatgboot(:),muhatg2(:)
+muhatgboot(:),muhatg2(:),fpar(:,:),fpar_est(:),Xaux(:)
 
 
-allocate (errg(n),muhatg(n),Yboot(n),errgboot(n),muhatgboot(n),muhatg2(n))
-
+allocate (errg(n),muhatg(n),Yboot(n),errgboot(n),muhatgboot(n),muhatg2(n),&
+  fpar_est(n),fpar(n,nf),Xaux(n))
 
+h0i = h0
+hi = h
+hgi = h0
 
 Xb=-1
 Pb=-1
@@ -754,6 +997,9 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
 call rfast_h(X,Y,W,n,h0,p,Xb,Pb,kbin,kernel,nh)
 call Interpola (Xb,Pb(1,1,1),kbin,X,muhatg,n)
 
+!print *, h0
+
+!print *, muhatg
 
 
 
@@ -761,9 +1007,9 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
  errg(i)=Y(i)-muhatg(i)
 end do
 
-do i=1,kbin
- pred0(i)=Pb(i,r+1,1)
-end do
+!do i=1,kbin
+! pred0(i)=Pb(i,r+1,1)
+!end do
 
 
 
@@ -777,10 +1023,40 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
  do i=1,kbin
   pred1(i,j)=Pb(i,r+1,1)
  end do
+ call Interpola (Xb,Pb(1,1,1),kbin,X,fpar(1,j),n)
+end do
+
+
+do i=1,n
+ do j=1,nf
+  if(F(i).eq.fact(j)) fpar_est(i)=fpar(i,j)
+ end do
 end do
 
+! !interpolamos efectos parciales para RSS
+! icont=0
+!   do i=1,n
+!    if (F(i).eq.fact(j)) then
+!     icont=icont+1
+!     Xaux(icont)=X(i)
+!    end if
+!   end do
+!   call Interpola (Xb,Pb(1,1,1),kbin,Xaux,fpar,icont)
+  
+!   ii=0
+!   do i=1,n
+!   if (F(i).eq.fact(j)) then
+!    ii=ii+1
+!    fpar_est(i)=fpar(ii)
+!   end if
+!   end do
+
+! end do
+
+
 
 
+!print *, h(1), h(2)
 
 
 ! estimo polinomios
@@ -813,21 +1089,56 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
 errg(1:n)=Y(1:n)-muhatg2(1:n)
 !**********************************
 
+!centro errores
+meanerr=sum(errg(1:n))/n
+do i=1,n
+ errg(i)=errg(i)-meanerr
+end do
 
 
 
 !Estadistico
 
-T=0
+T(1)=0
 do j=1,nf
 do i=1,kbin
+!do i=1,n
 ! T=T+abs(pred0(i)-pred1(i,j))
-T=T+abs(pred1(i,j))
+if(Xb(i).ge.-1.and.Xb(i).le.1) T(1)=T(1)+abs(pred1(i,j))
+!if(X(i).ge.-1.5.and.X(i).le.1.5) T(1)=T(1)+abs(fpar(i,j))
+end do
+end do
+
+
+
+! para la g
+T(2)=0
+do j=1,nf
+ Waux=0
+ do i=1,n
+  if (F(i).eq.fact(j)) Waux(i)=W(i)
+ end do
+ call rfast_h(X,errg,Waux,n,hg(j),p,Xb,Pb,kbin,kernel,nh)
+ do i=1,kbin
+   if(Xb(i).ge.-2.and.Xb(i).le.2) T(2)=T(2)+abs(Pb(i,1,1))
 end do
 end do
 
 
 
+RSS0=0
+RSS1=0
+do i=1,n
+ RSS0=RSS0+(Y(i)-muhatg2(i))**2
+ RSS1=RSS1+(Y(i)- muhatg(i) - fpar_est(i) )**2
+end do
+T(3)=RSS0-RSS1
+T(4)=T(3)/RSS1
+
+
+
+
+
 
 
 
@@ -847,6 +1158,9 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
   end if
 end do
 
+!h0 = h0i
+!h = hi
+!hg = h0i
 
 call rfast_h(X,Yboot,W,n,h0,p,Xb,Pb,kbin,kernel,nh)
 call Interpola (Xb,Pb(1,1,1),kbin,X,muhatgboot,n)
@@ -862,30 +1176,112 @@ subroutine globaltest_(F,X,Y,W,n,h0,h,nh,p,kbin,fact,nf,kernel,nboot,r,T,&
 end do
 
 
-
+!efectos parciales
 do j=1,nf
 Waux=0
 do i=1,n
 if (F(i).eq.fact(j)) Waux(i)=W(i)
 end do
 call rfast_h(X,errgboot,Waux,n,h(j),p,Xb,Pb,kbin,kernel,nh)
-
 do i=1,kbin
 pred1(i,j)=Pb(i,r+1,1)
 end do
+call Interpola (Xb,Pb(1,1,1),kbin,X,fpar(1,j),n) !interpolamos efectos parciales para RSS
+end do
+
+do i=1,n
+ do j=1,nf
+  if(F(i).eq.fact(j)) fpar_est(i)=fpar(i,j)
+ end do
+end do
+
+
+
+
+
+
+!polinomios
+
+
+! estimo polinomios
+hp=0 !ventana para pol
+
+
+if(r.eq.1) pp=0 !grado pol, solo calcula medias
+if(r.eq.2) pp=1
+
+
+do j=1,nf
+ Waux=0
+ do i=1,n
+  if (F(i).eq.fact(j)) Waux(i)=W(i)
+ end do
+ call rfast_h(X,errgboot,Waux,n,hp(j),pp,Xb,Pb,kbin,kernel,nh)
+ call Interpola (Xb,Pb(1,1,1),kbin,X,pol(1,j),n)
 end do
 
+if(r.eq.0) pol=0
 
 
-Tboot=0
+!**********************************
+do i=1,n
+ do j=1,nf
+  if(F(i).eq.fact(j)) muhatg2(i)=muhatgboot(i)+pol(i,j)
+ end do
+end do
+
+errgboot(1:n)=Yboot(1:n)-muhatg2(1:n)
+!**********************************
+
+
+
+
+
+
+
+Tboot(1)=0
 do k=1,nf
  do z=1,kbin
+!  do z=1,n
  ! Tboot=Tboot+abs(pred0(z)-pred1(z,k))
-  Tboot=Tboot+abs(pred1(z,k))
+  if(Xb(z).ge.-1.and.Xb(z).le.1) Tboot(1)=Tboot(1)+abs(pred1(z,k))
+!   if(X(z).ge.-1.and.X(z).le.1) Tboot(1)=Tboot(1)+abs(fpar(z,k))
  end do
 end do
 
-if(Tboot.gt.T) pvalor=pvalor+1
+
+
+
+
+! para la g
+Tboot(2)=0
+do j=1,nf
+ Waux=0
+ do i=1,n
+  if (F(i).eq.fact(j)) Waux(i)=W(i)
+ end do
+ call rfast_h(X,errgboot,Waux,n,hg(j),p,Xb,Pb,kbin,kernel,nh)
+ do i=1,kbin
+   if(Xb(i).ge.-2.and.Xb(i).le.2) Tboot(2)=Tboot(2)+abs(Pb(i,1,1))
+     ! if(Xb(i).ge.-2.and.Xb(i).le.2) T(2)=T(2)+abs(Pb(i,1,1))
+end do
+end do
+
+RSS0=0
+RSS1=0
+do i=1,n
+ RSS0=RSS0+(Yboot(i)-muhatg2(i))**2
+ RSS1=RSS1+(Yboot(i)- muhatgboot(i)-fpar_est(i) )**2
+end do
+Tboot(3)=RSS0-RSS1
+Tboot(4)=Tboot(3)/RSS1
+
+
+
+do j=1,4
+if(Tboot(j).gt.T(j)) pvalor(j)=pvalor(j)+1
+end do
+
 
 end do