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

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

ll_cv<-function(x,y,bandwidth){

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


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

nw_cv<-function(x,y,bandwidth){
apply(as.matrix(bandwidth),1,function(h){(1/length(x))*sum(((y-nw(x,x,y,h))/(1-(K(0)/apply(K(outer(x,x,"-")/h),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)

##--Save the values of the ross-validation function aplicated to no-parametric estimators:
#-- local-linear and Nadaraya-Watson.


##--Save the values of the ross-validation function aplicated to no-parametric estimators:
#-- local-linear and Nadaraya-Watson.
CV_ll<-numeric()
CV_nw<-numeric()


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


##-- Graphic representacion of of the cross-validation function aplicated to no-parametric estimators:
#-- local-linear and Nadaraya-Watson.
windows()
plot(CV_ll~h,type="l",main="cross validaton Linear-local")
windows()
plot(CV_nw~h,type="l",main="cross validation Nadaraya-Watson")


##-- 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))]
##-- Estimator of Nadaraya-Watson regresion with the cross-validation bandwidth.
nadaraya_watson<-nw(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],h[which(CV_nw==min(CV_nw,na.rm=T))])
h[which(CV_nw==min(CV_nw,na.rm=T))]
