rugarch - question w.r.t. (robust) SE and bias in estimates of a particular MA(1)-eGARCH(1, 1)-STD model

Dear all,

I am currently experiencing troubles with a particular 
ARMA(0,1)-eGARCH(1,1)-STD (i.e. Student's t distribution) model that 
should be fitted to a series of business-daily log-returns of the DACX 
Performance index (data downloaded from bloomberg).
Since I have done similar analyses in the past without any problems I 
think that there might be a problem with the particular model and the 
particular data.

The sample comprises 390 observations and spans the period (1992-01-02 , 
1996-10-01) with the convention yyyy-mm-dd.
Estimation is done by means of the rugarch-package.
I have uploaded the data and a plot at 
http://www.file-upload.net/download-9194647/eGARCH.zip.html

Thanks a lot for any ideas or suggestions in advance!


---------------------------
Problems:

First of all I can't obtain robust standard errors, even when setting 
the option "fit.control=list(scale=1)" as recommended in related threads.
No parameter is at its limit, but persistence of the estimated model is 
very high (0.9892).

Secondly the t values and p-values according to the non-robust estimates 
seem to be too good to be true for such a small sample.

Conditional Variance Dynamics 	
-----------------------------------
GARCH Model	: eGARCH(1,1)
Mean Model	: ARFIMA(0,0,1)
Distribution	: std

Optimal Parameters
------------------------------------
         Estimate  Std. Error  t value Pr(>|t|)
mu      0.000481    0.000000  2915.43        0
ma1     0.084402    0.000724   116.58        0
omega  -0.107101    0.000064 -1670.67        0
alpha1 -0.090824    0.000059 -1549.21        0
beta1   0.989236    0.000231  4290.18        0
gamma1 -0.073881    0.000131  -565.25        0
shape   6.738614    0.010256   657.02        0


Robust Standard Errors:
         Estimate  Std. Error  t value Pr(>|t|)
mu      0.000481         NaN      NaN      NaN
ma1     0.084402         NaN      NaN      NaN
omega  -0.107101         NaN      NaN      NaN
alpha1 -0.090824         NaN      NaN      NaN
beta1   0.989236         NaN      NaN      NaN
gamma1 -0.073881         NaN      NaN      NaN
shape   6.738614         NaN      NaN      NaN


Thirdly  when investigating the small sample properties via the 
ugarchdistribution function and according to the plots as suggested on 
this website 
http://unstarched.net/2012/12/26/garch-parameter-uncertainty-and-data-size/
there are weird results for the parameter "gamma1". There seems to ba a 
bias that emerges as the sample size increases from 190 to 2990.
(Results schown here are based on m.sim=100. I have repeated the 
exercise a couple of times, also with m.sim=1000 but results did not 
change substantially.)

*------------------------------------*
*    GARCH Parameter Distribution    *
*------------------------------------*
Model : eGARCH
No. Paths (m.sim) : 100
Length of Paths (n.sim) : 190
Recursive : TRUE
Recursive Length : 2990
Recursive Window : 200

Coefficients: True vs Simulation Mean (Window-n)
                     mu      ma1    omega    alpha1   beta1     gamma1  shape
true-coef   0.00048106 0.084402 -0.10710 -0.090824 0.98924 -0.0738815 6.7386
window-190  0.00028241 0.075118 -0.81085 -0.126335 0.91744 -0.0953088 8.1002
window-390  0.00021380 0.083615 -0.21620 -0.112777 0.97757 -0.0385150 6.0713
window-590  0.00030893 0.083906 -0.23531 -0.103314 0.97612 -0.0117405 6.3599
window-790  0.00033639 0.081546 -0.12463 -0.094743 0.98720 -0.0093634 6.1996
window-990  0.00035573 0.082678 -0.12063 -0.093973 0.98767  0.0021929 6.1567
window-1190 0.00037253 0.081575 -0.11090 -0.091895 0.98865  0.0010467 5.9618
window-1390 0.00040460 0.080915 -0.11374 -0.093049 0.98840  0.0044305 5.8192
window-1590 0.00039122 0.080130 -0.11148 -0.090192 0.98862  0.0039467 5.8280
window-1790 0.00041750 0.083254 -0.10783 -0.090696 0.98902  0.0049624 5.6801
window-1990 0.00042600 0.084380 -0.10752 -0.090184 0.98906  0.0077372 5.7210
window-2190 0.00042967 0.082424 -0.10863 -0.090388 0.98894  0.0113595 5.6812
window-2390 0.00043955 0.082570 -0.10982 -0.091506 0.98884  0.0093621 5.7262
window-2590 0.00043285 0.083678 -0.10853 -0.091099 0.98896  0.0097338 5.7423
window-2790 0.00044481 0.082849 -0.10450 -0.089724 0.98938  0.0114610 5.6936
window-2990 0.00045002 0.081756 -0.10647 -0.090193 0.98919  0.0119148 5.7327

window-390 no. of non-converged fits: 3

window-990 no. of non-converged fits: 1


One can see that the mean of the gamma1 parameter estimates is actually 
moving away from the true value as the sample size increases - which is 
very counter-intuitive.
The pattern is even more clear in the uploaded plot.
Results can be replicated in R so it is not due to a particular Monte 
Carlo run.
Maybe estimation is very biased and inconsistent here.
But it would be the first time that I encounter such big problems.

I have tried to reproduce the results using STATA (version 11.2) but the 
model did not converge (message: flat log likelihood encountered, cannot 
find uphill direction).

R and rugarch are up to date. The R syntax is given as:


############################################################
#
#
require('rugarch')
require('parallel')

setwd("D:/MA")
mydata = read.table("mydata.txt" , header=T)

spec1 = ugarchspec(
      variance.model=list( model="eGARCH", garchOrder=c(1,1), 
submodel=NA ),
      mean.model=list( armaOrder = c(0,1) ), distribution.model="std" )

fit1 = ugarchfit( spec=spec1, data=mydata , solver="solnp" , 
fit.control=list(scale=1) )
fit1
persistence(fit1)

cl = makePSOCKcluster(4)       # please adjust the number of cores to be 
used for parallel processing
mod = ugarchdistribution( fitORspec=fit1 , n.sim=190, n.start=250, 
m.sim=100, recursive=TRUE,
      recursive.length=2990, recursive.window=200, rseed=654911, 
solver="solnp", cluster=cl )
stopCluster(cl)
mod

parvals = 
t(as.data.frame(mod at model$pars)$Level[as.data.frame(mod at model$pars)$Include==1])
colnames(parvals) = 
rownames(mod at model$pars)[as.data.frame(mod at model$pars)$Include==1]
rownames(parvals) = ""
parvals                                  # names and values for the 
parameter estimates

mu = sapply(mod at dist, FUN = function(x) x$simcoef[, 1])
mu$details = NULL
ma1 = sapply(mod at dist, FUN = function(x) x$simcoef[, 2])
ma1$details = NULL
omega = sapply(mod at dist, FUN = function(x) x$simcoef[, 3])
omega$details = NULL
alpha1 = sapply(mod at dist, FUN = function(x) x$simcoef[, 4])
alpha1$details = NULL
beta1 = sapply(mod at dist, FUN = function(x) x$simcoef[, 5])
beta1$details = NULL
gamma1 = sapply(mod at dist, FUN = function(x) x$simcoef[, 6])
gamma1$details = NULL
shape = sapply(mod at dist, FUN = function(x) x$simcoef[, 7])
shape$details = NULL

n = length(mod at dist) - 1
clr <- rainbow(n, alpha = 1, start = 0.55, end = 0.8)

par(mfrow = c(2, 4))
boxplot(na.omit(mu), names =  (c(1:n)*200-10) , col = clr)
abline(h = parvals[1], col = 2)
title('Mu')
boxplot(na.omit(ma1), names = (c(1:n)*200-10), col = clr)
abline(h = parvals[2], col = 2)
title('MA 1')
boxplot(na.omit(omega), names = (c(1:n)*200-10), col = clr)
abline(h = parvals[3], col = 2)
title('Omega')
boxplot(na.omit(alpha1), names = (c(1:n)*200-10), col = clr)
abline(h = parvals[4], col = 2)
title('Alpha 1')
boxplot(na.omit(beta1), names = (c(1:n)*200-10), col = clr)
abline(h = parvals[5], col = 2)
title('Beta 1')
boxplot(na.omit(gamma1), names = (c(1:n)*200-10), col = clr)
abline(h = parvals[6], col = 2)
title('Gamma 1')
boxplot(na.omit(shape), names = (c(1:n)*200-10), col = clr)
abline(h = parvals[7], col = 2)
title('Shape')
#
#
################################################################




	[[alternative HTML version deleted]]