R/arfima-main.R

Defines functions fitandextractarfima arfimadistribution arfimaroll arfimapath arfimasim arfimaforecast2 arfimaforecast arfimafilter arfimaLLH arfimafit arfimaspec

Documented in arfimadistribution arfimafilter arfimafit arfimaforecast arfimapath arfimaroll arfimasim arfimaspec

#################################################################################
##
##   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.
##
#################################################################################

#---------------------------------------------------------------------------------
# SECTION ARFIMA spec
#---------------------------------------------------------------------------------
.arfimaspec = function(mmodel, dmodel, start.pars, fixed.pars)
{
	# extract values:
	# mean model
	include.mean = mmodel$include.mean
	armap = mmodel$armaOrder[1]
	armaq = mmodel$armaOrder[2]
	arfima = mmodel$arfima
	if(!is.null(mmodel$external.regressors)){
		include.mex = TRUE
		mexdata = as.matrix(mmodel$external.regressors)
		mxn = dim(mexdata)[2]
	} else{
		include.mex = FALSE
		mexdata = NULL
		mxn = 0
	}
	
	# dmodel/skew/shape
	include.dlambda = dmodel$include.dlambda
	include.skew = dmodel$include.skew
	include.shape = dmodel$include.shape
	maxOrder = max(c(armap, armaq))
	
	mmodel = list(include.mean = include.mean, armaOrder = mmodel$armaOrder, arfima = arfima,
			include.mex = include.mex, mexdata = mexdata, mxn = mxn)
	
	#specnames = c("mu","ar","ma","darfima","mxreg")

	modelnames = c(
			if(include.mean) "mu",
			if(armap > 0) paste("ar", 1:armap, sep = ""),
			if(armaq > 0) paste("ma", 1:armaq, sep = ""),
			if(arfima) paste("darfima"),
			if(include.mex)  paste("mxreg", 1:mxn, sep = ""),
			paste("sigma"),
			if(include.dlambda) "dlambda",
			if(include.skew) "skew",
			if(include.shape) "shape")
	
	pos = 1
	pos.matrix = matrix(0, ncol = 3, nrow = 9)
	colnames(pos.matrix) = c("start", "stop", "include")
	rownames(pos.matrix) = c("mu", "ar", "ma", "darfima", "mxreg", "sigma", "dlambda", "skew", "shape")
	if(include.mean){ pos.matrix[1,1:3] = 1
		pos = pos.matrix[1,2]+1}
	if(armap>0){pos.matrix[2,1:3] = c(pos, pos+armap-1, 1)
		pos = max(pos.matrix[1:2,2])+1}
	if(armaq>0){pos.matrix[3,1:3] = c(pos, pos+armaq-1, 1)
		pos = max(pos.matrix[1:3,2])+1}
	if(arfima){pos.matrix[4,1:3] = c(pos, pos, 1)
		pos = max(pos.matrix[1:4,2])+1}
	if(include.mex){pos.matrix[5,1:3] = c(pos, pos+mxn-1, 1)
		pos = max(pos.matrix[1:5,2])+1}
	pos.matrix[6,1:3] = c(pos,pos,1)
	pos = max(pos.matrix[1:6,2])+1
	if(include.dlambda){pos.matrix[7,1:3] = c(pos, pos, 1)
		pos = max(pos.matrix[1:7,2])+1}
	if(include.skew){pos.matrix[8,1:3] = c(pos, pos, 1)
		pos = max(pos.matrix[1:8,2])+1}
	if(include.shape){pos.matrix[9,1:3] = c(pos, pos, 1)}
	
	nn = length(modelnames)
	modelmatrix = matrix(0, ncol = 3, nrow = nn)
	rownames(modelmatrix) = modelnames
	colnames(modelmatrix) = c("opt", "fixed", "start")
	fixed.names = names(fixed.pars)
	fp = charmatch(fixed.names, modelnames)
	
	if(!is.null(fixed.names) && any(!is.na(fp))){
		fixed = fp[!is.na(fp)]
		modelmatrix[fixed,2] = 1
		fz = charmatch(modelnames, fixed.names)
		fz = fz[!is.na(fz)]
		fixed.pars = fixed.pars[fz]
		names(fixed.pars) = fixed.names[fz]
	} else{
		fixed.pars = NULL
	}
	modelmatrix[,1] = 1 - modelmatrix[,2]
	start.names = names(start.pars)
	sp = charmatch(start.names, modelnames)
	if(!is.null(start.names) && any(!is.na(sp))){
		start = sp[!is.na(sp)]
		modelmatrix[start,3] = 1
		sz = charmatch(modelnames, start.names)
		sz = sz[!is.na(sz)]
		start.pars = start.pars[sz]
	} else{
		start.pars = NULL
	}
	# the way we handle fixed pars will be through upper and lower bound constraints
	opt.pars = vector(mode = "numeric", length = length(modelnames))
	names(opt.pars) = modelnames
	# optimization model
	omodel = list(modelnames = modelnames, opt.pars = opt.pars, start.pars = start.pars, 
			fixed.pars = fixed.pars, modelmatrix = modelmatrix, maxOrder = maxOrder, 
			pos.matrix = pos.matrix)
	
	ans = new("ARFIMAspec",
			mean.model = mmodel,
			distribution.model = dmodel,
			optimization.model = omodel)
	return(ans)
}


