K<-function(x){
1/sqrt(2*pi)*exp(-x^2/2)
}
S<-function(t,x,bandwidth,r){
diag(t(K(outer(x,t,"-")/bandwidth))%*%(outer(x,t,"-")^r))
}

EstLL<-function(t,x,y,bandwidth){        
return(apply(as.matrix(bandwidth),1,
function(h){(t(K(outer(x,t,"-")/h))*(S(t,x,h,2)-t(outer(x,t,"-"))*S(t,x,h,1))
/apply( t(K(outer(x,t,"-")/h))*(S(t,x,h,2) - t(outer(x,t,"-")) * S(t,x,h,1)),1,sum ) )%*%y}))
}

ll_cv<-function(x,y,bandwidth){
return(apply(as.matrix(bandwidth),1,
function(h){1/length(x)*sum(((as.numeric(y -EstLL(x,x,y,h)))
/(1-( K(0)*S(x,x,h,2)/
    apply(t(K(outer(x,x,"-")/h))*(S(x,x,h,2)-t(outer(x,x,"-"))*S(x,x,h,1)),1,sum))))^2)}))
}

data(airquality)
names(airquality)

##-- Positions where Ozone value Na.
ind<-which(is.na(airquality$Ozone))

##- Windows sequency
h<-seq(0,20,by=0.02)

CV_ll<-numeric()
CV_ll<-ll_cv(airquality$Temp[-ind][order(airquality$Temp[-ind])],
airquality$Ozone[-ind][order(airquality$Temp[-ind])],h)

windows()
plot(CV_ll~h,type="l",main="cross validaton Linear-local")

##-- Estimator of local-linear regresion with the cross-validation bandwidth.
local_linear<-EstLL(airquality$Temp[-ind][order(airquality$Temp[-ind])],
airquality$Temp[-ind][order(airquality$Temp[-ind])],
airquality$Ozone[-ind][order(airquality$Temp[-ind])],h[which(CV_ll==min(CV_ll,na.rm=T))])
h[which(CV_ll==min(CV_ll,na.rm=T))]
