####Lab 2: Nonlinear regression####

###Libraries##################################
library(readxl)
library(DIMORA)
##############################################


#############################################################
#############################################################
####Nonlinear models for new product growth (diffusion models)
#############################################################

###Music data (RIAA)
music<- read_excel("music.xlsx")
str(music)

##create the variable cassette
cassette<- music$cassette[1:36]

###some simple plots
plot(cassette, type="b")
plot(cumsum(cassette), type="b")

###a better plot of the yearly time series
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)])

###we estimate a simple Bass Model 
bm_cass<-BM(cassette,display = T)
summary(bm_cass)

###prediction (out-of-sample)
pred_bmcas<- predict(bm_cass, newx=c(1:50))
pred.instcas<- make.instantaneous(pred_bmcas)

###plot of fitted model 
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)


###we estimate the model with 50% of the data

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)


###we estimate the model with 25% of the data
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)


###Comparison between models (instantaneous)
###instantaneous
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)


###Comparison between models (cumulative)
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)


###exercise: try the same with the CD time series


###Twitter (revenues)

twitter<- read_excel("twitter.xlsx")
length(twitter$twitter)


plot(twitter$twitter, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37,46), labels=twitter$quarter[c(1,10,19,28,37,46)])



###BM
tw<- (twitter$twitter)
Acf(tw)


bm_tw<-BM(tw,display = T)
summary(bm_tw)


pred_bmtw<- predict(bm_tw, newx=c(1:60))
pred.insttw<- make.instantaneous(pred_bmtw)

plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,60), ylim=c(0,40000))
lines(pred_bmtw, lwd=2, col=2)

plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,60))
lines(pred.insttw, lwd=2, col=2)


###GBMr1
GBMr1tw<- GBM(tw,shock = "rett",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 24,38,-0.1))


######GBMe1

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)

plot(cumsum(tw), type= "b",xlab="Quarter", ylab="Cumulative revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,60), ylim=c(0,50000))
lines(pred_GBMe1tw, lwd=2, col=2)


plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,60))
lines(pred_GBMe1tw.inst, lwd=2, col=2)


######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)

plot(tw, type= "b",xlab="Quarter", ylab="Quarterly revenues",  pch=16, lty=3, cex=0.6, xlim=c(1,60))
lines(pred_GGM_tw.inst, lwd=2, col=2)


###Analysis of residuals
res_GGMtw<- residuals(GGM_tw)
acf<- acf(residuals(GGM_tw))

plot(res_GGMtw)





#####################################################################################
###Let us consider the Germany energy transition case. 
#We study consumption of coal, nuclear, and renewables, until year 2019. 

##read the data
bp<- read_excel("BP1.xlsx")
str(bp)
 
###we consider coal, nuclear and renewables
GC<- bp$GermanyC[1:55]
GN<- bp$GermanyN[1:55]
GR<- bp$GermanyR[28:55] 
###eliminate some NA in the data

####Plots

plot(GC, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37,46,55), labels=bp$year[c(1,10,19,28,37,46,55)])


plot(GN, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37,46,55), labels=bp$year[c(1,10,19,28,37,46,55)])


plot(GR, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37), labels=bp$year[28:56][c(1,10,19,28,37)])



#####Coal modelling

##BM
bm_GC<-BM(GC,display = T)
summary(bm_GC)


pred_bmGC<- predict(bm_GC, newx=c(1:60))
pred.instGC<- make.instantaneous(pred_bmGC)


plot(GC, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37,46,55), labels=bp$year[c(1,10,19,28,37,46,55)])
lines(pred.instGC, lwd=2, col=2)



####Nuclear modelling


bm_GN<-BM(GN,display = T)
summary(bm_GN)

pred_bmGN<- predict(bm_GN, newx=c(1:60))
pred.instGN<- make.instantaneous(pred_bmGN)



plot(GN, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6, ylim=c(0,2))
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=2)

### we fit a GGM (better model...but we need to interpret well the parameters)

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:60))
pred.instGGMGN<- make.instantaneous(pred_GGMGN)


plot(GN, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
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=2)



########Renewables modelling

bm_GR<-BM(GR,display = T)
summary(bm_GR)

pred_bmGR<- predict(bm_GR, newx=c(1:60))
pred.instGR<- make.instantaneous(pred_bmGR)


plot(GR, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37), labels=bp$year[28:56][c(1,10,19,28,37)])
lines(pred.instGR, lwd=2, col=1)


###we fit a GBM with one exponential shock

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:60))
pred.instGBMe1GR<- make.instantaneous(pred_GBMe1GR)


plot(GR, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, xaxt="n", cex=0.6)
axis(1, at=c(1,10,19,28,37), labels=bp$year[28:56][c(1,10,19,28,37)])
lines(pred.instGBMe1GR, lwd=2, col=1)


####Comparison between BM and GBMe1

plot(GR, type= "b",xlab="Year", ylab="Annual consumption",  pch=16, lty=3, cex=0.6, xlim=c(1,60), ylim=c(0,3))

lines(pred.instGBMe1GR, lwd=2, col=2)
lines(pred.instGR, lwd=2, col=3)


### We Perform a competition study between nuclear and renewables

year<- bp$year[1:55]

###make a plot of the two series
plot(year,GN,xlab="year", ylab="Consumption (Mtoe)",   type= "b", pch=16, lty=3, cex=0.7)
points(year[28:55],GR, type= "b", pch=16, lty=3, cex=0.7, col=2)

###estimate the UCRCD (with delta and gamma)
ucrcdNR<- UCRCD(GN,GR, display=T)
summary(ucrcdNR)


##we make a plot of the UCRCD model 

plot(year,GN,xlab="year", ylab="Consumption (Mtoe)",   type= "b", pch=16, lty=3, cex=0.7)
points(year[28:55],GR, type= "b", pch=16, lty=3, cex=0.7, col=2)

lines(year[1:55], (ucrcdNR$fitted.i[[1]]), lwd=2, col=1)
lines(year[28:55], (ucrcdNR$fitted.i[[2]]), lwd=2, col=2)

###we also add a line and some text within the plot
abline(v=1991, lty=2)
text(1991, 0.8, pos=4, "renewables enter in 1992")