##### SIMULAZIONE SERIE CON TREND E VERIFICA ACF E PACF #####
# scelgo i parametri
alpha = 2.0
beta  = 0.3

# Simulo la serie di 200 osservazioni
Sim = c()
for (i in 1:200){
  Sim[i] = alpha + beta * i + rnorm(1, 0, 2)  
}

plot.ts(Sim)                         # analisi grafica
abline(alpha, beta, col = "red")     # la retta "vera"

acf(Sim)      # AutoCorrelation Function
pacf(Sim)     # Partial AutoCorrelation Function

ESERCIZIO PER CASA: provare a stimare un modello di regressione e fare l'analisi dei residui
**********************************************************
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Possiamo simulare ogni processo nella classe dei processi ARIMA(p, i, q)
# usando la funzione arima.sim()

##### I PROCESSI A MEDIA MOBILE #####
# Iniziamo con i processi MA(1)
# simuliamo un processo invertibile con coefficiente positivo
Y1 = arima.sim(n=1000, list(ma= 0.8), sd = sqrt(1))
# e con coefficiente negativo
Y2 = arima.sim(n=1000, list(ma= -0.8), sd = sqrt(1))

par(mfrow=c(2,1))
ts.plot(Y1)
ts.plot(Y2)

# confrontiamo le acf dei due processi
acf(Y1)
acf(Y2)

# confrontiamo le pacf dei due processi
pacf(Y1)
pacf(Y2)

# possiamo ottenere lo stesso processo simulandolo senza la funzione arima
Ysim1 = c(); Ysim1[1] = 0
Ysim2 = c(); Ysim2[1] = 0
shock = rnorm(1000, 0, 1)

for (i in 2:1000){
  Ysim1[i] =  0.8 * shock[i-1] + shock[i]
  Ysim2[i] =  -0.8 * shock[i-1] + shock[i]
}

par(mfrow=c(2,1))
ts.plot(Ysim1)
ts.plot(Ysim2)

# confrontiamo le acf dei due processi
acf(Ysim1)
acf(Ysim2)

# confrontiamo le pacf dei due processi
pacf(Ysim1)
pacf(Ysim2)

# Per ogni rappresentazione MA(1) invertibile esiste una rappresentazione
# MA(1) non invertibile, di parametro > 1, che possiede la stessa ACF.
theta = 1/0.8
Y3 = arima.sim(n=1000, list(ma= theta), sd = sqrt(1/(theta^2)))
# e con coefficiente negativo
Y4 = arima.sim(n=1000, list(ma= -theta), sd = sqrt(1/(theta^2)))

par(mfrow=c(2,1))
ts.plot(Y3) # graficamente non si evidenzia nessuna anomalia
ts.plot(Y4)

# confrontiamo le acf dei due processi
acf(Y1, main='Processo invertibile')             # processo invertivbile
acf(Y3, main='Processo analogo non invertibile') # processo non invertibile

# confrontiamo le pacf dei due processi
pacf(Y2, main='Processo invertibile')             # processo invertivbile
pacf(Y4, main='Processo analogo non invertibile') # processo non invertibile

# I due processi hanno identiche proprietaŐ e dunque sarebbe impossibile 
# discriminarli a partire da una serie storica. Tale problema di identificazione 
# viene risolto appunto vincolando il parametro nell'intervallo (-1,+1).
# Il vincolo permette anche di riscrivere il processo come un AR(ě) stazionario:

# Passiamo ai processi MA(2)... possiamo avere quattro casi:
# 1. coefficienti entrambi positivi
Y5 = arima.sim(n=1000, list(ma= c(0.6, 0.3)), sd = sqrt(1))
# 2. coefficienti entrambi negativi
Y6 = arima.sim(n=1000, list(ma= c(-0.6, -0.3)), sd = sqrt(1))
# 3. coefficienti con segni alterni
Y7 = arima.sim(n=1000, list(ma= c(0.6, -0.3)), sd = sqrt(1))
# 4. coefficienti con segni alterni
Y8 = arima.sim(n=1000, list(ma= c(-0.6, 0.3)), sd = sqrt(1))
               
