"Igarch" <- function(rtn,include.mean=F,volcnt=F){ # Estimation of a Gaussian IGARCH(1,1) model. # rtn: return series # include.mean: flag for the constant in the mean equation. # volcnt: flag for the constant term of the volatility equation. #### default is the RiskMetrics model # Idata <<- rtn Flag <<- c(include.mean,volcnt) # Mean=mean(Idata); Var = var(Idata); S = 1e-6 if((volcnt)&&(include.mean)){ params=c(mu = Mean,omega=0.1*Var,beta=0.85) lowerBounds = c(mu = -10*abs(Mean), omega= S^2, beta= S) upperBounds = c(mu = 10*abs(Mean), omega = 100*Var, beta = 1-S) } if((volcnt)&&(!include.mean)){ params=c(omega=0.1*Var, beta=0.85) lowerBounds=c(omega=S^2,beta=S) upperBounds=c(omega=100*Var,beta=1-S) } # if((!volcnt)&&(include.mean)){ params=c(mu = Mean, beta= 0.8) lowerBounds = c(mu = -10*abs(Mean), beta= S) upperBounds = c(mu = 10*abs(Mean), beta = 1-S) } if((!volcnt)&&(!include.mean)){ params=c(beta=0.85) 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){ include.mean=Flag[1] volcnt=Flag[2] mu=0; omega = 0 if((include.mean)&&(volcnt)){ my=parm[1]; omega=parm[2]; beta=parm[3]} if((!include.mean)&&(volcnt)){ omega=parm[1];beta=parm[2]} if((!include.mean)&&(!volcnt))beta=parm[1] if((include.mean)&&(!volcnt)){mu=parm[1]; beta=parm[2]} # z = (Idata - mu); Meanz = mean(z^2) e= omega + (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 } # Step 5: Estimate Parameters and Compute Numerically Hessian: fit = nlminb(start = params, objective = igarchLLH, lower = lowerBounds, upper = upperBounds) ##lower = lowerBounds, upper = upperBounds, control = list(trace=3)) epsilon = 0.0001 * fit\$par cat("Estimates: ",fit\$par,"\n") 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]) } } cat("Maximized log-likehood: ",igarchLLH(fit\$par),"\n") # 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) if((include.mean)&&(volcnt)){ mu=fit\$par[1]; omega=fit\$par[2]; beta = fit\$par[3] } if((include.mean)&&(!volcnt)){ mu = fit\$par[1]; beta = fit\$par[2]; omega = 0 } if((!include.mean)&&(volcnt)){ mu=0; omega=fit\$par[1]; beta=fit\$par[2] } if((!include.mean)&&(!volcnt)){ mu=0; omega=0; beta=fit\$par[1] } z=Idata-mu; Mz = mean(z^2) e= omega + (1-beta)*c(Mz,z[-length(z)]^2) h = filter(e,beta,"r",init=Mz) vol = sqrt(abs(h)) Igarch <- list(par=fit\$par,volatility = vol) }