R/dlsem.r

Defines functions causalEff cumulCalc findlag2sum dpathFind edgeMat isIndep morgraph angraph topOrder nodeDescen nodeAnces nodeMB chldsets plot.dlsem residualPlot edgeCoeff as.graphNEL makeGraph drawSample getStruct predict.dlsem residuals.dlsem fitted.dlsem formatFit compareModels logLik.dlsem confint.dlsem vcov.dlsem coef.dlsem nobs.dlsem print.summary.dlsem summary.dlsem formatSumm print.dlsem auto.lagPlot dlsem collCheck preProcess impoptAdj diffoptAdj autoCode findOp checkName checkNA lconAdj gconAdj intCheck lagPlot lagShapes makeShape EM.imputation addLags arFind splRecons linImp impuPred impuFit print.unirootTest doURT adft urtFun unirootTest isTimeVar applyDiff lagEff confint.hac summary.hac vcov.hac doHAC lmHAC W_hac dlaglm doLS deSeas getLev extrName creatForm searchGrid gammaParGen gamdefault gamlead scanForm getWBrk isQuant uncons gam ld qd ecq findLagLim Zmat genLag lagwei

Documented in as.graphNEL autoCode auto.lagPlot causalEff compareModels dlsem drawSample ecq gam isIndep lagPlot lagShapes ld lmHAC qd residualPlot unirootTest

defpar <- par()[setdiff(names(par()),c("cin","cra","csi","cxy","din","page"))]

#########################################
#
###  Available constrained lag shapes ###
#
# "ecq": endpoint-constrained quadratic
#  "qd": quadratic decreasing
#  "ld": linearly decreasing
# "gam": gamma
#
#########################################


# lag weights (internal use only)
lagwei <- function(theta,lag,type) {
  res <- c()
  xa <- 0  #####  
  for(i in 1:length(lag)) {
    if(type=="ecq") {
      if(lag[i]<theta[1] | lag[i]>theta[2]) {
        res[i] <- 0
        } else {
        res[i] <- -4/(theta[2]-theta[1]+2)^2*(lag[i]^2-(theta[1]+theta[2])*lag[i]+(theta[1]-1)*(theta[2]+1))
        }
      } else if(type=="qd") {
      if(lag[i]<theta[1] | lag[i]>theta[2]) {
        res[i] <- 0
        } else {
        res[i] <- (lag[i]^2-2*(theta[2]+1)*lag[i]+(theta[2]+1)^2)/(theta[2]-theta[1]+1)^2
        }
      } else if(type=="ld") {
      if(lag[i]<theta[1] | lag[i]>theta[2]) {
        res[i] <- 0
        } else {
        res[i] <- (theta[2]+1-lag[i])/(theta[2]+1-theta[1])
        }
      } else if(type %in% c(4,"gam")) {
      if(lag[i]>=xa) {
        bnum <- (lag[i]-xa+1)^(theta[1]/(1-theta[1]))*theta[2]^(lag[i]-xa)
        xM <- (theta[1]/(theta[1]-1))/log(theta[2])+xa-1
        bden <- (xM-xa+1)^(theta[1]/(1-theta[1]))*theta[2]^(xM-xa)
        res[i] <- bnum/bden
        } else {
        res[i] <- 0   
        }
      }
    }
  names(res) <- lag
  if(type %in% c(4,"gam")) {
    if(sum(res,na.rm=T)==0) res[1] <- 1
    }
  res
  }

# generate lagged instances of a variable (internal use only)
genLag <- function(x,maxlag,past=F) {
  if(maxlag>0) {
    if(past==T) x[which(is.na(x))] <- mean(x,na.rm=T)
    out <- x
    for(w in 1:maxlag) {
      if(w<length(x)) {
        if(past==F) {
          wx <- c(rep(NA,w),x[1:(length(x)-w)])
          } else {
          wx <- c(rep(mean(x,na.rm=T),w),x[1:(length(x)-w)])       
          }
        } else {
        if(past==F) {
          wx <- rep(NA,length(x))
          } else {
          wx <- rep(mean(x,na.rm=T),length(x))          
          }
        }
      out <- cbind(out,wx)        
      }
    colnames(out) <- NULL
    out
    } else {
    x
    }
  }

# distributed-lag transformation (internal use only)
Zmat <- function(x,type,theta,nlag) {
  if(type=="none") {
    matrix(x,ncol=1)
    } else {
    if(type=="gam") {
      xa <- 0  #####
      bhat <- gamlead(theta[1],theta[2])
      if(is.null(nlag)) {
        laglim <- min(bhat,trunc(2/3*length(x)))
        } else {
        laglim <- min(nlag,trunc(2/3*length(x)))
        }
      } else {
      if(is.null(nlag)) {
        laglim <- min(theta[2],trunc(2/3*length(x)))
        } else {
        laglim <- min(nlag,trunc(2/3*length(x)))
        }
      }
    H <- lagwei(theta,0:laglim,type)
    as.numeric(genLag(x,laglim)%*%matrix(H,ncol=1))
    }
  }

# compute the lag limit (internal use only)
findLagLim <- function(data,group=NULL) {
  if(is.null(group)) {
    trunc(nrow(na.omit(data))*2/3)
    } else {
    auxlaglim <- c()
    gruppi <- levels(factor(data[,group]))
    for(i in 1:length(gruppi)) {
      auxlaglim[i] <- nrow(na.omit(data[which(data[,group]==gruppi[i]),]))
      }
    trunc(min(auxlaglim,na.rm=T)*2/3)
    }
  }

# ECQ constructor
ecq <- function(x,a,b,x.group=NULL,nlag=NULL) {
  if(missing(a)==F & (length(a)!=1 || !is.numeric(a) || a<0 || a!=round(a))) stop("Argument 'a' must be a non-negative integer value",call.=F)
  if(missing(b)==F & (length(b)!=1 || !is.numeric(b) || b<0 || b!=round(b))) stop("Argument 'b' must be a non-negative integer value",call.=F)
  if(a>b) stop("Argument 'a' must be no greater than argument 'b'",call.=F)
  if(!is.null(nlag) & (length(nlag)!=1 || !is.numeric(nlag) || nlag<0 || nlag!=round(nlag))) stop("Argument 'nlag' must be a non-negative integer value",call.=F)
  if(is.null(x.group)) {
    res <- Zmat(x,"ecq",c(a,b),nlag)
    } else {
    res <- c()
    gruppi <- levels(factor(x.group))
    for(i in 1:length(gruppi)) {
      auxind <- which(x.group==gruppi[i])
      ires <- Zmat(x[auxind],"ecq",c(a,b),nlag)
      res[auxind] <- ires
      }
    }
  res
  }

# QD constructor
qd <- function(x,a,b,x.group=NULL,nlag=NULL) {
  if(missing(a)==F & (length(a)!=1 || !is.numeric(a) || a<0 || a!=round(a))) stop("Argument 'a' must be a non-negative integer value",call.=F)
  if(missing(b)==F & (length(b)!=1 || !is.numeric(b) || b<0 || b!=round(b))) stop("Argument 'b' must be a non-negative integer value",call.=F)
  if(a>b) stop("Argument 'a' must be no greater than argument 'b'",call.=F)
  if(!is.null(nlag) & (length(nlag)!=1 || !is.numeric(nlag) || nlag<0 || nlag!=round(nlag))) stop("Argument 'nlag' must be a non-negative integer value",call.=F)
  if(is.null(x.group)) {
    res <- Zmat(x,"qd",c(a,b),nlag)
    } else {
    res <- c()  
    gruppi <- levels(factor(x.group))
    for(i in 1:length(gruppi)) {
      auxind <- which(x.group==gruppi[i])
      ires <- Zmat(x[auxind],"qd",c(a,b),nlag)
      res[auxind] <- ires
      }
    }
  res
  }

# LD constructor
ld <- function(x,a,b,x.group=NULL,nlag=NULL) {
  if(missing(a)==F & (length(a)!=1 || !is.numeric(a) || a<0 || a!=round(a))) stop("Argument 'a' must be a non-negative integer value",call.=F)
  if(missing(b)==F & (length(b)!=1 || !is.numeric(b) || b<0 || b!=round(b))) stop("Argument 'b' must be a non-negative integer value",call.=F)
  if(a>b) stop("Argument 'a' must be no greater than argument 'b'",call.=F)
  if(!is.null(nlag) & (length(nlag)!=1 || !is.numeric(nlag) || nlag<0 || nlag!=round(nlag))) stop("Argument 'nlag' must be a non-negative integer value",call.=F)
  if(is.null(x.group)) {
    res <- Zmat(x,"ld",c(a,b),nlag)
    } else {
    res <- c()  
    gruppi <- levels(factor(x.group))
    for(i in 1:length(gruppi)) {
      auxind <- which(x.group==gruppi[i])
      ires <- Zmat(x[auxind],"ld",c(a,b),nlag)
      res[auxind] <- ires
      }
    }
  res
  }

# GAM constructor
gam <- function(x,a,b,x.group=NULL,nlag=NULL) {
  if(missing(a)==F & (length(a)!=1 || !is.numeric(a) || a<=0 || a>=1)) stop("Argument 'a' must be a value in the interval (0,1)",call.=F)
  if(missing(b)==F & (length(b)!=1 || !is.numeric(b) || b<=0 || b>=1)) stop("Argument 'b' must be a value in the interval (0,1)",call.=F)
  if(!is.null(nlag) & (length(nlag)!=1 || !is.numeric(nlag) || nlag<0 || nlag!=round(nlag))) stop("Argument 'nlag' must be a non-negative integer value",call.=F)
  if(is.null(x.group)) {
    res <- Zmat(x,"gam",c(a,b),nlag)
    } else {
    res <- c()
    gruppi <- levels(factor(x.group))
    for(i in 1:length(gruppi)) {
      auxind <- which(x.group==gruppi[i])
      ires <- Zmat(x[auxind],"gam",c(a,b),nlag)
      res[auxind] <- ires
      }
    }
  res
  }

# unconstrained constructor (internal use only)
uncons <- function(x,a,b,x.group=NULL) {
  if(missing(a)==F & (length(a)!=1 || !is.numeric(a) || a<0 || a!=round(a))) stop("Argument 'a' must be a non-negative integer value",call.=F)
  if(missing(b)==F & (length(b)!=1 || !is.numeric(b) || b<0 || b!=round(b))) stop("Argument 'b' must be a non-negative integer value",call.=F)
  if(a>b) stop("Argument 'a' must be no greater than argument 'b'",call.=F)
  if(is.null(x.group)) {
    if(b>0) {
      res <- x
      for(i in 1:b) {
        res <- cbind(res,c(rep(NA,i),x[1:(length(x)-i)]))
        }
      colnames(res) <- 0:b
      } else {
      res <- x
      }
    } else {
    res <- c()  
    gruppi <- levels(factor(x.group))
    for(i in 1:length(gruppi)) {
      auxind <- which(x.group==gruppi[i])
      ires <- idat <- x[auxind]
      for(j in 1:b) {
        ires <- cbind(ires,c(rep(NA,j),idat[1:(length(auxind)-j)]))
        }
      res <- rbind(res,ires)
      }
    colnames(res) <- 0:b
    }
  res[,(a+1):(b+1)]
  }

# check if a variable is quantitative (internal use only)
isQuant <- function(x) {
  if(is.numeric(x)) {
    if(identical(sort(unique(na.omit(x))),c(0,1))) F else T
    } else {
    F
    }
  }

# get info within brackets (internal use only)
getWBrk <- function(x) {
  regmatches(x, gregexpr("(?<=\\().*?(?=\\))",x,perl=T))[[1]]
  #xstr <- strsplit(x,"\\(")[[1]]
  #if(length(xstr)==2) {
  #  gsub("\\)","",xstr[2])
  #  } else {
  #  paste(paste(gsub("\\)","",xstr[-1]),collapse="("),paste(rep(")",length(xstr)-2),collapse=""),sep="")
  #  }
  }

