# carico libreria "sse"
source("c:/discod/didattica/serie storiche/sse_v5.txt")

##########################################################
#  LA SERIE MENSILE DELL'INDICE DELLA PRODUZIONE INDUSTRIALE
##########################################################

setwd("c:/discod/didattica/serie storiche/ss2627/Esercitazioni_R/esercitazione_2/")

# carico i dati considerando solo il periodo 1990:01 - 2020:08
y0=ts(scan("indprod.txt"),start=c(1990,1),end=c(2020,8),frequency=12)
# seleziono periodo  1990:01 - 2008:08
y=window(y0, start=c(1990,1), end=c(2008,8))

plot(y0)
plot(y)

### Ipotizziamo un modello additivo con stagionalità di ordine 12 

# stima del trend, utilizzo una m.m. a 13 termini ottenuta come centratura
# di due m.m. a 12 termini. Perdiamo 6 osservazioni iniziali e 6 finali

T0.hat=filter(y, filter=c(1/24,rep(2/24,11),1/24))
T0.hat=ts(T0.hat, start=c(1990,1),end=c(2020,8),frequency=12)
plot(y)
lines(T0.hat, col="red")

# calcolo dei coefficienti specifici di stagionalità
CS=y-T0.hat
plot(CS,type="l")
plot(CS,type="b",cex=0.7)
CS

# metto i dati in una matrice con gli anni nelle righe e i
# mesi nelle colonne. Dato che mancano gli ultimi 4 mesi del 2008
# aggiungo 4 NA finali.
CS.mat=CS.matrice(CS,n.periodi=19, periodo=12, 
                   na.inizio=0, na.fine=4)

# boxplot delle distribuzioni dei mesi
boxplot(CS.mat,col=rainbow(12))

# plotto tutti i coeff specifi
nr=nrow(CS.mat)
plot(rep(1,nr),CS.mat[,1], xlim=c(0.5,12.5),ylim=c(-65,23),xlab="mesi", ylab="Coeff. specifici")
for (i in 2:12) points(rep(i,nr),CS.mat[,i])
      
# calcolo dei coefficienti grezzi
CG=apply(CS.mat,2,mean,na.rm=T)
CG

# inserisco i CG nel grafico dei CS
for (i in 1:12) points(i,CG[i], cex=1.2, pch=19, col="red")
       
# controllo la somma dei coefficienti grezzi
sum(CG)

# calcolo dei coefficienti ideali
CI=CG-mean(CG)   # rappresentano la compon. stag. S_t
plot(CI, type="b")

# inserisco CI in un vettore della lunghezza di y
CI.vettore=rep(CI,20)
CI.vettore=CI.vettore[1:length(y)]
plot(CI.vettore,type="l")

# serie destagionalizzata
yd.a=y-CI.vettore
plot(y,type="l")
lines(yd.a, type="l", col="red")

###### solo comandi strettamente necessari ######
#T0.hat=filter(y, filter=c(1/24,rep(2/24,11),1/24))
#T0.hat=ts(T0.hat, start=c(1990,1),end=c(2020,8),frequency=12)
#CS=y-T0.hat
#CS.mat=CS.matrice(CS,n.periodi=19, periodo=12, na.inizio=0, na.fine=4)
#CG=apply(CS.mat,2,mean,na.rm=T)
#CI=CG-mean(CG)   
#CI.vettore=rep(CI,20)
#CI.vettore=CI.vettore[1:length(y)]
#yd.a=y-CI.vettore
##################################################


### Ipotizziamo ora un modello moltiplicativo con stagionalit? di ordine 12 

# per la stima preliminare del ciclo-trend, utilizzo la stessa m.m. a 13 termini considerata
# nel modello additivo
T0.hat=filter(y,filter=c(1/24,rep(2/24,11),1/24))
T0.hat=ts(T0.hat, start=c(1990,1),end=c(2008,8),frequency=12)
plot(y)
lines(T0.hat, col="red")

# calcolo dei coefficienti specifici di stagionalità;
# ora dobbiamo dividere per T0.hat
CS=y/T0.hat
plot(CS,type="l")
plot(CS,type="b")
CS  # notare valore dei CS

# metto coeff specifici in una matrice
CS.mat=CS.matrice(CS,n.periodi=19, periodo=12, 
                  na.inizio=0, na.fine=4)

# boxplot delle distribuzioni dei mesi
boxplot(CS.mat, col=rainbow(12))

# plotto tutti i coeff specifi
nr=nrow(CS.mat)
plot(rep(1,nr),CS.mat[,1], xlim=c(0.5,12.5),ylim=c(0.4,1.2),xlab="mesi", ylab="Coeff. specifici")
for (i in 2:12) points(rep(i,nr),CS.mat[,i])
           
# calcolo dei coefficienti grezzi; tutto come per modello additivo
CG2=apply(CS.mat,2,mean,na.rm=T)
CG2

# inserisco nel grafico dei CS
for (i in 1:12) points(i,CG2[i], cex=1.5, pch=19, col="red")
                
# ora dobbiamo controllare se il prodotto dei coefficienti grezzi ? 1
prod(CG2)  # prodotto diverso da 1

# calcolo dei coefficienti ideali
CG2.m.geom=prod(CG2)^(1/12)    # ora devo usare media geometrica
CI=CG2/CG2.m.geom           # rappresentano la compon. stag. S_t
plot(CI, type="b")