#---------------------------------------------------------------------------------
# SECTION sGARCH fit
#---------------------------------------------------------------------------------
.arfimafit = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(), 
		fit.control = list(fixed.se = 0, scale = 0))
{
	tic = Sys.time()
	if(is.null(fit.control$fixed.se)) fit.control$fixed.se = 0
	if(is.null(fit.control$scale)) fit.control$scale = 0
	# at the moment we do not allow scaling with fixed parameters
	if(!is.null(spec@optimization.model$fixed.pars)) fit.control$scale = 0
	xdata = .extractdata(data)
	if(!is.numeric(out.sample)) stop("\narfimafit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample)<0) stop("\narfimafit-->error: out.sample must be positive\n")
	n.start = round(out.sample,0)
	n = length(xdata$data)
	if((n-n.start)<100) stop("\narfimafit-->error: function requires at least 100 data\n points to run\n")
	data = xdata$data[1:(n-n.start)]
	dates = xdata$pos[1:(n-n.start)]
	origdata = xdata$data
	origdates = xdata$pos
	dformat = xdata$dformat
	# create a temporary environment to store values (deleted at end of function)
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	# expand the spec object and assign spec lists
	mmodel = spec@mean.model
	dmodel = spec@distribution.model
	omodel = spec@optimization.model
	if(mmodel$include.mex) {
		mmodel$origmexdata = mmodel$mexdata
		mmodel$mexdata = mmodel$mexdata[1:(n-n.start),]
		colnames(mmodel$mexdata) = paste("mexdata", 1:dim(mmodel$mexdata)[2], sep = "")
	} else{
		mmodel$origmexdata = NULL
	}
	
	assign("mmodel", mmodel, envir = tempenvir)
	assign("dmodel", dmodel, envir = tempenvir)
	assign("omodel", omodel, envir = tempenvir)
	assign("dates", dates, envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	include.mean = mmodel$include.mean
	m =  omodel$maxOrder
	T = length(as.numeric(data))
	dist = dmodel$distribution
	if(fit.control$scale) dscale = sd(data) else dscale = 1
	zdata = data/dscale
	assign("dscale", dscale, envir = tempenvir)
	# Optimization Starting Parameters Vector & Bounds
	ipars = .arfimastart(data = zdata, garchenv = tempenvir)
	pars = ipars$pars
	LB = ipars$modelLB
	UB = ipars$modelUB
	assign("npar", length(pars), envir = tempenvir)
	# we now split out any fixed parameters
	if( any(LB == UB) ){
		fixid = which( LB == UB )
		if( length(fixid) == length(pars) ){
			if( fit.control$fixed.se == 0 ) {
				# if all parameters are fixed an no standard erros are to
				# be calculated then we return a ugarchfilter object
				cat("\narfimafit-->warning: all parameters fixed...returning arfimafilter object instead\n")
				return( arfimafilter(data = data, spec = spec) )
			} else{
				# if all parameters are fixed but we require standard errors, we
				# skip the solver
				use.solver = 0
				fixed = pars[fixid]
				pars = pars[-fixid]
			}
		} else{
			# with some parameters fixed we extract them (to be rejoined at end)
			# so that they do not enter the solver
			use.solver = 1
			fixed = pars[fixid]
			pars = pars[-fixid]
			LB = LB[-fixid]
			UB = UB[-fixid]
		}
	} else{
		# no fixed parameters, all normal
		fixed = NULL
		fixid = NULL
		use.solver = 1
	}
	assign("fixed", fixed, envir = tempenvir)
	assign("fixid", fixid, envir = tempenvir)
	# start counter
	assign(".llh", 1, envir = tempenvir)
	# scale parameters (for non-solnp solvers)
	parscale = rep( 1, length = length(pars) )
	names(parscale) = names(pars)
	if(include.mean) parscale["mu"] = abs(mean(zdata))
	parscale["sigma"] = sd(zdata)
	Ifn = function(pars, data, returnType, garchenv){
		arx = 0.5
		omodel = get("omodel", garchenv)
		pos = omodel$pos.matrix		
		if(pos[2,3] == 1 && pos[4,3] == 0) arx = pars[pos[2,1]:pos[2,2]]
		min(Mod(polyroot(c(1, -arx))))
	}
	ILB = 1.01
	IUB = 100
	
	if(use.solver){
		solution = .garchsolver(solver, pars, fun = .arfimaLLH, Ifn, ILB, IUB, 
				gr = NULL, hessian = NULL, parscale = parscale, 
				control = solver.control, LB = LB, UB = UB, ux = NULL, ci = NULL, 
				mu = NULL, data = zdata, returnType = "llh", garchenv = tempenvir)
		sol = solution$sol
		hess = solution$hess
		timer = Sys.time()-tic
		opt.pars = sol$par
		if(!is.null(opt.pars)) names(opt.pars) = names(pars)
		if(is.null(fixed)){
			if(include.mean) opt.pars["mu"] = opt.pars["mu"] * dscale
			if(mmodel$include.mex){
				opt.pars[omodel$pos[5,1]:omodel$pos[5,2]] = opt.pars[omodel$pos[5,1]:omodel$pos[5,2]] * dscale
			}
			opt.pars["sigma"] = (opt.pars["sigma"] * dscale)
		}
		convergence = sol$convergence
	} else{
		hess = NULL
		timer = Sys.time()-tic
		opt.pars = NULL
		convergence = 0
		sol = list()
		sol$message = "all parameters fixed"
	}
	fit = list()
	model = list()
	# check convergence else write message/return
	if( convergence == 0 ){
		if( !is.null(fixed) && fit.control$fixed.se == 1 ){
			fixed =  get("fixed", tempenvir)
			npar = get("npar", tempenvir)
			pall = vector(mode = "numeric", length = npar)
			pall[fixid] = fixed
			pall[-fixid] = opt.pars
			opt.pars = pall
			Names = omodel$modelnames
			names(opt.pars) = Names
			assign("fixed", NULL, envir = tempenvir)
			assign("fixid", NULL, envir = tempenvir)
		}
		fit = .makearfimafitmodel(f = .arfimaLLH, opt = opt.pars, data = data,  T = T, m = m, timer = timer, 
				convergence = convergence, message = sol$message, hess, garchenv = tempenvir)
		fit$dates = dates
		fit$date.format = dformat
		fit$origdata = origdata
		fit$origdates = origdates
	} else{
		fit$message = sol$message
		fit$convergence = 1
	}
	
	# make model list to return some usefule information which
	# will be called by other functions (show, plot, sim etc)
	model = .makearfimamodel(tempenvir)
	model$n.start = n.start
	# return object ARFIMAfit
	ans = new("ARFIMAfit",
			fit = fit,
			model = model)
	rm(tempenvir)
	return(ans)
}

#---------------------------------------------------------------------------------
# SECTION sGARCH LLH
#---------------------------------------------------------------------------------
.arfimaLLH = function(pars, data, returnType = "llh", garchenv)
{
	# prepare inputs
	# rejoin fixed and pars
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	npar = get("npar", garchenv)
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		pall[fixid] = fixed
		pall[-fixid] = pars
		pars = pall
	}
	T = length(data)
	mmodel = get("mmodel", garchenv)
	dmodel = get("dmodel", garchenv)
	omodel = get("omodel", garchenv)
	fit.control = get("fit.control", garchenv)
	armaOrder = c(mmodel$armaOrder,
			0,
			mmodel$arfima,
			mmodel$include.mex)
	m = omodel$maxOrder
	N = c(m,T)
	mexdata = mmodel$mexdata
	distribution = dmodel$distribution
	# distribution number
	# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
	dist = dmodel$distno
	pos = omodel$pos.matrix
	
	mu = ar = ma = darfima = mxreg = skew = shape = dlambda = 0
	if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
	if(pos[2,3] == 1) ar = pars[pos[2,1]:pos[2,2]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	if(pos[4,3] == 1) darfima = pars[pos[4,1]:pos[4,2]]
	if(pos[5,3] == 1) mxreg = pars[pos[5,1]:pos[5,2]]
	if(pos[6,3] == 1) sigma = pars[pos[6,1]:pos[6,2]]
	if(pos[7,3] == 1) dlambda = pars[pos[7,1]:pos[7,2]]
	if(pos[8,3] == 1) skew = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) shape = pars[pos[9,1]:pos[9,2]]
	
	rx = .arfimaxfilterx(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, darfima = darfima, 
			mxreg = mxreg, mexdata = mexdata, data = data, N = N, garchenv)
	
	res = rx$res
	zrf = rx$zrf
	
	mvar = mean(res*res)
	
	# C Code [sgarch.c] via wrapper
	
	ans = .arfimafitC(armaOrder, mu, ar, ma, mxreg, mexdata, sigma, dlambda, skew, shape, 
			dist, data, zrf, N, res, garchenv)
	
	if(get(".csol", garchenv) == 1){
		cat(paste("\narfimafit-->warning: ", get(".filtermessage",garchenv),"\n", sep=""))
		return(llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv)))))
	}
	
	z = ans$z
	sigma = ans$sigma
	epsx = ans$res
	llh = ans$llh
	
	if(is.finite(llh) | !is.na(llh) | !is.nan(llh)){
		assign(".llh", llh, envir = garchenv)
	} else {
		llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	# LHT = raw scores
	LHT = -ans$LHT
	ans = switch(returnType,
			llh = llh,
			LHT = LHT,
			all = list(llh = llh, sigma = sigma, epsx = epsx, z = z, LHT = LHT))
	
	return(ans)
}