# scan formula (internal use only)
scanForm <- function(x,warn=F) {
  auxform <- gsub(" ","",formula(x))[-1]
  ynam <- gsub(" ","",auxform[1])
  auX <- gsub(" ","",strsplit(gsub("-","+",auxform[2]),"\\+")[[1]])
  auX <- setdiff(auX,c("1",""))
  if(length(auX)>0) {
    lnames <- ltype <- rep(NA,length(auX))
    lpar <- list()
    for(i in 1:length(auX)) {        
      if(nchar(auX[i])>0) {
        # check deprecated constructors  <-------------------- !!
        if(identical("quec.lag(",substr(auX[i],1,9))) {
          auX[i] <- gsub("^quec.lag\\(","ecq\\(",auX[i])
          if(warn==T) warning("The constructor 'quec.lag()' is deprecated, 'ecq' was used instead",call.=F)
          }
        if(identical("qdec.lag(",substr(auX[i],1,9))) {
          auX[i] <- gsub("^qdec.lag\\(","qd\\(",auX[i])
          if(warn==T) warning("The constructor 'qdec.lag()' is deprecated, 'qd' was used instead",call.=F)
          }
        if(identical("gamma.lag(",substr(auX[i],1,10))) {
          auX[i] <- gsub("^gamma.lag\\(","gam\\(",auX[i])
          if(warn==T) warning("The constructor 'gamma.lag()' is deprecated, 'gam' was used instead",call.=F)
          }
        if(identical("gamm.lag(",substr(auX[i],1,9))) {
          auX[i] <- gsub("^gamm.lag\\(","gam\\(",auX[i])
          if(warn==T) warning("The constructor 'gamm.lag()' is deprecated, 'gam' was used instead",call.=F)
          }
        # find parameters in lag shape constructors
        if(strsplit(auX[i],"\\(")[[1]][1] %in% c("ecq","qd","ld","gam")) {
          istr <- strsplit(getWBrk(auX[i]),",")[[1]]
          ltype[i] <- strsplit(auX[i],"\\(")[[1]][1]
          lnames[i] <- istr[1]
          #auX[i] <- istr[1]
          lpar[[i]] <- as.numeric(istr[2:3])
          } else {
          lnames[i] <- auX[i]
          ltype[i] <- "none"
          lpar[[i]] <- NA
          }
        }
      }
    names(lpar) <- names(ltype) <- lnames
    } else {
    lpar <- ltype <- c()
    }
  list(y=ynam,X=auX,ltype=ltype,lpar=lpar)
  }

# lead lag of a gamma lag shape (internal use only)
gamlead <- function(delta,lambda,tol=1e-4) {
  m <- ceiling(delta/((delta-1)*log(lambda)))
  res <- m
  b <- lagwei(c(delta,lambda),res,"gam")
  ind <- 1
  while(b>tol && ind<10*m) {
    res <- res+1
    b <- lagwei(c(delta,lambda),res,"gam")
    ind <- ind+1
    }
  res
  }

# default gamma lag shape (internal use only)
gamdefault <- function(maxlag) {
  G <- matrix(
  c(0.115, 0.005,
    0.05, 0.02,
    0.46, 0.04,
    0.42, 0.08,
    0.18, 0.13,
    0.32, 0.17,
    0.35, 0.21,
    0.35, 0.25,
    0.47, 0.28,
    0.51, 0.31,
    0.38, 0.35,
    0.47, 0.37,
    0.48, 0.40,
    0.48, 0.42,
    0.43, 0.45,
    0.53, 0.46,
    0.51, 0.48),byrow=T,ncol=2)
  if(maxlag<18) {
    G[max(1,maxlag),]
    } else {
    c(0.5,0.5)
    }
  }

# generate gamma parameters (internal use only)
gammaParGen <- function(by) {
  xseq <- seq(0+by,1-by,by=by)
  auxmat <- as.matrix(expand.grid(xseq,xseq))
  limmat <- c()
  for(i in 1:nrow(auxmat)) {
    limmat <- c(limmat,gamlead(auxmat[i,1],auxmat[i,2]))
    }
  res <- cbind(auxmat,limmat)
  colnames(res) <- c("delta","lambda","lead_lag")
  res
  }

# generate search grid (internal use only)
searchGrid <- function(mings,maxgs,minwd,maxld,lag.type,gammaMat) {
  if(lag.type=="gam") {
    gammaMat[which(gammaMat[,3]<=maxld & gammaMat[,3]>=minwd),1:2]
    } else {
    auxmat <- c()
    for(i in 0:maxld) {
      for(j in i:maxld) { 
        if(i>=mings & i<=maxgs & j-i>=minwd) auxmat <- rbind(auxmat,c(i,j))
        }   
      }
    auxmat
    }
  }

# create lm formula (internal use only)
creatForm <- function(y,X,group,type,theta,nlag) {
  xnam <- c()
  if(length(X)>0) {
    nomi <- setdiff(names(theta),group)
    for(i in 1:length(nomi)) {
      ilab <- nomi[i]
      if(type[ilab]=="none") {
        xnam[i] <- ilab
        } else {
        ixnam <- paste(type[ilab],"(",ilab,",",theta[[ilab]][1],",",theta[[ilab]][2],sep="")
        if(!is.null(group)) {
          if(!is.null(nlag)) {
            ixnam <- paste(ixnam,",",group,",",nlag,sep="")
            } else {
            ixnam <- paste(ixnam,",",group,sep="")  
            }
          } else {
          if(!is.null(nlag)) ixnam <- paste(ixnam,",,",nlag,sep="")
          }
        xnam[i] <- paste(ixnam,")",sep="")
        }
      }
    }
  if(is.null(group)) {
    if(length(X)>0) {
      res <- paste(y,"~",paste(xnam,collapse="+"),sep="")    
      } else {
      res <- paste(y,"~1",sep="")
      }
    } else {
    if(length(X)>0) {
      res <- paste(y,"~-1+",group,"+",paste(xnam,collapse="+"),sep="")    
      } else {
      res <- paste(y,"~-1+",group,sep="")
      }
    }
  formula(res)
  }

# extract the name from a lag shape (internal use only)
extrName <- function(x) {
  auxstr <- strsplit(x,"\\(")[[1]]
  if(auxstr[1] %in% c("ecq","qd","ld","gam")) {
    gsub("\\)","",strsplit(auxstr[2],",")[[1]][1])
    } else {
    x  
    }
  }

# get levels of variables in data (internal use only)
getLev <- function(data) {
  auxq <- sapply(data,isQuant)
  res <- list()
  for(i in 1:length(auxq)) {
    if(auxq[i]==T) {
      res[[i]] <- NA
      } else {
      res[[i]] <- paste(colnames(data)[i],levels(factor(data[,i])),sep="")
      }
    }
  names(res) <- colnames(data)
  res
  }

# deseasonalization (internal use only)
deSeas <- function(x,seas,group,data) {
  res <- data
  if(is.null(group)) {
    xseas <- factor(rep(1:seas,length=nrow(data)))
    if(nlevels(xseas)>=2 && min(table(xseas))>=2) {
      for(i in 1:length(x)) {
        iform <- paste(x[i],"~xseas",sep="")
        imod <- lm(formula(iform),data=data)
        res[,x[i]] <- mean(data[,x[i]],na.rm=T)+residuals(imod)
        }
      }
    } else {
    gruppi <- levels(factor(data[,group]))
    for(w in 1:length(gruppi)) {
      auxind <- which(data[,group]==gruppi[w])
      xseas <- factor(rep(1:seas,length=length(auxind)))
      if(nlevels(xseas)>=2 && min(table(xseas))>=2) {
        for(i in 1:length(x)) {
          iform <- paste(x[i],"~xseas",sep="")
          imod <- lm(formula(iform),data=data[auxind,])
          res[auxind,x[i]] <- mean(data[auxind,x[i]],na.rm=T)+residuals(imod)
          }
        }
      }
    }
  res
  }

# perform ols estimation (internal use only)
doLS <- function(formula,group,data) {
  formOK <- formula
  Xm0 <- model.matrix(formOK,data=data)
  auxdel <- names(which(apply(Xm0,2,var)==0))
  if(length(auxdel)>0) {
    auxlev <- getLev(data)
    x2del <- c()
    for(i in 1:length(auxdel)) {
      if(auxdel[i] %in% colnames(data)) {
        x2del <- c(x2del,auxdel[i])
        } else {
        auxf <- sapply(auxlev,function(z){auxdel[i] %in% z})
        x2del <- c(x2del,names(which(auxf==T)))
        }
      }
    x2del <- setdiff(x2del,group)
    if(length(x2del)>0) {
      auxform <- scanForm(formula)
      if(is.null(group)) auxsep <- "" else auxsep <- "-1+"
      formOK <- formula(paste(auxform$y,"~",auxsep,paste(setdiff(auxform$X,x2del),collapse="+"),sep=""))
      }
    }
  res <- lm(formOK,data=data)
  res$call$formula <- formOK
  res
  }

# fit a distributed-lag linear regression (internal use only)
dlaglm <- function(formula,group,data,adapt,no.select,min.gestation,max.gestation,min.width,max.lead,sign,ndiff,gamma.by,mess,nblank) {
  auxscan <- scanForm(formula)
  y <- auxscan$y
  if(length(auxscan$X)>0) {
    lagPar <- auxscan$lpar
    lagType <- auxscan$ltype
    lagNam <- names(lagPar)
    no.lag <- names(which(lagType=="none"))
    xOK <- c(no.select,no.lag)
    xtest <- setdiff(lagNam,xOK)
    if(length(xtest)==0) adapt <- F
    if(adapt==F) {
      warn1 <- warn2 <- 0
      for(i in 1:length(lagPar)) {
        ilimit <- findLagLim(data[,c(y,lagNam[i],group)],group=group)
        if(lagType[i] %in% c("ecq","qd","ld")) {
          if(lagPar[[i]][1]!=round(lagPar[[i]][1])) {
            lagPar[[i]][1] <- round(lagPar[[i]][1])
            warn1 <- warn1+1
            }
          if(lagPar[[i]][1]<0) {
            lagPar[[i]][1] <- 0
            warn1 <- warn1+1
            }
          if(lagPar[[i]][2]!=round(lagPar[[i]][2])) {
            lagPar[[i]][2] <- round(lagPar[[i]][2])
            warn1 <- warn1+1
            }
          if(lagPar[[i]][2]<0) {
            lagPar[[i]][2] <- 0
            warn1 <- warn1+1
            }
          if(lagPar[[i]][1]>lagPar[[i]][2]) {
            lagPar[[i]][2] <- lagPar[[i]][1]
            warn1 <- warn1+1
            }
          if(lagPar[[i]][1]>ilimit) {
            lagPar[[i]][1] <- ilimit
            warn2 <- warn2+1
            }
          if(lagPar[[i]][2]>ilimit) {
            lagPar[[i]][2] <- ilimit
            warn2 <- warn2+1
            }
          } else if(lagType[i]=="gam") {
          if(lagPar[[i]][1]<=0) {
            lagPar[[i]][1] <- 0.01
            warn1 <- warn1+1
            }
          if(lagPar[[i]][1]>=1) {
            lagPar[[i]][1] <- 0.99
            warn1 <- warn1+1
            }
          if(lagPar[[i]][2]<=0) {
            lagPar[[i]][2] <- 0.01
            warn1 <- warn1+1
            }
          if(lagPar[[i]][2]>=1) {
            lagPar[[i]][2] <- 0.99
            warn1 <- warn1+1
            }
          iglim <- gamlead(lagPar[[i]][1],lagPar[[i]][2])
          if(iglim>ilimit) {
            lagPar[[i]] <- gamdefault(ilimit)
            warn2 <- warn2+1
            }
          }
        }
      if(warn1>0) warning("Invalid lag shapes in the regression of '",y,"' replaced with the nearest valid ones",call.=F)
      if(warn2>0) warning("Too large lead lags in the regression of '",y,"' replaced with the maximum possible ones",call.=F)
      if(!is.null(mess)) {
        iblchar <- ""
        if(nblank>0) iblchar <- paste(rep(" ",nblank),collapse="")
        cat('\r')
        cat(paste(mess,iblchar,sep=""))
        flush.console()
        }
      formOK <- creatForm(y,names(lagPar),group,lagType,lagPar,NULL)
      modOK <- doLS(formula=formOK,group=group,data=data)
      } else {
      if(sum(lagType=="gam")>0) {
        gammaMat <- gammaParGen(gamma.by)
        } else {
        gammaMat <- NULL  
        }
      bestPar <- vector("list",length=length(lagNam))
      names(bestPar) <- lagNam
      consList <- list()
      for(i in 1:length(xtest)) {
        imings <- min.gestation[xtest[i]]
        imaxgs <- max.gestation[xtest[i]]
        iminwd <- min.width[xtest[i]]
        imaxld <- max.lead[xtest[i]]
        icons <- searchGrid(imings,imaxgs,iminwd,imaxld,lagType[xtest[i]],gammaMat)
        if(lagType[[xtest[i]]]=="gam") {
          igamdef <- gamdefault(imaxld)
          if(nrow(icons)==0) icons <- igamdef
          bestPar[[xtest[i]]] <- igamdef
          } else {
          bestPar[[xtest[i]]] <- c(imings,imaxld)
          }
        consList[[i]] <- icons
        }
      names(consList) <- xtest
      nittVet <- sapply(consList,nrow)
      nittTot <- 0
      for(i in 1:length(nittVet)) {
        nittTot <- nittTot+sum(nittVet[i:length(nittVet)])
        }
      fine <- nitt <- 0
      while(fine==0) {
        xtest <- setdiff(lagNam,xOK)
        ntest <- length(xtest)
        if(ntest>0) {
          currentBIC <- rep(NA,ntest) 
          currentPar <- vector("list",length=length(xtest))
          names(currentBIC) <- names(currentPar) <- xtest
          for(i in 1:ntest) {
            auxcons <- consList[[xtest[i]]]
            if(nrow(auxcons)==0) auxcons <- matrix(lagPar[[xtest[i]]],nrow=1)
            mse0 <- bhat0 <- c()
            for(j in 1:nrow(auxcons)) {
              nitt <- nitt+1
              iperc <- round(100*nitt/nittTot)
              if(!is.null(mess) && iperc%%5==0) {
                iblchar <- ""
                if(nblank>0) iblchar <- paste(rep(" ",nblank),collapse="")
                imess <- paste(mess," ... ",iperc,"%",iblchar,sep="")
                cat('\r')
                cat(imess)
                flush.console()
                }
              #
              testType <- lagType
              #testType[setdiff(xtest,xtest[i])] <- "none"  #####
              testPar <- bestPar
              testPar[[xtest[i]]] <- auxcons[j,]
              #
              form0 <- creatForm(y,names(testPar),group,testType,testPar,max.lead[xtest[i]])  #####
              #form0 <- creatForm(y,names(testPar),group,testType,testPar,NULL)
              #
              mod0 <- doLS(formula=form0,group=group,data=data)
              est0 <- mod0$coefficients
              ixall <- names(est0)
              iauxlab <- sapply(ixall,extrName)
              ixlab <- ixall[which(iauxlab==xtest[i])]
              #
              if(length(ixlab)>0 && ixlab %in% ixall) { 
                bhat0[j] <- est0[ixlab]
                mse0[j] <- mean(residuals(mod0)^2)
                } else {
                bhat0[j] <- NA
                mse0[j] <- Inf
                }
              }
            isign <- sign[xtest[i]]
            #if(!is.null(isign)) {
            if(isign!=F) {
              if(isign=="+") {
                auxsign <- which(bhat0>0)
                } else {
                auxsign <- which(bhat0<0)                  
                }
              if(length(auxsign)>0) {
                auxbest <- auxsign[which.min(mse0[auxsign])]
                } else {
                auxbest <- which.min(mse0)
                }
              } else {
              auxbest <- which.min(mse0)
              }
            currentPar[[xtest[i]]] <- auxcons[auxbest,]
            currentBIC[xtest[i]] <- mse0[auxbest]
            }
          xnew <- names(currentBIC)[which.min(currentBIC)]
          bestPar[[xnew]] <- currentPar[[xnew]]
          xOK <- c(xOK,xnew)
          } else {
          fine <- 1
          }
        }
      formOK <- creatForm(y,names(bestPar),group,lagType,bestPar,NULL)  #####
      modOK <- doLS(formula=formOK,group=group,data=data)
      }
    } else {
    if(is.null(group)) {
      formOK <- formula(paste(y,"~1",sep=""))
      } else {
      formOK <- formula(paste(y,"~-1+",group,sep=""))        
      }
    modOK <- doLS(formula=formOK,group=group,data=data)
    }
  modOK
  }