# inserisco CI in un vettore della lunghezza di y
CI.vettore=rep(CI,20)
CI.vettore=CI.vettore[1:length(y)]

# serie destagionalizzata
yd.b=y/CI.vettore
plot(y,type="l",col="white",ylim=c(80,140))
lines(yd.b, type="l", col="red")

lines(yd.a,col="blue")


###### Replichiamo ora sull'intera serie il modello additivo
y=window(y0, start=c(1990,1), end=c(2020,8))
plot(y)
T0.hat=filter(y, filter=c(1/24,rep(2/24,11),1/24))
T0.hat=ts(T0.hat, start=c(1990,1),end=c(2020,8),frequency=12)
CS=y0-T0.hat
CS.mat=CS.matrice(CS,n.periodi=31, periodo=12, na.inizio=0, na.fine=)
CG=apply(CS.mat,2,mean,na.rm=T)
CI=CG-mean(CG)   
CI.vettore=rep(CI,31)
CI.vettore=CI.vettore[1:length(y)]
yd.a=y-CI.vettore
lines(yd.a, col="red")

plot(yd.a,col="red")
## NB - attenzione comunque ai cambi di regime!


########################################?
# Effetto Slutsky-Yule
#
# L'applicazione di una media mobile su una serie White Noise pu? creare autocorrelazioni 
# spurie. Per verificare empiricamente cosa succede, generiamo un dataset che contenga 
# una serie di 100 dati casuali con distribuzione White Noise Gaussiana di media 1 e 
# varianza 2.
  y=rnorm(5000,mean=0,sd=sqrt(2))
  plot(y, type="l")

# Proviamo ora ad applicare a tale serie una media mobile 
#centrata di ordine 5 con pesi rispettivamente 
# (1/8, 1/4, 1/4, 1/4, 1/8).
# mmy=(y(-2)+2*y(-1)+2*y+2*y(1)+y(2))/8

mmy=filter(y, filter=c(1/8,1/4,1/4,1/4,1/8))
mmy[1:5]
mmy[95:100]

theta=c(1/8,1/4,1/4,1/4,1/8)
RRV=sum(theta^2)
RRV
var(y)
var(mmy,na.rm=T)
var(mmy,na.rm=T)/var(y)

#####################################################################
##  DESTAGIONALIZZAZIONE: 
##  LA SERIE TRIMESTRALE DEI CONSUMI FINALI INTERNI DELLE FAMIGLIE
####################################################################

# Consideraimo la serie storicha dei consumi finali interni delle famiglie 
# contenuti nel dal file cfi.txt (Es. pag. 103 del testo). Sono dati trimestrali che coprono il periodo 1997:1 - 2001:4.

## la serie trimestrale dei consumi finali interni delle famiglie 
y=ts(scan("cfi.txt"),start=c(1997,1),end=c(2001,4),frequency=4)
plot(y)

# La serie ha periodicità 4, quindi dobbiamo utilizzare la
# m.m. a 5 che si ottiene dall'operazione di centratura di due
# m.m. non centrate e a 4 termini (pesi 1/8,1/4,1/4,1/4,1/8)

# Consideriamo dapprima un MODELLO ADDITIVO

# calcolo della media mobile 
T0.hat=filter(y, filter=c(1/8,rep(1/4,3),1/8))
T0.hat=ts(T0.hat, start=c(1997,1),end=c(2001,4),frequency=4)
plot(y)
lines(T0.hat,col="red")

# Coefficienti specifici di stagionalità
CS=y-T0.hat

# Coefficienti grezzi di stagionalità
CS.mat=CS.matrice(CS,n.periodi=5, periodo=4, na.inizio=0, na.fine=0)
CG=apply(CS.mat,2,mean,na.rm=T)

# controllo media dei CG
mean(CG)  # diversa da zero, anche se di poco

# Coefficienti ideali di stagionalità
CI=CG-mean(CG)   

CI.vettore=rep(CI,20)
CI.vettore=CI.vettore[1:length(y)]

# serie destagionalizzata

yd.a=y-CI.vettore

#grafico
plot(y)
lines(T0.hat,col="red")
lines(yd.a,col="blue")

# Consideriamo dapprima un MODELLO MOLTIPLICATIVO

# calcolo della media mobile 
T0.hat=filter(y, filter=c(1/8,rep(1/4,3),1/8))
T0.hat=ts(T0.hat, start=c(1997,1),end=c(2001,4),frequency=4)
plot(y)
lines(T0.hat,col="red")

# Coefficienti specifici di stagionalità
CS=y/T0.hat

# Coefficienti grezzi di stagionalità
CS.mat=CS.matrice(CS, n.periodi=5, periodo=4, na.inizio=0, na.fine=0)
CG=apply(CS.mat,2,mean,na.rm=T)

# controllo media dei CG
prod(CG)  # diversa da uno, anche se di poco

# Coefficienti ideali di stagionalità
CG.m.geom=prod(CG)^(1/4)    # ora devo usare media geometrica
CI=CG/CG.m.geom           # rappresentano la compon. stag. S_t
plot(CI, type="b")

CI.vettore=rep(CI,20)
CI.vettore=CI.vettore[1:length(y)]

# serie destagionalizzata

yd.b=y/CI.vettore

#grafico
plot(y)
lines(yd.b,col="red")
lines(yd.a,col="blue")