R/rgarch-misc.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.
##
#################################################################################
# CONSTANTS used in package:
eps<-.Machine$double.eps
TinY = 1.0e-8

##################################################################################
# Helper Functions
##################################################################################

.makemodel = function(garchenv)
{
	model = vector(mode="list")
	mmodel = get("mmodel", garchenv)
	vmodel = get("vmodel", garchenv)
	omodel = get("omodel", garchenv)
	dmodel = get("dmodel", garchenv)
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	model$fixed = fixed
	model$fixid = fixid
	model$model = vmodel$model
	model$submodel = vmodel$submodel
	model$garchOrder = vmodel$garchOrder
	model$include.mex = mmodel$include.mex
	model$mexdata = mmodel$origmexdata
	model$vexdata = vmodel$origvexdata
	model$mxn = mmodel$mxn
	model$vxn = vmodel$vxn
	model$include.vex = vmodel$include.vex
	model$armaOrder = mmodel$armaOrder
	model$garchInMean = mmodel$garchInMean
	model$arfima = mmodel$arfima
	model$include.mex = mmodel$include.mex
	model$inMeanType = mmodel$inMeanType
	model$im = mmodel$im
	model$include.mean = mmodel$include.mean
	model$include.omega = dmodel$include.omega
	model$distribution = dmodel$distribution
	model$include.dlambda = dmodel$include.dlambda
	model$include.skew = dmodel$include.skew
	model$include.shape = dmodel$include.shape
	model$modelnames = omodel$modelnames
	model$start.pars = omodel$start.pars
	model$fixed.pars = omodel$fixed.pars
	model$distribution = dmodel$distribution
	model$modelmatrix = omodel$modelmatrix
	model$maxOrder = omodel$maxOrder
	model$maxOrder2 = omodel$maxOrder2
	model$pos = omodel$pos.matrix
	return(model)
}
.makearfimamodel = function(garchenv)
{
	model = vector(mode="list")
	mmodel = get("mmodel", garchenv)
	omodel = get("omodel", garchenv)
	dmodel = get("dmodel", garchenv)
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	model$fixed = fixed
	model$fixid = fixid
	model$include.mex = mmodel$include.mex
	model$mexdata = mmodel$origmexdata
	model$mxn = mmodel$mxn
	model$armaOrder = mmodel$armaOrder
	model$arfima = mmodel$arfima
	model$include.mex = mmodel$include.mex
	model$include.mean = mmodel$include.mean
	model$distribution = dmodel$distribution
	model$include.dlambda = dmodel$include.dlambda
	model$include.skew = dmodel$include.skew
	model$include.shape = dmodel$include.shape
	model$modelnames = omodel$modelnames
	model$start.pars = omodel$start.pars
	model$fixed.pars = omodel$fixed.pars
	model$distribution = dmodel$distribution
	model$modelmatrix = omodel$modelmatrix
	model$maxOrder = omodel$maxOrder
	model$pos = omodel$pos.matrix
	return(model)
}