#---------------------------------------------------------------------------------
# SECTION sGARCH filter
#---------------------------------------------------------------------------------
.arfimafilter = function(spec, data, out.sample = 0, n.old = NULL)
{
	# n.old is optional and indicates the length of the original dataseries (in
	# cases when this represents a dataseries augmented by newer data). The reason
	# for using this is so that the old and new datasets agree since the original
	# recursion uses the sum of the residuals to start the recursion and therefore
	# is influenced by new data. For a small augmentation the values converge after
	# x periods, but it is sometimes preferable to have this option so that there is
	# no forward looking information contaminating the study.
	.garchenv = environment()
	if(!is.null(n.old)) Nx = n.old else Nx = length(data)
	xdata = .extractdata(data)
	data = xdata$data
	dates = xdata$pos
	dformat = xdata$dformat
	origdata = data
	origdates = dates
	T = length(origdata)  - out.sample
	data = origdata[1:T]
	dates = origdates[1:T]
	
	mmodel = spec@mean.model
	dmodel = spec@distribution.model
	omodel = spec@optimization.model
	include.mean = mmodel$include.mean
	modelnames = omodel$modelnames
	pars = unlist(spec@optimization.model$fixed.pars)
	parnames = names(pars)
	if(length(modelnames) != length(pars))
		stop("\narfimafilter-->error: parameters do not match expected specification length")
	if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
		stop("\narfimafilter-->error: parameters names do not match specification")
	} else{
		pars = pars[modelnames]
	}
	m =  omodel$maxOrder
	T = length(as.numeric(data))
	dist = dmodel$distribution
	armaOrder = c(mmodel$armaOrder,
			0,
			mmodel$arfima,
			mmodel$include.mex)

	mexdata = mmodel$mexdata
	distribution = dmodel$distribution
	# distribution number
	# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
	dist = dmodel$distno
	N = c(m,T)
	pos = omodel$pos.matrix
	
	mu = ar = ma = darfima = mxreg = skew = shape = dlambda = 0
	if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
	if(pos[2,3] == 1) ar = pars[pos[2,1]:pos[2,2]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	if(pos[4,3] == 1) darfima = pars[pos[4,1]:pos[4,2]]
	if(pos[5,3] == 1) mxreg = pars[pos[5,1]:pos[5,2]]
	if(pos[6,3] == 1) sigma = pars[pos[6,1]:pos[6,2]]
	if(pos[7,3] == 1) dlambda = pars[pos[7,1]:pos[7,2]]
	if(pos[8,3] == 1) skew = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) shape = pars[pos[9,1]:pos[9,2]]
	
	
	rx = .arfimaxfilterx(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, darfima = darfima, 
			mxreg = mxreg, mexdata = mexdata, data = data, N = N, .garchenv)
	res = rx$res
	zrf = rx$zrf
	
	if(!is.null(n.old)){
		rx2 = .arfimaxfilterx(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, darfima = darfima, 
				mxreg = mxreg, mexdata = mexdata[1:Nx, ,drop = FALSE], data = data[1:Nx], N = c(m, Nx), .garchenv)
		res2 = rx2$res
		# unconditional sigma value
		mvar = mean(res2*res2)
	} else{
		mvar = mean(res*res)
	}
	
	ans = .arfimafitC(armaOrder, mu, ar, ma, mxreg, mexdata, sigma, dlambda, skew, shape, 
			dist, data, zrf, N, res, .garchenv)
	
	filter = list()
	filter$z = ans$z
	filter$residuals = ans$res
	filter$LLH = -ans$llh
	filter$log.likelihoods = ans$LHT
	filter$data = data
	filter$dates = dates
	filter$origdata = origdata
	filter$origdates = origdates
	filter$date.format = dformat
	filter$distribution = distribution
	filter$pars = pars
	model = c(mmodel, dmodel)
	rm(.garchenv)
	sol = new("ARFIMAfilter",
			filter = filter,
			model = model)
	return(sol)
}

#---------------------------------------------------------------------------------
# SECTION sGARCH forecast
#---------------------------------------------------------------------------------
.arfimaforecast = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0, 
		external.forecasts = list(mregfor = NULL), ...)
{
	fit = fitORspec
	data = fit@fit$data
	N = length(as.numeric(data))
	origData = fit@fit$origdata
	origDates = fit@fit$origdates
	dformat = fit@fit$date.format
	Nor = length(as.numeric(origData))
	ns = fit@model$n.start
	if( n.roll > ns ) stop("\narfimaforecast-->error: n.roll must not be greater than out.sample!")
	
	include.mean = fit@model$include.mean
	armaOrder = fit@model$armaOrder
	arfima = fit@model$arfima
	include.mex = fit@model$include.mex
	mxn = fit@model$mxn
	if( mxn > 0 ) mexdata = matrix(fit@model$mexdata, ncol = mxn) else mexdata = NULL
	modelnames = fit@model$modelnames
	pars = fit@fit$coef
	distribution = fit@model$distribution
	include.dlambda = fit@model$include.dlambda
	include.skew = fit@model$include.skew
	include.shape = fit@model$include.shape
	
	# check if necessary the external regressor forecasts provided first
	xreg = .forcregressors(include.mex, mxn, fit@model$mexdata, mregfor = external.forecasts$mregfor, 
			FALSE, NULL, NULL, NULL, pars, n.ahead, Nor, out.sample = ns, n.roll)
	
	mxreg = xreg$mxreg
	mxf = xreg$mxf
	
	pos = fit@model$pos
	
	ar = armap = ma = armaq = darfima = 0
	skew = shape = dlambda = 0
	
	if(pos[2,3] == 1){
		ar = pars[pos[2,1]:pos[2,2]]
		armap = length(ar)
	}
	
	if(pos[3,3] == 1){ 
		ma = pars[pos[3,1]:pos[3,2]]
		armaq = length(ma)
	}
	
	if(pos[4,3] == 1) darfima = pars[pos[4,1]:pos[4,2]]
	
	sigma = pars[pos[6,1]:pos[6,2]]
		
	if(pos[7,3] == 1) dlambda = pars[pos[7,1]:pos[7,2]]
	if(pos[8,3] == 1) skew = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) shape = pars[pos[9,1]:pos[9,2]]
	
	
	# filter data (check external regressor data - must equal length of origData)
	fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead + n.roll, ns)
	fspec = arfimaspec(
			mean.model = list(armaOrder = c(armap, armaq),
					include.mean = include.mean, arfima = arfima, 
					external.regressors = mxf[1:(N + fcreq), , drop = FALSE]), 
			distribution.model = distribution, fixed.pars = as.list(pars))
	tmp =  data.frame(origData[1:(N + fcreq)])
	rownames(tmp) = as.character(origDates[1:(N + fcreq)])
	flt = .arfimafilter(data = tmp[,1,drop=FALSE], spec = fspec, n.old = length(fit@fit$data))
	resfilter = flt@filter$residuals
	zfilter = flt@filter$z
	# forecast GARCH process
	forecasts = vector(mode = "list", length = n.roll+1)
	fwdd = vector(mode="list", length = n.roll+1)
	for(i in 1:(n.roll+1)){
		np = N + i - 1
		if(pos[1,3] == 1){
			mu = rep(pars[pos[1,1]:pos[1,2]], N+i+n.ahead-1)
		} else{
			mu = rep(0, N+i+n.ahead-1)
		}
		epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
		x = c(origData[1:(N+i-1)], rep(0, n.ahead))
		z = c(zfilter[1:(N+i-1)], rep(0, n.ahead))
		# forecast of externals is provided outside the system
		mxfi = mxf[1:(N+i-1+n.ahead), , drop = FALSE]
		if(arfima){
			ans = arfimaf(x[1:np], h = 0, armaOrder, mu, ar, ma, darfima, inmean = 0, im = 0, 
					epsx, include.mex, mxreg, mxfi, n.ahead)
		} else{
			ans = armaf(x[1:np], h = 0, armaOrder, mu, ar, ma, inmean = 0, im = 0,
					epsx, include.mex, mxreg, mxfi, n.ahead)
		}
		ans = ans[(np + 1):(np + n.ahead)]
		fdf = data.frame(series = ans)
		fwdd[[i]] = .forcdates( origDates, n.ahead, N, i, ns, dformat )
		rownames(fdf) = as.character(fwdd[[i]])
		forecasts[[i]] = fdf
	}
	fcst = list()
	fcst$n.ahead = n.ahead
	fcst$N = N + ns
	fcst$n.start = ns
	fcst$n.roll = n.roll
	fcst$series = fit@fit$data
	fcst$forecasts = forecasts
	fcst$fdates = fwdd
	model = fit@model
	model$dates = fit@fit$dates
	model$origdata = fit@fit$origdata
	model$origdates = fit@fit$origdates
	ans = new("ARFIMAforecast",
			forecast = fcst,
			model = model,
			flt)
	return(ans)
}

