Nothing
R/rugarch-figarch.R
In rugarch: Univariate GARCH Models
Defines functions
.figarchpath1
.figarchpath
.figarchsim1
.figarchsim
.figarchforecast2
.nfigarchforecast
.figarchforecast
.figarchfilter
.figarchLLH
.figarchfit
#################################################################################
##
## R package rugarch by Alexios Galanos Copyright (C) 2008-2022.
## This file is part of the R package rugarch.
##
## The R package rugarch 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 rugarch 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 fiGARCH fit
#---------------------------------------------------------------------------------
# [mu ar ma arfima im mxreg omega alpha beta gamma gamma11 gamma21 delta lambda vxreg skew shape dlamda aux aux aux aux]
.figarchfit = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(stationarity = 1, fixed.se = 0, scale = 0, rec.init = 'all',trunclag = 1000),
numderiv.control = list(grad.eps=1e-4, grad.d=0.0001, grad.zero.tol=sqrt(.Machine$double.eps/7e-7),
hess.eps=1e-4, hess.d=0.1, hess.zero.tol=sqrt(.Machine$double.eps/7e-7), r=4, v=2))
{
tic = Sys.time()
if(is.null(solver.control$trace)) trace = 0 else trace = solver.control$trace
if(is.null(fit.control$stationarity)) fit.control$stationarity = TRUE
if(is.null(fit.control$fixed.se)) fit.control$fixed.se = FALSE
if(is.null(fit.control$scale)) fit.control$scale = FALSE
if(is.null(fit.control$rec.init)) fit.control$rec.init = 'all'
if(is.null(fit.control$trunclag)) truncLag = 1000 else truncLag = fit.control$trunclag
mm = match(names(fit.control), c("stationarity", "fixed.se", "scale", "rec.init", "trunclag"))
if(any(is.na(mm))){
idx = which(is.na(mm))
enx = NULL
for(i in 1:length(idx)) enx = c(enx, names(fit.control)[idx[i]])
warning(paste(c("unidentified option(s) in fit.control:\n", enx), sep="", collapse=" "), call. = FALSE, domain = NULL)
}
# if we have external regressors in variance turn off scaling
if(spec@model$modelinc[15] > 0) fit.control$scale = FALSE
# if we have arch-in-mean turn off scaling
if(spec@model$modelinc[5] > 0) fit.control$scale = FALSE
# if there are fixed pars we do no allow scaling as there would be no way of mixing scaled
# amd non scaled parameters
if(sum(spec@model$pars[,2]) > 0) fit.control$scale = FALSE
xdata = .extractdata(data)
if(!is.numeric(out.sample)) stop("\nugarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample)<0) stop("\nugarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample,0)
n = length(xdata$data)
#if((n-n.start)<100) stop("\nugarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start)]
index = xdata$index[1:(n-n.start)]
origdata = xdata$data
origindex = xdata$index
period = xdata$period
# arglist replaces the use of custom environments (resolves the problem of
# non-shared environment in parallel estimation, particularly windows)
garchenv = new.env(hash = TRUE)
arglist = list()
arglist$garchenv <- garchenv
arglist$pmode = 0
arglist$truncLag = truncLag
model = spec@model
modelinc = model$modelinc
pidx = model$pidx
# expand the spec object and assign spec lists
if(modelinc[6] > 0){
mexdata = model$modeldata$mexdata[1:(n-n.start), , drop = FALSE]
} else{
mexdata = NULL
}
if(modelinc[15] > 0){
vexdata = model$modeldata$vexdata[1:(n-n.start), ,drop = FALSE]
} else{
vexdata = NULL
}
arglist$index = index
arglist$trace = trace
m = model$maxOrder
# store length of data for easy retrieval
model$modeldata$T = T = length(as.numeric(data))
dist = model$modeldesc$distribution
if(fit.control$scale) dscale = sd(data) else dscale = 1
zdata = data/dscale
recinit = .checkrec(fit.control$rec.init, T)
arglist$data = zdata
arglist$recinit = recinit
arglist$dscale = dscale
arglist$model = model
ipars = model$pars
# Optimization Starting Parameters Vector & Bounds
tmp = .garchstart(ipars, arglist)
####################################################
# Scale Omega in order to avoid problems: (25-01-2013)
#arglist$omegascale = 1/1e8
#tmp$ipars["omega",c(1,5,6)] = tmp$ipars["omega",c(1,5,6)] *1e8
####################################################
ipars = arglist$ipars = tmp$pars
arglist$tmph = tmp$tmph
# we now split out any fixed parameters
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
arglist$fit.control = fit.control
npars = sum(estidx)
if(any(ipars[,2]==1)){
if(npars == 0){
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
warning("\nugarchfit-->warning: all parameters fixed...returning ugarchfilter object instead\n")
return(ugarchfilter(data = xts(origdata, origindex), spec = spec, out.sample = out.sample))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = 0
ipars[ipars[,2]==1, 4] = 1
ipars[ipars[,2]==1, 2] = 0
arglist$ipars = ipars
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = 1
}
} else{
use.solver = 1
}
# start counter
assign("rugarch_llh", 1, envir = garchenv)
# assign solver constraints (solnp directly else exterior type penalty
# for other solvers)
if(fit.control$stationarity == 1 && modelinc[15] == 0){
Ifn = .figarchcon
ILB = rep(0,2)
IUB = rep(1,2)
if(solver == "solnp" | solver == "gosolnp" | solver == "hybrid") fit.control$stationarity = 0
} else{
Ifn = NULL
ILB = NULL
IUB = NULL
}
# conditions controls the non-solnp solver penalty
arglist$fit.control = fit.control
if(use.solver){
parscale = rep(1, length = npars)
names(parscale) = rownames(ipars[estidx,])
if(modelinc[1] > 0) parscale["mu"] = abs(mean(zdata))
if(modelinc[7] > 0) parscale["omega"] = var(zdata)
arglist$returnType = "llh"
solution = .garchsolver(solver, pars = ipars[estidx, 1], fun = .figarchLLH,
Ifn, ILB, IUB, gr = NULL, hessian = NULL, parscale = parscale,
control = solver.control, LB = ipars[estidx, 5],
UB = ipars[estidx, 6], ux = NULL, ci = NULL, mu = NULL, arglist)
sol = solution$sol
hess = solution$hess
timer = Sys.time()-tic
if(!is.null(sol$par)){
ipars[estidx, 1] = sol$par
if(modelinc[7]==0){
# call it once more to get omega
tmpx = .figarchLLH(sol$par, arglist)
ipars[pidx["omega",1], 1] = get("omega", garchenv)
}
if(sum(ipars[,2]) == 0){
if(modelinc[1] > 0) ipars[pidx["mu",1]:pidx["mu",2], 1] = ipars[pidx["mu",1]:pidx["mu",2], 1] * dscale
if(modelinc[6] > 0){
ipars[pidx["mxreg", 1]:pidx["mxreg", 2], 1] = ipars[pidx["mxreg", 1]:pidx["mxreg", 2], 1] * dscale
}
ipars[pidx["omega",1], 1] = ipars[pidx["omega",1],1] * dscale^2
}
} else{
ipars[estidx, 1] = NA
}
arglist$ipars = ipars
convergence = sol$convergence
if(convergence != 0) warning("\nugarchfit-->warning: solver failer to converge.")
} else{
solution = NULL
hess = NULL
timer = Sys.time()-tic
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
# check convergence else write message/return
# create a copy of ipars in case we need to change it below to calculate standard errors
# which we will need to reset later (because for example, infocriteria uses estimated
# parameters, not fixed.
ipars2 = ipars
if(convergence == 0){
arglist$dscale = 1
if(sum(ipars[,2]) > 0 && fit.control$fixed.se == 1){
ipars[ipars[,2]==1, 4] = 1
ipars[ipars[,2]==1, 2] = 0
arglist$ipars = ipars
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
}
arglist$data = data
fit = .makefitmodel(garchmodel = "fiGARCH", f = .figarchLLH, T = T, m = m,
timer = timer, convergence = convergence, message = sol$message,
hess = hess, arglist = arglist, numderiv.control = numderiv.control)
model$modelinc[7] = modelinc[7]
model$modeldata$data = origdata
model$modeldata$index = origindex
model$modeldata$period = period
model$pars[, 1] = fit$ipars[,1]
model$pars[, 5:6] = ipars2[,5:6]
fit$ipars[, 4] = ipars2[, 4]
fit$ipars[, 2] = ipars2[, 2]
fit$ipars[, 5:6] = ipars2[,5:6]
# make sure omega is now included (for working with object post-estimation)
fit$ipars["omega", 3] = 1
model$pars["omega", 3] = 1
} else{
fit$message = sol$message
fit$convergence = 1
model$modeldata$data = origdata
model$modeldata$index = origindex
model$modeldata$period = period
}
# make model list to return some usefule information which
# will be called by other functions (show, plot, sim etc)
model$n.start = n.start
fit$solver = solution
model$trunclag = truncLag
ans = new("uGARCHfit",
fit = fit,
model = model)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION fiGARCH LLH
#---------------------------------------------------------------------------------
.figarchLLH = function(pars, arglist)
{
# prepare inputs
data = arglist$data
returnType = arglist$returnType
garchenv = arglist$garchenv
# rejoin fixed and pars
model = arglist$model
estidx = arglist$estidx
idx = model$pidx
ipars = arglist$ipars
ipars[estidx, 1] = pars
trace = arglist$trace
T = length(data)
dscale = arglist$dscale
recinit = arglist$recinit
fit.control = arglist$fit.control
m = model$maxOrder
N = c(m,T)
truncLag = arglist$truncLag
mexdata = coredata(model$modeldata$mexdata[1:T,, drop = FALSE])
vexdata = coredata(model$modeldata$vexdata[1:T,, drop = FALSE])
distribution = model$modeldesc$distribution
modelinc = model$modelinc
# persist used equivalently here for the positivity constraints
persist = .figarchcon(pars, arglist)
# distribution number
# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
dist = model$modeldesc$distno
hm = arglist$tmph
rx = .arfimaxfilter(modelinc[1:21], pars = ipars[,1], idx = idx, mexdata = mexdata, h = hm, data = data, N = N)
res = rx$res
zrf = rx$zrf
res[is.na(res) | !is.finite(res) | is.nan(res)] = 0
# recursion initialization
mvar = ifelse(recinit$type==1, mean(res[1:recinit$n]*res[1:recinit$n]), backcastv(res, T, recinit$n))
if(modelinc[15]>0) {
mv = sum(apply(matrix(vexdata, ncol = modelinc[15]), 2, "mean")*ipars[idx["vxreg",1]:idx["vxreg",2],1])
} else{
mv = 0
}
hEst = mvar
# variance targeting not used
ipars[idx["omega",1],1] = max(eps, ipars[idx["omega",1],1])
if(is.na(hEst) | !is.finite(hEst) | is.nan(hEst)) hEst = var(data)
delta=ipars[idx["delta",1],1]
# figarch setup
NL=T+truncLag
eps=matrix(mvar, ncol=1, nrow=NL)
ebar = rep(0, T)
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
# if we have external regressors in variance equation we cannot have
# stationarity checks in likelihood. Stationarity conditions not valid for
# solnp solver which implements them internally as constraints.
if(fit.control$stationarity == 1 && modelinc[15] == 0){
if(any(is.na(persist)) || any(persist >= 1) || any(persist <= 0)){
if(arglist$pmode!=1){
return(llh = get("rugarch_llh", garchenv) + 0.1*(abs(get("rugarch_llh", garchenv))))
} else{
return(llh = 1e10)
}
}
}
ebar[1]=t(be)%*%eps[1:truncLag]
if(modelinc[6]>0) mexdata = as.double(as.vector(mexdata)) else mexdata = double(1)
if(modelinc[15]>0) vexdata = as.double(as.vector(vexdata)) else vexdata = double(1)
# modelinc (1:21) since 22 is either NA or numeric and no longer needed (paased to ipars if used)
ans = try( .C("figarchfilterC", model = as.integer(modelinc[1:21]),
pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
hEst = as.double(hEst), x = as.double(data), res = as.double(res),
e = double(T), ebar = as.double(ebar), eps=as.double(eps), be=as.double(be),
mexdata = mexdata, vexdata = vexdata,
zrf = as.double(zrf), constm = double(T), condm = double(T),
m = as.integer(m), T = as.integer(T), N = as.integer(truncLag),
h = double(T), z = double(T), llh = double(1), LHT = double(T),
PACKAGE = "rugarch"), silent = TRUE )
if(inherits(ans, "try-error")){
if(arglist$pmode!=1){
assign("rugarch_csol", 1, envir = garchenv)
assign("rugarch_filtermessage", ans, envir = garchenv)
if( trace > 0 ) cat(paste("\narfimafit-->warning: ", get("rugarch_filtermessage", garchenv),"\n", sep=""))
return(llh = (get("rugarch_llh", garchenv) + 0.1*(abs(get("rugarch_llh", garchenv)))))
} else{
return(llh = 1e10)
}
} else{
if(arglist$pmode!=1){
assign("rugarch_csol", 0, envir = garchenv)
}
}
z = ans$z
h = ans$h
epsx = ans$res
eps = ans$eps
ebar = ans$ebar
llh = ans$llh
if(is.finite(llh) && !is.na(llh) && !is.nan(llh)){
if(arglist$pmode!=1) assign("rugarch_llh", llh, envir = garchenv)
} else {
if(arglist$pmode!=1) llh = (get("rugarch_llh", garchenv) + 0.1*(abs(get("rugarch_llh", garchenv)))) else llh = 1e10
}
# LHT = raw scores
# ? -ans$LHT[(m+1):T]
LHT = -ans$LHT
ans = switch(returnType,
llh = llh,
LHT = LHT,
all = list(llh = llh, h = h, epsx = epsx, z = z, eps=eps, ebar=ebar, be=be, kappa = kappa, LHT = LHT, persistence = persist))
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION sGARCH filter
#---------------------------------------------------------------------------------
.figarchfilter = function(spec, data, out.sample = 0, n.old = NULL, rec.init = 'all', trunclag = 1000)
{
# 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.
xdata = .extractdata(data)
data = xdata$data
index = xdata$index
period = xdata$period
origdata = data
origindex = index
T = length(origdata) - out.sample
data = origdata[1:T]
index = origindex[1:T]
if(!is.null(n.old)) Nx = n.old else Nx = length(data)
recinit = .checkrec(rec.init, Nx)
model = spec@model
ipars = model$pars
pars = unlist(model$fixed.pars)
parnames = names(pars)
modelnames = .checkallfixed(spec)
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
cat("\nugarchfilter-->error: parameters names do not match specification\n")
cat("Expected Parameters are: ")
cat(paste(modelnames))
cat("\n")
stop("Exiting", call. = FALSE)
}
# once more into the spec
# NB Any changes made to the spec are not preserved once we apply set fixed
setfixed(spec)<-as.list(pars)
model = spec@model
model$modeldata$T = T
ipars = model$pars
idx = model$pidx
modelinc = model$modelinc
m = model$maxOrder
N = c(m,T)
mexdata = model$modeldata$mexdata[1:T, , drop = FALSE]
vexdata = model$modeldata$vexdata[1:T, , drop = FALSE]
distribution = model$modeldesc$distribution
dist = model$modeldesc$distno
kappa = 1
#persist = (sum(ipars[idx["alpha",1]:idx["alpha",2],1]) + sum(ipars[idx["beta",1]:idx["beta",2],1]))
rx = .arfimaxfilter(modelinc[1:21], ipars[,1], idx, mexdata = mexdata, h = 0, data = data, N = N)
res = rx$res
zrf = rx$zrf
if(!is.null(n.old)){
rx2 = .arfimaxfilter(modelinc[1:21], ipars[,1], idx, mexdata = mexdata[1:Nx, , drop = FALSE], h = 0, data = origdata[1:Nx], N = c(m, Nx))
res2 = rx2$res
mvar = ifelse(recinit$type==1, mean(res2[1:recinit$n]*res2[1:recinit$n]), backcastv(res2, Nx, recinit$n))
} else{
mvar = ifelse(recinit$type==1, mean(res[1:recinit$n]*res[1:recinit$n]), backcastv(res, T, recinit$n))
}
hEst = mvar
if(modelinc[15]>0) {
mv = sum(apply(matrix(vexdata, ncol = modelinc[15]), 2, "mean")*ipars[idx["vxreg",1]:idx["vxreg",2],1])
} else{
mv = 0
}
ipars[idx["omega",1],1] = max(eps, ipars[idx["omega",1],1])
if(is.na(hEst) | !is.finite(hEst) | is.nan(hEst)) hEst = var(data)
delta=ipars[idx["delta",1],1]
# figarch setup
truncLag=trunclag
NL=T+truncLag
eps=matrix(mvar, ncol=1, nrow=NL)
ebar = rep(0, T)
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
ebar[1]=t(be)%*%eps[1:truncLag]
arglist=list()
arglist$ipars=ipars
arglist$truncLag = truncLag
persist=.figarchcon(pars, arglist)
if(any(persist<=0) | any(persist>=1)) warning("\nPositivity Contraints NOT met (check parameters)!")
if(modelinc[6]>0) mexdata = as.double(as.vector(mexdata)) else mexdata = double(1)
if(modelinc[15]>0) vexdata = as.double(as.vector(vexdata)) else vexdata = double(1)
ans = try( .C("figarchfilterC", model = as.integer(modelinc[1:21]),
pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
hEst = as.double(hEst), x = as.double(data), res = as.double(res),
e = double(T), ebar = as.double(ebar), eps=as.double(eps), be=as.double(be),
mexdata = mexdata, vexdata = vexdata,
zrf = as.double(zrf), constm = double(T), condm = double(T),
m = as.integer(m), T = as.integer(T), N = as.integer(truncLag),
h = double(T), z = double(T), llh = double(1), LHT = double(T),
PACKAGE = "rugarch"), silent = TRUE )
filter = list()
filter$z = ans$z
filter$ebar = ans$ebar
filter$eps = ans$eps
filter$be = ans$be
filter$sigma = sqrt(ans$h)
filter$residuals = ans$res
filter$LLH = -ans$llh
filter$log.likelihoods = ans$LHT
filter$persistence = persist
filter$distribution = distribution
filter$ipars = ipars
model$modeldata$data = origdata
model$modeldata$index = origindex
model$modeldata$period = period
model$n.start = out.sample
model$trunclag = trunclag
sol = new("uGARCHfilter",
filter = filter,
model = model)
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION fiGARCH forecast
#---------------------------------------------------------------------------------
# Recipe for a GARCH forecast:
# 1. split forecast into arma and garch forecasts.
# 2. take into account garch-in-mean, arfima, exogenous regressors
# 3. check kappa and persistence
# 4. allow for rolling forecast based on filtering out.sample data from fit stage
# n.roll signifies whether to filter/roll the sigma else will use unconditional
# i.e. n.ahead=10 with n.roll=1 means that the first forecast is based on the
# previous value whereas all subsequent forecasts are based on the unconditional
# expectation formula
# with n.roll = n means that while n.ahead < n.roll we filter the data using the coef
# of the garch model to obtain filtered sigma estimates which are used to roll
# the forecast. This requires that during the fit process the out.sample option
# was used and that (n.roll) < out.sample (otherwise we revert to the unconditional
# expectation formula for the long run sigma.
.figarchforecast = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
fit = fitORspec
data = fit@model$modeldata$data
Nor = length(as.numeric(data))
index = fit@model$modeldata$index
period = fit@model$modeldata$period
ns = fit@model$n.start
N = Nor - ns
model = fit@model
ipars = fit@fit$ipars
modelinc = model$modelinc
idx = model$pidx
if( n.roll > ns ) stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
pars = fit@fit$coef
ipars = fit@fit$ipars
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(model, external.forecasts$mregfor, external.forecasts$vregfor, n.ahead, Nor, out.sample = ns, n.roll)
mxf = xreg$mxf
vxf = xreg$vxf
# filter data (check external regressor data - must equal length of origData)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "fiGARCH",
garchOrder = c(modelinc[8], modelinc[9]), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(modelinc[2], modelinc[3]),
include.mean = modelinc[1],
archm = ifelse(modelinc[5]>0,TRUE,FALSE), archpow = modelinc[5], arfima = modelinc[4],
external.regressors = mxf[1:(N + fcreq), , drop = FALSE], archex = modelinc[20]),
distribution.model = model$modeldesc$distribution, fixed.pars = as.list(pars))
tmp = xts(data[1:(N + fcreq)], index[1:(N + fcreq)])
flt = .figarchfilter(data = tmp, spec = fspec, n.old = N, trunclag = fit@model$trunclag)
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
truncLag = fit@model$trunclag
# ebar,eps and be from the filtered is now available
ebar =c(flt@filter$ebar,rep(0,n.ahead))
eps = flt@filter$eps
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(eps[-c(1:truncLag)], rep(0,n.ahead))
# forecast GARCH process
seriesFor = sigmaFor = matrix(NA, ncol = n.roll+1, nrow = n.ahead)
#seriesFor[1,] = fitted(flt)[(N+1):(N+n.roll+1)]
#sigmaFor[1,] = sigma(flt)[(N+1):(N+n.roll+1)]
# n.ahead x n.roll +1 (n.ahead=1 generted by model)
colnames(seriesFor) = colnames(sigmaFor) = as.character(index[N:(N+n.roll)])
rownames(seriesFor) = rownames(sigmaFor) = paste("T+", 1:n.ahead, sep="")
for(i in 1:(n.roll+1)){
np = N + i - 1
if(modelinc[1] > 0){
mu = rep(ipars[idx["mu",1]:idx["mu",2], 1], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
omega = rep(ipars[idx["omega",1]:idx["omega",2], 1], N+i+n.ahead-1)
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(data[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]
vxfi = vxf[1:(N+i-1+n.ahead), , drop = FALSE]
# CONTINUE Here
# need to append ebar/eps/ and indicate nlagbin
ans = .nfigarchforecast(ipars, modelinc, idx, mu, omega, mxfi, vxfi, h, epsx, z, eps=eps, ebar=ebar, be=be, truncLag=truncLag, data = x, N = np, n.ahead)
sigmaFor[,i] = ans$h
seriesFor[,i] = ans$x
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigmaFor = sigmaFor
fcst$seriesFor = seriesFor
model$modeldata$sigma = flt@filter$sigma
model$modeldata$residuals = flt@filter$residuals
ans = new("uGARCHforecast",
forecast = fcst,
model = model)
return(ans)
}
.nfigarchforecast = function(ipars, modelinc, idx, mu, omega, mxfi, vxfi, h, epsx, z, eps, ebar, be, truncLag, data, N, n.ahead)
{
if(modelinc[15]>0){
omega = omega + vxfi%*%t(matrix(ipars[idx["vxreg",1]:idx["vxreg",2],1], ncol = modelinc[15]))
}
for(i in 1:n.ahead){
ebar[N+i] = be%*%eps[(N-truncLag+i):(N+i-1)]
h[N+i] = omega[N+i] - ebar[N+i]
if(modelinc[9]>=i){
h[N+i] = h[N+i] + sum(ipars[idx["beta",1]:idx["beta",2],1]*(h[N+i-(1:modelinc[9])]^2-eps[(N+i-(1:modelinc[9]))]))
}
if(modelinc[8]>0){
for (j in 1:modelinc[8]){
if ((i-j) > 0){
s = h[N + i - j]^2
} else{
s = epsx[N + i - j]^2
}
h[N+i] = h[N+i] + ipars[idx["alpha",1]+j-1,1] * (s + ebar[N + i - j])
}
}
eps[N+i] = h[N+i]
h[N+i] = sqrt(h[N+i])
}
if(modelinc[4]>0){
res = arfimaf(ipars, modelinc[1:21], idx, mu, mxfi, h, epsx, z, data, N, n.ahead)
} else{
res = armaf(ipars, modelinc[1:21], idx, mu, mxfi, h, epsx, z, data, N, n.ahead)
}
return(list(h = h[(N+1):(N+n.ahead)], x = res[(N+1):(N+n.ahead)], ebar=ebar[(N+1):(N+n.ahead)], eps=eps[(N+1):(N+n.ahead)]))
}
#---------------------------------------------------------------------------------
# 2nd dispatch method for forecast
.figarchforecast2 = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), trunclag = 1000, ...)
{
# 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
N = length(as.numeric(data))
index = xdata$index
period = xdata$period
ns = out.sample
if( n.roll > ns ) stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
N = Nor - ns
model = spec@model
ipars = model$pars
pars = unlist(model$fixed.pars)
parnames = names(pars)
modelnames = .checkallfixed(spec)
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
cat("\nugarchforecast-->error: parameters names do not match specification\n")
cat("Expected Parameters are: ")
cat(paste(modelnames))
cat("\n")
stop("Exiting", call. = FALSE)
}
# once more into the spec
setfixed(spec)<-as.list(pars)
model = spec@model
idx = model$pidx
ipars = model$pars
modelinc = model$modelinc
model$modeldata$data = data
model$modeldata$index = index
model$modeldata$period = period
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(model, external.forecasts$mregfor, external.forecasts$vregfor, n.ahead, Nor, out.sample = ns, n.roll)
mxf = xreg$mxf
vxf = xreg$vxf
# filter data (check external regressor data - must equal length of origData)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "fiGARCH",
garchOrder = c(modelinc[8], modelinc[9]), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(modelinc[2], modelinc[3]),
include.mean = modelinc[1],
archm = ifelse(modelinc[5]>0,TRUE,FALSE), archpow = modelinc[5], arfima = modelinc[4],
external.regressors = mxf[1:(N + fcreq), , drop = FALSE], archex = modelinc[20]),
distribution.model = model$modeldesc$distribution, fixed.pars = as.list(pars))
tmp = xts(data[1:(N + fcreq)], index[1:(N + fcreq)])
flt = .figarchfilter(data = tmp, spec = fspec, n.old = N, trunclag = trunclag)
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
truncLag = trunclag
# ebar,eps and be from the filtered is now available
ebar =c(flt@filter$ebar,rep(0,n.ahead))
eps = flt@filter$eps
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(eps[-c(1:truncLag)], rep(0,n.ahead))
# forecast GARCH process
seriesFor = sigmaFor = matrix(NA, ncol = n.roll+1, nrow = n.ahead)
#seriesFor[1,] = fitted(flt)[(N+1):(N+n.roll+1)]
#sigmaFor[1,] = sigma(flt)[(N+1):(N+n.roll+1)]
# n.ahead x n.roll +1 (n.ahead=1 generted by model)
colnames(seriesFor) = colnames(sigmaFor) = as.character(index[N:(N+n.roll)])
rownames(seriesFor) = rownames(sigmaFor) = paste("T+", 1:n.ahead, sep="")
for(i in 1:(n.roll+1)){
np = N + i - 1
if(modelinc[1] > 0){
mu = rep(ipars[idx["mu",1]:idx["mu",2], 1], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
omega = rep(ipars[idx["omega",1]:idx["omega",2], 1], N+i+n.ahead-1)
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(data[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]
vxfi = vxf[1:(N+i-1+n.ahead), , drop = FALSE]
ans = .nfigarchforecast(ipars, modelinc, idx, mu, omega, mxfi, vxfi, h, epsx, z, eps=eps, ebar=ebar, be=be, truncLag=truncLag, data = x, N = np, n.ahead)
sigmaFor[,i] = ans$h
seriesFor[,i] = ans$x
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigmaFor = sigmaFor
fcst$seriesFor = seriesFor
model$modeldata$sigma = flt@filter$sigma
model$modeldata$residuals = flt@filter$residuals
ans = new("uGARCHforecast",
forecast = fcst,
model = model)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION sGARCH simulate
#---------------------------------------------------------------------------------
.figarchsim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL, ...)
{
# if( (n.sim+n.start) < 1000 && m.sim > 100 ){
# ans = .figarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
# startMethod = startMethod, presigma = presigma, prereturns = prereturns,
# preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
# mexsimdata = mexsimdata, vexsimdata = vexsimdata)
# } else{
ans = .figarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
startMethod = startMethod, presigma = presigma, prereturns = prereturns,
preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
mexsimdata = mexsimdata, vexsimdata = vexsimdata)
#}
return( ans )
}
.figarchsim1 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
if(fit@model$modelinc[4]>0){
if(n.start<fit@model$modelinc[3]){
warning("\nugarchsim-->warning: n.start>=MA order for arfima model...automatically setting.")
n.start = fit@model$modelinc[3]
}
}
# some checks
if(is.na(rseed[1])){
sseed = NA
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
# need to allow for arfima case:
n = n.sim + n.start
startMethod = startMethod[1]
data = fit@model$modeldata$data
N = length(as.numeric(data))
data = data[1:(N - fit@model$n.start)]
N = length(as.numeric(data))
m = fit@model$maxOrder
resids = fit@fit$residuals
sigma = fit@fit$sigma
model = fit@model
modelinc = model$modelinc
idx = model$pidx
ipars = fit@fit$ipars
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(model, mexsimdata, vexsimdata, N, n, m.sim, m)
mexsim = xreg$mexsimlist
vexsim = xreg$vexsimlist
if(N < n.start){
startmethod[1] = "unconditional"
warning("\nugarchsim-->warning: n.start greater than length of data...using unconditional start method...\n")
}
# Random Samples from the Distribution
if(length(sseed) == 1){
zmatrix = data.frame(dist = model$modeldesc$distribution, lambda = ipars[idx["ghlambda",1], 1],
skew = ipars[idx["skew",1], 1], shape = ipars[idx["shape",1], 1], n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(model$modeldesc$distribution, m.sim), lambda = rep(ipars[idx["ghlambda",1], 1], m.sim),
skew = rep(ipars[idx["skew",1], 1], m.sim), shape = rep(ipars[idx["shape",1], 1], m.sim),
n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
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)
}
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
if(startMethod[1] == "unconditional"){
#warning("\nunconditional start method not available for FIGARCH model. Using sample method.")
presigma = tail(sigma, m)
} else{
presigma = tail(sigma, m)
}
}
if(is.na(prereturns[1])){
if(startMethod[1] == "unconditional"){
prereturns = as.numeric(rep(uncmean(fit), m))
}
else{
prereturns = tail(data, m)
}
}
truncLag = fit@model$trunclag
# ebar,eps and be from the filtered is now available
ebar =c(tail(fit@fit$ebar, m),rep(0,n))
eps = fit@fit$eps
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(tail(eps,truncLag+m), rep(0,n))
# input vectors/matrices
h = c(presigma^2, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(ipars[idx["mu",1]:idx["mu",2], 1], ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
ebarSim =matrix(0, ncol = m.sim, nrow = n.sim)
epsSim = matrix(0, ncol = m.sim, nrow = n.sim)
z[is.na(z) | is.nan(z) | !is.finite(z)] = 0
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
if(startMethod[1] == "unconditional"){
preres = as.numeric(z[1:m, i])*presigma
} else{
preres = tail(resids, m)
}
}
res = c(preres, rep(0, n))
ans1 = try(.C("figarchsimC", model = as.integer(modelinc[1:21]), pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
h = as.double(h), z = as.double(z[,i]), res = as.double(res), e = as.double(res*res),
ebar = as.double(ebar), eps = as.double(eps), be = as.double(be),
vexdata = as.double(vexsim[[i]]), T = as.integer(n+m), N = as.integer(truncLag), m = as.integer(m), PACKAGE = "rugarch"), silent = TRUE)
if(inherits(ans1, "try-error")) stop("\nugarchsim-->error: error in calling C function....\n")
sigmaSim[,i] = ans1$h[(n.start + m + 1):(n+m)]^(1/2)
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
ebarSim[,i] = ans1$ebar[(n.start + m + 1):(n+m)]
epsSim[,i] = ans1$eps[(n.start + m + 1):(n+m)]
# ToDo: change to accomodate modelinc[20]
if(modelinc[6]>0){
mxreg = matrix( ipars[idx["mxreg",1]:idx["mxreg",2], 1], ncol = modelinc[6] )
if(modelinc[20]==0){
constm[,i] = constm[,i] + mxreg %*%t( matrix( mexsim[[i]], ncol = modelinc[6] ) )
} else{
if(modelinc[20] == modelinc[6]){
constm[,i] = constm[,i] + mxreg %*%t( matrix( mexsim[[i]]*sqrt(ans1$h), ncol = modelinc[6] ) )
} else{
constm[,i] = constm[,i] + mxreg[,1:(modelinc[6]-modelinc[20]),drop=FALSE] %*%t( matrix( mexsim[[i]][,1:(modelinc[6]-modelinc[20]),drop=FALSE], ncol = modelinc[6]-modelinc[20] ) )
constm[,i] = constm[,i] + mxreg[,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE] %*%t( matrix( mexsim[[i]][,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE]*sqrt(ans1$h), ncol = modelinc[20] ) )
}
}
}
if(modelinc[5]>0) constm[,i] = constm[,i] + ipars[idx["archm",1]:idx["archm",2], 1]*(sqrt(ans1$h)^modelinc[5])
if(modelinc[4]>0){
fres = c(ans1$res[(m+1):(n+m)], if(modelinc[3]>0) rep(0, modelinc[3]) else NULL)
ans2 = .arfimaxsim(modelinc[1:21], ipars, idx, constm[1:n, i], fres, T = n)
seriesSim[,i] = head(ans2$series, n.sim)
} else{
ans2 = .armaxsim(modelinc[1:21], ipars, idx, constm[,i], x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim, ebarSim = ebarSim, epsSim = epsSim, be = be)
model$modeldata$sigma = sigma
sol = new("uGARCHsim",
simulation = sim,
model = model,
seed = as.integer(sseed))
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION sGARCH path
#---------------------------------------------------------------------------------
.figarchpath = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL, trunclag=1000, ...)
{
ans = .figarchpath1(spec = spec, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
presigma = presigma, prereturns = prereturns, preresiduals = preresiduals,
rseed = rseed, custom.dist = custom.dist, mexsimdata = mexsimdata,
vexsimdata = vexsimdata, trunclag=trunclag)
return( ans )
}
.figarchpath1 = function(spec, n.sim = 1000, n.start = 0, m.sim = 1, presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL, trunclag=1000)
{
if(spec@model$modelinc[4]>0){
if(n.start<spec@model$modelinc[3]){
warning("\nugarchpath-->warning: n.start>=MA order for arfima model...automatically setting.")
n.start = spec@model$modelinc[3]
}
}
# some checks
if(is.na(rseed[1])){
sseed = NA
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
model = spec@model
ipars = model$pars
pars = unlist(model$fixed.pars)
parnames = names(pars)
modelnames = .checkallfixed(spec)
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
cat("\nugarchpath-->error: parameters names do not match specification\n")
cat("Expected Parameters are: ")
cat(paste(modelnames))
cat("\n")
stop("Exiting", call. = FALSE)
}
# once more into the spec
setfixed(spec)<-as.list(pars)
model = spec@model
ipars = model$pars
idx = model$pidx
modelinc = model$modelinc
# Enlarge Series:
n = n.sim + n.start
m = model$maxOrder
N = 0
if(modelinc[6]>0) {
mexdata = matrix(model$modeldata$mexdata, ncol = modelinc[6])
N = dim(mexdata)[1]
} else { mexdata = NULL }
if(modelinc[15]>0) {
vexdata = matrix(model$modeldata$vexdata, ncol = modelinc[15])
N = dim(vexdata)[1]
} else { vexdata = NULL }
distribution = model$modeldesc$distribution
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(model, mexsimdata, vexsimdata, N, n, m.sim, m)
mexsim = xreg$mexsimlist
vexsim = xreg$vexsimlist
# check positivity of coefficients
# Random Samples from the Distribution
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, lambda = ipars[idx["ghlambda",1], 1],
skew = ipars[idx["skew",1], 1], shape = ipars[idx["shape",1], 1],
n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), lambda = rep(ipars[idx["ghlambda",1], 1], m.sim),
skew = rep(ipars[idx["skew",1], 1], m.sim), shape = rep(ipars[idx["shape",1], 1], 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)
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchpath-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchpath-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchpath-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
hEst = ipars["omega",1]/(1-ipars["beta1",1]-ipars["alpha1",1])
# we need some estimate to start it off (use the wrong one but not unreasonable)
presigma = as.numeric(rep(sqrt(hEst), times = m))
}
if(is.na(prereturns[1])){
prereturns = as.numeric(rep(uncmean(spec), times = m))
}
truncLag = trunclag
ebar =c(rep(0, m),rep(0,n))
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(rep(presigma^2,truncLag+m), rep(0,n))
# input vectors/matrices
h = c(presigma^2, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(ipars[idx["mu",1]:idx["mu",2],1], ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
ebarSim =matrix(0, ncol = m.sim, nrow = n.sim)
epsSim = matrix(0, ncol = m.sim, nrow = n.sim)
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
preres = as.numeric(z[1:m,i])*presigma
}
z[1:m, 1:m.sim] = preres[1:m]/presigma[1:m]
z[is.na(z) | is.nan(z) | !is.finite(z)] = 0
res = c(preres, rep(0, n))
ans1 = try(.C("figarchsimC", model = as.integer(modelinc[1:21]), pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
h = as.double(h), z = as.double(z[,i]), res = as.double(res), e = as.double(res*res),
ebar = as.double(ebar), eps = as.double(eps), be = as.double(be),
vexdata = as.double(vexsim[[i]]), T = as.integer(n+m), N = as.integer(truncLag),
m = as.integer(m), PACKAGE = "rugarch"), silent = TRUE)
if(inherits(ans1, "try-error")) stop("\nugarchpath-->error: error in calling C function....\n")
sigmaSim[,i] = ans1$h[(n.start + m + 1):(n+m)]^(1/2)
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
ebarSim[,i] = ans1$ebar[(n.start + m + 1):(n+m)]
epsSim[,i] = ans1$eps[(n.start + m + 1):(n+m)]
if(modelinc[6]>0){
mxreg = matrix( ipars[idx["mxreg",1]:idx["mxreg",2], 1], ncol = modelinc[6] )
if(modelinc[20]==0){
constm[,i] = constm[,i] + as.numeric( mxreg %*%t( matrix( mexsim[[i]], ncol = modelinc[6] ) ) )
} else{
if(modelinc[20] == modelinc[6]){
constm[,i] = constm[,i] + as.numeric( mxreg %*%t( matrix( mexsim[[i]]*sqrt(ans1$h), ncol = modelinc[6] ) ) )
} else{
constm[,i] = constm[,i] + as.numeric( mxreg[,1:(modelinc[6]-modelinc[20]),drop=FALSE] %*%t( matrix( mexsim[[i]][,1:(modelinc[6]-modelinc[20]),drop=FALSE], ncol = modelinc[6]-modelinc[20] ) ) )
constm[,i] = constm[,i] + as.numeric( mxreg[,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE] %*%t( matrix( mexsim[[i]][,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE]*sqrt(ans1$h), ncol = modelinc[20] ) ) )
}
}
}
if(modelinc[5]>0) constm[,i] = constm[,i] + ipars[idx["archm",1]:idx["archm",2], 1]*(sqrt(ans1$h)^modelinc[5])
if(modelinc[4]>0){
fres = c(ans1$res[(m+1):(n+m)], if(modelinc[3]>0) rep(0, modelinc[3]) else NULL)
ans2 = .arfimaxsim(modelinc[1:21], ipars, idx, constm[1:n, i], fres, T = n)
seriesSim[,i] = head(ans2$series, n.sim)
} else{
#if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .armaxsim(modelinc[1:21], ipars, idx, constm[,i], x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
path = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sol = new("uGARCHpath",
path = path,
model = model,
seed = as.integer(sseed))
return(sol)
}
Try the rugarch package in your browser
Any scripts or data that you put into this service are public.
rugarch documentation built on Sept. 20, 2023, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .figarchpath1 .figarchpath .figarchsim1 .figarchsim .figarchforecast2 .nfigarchforecast .figarchforecast .figarchfilter .figarchLLH .figarchfit
#################################################################################
##
## R package rugarch by Alexios Galanos Copyright (C) 2008-2022.
## This file is part of the R package rugarch.
##
## The R package rugarch 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 rugarch 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 fiGARCH fit
#---------------------------------------------------------------------------------
# [mu ar ma arfima im mxreg omega alpha beta gamma gamma11 gamma21 delta lambda vxreg skew shape dlamda aux aux aux aux]
.figarchfit = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(stationarity = 1, fixed.se = 0, scale = 0, rec.init = 'all',trunclag = 1000),
numderiv.control = list(grad.eps=1e-4, grad.d=0.0001, grad.zero.tol=sqrt(.Machine$double.eps/7e-7),
hess.eps=1e-4, hess.d=0.1, hess.zero.tol=sqrt(.Machine$double.eps/7e-7), r=4, v=2))
{
tic = Sys.time()
if(is.null(solver.control$trace)) trace = 0 else trace = solver.control$trace
if(is.null(fit.control$stationarity)) fit.control$stationarity = TRUE
if(is.null(fit.control$fixed.se)) fit.control$fixed.se = FALSE
if(is.null(fit.control$scale)) fit.control$scale = FALSE
if(is.null(fit.control$rec.init)) fit.control$rec.init = 'all'
if(is.null(fit.control$trunclag)) truncLag = 1000 else truncLag = fit.control$trunclag
mm = match(names(fit.control), c("stationarity", "fixed.se", "scale", "rec.init", "trunclag"))
if(any(is.na(mm))){
idx = which(is.na(mm))
enx = NULL
for(i in 1:length(idx)) enx = c(enx, names(fit.control)[idx[i]])
warning(paste(c("unidentified option(s) in fit.control:\n", enx), sep="", collapse=" "), call. = FALSE, domain = NULL)
}
# if we have external regressors in variance turn off scaling
if(spec@model$modelinc[15] > 0) fit.control$scale = FALSE
# if we have arch-in-mean turn off scaling
if(spec@model$modelinc[5] > 0) fit.control$scale = FALSE
# if there are fixed pars we do no allow scaling as there would be no way of mixing scaled
# amd non scaled parameters
if(sum(spec@model$pars[,2]) > 0) fit.control$scale = FALSE
xdata = .extractdata(data)
if(!is.numeric(out.sample)) stop("\nugarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample)<0) stop("\nugarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample,0)
n = length(xdata$data)
#if((n-n.start)<100) stop("\nugarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start)]
index = xdata$index[1:(n-n.start)]
origdata = xdata$data
origindex = xdata$index
period = xdata$period
# arglist replaces the use of custom environments (resolves the problem of
# non-shared environment in parallel estimation, particularly windows)
garchenv = new.env(hash = TRUE)
arglist = list()
arglist$garchenv <- garchenv
arglist$pmode = 0
arglist$truncLag = truncLag
model = spec@model
modelinc = model$modelinc
pidx = model$pidx
# expand the spec object and assign spec lists
if(modelinc[6] > 0){
mexdata = model$modeldata$mexdata[1:(n-n.start), , drop = FALSE]
} else{
mexdata = NULL
}
if(modelinc[15] > 0){
vexdata = model$modeldata$vexdata[1:(n-n.start), ,drop = FALSE]
} else{
vexdata = NULL
}
arglist$index = index
arglist$trace = trace
m = model$maxOrder
# store length of data for easy retrieval
model$modeldata$T = T = length(as.numeric(data))
dist = model$modeldesc$distribution
if(fit.control$scale) dscale = sd(data) else dscale = 1
zdata = data/dscale
recinit = .checkrec(fit.control$rec.init, T)
arglist$data = zdata
arglist$recinit = recinit
arglist$dscale = dscale
arglist$model = model
ipars = model$pars
# Optimization Starting Parameters Vector & Bounds
tmp = .garchstart(ipars, arglist)
####################################################
# Scale Omega in order to avoid problems: (25-01-2013)
#arglist$omegascale = 1/1e8
#tmp$ipars["omega",c(1,5,6)] = tmp$ipars["omega",c(1,5,6)] *1e8
####################################################
ipars = arglist$ipars = tmp$pars
arglist$tmph = tmp$tmph
# we now split out any fixed parameters
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
arglist$fit.control = fit.control
npars = sum(estidx)
if(any(ipars[,2]==1)){
if(npars == 0){
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
warning("\nugarchfit-->warning: all parameters fixed...returning ugarchfilter object instead\n")
return(ugarchfilter(data = xts(origdata, origindex), spec = spec, out.sample = out.sample))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = 0
ipars[ipars[,2]==1, 4] = 1
ipars[ipars[,2]==1, 2] = 0
arglist$ipars = ipars
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = 1
}
} else{
use.solver = 1
}
# start counter
assign("rugarch_llh", 1, envir = garchenv)
# assign solver constraints (solnp directly else exterior type penalty
# for other solvers)
if(fit.control$stationarity == 1 && modelinc[15] == 0){
Ifn = .figarchcon
ILB = rep(0,2)
IUB = rep(1,2)
if(solver == "solnp" | solver == "gosolnp" | solver == "hybrid") fit.control$stationarity = 0
} else{
Ifn = NULL
ILB = NULL
IUB = NULL
}
# conditions controls the non-solnp solver penalty
arglist$fit.control = fit.control
if(use.solver){
parscale = rep(1, length = npars)
names(parscale) = rownames(ipars[estidx,])
if(modelinc[1] > 0) parscale["mu"] = abs(mean(zdata))
if(modelinc[7] > 0) parscale["omega"] = var(zdata)
arglist$returnType = "llh"
solution = .garchsolver(solver, pars = ipars[estidx, 1], fun = .figarchLLH,
Ifn, ILB, IUB, gr = NULL, hessian = NULL, parscale = parscale,
control = solver.control, LB = ipars[estidx, 5],
UB = ipars[estidx, 6], ux = NULL, ci = NULL, mu = NULL, arglist)
sol = solution$sol
hess = solution$hess
timer = Sys.time()-tic
if(!is.null(sol$par)){
ipars[estidx, 1] = sol$par
if(modelinc[7]==0){
# call it once more to get omega
tmpx = .figarchLLH(sol$par, arglist)
ipars[pidx["omega",1], 1] = get("omega", garchenv)
}
if(sum(ipars[,2]) == 0){
if(modelinc[1] > 0) ipars[pidx["mu",1]:pidx["mu",2], 1] = ipars[pidx["mu",1]:pidx["mu",2], 1] * dscale
if(modelinc[6] > 0){
ipars[pidx["mxreg", 1]:pidx["mxreg", 2], 1] = ipars[pidx["mxreg", 1]:pidx["mxreg", 2], 1] * dscale
}
ipars[pidx["omega",1], 1] = ipars[pidx["omega",1],1] * dscale^2
}
} else{
ipars[estidx, 1] = NA
}
arglist$ipars = ipars
convergence = sol$convergence
if(convergence != 0) warning("\nugarchfit-->warning: solver failer to converge.")
} else{
solution = NULL
hess = NULL
timer = Sys.time()-tic
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
# check convergence else write message/return
# create a copy of ipars in case we need to change it below to calculate standard errors
# which we will need to reset later (because for example, infocriteria uses estimated
# parameters, not fixed.
ipars2 = ipars
if(convergence == 0){
arglist$dscale = 1
if(sum(ipars[,2]) > 0 && fit.control$fixed.se == 1){
ipars[ipars[,2]==1, 4] = 1
ipars[ipars[,2]==1, 2] = 0
arglist$ipars = ipars
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
}
arglist$data = data
fit = .makefitmodel(garchmodel = "fiGARCH", f = .figarchLLH, T = T, m = m,
timer = timer, convergence = convergence, message = sol$message,
hess = hess, arglist = arglist, numderiv.control = numderiv.control)
model$modelinc[7] = modelinc[7]
model$modeldata$data = origdata
model$modeldata$index = origindex
model$modeldata$period = period
model$pars[, 1] = fit$ipars[,1]
model$pars[, 5:6] = ipars2[,5:6]
fit$ipars[, 4] = ipars2[, 4]
fit$ipars[, 2] = ipars2[, 2]
fit$ipars[, 5:6] = ipars2[,5:6]
# make sure omega is now included (for working with object post-estimation)
fit$ipars["omega", 3] = 1
model$pars["omega", 3] = 1
} else{
fit$message = sol$message
fit$convergence = 1
model$modeldata$data = origdata
model$modeldata$index = origindex
model$modeldata$period = period
}
# make model list to return some usefule information which
# will be called by other functions (show, plot, sim etc)
model$n.start = n.start
fit$solver = solution
model$trunclag = truncLag
ans = new("uGARCHfit",
fit = fit,
model = model)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION fiGARCH LLH
#---------------------------------------------------------------------------------
.figarchLLH = function(pars, arglist)
{
# prepare inputs
data = arglist$data
returnType = arglist$returnType
garchenv = arglist$garchenv
# rejoin fixed and pars
model = arglist$model
estidx = arglist$estidx
idx = model$pidx
ipars = arglist$ipars
ipars[estidx, 1] = pars
trace = arglist$trace
T = length(data)
dscale = arglist$dscale
recinit = arglist$recinit
fit.control = arglist$fit.control
m = model$maxOrder
N = c(m,T)
truncLag = arglist$truncLag
mexdata = coredata(model$modeldata$mexdata[1:T,, drop = FALSE])
vexdata = coredata(model$modeldata$vexdata[1:T,, drop = FALSE])
distribution = model$modeldesc$distribution
modelinc = model$modelinc
# persist used equivalently here for the positivity constraints
persist = .figarchcon(pars, arglist)
# distribution number
# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
dist = model$modeldesc$distno
hm = arglist$tmph
rx = .arfimaxfilter(modelinc[1:21], pars = ipars[,1], idx = idx, mexdata = mexdata, h = hm, data = data, N = N)
res = rx$res
zrf = rx$zrf
res[is.na(res) | !is.finite(res) | is.nan(res)] = 0
# recursion initialization
mvar = ifelse(recinit$type==1, mean(res[1:recinit$n]*res[1:recinit$n]), backcastv(res, T, recinit$n))
if(modelinc[15]>0) {
mv = sum(apply(matrix(vexdata, ncol = modelinc[15]), 2, "mean")*ipars[idx["vxreg",1]:idx["vxreg",2],1])
} else{
mv = 0
}
hEst = mvar
# variance targeting not used
ipars[idx["omega",1],1] = max(eps, ipars[idx["omega",1],1])
if(is.na(hEst) | !is.finite(hEst) | is.nan(hEst)) hEst = var(data)
delta=ipars[idx["delta",1],1]
# figarch setup
NL=T+truncLag
eps=matrix(mvar, ncol=1, nrow=NL)
ebar = rep(0, T)
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
# if we have external regressors in variance equation we cannot have
# stationarity checks in likelihood. Stationarity conditions not valid for
# solnp solver which implements them internally as constraints.
if(fit.control$stationarity == 1 && modelinc[15] == 0){
if(any(is.na(persist)) || any(persist >= 1) || any(persist <= 0)){
if(arglist$pmode!=1){
return(llh = get("rugarch_llh", garchenv) + 0.1*(abs(get("rugarch_llh", garchenv))))
} else{
return(llh = 1e10)
}
}
}
ebar[1]=t(be)%*%eps[1:truncLag]
if(modelinc[6]>0) mexdata = as.double(as.vector(mexdata)) else mexdata = double(1)
if(modelinc[15]>0) vexdata = as.double(as.vector(vexdata)) else vexdata = double(1)
# modelinc (1:21) since 22 is either NA or numeric and no longer needed (paased to ipars if used)
ans = try( .C("figarchfilterC", model = as.integer(modelinc[1:21]),
pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
hEst = as.double(hEst), x = as.double(data), res = as.double(res),
e = double(T), ebar = as.double(ebar), eps=as.double(eps), be=as.double(be),
mexdata = mexdata, vexdata = vexdata,
zrf = as.double(zrf), constm = double(T), condm = double(T),
m = as.integer(m), T = as.integer(T), N = as.integer(truncLag),
h = double(T), z = double(T), llh = double(1), LHT = double(T),
PACKAGE = "rugarch"), silent = TRUE )
if(inherits(ans, "try-error")){
if(arglist$pmode!=1){
assign("rugarch_csol", 1, envir = garchenv)
assign("rugarch_filtermessage", ans, envir = garchenv)
if( trace > 0 ) cat(paste("\narfimafit-->warning: ", get("rugarch_filtermessage", garchenv),"\n", sep=""))
return(llh = (get("rugarch_llh", garchenv) + 0.1*(abs(get("rugarch_llh", garchenv)))))
} else{
return(llh = 1e10)
}
} else{
if(arglist$pmode!=1){
assign("rugarch_csol", 0, envir = garchenv)
}
}
z = ans$z
h = ans$h
epsx = ans$res
eps = ans$eps
ebar = ans$ebar
llh = ans$llh
if(is.finite(llh) && !is.na(llh) && !is.nan(llh)){
if(arglist$pmode!=1) assign("rugarch_llh", llh, envir = garchenv)
} else {
if(arglist$pmode!=1) llh = (get("rugarch_llh", garchenv) + 0.1*(abs(get("rugarch_llh", garchenv)))) else llh = 1e10
}
# LHT = raw scores
# ? -ans$LHT[(m+1):T]
LHT = -ans$LHT
ans = switch(returnType,
llh = llh,
LHT = LHT,
all = list(llh = llh, h = h, epsx = epsx, z = z, eps=eps, ebar=ebar, be=be, kappa = kappa, LHT = LHT, persistence = persist))
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION sGARCH filter
#---------------------------------------------------------------------------------
.figarchfilter = function(spec, data, out.sample = 0, n.old = NULL, rec.init = 'all', trunclag = 1000)
{
# 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.
xdata = .extractdata(data)
data = xdata$data
index = xdata$index
period = xdata$period
origdata = data
origindex = index
T = length(origdata) - out.sample
data = origdata[1:T]
index = origindex[1:T]
if(!is.null(n.old)) Nx = n.old else Nx = length(data)
recinit = .checkrec(rec.init, Nx)
model = spec@model
ipars = model$pars
pars = unlist(model$fixed.pars)
parnames = names(pars)
modelnames = .checkallfixed(spec)
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
cat("\nugarchfilter-->error: parameters names do not match specification\n")
cat("Expected Parameters are: ")
cat(paste(modelnames))
cat("\n")
stop("Exiting", call. = FALSE)
}
# once more into the spec
# NB Any changes made to the spec are not preserved once we apply set fixed
setfixed(spec)<-as.list(pars)
model = spec@model
model$modeldata$T = T
ipars = model$pars
idx = model$pidx
modelinc = model$modelinc
m = model$maxOrder
N = c(m,T)
mexdata = model$modeldata$mexdata[1:T, , drop = FALSE]
vexdata = model$modeldata$vexdata[1:T, , drop = FALSE]
distribution = model$modeldesc$distribution
dist = model$modeldesc$distno
kappa = 1
#persist = (sum(ipars[idx["alpha",1]:idx["alpha",2],1]) + sum(ipars[idx["beta",1]:idx["beta",2],1]))
rx = .arfimaxfilter(modelinc[1:21], ipars[,1], idx, mexdata = mexdata, h = 0, data = data, N = N)
res = rx$res
zrf = rx$zrf
if(!is.null(n.old)){
rx2 = .arfimaxfilter(modelinc[1:21], ipars[,1], idx, mexdata = mexdata[1:Nx, , drop = FALSE], h = 0, data = origdata[1:Nx], N = c(m, Nx))
res2 = rx2$res
mvar = ifelse(recinit$type==1, mean(res2[1:recinit$n]*res2[1:recinit$n]), backcastv(res2, Nx, recinit$n))
} else{
mvar = ifelse(recinit$type==1, mean(res[1:recinit$n]*res[1:recinit$n]), backcastv(res, T, recinit$n))
}
hEst = mvar
if(modelinc[15]>0) {
mv = sum(apply(matrix(vexdata, ncol = modelinc[15]), 2, "mean")*ipars[idx["vxreg",1]:idx["vxreg",2],1])
} else{
mv = 0
}
ipars[idx["omega",1],1] = max(eps, ipars[idx["omega",1],1])
if(is.na(hEst) | !is.finite(hEst) | is.nan(hEst)) hEst = var(data)
delta=ipars[idx["delta",1],1]
# figarch setup
truncLag=trunclag
NL=T+truncLag
eps=matrix(mvar, ncol=1, nrow=NL)
ebar = rep(0, T)
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
ebar[1]=t(be)%*%eps[1:truncLag]
arglist=list()
arglist$ipars=ipars
arglist$truncLag = truncLag
persist=.figarchcon(pars, arglist)
if(any(persist<=0) | any(persist>=1)) warning("\nPositivity Contraints NOT met (check parameters)!")
if(modelinc[6]>0) mexdata = as.double(as.vector(mexdata)) else mexdata = double(1)
if(modelinc[15]>0) vexdata = as.double(as.vector(vexdata)) else vexdata = double(1)
ans = try( .C("figarchfilterC", model = as.integer(modelinc[1:21]),
pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
hEst = as.double(hEst), x = as.double(data), res = as.double(res),
e = double(T), ebar = as.double(ebar), eps=as.double(eps), be=as.double(be),
mexdata = mexdata, vexdata = vexdata,
zrf = as.double(zrf), constm = double(T), condm = double(T),
m = as.integer(m), T = as.integer(T), N = as.integer(truncLag),
h = double(T), z = double(T), llh = double(1), LHT = double(T),
PACKAGE = "rugarch"), silent = TRUE )
filter = list()
filter$z = ans$z
filter$ebar = ans$ebar
filter$eps = ans$eps
filter$be = ans$be
filter$sigma = sqrt(ans$h)
filter$residuals = ans$res
filter$LLH = -ans$llh
filter$log.likelihoods = ans$LHT
filter$persistence = persist
filter$distribution = distribution
filter$ipars = ipars
model$modeldata$data = origdata
model$modeldata$index = origindex
model$modeldata$period = period
model$n.start = out.sample
model$trunclag = trunclag
sol = new("uGARCHfilter",
filter = filter,
model = model)
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION fiGARCH forecast
#---------------------------------------------------------------------------------
# Recipe for a GARCH forecast:
# 1. split forecast into arma and garch forecasts.
# 2. take into account garch-in-mean, arfima, exogenous regressors
# 3. check kappa and persistence
# 4. allow for rolling forecast based on filtering out.sample data from fit stage
# n.roll signifies whether to filter/roll the sigma else will use unconditional
# i.e. n.ahead=10 with n.roll=1 means that the first forecast is based on the
# previous value whereas all subsequent forecasts are based on the unconditional
# expectation formula
# with n.roll = n means that while n.ahead < n.roll we filter the data using the coef
# of the garch model to obtain filtered sigma estimates which are used to roll
# the forecast. This requires that during the fit process the out.sample option
# was used and that (n.roll) < out.sample (otherwise we revert to the unconditional
# expectation formula for the long run sigma.
.figarchforecast = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
fit = fitORspec
data = fit@model$modeldata$data
Nor = length(as.numeric(data))
index = fit@model$modeldata$index
period = fit@model$modeldata$period
ns = fit@model$n.start
N = Nor - ns
model = fit@model
ipars = fit@fit$ipars
modelinc = model$modelinc
idx = model$pidx
if( n.roll > ns ) stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
pars = fit@fit$coef
ipars = fit@fit$ipars
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(model, external.forecasts$mregfor, external.forecasts$vregfor, n.ahead, Nor, out.sample = ns, n.roll)
mxf = xreg$mxf
vxf = xreg$vxf
# filter data (check external regressor data - must equal length of origData)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "fiGARCH",
garchOrder = c(modelinc[8], modelinc[9]), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(modelinc[2], modelinc[3]),
include.mean = modelinc[1],
archm = ifelse(modelinc[5]>0,TRUE,FALSE), archpow = modelinc[5], arfima = modelinc[4],
external.regressors = mxf[1:(N + fcreq), , drop = FALSE], archex = modelinc[20]),
distribution.model = model$modeldesc$distribution, fixed.pars = as.list(pars))
tmp = xts(data[1:(N + fcreq)], index[1:(N + fcreq)])
flt = .figarchfilter(data = tmp, spec = fspec, n.old = N, trunclag = fit@model$trunclag)
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
truncLag = fit@model$trunclag
# ebar,eps and be from the filtered is now available
ebar =c(flt@filter$ebar,rep(0,n.ahead))
eps = flt@filter$eps
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(eps[-c(1:truncLag)], rep(0,n.ahead))
# forecast GARCH process
seriesFor = sigmaFor = matrix(NA, ncol = n.roll+1, nrow = n.ahead)
#seriesFor[1,] = fitted(flt)[(N+1):(N+n.roll+1)]
#sigmaFor[1,] = sigma(flt)[(N+1):(N+n.roll+1)]
# n.ahead x n.roll +1 (n.ahead=1 generted by model)
colnames(seriesFor) = colnames(sigmaFor) = as.character(index[N:(N+n.roll)])
rownames(seriesFor) = rownames(sigmaFor) = paste("T+", 1:n.ahead, sep="")
for(i in 1:(n.roll+1)){
np = N + i - 1
if(modelinc[1] > 0){
mu = rep(ipars[idx["mu",1]:idx["mu",2], 1], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
omega = rep(ipars[idx["omega",1]:idx["omega",2], 1], N+i+n.ahead-1)
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(data[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]
vxfi = vxf[1:(N+i-1+n.ahead), , drop = FALSE]
# CONTINUE Here
# need to append ebar/eps/ and indicate nlagbin
ans = .nfigarchforecast(ipars, modelinc, idx, mu, omega, mxfi, vxfi, h, epsx, z, eps=eps, ebar=ebar, be=be, truncLag=truncLag, data = x, N = np, n.ahead)
sigmaFor[,i] = ans$h
seriesFor[,i] = ans$x
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigmaFor = sigmaFor
fcst$seriesFor = seriesFor
model$modeldata$sigma = flt@filter$sigma
model$modeldata$residuals = flt@filter$residuals
ans = new("uGARCHforecast",
forecast = fcst,
model = model)
return(ans)
}
.nfigarchforecast = function(ipars, modelinc, idx, mu, omega, mxfi, vxfi, h, epsx, z, eps, ebar, be, truncLag, data, N, n.ahead)
{
if(modelinc[15]>0){
omega = omega + vxfi%*%t(matrix(ipars[idx["vxreg",1]:idx["vxreg",2],1], ncol = modelinc[15]))
}
for(i in 1:n.ahead){
ebar[N+i] = be%*%eps[(N-truncLag+i):(N+i-1)]
h[N+i] = omega[N+i] - ebar[N+i]
if(modelinc[9]>=i){
h[N+i] = h[N+i] + sum(ipars[idx["beta",1]:idx["beta",2],1]*(h[N+i-(1:modelinc[9])]^2-eps[(N+i-(1:modelinc[9]))]))
}
if(modelinc[8]>0){
for (j in 1:modelinc[8]){
if ((i-j) > 0){
s = h[N + i - j]^2
} else{
s = epsx[N + i - j]^2
}
h[N+i] = h[N+i] + ipars[idx["alpha",1]+j-1,1] * (s + ebar[N + i - j])
}
}
eps[N+i] = h[N+i]
h[N+i] = sqrt(h[N+i])
}
if(modelinc[4]>0){
res = arfimaf(ipars, modelinc[1:21], idx, mu, mxfi, h, epsx, z, data, N, n.ahead)
} else{
res = armaf(ipars, modelinc[1:21], idx, mu, mxfi, h, epsx, z, data, N, n.ahead)
}
return(list(h = h[(N+1):(N+n.ahead)], x = res[(N+1):(N+n.ahead)], ebar=ebar[(N+1):(N+n.ahead)], eps=eps[(N+1):(N+n.ahead)]))
}
#---------------------------------------------------------------------------------
# 2nd dispatch method for forecast
.figarchforecast2 = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), trunclag = 1000, ...)
{
# 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
N = length(as.numeric(data))
index = xdata$index
period = xdata$period
ns = out.sample
if( n.roll > ns ) stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
N = Nor - ns
model = spec@model
ipars = model$pars
pars = unlist(model$fixed.pars)
parnames = names(pars)
modelnames = .checkallfixed(spec)
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
cat("\nugarchforecast-->error: parameters names do not match specification\n")
cat("Expected Parameters are: ")
cat(paste(modelnames))
cat("\n")
stop("Exiting", call. = FALSE)
}
# once more into the spec
setfixed(spec)<-as.list(pars)
model = spec@model
idx = model$pidx
ipars = model$pars
modelinc = model$modelinc
model$modeldata$data = data
model$modeldata$index = index
model$modeldata$period = period
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(model, external.forecasts$mregfor, external.forecasts$vregfor, n.ahead, Nor, out.sample = ns, n.roll)
mxf = xreg$mxf
vxf = xreg$vxf
# filter data (check external regressor data - must equal length of origData)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "fiGARCH",
garchOrder = c(modelinc[8], modelinc[9]), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(modelinc[2], modelinc[3]),
include.mean = modelinc[1],
archm = ifelse(modelinc[5]>0,TRUE,FALSE), archpow = modelinc[5], arfima = modelinc[4],
external.regressors = mxf[1:(N + fcreq), , drop = FALSE], archex = modelinc[20]),
distribution.model = model$modeldesc$distribution, fixed.pars = as.list(pars))
tmp = xts(data[1:(N + fcreq)], index[1:(N + fcreq)])
flt = .figarchfilter(data = tmp, spec = fspec, n.old = N, trunclag = trunclag)
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
truncLag = trunclag
# ebar,eps and be from the filtered is now available
ebar =c(flt@filter$ebar,rep(0,n.ahead))
eps = flt@filter$eps
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(eps[-c(1:truncLag)], rep(0,n.ahead))
# forecast GARCH process
seriesFor = sigmaFor = matrix(NA, ncol = n.roll+1, nrow = n.ahead)
#seriesFor[1,] = fitted(flt)[(N+1):(N+n.roll+1)]
#sigmaFor[1,] = sigma(flt)[(N+1):(N+n.roll+1)]
# n.ahead x n.roll +1 (n.ahead=1 generted by model)
colnames(seriesFor) = colnames(sigmaFor) = as.character(index[N:(N+n.roll)])
rownames(seriesFor) = rownames(sigmaFor) = paste("T+", 1:n.ahead, sep="")
for(i in 1:(n.roll+1)){
np = N + i - 1
if(modelinc[1] > 0){
mu = rep(ipars[idx["mu",1]:idx["mu",2], 1], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
omega = rep(ipars[idx["omega",1]:idx["omega",2], 1], N+i+n.ahead-1)
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(data[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]
vxfi = vxf[1:(N+i-1+n.ahead), , drop = FALSE]
ans = .nfigarchforecast(ipars, modelinc, idx, mu, omega, mxfi, vxfi, h, epsx, z, eps=eps, ebar=ebar, be=be, truncLag=truncLag, data = x, N = np, n.ahead)
sigmaFor[,i] = ans$h
seriesFor[,i] = ans$x
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigmaFor = sigmaFor
fcst$seriesFor = seriesFor
model$modeldata$sigma = flt@filter$sigma
model$modeldata$residuals = flt@filter$residuals
ans = new("uGARCHforecast",
forecast = fcst,
model = model)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION sGARCH simulate
#---------------------------------------------------------------------------------
.figarchsim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL, ...)
{
# if( (n.sim+n.start) < 1000 && m.sim > 100 ){
# ans = .figarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
# startMethod = startMethod, presigma = presigma, prereturns = prereturns,
# preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
# mexsimdata = mexsimdata, vexsimdata = vexsimdata)
# } else{
ans = .figarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
startMethod = startMethod, presigma = presigma, prereturns = prereturns,
preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
mexsimdata = mexsimdata, vexsimdata = vexsimdata)
#}
return( ans )
}
.figarchsim1 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
if(fit@model$modelinc[4]>0){
if(n.start<fit@model$modelinc[3]){
warning("\nugarchsim-->warning: n.start>=MA order for arfima model...automatically setting.")
n.start = fit@model$modelinc[3]
}
}
# some checks
if(is.na(rseed[1])){
sseed = NA
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
# need to allow for arfima case:
n = n.sim + n.start
startMethod = startMethod[1]
data = fit@model$modeldata$data
N = length(as.numeric(data))
data = data[1:(N - fit@model$n.start)]
N = length(as.numeric(data))
m = fit@model$maxOrder
resids = fit@fit$residuals
sigma = fit@fit$sigma
model = fit@model
modelinc = model$modelinc
idx = model$pidx
ipars = fit@fit$ipars
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(model, mexsimdata, vexsimdata, N, n, m.sim, m)
mexsim = xreg$mexsimlist
vexsim = xreg$vexsimlist
if(N < n.start){
startmethod[1] = "unconditional"
warning("\nugarchsim-->warning: n.start greater than length of data...using unconditional start method...\n")
}
# Random Samples from the Distribution
if(length(sseed) == 1){
zmatrix = data.frame(dist = model$modeldesc$distribution, lambda = ipars[idx["ghlambda",1], 1],
skew = ipars[idx["skew",1], 1], shape = ipars[idx["shape",1], 1], n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(model$modeldesc$distribution, m.sim), lambda = rep(ipars[idx["ghlambda",1], 1], m.sim),
skew = rep(ipars[idx["skew",1], 1], m.sim), shape = rep(ipars[idx["shape",1], 1], m.sim),
n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
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)
}
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
if(startMethod[1] == "unconditional"){
#warning("\nunconditional start method not available for FIGARCH model. Using sample method.")
presigma = tail(sigma, m)
} else{
presigma = tail(sigma, m)
}
}
if(is.na(prereturns[1])){
if(startMethod[1] == "unconditional"){
prereturns = as.numeric(rep(uncmean(fit), m))
}
else{
prereturns = tail(data, m)
}
}
truncLag = fit@model$trunclag
# ebar,eps and be from the filtered is now available
ebar =c(tail(fit@fit$ebar, m),rep(0,n))
eps = fit@fit$eps
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(tail(eps,truncLag+m), rep(0,n))
# input vectors/matrices
h = c(presigma^2, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(ipars[idx["mu",1]:idx["mu",2], 1], ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
ebarSim =matrix(0, ncol = m.sim, nrow = n.sim)
epsSim = matrix(0, ncol = m.sim, nrow = n.sim)
z[is.na(z) | is.nan(z) | !is.finite(z)] = 0
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
if(startMethod[1] == "unconditional"){
preres = as.numeric(z[1:m, i])*presigma
} else{
preres = tail(resids, m)
}
}
res = c(preres, rep(0, n))
ans1 = try(.C("figarchsimC", model = as.integer(modelinc[1:21]), pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
h = as.double(h), z = as.double(z[,i]), res = as.double(res), e = as.double(res*res),
ebar = as.double(ebar), eps = as.double(eps), be = as.double(be),
vexdata = as.double(vexsim[[i]]), T = as.integer(n+m), N = as.integer(truncLag), m = as.integer(m), PACKAGE = "rugarch"), silent = TRUE)
if(inherits(ans1, "try-error")) stop("\nugarchsim-->error: error in calling C function....\n")
sigmaSim[,i] = ans1$h[(n.start + m + 1):(n+m)]^(1/2)
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
ebarSim[,i] = ans1$ebar[(n.start + m + 1):(n+m)]
epsSim[,i] = ans1$eps[(n.start + m + 1):(n+m)]
# ToDo: change to accomodate modelinc[20]
if(modelinc[6]>0){
mxreg = matrix( ipars[idx["mxreg",1]:idx["mxreg",2], 1], ncol = modelinc[6] )
if(modelinc[20]==0){
constm[,i] = constm[,i] + mxreg %*%t( matrix( mexsim[[i]], ncol = modelinc[6] ) )
} else{
if(modelinc[20] == modelinc[6]){
constm[,i] = constm[,i] + mxreg %*%t( matrix( mexsim[[i]]*sqrt(ans1$h), ncol = modelinc[6] ) )
} else{
constm[,i] = constm[,i] + mxreg[,1:(modelinc[6]-modelinc[20]),drop=FALSE] %*%t( matrix( mexsim[[i]][,1:(modelinc[6]-modelinc[20]),drop=FALSE], ncol = modelinc[6]-modelinc[20] ) )
constm[,i] = constm[,i] + mxreg[,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE] %*%t( matrix( mexsim[[i]][,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE]*sqrt(ans1$h), ncol = modelinc[20] ) )
}
}
}
if(modelinc[5]>0) constm[,i] = constm[,i] + ipars[idx["archm",1]:idx["archm",2], 1]*(sqrt(ans1$h)^modelinc[5])
if(modelinc[4]>0){
fres = c(ans1$res[(m+1):(n+m)], if(modelinc[3]>0) rep(0, modelinc[3]) else NULL)
ans2 = .arfimaxsim(modelinc[1:21], ipars, idx, constm[1:n, i], fres, T = n)
seriesSim[,i] = head(ans2$series, n.sim)
} else{
ans2 = .armaxsim(modelinc[1:21], ipars, idx, constm[,i], x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim, ebarSim = ebarSim, epsSim = epsSim, be = be)
model$modeldata$sigma = sigma
sol = new("uGARCHsim",
simulation = sim,
model = model,
seed = as.integer(sseed))
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION sGARCH path
#---------------------------------------------------------------------------------
.figarchpath = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL, trunclag=1000, ...)
{
ans = .figarchpath1(spec = spec, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
presigma = presigma, prereturns = prereturns, preresiduals = preresiduals,
rseed = rseed, custom.dist = custom.dist, mexsimdata = mexsimdata,
vexsimdata = vexsimdata, trunclag=trunclag)
return( ans )
}
.figarchpath1 = function(spec, n.sim = 1000, n.start = 0, m.sim = 1, presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL, trunclag=1000)
{
if(spec@model$modelinc[4]>0){
if(n.start<spec@model$modelinc[3]){
warning("\nugarchpath-->warning: n.start>=MA order for arfima model...automatically setting.")
n.start = spec@model$modelinc[3]
}
}
# some checks
if(is.na(rseed[1])){
sseed = NA
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
model = spec@model
ipars = model$pars
pars = unlist(model$fixed.pars)
parnames = names(pars)
modelnames = .checkallfixed(spec)
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
cat("\nugarchpath-->error: parameters names do not match specification\n")
cat("Expected Parameters are: ")
cat(paste(modelnames))
cat("\n")
stop("Exiting", call. = FALSE)
}
# once more into the spec
setfixed(spec)<-as.list(pars)
model = spec@model
ipars = model$pars
idx = model$pidx
modelinc = model$modelinc
# Enlarge Series:
n = n.sim + n.start
m = model$maxOrder
N = 0
if(modelinc[6]>0) {
mexdata = matrix(model$modeldata$mexdata, ncol = modelinc[6])
N = dim(mexdata)[1]
} else { mexdata = NULL }
if(modelinc[15]>0) {
vexdata = matrix(model$modeldata$vexdata, ncol = modelinc[15])
N = dim(vexdata)[1]
} else { vexdata = NULL }
distribution = model$modeldesc$distribution
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(model, mexsimdata, vexsimdata, N, n, m.sim, m)
mexsim = xreg$mexsimlist
vexsim = xreg$vexsimlist
# check positivity of coefficients
# Random Samples from the Distribution
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, lambda = ipars[idx["ghlambda",1], 1],
skew = ipars[idx["skew",1], 1], shape = ipars[idx["shape",1], 1],
n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), lambda = rep(ipars[idx["ghlambda",1], 1], m.sim),
skew = rep(ipars[idx["skew",1], 1], m.sim), shape = rep(ipars[idx["shape",1], 1], 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)
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchpath-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchpath-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchpath-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
hEst = ipars["omega",1]/(1-ipars["beta1",1]-ipars["alpha1",1])
# we need some estimate to start it off (use the wrong one but not unreasonable)
presigma = as.numeric(rep(sqrt(hEst), times = m))
}
if(is.na(prereturns[1])){
prereturns = as.numeric(rep(uncmean(spec), times = m))
}
truncLag = trunclag
ebar =c(rep(0, m),rep(0,n))
delta = ipars["delta",1]
k=seq(1,truncLag)
be=(k-1-delta)/k
be=rev(cumprod(be))
eps = c(rep(presigma^2,truncLag+m), rep(0,n))
# input vectors/matrices
h = c(presigma^2, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(ipars[idx["mu",1]:idx["mu",2],1], ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
ebarSim =matrix(0, ncol = m.sim, nrow = n.sim)
epsSim = matrix(0, ncol = m.sim, nrow = n.sim)
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
preres = as.numeric(z[1:m,i])*presigma
}
z[1:m, 1:m.sim] = preres[1:m]/presigma[1:m]
z[is.na(z) | is.nan(z) | !is.finite(z)] = 0
res = c(preres, rep(0, n))
ans1 = try(.C("figarchsimC", model = as.integer(modelinc[1:21]), pars = as.double(ipars[,1]), idx = as.integer(idx[,1]-1),
h = as.double(h), z = as.double(z[,i]), res = as.double(res), e = as.double(res*res),
ebar = as.double(ebar), eps = as.double(eps), be = as.double(be),
vexdata = as.double(vexsim[[i]]), T = as.integer(n+m), N = as.integer(truncLag),
m = as.integer(m), PACKAGE = "rugarch"), silent = TRUE)
if(inherits(ans1, "try-error")) stop("\nugarchpath-->error: error in calling C function....\n")
sigmaSim[,i] = ans1$h[(n.start + m + 1):(n+m)]^(1/2)
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
ebarSim[,i] = ans1$ebar[(n.start + m + 1):(n+m)]
epsSim[,i] = ans1$eps[(n.start + m + 1):(n+m)]
if(modelinc[6]>0){
mxreg = matrix( ipars[idx["mxreg",1]:idx["mxreg",2], 1], ncol = modelinc[6] )
if(modelinc[20]==0){
constm[,i] = constm[,i] + as.numeric( mxreg %*%t( matrix( mexsim[[i]], ncol = modelinc[6] ) ) )
} else{
if(modelinc[20] == modelinc[6]){
constm[,i] = constm[,i] + as.numeric( mxreg %*%t( matrix( mexsim[[i]]*sqrt(ans1$h), ncol = modelinc[6] ) ) )
} else{
constm[,i] = constm[,i] + as.numeric( mxreg[,1:(modelinc[6]-modelinc[20]),drop=FALSE] %*%t( matrix( mexsim[[i]][,1:(modelinc[6]-modelinc[20]),drop=FALSE], ncol = modelinc[6]-modelinc[20] ) ) )
constm[,i] = constm[,i] + as.numeric( mxreg[,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE] %*%t( matrix( mexsim[[i]][,(modelinc[6]-modelinc[20]+1):modelinc[6],drop=FALSE]*sqrt(ans1$h), ncol = modelinc[20] ) ) )
}
}
}
if(modelinc[5]>0) constm[,i] = constm[,i] + ipars[idx["archm",1]:idx["archm",2], 1]*(sqrt(ans1$h)^modelinc[5])
if(modelinc[4]>0){
fres = c(ans1$res[(m+1):(n+m)], if(modelinc[3]>0) rep(0, modelinc[3]) else NULL)
ans2 = .arfimaxsim(modelinc[1:21], ipars, idx, constm[1:n, i], fres, T = n)
seriesSim[,i] = head(ans2$series, n.sim)
} else{
#if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .armaxsim(modelinc[1:21], ipars, idx, constm[,i], x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
path = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sol = new("uGARCHpath",
path = path,
model = model,
seed = as.integer(sseed))
return(sol)
}
Try the rugarch package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.