Discussion :

[geofnich]

Bonjour,

Je voudrais obtenir les intervalles de confiance des prévisions effectué à l'aide de processus ARMA, est ce possible?


Merci d'avance

[geofnich]

Aucune idée?

J'ai peut-etre mal posé ma question...


Merci de votre aide

[Hydro999]

Voici le code pour l'étude d'un processus ARMA.
Biensur vous etes invité à l'améliorer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
 
## Etude de série chronologique sur le pétrole
 
# source("R/TimeSerie/EtudeOil4.R")
 
# lien pour télécharger les données
lienOil ="http://research.stlouisfed.org/fred2/data/OILPRICE.csv"
# les données sont placés dans oil
oil <- read.table(lienOil,sep=",",dec=".",header=T)
 
# On prend que les valeurs
oil=oil$VALUE
 
# Calcul de Autocorrelation (ACF) et de l'Autocorrelation Partiel (PACF) 
# On fix l'interval de confiance à 0.95
IC=0.95
#acf(oil,ci=IC ,50)
 
# On differencie gnp
oilgr=diff(log(oil))
Loil=log(oil)
 
par(mfrow=c(2,1))
acf(oilgr,ci=IC ,30)
pacf(oilgr,ci=IC,30)
# on en déduit un ma(1) et un ar(1)
 
# Modèle
oilgr.ar=arima(oilgr,order=c(1,0,0))
oilgr.ma=arima(oilgr,order=c(0,0,1))
 
Loil.arma=arima(Loil,order=c(1,1,1))
Loil.arma2=arima(Loil,order=c(6,1,0))
 
# etude de normalité
quartz() # On ouvre une nouvelle fenetre
tsdiag(Loil.arma2,gof.lag=20)
 
quartz() # On ouvre une nouvelle fenetre
hist(Loil.arma2$resid, br=12)
quartz() # On ouvre une nouvelle fenetre
qqnorm(Loil.arma2$resid)
shapirotest=shapiro.test(Loil.arma2$resid)
 
# Prédiction
oil.pr=predict(Loil.arma2,n.ahead=10,se.value=.60)
IP = 0.70 # intervalle de prévision
U= oil.pr$pred + qnorm(IP)*oil.pr$se
L= oil.pr$pred - qnorm(IP)*oil.pr$se
 
#exp
oilPred=exp(oil.pr$pred)
U=exp(U)
L=exp(L)
 
 
month=(length(oil)-100):length(oil)
 
quartz() # On ouvre une nouvelle fenetre
#plot(month,oil[month],type="o",xlim=c(679,779),ylab="Oil",main="Spot Oil Price: West Texas Intermediate")
 
xlimD=length(oil)-100
 
xlimF=length(oil)+10
 
plot(month,oil[month],type="o",xlim=c(xlimD, xlimF),ylab="Oil",main="Spot Oil Price: West Texas Intermediate")
 
lines(oilPred,col="red",type="o")
lines(U,col="blue",lty="dashed")
lines(L,col="blue",lty="dashed")
V=length(oil)+1/2
abline(v=V, lty="dotted")

Il suffit de faire un copier coller dans un fichier : EtudeOil.R dans le répertoire qui va bien et de taper la commande

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

source("R/TimeSerie/EtudeOil.R")

Cdt

[geofnich]

Merci beaucoup pour les intervalles de confiance.

J'avais déjà un code pour la modélisation ARMA, je le poste si cela vous intérèsse.

[Hydro999]

Oui ça m'intéresse .
Je commence aussi à regarder les processus ARCH et Garch.

[geofnich]

Désolé pour le retard!
As tu un déjà commencé à programmer les ARCH/GARCH sous R?(je l'ai fais sous SAS masi cela m'intéresserait de le faire sous R...)

Voici le code pour un ARMA :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
 
#On efface tout!!!!! :
rm(list= ls())
library(tseries)
setwd("C:/temp")
 
 
Differenciation=("Ordre d de différenciation")
Ordre=("Ordre")
#############################################################################
# PARTIE 1 : CREATION DE LA SERIE : 
S1=rbind(63.49,49.82,87.10,95.11,94.61,86.35,87.49,75.10,68.36,34.4)
#############################################################################
#PARTIE 2 : TRANSFORMATION DE LA TABLE EN SERIE CHRONOLOGIQUE : 
S1 = ts(S1)
 
#############################################################################
#PARTIE 3 :STATIONNARISATION DE LA SERIE : 
#On passe au log et on différencie :
serie=log(S1)
for(i in 1:length(serie)){
	if (is.infinite(serie[i])){
		serie[i]=0}
}
 
