想用R软件实现NGARCH（1，1）模型的估计，不知道R有什么程序包可以实现，看了rugarch包中的NGARCH，仔仔细细浏览了一遍rugarch包的Manual，发现此NGARCH非我所要估计的NGARCH，我所要估计的具体模型如下图，想请哪位高手不吝指教。此外，本人也尝试了用maxLik函数自编MLE程序来实现，开始我先从简单的GARCH（1，1）（也就是theta=0的情况）来尝试编程的效果，因为GARCH（1，1）在R中可以使用rugarch包实现，但是发现总是跟用rugarch包得到的结果相差很大，想了好久都不知道为何？也想请教哪位高手指正编程中的错误和不当之处，万分感谢！！！    待估模型、rugarch包实现的GARCH（1，1）以及自编MLE程序如下：
rugarch包实现的GARCH（1，1）估计结果
*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics         
-----------------------------------
GARCH Model        : sGARCH(1,1)
Mean Model        : ARFIMA(0,0,0)
Distribution        : norm

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.019888    0.021256  0.93562  0.34947
omega   0.025482    0.027440  0.92864  0.35308
alpha1  0.030516    0.020506  1.48817  0.13671
beta1   0.896523    0.092633  9.67825  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.019888    0.018322  1.08545  0.27772
omega   0.025482    0.046832  0.54412  0.58636
alpha1  0.030516    0.035093  0.86957  0.38454
beta1   0.896523    0.160244  5.59473  0.00000

LogLikelihood : -672.6365


本人自编的GARCH（1，1）MLE估计程序：

library(maxLik)  #GARCH(1,1)的MLE估计

T<-length(r1)

loglik<-function(param)

      {

         mu<-param[1]

         omega<-param[2]

         alpha<-param[3]         

         beta<-param[4]

         theta<-0

         r10<-mu

         sigma_sq0<-omega/(1-alpha*(1+theta^2)-beta)         

         sigma_sq<-rep(0,T)

         for(i in 1:T)

             {

                     if(i==1)

                     sigma_sq<-omega+alpha*(r10-mu-theta*sqrt(sigma_sq0))^2+beta*sigma_sq0

                     else

                     sigma_sq<-omega+alpha*(r1[i-1]-mu-theta*sqrt(sigma_sq[i-1]))^2+beta*sigma_sq[i-1]                  

             }

         z<-rep(0,T)

         for(j in 1:T)

             {

                     z[j]<-(r1[j]-mu)/sqrt(sigma_sq[j])

             }

         ll<-T*(-1/2*log(2*pi))-sum(z^2/2)

         return(ll)           

       }

res<-maxLik(loglik,start=c(0.2,0.05,0.4,0.5))

summary(res)