K<-function(x){
1/sqrt(2*pi)*exp(-x^2/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.



CV_nw<-numeric()



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


windows()
plot(CV_nw~h,type="l",main="cross validation Nadaraya-Watson")



##-- 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))]