# HAC weights (internal use only)
W_hac <- function(Xmat,res,maxlag) {
  n <- nrow(Xmat)
  p <- ncol(Xmat)
  W <- matrix(0,nrow=p,ncol=p)
  for(i in 1:n) {
    W <- W+res[i]^2*Xmat[i,]%*%t(Xmat[i,])
    }
  if(maxlag>0) {
    for(j in 1:maxlag) {
      wi <- 0
      for(i in (j+1):n) {
        wi <- wi+res[i]*res[i-j]*(Xmat[i,]%*%t(Xmat[i-j,])+Xmat[i-j,]%*%t(Xmat[i,]))
        }
      W <- W+(1-j/(1+maxlag))*wi
      }
    }
  W
  }

# apply HAC covariance matrix to a lm object
lmHAC <- function(x,group=NULL) {
  if(("lm" %in% class(x))==F & ("dlsem" %in% class(x))==F) stop("Argument 'x' must be an object of class 'lm' or 'dlsem'",call.=F)
  if("lm" %in% class(x)) {
    x$vcov <- doHAC(x=x,group=group)
    class(x) <- c("hac","lm")
    } else {
    vcovL <- lapply(x$estimate,doHAC,group=group)
    for(i in 1:length(x$estimate)) {
      x$estimate[[i]] <- vcovL[[i]]
      class(x$estimate[[i]]) <- c("hac","lm")
      }
    }
  x
  }

# HAC for class 'lm' (internal use only)
doHAC <- function(x,group) {
  Xmat <- model.matrix(x)
  Smat <- summary(x)$cov.unscaled
  Xmat <- Xmat[,colnames(Smat),drop=F]  ## delete collinear terms
  n <- nrow(Xmat)
  p <- ncol(Xmat)
  res <- x$residuals
  if(!is.null(group) && is.na(group)) group <- NULL
  if(is.null(group)) {
    maxlag <- ar(res)$order
    W <- W_hac(Xmat,res,maxlag)
    } else {
    glev <- x$xlevels[[group]]
    gnam <- paste(group,glev,sep="")
    Wsum <- matrix(0,nrow=ncol(Xmat),ncol=ncol(Xmat))
    W <- matrix(0,nrow=p,ncol=p)
    maxlag <- c()
    for(i in 1:length(gnam)) {
      if(gnam[i] %in% colnames(Xmat)) {
        iind <- names(which(Xmat[,gnam[i]]==1))
        } else {
        iind <- names(which(apply(Xmat[,gnam[setdiff(1:length(gnam),i)]],1,sum)==0))
        }
      maxlag[i] <- ar(res[iind])$order
      W <- W+W_hac(Xmat[iind,],res[iind],maxlag[i])
      }
    names(maxlag) <- gnam
    }
  out <- n/(n-p)*Smat%*%W%*%Smat
  attr(out,"max.lag") <- maxlag
  out
  }

# vcov method for class 'hac'
vcov.hac <- function(object,...)  {
  object$vcov
  }

# summary method for class 'hac'
summary.hac <- function(object,...)  {
  res <- summary.lm(object)
  res$coefficients[colnames(object$vcov),2] <- sqrt(diag(object$vcov))
  res$coefficients[,3] <- res$coefficients[,1]/res$coefficients[,2]
  res$coefficients[,4] <- 2*pt(-abs(res$coefficients[,3]),object$df.residual)
  res
  }

# confint method for class 'hac'
confint.hac <- function(object,parm,level=0.95,...) {
  summ <- summary(object)$coefficients
  quan <- qt((1+level)/2,object$df.residual)
  res <- cbind(summ[,1]-quan*summ[,2],summ[,1]+quan*summ[,2])
  colnames(res) <- paste(100*c(1-level,1+level)/2," %",sep="")
  res
  }
  
# compute lag effects of a covariate (internal use only)
lagEff <- function(model,x,cumul,lag) {
  formstr <- strsplit(gsub(" ","",as.character(model$call$formula)[3]),"\\+")[[1]]
  xall <- names(model$coefficients)
  auxlab <- sapply(xall,extrName)
  xlab <- xall[which(auxlab==x)]
  auxscan <- scanForm(model$call)
  if(auxscan$ltype[x] %in% c("ecq","qd","ld")) {
    sx <- auxscan$lpar[[x]][1]
    dx <- auxscan$lpar[[x]][2]
    imu <- model$coefficients[xlab]
    icov <- vcov(model)[xlab,xlab]
    if(is.null(lag)) {
      xgrid <- 0:dx
      } else {
      xgrid <- lag
      }
    iH <- matrix(lagwei(c(sx,dx),xgrid,auxscan$ltype[x]),ncol=1)
    } else if(auxscan$ltype[x]=="gam") {
    delta <- auxscan$lpar[[x]][1]
    lambda <- auxscan$lpar[[x]][2]
    ilim <- c(0,gamlead(delta,lambda))  #####
    imu <- model$coefficients[xlab]
    icov <- vcov(model)[xlab,xlab]
    if(is.null(lag)) {
      xa <- 0  #####
      xgrid <- 0:ilim[2]
      } else {
      xgrid <- lag
      }
    iH <- matrix(lagwei(c(delta,lambda),xgrid,"gam"),ncol=1)
    idel <- setdiff(xgrid,ilim[1]:ilim[2])
    if(length(idel)>0) iH[sapply(idel,function(z){which(xgrid==z)}),] <- 0
    } else {  
    xgrid <- 0  
    imu <- model$coefficients[x]
    icov <- vcov(model)[x,x]
    iH <- matrix(1,nrow=1,ncol=1)
    }
  ibhat <- iH%*%imu
  ibse <- sqrt(diag(iH%*%icov%*%t(iH)))
  #
  out <- cbind(ibhat,ibse)
  rownames(out) <- xgrid
  colnames(out) <- c("estimate","std. err.")
  if(cumul==T) out <- cumulCalc(out)
  out
  }

# apply differentiation (internal use only)
applyDiff <- function(x,group,data,k) {
  #
  diffFun <- function(z,k) {
    if(k>0 & k<length(z)) {
      res <- z
      for(i in 1:k) {
        res <- res-c(NA,res[1:(length(res)-1)])
        }
      res
      } else if(k<=0) {
      z
      } else {
      rep(NA,length(z))
      }
    }
  #
  diffdat <- data
  if(is.null(group)) {
    for(w in 1:length(x)) {
      if(isQuant(data[,x[w]])) {
        wdat <- data[,x[w]]
        diffdat[,x[w]] <- diffFun(wdat,k[w])      
        }
      }
    } else {
    data[,group] <- factor(data[,group])
    gruppi <- levels(factor(data[,group]))
    for(i in 1:length(gruppi)) {
      auxind <- which(data[,group]==gruppi[i])
      for(w in 1:length(x)) {
        if(isQuant(data[,x[w]])) {
          wdat <- data[auxind,x[w]]
          diffdat[auxind,x[w]] <- diffFun(wdat,k[w])
          }
        }
      } 
    }
  diffdat
  }

# check the time variable (internal use only)
isTimeVar <- function(x) {
  #
  standarDates <- function(string) {
    patterns = c('[0-9][0-9][0-9][0-9]/[0-9][0-9]/[0-9][0-9]','[0-9][0-9]/[0-9][0-9]/[0-9][0-9][0-9][0-9]','[0-9][0-9][0-9][0-9]-[0-9][0-9]-[0-9][0-9]')
    formatdates = c('%Y/%m/%d','%d/%m/%Y','%Y-%m-%d')
    standardformat='%d/%m/%Y'
    for(i in 1:3){
      if(grepl(patterns[i], string)){
        aux=as.Date(string,format=formatdates[i])
        if(!is.na(aux)){
          format(aux, standardformat)
          }
        }
      }
    F
    }
  #
  res <- T
  if(!is.numeric(x)) {
    if(F %in% sapply(x,standarDates)) res <- F
    }
  res
  }

