#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Laboratorio 2- MODELLI PER CICLO DI VITA DEL PRODOTTO - Marzo 2025 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##################################################################### ##libraries --------------------------------------------------------------- library("readxl") library(DIMORA) #------------------------------------------------------------------------ ######################### ###Music data (RIAA)#### ######################## music<- read_excel("music.xlsx") str(music) ##creazione della variabile 'cassette' cassette<- music$cassette[1:36] ###alcuni grafici plot(cassette, type="b") plot(cumsum(cassette), type="b") ###un grafico migliore (dati istantanei) plot(cassette, type= "b",xlab="Year", ylab="Annual sales", pch=16, lty=3, xaxt="n", cex=0.6) axis(1, at=c(1,10,19,28,37), labels=music$year[c(1,10,19,28,37)]) ###Bass Model bm_cass<-BM(cassette,display = T) summary(bm_cass) ###previsione (out-of-sample) pred_bmcas<- predict(bm_cass, newx=c(1:50)) pred.instcas<- make.instantaneous(pred_bmcas) ###grafico delle previsioni plot(cassette, type= "b",xlab="Year", ylab="Annual sales", pch=16, lty=3, xaxt="n", cex=0.6) axis(1, at=c(1,10,19,28,37), labels=music$year[c(1,10,19,28,37)]) lines(pred.instcas, lwd=2, col=2) ###stimiamo il modello con il 50% dei dati bm_cass50<-BM(cassette[1:18],display = T) summary(bm_cass50) pred_bmcas50<- predict(bm_cass50, newx=c(1:50)) pred.instcas50<- make.instantaneous(pred_bmcas50) plot(cassette, type= "b",xlab="Year", ylab="Annual sales", pch=16, lty=3, xaxt="n", cex=0.6) axis(1, at=c(1,10,19,28,37), labels=music$year[c(1,10,19,28,37)]) lines(pred.instcas50, lwd=2, col=2) ###stimiamo il modello con il 25% dei dati bm_cass75<-BM(cassette[1:9],display = T) summary(bm_cass75) pred_bmcas75<- predict(bm_cass75, newx=c(1:50)) pred.instcas75<- make.instantaneous(pred_bmcas75) ###Confronto tra modelli ###dati istantanei plot(cassette, type= "b",xlab="Year", ylab="Annual sales", pch=16, lty=3, xaxt="n", cex=0.6) axis(1, at=c(1,10,19,28,37), labels=music$year[c(1,10,19,28,37)]) lines(pred.instcas75, lwd=2, col=2) lines(pred.instcas50, lwd=2, col=3) lines(pred.instcas, lwd=2, col=4) ###Confronto tra modelli ###dati cumulati plot(cumsum(cassette), type= "b",xlab="Year", ylab="Annual sales", pch=16, lty=3, xaxt="n", cex=0.6) axis(1, at=c(1,10,19,28,37), labels=music$year[c(1,10,19,28,37)]) lines(pred_bmcas75, lwd=2, col=2) lines(pred_bmcas50, lwd=2, col=3) lines(pred_bmcas, lwd=2, col=4) ################################# ####Ricavi trimestrali di Twitter ################################## twitter<- read_excel("twitter.xlsx") str(twitter) tw<- (twitter$twitter) ###Grafici preliminari ##ricavi istantanei plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6) ##ricavi cumulati plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues", pch=16, lty=3, cex=0.6) ##studiamo la funzione di autocorrelazione Acf(tw) ####Modelli#### ###Modello di Bass standard bm_tw<-BM(tw,display = T) summary(bm_tw) ###Previsioni con il modello di Bass standard pred_bmtw<- predict(bm_tw, newx=c(1:60)) pred.insttw<- make.instantaneous(pred_bmtw) ###Grafici delle previsioni ##istantaneee plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, col=1, xlim=c(1,60)) lines(pred.insttw, lwd=2, col=1) ##cumulate plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues", pch=16, lty=3, xaxt="n", cex=0.6, col=1, xlim=c(1,60), ylim=c(1,40000)) lines(pred_bmtw, lwd=2, col=1) ###GBM con uno shock rettangolare GBMr1tw<- GBM(tw,shock = "rett",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 24,38,-0.1)) summary(GBMr1tw) ###Previsioni con GBM con uno shock rettangolare pred_GBMr1tw<- predict(GBMr1tw, newx=c(1:60)) pred_GBMr1tw.inst<- make.instantaneous(pred_GBMr1tw) ###Grafici delle previsioni ##istantaneee plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, col=1, xlim=c(1,60)) lines(pred_GBMr1tw.inst, lwd=2, col=1) ##cumulate plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues", pch=16, lty=3, xaxt="n", cex=0.6, col=1, xlim=c(1,60), ylim=c(1,40000)) lines(pred_GBMr1tw, lwd=2, col=1) ######GBM con uno shock esponenziale GBMe1tw<- GBM(tw,shock = "exp",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 12,-0.1,0.1)) summary(GBMe1tw) pred_GBMe1tw<- predict(GBMe1tw, newx=c(1:60)) pred_GBMe1tw.inst<- make.instantaneous(pred_GBMe1tw) ###Grafici delle previsioni ##istantaneee plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, col=1, xlim=c(1,60)) lines(pred_GBMe1tw.inst, lwd=2, col=1) ##cumulate plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues", pch=16, lty=3, xaxt="n", cex=0.6, col=1, xlim=c(1,60), ylim=c(1,40000)) lines(pred_GBMe1tw, lwd=2, col=1) ######GGM GGM_tw<- GGM(tw, prelimestimates=c(4.463368e+04, 0.001, 0.01, 1.923560e-03, 9.142022e-02)) summary(GGM_tw) pred_GGM_tw<- predict(GGM_tw, newx=c(1:60)) pred_GGM_tw.inst<- make.instantaneous(pred_GGM_tw) ###Grafici delle previsioni ##istantaneee plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, col=1, xlim=c(1,60)) lines(pred_GGM_tw.inst, lwd=2, col=1) ##cumulate plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues", pch=16, lty=3, xaxt="n", cex=0.6, col=1, xlim=c(1,60), ylim=c(1,40000)) lines(pred_GGM_tw, lwd=2, col=1) ###Facciamo un grafico complessivo per confrontare i 3 modelli plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, col=1, xlim=c(1,60)) lines(pred_GGM_tw.inst, lwd=2, col=1) lines(pred_GBMe1tw.inst, lty=2, lwd=2) lines(pred_GBMr1tw.inst, lty=3, lwd=2) legend("topleft", legend=c("GBM exponential shock", "GBM rectangular shock", "GGM"), lty=c(2,3,1), lwd=c(2,3,2)) ###Analisi dei residui res_GGMtw<- residuals(GGM_tw) acftw<- acf(residuals(GGM_tw)) fit_GGMtw<- fitted(GGM_tw) fit_GGMtw_inst<- make.instantaneous(fit_GGMtw) ####Affinamento ARMAX library(forecast) s2 <- Arima(cumsum(tw), order = c(2,0,2), xreg = fit_GGMtw) summary(s2) pres2 <- make.instantaneous(fitted(s2)) ###Grafico con affinamento ARMAX plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, col=1, xlim=c(1,60)) lines(pred_GGM_tw.inst, lwd=2, col=1) lines(pres2, lty=1,lwd=1) ################################### ####Ciclo di vita di Apple iPod#### ################################### ##Data data<- read.csv("ipod.csv") ipod<- data$ipod # ###Grafici preliminari plot(ipod, type= "b",xlab="Quarter", ylab="Quarterly sales", pch=16, lty=3, xaxt="n", cex=0.6) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) # # ##BM bm.fit<- BM(ipod,display = T) summary(bm.fit) # pred.bm<- predict(bm.fit, newx=c(1:60)) pred.bmi<- make.instantaneous(pred.bm) # # ###Dati osservati e previsti #istantanei plot(ipod, type= "b",xlab="Quarter", ylab="Quarterly sales", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) lines(pred.bmi, lwd=2, col=2) #cumulati plot(cumsum(ipod), type= "b",xlab="Quarter", ylab="Cumulative sales", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) lines(pred.bm, lwd=2, col=2) # # ###GGM ggm.fit<-GGM(ipod, prelimestimates=c(4.049534e+03, 0.001, 0.1,1.634853e-03,1.493774e-01)) summary(ggm.fit) # ###Previsione con il GGM pred.ggm<- predict(ggm.fit, newx=c(1:60)) pred.ggmi<- make.instantaneous(pred.ggm) # plot(ipod, type= "b",xlab="Quarter", ylab="Quarterly sales", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) lines(pred.ggmi, lwd=2, col=2) # plot(cumsum(ipod), type= "b",xlab="Quarter", ylab="Cumulative sales", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) lines(pred.ggm, lwd=2, col=2) # # ###Affinamento SARMAX # ####Grafico dei residui e autocorrelazione dei residui res.ggm<- residuals(ggm.fit) #Autocorrelazione acf<- acf(residuals(ggm.fit)) # ##Grafico plot(res.ggm, type= "b",xlab="Quarter", ylab="Residuals", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) # ##Calcolo le previsioni in-sample del GGM fitted.ggm<- fitted(ggm.fit) #cumulate fitted.ggmi<- make.instantaneous(fitted.ggm) #istantanee # # ####SARMAX con regressore esterno 'fitted.GGM' s2 <- Arima(cumsum(ipod), order = c(2,0,2), seasonal=list(order=c(1,0,1), period=4),xreg = fitted.ggm) summary(s2) ##previsioni istantanee del modello s2 pres2 <- make.instantaneous(fitted(s2)) # # #####Grafico finale con previsioni secondo il GGM e con affinamento SARMAX plot(ipod, type= "b",xlab="Quarter", ylab="Quarterly sales", pch=16, lty=3, xaxt="n", cex=0.6, col=1) axis(1, at=c(1,10,19,28,37, 46), labels=data$quarter[c(1,10, 19,28,37, 46)]) lines(fitted.ggmi, lwd=1, lty=2) lines(pres2, lty=1,lwd=1) legend("topleft", legend=c("GGM", "GGM+SARMAX "), lty=c(2,1), lwd=1) # ##ACF dei residui del modello s2 resAR<- residuals(s2) AR.acf<- Acf(resAR) # ###################################### ####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)