#  Importation des données# 
demo=read.table("C:/Documents and Settings/CHEIK/Bureau/ETE 2014/plage.csv",sep=",",header=TRUE)
tail(demo)

#  Nettoyage des observations manquantes de l'année 2014 #
demo=demo[!is.na(demo[,2]),]
tail(demo)

# Transformation en données mensuelles
hebdo=ts(demo$plage,start=2004, end=c(2014,7),frequency=52) 
plot(hebdo) 

H2M=function(BASE){ 
X=BASE[,2] 
Y=BASE[,1] 
date1=substr(as.character(Y),1,10) 
date2=substr(as.character(Y),14,23) 
D1=as.Date(date1,"%Y-%m-%d") 
D2=as.Date(date2,"%Y-%m-%d") 
vm=vy=N=NA 
for(t in 1:length(D1)){ 
  mois=seq(D1[t],D2[t],length=7) 
  vm=c(vm,as.POSIXlt(mois)$mon+1) 
  vy=c(vy,as.POSIXlt(mois)$year+1900) 
  N=c(N,rep(X[t],7))} 
N=N[-1]; 
vm=vm[-1]; 
vy=vy[-1] 
YM=vy*100+vm 
Z=tapply(N,as.factor(YM),mean)
Zts=ts(as.numeric(Z),start=c(2004,1),frequency=12) 
return(Zts)} 

# Application aux données

plage=H2M(demo)
plage

# Extraction d'un échantillon pour la prévision
b=7
n=length(plage)
base=ts(plage[-((n-b+1):n)],start=c(2004,1),end=c(2013,12),frequency=12)
base

base.test=ts(plage[(n-b+1):n],start=c(2013,12),frequency=12) 
base.test

# Pour simplifier

Y=ts(base, start=c(2004,1),frequency=12)

# Stationnarité et saisonnalité

plot(Y,type="l",main="Représentation de la série Plage")

# test tendance

trend=1:length(Y)
Base=data.frame(Y,trend)
reg0=lm(Y~trend,data=Base)
summary(reg0)

#Coefficients:
#  Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  18.9651     2.5932   7.313 3.40e-11 ***
#  trend         0.1761     0.0372   4.734 6.18e-06 ***

# test saisonnalité

Mois=as.numeric(cycle(Y))
Base$Mois=as.numeric(cycle(Y))
head(Base)

reg0m=lm(Y~0+as.factor(Mois),data=Base)
summary(reg0m)

# graphiquement

win.graph(width=8, height=7,pointsize=7.25)
plot(Y,type="l",main="Représentation de la série Plage")
Month=c("J","F","M","A","M","J","J","A","S","O","N","D")
win.graph(width=8, height=7,pointsize=7.25)
plot(Y,type="l")
points(Y,pch=Month)

# Elimination de la tendance

tY=Y-predict(reg0) # correspond aux résidus de reg0
plot(tY)


# Elimination de la saisonnalité, par passage en différence saisonnière

dY=diff(tY,12)
plot(dY,type="l", main="Differentiation première saisonnière")

# Test de stationnarité de Dick-Fuller Augmenté

library(forecast)
library(tseries)

adf.test(dY)

# ACF et PACF de dY 

par(mfrow=c(1,2))
acf(dY,lag=50,main="dY",col="red",lwd=4)
pacf(dY,lag=50,main="dY",col="blue",lwd=4)

# Modèles

model1=arima(Y,order=c(1,1,0),seasonal=list(order=c(1,1,1),period=12)) 
acf(as.numeric(residuals(model1)),lwd=6,lag=36) 
pacf(as.numeric(residuals(model1)),lwd=6,lag=36)

# Test de bruit blanc

# 1)- test d'autocorrelation des résidus

Box.test(residuals(model1),type="Ljung-Box")
Box.test(residuals(model1),type="Box-Pierce")

# 2)- test de nulité de la moyenne des résidus

t.test(residuals(model1),mu=0,conf.level=0.95) 

# test de normalité des résidus

qqnorm(residuals(model1))
plot(density(residuals(model1)))
shapiro.test(residuals(model1))

# Prédiction et intervalle de confiance

Tnew=matrix(0,7,1)
Tnew[1]=2013+(12/12)
for(i in 2:7){Tnew[i]=2014+(i-2)/12}
Tnew

Ynew=matrix(0,7,1)
Ynew
for(i in 1:7){Ynew[i]=base.test[i]}
Ynew

#  Prédiction avec le modèle pour h=7, 7 mois en avant (année 2014)
# prévision
prevision=predict(model1,7)$pred
prevision


# Prévisions et observée

plot(Tnew,prevision,type="l",lwd=2,col="blue") 
lines(Tnew,Ynew, col="red") 
legend("topleft",c("prevision","valeurs exactes"),lty=1,col=c("blue","red"),bty="n")

# Lissage exponentiel double à deux paramètres Holt-Winter

HW=HoltWinters(Y)
Yhwpred=as.numeric(predict(HW,7)) 
plot(Tnew,Yhwpred,type="l",lwd=2,col="blue")
lines(Tnew,Ynew, col="red")
legend("topleft",c("hwpred","valeurs exactes"),lty=1,col=c("blue","red"),bty="n")