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


########################
## esercizio 1
########################

# genero dati da un WNG(0,1.5)
set.seed(652)
wn100=rnorm(100,mean=0,sd=sqrt(1.5))
wn1000=rnorm(1000,mean=0,sd=sqrt(1.5))
plot(wn100,type="l")
plot(wn1000,type="l")
acf(wn100,lag.max=25)
limiti.100=qnorm(0.025)/sqrt(100)
limiti.100
pacf(wn100,25)
cfs(wn100,25,new=T)  # comando di sse
acf(wn1000,25,ylim=c(-0.1,0.1))
limiti.1000=qnorm(0.025)/sqrt(1000)
limiti.1000
pacf(wn1000,25, ylim=c(-0.1,0.1))
cfs(wn1000,25,new=T)  # comando di sse

#################
## esercizio 2
#################

# considero media mobile aritmetica a 5 termini
wn100.mm5=filter(wn100, filter=rep(1/5,5))
acf(wn100,25)
acf(wn100.mm5,lag.max=25, na.action=na.pass)
# la m.m. di dati incorrelati è autocorrelata

#####################
### Esercizio 3.
#####################

n=500
burn.out=50
y=rep(0,n+burn.out)
y[1]=0.1
y[2]=0.1
phi1=0.5
phi2=0.3
#phi1=0.25
#phi2=0.15
phi1=0.65
phi2=0.35
set.seed(123)
eps=rnorm(n+burn.out)
for (i in 3:(n+burn.out))
{
  y[i]=phi1*y[i-1]+phi2*y[i-2]+eps[i]
}
y=y[(burn.out+1):(n+burn.out)]

plot(y,type="l")
acf(y,20)
pacf(y,20,ylim=c(0,1))
cfs(y,20,new=T)
upper.limit=1.96/sqrt(n)
lower.limit=-1.96/sqrt(n)
upper.limit
lower.limit



mm=filter(y,rep(1/11,11))
plot(mm,type="l")
abline(h=mean(y))


###################
## Esercizio 4
###################

# autocorrelazioni teoriche
#ARMAacf(ar = c(phi1,phi2), ma =NULL, lag.max =lmax, pacf = FALSE)
# generazione da modelli ARMA
#arima.sim(n,model=list(order=c(2,0,1),ar=c(0.5,0.3), ma=c(0.4)), sd = 1)

###################################
# AR(2) con phi1=0.5, phi2=0.3
###################################
# funzioni di autocorrelazione teoriche
phi1=0.5
phi2=0.3
lmax=25
rho=ARMAacf(ar = c(phi1,phi2), ma =NULL, lag.max =lmax, pacf = FALSE)
plot(rho[1:(lmax+1)],type="h")
rho.part=ARMAacf(ar = c(phi1,phi2), ma =NULL, lag.max =lmax, pacf = TRUE)
plot(rho.part[1:lmax],type="h")

# generazione e autocorrelazioni campionarie
n=100
set.seed(1536)
ar2=arima.sim(n=n,model=list(order=c(2,0,0),
                             ar=c(0.5,0.3), 
                             ma=NULL), 
                             sd = 0.8)
acf(ar2,25)
pacf(ar2,25)
cfs(ar2,25,new=T)

## NOTA BENE: nella versione di R gli epsilon della 
## parte MA hanno segni positivi. Di conseguenza tutti
## i risultati delle autocorrelazioni degli MA devono
## essere correttamente interpretati
## Questo  vale sia per la simulazione che per la stima


###################################
# MA(1) con theta=0.6
###################################
# funzioni di autocorrelazione teoriche
theta=0.6
lmax=25
# NB - rho[1]=-theta/(1+theta^2)=-0.44, quindi ci aspetteremmo un 
# coefficiente negativo, invece...
rho=ARMAacf(ar = NULL, ma =theta, lag.max =lmax, pacf = FALSE)
plot(rho[1:lmax],type="h")
rho.part=ARMAacf(ar = NULL, ma =0.6, lag.max =lmax, pacf = TRUE)
plot(rho.part[1:lmax],type="h")
# Per avere lo stesso risultato delle slide bisogna mettere -0.6


# lo stesso per le autocorrelazioni campionarie
n=100
set.seed(1415)
ma1=arima.sim(n=n,model=list(order=c(0,0,1),ar=NULL, ma=theta), sd = 0.8)
acf(ma1,25)
pacf(ma1,25)
cfs(ma1,25,new=T)