.makearfimafitmodel = function(f, opt, data, T, m, timer, convergence, message, hess, garchenv)
{
	fit = vector(mode = "list")
	model = vector(mode="list")
	mmodel = get("mmodel", garchenv)
	omodel = get("omodel", garchenv)
	dmodel = get("dmodel", garchenv)
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	npar = get("npar", garchenv)
	assign(".llh",1, envir = garchenv)
	
	if(is.null(hess)){
		#fit$hessian = .hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)
		fit$hessian = hessian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "llh", garchenv = garchenv)		
	} else{
		fit$hessian = hess
	}
	fit$cvar = try(solve(fit$hessian), silent = TRUE)
	# (might also fail in which case user will see an error about the inversion failure)
	if(inherits(fit$cvar, "try-error")){
		zz = try(solve(.hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)), silent=TRUE)
		if(inherits(zz, "try-error")) {
			fit$cvar = NULL
			cat("\nrgarch-->warning: failed to invert hessian\n")
		} else{
			fit$cvar = zz
		}
	}
	#assign(".llh", 1, envir = garchenv)
	#		garchenv)
	# call LLH once more to get optimal parameters not report with returnType "all"
	temp = f(pars = opt, data = data, returnType="all", garchenv)
	fit$z = temp$z
	fit$data = data
	fit$dates = get("dates", garchenv)
	fit$LLH = -temp$llh
	fit$log.likelihoods = temp$LHT
	fit$residuals = temp$epsx
	assign(".llh",1, envir = garchenv)
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		nfix = length(fixed)
		pall[fixid] = fixed
		pall[-fixid] = opt
		fit$coef = pall
		Names = omodel$modelnames
		names(fit$coef) = Names
		complNULL = rep(NA, nfix)
		if(is.null(fit$cvar)){
			fit$se.coef = rep(NA,length(opt))
			fit$tval = rep(NA,length(opt))
			# change here
			fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
					fit$tval, rep(NA,length(opt)))
			fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
					complNULL, complNULL)
			fit$robust.se.coef = rep(NA,length(opt))
			fit$robust.tval = rep(NA,length(opt))
			# change here
			fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef, 
					fit$robust.tval, rep(NA,length(opt)))
			fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL, 
					complNULL, complNULL)
			fit$hessian.message = "failed to invert hessian"
		} else{
			tmp = robustvcv(fun = f, pars = opt, nlag = 0, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
			fit$robust.cvar = tmp$vcv
			fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv) 
			fit$se.coef = sqrt(diag(abs(fit$cvar)))
			fit$tval = fit$coef[-fixid]/fit$se.coef
			# change here
			fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
					fit$tval, 2*(1-pnorm(abs(fit$tval))))
			fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
					complNULL, complNULL)
			fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
			fit$robust.tval = fit$coef[-fixid]/fit$robust.se.coef
			# change here
			fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef,
					fit$robust.tval, 2*(1-pnorm(abs(fit$robust.tval))))
			fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL, 
					complNULL, complNULL)
			fit$hessian.message = NULL
		}
	} else{
		fit$coef = opt
		names(fit$coef) = omodel$modelnames
		if(is.null(fit$cvar)){
			fit$se.coef = rep(NA,length(opt))
			fit$tval = rep(NA,length(opt))
			# change here
			fit$matcoef = cbind(fit$coef, fit$se.coef,
					fit$tval, rep(NA,length(opt)))
			fit$robust.se.coef = rep(NA,length(opt))
			fit$robust.tval = rep(NA,length(opt))
			# change here
			fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval, 
					rep(NA,length(opt)))
			fit$hessian.message = "failed to invert hessian"
		} else{
			nlag=min(floor(1.2*(T)^(1/3)),(T))
			tmp = robustvcv(fun = f, pars = opt, nlag = nlag, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
			fit$robust.cvar = tmp$vcv
			fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv) 
			fit$se.coef = sqrt(diag(abs(fit$cvar)))
			fit$tval = fit$coef/fit$se.coef
			# change here
			fit$matcoef = cbind(fit$coef, fit$se.coef,
					fit$tval, 2*(1-pnorm(abs(fit$tval))))
			fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
			fit$robust.tval = fit$coef/fit$robust.se.coef
			# change here
			fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval, 
					2*(1-pnorm(abs(fit$robust.tval))))
			fit$hessian.message = NULL
		}
	}
	dimnames(fit$matcoef) = list(names(fit$coef), c(" Estimate",
					" Std. Error", " t value", "Pr(>|t|)"))
	dimnames(fit$robust.matcoef) = list(names(fit$coef), c(" Estimate",
					" Std. Error", " t value", "Pr(>|t|)"))
	fit$fitted.values = data - fit$residuals
	fit$convergence = convergence
	fit$message = message
	fit$timer = timer
	return(fit)
}


