%\VignetteIndexEntry{Annexe Chapitre 7}
%\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{Simulation (compl\'ements du Chapitre 7)}
\author{Yves Aragon\footnote{yves.aragon@gmail.com} \cr
{\normalsize Universit\'e Toulouse Capitole} }
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\setcounter{section}{7}
\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/annexe_simul-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{Exercices}

\begin{Exercice}[Simulation d'un SARMA]
On veut simuler une s\'erie ob\'eissant \`a  (1.2).

\begin{itemize}
  \item Tirer d'abord 290 observations i.i.d.   suivant la loi de $z_t$.
  \item Simuler d'apr\`es cette s\'erie, une s\'erie de 240 valeurs ob\'eissant \`a (1.2) 
\end{itemize}

\end{Exercice}

\textbf{R\'eponse. }

Cette simulation est effectu\'ee par 

{\small
<<label=simusar.ini>>=
# graine
set.seed(2761)
innov1 <- rnorm(290, sd = 4.18)
y <- arima.sim(list(order = c(12, 0, 1), 
                    ma = -.7, ar = c(rep(0, 11), .9)), 
               innov = innov1, n.start = 50, n = 240) + 50
y.ts <- ts(y, frequency = 12, start = c(1920, 1))
ytr <- cbind(y.ts, nottem)
colnames(ytr) <- c("serie simulee", "temperature")
@
}

Notons simplement que \code{innov1} correspond \`a $z_t$, \code{y} \`a $y_t - 50$ et qu'on a d\^u expliciter que
les coefficients de la r\'egression de $y_t$ sur ses valeurs pass\'ees jusqu'en $t-11$, sont nuls.



\begin{Exercice}[ARIMA]
On veut simuler une s\'erie de 200 valeurs d'une autor\'egression dont le polyn\^ome a deux racines strictement sup\'erieurs \`a 1 et une racine \'egale \`a 1 :
\[
(1 - \frac{\bm}{1.4}) (1  - \bm)(1 - \frac{\bm}{1.9}).
\]
et la  variance du bruit est \'egale \`a 1. 
\begin{itemize}
   \item Calculer le polyn\^ome d'autor\'egression.
   \item Si on essaie de simuler cette s\'erie directement  \`a l'aide de \code{arima.sim()}, qu'observe-t-on ?
La s\'erie ob\'eit \`a un ARIMA$(2,1,0)$.   Apr\`es avoir consult\'e l'aide en ligne de  cette fonction,
reformuler la simulation pour pouvoir utiliser \code{arima.sim()}.
 
  \item Simuler la s\'erie \`a l'aide de \code{simulate()}.
\end{itemize}
\end{Exercice}

\textbf{R\'eponse.}

\begin{itemize}
   \item Calcul du polyn\^ome d'autor\'egression.

   {\small
<<label= autop>>=
require("polynom")
autop <- polynomial(c(1, -1/1.4)) * polynomial(c(1, -1)) * 
  polynomial(c(1, -1/1.9))
@
}
\noindent

   \item
Simulation directe.
On obtient une erreur car
  la partie autor\'egressive du mod\`ele n'est pas stationnaire. Le processus \`a simuler \'etant un ARIMA(2,1,0), on peut exprimer
le facteur du terme en $(1 - \bm)$ : $(1 - \bm/1.4) (1 - \bm/1.9)$ et simuler l'ARIMA :

<<label= autop2>>=
autop1 <- polynomial(c(1, -1/1.4)) * 
  polynomial(c(1, -1/1.9))
asim8b <- arima.sim(n = 60, list(ar = -autop1[-1], 
                                 order = c(2, 1, 0)))
@



\item  Pour simuler la s\'erie \`a l'aide de \code{simulate()}, on construit le mod\`ele via \code{ARMA()}
 puis on le simule

{\small
<<label=construct3>>=
require(dse)
AR <- array(autop1, c(length(autop1), 1, 1))
MA <- array(1, c(1, 1, 1))
mod2 <- ARMA(A = AR, B = MA)
asim8c <- simulate(mod2, sampleT = 60, sd = 1.5)
@
}
\noindent
ainsi, alors que \code{arima.sim()} ne peut simuler que des ARMA ou des ARIMA explicites, \code{simulate()}, comme \code{filter()} peut  simuler toute autor\'egression.
\end{itemize}



\subsection{Intervention}

\begin{Exercice}
On dispose d'une s\'erie de 100 observations. On sait qu'\`a la date $t1 = 10$, une intervention a provoqu\'e une hausse brutale du
niveau moyen de la s\'erie qui est progressivement revenue \`a son niveau ant\'erieur \`a la date t1. D'autre part, en t2 = 25, une autre intervention
a provoqu\'e une baisse progressive, avec des oscillations, du niveau moyen vers un niveau durablement inf\'erieur.
\begin{enumerate}
  \item Ecrire formellement ce m\'ecanisme. On notera $\omega_i,\: \delta_i, \;i=1,2$ les param\`etres des deux interventions.
  \item Ecrire le code \proglang{R} pour calculer cet effet. Choisir des valeurs sens\'ees pour les param\`etres.
\end{enumerate}
\end{Exercice}


\textbf{R\'eponse.}

L'intervention en $t1$ est associ\'ee \`a une impulsion $P_t^{t1}$ et un amortissement du type
\[
\frac{\omega_1}{1 - \delta_1 \bm},\hspace{1cm} \omega_1,\; \delta_1 > 0
\]
 et celle en $t2$, qui dure, est associ\'ee \`a un \'echelon $S_t^{t2}$ et l'amortissement est du m\^eme type avec maintenant
 $ \omega_2 >0,\; \delta_2 <0$.
 Sans autres pr\'ecisions, l'intervention en $t2$ va provoquer un saut de $\omega_2$. On peut l'att\'enuer en introduisant une intervention ponctuelle
 $P_t^{t2}$
 de coefficient $\omega_3 < 0$.
  
 
A la date $t2$ la s\'erie est n\'ecessairement en train de revenir \`a son niveau moyen initial quand survient l'\'ev\'enement.



\end{document}
