######################################################
#            REGRESIÓN Y ANÁLISIS DE DATOS           #
######################################################

##--Data.
data(airquality)
names(airquality)

airquality$Ozone
airquality$Temp

##-- Positions where Ozone value Na.
ind<-which(is.na(airquality$Ozone))



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

##-- Ozone and Temperature time series.

windows()
par(mfrow=c(1,2))
plot(ts(airquality$Ozone,frequency=30), type="o", ylab="Ozone")
plot(ts(airquality$Temp,frequency=30), type="o", ylab="Temperature")


##-- Regression model.
windows()
lm(airquality$Ozone[-ind]~airquality$Temp[-ind])
plot(lm(airquality$Ozone[-ind]~airquality$Temp[-ind]))
windows()
plot(airquality$Temp[-ind],airquality$Ozone[-ind],main="Regression model")
abline(lm(airquality$Ozone[-ind]~airquality$Temp[-ind]),lwd=2,col="blue")


##-- Ozone and Temperature empirical distribution functions 
#- without positions before were difined (ind).

windows()
par(mfrow=c(1,2))
plot(ecdf(airquality$Ozone[-ind]),main="Ozone")
plot(ecdf(airquality$Temp[-ind]),main="Temperature")

 
##-- Ozone and Temperature density functions 
#- without positions before were difined (ind).
 
par(mfrow=c(1,2))
plot(density(airquality$Ozone[-ind]),main="Ozono")
plot(density(airquality$Temp[-ind]),main="Temperatura")
 

##-- Shapiro-Wilks test is applicated to Ozone and Temperature.
shapiro.test(airquality$Ozone[-ind])
shapiro.test(airquality$Temp[-ind])


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


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

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

K<-function(mas){
#--Gausiann kernel

return(1/sqrt(2*pi)*exp(-mas^2/2))
}


S<-function(t,mas,bandwidth,exponent){
#-- t vector of points where S will be evaluated.
#-- mas values of the independ variable.
#-- bandwidth window 
#-- r exponent.

    return(diag((t(K(outer(mas,t,"-")/bandwidth)))%*%(outer(mas,t,"-")^exponent)))
}


EstLL<-function(t,mas,y,bandwidth){
#-- This function returns the estimation which was caluted by no-parametric estimator
#-- of local-linear regresion.
#-- t vector of points where S will be evaluated.
#-- x values of the independ variable.
#-- y values of the depend variable.
#-- bandwidth is the window
        
    return( ( t( K( outer(mas,t,"-") /bandwidth) ) * ( S(t,mas,bandwidth,2) - t(outer(mas,t,"-")) * S(t,mas,bandwidth,1) ) /apply( t( K( outer(mas,t,"-") /bandwidth) ) * ( S(t,mas,bandwidth,2) - t(outer(mas,t,"-")) * S(t,mas,bandwidth,1) ),1,sum ) )%*%y)
}


###########################

##-- CROSS-VALIDATION BANDWIDTH

ll_cv<-function(mas,y,bandwidth){
#--It's the cross-validation function aplicated to no-parametric estimator
#-- of local-linear regresion.

    return( 1/length(mas) * sum( ( as.numeric(y - EstLL(mas,mas,y,bandwidth)) / (1 - ( K(0)*S(mas,mas,bandwidth,2) / apply( t( K( outer(mas,mas,"-") /bandwidth) ) * ( S(mas,mas,bandwidth,2) - t(outer(mas,mas,"-")) * S(mas,mas,bandwidth,1) ), 1, sum) )  ) )^2 ) )
}


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

##--Save the values of the ross-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])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],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<-EstLL(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],h[which(CV_ll==min(CV_ll,na.rm=T))])

########################################################
#  3.- Other window selection                          #
########################################################

library(KernSmooth)
##--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<-EstLL(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],dpill(airquality$Temp[-ind],airquality$Ozone[-ind]))

#########################################################
# 4.- Nadaraya-Watson estimator                         #
#########################################################

nw<-function(t,mas,y,bandwidth){
#-- This function returns the estimation which was caluted by no-parametric estimator
#-- of Nayadara-Watson regresion.
#-- t vector of points where S will be evaluated.
#-- x values of the independ variable.
#-- y values of the depend variable.
#-- bandwidth is the window

return((K(outer(t,mas,"-")/bandwidth)/apply(K(outer(t,mas,"-")/bandwidth),1,sum))%*%y)

}



cv_nw<-function(mas,y,bandwidth){
#--It's the cross-validation function aplicated to no-parametric estimator
#-- of Nadaraya-Watson regresion.

return(1/length(mas) * sum( ((y - nw(mas,mas,y,bandwidth)) / (1 - (K(0)/ apply( K(outer(mas,mas,"-")/bandwidth) ,1,sum) ) ) )^2 ))

}


##-- Sequency of windows
h<-seq(0,20,by=0.2)
##--Save the values of the ross-validation function aplicated to Nadaraya-Watson no-parametric estimator
#-- of regresion.

##--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])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],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<-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))])

##################################################################
#           Graphic Representation                               #
##################################################################

windows()
par(mfrow=c(2,2))
plot(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],xlab="Temperatura", ylab="Ozone",main="Local linear CV")
lines(local_linear~airquality$Temp[-ind][order(airquality$Temp[-ind])],lwd=2)


plot(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],xlab="Temperatura", ylab="Ozone",main="Local linear plug-in")
lines(local_linear_plugin~airquality$Temp[-ind][order(airquality$Temp[-ind])],col="red",lwd=2)


plot(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],xlab="Temperatura", ylab="Ozone",main="Nadaraya-Watson estimator")
lines(nadaraya_watson~airquality$Temp[-ind][order(airquality$Temp[-ind])],col="blue",lwd=2)

plot(airquality$Temp[-ind][order(airquality$Temp[-ind])],airquality$Ozone[-ind][order(airquality$Temp[-ind])],xlab="Temperatura", ylab="Ozone",main="No parametric Estimators")
lines(local_linear~airquality$Temp[-ind][order(airquality$Temp[-ind])],lwd=2)
lines(local_linear_plugin~airquality$Temp[-ind][order(airquality$Temp[-ind])],col="red",lwd=2,lty=2)
lines(nadaraya_watson~airquality$Temp[-ind][order(airquality$Temp[-ind])],col="blue",lwd=2,lty=3)




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