load("/Users/utente/Documents/prova r pkg/lucidi del corso/RMatlab.RData")
tsgtm <- ts(datigtm[,2:3], frequency = 12, start = c(2004, 1)) # frequency 12 => Monthly Data
par(mar=c(4.5,4.5,1,1))
plot(tsgtm,cex.lab=2,main="",plot.type = "single",lty=1,lwd=c(2,1),ylab="")
legend(2012,90, c("R","Matlab"), lwd=c(2,1),cex=1.5)
library(lubridate)
tsgtw <- ts(datigtw[,2:3], freq=365.25/7, start=decimal_date(ymd("2014-06-29")))
par(mar=c(4.5,4.5,1,1))
plot(tsgtw,cex.lab=2,main="",plot.type = "single",lty=1,lwd=c(2,1),ylab="")
legend(2015,50, c("R","Matlab"), lwd=c(2,1),cex=1.5)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/percapitagdp.RData")
tsgdp <- ts(gdppc[6:189,2:183], frequency = 1, start = c(1820))
plot(tsgdp[,c(5,7)],cex.lab=2,main="",plot.type = "single",lty=c(1,1),lwd=c(2,1),ylab="")
legend(1850,15000, c("Francia","Italia"), lty=c(1,1),lwd=c(2,1),cex=1.5)
tsgdp_it_bi <- ts(gdppc_it_bi[,2:4], frequency = 1, start = c(1861))
plot(tsgdp_it_bi[,3],main="",lwd=2,ylab="")

data("AirPassengers")
AP <- AirPassengers
plot(AP, ylab="Passengers (1000s)", type="o", pch =19, cex=0.7)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/prezzistat.RData")
prezzi_istat <- ts(prezzi[,2:17], frequency = 1, start = c(1861))
par(mar=c(4.5,1,1,1))
plot(prezzi_istat[,c(1,5,16)],cex.lab=2,main="",plot.type = "single",lty=c(1,2,2),lwd=c(2,2,1),ylab="")
legend(1865,15, c("Pane","Carne Bov","Zucchero"), lty=c(1,2,2),lwd=c(2,2,1),cex=1)

load("/Users/utente/Documents/prova r pkg/lucidi del corso/apple.RData")
apple$Date <- as.Date(apple$Date)
par(mar=c(4.5,4.5,1,1))
plot(apple$Date,apple$Close, main="",lwd=1,type="l",xlab="Tempo",ylab="Quotazioni Apple")

load("/Users/utente/Documents/prova r pkg/lucidi del corso/oilgasoline.RData")
tscoil <- ts(coil[,2:3], start = c(1999, 6),frequency = 12)
par(mar=c(1,1,1,1))
plot(tscoil, main="",lwd=1)

decAP <- decompose(AP, type = "multiplicative")
plot(decAP)
decAP <- decompose(AP, type = "additive")
plot(decAP)

decMatlab <- decompose(tsgtm[,2], type = "multiplicative")
plot(decMatlab)