#---------------------------------------------------------------------------------
# 2nd dispatch method for forecast
.arfimaforecast2 = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0, 
		external.forecasts = list(mregfor = NULL), ...)
{
	# first we filter the data to get the results:
	spec = fitORspec
	if(is.null(data))
		stop("\nugarchforecast-->error: data must not be NULL when using a specification!")
	# we then extract the data/coefs etc
	xdata = .extractdata(data)
	Nor = length(as.numeric(xdata$data))
	data = xdata$data[1:(Nor - out.sample)]
	N = length(as.numeric(data))
	dates = xdata$pos[1:(Nor - out.sample)]
	origData =  xdata$data
	origDates =  xdata$pos
	dformat = xdata$dformat
	mmodel = spec@mean.model
	dmodel = spec@distribution.model
	omodel = spec@optimization.model
	
	ns = out.sample
	if( n.roll > ns )
		stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
	
	include.mean = mmodel$include.mean
	armaOrder = mmodel$armaOrder
	arfima = mmodel$arfima
	include.mex = mmodel$include.mex
	mxn = mmodel$mxn
	if(mxn>0) mexdata = matrix(mmodel$mexdata, ncol = mxn) else mexdata = NULL
	pow = 2
	modelnames = omodel$modelnames
	pars = unlist(omodel$fixed.pars)
	distribution = dmodel$distribution
	include.dlambda = dmodel$include.dlambda
	include.skew = dmodel$include.skew
	include.shape = dmodel$include.shape
	
	# check if necessary the external regressor forecasts provided first
	xreg = .forcregressors(include.mex, mxn, mexdata, mregfor = external.forecasts$mregfor, 
			FALSE, NULL, NULL, NULL, pars, n.ahead, Nor, out.sample = ns, n.roll)
	
	mxreg = xreg$mxreg
	mxf = xreg$mxf
	pos = omodel$pos
	
	ar = armap = ma = armaq = darfima = 0
	skew = shape = dlambda = 0
	
	if(pos[2,3] == 1){
		ar = pars[pos[2,1]:pos[2,2]]
		armap = length(ar)
	}
	
	if(pos[3,3] == 1){ 
		ma = pars[pos[3,1]:pos[3,2]]
		armaq = length(ma)
	}
	
	if(pos[4,3] == 1) darfima = pars[pos[4,1]:pos[4,2]]
	
	sigma = pars[pos[6,1]:pos[6,2]]
	
	if(pos[7,3] == 1) dlambda = pars[pos[7,1]:pos[7,2]]
	if(pos[8,3] == 1) skew = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) shape = pars[pos[9,1]:pos[9,2]]
	
	
	# filter data (check external regressor data - must equal length of origData)
	fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead + n.roll, ns)
	fspec = arfimaspec(
			mean.model = list(armaOrder = c(armap, armaq),
					include.mean = include.mean, arfima = arfima, 
					external.regressors = mxf[1:(N + fcreq), ,drop = FALSE]), 
			distribution.model = distribution, fixed.pars = as.list(pars))
	tmp =  data.frame(origData[1:(N + fcreq)])
	rownames(tmp) = as.character(origDates[1:(N + fcreq)])
	flt = .arfimafilter(data = tmp[,1,drop=FALSE], spec = fspec, n.old = N)
	resfilter = flt@filter$residuals
	zfilter = flt@filter$z
	# forecast GARCH process
	forecasts = vector(mode = "list", length = n.roll+1)
	fwdd = vector(mode="list", length = n.roll+1)
	for(i in 1:(n.roll+1)){
		np = N + i - 1
		if(pos[1,3] == 1){
			mu = rep(pars[pos[1,1]:pos[1,2]], N+i+n.ahead-1)
		} else{
			mu = rep(0, N+i+n.ahead-1)
		}
		epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
		x = c(origData[1:(N+i-1)], rep(0, n.ahead))
		z = c(zfilter[1:(N+i-1)], rep(0, n.ahead))
		# forecast of externals is provided outside the system
		mxfi = mxf[1:(N+i-1+n.ahead), , drop = FALSE]
		
		if(arfima){
			ans = arfimaf(x[1:np], h = 0, armaOrder, mu, ar, ma, darfima, inmean = 0, im = 0, 
					epsx, include.mex, mxreg, mxf, n.ahead)
		} else{
			ans = armaf(x[1:np], h = 0, armaOrder, mu, ar, ma, inmean = 0, im = 0,
					epsx, include.mex, mxreg, mxf, n.ahead)
		}
		ans = ans[(np + 1):(np + n.ahead)]
		fdf = data.frame(series = ans)
		fwdd[[i]] = .forcdates( origDates, n.ahead, N, i, ns, dformat )
		rownames(fdf) = as.character(fwdd[[i]])
		forecasts[[i]] = fdf
	}
	fcst = list()
	fcst$n.ahead = n.ahead
	fcst$N = N + ns
	fcst$n.start = ns
	fcst$n.roll = n.roll
	fcst$series = data
	fcst$forecasts = forecasts
	fcst$fdates = fwdd
	model = c(mmodel, dmodel)
	model$dates = dates
	model$origdata = origData
	model$origdates = origDates
	ans = new("ARFIMAforecast",
			forecast = fcst,
			model = model,
			flt)
	return(ans)
	
}
#---------------------------------------------------------------------------------
# SECTION sGARCH simulate
#---------------------------------------------------------------------------------

