Hi, I am trying to write the Realized GARCH model with order (1,1)

The model can be describe bellow:

r_t = sqrt( h_t) * z_t

logh_t = w + b*logh_(t-1) + r*logx_(t-1)

logx_t = c + q*logh_t + t1*z_t +t2*(z_t ^2 -1) + u_t

and z follow N(0,1) , u follow N(0, sigma.u^2)

But I'm troubled with the simulation check for my code.

After I simulate data from the model and estimate the data,

I can't get precise estimation for my setting parameters.

This is my simulation code:
=======================
sim<-function(theta)
{
 omega = theta[1]
 bet = theta[2]
 gam = theta[3]
 xi = theta[4]
 phi = theta[5]
 tau1 = theta[6]
 tau2 = theta[7]
 sigma.u = theta[8]
 n = theta[9]
 n.warm = 500
 n = n+n.warm
 z = rnorm(n+1)
 u = rnorm((n+1),0,sigma.u)
 logh.pre = 1
 logx.pre = xi+phi*logh.pre+tau1*z[1]+tau2*(z[1]^2-1)+u[1]
 logh = c(logh.pre,rep(0,n))
 r = c(rep(0,(n+1)))
 logx = c(logx.pre,rep(0,n))
 for(i in 2:(n+1))
 {
  logh[i] = omega+bet*logh[(i-1)]+gam*logx[(i-1)]
  logx[i] = xi+phi*logh[i]+tau1*z[i]+tau2*(z[i]^2-1)+u[i]
  r[i] = sqrt(exp(logh[i]))*z[i]
 }
 gdata = rbind(r,exp(logx))
 ans = gdata[,-(1:(n.warm+1))]
 ans
}
===========================

This is my estimation code:
============================
 LLH<-function(data,theta)
 {
  r = data[1,]
  x = data[2,]
  n = dim(data)[2]
  logx = log(x)
  logh = rep(0,n)
  u = rep(0,n)
  garchDist = function(r, hh) { dnorm(x = r/hh)/hh }
  measureDist = function(u,sigma){dnorm(x = u/sigma)/sigma}
  omega = theta[1]
  bet = theta[2]
  gam = theta[3]
  xi = theta[4]
  phi = theta[5]
  tau1 = theta[6]
  tau2 = theta[7]
  sigma.u = theta[8]
  logh.pre = 0.1
  logx.pre = 0.1
  logh[1]<-omega+bet*logh.pre+gam*logx.pre

u[1]<-logx[1]-xi-phi*logh[1]-tau1*(r[1]/sqrt(exp(logh[1])))-tau2*((r[1]/sqrt(exp(logh[1])))^2-1)
  for(i in 2:n)
  {
   logh[i] = omega+bet*logh[(i-1)]+gam*logx[(i-1)]
   u[i] =
logx[i]-xi-phi*logh[i]-tau1*(r[i]/sqrt(exp(logh[i])))-tau2*((r[i]/sqrt(exp(logh[i])))^2-1)
  }
  h = exp(logh)
  hh = sqrt(h)
  llh = -sum(log(garchDist(r, hh)))-sum(log(measureDist(u,sigma.u)))
  llh
 }

fitting<-function(data)
{
 LLH1<-function(theta) {LLH(data,theta)}
 ini = c(0.1,0.6,0.3,0.07,1,-0.03,0.1,0.4)
 for(i in 1:10)
 {
  fit = optim(ini,LLH1,control=list(trace=2),hessian=T)
  ini = fit$par
 }
 hessian = fit$hessian
 std = sqrt(diag(solve(hessian)))
 A = rbind(fit$par,std)
 A
}
=====================

Does anyone can point out that where I am wrong?  Thank you very much!

	[[alternative HTML version deleted]]

