"RMfit" <- function(rtn,estim=TRUE){
# Estimation of the RiskMetrics special IGARCH(1,1) model.
# rtn: return series 
#
Idata <<- rtn
Var = var(Idata); S = 1e-6
beta=0.96
if(estim){
params=c(beta = 0.9)
lowerBounds = c(beta= S)
upperBounds = c(beta = 1-S)
# Step 3: set conditional distribution function:
igarchDist = function(z,hh){dnorm(x = z/hh)/hh}
# Step 4: Compose log-likelihood function:
igarchLLH = function(parm){
beta=parm[1]
z = Idata; Meanz=mean(z^2)
e= (1-beta) * c(Meanz, z[-length(Idata)]^2)
h = filter(e, beta, "r", init=Meanz)
hh = sqrt(abs(h))
llh = -sum(log(igarchDist(z, hh)))
llh
}
###print(igarchLLH(params))
# Step 5: Estimate Parameters and Compute Numerically Hessian:
fit = nlminb(start = params, objective = igarchLLH,
lower = lowerBounds, upper = upperBounds) #, control = list(trace=3))
epsilon = 0.0001 * fit$par
npar=length(params)
Hessian = matrix(0, ncol = npar, nrow = npar)
for (i in 1:npar) {
for (j in 1:npar) {
x1 = x2 = x3 = x4  = fit$par
x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (igarchLLH(x1)-igarchLLH(x2)-igarchLLH(x3)+igarchLLH(x4))/
(4*epsilon[i]*epsilon[j])
}
}
# Step 6: Create and Print Summary Report:
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
beta=fit$par[1]
}
else{cat("Default beta = 0.96 is used.","\n")}
###
z=Idata; Mz = mean(z^2)
e= (1-beta)*c(Mz,z[-length(z)]^2)
h = filter(e,beta, "r", init=Mz)
vol = sqrt(abs(h))

# Step 7: compute VaR and ES based on the estimated results
T=length(h)
Sigma=beta*h[T]+(1-beta)*z[T]^2
Vpred=sqrt(Sigma)
tt=matrix(c(T,Vpred),1,2)
colnames(tt) <- c("Orig","Vpred")
cat("\n Volatility prediction:\n")
print(tt)
prob=c(0.95,0.99,0.999)
q1=qnorm(prob)
d1=dnorm(q1); d11=d1/(1-prob)
VaR=q1*Vpred
ES=d11*Vpred
tbl=cbind(prob,VaR,ES)
cat("\n Risk measure based on RiskMetrics:\n")
print(tbl)

RMfit <- list(par=beta,volatility = vol)
}