K<-function(x){
(1/sqrt(2*pi))*exp(-x^2/2)
}


L<-function(t,x,bandwidth){
K(outer(t,x,"-")/bandwidth)/apply(K(outer(t,x,"-")/bandwidth),1,sum)
}

nw<-function(t,x,y,bandwidth){
L(t,x,bandwidth)%*%y
}


CV<-function(x,y,bandwidth){
mean((y-nw(x,x,y,bandwidth)/(1-diag(L(x,x,bandwidth))))^2)
}

##DATOS
library(MASS)
data(mcycle)
plot(mcycle)
names(mcycle)
x<-mcycle$times
y<-mcycle$accel

##SELECCIÓN DEL PARÁMETRO
bandwidth=seq(0,8,by=0.02)
CV_nw<-apply(as.matrix(bandwidth),1,function(h) CV(x,y,h))

##ESTIMACIÓN DE LA REGRESIÓN
h_opt<-bandwidth[which(CV_nw==min(CV_nw,na.rm=T))]
regresion<-nw(x,x,y,bandwidth[which(CV_nw==min(CV_nw,na.rm=T))])
points(x,regresion,lwd=2,type="l", col="blue")

##CONSTRUCCIÓN BANDAS DE CONFIANZA
#--Cálculo de sigma
y0<-sort(y)
y1<-y0[-1]
y2<-y0[-length(y)]
sigma<-(1/(2*(length(y)-1)))*sum((y1-y2)^2)

alpha=0.05
inf<-regresion-qnorm(1-(alpha/2))*sqrt(sum(L(x,x,h_opt)^2)*sigma)
sup<-regresion+qnorm(1-(alpha/2))*sqrt(sum(L(x,x,h_opt)^2)*sigma)
lines(x,inf,col="green")
lines(x,sup,col="green")