.makefitmodel = function(garchmodel, f, opt, data, T, m, timer, convergence, message, hess,
		garchenv)
{
	fit = vector(mode = "list")
	model = vector(mode="list")
	mmodel = get("mmodel", garchenv)
	vmodel = get("vmodel", garchenv)
	omodel = get("omodel", garchenv)
	dmodel = get("dmodel", garchenv)
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	npar = get("npar", garchenv)
	assign(".llh",1, envir = garchenv)
	
	if(is.null(hess)){
		#fit$hessian = .hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)
		fit$hessian = hessian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "llh", garchenv = garchenv)
		
	} else{
		fit$hessian = hess
	}
	fit$cvar = try(solve(fit$hessian), silent = TRUE)
	# (might also fail in which case user will see an error about the inversion failure)
	if(inherits(fit$cvar, "try-error")){
		zz = try(solve(.hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)), silent=TRUE)
		if(inherits(zz, "try-error")) {
			fit$cvar = NULL
			cat("\nrgarch-->warning: failed to invert hessian\n")
		} else{
			fit$cvar = zz
		}
	}
	#assign(".llh", 1, envir = garchenv)
	#		garchenv)
	# call LLH once more to get optimal parameters not report with returnType "all"
	temp = f(pars = opt, data = data, returnType="all", garchenv)
	if(garchmodel == "fGARCH" || garchmodel == "apARCH"){
		# apARCH and fGARCH are calculated in sigma rather than sigma^d (since d is a variable
		# to be calculated, unlike sGARCH where it is fixed at 2)
		fit$var = temp$h^2
		fit$sigma = temp$h
	} else{
		fit$var = abs(temp$h)
		fit$sigma = sqrt(abs(temp$h))
	}
	fit$z = temp$z
	fit$data = data
	fit$dates = get("dates", garchenv)
	fit$LLH = -temp$llh
	fit$log.likelihoods = temp$LHT
	fit$residuals = temp$epsx
	assign(".llh",1, envir = garchenv)
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		nfix = length(fixed)
		pall[fixid] = fixed
		pall[-fixid] = opt
		fit$coef = pall
		Names = omodel$modelnames
		names(fit$coef) = Names
		complNULL = rep(NA, nfix)
		if(is.null(fit$cvar)){
			fit$se.coef = rep(NA,length(opt))
			fit$tval = rep(NA,length(opt))
			# change here
			fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
					fit$tval, rep(NA,length(opt)))
			fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
					complNULL, complNULL)
			fit$robust.se.coef = rep(NA,length(opt))
			fit$robust.tval = rep(NA,length(opt))
			# change here
			fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef, 
					fit$robust.tval, rep(NA,length(opt)))
			fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL, 
					complNULL, complNULL)
			fit$hessian.message = "failed to invert hessian"
		} else{
			tmp = robustvcv(fun = f, pars = opt, nlag = 0, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
			fit$robust.cvar = tmp$vcv
			fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv) 
			fit$se.coef = sqrt(diag(abs(fit$cvar)))
			fit$tval = fit$coef[-fixid]/fit$se.coef
			# change here
			fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
					fit$tval, 2*(1-pnorm(abs(fit$tval))))
			fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
					complNULL, complNULL)
			fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
			fit$robust.tval = fit$coef[-fixid]/fit$robust.se.coef
			# change here
			fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
			fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef,
					fit$robust.tval, 2*(1-pnorm(abs(fit$robust.tval))))
			fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL, 
					complNULL, complNULL)
			fit$hessian.message = NULL
		}
	} else{
		fit$coef = opt
		names(fit$coef) = omodel$modelnames
		if(is.null(fit$cvar)){
			fit$se.coef = rep(NA,length(opt))
			fit$tval = rep(NA,length(opt))
			# change here
			fit$matcoef = cbind(fit$coef, fit$se.coef,
					fit$tval, rep(NA,length(opt)))
			fit$robust.se.coef = rep(NA,length(opt))
			fit$robust.tval = rep(NA,length(opt))
			# change here
			fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval, 
					rep(NA,length(opt)))
			fit$hessian.message = "failed to invert hessian"
		} else{
			nlag=min(floor(1.2*(T)^(1/3)),(T))
			# tmp = .robustSE(fun = f, pars = opt, hess = NULL, n = T-m, data = data, returnType = "LHT", garchenv)
			tmp = robustvcv(fun = f, pars = opt, nlag = nlag, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
			fit$robust.cvar = tmp$vcv
			fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv) 
			fit$se.coef = sqrt(diag(abs(fit$cvar)))
			fit$tval = fit$coef/fit$se.coef
			# change here
			fit$matcoef = cbind(fit$coef, fit$se.coef,
					fit$tval, 2*(1-pnorm(abs(fit$tval))))
			fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
			fit$robust.tval = fit$coef/fit$robust.se.coef
			# change here
			fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval, 
					2*(1-pnorm(abs(fit$robust.tval))))
			fit$hessian.message = NULL
		}
	}
	
	dimnames(fit$matcoef) = list(names(fit$coef), c(" Estimate",
					" Std. Error", " t value", "Pr(>|t|)"))
	dimnames(fit$robust.matcoef) = list(names(fit$coef), c(" Estimate",
					" Std. Error", " t value", "Pr(>|t|)"))
	fit$fitted.values = data-fit$residuals
	fit$convergence = convergence
	fit$message = message
	fit$kappa = temp$kappa
	fit$persistence = temp$persistence
	fit$timer = timer
	return(fit)
}

