Nothing
R/arfima-main.R
In rgarch: Flexible GARCH modelling in R
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
## This file is part of the R package rgarch.
##
## The R package rgarch is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## The R package rgarch is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
#################################################################################
#---------------------------------------------------------------------------------
# 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 2, 2019, 5:22 p.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.
#################################################################################
##
## 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.
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.