Dans ce TP, on utilise les librairies :

library(caschrono)
library(forecast)
library(astsa)

EXERCICE 1

Question 1

L’équation est celle d’un processus AR(2)

AR2 <- arima.sim(list(ar=c(1.5,-0.75)), n=200); plot(AR2)

Ce processus est stationnaire car les racines \(r_1\) et \(\bar{r_1}\) du polynome associé \(P(z)=0.75z^2-1.5z+1\) sont de module différent de 1 :

D <- 1.5^2-4*0.75; r1 <- (1.5+sqrt(-D)*1i)/(2*0.75) ; abs(r1)
## [1] 1.154701

Question 2

#acf(AR2,10)
#pacf(AR2,10)
acf2y(AR2,10)

##       LAG        ACF1         PACF
##  [1,]   1  0.85603546  0.856035455
##  [2,]   2  0.52616947 -0.773295960
##  [3,]   3  0.13277324 -0.037102045
##  [4,]   4 -0.21438653 -0.029052612
##  [5,]   5 -0.44631054 -0.071530177
##  [6,]   6 -0.53575016 -0.036093056
##  [7,]   7 -0.48730469  0.027126412
##  [8,]   8 -0.35079430 -0.140712374
##  [9,]   9 -0.17891370 -0.006098879
## [10,]  10 -0.00757228  0.078017548

L’ACF ne semble pas être celle d’un MA, et la PACF ayant deux “piques” correspond à celle d’un AR(2).

Question 3

On teste un AR(4) :

ARtest <- arima.sim(list(ar=c(0.4,-0.75,0.5,-0.7)),n=1000)
acf2y(ARtest,10)

##       LAG        ACF1         PACF
##  [1,]   1 -0.06667859 -0.066678594
##  [2,]   2 -0.44566190 -0.452118064
##  [3,]   3  0.37820741  0.384961468
##  [4,]   4 -0.25902387 -0.663197174
##  [5,]   5 -0.52566172 -0.075178763
##  [6,]   6  0.46187512  0.014527176
##  [7,]   7  0.18316606  0.034646558
##  [8,]   8 -0.31032817  0.015894247
##  [9,]   9  0.31466680 -0.011568385
## [10,]  10  0.10546629  0.005375511

L’ACF nous montre clairement que ce n’est pas un MA, et la PACF indique un AR(4), car 4 piques en dehors de l’intervalle de confiance puis plus rien.

Question 4

On commence toujours par une rapide étude graphique :

ARMA21 <- arima.sim(list(ar = c(0.89, -0.5), ma = -0.23), n=1000) ; plot(ARMA21)

acf2y(ARMA21,10)

##       LAG        ACF1          PACF
##  [1,]   1  0.47055203  0.4705520343
##  [2,]   2 -0.11506388 -0.4321749355
##  [3,]   3 -0.33611263 -0.0871954998
##  [4,]   4 -0.18849998  0.0372538536
##  [5,]   5  0.02998299 -0.0009542228
##  [6,]   6  0.12932678  0.0155905941
##  [7,]   7  0.07389693 -0.0201419284
##  [8,]   8 -0.02455407 -0.0120333818
##  [9,]   9 -0.11643660 -0.0839408176
## [10,]  10 -0.10619936 -0.0076294751

On estime l’ordre avec auto.arima et puis on étudie résidus :

auto.arima(ARMA21)
## Series: ARMA21 
## ARIMA(2,0,2) with zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       0.7481  -0.5290  -0.1084  0.0992
## s.e.  0.0721   0.0405   0.0797  0.0601
## 
## sigma^2 estimated as 0.9793:  log likelihood=-1406.81
## AIC=2823.62   AICc=2823.68   BIC=2848.16
Res <- auto.arima(ARMA21)$residuals
acf2y(Res,10) 

##       LAG          ACF1          PACF
##  [1,]   1 -0.0002302169 -0.0002302169
##  [2,]   2  0.0010053775  0.0010053245
##  [3,]   3 -0.0067278785 -0.0067274230
##  [4,]   4  0.0133738555  0.0133704151
##  [5,]   5 -0.0181590727 -0.0181448727
##  [6,]   6  0.0006331820  0.0005611278
##  [7,]   7 -0.0367027314 -0.0365094739
##  [8,]   8  0.0167681818  0.0163849149
##  [9,]   9 -0.0536310642 -0.0532279039
## [10,]  10 -0.0409455187 -0.0418541688
Box.test(Res)
## 
##  Box-Pierce test
## 
## data:  Res
## X-squared = 5.3e-05, df = 1, p-value = 0.9942

Conclusion : Le test de Box-Pierce ne permet pas de rejeter l’hypothèse d’un bruit blanc. La modélisation est donc pertinente.

Exercice 2

help(oil)
plot(oil)

Ce chronogramme semble avoir une tendance : ça ne ressemble pas à une réalisation d’un processus stationnaire.

# Calcul des log-rendements 
r.oil <- diff(log(oil))
plot(r.oil)

Ce chronogramme ressemble plus à la réalisation d’un processus stationnaire que le précédent. La tendance a été lisée.

acf2y(r.oil,20)

