Tag Archives: M2

Devoir, Méthodes de Prévisions (consommation électrique)

Pour le projet de l’UE de prévision, je demande une prévision de la consommation d’électricité. Pour construire les modèles, je donne accès à un gros fichier qui donne, pour toute la France, la consommation d’électricité, entre 2012 et 2014, heure par heure.

  • comme en cours, la première étape sera de proposer une série de 52 prévisions pour chaque semaine de 2015,
  • la seconde étape sera de proposer deux séries de 31 prévisions, pour chaque jour d’août 2015 et de décembre 2015
  • enfin, la dernière étape sera une prévision de 168 prévisions, pour chaque heure de la semaine du 19 au 25 octobre 2015

Pour aider (?), je mets un lien vers la température en France  dans quelques villes, paris, toulouse, grenoble, nantes, dunkerque (via ftp.ncdc.noaa.gov). Il sera possible d’utiliser la vraie température de 2015 pour faire les prévisions de consommation électrique.  Idéalement, je veux un fichier html (markdown) et 4 fichiers csv avec les 4 séries de prévision.

Moyennes de Prédicteurs

La dernière fois, nous avions vu 4 modèles de séries temporelles, pour modéliser la série hebdomadaire de températures. Le premier était obtenu avec auto.arima, les deux étaient ensuite ses modèles SARIMA (le second était une forte autocorélation à l’ordre 52, et le troisième avait une racine unité saisonnière), et le dernier était obtenu par la méthode de Buys-Ballot (cyclique, avec un bruit autocorrélé).

Au lieu de se demander lequel des quatre est le “meilleur”, on peut tenter de trouver une moyenne pondérée, \hat{y}_t = \sum_{i}\omega_i \hat{y}_t^{(i)}\omega=(\omega_i) est un point du simplexe, au sens où \omega_i\geq 0 et \omega_1+...+\omega_k=1. On veut ensuite trouver les “meilleurs” pondérations, par exemple au sens de la norme L^2. On va calibrer nos quatre modèles sur les m premières valeurs, et on compare ensuite la combinaison des prévions sur les 111 prochaines valeurs (et les vraies valeurs), \omega^\star=\text{argmin}\{\sum_{h=1}^{111}\left(y_{t+h}- \sum_{i}\omega_i \hat{y}_{t+h}^{(i)}\right)^2\}

Par exemple, on peut utiliser les 200 premières valeurs.

m=200
Ym=ts(y[1:m],start=1996,frequency=52)

On va ensuite estimer les 4 modèles sur ces 200 valeurs

model1=arima(Ym,order = c(4,0,0), 
seasonal = list(order = c(1,0,0))) 
model2=arima(Ym,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
model3=arima(Ym,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
regL=lm(Temp~0+as.factor(NumWeek),
data=electricite[1:m,])
E=residuals(regL)
model4=arima(ts(E,start=1996,
frequency=52),order = c(1,0,0))

On peut ensuite collecter les prédictions

DONNEES = data.frame(
  y1=predict(model1,n.ahead = 111)$pred,
  y2=predict(model2,n.ahead = 111)$pred,
  y3=predict(model3,n.ahead = 111)$pred,
  y4=predict(model4,n.ahead = 111)$pred+
     predict(regL,newdata=electricite[m+(1:111),]),
  vraie=electricite[m+(1:111),"Temp"])

Pour créer une combinaison linéaire des prévisions, on utilise

a=rep(1/4,4)
y_prev = as.matrix(DONNEES[,1:4]) %*% a

On peut ainsi visualiser les 4 prévisions, leur combinaison linéaire (ici avec des équi-pondérations), et leurs observations

plot(1:111,DONNEES[,5],type="l",lwd=2)
for(i in 1:4) lines(1:111,DONNEES[,i],col="blue")
lines(1:111,yp,col="red")

Pour trouver la valeur “optimale” des pondérations, qui minimisera la somme des carrés des erreurs, on utilise le code suivant (on fera une transformation des données, ici sur le simplexe, pour s’assurer que l’optimisation donnera un élément du simplexe)

erreur=function(a) sum((DONNEES[,1:4] 
%*% a-DONNEES[,5])^2)
erreur2=function(b){
  a=rep(NA,length(b)+1)
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  erreur(a)
}

On obtient l’optimum, et on peut obtenir une prévision

b=optim(par=c(0,0,0),erreur2)$par
a=rep(NA,length(b)+1)
a[1:length(b)]=exp(b)/(1+sum(exp(b)))
a[length(b)+1]=1/(1+sum(exp(b)))
yp=DONNEES[,1:4] %*% a
lines(1:111,yp,col="red",lwd=4)

On peut alors faire changer le nombre d’observations utilisées

Techniquement, ca pose un petit soucis, car il faut s’assurer la convergence de deux algorithmes : celui d’estimation des paramètres des SARIMA, celui de recherche du paramètre de poids.

poids_optimaux=function(m){
  Ym=ts(y[1:m],start=1996,frequency=52)
  TRY=try(model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0))),FALSE)
  if(inherits(TRY, "try-error")) model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0)),method="CSS")
  TRY=try(model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0))),FALSE)
  if(inherits(TRY, "try-error")) model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0)),method="CSS")
  TRY=try(model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1))),FALSE)
  if(inherits(TRY, "try-error")) model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1)),method="CSS")