.forcregressors = function(include.mex, mxn, mexdata, mregfor, include.vex, 
		vxn, vexdata, vregfor, pars, n.ahead, N, out.sample, n.roll)
{
	# N is the original length
	treq = n.ahead + n.roll
	if(include.mex){
		
		if(!is.null(mregfor)){
			nmex = dim(as.matrix(mregfor))[1]
			mmex = dim(as.matrix(mregfor))[2]
		} else{
			nmex = 0
			mmex = 0
		}
		
		if(!is.null(mregfor) && mmex != mxn)
		{
			cat("\nugarchforecast-->error: Column dimension of external mean forecast matrix is wrong.")
			cat(paste("\nModel has ", mxn, " external regressors but forecast matrix has", mmex, sep = ""))
			stop("\n...exiting\n")
		}
		
		if(!is.null(mregfor) && nmex < treq)
		{
			cat(paste("\nugarchforecast-->error: You requested ", treq ," actual forecasts (including 
the rolling periods) but external mean forecasts provided have only",
nmex, " rows",sep=""))
			cat(paste("\nA minimum of ", treq," rows are required", sep=""))
			stop("\n...exiting")
		}
		if(is.null(mregfor)){
			mxf = rbind(as.matrix(mexdata)[1:(N-out.sample), ,drop = FALSE], matrix(0, ncol = mxn, nrow = treq))
		} else {
			mxf = rbind(as.matrix(mexdata)[1:(N-out.sample), ,drop = FALSE], as.matrix(mregfor)[1:treq,,drop = FALSE])
		}
		mxreg = pars[paste("mxreg", 1:mxn, sep = "")]
	} else{
		mxreg = NULL
		mxf = NULL
	}
	if(include.vex){
		
		if(!is.null(vregfor)){
			nvex = dim(as.matrix(vregfor))[1]
			mvex = dim(as.matrix(vregfor))[2]
		} else{
			nvex = 0
			mvex = 0
		}
		if(!is.null(vregfor) && mvex != vxn)
		{
			cat("\nugarchforecast-->error: Column dimension of external variance forecast matrix is wrong.")
			cat(paste("\nModel has ",vxn," external regressors but forecast matrix has", mvex, sep = ""))
			stop("\n...exiting\n")
		}
		# N is the original length
		if(!is.null(vregfor) && nvex < treq)
		{
			cat(paste("\nugarchforecast-->error: You requested ", treq ," actual forecasts (including 
									the rolling periods) but external variance forecasts provided have only",
							nvex, " rows",sep=""))
			cat(paste("\nA minimum of ", treq," rows are required", sep=""))
			stop("\n...exiting")
		}
		
		if(is.null(vregfor)){
			vxf = rbind(as.matrix(vexdata)[1:(N-out.sample), ,drop = FALSE], matrix(0, ncol = vxn, nrow = treq))
		} else {
			vxf = rbind(as.matrix(vexdata)[1:(N-out.sample), ,drop = FALSE], as.matrix(vregfor)[1:treq,,drop = FALSE])
		}
		vxreg = pars[paste("vxreg", 1:vxn, sep = "")]
	} else{
		vxreg = NULL
		vxf = NULL
	}
	return(list(mxreg = mxreg, mxf = mxf, vxreg = vxreg, vxf = vxf))
}

