# figura2-1.txt

# setwd("c:/discod/didattica/serie storiche/ss2627/esercitazioni_R/esercitazione_1/")

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

############ ESEMPIO 1 - tasso di disoccupazione ################

y=scan("figura2-1.txt")
#win.graph()
plot(y, type="l")

n=length(y)
tempo=seq(1,n,1)  # t
tempo2=tempo^2    # t^2
tempo3=tempo^3    # t^3
tempo4=tempo^4    # t^4
# apro finestra grafica doppia
#x11(width=8,height=8)
#par(mfrow=c(2,1),mai=c(0.8,0.8,0.6,0.6))

### polinomio di grado zero

fit0=lm(y~1)
summary(fit0)

plot(y,type="l")
lines(fit0$fitted, col="red")  # valori stimati

plot(fit0$resid, type="l")     # residui
abline(h=0, col="black",lty=2)

# per avere i valori esatti  di R^2 e R^2 corretto
R2.0=summary(fit0)$r.squared
R2cor.0=summary(fit0)$adj.r.squared
R2.0
R2cor.0


### polinomio di grado uno

fit1=lm(y~1+tempo)
summary(fit1)

plot(y,type="l")
lines(fit1$fitted, col="red")  # valori stimati
plot(fit1$resid, type="l")     # residui
abline(h=0, col="black",lty=2) # linea dello zero

# per avere i valori esatti  di R^2 e R^2 corretto
R2.1=summary(fit1)$r.squared
R2cor.1=summary(fit1)$adj.r.squared
R2.1
R2cor.1


### polinomio di grado due

fit2=lm(y~1+tempo+tempo2)
summary(fit2)

plot(y,type="l")
lines(fit2$fitted, col="red")  # valori stimati
plot(fit2$resid, type="l")     # residui
abline(h=0, col="black", lty=2)

# per avere i valori esatti  di R^2 e R^2 corretto
R2.2=summary(fit2)$r.squared
R2cor.2=summary(fit2)$adj.r.squared
R2.2
R2cor.2


## polinomio di grado tre

fit3=lm(y~1+tempo+tempo2+ tempo3)  # provare a togliere tempo2
summary(fit3)

plot(y,type="l")
lines(fit3$fitted, col="red")  # valori stimati
plot(fit3$resid, type="l")     # residui
abline(h=0, col="black", lty=2)

# per avere i valori esatti  di R^2 e R^2 corretto
R2.3=summary(fit3)$r.squared
R2cor.3=summary(fit3)$adj.r.squared
R2.3
R2cor.3

## polinomio di grado quattro

fit4=lm(y~1+tempo+tempo2+tempo3+tempo4)
summary(fit4)

plot(y,type="l")
lines(fit4$fitted, col="red")  # valori stimati
plot(fit4$resid, type="l")     # residui
abline(h=0, col="black", lty=2)

# per avere i valori esatti  di R^2 e R^2 corretto
R2.4=summary(fit4)$r.squared
R2cor.4=summary(fit4)$adj.r.squared
R2.4
R2cor.4

plot(fit2$resid, type="l")     # residui
lines(fit4$resid, col="red")   # residui
abline(h=0)


R2=c(R2.0,R2.1,R2.2,R2.3,R2.4)
plot(R2,type="b")
R2

R2cor=c(R2cor.0,R2cor.1,R2cor.2,R2cor.3,R2cor.4)
plot(R2cor,type="b")
R2cor

#....non dimentichiamo che la serie di partenza
# contiene anche la componente stagionale che noi
# non abbiamo modellato...


# Proviamo ad invididuare il grado del polinomio metodo delle differenze successive
d1=diff(y) #diff(y,differences=1)
d2=diff(d1) #diff(y,differences=2)
d3=diff(d2) #diff(y,differences=3)
d4=diff(d3) #diff(y,differences=4)
plot(d1,type="l")
plot(d2,type="l")
plot(d3,type="l")
plot(d4,type="l")


########### ora modello con la sola stagionalità

n=length(y)
# creo le  variabili dummy
# start="stagione" da cui parte la serie
# comando make.dummy fa parte della libreria
# sse_v5.txt; non è comando standard di R
dv=make.dummy(n=n, freq=4, start=4)
colnames(dv)=c("d1t", "d2t", "d3t", "d4t")
dv


