r - R中的多元GARCH-M

Question

我想知道是否有可以在 R 中实现多元 GARCH-M 模型的 R 包。我知道有一些包可以处理多元 GARCH 模型（例如 BEKK、DCC、CCC），但我还没有找到方法估计一个 MGARCH-M 模型。特别是，我对具有 VAR 结构并包括波动率项和滞后项的平均方程感兴趣，如 Elder&Serletis (2010) 中所建模的

任何意见，将不胜感激。

谢谢你，拉斐尔

Answer 1

这是来自Ruey Tsay的 R 脚本

要使用：

garchM(data, type= ). data是您的 y 或矩阵（将仅使用 y 值）。type有三个选项，type = 1用于均值方差、均值type = 2波动率和均值type = 3对数（方差）。

"garchM" <- function(rtn,type=1){
# Estimation of a Gaussian GARCH(1,1)-M model.
##### The program uses GARCH(1,1) results as initial values.
# rtn: return series 
# type = 1 for Variance-in-mean
#      = 2 for volatility-in-mean
#      = 3 for log(variance)-in-mean
#
if(is.matrix(rtn))rtn=c(rtn[,1])
garchMdata <<- rtn
# obtain initial estimates
m1=garch11FIT(garchMdata)
est=as.numeric(m1$par); v1=m1$ht  ## v1 is sigma.t-square
Mean=est[1]; cc=est[2]; ar=est[3]; ma=est[4]; S=1e-6
if(type==2)v1=sqrt(v1)
if(type==3)v1=log(v1)
#### Obtain initial estimate of the parameters for the mean equation
m2=lm(rtn~v1)
Cnst=as.numeric(m2$coefficients[1])
gam=as.numeric(m2$coefficients[2])
params=c(mu=Cnst,gamma=gam, omega=cc, alpha=ar,beta=ma)
lowBounds=c(mu=-5*abs(Mean),gamma=-20*abs(gam), omega=S, alpha=S, beta=ma*0.6)
uppBounds=c(mu=5*abs(Mean),gamma=100*abs(gam), omega=cc*5 ,alpha=3*ar,beta=1-S)
### Pass model information via defining global variable
Vtmp <<- c(type,v1[1])
#
fit=nlminb(start = params, objective= glkM, lower=lowBounds, upper=uppBounds)
##,control=list(trace=3,rel.tol=1e-5))
epsilon = 0.0001 * fit$par
npar=length(params)
Hessian = matrix(0, ncol = npar, nrow = npar)
for (i in 1:npar) {
for (j in 1:npar) {
x1 = x2 = x3 = x4  = fit$par
x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (glkM(x1)-glkM(x2)-glkM(x3)+glkM(x4))/
(4*epsilon[i]*epsilon[j])
}
}
cat("Maximized log-likehood: ",-glkM(fit$par),"\n")
# Step 6: Create and Print Summary Report:
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)

m3=ResiVol(fit$par)

garchM <- list(residuals=m3$residuals,sigma.t=m3$sigma.t)
}

glkM = function(pars){
rtn <- garchMdata
mu=pars[1]; gamma=pars[2]; omega=pars[3]; alpha=pars[4]; beta=pars[5]
type=Vtmp[1]
nT=length(rtn)
# use conditional variance
if(type==1){
ht=Vtmp[2]
et=rtn[1]-mu-gamma*ht
at=c(et)
for (i in 2:nT){
sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
ept = rtn[i]-mu-gamma*sig2t
at=c(at,ept)
ht=c(ht,sig2t)
}
}
# use volatility
if(type==2){
ht=Vtmp[2]^2
et=rtn[1]-mu-gamma*Vtmp[2]
at=c(et)
for (i in 2:nT){
sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
ept=rtn[i]-mu-gamma*sqrt(sig2t)
at=c(at,ept)
ht=c(ht,sig2t)
}
}
# use log(variance)
if(type==3){
ht=exp(Vtmp[2])
et=rtn[1]-mu-gamma*Vtmp[2]
at=c(et)
for (i in 2:nT){
sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
ept=rtn[i]-mu-gamma*log(abs(sig2t))
at=c(at,ept)
ht=c(ht,sig2t)
}
}
#
hh=sqrt(abs(ht))
glk=-sum(log(dnorm(x=at/hh)/hh))

glk
}


ResiVol = function(pars){
rtn <- garchMdata
mu=pars[1]; gamma=pars[2]; omega=pars[3]; alpha=pars[4]; beta=pars[5]
type=Vtmp[1]
nT=length(rtn)
# use conditional variance
if(type==1){
ht=Vtmp[2]
et=rtn[1]-mu-gamma*ht
at=c(et)
for (i in 2:nT){
sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
ept = rtn[i]-mu-gamma*sig2t
at=c(at,ept)
ht=c(ht,sig2t)
}
}
# use volatility
if(type==2){
ht=Vtmp[2]^2
et=rtn[1]-mu-gamma*Vtmp[2]
at=c(et)
for (i in 2:nT){
sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
ept=rtn[i]-mu-gamma*sqrt(sig2t)
at=c(at,ept)
ht=c(ht,sig2t)
}
}
# use log(variance)
if(type==3){
ht=exp(Vtmp[2])
et=rtn[1]-mu-gamma*Vtmp[2]
at=c(et)
for (i in 2:nT){
sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
ept=rtn[i]-mu-gamma*log(abs(sig2t))
at=c(at,ept)
ht=c(ht,sig2t)
}
}
#

ResiVol <- list(residuals=at,sigma.t=sqrt(ht))
}

garch11FIT = function(x){
# Step 1: Initialize Time Series Globally:
tx <<- x
# Step 2: Initialize Model Parameters and Bounds:
Mean = mean(tx); Var = var(tx); S = 1e-6
params = c(mu = Mean, omega = 0.1*Var, alpha = 0.1, beta = 0.8)
lowerBounds = c(mu = -10*abs(Mean), omega = S^2, alpha = S, beta = S)
upperBounds = c(mu = 10*abs(Mean), omega = 100*Var, alpha = 1-S, beta = 1-S)
# Step 3: Set Conditional Distribution Function:
garchDist = function(z, hh) { dnorm(x = z/hh)/hh }
# Step 4: Compose log-Likelihood Function:
garchLLH = function(parm) {
mu = parm[1]; omega = parm[2]; alpha = parm[3]; beta = parm[4]
z = tx-mu; Mean = mean(z^2)
# Use Filter Representation:
e = omega + alpha * c(Mean, z[-length(tx)]^2)
h = filter(e, beta, "r", init = Mean)
hh = sqrt(abs(h))
llh = -sum(log(garchDist(z, hh)))
llh }
#####print(garchLLH(params))
# Step 5: Estimate Parameters and Compute Numerically Hessian:
fit = nlminb(start = params, objective = garchLLH,
lower = lowerBounds, upper = upperBounds)
#
est=fit$par
# compute the sigma.t^2 series
z=tx-est[1]; Mean=mean(z^2)
e=est[2]+est[3]*c(Mean,z[-length(tx)]^2)
h=filter(e,est[4],"r",init=Mean)

garch11Fit <- list(par=est,ht=h)
}