.simregressors = function(include.mex, mxn, mexsimdata, mexdata, include.vex, vxn, 
		vexsimdata, vexdata, N, pars, n, m.sim, m)
{
	if(include.mex){
		if(is.null(mexsimdata)){
			mexsimdata = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = n)
		}
		if(!is.null(mexsimdata))
		{
			if(!is.list(mexsimdata)) 
				stop("\nugarchsim-->error: mexsimdata should be a list of length m.sim")
			if(length(mexsimdata) != m.sim){
				msd = vector(mode = "list", length = m.sim)
				for(i in 1:m.sim) msd[[i]] = as.matrix(mexsimdata[[1]])
				mexsimdata = msd
				cat("\nugarchsim-->warning: length of mexsimdata list not equal to m.sim...
replicating first list element m.sim times.\n")
			}
			for(i in 1:m.sim){
				if(dim(as.matrix(mexsimdata[[i]]))[2] != mxn ) 
					stop(paste("\nugarchsim-->error: mexsimdata ", i," has wrong no. of column", sep=""))
				if(dim(as.matrix(mexsimdata[[i]]))[1] != n )
					stop(paste("\nugarchsim-->error: mexsimdata ", i," has wrong no. of rows", sep=""))
			}		
		}
		
		mxreg = pars[paste("mxreg",1:mxn,sep="")]
		mexsimlist = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim){
			premexdata = mexdata[(N-m+1):N,,drop=FALSE]
			mexsimlist[[i]] = matrix(rbind(premexdata, mexsimdata[[i]]), ncol = mxn)
		}
	} else{
		mxreg = 0
		mexsimlist = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) mexsimlist[[i]]=0
	}
	if(include.vex){
		if(is.null(vexsimdata)){
			vexsimdata = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) vexsimdata[[i]] = matrix(0, ncol = vxn, nrow = n)
		}
		if(!is.null(vexsimdata))
		{
			if(!is.list(vexsimdata)) 
				stop("\nugarchsim-->error: mexsimdata should be a list of length m.sim")
			if(length(vexsimdata) != m.sim){
				vsd = vector(mode = "list", length = m.sim)
				for(i in 1:m.sim) vsd[[i]] = as.matrix(vexsimdata[[1]])
				vexsimdata = vsd
				cat("\nugarchsim-->warning: length of vexsimdata list not equal to m.sim...
replicating first list element m.sim times.\n")
			}
			for(i in 1:m.sim){
				if(dim(as.matrix(vexsimdata[[i]]))[2] != vxn ) 
					stop(paste("\nugarchsim-->error: vexsimdata ", i," has wrong no. of column", sep=""))
				if(dim(as.matrix(vexsimdata[[i]]))[1] != n )
					stop(paste("\nugarchsim-->error: vexsimdata ", i," has wrong no. of rows", sep=""))
			}		
		}
		
		vxreg = pars[paste("vxreg",1:vxn,sep="")]
		vexsimlist = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim){
			prevexdata = vexdata[(N-m+1):N,,drop=FALSE]
			vexsimlist[[i]] = matrix(rbind(prevexdata, vexsimdata[[i]]), ncol = vxn)
		}
	} else{
		vxreg = 0
		vexsimlist = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) vexsimlist[[i]]=0
	}
	return(list(mxreg = mxreg, mexsimlist = mexsimlist, vxreg = vxreg, 
					vexsimlist = vexsimlist))
}

