R/gogarch-main.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.
##
#################################################################################

# The next release should include the multivariate affine skew-student of Fernanfez and Steel,
# and perhaps some of their other skewed distribution which admit an affine representation of margins.

.gogarchspec = function(mean.model = list(demean = c("constant", "AR", "VAR", "robVAR"), lag = 1, lag.max = NULL, 
				lag.criterion = c("AIC", "HQ", "SC", "FPE"), external.regressors = NULL,
				robust.control = list(gamma = 0.25, delta = 0.01, nc = 10, ns = 500)), 
		variance.model = list(model = "sGARCH", garchOrder = c(1,1), submodel = NULL, variance.targeting = FALSE), 
		distribution.model = list(distribution = c("mvnorm", "manig", "magh")),
		ica = c("fastica", "pearson", "jade", "radical"), ica.fix = list(A = NULL, K = NULL), ...)
{
	mmodel = list(demean = "constant", lag = 1, lag.max = NULL, lag.criterion = "AIC", external.regressors = NULL,
			robust.control = list(gamma = 0.25, delta = 0.01, nc = 10, ns = 500))
	mmatch = match(names(mean.model), c("demean", "lag", "lag.max", "lag.criterion", "external.regressors", "robust.control"))
	if(length(mmatch[!is.na(mmatch)]) > 0){
		mx = which(!is.na(mmatch))
		mmodel[mmatch[!is.na(mmatch)]] = mean.model[mx]
	}
	valid.demean = c("constant", "AR", "VAR", "robVAR")
	if(!any(mmodel$demean[1] == valid.demean)) stop("\ngogarchspec-->error: Invalid demean choice\n")
	valid.criterion = c("AIC", "HQ", "SC", "FPE")
	if(!any(mmodel$lag.criterion[1] == valid.criterion)) stop("\ngogarchspec-->error: Invalid VAR criterion choice\n")
	mmodel$model = mmodel$demean
	if(mmodel$model == "robVAR"){
		mmodel$model = "VAR"
		mmodel$demean = "VAR"
		mmodel$robust = TRUE
	} else{
		mmodel$robust = FALSE
	}
	
	if(!is.null(mean.model$external.regressors)){
		mmodel$external.regressors = as.matrix(mean.model$external.regressors)
		mmodel$mxn = dim(mmodel$external.regressors)[2]
	} else{
		mmodel$external.regressors = NULL
		mmodel$mxn = 0
	}
	if( mmodel$model == "constant" ) mmodel$lag = 0
	
	vmodel =  list(model = "sGARCH", garchOrder = c(1,1), submodel = NULL, variance.targeting = FALSE)
	vmatch = match(names(variance.model), c("model", "garchOrder", "submodel", "variance.targeting"))
	if(length(vmatch[!is.na(vmatch)]) > 0){
		vx = which(!is.na(vmatch))
		vmodel[vmatch[!is.na(vmatch)]] = variance.model[vx]
	}
	# we let the univariate GARCH specification take care of the other checks
	
	dmodel = list(distribution = "mvnorm")
	dmatch = match(names(distribution.model), c("distribution"))
	if(length(dmatch[!is.na(dmatch)]) > 0){
		dx = which(!is.na(dmatch))
		dmodel[dmatch[!is.na(dmatch)]] = distribution.model[dx]
	}
	dmodel$distribution = dmodel$distribution[1]
	
	valid.dist = c("mvnorm", "manig", "magh")
	if(!any(dmodel$distribution == valid.dist)) stop("\ngogarchspec-->error: Invalid distribution choice\n")
	
	dist = switch(dmodel$distribution,
			manig = "nig",
			magh  = "ghyp",
			mvnorm = "norm")
	
	# this is the indepedent factor specification
	spec = ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE),
				variance.model = list(model = vmodel$model, garchOrder = vmodel$garchOrder, 
						submodel = vmodel$submodel, variance.targeting = vmodel$variance.targeting),
				distribution.model = dist)
	
	mspec = list()
	
	if(is.null(ica)){
		mspec$ica = "fastica"
	} else{
		valid.models = c("fastica", "pearson", "jade", "radical")
		if(!any(tolower(ica[1]) == valid.models)) stop("\ngogarchspec-->error: Invalid ICA model chosen\n", call. = FALSE)
		mspec$ica = tolower(ica[1])
	}
	
	if(is.null(ica.fix)) mspec$ica.fix = list() else mspec$ica.fix = ica.fix
	
	mspec$mmodel = mmodel
	
	mspec$distribution = dmodel$distribution

	ans = new("goGARCHspec",
			mspec = mspec,
			uspec = spec)
	return(ans)
}


