write("SERIES CHRONOLOGIQUES : TP PRESSE",file="Resul.data")
write(" ",file="Resul.data",append=T)#saut de ligne
write("N.B. Les reponses aux questions concernant les representations graphiques sont a consulter dans le fichier Fig.pdf",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

write("1 MODELE DE BUYS-BALLOT",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write("1.1 Approche descriptive",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write("1.1.1 Donnees et constantes",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write("Chiffre d'affaires de la presse parisienne dans une petite ville de province de 1981 a 1985 (unite=1KF)",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Titre
Titre<-"Chiffre d'affaires de la presse parisienne dans une petite ville de province de 1981 a 1985"

#Presentation des donnees
Donij<-read.table("Presse.txt",header=T,sep="\t",row.names=1,skip=1)
n<-nrow(Donij);p<-ncol(Donij);T=n*p
mois<-row.names(t(Donij));annees<-row.names(Donij)
write("Tableau des donnees",file="Resul.data",append=T)
write.table(format(Donij),file="Resul.data",quote=F,sep="\t",col.names=NA,append=T)
write(" ",file="Resul.data",append=T)
write(c("Parametres : n = ",n,"  p = ",p,"  T = ",T),ncolumns=6,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

pdf("Fig.pdf",width=9, height = 6) #creation du fichier pdf

#Representation directe de la chronique"
Donm<-as.matrix(Donij)
Dont<-as.vector(t(Donm))
Temps<-sequence(T)
rx<-c(0,T);ex<-seq(0,T,p)
ry<-c(75,225);ey<-seq(75,225,25)
par(bg="yellow")
plot(Temps,Dont,type="l",xlim=rx, ylim=ry, xlab="Temps en mois",ylab="Chiffre d'affaires en KF",main=Titre,axes=F,col="red")
axis(1,ex,pos=75);axis(2,ey,pos=0)
title("Representation de la chronique",line=0,font.main=3)

#Representation de la chronique avec la fonction plot.ts
Donts<-ts(Dont,frequency=p,start=c(1981,1))
plot.ts(Donts,xlab="Temps en mois",ylab="Chiffre d'affaires en KF",main=Titre,col="red")
title("Representation de la chronique",line=-1,font.main=3)

#Representation du mouvement saisonnier
xmvt<-matrix(1:(p+2),nr=p+2,nc=n)
ymvt<-cbind(as.matrix(c(NA,Donm[1:n-1,p])),Donm,as.matrix(c(Donm[2:n,1],NA)))
matplot(xmvt,t(ymvt),type="l",axes=F,xlab="mois",ylab="Chiffre d'affaires en KF",main=Titre,xlim=c(0,p+2),ylim=ry)
axis(1,at=0:(p+2),labels=c(NA,mois[12],mois,mois[1]),pos=75);axis(2,ey,pos=0)
text(rep(2,n),Donij[,1],row.names(Donij))
title("Representation du mouvement saisonnier",line=0,font.main=3)

write("1.1.3 Moyennes, tendance et mouvement saisonnier",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
#Moyennes mensuelles
write("Moyennes mensuelles ",file="Resul.data",append=T)
mmens<-as.data.frame(t(mean(Donij)))
write.table(format(round(mmens)),file="Resul.data",quote=F,col.names=TRUE,row.names=F,sep="\t",append=T)
write(" ",file="Resul.data",append=T)

#Moyennes annuelles
write("Moyennes annuelles",file="Resul.data",append=T)
mannu<-as.data.frame(t(mean(as.data.frame(t(Donij)))))
write.table(format(round(mannu)),file="Resul.data",quote=F,col.names=TRUE,row.names=F,sep="\t",append=T)
write(" ",file="Resul.data",append=T)

#Moyenne generale
mge<-mean(Donm)
write(c("Moyenne generale : ",round(mge)),ncolumns=2,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Calcul de la tendance
A<-12/(n*p*(n^2-1))*(c(1:5)%*%t(mannu)-n*(n+1)/2*mge)
B<-mge-A*(n*p+1)/2
write(c("Tendance par calcul direct : pente = ",round(A,2),"  ordonnee a l'origine = ",round(B)),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Calcul direct des coefficients saisonniers
coefs<-mmens-mge-A*(c(1:12)-(p+1)/2)
write("Coefficients saisonniers obtenus directement",file="Resul.data",append=T)
write.table(format(round(coefs)),quote=F,sep="\t",row.names=F,col.names=TRUE,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

write("1.1.4 Tendance et mouvement saisonnier en utilisant la regression lineaire",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write("Difference entre le calcul direct et l'utilisation du modele lineaire pour la tendance et des moyennes mensuelles de la serie initiale privee de tendance pour l'effet saisonnier, unite = 1x10^-15",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
#Utilisation du modele lineaire pour la tendance
mlt<-lm(t(mannu)~seq((p+1)/2,T,p))
write(c("Tendance : pente : ",round(10^15*(A-mlt$coef[2])),"  ordonnee a l'origine : ",round(10^15*(B-mlt$coef[1]))),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Coefficients saisonniers comme moyennes mensuelles de la serie initiale privee de la tendance estimee
write("Coefficients saisonniers :",file="Resul.data",append=T)
coefsm<-mean(as.data.frame(t(matrix(Dont-mlt$coef[2]*Temps-mlt$coef[1],nrow=p,ncol=n))))
write.table(format(round(10^15*(coefs-coefsm))),quote=F,sep="\t",row.names=F,col.names=TRUE,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

write("1.1.5 Serie ajustee et prevision",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Serie ajustee
ajust<-A*Temps+B+rep(as.numeric(coefs),n)
ajustij<-as.data.frame(t(matrix(ajust,nrow=p,ncol=n)))
row.names(ajustij)<-annees
ajustji<-as.data.frame(t(ajustij))
row.names(ajustji)<-mois
ajustij<-as.data.frame(t(ajustji))
write("Serie ajustee",file="Resul.data",append=T)
write.table(format(round(ajustij)),file="Resul.data",quote=F,sep="\t",col.names=NA,append=T)
write(" ",file="Resul.data",append=T)

#Prevision
Tprev<-((T+1):(T+p))
prev<-round(A*Tprev+B+as.numeric(coefs))
previj<-as.data.frame(t(matrix(prev,nrow=p,ncol=1)))
row.names(previj)<-"1986"
prevji<-as.data.frame(t(previj))
row.names(prevji)<-mois
previj<-as.data.frame(t(prevji))
write("Prevision",file="Resul.data",append=T)
write.table(format(previj),file="Resul.data",quote=F,sep="\t",col.names=NA,append=T)
write(" ",file="Resul.data",append=T)

#Representation de la chronique, de la tendance, de la serie ajustee et de la prevision
tendts<-ts(A*Temps+B,frequency=p,start=c(1981,1))
ajustts<-ts(ajust,frequency=p,start=c(1981,1))
seteajts<-cbind(Donts,tendts,ajustts)
ts.plot(seteajts,xlab="Temps en mois",ylab="Chiffre d'affaires en KF",main=Titre,col=c("red","black","blue"))
title("Chronique, tendance et serie ajustee",line=-1,font.main=3)
prevts<-ts(prev,frequency=p,start=c(1986,1))
seteajprts<-cbind(seteajts,prevts)
ts.plot(seteajprts,xlab="Temps en mois",ylab="Chiffre d'affaires en KF",main=Titre,col=c("red","black","blue","green"))
title("Chronique, tendance, serie ajustee et prevision",line=-1,font.main=3)

#Residus
resij<-Donij-ajustij
write("1.1.7 Residus",file="Resul.data",append=T)
write.table(format(round(resij)),file="Resul.data",quote=F,sep="\t",col.names=NA,append=T)
write(" ",file="Resul.data",append=T)

#Representation des residus
rests<-Donts-ajustts
plot.ts(rests,xlab="Temps en mois",ylab="Chiffre d'affaires en KF",main=Titre,col="blue")
title("Residus",line=-1,font.main=3)

write("1.2 Approche inductive",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Estimation de la variance de l'erreur
sigma2<-sum(resij^2)/(T-p-1);sigma<-sqrt(sigma2)
write("1.2.1 Estimation de la variance et de l'ecart-type de l'erreur",file="Resul.data",append=T)
write(c("Variance = ",round(sigma2,2),"  ecart-type = ",round(sigma,1)),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Calcul des residus studentises
resSm<-sqrt(T)*diag(1/(sqrt(p*(n-1)-12*((1:n)-(n+1)/2)^2/(n^2-1))))%*%as.matrix(resij)/sigma
resSij<-as.data.frame(resSm)
row.names(resSij)<-annees
write("1.2.2 Residus studentises",file="Resul.data",append=T)
write.table(format(round(resSij,2)),file="Resul.data",quote=F,sep="\t",col.names=NA,append=T)
write(" ",file="Resul.data",append=T)
write(c("Seuil de l'intervalle de confiance bilateral a 95% pour la loi de Student a ",T-p-1," degres de liberte : ",round(qt(.975,T-p-1),2)),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Representation des residus studentises
resSt<-as.vector(t(resSm))
resSts<-ts(resSt,frequency=p,start=c(1981,1))
plot.ts(resSts,xlab="Temps en mois",ylab="Residus studentises",main=Titre,col="blue")
abline(h=qt(.975,T-p-1),col="red");abline(h=-qt(.975,T-p-1),col="red")
title("Residus studentises",line=-1,font.main=3)

write("1.2.3 Version studentisee des estimations",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)

#Version studentisee des estimations des parametres

#Tendance
AS<-sqrt(T*p^2*(n^2-1)/12)*A/sigma;pAS<-2*(1-pt(abs(AS),T-p-1))
BS<-sqrt(T/(1+3*(T+1)^2/p^2/(n^2-1)))*B/sigma;pBS<-2*(1-pt(abs(BS),T-p-1))
tendS<-as.data.frame(cbind(c("parametre       ","pente             ","ordonnee a l'origine"),c("reel",round(A,2),round(B)),c(" student",round(c(AS,BS),2)),c("p-valeur",format.pval(c(pAS,pBS)))))
write("Estimation de la tendance",file="Resul.data",append=T)
write.table(tendS,file="Resul.data",quote=F,sep="     ",col.names=F,row.names=F,append=T)
write(" ",file="Resul.data",append=T)

#Coefficients saisonniers
write("Estimation des coefficients saisonniers",file="Resul.data",append=T)
coefsS<-sqrt(T)*diag(1/(sqrt(p-1+12*((1:p)-(p+1)/2)^2/p^2/(n^2-1))))%*%t(coefs)/sigma
coefsSd<-as.data.frame(cbind(c("mois",mois),c("reel",format(round(t(as.vector(coefs)),1),nsmall=1)),c("student",format(round(as.vector(coefsS),1),justify="right")),c("p-valeur",format.pval(2*(1-pt(abs(coefsS),T-p-1))))))
write.table(coefsSd,file="Resul.data",quote=F,sep="     ",col.names=F,row.names=F,append=T)
write(" ",file="Resul.data",append=T)

#Test de la necessite du mouvement saisonnier
write("1.2.4 Test de la necessite du mouvement saisonnier",file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write("Tendance avec ou sans mouvement saisonnier",file="Resul.data",append=T)
write(c("Avec : pente : ",format(round(A,3),nsmall=3),"  ordonnee a l'origine : ",round(B,1)),ncolumns=4,file="Resul.data",append=T)
ml0<-lm(Dont~Temps)
write(c("Sans : pente : ",round(ml0$coef[2],3),"  ordonnee a l'origine : ",round(ml0$coef[1],1)),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write("Estimation de la variance de l'erreur avec ou sans mouvement saisonnier",file="Resul.data",append=T)
write(c("Avec : variance = ",round(sigma2,2),"  ecart-type = ",round(sigma,1)),ncolumns=4,file="Resul.data",append=T)
sigma0<-as.numeric(summary(ml0)["sigma"])
write(c("Sans : variance = ",round(sigma0^2,2),"  ecart-type = ",round(sigma0,1)),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)
write(c("Test de Fisher a ",p-1," et ",T-p-1," degres de libertes"),ncolumns=5,file="Resul.data",append=T)
Fobs<-((T-2)* sigma0^2-(T-p-1)*sigma2)/((p-1)*sigma2)
write(c("Statistique observee : ",round(Fobs,2),"    p-valeur : ",format.pval(1-pf(Fobs,p-1,T-p-1))),ncolumns=4,file="Resul.data",append=T)
write(" ",file="Resul.data",append=T)


dev.off(2) #fermeture du fichier pdf