Nothing
#################################################################################
##
## 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.
##
#################################################################################
# CONSTANTS used in package:
eps<-.Machine$double.eps
TinY = 1.0e-8
##################################################################################
# Helper Functions
##################################################################################
.makemodel = function(garchenv)
{
model = vector(mode="list")
mmodel = get("mmodel", garchenv)
vmodel = get("vmodel", garchenv)
omodel = get("omodel", garchenv)
dmodel = get("dmodel", garchenv)
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
model$fixed = fixed
model$fixid = fixid
model$model = vmodel$model
model$submodel = vmodel$submodel
model$garchOrder = vmodel$garchOrder
model$include.mex = mmodel$include.mex
model$mexdata = mmodel$origmexdata
model$vexdata = vmodel$origvexdata
model$mxn = mmodel$mxn
model$vxn = vmodel$vxn
model$include.vex = vmodel$include.vex
model$armaOrder = mmodel$armaOrder
model$garchInMean = mmodel$garchInMean
model$arfima = mmodel$arfima
model$include.mex = mmodel$include.mex
model$inMeanType = mmodel$inMeanType
model$im = mmodel$im
model$include.mean = mmodel$include.mean
model$include.omega = dmodel$include.omega
model$distribution = dmodel$distribution
model$include.dlambda = dmodel$include.dlambda
model$include.skew = dmodel$include.skew
model$include.shape = dmodel$include.shape
model$modelnames = omodel$modelnames
model$start.pars = omodel$start.pars
model$fixed.pars = omodel$fixed.pars
model$distribution = dmodel$distribution
model$modelmatrix = omodel$modelmatrix
model$maxOrder = omodel$maxOrder
model$maxOrder2 = omodel$maxOrder2
model$pos = omodel$pos.matrix
return(model)
}
.makearfimamodel = function(garchenv)
{
model = vector(mode="list")
mmodel = get("mmodel", garchenv)
omodel = get("omodel", garchenv)
dmodel = get("dmodel", garchenv)
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
model$fixed = fixed
model$fixid = fixid
model$include.mex = mmodel$include.mex
model$mexdata = mmodel$origmexdata
model$mxn = mmodel$mxn
model$armaOrder = mmodel$armaOrder
model$arfima = mmodel$arfima
model$include.mex = mmodel$include.mex
model$include.mean = mmodel$include.mean
model$distribution = dmodel$distribution
model$include.dlambda = dmodel$include.dlambda
model$include.skew = dmodel$include.skew
model$include.shape = dmodel$include.shape
model$modelnames = omodel$modelnames
model$start.pars = omodel$start.pars
model$fixed.pars = omodel$fixed.pars
model$distribution = dmodel$distribution
model$modelmatrix = omodel$modelmatrix
model$maxOrder = omodel$maxOrder
model$pos = omodel$pos.matrix
return(model)
}
.makearfimafitmodel = function(f, opt, data, T, m, timer, convergence, message, hess, garchenv)
{
fit = vector(mode = "list")
model = vector(mode="list")
mmodel = get("mmodel", garchenv)
omodel = get("omodel", garchenv)
dmodel = get("dmodel", garchenv)
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
npar = get("npar", garchenv)
assign(".llh",1, envir = garchenv)
if(is.null(hess)){
#fit$hessian = .hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)
fit$hessian = hessian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "llh", garchenv = garchenv)
} else{
fit$hessian = hess
}
fit$cvar = try(solve(fit$hessian), silent = TRUE)
# (might also fail in which case user will see an error about the inversion failure)
if(inherits(fit$cvar, "try-error")){
zz = try(solve(.hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)), silent=TRUE)
if(inherits(zz, "try-error")) {
fit$cvar = NULL
cat("\nrgarch-->warning: failed to invert hessian\n")
} else{
fit$cvar = zz
}
}
#assign(".llh", 1, envir = garchenv)
# garchenv)
# call LLH once more to get optimal parameters not report with returnType "all"
temp = f(pars = opt, data = data, returnType="all", garchenv)
fit$z = temp$z
fit$data = data
fit$dates = get("dates", garchenv)
fit$LLH = -temp$llh
fit$log.likelihoods = temp$LHT
fit$residuals = temp$epsx
assign(".llh",1, envir = garchenv)
if(!is.null(fixed)){
pall = vector(mode = "numeric", length = npar)
nfix = length(fixed)
pall[fixid] = fixed
pall[-fixid] = opt
fit$coef = pall
Names = omodel$modelnames
names(fit$coef) = Names
complNULL = rep(NA, nfix)
if(is.null(fit$cvar)){
fit$se.coef = rep(NA,length(opt))
fit$tval = rep(NA,length(opt))
# change here
fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
fit$tval, rep(NA,length(opt)))
fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$robust.se.coef = rep(NA,length(opt))
fit$robust.tval = rep(NA,length(opt))
# change here
fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef,
fit$robust.tval, rep(NA,length(opt)))
fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$hessian.message = "failed to invert hessian"
} else{
tmp = robustvcv(fun = f, pars = opt, nlag = 0, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
fit$robust.cvar = tmp$vcv
fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv)
fit$se.coef = sqrt(diag(abs(fit$cvar)))
fit$tval = fit$coef[-fixid]/fit$se.coef
# change here
fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
fit$tval, 2*(1-pnorm(abs(fit$tval))))
fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
fit$robust.tval = fit$coef[-fixid]/fit$robust.se.coef
# change here
fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef,
fit$robust.tval, 2*(1-pnorm(abs(fit$robust.tval))))
fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$hessian.message = NULL
}
} else{
fit$coef = opt
names(fit$coef) = omodel$modelnames
if(is.null(fit$cvar)){
fit$se.coef = rep(NA,length(opt))
fit$tval = rep(NA,length(opt))
# change here
fit$matcoef = cbind(fit$coef, fit$se.coef,
fit$tval, rep(NA,length(opt)))
fit$robust.se.coef = rep(NA,length(opt))
fit$robust.tval = rep(NA,length(opt))
# change here
fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval,
rep(NA,length(opt)))
fit$hessian.message = "failed to invert hessian"
} else{
nlag=min(floor(1.2*(T)^(1/3)),(T))
tmp = robustvcv(fun = f, pars = opt, nlag = nlag, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
fit$robust.cvar = tmp$vcv
fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv)
fit$se.coef = sqrt(diag(abs(fit$cvar)))
fit$tval = fit$coef/fit$se.coef
# change here
fit$matcoef = cbind(fit$coef, fit$se.coef,
fit$tval, 2*(1-pnorm(abs(fit$tval))))
fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
fit$robust.tval = fit$coef/fit$robust.se.coef
# change here
fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval,
2*(1-pnorm(abs(fit$robust.tval))))
fit$hessian.message = NULL
}
}
dimnames(fit$matcoef) = list(names(fit$coef), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
dimnames(fit$robust.matcoef) = list(names(fit$coef), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
fit$fitted.values = data - fit$residuals
fit$convergence = convergence
fit$message = message
fit$timer = timer
return(fit)
}
.makefitmodel = function(garchmodel, f, opt, data, T, m, timer, convergence, message, hess,
garchenv)
{
fit = vector(mode = "list")
model = vector(mode="list")
mmodel = get("mmodel", garchenv)
vmodel = get("vmodel", garchenv)
omodel = get("omodel", garchenv)
dmodel = get("dmodel", garchenv)
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
npar = get("npar", garchenv)
assign(".llh",1, envir = garchenv)
if(is.null(hess)){
#fit$hessian = .hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)
fit$hessian = hessian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "llh", garchenv = garchenv)
} else{
fit$hessian = hess
}
fit$cvar = try(solve(fit$hessian), silent = TRUE)
# (might also fail in which case user will see an error about the inversion failure)
if(inherits(fit$cvar, "try-error")){
zz = try(solve(.hessian2sided(f, opt, data = data, returnType = "llh", garchenv = garchenv)), silent=TRUE)
if(inherits(zz, "try-error")) {
fit$cvar = NULL
cat("\nrgarch-->warning: failed to invert hessian\n")
} else{
fit$cvar = zz
}
}
#assign(".llh", 1, envir = garchenv)
# garchenv)
# call LLH once more to get optimal parameters not report with returnType "all"
temp = f(pars = opt, data = data, returnType="all", garchenv)
if(garchmodel == "fGARCH" || garchmodel == "apARCH"){
# apARCH and fGARCH are calculated in sigma rather than sigma^d (since d is a variable
# to be calculated, unlike sGARCH where it is fixed at 2)
fit$var = temp$h^2
fit$sigma = temp$h
} else{
fit$var = abs(temp$h)
fit$sigma = sqrt(abs(temp$h))
}
fit$z = temp$z
fit$data = data
fit$dates = get("dates", garchenv)
fit$LLH = -temp$llh
fit$log.likelihoods = temp$LHT
fit$residuals = temp$epsx
assign(".llh",1, envir = garchenv)
if(!is.null(fixed)){
pall = vector(mode = "numeric", length = npar)
nfix = length(fixed)
pall[fixid] = fixed
pall[-fixid] = opt
fit$coef = pall
Names = omodel$modelnames
names(fit$coef) = Names
complNULL = rep(NA, nfix)
if(is.null(fit$cvar)){
fit$se.coef = rep(NA,length(opt))
fit$tval = rep(NA,length(opt))
# change here
fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
fit$tval, rep(NA,length(opt)))
fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$robust.se.coef = rep(NA,length(opt))
fit$robust.tval = rep(NA,length(opt))
# change here
fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef,
fit$robust.tval, rep(NA,length(opt)))
fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$hessian.message = "failed to invert hessian"
} else{
tmp = robustvcv(fun = f, pars = opt, nlag = 0, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
fit$robust.cvar = tmp$vcv
fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv)
fit$se.coef = sqrt(diag(abs(fit$cvar)))
fit$tval = fit$coef[-fixid]/fit$se.coef
# change here
fit$matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$se.coef,
fit$tval, 2*(1-pnorm(abs(fit$tval))))
fit$matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
fit$robust.tval = fit$coef[-fixid]/fit$robust.se.coef
# change here
fit$robust.matcoef = matrix(NA, ncol = 4, nrow = length(fit$coef))
fit$robust.matcoef[-fixid,] = cbind(fit$coef[-fixid], fit$robust.se.coef,
fit$robust.tval, 2*(1-pnorm(abs(fit$robust.tval))))
fit$robust.matcoef[fixid,] = cbind(fit$coef[fixid], complNULL,
complNULL, complNULL)
fit$hessian.message = NULL
}
} else{
fit$coef = opt
names(fit$coef) = omodel$modelnames
if(is.null(fit$cvar)){
fit$se.coef = rep(NA,length(opt))
fit$tval = rep(NA,length(opt))
# change here
fit$matcoef = cbind(fit$coef, fit$se.coef,
fit$tval, rep(NA,length(opt)))
fit$robust.se.coef = rep(NA,length(opt))
fit$robust.tval = rep(NA,length(opt))
# change here
fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval,
rep(NA,length(opt)))
fit$hessian.message = "failed to invert hessian"
} else{
nlag=min(floor(1.2*(T)^(1/3)),(T))
# tmp = .robustSE(fun = f, pars = opt, hess = NULL, n = T-m, data = data, returnType = "LHT", garchenv)
tmp = robustvcv(fun = f, pars = opt, nlag = nlag, hess = fit$hessian, n = T, data = data, returnType = "LHT", garchenv)
fit$robust.cvar = tmp$vcv
fit$scores = jacobian(func = f, x = opt, method="Richardson", method.args=list(), data = data, returnType = "LHT", garchenv)
fit$se.coef = sqrt(diag(abs(fit$cvar)))
fit$tval = fit$coef/fit$se.coef
# change here
fit$matcoef = cbind(fit$coef, fit$se.coef,
fit$tval, 2*(1-pnorm(abs(fit$tval))))
fit$robust.se.coef = sqrt(diag(fit$robust.cvar))
fit$robust.tval = fit$coef/fit$robust.se.coef
# change here
fit$robust.matcoef = cbind(fit$coef, fit$robust.se.coef,fit$robust.tval,
2*(1-pnorm(abs(fit$robust.tval))))
fit$hessian.message = NULL
}
}
dimnames(fit$matcoef) = list(names(fit$coef), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
dimnames(fit$robust.matcoef) = list(names(fit$coef), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
fit$fitted.values = data-fit$residuals
fit$convergence = convergence
fit$message = message
fit$kappa = temp$kappa
fit$persistence = temp$persistence
fit$timer = timer
return(fit)
}
.forcregressors = function(include.mex, mxn, mexdata, mregfor, include.vex,
vxn, vexdata, vregfor, pars, n.ahead, N, out.sample, n.roll)
{
# N is the original length
treq = n.ahead + n.roll
if(include.mex){
if(!is.null(mregfor)){
nmex = dim(as.matrix(mregfor))[1]
mmex = dim(as.matrix(mregfor))[2]
} else{
nmex = 0
mmex = 0
}
if(!is.null(mregfor) && mmex != mxn)
{
cat("\nugarchforecast-->error: Column dimension of external mean forecast matrix is wrong.")
cat(paste("\nModel has ", mxn, " external regressors but forecast matrix has", mmex, sep = ""))
stop("\n...exiting\n")
}
if(!is.null(mregfor) && nmex < treq)
{
cat(paste("\nugarchforecast-->error: You requested ", treq ," actual forecasts (including
the rolling periods) but external mean forecasts provided have only",
nmex, " rows",sep=""))
cat(paste("\nA minimum of ", treq," rows are required", sep=""))
stop("\n...exiting")
}
if(is.null(mregfor)){
mxf = rbind(as.matrix(mexdata)[1:(N-out.sample), ,drop = FALSE], matrix(0, ncol = mxn, nrow = treq))
} else {
mxf = rbind(as.matrix(mexdata)[1:(N-out.sample), ,drop = FALSE], as.matrix(mregfor)[1:treq,,drop = FALSE])
}
mxreg = pars[paste("mxreg", 1:mxn, sep = "")]
} else{
mxreg = NULL
mxf = NULL
}
if(include.vex){
if(!is.null(vregfor)){
nvex = dim(as.matrix(vregfor))[1]
mvex = dim(as.matrix(vregfor))[2]
} else{
nvex = 0
mvex = 0
}
if(!is.null(vregfor) && mvex != vxn)
{
cat("\nugarchforecast-->error: Column dimension of external variance forecast matrix is wrong.")
cat(paste("\nModel has ",vxn," external regressors but forecast matrix has", mvex, sep = ""))
stop("\n...exiting\n")
}
# N is the original length
if(!is.null(vregfor) && nvex < treq)
{
cat(paste("\nugarchforecast-->error: You requested ", treq ," actual forecasts (including
the rolling periods) but external variance forecasts provided have only",
nvex, " rows",sep=""))
cat(paste("\nA minimum of ", treq," rows are required", sep=""))
stop("\n...exiting")
}
if(is.null(vregfor)){
vxf = rbind(as.matrix(vexdata)[1:(N-out.sample), ,drop = FALSE], matrix(0, ncol = vxn, nrow = treq))
} else {
vxf = rbind(as.matrix(vexdata)[1:(N-out.sample), ,drop = FALSE], as.matrix(vregfor)[1:treq,,drop = FALSE])
}
vxreg = pars[paste("vxreg", 1:vxn, sep = "")]
} else{
vxreg = NULL
vxf = NULL
}
return(list(mxreg = mxreg, mxf = mxf, vxreg = vxreg, vxf = vxf))
}
.simregressors = function(include.mex, mxn, mexsimdata, mexdata, include.vex, vxn,
vexsimdata, vexdata, N, pars, n, m.sim, m)
{
if(include.mex){
if(is.null(mexsimdata)){
mexsimdata = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = n)
}
if(!is.null(mexsimdata))
{
if(!is.list(mexsimdata))
stop("\nugarchsim-->error: mexsimdata should be a list of length m.sim")
if(length(mexsimdata) != m.sim){
msd = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) msd[[i]] = as.matrix(mexsimdata[[1]])
mexsimdata = msd
cat("\nugarchsim-->warning: length of mexsimdata list not equal to m.sim...
replicating first list element m.sim times.\n")
}
for(i in 1:m.sim){
if(dim(as.matrix(mexsimdata[[i]]))[2] != mxn )
stop(paste("\nugarchsim-->error: mexsimdata ", i," has wrong no. of column", sep=""))
if(dim(as.matrix(mexsimdata[[i]]))[1] != n )
stop(paste("\nugarchsim-->error: mexsimdata ", i," has wrong no. of rows", sep=""))
}
}
mxreg = pars[paste("mxreg",1:mxn,sep="")]
mexsimlist = vector(mode = "list", length = m.sim)
for(i in 1:m.sim){
premexdata = mexdata[(N-m+1):N,,drop=FALSE]
mexsimlist[[i]] = matrix(rbind(premexdata, mexsimdata[[i]]), ncol = mxn)
}
} else{
mxreg = 0
mexsimlist = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) mexsimlist[[i]]=0
}
if(include.vex){
if(is.null(vexsimdata)){
vexsimdata = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) vexsimdata[[i]] = matrix(0, ncol = vxn, nrow = n)
}
if(!is.null(vexsimdata))
{
if(!is.list(vexsimdata))
stop("\nugarchsim-->error: mexsimdata should be a list of length m.sim")
if(length(vexsimdata) != m.sim){
vsd = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) vsd[[i]] = as.matrix(vexsimdata[[1]])
vexsimdata = vsd
cat("\nugarchsim-->warning: length of vexsimdata list not equal to m.sim...
replicating first list element m.sim times.\n")
}
for(i in 1:m.sim){
if(dim(as.matrix(vexsimdata[[i]]))[2] != vxn )
stop(paste("\nugarchsim-->error: vexsimdata ", i," has wrong no. of column", sep=""))
if(dim(as.matrix(vexsimdata[[i]]))[1] != n )
stop(paste("\nugarchsim-->error: vexsimdata ", i," has wrong no. of rows", sep=""))
}
}
vxreg = pars[paste("vxreg",1:vxn,sep="")]
vexsimlist = vector(mode = "list", length = m.sim)
for(i in 1:m.sim){
prevexdata = vexdata[(N-m+1):N,,drop=FALSE]
vexsimlist[[i]] = matrix(rbind(prevexdata, vexsimdata[[i]]), ncol = vxn)
}
} else{
vxreg = 0
vexsimlist = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) vexsimlist[[i]]=0
}
return(list(mxreg = mxreg, mexsimlist = mexsimlist, vxreg = vxreg,
vexsimlist = vexsimlist))
}
.custzdist = function(custom.dist, zmatrix, m.sim, n)
{
if(is.na(custom.dist$name) | is.na(custom.dist$name)){
z = matrix(apply(zmatrix, 1, FUN = function(x) .makeSample(as.character(x[1]), lambda = as.numeric(x[2]),
skew = as.numeric(x[3]), shape = as.numeric(x[4]), n = as.numeric(x[5]),
seed = as.integer(x[6]))), n, m.sim)
#nx = dim(zmatrix)[1]
#tmp = apply(as.data.frame(1:nx), 1, FUN = function(i){
# set.seed(zmatrix[i,6]);
# .C("distributionsample",
# shape = as.double(zmatrix[i,4]), skew = as.double(zmatrix[i,3]),
# lambda = as.double(zmatrix[i,2]),
# ndis = as.integer(zmatrix[i,1]), n = as.integer(zmatrix[i,5]),
# rvec = double(zmatrix[i,5]), DUP = FALSE, PACKAGE ="rgarch")$rvec})
#z = matrix(tmp, nrow = n, ncol = m.sim)
}
if(!is.na(custom.dist$name) && !is.na(custom.dist$distfit)){
if(is.matrix(custom.dist$distfit))
{
if(dim(custom.dist$distfit)[2]!=m.sim) stop("column dimension of custom
innovations\n matrix must be equal to m.sim")
if(dim(custom.dist$distfit)[1]!=n) stop("row dimension
of custom innovations\n matrix must be equal to n.sim+n.start")
z = custom.dist$distfit
} else{
if(!is.character(custom.dist$name)) stop("custom distribution must be a
character string")
temp = paste("r", custom.dist$name, sep="")
if(is.null(custom.dist$distfit)) stop("custom distribution missing a
distfit object")
# if this is often used we might consider using apply to
# use a new seed for each m.sim
.rdist = eval(parse(text=paste(temp)))
tmp = .rdist(n*m.sim, custom.dist$distfit)
z = matrix(as.numeric(tmp), ncol = m.sim, nrow = n, byrow=TRUE)
}
}
return(z)
}
.stars = function(testvector,levels=c(0.01, 0.05, 0.1))
{
N = length(testvector)
ans = vector(mode="character",length=N)
#recursive replacement
z = which(testvector<levels[3])
ans[z] = c("*")
z = which(testvector<levels[2])
ans[z] = c("**")
z = which(testvector<levels[1])
ans[z] = c("***")
ans
}
# method from a spec and data object
.safefit = function(spec, data, out.sample, solver, fit.control, solver.control)
{
ans = try(ugarchfit(spec = spec, data = data, out.sample = out.sample, fit.control = fit.control,
solver = solver, solver.control = solver.control), silent = TRUE)
if(inherits(ans, "try-error")) ans = NULL
return(ans)
}
.safefitarfima = function(spec, data, out.sample, solver, fit.control, solver.control)
{
ans = try(arfimafit(spec = spec, data = data, out.sample = out.sample, fit.control = fit.control,
solver = solver, solver.control = solver.control), silent = TRUE)
if(inherits(ans, "try-error")) ans = NULL
return(ans)
}
.boxcoxtransform = function(x, lambda)
{
if(lambda!=0) ret = (x^lambda - 1)/lambda else ret = log(x)
return(ret)
}
repmat = function(a, n, m)
{
kronecker(matrix(1, n, m), a)
}
size = function(x, n = NULL)
{
x = as.matrix(x)
if(missing(n)) sol = c(n = dim(x)[1], m = dim(x)[2]) else sol = dim(x)[n]
return(sol)
}
zeros = function(n = 1, m = 1)
{
if(missing(m)) m = n
sol = matrix(0, nrow = n, ncol = m)
return(sol)
}
ones = function(n = 1, m = 1)
{
if(missing(m)) m = n
sol = matrix(1, nrow = n, ncol = m)
return(sol)
}
newlagmatrix = function(x,nlags,xc)
{
nlags = nlags+1
xt = size(x, 1);
newX = rbind(x, zeros(nlags, 1))
lagmatrix = repmat(newX, nlags, 1)
lagmatrix = matrix(lagmatrix[1:(size(lagmatrix,1)-nlags)], nrow = (xt+nlags-1), ncol = nlags)
lagmatrix = lagmatrix[nlags:xt,]
y = lagmatrix[,1]
x = lagmatrix[,2:nlags]
if(xc == 1) x = cbind(ones(size(x,1), 1), x)
return(data.frame(y = y, x = x))
}
.colorgradient = function(n = 50, low.col = 0.6, high.col=0.7, saturation = 0.8) {
if (n < 2) stop("n must be greater than 2")
n1 = n%/%2
n2 = - n - n1
c(hsv(low.col, saturation, seq(1, 0.5, length = n1)),
hsv(high.col, saturation, seq(0.5, 1, length = n2)))
}
.simlayout = function(m)
{
if(m == 1){
nf = c(1, 1, 2, 2)
nf = layout(matrix(nf, 2, 2, byrow = TRUE), respect = TRUE)
middle.plot = 1
}
if(m == 2){
nf = c(1, 1, 1, 1, 0, 2, 2, 0, 3, 3, 3, 3)
nf = layout(matrix(nf, 3, 4, byrow = TRUE), respect = TRUE)
middle.plot = 2
}
if(m == 3){
nf = c(1, 0, 0, 2, 0, 3, 3, 0, 4, 0, 0, 5)
nf = layout(matrix(nf, 3, 4, byrow = TRUE), respect = TRUE)
middle.plot = 3
}
if(m == 12){
nf = c(1, 2, 3, 4,
5, 6, 6, 7,
8, 6, 6, 9,
10, 11, 12, 13)
nf = layout(matrix(nf, 4, 4, byrow = TRUE), respect = TRUE)
}
}
.sdigit = function(x){
sid = as.numeric(strsplit(format(as.numeric(x), scientific=TRUE), "e")[[1]])[2]
10^(-sid)
}
.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]
data = matrix(data,ncol=1)
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))
}
.lagx = function(data, n.lag=1, removeNA = FALSE, pad = NA)
{
# has NAs
data = as.matrix(data)
n = dim(data)[1]
d = dim(data)[2]
if(dim(data)[2]==1) data=matrix(data,ncol=1)
z = apply(data,2,FUN=function(x) .embed(x,n.lag+1)[,n.lag+1])
if(!removeNA) z = rbind(matrix(pad,ncol=d,nrow=n.lag),z)
return(z)
}
.lagmatrix = function(data, n.lag = 1, pad = 0)
{
n = length(as.numeric(data))
z = matrix(NA, ncol = n.lag, nrow = n)
for(i in 1:n.lag) z[,i] = .lagx(as.numeric(data), i, removeNA = FALSE, pad = pad)
return(z)
}
.abind = function(x, y)
{
m = dim(x)[1]
if( is.matrix(x) ) {
n1 = 1
x = array(x, dim = c(m, m, 1))
} else{
n1 = dim(x)[3]
}
if( is.matrix(y) ) {
n2 = 1
y = array(y, dim = c(m, m, 1))
} else{
n2 = dim(y)[3]
}
nw = array(NA, dim = c(m, m, n1 + n2))
nw[ , , 1:n1] = x[, , 1:n1]
nw[ , , (n1+1):(n1+n2)] = y[, , 1:n2]
nw
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.