Codes/Figures                        Données                         Compléments                     Erratum

Les compléments R par chapitre

• Comp1. Démarche de base en ST
• Comp2. R pour les ST
• Comp3. Rappels de statistique
• Comp4. Modèles de base en ST
• Comp5. ST non stationnaires
• Comp6. Lissage exponentiel
• Comp7. Simulation
• Comp8. Aéroport Toulouse-Blagnac
• Comp9. Température à Nottingham
• Comp10. Consommation d'éléctricité
• Comp11. Production de lait
• Comp12. Modèle ARCH/GARCH

# 12.3 ARCH(1) et test d'hétéroscédasticité conditionnelle
# Exemple 12.1

require("fGarch")
spec.1 = garchSpec(model = list(mu=5,omega= 0.05, alpha = 0.9, beta = 0),rseed = 397)
archsim.1 =   garchSim(extended = TRUE,spec.1, n = 300)
mar.p = c(4.1,4,1,1)
op<-par(oma=rep(0,4),cex.lab=1.2)
par(mfrow=c(2,2))
par(mar=mar.p)
plot.ts(archsim.1[,1], xlab="temps", ylab="série")
#
par(mar=mar.p)
plot.ts(archsim.1[,2], xlab="temps", ylab="écart-type conditionnel")
#
par(mar=mar.p)
acf(archsim.1[,1],lag.max=30,xlab="retard", ylab="acf de y", main="")
v = (archsim.1[,1]-5)^2
#
par(mar=mar.p)
acf(v, lag.max=30, main="",xlab="retard", ylab=expression(paste("acf de ",(y-5)^2,sep="")))
par(op)


# simulation de 100 trajectoires de 20 valeurs du GARCH(2,1) 
nsim=100
spec=garchSpec(model=list(mu=2,omega=0.09,alpha =c(0.15, 0.3),
  beta = 0.4), rseed=9647)
var.margi = 0.09/(1 - 0.15 - 0.3-0.4)
y = garchSim(spec, n = 420, extended = TRUE)
y1 = y[21:420,1]
m1 = mean(y1)
(q5 = quantile(abs(y1-m1), probs=c(.975,.98,.985,.99,.995)))
extrem985 = which(abs(y1-m1) > q5[3])
cat("nombre de points : ",length(extrem985),"\n")
y1b = y1
y1b[extrem985]  = q5[3]* sign(y1[extrem985]-m1)
mod1=garchFit(~garch(2,1),data=y1b,trace=FALSE,include.mean=TRUE)
summary(mod1)



interv = 398:400
presamp = cbind(residuals(mod1)[interv], mod1@h.t[interv],y1b[interv,1])

# spécification du modèle à partir des résultats d'estimation
spec  = garchSpec(model = list(mu = mod1@fit$par[1], omega= mod1@fit$par[2], alpha =mod1@fit$par[c(3,4)] , beta =  mod1@fit$par[5]), presample= presamp)
# tableau des simulations
ysim = matrix(0,nrow=20,ncol=nsim)
for(isim in 1:nsim)
{
ysim[,isim] = garchSim(spec, n = 20)
}
# bande à 95%
# repérage pour chaque horizon des quantiles empiriques d'ordres 2.5% et 97.5%
p.cond=predict(mod1, mse="cond", n.ahead=20, plot=FALSE)
b.inf= p.cond$meanForecast -1.96*p.cond$standardDeviation
b.sup= p.cond$meanForecast +1.96*p.cond$standardDeviation
matpr= cbind(y[401:420,1],b.inf,b.sup)
bornes = apply(ysim,1,quantile, probs = c(0.025,0.975))
matplot(1:20,cbind(matpr,t(bornes)),type='l', lty=c(1,2,2,3,3), col=c(1,2,2,3,3))


# Représentation de la Fig. 12.7
density.plot=function(x,legende=FALSE,...)
{
H<-hist(x, sub = NULL, ylab = "densité", freq=FALSE, ...)
abline(v=0,lwd=2)
rug(x, ticksize = 0.01)
xmin=par()$usr[1];xmax=par()$usr[2]
tab<-seq(xmin,xmax,0.1)
lines(tab, dnorm(tab, mean(x),sd(x)), col="red", lty=2, lwd=2)
lines(density(x),lwd=2, col="orange")
if(legende)
 {lg0 =c("estimation n.p. de la densité","estimation d'une   gaussienne")
 legend("topright",legend=lg0,lty=c(1,2),lwd=2,
 col=c("orange","red"),cex=0.9)
 }
}

mar.p = c(4.1,4,2.5,1)
op<-par(oma=rep(0,4),cex.lab=1.2,mar=mar.p)
par(mfcol=c(2,2))
# ou alors
rr= rendav.12[,4][complete.cases(rendav.12) == TRUE]
density.plot(rr, main="L'Oréal avant", xlab = "rendement", ylim=c(0,0.33),
nclass=12, xlim=c(-10,10))
#
rr= rendapr.12[,4][complete.cases(rendapr.12) == TRUE]
density.plot(rr, main="L'Oréal pendant", xlab = "rendement", ylim=c(0,0.33),
nclass=18, xlim=c(-10,10))
#
rr= rendav.06[,2][complete.cases(rendav.06) == TRUE]
density.plot(rr, main="Société Générale avant", xlab = "rendement", ylim=c(0,0.35),
nclass=12, xlim=c(-12,12))
#
rr= rendapr.06[,2][complete.cases(rendapr.06) == TRUE]
density.plot(rr, main="Société Générale pendant", xlab = "rendement", ylim=c(0,0.35),
nclass=30, xlim=c(-12,12))
par(op)


# Tableau 12.4
xtable(hetero,caption="Test d'hétéroscédasticité conditionnelle avant et pendant la crise : p-value",label="hetero.avap",digits=4)