首先得看下Garch(1,1)模型长什么样
(1)首先生成一个Garch(1,1)序列
n <- 10000a <- c(0.1, 0.3, 0.6) # GARCH(1,1) coefficientssigma_0 <- sqrt(a[1]/(1 - a[2] - a[3]))r_0 <- sigma_0 * rnorm(1)r <- double(n)sigma <- double(n)sigma[1] <- sqrt(a[1] + a[2] * r_0^2 + a[3] * sigma_0^2)r[1] <- sigma[1] * rnorm(1)for(i in 2:n) # Generate GARCH(1,1) process{ sigma[i] <- sqrt(a[1] + a[2]*r[i-1]^2 + a[3]*sigma[i-1]^2) r[i] <- sigma[i] * rnorm(1)}(2)用R的rugarch包来求下参数
library(rugarch)r.garch <- ugarchspec(mean.model = list(armaOrder = c(0,0),include.mean = FALSE),variance.model = list(garchOrder = c(1, 1))) # Fit garch(1,1) fit = ugarchfit(data = r, spec = r.garch) # Diagnostic tests(3)写出似然函数,用optim求得参数
loglik <- function(param){ a0 <- param[1]; a1 <- param[2] ; a2 <- param[3] sigma_0 <- sqrt(a[1]/(1 - a[2] - a[3])) r_0 <- sigma_0 * rnorm(1) sigma <- double(n) sigma[1] <- sqrt(a[1] + a[2] * r_0^2 + a[3] * sigma_0^2) d <- 0 for(i in 2:n) { sigma[i] <- sqrt(a[1] + a[2]*r[i-1]^2 + a[3]*sigma[i-1]^2) } for(i in 1:n){ d <- d + log(dnorm(r[i],0,sigma[i])) } return(-d)}loglik(c(0.1,0.1,0.3))optim(c(0.1,0.1,0.3), loglik)