par(mfrow=c(2,2))
ts.plot(Y5)
ts.plot(Y6)
ts.plot(Y7)
ts.plot(Y8)

# confrontiamo le acf dei quattro processi
acf(Y5, main = 'coeff. positivi')
acf(Y6, main = 'coeff. negativi')
acf(Y7, main = 'coeff. alterni (+, -)')
acf(Y8, main = 'coeff. alterni (-, +)')

# confrontiamo le pacf dei quattro processi
pacf(Y5, main = 'coeff. positivi')
pacf(Y6, main = 'coeff. negativi')
pacf(Y7, main = 'coeff. alterni (+, -)')
pacf(Y8, main = 'coeff. alterni (-, +)')

# possiamo ottenere gli stessi processi simulandoli senza la funzione arima
Ysim5 = c(); Ysim5[1] = 0; Ysim5[2] = 0; 
Ysim6 = c(); Ysim6[1] = 0; Ysim6[2] = 0; 
Ysim7 = c(); Ysim7[1] = 0; Ysim7[2] = 0; 
Ysim8 = c(); Ysim8[1] = 0; Ysim8[2] = 0; 

for (i in 3:1000){
  Ysim5[i] =   0.6 * shock[i-1] + 0.3 * shock[i-2] + shock[i]
  Ysim6[i] =  -0.6 * shock[i-1] - 0.3 * shock[i-2] + shock[i]
  Ysim7[i] =   0.6 * shock[i-1] - 0.3 * shock[i-2] + shock[i]
  Ysim8[i] =  -0.6 * shock[i-1] + 0.3 * shock[i-2] + shock[i]
}

par(mfrow=c(2,2))
ts.plot(Ysim5)
ts.plot(Ysim6)
ts.plot(Ysim7)
ts.plot(Ysim8)

*******************************************************************
##### I PROCESSI AUTO REGRESSIVI #####
# Iniziamo con i processi AR(1)
*******************************************************************
Y1 = arima.sim(n=1000, list(ar= 0.9), sd = sqrt(1))
Y1bis = arima.sim(n=1000, list(ar= 0.2), sd = sqrt(1))
# e con coefficiente negativo
Y2 = arima.sim(n=1000, list(ar= -0.8), sd = sqrt(1))

par(mfrow=c(2,1))
ts.plot(Y1)
ts.plot(Y1bis)    # vediamo l'effetto del coefficiente sulla serie

ts.plot(Y1)
ts.plot(Y2)

# confrontiamo le acf dei due processi
acf(Y1)
acf(Y2)

# confrontiamo le pacf dei due processi
pacf(Y1)
pacf(Y2)

# possiamo ottenere lo stesso processo simulandolo senza la funzione arima
Ysim1 = c(); Ysim1[1] = 0
Ysim2 = c(); Ysim2[1] = 0

for (i in 2:100){
  Ysim1[i] =  0.8 * Ysim1[i-1] + rnorm(1, 0, 1)
  Ysim2[i] =  -0.8 * Ysim2[i-1] + rnorm(1, 0, 1)
}

par(mfrow=c(2,1))
ts.plot(Ysim1)
ts.plot(Ysim2)

# confrontiamo le acf dei due processi
acf(Ysim1)
acf(Ysim2)

# confrontiamo le pacf dei due processi
pacf(Ysim1)
pacf(Ysim2)

# La stazionariet… Š una caratteristica fondamentale che desideriamo da un processo AR
# se il processo Š non stazionario diventa impossibile studiarne le propriet…
Y3 = c(); Y3[1]=0
for (i in 2:30){
  Y3[i] =  1.1 * Y3[i-1] + rnorm(1, 0, 1)
}

par(mfrow=c(1,1))
ts.plot(Y3)

# esiste per• un particolare tipo di processo non stazionario chiamato Random Walk in cui
# il coefficiente che determina l'impatto di ogni osservazione sulla successiva Š esattamente 1
Y4 = c(); Y4[1]=0
for (i in 2:1000){
  Y4[i] =  1.0 * Y4[i-1] + rnorm(1, 0, 1)
}

