Monter un GARCH (1,1) - modèle avec covariables en R

J'ai quelques expériences avec la modélisation de séries chronologiques, sous forme de modèles ARIMA simples et ainsi de suite. Maintenant, j'ai quelques données qui présentent des regroupements de volatilité, et je voudrais essayer de commencer par ajuster un modèle GARCH (1,1) sur les données.

J'ai une série de données et un certain nombre de variables, je pense, l'influencent. Donc, en termes de régression de base, cela ressemble à:

$$y_t = \alpha + \beta_1 x_{t1} + \beta_2 x_{t2} + \epsilon_t .$$

Mais je ne sais vraiment pas comment implémenter cela dans un modèle GARCH (1,1)? J'ai regardé le rugarch-package et le fGarch-package dedans R, mais je n'ai rien pu faire de plus significatif que les exemples que l'on peut trouver sur Internet.

Voici un exemple d'implémentation utilisant le rugarchpackage et avec de fausses données. La fonction ugarchfitpermet d'inclure des régresseurs externes dans l'équation moyenne (notez l'utilisation de external.regressorsin fit.specdans le code ci-dessous).

$$\begin{align*} y_t &= \lambda_0 + \lambda_1 x_{t,1} + \lambda_2 x_{t,2} + \epsilon_t, \\ \epsilon_t &= \sigma_t Z_t , \\ \sigma_t^2 &= \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma_{t-1}^2 , \end{align*}$$ $x_{t,1}$$x_{t,2}$$t$$Z_t$

Les valeurs des paramètres utilisées dans l'exemple sont les suivantes.

## Model parameters
nb.period  <- 1000
omega      <- 0.00001
alpha      <- 0.12
beta       <- 0.87
lambda     <- c(0.001, 0.4, 0.2)

$x_{t,1}$$x_{t,2}$$y_t$R

## Dependencies
library(rugarch)

## Generate some covariates
set.seed(234)
ext.reg.1 <- 0.01 * (sin(2*pi*(1:nb.period)/nb.period))/2 + rnorm(nb.period, 0, 0.0001)
ext.reg.2 <- 0.05 * (sin(6*pi*(1:nb.period)/nb.period))/2 + rnorm(nb.period, 0, 0.001)
ext.reg   <- cbind(ext.reg.1, ext.reg.2)

## Generate some GARCH innovations
sim.spec    <- ugarchspec(variance.model     = list(model = "sGARCH", garchOrder = c(1,1)), 
                          mean.model         = list(armaOrder = c(0,0), include.mean = FALSE),
                          distribution.model = "norm", 
                          fixed.pars         = list(omega = omega, alpha1 = alpha, beta1 = beta))
path.sgarch <- ugarchpath(sim.spec, n.sim = nb.period, n.start = 1)
epsilon     <- as.vector(fitted(path.sgarch))

## Create the time series
y <- lambda[1] + lambda[2] * ext.reg[, 1] + lambda[3] * ext.reg[, 2] + epsilon

## Data visualization
par(mfrow = c(3,1))
plot(ext.reg[, 1], type = "l", xlab = "Time", ylab = "Covariate 1")
plot(ext.reg[, 2], type = "l", xlab = "Time", ylab = "Covariate 2")
plot(y, type = "h", xlab = "Time")
par(mfrow = c(1,1))

Un ajustement est effectué ugarchfitcomme suit.

## Fit
fit.spec <- ugarchspec(variance.model     = list(model = "sGARCH",
                                                 garchOrder = c(1, 1)), 
                       mean.model         = list(armaOrder = c(0, 0),
                                                 include.mean = TRUE,
                                                 external.regressors = ext.reg), 
                       distribution.model = "norm")
fit      <- ugarchfit(data = y, spec = fit.spec)

Les estimations des paramètres sont

## Results review
fit.val     <- coef(fit)
fit.sd      <- diag(vcov(fit))
true.val    <- c(lambda, omega, alpha, beta)
fit.conf.lb <- fit.val + qnorm(0.025) * fit.sd
fit.conf.ub <- fit.val + qnorm(0.975) * fit.sd
> print(fit.val)
#     mu       mxreg1       mxreg2        omega       alpha1        beta1 
#1.724885e-03 3.942020e-01 7.342743e-02 1.451739e-05 1.022208e-01 8.769060e-01 
> print(fit.sd)
#[1] 4.635344e-07 3.255819e-02 1.504019e-03 1.195897e-10 8.312088e-04 3.375684e-04

Et les vraies valeurs correspondantes sont

> print(true.val)
#[1] 0.00100 0.40000 0.20000 0.00001 0.12000 0.87000

La figure suivante montre une estimation des paramètres avec des intervalles de confiance à 95% et les valeurs réelles. Le Rcode utilisé pour le générer est fourni ci-dessous.

plot(c(lambda, omega, alpha, beta), pch = 1, col = "red",
     ylim = range(c(fit.conf.lb, fit.conf.ub, true.val)),
     xlab = "", ylab = "", axes = FALSE)
box(); axis(1, at = 1:length(fit.val), labels = names(fit.val)); axis(2)
points(coef(fit), col = "blue", pch = 4)
for (i in 1:length(fit.val)) {
    lines(c(i,i), c(fit.conf.lb[i], fit.conf.ub[i]))
}
legend( "topleft", legend = c("true value", "estimate", "confidence interval"),
        col = c("red", "blue", 1), pch = c(1, 4, NA), lty = c(NA, NA, 1), inset = 0.01)