
######################################################
######################################################
#            REGRESSION KERNEL ESTIMATION            #
######################################################
######################################################

##-- Attaching the data.

data(airquality)
names(airquality)



##-- Positions where Ozone´s value is NA.

ind_1<-which(is.na(airquality$Ozone))
ind_2<-order(airquality$Temp[-ind_1])

Ozn<-airquality$Ozone[-ind_1]
Tmp<-airquality$Temp[-ind_1]


######################################################
# 1.- Analisys of univariant distribution            #
######################################################

##-- Temperature time series.

windows()
temp.ts<-ts(airquality$Temp)
plot(temp.ts, type="o", ylab="Temperatura",xlab="Día",main="Serie de tiempo de la temperatura")


##-- Ozone and Temperature Empirical Distribution Functions 

windows()
par(mfrow=c(1,2))
plot(ecdf(Ozn),main="Ozono")
plot(ecdf(Tmp),main="Temperatura")

##-- Ozone and Temperature Empirical Distribution Functions 

windows()
par(mfrow=c(1,2))
plot(ecdf(Ozn),main="Ozone")
plot(ecdf(Tmp),main="Temperature")


##-- Ozone and Temperature density functions 


kern_est<-function(x,h,mas){
    return(1/(length(mas)*h)*t(matrix(apply(dnorm(apply(outer(x,mas,"-"),2,function(x)outer(x,h,"/"))),1,sum),ncol=length(h))))
}

LCV_sub<-function(h,mas,tri_matrix){
return(sum(dnorm(outer(mas,mas,"-")/h,sd=sqrt(2)))/(length(mas)^2*h)
    - (4/(length(mas)*(length(mas)-1)*h)) * sum(dnorm((outer(mas,mas,"-")/h)[tri_matrix])))
}

main<-function(mas,subint,grid,tri_matrix){

opt<- apply(subint,1,function(y) as.numeric(optimize(LCV_sub,interval=c(y[1],y[2]),mas=mas,tri_matrix=tri_matrix_sub)))
h<-opt[1,][which.min(opt[2,])]
est<- kern_est(grid,h,mas)

return(as.numeric(est))

}

s<-seq(0,8,by=0.5)
grid<-seq(min(Tmp),max(Tmp),by=0.5)
subint_8<-cbind(c(0,s),c(s,0))[-1,][-length(s),]
tri_matrix_sub<-upper.tri(matrix(0,nr=length(Tmp),nc=length(Tmp)))

mas<-Tmp
estimador<-main(mas,subint_8,grid,tri_matrix)

windows() 
par(mfrow=c(1,2))
hist(Ozn,freq=F,main="Densidad del ozono")
lines(density(Ozn),main="Ozono",col="blue")
hist(Tmp,freq=F,main="Densidad de la temperatura")
lines(density(Tmp),main="Temperatura",col="blue")


##-- Shapiro-Wilks tests 
shapiro.test(Ozn)
shapiro.test(Tmp)


##-- Regression model.

windows()
lm(Ozn~Tmp)
plot(lm(Ozn~Tmp))
windows()
plot(Tmp,Ozn,main="Regression model")
abline(lm(Ozn~Tmp),lwd=2,col="blue")
 

##-- Ozone and Temperature box-plots
windows()
par(mfrow=c(1,2))
boxplot(airquality$Ozone[-ind1],main="Ozone")
boxplot(airquality$Temp[-ind1],main="Temperature")


##-- Ozone and Temperature descriptions
summary(airquality$Ozone[-ind_1])
summary(airquality$Temp[-ind_1])


#######################################################
# 2.- Local-linear estimator                          #
#######################################################


evals<-function(x,t,h){

# Description
#

# Arguments
#

# Output
#ev_1=x-t
#ev_2=(x-t)/h

    ev_1<-outer(x,t,"-")
    ev_2<-outer(x,t,"-")/h
    return(list(ev_1,ev_2))
}



S<-function(r,mat,mat2){

# Description
#

# Arguments
#

# Output
#

#-- x values of the independ variable.
#-- bandwidth window 
#-- r exponent.
 return(t(mat)%*%((mat2)^r))
 #return(diag((t())%*%(outer(mas,t,"-")^exponent)))

}


l<-function(opt_m,opt_cv,ev){

# Description
#

# Arguments
#

# Output
#

   
        b<-t( dnorm(ev[[2]]) )
    if(opt_m==1){
        

        b<-b * ( S(2,ev[[2]],ev[[1]]) - t(ev[[1]]) * S(1,ev[[2]],ev[[1]]) )
        l_ll<-b/apply(b,1,sum) 

        if(opt_cv==1){

            
            l_diag<-diag(l_ll)

            return(list(l_ll,l_diag))
        }

        else{
            return(l_ll)
        }

    }
    else{
    
        l_nw<-b/apply(b,1,sum)

        if(opt_cv==1){

            l_diag<-diag( l_nw)

            return(list(l_nw,l_diag))
        }
        else{
        return(l_nw)
        }
    
    }
        
}


