[R-SIG-Finance] Hessian in GARCH-Models (rugarch-package)
alexios
alexios at 4dscape.com
Thu Feb 9 13:28:37 CET 2012
Dear Linus,
If you want to calculate the robust covariance matrix of the parameters
you should look at the 'robustvcv' function in the 'rugarch-numderiv.R'
file.
It is based on White's (1982) method and descibed in the Vignette
(original coding is due to K.Sheppard and translated from matlab into
R). You will need to enable the returning of the "non-summed" vector of
likelihoods from your function (just use a switch as I do in my code).
To answer your other question, "x" is the vector of parameters at the
optimum, "f" the likelihood function (returning the scalar likelihood),
"..." any additional parameters required by your likelihood function.
HTH.
Regards,
Alexios
On 08/02/2012 22:30, Lin23 wrote:
> Hello,
>
> i wrote a code for a GJR-GARCH-Model. With my current data and model there
> is something wrong with the robust standard errors and the p-value. I use
> maxLik-package with BHHH-algorithm. The coefficients in my model are equal
> to those in the rugarch-package or EViews. But the p-value for the intercept
> in the variance-equation is way too high in my model. When i estimate my
> model with non-robust standard errors (-crossprod(fit$gradientObs))
> everything is all right. So I think there is something wrong with the
> hessian that is computed by the maxLik-package. I want to compute my own
> hessian. I used the code for the hessian from the rugarch-package
> (rugarch-numderiv.R). But I failed. Maybe someone can help me with my
> problem. I would appreciate it very much.
>
> My GARCH-code:
>
> garch<-function(y,Di,Mi,Do,Fr,x1,x2,x3,x4,Dum,backcast){
>
> n<-length(y)
>
> reg<- lm(y ~ Di + Mi + Do + Fr + x1 + x2 + x3 + x4+(Dum:x3)+(Dum:x4))
> #Start Values
>
> exp.smooth<- function(lambda){
> smooth<- as.vector(1:n-1)
> init.hh<- lambda^n*mean(reg$res^2)+(1-lambda)*sum(lambda^smooth*reg$res^2)
> # Backcast initial variance
> return(init.hh)
> }
> init.hh<- exp.smooth(lambda=backcast)
>
>
> LLH<- function(par){ # Log.Likelihood
> a<-par[1]
> di<- par[2]
> mi<- par[3]
> do<- par[4]
> fr<- par[5]
> b1<- par[6]
> b2<- par[7]
> b3<- par[8]
> b4<- par[9]
> b5<- par[10]
> b6<- par[11]
> omega<- par[12]
> alpha<- par[13]
> beta<- par[14]
> gjr<- par[15]
> dummy<- par[16]
>
> res<-array(n)
> hh<-array(n)
> ll<-numeric(n)
>
> for (i in 1:n){
>
> res[i]<-y[i]-a-di*Di[i]-mi*Mi[i]-do*Do[i]-fr*Fr[i]-b1*x1[i]-b2*x2[i]-b3*x3[i]-b4*x4[i]-b5*Dum[i]*x3[i]-b6*Dum[i]*x4[i]
> #Mean-Equation
> }
>
> hh[1]<- init.hh
> for(i in 2:n){
> hh[i]<-
> (omega+alpha*res[i-1]^2+beta*(hh[i-1]-omega)+ifelse(res[i-1]<0,
> gjr*res[i-1]^2, 0))*(1+dummy*Dum[i]) #Variance-Equation
> }
>
> for (i in 1:n){
> ll[i]<- -1/2*log(2*pi*hh[i]) - 1/2*res[i]^2/hh[i] #
> Log.Likelihood
> }
> ll
>
> }
>
> param<- c(reg$coef, omega = 0.1*mean(reg$res^2),alpha = 0.15, beta =
> 0.6,gjr=0.05,dummy=0)
> fit<-maxLik(start=param,
> logLik=LLH,method="BHHH",finalHessian=TRUE,iterlim=2500) # fit with
> maxLik-package
>
>
> grad<- fit$gradientObs
> robust.vcov<-
> solve(fit$hessian)%*%t(grad)%*%grad%*%solve(fit$hessian)
> se.coef<- sqrt(diag(robust.vcov)) # robust standard errors
> tval<- fit$estimate/se.coef
> matcoef<- cbind(fit$estimate, se.coef, tval, 2*(1-pnorm(abs(tval))))
> dimnames(matcoef)<- list(names(tval), c("Estimate", "Std. Error","t
> value", "Pr(>|t|)"))
>
>
> cat("\n---------------------------------------------------------------------\n")
> #Output
> cat("\nKoeffizienten\n")
>
> cat("\n---------------------------------------------------------------------\n")
> printCoefmat(matcoef, digits = 6, signif.stars=TRUE)
> }
> EST1<-
> garch(y=dat2$r_csi,Di=dat2$Di,Mi=dat2$Mi,Do=dat2$Do,Fr=dat2$Fr,x1=dat2$r_t,x2=dat2$r_sp_l,x3=dat2$r_csi_l,x4=dat2$r_csi_l_3,Dum=dat2$f,backcast=0.7)
>
> See below for the hessian-code from the rugarch-package. I dont know what I
> to enter for x and for ... in my case. I hope someone can help.
>
> hessian-code from rugarch-package:
>
> .hessian2sided = function(f, x, ...)
> {
> n = length(x)
> fx = f(x, ...)
> eps = .Machine$double.eps
>
> # Compute the stepsize (h)
> h = eps^(1/3)*apply(as.data.frame(x), 1,FUN = function(z) max(abs(z),
> 1e-2))
> xh = x+h
> h = xh-x
> if(length(h) == 1) ee = matrix(h, ncol = 1, nrow = 1) else ee =
> as.matrix(diag(h))
>
> # Compute forward and backward steps
> gp = vector(mode = "numeric", length = n)
> gp = apply(ee, 2, FUN = function(z) f(x+z, ...))
> gm = vector(mode="numeric",length=n)
> gm = apply(ee, 2, FUN = function(z) f(x-z, ...))
> H = h%*%t(h)
> Hm = H
> Hp = H
> # Compute "double" forward and backward steps
> for(i in 1:n){
> for(j in i:n){
> Hp[i,j] = f(x+ee[,i]+ee[,j], ...)
> Hp[j,i] = Hp[i,j]
> Hm[i,j] = f(x-ee[,i]-ee[,j], ...)
> Hm[j,i] = Hm[i,j]
> }
> }
> #Compute the hessian
> for(i in 1:n){
> for(j in i:n){
> H[i,j] = ( (Hp[i,j]-gp[i]-gp[j]+fx+fx-gm[i]-gm[j]+Hm[i,j]) /H[i,j] )/2
> H[j,i] = H[i,j]
> }
> }
> return(H)
> }
>
> Thank you all!!!
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Hessian-in-GARCH-Models-rugarch-package-tp4371133p4371133.html
> Sent from the Rmetrics mailing list archive at Nabble.com.
>
> _______________________________________________
> R-SIG-Finance at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
More information about the R-SIG-Finance
mailing list