.arfimasim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = 
				c("unconditional","sample"), prereturns = NA, 
		preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA, type = "z"), 
		mexsimdata = NULL)
{
	# some checks
	if (is.na(rseed[1])){
		sseed = as.integer(runif(1*m.sim,0,65000))
	} else{
		if(length(rseed) < m.sim){
			stop("\nugarchsim-->error: seed length must equal m.sim\n")
		} else {
			sseed = rseed[1:m.sim]
		}
	}
	# Enlarge Series:
	# need to allow for arfima case:
	arfima = fit@model$arfima
	armap = fit@model$armaOrder[1]
	armaq = fit@model$armaOrder[2]
	if(arfima) n.start = max(armap+armaq, n.start)
	
	n = n.sim + n.start
	startMethod = startMethod[1]
	data = fit@fit$data
	N = length(as.numeric(data))
	m = fit@model$maxOrder
	resids = fit@fit$residuals
	include.mean = fit@model$include.mean
	armaOrder = fit@model$armaOrder
	arfima = fit@model$arfima
	include.mex = fit@model$include.mex
	mxn = fit@model$mxn
	if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol = mxn) else mexdata = NULL
	modelnames = fit@model$modelnames
	pars = fit@fit$coef
	distribution = fit@model$distribution
	include.dlambda = fit@model$include.dlambda
	include.skew = fit@model$include.skew
	include.shape = fit@model$include.shape
	
	# check if necessary the external regressor forecasts provided first
	xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex = FALSE, 
			vxn = NULL, vexsimdata = NULL, vexdata = NULL, N, pars, n, m.sim, m)
	
	if(mxn > 0) mxreg = matrix(xreg$mxreg, ncol = mxn)
	mexsim = xreg$mexsimlist
	
	pos = fit@model$pos
	
	ar = armap = ma = armaq = inmean = darfima = mu = 0
	skew = shape = dlambda = 0
	
	if(pos[1,3] == 1){
		mu = pars[pos[1,1]:pos[1,2]]
	}
	
	if(pos[2,3] == 1){
		ar = pars[pos[2,1]:pos[2,2]]
		armap = length(ar)
	}
	
	if(pos[3,3] == 1){ 
		ma = pars[pos[3,1]:pos[3,2]]
		armaq = length(ma)
	}
	
	if(pos[4,3] == 1) darfima = pars[pos[4,1]:pos[4,2]]
	
	sigma = pars[pos[6,1]:pos[6,2]]
	
	if(pos[7,3] == 1) dlambda = pars[pos[7,1]:pos[7,2]]
	if(pos[8,3] == 1) skew = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) shape = pars[pos[9,1]:pos[9,2]]
	
	
	if(N < n.start){
		startmethod[1] = "unconditional"
		cat("\narfimasim-->warning: n.start greater than length of data...using unconditional 
						start method...\n")
	}
	
	# Random Samples from the Distribution
	# For the arfima method this is the residuals, NOT the standardized residuals
	z = matrix(0, ncol = m.sim, nrow = n)
	zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim), 
			skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
	z = matrix(.custzdist(custom.dist, zmatrix, m.sim, n), ncol = m.sim)
	if( is.matrix( custom.dist$distfit ) && custom.dist$type != "z") cres = TRUE else cres = FALSE
	if(startMethod == "unconditional"){
		z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
	} else{
		z = rbind(matrix(tail(fit@fit$z, m), nrow = m, ncol = m.sim), z)
	}
	
	if(!any(is.na(prereturns))){
		prereturns = as.vector(prereturns)
		if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
	}
	if(!any(is.na(preresiduals))){
		preresiduals = as.vector(preresiduals)
		preres = preresiduals
		if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
	}
	if(any(is.na(prereturns))){
		if(startMethod[1] == "unconditional"){
			prereturns = as.numeric(rep(uncmean(fit), m))
		}
		else{
			prereturns = tail(data, m)
		}
	}
	
	# input vectors/matrices
	x = c(prereturns, rep(0, n))
	constm = matrix(mu, ncol = m.sim, nrow = n + m)
	
	# outpus matrices
	seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
	residSim =  matrix(0, ncol = m.sim, nrow = n.sim)
	
	for(i in 1:m.sim){
		if(any(is.na(preresiduals))){
			if(startMethod[1] == "sample"){
				preres 	= tail(resids, m)
				res 	= c(preres, if(cres) tail(z[(m+1):(n+m), i], n) else tail(z[(m+1):(n+m), i], n) * sigma)
				residSim[,i] = tail(res, n.sim)
			} else{
				res 	= if(cres) as.numeric(z[,i]) else as.numeric(z[,i])*sigma
				residSim[,i] = tail(res, n.sim)
			}
		}
		if(include.mex){
			constm[,i] = constm[,i] + mxreg%*%t(matrix(mexsim[[i]], ncol = mxn))
		}
		if(arfima){
			#if(constant) constm[,i] = constm[,i]*(1-sum(ar))
			x = c(prereturns, rep(0, n + armaq))
			res = c(if(armaq>0) rep(0, armaq) else NULL, preres[1:m], if(cres) as.numeric(z[(m+1):(n+m),i]) else as.numeric(z[(m+1):(n+m),i]) * sigma)
			ans2 = .arfimaxsim(armap, armaq, constm[1:(n + m), i], ar, ma, darfima, x[1:(n + m + armaq)], res[1:(n + m + armaq)], T = n + m)
			seriesSim[,i] = ans2$s[(m+1):(n+m)]
		} else{
			#if(constant) constm[,i] = constm[,i]*(1-sum(ar))
			ans2 = .armaxsim(armap, armaq, constm[,i], ar, ma, x, res, T = n + m, m)
			seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
			# constm[m+1,1] + sum(ar * rev(x[1:2] - constm[1:2,1])) + sum(ma * rev(res[1:2]))
		}
	}
	smodel = fit@model
	smodel$dates = fit@fit$dates
	sim = list(seriesSim = seriesSim, residSim = residSim)
	sim$series = data
	sol = new("ARFIMAsim",
			simulation = sim,
			model = smodel,
			seed = as.integer(sseed))
	return(sol)
}


#---------------------------------------------------------------------------------
# SECTION sGARCH path
#---------------------------------------------------------------------------------

.arfimapath = function(spec, n.sim = 1000, n.start = 0, m.sim = 1, prereturns = NA, preresiduals = NA, 
		rseed = NA, custom.dist = list(name = NA, distfit = NA, type = "z"), mexsimdata = NULL)
{
	# some checks
	if (is.na(rseed[1])){
		sseed = as.integer(runif(1*m.sim,0,65000))
	} else{
		if(length(rseed) < m.sim){
			stop("\nugarchsim-->error: seed length must equal m.sim\n")
		} else {
			sseed = rseed[1:m.sim]
		}
	}
	# Enlarge Series:
	n = n.sim + n.start
	
	armaOrder = spec@mean.model$armaOrder
	armap = armaOrder[1]
	armaq = armaOrder[2]
	include.mean = spec@mean.model$include.mean
	include.mex = spec@mean.model$include.mex
	mxn = spec@mean.model$mxn
	N = 0
	if(include.mex) {
		mexdata = matrix(spec@mean.model$mexdata, ncol=mxn)
		N = dim(mexdata)[1]
	} else { mexdata = NULL }
	distribution = spec@distribution.model$distribution
	pos = spec@optimization.model$pos.matrix
	pars = unlist(spec@optimization.model$fixed.pars)
	# perform checks on supplied parameters
	exp.pars = rownames(spec@optimization.model$modelmatrix)
	spec.pars = names(pars)
	tmatch = match(sort(spec.pars), sort(exp.pars))
	if(any(is.na(tmatch))){
		cat("\narfimapath-->error: supplied parameters do not match expected parameters\n")
		cat("expected parameters: \n")
		cat(paste(sort(exp.pars), " ", sep="")) 
		cat("\nsupplied parameters: \n")
		cat(paste(sort(spec.pars), " ", sep="")) 
		stop("\n", call. = FALSE)
	}
	
	mu = ar = ma = darfima = mxreg = skew = shape = dlambda = 0
	if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
	if(pos[2,3] == 1) ar = pars[pos[2,1]:pos[2,2]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	if(pos[4,3] == 1) darfima = pars[pos[4,1]:pos[4,2]]
	if(pos[5,3] == 1) mxreg = pars[pos[5,1]:pos[5,2]]
	if(pos[6,3] == 1) sigma = pars[pos[6,1]:pos[6,2]]
	if(pos[7,3] == 1) dlambda = pars[pos[7,1]:pos[7,2]]
	if(pos[8,3] == 1) skew = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) shape = pars[pos[9,1]:pos[9,2]]
	
	arfima = spec@mean.model$arfima
	m = spec@optimization.model$maxOrder
	modelnames = spec@optimization.model$modelnames
	
	include.dlambda = spec@distribution.model$include.dlambda
	include.skew = spec@distribution.model$include.skew
	include.shape = spec@distribution.model$include.shape
	
	# check if necessary the external regressor forecasts provided first
	xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex = FALSE, vxn = NULL, 
			vexsimdata = NULL, vexdata = NULL, N, pars, n, m.sim, m)
	
	if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
	mexsim = xreg$mexsimlist
	
	# Random Samples from the Distribution
	z = matrix(0, ncol = m.sim, nrow = n)
	zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim), 
			skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
	z = .custzdist(custom.dist, zmatrix, m.sim, n)
	z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
	if( is.matrix( custom.dist$distfit ) && custom.dist$type != "z") cres = TRUE else cres = FALSE
	# create the presample information
	if(!any(is.na(prereturns))){
		prereturns = as.vector(prereturns)
		if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
	}
	
	if(!any(is.na(preresiduals))){
		preresiduals = as.vector(preresiduals)
		preres = preresiduals
		if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
	}
	
	if(any(is.na(prereturns))){
		prereturns = rep( as.numeric(mu), m)
	}
	
	# input vectors/matrices
	x = c(prereturns, rep(0, n))
	constm = matrix(mu, ncol = m.sim, nrow = n + m)
	
	# outpus matrices
	seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
	residSim =  matrix(0, ncol = m.sim, nrow = n.sim)
	
	for(i in 1:m.sim){
		if(any(is.na(preresiduals))){
			preres = if(cres) as.numeric(z[1:m,i]) else as.numeric(z[1:m,i])*sigma
		} 
		res = c(preres, if(cres) as.numeric(z[(m+1):(n+m),i]) else as.numeric(z[(m+1):(n+m),i]) * sigma)
		residSim[,i] = tail(res, n.sim)
		

		if(include.mex){
			constm[,i] = constm[,i] + mxreg%*%t(matrix(mexsim[[i]], ncol = mxn))
		}
		if(arfima){
			#if(constant) constm[,i] = constm[,i]*(1-sum(ar))
			x = c(prereturns, rep(0, n + armaq))
			res = c(if(armaq>0) rep(0, armaq) else NULL, preres[1:m], if(cres) as.numeric(z[(m+1):(n+m),i]) else as.numeric(z[(m+1):(n+m),i]) * sigma)
			ans2 = .arfimaxsim(armap, armaq, constm[1:(n + m), i], ar, ma, darfima, x[1:(n + m + armaq)], res[1:(n + m + armaq)], T = n + m)
			seriesSim[,i] = ans2$s[(m+1):(n+m)]
		} else{
			#if(constant) constm[,i] = constm[,i]*(1-sum(ar))
			ans2 = .armaxsim(armap, armaq, constm[,i], ar, ma, x, res, T = n + m, m)
			seriesSim[,i] = tail(ans2$x, n.sim)
		}
	}
	
	path = list(seriesSim = seriesSim, residSim = residSim)
	sol = new("ARFIMApath",
			path = path,
			seed = as.integer(sseed))
	return(sol)
}

