%\VignetteIndexEntry{Annexe Chapitre 12}
%\VignetteDepends{}
%\VignetteKeywords{ts}
%\VignettePackage{caschrono}
\documentclass{article}
\usepackage{Sweave}
\usepackage[latin1]{inputenc}
\usepackage{times}
\usepackage{mathptm}
\usepackage{hyperref}
\usepackage{natbib}
\usepackage{pratiquer}
\setkeys{Gin}{width=0.95\textwidth}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}

\DeclareMathOperator*{\argmax}{argmax}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\BB}{BB}
\DeclareMathOperator{\BBN}{BBN}
\DeclareMathOperator{\rang}{rang}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\V}{V}
\DeclareMathOperator{\C}{Cov}
\DeclareMathOperator{\Prob}{Pr}
\DeclareMathOperator{\diag}{diag}
\def\agra{\boldsymbol{a}} % a grave
%
\def\alg{\boldsymbol{\alpha}}
\def\b{\mbox{\bf b}}
%
\def\B{\mbox{\bf B}}
\def\BB#1{\mbox{BB}(0,\sigma_{#1}^2)}
\def\bm{\mbox{B}}
\def\betg{\boldsymbol{\beta}}
\def\bar#1#2{\overline{#1}_{#2}} % surlign\'e indic\'e
\def\bgr{\mbox{\bf b}}
\def\bla{\mbox{~}}
\def\blb{\mbox{~~}}
%
\def\card{\mbox{card}}
\def\corr{\mbox{corr}}
\def\cov{\mbox{\sf cov}}
\def\cv{\mbox{CV}}
%
\def\D{\mbox{\bf D}}
\def\d{\bf d}
\def\degr{^{\circ}}
\def\deltag{\boldsymbol{\delta}}
\def\Deltag{\boldsymbol{\Delta}} % ligne 64
\def\diag{\mbox{diag}}
\def\dir{\mbox{\tiny{DIR}}}
\def\dsp{\displaystyle}
%
\def\e{\mbox{\bf e}}
\def\eps{{\bf \epsilon}}
\def\esp{\mbox{\sf E}}
\def\eti{\tilde{\epsilon}}

\def\F{\mbox{\bf F}}
\def\g{\mbox{\bf g}}
\def\gamg{\boldsymbol{\Gamma}}
\def\gamgp{\boldsymbol{\gamma}}
\def\hors{\mathbin{\in\mkern-12mu/}}

