|
require(caschrono)
|
|
require(chron)
|
|
require(urca)
|
|
|
|
#--------------------------------------------------
|
|
# Entree donnees
|
|
#--------------------------------------------------
|
|
qj<-read.table("q_Nosiarivo_jour_8199.txt",,dec=".") #débit journalier
|
|
qm<-read.table("q_Nosiarivo_mois_8199.txt",,dec=".") #débit mensuel
|
|
p<-read.table("p_Toliara_6014.txt",,dec=".") #pluie
|
|
|
|
qm8183<-read.table("q_Nosiarivo_mois_8183.txt",,dec=".") #Années 80
|
|
qm9499<-read.table("q_Nosiarivo_mois_9499.txt",,dec=".") #Années 90
|
|
qm9499_fill<-read.table("q_Nosiarivo_mois_9499_fill.txt",,dec=".") #Années 90 _ on a rempli les données manquantes par trend + seasonality
|
|
qm9499_fill2<-read.table("q_Nosiarivo_mois_9499_fill2.txt",,dec=".") #Années 90 _ on a rempli les données manquantes par mean + seasonality
|
|
|
|
|
|
Tmin<-read.table("T(min)_Toliara6012.txt",,dec=".") #Températures min et max
|
|
Tmax<-read.table("T(max)_Toliara6012.txt",,dec=".") #
|
|
|
|
# Conversion en séries temporelles
|
|
p.ts=ts(p,start=c(1960,1),frequency=12)
|
|
qm.ts=ts(qm,start=c(1981,8),frequency=12)
|
|
qj.ts=ts(qj,start=c(1981,213),frequency=365)
|
|
|
|
qm8.ts=ts(qm8183,start=c(1981,9),frequency=12)
|
|
qm9.ts=ts(qm9499,start=c(1994,1),frequency=12)
|
|
qm9_fill.ts=ts(qm9499_fill,start=c(1994,1),frequency=12)
|
|
qm9_fill2.ts=ts(qm9499_fill2,start=c(1994,1),frequency=12)
|
|
|
|
Tmin.ts=ts(Tmin,start=c(1960,1),frequency=12)
|
|
Tmax.ts=ts(Tmax,start=c(1960,1),frequency=12)
|
|
|
|
|
|
#--------------------------------------------------
|
|
# Visualisation
|
|
#--------------------------------------------------
|
|
|
|
|
|
# Vue d'ensemble : Chronogramme
|
|
plot.ts(p.ts,xlab='année',ylab='pluviométrie(mm)')
|
|
plot.ts(qm.ts,xlab='année',ylab='Débit(m3/s)')
|
|
plot.ts(qj.ts,xlab='année',ylab='Débit(m3/s)')
|
|
|
|
plot.ts(qm8.ts,xlab='année',ylab='Débit(m3/s)')
|
|
plot.ts(qm9.ts,xlab='année',ylab='Débit(m3/s)')
|
|
|
|
|
|
# Aperçu de la saisonalité : monthplot
|
|
monthplot(p.ts,xlab='année',ylab='pluviométrie(mm)')
|
|
monthplot(qm.ts,xlab='année',ylab='Débit(m3/s)')
|
|
|
|
monthplot(qm8.ts,xlab='années 80',ylab='Débit(m3/s)')
|
|
monthplot(qm9.ts,xlab='années 90',ylab='Débit(m3/s)')
|
|
|
|
# Diagramme retardée - lagpplot
|
|
|
|
lag.plot(rev(qm8.ts),12,layout=c(4,3),do.lines=FALSE,diag.col="red",col.main="blue")
|
|
|
|
#Décomposition saisonnière de la série
|
|
|
|
decompose(p.ts)
|
|
|
|
#---------------------------------------------------------
|
|
# premiers mod?les simples (lissage exponentiel)
|
|
#---------------------------------------------------------
|
|
|
|
#Prévision par lissage exponentiel
|
|
|
|
p.ets=ets(p.ts,model="ANA") #prédiction à l'horizon 12
|
|
p.predict=predict(p.ets,12)
|
|
|
|
xlisse <- HoltWinters(qm9_fill.ts,seasonal="add")
|
|
xlisse_predict<-predict(xlisse,n.ahead=12)
|
|
|
|
plot((xlisse, n.ahead = 12, col = "black", ylab = object$series, lty = 2, n1, newxreg, transform, Plot=TRUE, ...))
|
|
|
|
#Prévision par arma
|
|
|
|
fit <- Arima(qm9_fill.ts,order=c(3,1,0))
|
|
|
|
plot(fit$x,col="blue")
|
|
lines(fitted(fit),col="red")
|
|
|
|
|
|
#--------------------------------------------------------------
|
|
# Analyse des s?ries et estimation des mod?les
|
|
#--------------------------------------------------------------
|
|
|
|
# test de stationnarit?
|
|
# KPSS test
|
|
|
|
# test pour savoir si y est stationnaire de moyenne constante, option type = "mu" dans ur.kpss()
|
|
|
|
summary(ur.kpss(p.ts,type="mu"))
|
|
summary(ur.kpss(qm9_fill.ts,type="mu"))
|
|
|
|
|
|
# test pour savoir si y est stationnaire ? une tendance lin?aire pr?s, option type = "tau" dans ur.kpss()
|
|
|
|
summary(ur.kpss(p.ts,type="tau"))
|
|
summary(ur.kpss(qm9_fill.ts,type="tau"))
|
|
|
|
|
|
# ACF de la pluie et ACF de la s?rie diff?renci?e de la pluie
|
|
|
|
plot2acf(p.ts,diff(p.ts,12),lag.max=40,main=c("pluie",expression(paste("(1-",B^{12},") pluie", sep =""))))
|
|
plot2acf(p.ts,diff(p.ts,12),lag.max=80,main=c("pluie Toliara",expression(paste("(1-",B^{12},") pluie", sep =""))))
|
|
|
|
|
|
# ACF des d?bits et de sa s?rie diff?renci? ? l'ordre 12
|
|
|
|
plot2acf(qm9_fill.ts,diff(p.ts,12),lag.max=40,main=c("qm9_fill",expression(paste("(1-",B^{12},") qm9_fill", sep =""))))
|
|
|
|
# ACF des d?bits et de la pluie
|
|
|
|
plot2acf(p.ts ,qm9_fill.ts, lag.max=40,main=c("ACF_pluie_toliara", "ACF_qm9_fill_Nosiarivo",sep =""))
|
|
|
|
# Identification du mod?le ARMA
|
|
|
|
armaselect(p.ts,nbmod=4) # s?lection des ordres de l'ARMA, pas adapt? pour les chroniques saisonniers.
|
|
|
|
mod0=auto.arima(p.ts)
|
|
|
|
ret=c(6,12,18,24,30)
|
|
Box.test.2(residuals(mod0),nlag=ret,type="Ljung-Box",decim=2) # test de blancheur du r?sidu
|
|
|
|
modp1=Arima(p.ts,order=c(1,0,1),seasonal=list(order=c(1,0,1),period=12))
|
|
Box.test.2(residuals(modp1),nlag=ret,type="Ljung-Box",decim=2)
|
|
|
|
t_stat(modp1) # test de significativit? des coefficients
|
|
|
|
modp2=auto.arima(p.ts, D=1)
|
|
Box.test.2(residuals(modp2),nlag=ret,type="Ljung-Box",decim=2)
|
|
t_stat(modp2)
|
|
|
|
modq0=auto.arima(qm9_fill.ts)
|
|
Box.test.2(residuals(modq0),nlag=ret,type="Ljung-Box",decim=2)
|
|
t_stat(modq0)
|
|
|
|
modq1=Arima(qm9_fill.ts,seasonal=list(order=c(1,1,0),period=12),include.drift=TRUE)
|
|
|
|
# voir le graphe du mod?le
|
|
|
|
plot(modq0$x,col="blue")
|
|
lines(fitted(modq0),col="red")
|
|
|
|
fitnot=Arima(p,order=c(1,0,0),list(order=c(2,1,0),period=12))
|
|
|
|
summary(fitnot)
|
|
|
|
help(Arima)
|
|
|
|
|
|
|
|
#-------------------
|
|
# pr?vision
|
|
#-------------------
|
|
|
|
pred.modq1=forecast(modq1,h=36)
|
|
|
|
|
|
|
|
|
|
|
|
|