set.seed(12345)
par(mfrow=c(3,1),mar=c(4.5,4.5,1,1))
sersim <- arima.sim(n = 300, list(ar = c(0.9, 0), ma = c(0, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
par(mfrow=c(3,1),mar=c(4.5,4.5,1,1))
sersim <- arima.sim(n = 300, list(ar = c(-0.9, 0), ma = c(0, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
par(mfrow=c(2,3),mar=c(4.5,4.5,1,1))
sersim <- arima.sim(n = 300, list(ar = c(0.6, 0.3), ma = c(0, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
sersim <- arima.sim(n = 300, list(ar = c(-0.6, -0.3), ma = c(0, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
par(mfrow=c(2,3),mar=c(4.5,4.5,1,1))
sersim <- arima.sim(n = 300, list(ar = c(0.6, -0.3), ma = c(0, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
sersim <- arima.sim(n = 300, list(ar = c(-0.6, 0.3), ma = c(0, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
par(mfrow=c(3,1),mar=c(4.5,4.5,1,1))
sersim <- arima.sim(n = 300, list(ar = c(0, 0), ma = c(0.9, 0)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
par(mfrow=c(2,3),mar=c(4.5,4.5,1,1))
sersim <- arima.sim(n = 300, list(ma = c(0.6, 0.3)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)
set.seed(12345)
sersim <- arima.sim(n = 300, list(ma = c(-0.6, -0.3)), sd = sqrt(2))
plot(sersim)
acfsim <- acf(sersim,20,type="correlation")
pacfsim <- pacf(sersim,20)

par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(LakeHuron,main="",lwd=1,ylab="Livelli Annuali del Lake Huron dal 1875 al 1972")
par(mfrow=c(2,1),mar=c(4.5,4.5,1,1))
acflak <- acf(LakeHuron,20,type="correlation")
pacflak <- pacf(LakeHuron,20)
stilak <- arima(LakeHuron,order=c(2,0,0))
library(lmtest)
coeftest(stilak)
par(mfrow=c(2,1),mar=c(4.5,4.5,1,1))
acf(stilak$residuals,20,type="correlation")
hist(stilak$residuals,col="gray",breaks = 15)

par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(nhtemp,main="",lwd=1,ylab="Temperature medie annuali in New Haven dal 1912 al 1971")
par(mfrow=c(2,1),mar=c(4.5,4.5,1,1))
acfnh <- acf(nhtemp,20,type="correlation")
pacfnh <- pacf(nhtemp,20)
stinh <- arima(nhtemp,order=c(1,0,2))
stinh
coeftest(stinh)
stinh <- arima(nhtemp,order=c(1,0,1))
stinh
coeftest(stinh)
par(mfrow=c(2,1),mar=c(4.5,4.5,1,1))
acf(stinh$residuals,20,type="correlation")
hist(stinh$residuals,col="gray",breaks = 15)

library(strucchange)
bptsgdp_it_bi <- breakpoints(tsgdp_it_bi[,3] ~ 1)
summary(bptsgdp_it_bi) 
ci_bptsgdp_it_bi <- confint(bptsgdp_it_bi)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(tsgdp_it_bi[,3],main="",ylab="GDP procapite, ricostruzione BI")
lines(bptsgdp_it_bi)
lines(ci_bptsgdp_it_bi)

library(tseries)
library(forecast)
decAP <- decompose(AP, type = "additive")
plot(decAP)
adf.test(AP, alternative ="stationary", k=12)
findbest <- auto.arima(AP)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(forecast(findbest,h=20))
fitAP <- arima(AP, order=c(0, 1, 1), list(order=c(0, 1, 0), period = 12))
fore <- predict(fitAP, n.ahead=24)
U <- fore$pred + 2*fore$se
L <- fore$pred - 2*fore$se
ts.plot(AP, fore$pred, U, L, col=c(1, 2, 4, 4), lty=c(1, 1, 2, 2))
legend("topleft", c("Actual", "Forecast", "Error Bounds (95% prediction interval)"),col=c(1,2,4),lty=c(1,1,2))

decR <- decompose(tsgtm[,1], type = "additive")
plot(decR)
adf.test(tsgtm[,1], alternative ="stationary", k=12)
findbest <- auto.arima(tsgtm[,1])
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(forecast(findbest,h=20))
fitR <- arima(tsgtm[,1], order=c(1, 1, 1), list(order=c(2, 0, 0), period = 12))
fore <- predict(fitR, n.ahead=24)
U <- fore$pred + 2*fore$se
L <- fore$pred - 2*fore$se
ts.plot(tsgtm[,1], fore$pred, U, L, col=c(1, 2, 4, 4), lty=c(1, 1, 2, 2),ylab="Google Trend - R")
legend("topleft", c("Actual", "Forecast", "Error Bounds (95%)"),col=c(1,2,4),lty=c(1,1,2))

bptsgtm <- breakpoints(tsgtm[,1] ~ 1)
summary(bptsgtm) 
ci_bptsgtm <- confint(bptsgtm)
par(mfrow=c(1,1),mar=c(4.5,4.5,1,1))
plot(tsgtm[,1],main="",ylab="Google Trend - R")
lines(bptsgtm)
lines(ci_bptsgtm)

library(forecast)
par(mar=c(4,4,1,1))
fcast <- hw(AP,h=36,seasonal = "multiplicative",exponential = T)
plot(fcast,main="",xlab="Tempo",ylab="Passeggeri",cex.lab=2)
fcast <- hw(AP,h=36,seasonal = "additive",dumped = T)
plot(fcast,main="",xlab="Tempo",ylab="Passeggeri",cex.lab=2)