
k<-function(x){     #kernel
return((1/sqrt(2*pi))*exp(-x^2/2))
}

nw<-function(t,x,y,h){#Nadaraya-Watson
    txh<-outer(t,x,"-")/h     #Da todas las restas t-x
    ktxh<-k(txh)              #Hace k de la matriz
    L<-ktxh/apply(ktxh,1,sum) #da una matriz donde las filas suman uno
                        #divido cada fila por cada fila 
    return(L%*%y)}
    


matrizLNW<-function(x,h){
    #usamos nw para calcukar matriz L
    txh<-outer(x,x,"-")/h     #Da todas las restas t-x
    ktxh<-k(txh)              #Hace k de la matriz
    L<-ktxh/apply(ktxh,1,sum) 
    return(L)
}


CV<-function(x,y,h){
    L=matrizLNW(x,h)
    mgorro=L%*%y
    mean((y-mgorro)^2/(1-diag(L))^2)
}

library(MASS)
plot(mcycle)
x=mcycle$times
y=mcycle$accel
y0<-y[order(x)]
y1<-y0[-1]
y2<-y0[-length(y)]
sigma<-(1/2*(length(y)-1))*sum((y1-y2)^2)
bandwith<-seq(0,8,by=0.01)
CV_NW<-apply(as.matrix(bandwith),1,function(h){CV(x,y,h)})
bandwith[which(CV_NW==min(CV_NW,na.rm=T))]
mgorro=nw(x,x,y,bandwith[which(CV_NW==min(CV_NW,na.rm=T))])#h=3
points(x,mgorro,type="l",col="red",lwd=5)

banda_inf=mgorro-qnorm(1-(alpha/2))*sd(mgorro)
banda_sup=mgorro+qnorm(1-(alpha/2))*sd(mgorro)
lines(x,banda_inf,col="red")
lines(x,banda_sup,col="red")


banda_inf1=mgorro-qnorm(1-(alpha/2))
*sqrt(sum(matrizLNW(x,bandwith[which(CV_NW==min(CV_NW,na.rm=T))])^2)*sigma)
banda_sup1=mgorro+qnorm(1-(alpha/2))
*sqrt(sum(matrizLNW(x,bandwith[which(CV_NW==min(CV_NW,na.rm=T))])^2)*sigma)
lines(x,banda_inf1,col="blue")
lines(x,banda_sup1,col="blue")