.gogarchfit = function(spec, data, out.sample = 0, solver = "solnp", fit.control = list(stationarity = 1), 
		solver.control = list(), parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	if(is.null(data)) stop("\ngogarchfit-->error: gogarchfit requires a data object", call. = FALSE)
	if(!is.matrix(data) & !is.data.frame(data)) stop("\ngogarchfit-->error: gogarchfit only supports matrix or data.frame objects for the data", call. = FALSE)
	md = spec@mspec$distribution
	ans = switch(md,
			manig  = .gogarchfit.nig(spec = spec, data = data, out.sample = out.sample, solver = solver, 
					fit.control = fit.control, solver.control = solver.control, 
					parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
			magh   = .gogarchfit.ghyp(spec = spec, data = data, out.sample = out.sample, solver = solver, 
					fit.control = fit.control, solver.control = solver.control, 
					parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
			mvnorm = .gogarchfit.norm(spec = spec, data = data, out.sample = out.sample, solver = solver, 
					fit.control = fit.control, solver.control = solver.control, 
					parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit, ...))
	return(ans)
}

.gogarchfilter = function(fit, data, out.sample = 0, n.old = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	if(is.null(data)) stop("\ngogarchfilter-->error: gogarchfilter requires a data object", call. = FALSE)
	if(!is.matrix(data) & !is.data.frame(data)) stop("\ngogarchfilter-->error: gogarchfilter only supports matrix or data.frame objects for the data", call. = FALSE)
	md = fit@mfit$mdist
	ans = switch(md,
			manig   = .gogarchfilter.nig(fit = fit, data = data, out.sample = out.sample, 
					n.old = n.old, parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
			magh    = .gogarchfilter.ghyp(fit = fit, data = data, out.sample = out.sample, 
					n.old = n.old, parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
			mvnorm  = .gogarchfilter.norm(fit = fit, data = data, out.sample = out.sample, 
					n.old = n.old, parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit, ...))
	return(ans)
}

.gogarchforecast = function(fit, n.ahead = 10, n.roll = 0, external.forecasts = list(mregfor = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	md = fit@mfit$mdist
	ans = switch(md,
			manig	= .gogarchforecast.nig(fit = fit, n.ahead = n.ahead, n.roll = n.roll, 
					external.forecasts = external.forecasts, parallel = parallel, parallel.control = parallel.control, ...),
			magh	= .gogarchforecast.ghyp(fit = fit, n.ahead = n.ahead, n.roll = n.roll, 
					external.forecasts = external.forecasts, parallel = parallel, parallel.control = parallel.control, ...),
			mvnorm	= .gogarchforecast.norm(fit = fit, n.ahead = n.ahead, n.roll = n.roll, 
					external.forecasts = external.forecasts, parallel = parallel, parallel.control = parallel.control, ...))
	return(ans)
}

.gogarchsim = function(fit, n.sim = 1, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"), 
		prereturns = NA, mexsimdata = NULL, rseed = NULL, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	md = fit@mfit$mdist
	ans = switch(md,
			manig	= .gogarchsim.nig(fit,  n.sim, n.start, m.sim, startMethod, prereturns, mexsimdata, rseed, parallel,
					parallel.control, ...),
			magh	= .gogarchsim.ghyp(fit, n.sim, n.start, m.sim, startMethod, prereturns, mexsimdata, rseed, parallel,
					parallel.control, ...),
			mvnorm	= .gogarchsim.norm(fit,  n.sim, n.start, m.sim, startMethod, prereturns, mexsimdata, rseed, parallel,
					parallel.control, ...))
	return(ans)
}

.rollgogarch = function(spec, data, n.ahead = 1, forecast.length = 50, 
		refit.every = 25, refit.window = c("recursive", "moving"), parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
		solver = "solnp", fit.control = list(), solver.control = list(), rseed = NULL, 
		save.fit = FALSE, save.wdir = NULL, ...)
{
	if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
	forecast.length = floor(forecast.length/refit.every) * refit.every
	
	if(is.null(rseed)) rseed = as.numeric(Sys.time()) else rseed = as.integer(rseed)
	WD <- getwd()
	if(is.null(save.wdir)){
		if (!is.null(WD)) setwd(WD)
	} else{
		ND = save.wdir
		if (!is.null(ND)) setwd(ND)
	}
	# steps:
	# 1. Fit using out.sample (gofit)
	# 2. Forecast with rolling (goforecast)
	# 3. Convolution (convolution)
	# 4. VaR (makeQ)
	# Setup Indexed Data for Fitting
	asset.names = colnames(data)
	xdata = .extractdata(data[,1,drop = FALSE])
	data = as.matrix(data)
	dates = xdata$pos
	if(is.null(fit.control$stationarity)) fit.control$stationarity = 1
	if(is.null(fit.control$fixed.se)) fit.control$fixed.se = 0
	T = dim(data)[1]
	startT = T - forecast.length
	include.skew  = spec@uspec@distribution.model$include.skew
	include.shape = spec@uspec@distribution.model$include.shape
	# forecast window index
	fwindex = t(.embed((T-forecast.length - n.ahead + 2):T, refit.every, by = refit.every, TRUE ) )
	fitindx.start = c(fwindex[1,] - 1)
	fwindex.end = fwindex[refit.every,]
	nf = length(fwindex.end)
	fitdata = vector(mode = "list", length = nf)
	gofit = vector(mode = "list", length = nf)
	
	#gofit = vector(mode = "list", length = nf)
	outdata = vector(mode = "list", length = nf)
	distribution = spec@mspec$distribution
	forecastlist = vector(mode = "list", length = nf)
	
	mi = 0
	for(i in 1:nf){
		if(refit.window[1]=="recursive"){
			fitdata = data[1:(fitindx.start[i]+refit.every),]
		} else{
			fitdata = data[(i*refit.every):(fitindx.start[i]+refit.every),]
		}
		if(i == 1){
			set.seed(rseed)
			gofit = gogarchfit(spec, data = fitdata, out.sample = refit.every, solver = solver,
					fit.control = fit.control, solver.control = solver.control, parallel = parallel, parallel.control = parallel.control, ...)
			A = gofit@mfit$A
			forecastlist[[i]] = gogarchforecast(gofit, n.ahead = n.ahead, n.roll = refit.every-1, parallel = parallel, parallel.control = parallel.control)
			spec@mspec$ica.fix = list(A = NULL, K = NULL)
		} else{
			
			gofit = gogarchfit(spec, data = fitdata, out.sample = refit.every, solver = solver,
					fit.control = fit.control, solver.control = solver.control, parallel = parallel, parallel.control = parallel.control, 
					A.init = A, ...)
			
			A = gofit@mfit$A
			forecastlist[[i]] = gogarchforecast(gofit, n.ahead = n.ahead, n.roll = refit.every-1, parallel = parallel, parallel.control = parallel.control)
		}
		if(save.fit){
			eval(parse(text = paste("gofitroll",i,"=gofit",sep = "")))
			eval(parse(text = paste("save(gofitroll",i,",file='gofitroll",i,".rda')",sep = "")))
		}
		if(trace) print(i)
		mi = mi+1
	}
	mspec = list()
	mspec$refit.every = refit.every
	mspec$forecast.length = forecast.length
	mspec$asset.names = asset.names
	mspec$gospec = spec
	ans = new("goGARCHroll",
			forecast = forecastlist,
			spec = mspec)
	# return to the current directory
	setwd(WD)
	return(ans)
}