###################################
# MA(2) con theta1=0.5, theta2=0.4
###################################
# funzioni di autocorrelazione teoriche
theta1=0.5
theta2=0.4
lmax=25
rho=ARMAacf(ar = NULL, ma =c(theta1,theta2), lag.max =lmax, pacf = FALSE)
plot(rho[1:lmax],type="h")
rho.part=ARMAacf(ar = NULL, ma =c(theta1,theta2), lag.max =lmax, pacf = TRUE)
plot(rho.part[1:lmax],type="h")

# lo stesso per le autocorrelazioni campionarie
n=150 # provare n=1000
set.seed(1453)  #1453
ma2=arima.sim(n=n,model=list(order=c(0,0,2),ar=NULL, ma=c(theta1,theta2)), sd = 0.8)
acf(ma2,25)
pacf(ma2,25)
cfs(ma2,25,new=F)

###################################
# ARMA(1,1) 
###################################
# funzioni di autocorrelazione teoriche
phi1=0.5
theta1=0.5
lmax=25
rho=ARMAacf(ar = phi1, ma =theta1, lag.max =lmax, pacf = FALSE)
plot(rho[1:lmax],type="h")
rho.part=ARMAacf(ar = phi1, ma =c(theta1), lag.max =lmax, pacf = TRUE)
plot(rho.part[1:lmax],type="h")


phi1=0.5
theta1=-0.5  
lmax=25
rho=ARMAacf(ar = phi1, ma =c(theta1), lag.max =lmax, pacf = FALSE)
plot(rho[1:lmax],type="h")
rho.part=ARMAacf(ar = phi1, ma =c(theta1), lag.max =lmax, pacf = TRUE)
plot(rho.part[1:lmax],type="h")
# in questo  caso l'arma(1,1) si può semplificare e coincide con un wn


# lo stesso per le autocorrelazioni campionarie
phi1=0.5
theta1=0.5
n=120 # provare n=1000
set.seed(28) 
arma11=arima.sim(n=n,model=list(order=c(1,0,1), 
                  ar=phi1, ma=theta1), sd = 0.8)
#arma11=arima.sim(n=n,model=list(ar=phi1, ma=theta1), sd = 0.8)
acf(arma11,25)
pacf(arma11,25)
cfs(arma11,25,new=F)


######################################
## Esercizio 5 - proprietà delle stime
######################################

phi1=0.5
theta1=-0.5
sigma2=1.2
n=80
var.teor=sigma2*((1+(theta1^2)-2*phi1*theta1)/(1-phi1^2))
var.teor
rho1.teor=((phi1-theta1)*(1-theta1*phi1))/
                  (1+(theta1^2)-2*phi1*theta1)
rho1.teor
#rho1.teor.2=ARMAacf(ar = phi1, ma =c(-theta1), lag.max =2, pacf = FALSE)[2]
#rho1.teor.2
M=100
n=200
var.stim=rep(0,M)
rho1.stim=rep(0,M)
for (i in 1:M)
{
  # si noti che, poiché R usa la verisone MA con il segno più, devo considerare
  # un valore di theta1 opposto a quello ipotizzato, cioè -theta1
  arma11=arima.sim(n=n,model=list(order=c(1,0,1),ar=phi1, ma=c(-theta1)), sd = sqrt(sigma2))
  var.stim[i]=var(arma11)
  rho1.stim[i]=acf(arma11,1)$acf[2]
}
var.stim
rho1.stim
# confronto tra varianze
var.teor
mean(var.stim)
abs(mean(var.stim)-var.teor)
# confronto tra rho1
rho1.teor
mean((rho1.stim))
abs(mean(rho1.stim)-rho1.teor)


#####################################
## Esercizio 6 -  Calcolo della PACF
#####################################


phi1=0.6
theta1=0.5
n=500 # provare n=1000
set.seed(38) #28 
arma11=arima.sim(n=n,model=list(order=c(1,0,1), 
                                ar=phi1, ma=theta1), sd = 0.8)
y=arma11[4:n]
y1=arma11[3:(n-1)]
y2=arma11[2:(n-2)]
y3=arma11[1:(n-3)]
kk=pacf(y,5)

P1=kk$acf[1]
P2=kk$acf[2]
P3=kk$acf[3]
c(P1,P2,P3)


mod1=lm(y~-1+y1)
phi11=mod1$coef

mod2=lm(y~-1+y1+y2)
phi22=mod2$coef[2]

mod3=lm(y~-1+y1+y2+y3)
phi33=mod3$coef[3]

c(phi11,phi22,phi33)
c(P1,P2,P3)