regL=lm(Temp~0+as.factor(NumWeek),data=electricite[1:m,])
  E=residuals(regL)
model4=arima(ts(E,start=1996,frequency=52),order = c(1,0,0))
  DONNEES = data.frame(
    y1=predict(model1,n.ahead = 111)$pred,
    y2=predict(model2,n.ahead = 111)$pred,
    y3=predict(model3,n.ahead = 111)$pred,
    y4=predict(model4,n.ahead = 111)$pred+
      predict(regL,newdata=electricite[m+(1:111),]),
    vraie=electricite[m+(1:111),"Temp"])
  erreur=function(a) sum((as.matrix(DONNEES[,1:4]) %*% a-DONNEES[,5])^2)
  erreur2=function(b){
    a=rep(NA,length(b)+1)
    a[1:length(b)]=exp(b)/(1+sum(exp(b)))
    a[length(b)+1]=1/(1+sum(exp(b)))
    erreur(a)
  }
  TRY=try(b <- optim(par=c(0,0,0),erreur2)$par, FALSE)
  a=rep(NA,length(b)+1)
  if(!inherits(TRY, "try-error")) {
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  }
  return(a)
}

On peut alors regarder l’évolution des pondérations au cours du temps. On peut vectoriser, mais si l’algorithme aboutit sur une erreur, c’est plus simple d’identifier

VM=seq(200,480,by=10)
VA=matrix(NA,length(VM),4)
for(i in 1:length(VM)) VA[i,]=poids_optimaux(VM[i])
colnames(VA)=c("auto","SARMA","SARIMA","Buys Ballot")

(pour débugger, c’est plus simple les boucles). On peut alors obtenir le graphique suivant, avec en rose Buys-Ballot, et en vert le modèle SARIMA, avec une racine unité saisonnière,

tva=t(VA)
barplot(tva,legend = rownames(counts),
 col=c("yellow","orange","green","pink"))

On se rend compte qu’assez vite, le modèle le plus intéressant est le modèle de Buys Ballot.

Il est possible de changer de fonction de perte, par exemple si on s’intéresse à un quantile à 90%,

tau=.9
loss=function(e) e*(tau-(e<=0)*1)

Dans la fonction, on utilise alors

erreur=function(a) sum(loss(as.matrix(
DONNEES[,1:4]) %*% a-DONNEES[,5]))

C’est, cette fois, beaucoup plus incertain, même si les deux modèles les plus intéressants sont toujours le SARIMA et Buys-Ballot.

Consommation et production électrique

Une version temporaire des slides est en ligne. Pour le TP, on continuera à utiliser les données suivante

url="http://freakonometrics.free.fr/elec.RData"
download.file(url,"elec.RData")
load("elec.RData")
k=620
n=nrow(electricite)
passe=1:k
futur=(k+1):n
y=electricite$Load[passe]
plot(y,type="l")

La semaine passée, nous avions commencé de modéliser la série de température (cette dernière ayant un impact important comme déterminant de la consommation électrique)

y=electricite$Temp[passe]
plot(y,type="l")
Y=ts(y,start=1996,frequency=52)
library(forecast)

plot(y[(53:k)-52],y[53:k])
abline(lm(y[53:k]~y[(53:k)-52]),col="red")
abline(a=0,b=1,lty=2)

La série temporelle est relativement autocorrélée, à l’ordre 52

acf(y,lag=120)

model1=auto.arima(Y)
acf(residuals(model1),120)

On va garder ce modèle en mémoire, et regarder par la suite ses prévisions. Tentons ici un SARIMA, à la main

