Tab7.1 <- read.table(file="Tab7.1.txt", h=T, sep="\t", dec=",")

## esempi di visualizzazione

## creiamo un oggetto 'ts' (dati+etichette temporali, o "timestamps")
vend <- ts(Tab7.1$vendite, start=c(2005, 1), freq=12)

## grafico standard
plot(vend)

## monthplot
monthplot(vend)
monthplot(stl(vend, "periodic"), choice="seasonal")
monthplot(stl(vend, "periodic"), choice="trend")

plot(decompose(vend))

## correlogramma
plot(acf(vend))

## es. autocorrelazione al secondo ritardo
N <- length(vend)
cov(vend[1:(N-2)], vend[3:N])/var(vend)

###############################################

### Esempio da Tab. 7.3
Tab7.3 <- read.table(file="Tab7.3.txt", h=T, sep="\t", dec=",")

plot(Tab7.3$vendite, type="l", lwd=2)


## seasonal plot:
anni <- unique(Tab7.3$anno)
plot(Tab7.3[Tab7.3$anno == anni[1], "vendite"], type="l",
     ylim=c(min(Tab7.3$vendite), max(Tab7.3$vendite)))
for(i in 2:12) lines(Tab7.3[Tab7.3$anno == anni[i], "vendite"], col=i)


## Stima del trend con il metodo "naif": regressione contro il tempo
## crea i dati annuali sommando i mensili per ogni anno
vend.anno <- tapply(Tab7.3$vendite, Tab7.3$anno, sum)

## OLS vend.anno contro t
trendmod <- lm(vend.anno ~ I(1:length(vend.anno)))
summary(trendmod)
plot(vend.anno, type="l", lwd=2)

## ispezioniamo l'adattamento
lines(fitted(trendmod), col="red", lwd=2)

## meglio una quadratica?
trendmod2 <- lm(vend.anno ~ I(1:length(vend.anno))+I((1:length(vend.anno))^2))
summary(trendmod2)
lines(fitted(trendmod2), col="green3", lty=2, lwd=2)


## alternativamente: stima del trend sui dati mensili con dummies-mese
plot(Tab7.3$vendite, type="l", lwd=2)

## OLS vend contro t, trend lineare
trendmod.m <- lm(vendite ~ I(1:length(vendite))+as.factor(mese), data=Tab7.3)
summary(trendmod.m)

## ispezioniamo l'adattamento
lines(fitted(trendmod.m), col="red", lwd=2)

## meglio una quadratica?
trendmod2.m <- lm(vendite ~ I(1:length(vendite))+I((1:length(vendite))^2)+
                as.factor(mese), data=Tab7.3)
summary(trendmod2.m)
lines(fitted(trendmod2.m), col="green3", lwd=2)

## residui:
par(mfrow=c(2,1))
plot(resid(trendmod.m), col="blue", lwd=2, type="l")
abline(h=0, col="orange")
plot(resid(trendmod2.m), col="blue", lwd=2, type="l")
abline(h=0, col="orange")

## alternativamente: detrending con medie mobili

## una semplice funzione "media mobile" autoprodotta:
mm12 <- function(x) {
  mx <- rep(NA, length(x))
  for(i in 7:(length(x)-6)) {
    mx[i] <- (sum(x[(i-5):(i+5)])+0.5*sum(x[c(i-6, i+6)]))/12
  }
  return(mx)
}

vend.t <- mm12(Tab7.3$vendite)

plot(Tab7.3$vendite, type="l")
lines(vend.t, col="red", lwd=2)

## detrend
vend.dt1 <- Tab7.3$vendite - vend.t

plot(vend.dt1, type="l", col="blue")

## calcola coefficienti stagionali:
Shat <- tapply(vend.dt1, Tab7.3$mese, mean, na.rm=TRUE)

## verifica se la media è zero
mean(Shat)
S2hat <- Shat-mean(Shat)


####### metodi per serie storiche (oggetti 'ts')

t73<-ts(Tab7.3$vendite, start=c(2000,1), freq=12)
plot(t73)


## scomposizione della serie 
plot(decompose(t73))

## singole componenti:
t73.t <- decompose(t73)$trend
plot(t73.t)

## serie grezza e componente di trend (rosso)
plot(t73)
lines(t73.t, col="red", lwd=2)

## componente stagionale
t73.s <- decompose(t73)$seasonal
plot(t73.s)

##
par(mfrow=c(2,1))

## previsione "automatica" delle serie storiche
## (Hyndman and Khandakar, JSS 27/2008)
library(forecast)

## livellamento esponenziale, specificazione "automatica"
myets <- ets(t73)
summary(myets)

plot(t73)
lines(fitted(myets), col="red")

forecast(myets, 12)
plot(forecast(myets, 12))


## "postmortem assessment"

## dimentica l'ultimo anno
plot(t73, col="white")
lines(window(t73, start=c(2000, 1), end=c(2010, 12)))

myets0 <- ets(window(t73, start=c(2000, 1), end=c(2010, 12)))
## previsioni:
lines(forecast(myets0, 12)$mean, col="blue", lwd=2)
## dati effettivamente osservati:
lines(window(t73, start=c(2011, 1), end=c(2011, 12)), col="red", lwd=2)

## idem: ARIMA
mya0 <- auto.arima(window(t73, start=c(2000, 1), end=c(2010, 12)))
summary(mya0)
lines(forecast(mya0, 12)$mean, col="green3", lwd=2)


#### Tabella 7.4 (pezzi del ricambio "C" per settimana)

Tab7.4 <- read.table(file="Tab7.4.txt", h=T, sep="\t", dec=",")

pezzi <- ts(Tab7.4$pezzi)
plot(pezzi)

## previsione con ETS "automatico"
pets <- ets(pezzi)
lines(fitted(pets), col="red")

plot(forecast(pets), 10)

plot(c(pezzi, rep(NA, 10)), type="l", ylim=c(900, 1800))
lines(c(rep(NA, 100), forecast(pets)$mean), col="red")

## interpolazione "naif" del trend
tmod <- lm(Tab7.4[55:100, "pezzi"] ~ I(1:46))
lines(c(rep(NA, 54), fitted(tmod)), col="green3", lty=2, lwd=2)
## previsione "naif"
lines(c(rep(NA, 100), 
        coef(tmod)[[1]] + coef(tmod)[[2]]*47:56),
      col="blue", lwd=2)