R/rdcc-postestimation.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.
##
#################################################################################
.eps = .Machine$double.eps
	
.dcchessian = function(f, pars, garchpars, N, data, garchenv, fname)
{
	parallel = get("parallel", garchenv)
	parallel.control = get("parallel.control", garchenv)
	
	x = c(garchpars, pars)
	n = length(x)
	dccn = length(pars)
	garchn = length(garchpars)
	fx = f(pars, garchpars, data, returnType = "llh", garchenv = garchenv)
	.eps = .Machine$double.eps
	# Compute the stepsize (h)
	h = .eps ^ (1/3) * pmax( abs( x ), 1 )
	xh = x + h
	h  = xh - x
	ee = Matrix( diag( h ), sparse = TRUE )
	
	# Compute forward and backward steps
	g = vector(mode = "numeric", length = n)
	if( parallel ){
		if( parallel.control$pkg == "multicore" ){
			tmp = mclapply(1:n, FUN = function(i){
						cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
						tmpg =  f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i], 
								garchpars = x[1:garchn] + ee[1:garchn, i], 
					data = data, returnType = "llh", garchenv = garchenv)
					return( tmpg )
					}, mc.cores = parallel.control$cores)
			g = as.numeric( unlist(tmp) )
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("x", "garchn", "ee", "data", "garchenv", "n", "fname", local = TRUE)
			tmp = sfLapply(as.list(1:n), fun = function(i){
						cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
						tmpg =  eval(parse(text = paste("rgarch:::", fname, 
						"( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i], garchpars = x[1:garchn] + ee[1:garchn, i], 
								data = data, returnType = 'llh', garchenv = garchenv)", sep = "")))
						return( tmpg )
					})
			# sfStop() we'll stop it further down
			g = as.numeric( unlist(tmp) )
		}
	} else{
		tmp = lapply(as.list(1:n), FUN = function(i){
					cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
					tmpg =  f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i], 
							garchpars = x[1:garchn] + ee[1:garchn, i], 
							data = data, returnType = "llh", garchenv = garchenv)
					return( tmpg )
				})
		g = as.numeric( unlist(tmp) )
	}
	
	H = h %*% t( h )
	#idx = 1
	# Compute "double" forward and backward steps
	#cat("\n")
	#for(i in 1:n){
	#	for(j in  (n - N + 1):n){
	#		cat(paste("Evaluating Function ",idx," out of ",n*N,"\n",sep=""))
	#		if(i <= j){
	#			H[i, j] = (f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j], 
	#						garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j], 
	#						data = data, returnType = "llh", garchenv = garchenv) - g[i] - g[j] + fx) / H[i, j]
	#			H[j, i] = H[i, j]
	#		}
	#		idx = idx + 1
	#	}
	#}
	
	if( parallel ){
		if( parallel.control$pkg == "multicore" ){
			tmp = mclapply(1:n, FUN = function(i){
					Htmp = H
					for(j in  (n - N + 1):n){
						if(i <= j){
							Htmp[i, j] = (f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j], 
												garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j], 
												data = data, returnType = "llh", garchenv = garchenv) - g[i] - g[j] + fx) / Htmp[i, j]
							Htmp[j, i] = Htmp[i, j]
						}
					}
					return(Htmp)
				})
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("x", "H", "garchn", "ee", "data", "garchenv", "n", "N", "g", "fx", "fname", local = TRUE)
			tmp = sfLapply(as.list(1:n), fun = function(i){
						Htmp = H
						for(j in  (n - N + 1):n){
							if(i <= j){
								Htmp[i, j] = eval(parse(text = paste("(rgarch:::", fname, "( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j], 
													garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j], 
													data = data, returnType = 'llh', garchenv = garchenv) - g[i] - g[j] + fx) / Htmp[i, j]", sep = "")))
								Htmp[j, i] = Htmp[i, j]
							}
						}
						return(Htmp)
					})
			sfStop()
		}
		for(i in 1:n){
			for(j in  (n - N + 1):n){
				if(i <= j){
					H[i, j] = tmp[[i]][i, j]
					H[j, i] = tmp[[i]][j, i]
				}
			}
		}
	} else{
		tmp = lapply(as.list(1:n), FUN = function(i){
					Htmp = H
					for(j in  (n - N + 1):n){
						if(i <= j){
							Htmp[i, j] = (f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j], 
												garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j], 
												data = data, returnType = "llh", garchenv = garchenv) - g[i] - g[j] + fx) / Htmp[i, j]
							Htmp[j, i] = Htmp[i, j]
						}
					}
					return(Htmp)
				})
		for(i in 1:n){
			for(j in  (n - N + 1):n){
				if(i <= j){
					H[i, j] = tmp[[i]][i, j]
					H[j, i] = H[i, j]
				}
			}
		}
	}
	
	newH = H[(n - N + 1):n, ]
	H = newH
	return(H)
}


