TP 5
Adrien Hardy
library(tseries)
library(forecast)
library(caschrono)
library(fGarch)
# set.seed(123) # je fixe les réalisations aléatoiresEXERCICE I
?garchSim Modèle ARCH(1)
spec = garchSpec(model = list(omega = 1, alpha = c(0.8), beta = 0))
X <- garchSim(spec, n = 500)
plot(X)
acf(X)
Box.test.2(X,1:20)## Retard p-value
## [1,] 1 0.37479047
## [2,] 2 0.05804548
## [3,] 3 0.00527414
## [4,] 4 0.00907868
## [5,] 5 0.01663220
## [6,] 6 0.00218158
## [7,] 7 0.00428588
## [8,] 8 0.00389009
## [9,] 9 0.00701217
## [10,] 10 0.00851219
## [11,] 11 0.01401100
## [12,] 12 0.02148134
## [13,] 13 0.03023286
## [14,] 14 0.01787048
## [15,] 15 0.02352157
## [16,] 16 0.01000353
## [17,] 17 0.01132053
## [18,] 18 0.00830769
## [19,] 19 0.01224506
## [20,] 20 0.01465810Residus <- residuals(arima(X^2, order=c(1,0,0)))
acf(Residus)
Box.test.2(Residus,1:20)## Retard p-value
## [1,] 1 0.25364361
## [2,] 2 0.00838127
## [3,] 3 0.01475421
## [4,] 4 0.01760368
## [5,] 5 0.01954620
## [6,] 6 0.03614874
## [7,] 7 0.05977032
## [8,] 8 0.05980713
## [9,] 9 0.08535465
## [10,] 10 0.12450461
## [11,] 11 0.17193081
## [12,] 12 0.22437568
## [13,] 13 0.27086284
## [14,] 14 0.32706449
## [15,] 15 0.39487266
## [16,] 16 0.39987510
## [17,] 17 0.46561077
## [18,] 18 0.53364359
## [19,] 19 0.58606819
## [20,] 20 0.63303869Modèle ARCH(3)
spec = garchSpec(model = list(omega = 1, alpha = c(0.5,0.1,0.3), beta = 0))
X <- garchSim(spec, n = 500)
plot(X)
acf(X)
Residus <- residuals(arima(X^2, order=c(3,0,0)))
acf(Residus)
Box.test.2(Residus,1:20)## Retard p-value
## [1,] 1 0.9644953
## [2,] 2 0.9870417
## [3,] 3 0.9986589
## [4,] 4 0.9986184
## [5,] 5 0.9973485
## [6,] 6 0.9973354
## [7,] 7 0.9957780
## [8,] 8 0.7728837
## [9,] 9 0.7418514
## [10,] 10 0.6389294
## [11,] 11 0.4965860
## [12,] 12 0.5647410
## [13,] 13 0.6434901
## [14,] 14 0.4620679
## [15,] 15 0.4909645
## [16,] 16 0.4100683
## [17,] 17 0.4532260
## [18,] 18 0.5219448
## [19,] 19 0.5702670
## [20,] 20 0.6331217Modèle GARCH(1,1)
spec = garchSpec(model = list(omega = 1, alpha = c(0.5), beta = c(0.4)))
X <- garchSim(spec, n = 500)
plot(X)
acf(X)
Residus <- residuals(arima(X^2, order=c(1,0,1)))
acf(Residus)
Box.test.2(Residus,1:20)## Retard p-value
## [1,] 1 0.28732449
## [2,] 2 0.00140039
## [3,] 3 0.00000000
## [4,] 4 0.00000000
## [5,] 5 0.00000000
## [6,] 6 0.00000000
## [7,] 7 0.00000000
## [8,] 8 0.00000000
## [9,] 9 0.00000000
## [10,] 10 0.00000000
## [11,] 11 0.00000000
## [12,] 12 0.00000000
## [13,] 13 0.00000000
## [14,] 14 0.00000000
## [15,] 15 0.00000000
## [16,] 16 0.00000000
## [17,] 17 0.00000000
## [18,] 18 0.00000000
## [19,] 19 0.00000000
## [20,] 20 0.00000000Modèle GARCH(2,3)
spec = garchSpec(model = list(omega = 1, alpha = c(0.06,0.1), beta = c(0.2,0.15,0.1)))
X <- garchSim(spec, n = 500)
plot(X)
acf(X)
Residus <- residuals(arima(X^2, order=c(3,0,3)))
acf(Residus)
Box.test.2(Residus,1:20)## Retard p-value
## [1,] 1 0.5592017
## [2,] 2 0.8378062
## [3,] 3 0.6804462
## [4,] 4 0.8240519
## [5,] 5 0.6955286
## [6,] 6 0.8044649
## [7,] 7 0.8668329
## [8,] 8 0.9143899
## [9,] 9 0.9237524
## [10,] 10 0.9550390
## [11,] 11 0.9689032
## [12,] 12 0.9285905
## [13,] 13 0.9537861
## [14,] 14 0.9712432
## [15,] 15 0.9778869
## [16,] 16 0.9861581
## [17,] 17 0.9462377
## [18,] 18 0.9632284
## [19,] 19 0.9664228
## [20,] 20 0.9727337EXERCICE 2
On récupère le jeu de données INTC (en specifiant d’abord le répertoire où R doit chercher) :
# setwd("/Users/adrien/R/Data")
Dat0 <- read.table("/Users/adrien/R/Data/INTC", header=T)
INTC <- ts(data = Dat0, start=c(1973,1),end=c(2003,12), frequency=12)
plot(INTC)
On étudie la série INTC :
acf(INTC,lag=100)
pacf(INTC,lag=100)
Box.test.2(INTC,1:50)## Retard p-value
## [1,] 1 0.89386607
## [2,] 2 0.91427190
## [3,] 3 0.37222884
## [4,] 4 0.32435131
## [5,] 5 0.43077263
## [6,] 6 0.48746844
## [7,] 7 0.08921509
## [8,] 8 0.06034335
## [9,] 9 0.08975699
## [10,] 10 0.11540033
## [11,] 11 0.10113237
## [12,] 12 0.11191248
## [13,] 13 0.11873731
## [14,] 14 0.03322213
## [15,] 15 0.04178934
## [16,] 16 0.05883222
## [17,] 17 0.08055149
## [18,] 18 0.10015606
## [19,] 19 0.11179859
## [20,] 20 0.13825313
## [21,] 21 0.11193471
## [22,] 22 0.12569008
## [23,] 23 0.10902431
## [24,] 24 0.13738412
## [25,] 25 0.13363336
## [26,] 26 0.13808178
## [27,] 27 0.16901174
## [28,] 28 0.19914106
## [29,] 29 0.19749051
## [30,] 30 0.23010734
## [31,] 31 0.26453123
## [32,] 32 0.30529418
## [33,] 33 0.31421034
## [34,] 34 0.28614257
## [35,] 35 0.32761383
## [36,] 36 0.36553451
## [37,] 37 0.40743047
## [38,] 38 0.44907218
## [39,] 39 0.48802732
## [40,] 40 0.51625613
## [41,] 41 0.55752500
## [42,] 42 0.59868650
## [43,] 43 0.62994721
## [44,] 44 0.66410356
## [45,] 45 0.66202762
## [46,] 46 0.67220471
## [47,] 47 0.70436039
## [48,] 48 0.73974665
## [49,] 49 0.77192588
## [50,] 50 0.79402710D’après le test de Box-Pierce, l’hypothèse d’un bruit blanc ne peut pas être rejetée.
On étudie maintenant cette série élevée au carré pour débusquer un effet ARCH :
acf(INTC^2,lag=100) 
pacf(INTC^2,lag=100)
Box.test(INTC^2,3)##
## Box-Pierce test
##
## data: INTC^2
## X-squared = 34.233, df = 3, p-value = 1.769e-07En effet, si une série est un bruit blanc, alors la série élevée au carré doit l’être aussi, et ce n’est ici manifestement pas la cas en vertu du test de Box-Pierce. On tente de modéliser INTC^2 avec un ARMA :
ARMA <- auto.arima(INTC^2,stationary=TRUE, seasonal=FALSE)
ARMA## Series: INTC^2
## ARIMA(1,0,2) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 mean
## 0.8746 -0.8027 0.0725 0.0181
## s.e. 0.0758 0.0903 0.0610 0.0041
##
## sigma^2 estimated as 0.001384: log likelihood=698.47
## AIC=-1386.95 AICc=-1386.78 BIC=-1367.35Box.test.2(residuals(ARMA),1:50)## Retard p-value
## [1,] 1 0.9974349
## [2,] 2 0.9999866
## [3,] 3 0.9302030
## [4,] 4 0.9617276
## [5,] 5 0.9784405
## [6,] 6 0.9841967
## [7,] 7 0.9807108
## [8,] 8 0.9886544
## [9,] 9 0.9950860
## [10,] 10 0.9929692
## [11,] 11 0.9685373
## [12,] 12 0.4814540
## [13,] 13 0.5566356
## [14,] 14 0.3151310
## [15,] 15 0.3592240
## [16,] 16 0.3676372
## [17,] 17 0.3887653
## [18,] 18 0.4180108
## [19,] 19 0.4644522
## [20,] 20 0.5245124
## [21,] 21 0.4692977
## [22,] 22 0.5114473
## [23,] 23 0.5075649
## [24,] 24 0.5663583
## [25,] 25 0.6230293
## [26,] 26 0.6686375
## [27,] 27 0.6876068
## [28,] 28 0.7207481
## [29,] 29 0.7621122
## [30,] 30 0.7980296
## [31,] 31 0.8231999
## [32,] 32 0.8522068
## [33,] 33 0.8764502
## [34,] 34 0.9005882
## [35,] 35 0.9211340
## [36,] 36 0.9322147
## [37,] 37 0.9470294
## [38,] 38 0.9574408
## [39,] 39 0.9667697
## [40,] 40 0.9743623
## [41,] 41 0.9782832
## [42,] 42 0.9830756
## [43,] 43 0.9797898
## [44,] 44 0.9848187
## [45,] 45 0.9883363
## [46,] 46 0.9913175
## [47,] 47 0.9936642
## [48,] 48 0.9954219
## [49,] 49 0.9966068
## [50,] 50 0.9971864Une modélisation ARMA(1,2) semble raisonnable. Cependant, si INTC était un GARCH(p,q) et, supposant que la série à un moment d’ordre 4, INTC^2 serait un ARMA(max(p,q),q), qui ne contient pas la classe ARMA(1,2). Cepedant, on voit que parmi les coeffcients :
coef(ARMA)## ar1 ma1 ma2 intercept
## 0.87457353 -0.80266059 0.07251611 0.01809983le coefficient ma2 est negligeable devant ar1 et ma1. On peut tenter de négliger ce coefficient et modéliser avec un ARMA(1,1) :
arima(INTC^2,order=c(1,0,1)) ##
## Call:
## arima(x = INTC^2, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.9122 -0.7971 0.0182
## s.e. 0.0469 0.0685 0.0044
##
## sigma^2 estimated as 0.001374: log likelihood = 697.75, aic = -1387.5Box.test.2(residuals(arima(INTC^2,order=c(1,0,1))),1:50)## Retard p-value
## [1,] 1 0.4062199
## [2,] 2 0.5891072
## [3,] 3 0.5198242
## [4,] 4 0.6824185
## [5,] 5 0.8014312
## [6,] 6 0.8586763
## [7,] 7 0.8770328
## [8,] 8 0.9094476
## [9,] 9 0.9472803
## [10,] 10 0.9523248
## [11,] 11 0.8707045
## [12,] 12 0.3483931
## [13,] 13 0.4247163
## [14,] 14 0.2560471
## [15,] 15 0.3053455
## [16,] 16 0.3159037
## [17,] 17 0.3461926
## [18,] 18 0.3710159
## [19,] 19 0.4182775
## [20,] 20 0.4659493
## [21,] 21 0.4255396
## [22,] 22 0.4608356
## [23,] 23 0.4579136
## [24,] 24 0.5166969
## [25,] 25 0.5744912
## [26,] 26 0.6148837
## [27,] 27 0.6383801
## [28,] 28 0.6674724
## [29,] 29 0.7124502
## [30,] 30 0.7521049
## [31,] 31 0.7777328
## [32,] 32 0.8116283
## [33,] 33 0.8380984
## [34,] 34 0.8670140
## [35,] 35 0.8924016
## [36,] 36 0.9058884
## [37,] 37 0.9250510
## [38,] 38 0.9391823
## [39,] 39 0.9514647
## [40,] 40 0.9621827
## [41,] 41 0.9686885
## [42,] 42 0.9752068
## [43,] 43 0.9707715
## [44,] 44 0.9777078
## [45,] 45 0.9824633
## [46,] 46 0.9867298
## [47,] 47 0.9901243
## [48,] 48 0.9927310
## [49,] 49 0.9945922
## [50,] 50 0.9955474Cette modélisation est pertinente.
On entrevoit ainsi la possibilité que INTC provienne d’un GARCH(1,1).
L’analyse suivante nous montre que c’est une hypothèse raisonnable :
Fit <- garch(INTC,order=c(1,1)) ##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.609728e-02 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -5.619e+02
## 1 6 -5.619e+02 6.38e-05 1.21e-04 3.2e-03 6.0e+05 3.4e-04 3.64e+01
## 2 9 -5.620e+02 1.53e-04 1.58e-04 2.6e-02 2.1e+00 2.7e-03 3.66e+00
## 3 13 -5.654e+02 5.99e-03 9.72e-03 4.5e-01 2.0e+00 8.9e-02 3.43e+00
## 4 14 -5.664e+02 1.84e-03 4.64e-03 2.9e-01 2.0e+00 8.9e-02 2.74e-01
## 5 15 -5.686e+02 3.79e-03 4.17e-03 2.3e-01 2.0e+00 8.9e-02 2.74e-01
## 6 17 -5.711e+02 4.47e-03 6.83e-03 3.8e-01 1.9e+00 3.0e-01 1.15e-01
## 7 22 -5.713e+02 2.79e-04 6.59e-04 3.3e-04 7.1e+00 3.5e-04 9.62e-03
## 8 23 -5.713e+02 1.46e-05 1.28e-05 3.2e-04 2.0e+00 3.5e-04 8.19e-03
## 9 29 -5.721e+02 1.42e-03 1.12e-03 7.8e-02 9.1e-01 9.0e-02 8.60e-03
## 10 30 -5.733e+02 2.11e-03 3.31e-03 1.3e-01 2.0e+00 1.8e-01 1.89e-01
## 11 31 -5.738e+02 8.27e-04 6.10e-03 7.3e-02 0.0e+00 1.3e-01 6.10e-03
## 12 33 -5.746e+02 1.44e-03 5.87e-03 2.5e-02 8.9e-01 5.7e-02 1.06e-02
## 13 35 -5.752e+02 1.02e-03 2.38e-03 7.0e-03 1.4e+00 1.1e-02 3.18e-03
## 14 36 -5.756e+02 6.79e-04 6.81e-04 6.9e-03 1.8e+00 1.1e-02 1.70e-03
## 15 38 -5.760e+02 6.37e-04 9.14e-04 2.7e-02 9.1e-01 4.5e-02 1.93e-03
## 16 40 -5.760e+02 1.24e-05 3.49e-05 2.3e-03 6.7e-01 3.8e-03 4.86e-05
## 17 41 -5.760e+02 3.64e-06 9.07e-06 1.1e-03 0.0e+00 2.0e-03 9.07e-06
## 18 43 -5.760e+02 1.33e-07 3.40e-07 1.8e-04 8.2e-01 4.0e-04 4.59e-07
## 19 45 -5.760e+02 4.05e-07 8.55e-08 1.8e-04 0.0e+00 4.4e-04 8.55e-08
## 20 46 -5.760e+02 4.99e-08 8.92e-10 6.6e-06 0.0e+00 1.6e-05 8.92e-10
## 21 47 -5.760e+02 1.10e-08 1.59e-11 8.8e-07 0.0e+00 1.6e-06 1.59e-11
## 22 48 -5.760e+02 1.74e-09 1.32e-13 3.0e-07 0.0e+00 5.6e-07 1.32e-13
## 23 49 -5.760e+02 4.43e-11 1.07e-16 8.6e-09 0.0e+00 1.5e-08 1.07e-16
## 24 50 -5.760e+02 4.35e-13 3.29e-20 1.2e-10 0.0e+00 2.1e-10 3.29e-20
## 25 51 -5.760e+02 -3.36e-15 1.25e-24 6.5e-13 0.0e+00 1.4e-12 1.25e-24
##
## ***** X- AND RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -5.759765e+02 RELDX 6.479e-13
## FUNC. EVALS 51 GRAD. EVALS 25
## PRELDF 1.254e-24 NPRELDF 1.254e-24
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.158365e-03 1.000e+00 -1.605e-07
## 2 8.739503e-02 1.000e+00 -2.417e-09
## 3 8.456534e-01 1.000e+00 -2.211e-09coef(Fit)## a0 a1 b1
## 0.001158365 0.087395034 0.845653382Er <- residuals(Fit)
plot.ts(Er) ; Box.test.2(Er,1:100) 
## Retard p-value
## [1,] 1 0.7192460
## [2,] 2 0.8464215
## [3,] 3 0.5421388
## [4,] 4 0.4955916
## [5,] 5 0.6374841
## [6,] 6 0.6352119
## [7,] 7 0.3365955
## [8,] 8 0.3140770
## [9,] 9 0.4056276
## [10,] 10 0.4827762
## [11,] 11 0.5419189
## [12,] 12 0.5709801
## [13,] 13 0.5042275
## [14,] 14 0.3256200
## [15,] 15 0.3803232
## [16,] 16 0.4488597
## [17,] 17 0.5138924
## [18,] 18 0.5221752
## [19,] 19 0.5772624
## [20,] 20 0.6405308
## [21,] 21 0.5930960
## [22,] 22 0.5013889
## [23,] 23 0.4377090
## [24,] 24 0.4960664
## [25,] 25 0.4425148
## [26,] 26 0.4799296
## [27,] 27 0.5353108
## [28,] 28 0.5851709
## [29,] 29 0.5898103
## [30,] 30 0.6360464
## [31,] 31 0.6813499
## [32,] 32 0.7261372
## [33,] 33 0.7213385
## [34,] 34 0.6941055
## [35,] 35 0.7335554
## [36,] 36 0.7676821
## [37,] 37 0.7962666
## [38,] 38 0.8285838
## [39,] 39 0.8537805
## [40,] 40 0.8669364
## [41,] 41 0.8868791
## [42,] 42 0.9070094
## [43,] 43 0.9218327
## [44,] 44 0.9368688
## [45,] 45 0.9247388
## [46,] 46 0.9251788
## [47,] 47 0.9356199
## [48,] 48 0.9478785
## [49,] 49 0.9580580
## [50,] 50 0.9653318
## [51,] 51 0.9682192
## [52,] 52 0.9746951
## [53,] 53 0.9718667
## [54,] 54 0.9701964
## [55,] 55 0.9751967
## [56,] 56 0.9791780
## [57,] 57 0.9821518
## [58,] 58 0.9820760
## [59,] 59 0.9859155
## [60,] 60 0.9855846
## [61,] 61 0.9876990
## [62,] 62 0.9757557
## [63,] 63 0.9796232
## [64,] 64 0.9821014
## [65,] 65 0.9817073
## [66,] 66 0.9849550
## [67,] 67 0.9768166
## [68,] 68 0.9792963
## [69,] 69 0.9826463
## [70,] 70 0.9764412
## [71,] 71 0.9799386
## [72,] 72 0.9821583
## [73,] 73 0.9857106
## [74,] 74 0.9880569
## [75,] 75 0.9904210
## [76,] 76 0.9902713
## [77,] 77 0.9922968
## [78,] 78 0.9917486
## [79,] 79 0.9932373
## [80,] 80 0.9894478
## [81,] 81 0.9872160
## [82,] 82 0.9884680
## [83,] 83 0.9882724
## [84,] 84 0.9744134
## [85,] 85 0.9787206
## [86,] 86 0.9804616
## [87,] 87 0.9791726
## [88,] 88 0.9811699
## [89,] 89 0.9839301
## [90,] 90 0.9860001
## [91,] 91 0.9885858
## [92,] 92 0.9905897
## [93,] 93 0.9923090
## [94,] 94 0.9932638
## [95,] 95 0.9937527
## [96,] 96 0.9948069
## [97,] 97 0.9953887
## [98,] 98 0.9962485
## [99,] 99 0.9966755
## [100,] 100 0.9972576On peut alors tenter de prédire la volatilité \(\sigma_t^2\) de la série INTC \(X_t=\sigma_t\varepsilon_t\) à l’aide des formules du cours. En effet, on a pour un GARCH(1,1) :
\[ X_t=\sigma_t\varepsilon_t,\qquad \sigma_t^2=\omega+ \alpha X_{t-1}^2+\beta\sigma_{t-1}^2. \]
On nommes les constantes :
omega = coef(Fit)[1]
alpha = coef(Fit)[2]
beta = coef(Fit)[3]
n = length(INTC) ; n## [1] 372On commence par estimer \(\sigma_n^2\) par récurrence (en supposant que \(\sigma_0=0\)) :
sigma2 <- vector()
sigma2[1]=0
for (i in 1:n){sigma2[i+1]= omega +alpha*INTC[i]^2 + beta*sigma2[i]}
sigma2_n=sigma2[n+1]Ensuite, on crée une fonction \(\verb'sigma2_est'\) de \(\ell\) qui nous renvoie \(\hat\sigma_n^2(\ell)\) comme suite :
sigma2_est <- function(l){
sig_est1 <- unname(omega + alpha*INTC[n]^2 + beta*sigma2_n)
if (l==1){return(sig_est1)}
else {
Est <- vector()
Est[1]=sig_est1
for (i in 2:l){Est[i]=omega*(1-(alpha+beta)^l)/(1-alpha-beta)+(alpha+beta)^l*Est[i-1]}
return(Est[l])
}
}Par exemple, on prédit à l’horizon \(\ell=1,2,3,4\) :
sigma2_est(1) ; sigma2_est(2) ; sigma2_est(3) ; sigma2_est(4)## [1] 0.01503668## [1] 0.0153298## [1] 0.01580714## [1] 0.0163155