main<-function(t,x,y,h,op){

# Description
#

# Arguments
#

# Output
#
    if(op==1){
    ## Local - linear

    ev<-evals(x,t,h)

    ll_cv<-l(1,1,ev)

    m_ll<-ll_cv[[1]]%*%y
    



    cv_ll<-1/length(x) * (y  - m_ll )/ (1 - ll_cv[[2]]) )^2 
    h_cv_ll<-h[which.min(cv_ll)]
    
    #windows()
    #plot(cv_ll~h,type="l",lwd=2,main="CV para E.L.L.",ylab="CV(h)")

    windows()
    plot(Tmp[ind_2],Ozn[ind_2],xlab="Temperatura", ylab="Ozono",main="Estimación por el E.L.L")
    lines(m_ll[,which.min(cv_ll)][ind_2]~Tmp[ind_2],lwd=2)

    ## Nadaraya - Watson 

    nw<-l(2,1,ev)

    m_nw<-nw[[1]]%*%y
    cv_nw<-1/length(mas) * sum( ((y - m_nw) / (1 - nw[[2]]) ) )^2 
    h_cv_nw<-h[which.min(cv_nw)]

    #windows()
    #plot(cv_nw~h,type="l",main="Cross Validaton Nadaraya-Watson")

    windows()
    plot(Tmp[ind_2],Ozn[ind_2],xlab="Temperatura", ylab="Ozono",main="Comparación L.L - N.W")
    lines(m_ll[,which.min(cv_ll)][ind_2]~Tmp[ind_2],lwd=2)
    lines(m_nw[,which.min(cv_nw)][ind_2]~Tmp[ind_2],lwd=2,col="red")
    legend(60,150,legend=c("L.L.","N.W"),lwd=c(2,2),col=c("black","red"))
    }
    ## Plug - in
    if(op==2){
    library(KernSmooth)

    h_pi<- dpill(airquality$Temp[-ind_1],airquality$Ozone[-ind_1])
    ev<-evals(x,t,h_pi)
    ll_dpill<-l(1,0,ev)

    m_pi<-ll_dpill%*%y
    
    ev<-evals(x,t,h)

    ll_cv<-l(1,1,ev)

    m_ll<-ll_cv[[1]]%*%y
    cv_ll<-1/length(x) * sum( ( as.numeric(y  - m_ll ) )/ (1 - ll_cv[[2]]) )^2 
    h_cv_ll<-h[which.min(cv_ll)]
    
    
    windows()
    plot(Tmp[ind_2],Ozn[ind_2],xlab="Temperatura", ylab="Ozone",main="Comparación L.L. con ventanas C.V. y P.I.")
    lines(m_ll[,which.min(cv_ll)][ind_2]~Tmp[ind_2],lwd=2)
    lines(m_pi[ind_2]~Tmp[ind_2],lwd=2,col="blue")
    legend(60,150,legend=c("C.V.","P.I."),lwd=c(2,2),col=c("black","blue"))
}

}
## Results

h<-seq(0.5,5,by=0.02)
results1<-apply(as.matrix(h),1,function(x) main(Tmp,Tmp,Ozn,x,1))
    
results2<-apply(as.matrix(h),1,function(x) main(Tmp,Tmp,Ozn,x,2))



##-- Without outliers: 
# Observations 57,70,83. Sample in order
plot(lm(airquality$Temp[-ind][order(airquality$Temp[-ind])]~airquality$Ozone[-ind][order(airquality$Temp[-ind])]))

windows()
par(mfrow=c(1,2))


## Local-linear estimator using cross-validation ##

##-- Sequency of windows
h<-seq(0,20,by=0.02)

##--Save the values of the cross-validation function aplicated to no-parametric estimator
#-- of local-linear regresion.
CV_ll<-numeric()
for(i in 1:length(h)){CV_ll[i]<-ll_cv(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],h[i])}

##-- Graphic representacion of of the ross-validation function aplicated to no-parametric estimator
#-- of local-linear regresion respect to the sequuency of windows.
plot(CV_ll~h,type="l",main="cross validaton Linear-local")

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


## Local-linear estimator using plug-in ##

##--Local linear Gaussian kernel regression estimate.
##--Use direct plug-in methodology to select the bandwidth of a local linear Gaussian kernel regression estimate.
local_linear_plugin_without<-EstLL(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],dpill(airquality$Temp[-ind],airquality$Ozone[-ind]))


## Nadaraya-Watson estimation using cross-validation ##

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

##--Save the values of the ross-validation function aplicated to no-parametric estimator
#-- of Nadaraya-Watson regresion.
CV_nw<-numeric()
for(i in 1:length(h)){CV_nw[i]<-cv_nw(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],h[i])}

