%\VignetteIndexEntry{Annexe Chapitre 4}
%\VignetteDepends{}
%\VignetteKeywords{ts}
%\VignettePackage{caschrono}
\documentclass{article}
\usepackage{Sweave}
\usepackage[latin1]{inputenc}
\usepackage{times}
\usepackage{enumerate}
\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{Mod\`eles de base en s\'eries temporelles (compl\'ements du Chapitre 4)}
\author{Yves Aragon\footnote{yves.aragon@gmail.com} \cr
{\normalsize Universit\'e Toulouse Capitole} }
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\setcounter{section}{4}
\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 = "anx4-bitmap-"
@
% 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")
@

<<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="")
@ 

\subsection{Stationnarit\'e}
% ###################################################################
\begin{Exercice} [Danone - test de blancheur] Test de la blancheur du rendement de l'action Danone et celle de son carr\'e. On suivra les \'etapes suivantes :
\begin{enumerate}
\item
Calculer le  carr\'e du rendement centr\'e.
\item
Tester la blancheur du rendement sur toute la s\'erie (on pourra tester la blancheur aux retards 3, 6, 9 et 12). On obtient un r\'esultat inattendu. Lequel ?
\item
Apr\`es examen du chronogramme du rendement (fig. 1.4, chap. 1) on d\'ecide de se limiter \`a l'\'etude de la s\'erie des 600 premi\`eres valeurs. Pourquoi ce choix ?
Qu'a-t-on observ\'e sur le chronogramme ?
\item
Tester la blancheur du rendement et du rendement centr\'e au carr\'e sur la s\'erie de ces 600 valeurs. Conclusion ?
\item
Au vu de ces r\'esultats, le rendement peut-il \^etre un bruit blanc gaussien ?
\end{enumerate}
\end{Exercice} 

\textbf{R\'eponse.} 
\begin{enumerate}[{Q}-1.]
\item Le calcul du rendement s'effectue sans difficult\'e, voir le code ci-dessous.
\item
Si nous prenons toute la s\'erie,
 nous rejetons   la blancheur du
rendement, colonne \code{p-val. s\'erie compl.} du Tableau~\ref{danoblanc}, or il est tr\`es rare qu'un rendement ne soit pas un bruit blanc. 
Nous examinons donc  la s\'erie et constatons,
fig. 1.4, chap. 1,
qu'elle fluctue beaucoup autour de la 700{\ieme}
observation. Aussi nous effectuons maintenant le test sur la sous-s\'erie des 600 premi\`eres
valeurs.

{\small
<<label=lb2>>=
require("caschrono")
require("timeSeries")
data("csdl")
aa <- returns(csdl, percentage = TRUE)
aab <- aa[complete.cases(aa) == TRUE,]
# in previous version we use package its which will not be maintained anymo
# r.csdl = its(aab, as.POSIXct(row.names(aab)))
r.csdl <- zoo(aab, as.POSIXct(row.names(aab)))
r.danone <- r.csdl[, 3]
rdt2 <- (r.danone-mean(r.danone))^2
a0 <- Box.test.2(r.danone, nlag = c(3, 6, 9, 12),
                 type = "Ljung-Box", decim = 4)
a1 <- Box.test.2(r.danone[1:600], nlag = c(3, 6, 9, 12),
                 type = "Ljung-Box", decim = 4)
a2 <- Box.test.2(rdt2[1:600], nlag = c(3, 6, 9, 12),
                 type = "Ljung-Box", decim = 4)
a12 <- cbind(a0, a1[, 2], a2[, 2])
colnames(a12) <- c("Retard", "p-val. serie compl.",
                   "p-val. 600 obs.", "p-val. rdt carre")
@
}

<<label=lb2.tex, results=tex>>=
require(xtable)
xtable(a12, caption = "Table a12 : test de blancheur du rendement 
       de Danone et de son carre.", label = "danoblanc",
       digits = 4)
@


\noindent
\item On constate que le rendement est un bruit blanc, mais que son carr\'e, derni\`ere colonne de la table (\ref{danoblanc}), n'en est pas un.
Notons qu'on a  centr\'e le rendement avant de l'\'elever au carr\'e, au cas o\`u il serait de moyenne non nulle.
\item Un bruit blanc gaussien ($\BBN$) est une suite de v.a. \textit{ind\'ependantes} identiquement distribu\'ees. Si le rendement de Danone \'etait un $\BBN$ alors son 
carr\'e serait une suite de v.a. ind\'ependantes,
identiquement distribu\'ees,
et donc un bruit blanc, ce qu'il n'est pas, donc le rendement de Danone ne peut pas \^etre un $\BBN$.
\end{enumerate}


\subsection{S\'eries lin\'eaires}
Repr\'esentation MA$(\infty)$ des s\'eries :

\begin{subequations} \label{portem.ex}
\begin{eqnarray}
y_t &= -0.7 y_{t-1} + z_t, \;\; \mbox{s\'erie \texttt{y1}}  \notag\\
y_t &= -0.7 y_{t-12} + z_t, \;\; \mbox{s\'erie \texttt{y2}}\notag
\end{eqnarray}
\end{subequations}

Par \code{ARMAtoMA()} :

{\small
<<>>=
ARMAtoMA(-.7, 0, 10)
ARMAtoMA(c(rep(0, 11), -.7), 0, 25)
@
}

Par \code{ImpulseCoefficientsARMA()} de \pkg{FitARMA} :

{\small
<<>>=
require(FitARMA)
ImpulseCoefficientsARMA(-.7, 0, lag.max = 10)
ImpulseCoefficientsARMA(c(rep(0, 11), -.7), 0, lag.max = 25)
@
}
On appelle \'egalement fonction de r\'eponse impulsionnelle la suite des $\psi_i$ dans la repr\'esentation MA$(\infty)$
car $\psi_i$ est la r\'eponse de $y_t$ \`a un choc de "1" sur le bruit en $t-i$.




%\section{Construction d'un  ARMA ou d'un SARMA }
%Dans la d\'efinition d'un ARMA on a pr\'ecis\'e que les polyn\^omes d'autor\'egression et de moyenne mobile ne doivent pas avoir de racine commune.
%Illustrons les cons\'equences de la pr\'esence d'une telle racine. Pour ca nous simulons un ARMA(2,2) avec une racine commune et estimons le mod\`ele sans pr\'ecaution.
%



\subsection{Comparaison ACF et PACF th\'eoriques et empiriques -- exemples}
Dans ces   exemples   on simule des trajectoires de 200 observations   suivant des mod\`eles particuliers.
On trace l'ACF et la PACF th\'eoriques, ainsi que  leurs versions empiriques pour la trajectoire simul\'ee. La comparaison
des versions th\'eoriques et empiriques permet de
se convaincre que l'identification d'un mod\`ele par les ACF et PACF empiriques n'est pas toujours chose ais\'ee.
\\
$ \blacktriangleright$ Mod\`ele MA(2)
Consid\'erons le mod\`ele
\begin{eqnarray}
y_t =  (1 - 0.3  \bm +0.6 \bm^2 )  z_t \hspace{1cm} z_t \sim \BBN(0,1.5). \label{traj.ma2}
\end{eqnarray}
D'une part, \code{ARMAacf()} permet de calculer les ACF et PACF th\'eoriques de cette s\'erie.
D'autre part, nous en simulons une trajectoire  par \code{arima.sim()} enfin
 par \code{acf()} appliqu\'ee \`a la trajectoire simul\'ee, nous calculons
 des versions empiriques de l'ACF et de la PACF.
Nous pouvons donc comparer les versions empiriques et th\'eoriques obtenues.


{\small
<<label=prep_ma2,echo=TRUE>>=
set.seed(951)
ya <- arima.sim(n = 200, list(ma = c(-0.3, 0.6)), sd = sqrt(1.5))
@
}
\noindent
Les diff\'erents r\'esultats sont organis\'es en liste pour la commodit\'e de la repr\'esentation graphique des quatre fonctions.
Les bandes autour de 0  sont d'ordonn\'ees $ \pm 1.96/\sqrt{200}$ (fig. \ref{plot.ma2}).
Les deux ACF sont tr\`es ressemblantes et typiques d'un MA(2). On voit sur la figure que l'ACF empirique prend des valeurs  non significativement diff\'erentes de 0 d\`es le retard 3 et, au  retard 11,
la valeur semble un peu significative.
Les PACF ne sont pas simples \`a interpr\'eter ici, mais leur examen suffit \`a comprendre que le mod\`ele n'est pas un AR pur.
Les graphiques de ces quatre fonctions, sans marges entre eux, sont obtenus  par \code{plotacfthemp()} de \pkg{caschrono}.


<<label=plot.ma2,fig=TRUE,echo=FALSE,eval=FALSE,height=6>>=
<<prep_ma2>>
titre <- "MA(2)"
plotacfthemp(ya, ma = c(-0.3, 0.6), lag.max = 20)
@
\begin{figure}[htbp]
\begin{center}
<<results=tex,echo=FALSE>>=
<<bfig4>>
<<plot.ma2>>
<<zfig2>>
<<bfigps4>>
<<plot.ma2>>
<<zfig2>>
<<zfiginclude>>
@
\end{center}
\caption{Fonctions d'autocorr\'elation d'un MA(2).}
\label{plot.ma2}
\end{figure}
Si l'on veut estimer un MA(2) sur la s\'erie \code{ya}, on \'ecrit la commande

{\small
<<label=ma2.estim>>=
require("forecast")
(mod.ma2 = Arima(ya, order = c(0, 0,2), include.mean = FALSE))
@
}

~\\
$ \blacktriangleright$ Mod\`ele AR(1) :
\begin{eqnarray}
y_t =  -10 + \frac{1}{1+0.8 \bm}  z_t \hspace{1cm} z_t \sim \BBN(0,1.5). \label{acf.ar1.a}
\end{eqnarray}
On en simule une trajectoire de 200 observations et on trace, fig. \ref{plot.acf.ar1}, les ACF et PACF th\'eoriques et empiriques.


{\small
<<label=prep_ar1, echo=TRUE>>=
set.seed(5419)
n2 <- 210
yb <- arima.sim(n = 200, list(ar = -0.8), sd = sqrt(1.5))
yb <- yb - 10
@
}

<<label=plot.ar1,fig=TRUE,echo=FALSE,eval=FALSE>>=
<<prep_ar1>>
plotacfthemp(yb, ar = -0.8, lag.max = 20)
@
\begin{figure}[htbp]
\begin{center}
<<results=tex,echo=FALSE>>=
<<bfig>>
<<plot.ar1>>
<<zfig2>>
<<bfigps>>
<<plot.ar1>>
<<zfig2>>
<<zfiginclude>>
@
\end{center}
\caption{Fonctions d'autocorr\'elation d'un AR(1).}
\label{plot.acf.ar1}
\end{figure}
Les PACF th\'eorique et empirique sont proches. Sur le graphique de l'ACF th\'eorique on v\'erifie qu'aucune autocorr\'elation n'est non nulle, alors que
sur le graphique de l'ACF empirique, si on ne dispose pas simultan\'ement de celui de la PACF empirique, on ne sait pas si les valeurs  dans la bande autour de 0
sont des estimations de 0 ou bien des estimations de quantit\'es faibles mais non nulles. Dans la pratique, il faut donc  envisager les deux possibilit\'es, estimation de 0 ou
estimation de quantit\'es faibles, tant que l'une ou l'autre de ces deux situations n'est pas tr\`es improbable.

L'estimation de (\ref{acf.ar1.a}), par exemple sur \code{yb}, s'obtient par

{\small
<<label=estim.12>>=
(mod12 = Arima(yb, order = c(1, 0, 0)))
@
}
\noindent
On voit que \code{Arima()} note \code{intercept}, ce qui est ici, la moyenne de la s\'erie. La notation est ambigue ; nous y reviendrons au chapitre 5 et dans ses exercices.
 

%

\subsubsection{Fonction d'autocorr\'elation partielle}
% ###################################################################
\begin{Exercice}
Simuler une trajectoire  de 200 valeurs d'un processus autor\'egressif ob\'eissant \`a :
\begin{equation} \label{ar2pacf}
y_t = -0.7 y_{t-1} + 0.2 y_{t-2}+ z_t
\end{equation}
et calculer la PACF empirique jusqu'au retard 4. Comparer avec l'exemple pr\'ec\'edent.
\end{Exercice}








\textbf{R\'eponse.} Il nous faut fixer la variance du bruit. Choisissons $\sigma^2_z = 2$.
D'autre part nous effectuons \code{nsim} simulations pour comparer la pr\'ecision de l'estimateur et nous superposons  valeur th\'eorique
et quantile empirique d'ordres .25 et .75 pour chaque retard :


{\small
<<label=simar2>>=
require(FitARMA)
set.seed(51)
nsim <- 100
nobs <- 200
nsim <- 50
nlag <- 20
y.mat <- matrix(0, nrow = nobs, ncol = nsim)
facp.mat <- matrix(0, nrow = nlag, ncol = nsim)
y.mat <- matrix(0, nrow = 200, ncol = nsim)
facp.mat <- matrix(0, nrow = nlag, ncol = nsim)
for(isim in 1:nsim){
 y.mat[, isim] <- arima.sim(n = nobs, 
                            list(ar = c(-0.7, 0.2)), 
                            sd = sqrt(2))
 facp.mat[, isim] = pacf(y.mat[,isim], 20, plot = FALSE)$acf
}
aa <- t(apply(facp.mat,1, "quantile", probs = c(0.25, .75)))
# pacf theo
theo <- TacvfARMA(phi = c(-.7, .2), lag.max = 20)
pacf.theo <- PacfDL(theo/theo[1], LinearPredictor = TRUE)$Pacf
# intervalle autour de 0 \`a 50%
binf <- qnorm(.25)/nobs^.5
bsup <- qnorm(.75)/nobs^.5
aaa <- cbind(aa, pacf.theo, binf, bsup)
@
}
\code{binf} et \code{bsup} d\'efinissent la bande de confiance asymptotique, valable \'evidemment \`a partir du retard 3, voir la propri\'et\'e 4.6, alors que les quantiles
d'ordres .25 et .75, matrice \code{aa}, donnent les quantiles empiriques d'apr\`es les \Sexpr{nsim} simulations. 
On obtient le graphique  
(fig. \ref{matplot}) pour 50 simulations par

{\small
<<eval=FALSE>>=
matplot(1:20, aaa, type = "l", ylab = "PACF", 
        xlab = "retard", col = "black")
legend("topright", paste("nombre de simulations : ", 
                         as.character(nsim)))
@
}

Il est instructif de faire le graphique pour diff\'erents nombres de simulations ; on se fait ainsi une id\'ee de la vitesse de convergence 
vers la loi asymptotique de la propri\'et\'e 4.6.


<<label=matplot0,fig=TRUE,echo=FALSE,eval=FALSE,height=6>>=
matplot(1:20, aaa, type = "l", ylab = "PACF", 
        xlab = "retard", col = "black")
legend("topright", paste("nombre de simulations : ", 
                         as.character(nsim)))
@

\begin{figure}[htbp]
\begin{center}
<<results=tex,echo=FALSE>>=
<<bfig4>>
<<matplot0>>
<<zfig2>>
<<bfigps4>>
<<matplot0>>
<<zfig2>>
<<zfiginclude>>
@
\end{center}
\caption{Estimations de PACF.}
\label{matplot}
\end{figure}



% ###################################################################
\begin{Exercice} Ajuster un MA$(1)$ \`a $y_t$ simul\'e suivant :

\begin{subequations} \label{portem.ex}
\begin{eqnarray}
y_t &= -0.7 y_{t-1} + z_t, \;\; \mbox{s\'erie \texttt{y1}}  \label{portem.ex1}\\
y_t &= -0.7 y_{t-12} + z_t, \;\; \mbox{s\'erie \texttt{y2}}\label{portem.ex2}
\end{eqnarray}
\end{subequations}
o\`u $z_t$ est un bruit blanc gaussien de variance 4,
 {\ie} un mod\`ele incorrect. Effectuer un test montrant que ce mod\`ele ne convient pas.
\end{Exercice}

\textbf{R\'eponse.}
Nous ajustons  un MA(1) \`a \code{y1} qui suit en r\'ealit\'e un AR(1).

<<label=prep.y1, echo=FALSE>>=
set.seed(123)
y1 <- arima.sim(n = 100, list(ar = -.7), 
                sd = sqrt(4))
y2 <- arima.sim(n = 100, list(ar = c(rep(0, 11), -.7)), 
                sd = sqrt(4))
@
{\small
<<label=estim.ma1>>=
(mod2 <- Arima(y1, order = c(0, 0, 1), include.mean = FALSE))
@
}
\noindent
On n'observe rien de particulier, il n'y  a pas d'avertissement sur une \'eventuelle mauvaise convergence de l'algorithme d'estimation.
Testons maintenant la blancheur du r\'esidu.

{\small
<<label=test.ma1>>=
aa <- Box.test.2(residuals(mod2), nlag = c(3, 6, 9, 12), 
                type = "Ljung-Box",
                decim = 4, fitdf = 1)
colnames(aa) <- c("Retard", "p-val.")
t(aa)
@
}
\noindent
On voit que le r\'esidu de l'ajustement n'est pas un bruit blanc car, au moins  parmi les autocorr\'elations d'ordre 1 \`a 6 certaines ne peuvent \^etre consid\'er\'ees comme nulles. L'ajustement est mauvais et il
faut chercher autre chose. C'est l'examen de l'ACF et de la PACF  de la s\'erie
qui guide vers un bon mod\`ele. Cette d\'emarche, appel\'ee \textit{identification} de la s\'erie,  est d\'etaill\'ee \`a la section 4.6.





% ###################################################################
\begin{Exercice}[Mod\`ele MA$(2)$] On consid\`ere le mod\`ele MA$(2)$
\begin{eqnarray}
y_t =  (1 - 0.3  \bm +0.6 \bm^2 )  z_t \hspace{1cm} z_t \sim \BBN(0,1.5). \label{traj.ma2}
\end{eqnarray}
\begin{enumerate}
  \item Calculer les ACF et PACF th\'eoriques de cette s\'erie par \code{ARMAacf()}.
  \item Simuler une trajectoire de 200 observations  de (\ref{traj.ma2})
par \code{arima.sim()}.
  \item Par \code{acf()} appliqu\'ee \`a la trajectoire simul\'ee, calculer
 des versions empiriques de l'ACF et de la PACF.
\item  Comparer graphiquement les versions th\'eorique et empirique de chaque fonction.
\end{enumerate}
\end{Exercice}

\textbf{R\'eponse.}

{\small
<<>>=
theo <- TacvfARMA(theta = c(.3, -.6), lag.max = 20)
(acf.theo <- theo[-1]/theo[1])
(pacf.theo <- PacfDL(theo/theo[1], LinearPredictor = TRUE)$Pacf)
set.seed(12)
y <- arima.sim(n = 200, list(ma = c(-0.3, .6)), sd = sqrt(1.5))
(acf.emp <- acf(y, 20, plot = FALSE)$acf[-1])
(pacf.emp <- pacf(y, 20, plot = FALSE)$acf)
@
}
\code{acf()} fournit les autocorr\'elations \`a partir du retard 0, tandis que \code{pacf()} fournit les autocorr\'elations partielles \`a partir du retard 1.
Noter \'egalement les changements de convention pour les param\`etres $\phi$ et $\theta$ entre \code{arima()} et les fonctions de \pkg{FitAR} ou \pkg{FitARMA}

% ###################################################################
\begin{Exercice}[Mod\`ele AR$(1)$]
On consid\`ere le mod\`ele AR$(1)$ :
\begin{eqnarray}
y_t =  -10 + \frac{1}{1+0.8 \bm}  z_t \hspace{1cm} z_t \sim \BBN(0,1.5). \label{acf.ar1}
\end{eqnarray}
R\'epondre pour ce mod\`ele, aux questions de l'exercice pr\'ec\'edent.
\end{Exercice}

\textbf{R\'eponse.}

{\small
<<>>=
theo <- TacvfARMA(phi = -.8, lag.max = 20)
(acf.theo <- theo[-1]/theo[1])
(pacf.theo <- PacfDL(theo/theo[1], LinearPredictor = TRUE)$Pacf)
set.seed(23)
y <- arima.sim(n = 200, list(ar = -.8),
               sd = sqrt(1.5)) - 10
(acf.emp <- acf(y, 20, plot = FALSE)$acf[-1])
(pacf.emp <- pacf(y, 20, plot = FALSE)$acf)
@
}



% ###################################################################
\begin{Exercice}
Simuler une trajectoire de 200 observations du mod\`ele  ARMA$(1,2)$ combinant les composantes autor\'egressive et moyenne mobile des mod\`eles pr\'ec\'edents (\ref{acf.ar1} et \ref{traj.ma2} )
et comparer les ACF et PACF th\'eoriques et empiriques.
\end{Exercice}
\textbf{R\'eponse.} (Simuler une trajectoire de 200 ...)
Le mod\`ele \`a simuler est
\begin{eqnarray}
y_t = -10  + \frac{1 - 0.3  \bm +0.6 \bm^2 }{1+0.8  \bm}  z_t, \hspace{1cm} z_t \sim \BBN(0,1.5). \label{arma12.a}
\end{eqnarray}
La simulation ne pose pas de difficult\'es :

{\small
<<label=prep_arma1_21a, echo=TRUE>>=
set.seed(4123)
yc <- arima.sim(n = 200, list(ar = -0.8, ma = c(-0.3, 0.6)), 
                sd = sqrt(1.5)) - 10
@
}
On calcule d'autre part les ACF et PACF th\'eroriques

{\small
<<label=prep_arma1_21a0, echo=TRUE>>=
acf.th <- ARMAacf(ar = -0.8, ma = c(-0.3, 0.6), 
                  lag.max = 20, pacf = FALSE)
pacf.th <- ARMAacf(ar = -0.8, ma = c(-0.3, 0.6), 
                   lag.max = 20, pacf = TRUE)
@
}
\noindent
Enfin on repr\'esente ces quatre fonctions sur un m\^eme graphique, fig. \ref{acf.arma12} :



{\small
<<label=plot.arma1,fig=TRUE,echo=TRUE,eval=FALSE>>=
plotacfthemp(yc, ar = -0.8, ma = c(-0.3, 0.6), 
             lag.max = 20)
@
}
\noindent
\begin{figure}
\begin{center}
<<results=tex,echo=FALSE>>=
<<bfig>>
<<plot.arma1>>
<<zfig2>>
<<bfigps>>
<<plot.arma1>>
<<zfig2>>
<<zfiginclude>>
@
\end{center}
\caption{Fonctions d'autocorr\'elation d'un ARMA(1,2).}
\label{acf.arma12}
\end{figure}
L'estimation de (\ref{arma12.a}), par exemple sur \code{yc}, s'obtient par

{\small
<<label=estim.12>>=
(mod12 = Arima(yc, order=c(1, 0, 2)))
@
}
\noindent
On voit que \code{Arima()} note \code{intercept}, ce qui est ici, la moyenne de la s\'erie.



% ###################################################################
\begin{Exercice}
\code{armasubsets()} de  \pkg{TSA} aide \`a l'identification des mod\`eles ARMA suivant le m\^eme principe que \code{armaselect()} mais reconna\^{\i}t les trous (les plages de coefficients nuls)
dans les ordres de r\'egression.
Utiliser \code{armasubsets()} pour identifier le mod\`ele de \code{yd} simul\'e suivant 
\begin{eqnarray}
y_t = 4  + \frac{1 +0.6 \bm^2 }{(1+0.8  \bm)(1 - 0.7\bm^4)}  z_t, \hspace{1cm} z_t \sim \BBN(0,1.5). \label{ex.sarma}
\end{eqnarray}

\end{Exercice}

\textbf{R\'eponse.}

{\small
<<label=prep_arma1_21b>>=
set.seed(951)
ya <- arima.sim(n = 200, list(ma = c(-0.3, 0.6)), sd = sqrt(1.5))
set.seed(7392)
require("polynom")
autopol <- polynomial(c(1, 0.8)) * polynomial(c(1, 0, 0, 0, -0.7))
yd <- arima.sim(n = 200, list(ar = -autopol[-1], 
                              ma = c(0, 0.6)), sd = sqrt(1.5))
yd <- yd + 4
@
}


{\small
<<label=id.saiso,fig=TRUE,echo=TRUE,eval=FALSE,height=6>>=
require("TSA")
res <- armasubsets(y = yd, nar = 10, nma = 10, y.name = "yd", 
                   ar.method = "ols")
plot(res)
@
}

\begin{figure}[htbp]
\begin{center}
<<results=tex,echo=FALSE>>=
<<bfig>>
<<id.saiso>>
<<zfig2>>
<<bfigps>>
<<id.saiso>>
<<zfig2>>
<<zfiginclude>>
@
\end{center}
\caption{Identification de yd par armasubsets().}
\label{id.saiso}
\end{figure}

Les zones  les plus sombres du graphique (fig. \ref{id.saiso}) sugg\`erent les ordres \`a retenir mais l'interpr\'etation n'est pas simple.

\end{document}



