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

##-- Attaching the data.

data(airquality)
names(airquality)


##-- Positions where Ozone´s value is NA.
ind_1<-which(is.na(airquality$Ozone))

##-- Positions where Temperatures's values, whithout positions where Ozone´s value is NA,are ordered.
ind_2<-order(airquality$Temp[-ind_1])

##--Ozone and Temperature values whitout positions where Ozone´s value is NA.
Ozn<-airquality$Ozone[-ind_1]
Tmp<-airquality$Temp[-ind_1]


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

##-- Ozone and Temperature descriptions
summary(Ozn)
sd(Ozn)
summary(Tmp)
sd(Tmp)

##-- Ozone and Temperature box-plots
windows()
boxplot(Ozn,Tmp,main="Gráficos boxplot",names=c("Ozono","Temperatura"))
##--Positions where Ozone´s value is outlier. 
#--Temperature´s value don´t have outliers
outliers=c(which(Ozn[ind_2]==boxplot(Ozn)$out[1]),which(Ozn[ind_2]==boxplot(Ozn)$out[2]))


##-- 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 density functions 


kern_est<-function(x,h,mas){
# Description:
#   This function implements the kernel estimator expression:
#       f.hat(x)=1/(n*h) sum_{i=1}^{n}{K((x-X_i))/h},
#   for a Normal kernel.

# Arguments:
#   x<- The point(s) where to evaluate the estimator.
#   h<- The bandwidth's sequency to use.
#   mas<- The sample X_1,...,X_n to use.    

    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){
# Description:
#   This function implements the the least cross validation expression:

#       LCV(h)= 1/(n^{2} * h)*sum_{i,j}{K*K ((X_i-X_j)/h)} 
#             - 4/(n*(n-1)*h) *sum_{i=1}^{n} * sum_{j>i}{K((X_i-X_j)/h)},

#   for a Normal kernel (so the convolution K*K corresponds to a Normal density whith mu=0 and sd=sqrt(2)).
#   When minimazing this function over h, the cross validation bandwidth is obtained.

# Arguments:
#   h<- The bandwidth to use.
#   mas <- The sample X_1,...,X_n to use.
#   tri_matrix<- An upper triangular matrix (logical) of dimension (n x n).

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

estim<-function(mas,subint,grid,tri_matrix){
# Description:
#   This function computes cross-validation and bandwiths for a sample. 
#   In addition, the kernel estimator is performed whith that bandwidth.

# Arguments:
#   mas<- The sample.
#   grid<- A grid where to evaluate the estimator.
#   sub_int<- Subintervals where the optimization of LCV is carried out.
#   tri_matrix<- An upper triangular matrix (logical) of dimension (n x n).


# Output:
#   This function returns the bandwidth.
#   It also returns the approximated density (for cross validation ) over a grid.

    #Compute the cross - validation bandwidth by minimazind the LCV function (in subintervals).
    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)  ## Kernel estimation over grid with h_cv

    return(list(as.numeric(est),h))

}

s<-seq(0.5,8,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)))

#--Temperature
gridTmp<-seq(min(Tmp),max(Tmp),by=0.5)
resultsTmp<-estim(Tmp,subint_8,gridTmp,tri_matrix)
estimatorTmp<-resultsTmp[[1]]
bandwidthTmp<-resultsTmp[[2]]

#--Ozone
gridOzn<-seq(min(Ozn),max(Ozn),by=0.5)
resultsOzn<-estim(Ozn,subint_8,gridOzn,tri_matrix)
estimatorOzn<-resultsOzn[[1]]
bandwithOzn<-resultsOzn[[2]]

##-- Graffic Representation
windows()
par(mfrow=c(1,2))
hist(Tmp,freq=F,main="Densidad de la Temperatura",xlab="Temperatura")
lines(estimatorTmp~gridTmp,main="Ventana CV; h=4.728413",col="red")
hist(Ozn,freq=F,main="Densidad del Ozono",xlab="Ozono",ylim=c(0,0.02),xlab="Ozono")
lines(estimatorOzn~gridOzn,main="Ventana CV; h= 6.428435",col="blue")


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


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


##-- Scatterplot of ozone~temperature
windows()
plot(Ozn~Tmp,xlab="Temperatura",ylab="Ozono",main="Gráfico de dispersión")


##-- Regression model.
windows()
Modelo<-lm(Ozn~Tmp)
plot(Modelo)
summary(Modelo)
windows()
plot(Tmp,Ozn,main="Modelo de regresión simple",xlab="temperatura",ylab="ozono")
abline(Modelo,lwd=2,col="blue")


