%\VignetteIndexEntry{Annexe Chapitre 5}
%\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{S\'eries temporelles non stationnaires (compl\'ements du Chapitre 5)}
\author{Yves Aragon\footnote{aragon@cict.fr} \cr
{\normalsize Universit\'e Toulouse 1 Capitole} }
\begin{document}
\maketitle
\setcounter{section}{5}
\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 = "./Figures/anx5-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")
@

<<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{S\'eries int\'egr\'ees - Mod\`eles ARIMA et SARIMA \label{eth1:arisari}}

\begin{exercice}[Diff\'erenciation  saisonni\`ere]
V\'erifier empiriquement l'effet d'une diff\'erenciation saisonni\`ere sur 

\begin{equation}\label{trend.saiso}
y_t=  a + b \cos(2 \pi t/4 )+ c \sin( 2 \pi t/4) + u_t \mbox{~avec~~} u_t  =  \frac{1 }{1 -0.9\: \bm}z_t.
\end{equation}

Pour cela on pourra d\'efinir les s\'eries 
 $\cos(2 \pi t/4 )$ et $ \sin( 2 \pi t/4)$, $t=1,\cdots,48$ et calculer leurs diff\'erences saisonni\`eres.
 
Pour des compl\'ements th\'eoriques on peut consulter  Gourieroux et Monfort (1995),  chap. 3,
qui pr\'esentent les propri\'et\'es alg\'ebriques des 
 filtres
de moyenne mobile, appel\'es aussi filtres de moyenne glissante (\textit{running mean}).  Ladiray et Quenneville (2000) expliquent en d\'etail l'usage de ces filtres en macro-\'economie.

\end{exercice}

\textbf{R\'eponse.}

{\small
<<label=period>>=
# serie de 48 points
# x1 et x2 sont les fonctions periodiques
# et on verifie que leurs differences saisonnieres sont nulles.
temps = 1:48
x1 = cos(2*pi*temps/4)
diff(x1,4)
x2 = sin(2*pi*temps/4)
diff(x2,4)
@
}

On voit que les diff\'erenciations saisonni\`eres de \code{x1} et \code{x2} sont nulles.

\begin{exercice}[Estimation d'un SARIMA avec d\'erive]
On a indiqu\'e \`a la section 5.1, comment \proglang{R} estime les mod\`eles int\'egr\'es. Ainsi, s'il y avait diff\'erenciation aux ordres 1 et 12, il faudrait introduire le
r\'egresseur $t^2$ (la constante et le r\'egresseur $t$ sont \'elimin\'es dans les diff\'erenciations) :
\[
Y_t = c \; t^2 + e_t.
\]
La diff\'erenciation aux ordres 1 et 12 donne
\[
  (1 - \bm)(1 - \bm^{12}) Y_t = 24 \;c    +  (1 - \bm)(1 - \bm^{12}) e_t,
\]
 et \proglang{R} fournit donc une estimation de $c$ et non de $24 \;c$.
 
 V\'erifier cette assertion en simulant un SARIMA$(1,1,0)(0,1,1)_{12}$ puis en l'estimant.
 
 \end{exercice}
 
\textbf{ R\'eponse.}
 Il faut \'evidemment simuler un SARIMA avec d\'erive. Nous allons faire cette simulation de deux fa\c{c}ons.
 D'abord, nous choisissons un mod\`ele, un \\
  SARIMA$(1,1,0)(0,1,1)_{12}$  avec d\'erive :
 \[
 (1 - \bm)(1-\bm^{12} ) y_t  = \texttt{moy} +\frac{1 +0.5\: \bm^{12}}{1+0.8\: \bm} z_t,\;z_t \sim \BBN(0,1)
 \]
\code{moy} est la moyenne du processus diff\'erenci\'e aux ordres 1 et 12 et donne une d\'erive sur le processus $y_t$.
Nous \'ecrivons les diff\'erents polyn\^omes qui seront utilis\'es.

{\small
<<>>=
require(dse)
require(polynom)
require(forecast)
nobs=200
set.seed(234)
ar.1=c(1,.8)
masaiso = polynomial(c(1,rep(0,11),.5))
ar.13=  polynomial(c(1,rep(0,11),-1))*polynomial(ar.1)
ar.2 = polynomial(c(1,-1))*polynomial(ar.1)
moy=5
@
}
\begin{itemize}
  \item \textbf{Simulation par} \code{arima.sim}\\
On commence par  simuler $ (1 - \bm)(1-\bm^{12} ) y_t $ :

{\small
<<label=arimasim>>=
ya = arima.sim(n=nobs,list(order=c(1,0,12),ar=-.8,
ma=c(rep(0,11),.5))) + moy
mean(ya)
Arima(ya,order=c(1,0,0),seasonal=list(order=c(0,0,1),period=12),
include.mean=TRUE)
@
}

Comme v\'erification, nous avons imprim\'e  la moyenne de la s\'erie simul\'ee et estim\'e le mod\`ele attendu.
On int\`egre maintenant la s\'erie en deux \'etapes et \`a chacune, on estime le mod\`ele de la s\'erie int\'egr\'ee.
On proc\`ede de deux fa\c{c}ons, ou bien en utilisant l'option \code{include.drift=TRUE} ou bien en r\'egressant sur la
s\'erie $1, 2,\cdots$

{\small
<<label=arimasim2>>=
# I(1)
yd0=diffinv(ya,1)[-1]
x1 = as.matrix(seq(1,length(yd0)))
Arima(yd0,order=c(1,1,0),
seasonal=list(order=c(0,0,1), period=12),xreg=x1)
Arima(yd0,order=c(1,1,0),
seasonal=list(order=c(0,0,1), period=12),
include.drift=TRUE)
@
}
On obtient les m\^emes r\'esultats par les deux m\'ethodes.
Enfin pour simuler la s\'erie $y_t$, on int\`egre saisonni\`erement \code{yd0}

{\small
<<label=arimasim3>>=
# I(12)
yd2=diffinv(yd0,12)[-(1:12)]
x2 = as.matrix(seq(1,length(yd0))^2)
(m1=Arima(yd2,order=c(1,1,0),seasonal=list(order=c(0,1,1),
         frequency=12),xreg=x2))
(m1b=Arima(yd2,order=c(1,1,0),seasonal=list(order=c(0,1,1),
       frequency=12),include.drift=TRUE))
@
}
\code{m1} estime le coefficient de \code{x2} et la d\'erive 24 fois ce coefficient.
Dans l'estimation utilisant l'option \code{include.drift=TRUE}
on note que la d\'erive n'est pas significative, alors qu'elle vaut 5 dans le mod\`ele simul\'e.

 

   
   \item \textbf{Simulation par} \code{simulate} \\Il nous faut   d\'evelopper le polyn\^ome d'autor\'egression 
   dans
   \[
    (1 - \bm)(1-\bm^{12} ) (1+0.8 \:\bm) y_t= (1+0.8 \:\bm) \texttt{moy} + (1 +0.5\: \bm^{12}) z_t
   \]

On d\'efinit donc les trois polyn\^omes qui d\'efinissent l'autor\'egression
    
{\small
<<label=simulate1>>=    
(ar.1=polynomial(c(1,.8)))
(ar.2 = polynomial(c(1,-1))*ar.1)
# terme autoregressif complet jusqu'au retard 14
(ar.14 =  polynomial(c(1,rep(0,11),-1))*ar.2)
moy=5   
(masaiso = polynomial(c(1,rep(0,11),.5)))
(MA2= array(masaiso,c(13,1,1)))
# 
AR2 = array(ar.14,c(15,1,1))
MA2= array(masaiso,c(13,1,1))
modsarima=   ARMA(A=AR2, B=MA2,C=1)
derive= moy
cte =  derive*predict(polynomial(ar.1),1)
u.t = cte*matrix(rep(1,nobs))
y2 = simulate(modsarima,input=u.t,sampleT=nobs)$output
@
}

Pour v\'erifier, on calcule la moyenne de la s\'erie diff\'erenci\'ee une fois simplement et une fois saisonni\`erement :

{\small
<<>>=
mean(diff(diff(y2,1),12))
@
}
puis estimons le mod\`ele dont on a simul\'e une trajectoire :

{\small
<<>>=
x1 = as.matrix(seq(1,length(y2)))
x2 = as.matrix(seq(1,length(y2))^2)
(m1=Arima(y2,order=c(1,1,0),
seasonal=list(order=c(0,1,1), period=12), 
xreg=x1))
(m2=Arima(y2,order=c(1,1,0),
seasonal=list(order=c(0,1,1),period=12),
include.drift=TRUE))
(m3=Arima(y2,order=c(1,1,0),
seasonal=list(order=c(0,1,1),period=12),
xreg=x2))
@
}
On peut observer que l'estimation n'a pas fonctionn\'e correctement aussi bien en introduisant le r\'egresseur $1,2,\cdots$
<<echo=FALSE>>=
cf3 = m3$coef
@
qu'une d\'erive. On peut enfin v\'erifier que 24 $\times$ \code{cf3["x2"]} = \Sexpr{round(24*cf3["x2"],digits=2)} qui est tr\`es proche de la moyenne empirique de la s\'erie diff\'erenci\'ee deux fois.


\end{itemize}


\begin{exercice}[Lag plot  d'une s\'erie avec d\'erive]
Simuler des s\'eries de 200 points suivant 
\[
y_t = c+ y_{t-4}+u_t,\mbox{~~avec~} u_t = \frac{1}{1 - 0.9 \:\bm }z_t
\]
 avec $\sigma_z^2 = 1$ et 
les valeurs initiales \'egales \`a 0. D'abord avec $c=-.2$ puis $c=.2$.
Dessiner les lag plots de ces s\'eries jusqu'au retard 6 et observer la d\'erive sur ces graphes. 
Commenter.
\end{exercice}

\textbf{R\'eponse.}
Le plus simple est de simuler l'ARMA de la s\'erie diff\'erenci\'ee puis de l'int\'egrer.

{\small
<<label=derive,eval=FALSE>>=
set.seed(51)
c=-.2
y0=arima.sim(list(order=c(1,0,0),ar=.9),n=200)+c
y1=diffinv(y0,4)[-(1:4)]
lag.plot(rev(y1),4,layout=c(2,2),do.lines=FALSE,
   main="Derive negative",diag.col="red")
c=.2
y0=arima.sim(list(order=c(1,0,0),ar=.9),n=200)+c
y1=diffinv(y0,4)[-(1:4)]
lag.plot(rev(y1),4,layout=c(2,2),ask=NA,do.lines=FALSE,
                     main="Derive positive",diag.col="red")
@
}
On observe dans la plupart des simulations qu'au retard 4 il y a plus de points au-dessus de la droite $y=x$
quand la d\'erive est positive et moins de points au-dessus de la droite $y=x$ quand elle est n\'egative.
L'\'ecart vertical \`a la diagonale est de l'ordre de $c$.





\begin{exercice}[R\'egression d'une marche al\'eatoire]

Simuler deux marches al\'eatoires \code{x} et \code{y} ind\'ependantes, par exemple \`a l'aide du code :

{\small
<<label=spur1>>=
set.seed(514) ; nobs=300
y0 = 2+rnorm(nobs)
y = diffinv(y0)
x0 = 2-1.5*rnorm(nobs)
x = diffinv(x0)
@
}

\begin{enumerate}
  \item Dessiner le diagramme de dispersion de (\code{x},\code{y}). 
  \item Superposer les chronogrammes des deux s\'eries. 
  \item Que sugg\`erent ces graphiques ?
  \item Effectuer la r\'egression lin\'eaire de \code{y} sur \code{x}. 
  \item Etudier la stationnarit\'e du r\'esidu  et conclure sur la pertinence de cette r\'egression.
\end{enumerate}

\end{exercice}


\textbf{R\'eponse}

\begin{enumerate}
  \item 
  
  {\small
<<label=spur20,eval=FALSE>>=
plot(x,y,type=l)
@
}

  \item 
 {\small
<<label=spur2,eval=FALSE>>=
plot.ts(cbind(x,y), plot.type = "single",lty=1:2)
legend("topleft",c("x","y"),lty=1:2)
@
}

Le diagramme de dispersion (fig. \ref{spur1c}  haut) sugg\`ere 
contre tout bon sens 
une forte corr\'elation, et les chronogrammes (fig. \ref{spur1c}  bas)  \'evoluent parall\`element si on n\'eglige le
d\'ecrochement entre les deux s\'eries.



<<label=spur1c,fig=TRUE,echo=FALSE,eval=FALSE,height=5>>=
op <- par(mfrow=c(2,1), mar=c(4,4,1,1), oma=c(0,0,0,0))
plot(x,y, xlab='x', ylab='y',pch='+',cex=.6)
plot.ts(cbind(x,y), plot.type='single', ylab='x, y', xlab='temps')
par(op)
@


\begin{figure}[h]
\begin{center}
<<label=spur2.a,  echo=FALSE,results=tex>>=
<<bfig2>>
<<spur1c>>
<<zfig2>>
<<bfigps2>>
<<spur1c>>
<<zfig2>>
<<zfiginclude>>
@
\caption{Diagramme de dispersion (haut) et chronogrammes (bas) de deux marches al\'eatoires, \texttt{x} et \texttt{y}, ind\'ependantes.}\label{spur1c}
\end{center}
\end{figure} 
  \item Les graphiques sugg\`erent une \'evolution parall\`ele des deux s\'eries.
  \item Nous  r\'egressons \code{y} sur \code{x}.

{\small
<<label=spur1>>=
mod4 = lm(y~x)
(aa=summary(mod4))
@
}  

Nous obtenons  \code{R2} $= \Sexpr{round(aa$r.squared,4)}$, valeur  \'elev\'ee, et la r\'egression semble tr\`es significative. Pour aller plus loin,
nous dessinons le chronogramme du r\'esidu et son ACF (fig. \ref{spur2d}). Le chronogramme du r\'esidu montre de longues s\'eries de valeurs de m\^eme signe ; typiquement,
ce r\'esidu n'est pas stationnaire. L'accumulation de valeurs cons\'ecutives de m\^eme signe entra\^{\i}ne une forte valeur de l'autocorr\'elation d'ordre 1. Evidemment, cette autocorr\'elation
est purement num\'erique et n'est pas l'estimation d'une autocorr\'elation th\'eorique, non d\'efinie dans cet exemple.
  \item 
  
{\small
<<label=spur4,eval=FALSE>>=
acf(resid)
@
}    
  
{\small
<<label=spur2d,fig=TRUE,echo=TRUE,eval=FALSE,height=6>>=
op <- par(mfrow=c(2,1), mar=c(4,3,1,0), oma=c(0,0,0,0))
plot.ts(mod4$residuals, xlab='temps',ylab='residu MCO')
abline(h=0)
acf(mod4$residuals,main="")
par(op)
@
}

\begin{figure}[h]
\begin{center}
<<label=spur2.b,  echo=FALSE,results=tex>>=
<<bfig>>
<<spur2d>>
<<zfig2>>
<<bfigps>>
<<spur2d>>
<<zfig2>>
<<zfiginclude>>
@
\caption{R\'esidu de la r\'egression de y sur x - Chronogramme  et ACF.}\label{spur2.b}
\end{center}
\end{figure}

\end{enumerate}






\noindent


R\'esumons. Cet exemple nous a conduit \`a effectuer une r\'egression ayant une significativit\'e illusoire : un R2 \'elev\'e et une 
r\'egression apparemment tr\`es significative, 
le diagramme de dispersion des
 deux s\'eries est trompeur mais
les chronogrammes montrent l'absence de lien entre elles. 
 

Pour compl\'eter le traitement de cet exemple, voyons ce que sugg\`ere un test de racine unit\'e. On sait que le r\'esidu est de moyenne nulle,
d'o\`u, \code{type=none}, indiqu\'e dans l'appel de \code{ur.df()}, ci-dessous.
 Des essais avec diff\'erentes valeurs de \code{lags} nous conduisent \`a \code{lags=0}
%http://en.wikibooks.org/wiki/LaTeX/Floats,_Figures_and_Captions

{\small
\begin{verbatim}
require(urca)
ur.mod4=ur.df(mod4$residuals, type = "none",lags = 0) 
summary(ur.mod4)
\end{verbatim}
}

\noindent
On obtient une p-value assez \'elev\'ee, proche de 10\%. On conserve donc l'hypoth\`ese nulle : la s\'erie du r\'esidu est non stationnaire.

\subsection*{Bibliographie}
Gourieroux C. et Monfort A. (1995). S\'eries temporelles et mod\`eles dynamiques.
Economica, 2 edn.. \newline

Ladiray D. et Quenneville B. (2000). D\'esaisonnaliser avec la m\'ethode x-11. Tech.
Rep.. \url{http://www.census.gov/ts/papers/x11french.pdf}.
\end{document}