# unit root test
unirootTest <- function(x=NULL,group=NULL,time=NULL,data,test=NULL,log=FALSE) {
  if(missing(data)==F & !identical(class(data),"data.frame")) stop("Argument 'data' must be a data.frame",call.=F)
  if(!is.null(x)) { 
    x2del <- c()
    for(i in 1:length(x)) {
      if((x[i] %in% colnames(data))==F) {
        warning("Variable '",x[i],"' not found in data and ignored",call.=F)
        x2del <- c(x2del,x[i])
        } else {
        if(isQuant(data[,x[i]])==F) {
          warning("'",x[i],"' is not a quantitative variable and was ignored",call.=F)
          x2del <- c(x2del,x[i])
          }
        }
      }
    x <- setdiff(x,x2del)
    if(length(x)<1) stop("No quantitative variables provided to argument 'x'",call.=F)
    } else {
    allnam <- setdiff(colnames(data),c(group,time))
    for(i in 1:length(allnam)) {
      if(isQuant(data[,allnam[i]])) x <- c(x,allnam[i])
      }
    }
  if(!is.null(group)) {
    if(is.na(group)) group <- NULL
    if(length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
    if((group %in% colnames(data))==F) stop("Variable '",group,"' provided to argument 'group' not found in data",call.=F)
    if(group %in% x) stop("Variable '",group,"' is provided to both arguments 'x' and 'group'",call.=F)
    if(group %in% time) stop("Variable '",group,"' is provided to both arguments 'group' and 'time'",call.=F)
    data[,group] <- factor(data[,group])
    gruppi <- levels(data[,group])
    if(length(gruppi)<2) stop("The group factor must have at least 2 unique values",call.=F)
    n <- min(table(data[,group]))
    if(n<3) stop("There must be at least 3 observations per group",call.=F)
    } else {
    n <- nrow(data)
    if(n<3) stop("There must be at least 3 observations",call.=F)  
    }
  if(!is.null(time)) {
    if(is.na(time)) time <- NULL
    if(length(time)!=1) stop("Argument 'time' must be of length 1",call.=F)
    if((time %in% colnames(data))==F) stop("Variable '",time,"' provided to argument 'time' not found in data",call.=F)
    if(time %in% x) stop("Variable '",time,"' is provided to both arguments 'x' and 'time'",call.=F)
    if(time %in% group) stop("Variable '",time,"' is provided to both arguments 'group' and 'time'",call.=F)
    if(isTimeVar(data[,time])==F) stop("The time variable is neither numeric nor a date",call.=F)
    if(is.null(group)) {
      if(sum(duplicated(data[,time]))>0) stop("The time variable has duplicated values",call.=F)
      } else {
      timesplit <- split(data[,time],data[,group])
      if(sum(sapply(timesplit,function(z){sum(duplicated(z))}))>0) stop("The time variable has duplicated values",call.=F)  
      }
    }
  if(is.null(test)) {
    test <- ifelse(n<100,"kpss","adf")
    } else {
    if((test %in% c("kpss","adf"))==F) stop("Argument 'test' must be one among 'kpss' or 'adf'",call.=F)  
    }
  if(identical(log,T)) {
    for(i in 1:length(x)) {
      if(sum(data[,x[i]]<=0,na.rm=T)>0) {
        warning("Logarithmic transformation not applied to variable '",x[i],"'",call.=F)
        } else {
        data[,x[i]] <- log(data[,x[i]])
        }
      }
    } else if(!identical(log,F)) {
    for(i in 1:length(log)) {
      if((log[i] %in% colnames(data))==F) {
        warning("Variable '",log[i],"' provided to argument 'log' not found in data",call.=F)
        } else if(log[i] %in% c(group,time)) {
        warning("Logarithmic transformation not applied to variable '",log[i],"'",call.=F)
        } else {
        if(sum(data[,log[i]]<=0,na.rm=T)>0) {
          warning("Logarithmic transformation not applied to variable '",log[i],"'",call.=F)
          } else {
          data[,log[i]] <- log(data[,log[i]])  
          }
        }
      }
    }
  urtFun(x,group,time,data,test,log)
  }
  
# interface for unit root test (internal use only)
urtFun <- function(x,group,time,data,test,log) {
  if(!is.null(group)) {
    g.id <- as.numeric(data[,group])
    glab <- levels(data[,group])
    } else {
    g.id <- glab <- rep(1,nrow(data))
    }
  data[which(abs(data[,x])==Inf),x] <- NA
  if(!is.null(time)) {
    for(i in 1:length(g.id)) {
      auxind <- which(g.id==glab[i])
      idat <- data[auxind,]
      data[auxind,] <- idat[order(idat[,time]),]
      }
    }
  res <- vector("list",length=length(x))
  for(i in 1:length(x)) {
    if(test=="adf") {
      if(is.null(group)) {
        ikmax <- trunc((length(na.omit(data[,x[i]]))-1)^(1/3))
        } else {
        ik <- table(na.omit(data[,c(x[i],group)])[,group])
        ikmax <- sapply(ik,function(z){trunc((z-1)^(1/3))})
        }
      res[[i]] <- doURT(data[,x[i]],g.id=g.id,test="adf",glab=glab,par=ikmax)
      } else {
      res[[i]] <- doURT(data[,x[i]],g.id=g.id,test="kpss",glab=glab,par=c(F,T))  #####
      }
    }
  names(res) <- x
  if(!is.null(group)) attr(res,"group") <- group
  attr(res,"test") <- test
  class(res) <- "unirootTest"
  res
  }

# adf test (internal use only)
adft <- function(x,kmax) {
  #
  doADF <- function(k) {
    y <- diff(x)
    n <- length(y)
    k <- k+1
    z <- embed(y,k)
    yt <- z[,1]
    xt1 <- x[k:n]
    tt <- k:n
    if(k>1) {
      yt1 <- z[,2:k,drop=F]
      res <- lm(yt~xt1+tt+yt1)
      } else {
      res <- lm(yt~xt1+tt)
      }
    res.sum <- summary(res)$coefficients
    if(nrow(res.sum)>=2) {
      STAT <- res.sum[2,1]/res.sum[2,2]
      table <- -1*cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96),
                        c(3.95, 3.8, 3.73, 3.69, 3.68, 3.66),
                        c(3.6, 3.5, 3.45, 3.43, 3.42, 3.41),
                        c(3.24, 3.18, 3.15, 3.13, 3.13, 3.12),
                        c(1.14, 1.19, 1.22, 1.23, 1.24, 1.25),
                        c(0.8, 0.87, 0.9, 0.92, 0.93, 0.94),
                        c(0.5, 0.58, 0.62, 0.64, 0.65, 0.66),
                        c(0.15, 0.24, 0.28, 0.31, 0.32, 0.33))
      tablen <- dim(table)[2]
      tableT <- c(25, 50, 100, 250, 500, 1e+05)
      tablep <- c(0.01, 0.025, 0.05, 0.1, 0.9, 0.95, 0.975, 0.99)
      tableipl <- numeric(tablen)
      for(i in (1:tablen)) {
        tableipl[i] <- approx(tableT,table[,i],n,rule=2)$y
        }
      PVAL <- approx(tableipl,tablep,STAT,rule=2)$y
      } else {
      STAT <- PVAL <- NA
      }
    c(STAT,PVAL)
    }
  #
  k <- kmax
  res <- doADF(k)
  while(is.na(res[1])||(abs(res[1])>1.6 & k>0)) {
    k <- k-1
    res <- doADF(k)
    }
  list(statistic=res[1],lag.order=k,p.value=res[2])
  }

# kpss test (internal use only)
kpsst <- function (x,trend,lshort) {
  n <- length(x)
  if(trend==T) {
    t <- 1:n
    e <- residuals(lm(x ~ t))
    table <- c(0.216, 0.176, 0.146, 0.119)
    } else {
    e <- residuals(lm(x ~ 1))
    table <- c(0.739, 0.574, 0.463, 0.347)
    }
  tablep <- c(0.01, 0.025, 0.05, 0.1)
  s <- cumsum(e)
  eta <- sum(s^2)/(n^2)
  s2 <- sum(e^2)/n
  if(lshort==T) {
    l <- trunc(4*(n/100)^0.25)
    } else {
    l <- trunc(12*(n/100)^0.25)
    }
  k <- 0
  for(i in 1:l) {
    ik <- 0
    for(j in (i+1):n) {
      ik <- ik+e[j]*e[j-i]
      }
    k <- k+(1-i/(l+1))*ik
    }
  STAT <- eta/(s2+2*k/n)
  PVAL <- approx(table,tablep,STAT,rule=2)$y
  list(statistic=STAT,lag.order=l,p.value=PVAL)
  }

# function for uniroot test (internal use only)  
doURT <- function(x,g.id,test,glab,par) {
  gruppi <- sort(unique(g.id))
  auxord <- auxstat <- auxp <- nwm <- c()
  h0 <- ifelse(test=="adf","unit root","stationarity")
  auxwarn <- options()$warn
  options(warn=-1)
  for(i in 1:length(gruppi)) {
    auxind <- which(g.id==gruppi[i])
    auxdat <- na.omit(splRecons(x[auxind]))
    nwm[i] <- length(auxdat)
    if(length(auxdat)>4 && var(auxdat)>0) {
      if(test=="adf") {
        auxurt <- adft(auxdat,kmax=par[i])
        } else {
        auxurt <- kpsst(auxdat,trend=par[1],lshort=par[2])  
        }
      auxstat[i] <- auxurt$statistic
      auxp[i] <- auxurt$p.value
      auxord[i] <- auxurt$lag.order
      } else {
      auxstat[i] <- auxp[i] <- auxord[i] <- NA
      }
    }
  if(length(auxstat)>1) names(auxstat) <- names(nwm) <- names(auxord) <- glab
  auxp[which(auxp>1)] <- 1
  auxp[which(auxp<0)] <- 0
  if(length(auxp)==1) {
    res <- list(statistic=auxstat,lag.order=auxord,n=nwm,null=h0,z.value=auxstat,p.value=auxp)
    } else {
    m <- length(auxp)
    logp <- qnorm(auxp)
    rhat <- 1-var(logp)
    rstar <- max(rhat,-1/(m-1))
    auxz <- sum(logp)/sqrt(m*(1+(m-1)*(rstar+0.2*sqrt(2/(m+1))*(1-rstar))))
    #auxz <- sum(logp)/sqrt(m)
    auxpval <- 2*pnorm(-abs(auxz))
    res <- list(statistic=auxstat,lag.order=auxord,n=nwm,null=h0,z.value=auxz,p.value=auxpval)
    }
  options(warn=auxwarn)
  res
  }

# print method for class 'unirootTest'
print.unirootTest <- function(x,...) {
  if(attr(x,"test")=="adf") {
    cat("ADF test (null hypothesis is unit root)","\n")
    } else {
    cat("KPSS test (null hypothesis is stationarity)","\n")
    }
  res <- sapply(x,function(z){z$p.value})
  print(round(res,4))
  }

# estimate parameters of the imputation model (internal use only)
impuFit <- function(xcont,xqual,group,data) {
  res <- G <- list()
  logL <- 0
  z.names <- c(group,xqual)
  for(i in 1:length(xcont)) {
    if(i==1) {
      if(is.null(z.names)) ipar <- "1" else ipar <- z.names  
      } else {
      ipar <- c(xcont[1:(i-1)],z.names)        
      }
    iform <- formula(paste(xcont[i],"~",paste(ipar,collapse="+"),sep=""))
    imod <- lm(iform,data=data)
    imod$call$formula <- iform
    logL <- logL-0.5*AIC(imod,k=0)
    res[[i]] <- imod
    G[[i]] <- ipar
    }
  names(res) <- names(G) <- xcont
  list(hat=res,G=G,logL=logL)
  }

# predict missing values from the imputation model (internal use only)
impuPred <- function(mod,data) {
  est <- mod$hat
  G <- mod$G
  res <- data
  auxwarn <- options()$warn
  options(warn=-1)
  for(i in 1:length(est)) {
    iest <- est[[i]]
    inam <- names(est)[i]
    ipar <- G[[inam]]
    ina <- which(is.na(data[,inam]))
    if(length(ina)>0) {res[ina,inam] <- predict(iest,res[ina,ipar,drop=F])}
    }
  options(warn=auxwarn)
  res
  }

# linear interpolation (internal use only)
linImp <- function(x) {
  res <- x
  auxNA <- which(is.na(x))
  if(length(auxNA)>0) {
    naL <- split(auxNA,cumsum(c(1,diff(auxNA)!=1)))
    for(i in 1:length(naL)) {
      ina <- naL[[i]]
      x1 <- min(ina)-1
      x2 <- max(ina)+1
      y1 <- x[x1]
      y2 <- x[x2]
      b <- (y2-y1)/(x2-x1)
      a <- y1-b*x1
      res[ina] <- a+b*ina
      }
    }
  res
  }

# spline reconstruction (internal use only)
splRecons <- function(x) {
  res <- x
  auxNA <- which(is.na(x))
  if(length(auxNA)>0&length(auxNA)<length(x)) {
    auxO <- which(!is.na(x))
    auxI <- intersect(auxNA,min(auxO):max(auxO))
    yI <- spline(1:length(x),x,xout=1:length(x))
    res[auxI] <- yI$y[auxI]
    }
  res
  }

# lag order determination (internal use only)
arFind <- function(x,group,data) {
  if(!is.null(group)) {
    data[,group] <- factor(data[,group])
    gruppi <- levels(data[,group])    
    res <- rep(0,length(x))
    for(w in 1:length(gruppi)) {
      auxind <- which(data[,group]==gruppi[w])
      for(i in 1:length(x)) {
        ix <- na.omit(splRecons(data[auxind,x[i]]))
        if(var(ix)>0) res[i] <- ar(na.omit(ix))$order
        }
      }    
    } else {
    res <- c()
    for(i in 1:length(x)) {
      ix <- na.omit(splRecons(data[,x[i]]))
      if(var(ix)>0) res[i] <- ar(na.omit(ix))$order
      }
    }
  res
  }

# add lagged instances (internal use only)
addLags <- function(x,group=NULL,data,k=0) {
  if(k>0) {
    if(is.null(group)) {
      for(i in 1:length(x)) {
        for(j in 1:k) {
          data[,paste(x[i],"_lag",j,sep="")] <- c(rep(NA,j),data[1:(nrow(data)-j),x[i]])
          }
        }
      } else {
      gruppi <- levels(factor(data[,group]))
      for(w in 1:length(gruppi)) {
        auxind <- which(data[,group]==gruppi[w])
        for(i in 1:length(x)) {
          for(j in 1:k) {
            data[auxind,paste(x[i],"_lag",j,sep="")] <- c(rep(NA,j),data[1:(length(auxind)-j),x[i]])
            }
          }
        }      
      }
    }
  data
  }