\def\I{\mbox{\bf I}}
\def\indic{\mathrm{I\mkern-8muI}}
\def\ie{c'est-\`a-dire}
\def\iid{\sim_{\mbox{i.i.d.}}}
\def\lm{\mbox{L}}

\def\mapsous#1{\smash{ \mathop{\longrightarrow}\limits_{#1}}} %
\def\mapsur#1{\smash{ \mathop{\longrightarrow}\limits^{#1}}}  %
\def\M{\mbox{\bf M}}
\def\mug{\boldsymbol{\mu}}
\def\nor{\mathcal{N}}
\def\nug{\boldsymbol{\nu}}
\def\Nx{\mbox{\bf N}}
\def\P{\mbox{\bf P}}
\def\Phig{\boldsymbol{\Phi}}
\def\phig{\boldsymbol{\phi}}
\def\poi{\mbox{\cal P}}
\def\pr{\mbox{Pr}}
\def\prim{^{\boldsymbol{\prime}}}
\def\px{\mbox{\bf x}}

\def\Rho{\boldsymbol{\rho}}
\def\rl{\mathbin{I\mkern-8muR}} % Reel
\def\RR{\textsf{R}\/}

\newcommand{\sig}[2]{\Sigma_{{#1},{#2}}}
\def \sitst{\textsf{SiteST}\/}
\def\SP{\texttt{S-PLUS}\/}
\def\T{\mbox{\bf T}}
\def\tr{\triangle}
\def\t{\mbox{\bf t}}
\def\tra{\mbox{tr}}

\def\U{\mbox{\bf U}}
\def\Unif{\emph{Unif}}
\def\u{\mbox{\bf u}}
\def\v{\mbox{\bf v}}

\def\w{\mbox{\bf w}}
\def\var{\mbox{\sf var}}
\def\W{\mbox{\bf W}}
\def\X{\textbf{X}}
\def\x{\textbf{x}}

\def\Y{\textbf{Y}}
\def\yg{\textbf{y}}
\def\y0{y^0}

\def\Z{\textbf{Z}}
\def\z{ \textbf{z}}
\def\zer{\large{0}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\urlstyle{sf}
\def\rmdefault{cmr}

%%%% pour les chiffres des \'equations
\makeatletter
\renewcommand\theequation{\thesection.\arabic{equation}}
\@addtoreset{equation}{section}
\makeatother

\title{H\'et\'erosc\'edasticit\'e conditionnelle (compl\'ements du Chapitre 12)}
\author{Yves Aragon\footnote{yves.aragon@gmail.com} \cr
{\normalsize Universit\'e Toulouse Capitole} }
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\setcounter{section}{12}
\renewcommand{\thefootnote}{\arabic{footnote}}

\SweaveOpts{keep.source=TRUE}
<<echo=FALSE>>=
owidth <- getOption("width") # largeur des sorties
options(width=60, continue="+ ","warn"=-1 )
.PngNo <- 0
nom.fich = "anx12-bitmap-"
@
%
@
%
\SweaveOpts{keep.source=TRUE}
%
% les diff\'erents types de graphiques
%  ancien
<<label=bfig,echo=FALSE,eval=FALSE>>=
.PngNo <- .PngNo + 1; file = paste(nom.fich, .PngNo, sep="")
pdf(file=paste(file,".pdf",sep=""), width = 7, height = 7, pointsize = 12, bg = "white")
@

<<label=bfigps,echo=FALSE,eval=FALSE>>=
postscript(file=paste(file,".ps",sep=""), width = 7, height = 7, pointsize = 12, bg = "white",horizontal= FALSE,paper="special")
@


% 1111111111111
<<label=bfig1,echo=FALSE,eval=FALSE>>=
.PngNo <- .PngNo + 1; file = paste(nom.fich, .PngNo, sep="")
pdf(file=paste(file,".pdf",sep=""), width = 5, height = 2, pointsize = 10, bg = "white")
@

<<label=bfigps1,echo=FALSE,eval=FALSE>>=
postscript(file=paste(file,".ps",sep=""),  width = 5, height =2, pointsize = 10, bg = "white",horizontal= FALSE,paper="special")
@
% 222222222222222
<<label=bfig2,echo=FALSE,eval=FALSE>>=
.PngNo <- .PngNo + 1; file = paste(nom.fich, .PngNo, sep="")
pdf(file=paste(file,".pdf",sep=""), width = 3.9, height = 3.1, pointsize = 10, bg = "white")
@

<<label=bfigps2,echo=FALSE,eval=FALSE>>=
postscript(file=paste(file,".ps",sep=""), width = 3.9, height = 3.1,   pointsize = 10, bg = "white",horizontal= FALSE,paper="special")
@
%   3333333333333333333333333333

<<label=bfig3,echo=FALSE,eval=FALSE>>=
.PngNo <- .PngNo + 1; file = paste(nom.fich, .PngNo, sep="")
pdf(file=paste(file,".pdf",sep=""), width = 5.92, height = 6.74, pointsize = 10, bg = "white")
@
<<label=bfigps3,echo=FALSE,eval=FALSE>>=
postscript(file=paste(file,".ps",sep=""), width = 5.92, height = 6.74, pointsize = 10, bg = "white",horizontal= FALSE,paper="special")
@

<<label=bfig4,echo=FALSE,eval=FALSE>>=
.PngNo <- .PngNo + 1; file = paste(nom.fich, .PngNo, sep="")
pdf(file=paste(file,".pdf",sep=""), width = 6, height = 6, pointsize = 10, bg = "white")
@
<<label=bfigps4,echo=FALSE,eval=FALSE>>=
postscript(file=paste(file,".ps",sep=""), width = 6, height = 6, pointsize = 10, bg = "white",horizontal= FALSE,paper="special")
@

%5555555555555555
<<label=bfig5,echo=FALSE,eval=FALSE>>=
.PngNo <- .PngNo + 1; file = paste(nom.fich, .PngNo, sep="")
pdf(file=paste(file,".pdf",sep=""), width = 7, height = 4, pointsize = 10, bg = "white")
@

<<label=bfigps5,echo=FALSE,eval=FALSE>>=
postscript(file=paste(file,".ps",sep=""),  width = 7, height =4, pointsize = 10, bg = "white",horizontal= FALSE,paper="special")
@


<<label=zfig2,echo=FALSE,eval=FALSE>>=
dev.null <- dev.off()
@

<<label=zfiginclude,echo=FALSE,eval=FALSE>>=
cat("\\includegraphics[width=0.9\\textwidth]{", file, "}\n\n", sep="")
@


\begin{Exercice}
Faire appara\^itre A partir du  rendement  compos\'e quotidien $r_t = \log(x_t) - \log(x_{t-1})$ calculer le rendement compos\'e sur une semaine
de 5 jours.
\end{Exercice}

\textbf{R\'eponse.}
Le rendement compos\'e sur 5 jours est
\begin{multline}
 r_{5,t} = \log(x_t) - \log(x_{t-5}) \\
 = (\log(x_t) - \log(x_{t-1})) + ( \log(x_{t-1})-   \log(x_{t-2})+\cdots +(\log(x_{t-4})-   \log(x_{t-5}))
 \end{multline}
il est donc la somme des rendements quotidiens des jours de la p\'eriode.



\begin{Exercice}
\begin{itemize}
  \item Simuler 300 observations de $y_t $ ob\'eissant  \`a  :
  \begin{eqnarray}\label{sim.arch1}
    y_t = 5 +  \epsilon_t   \\
z^2   \epsilon_t = \sigma_t z_t,\; z_t \thicksim \BBN(0,\sigma_z^2) \\
    \sigma^2_t = 0.1 + 0.9 \: \epsilon_{t-1}^2
  \end{eqnarray}
(On simulera 10 valeurs
avant de commencer \`a stocker la s\'erie et on initialisera  \`a 0 la premi\`ere
variance conditionnelle.) On reproduira 
explicitement les \'equations. Visualiser la s\'erie et recommencer les simulations avec diff\'erentes graines. 
  \item Tester la blancheur du carr\'e centr\'e de la s\'erie. Commenter.
  
\end{itemize}

\end{Exercice} 

\textbf{R\'eponse. } 
Il nous faut essentiellement : (1) tirer les valeurs de $z_t$ i.i.d. $\nor(0,1)$, (2) calculer les variances conditionnelles par r\'ecurrence. Apr\`es quoi on peut   calculer $y_t$.

{\small
<<label=sim.manu>>=
require(caschrono)
set.seed(1618) ; nobs = 300
zt = rnorm(nobs+10)
# initialisation
sig2 = rep(0, nobs+10)
epsil = sig2
y = sig2
for(i in 2:(nobs+10))
 {sig2[i] <- 0.1 + 0.9*epsil[i-1]^2
  epsil[i] <- sig2[i]^.5 *zt[i]
  y[i] <- 5+epsil[i] 
 }
y <- y[-(1:10)]
@
}
Pour obtenir les graphiques on lancera le code

{\small
\begin{verbatim}
plot.ts(y)
acf(y)
\end{verbatim}
}

\begin{Exercice} [Test d'h\'et\'erosc\'edasticit\'e conditionnelle bas\'e sur un test de blancheur]
Tester   l'h\'et\'erosc\'edasticit\'e conditionnelle de ces deux s\'eries par un test de blancheur de Box-Pierce.

\end{Exercice}

\textbf{R\'eponse.}
Les deux s\'eries en question sont le processus ARCH(1) simul\'e et le bruit blanc \`a la base, colonnes 1 et  3 dans
\code{archsim.1} simul\'e section 12.3.1 suivant le mod\`ele \ref{sim.arch1}. On en refait la simulation ici :


{\small
<<>>=
require(fGarch)
spec.1 <- garchSpec(model = list(mu = 5, omega = 0.1, 
                                 alpha = 0.9, beta = 0),
                    rseed = 397)
archsim.1 <- garchSim(extended = TRUE, spec.1, n = 300, n.start=10)
head(archsim.1,2)
@
}


On a not\'e que l'h\'et\'erosc\'edasticit\'e conditionnelle 
d'un bruit blanc $x_t$,  se traduit par la non blancheur de son carr\'e,  un test de blancheur du carr\'e de la s\'erie est
donc  un test d'homosc\'edasticit\'e de la s\'erie. La blancheur de l'ARCH(1) simul\'e dans le livre et celle de la s\'erie $z_t$ sous-jacente
s'obtiennent par :

{\small
<<label=test.hetero.2>>=
moyenne <- mean(archsim.1[,1])
ret <- c(6, 12, 18)
a1 <- Box.test.2(archsim.1[,1], ret, 
                 type = "Box-Pierce", decim = 8)
a2 <- Box.test.2((archsim.1[,1]-moyenne)^2, ret, 
                 type = "Box-Pierce", decim = 8)
a3 <- Box.test.2(archsim.1[,3]^2, ret, 
                 type = "Box-Pierce", decim = 8)
a123 <- cbind(a1, a2[,2], a3[,2])
colnames(a123) <-
   c("Retard", "p.val y1", "p.val y1 centre carre", "p.val z_t")
@
}
On centr\'e la s\'erie simul\'ee car elle n'est pas de moyenne nulle.  
{\small
<<label=hetero.tex, results=tex, echo=FALSE>>=
require(xtable)
xtable(a123 , caption = "Test de blancheur, p-values", label = "blanc.y", 
 digits = 4)
@
}
\noindent
On a test\'e que les coefficients d'autocorr\'elation  jusqu'aux ordres 6, 12, 18 sont nuls, pour la s\'erie simul\'ee \code{y1}, son carr\'e centr\'e,
et pour la s\'erie du bruit blanc gaussien servant \`a engendrer l'ARCH. On voit (table \ref{blanc.y}) qu'on accepte
l'hypoth\`ese de blancheur pour la s\'erie simul\'ee et  pour le bruit
blanc gaussien sous-jacent, mais qu'on la rejette \'evidemment pour le carr\'e de la s\'erie centr\'ee.


\begin{exercice}[Tests sur \code{mod1}]
L'estimation de (\ref{arch1.sim}) a \'et\'e obtenue dans \code{mod1} :
\begin{eqnarray}  \label{arch1.sim}
    y_t &= & 2+ \sigma_t z_t \\
       \sigma^2_t& =& 0.09 + 0.15 \epsilon_{t-1}^2 +  0.3 \epsilon_{t-2}^2  + 0.4 \sigma^2_{t-1}. \notag
 \end{eqnarray}
 Faisons l'hypoth\`ese que les estimateurs des param\`etres sont approximativement 
normalement distribu\'es. 
\begin{itemize}
\item  Extraire de la sortie de l'estimation,  la matrice des variances-covariances estim\'ees des param\`etres.
  \item Tester l'hypoth\`ese que $\alpha_2 = 0.3$.
  \item Tester l'hypoth\`ese que $(\alpha_2, \beta_1)  = (0.3, 0.4)$
 \item Conclure.
\end{itemize}

\end{exercice}

\textbf{R\'eponse.}
D'abord nous simulons la m\^eme s\'erie que dans le chapitre :


{\small
<<label=sim.garch>>=
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]
@
}
\noindent
La s\'erie \`a \'etudier est form\'ee des 400 derni\`eres   observations   contenues dans la premi\`ere colonne de \verb"y".
Avant d'ajuster  un  mod\`ele  \`a cette s\'erie nous avons examin\'e son chronogramme :

\begin{rcode}
plot.ts(y1, xlab='temps')
\end{rcode}


Il y a plusieurs valeurs extr\^emes que  nous d\'ecidons de les remplacer par des valeurs plus raisonnable. 

{\small
<<label=y1.ext1, eval=TRUE>>=
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)
@
}
\noindent
Nous travaillons sur la s\'erie corrig\'ee \code{y1b}.

\begin{itemize}
\item  Extraction de la matrice des variances-covariances estim\'ees des param\`etres.
L'estimation est conduite par :

{\small
<<label=y1.ext2, eval=TRUE>>=
mod1 <- garchFit(~garch(2,1), data = y1b, trace = FALSE, 
                 include.mean = TRUE)
@
}
On situe les statistiques qui nous int\'eressent par \code{str(mod1)}. On y r\'ecup\`ere aussi les noms pr\'ecis des param\`etres. 
L'estimation des param\`etres alpha2 et beta1 est extraite de \code{mod1} par : 

{\small
<<>>=
param = c("alpha2", "beta1")
estim = as.matrix(coef(mod1)[param])
@
}
et celle de leur matrice des covariances estim\'ee par :

{\small
<<>>=
# matrice des covariances des estimateurs
vcov.par = mod1@fit$cvar[param,param]
@
}
On peut noter \'egalement que \verb"mod1@fit$se.coef" contient les \'ecarts-types d'estimation, {\ie}~la racine carr\'ee de la diagonale de
la matrice des covariances.

 \item Tester l'hypoth\`ese H$_0$ : $\alpha_2 = 0.3$.

 Nous devons former la statistique :
$T = \frac{\widehat{\alpha}_2 - 0.3}{s_{\widehat{\alpha}_2}}$ (cf. chap.3). 
Sous H$_0$, elle  est approximativement normalement distribu\'ee $\nor(0,1)$.
 
{\small
<<label=test.par.2>>=
t.alpha = (coef(mod1)["alpha2"]-0.3)/mod1@fit$se.coef["alpha2"]
p.val = 2* pnorm(-abs(t.alpha))
cat("p-value : ", p.val, "\n")
@
}
\noindent
 La p-value est   \'elev\'ee,
et donc l'estimation est coh\'erente avec le mod\`ele simul\'e.


\item Tester l'hypoth\`ese que $(\alpha_2, \beta_1)  = (0.3, 0.4)$.

Sous l'hypoth\`ese nulle :
\[
\mbox{H}_0 : \begin{bmatrix}
\alpha_2 \\
\beta_1
\end{bmatrix} = \begin{bmatrix}
0.3 \\
0.4
\end{bmatrix},\quad \begin{bmatrix}
\widehat{\alpha}_2 \\
\widehat{\beta_1}
\end{bmatrix}  \sim \mbox{~AN}(\begin{bmatrix}
0.3 \\
0.4
\end{bmatrix},\widehat{V}),
\]
o\`u  $\widehat{V}$ est la matrice des covariances estim\'ee pour ce couple d'estimateurs et
\[
[\widehat{\alpha}_2 - 0.3 , \widehat{\beta_1} - 0.4] \widehat{V}^{-1} \begin{bmatrix}
\widehat{\alpha}_2 -0.3 \\
\widehat{\beta_1} -0.4
\end{bmatrix}  \sim \mbox{~A}\chi^2_2
\]
On rejette l'hypoth\`ese nulle pour de grandes valeurs de cette statistique. Pour effectuer ce test,
il nous faut   : (1) r\'ecup\'erer les estimations de
ces param\`etres ainsi que  la matrice des covariances estim\'ee correspondante, ce que nous avons fait au d\'ebut de l'exercice, (2)
calculer la distance du $\chi^2$ de la moyenne estim\'ee \`a la moyenne th\'eorique, et enfin (3) la p-value.

{\small
<<label=test.par.1>>=
# moyenne th\'eorique
theo = as.matrix(c(0.3, 0.4))
# difference
differ = estim-theo
# stat de test
khi = t(differ)%*%solve(vcov.par)%*%differ;
# statistique de test et p-value
cat("stat de test et p-value : ", 
    c(khi, 1-pchisq(khi, df=2)), "\n")
@
}
\noindent

  item En r\'esum\'e, la simplification d'estimation envisag\'ee pour ces deux param\`etres est coh\'erente avec le mod\`ele simul\'e.
Il faut cependant  garder \`a l'esprit que ces estimations de covariances sont bas\'ees sur des approximations grossi\`eres et donc les 
conclusions des tests sont d'abord \`a comprendre comme des indications.

\end{itemize}


\begin{Exercice}[L'Or\'eal]

Dans le traitement de L'Or\'eal, on a encha\^in\'e la mod\'elisation ARMA de la s\'erie et la mod\'elisation GARCH du r\'esidu de cet ajustement.
\begin{itemize}
  \item Effectuer une mod\'elisation GARCH du rendement sans passer par l'\'etape de mod\'elisation par un ARMA.
  \item Superposer sur un m\^eme graphique : la s\'erie \`a pr\'edire et les intervalles de pr\'evision (\`a 90\%) par un ARMA  et  par un GARCH. Commenter.
\end{itemize}

\end{Exercice}

\textbf{R\'eponse.}

\begin{itemize}
  \item Effectuer une mod\'elisation GARCH du rendement sans passer par l'\'etape de mod\'elisation par un ARMA.
  
  
Dans le chapitre, nous avons obtenu le mod\`ele 12.13 en encha\^inant une mod\'elisation ARMA du rendement puis une mod\'elisation GARCH 
du r\'esidu du premier ajustement. Tout d'abord, nous r\'eestimons le mod\`ele ARMA identifi\'e pour L'Or\'eal :

{\small
<<label=lor.1a, echo=TRUE>>=
require(fBasics)
data(csdl)
aa = returns(csdl, percentage = TRUE)
aab = aa[complete.cases(aa) == TRUE,]
# ERRATUM: its package not maintained anymore
# r.csdl = its(aab, as.POSIXct(row.names(aab)))
# r.lor <- rangeIts(r.csdl[,"L_Oreal"], start = "2007-12-28")
# its replaced by zoo package
r.csdl = zoo(aab, as.POSIXct(row.names(aab)))
r.lor <- window(r.csdl[,"L_Oreal"], start = "2007-12-28")
r.lor.0 = r.lor[1:(length(r.lor)-51)]
r.lor.1 = r.lor[(length(r.lor)-50):length(r.lor)]
require("forecast")
mod.r.lor = Arima(as.numeric(r.lor.0), order = c(0, 0, 4), 
                  include.mean = FALSE,
                  fixed = c(NA, NA, 0, NA))
@
}

D'autre part, nous estimons sur le rendement, le mod\`ele GARCH identifi\'e sur les r\'esidus 
de l'ajustement \`a un ARMA :

{\small
<<label=lor.12, echo=TRUE>>=
<<label=esti.var>>=
var.marg.est <- function(mod)
 {param.estim <- mod@fit$par
  std.estim <- mod@fit$se.coef
  k <- which(names(param.estim)=="omega")
  value <- param.estim[k]/(1 - 
        sum(param.estim[(k + 1):length(param.estim)]))
  cat("variance marginale : ", value, "\n")
 }
mod.garch.lor <- garchFit(~garch(1,1), data = as.numeric(r.lor.0), 
                          trace = FALSE, include.mean = TRUE, 
                          na.action = na.pass)
summary(mod.garch.lor)
var.marg.est(mod.garch.lor)
@
}
Nous pouvons voir que le mod\`ele capte bien l'h\'et\'erosc\'edasticit\'e du rendement et que les param\`etres sont significatifs.
Mais le mod\`ele fournit un r\'esidu qui  n'est pas blanc, au contraire de ce qu'avait fourni l'ajustement d'un ARMA/GARCH
section 12.8.1.



  \item Superposer sur un m\^eme graphique : la s\'erie \`a pr\'edire et les intervalles de pr\'evision (\`a 90\%) par un ARMA  et  par un GARCH. Commenter.

Il faut donc pr\'edire le rendement de L'Or\'eal par les deux m\'ethodes.


{\small
<<label=stock.lor>>=
npred = 50
pred.arima.lor = predict(mod.r.lor,n.ahead=npred)
pred.garch.lor = predict(mod.garch.lor,n.ahead=npred)

str(pred.arima.lor)
demi = qnorm(0.95)*pred.arima.lor$se
binf.arima.lor =  pred.arima.lor$pred-demi
bsup.arima.lor =  pred.arima.lor$pred+demi

demi = qnorm(0.95)*pred.garch.lor$standardDeviation
binf.garch.lor = pred.garch.lor$meanForecast-demi
bsup.garch.lor = pred.garch.lor$meanForecast+demi
@
}
 \code{mat.lor} rassemble  les diff\'erentes s\'eries dont on veut superposer les chronogrammes :
l'estimation du rendement moyen, le rendement observ\'e, les bornes inf\'erieures et sup\'erieures obtenues par la mod\'elisation GARCH
puis par la mod\'elisation ARMA.

{\small
<<label=ajust,height=10,fig=TRUE,echo=TRUE,eval=FALSE>>=
par(oma=rep(0.5,4))
mat.lor = cbind(binf.arima.lor, bsup.arima.lor, 
 binf.garch.lor, bsup.garch.lor, as.numeric(r.lor.1[1:npred]))
matplot(1:npred,mat.lor, type = "l", col = "black", 
        lty = c(1, 1, 2, 2, 3), lwd = 2, 
        xlab = "horizon 50", ylab = "rendement")
leg.txt = c("GARCH", "AR", "realisation")
legend(14, 3, leg.txt, lty = c(1, 2, 3))
@
}
\noindent

On note sur le graphique superposant les intervalles de pr\'evision par les deux m\'ethodes (fig. \ref{previreal.Loreal}) que la seule mod\'elisation
par un ARMA fournit des intervalles de pr\'ediction trop amples.


\noindent
\begin{figure}
\begin{center}
<<label=previreal.Loreal, echo=FALSE,results=tex>>=
<<bfig5>>
<<ajust>>
<<zfig2>>
<<bfigps5>>
<<ajust>>
<<zfig2>>
<<zfiginclude>>
@
\end{center}
\caption{L'Or\'eal : intervalles de pr\'evision (ARMA et GARCH) \`a 80\% et r\'ealisation du rendement.}\label{previreal.Loreal}
\end{figure}
%##############################  Fin L'Oreeal #############################################

\end{itemize}


\begin{Exercice}[Danone]

Dans le traitement de Danone, on a encha\^in\'e la mod\'elisation ARMA de la s\'erie et la mod\'elisation GARCH du r\'esidu de cet ajustement.
\begin{itemize}
  \item Effectuer une mod\'elisation GARCH du rendement sans passer par l'\'etape de mod\'elisation par un ARMA.
  \item Superposer sur un m\^eme graphique : la s\'erie \`a pr\'edire et les intervalles de pr\'evision (\`a 90\%) par un ARMA, un ARMA-GARCH et   un GARCH. Commenter.
\end{itemize}
\end{Exercice}

\textbf{R\'eponse.}

\begin{itemize}
  \item Effectuer une mod\'elisation GARCH du rendement sans passer par l'\'etape de mod\'elisation par un ARMA.
Dans le chapitre, nous avons obtenu le mod\`ele 12.16 en encha\^inant une mod\'elisation ARMA du rendement puis une mod\'elisation GARCH 
du r\'esidu du premier ajustement.

{\small
<<label=rdan.0>>=
# ERRATUM
# r.dan = rangeIts(r.csdl[,"Danone"], start= "2007-12-28")
# replaced by :
r.dan <- window(r.csdl[, "Danone"], start = "2007-12-28")
r.dan.0 = r.dan[1:(length(r.dan)-50)]
r.dan.1 = r.dan[(length(r.dan)-49):length(r.dan)]
@
}

{\small
<<label=armagarch.dan,echo=FALSE>>=
mod.arma.dan = Arima(as.numeric(r.dan.0), 
                     order = c(2, 0, 0), 
                     include.mean = FALSE)
mod.garch.dan = garchFit(~garch(1, 1), 
                         data = as.numeric(r.dan.0), 
                         trace = FALSE,
                         include.mean = FALSE, 
                         na.action = na.pass)
@
}
 
  \item Superposer sur un m\^eme graphique : la s\'erie \`a pr\'edire et les intervalles de pr\'evision (\`a 90\%) par un ARMA  et  par un GARCH. Commenter.


Nous calculons les pr\'edictions par les deux mod\`eles.


{\small
<<label=stock.dan>>=
npred = 50
pred.arima.dan = predict(mod.arma.dan, n.ahead = npred)
pred.garch.dan = predict(mod.garch.dan, n.ahead = npred)
str(pred.arima.dan)
demi = qnorm(0.95)*pred.arima.dan$se
binf.arima.dan = pred.arima.dan$pred - demi
bsup.arima.dan = pred.arima.dan$pred + demi
#
demi = qnorm(0.95)*pred.garch.dan$standardDeviation
binf.garch.dan = pred.garch.dan$meanForecast - demi
bsup.garch.dan = pred.garch.dan$meanForecast + demi
@
}
 \code{mat.dan} rassemble  les diff\'erentes s\'eries dont on veut superposer les chronogrammes :
le rendement observ\'e et les bornes inf\'erieures et sup\'erieures obtenues par la mod\'elisation GARCH
puis par la mod\'elisation ARMA.

{\small
<<label=pred.dan,height=10,fig=TRUE,echo=TRUE,eval=FALSE>>=
par(oma=rep(0.5, 4))
mat.dan = cbind(binf.arima.dan, bsup.arima.dan, 
 binf.garch.dan, bsup.garch.dan, as.numeric(r.dan.1[1:npred]))
matplot(1:npred, mat.dan, type = "l", col = "black", 
        lty = c(1, 1, 2, 2, 3), lwd = 2, 
        xlab = "horizon 50", ylab = "rendement")
leg.txt = c("GARCH", "AR", "realisation")
legend(14, 3, leg.txt, lty = c(1, 2, 3))
@
}
\noindent
Nous observons que pour Danone  il y a peu de diff\'erence dans les pr\'evisions GARCH et ARMA (fig. \ref{pred.dan}).

\begin{figure}
\begin{center}
<<label=pred.dan.plot,  echo=FALSE,results=tex>>=
<<bfig5>>
<<pred.dan>>
<<zfig2>>
<<bfigps5>>
<<pred.dan>>
<<zfig2>>
<<zfiginclude>>
@
\caption{Danone : pr\'evision ARMA et GARCH (90\%) et r\'ealisation.}\label{pred.dan}
\end{center}
\end{figure}
\end{itemize}



%% %%%%%%%%%%%%%%%%%%%%%%%%%%% section essilor
%\section{Effets d'une valeur aberrante \label{sec.aber}}
%Nous avons vu au chapitre 2 que le cours d'Essilor pr\'esente un point aberrant et en avons envisag\'e diff\'erentes corrections.
%Ici, nous le corrigeons en le rempla\c{c}ant par la moyenne de ses voisins.
%
%{\small
%<<label=essil.aber>>=
%data(essil)
%i.max = which(essil == max(essil))
%essil[i.max]
%essil1 = essil
%essil1@.Data[i.max] = 0.5*(essil@.Data[i.max-1]+essil@.Data[i.max+1])
%@
%}
%\noindent
%Examinons maintenant l'effet de ce point sur diff\'erentes statistiques en les calculant pour la s\'erie originale et la s\'erie corrig\'ee.
%Nous faisons finir la  premi\`ere p\'eriode au 30 juin 2007.
%
%{\small
%<<label=essil.aber2>>=
%resume = function(cours, titre="")
%{
% avan.06 = rangeIts(cours, end= "2007-06-30")
% apres.06 = rangeIts(cours, start= "2007-07-01")
% rendav.06 = returns(avan.06,percentage= TRUE)
% rendapr.06 = returns(apres.06,percentage= TRUE)
% colnames(rendapr.06) = colnames(cours)
% sk.av06   =  apply(rendav.06,2,skewness, na.rm=TRUE)
% kurt.av06 =  apply(rendav.06,2,kurtosis, na.rm=TRUE)
% sk.apr06  =  apply(rendapr.06,2,skewness, na.rm=TRUE)
% kurt.apr06=  apply(rendapr.06,2,kurtosis, na.rm=TRUE)
% sk06 =  rbind(sk.av06,  sk.apr06, kurt.av06 , kurt.apr06)
% colnames(sk06) =  titre
% rownames(sk06) = c("asym.av","asym.apr","aplat.av","aplat.apr")
% sk06
%}
%non=resume(essil, "s\'erie non corrig\'ee")
%oui=resume(essil1, "s\'erie corrig\'ee")
%@
%}
%
%<<label=essil.aber3,results=tex, echo=FALSE>>=
%res.comp= matrix(NA, ncol=3, nrow=2)
%res.comp[,1]= c(non[c(1,3)])
%res.comp[,2]= c(oui[c(1,3)])
%res.comp[,3]=c(oui[c(2,4)])
%rownames(res.comp)=c("Asym.","Aplat.")
%colnames(res.comp)=c("Non cor. AV","Cor. AV","Apres")
%xtable(res.comp, caption="Asym\'etrie et aplatissement, avant (sans et avec correction)et apr\`es"  , label="avapr", digits=4)
%@
%
%Nous  constatons,  table (\ref{avapr}), que la valeur aberrante modifie tr\`es sensiblement le coefficient d'aplatissement.
%
\end{document}
