R/rgarch-numderiv.R

In rgarch: Flexible GARCH modelling in R

#################################################################################
##
##   R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
##   This file is part of the R package rgarch.
##
##   The R package rgarch is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package rgarch is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
#################################################################################

# Fitting procedure required functions:
.hessian2sided = function(f, x, ...)
{
	n = length(x)
	fx = f(x, ...)
	eps = .Machine$double.eps
	
	# Compute the stepsize (h)
	h = eps^(1/3)*apply(as.data.frame(x), 1,FUN = function(z) max(abs(z), 1e-2))
	xh = x+h
	h = xh-x
	if(length(h) == 1) ee = matrix(h, ncol = 1, nrow = 1) else ee = Matrix(diag(h), sparse = TRUE)
	
	# Compute forward and backward steps
	gp = vector(mode = "numeric", length = n)
	gp = apply(ee, 2, FUN = function(z) f(x+z, ...))
	gm = vector(mode="numeric",length=n)
	gm = apply(ee, 2, FUN = function(z) f(x-z, ...))
	H = h%*%t(h)
	Hm = H
	Hp = H
	# Compute "double" forward and backward steps
	for(i in 1:n){
		for(j in  i:n){
			Hp[i,j] = f(x+ee[,i]+ee[,j], ...)
			Hp[j,i] = Hp[i,j]
			Hm[i,j] = f(x-ee[,i]-ee[,j], ...)
			Hm[j,i] = Hm[i,j]
		}
	}
	#Compute the hessian
	for(i in 1:n){
		for(j in  i:n){
			H[i,j] = ( (Hp[i,j]-gp[i]-gp[j]+fx+fx-gm[i]-gm[j]+Hm[i,j]) /H[i,j] )/2
			H[j,i] = H[i,j]
		}
	}
	return(H)
}
# robust standard errors
.robustSE = function(fun, pars, hess, n, ...)
{
	stderrors = solve(hess)
	np = length(pars)	
	h = apply(as.data.frame(pars), 1, FUN = function(x) min(abs(x/2)+1e-4, max(x,1e-2))*eps^(1/3))
	hplus =  pars+h
	hminus = pars-h
	likelihoodsminus = likelihoodsplus = zeros(n, np)
	hparsminus = hparsplus = matrix(pars, ncol = np, nrow = np, byrow = T)
	diag(hparsplus) = hplus
	diag(hparsminus) = hminus
	likelihoodsplus =  apply(hparsplus, 1,  FUN = function(x) -fun(x, ...))
	likelihoodsminus = apply(hparsminus, 1, FUN = function(x) -fun(x, ...))
	scores = matrix(0, ncol = np, nrow = n)
	likpm = likelihoodsplus-likelihoodsminus
	scores = likpm/(2*repmat(t(h), n, 1))
	scores = apply(scores, 2, FUN = function(x) scale(x, center = TRUE, scale = FALSE))
	B = (t(scores)%*%scores)
	robustSE = stderrors%*%B%*%stderrors
	return(list(scores = scores, robustSE = robustSE))
}




# The following functions are based on on Kevin Sheppard's MFE toolbox
#---------------------------------------------------------------------
neweywestcv = function(data, nlag = NULL, center = TRUE)
{
	# Long-run covariance estimation using Newey-West (Bartlett) weights
	#  if nlag empty=NULL NLAG=min(floor(1.2*T^(1/3)),T)
	N = dim(as.matrix(data))[1]
	if(is.null(nlag)) nlag=min(floor(1.2*N^(1/3)),N)
	if(center) data = apply(data, 2, FUN = function(x) scale(x, center = TRUE, scale = FALSE))
	# weights
	bw = (nlag+1-(seq(0,nlag,by=1)))/(nlag+1)
	cv = 1/N * t(data)%*%data
	for(i in 1:nlag){
		gmi = 1/N * (t(data[(i+1):N,])%*%data[1:(N-i),])
		gpp = gmi + t(gmi)
		cv = cv + bw[i+1]*gpp
	}
	return(cv)
}

robustvcv = function(fun, pars, nlag = 0, hess, n, ...)
{
	k = length(pars)
	h = apply(as.data.frame(pars), 1, FUN = function(x) max(abs(x*eps^(1/3)), 1e-7))
	
	hplus =  pars+h
	hminus = pars-h
	hparsminus = hparsplus = matrix(pars, ncol = k, nrow = k, byrow = T)
	diag(hparsplus) = hplus
	diag(hparsminus) = hminus
	
	likelihoodsminus = likelihoodsplus = zeros(n, k)
	likelihoodsplus =  apply(hparsplus, 1,  FUN = function(x) fun(x, ...))
	likelihoodsminus = apply(hparsminus, 1, FUN = function(x) fun(x, ...))
	
	scores = zeros(n,k)
	likpm = likelihoodsplus-likelihoodsminus
	scores = likpm/(2*repmat(t(h), n, 1))
	A = hess/n
	hess = A
	Ainv = try( solve(A, tol = 1e-24), silent = TRUE )
	if( inherits(Ainv, "try-error")){
		info = 1
		vcv = NA
	} else{
		if(nlag>0){
			B = neweywestcv(scores, nlag)
			vcv = (Ainv%*%B%*%Ainv)/n
		} else{
			B = cov(scores)
			vcv = (Ainv%*%B%*%Ainv)/n
		}
		info = 0
	}
	return(list(vcv = vcv, scores = scores, info = info))
}