%\VignetteIndexEntry{Annexe Chapitre 10}
%\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{Consommation d'\'electricit\'e (compl\'ements du Chapitre 10)}
\author{Yves Aragon\footnote{yves.aragon@gmail.com} \cr
{\normalsize Universit\'e Toulouse 1 Capitole} }
\begin{document}
\maketitle
\setcounter{section}{10}
\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 = "anx10-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")
@
%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}
Superposer les chronogrammes de   \code{u}, d\'efini section 1 et   \code{u.3c}. Commenter.
\end{exercice}

\textbf{R\'eponse.}

{\small
<<label=init.consel>>=
require(caschrono)
data(khct)
@
}
\noindent
Nous   formons le data frame des variables de la p\'eriode d'apprentissage  :

{\small
<<label=prep.df>>=
khct.df <- as.data.frame(window(cbind(khct, time(khct),
 (time(khct)-1977)^2), end = c(1983, 12)))
colnames(khct.df) <- c("kwh", "htdd", "cldd", "t1", "t1.2")
@
}
\noindent
et  r\'eestimer (\ref{kwh.mco.2}).

{\small
<<label=prep.df2>>=
mod2 = lm(sqrt(kwh) ~ htdd + cldd + t1 + t1.2, data = khct.df)
u = ts(residuals(mod2), start = c(1970, 1), frequency = 12)
@
}



\begin{multline} 
\sqrt{\texttt{kwh}}_t \label{kwh.mco.2}
=
\\
-673.06 + 0.00066 \;\texttt{htdd}_t + 0.01\; \texttt{cldd}_t +
 0.3456\;\texttt{temps}_t -0.0059 \; (\texttt{temps} - 1977)_t^2 + u_t.
\end{multline}

Nous avons besoin du d\'etail de l'estimation conduisant \`a \code{u.3c} :

{\small
<<label=prep.conselec2>>=
kwh1rc <- window(sqrt(khct[,"kwh"]), end = c(1983,12))
xreg1 <- khct.df[ ,c("htdd", "cldd", "t1", "t1.2")]
xreg2 <- xreg1[,-4]
require("forecast")
(mdarx3c <- Arima(kwh1rc, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 1)),
  xreg = xreg2))
u.3c <- kwh1rc - as.matrix(xreg2)%*%as.matrix(mdarx3c$coef[5:7]) - mdarx3c$coef[4]
@
}
Enfin nous superposons les deux r\'esidus :

{\small
<<label=conselec.super,fig=TRUE,eval=FALSE>>=
plot.ts(cbind(u,u.3c), plot.type = "single", lty=1:2)
abline(h=0)
@
}

\begin{figure}[h]
\centering
<<label=conselec0,  echo=FALSE,results=tex>>=
<<bfig5>>
<<conselec.super>>
<<zfig2>>
<<bfigps5>>
<<conselec.plot1>>
<<zfig2>>
<<zfiginclude>>
@
\caption{Chronogrammes superpos\'es des r\'esidus u et u.3c.}\label{conselec.super}
\end{figure}

On observe que les r\'esidus \code{u} et  \code{u.3c} sont  tr\`es proches (fig. \ref{conselec.super}) alors qu'il y a une 
variable explicative de la r\'egression MCO \code{t1.2} absente de la r\'egression MCG avec erreur SARMA$(1,0)(1,1)_{12}$.

\begin{exercice}
Tester que le coefficient de  \verb"cldd" est 10 fois plus grand que celui de \verb"htdd", dans (\ref{conselec_mod1}).

\end{exercice}

\textbf{R\'eponse.}
Le mod\`ele estim\'e :
 \begin{eqnarray}  \label{conselec_mod1}
\begin{split} \label{conselec_mod1}
  \sqrt{\texttt{kwh}}_ t &= -680.9117  + 0.00057 \;\texttt{htdd}_t+  0.0073\;\texttt{cldd}_t+
  \\   & \hspace{1cm}~
  0.3496\; \texttt{temps}_t + u_t, \; t = 1,\cdots,168 \\
  u_t &= \frac{1 - 0.7766 \: \bm^{12}}{(1 - 0.6323\: \bm)(1 -0.984\:\bm^{12})} z_t,
  \hspace{1cm} \widehat{\var}(z_t) = 0.03524
\end{split}
 \end{eqnarray}
  est contenu dans l'objet \code{mdarx3c}. Cet objet contient \'egalement la matrice des covariances estim\'ee des estimateurs des param\`etres.
  Notons $\Sigma$ la sous-matrice des covariances des estimateurs de  $\beta_{cldd}$ et $\beta_{htdd}$. 
  
 Sous l'hypoth\`ese  : $\beta_{cldd} = 10 \;\beta_{htdd}$, 
$\widehat{ \beta}_{cldd} - 10 \;\widehat{\beta}_{htdd} \sim AN(0, \sigma^2)$
o\`u 
\[\sigma^2= \begin{bmatrix}1 & - 10 \end{bmatrix} \Sigma  \begin{bmatrix}1 \\- 10   \end{bmatrix}\]

D'o\`u la statistique de test :
\[
t = \frac{\widehat{ \beta}_{cldd} - 10 \;\widehat{\beta}_{htdd}}{\widehat{\sigma}}
\]
sous $H_0$ elle suit approximativement une loi $\nor(0,1)$. On  remplace  dans $\sigma^2$, tous les param\`etres par leurs estimations.
On rejettera l'hypoth\`ese nulle pour de grandes valeurs absolues de $t$.
Nous sommes maintenant en mesure de calculer $t$.
Par \code{str(mdarx3c)} nous rep\'erons   les estimations des param\`etres, \'el\'ement \verb"mdarx3c$coef",  et leur matrice des covariances estim\'ees,  \'el\'ement 
\verb"mdarx3c$var.coef".

{\small
<<>>=
vcovar = mdarx3c$var.coef[c("cldd", "htdd"), c("cldd", "htdd")]
delta0 = mdarx3c$coef["cldd"] - 10*mdarx3c$coef["htdd"]
sig2 = t(matrix(c(1, -10)))%*%vcovar%*%matrix(c(1, -10))
(t0 = delta0/sig2^.5)
(p.val = 2*(1-pnorm(t0)))
@
}
La p.value est sup\'erieure \`a 18\% ; on peut conserver cette hypoth\`ese.


\end{document}
