Fit a GARCH (1,1) - model dengan kovariat di R

Saya memiliki beberapa pengalaman dengan pemodelan deret waktu, dalam bentuk model ARIMA sederhana dan sebagainya. Sekarang saya memiliki beberapa data yang menunjukkan pengelompokan volatilitas, dan saya ingin mencoba memulai dengan memasang model GARCH (1,1) pada data.

Saya memiliki serangkaian data dan sejumlah variabel yang saya pikir mempengaruhinya. Jadi dalam istilah regresi dasar, sepertinya:

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

Tapi saya benar-benar bingung bagaimana menerapkan ini menjadi model GARCH (1,1)? Saya telah melihat rugarch-paket dan fGarch-paket dalam R, tapi saya belum bisa melakukan sesuatu yang berarti selain contoh yang dapat ditemukan di internet.

Berikut adalah contoh implementasi menggunakan rugarchpaket dan dengan beberapa data palsu. Fungsi ini ugarchfitmemungkinkan untuk dimasukkannya regressor eksternal dalam persamaan rata-rata (perhatikan penggunaan external.regressorsdi fit.specdalam kode di bawah).

Untuk memperbaiki notasi, modelnya adalah mana dan menunjukkan kovariat pada waktu , dan dengan asumsi / persyaratan "biasa" pada parameter dan proses inovasi .

$$\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$

Nilai parameter yang digunakan dalam contoh adalah sebagai berikut.

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

Gambar di bawah ini menunjukkan seri kovariat dan serta seri . The kode yang digunakan untuk menghasilkan mereka disediakan di bawah ini.$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))

Penyesuaian dilakukan dengan ugarchfitsebagai berikut.

## 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)

Estimasi parameter adalah

## 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

Dan nilai sebenarnya yang sesuai adalah

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

Gambar berikut menunjukkan estimasi parameter dengan interval kepercayaan 95%, dan nilai sebenarnya. The Rkode yang digunakan untuk menghasilkan itu disediakan di bawah.

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)