model2=arima(Y,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
summary(model2)
model2=arima(Y,order = c(0,0,0), 
seasonal = list(order = c(2,0,0)))
summary(model2)
acf(residuals(model2),120)
model2=arima(Y,order = c(2,0,0), 
seasonal = list(order = c(1,0,0)))            
model2=arima(Y,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
summary(model2)
acf(residuals(model2),120)

Tentons ensuite avec une racine unité saisonnière

Z=diff(Y,52)
model3=arima(Z,order = c(0,0,1), 
seasonal = list(order = c(0,0,1)))
summary(model3)
model3=arima(Z,order = c(0,0,2), 
seasonal = list(order = c(0,0,1)))
summary(model3)
acf(residuals(model3),120)
model3=arima(Y,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
summary(model3)

On peut ensuite tenter un modèle à la Buys-Ballot

reg=lm(Temp~0+as.factor(NumWeek),
data=electricite[passe,])
cycle=predict(reg,newdata=electricite[passe,])
plot(y,type="l")
lines(cycle,col="red")
E=residuals(reg)
plot(E,type="l")
acf(E,120)
model4=arima(E,order = c(1,0,0))
summary(model4)

horizon=n-k
prev4E=predict(model4,n.ahead = horizon)$pred
prev4C=predict(reg,newdata=electricite[futur,])
prev4 =prev4E+prev4C
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(prev4,col="blue")

prev1=predict(model1,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev1,col="purple")

prev2=predict(model2,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev2,col="green")

prev3=predict(model3,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev3,col="orange")

MODELES = data.frame(y1=prev1,
y2=prev2,y3=prev3,y4=prev4, 
y =electricite$Temp[futur])
plot(MODELES$y,type="l")
lines(1:111,MODELES$y1,col="purple")
lines(1:111,MODELES$y2,col="green")
lines(1:111,MODELES$y3,col="orange")
lines(1:111,MODELES$y4,col="blue")

 

Material for the Course on Inequality

Here is the video that we discussed during the first course,

References for the course are the following,

For similar course lectures notes, see Emmanuel Flachaire’s ECON-473 webpage, as well as Michel Lubrano’s notes.

For the project, some references in French,

and in English

There is also the following report (one can focus only on a chapter of the report)

 

Exposés de Méthodologie Statistique (M2)

Tout au long de l’année, je vais animer un cours d’exposés de méthodologie statistique, pour les étudiants de dernière année de Master. L’objectif de ce travail est d’apprendre à présenter des outils ou techniques statistiques avancées. Je reviendrais jeudi sur l’idée du cours, et sur ce qui est demandé (cf doc).

Le travail sera réalisé par groupe de 2 étudiants. Des créneaux seront réservés sur l’emploi du temps pour le suivi des projets. Des soutenances seront organisées  auxquelles tous les étudiants devront assister.

D’un point de vue plus pratique, il s’agit de préparer un exposé de 45 minutes présentant une méthode statistique choisie dans la liste ci-dessous. L’exposé s’adressera à un public de scientifiques et/ou d’économistes ayant déjà un bagage mathématique et statistique et devra leur permettre de comprendre la méthode, son utilisation et son application. La présentation de composera

  • d’une introduction permettent notamment de se faire une idée des applications possibles de la méthode présentée ; on situera aussi la méthode par rapport à d’autres méthodes existantes.
  • d’une partie théorique.
  • d’un exemple illustratif (les étudiants choisiront eux-mêmes le problème, avec des données associées et l’exemple sera traité avec un logiciel statistique). L’application doit bien montrer en quoi les résultats obtenus apportent des informations importantes et utilisables dans le cadre d’une étude.

De plus, il sera important de dominer la portée interprétative des résultats statistiques en relation avec la réalité étudiée et en rapport avec les raisons pour lesquelles cette réalité est étudiée comme on le fait dans un contexte professionnel.

  1. Régressions Ridge et Lasso
  2. Régression PLS
  3. Arbres de décision et méthode CART
  4. Bagging et Boosting
  5. Réseaux de neurones
  6. Perceptron et SVM
  7. Processus AR à seuil
  8. Processus AR fractionnaire
  9. Modèles de chaînes de Markov cachées
  10. Quantiles et expectiles
  11. Algorithme EM
  12. Estimation de mélanges de lois
  13. Théorie des Valeurs Extrêmes
  14. Interpolation spatiale et Krigeage
  15. ICA, Analyse en Composantes Indépendantes
  16. Fonction de profondeur et quantiles multivariés