# imputation of missing data (internal use only)
EM.imputation <- function(xcont,xqual,group,data,tol,maxiter,quiet=F) {
  nmiss <- apply(data[,xcont],2,function(v){sum(is.na(v))})
  xcont <- xcont[order(nmiss)]
  currentDat <- data
  for(i in 1:length(xcont)) {
    currentDat[which(is.na(currentDat[,xcont[i]])),xcont[i]] <- mean(data[,xcont[i]],na.rm=T)
    }
  currentFit <- impuFit(xcont=xcont,xqual=xqual,group=group,data=currentDat)
  currentLik <- -Inf
  fine <- forcend <- 0
  count <- 1
  if(quiet==F) {
    cat("Starting EM...")
    flush.console() 
    }
  while(fine==0) {
    newDat <- impuPred(currentFit,data)
    newFit <- impuFit(xcont=xcont,xqual=xqual,group=group,data=newDat)
    newLik <- newFit$logL
    if(quiet==F) {
      cat('\r',"EM iteration ",count,". Log-likelihood: ",newLik,sep="")
      flush.console() 
      }
    if(newLik<currentLik) {
      newLik <- currentLik
      newDat <- currentDat
      #warning("Forced stop of EM algorithm because likelihood has decreased",call.=F)
      fine <- 1
      } else {
      if(newLik-currentLik>tol & count<maxiter) {
        currentFit <- newFit
        currentLik <- newLik
        currentDat <- newDat
        count <- count+1
        } else {
        fine <- 1
        if(count>=maxiter) forcend <- 1
        }
      }
    }
  if(quiet==F) {
    if(forcend==0) {
      cat('\r',"EM converged after ",count," iterations. Log-likelihood: ",newLik,sep="","\n")
      } else {
      cat('\r',"EM stopped after ",maxiter," iterations. Log-likelihood: ",newLik,sep="","\n")      
      }
    }
  newDat
  }

# function to plot a lag shape (internal use only)
makeShape <- function(bmat,maxlag,cumul,bcum,conf,ylim,title) {
  if(!is.null(maxlag)) {
    maxlag <- maxlag-1  #####
    ymlag <- max(as.numeric(rownames(bmat)))
    if(maxlag>=ymlag) {
      if(cumul==F) {
        addmat <- matrix(0,nrow=maxlag-ymlag+1,ncol=3)
        } else {
        addmat <- matrix(bmat[nrow(bmat),],nrow=maxlag-ymlag+1,ncol=3,byrow=T)
        }
      bmat <- rbind(bmat,addmat)
      rownames(bmat) <- -1:(maxlag+1)
      } else {                
      bmat <- bmat[1:(which(rownames(bmat)==as.character(maxlag))+1),]
      }
    auxNZ <- which(bmat[,1]!=0)
    } else {
    auxNZ <- which(bmat[,1]!=0)
    if(cumul==F) {
      if(nrow(bmat)>max(auxNZ)+1) {
        bmat <- bmat[1:(max(auxNZ)+1),]
        }
      } else {                          
      auxCm <- min(intersect(which(diff(bmat[,1])==0)+1,auxNZ))
      if(nrow(bmat)>auxCm) {
        bmat <- bmat[1:auxCm,]
        }
      auxNZ <- which(bmat[,1]!=0)
      }
    }
  xaux <- as.numeric(rownames(bmat))  #####
  xaux <- c(xaux,max(xaux)+1)
  #
  if(is.null(ylim)) {
    upLim <- 1.05*max(bmat)
    lowLim <- 1.05*min(bmat)
    upLim <- max(abs(c(upLim,lowLim)))
    lowLim <- -max(abs(c(upLim,lowLim)))
    } else {
    lowLim <- ylim[1]
    upLim <- ylim[2]
    }
  auxs <- which(bmat[,1]!=0)
  bmat_s <- bmat[c(max(1,min(auxs)-1):min(nrow(bmat),max(auxs)+1)),]
  bval <- as.numeric(rownames(bmat_s))
  #
  m0 <- lm(bmat_s[which(bmat_s[,1]!=0),1]~bval[which(bmat_s[,1]!=0)])
  smooth <- ifelse(sum(residuals(m0)^2)<0.001,F,T)
  #
  if(smooth==T) {
    xgrid <- sort(unique(c(bval,seq(min(bval),max(bval),length=100))))
    ygrid <- cbind(spline(bval,bmat_s[,1],xout=xgrid)$y,spline(bval,bmat_s[,2],xout=xgrid)$y,spline(bval,bmat_s[,3],xout=xgrid)$y)
    } else {
    xgrid <- bval
    ygrid <- bmat_s
    }
  #
  auxsgn <- sign(bmat[auxs,1])
  if(sum(auxsgn==-1)==0) {
    auxdel <- which(ygrid[,1]<0)
    } else if(sum(auxsgn==1)==0) {
    auxdel <- which(ygrid[,1]>0)
    } else {
    auxdel <- c()  
    }
  if(length(auxdel)>0) {
    auxInt <- c()
    for(i in 1:length(bval)) {
      auxInt[i] <- which(xgrid==bval[i])
      }
    for(i in 1:length(auxdel)) {
      isx <- auxInt[max(which(auxInt<=auxdel[i]))]
      idx <- auxInt[min(which(auxInt>=auxdel[i]))]
      ygrid[isx:idx,] <- NA
      }
    ygrid[c(1,nrow(ygrid)),] <- 0
    for(j in 1:ncol(ygrid)) {
      ygrid[,j] <- linImp(ygrid[,j])
      }
    }
  #
  plot(0,type="n",xlim=c(min(xaux),max(xaux)),ylim=c(lowLim,upLim),yaxs="i",xaxs="i",cex.lab=1.2,
    lwd=2,xaxt="n",yaxt="n",xlab="Lag",ylab="Coefficient",main=title,cex.main=1.2) 
  if(cumul==T) mtext("cumulative lag shape",cex=0.9)
  polygon(c(xgrid,rev(xgrid)),c(ygrid[,2],rev(ygrid[,3])),border=NA,col="grey80")
  yaxaux <- seq(lowLim,upLim,length=21)
  ylabaux <- signif(yaxaux,3)
  ylabaux[11] <- 0
  xaxaux <- seq(min(xaux),max(xaux))
  auxby <- max(1,round((max(xaux)-min(xaux)+1)/30))
  xlabaux1 <- xlabaux2 <- seq(min(xaux),max(xaux),by=auxby)
  xlabaux2[c(1,length(xlabaux1))] <- NA
  abline(h=yaxaux,v=seq(min(xaux),max(xaux),by=auxby),col="grey75",lty=2)
  abline(h=0,lty=2,col="grey35")                                        
  lines(ygrid[,1]~xgrid,col="grey40",lty=2)
  #
  auxpoi <- which(bmat[,1]!=0)
  auxpoiOK <- max(1,(min(auxpoi)-1)):min(nrow(bmat),(max(auxpoi)+1))
  points(bmat[auxpoiOK,1]~as.numeric(names(bmat[,1]))[auxpoiOK],col="grey35",lty=2,cex=0.6)
  #
  axis(1,at=xlabaux1,labels=xlabaux2,cex.axis=1.1)
  axis(2,at=yaxaux,labels=ylabaux,cex.axis=1.1)
  confLeg <- paste("   ",conf*100,"% CI: (",bcum[2],", ",bcum[3],")",sep="")      
  if(max(bmat[,1])>0) {
    legpos <- "bottomright"
    } else {
    legpos <- "topright"
    }                                       
  est <- bmat[,1]
  if(cumul==T) {
    newest <- est[1]
    for(i in 2:length(est)) {
      newest[i] <- est[i]-est[i-1]
      } 
    est <- newest
    }
  minlag <- min(as.numeric(rownames(bmat)[which(est!=0)]))
  maxlag <- max(as.numeric(rownames(bmat)[which(est!=0)]))             
  legend(legpos,legend=c(paste("Relevant lags: ",minlag," to ",maxlag,sep=""),paste("Cumulative coefficient: ",bcum[1],sep=""),confLeg),bty="n",cex=1.1)
  box()
  }

# estimated lag shapes
lagShapes <- function(x,cumul=FALSE) {
  if(("dlsem" %in% class(x))==F) stop("Argument 'x' must be an object of class 'dlsem'",call.=F)
  if(length(cumul)!=1 || !is.logical(cumul)) stop("Argument 'cumul' must be a logical value",call.=F)
  G <- makeGraph(x)$full.graph
  est <- x$estimate
  nomi <- names(x$estimate)
  pset <- getStruct(x)
  res <- vector("list",length=length(nomi))
  names(res) <- nomi
  for(i in 1:length(nomi)) {
    ires <- list()
    ipar <- pset[[nomi[i]]]
    if(length(ipar)>0) {
      for(j in 1:length(ipar)) {
        ijtab <- lagEff(est[[nomi[i]]],x=ipar[j],lag=NULL,cumul=cumul)
        if(cumul==F) {
          ijtabOK <- rbind(ijtab,c(0,0))
          } else {
          ijtabOK <- rbind(ijtab,ijtab[nrow(ijtab),])
          }
        rownames(ijtabOK)[nrow(ijtabOK)] <- nrow(ijtabOK)-1
        ires[[j]] <- ijtabOK
        }
      names(ires) <- ipar
      }
    res[[i]] <- ires
    }
  res
  }

# plot the lag shape associated to an overall causal effect or a path
lagPlot <- function(x,from=NULL,to=NULL,path=NULL,maxlag=NULL,cumul=FALSE,conf=0.95,use.ns=FALSE,ylim=NULL,title=NULL) {
  if(("dlsem" %in% class(x))==F) stop("Argument 'x' must be an object of class 'dlsem'",call.=F)
  if(!is.null(maxlag) && (length(maxlag)!=1 || !is.numeric(maxlag) || maxlag<=0 || maxlag!=round(maxlag))) stop("Argument 'maxlag' must be a positive integer number",call.=F)
  if(length(cumul)!=1 || !is.logical(cumul)) stop("Argument 'cumul' must be a logical value",call.=F)
  if(length(conf)!=1 || !is.numeric(conf) || conf<=0 || conf>=1) stop("Argument 'conf' must be a real number in the interval (0,1)",call.=F)
  if(length(use.ns)!=1 || !is.logical(use.ns)) stop("Argument 'use.ns' must be a logical value",call.=F)
  if(!is.null(ylim) && (length(ylim)!=2 || ylim[1]>=ylim[2])) stop("Invalid argument 'ylim'",call.=F)
  #
  if(!is.null(from) && !is.na(from) && (is.null(to) & is.null(path))) {
    path <- from
    auxstr <- strsplit(path,"\\*")[[1]]
    if(length(auxstr)<2) {
      stop("Argument 'to' is missing",call.=F)
      } else {  
      from <- NULL
      }
    } else {
    if(is.null(from)||is.null(to)) {
      auxstr <- strsplit(path,"\\*")[[1]]
      if(length(auxstr)<2) stop("Invalid path length",call.=F)
      from <- to <- NULL
      } else {
      path <- NULL
      }
    }
  #
  if(is.null(path) && (is.null(from) || is.na(from))) stop("Argument 'from' is missing",call.=F)
  if(is.null(path) && (is.null(to) || is.na(to))) stop("Argument 'to' is missing",call.=F)
  xedgF <- edgeMat(x,conf=conf,full=T)
  xedg <- edgeMat(x,conf=conf,full=F)
  if(is.null(path)) {                        
    auxpa <- causalEff(x,from=from,to=to,lag=NULL,cumul=cumul,conf=conf,use.ns=use.ns)$overall
    } else {
    path <- gsub(" ","",path)
    auxstr <- strsplit(path,"\\*")[[1]]
    pathchk <- setdiff(auxstr,names(x$estimate))
    if(length(pathchk)>0) stop("Unknown variable '",pathchk[1],"' in the path",call.=F)
    from <- auxstr[1]
    to <- rev(auxstr)[1] 
    isIn <- isInF <- 1
    for(i in 2:length(auxstr)) {
      isIn <- isIn*length(which(xedg[,1]==auxstr[i-1] & xedg[,2]==auxstr[i]))
      isInF <- isInF*length(which(xedgF[,1]==auxstr[i-1] & xedgF[,2]==auxstr[i]))
      }
    if((use.ns==T & isInF>0) | isIn>0) {
      auxpa <- causalEff(x,from=from,to=to,cumul=cumul,conf=conf,use.ns=use.ns)
      auxpa <- auxpa[[path]] 
      } else {
      #if(isInF>0) {
      #  stop("Path not found. Try to reduce 'conf' or to set 'use.ns' to TRUE",call.=F)
      #  } else {
      #  stop("Inexistent path",call.=F)
      #  }
      auxpa <- NULL
      }
    }
  if(!is.null(auxpa)) {
    yaux <- rbind(rep(0,ncol(auxpa)),auxpa)
    rownames(yaux) <- c(-1:(nrow(yaux)-2))
    if(is.null(title)) {
      if(is.null(path)) {
        title <- paste(to," ~ ",from,sep="")
        } else {
        title <- paste(auxstr,collapse=" * ")
        }
      }
    bmat <- yaux[,c(1,3,4)]
    if(cumul==F) {
      auxbcum <- causalEff(x,from=from,to=to,cumul=T,conf=conf,use.ns=use.ns)$overall
      bcum <- signif(auxbcum[nrow(auxbcum),c(1,3,4)],5)
      } else {
      bcum <- signif(bmat[nrow(bmat),],5)
      }
    makeShape(bmat=bmat,maxlag=maxlag,cumul=cumul,bcum=bcum,conf=conf,ylim=ylim,title=title)
    } else {
    NULL  
    }
  }