##-- Graphic representacion of of the ross-validation function aplicated to no-parametric estimator
#- of Nadaraya-Watson regresion respect to the sequency of windows.
plot(CV_nw~h,type="l",main="cross validation Nadaraya-Watson")

##-- Estimator of Nadaraya-Watson regresion with the cross-validation bandwidth.
nadaraya_watson_without<-nw(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],h[which(CV_nw==min(CV_nw,na.rm=T))])



windows()
par(mfrow=c(2,2))

plot(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],xlab="Temperatura without outliers", ylab="Ozone without outliers",main="Local linear CV")
lines(local_linear_without~airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],lwd=2)

plot(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],xlab="Temperature without outliers", ylab="Ozone without outliers",main="Local linear plug-in")
lines(local_linear_plugin_without~airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],col="red",lwd=2)

plot(airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],airquality$Ozone[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],xlab="Temperature without outliers", ylab="Ozone without outliers",main="Nadaraya-Watson estimator")
lines(nadaraya_watson_without~airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],col="blue",lwd=2)


##########################################################
#   5.- Sample model: Piecewise Linear Regression        #
##########################################################

windows()
par(mfrow=c(2,2))
plot(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],xlab="Temperature",ylab="Ozone", main="Local-lineal estimator")
legend(x=60, y=100, legend = c("Plug-in","Cross-validation"), col=c("blue", "green"), pch=1, bty="n")

lines(airquality$Temp[-ind][order(airquality$Temp[-ind])][1:max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==75))],lm(airquality$Ozone[-ind][order(airquality$Temp[-ind])][1:max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==75))]~airquality$Temp[-ind][order(airquality$Temp[-ind])][1:max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==75))])$fitted.values,lwd=2,col="blue")

lines(airquality$Temp[-ind][order(airquality$Temp[-ind])][max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==75)):max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==91))],lm(airquality$Ozone[-ind][order(airquality$Temp[-ind])][max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==75)):max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==91))]~airquality$Temp[-ind][order(airquality$Temp[-ind])][max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==75)):max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==91))])$fitted.values,lwd=2,col="blue")

lines(airquality$Temp[-ind][order(airquality$Temp[-ind])][max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==91)):max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==96))],lm(airquality$Ozone[-ind][order(airquality$Temp[-ind])][max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==91)):max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==96))]~airquality$Temp[-ind][order(airquality$Temp[-ind])][max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==91)):max(which(airquality$Temp[-ind][order(airquality$Temp[-ind])]==96))])$fitted.values,lwd=2,col="blue")

lines(local_linear_without~airquality$Temp[-ind][order(airquality$Temp[-ind])][-c(57,70,83)],lwd=2,col="green")


##-Whith splines


library(splines)
##-- This package provides functions for working with regression splines using the B-spline basis, bs, and the natural cubic spline basis, ns. 

##-- function bs: Generate the B-spline basis matrix for a polynomial spline
##-- ARGUMENTS
#- degree: degree of the piecewise polynomial?default is 3 for cubic splines.
#- knots:  the internal breakpoints that define the spline. The default is NULL, which results in a basis for ordinary polynomial regression. 
#-         Typical values are the mean or median for one knot, quantiles for more knots. 


lm(Ozone~bs(Temp,knots=c(75,91,96),degree=1),data=airquality)
summary(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=1),data=airquality))
plot(airquality$Temp,airquality$Ozone, xlab="Temperature", ylab="Ozone",
main="Piecewise Linear Regression 1")
new <- data.frame(Temp = seq(min(airquality$Temp),max(airquality$Temp),1))
lines(new$Temp[-length(new$Temp)],predict(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=1,df=3),data=airquality),new)[-length(predict(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=1,df=3),data=airquality),new))],lty=1,col="red",lwd=3)


lm(Ozone~bs(Temp,knots=c(75,91,96),degree=2),data=airquality)
summary(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=2),data=airquality))
plot(airquality$Temp,airquality$Ozone, xlab="Temperature", ylab="Ozone",
main="Piecewise Linear Regression 2")
new <- data.frame(Temp = seq(min(airquality$Temp),max(airquality$Temp),1))
lines(new$Temp[-length(new$Temp)],predict(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=2),data=airquality),new)[-length(predict(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=1,df=3),data=airquality),new))],lty=1,col="red",lwd=3)


lm(Ozone~bs(Temp,knots=c(75,91,96),degree=3),data=airquality)
summary(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=3),data=airquality))
plot(airquality$Temp,airquality$Ozone, xlab="Temperature", ylab="Ozone",
main="Piecewise Linear Regression 3")
new <- data.frame(Temp = seq(min(airquality$Temp),max(airquality$Temp),1))
lines(new$Temp[-length(new$Temp)],predict(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=3),data=airquality),new)[-length(predict(lm(Ozone~bs(Temp,knots=c(75,91,96),degree=1,df=3),data=airquality),new))],lty=1,col="red",lwd=3)