.custzdist = function(custom.dist, zmatrix, m.sim, n)
{
	if(is.na(custom.dist$name) | is.na(custom.dist$name)){
		
		z = matrix(apply(zmatrix, 1, FUN = function(x) .makeSample(as.character(x[1]), lambda = as.numeric(x[2]), 
							skew = as.numeric(x[3]), shape = as.numeric(x[4]), n = as.numeric(x[5]), 
							seed = as.integer(x[6]))), n, m.sim)
		
		#nx = dim(zmatrix)[1]
		#tmp = apply(as.data.frame(1:nx), 1, FUN = function(i){
		#			set.seed(zmatrix[i,6]);
		#			.C("distributionsample", 
		#			shape = as.double(zmatrix[i,4]), skew = as.double(zmatrix[i,3]), 
		#			lambda = as.double(zmatrix[i,2]), 
		#			ndis = as.integer(zmatrix[i,1]), n = as.integer(zmatrix[i,5]), 
		#			rvec = double(zmatrix[i,5]), DUP = FALSE, PACKAGE ="rgarch")$rvec})
		#z = matrix(tmp, nrow = n, ncol = m.sim)
	}
	if(!is.na(custom.dist$name) && !is.na(custom.dist$distfit)){
		
		if(is.matrix(custom.dist$distfit))
		{
			if(dim(custom.dist$distfit)[2]!=m.sim) stop("column dimension of custom 
				innovations\n matrix must be equal to m.sim")
			if(dim(custom.dist$distfit)[1]!=n) stop("row dimension 
				of custom innovations\n matrix must be equal to n.sim+n.start")
			z = custom.dist$distfit
		} else{
			if(!is.character(custom.dist$name)) stop("custom distribution must be a 
				character string")
			temp = paste("r", custom.dist$name, sep="")
			if(is.null(custom.dist$distfit)) stop("custom distribution missing a 
				distfit object")
			# if this is often used we might consider using apply to
			# use a new seed for each m.sim
			.rdist = eval(parse(text=paste(temp)))
			tmp = .rdist(n*m.sim, custom.dist$distfit)
			z = matrix(as.numeric(tmp), ncol = m.sim, nrow = n, byrow=TRUE)
		}
	}
	return(z)
}
.stars = function(testvector,levels=c(0.01, 0.05, 0.1))
{
	N = length(testvector)
	ans = vector(mode="character",length=N)
	#recursive replacement
	z = which(testvector<levels[3])
	ans[z] = c("*")
	z = which(testvector<levels[2])
	ans[z] = c("**")
	z = which(testvector<levels[1])
	ans[z] = c("***")
	ans
}

# method from a spec and data object

.safefit = function(spec, data, out.sample, solver, fit.control, solver.control)
{
	ans = try(ugarchfit(spec = spec, data = data, out.sample = out.sample, fit.control = fit.control,
					solver = solver, solver.control = solver.control), silent = TRUE)
	if(inherits(ans, "try-error")) ans = NULL
	return(ans)
}

.safefitarfima = function(spec, data, out.sample, solver, fit.control, solver.control)
{
	ans = try(arfimafit(spec = spec, data = data, out.sample = out.sample, fit.control = fit.control,
					solver = solver, solver.control = solver.control), silent = TRUE)
	if(inherits(ans, "try-error")) ans = NULL
	return(ans)
}


.boxcoxtransform = function(x, lambda)
{
	if(lambda!=0) ret = (x^lambda - 1)/lambda else ret = log(x)
	return(ret)
}

repmat = function(a, n, m)
{
	kronecker(matrix(1, n, m), a)
}
size = function(x, n = NULL)
{
	x = as.matrix(x)
	if(missing(n)) sol = c(n = dim(x)[1], m = dim(x)[2]) else sol = dim(x)[n]
	return(sol)
}

zeros = function(n = 1, m = 1)
{
	if(missing(m)) m = n
	sol = matrix(0, nrow = n, ncol = m)
	return(sol)
}

ones = function(n = 1, m = 1)
{
	if(missing(m)) m = n
	sol = matrix(1, nrow = n, ncol = m)
	return(sol)
}

newlagmatrix = function(x,nlags,xc)
{
	nlags = nlags+1
	xt = size(x, 1);
	newX = rbind(x, zeros(nlags, 1))
	lagmatrix = repmat(newX, nlags, 1)
	lagmatrix = matrix(lagmatrix[1:(size(lagmatrix,1)-nlags)], nrow = (xt+nlags-1), ncol = nlags)
	lagmatrix = lagmatrix[nlags:xt,]
	y = lagmatrix[,1]
	x = lagmatrix[,2:nlags]
	if(xc == 1) x = cbind(ones(size(x,1), 1), x)
	return(data.frame(y = y, x = x))
}

.colorgradient = function(n = 50, low.col = 0.6, high.col=0.7, saturation = 0.8) {
	if (n < 2) stop("n must be greater than 2")
	n1 = n%/%2
	n2 = - n - n1
	c(hsv(low.col, saturation, seq(1, 0.5, length = n1)),
			hsv(high.col, saturation, seq(0.5, 1, length = n2)))
}

.simlayout = function(m)
{
	if(m == 1){
		nf = c(1, 1, 2, 2)
		nf = layout(matrix(nf, 2, 2, byrow = TRUE), respect = TRUE)
		middle.plot = 1
	}
	if(m == 2){
		nf = c(1, 1, 1, 1, 0, 2, 2, 0, 3, 3, 3, 3)
		nf = layout(matrix(nf, 3, 4, byrow = TRUE), respect = TRUE)
		middle.plot = 2
	}
	if(m == 3){
		nf = c(1, 0, 0, 2, 0, 3, 3, 0, 4, 0, 0, 5)
		nf = layout(matrix(nf, 3, 4, byrow = TRUE), respect = TRUE)
		middle.plot = 3
	}
	if(m == 12){
	nf = c(1, 2, 3, 4,
		   5, 6, 6, 7,
		   8, 6, 6, 9,
		  10, 11, 12, 13)		
	nf = layout(matrix(nf, 4, 4, byrow = TRUE), respect = TRUE)
	}
}

.sdigit = function(x){
	sid = as.numeric(strsplit(format(as.numeric(x), scientific=TRUE), "e")[[1]])[2]
	10^(-sid)
}

.embed<-function(data, k, by = 1, ascending = FALSE)
{
	# n = no of time points, k = number of columns
	# by = increment. normally =1 but if =b calc every b-th point
	# ascending If TRUE, points passed in ascending order else descending.
	# Note that embed(1:n,k) corresponds to embedX(n,k,by=1,rev=TRUE)
	# e.g. embedX(10,3)
	#if(is.null(dim(data)[1])) n<-length(data) else n<-dim(data)[1]
	data = matrix(data,ncol=1)
	n<-dim(data)[1]
	s <- seq(1,n-k+1,by)
	lens <- length(s)
	cols <- if (ascending) 1:k else k:1
	return(matrix(data[s + rep(cols,rep(lens,k))-1],lens))
}

.lagx = function(data, n.lag=1, removeNA = FALSE, pad = NA)
{
	# has NAs
	data = as.matrix(data)
	n = dim(data)[1]
	d = dim(data)[2]
	if(dim(data)[2]==1) data=matrix(data,ncol=1)
	z = apply(data,2,FUN=function(x) .embed(x,n.lag+1)[,n.lag+1])
	if(!removeNA) z = rbind(matrix(pad,ncol=d,nrow=n.lag),z)
	return(z)
}

.lagmatrix = function(data, n.lag = 1, pad = 0)
{
	n = length(as.numeric(data))
	z = matrix(NA, ncol = n.lag, nrow = n)
	for(i in 1:n.lag) z[,i] = .lagx(as.numeric(data), i, removeNA = FALSE, pad = pad)
	return(z)
}


.abind = function(x, y)
{
	m = dim(x)[1]
	if( is.matrix(x) ) {
		n1 = 1
		x = array(x, dim = c(m, m, 1))
	} else{
		n1 = dim(x)[3]
	}
	if( is.matrix(y) ) {
		n2 = 1
		y = array(y, dim = c(m, m, 1))
	} else{
		n2 = dim(y)[3]
	}
	nw = array(NA, dim = c(m, m, n1 + n2))
	nw[ , , 1:n1] = x[, , 1:n1]
	nw[ , , (n1+1):(n1+n2)] = y[, , 1:n2]
	nw
}

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.