.arfimaroll = function(spec,  data, n.ahead = 1, forecast.length = 500, 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(), calculate.VaR = TRUE, VaR.alpha = c(0.01, 0.05), ...)
{

	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	dataname = names(data)
	xdata = .extractdata(data)
	data = xdata$data
	dates = xdata$pos
	if( is.null(fit.control$fixed.se) ) fit.control$fixed.se = 0
	T = length(data)
	startT = T-forecast.length
	include.dlambda = spec@distribution.model$include.dlambda
	include.skew 	= spec@distribution.model$include.skew
	include.shape 	= spec@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)
	mexdata = vector(mode="list",length = nf)
	fitlist = vector(mode="list",length = nf)
	outdata = vector(mode="list",length = nf)
	coefv = vector(mode="list", length = nf)
	distribution = spec@distribution.model$distribution
	forecastlist = vector(mode="list", length = nf)
	forecast.length = floor(forecast.length/refit.every) * refit.every
	
	for(i in 1:nf){
		if(refit.window[1]=="recursive"){
			fitdata[[i]] = data[1:(fitindx.start[i]+refit.every)]
			outdata[[i]] = data[(fitindx.start[i]+1):(fitindx.start[i]+refit.every)]
			if( spec@mean.model$include.mex ){
				mexdata[[i]] = spec@mean.model$mexdata[1:(fitindx.start[i]+refit.every), , drop = FALSE]
			} else{
				mexdata[[i]] = NULL
			}
		} else{
			fitdata[[i]] = data[(i*refit.every):(fitindx.start[i]+refit.every)]
			outdata[[i]] = data[(fitindx.start[i]+1):(fitindx.start[i]+refit.every)]
			if( spec@mean.model$include.mex ){
				mexdata[[i]] = spec@mean.model$mexdata[(i*refit.every):(fitindx.start[i]+refit.every), , drop = FALSE]
			} else{
				mexdata[[i]] = NULL
			}
		}
	}
	
	cat("\n...estimating refit windows...\n")
	specx = vector(mode = "list", length = nf)
	for(i in 1:nf){
		specx[[i]] = spec
		specx[[i]]@mean.model$mexdata = if(spec@mean.model$include.mex) mexdata[[i]] else NULL
	}
	if( parallel ){
		if( parallel.control$pkg == "multicore" ){
			if(!exists("mclapply")){
				require('multicore')
			}
			fitlist = mclapply(1:nf, FUN = function(i) arfimafit(spec = specx[[i]], data = fitdata[[i]], 
								solver = solver, out.sample = refit.every, solver.control = solver.control,
								fit.control = fit.control), mc.cores = parallel.control$cores)
		} else{
			if(!exists("sfLapply")){
				require('snowfall')
			}
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("specx", "fitdata", "refit.every", "solver", "solver.control", "fit.control",local = TRUE)
			fitlist = sfLapply(as.list(1:nf), fun = function(i) rgarch::ugarchfit(spec = specx[[i]], 
								data = fitdata[[i]], solver = solver, out.sample = refit.every, 
								solver.control = solver.control, fit.control = fit.control))
			sfStop()
		}
	} else{
		for(i in 1:nf){
			fitlist[[i]] = arfimafit(data = fitdata[[i]], spec = specx[[i]], solver = solver, 
					out.sample = refit.every, solver.control = solver.control,
					fit.control = fit.control)
			if(fitlist[[i]]@fit$convergence!=0){
				if(i>1){
					specx[[i]]@optimization.model$start.pars = as.list(coef(fitlist[[i-1]]))
					fitlist[[i]] = arfimafit(data = fitdata[[i]], spec = specx[[i]], solver = solver, 
							out.sample = refit.every, solver.control = solver.control,
							fit.control = fit.control)
					specx[[i]]@optimization.model$start.pars = NULL
				}
			}
		}
	}
	converge = sapply(fitlist, FUN=function(x) x@fit$convergence)
	if(any(converge!=0)){
		ncon = which(converge!=0)
		stop(paste("\nno convergence in the following fits: ", ncon, "...exiting", sep=""), call.=FALSE)
	}
	cat("\nDone!...all converged.\n")
	coefv = t(sapply(fitlist,FUN=function(x) x@fit$coef))
	coefmat = lapply(fitlist,FUN=function(x) x@fit$robust.matcoef)
	LLH = sapply(fitlist, FUN=function(x) x@fit$LLH)
	# filtered sigma/series(actual) estimates
	filter.sigma = vector(mode = "numeric", length = forecast.length)
	filter.series = unlist(outdata)
	
	if( parallel ){
		if( parallel.control$pkg == "multicore" ){
			forecastlist = mclapply(fitlist, FUN = function(x) arfimaforecast(x, 
								n.ahead = n.ahead, n.roll = refit.every-1), mc.cores = parallel.control$cores)
		} else{
			nx = length(fitlist)
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("fitlist", "n.ahead", "refit.every", local = TRUE)
			forecastlist = sfLapply(as.list(1:nx), fun = function(i) rgarch::arfimaforecast(fitlist[[i]], 
								n.ahead = n.ahead, n.roll = refit.every-1))
			sfStop()
		}
	} else{
		forecastlist = lapply(fitlist, FUN = function(x) arfimaforecast(x, n.ahead = n.ahead, n.roll = refit.every-1))
	}
	
	eindex = t(.embed(1:forecast.length, refit.every, refit.every, TRUE))
	# collect forecasts [mu sigma (skew shape)]
	# 2 cases : n.ahead = 1 we use a matrix, else
	# n.ahead > 1 we use a list object
	if(n.ahead == 1){
		fpmlist = vector(mode="list", length = 1)
		fser = as.numeric(sapply(forecastlist,FUN = function(x) sapply(x@forecast$forecasts, FUN=function(x) x)))
		fdates = dates[as.numeric(fwindex)]
		rdat = data[(fitindx.start[1]+1):(fitindx.start[1] + forecast.length + 1 - 1)]
		fpmlist[[1]] = .forctestsarfima3(fser, fitindx.start[1]+1, n.roll = forecast.length, 1, rdat, fdates, type = 2)
		eindex = t(.embed(1:forecast.length, refit.every, refit.every, TRUE))
		if(include.dlambda) dlambda.f = sapply(fitlist,FUN=function(x) coef(x)["dlambda"])
		if(include.skew) 	skew.f 	= sapply(fitlist,FUN=function(x) coef(x)["skew"])
		if(include.shape) 	shape.f = sapply(fitlist,FUN=function(x) coef(x)["shape"])
		sigma.f = sapply(fitlist,FUN=function(x) coef(x)["sigma"])
		fmatrix = matrix(NA, ncol = 5, nrow = forecast.length)
		colnames(fmatrix) = c("muf", "sigmaf", "dlambdaf", "skewf", "shapef")
		fmatrix[,1] = fser
		for(i in 1:dim(eindex)[2]){
			fmatrix[eindex[,i], 2] = rep(sigma.f[i], refit.every)
			fmatrix[eindex[,i], 3] = rep(if(include.dlambda) dlambda.f[i] else 0, refit.every)
			fmatrix[eindex[,i], 4] = rep(if(include.skew) skew.f[i] else 0, refit.every)
			fmatrix[eindex[,i], 5] = rep(if(include.shape) shape.f[i] else 0, refit.every)
		}
		#re - scale the fmatrix to returns-based density
		smatrix = .scaledist(dist = distribution, fmatrix[,1], fmatrix[,2], fmatrix[,3], fmatrix[,4], fmatrix[,5])
		#fdates = dates[as.numeric(fwindex)]
		fdensity = vector(mode = "list", length = 1)
		f01density = vector(mode = "list", length = 1)
		fdensity[[1]] = data.frame(fdate = fdates, fmu=smatrix[,1], fsigma=smatrix[,2], fdlambda = fmatrix[,3], fskew=smatrix[,3], fshape=smatrix[,4])
		f01density[[1]] = data.frame(f01date=fdates, f01mu=fmatrix[,1], f01sigma=fmatrix[,2], f01dlambda = fmatrix[,3], f01skew=fmatrix[,3], f01shape=fmatrix[,4])
		if(calculate.VaR){
			VaR.list = vector(mode="list", length = 1)
			cat("\n...calculating VaR...\n")
			if(is.null(VaR.alpha)) VaR.alpha = c(0.01, 0.05)
			n.v = dim(fdensity[[1]])[1]
			m.v = length(VaR.alpha)
			VaR.matrix = matrix(NA, ncol = m.v + 1, nrow = n.v)
			for(i in 1:m.v){
				VaR.matrix[,i] = .qdensity(rep(VaR.alpha[i], n.v) , mu = smatrix[,1], sigma = smatrix[,2], lambda = fmatrix[,3], 
						skew = smatrix[,3], shape = smatrix[,4], distribution = distribution)
			}
			VaR.matrix[,m.v+1] = unlist(outdata)
			colnames(VaR.matrix) = c(paste("alpha(", round(VaR.alpha,2)*100, "%)",sep=""), "actual")
			rownames(VaR.matrix) = as.character(fdates)
			VaR.list[[1]] = VaR.matrix
			cat("\nDone!\n")
		} else{
			VaR.matrix = NULL
			VaR.alpha = NULL
		}
		rolllist = list()
		rolllist$n.refit = nf
		rolllist$refit.every = refit.every
		rolllist$distribution = distribution
		rolllist$fdensity = fdensity
		rolllist$f01density = f01density
		rolllist$coefs = coefv
		rolllist$coefmat = coefmat
		rolllist$LLH = LLH
		rolllist$VaR.out = VaR.list
		rolllist$n.ahead = n.ahead
		rolllist$fpm = fpmlist
		rolllist$filterseries = filter.series
		rolllist$forecast.length = forecast.length
		rolllist$VaR.alpha =VaR.alpha
		rolllist$dataname = dataname
		ans = new("ARFIMAroll",
				roll = rolllist,
				forecast = forecastlist,
				spec)
	} else{
		fpmlist = vector(mode="list", length = nf)
		eindex = t(.embed(1:forecast.length, refit.every, refit.every, TRUE))
		sigma.f = sapply(fitlist,FUN=function(x) coef(x)["sigma"])
		if(include.dlambda) dlambda.f = sapply(fitlist,FUN=function(x) coef(x)["dlambda"])
		if(include.skew) skew.f = sapply(fitlist,FUN=function(x) coef(x)["skew"])
		if(include.shape) shape.f = sapply(fitlist,FUN=function(x) coef(x)["shape"])
		# split the array into 1:n-day ahead forecasts
		flist = vector(mode="list", length = n.ahead)
		f01density = vector(mode="list", length = n.ahead)
		fdensity = vector(mode="list", length = n.ahead)
		#(NA, dim = c(forecast.length, 4, n.ahead))
		for(i in 1:n.ahead){
			fser =  unlist(lapply(forecastlist, FUN=function(x) sapply(x@forecast$forecast, FUN=function(x) x[i])))
			rdat = data[(fitindx.start[1]+i):(fitindx.start[1] + forecast.length + i - 1)]
			fdates = as.character(dates[(fitindx.start[1]+i):(fitindx.start[1] + forecast.length + i - 1)])
			fpmlist[[i]] = .forctestsarfima3(fser, fitindx.start[1]+i, n.roll = forecast.length, i, rdat, fdates, type = 1)
			flist[[i]] = fdensity[[i]] = f01density[[i]] = matrix(NA, ncol = 5, nrow = forecast.length)
			f01density[[i]][,1] = fser
			for(j in 1:dim(eindex)[2]){
				f01density[[i]][eindex[,j],2] = rep(sigma.f[j], refit.every)
				f01density[[i]][eindex[,j],3] = rep(if(include.dlambda) dlambda.f[j] else 0, refit.every)
				f01density[[i]][eindex[,j],4] = rep(if(include.skew) skew.f[j] else 0, refit.every)
				f01density[[i]][eindex[,j],5] = rep(if(include.shape) shape.f[j] else 0, refit.every)
			}
			fdensity[[i]] = .scaledist(dist = distribution, f01density[[i]][,1], f01density[[i]][,2], f01density[[i]][,3], f01density[[i]][,4], f01density[[i]][,5])
			fdensity[[i]] = as.data.frame(fdensity[[i]])
			f01density[[i]] = as.data.frame(f01density[[i]])
			
			fdensity[[i]] = cbind(fdates, fdensity[[i]])
			
			f01density[[i]] = cbind(fdates, f01density[[i]])
			
			colnames(fdensity[[i]])  = c("fdate", "fmu", "fsigma", "fskew", "fshape")
			colnames(f01density[[i]]) = c("f01date", "f01mu", "f01sigma", "f01dlambda", "f01skew", "f01shape")
		}
		if(calculate.VaR){
			# should consider implementing mclapply here as well
			cat("\n...calculating VaR...\n")
			if(is.null(VaR.alpha)) VaR.alpha = c(0.01, 0.05)
			n.v = forecast.length
			m.v = length(VaR.alpha)
			VaR.list = vector(mode="list", length = n.ahead)
			for(j in 1:n.ahead){
				VaR.list[[j]] = matrix(NA, ncol = m.v+1, nrow = n.v)
				for(i in 1:m.v){
					VaR.list[[j]][,i] = .qdensity(rep(VaR.alpha[i], n.v) , mu = fdensity[[j]][,2],
							sigma = fdensity[[j]][,3], lambda = f01density[[i]][,3], skew = fdensity[[j]][,4], 
							shape = fdensity[[j]][,5],distribution = distribution)
				}
				VaR.list[[j]][,m.v+1] = data[(fitindx.start[1]+j):(fitindx.start[1] + forecast.length + j - 1)]
				colnames(VaR.list[[j]]) = c(paste("alpha(", round(VaR.alpha,2)*100, "%)",sep=""), "actual")
				rownames(VaR.list[[j]]) = as.character(dates[(fitindx.start[1]+j):(fitindx.start[1] + forecast.length + j - 1)])
			}
			cat("\nDone!\n")
		} else{
			VaR.list = NULL
			VaR.alpha = NULL
		}
		rolllist = list()
		rolllist$n.refit = nf
		rolllist$fpm  = fpmlist
		rolllist$refit.every = refit.every
		rolllist$distribution = distribution
		rolllist$fdensity = fdensity
		rolllist$f01density = f01density
		rolllist$coefs = coefv
		rolllist$coefmat = coefmat
		rolllist$LLH = LLH
		rolllist$VaR.out = VaR.list
		rolllist$n.ahead = n.ahead
		rolllist$filterseries = filter.series
		rolllist$forecast.length = forecast.length
		rolllist$VaR.alpha =VaR.alpha
		rolllist$dataname = dataname
		ans = new("ARFIMAroll",
				roll = rolllist,
				forecast = forecastlist,
				spec)
	}
	return(ans)
}