##       LAG         ACF1          PACF
##  [1,]   1  0.131303331  0.1313033315
##  [2,]   2 -0.067151512 -0.0858725683
##  [3,]   3  0.134740462  0.1594958503
##  [4,]   4 -0.007416312 -0.0596497008
##  [5,]   5  0.015468808  0.0542317276
##  [6,]   6 -0.033768496 -0.0768873667
##  [7,]   7 -0.034713442 -0.0003380219
##  [8,]   8  0.129331408  0.1206751658
##  [9,]   9  0.075972545  0.0517216606
## [10,]  10  0.021104848  0.0317191090
## [11,]  11  0.013195204 -0.0231050503
## [12,]  12  0.003433059 -0.0014644320
## [13,]  13 -0.016836043 -0.0336696012
## [14,]  14  0.061709186  0.0871322846
## [15,]  15 -0.052122650 -0.0749460361
## [16,]  16 -0.085654021 -0.0566749318
## [17,]  17  0.032952206  0.0060030598
## [18,]  18  0.047620321  0.0448886395
## [19,]  19 -0.045310214 -0.0472274611
## [20,]  20 -0.067042969 -0.0538996946
# Modélisation avec un ARMA
auto.arima(r.oil, stationary=TRUE, seasonal=FALSE)
## Series: r.oil 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##           ar1     ma1
##       -0.5253  0.7142
## s.e.   0.0872  0.0683
## 
## sigma^2 estimated as 0.002112:  log likelihood=904.58
## AIC=-1803.15   AICc=-1803.11   BIC=-1790.25
# Etude des résidus de la modélisation
ResOil <- auto.arima(r.oil, stationary=TRUE, seasonal=FALSE)$residuals
acf2y(ResOil) 

##       LAG         ACF1          PACF
##  [1,]   1 -0.018080995 -0.0180809946
##  [2,]   2  0.005879764  0.0055546576
##  [3,]   3  0.086546162  0.0867833650
##  [4,]   4  0.006777385  0.0099812647
##  [5,]   5  0.001092699  0.0002854634
##  [6,]   6 -0.011439596 -0.0191994585
##  [7,]   7 -0.055006498 -0.0575393113
##  [8,]   8  0.135755564  0.1351362760
##  [9,]   9  0.041865787  0.0513482077
## [10,]  10  0.034423461  0.0451948968
## [11,]  11 -0.006297216 -0.0294973887
## [12,]  12  0.022204374  0.0090662076
## [13,]  13 -0.037243228 -0.0467053757
## [14,]  14  0.074619955  0.0805743501
## [15,]  15 -0.058910616 -0.0424396009
## [16,]  16 -0.066082718 -0.0763687593
## [17,]  17  0.023485310 -0.0007418828
## [18,]  18  0.049753157  0.0469799741
## [19,]  19 -0.048030371 -0.0313131285
## [20,]  20 -0.047934253 -0.0600132292
## [21,]  21  0.017451110  0.0238645761
## [22,]  22  0.102294838  0.0884126905
## [23,]  23 -0.070880517 -0.0557195132
## [24,]  24 -0.043816919 -0.0352702550
## [25,]  25 -0.089168216 -0.0997125787
## [26,]  26  0.033080643  0.0240684600
## [27,]  27 -0.111645902 -0.0988630381
## [28,]  28 -0.046490583 -0.0201703596
## [29,]  29  0.007408309  0.0152872635
## [30,]  30 -0.006171383 -0.0116656792
## [31,]  31 -0.039181566 -0.0437828405
## [32,]  32 -0.026896855 -0.0454250556
## [33,]  33 -0.045730018 -0.0019079711
## [34,]  34  0.030869790  0.0407717776
## [35,]  35 -0.104857955 -0.0832226429
## [36,]  36  0.019939153  0.0114336556
## [37,]  37 -0.059910906 -0.0352408758
## [38,]  38  0.005485258  0.0282566619
## [39,]  39  0.001750034 -0.0028178992
## [40,]  40  0.001855440 -0.0105790238
Box.test(ResOil)
## 
##  Box-Pierce test
## 
## data:  ResOil
## X-squared = 0.17785, df = 1, p-value = 0.6732

Conclusion : une modélisation des log-rendements avec un ARMA(1,1) nous rend des résidus qui passent le test de blancheur de Box-Pierce. La modélisation est donc raisonnable.

Exercice 3 (algorithme Durbin-Levinson)

Voici une façon de coder l’algorithme.

DLalgo <- function(x,p){
    ## on nomme la fonction d'autocovariance empirique
 gamma <- function(h){acf(x,h,plot=F,type="covariance")$acf[h+1]}
 
  ## initiatilisation
 k <- 1
 kappa = gamma(1)/gamma(0)
 sigma2 = gamma(0)*(1-kappa^2)
 phi <- 1:p ; phi[1] <- kappa
 
  ## algorithme
 while (k <= p){
   
   ### calcul annexe de la somme
   S = 0 
   for (j in 1:k){ S <- S + phi[j] * gamma(k-j+1) }
   
   ### kappa 
   kappa <- (gamma(k+1)-S)/sigma2
   
   ### phi
   for (j in 1:k){ phi[j] <- phi[j] - kappa * phi[k-j+1] }
   phi[k+1] <-kappa
   
   ### sigma2
   sigma2 <- sigma2 * (1-kappa^2)
   
   ### on incremente la boucle
  k <- k+1
 }
   ## resultats
   Result = list(phi,sigma2)
   return(Result)
 }

Test

On retrouve les paramètres du modèle AR(2) de la Exercice 1, question 1 :

DLalgo(AR2,1)
## [[1]]
## [1]  1.518004 -0.773296
## 
## [[2]]
## [1] 0.8935402