##-- Comparations between boxplot 
Comp<-Tmp
Comp[which(Comp<=75)]<-1
Comp[which(Comp>75)]<-2
windows()
boxplot(Ozn~Comp,names=c("temp<=75","temp>75"),main="Gráficos boxplot del ozono")
 


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


evals<-function(x,t,h){
# Description:
#   This function implements two auxiliar functions:
#       f(x,t)=(x-t)
#       g(x,t)=f(x,t)/h

#The following expresion will be calculated with this function:
#             (Xj-ti)/h
#             K[(Xj-ti)/h]
#   where K is kernel estimator.


# Arguments:
#   x<- The sample.
#   t<- The point(s) where to evaluate the estimator.
#   h<- The bandwidth to use.
 

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


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

# Description
#    This function implenent an auxilar function. This function will be used to caltulate the Local Linear estimator 
#    and its function of cross- validation.


# Arguments
#    x<-The sample
#    r<- The exponent
#    mat and mat2 two matrix

    return(apply(matrix(mat*mat2^r,nr=length(x)),2,sum))

}


l<-function(x,t,h,m,opt_m,opt_cv,ev){

# Description
#      This function calcule the matrix L of Local-Linear or Nadaraya-Watson estimator and its diagonal.
#      Lii will be used in cross-validation function of Local-linear estimation.


# Arguments
#    x<-      The sample
#    t<-      The point(s) where to evaluate the estimator.
#    m<-      A matrix
#    opt_m<-  This argument allows to choose an estimator: local linear.
#                                                         Nadaraya-Watson. 
#
#    opt_cv<- This argument allows to choose between to do cross-validation or not.
#    ev<-     In this argument are saved the results given by evals function.

# Output
#     This function return only the matrix L of Nadaraya- Watson or Local-Linear estimator or this matrix and its diagonal.

    
    mat_1<-matrix(t(dnorm(ev[[2]])),nc=length(t),byrow=T) ## Evalue K[(Xj-ti)/h]
    b<- t(matrix(mat_1,nr=length(x)))

    if(opt_m==1){
       #-- It was choosen Local-Linear estimator
        mat_2<- matrix(t(ev[[1]]),nc=length(t),nr=length(x)*length(h),byrow=T)
        b<-b*(S(x,2,mat_1,mat_2) - t(matrix(mat_2,nr=length(x)))* S(x,1,mat_1,mat_2))
        l_ll<-b/apply(b,1,sum)  #-- L

        if(opt_cv==1){
            #-- It was choosen to do cross-validation.
            b_diag<-matrix(t(l_ll),nr=length(x)*length(h))
            l_diag<-matrix(b_diag[m==T],nc=length(h),byrow=T) #--Lii

            return(list(l_ll,l_diag))  #-- L and Lii
        }

        else{
            #-- It wasn´t choosen to do cross-validation.
            return(l_ll) #-- L
        }

    }
    else{
        #-- It was choosen Local-Linear estimator
        l_nw<-b/apply(b,1,sum) #-- L

        if(opt_cv==1){
            #-- It was choosen to do cross-validation.
            b_diag<-matrix(t(l_nw),nr=length(x)*length(h))
            l_diag<-matrix(b_diag[m==T],nc=length(h),byrow=T) ##--Lii

            return(list(l_nw,l_diag)) #-- L and Lii
        }
        else{
        #-- It wasn´t choosen to do cross-validation.
        return(l_nw) #-- L
        }
    
    }
        
}

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

# Description
#    This function evaluate the Local Linear and Nadaraya-Watson estimator in a sequency, using cross-validation or plug-in bandwith


# Arguments
#    t<-   The point(s) where to evaluate the estimator.
#    x<-   The sample's  dependent variable
#    y<-   The sample's  independent variable
#    h<-   The bandwith's sequency
#    mat<- The matrix