# check if integer (internal use only)
intCheck <- function(x) {
  if(is.numeric(x) && x[1]==round(x[1]) && x[1]>=0) T else F
  }

# adjust global control options (internal use only)
gconAdj <- function(x) {
  nomi <- names(x)
  warn1 <- warn2 <- 0
  if(!is.null(x) && !is.list(x)) warn2 <- warn2+1
  unknam <- setdiff(nomi,c("adapt","min.gestation","max.gestation","min.width","max.lead","sign"))
  if(length(unknam)>0) {
    x <- x[setdiff(nomi,unknam)]
    warning("Some components with unknown names in argument 'global.control' were ignored",call.=F)
    }
  #
  if("adapt" %in% nomi) {
    xad <- x$adapt
    if(is.null(xad)) xad <- F
    if(is.logical(xad)) {
      if(length(xad)>1) warn1 <- warn1+1
      xad <- xad[1]
      } else {
      xad <- F
      warn2 <- warn2+1
      }
    } else {
    xad <- F
    }
  #
  if("min.gestation" %in% nomi) {
    xming <- x$min.gestation
    if(is.null(xming)) xming <- 0
    if(length(xming)>1) warn1 <- warn1+1
    xming <- xming[1]
    if(intCheck(xming)==F) {
      xming <- 0
      warn2 <- warn2+1
      }
    } else {
    xming <- 0
    }
  #
  if("max.gestation" %in% nomi) {
    xmaxg <- x$max.gestation
    if(is.null(xmaxg)) xmaxg <- Inf
    if(length(xmaxg)>1 | xmaxg<xming) warn1 <- warn1+1
    xmaxg <- max(xming,xmaxg[1])
    if(intCheck(xmaxg)==F) {
      xmaxg <- Inf
      warn2 <- warn2+1
      }
    } else {
    xmaxg <- Inf
    }
  #
  if("max.lead" %in% nomi) {
    xml <- x$max.lead
    if(is.null(xml)) xml <- Inf
    if(length(xml)>1) warn1 <- warn1+1
    xml <- xml[1]
    if(intCheck(xml)==F) {
      xml <- Inf
      warn2 <- warn2+1
      }
    } else {
    xml <- Inf
    }
  if(xml<Inf & xmaxg==Inf) xmaxg <- xml
  #
  if("min.width" %in% nomi) {
    xminw <- x$min.width
    if(is.null(xminw)) xminw <- 0
    if(length(xminw)>1) warn1 <- warn1+1
    xminw <- xminw[1]
    if(intCheck(xminw)==F) {
      xminw <- 0
      warn2 <- warn2+1
      }
    } else {
    xminw <- 0
    }
  #
  if("sign" %in% nomi) {
    xsg <- x$sign
    if(is.null(xsg)) xsg <- F
    if(xsg[1] %in% c("+","-")) {
      if(length(xsg)>1) warn1 <- warn1+1
      xsg <- xsg[1]
      } else {
      xsg <- F
      warn2 <- warn2+1
      }
    } else {
    xsg <- F
    }
  #
  if(xming>xml) {
    warn2 <- warn2+1
    xming <- xml
    }
  if(xmaxg>xml) {
    warn2 <- warn2+1
    xmaxg <- xml
    }
  if(xminw>xml) {
    warn2 <- warn2+1
    xminw <- xml
    }
  #
  if(warn1>0) warning("Some components in argument 'global.control' had length >1 and only the first element was used",call.=F)
  if(warn2>0) warning("Some invalid or uncoherent components in argument 'global.control' were adjusted for coherence",call.=F)
  list(adapt=xad,min.gestation=xming,max.gestation=xmaxg,min.width=xminw,max.lead=xml,sign=xsg)
  }

# adjust local control options (internal use only)
lconAdj <- function(x,gcon,pset) {
  warn1 <- warn2 <- 0
  nomi <- names(pset)
  if(!is.null(x) && !is.list(x)) warn2 <- warn2+1
  xad <- vector("list",length=length(nomi))
  names(xad) <- nomi
  for(j in 1:length(nomi)) {
    if("adapt" %in% names(x) && nomi[j] %in% names(x[["adapt"]])) {
      xad[[j]] <- x$adapt[[nomi[j]]]
      if(is.null(xad[[j]])) xad[[j]] <- gcon$adapt
      if(is.logical(xad[[j]])) {
        if(length(xad[[j]])>1) warn1 <- warn1+1
        xad[[j]] <- xad[[j]][1]
        } else {
        xad[[j]] <- gcon$adapt
        warn2 <- warn2+1
        }
      } else {
      xad[[j]] <- gcon$adapt
      }
    }
  #
  comp <- c("min.gestation","max.gestation","min.width","max.lead","sign")
  res <- vector("list",length=length(comp))
  names(res) <- comp
  for(i in 1:length(comp)) {
    ires <- vector("list",length=length(nomi))
    names(ires) <- nomi
    for(j in 1:length(nomi)) {
      ijpset <- pset[[nomi[[j]]]]
      ijres <- c()
      if(length(ijpset)>0) {
        for(k in 1:length(ijpset)) {
          if(comp[i] %in% names(x) && nomi[j] %in% names(x[[comp[i]]])) {
            ijval <- x[[comp[i]]][[nomi[[j]]]]
            if(ijpset[k] %in% names(ijval)) {
              ijres[k] <- ijval[ijpset[k]]
              if(is.null(ijres[k])) ijres[k] <- gcon[[comp[i]]]
              if(comp[i]=="sign") {
                ijkch <- ijres[k] %in% c("+","-")
                } else {
                ijkch <- intCheck(ijres[k]) 
                }
              if(ijkch==F) {
                ijres[k] <- gcon[[comp[i]]]
                warn2 <- warn2+1
                }
              } else {
              ijres[k] <- gcon[[comp[i]]]  #####
              }
            } else {
            ijres[k] <- gcon[[comp[i]]]
            }
          }
        names(ijres) <- ijpset
        ires[[j]] <- ijres
        }
      }
    res[[i]] <- ires
    }
  if(warn1>0) warning("Some components in argument 'local.control' had length >1 and only the first element was used",call.=F)
  if(warn2>0) warning("Some invalid or uncoherent components in argument 'local.control' were adjusted for coherence",call.=F)
  c(adapt=list(xad),res)
  }

# check missing values (internal use only)
checkNA <- function(x,group,data) {
  if(sum(!is.na(data[,x]))<3) stop("Variable '",x,"' has less than 3 observed values",call.=F)
  if(!is.null(group)) {
    gruppi <- levels(factor(data[,group]))
    for(i in 1:length(gruppi)) {
      auxind <- which(data[,group]==gruppi[i])
      if(sum(!is.na(data[auxind,x]))<1) {
        stop("Variable '",x,"' has no observed values in group '",gruppi[i],"'",call.=F)
        }
      }
    }
  }

# check variable name (internal use only)
checkName <- function(x) {
  res <- T
  if((substr(x,1,1) %in% c(letters,toupper(letters)))==F) res <- F
  if(length(grep("[-\\+\\*\\/]",x))>0) res <- F
  res
  }

# find operator in a variable name (internal use only)
findOp <- function(x) {
  if(gregexpr("\\(",x)[[1]][1]>0) {
    strsplit(x,"\\(")[[1]][1]
    } else {
    NULL 
    }
  }

# automated full model code
autoCode <- function(var.names,lag.type="ecq") {
  if(missing(var.names)==F & length(var.names)<2) stop("Argument 'var.names' must be at least of length 2",call.=F)
  if(length(lag.type)!=1) stop("Argument 'lag.type' must be of length 1",call.=F)
  if((lag.type %in% c("ecq","qd","ld","gam"))==F) stop("Argument 'lag.type' must be one among 'ecq', 'qd', 'ld' and 'gam'",call.=F)
  res <- list()
  if(checkName(var.names[1])==F) stop("'",var.names[1]," 'is not a valid variable name",call.=F)
  res[[1]] <- formula(paste(var.names[1],"~1",sep=""))
  for(i in 2:length(var.names)) {
    iy <- var.names[i]
    if(checkName(iy)==F) stop("'",iy," 'is not a valid variable name",call.=F)
    iX <- var.names[1:(i-1)]
    res[[i]] <- formula(paste(iy,"~",paste(paste(lag.type,"(",iX,",,)",sep=""),collapse="+"),sep=""))    
    }
  res
  }

# adjust diff options (internal use only)
diffoptAdj <- function(x) {
  warn1 <- warn2 <- 0
  if(!is.null(x) && !is.list(x)) warn2 <- warn2+1
  nomi <- names(x)
  unknam <- setdiff(nomi,c("test","maxdiff","ndiff"))
  if(length(unknam)>0) {
    x <- x[setdiff(nomi,unknam)]
    warning("Some components with unknown names in argument 'diff.options' were ignored",call.=F)
    }
  # 
  if("test" %in% nomi) {
    test <- x$test
    if(is.null(test)) test <- NULL
    if(length(test)>1) warn1 <- warn1+1
    test <- test[1]
    if(!is.null(test) && (test %in% c("kpss","adf"))==F) {
      test <- NULL
      warn2 <- warn2+1
      }
    } else {
    test <- NULL  
    }
  #
  if("maxdiff" %in% nomi) {
    maxdiff <- x$maxdiff
    if(is.null(maxdiff)) maxdiff <- 2
    if(length(maxdiff)>1) warn1 <- warn1+1
    maxdiff <- maxdiff[1]
    if(intCheck(maxdiff)==F) {
      maxdiff <- 2
      warn2 <- warn2+1
      }
    } else {
    maxdiff <- 2
    }
  #
  if("ndiff" %in% nomi) {
    ndiff <- x$ndiff
    if(length(ndiff)>1) warn1 <- warn1+1
    ndiff <- ndiff[1]
    if(!is.null(ndiff) && intCheck(ndiff)==F) {
      ndiff <- NULL
      warn2 <- warn2+1
      }
    } else {
    ndiff <- NULL
    }
  #
  if(warn1>0) warning("Some components in argument 'diff.options' had length >1 and only the first element was used",call.=F)
  if(warn2>0) warning("Some invalid or uncoherent components in argument 'diff.options' were adjusted",call.=F)
  list(test=test,maxdiff=maxdiff,ndiff=ndiff)
  }

# adjust imput options (internal use only)
impoptAdj <- function(x) {
  warn1 <- warn2 <- 0
  if(!is.null(x) && !is.list(x)) warn2 <- warn2+1
  nomi <- names(x)
  unknam <- setdiff(nomi,c("tol","maxiter","maxlag","no.imput"))
  if(length(unknam)>0) {
    x <- x[setdiff(nomi,unknam)]
    warning("Some components with unknown names in argument 'imput.options' were ignored",call.=F)
    }
  # 
  if("no.imput" %in% nomi) {
    noimp <- x$no.imput
    } else {
    noimp <- NULL  
    }
  #
  if("tol" %in% nomi) {
    tol <- x$tol
    if(is.null(tol)) tol <- 0.0001
    if(length(tol)>1) warn1 <- warn1+1
    tol <- tol[1]
    if(tol<=0) {
      tol <- 0.0001
      warn2 <- warn2+1
      }
    } else {
    tol <- 0.0001
    }
  #
  if("maxiter" %in% nomi) {
    maxiter <- x$maxiter
    if(is.null(maxiter)) maxiter <- 500
    if(length(maxiter)>1) warn1 <- warn1+1
    maxiter <- maxiter[1]
    if(intCheck(maxiter)==F) {
      maxiter <- 500
      warn2 <- warn2+1
      }
    } else {
    maxiter <- 500
    }
  #
  if("maxlag" %in% nomi) {
    maxlag <- x$maxlag
    if(is.null(maxlag)) maxlag <- 2
    if(length(maxlag)>1) warn1 <- warn1+1
    maxlag <- maxlag[1]
    if(intCheck(maxlag)==F) {
      maxlag <- 2
      warn2 <- warn2+1
      }
    } else {
    maxlag <- 2
    }
  #
  if(warn1>0) warning("Some components in argument 'imput.options' had length >1 and only the first element was used",call.=F)
  if(warn2>0) warning("Some invalid components in argument 'imput.options' were adjusted",call.=F)
  list(tol=tol,maxiter=maxiter,maxlag=maxlag,no.imput=noimp)
  }