# 2-stage partitioned standard errors for DCC type models

.dccmakefitmodel = function(garchmodel, f, pars, garchpars, data, garchenv, timer, message, fname)
{
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	mspec = get("mspec", garchenv)
	npar = get("npar", garchenv)
	
	parallel = get("parallel", garchenv)
	parallel.control = get("parallel.control", garchenv)
	
	eval.se = get("eval.se", garchenv)
	
	fitlist = get("fitlist", garchenv)
	paridx = get("paridx", garchenv)
	paridxa = get("paridxa", garchenv)
	dccnames = mspec$optimization.model$modelnames
	resids = residuals(fitlist)
	sigmas = sigma(fitlist)
	dccOrder = mspec$dccOrder
	dccP = dccOrder[1]
	dccQ = dccOrder[2]
	m = length(fitlist@fit)
	datanames = colnames(data)
	mo = mspec$optimization.model$maxOrder
	garchparameters = vector(mode = "numeric")
	garchNames = vector(mode = "character")
	
	for(i in 1:m){
		cf = coef(fitlist@fit[[i]])
		cfnames = names(cf)
		garchNames = c(garchNames, paste(datanames[i], cfnames, sep = ".") )
	}
	
	np = length(garchpars) + length(pars)
	ns = np - length(pars) + 1
	
	sol = f(pars, garchpars, data = data, returnType = "all", garchenv = garchenv)
	likelihoods 	= sol$lik
	loglikelihood 	= sol$llh
	Rtout = sol$Rt
	Qtout = sol$Qt
	
	N = dim(data)[1]
	Ht = array( 0, dim = c(m, m, N) )
	stdresid = matrix(0, nrow = N, ncol = m)
	
	if( parallel ){
		if( parallel.control$pkg == "multicore" ){
			tmp = mclapply(1:N, FUN = function(i){
						tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
						zz = eigen( tmph )
						sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
						tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
						return( list( H = tmph, Z = tmpz ) )
					})
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("sigmas", "Rtout", "resids", local = TRUE)
			tmp = sfLapply(as.list(1:N), fun = function(i){
						tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
						zz = eigen( tmph )
						sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
						tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
						return( list( H = tmph, Z = tmpz ) )
					})
			if(!eval.se) sfStop()
			# else we'll stop the cluster when the standard errors are evaluated.
		}
		for(i in 1:N){
			Ht[,,i] = tmp[[i]]$H
			stdresid[i,] = tmp[[i]]$Z
		}
	} else{
		tmp = lapply(as.list(1:N), FUN = function(i){
					tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
					zz = eigen( tmph )
					sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
					tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
					return( list( H = tmph, Z = tmpz ) )
				})
		for(i in 1:N){
			Ht[,,i] = tmp[[i]]$H
			stdresid[i,] = tmp[[i]]$Z
		}
	}
	
	assign("stdresid", stdresid, envir = garchenv)
	assign("Ht", Ht, envir = garchenv)
	
	if(eval.se){
		A = zeros( np, np )
		tidx = 1
		for(i in 1:m){
			cvar = fitlist@fit[[i]]@fit$cvar
			workingsize = dim(cvar)[1]
			A[(tidx:(tidx + workingsize - 1)), (tidx:(tidx + workingsize - 1))] = solve(cvar)
			tidx = tidx + workingsize
		}
		
		cat("\n\nCalculating Standard Errors, this can take a while\n")
		otherA = .dcchessian(f = f, pars = pars, garchpars = garchpars, N = length(pars), data = data, 
				garchenv = garchenv, fname)
		A[(np - length(pars) + 1):np, ] = otherA
		jointscores = zeros(N, np)
		tidx = 1
		for(i in 1:m){
			cf = coef(fitlist@fit[[i]])
			workingsize = length(cf)
			# head(fitlist@fit[[i]]@fit$scores, 22)
			jointscores[,(tidx:(tidx + workingsize - 1))] = fitlist@fit[[i]]@fit$scores
			tidx = tidx + workingsize
		}
		#Now all we need to do is calculate the scores form the dcc estimator and we have everything
		h = pmax( abs( pars/2 ), 1e-2 ) * .eps^(1/3)
		hplus = pars + h
		hminus = pars - h
		likelihoodsplus =  zeros( N, length(pars) )
		likelihoodsminus = zeros( N, length(pars) )
		for(i in 1:length(pars)){
			hparameters1 = pars
			hparameters2 = pars
			hparameters1[i] = hplus[i]
			hparameters2[i] = hminus[i]
			LHT1 = f(pars = hparameters1, garchpars, data = data, returnType = "lik", garchenv = garchenv)
			LHT2 = f(pars = hparameters2, garchpars, data = data, returnType = "lik", garchenv = garchenv)
			likelihoodsplus[, i]  = LHT1
			likelihoodsminus[, i] = LHT2
		}
		
		sctemp = likelihoodsplus - likelihoodsminus
		DCCscores = matrix(NA, ncol = dim(sctemp)[2], nrow = dim(sctemp)[1])
		sdtemp = 2 * repmat( t( h ), N, 1 )
		for(i in 1:dim(sctemp)[2]){
			DCCscores[,i] = sctemp[,i] / sdtemp[,i]
		}
		jointscores[, ns:np] = DCCscores
		B = cov( jointscores )
		A = A/ (N) 
		dcccvar = ( solve( A ) %*% B %*% solve( A ) ) / N
		
		allpars = c(garchpars, pars)
		allnames = c(garchNames, dccnames)
		se.coef = sqrt(diag(abs(dcccvar)))
		tval = as.numeric( allpars/se.coef )
		pval = 2* ( 1 - pnorm( abs( tval ) ) )
		matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
		matcoef[, 1] = allpars
		matcoef[, 2] = se.coef
		matcoef[, 3] = tval
		matcoef[, 4] = pval
		dimnames(matcoef) = list(allnames, c(" Estimate",
						" Std. Error", " t value", "Pr(>|t|)"))
		
	} else{
		allpars = c(garchpars, pars)
		allnames = c(garchNames, dccnames)
		se.coef = rep(NA, length(allpars))
		tval = rep(NA, length(allpars))
		pval = rep(NA, length(allpars))
		matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
		matcoef[, 1] = allpars
		dimnames(matcoef) = list(allnames, c(" Estimate",
						" Std. Error", " t value", "Pr(>|t|)"))
		dcccvar = NULL
		jointscores = NULL
	}
	
	dccfit = list()
	dccfit$coef = allpars
	names(dccfit$coef) = allnames
	dccfit$dcccoef = pars
	names(dccfit$dcccoef) = dccnames
	dccfit$garchcoef = garchpars
	names(dccfit$garchcoef) = garchNames
	dccfit$matcoef = matcoef
	dccfit$cvar = dcccvar
	dccfit$scores = jointscores
	dccfit$R = Rtout
	dccfit$H = Ht
	dccfit$Q = Qtout
	dccfit$llh = loglikelihood
	dccfit$log.likelihoods = likelihoods
	dccfit$garchmodel = garchmodel
	dccfit$timer = timer
	dccfit$convergence = 0
	dccfit$message = message
	return( dccfit )
}


