# R/dlsem.r In dlsem: Distributed-Lag Linear Structural Equation Models

#### Documented in as.graphNELautoCodeauto.lagPlotcausalEffcompareModelsdlsemdrawSampleecqgamisIndeplagPlotlagShapesldlmHACqdresidualPlotunirootTest

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

#########################################
#
###  Available constrained lag shapes ###
#
#  "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  #####
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)
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)) {
}
res <- cbind(auxmat,limmat)
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)
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)
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
}
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]]
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,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]
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) {
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)) {
} 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) {
} 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(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))})
}
} 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)
#
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
while(is.na(res[1])||(abs(res[1])>1.6 & k>0)) {
k <- k-1
}
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()
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) {
} 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,...) {
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)
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) {
} else {
}
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) {
for(i in 2:length(est)) {
}
}
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)
nomi <- names(x)
warn1 <- warn2 <- 0
if(!is.null(x) && !is.list(x)) warn2 <- warn2+1
if(length(unknam)>0) {
x <- x[setdiff(nomi,unknam)]
warning("Some components with unknown names in argument 'global.control' were ignored",call.=F)
}
#
} else {
warn2 <- warn2+1
}
} else {
}
#
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(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)
}

# adjust local control options (internal use only)
warn1 <- warn2 <- 0
nomi <- names(pset)
if(!is.null(x) && !is.list(x)) warn2 <- warn2+1
for(j in 1:length(nomi)) {
} else {
warn2 <- warn2+1
}
} else {
}
}
#
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)
}

# 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)
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)
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(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) {
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)
#
if(is.null(diff.options\$test)) {
if(is.null(group)) {
n <- nrow(data)
} else {
n <- min(table(data[,group]))
}
}
ndiff <- diff.options\$ndiff
#
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) {
} 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) {
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
}
}
}
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(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)) {
} 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)) {
} 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)))
pset[[length(pset)+1]] <- character(0)
}
}
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) {
}
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]])) {
}
}
}
topG <- topOrder(G)
if(is.null(topG)) stop("The DAG contains directed cycles",call.=F)
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=""))
}
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)
}
#
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))
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)})
iauxsg <- local.control\$sign[[nomi[i]]]
isg[names(iauxsg)] <- iauxsg
#
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(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,