efficient code for nonlinear garch model

Hi,
I'm trying to think of efficient code for NGARCH model,
sigma[t]^2=omega + alpha*(R[t-1]-theta*sigma[t-1])^2 + beta*(sigma[t-1)]^2
from book (Elements of Financial Risk Menagment-Peter F Christoffersen page
77).

I already coded it the "usual way" but it's taking too long (around 30 s).
Now I'm trying somehow to implement -filter- function like in paper -
Parameter Estimation of ARMA Models with
GARCH/APARCH Errors An R and SPlus Software Implementation - page 6.
The problem I can't overcome is that in this model I can't figure out how
to present the middle part ( 2 * Alpha * R[t-1]*theta*sigma[t-1]) with
filter function or is it even possible to code it this way.

Any help would appreciated.Here is my slower code:

"Ngarch" <- function(rtn,value){
   # Estimation of a non-symmertic GARCH, NGARCH(1,1), model.
   # Assume normal innovations
   # rtn: return series
   #
   # The likelihood function "mxlk" can be modified to fit more general
NGARCH
   #  models.
   # obtain initial estimates
   rtn=as.matrix(rtn)
   mu=mean(rtn)
   par=c(0.0001,0.8,0.01,0.7)
   #
   mxlk <- function(par){
     mxlk=0
     ht=var(rtn)
     T=length(rtn)
     if(T > 40)ht=var(rtn[1:40])
     at=rtn-mu
     for (i in 2:T){

sig2t=par[1]+par[2]*ht[i-1]+par[3]*ht[i-1]*(at[i-1]/sqrt(ht[i-1])-par[4])^2
       ht[i]=sig2t
       mxlk=mxlk+0.5*(log(sig2t) + (at[i])^2/ht[i])
     }
     mxlk
   }

   low=c(1e-6,1e-6,1e-6,1e-6)
   upp=c(100*var(rtn),1-1e-6,1-1e-6,5)
   mm=optim(par,mxlk,method="Nelder-Mead",hessian=T)
   #mm=nlminb(par,mxlk,hessian=T,scale=1,lower=low,upper=upp)
   #mm=optimx(par,mxlk,method="L-BFGS-B",hessian=T,lower=low,upper=upp)
   ## Print the results
   par=mm$par
   H=mm$hessian
   Hi = solve(H)
   cat(" ","\n")
   cat("Estimation results of NGARCH(1,1) model:","\n")
   cat("estimates: ",par,"\n")
   se=sqrt(diag(Hi))
   cat("std.errors: ",se,"\n")
   tra=par/se
   cat("t-ratio: ",tra,"\n")
   # compute the volatility series and residuals

   ht=var(rtn)
   T=length(rtn)
   if(T > 40)ht=var(rtn[1:40])
   at=rtn-mu
   for (i in 2:T){
     sig2t=par[1]+par[2]*ht[i-1]+par[3]*(at[i-1]-par[4]*sqrt(ht[i-1]))^2
     ht=c(ht,sig2t)
   }
   sigma.t=sqrt(ht)
   Ngarch <- as.matrix(cbind(as.numeric(at),as.numeric(sigma.t)))
   colnames(Ngarch)=c("residuals","volatility")
   Ngarch
}




I think that Alexios has this function in his package (called NAGARCH) but
I need to code this myself for my master work.
I already implemented (read ported) his ugarch function in mathematica
wolfram and it's working great.

	[[alternative HTML version deleted]]