Nothing
R/rgarch-rolling.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.
##
#################################################################################
.rollfdensity = 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$stationarity)) fit.control$stationarity=1
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)
vexdata = 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
model = spec@variance.model$model
forecastlist = vector(mode="list", length = nf)
forecast.length = floor(forecast.length/refit.every) * refit.every
VaR.list = NULL
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
}
if( spec@variance.model$include.vex ){
vexdata[[i]] = spec@variance.model$vexdata[1:(fitindx.start[i]+refit.every), , drop = FALSE]
} else{
vexdata[[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
}
if( spec@variance.model$include.vex ){
vexdata[[i]] = spec@variance.model$vexdata[(i*refit.every):(fitindx.start[i]+refit.every), , drop = FALSE]
} else{
vexdata[[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
specx[[i]]@variance.model$vexdata = if(spec@variance.model$include.vex) vexdata[[i]] else NULL
}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
if(!exists("mclapply")){
require('multicore')
}
fitlist = mclapply(1:nf, FUN = function(i){
ans = ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solver,
out.sample = refit.every, solver.control = solver.control,
fit.control = fit.control)
if(ans@fit$convergence!=0){
if(solver == "solnp") solverx = "nlminb" else solverx = "gosolnp"
ans = ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solverx,
out.sample = refit.every, solver.control = list(),
fit.control = fit.control)
}
return(ans)
}, 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){
ans = rgarch::ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solver,
out.sample = refit.every, solver.control = solver.control,
fit.control = fit.control)
if(ans@fit$convergence!=0){
if(solver == "solnp") solverx = "nlminb" else solverx = "gosolnp"
ans = rgarch::ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solverx,
out.sample = refit.every, solver.control = list(),
fit.control = fit.control)
}
return(ans)
})
sfStop()
}
} else{
for(i in 1:nf){
fitlist[[i]] = ugarchfit(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]] = ugarchfit(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
}
}
}
}
# TODO: return fitlist and allow to resume
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) ugarchforecast(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::ugarchforecast(fitlist[[i]],
n.ahead = n.ahead, n.roll = refit.every-1))
sfStop()
}
} else{
forecastlist = lapply(fitlist, FUN = function(x) ugarchforecast(x, n.ahead = n.ahead, n.roll = refit.every-1))
}
eindex = t(.embed(1:forecast.length, refit.every, refit.every, TRUE))
for(i in 1:nf){
filter.sigma[eindex[,i]] = forecastlist[[i]]@filter$sigma[(fitindx.start[i]
+1):(fitindx.start[i]+refit.every)]
}
# 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[,2])))
fsig = as.numeric(sapply(forecastlist,FUN=function(x) sapply(x@forecast$forecasts, FUN=function(x) x[,1])))
fdates = dates[as.numeric(fwindex)]
rdat = data[(fitindx.start[1]+1):(fitindx.start[1] + forecast.length + 1 - 1)]
fpmlist[[1]] = .forctests3(fsig, fser, fitindx.start[1]+1, n.roll = forecast.length, 1, rdat, fdates,
spec@variance.model$model, spec@variance.model$submodel, 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"])
fmatrix = matrix(NA, ncol = 5, nrow = forecast.length)
colnames(fmatrix) = c("muf", "sigmaf", "dlambdaf", "skewf", "shapef")
fmatrix[,1] = fser
fmatrix[,2] = fsig
for(i in 1:dim(eindex)[2]){
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$garchmodel = model
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$filtersigma = filter.sigma
rolllist$forecast.length = forecast.length
rolllist$VaR.alpha =VaR.alpha
rolllist$dataname = dataname
ans = new(paste(spec@variance.model$model,"roll",sep=""),
roll = rolllist,
forecast = forecastlist,
spec)
} else{
fpmlist = vector(mode="list", length = nf)
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"])
# 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){
fsig = unlist(lapply(forecastlist, FUN=function(x) sapply(x@forecast$forecast, FUN=function(x) x[i,1])))
fser = unlist(lapply(forecastlist, FUN=function(x) sapply(x@forecast$forecast, FUN=function(x) x[i,2])))
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]] = .forctests3(fsig, fser, fitindx.start[1]+i, n.roll = forecast.length, i, rdat, fdates,
spec@variance.model$model, spec@variance.model$submodel, type = 1)
flist[[i]] = fdensity[[i]] = f01density[[i]] = matrix(NA, ncol = 5, nrow = forecast.length)
f01density[[i]][,1] = fser
f01density[[i]][,2] = fsig
for(j in 1:dim(eindex)[2]){
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$garchmodel = model
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$filtersigma = filter.sigma
rolllist$forecast.length = forecast.length
rolllist$VaR.alpha =VaR.alpha
rolllist$dataname = dataname
ans = new(paste(spec@variance.model$model,"roll",sep=""),
roll = rolllist,
forecast = forecastlist,
spec)
}
return(ans)
}
# forecast performance measures
.ugarchrollreport = function(object, type = "VaR", n.ahead = 1, VaR.alpha = 0.01, conf.level = 0.95)
{
switch(type,
VaR = .rollVaRreport(object, n.ahead, VaR.alpha, conf.level),
fpm = .rollfpmreport(object, n.ahead))
invisible(object)
}
.rollfpmreport = function(object, n.ahead = 1)
{
fpmobj = object@roll$fpm[[n.ahead]]
model = fpmobj@details$model
submodel = fpmobj@details$submodel
cat(paste("\nGARCH Roll Forecast Performance Measures", sep = ""))
cat(paste("\n----------------------------------------", sep = ""))
cat(paste("\nModel : ", model, sep = ""))
if(model == "fGARCH"){
cat(paste("\nSubModel : ", submodel, sep = ""))
}
cat(paste("\nno.refits : ", object@roll$n.refit, sep = ""))
cat(paste("\nn.ahead : ", fpmobj@details$n.ahead, sep = ""))
cat(paste("\nn.roll : ", fpmobj@details$n.roll, sep = ""))
#cat(paste("\nfpm type : ", fpmobj@details$type, sep = ""))
cat(paste("\nForecasts w/th Mean Loss Functions\n",n.ahead,"-ahead rolling\n", sep=""))
fc = cbind(fpmobj@forecasts$series, fpmobj@forecasts$sigma)
rownames(fc) = paste("roll-", seq(1, fpmobj@details$n.roll,by = 1), sep="")
fc = rbind(fc, cbind(fpmobj@seriesfpm$meanloss["mse"], fpmobj@sigmafpm$meanloss["mse"]))
fc = rbind(fc, cbind(fpmobj@seriesfpm$meanloss["mad"], fpmobj@sigmafpm$meanloss["mad"]))
fc = rbind(fc, cbind(fpmobj@seriesfpm$meanloss["dac"], NA))
fc = rbind(fc, cbind(NA, fpmobj@sigmafpm$meanloss["r2log"]))
fc = rbind(fc, cbind(NA, fpmobj@sigmafpm$meanloss["qlike"]))
colnames(fc) = c("series", "sigma")
cat("\n")
if(dim(fc)[1]>14){
print(head(round(fc,6),9), digits = 5)
cat(" ........ .......\n")
print(tail(round(fc,6),5), digits = 5)
} else{
print(round(fc,6), digits = 5)
}
cat("\n* performance measures for column 2 are on variance")
cat("\n\n")
}
.rollVaRreport = function(object, n.ahead = 1, VaR.alpha = 0.01, conf.level = 0.95)
{
n = object@roll$n.ahead
v.a = object@roll$VaR.alpha
if(is.null(object@roll$VaR.out)) stop("\nplot-->error: VaR was not calculated for this object\n", call.=FALSE)
if(n.ahead > n) stop("\nplot-->error: n.ahead chosen is not valid for object\n", call.=FALSE)
if(!is.null(v.a) && !any(v.a==VaR.alpha[1])) stop("\nplot-->error: VaR.alpha chosen is invalid for the object\n", call.=FALSE)
if(is.list(object@roll$VaR.out)){
dvar = object@roll$VaR.out[[n.ahead]]
m = dim(dvar)[2]
idx = which(colnames(dvar) == paste("alpha(", round(VaR.alpha,2)*100, "%)",sep=""))
.VaRreport(object@roll$dataname, object@roll$garchmodel, object@roll$distribution, p = VaR.alpha, actual=dvar[,m], VaR = dvar[, idx],
conf.level = conf.level)
} else{
dvar = object@roll$VaR.out
m = dim(dvar)[2]
idx = which(colnames(dvar) == paste("alpha(", round(VaR.alpha,2)*100, "%)",sep=""))
.VaRreport(object@roll$dataname, object@roll$garchmodel, object@roll$distribution, p = VaR.alpha, actual=dvar[,m], VaR = dvar[, idx],
conf.level = conf.level)
}
invisible(object)
}
.embed = function(data, k, by = 1, ascending = FALSE)
{
# n = no of time points, k = number of columns
# by = increment. normally =1 but if =b calc every b-th point
# ascending If TRUE, points passed in ascending order else descending.
# Note that embed(1:n,k) corresponds to embedX(n,k,by=1,rev=TRUE)
# e.g. embedX(10,3)
if(is.null(dim(data)[1])) n<-length(data) else n<-dim(data)[1]
s <- seq(1,n-k+1,by)
lens <- length(s)
cols <- if (ascending) 1:k else k:1
return(matrix(data[s + rep(cols,rep(lens,k))-1],lens))
}
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.
##
#################################################################################
.rollfdensity = 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$stationarity)) fit.control$stationarity=1
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)
vexdata = 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
model = spec@variance.model$model
forecastlist = vector(mode="list", length = nf)
forecast.length = floor(forecast.length/refit.every) * refit.every
VaR.list = NULL
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
}
if( spec@variance.model$include.vex ){
vexdata[[i]] = spec@variance.model$vexdata[1:(fitindx.start[i]+refit.every), , drop = FALSE]
} else{
vexdata[[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
}
if( spec@variance.model$include.vex ){
vexdata[[i]] = spec@variance.model$vexdata[(i*refit.every):(fitindx.start[i]+refit.every), , drop = FALSE]
} else{
vexdata[[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
specx[[i]]@variance.model$vexdata = if(spec@variance.model$include.vex) vexdata[[i]] else NULL
}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
if(!exists("mclapply")){
require('multicore')
}
fitlist = mclapply(1:nf, FUN = function(i){
ans = ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solver,
out.sample = refit.every, solver.control = solver.control,
fit.control = fit.control)
if(ans@fit$convergence!=0){
if(solver == "solnp") solverx = "nlminb" else solverx = "gosolnp"
ans = ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solverx,
out.sample = refit.every, solver.control = list(),
fit.control = fit.control)
}
return(ans)
}, 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){
ans = rgarch::ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solver,
out.sample = refit.every, solver.control = solver.control,
fit.control = fit.control)
if(ans@fit$convergence!=0){
if(solver == "solnp") solverx = "nlminb" else solverx = "gosolnp"
ans = rgarch::ugarchfit(spec = specx[[i]], data = fitdata[[i]], solver = solverx,
out.sample = refit.every, solver.control = list(),
fit.control = fit.control)
}
return(ans)
})
sfStop()
}
} else{
for(i in 1:nf){
fitlist[[i]] = ugarchfit(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]] = ugarchfit(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
}
}
}
}
# TODO: return fitlist and allow to resume
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) ugarchforecast(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::ugarchforecast(fitlist[[i]],
n.ahead = n.ahead, n.roll = refit.every-1))
sfStop()
}
} else{
forecastlist = lapply(fitlist, FUN = function(x) ugarchforecast(x, n.ahead = n.ahead, n.roll = refit.every-1))
}
eindex = t(.embed(1:forecast.length, refit.every, refit.every, TRUE))
for(i in 1:nf){
filter.sigma[eindex[,i]] = forecastlist[[i]]@filter$sigma[(fitindx.start[i]
+1):(fitindx.start[i]+refit.every)]
}
# 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[,2])))
fsig = as.numeric(sapply(forecastlist,FUN=function(x) sapply(x@forecast$forecasts, FUN=function(x) x[,1])))
fdates = dates[as.numeric(fwindex)]
rdat = data[(fitindx.start[1]+1):(fitindx.start[1] + forecast.length + 1 - 1)]
fpmlist[[1]] = .forctests3(fsig, fser, fitindx.start[1]+1, n.roll = forecast.length, 1, rdat, fdates,
spec@variance.model$model, spec@variance.model$submodel, 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"])
fmatrix = matrix(NA, ncol = 5, nrow = forecast.length)
colnames(fmatrix) = c("muf", "sigmaf", "dlambdaf", "skewf", "shapef")
fmatrix[,1] = fser
fmatrix[,2] = fsig
for(i in 1:dim(eindex)[2]){
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$garchmodel = model
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$filtersigma = filter.sigma
rolllist$forecast.length = forecast.length
rolllist$VaR.alpha =VaR.alpha
rolllist$dataname = dataname
ans = new(paste(spec@variance.model$model,"roll",sep=""),
roll = rolllist,
forecast = forecastlist,
spec)
} else{
fpmlist = vector(mode="list", length = nf)
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"])
# 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){
fsig = unlist(lapply(forecastlist, FUN=function(x) sapply(x@forecast$forecast, FUN=function(x) x[i,1])))
fser = unlist(lapply(forecastlist, FUN=function(x) sapply(x@forecast$forecast, FUN=function(x) x[i,2])))
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]] = .forctests3(fsig, fser, fitindx.start[1]+i, n.roll = forecast.length, i, rdat, fdates,
spec@variance.model$model, spec@variance.model$submodel, type = 1)
flist[[i]] = fdensity[[i]] = f01density[[i]] = matrix(NA, ncol = 5, nrow = forecast.length)
f01density[[i]][,1] = fser
f01density[[i]][,2] = fsig
for(j in 1:dim(eindex)[2]){
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$garchmodel = model
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$filtersigma = filter.sigma
rolllist$forecast.length = forecast.length
rolllist$VaR.alpha =VaR.alpha
rolllist$dataname = dataname
ans = new(paste(spec@variance.model$model,"roll",sep=""),
roll = rolllist,
forecast = forecastlist,
spec)
}
return(ans)
}
# forecast performance measures
.ugarchrollreport = function(object, type = "VaR", n.ahead = 1, VaR.alpha = 0.01, conf.level = 0.95)
{
switch(type,
VaR = .rollVaRreport(object, n.ahead, VaR.alpha, conf.level),
fpm = .rollfpmreport(object, n.ahead))
invisible(object)
}
.rollfpmreport = function(object, n.ahead = 1)
{
fpmobj = object@roll$fpm[[n.ahead]]
model = fpmobj@details$model
submodel = fpmobj@details$submodel
cat(paste("\nGARCH Roll Forecast Performance Measures", sep = ""))
cat(paste("\n----------------------------------------", sep = ""))
cat(paste("\nModel : ", model, sep = ""))
if(model == "fGARCH"){
cat(paste("\nSubModel : ", submodel, sep = ""))
}
cat(paste("\nno.refits : ", object@roll$n.refit, sep = ""))
cat(paste("\nn.ahead : ", fpmobj@details$n.ahead, sep = ""))
cat(paste("\nn.roll : ", fpmobj@details$n.roll, sep = ""))
#cat(paste("\nfpm type : ", fpmobj@details$type, sep = ""))
cat(paste("\nForecasts w/th Mean Loss Functions\n",n.ahead,"-ahead rolling\n", sep=""))
fc = cbind(fpmobj@forecasts$series, fpmobj@forecasts$sigma)
rownames(fc) = paste("roll-", seq(1, fpmobj@details$n.roll,by = 1), sep="")
fc = rbind(fc, cbind(fpmobj@seriesfpm$meanloss["mse"], fpmobj@sigmafpm$meanloss["mse"]))
fc = rbind(fc, cbind(fpmobj@seriesfpm$meanloss["mad"], fpmobj@sigmafpm$meanloss["mad"]))
fc = rbind(fc, cbind(fpmobj@seriesfpm$meanloss["dac"], NA))
fc = rbind(fc, cbind(NA, fpmobj@sigmafpm$meanloss["r2log"]))
fc = rbind(fc, cbind(NA, fpmobj@sigmafpm$meanloss["qlike"]))
colnames(fc) = c("series", "sigma")
cat("\n")
if(dim(fc)[1]>14){
print(head(round(fc,6),9), digits = 5)
cat(" ........ .......\n")
print(tail(round(fc,6),5), digits = 5)
} else{
print(round(fc,6), digits = 5)
}
cat("\n* performance measures for column 2 are on variance")
cat("\n\n")
}
.rollVaRreport = function(object, n.ahead = 1, VaR.alpha = 0.01, conf.level = 0.95)
{
n = object@roll$n.ahead
v.a = object@roll$VaR.alpha
if(is.null(object@roll$VaR.out)) stop("\nplot-->error: VaR was not calculated for this object\n", call.=FALSE)
if(n.ahead > n) stop("\nplot-->error: n.ahead chosen is not valid for object\n", call.=FALSE)
if(!is.null(v.a) && !any(v.a==VaR.alpha[1])) stop("\nplot-->error: VaR.alpha chosen is invalid for the object\n", call.=FALSE)
if(is.list(object@roll$VaR.out)){
dvar = object@roll$VaR.out[[n.ahead]]
m = dim(dvar)[2]
idx = which(colnames(dvar) == paste("alpha(", round(VaR.alpha,2)*100, "%)",sep=""))
.VaRreport(object@roll$dataname, object@roll$garchmodel, object@roll$distribution, p = VaR.alpha, actual=dvar[,m], VaR = dvar[, idx],
conf.level = conf.level)
} else{
dvar = object@roll$VaR.out
m = dim(dvar)[2]
idx = which(colnames(dvar) == paste("alpha(", round(VaR.alpha,2)*100, "%)",sep=""))
.VaRreport(object@roll$dataname, object@roll$garchmodel, object@roll$distribution, p = VaR.alpha, actual=dvar[,m], VaR = dvar[, idx],
conf.level = conf.level)
}
invisible(object)
}
.embed = function(data, k, by = 1, ascending = FALSE)
{
# n = no of time points, k = number of columns
# by = increment. normally =1 but if =b calc every b-th point
# ascending If TRUE, points passed in ascending order else descending.
# Note that embed(1:n,k) corresponds to embedX(n,k,by=1,rev=TRUE)
# e.g. embedX(10,3)
if(is.null(dim(data)[1])) n<-length(data) else n<-dim(data)[1]
s <- seq(1,n-k+1,by)
lens <- length(s)
cols <- if (ascending) 1:k else k:1
return(matrix(data[s + rep(cols,rep(lens,k))-1],lens))
}
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.