### TP5

### Exercice 3.1

# Question 1

# Simulation du modèle ARMA(1,1)

N=500
y=0*(1:N) #y=rep(0,N)
xi=rnorm(N)

y[1]=xi[1]
for (i in 2:N){
  y[i]=0.6*y[i-1] +0.5*xi[i-1]+ xi[i] 
}

# Question 2

# Représentation de la série simulée, des autocorrélations et des autocorrélations partielles sur une même fenêtre

par(mfrow=c(3,1))
plot(1:N,y, type = "l", col = "red") 
acf(y)
pacf(y)

# Question 3

# Ajustement d'un MA(1) sur la série simulée

est1<-arima(y,order=c(0,0,1)) 
param=est1$coef
param

# Question 4

# Ajustement d'un AR(1) sur la série simulée

est2<-arima(y,order=c(1,0,0)) 

# Question 5

# Ajustement d'un ARMA(1,1) sur la série simulée

est3<-arima(y,order=c(1,0,1))

# Question 6

# Comparaison des trois modèles 

AIC(est1,est2,est3)
BIC(est1,est2,est3)

# install.packages("forecast")
library(forecast)

# Choix du modèle ARIMA pour la série y par la fonction auto.sim

auto.arima(y)

### Exercice 3.2

# Question 2

# Simulation d'un AR(2)

y=arima.sim(n = 1000, list(ar = c(0.8897, -0.4858)), sd = 1) 
par(mfrow=c(3,1))
plot(y, type = "l", col = "red") 
acf(y)
pacf(y)

# Simulation d'un MA(2)

y=arima.sim(n = 1000, list( ma = c(-0.2279, 0.2488)), sd = 1)
par(mfrow=c(3,1))
plot(y, type = "l", col = "red") 
acf(y)
pacf(y)

# Simulation d'un ARMA(2,2)

y=arima.sim(n = 1000, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)), sd = 1)
par(mfrow=c(3,1))
plot(1:length(y),y, type = "l", col = "red")
acf(y)
pacf(y)

# Simulation d'un ARIMA(2,1,2)

y=arima.sim(n = 1000, list(order = c(2,1,2),ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)), sd = 1)
par(mfrow=c(3,1))
plot(1:length(y),y, type = "l", col = "red")
acf(y)
pacf(y)

plot(diff(y))

### Exercice 3.3

# Question 1

serie=scan("serie1.dat")
plot(serie,type = "l")

# La série est clairement non stationnaire, on observe une tendance 
# ou l'effet d'une marche aléatoire. On ne peut donc pas la modéliser 
# par un ARMA.

serie_diff=diff(serie)
plot(serie_diff,t="l")

par(mfrow=c(3,1))
plot(serie_diff,t="l")
acf(serie_diff)
pacf(serie_diff)

# Choix d'un modèle AR(p) avec un critère d'information.

est1<-arima(serie_diff,order=c(1,0,0)) 
est2<-arima(serie_diff,order=c(2,0,0)) 
est3<-arima(serie_diff,order=c(3,0,0))

AIC(est1,est2,est3)
BIC(est1,est2,est3)

model1=est2 

shapiro.test(model1$residuals)
Box.test(model1$residuals,lag=20)

# Une grande p-valeur, les résidus sont donc indépendants.

# Choix d'un modèle MA(q) avec un critère d'information.

est1<-arima(serie_diff,order=c(0,0,1)) 
est2<-arima(serie_diff,order=c(0,0,2)) 
est3<-arima(serie_diff,order=c(0,0,3))

AIC(est1,est2,est3)
BIC(est1,est2,est3)

model2=est3 

shapiro.test(model2$residuals)
Box.test(model2$residuals,lag=20)

# On rejette H0, les résidus sont donc dépendants.
# Au final entre le modèle AR(2) et le MA(3), on retient le AR(2) 
# qui est valide (la structure des résidus correspondants semble en adéquation
# avec la structure d'un bruit blanc).

auto.arima(serie)

# La fonction auto.arima confirme bien le choix précédent.