.arfimadistribution = function(fitORspec, data = NULL, n.sim = 2000, n.start = 1, 
		m.sim = 100, recursive = FALSE, recursive.length = 6000, recursive.window = 1000, 
		prereturns = NA, preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA, type = "z"), 
		mexsimdata = NULL, fit.control = list(), solver = "solnp", 
		solver.control = list(), parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	if(recursive){
		nwindows = 1 + round( (recursive.length - n.sim) / recursive.window )
		swindow = vector(mode = "list", length = nwindows)
		rwindow = vector(mode = "list", length = nwindows)
	} else{
		nwindows = 1
		swindow = vector(mode = "list", length = nwindows)
		rwindow = vector(mode = "list", length = nwindows)
		recursive.window = 0
	}
	if(is(fitORspec, "ARFIMAfit")){
		
		for(i in 1:nwindows){
			sim = arfimasim(fitORspec, n.sim = n.sim + (i-1)*recursive.window, 
					n.start = n.start, m.sim = m.sim, prereturns = prereturns, 
					preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist, 
					mexsimdata = mexsimdata)
			swindow[[i]]$path.df = as.data.frame(sim)
			swindow[[i]]$seed = sim@seed
		}
		fixpars = as.list(coef(fitORspec))
		truecoef = fitORspec@fit$robust.matcoef
		spec = getspec(fitORspec)
	}
	# simulate series paths
	if(is(fitORspec, "ARFIMAspec")){
		
		for(i in 1:nwindows){
			sim  = arfimapath(fitORspec, n.sim =  n.sim + (i-1)*recursive.window, 
					n.start = n.start, m.sim = m.sim, prereturns = prereturns, 
					preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist, 
					mexsimdata = mexsimdata)
			swindow[[i]]$path.df = as.data.frame(sim)
			swindow[[i]]$seed = sim@seed
			
		}
		spec = fitORspec
		spec@optimization.model$fixed.pars = NULL
		fixpars = fitORspec@optimization.model$fixed.pars
		truecoef = as.matrix(cbind(unlist(fitORspec@optimization.model$fixed.pars),rep(0,length(fixpars)),
						rep(10, length(fixpars)),rep(0,length(fixpars))))
	}
	# fit to simulated series (using starting parameters)
	spec@optimization.model$start.pars = fixpars
	
	fitlist = vector( mode = "list", length = m.sim )
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			if(!exists("mclapply")){
				require('multicore')
			}
			for(i in 1:nwindows){
				rwindow[[i]]$fitlist = mclapply(swindow[[i]]$path.df, FUN = function(x) .fitandextractarfima(spec, x, out.sample = 0, 
									solver = solver, fit.control = fit.control, solver.control = solver.control), mc.cores = parallel.control$cores)
			}
		} else{
			for(i in 1:nwindows){
				nx = dim(swindow[[i]]$path.df)[2]
				sfInit(parallel = TRUE, cpus = parallel.control$cores)
				sfExport("spec", "swindow", "solver", "fit.control", "solver.control", local = TRUE)
				rwindow[[i]]$fitlist = sfLapply(as.list(1:nx), fun = function(j) rgarch:::.fitandextractarfima(spec, swindow[[j]]$path.df, out.sample = 0, 
									solver = solver, fit.control = fit.control, solver.control = solver.control))
				sfStop()
			}
		}
	} else{
		for(i in 1:nwindows){
			rwindow[[i]]$fitlist = lapply(swindow[[i]]$path.df, FUN = function(x) .fitandextractarfima(spec, x, out.sample = 0, 
								solver = solver, fit.control = fit.control, solver.control = solver.control))
		}
	}
	
	reslist = vector(mode = "list", length = nwindows)
	for(j in 1:nwindows){
		reslist[[j]]$simcoef = 	matrix(NA, ncol = length(fixpars), nrow = m.sim)
		reslist[[j]]$rmse = 	rep(NA, length = length(fixpars))
		reslist[[j]]$simcoefse = matrix(NA, ncol = length(fixpars), nrow = m.sim)
		reslist[[j]]$likelist = rep(NA, length = m.sim)
		reslist[[j]]$mlongrun = rep(NA, length = m.sim)
		reslist[[j]]$simmaxdata  = matrix(NA, ncol = 2, nrow = m.sim)
		reslist[[j]]$simmindata  = matrix(NA, ncol = 2, nrow = m.sim)
		reslist[[j]]$simmeandata  = matrix(NA, ncol = 2, nrow = m.sim)
		reslist[[j]]$simmomdata = matrix(NA, ncol = 2, nrow = m.sim)
		reslist[[j]]$convergence = 	rep(1, length = m.sim)
		reslist[[j]]$seeds = rep(1, length = m.sim)
		for(i in 1:m.sim){
			if(rwindow[[j]]$fitlist[[i]]$convergence!=0) next()
			reslist[[j]]$simcoef[i, ] = rwindow[[j]]$fitlist[[i]]$simcoef
			reslist[[j]]$simcoefse[i,] = rwindow[[j]]$fitlist[[i]]$simcoefse
			reslist[[j]]$likelist[i] = 	rwindow[[j]]$fitlist[[i]]$llh
			reslist[[j]]$mlongrun[i] = 	rwindow[[j]]$fitlist[[i]]$mlongrun
			reslist[[j]]$simmaxdata[i, ] = rwindow[[j]]$fitlist[[i]]$maxdata
			reslist[[j]]$simmindata[i, ] = rwindow[[j]]$fitlist[[i]]$mindata
			reslist[[j]]$simmeandata[i, ] = rwindow[[j]]$fitlist[[i]]$meandata
			reslist[[j]]$simmomdata[i, ] = rwindow[[j]]$fitlist[[i]]$momdata
			reslist[[j]]$convergence[i] = rwindow[[j]]$fitlist[[i]]$convergence
		}
		reslist[[j]]$seed = swindow[[j]]$seed
		reslist[[j]]$rmse = .rmse(reslist[[j]]$simcoef, unlist(fixpars))
	}
	reslist$details = list(n.sim = n.sim, n.start = n.start, m.sim = m.sim,  
			recursive = recursive, recursive.length = recursive.length, recursive.window = recursive.window,
			nwindows = nwindows)
	ans = new("ARFIMAdistribution",
			dist = reslist,
			truecoef = truecoef,
			spec = spec)
	
	return(ans)
}

.fitandextractarfima = function(spec, x, out.sample = 0,  solver = "solnp", fit.control = list(), solver.control = list())
{
	dist = list()
	fit = .safefitarfima(spec, x, out.sample = 0, solver = solver, fit.control = fit.control, solver.control = solver.control)
	if( is.null(fit) || fit@fit$convergence == 1 || !is( fit, "ARFIMAfit" ) || any( is.na( coef( fit ) ) ) ){
		dist$convergence = 1
		return(dist)
	}
	dist = list()
	dist$simcoef = coef(fit)
	dist$simcoefse = fit@fit$robust.matcoef[, 2]
	dist$llh = likelihood(fit)
	dist$mlongrun = uncmean(fit)
	tmp = as.data.frame(fit)
	dist$maxdata = apply(tmp[, -2], 2, "max")
	dist$mindata = apply(tmp[, -2], 2, "min")
	dist$meandata = apply(tmp[, -2], 2, "mean")
	dist$momdata = c(.kurtosis(tmp[,1]), .skewness(tmp[,1]))
	dist$convergence = fit@fit$convergence
	return(dist)
}

Try the rgarch package in your browser

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

rgarch documentation built on May 31, 2017, 3:20 a.m.