par(mfrow=c(1,1))
ts.plot(Y4)

par(mfrow=c(2,1))
# confrontiamo le acf di un processo stazionario e un random walk
acf(Y1, main='Processo stazionario')
acf(Y4, main='Processo random walk')

# confrontiamo le pacf dei due processi
pacf(Y1, main='Processo stazionario')
pacf(Y4, main='Processo random walk')


# Passiamo ai processi AR(2)... possiamo avere quattro casi:
# Vedi triangolo stazionariet…
# 1. coefficienti con segni alterni (-, +)
Y5 = arima.sim(n=1000, list(ar= c(-0.6, 0.3)), sd = sqrt(1))
# 2. coefficienti entrambi positivi
Y6 = arima.sim(n=1000, list(ar= c(0.6, 0.3)), sd = sqrt(1))
# 3. coefficienti entrambi negativi e radici complesse
Y7 = arima.sim(n=1000, list(ar= c(-1.6, -0.8)), sd = sqrt(1))
# 4. coefficienti con segni alterni (+, -) e radici complesse
Y8 = arima.sim(n=1000, list(ar= c(1.6, -0.8)), sd = sqrt(1))
              
par(mfrow=c(2,2))
ts.plot(Y5)
ts.plot(Y6)
ts.plot(Y7)
ts.plot(Y8)

# confrontiamo le acf dei quattro processi
acf(Y5, main = 'segni alterni (-, +)')
acf(Y6, main = 'entrambi positivi')
acf(Y7, main = 'negativi e radici compl.')
acf(Y8, main = 'alterni (+, -) e radici compl.')

# confrontiamo le pacf dei quattro processi
pacf(Y5, main = 'segni alterni (-, +)')
pacf(Y6, main = 'entrambi positivi')
pacf(Y7, main = 'negativi e radici compl.')
pacf(Y8, main = 'alterni (+, -) e radici compl.')

# possiamo ottenere gli stessi processi simulandoli senza la funzione arima
Ysim5 = c(); Ysim5[1] = 0; Ysim5[2] = 0; 
Ysim6 = c(); Ysim6[1] = 0; Ysim6[2] = 0; 
Ysim7 = c(); Ysim7[1] = 0; Ysim7[2] = 0; 
Ysim8 = c(); Ysim8[1] = 0; Ysim8[2] = 0; 

for (i in 3:1000){
  Ysim5[i] =  -0.6 * Ysim5[i-1] + 0.3 * Ysim5[i-2] + rnorm(1, 0, 1)
  Ysim6[i] =   0.6 * Ysim6[i-1] + 0.3 * Ysim6[i-2] + rnorm(1, 0, 1)
  Ysim7[i] =  -1.6 * Ysim7[i-1] - 0.8 * Ysim7[i-2] + rnorm(1, 0, 1)
  Ysim8[i] =   1.6 * Ysim8[i-1] - 0.8 * Ysim8[i-2] + rnorm(1, 0, 1)
}

par(mfrow=c(2,2))
ts.plot(Ysim5)
ts.plot(Ysim6)
ts.plot(Ysim7)
ts.plot(Ysim8)

##### ESEMPI EMPIRICI DI PROCESSI AR #####
library(tseries)
data(NelPlo)

par(mfrow=c(3,1))
plot(cpi)
acf(cpi)    # Ci fa pensare ad un processo autoregressivo molto persistente
pacf(cpi)   # potrebbe essere di ordine 1

# come potrebbe essere l'analogo processo simulato?
YsimSpec = c(); YsimSpec[1] = 0

for (i in 2:100){
  YsimSpec[i] = 0.002 * YsimSpec[i-1] + 0.004*i^2 + rnorm(1, 0, 1)
}
ts.plot(YsimSpec)
acf(YsimSpec)
pacf(YsimSpec)

# Un altro esempio
plot(unemp)
acf(unemp)   # Potrebbe essere un AR(2) con radici complesse 
pacf(unemp)  # con segni alterni (+, -) come il caso Y8