# calcolo dei coefficienti grezzi
fit.stag=lm(y~-1+dv[,1]+dv[,2]+dv[,3]+dv[,4])
summary(fit.stag)
gamma.star=fit.stag$coef
# se usiamo l'espressione dello stimatore OLS 
# otteniamo lo stesso risultato
gamma.star.ols=solve(t(dv)%*%dv)%*%t(dv)%*%y

# se avessimo messo la costante...
# DV=cbind(rep(1,n),dv)
# gamma.star.ols2=solve(t(DV)%*%DV)%*%t(DV)%*%y
# fit.stag2=lm(y~1+dv[,1]+dv[,2]+dv[,3]+dv[,4])
# summary(fit.stag2)

# media dei coefficienti grezzi
gamma.bar=mean(gamma.star)
gamma.bar

# coefficienti ideali ("veri") di stagionalità
gamma=gamma.star-gamma.bar
gamma

#gamma=gamma.star-mean(gamma.star)

# componente stagionale stimata
# Devo ripetere i coefficienti stagionali 
# per avere una stagionalità  lunga quanto la serie
# Il primo periodo è il quarto trimestre....l'ultimo
# è il secondo trimestre. Quindi:
St=c(gamma[4],rep(gamma,8),gamma[1],gamma[2])

# serie destagonalizzata
#
yd=y-St
plot(yd,type="l")


# valori stimati
y.stimati=fit.stag$fitted

# grafico
plot(y,type="l", xlab="t")
lines(y.stimati, col="red")

#################################################
## stima simultanea di trend e stagionalità
#################################################

fit.sim=lm(y~-1+tempo+tempo2+dv[,1]+dv[,2]+dv[,3]+dv[,4])
summary(fit.sim)


# grafico valori stimati e residui
#x11(width=8,height=8)
#par(mfrow=c(2,1),mai=c(0.8,0.8,0.6,0.6))

y.stimati=fit.sim$fitted
plot(y,type="l", xlab="t")
lines(y.stimati, col="red")

plot(fit.sim$res, type="l", xlab="t")
abline(h=0,lty=2)


# trend di grado quattro
fit.sim2=lm(y~-1+tempo+tempo2+tempo3+tempo4+dv[,1]+dv[,2]+dv[,3]+dv[,4])
summary(fit.sim2)

# confronto tra modelli
plot(y,type="l", xlab="t")
lines(fit.sim$fitted, col="red")
plot(y,type="l", xlab="t")
lines(fit.sim2$fitted, col="red")

# confronto tra gli errori modelli
plot(fit.sim$resid,type="l", xlab="t") # nero grado 2
lines(fit.sim2$resid, col="red")  # rosso grado 4
abline(h=0,lty=2)

# meglio trend di grado 2 o di grado 4?

R2.2=summary(fit.sim)$r.squared
R2.4=summary(fit.sim2)$r.squared

R2.cor.2=summary(fit.sim)$adj.r.squared
R2.cor.4=summary(fit.sim2)$adj.r.squared

R2.2
R2.4

R2.cor.2
R2.cor.4
# praticamente equivalenti

######################################################################
# Consideriamo il modello con la costante al posto della prima dummy
######################################################################

fit.sim.cost=lm(y~1+dv[,2]+dv[,3]+dv[,4]+tempo+tempo2)
summary(fit.sim.cost)

# coeff stimati ("tilde")
gamma.tilde=fit.sim.cost$coef[1:4]

# coeff grezzi (gamma star)
gamma.star=rep(0,4)
gamma.star[1]=gamma.tilde[1]
gamma.star[2]=gamma.tilde[2]+gamma.tilde[1]
gamma.star[3]=gamma.tilde[3]+gamma.tilde[1]
gamma.star[4]=gamma.tilde[4]+gamma.tilde[1]

# i parametri gamma.star sommano a zero? 
sum(gamma.star)   # no
# quindi sono coefficienti grezzi

# per ottenre i coeff ideali gli sottraggo la loro media
gamma.hat=gamma.star-mean(gamma.star)
sum(gamma.hat)  
# per la destagionalizzazione devo avere il vettore con
# i coef stagionali ideali per ogni t
St=c(gamma.hat[4],rep(gamma.hat,8),gamma.hat[1],gamma.hat[2])


# serie destagionalizzata
yd=y-St 
plot(y, type="l")
lines(yd,col="red")