# preprocessing
preProcess <- function(x=NULL,group=NULL,time=NULL,seas=NULL,data,log=FALSE,
  diff.options=list(test=NULL,maxdiff=2,ndiff=NULL),
  imput.options=list(tol=0.0001,maxiter=500,maxlag=2,no.imput=NULL),quiet=FALSE) {
  #
  if(missing(data)==F & !identical(class(data),"data.frame")) stop("Argument 'data' must be a data.frame",call.=F)  
  if(!is.null(group) && length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
  if(!is.null(group) && is.na(group)) group <- NULL
  if(!is.null(seas) && (length(seas)!=1 || !is.numeric(seas) || seas<=0 || seas!=round(seas))) stop("Argument 'seas' must be a positive integer value",call.=F)
  if(!is.null(group)) {
    if(length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
    if((group %in% colnames(data))==F) stop("Variable '",group,"' provided to argument 'group' not found in data",call.=F)
    gruppi <- levels(factor(data[,group]))
    if(length(gruppi)<2) stop("The group factor must have at least 2 unique values",call.=F)
    data[,group] <- factor(data[,group])
    if(min(table(data[,group]))<3) stop("There must be at least 3 observations per group",call.=F)
    data <- data[order(data[,group]),]
    } else {
    if(nrow(data)<3) stop("There must be at least 3 observations",call.=F)  
    }
  if(!is.null(time)) {
    if(is.na(time)) time <- NULL
    if(length(time)!=1) stop("Argument 'time' must be of length 1",call.=F)
    if((time %in% colnames(data))==F) stop("Variable '",time,"' provided to argument 'time' not found in data",call.=F)
    if(time %in% group) stop("Variable '",time,"' is provided to both arguments 'group' and 'time'",call.=F)
    if(isTimeVar(data[,time])==F) stop("The time variable is neither numeric nor a date",call.=F)
    if(is.null(group)) {
      if(sum(duplicated(data[,time]))>0) stop("The time variable has duplicated values",call.=F)
      } else {
      timesplit <- split(data[,time],data[,group])
      if(sum(sapply(timesplit,function(z){sum(duplicated(z))}))>0) stop("The time variable has duplicated values",call.=F)  
      }
    #
    if(!is.null(group)) {
      for(i in 1:length(gruppi)) {
        iind <- which(data[,group]==gruppi[i])
        idat <- data[iind,]
        data[iind,] <- idat[order(idat[,time]),]
        }
      } else {
      data <- data[order(data[,time]),]
      }
    #
    }
  if(!is.null(group) && length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
  if(!is.null(group) && (group %in% colnames(data))==F) stop("Variable '",group,"' not found in data",sep="",call.=F)
  if(length(quiet)!=1 || !is.logical(quiet)) stop("Argument 'quiet' must be a logical value",call.=F)
  #if(!is.null(group)) {
  #  data[,group] <- factor(data[,group])
  #  gruppi <- levels(data[,group])
  #  }
  if(!is.null(x)) {
    x2del <- c()
    for(i in 1:length(x)) {
      if((x[i] %in% colnames(data))==F) {
        if(quiet==F) warning("Variable '",x[i],"' not found in data and ignored",call.=F)
        x2del <- c(x2del,x[i])
        }
      }
    x <- setdiff(x,x2del)
    if(length(x)<1) stop("No valid variables provided to argument 'x'",call.=F)
    } else {
    x <- setdiff(colnames(data),c(group,time))
    }
  if(length(quiet)!=1 || !is.logical(quiet)) stop("Argument 'quiet' must be a logical value",call.=F)
  #
  diff.options <- diffoptAdj(diff.options)
  if(is.null(diff.options$test)) {
    if(is.null(group)) {
      n <- nrow(data)
      } else {
      n <- min(table(data[,group]))  
      }
    diff.options$test <- ifelse(n>100,"adf","kpss")
    }
  ndiff <- diff.options$ndiff
  imput.options <- impoptAdj(imput.options)
  #
  if(!is.null(time)) {
    if(!is.null(group)) {
      for(i in 1:length(gruppi)) {
        iind <- which(data[,group]==gruppi[i])
        idat <- data[iind,]
        data[iind,] <- idat[order(idat[,time]),]
        }
      } else {
      data <- data[order(data[,time]),]
      }
    }
  nodenam <- x
  xfact <- c()
  for(i in 1:length(nodenam)) {
    if(isQuant(data[,nodenam[i]])==F) {
      if(sum(is.na(data[,nodenam[i]]))>0) stop("Variable ",nodenam[i]," is qualitative and contains missing values",sep="",call.=F)
      xfact <- c(xfact,nodenam[i])
      data[,nodenam[i]] <- factor(data[,nodenam[i]])
      }
    }
  # deseasonalization
  if(!is.null(seas) && seas>1) data[,nodenam] <- deSeas(setdiff(nodenam,xfact),seas=seas,group=group,data=data)
  # log transformation
  logOK <- c()
  if(identical(log,T)) {
    logtest <- setdiff(nodenam,xfact)            
    if(length(logtest)>0) {
      for(i in 1:length(logtest)) {
        if(sum(data[,logtest[i]]<=0,na.rm=T)>0) {
          if(quiet==F) warning("Logarithmic transformation not applied to variable '",logtest[i],"'",call.=F)      
          } else {
          data[,logtest[i]] <- log(data[,logtest[i]])
          logOK <- c(logOK,logtest[i])
          }
        }
      }
    } else if(!identical(log,F)) {
    for(i in 1:length(log)) {
      if((log[i] %in% colnames(data))==F) {
        if(quiet==F) warning("Variable '",log[i],"' provided to argument 'log' not found in data",call.=F)
        } else if(log[i] %in% c(group,time)) {
        if(quiet==F) warning("Logarithmic transformation not applied to variable '",log[i],"'",call.=F)
        } else {
        if(sum(data[,log[i]]<=0,na.rm=T)>0) {
          if(quiet==F) warning("Logarithmic transformation not applied to variable '",log[i],"'",call.=F)
          } else {
          data[,log[i]] <- log(data[,log[i]])  
          }
        }
      }
    }
  # imputation
  auxna <- apply(data[,nodenam],1,function(x){sum(is.na(x))})
  if(sum(auxna)>0) {
    auxOK <- unname(which(auxna<length(nodenam))) 
    if(sum(is.na(data[auxOK,]))>0 && imput.options$maxiter>0) {
      #
      auxch <- setdiff(imput.options$no.imput,colnames(data))
      if(length(auxch)>0) warning("Variable '",auxch[1],"' provided to component 'no.imput' in argument 'imput.options' not found in data and ignored",call.=F)
      #
      x2imp <- setdiff(nodenam,c(imput.options$no.imput,xfact))
      if(length(x2imp)>0) {
        for(i in 1:length(x2imp)) {
          nachk <- checkNA(x2imp[i],group,data[auxOK,])
          }
        nIm <- sum(is.na(data[,x2imp]))
        if(nIm>0) {
          emdat <- data[,c(group,xfact,x2imp),drop=F]
          for(i in 1:length(x2imp)) {
            iord <- min(arFind(x2imp[i],group,data),imput.options$maxlag)
            if(iord>0) {
              idx <- addLags(x=x2imp[i],data=data,k=iord)
              emdat <- cbind(emdat,idx) 
              }
            }
          filldat <- EM.imputation(xcont=x2imp,xqual=xfact,group=group,data=emdat[auxOK,],tol=imput.options$tol,maxiter=imput.options$maxiter,quiet=quiet)
          data[auxOK,c(group,xfact,x2imp)] <- filldat[,c(group,xfact,x2imp)]
          }
        }
      } else {
      if(quiet==F && imput.options$maxiter==0) cat("Imputation not performed","\n")  
      }
    }
  # differentiation
  difftest <- setdiff(nodenam,xfact)
  if(is.null(ndiff)) {
    ndiff <- 0
    if(length(difftest)>0 & diff.options$maxdiff>0) {
      fine <- 0
      if(quiet==F) cat("Checking stationarity ...")
      flush.console()
      data <- applyDiff(x=difftest,group=group,data=data,k=rep(0,length(difftest)))
      while(fine==0) {
        auxp <- c()
        urtList <- urtFun(x=difftest,group=group,time=NULL,data=data,test=diff.options$test,log=F)
        for(i in 1:length(difftest)) {
          ipvl <- urtList[[i]]$p.value
          if(is.null(ipvl)) {
            auxp[i] <- 0
            } else {
            if(is.na(ipvl)) {
              auxp[i] <- 0
              } else {
              auxp[i] <- ipvl
              }
            }
          }
        if(diff.options$test=="adf") {
          nUR <- length(which(auxp>0.05))
          } else {
          nUR <- length(which(auxp<0.05))
          }
        if(nUR>0) {
          if(ndiff<diff.options$maxdiff) {
            ndiff <- ndiff+1
            data <- applyDiff(x=difftest,group=group,data=data,k=rep(1,length(difftest)))                                        
            } else {
            fine <- 1
            }
          } else {
          fine <- 1
          }
        }
      } else {
      urtList <- NULL
      }
    } else {
    data <- applyDiff(x=difftest,group=group,data=data,k=rep(ndiff,length(difftest)))
    }
  if(quiet==F) {
    cat('\r')
    if(ndiff>0) {
      cat("Order",ndiff,"differentiation performed","\n")
      } else {
      cat("Differentiation not performed","\n")
      }
    }
  res <- data[,c(group,time,x),drop=F]
  attr(res,"log") <- logOK
  attr(res,"ndiff") <- ndiff
  res
  }

# check collinearity (internal use only)
collCheck <- function(y,x,group,data) {
  form <- paste(y,"~",paste(c(group,x),collapse="+"),sep="")
  m0 <- lm(formula(form),data=data)
  names(which(is.na(m0$coefficients)))
  }

# fit a dlsem
dlsem <- function(model.code,group=NULL,time=NULL,exogenous=NULL,data,
  hac=TRUE,gamma.by=0.05,global.control=NULL,local.control=NULL,log=FALSE,
  diff.options=list(test=NULL,maxdiff=2,ndiff=NULL),
  imput.options=list(tol=0.0001,maxiter=500,maxlag=2,no.imput=NULL),quiet=FALSE) {
  #
  if(missing(model.code)==F & (!is.list(model.code) || length(model.code)==0 || sum(sapply(model.code,class)!="formula")>0)) stop("Argument 'model code' must be a list of formulas",call.=F)
  if(missing(data)==F & !identical(class(data),"data.frame")) stop("Argument 'data' must be a data.frame",call.=F)
  #nameOfData <- deparse(substitute(data))
  if(!is.null(group) && length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
  if(!is.null(group) && is.na(group)) group <- NULL
  if(!is.null(group)) {
    if(length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
    if((group %in% colnames(data))==F) stop("Variable '",group,"' provided to argument 'group' not found in data",call.=F)
    if(group %in% exogenous) stop("Variable '",group,"' is provided to both arguments 'group' and 'exogenous'",call.=F)
    gruppi <- levels(factor(data[,group]))
    if(length(gruppi)<2) stop("The group factor must have at least 2 unique values",call.=F)
    data[,group] <- factor(data[,group])
    if(min(table(data[,group]))<3) stop("There must be at least 3 observations per group",call.=F)
    data <- data[order(data[,group]),]
    } else {
    if(nrow(data)<3) stop("There must be at least 3 observations",call.=F)  
    }
  if(!is.null(time)) {
    if(is.na(time)) time <- NULL
    if(length(time)!=1) stop("Argument 'time' must be of length 1",call.=F)
    if((time %in% colnames(data))==F) stop("Variable '",time,"' provided to argument 'time' not found in data",call.=F)
    if(time %in% group) stop("Variable '",time,"' is provided to both arguments 'group' and 'time'",call.=F)
    if(time %in% exogenous) stop("Variable '",time,"' is provided to both arguments 'time' and 'exogenous'",call.=F)
    if(isTimeVar(data[,time])==F) stop("The time variable is neither numeric nor a date",call.=F)
    if(is.null(group)) {
      if(sum(duplicated(data[,time]))>0) stop("The time variable has duplicated values",call.=F)
      } else {
      timesplit <- split(data[,time],data[,group])
      if(sum(sapply(timesplit,function(z){sum(duplicated(z))}))>0) stop("The time variable has duplicated values",call.=F)  
      }
    #
    if(!is.null(group)) {
      for(i in 1:length(gruppi)) {
        iind <- which(data[,group]==gruppi[i])
        idat <- data[iind,]
        data[iind,] <- idat[order(idat[,time]),]
        }
      } else {
      data <- data[order(data[,time]),]
      }
    }
  if(!is.null(exogenous) && identical(NA,exogenous)) exogenous <- NULL
  if(!is.null(group) && length(group)!=1) stop("Argument 'group' must be of length 1",call.=F)
  if(!is.null(group) && (group %in% colnames(data))==F) stop("Variable '",group,"' not found in data",sep="",call.=F)
  if(length(hac)!=1 || !is.logical(hac)) stop("Argument 'hac' must be a logical value",call.=F)
  if(length(gamma.by)!=1 || !is.numeric(gamma.by) || gamma.by<=0 || gamma.by>=1) stop("Argument 'gamma.by' must be a real number in the interval (0,1)",call.=F)
  if(length(quiet)!=1 || !is.logical(quiet)) stop("Argument 'quiet' must be a logical value",call.=F)
  rownames(data) <- 1:nrow(data)
  estL <- pset <- list()
  for(i in 1:length(model.code)) {
    if(sum(grepl("\\-",model.code[[i]]))>0) stop("Invalid character '-' in 'model.code', regression of '",model.code[[i]][2],"'",call.=F)
    if(sum(grepl("\\:",model.code[[i]]))>0) stop("Invalid character ':' in 'model.code', regression of '",model.code[[i]][2],"'",call.=F) #####
    if(sum(grepl("\\*",model.code[[i]]))>0) stop("Invalid character '*' in 'model.code', regression of '",model.code[[i]][2],"'",call.=F) #####
    }
  for(i in 1:length(model.code)) {
    iscan <- scanForm(model.code[[i]],warn=T)
    ynam <- iscan$y    
    if(checkName(ynam)==F) stop("'",ynam,"' is not a valid variable name",call.=F)
    #
    if((ynam %in% colnames(data))==F) {
      auxop <- findOp(ynam)
      if(is.null(auxop)) {
        stop("Variable '",ynam,"' not found in data",call.=F)
        } else {
        stop("Operator ",auxop,"() not allowed in 'model.code'",call.=F) 
        }
      }
    #
    ipar <- names(iscan$ltype)
    if(length(ipar)==0) ipar <- character(0)
    pset[[i]] <- ipar
    names(model.code)[i] <- names(pset)[i] <- ynam
    if(length(ipar)>0) {
      if(!is.null(group)) {
        if(group %in% ipar) stop("Variable '",group,"' is defined as a group factor and appears in 'model.code'",call.=F) 
        }
      #
      auxfun <- setdiff(ipar,colnames(data))
      if(length(auxfun)>0) {
        auxop <- findOp(auxfun[1])
        if(is.null(auxop)) {
          stop("Variable '",auxfun[1]," not found in data",call.=F)
          } else {
          stop("Operator ",auxop,"() not allowed in 'model.code'",call.=F) 
          }
        }
      #
      auxexo <- intersect(ipar,exogenous)
      if(length(auxexo)>0) stop("Variable '",auxexo[1],"' appears both in 'model.code' and in 'exogenous'",call.=F)
      auxdupl <- duplicated(ipar)
      if(sum(auxdupl)>0) stop("Duplicated covariate '",ipar[auxdupl][1],"' in 'model.code', regression of '",model.code[[i]][2],"'",call.=F)
      }
    auxcoll <- collCheck(ynam,c(ipar,exogenous),group,data)
    if(length(auxcoll)>0) stop("Collinearity problem with covariate '",auxcoll[1],"' in the regression of '",model.code[[i]][2],"'",call.=F)
    }
  auxdupl <- duplicated(names(model.code))
  if(sum(auxdupl)>0) stop("Duplicated response variable '",names(model.code)[auxdupl][1],"' in 'model.code'",call.=F)
  nodenam <- unique(c(names(pset),unlist(pset)))
  auxadd <- setdiff(nodenam,names(model.code))
  if(length(auxadd)>0) {
    for(i in length(auxadd):1) {
      model.code <- c(formula(paste(auxadd[i],"~1",sep="")),model.code)
      names(model.code)[1] <- auxadd[i]
      pset[[length(pset)+1]] <- character(0) 
      names(pset)[length(pset)] <- auxadd[i]
      }
    }
  if(length(nodenam)<2) stop("The model cannot contain less than 2 endogenous variables",call.=F)
  auxvar <- setdiff(c(nodenam,exogenous),colnames(data))
  if(length(auxvar)>0) {
    stop("Variable '",auxvar[1],"' not found in data",sep="",call.=F)
    }
  for(i in 1:length(nodenam)) {
    if(isQuant(data[,nodenam[i]])==F) stop("Qualitative variables cannot appear in 'model.code': ",nodenam[i],sep="",call.=F)
    }
  if(!is.null(exogenous)) {
    for(i in 1:length(exogenous)) {
      if(isQuant(data[,exogenous[i]])==F) {
        if(sum(is.na(data[,exogenous[i]]))>0) stop("Qualitative variables cannot contain missing values: ",exogenous[i],sep="",call.=F)
        }
      }                                                                                                                                                                          
    }
  G <- new("graphNEL",nodes=names(pset),edgemode="directed")    
  for(i in 1:length(pset)) {
    if(length(pset[[i]])>0) {
      for(j in 1:length(pset[[i]])) {
        G <- addEdge(pset[[i]][j],names(pset)[i],G,1)
        }
      }
    }          
  topG <- topOrder(G)
  if(is.null(topG)) stop("The DAG contains directed cycles",call.=F)  
  global.control <- gconAdj(global.control)
  local.control <- lconAdj(local.control,global.control,pset)
  nodenam <- c(exogenous,topG)
  #
  #ndiff <- attr(data,"ndiff")
  #if(is.null(ndiff)) {
  #  ndiff <- 0
  #  #if(quiet==F) {
  #  #  pre <- winDialog("yesno",message="Data have not been preprocessed. Do you want to apply preprocessing with default settings?")
  #  #  if(identical(pre,"YES")) data <- preProcess(nodenam,group=group,time=time,log=T,data=data)
  #  #  }
  #  warning("It seems that data have not been preprocessed",call.=F)
  #  }
  #
  # preprocessing
  data_orig <- data
  data <- preProcess(x=nodenam,group=group,time=time,seas=NULL,data=data,
    log=log,diff.options=diff.options,imput.options=imput.options,quiet=quiet)
  #
  nomi <- c()
  optList <- callList <- lparList <- codeList <- vector("list",length=length(model.code))
  if(quiet==F) cat("Starting estimation ...")
  flush.console()
  for(i in 1:length(model.code)) {
    nomi[i] <- as.character(model.code[[i]])[2]
    if(quiet==F) {
      if(i>1) {
        inbl0 <- nchar(imess)
        } else {
        inbl0 <- 0  
        }
      imess <- paste("Estimating regression ",i,"/",length(model.code)," (",nomi[i],")",sep="")
      inbl <- max(0,inbl0-nchar(imess)+2)
      } else {
      imess <- NULL
      inbl <- 0
      }
    if(is.null(exogenous)) {
      iform <- model.code[[i]]
      } else {
      iform <- formula(paste(as.character(model.code[[i]])[2],"~",paste(exogenous,collapse="+"),"+",as.character(model.code[[i]])[3],sep=""))
      }
    iad <- local.control[["adapt"]][nomi[i]]
    ipar <- pset[[nomi[i]]]
    iscan <- scanForm(model.code[[i]])
    ilimit <- Inf
    if(length(ipar)>0) {
      #
      ilagtype <- iscan$ltype
      ilagnam <- names(ilagtype[which(ilagtype!="none")])
      if(length(ilagnam)>0) {
        ilimit <- findLagLim(data[,c(nomi[i],ilagnam,group)],group=group)
        if(sum(sapply(iscan$lpar,function(z){sum(is.na(z))}))>0) iad <- T
        }
      #
      iming <- rep(min(ilimit,global.control$min.gestation),length(ipar))
      iges <- rep(min(ilimit,global.control$max.gestation),length(ipar))
      iwd <- rep(min(ilimit,global.control$min.width),length(ipar))
      ilead <- rep(min(ilimit,global.control$max.lead),length(ipar))
      isg <- rep(global.control$sign,length(ipar))
      names(iming) <- names(iges) <- names(iwd) <- names(ilead) <- names(isg) <- ipar
      #
      } else {
      iming <- iges <- iwd <- ilead <- isg <- c()
      }
    #
    iauxming <- local.control$min.gestation[[nomi[i]]]
    iming[names(iauxming)] <- sapply(iauxming,function(z){min(ilimit,z)})
    iauxges <- local.control$max.gestation[[nomi[i]]]
    iges[names(iauxges)] <- sapply(iauxges,function(z){min(ilimit,z)})
    iauxwd <- local.control$min.width[[nomi[i]]]
    iwd[names(iauxwd)] <- sapply(iauxwd,function(z){min(ilimit,z)})
    iauxlead <- local.control$max.lead[[nomi[i]]]
    ilead[names(iauxlead)] <- sapply(iauxlead,function(z){min(ilimit,z)})
    iauxsg <- local.control$sign[[nomi[i]]]
    isg[names(iauxsg)] <- iauxsg
    #
    if(iad==T) optList[[i]] <- list(adapt=iad,min.gestation=iming,max.gestation=iges,min.width=iwd,max.lead=ilead,sign=isg)
    imod <- dlaglm(iform,group=group,data=data,adapt=iad,no.select=exogenous,min.gestation=iming,max.gestation=iges,min.width=iwd,max.lead=ilead,sign=isg,ndiff=diff.options$ndiff,gamma.by=gamma.by,mess=imess,nblank=inbl)
    callList[[i]] <- icall <- imod$call$formula
    iscan <- scanForm(icall)
    ixnam <- setdiff(names(iscan$ltype),c(group,time,exogenous))
    codeList[[i]] <- creatForm(nomi[i],ixnam,NULL,iscan$ltype[ixnam],iscan$lpar[ixnam],NULL)
    if(length(ixnam)>0) {
      ilpar <- data.frame(iscan$ltype[ixnam],do.call(rbind,iscan$lpar[ixnam]))
      colnames(ilpar) <- c("type","par1","par2")
      lparList[[i]] <- ilpar
      }
    if(hac==T) {
      imod$vcov <- doHAC(imod,group=group)
      class(imod) <- c("hac","lm")
      } else {
      imod$vcov <- vcov(imod)
      }
    estL[[i]] <- imod
    }
  names(optList) <- names(lparList) <- nomi
  if(quiet==F) {
    blchar <- ""
    if(iad==F) {
      if(nchar(imess)>=20) blchar <- paste(rep(" ",nchar(imess)-20),collapse="")
      } else {
      if(nchar(imess)>=11) blchar <- paste(rep(" ",nchar(imess)-11),collapse="")
      }
    cat('\r')
    cat("Estimation completed",blchar,"\n")
    }
  names(estL) <- nomi
  # fitted and residuals
  epsL <- fitL <- list()
  for(i in 1:length(estL)) {
    epsL[[i]] <- residuals(estL[[i]])
    fitL[[i]] <- fitted(estL[[i]])
    }
  epsOK <- data.frame(do.call(cbind,lapply(epsL,formatFit,n=nrow(data))))
  fitOK <- data.frame(do.call(cbind,lapply(fitL,formatFit,n=nrow(data))))
  if(!is.null(time)) {
    epsOK <- cbind(data[,time],epsOK)
    fitOK <- cbind(data[,time],fitOK)
    } 
  if(!is.null(group)) {
    epsOK <- cbind(data[,group],epsOK)
    fitOK <- cbind(data[,group],fitOK)
    }
  colnames(epsOK) <- colnames(fitOK) <- c(group,time,names(estL))
  # autocorrelation order
  if(hac==T) {
    acOrder <- lapply(estL,function(z){attr(z$vcov,"max.lag")})
    } else {
    if(is.null(group)) {
      acOrder <- c()
      for(i in 1:length(estL)) {
        acOrder[i] <- ar(na.omit(epsOK[,names(estL)[i]]))$order
        }
      names(acOrder) <- names(estL)
      } else {
      acOrder <- list()
      for(i in 1:length(estL)) {
        iac <- c()
        for(j in 1:length(gruppi)) {
          iind <- which(epsOK[,group]==gruppi[j])
          iac[j] <- ar(na.omit(epsOK[iind,names(estL)[i]]))$order
          }
        names(iac) <- paste(group,gruppi,sep="")
        acOrder[[i]] <- iac
        }
      names(acOrder) <- names(estL)
      }
    }
  # output
  out <- list(estimate=estL,model.code=codeList,call=callList,lag.par=lparList,exogenous=exogenous,group=group,time=time,
    log=attr(data,"log"),ndiff=attr(data,"ndiff"),
    data=data[,c(group,time,nodenam)],data.orig=data_orig,
    fitted=fitOK,residuals=epsOK,autocorr=acOrder,hac=hac,adaptation=optList)
  class(out) <- "dlsem"
  out
  }

# automated plots of lag shapes
auto.lagPlot <- function(x,cumul=FALSE,conf=0.95,plotDir=NULL) {
  if(("dlsem" %in% class(x))==F) stop("Argument 'x' must be an object of class 'dlsem'",call.=F)
  if(length(cumul)!=1 || !is.logical(cumul)) stop("Argument 'cumul' must be a logical value",call.=F)
  if(length(conf)!=1 || !is.numeric(conf) || conf<=0 || conf>=1)