#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Laboratorio 3- MODELLI PER CICLO DI VITA DEL PRODOTTO, MODELLO PROPHET - #Aprile 2025 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##librerie --------------------------------------------------------------- library("readxl") library(DIMORA) colors <- grey(seq(0, 1, length = 5)) ##Modelli univariati (BM,GBM, GGM) e modello per competizione (UCRCD) ###################################### ####Transizione energetica in Germania ###################################### bp<- read_excel("BP1.xlsx") GC<- bp$GermanyC[1:55] ##carbone GN<- bp$GermanyN[1:55] ##nucleare GG<- bp$GermanyG[1:55] ##gas GR<- bp$GermanyR[28:55] ##rinnovabili year<- bp$year[1:55] ###Facciamo un grafico complessivo di tutte le serie plot(year,GN,xlab="Year", ylab="Annual consumption", type= "b", pch=16, lty=3, cex=0.7, col=1,ylim=c(0,7)) points(year[28:55],GR, type= "b", pch=2, lty=3, cex=0.7, col=2) points(year,GC, type= "b", pch=5, lty=3, cex=0.7, col=3) points(year,GG, type= "b", pch=6, lty=3, cex=0.7, col=4) ################ ###Nucleare#### ############### ###Modello di Bass bm_GN<-BM(GN,display = T) summary(bm_GN) pred_bmGN<- predict(bm_GN, newx=c(1:55)) pred.instGN<- make.instantaneous(pred_bmGN) ###Grafico delle previsioni con il modello di Bass plot(GN, type= "b",xlab="Year", ylab="Annual consumption", pch=16, lty=3, xaxt="n", cex=0.6, ylim=c(0,2), col=1) axis(1, at=c(1,10,19,28,37,46,55), labels=bp$year[c(1,10,19,28,37,46,55)]) lines(pred.instGN, lwd=2, col=1) #### ####Modello GGM (potenziale dinamico) GGM_GN<- GGM(GN, prelimestimates=c(56.339566617, 0.001, 0.01, 0.001481412, 0.129385437)) summary(GGM_GN) pred_GGMGN<- predict(GGM_GN, newx=c(1:55)) pred.instGGMGN<- make.instantaneous(pred_GGMGN) plot(GN, type= "b",xlab="Year", ylab="Annual consumption", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37,46,55), labels=bp$year[c(1,10,19,28,37,46,55)]) lines(pred.instGGMGN, lwd=2, col=1) ##################### ####Rinnovabili###### ###Modello di Bass bm_GR<-BM(GR,display = T) summary(bm_GR) pred_bmGR<- predict(bm_GR, newx=c(1:55)) pred.instGR<- make.instantaneous(pred_bmGR) ###Grafico delle previsioni con il modello di Bass plot(GR, type= "b",xlab="Year", ylab="Annual consumption", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37), labels=bp$year[28:55][c(1,10,19,28,37)]) lines(pred.instGR, lwd=2, col=1) ###Modello GBM con uno shock esponenziale GBMe1GR<- GBM(GR,shock = "exp",nshock = 1,prelimestimates = c(3.250461e+01, 5.708759e-04, 1.914512e-01, 15,-0.1,0.1)) summary(GBMe1GR) pred_GBMe1GR<- predict(GBMe1GR, newx=c(1:55)) pred.instGBMe1GR<- make.instantaneous(pred_GBMe1GR) ####Grafico di confronto tra BM e GBM plot(GR, type= "b",xlab="Year", ylab="Annual consumption", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37), labels=bp$year[28:55][c(1,10,19,28,37)]) lines(pred.instGBMe1GR, lwd=2, col=2) lines(pred.instGR, lty=2,lwd=2) ####################################################################### ######################################################################## ####Studiamo la competizione tra nucleare e rinnovabili. Unico modello di competizione che viene bene. ###Con il carbone non viene perché manca un pezzo di serie, e il gas ha problemi anche sulla modellazione univariata. ucrcdNR<- UCRCD(GN,GR, display=T) summary(ucrcdNR) plot(year,GN,xlab="Year", ylab="Annual consumption", type= "b", pch=16, lty=3, cex=0.7, col=1) points(year[28:55],GR, type= "b", pch=2, lty=3, cex=0.7, col=2) lines(year[1:55], (ucrcdNR$fitted.i[[1]]), lwd=1, col=3) lines(year[28:55], (ucrcdNR$fitted.i[[2]]), lwd=1, col=4) abline(v=1992, lty=3) legend("topleft", legend=c("Nuclear", "Renewables"), col=c(3,4), lty=1) ############################################################## ###Music data- studio della competizione tra cassette e CD#### ############################################################## music<- read_excel("music.xlsx") ###UCRCD: Competizione tra cassette e CD year<- music$year[1:36] cass<- music$cassette[1:36] cd<- music$cd[10:36] # ##Grafici esplorativi #istantanei plot(year[10:36],cass[10:36],xlab="Year", ylab="Annual sales", type= "b", pch=16, lty=3, cex=0.7, ylim=c(0,1000), col=colors[1]) points(year[10:36],cd, type= "b", pch=2, lty=3, cex=0.7, col=colors[3]) legend("topleft", legend=c("Cassette", "CD"), pch=c(16,2), col=c(colors[1], colors[3])) # #cumulati plot(year[10:36],sum(cass[1:9])+cumsum(cass[10:36]),xlab="Year", ylab="Cumulative sales", type= "b", pch=16, lty=3, cex=0.7, ylim=c(0,15000), col=colors[1]) points(year[10:36],cumsum(cd), type= "b", pch=2, lty=3, cex=0.7, col=colors[3]) legend("topleft", legend=c("Cassette", "CD"), pch=c(16,2), col=c(colors[1], colors[3])) # ####UCRCD con delta e gamma diversi (par = "double") ucrcd<- UCRCD(cass,cd, display=T) summary(ucrcd) # ##Grafico osservati vs previsti (dati istantanei, mi limito alla parte in competizione) plot(year[10:36],cass[10:36],xlab="Year", ylab="Annual sales", type= "b", pch=16, lty=3, cex=0.7, ylim=c(0,1000), col=colors[1]) points(year[10:36],cd, type= "b", pch=2, lty=3, cex=0.7, col=colors[3]) lines(year[10:36], (ucrcd$fitted.i[[1]])[10:36], lwd=1, col=colors[1]) lines(year[10:36], (ucrcd$fitted.i[[2]]), lwd=1, col=colors[3]) legend("topleft", legend=c("Cassette", "CD"), lty=1, lwd=2, col=colors[c(1,3)], pch=c(16,2), cex=0.7) # ############################################################################### ############################################################################### ############################################################################### ##Analisi con modello Prophet- Ciclo di vita di iPhone ############################################################################### #librerie library(prophet) ##dati iphone=read.csv("iphone.csv", sep=";") str(iphone) iphone$quarter=as.Date(iphone$quarter, format="%d/%m/%Y") str(iphone) colnames(iphone)=c("y","ds","cap") #grafico esplorativo plot(iphone$y, type="l", x=iphone$ds, ylab="", xlab="Anno") #trend non lineare #presenza di stagionalità #Modello Prophet # #modello senza stagionalità # m2=prophet(iphone, yearly.seasonality=F, daily.seasonality=F, weekly.seasonality = F, growth="logistic", n.changepoints=15) # #modello con solo stagionalità annuale # m2=prophet(iphone, yearly.seasonality=T, daily.seasonality=F, weekly.seasonality = F, growth="logistic", n.changepoints=15) ####################### #modello con selezione automatica delle stagionalità (il default è "auto" che la mette se ha senso e non la mette se non ha senso) #in questo caso inserisce solo la stagionalità annuale, e dà 2 warning per dire che se vuoi stagionalità giornaliera e settimanale la devi mettere con TRUE ##non specifico i changepoints, che quindi entrano in automatico m2=prophet(iphone, growth="logistic") #n.changepoints di default sono 25 summary(m2) m2$seasonalities future2 <- make_future_dataframe(m2, periods = 20, freq="quarter", include_history = T) tail(future2) future2$cap=100 forecast2 <- predict(m2, future2) tail(forecast2[c("ds", "yhat", "yhat_lower", "yhat_upper")]) #grafico delle previsioni plot(m2, forecast2) dyplot.prophet(m2, forecast2) #grafico con l'aggiunta dei punti di cambio plot(m2, forecast2)+add_changepoints_to_plot(m2, threshold=0) #si possono estrarre le date corrispondenti ai punti di cambio m2$changepoints #grafico della stagionalità plot_forecast_component(m2, forecast2, "yearly") ##modello Prophet specificando che non ci sono i punti di cambio e ponendo stagionalità moltiplicativa (in queso caso appare più ragionevole) m2=prophet(iphone, growth="logistic", n.changepoints=0, seasonality.mode='multiplicative') #summary(m2) #m2$seasonalities future2 <- make_future_dataframe(m2, periods = 20, freq="quarter", include_history = T) tail(future2) future2$cap=100 forecast2 <- predict(m2, future2) tail(forecast2[c("ds", "yhat", "yhat_lower", "yhat_upper")]) #grafico delle previsioni plot(m2, forecast2) dyplot.prophet(m2, forecast2)