v=-1
d=0
while(v==-1){
	kpss=kpss.test(serie)
	if (kpss$p.value<0.05) {serie=diff(serie,lag=1,differences=1)
					d=d+1}
	else {v=0}
	}
Differenciation=c(Differenciation,d)
 
#############################################################################
#PARTIE 4 : MODELISATION ARMA : B&J
 
 
	#Etape 1 : Identification du modèle
 
 
	op <- par(mfrow = c(2, 2))# 2 x 2 images
	#On dessine le graphe de la série : 
	plot(S1)
	title(main="Graphique de la série")
	#On plot la série( pour vérifier visuellement la stationnarité):
	plot(serie)
	title(main="Série stationnarisée en moyenne et en variance")
	#On trace la fonction d'autoco :
	corr=acf(serie,main="Fonction d'autocorrélation")
	#On trace la fonction d'autoco partielle : 
	autocorr=pacf(serie,main="Fonction d'autocorrélation partielle")
 
	#On détermine q à l'aide de la foncion d'autocorrélation :
	i=1
	q=-1
	while(q==-1){
		if (corr[i,]<0.2530344911 && corr[i,]>-0.2530344911)
			{q=i}
		i=i+1
	}
	#On détermine p à l'aide de la fonction d'autocorrélation partielle :
	A=ts(autocorr[[1]])
	print(A)
	paz=as.numeric(A)
	i=1
	p=-1
	while(p==-1){
		if (paz[i]<0.2530344911 && paz[i]>-0.2530344911)
			{p=i-1}
		i=i+1
	}
 
 
 
 
 
	#Etape 2 : Validation du modèle
 
 
		#Validation des paramètres : On effectue d'abord un test sur les paramètres pour detecter un superparamétrage du modèle : cela nous donnera plusieurs modèles
			#On commence par p :
			pp=p
			u=-1
			while(u==-1){
				test1=arima(serie,order= c(pp,0,0))
				tstat=abs(test1$coef[pp]/sqrt(abs(test1$var.coef[pp,pp])))
				if (abs(tstat)>1.97)
					{u=0}
				else {pp=(pp-1)}
				if (pp==0) {u=-2}
			}	
 
			#Puis q :
			u=-1
			qq=q
			while(u==-1){
				test1=arima(serie,order= c(0,0,qq))
				tstat=abs(test1$coef[qq]/sqrt(abs(test1$var.coef[qq,qq])))
				if (abs(tstat)>1.97)
					{u=0}
				else {qq=(qq-1)}
				if (qq==0) {u=-2}
			}
 
			#On test maintenant p et q
			ppp=p
			qqq=q
			w=-1
			while(w==-1)
			{	u=-1
				while(u==-1)
				{	test1=arima(serie,order= c(ppp,0,qqq))
					tstat1=abs(test1$coef[ppp]/sqrt(abs(test1$var.coef[ppp,ppp])))		
					tstat2=abs(test1$coef[qqq+ppp]/sqrt(abs(test1$var.coef[qqq+ppp,qqq+ppp])))
					if (abs(tstat1)>1.97 & abs(tstat2)>1.97) {u=0}
					else {qqq=(qqq-1)}
					if (qqq==0) {u=-2}
				}
				if (u==0) {w=0}
				if (w==-1) {ppp=(ppp-1)}
				if (w==-1) {qqq=q}
				if (ppp==0) {w=-2}
			}
			if (ppp==0) {qqq=0}
			if (qqq==0) {ppp=0}
 
			if(qqq<0){qqq=0}
			if(ppp<0){ppp=0}
 
 
 
 
 
		#Validation des résidus : 
			#Test d'autocorrélation de Box-Pierce :
				#Pour p :
				if (pp!=0){
				test1=arima(serie,order= c(pp,0,0))
				za=Box.test(test1$residuals, type = "Box-Pierce")
				if (za$p.value<0.05) {pp=0}
				}
				#Pour q
				if (qq!=0){
				test1=arima(serie,order= c(0,0,qq))
				za=Box.test(test1$residuals, type = "Box-Pierce")
				if (za$p.value<0.05) {qq=0}
				}
				#Pour ppp et qqq : 
				if (ppp!=0){
				test1=arima(serie,order= c(ppp,0,qqq))
				za=Box.test(test1$residuals, type = "Box-Pierce")
				if (za$p.value<0.05) {ppp=0 
							    qqq=0}
				}
 
 
		#Critère de choix : on minimise le BIC(=SIC) et le HQ(adoucit la pénalité du Bic : intermédiaire entre AIC et BIC);rq : HQ et SIC sont les seuls critères convergent!!!
			#Pour pp
			if (pp!=0){
			test1=arima(serie, order= c(pp,0,0))
			aic1 = log(test1$sigma2)+2*(pp)/length(S1)
			bic1 = log(test1$sigma2)+(pp)*log(length(S1))/length(S1)
			hq1=log(test1$sigma2)+(pp)*log(log(length(S1))/length(S1))
			print(aic1)
			print(bic1)
			print(hq1)
			}
			if (pp==0) {aic1=1000000
					  bic1=1000000
					  hq1=1000000}
 
			#Pour qq
			if (qq!=0){
			test1=arima(serie, order= c(0,0,qq))
			aic2 = log(test1$sigma2)+2*(qq)/length(S1)
			bic2 = log(test1$sigma2)+(qq)*log(length(S1))/length(S1)
			hq2=log(test1$sigma2)+(qq)*log(log(length(S1))/length(S1))
			print(aic2)
			print(bic2)
			print(hq2)
			}
			if (qq==0) {aic2=1000000
					  bic2=1000000
					  hq2=1000000}
			#Pour ppp et qqq
			if (ppp!=0){
			test1=arima(serie, order= c(ppp,0,qqq))
			aic3 = log(test1$sigma2)+2*(ppp+qqq)/length(S1)
			bic3 = log(test1$sigma2)+(ppp+qqq)*log(length(S1))/length(S1)
			hq3=log(test1$sigma2)+(ppp+qqq)*log(log(length(S1))/length(S1))
			print(aic3)
			print(bic3)
			print(hq3)
			}
			if (ppp==0) {aic3=1000000
					  bic3=1000000
					  hq3=1000000}
 
 
			#On choisit le min : 
			#min1=min(aic1,aic2,aic3) on s'en fout
			min2=min(bic1,bic2,bic3)#on regarde surtt celui la
			#min3=min(hq1,hq2,hq3) bien mais alpha choisit arbitrairement( c'est meme fau!!)à 1!!
			if (min2==bic1){p=pp
						q=0}
			if (min2==bic2){p=0 
						q=qq}
			if (min2==bic3){p=ppp 
						q=qqq}
 
 
 
	# Etape 3 : Estimation des paramètres et SSR
 
		#si p=0 et q=0 on fixe arbitrairement un ARMA(1,1) : normalement sa ne doit pas arriver!!
			if (p<=0 && q<=0){
				p=1
				q=1
			}
		test1=arima(serie, order= c(p,0,q))
		Ordre=c(Ordre,"p=",p,"q=",q)
 
	# Etape 4 : Prévisions et intervalles de confiances
 
 
		prev=predict(test1,n.ahead=5,se.value=.60)
		U= prev$pred + qnorm(0.95)*prev$se
		L= prev$pred - qnorm(0.95)*prev$se
		if(d==0){
			prevS1<-exp(prev$pred[1])
			IC1<-exp(U[1])
			IC2<-exp(L[1])
			for(i in (length(S1)+2):15){
				prevS1<-rbind(prevS1,exp(prev$pred[i-length(S1)]))
				IC1<-rbind(IC1,exp(U[i-length(S1)]))
				IC2<-rbind(IC2,exp(L[i-length(S1)]))
			}
		}
		if(d==1){
			prevS1<-S1[length(S1)]*exp(prev$pred[1])
			IC1<-S1[length(S1)]*exp(U[1])
			IC2<-S1[length(S1)]*exp(L[1])
			for(i in (length(S1)+2):15){
				prevS1<-rbind(prevS1,prevS1[i-(length(S1)+1)]*exp(prev$pred[i-length(S1)]))
				IC1<-rbind(IC1,IC1[i-(length(S1)+1)]*exp(U[i-length(S1)]))
				IC2<-rbind(IC2,IC2[i-(length(S1)+1)]*exp(L[i-length(S1)]))
			}
		}
		if(d==2){
			prevS1<-(S1[length(S1)]*S1[length(S1)]/S1[length(S1)-1])*exp(prev$pred[1])
			prevS1<-rbind(prevS1,prevS1[1]*prevS1[1]/S1[length(S1)])*exp(prev$pred[2])
			IC1<-(S1[length(S1)]*S1[length(S1)]/S1[length(S1)-1])*exp(U[1])
			IC1<-rbind(IC1,IC1[1]*IC1[1]/S1[length(S1)])*exp(U[2])
			IC2<-(S1[length(S1)]*S1[length(S1)]/S1[length(S1)-1])*exp(L[1])
			IC2<-rbind(IC2,IC2[1]*IC2[1]/S1[length(S1)])*exp(L[2])
			j=1
			for(i in (length(S1)+3):15){
				prevS1<-rbind(prevS1,(prevS1[j+1]*prevS1[j+1]/prevS1[j])*exp(prev$pred[i-length(S1)]))
				IC1<-rbind(IC1,(IC1[j+1]*IC1[j+1]/IC1[j])*exp(U[i-length(S1)]))
				IC2<-rbind(IC2,(IC2[j+1]*IC2[j+1]/IC2[j])*exp(L[i-length(S1)]))
				j=j+1
			}
		}
		if(d==3){
			j=1
			prevS1<-(S1[length(S1)]*S1[length(S1)]*S1[length(S1)]*S1[length(S1)-2]/(S1[length(S1)-1]*S1[length(S1)-1]*S1[length(S1)-1]))*exp(prev$pred[1])
			prevS1<-rbind(prevS1,prevS1[1]*prevS1[1]*prevS1[1]*S1[length(S1)-1]/(S1[length(S1)]*S1[length(S1)]*S1[length(S1)]))*exp(prev$pred[2])
			prevS1<-rbind(prevS1,prevS1[2]*prevS1[2]*prevS1[2]*S1[length(S1)]/(prevS1[1]*prevS1[1]*prevS1[1]))*exp(prev$pred[3])
			IC1<-(S1[length(S1)]*S1[length(S1)]*S1[length(S1)]*S1[length(S1)-2]/(S1[length(S1)-1]*S1[length(S1)-1]*S1[length(S1)-1]))*exp(U[1])
			IC1<-rbind(IC1,IC1[1]*IC1[1]*IC1[1]*S1[length(S1)-1]/(S1[length(S1)]*S1[length(S1)]*S1[length(S1)]))*exp(U[2])
			IC1<-rbind(IC1,IC1[2]*IC1[2]*IC1[2]*S1[length(S1)]/(IC1[1]*IC1[1]*IC1[1]))*exp(U[3])
			IC2<-(S1[length(S1)]*S1[length(S1)]*S1[length(S1)]*S1[length(S1)-2]/(S1[length(S1)-1]*S1[length(S1)-1]*S1[length(S1)-1]))*exp(L[1])
			IC2<-rbind(IC2,IC2[1]*IC2[1]*IC2[1]*S1[length(S1)-1]/(S1[length(S1)]*S1[length(S1)]*S1[length(S1)]))*exp(L[2])
			IC2<-rbind(IC2,IC2[2]*IC2[2]*IC2[2]*S1[length(S1)]/(IC2[1]*IC2[1]*IC2[1]))*exp(L[3])
			for(i in (length(S1)+4):130){
				prevS1<-rbind(prevS1,prevS1[j+2]*prevS1[j+2]*prevS1[j+2]*prevS1[j]/(prevS1[j+1]*prevS1[j+1]*prevS1[j+1])*exp(prev$pred[i-length(S1)]))
				IC1<-rbind(IC1,IC1[j+2]*IC1[j+2]*IC1[j+2]*IC1[j]/(IC1[j+1]*IC1[j+1]*IC1[j+1])*exp(U[i-length(S1)]))
				IC2<-rbind(IC2,IC2[j+2]*IC2[j+2]*IC2[j+2]*IC2[j]/(IC2[j+1]*IC2[j+1]*IC2[j+1])*exp(L[i-length(S1)]))
 
				j=j+1	
			}
		}	
	#On concatene les valeurs de la serie avec leur prévisions
	Prevision=rbind(S1,prevS1)
	IC1=rbind(S1,IC1)
	IC2=rbind(S1,IC2)
	Prevision=cbind(IC2,Prevision,IC1,Ordre,Differenciation)
 
 
	KKK=paste("Prevision",".csv",sep="")
	write.table(x = Prevision, file = KKK, sep = ";", row.names = FALSE)

[Hydro999]

Merci pour le code.
De mon coté j'ai trouvé cette fonction donnant le meilleur ARMA possible :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
get.best.arima <- function(x.ts, maxord = c(1,1,1,1,1,1))
	{
	best.aic <- 1e8
	n <- length(x.ts)
	for (p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
	for (P in 0:maxord[4]) for(D in 0:maxord[5]) for(Q in 0:maxord[6])
		{
		fit <- arima(x.ts, order = c(p,d,q),seas = list(order = c(P,D,Q),frequency(x.ts)), method = "CSS")
		fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef)
		if (fit.aic < best.aic)
			{
			best.aic <- fit.aic
			best.fit <- fit
			best.model <- c(p,d,q,P,D,Q)
			}
		}
		list(best.aic, best.fit, best.model)
	}

Tu utilise la fonction comme ça

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
best.arima.oil <- get.best.arima( log(Oil.ts),maxord = c(2,2,2,2,2,2))

Pour les arch garch j'utilise la librairie

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
library(tseries)

J'ai commencé mais je ne suis pas trés à l'aise pour le moment sur ce sujet.