.dccmakefiltermodel = function(garchmodel, f, pars, garchpars, data, garchenv, timer, message, fname)
{
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	mspec = get("mspec", garchenv)
	npar = get("npar", garchenv)
	
	parallel = get("parallel", garchenv)
	parallel.control = get("parallel.control", garchenv)
	eval.se = get("eval.se", garchenv)
	
	filterlist = get("filterlist", garchenv)
	paridx = get("paridx", garchenv)
	paridxa = get("paridxa", garchenv)
	dccnames = mspec$optimization.model$modelnames
	resids = residuals(filterlist)
	sigmas = sigma(filterlist)
	dccOrder = mspec$dccOrder
	dccP = dccOrder[1]
	dccQ = dccOrder[2]
	m = length(filterlist@filter)
	datanames = colnames(data)
	mo = mspec$optimization.model$maxOrder
	garchparameters = vector(mode = "numeric")
	garchNames = vector(mode = "character")
	
	for(i in 1:m){
		cf = coef(filterlist@filter[[i]])
		cfnames = names(cf)
		garchNames = c(garchNames, paste(datanames[i], cfnames, sep = ".") )
	}
	
	np = length(garchpars) + dccP + dccQ
	ns = np - dccP - dccQ + 1
	
	sol = f(pars, garchpars, data = data, returnType = "all", garchenv = garchenv)
	likelihoods 	= sol$lik
	loglikelihood 	= sol$llh
	Rtout = sol$Rt
	Qtout = sol$Qt
	
	N = dim(data)[1]
	Ht = array( 0, dim = c(m, m, N) )
	stdresid = matrix(0, nrow = N, ncol = m)
	
	
	if( parallel ){
		if( parallel.control$pkg == "multicore" ){
			tmp = mclapply(1:N, FUN = function(i){
						tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
						zz = eigen( tmph )
						sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
						tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
						return( list( H = tmph, Z = tmpz ) )
					})
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("sigmas", "Rtout", "resids", local = TRUE)
			tmp = sfLapply(as.list(1:N), fun = function(i){
						tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
						zz = eigen( tmph )
						sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
						tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
						return( list( H = tmph, Z = tmpz ) )
					})
			sfStop()
		}
		for(i in 1:N){
			Ht[,,i] = tmp[[i]]$H
			stdresid[i,] = tmp[[i]]$Z
		}
	} else{
		tmp = lapply(as.list(1:N), FUN = function(i){
					tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
					zz = eigen( tmph )
					sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
					tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
					return( list( H = tmph, Z = tmpz ) )
				})
		for(i in 1:N){
			Ht[,,i] = tmp[[i]]$H
			stdresid[i,] = tmp[[i]]$Z
		}
	}
	
	assign("stdresid", stdresid, envir = garchenv)
	assign("Ht", Ht, envir = garchenv)
	
	if( eval.se ){
		A = zeros( np, np )
		tidx = 1
		for(i in 1:m){
			cvar = filterlist@filter[[i]]@filter$cvar
			workingsize = dim(cvar)[1]
			A[(tidx:(tidx + workingsize - 1)), (tidx:(tidx + workingsize - 1))] = solve(cvar)
			tidx = tidx + workingsize
		}
		
		cat("\n\nCalculating Standard Errors, this can take a while\n")
		otherA = .dcchessian(f = f, pars = pars, garchpars = garchpars, N = dccP + dccQ, data = data, 
				garchenv = garchenv, fname)
		A[(np - dccP - dccQ + 1):np, ] = otherA
		jointscores = zeros(N, np)
		tidx = 1
		for(i in 1:m){
			cf = coef(filterlist@filter[[i]])
			workingsize = length(cf)
			jointscores[,(tidx:(tidx + workingsize - 1))] = filterlist@filter[[i]]@filter$scores
			tidx = tidx + workingsize
		}
		#Now all we need to do is calculate the scores form the dcc estimator and we have everything
		h = pmax( abs( pars/2 ), 1e-2 ) * .eps^(1/3)
		hplus = pars + h
		hminus = pars - h
		likelihoodsplus =  zeros( N, dccP + dccQ )
		likelihoodsminus = zeros( N, dccP + dccQ )
		for(i in 1:(dccP + dccQ)){
			hparameters1 = pars
			hparameters2 = pars
			hparameters1[i] = hplus[i]
			hparameters2[i] = hminus[i]
			LHT1 = f(pars = hparameters1, garchpars, data = data, returnType = "lik", garchenv = garchenv)
			LHT2 = f(pars = hparameters2, garchpars, data = data, returnType = "lik", garchenv = garchenv)
			likelihoodsplus[, i]  = LHT1
			likelihoodsminus[, i] = LHT2
		}
		
		sctemp = likelihoodsplus - likelihoodsminus
		DCCscores = matrix(NA, ncol = dim(sctemp)[2], nrow = dim(sctemp)[1])
		sdtemp = 2 * repmat( t( h ), N, 1 )
		for(i in 1:dim(sctemp)[2]){
			DCCscores[,i] = sctemp[,i] / sdtemp[,i]
		}
		jointscores[, ns:np] = DCCscores
		B = cov( jointscores )
		A = A/ N 
		dcccvar = ( solve( A ) %*% B %*% solve( A ) ) / N
		
		allpars = c(garchpars, pars)
		allnames = c(garchNames, dccnames)
		se.coef = sqrt(diag(abs(dcccvar)))
		tval = as.numeric( allpars/se.coef )
		pval = 2* ( 1 - pnorm( abs( tval ) ) )
		matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
		matcoef[, 1] = allpars
		matcoef[, 2] = se.coef
		matcoef[, 3] = tval
		matcoef[, 4] = pval
		dimnames(matcoef) = list(allnames, c(" Estimate",
						" Std. Error", " t value", "Pr(>|t|)"))
		
	} else{
		
		allpars = c(garchpars, pars)
		allnames = c(garchNames, dccnames)
		se.coef = rep(NA, length(allpars))
		tval = rep(NA, length(allpars))
		pval = rep(NA, length(allpars))
		matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
		matcoef[, 1] = allpars
		dimnames(matcoef) = list(allnames, c(" Estimate",
						" Std. Error", " t value", "Pr(>|t|)"))
		dcccvar = NULL
		jointscores = NULL
	}
	
	dccfilter = list()
	dccfilter$coef = allpars
	names(dccfilter$coef) = allnames
	dccfilter$dcccoef = pars
	names(dccfilter$dcccoef) = dccnames
	dccfilter$garchcoef = garchpars
	names(dccfilter$garchcoef) = garchNames
	dccfilter$matcoef = matcoef
	dccfilter$cvar = dcccvar
	dccfilter$scores = jointscores
	dccfilter$R = Rtout
	dccfilter$H = Ht
	dccfilter$Q = Qtout
	dccfilter$llh = loglikelihood
	dccfilter$log.likelihoods = likelihoods
	dccfilter$garchmodel = garchmodel
	dccfilter$timer = timer
	dccfilter$convergence = 0
	dccfilter$message = message
	return( dccfilter )
}


.sqrtsymmat = function( X )
{
	tmp = eigen( X )
	sqrttmp = ( tmp$vectors %*% diag( sqrt( tmp$values ) ) %*% solve( tmp$vectors ) )
	return( sqrttmp )
}

Try the rgarch package in your browser

Any scripts or data that you put into this service are public.

rgarch documentation built on May 2, 2019, 5:22 p.m.