######################################################################################################## #objectif: selectionner les variables et construire un modele GAM realisant la meilleur prevision #critere d'evaluation de la prevision: RMSE et MAPE ###########description des donnees ##########donnees de conso. francaise moyenne hebdomadaire de 1996 a 2009 ##########data0: de 1996 a 2008, donnees d'apprentissage ##########data1: 2009, donnees test #nom des variables: #Time: numero d'observation #"Day" "Month" "Year": la date #NumWeek: numero de la semaine dans l'annee /52 (va de 0 a 1 tous les ans) #Load: la consommation en MW #Load1 : la consommation en MW de la semaine precedente (lag 1) #Temp: la temperature moyenne hebdo. en ?C #Temp1: la temperature en ?C de la semaine precedente (lag 1) #IPI: Indice de production industriel (source INSEE) #IPI_CVS: Indice de production industriel corrige des variations saisonnieres ###############Initialisation rm(list=objects()) ###############packages library(mgcv) library(mgcViz) library(gridExtra) library(yarrr) ###############Import des donnees setwd("/Users/yannig/Documents/Enseignement/2017-2018/M2statML/TP/") #C:\Enseignement\2015-2016\Projet Data Mining\TP\Regression data0<-read.table("data_conso_hebdo0.txt", header=TRUE) data1<-read.table("data_conso_hebdo1.txt", header=TRUE) ###############evaluation criteria rmse<-function(eps) { return(round(sqrt(mean(eps^2,na.rm=TRUE)),digits=0)) } mape<-function(y,ychap) { return(round(100*mean(abs(y-ychap)/abs(y)),digits=2)) } #################################################### ##############qq graphiques pour l'analyse descriptive #################################################### attach(data0) plot(Load,type='l') plot(Temp, Load,pch=16,cex=0.5) plot(NumWeek, Load,pch=16,cex=0.5) plot(Load, Load1,pch=16,cex=0.5) acf(Load, lag.max=20) #################################################### ##############exercice 1: choice of k #################################################### plot(Temp, Load,pch=16,cex=0.5) g0 <- gam(Load~s(Temp, k=3, bs="cr"), data=data0) summary(g0) plot(g0) plot(data0$Temp, g0$residuals, pch=16) g1 <- gam(g0$residuals~ s(data0$Temp, k=5, bs="cr")) summary(g1) # # fit0 <- getViz(g0, nsim = 50) # plot(sm(fit0, 1), n = 400) + l_points() + l_fitLine() + l_ciLine() # ( check1D(fit0, "Temp") + l_gridCheck1D(gridFun = mean, stand = "sc", n=100) )$ggObj g1 <- gam(Load~s(Temp, k=10, bs="cr"), data=data0) summary(g1) sqrt(g0$gcv.ubre) sqrt(g1$gcv.ubre) fit1 <- getViz(g1, nsim = 50) summary(g1) plot(sm(fit1, 1), n = 400) + l_points() + l_fitLine() + l_ciLine() ( check1D(fit1, "Temp") + l_gridCheck1D(gridFun = mean, stand = "sc", n=100) )$ggObj ############################################################ ############block Cross Validation for choosing k ############################################################ univ<-function(k, block) { g<- gam(Load~s(Temp, k=k, bs="cr"), data=data0[-block,]) #g<- gam(Load~s(Temp, k=k, bs="cr") + s(Time, k=3) + s(NumWeek, k=20), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } Nblock<-10 borne_block<-seq(1, nrow(data0), length=Nblock+1)%>%floor block_list<-list() l<-length(borne_block) for(i in c(2:(l-1))) { block_list[[i-1]] <- c(borne_block[i-1]:(borne_block[i]-1)) } block_list[[l-1]]<-c(borne_block[l-1]:(borne_block[l])) K<-c(3:20) rmseK<-lapply(K, function(k){lapply(block_list, univ,k=k)%>%unlist%>%rmse} ) plot(K, rmseK, type='b', pch=20) ################################################################################ ############block Cross Validation for choosing k with different spline basis ################################################################################ univ<-function(k, block, bs) { g<- gam(Load~s(Temp, k=k, bs=bs), data=data0[-block,]) #g<- gam(Load~s(Temp, k=k, bs=bs) + s(Time, k=3), data=data0[-block,]) #g<- gam(Load~s(Temp, k=k, bs=bs) + s(Time, k=3) + s(NumWeek, k=20), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } K<-c(4:20) rmseKcr<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='cr')%>%unlist%>%rmse} ) rmseKtp<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='tp')%>%unlist%>%rmse} ) rmseKps<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='ps')%>%unlist%>%rmse} ) rmseKcs<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='cs')%>%unlist%>%rmse} ) col<-piratepal("basel", length.out = 9) plot(K, rmseKcr, type='b', pch=20, ylim= range(rmseKcr, rmseKps), col=col[1]) lines(K, rmseKtp, type='b', pch=20, col=col[2]) lines(K, rmseKps, type='b', pch=20, col=col[3]) lines(K, rmseKcs, type='b', pch=20, col=col[4]) legend("top", col=col, c("cr", "tp", "ps", "cs"), pch=20, ncol=2, bty='n') gcr <- gam(Load~s(Temp, k=5, bs="cr"), data=data0) gcr.fit <- predict(gcr, newdata=data0, type="terms") gtp <- gam(Load~s(Temp, k=5, bs="tp"), data=data0) gtp.fit <- predict(gtp, newdata=data0, type="terms") gps <- gam(Load~s(Temp, k=5, bs="ps"), data=data0) gps.fit <- predict(gps, newdata=data0, type="terms") gad <- gam(Load~s(Temp, k=10, bs="ad", m=5, xt=list(bs="cr")), data=data0) gad.fit <- predict(gad, newdata=data0, type="terms") par(mfrow=c(1,1)) o<-order(data0$Temp) plot(data0$Temp, data0$Load-gcr$coef[1], pch=16, cex=0.5) lines(data0$Temp[o], gcr.fit[o], col=col[1]) lines(data0$Temp[o], gtp.fit[o], col=col[2]) lines(data0$Temp[o], gps.fit[o], col=col[3]) lines(data0$Temp[o], gad.fit[o], col=col[4]) legend("top", col=col[1:4], c("cr", "tp", "ps", "ad"), pch=20, ncol=2, bty='n') ########################################################################################################################## #####plot of each splines scaled by their coef (colors) and the estimated effect (black) for cubic regression basis ########################################################################################################################## gcr.mat <- predict(gcr, newdata=data0, type="lpmatrix") gcr.mat%*%gcr$coefficients-gcr$fitted o<-order(data0$Temp) col<-piratepal("basel", length.out = 5) plot(data0$Temp, data0$Load-gcr$coef[1], pch=16, cex=0.5) for(i in c(2:ncol(gcr.mat))) { lines(data0$Temp[o], gcr.mat[o,i]*gcr$coefficients[i], col=col[i-1]) } lines(data0$Temp[o], gcr.fit[o], col="black", lwd=2) ########################################################################################## #####plot of the different splines basis ########################################################################################## par(mfrow=c(2,2)) col<-piratepal("basel", length.out = 5) listMod<-list(gcr, gtp, gps, gad) names(listMod)<-c('cr', 'tp', 'ps', 'ad') for(i in c(1:length(listMod))) { g.mat <- predict(listMod[[i]], newdata=data0, type="lpmatrix") print(dim(g.mat)) plot(data0$Temp[o],g.mat[o,2] , pch=16, cex=0.25, ylab='',xlab='', main=names(listMod)[i], col=col[1], type='l', ylim=range(g.mat[,-1])) for(j in c(3:ncol(g.mat))) { lines(data0$Temp[o], g.mat[o,j], col=col[j-1]) } } #############en choisissant les noeuds # gcr <- gam(Load~s(Temp, k=3, bs="cr"), data=data0, knots=list(Temp=c(5,10,20))) # g.mat <- predict(gcr, newdata=data0, type="lpmatrix") # # plot(data0$Temp[o],g.mat[o,2] , pch=16, cex=0.25, ylab='',xlab='', main=names(listMod)[i], col=col[1], type='l', ylim=range(g.mat[,-1])) # for(j in c(3:ncol(g.mat))) # { # lines(data0$Temp[o], g.mat[o,j], col=col[j-1]) # } # abline(v=c(5,10,20)) ########################################################################################## #####plot of each splines scaled by their coef (colors) and the estimated effect (black) ########################################################################################## par(mfrow=c(2,2)) col<-piratepal("basel", length.out = 5) listMod<-list(gcr, gtp, gps, gad) names(listMod)<-c('cr', 'tp', 'ps', 'ad') for(i in c(1:length(listMod))) { g.mat <- predict(listMod[[i]], newdata=data0, type="lpmatrix") plot(data0$Temp, data0$Load-listMod[[i]]$coef[1], pch=16, cex=0.25, ylab='',xlab='', main=names(listMod)[i], col='grey') for(j in c(2:ncol(g.mat))) { lines(data0$Temp[o], g.mat[o,j]*listMod[[i]]$coefficients[j], col=col[j-1]) } g.fit <- predict(listMod[[i]], newdata=data0, type="terms") lines(data0$Temp[o], g.fit[o], col="black", lwd=2) } ################################################################################# ############NumWeek effect ################################################################################# univ<-function(k, block, bs) { g<- gam(Load~s(NumWeek, k=k, bs=bs), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } K<-c(4:40) rmseKcr<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='cr')%>%unlist%>%rmse} ) rmseKtp<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='tp')%>%unlist%>%rmse} ) rmseKps<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='ps')%>%unlist%>%rmse} ) rmseKcs<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='cs')%>%unlist%>%rmse} ) par(mfrow=c(1,1)) col<-piratepal("basel", length.out = 9) plot(K, rmseKcr, type='b', pch=20, ylim= range(rmseKcr, rmseKps), col=col[1]) lines(K, rmseKtp, type='b', pch=20, col=col[2]) lines(K, rmseKps, type='b', pch=20, col=col[3]) lines(K, rmseKcs, type='b', pch=20, col=col[4]) legend("top", col=col, c("cr", "tp", "ps", "cs"), pch=20, ncol=2, bty='n') g<- gam(Load~s(NumWeek, k=30, bs='tp'), data=data0) summary(g) Kopt<-10 gcr <- gam(Load~s(NumWeek, k=Kopt, bs="cr"), data=data0) gcr.fit <- predict(gcr, newdata=data0, type="terms") gtp <- gam(Load~s(NumWeek, k=Kopt, bs="tp"), data=data0) gtp.fit <- predict(gtp, newdata=data0, type="terms") gps <- gam(Load~s(NumWeek, k=Kopt, bs="ps"), data=data0) gps.fit <- predict(gps, newdata=data0, type="terms") gad <- gam(Load~s(NumWeek, k=Kopt, bs="cr", m=10, xt=list(bs="cr")), data=data0) gad.fit <- predict(gad, newdata=data0, type="terms") par(mfrow=c(2,2)) o<-order(data0$NumWeek) col<-piratepal("basel", length.out = 5) listMod<-list(gcr, gtp, gps, gad) names(listMod)<-c('cr', 'tp', 'ps', 'ad') for(i in c(1:length(listMod))) { g.mat <- predict(listMod[[i]], newdata=data0, type="lpmatrix") plot(data0$NumWeek, data0$Load-listMod[[i]]$coef[1], pch=16, cex=0.25, ylab='',xlab='', main=names(listMod)[i], col='grey') for(j in c(2:ncol(g.mat))) { lines(data0$NumWeek[o], g.mat[o,j]*listMod[[i]]$coefficients[j], col=col[j-1]) } g.fit <- predict(listMod[[i]], newdata=data0, type="terms") lines(data0$NumWeek[o], g.fit[o], col="black", lwd=2) } ################################################################################# ############Trend effect ################################################################################# univ<-function(k, block, bs) { g<- gam(Load~s(Time, k=k, bs=bs), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } K<-c(3:20) rmseKcr<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='cr')%>%unlist%>%rmse} ) rmseKtp<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='tp')%>%unlist%>%rmse} ) par(mfrow=c(1,1)) col<-piratepal("basel", length.out = 9) plot(K, rmseKcr, type='b', pch=20, ylim= range(rmseKcr, rmseKps), col=col[1]) lines(K, rmseKtp, type='b', pch=20, col=col[2]) legend("top", col=col, c("cr", "tp"), pch=20, ncol=2, bty='n') g<- gam(Load~s(Time, k=4, bs='cr'), data=data0) summary(g) plot(g) ################################################################################# ############3 variables model ################################################################################# gam1<-gam(Load~s(Time, k=3)+s(NumWeek, k=30,bs='tp')+s(Temp,k=5)) summary(gam1) #####temp univ<-function(k, block, bs) { g<- gam(Load~s(Time,k=5)+s(NumWeek,k=30,bs='tp')+s(Temp,k=k, bs=bs), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } K<-c(3:15) rmseK<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='tp')%>%unlist%>%rmse} ) col<-piratepal("basel", length.out = 9) plot(K, rmseK, type='b', pch=20, ylim= range(rmseK), col=col[1]) #####NumWeek univ<-function(k, block, bs) { g<- gam(Load~s(Time,k=5)+s(NumWeek,k=k,bs=bs)+s(Temp,k=5, bs='tp'), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } K<-c(10:40) rmseK<-lapply(K, function(k){lapply(block_list, univ, k=k, bs='tp')%>%unlist%>%rmse} ) col<-piratepal("basel", length.out = 9) plot(K, rmseK, type='b', pch=20, ylim= range(rmseK), col=col[1]) ############################################################################ ##############exercice 2: fit a GAM model consedering all predictors ############################################################################ gam1<-gam(Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5), data=data0) summary(gam1) blockRMSE<-function(equation, block) { g<- gam(as.formula(equation), data=data0[-block,]) forecast<-predict(g, newdata=data0[block,]) return(data0[block,]$Load-forecast) } equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5) Block_residuals<-lapply(block_list, blockRMSE, equation=equation)%>%unlist rmseBloc1<-rmse(Block_residuals) boxplot(Block_residuals) plot(Block_residuals, type='l') plot(data0$Temp, Block_residuals, pch=16) plot(data0$NumWeek, Block_residuals, pch=16) plot(data0$Load1, Block_residuals, pch=16) gam_prov <- gam(Block_residuals~s(Load1), data=data0) summary(gam_prov) fit <- getViz(gam_prov, nsim = 50) plot(sm(fit, 1), n = 400) + l_points() + l_fitLine() + l_ciLine() equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5)+s(Load1, k=10) gam2<-gam(equation, data=data0) summary(gam2) Block_residuals2<-lapply(block_list, blockRMSE, equation=equation)%>%unlist rmseBloc2<-rmse(Block_residuals2) rmseBloc2 plot(data0$IPI, Block_residuals2, pch=16) gam_prov <- gam(Block_residuals2~s(data0$IPI)) summary(gam_prov) plot(data0$IPI_CVS, Block_residuals2, pch=16) gam_prov <- gam(Block_residuals2~s(data0$IPI_CVS)) summary(gam_prov) plot(gam_prov) equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5)+s(Load1, k=10)+s(IPI_CVS) gam3<-gam(equation, data=data0) summary(gam3) Block_residuals3<-lapply(block_list, blockRMSE, equation=equation)%>%unlist rmseBloc3<-rmse(Block_residuals3) rmseBloc3 #####change the IPI_CVS in linear effect equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5)+s(Load1, k=10)+IPI_CVS gam4<-gam(equation, data=data0) summary(gam4) Block_residuals4<-lapply(block_list, blockRMSE, equation=equation)%>%unlist rmseBloc4<-rmse(Block_residuals4) rmseBloc4 plot(data0$Temp1, Block_residuals4, pch=16) gam_prov <- gam(Block_residuals4~s(data0$Temp1)) summary(gam_prov) plot(gam_prov) equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5)+s(Load1, k=10)+IPI_CVS+s(Temp1) gam5<-gam(equation, data=data0) summary(gam5) Block_residuals5<-lapply(block_list, blockRMSE, equation=equation)%>%unlist rmseBloc5<-rmse(Block_residuals5) rmseBloc5 plot(gam5$residuals, type='l') noel = which(abs(data0$Day - 24) <= 3 & data0$Month == 12) consoNoel = vector("numeric", length(data0$Time)) consoNoel[noel] = 1 data0 <- data.frame(data0, consoNoel) plot(data0$Time, gam5$residuals, type='l') select<-which(data0$consoNoel==1) points(data0$Time[select], gam5$residuals[select], col='red', pch=20) equation <- Load~s(Time,k=3)+s(NumWeek,k=30, bs='cc')+s(Temp,k=5)+s(Load1, k=10)+IPI_CVS+s(Temp1) +consoNoel gam6<-gam(equation, data=data0) summary(gam6) Block_residuals6<-lapply(block_list, blockRMSE, equation=equation)%>%unlist rmseBloc6<-rmse(Block_residuals6) rmseBloc6 plot(data0$Time, gam6$residuals, type='l') select<-which(data0$consoNoel==1) points(data0$Time[select], gam6$residuals[select], col='red', pch=20) data0[which(abs(gam6$residuals)>3000), 1:3] noel = which(abs(data1$Day - 24) <= 3 & data1$Month == 12) consoNoel = vector("numeric", length(data1$Time)) consoNoel[noel] = 1 data1 <- data.frame(data1, consoNoel) ychap6 <- predict(gam6, newdata=data1) rmse(data1$Load-ychap) equation <- Load~s(Time,k=3)+s(NumWeek,k=30, bs='cc')+te(Time, Temp, k=c(3, 5))+s(Load1, k=10)+IPI_CVS+s(Temp1) gam7<-gam(equation, data=data0) summary(gam7) ychap <- predict(gam7, newdata=data1) rmse(data1$Load-ychap) par(mfrow=c(1,1)) plot(data1$Load, type='l') lines(ychap, col='red') lines(ychap6, col='blue') ######################################################################################################################## ###############################time weighted models ######################################################################################################################## ####best window windowForecast<-function(equation, Deb, NbYear) { selec0 <- which((data0$Year>=Deb)&(data0$Year<=Deb+NbYear-1)) selec1 <- which(data0$Year==Deb+NbYear) g<-gam(equation, data=data0[selec0,]) g$forecast<-predict(g, newdata=data0[selec1,]) residuals <-data0[selec1,]$Load-g$forecast return(residuals) } equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5)+s(Load1, k=10)+IPI_CVS+s(Temp1)+consoNoel rmseY <- NULL for(NbYear in c(2:9)) { YearList<-c(1996:(2008-NbYear)) residuals <- lapply(YearList, windowForecast, equation=equation, NbYear) rmseY <- c(rmseY, rmse(unlist(residuals))) } plot(c(2:9), rmseY, type='b', pch=20) min(rmseY) ####with exponentiallly decreasing weights NbYear <- 6 YearList <- c(1996:(2008-NbYear)) memory <- 1/100 rmseY <- NULL for(Deb in YearList) { selec0 <- which((data0$Year>=Deb)&(data0$Year<=Deb+NbYear-1)) selec1 <- which(data0$Year==Deb+NbYear) g<-gam(equation, data=data0[selec0,], weights=exp(memory*data0[selec0,]$Time)) g$forecast<-predict(g, newdata=data0[selec1,]) residuals[[Deb-YearList[1]+1]] <-data0[selec1,]$Load-g$forecast rmseY <- c(rmseY, rmse(residuals[[Deb-YearList[1]+1]])) } rmse(unlist(residuals)) ##########grille sur memory memory <- 1/100 weights=exp(memory*data0$Time) plot(data0$Time, weights, type='l') memory <- 1/200 weights=exp(memory*data0$Time) plot(data0$Time, weights, type='l') NbYear <- 6 YearList <- c(1996:(2008-NbYear)) memoryGrid <- 1/(c(1:10)*100) rmseMemo <- array(0, dim=length(memoryGrid)) for(i in c(1:length(memoryGrid))) { rmseY <- NULL memory <- memoryGrid[i] for(Deb in YearList) { selec0 <- which((data0$Year>=Deb)&(data0$Year<=Deb+NbYear-1)) selec1 <- which(data0$Year==Deb+NbYear) g<-gam(equation, data=data0[selec0,], weights=exp(memory*data0[selec0,]$Time)) g$forecast<-predict(g, newdata=data0[selec1,]) residuals[[Deb-YearList[1]+1]] <-data0[selec1,]$Load-g$forecast rmseY <- c(rmseY, rmse(residuals[[Deb-YearList[1]+1]])) } rmseMemo[i] <- rmse(unlist(residuals)) } plot(memoryGrid, rmseMemo, type='l') NbYear <- 6 YearList <- c(1996:(2008-NbYear)) memoryGrid <- 1/(c(1:10)*100) memoryGrid <- seq(1/200, 1/100, length=10) rmseMemo <- array(0, dim=length(memoryGrid)) for(i in c(1:length(memoryGrid))) { rmseY <- NULL memory <- memoryGrid[i] for(Deb in YearList) { selec0 <- which((data0$Year>=Deb)&(data0$Year<=Deb+NbYear-1)) selec1 <- which(data0$Year==Deb+NbYear) g<-gam(equation, data=data0[selec0,], weights=exp(memory*data0[selec0,]$Time)) g$forecast<-predict(g, newdata=data0[selec1,]) residuals[[Deb-YearList[1]+1]] <-data0[selec1,]$Load-g$forecast rmseY <- c(rmseY, rmse(residuals[[Deb-YearList[1]+1]])) } rmseMemo[i] <- rmse(unlist(residuals)) } plot(memoryGrid, rmseMemo, type='l') ######################################################################################################################## ###############################forecasting step ######################################################################################################################## noel = which(abs(data1$Day - 24) <= 3 & data1$Month == 12) consoNoel = vector("numeric", length(data1$Time)) consoNoel[noel] = 1 data1 <- data.frame(data1, consoNoel) gam6.forecast=predict(gam6,newdata=data1) rmse(data1$Load-gam6.forecast) mape(data1$Load,gam6.forecast) equation <- Load~s(Time,k=3)+s(NumWeek,k=30)+s(Temp,k=5)+s(Load1, k=10)+IPI_CVS+s(Temp1)+consoNoel memory <- memoryGrid[which.min(rmseMemo)] gam6.weighted<-gam(equation, data=data0, weights=exp(memory*data0$Time)) gam6.weighted.forecast=predict(gam6.weighted,newdata=data1) rmse(data1$Load-gam6.weighted.forecast) mape(data1$Load,gam6.weighted.forecast) par(mfrow=c(1,1)) plot(data1$Time, data1$Load, type='l') lines(data1$Time, gam6.forecast, col='red') lines(data1$Time, gam6.weighted.forecast, col='blue') ######################################################################################################################## ###############################rolling forecasts ######################################################################################################################## model <- gam6 ychap <- predict(model, newdata=data1[1,]) for(i in c(1:(nrow(data1)-1))) { data.update <- rbind(data0, data1[1:i,]) model <-gam(model$formula, data=data.update) ychap <- c(ychap, predict(model, newdata=data1[i+1,])) } rmse(data1$Load-ychap) mape(data1$Load,ychap) plot(data1$Load, type='l') lines(ychap, col='red') model <- gam6.weighted ychap.weighted <- predict(model, newdata=data1[1,]) for(i in c(1:(nrow(data1)-1))) { data.update <- rbind(data0, data1[1:i,]) model <-gam(model$formula, data=data.update, weights=exp(memory*data.update$Time)) ychap.weighted <- c(ychap.weighted, predict(model, newdata=data1[i+1,])) } rmse(data1$Load-ychap.weighted) mape(data1$Load,ychap.weighted) plot(data1$Load, type='l') lines(ychap.weighted, col='red') lines(ychap, col='blue') plot(cumsum(data1$Load-ychap), col='blue', type='l') lines(cumsum(data1$Load-ychap.weighted), col='red') ########################################################################################## ##########avec bam ########################################################################################## bam6 <- bam(equation, data=data0) model <- bam6 ychap.bam <- predict(model, newdata=data1[1,]) for(i in c(1:(nrow(data1)-1))) { model <- bam.update(model, data=data1[i,]) ychap.bam <- c(ychap.bam, predict(model, newdata=data1[i+1,])) } rmse(data1$Load-ychap.bam) mape(data1$Load,ychap.bam)