# Output
#   This function returns the bandwidths (h) which are obtained to cross-validation or plug-in.
#   It also returns  graphics associated to cross-valitation to Local-lineal and Nadaraya-Watson and a graphic asociated to estimation local linear.
#   Moreover, it returns a comparative graphic between local linear and Nadaraya-Watson estimation and 
#   finally, a comparative graphic between local linear estimation whith cross-validation bandwith and plug-in bandwith.





    ## Local - linear

    ev<-evals(x,t,h)

    ll_cv<-l(t,x,h,mat,1,1,ev)

    m_ll<-matrix(ll_cv[[1]]%*%y,nc=length(h),byrow=T)

    cv_ll<-1/length(x) * apply(( (matrix(y,nr=length(x),nc=length(h)) - m_ll ) / (1 - ll_cv[[2]]) )^2 ,2,sum)
    h_cv_ll<-h[which.min(cv_ll)]
    print(h_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(x,x,h,mat,2,1,ev)

    m_nw<-matrix(nw[[1]]%*%y,nc=length(h),byrow=T)

    cv_nw<-1/length(x) * apply( ( ( matrix(y,nr=length(x),nc=length(h)) - m_nw ) / (1 - nw[[2]]) )^2 ,2,sum)
    h_cv_nw<-h[which.min(cv_nw)]
    print(h_cv_nw)

    windows()
    plot(cv_nw~h,type="l",main="CV para E.N.W.",lwd=2)

    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

    h_pi<- dpill(x,y)
    print(h_pi)
    
    ev<-evals(x,t,h_pi)
    ll_dpill<-l(x,x,h_pi,mat,1,0,ev)

    m_pi<-ll_dpill%*%y
    
    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"))


}


library(KernSmooth)

## Results

h<-seq(0.5,5,by=0.01)
mat<-matrix(diag(T,nr=length(Tmp),nc=length(Tmp)),nc=length(Tmp),nr=length(Tmp)*length(h),byrow=T)
results<-main(Tmp,Tmp,Ozn,h,mat)
    

##-- Without outliers: 
# Observations 70,83. Sample in order

## Results without outliers
Tmp<-Tmp[-outliers]
ind_2<-order(Tmp)
mat<-matrix(diag(T,nr=length(Tmp),nc=length(Tmp)),nc=length(Tmp),nr=length(Tmp)*length(h),byrow=T)
Ozn<-Ozn[-outliers]
results<-main(Tmp,Tmp,Ozn,h,mat)
    
    





##########################################################
#   5.- Sample model: Piecewise Linear Regression        #
##########################################################
##-- Positions where Temperatures's values, whithout positions where Ozone´s value is NA,are ordered.
ind_2<-order(airquality$Temp[-ind_1])

##--Ozone and Temperature values whitout positions where Ozone´s value is NA.
Ozn<-airquality$Ozone[-ind_1]
Tmp<-airquality$Temp[-ind_1]


windows()
knots<-c(75,91,96)
piecese1<-1:max(which(Tmp[ind_2]==75))
piecese2<-max(which(Tmp[ind_2]==75)):max(which(Tmp[ind_2]==91))
piecese3<-max(which(Tmp[ind_2]==91)):max(which(Tmp[ind_2]==96))
plot(Tmp[ind_2],Ozn[ind_2],xlab="Temperatura",ylab="Ozono", main="Modelo alternativo")

lines(Tmp[ind_2][piecese1],lm(Ozn[ind_2][piecese1]~Tmp[ind_2][piecese1])$fitted.values,lwd=2,col="blue")

lines(Tmp[ind_2][piecese2],lm(Ozn[ind_2][piecese2]~Tmp[ind_2][piecese2])$fitted.values,lwd=2,col="yellow")

lines(Tmp[ind_2][piecese3],lm(Ozn[ind_2][piecese3]~Tmp[ind_2][piecese3])$fitted.values,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. 

windows()
par(mfrow=c(2,2))
model_spline<-lm(Ozn~bs(Tmp,knots=knots,degree=1),data=airquality)
summary(model_spline)
plot(Tmp,Ozn, xlab="Temperatura", ylab="Ozono",main="Regresión spline (grado 1)")
new <- data.frame(Tmp = seq(min(Tmp),max(Tmp),1))
lines(new$Tmp[-length(new$Tmp)],predict(model_spline,new)[-length(predict(model_spline,new))],lty=1,col="red",lwd=3)


model_spline<-lm(Ozn~bs(Tmp,knots=knots,degree=2),data=airquality)
summary(model_spline)
plot(Tmp,Ozn, xlab="Temperatura", ylab="Ozono",main="Regresión spline (grado 2)")
new <- data.frame(Tmp = seq(min(Tmp),max(Tmp),1))
lines(new$Tmp[-length(new$Tmp)],predict(model_spline,new)[-length(predict(model_spline,new))],lty=1,col="red",lwd=3)



model_spline<-lm(Ozn~bs(Tmp,knots=knots,degree=3),data=airquality)
summary(model_spline)
plot(Tmp,Ozn, xlab="Temperatura", ylab="Ozono",main="Regresión spline (grado 3)")
new <- data.frame(Tmp = seq(min(Tmp),max(Tmp),1))
lines(new$Tmp[-length(new$Tmp)],predict(model_spline,new)[-length(predict(model_spline,new))],lty=1,col="red",lwd=3)
