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.
##
#################################################################################
.eps = .Machine$double.eps
.dcchessian = function(f, pars, garchpars, N, data, garchenv, fname)
{
parallel = get("parallel", garchenv)
parallel.control = get("parallel.control", garchenv)
x = c(garchpars, pars)
n = length(x)
dccn = length(pars)
garchn = length(garchpars)
fx = f(pars, garchpars, data, returnType = "llh", garchenv = garchenv)
.eps = .Machine$double.eps
# Compute the stepsize (h)
h = .eps ^ (1/3) * pmax( abs( x ), 1 )
xh = x + h
h = xh - x
ee = Matrix( diag( h ), sparse = TRUE )
# Compute forward and backward steps
g = vector(mode = "numeric", length = n)
if( parallel ){
if( parallel.control$pkg == "multicore" ){
tmp = mclapply(1:n, FUN = function(i){
cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
tmpg = f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i],
garchpars = x[1:garchn] + ee[1:garchn, i],
data = data, returnType = "llh", garchenv = garchenv)
return( tmpg )
}, mc.cores = parallel.control$cores)
g = as.numeric( unlist(tmp) )
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("x", "garchn", "ee", "data", "garchenv", "n", "fname", local = TRUE)
tmp = sfLapply(as.list(1:n), fun = function(i){
cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
tmpg = eval(parse(text = paste("rgarch:::", fname,
"( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i], garchpars = x[1:garchn] + ee[1:garchn, i],
data = data, returnType = 'llh', garchenv = garchenv)", sep = "")))
return( tmpg )
})
# sfStop() we'll stop it further down
g = as.numeric( unlist(tmp) )
}
} else{
tmp = lapply(as.list(1:n), FUN = function(i){
cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
tmpg = f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i],
garchpars = x[1:garchn] + ee[1:garchn, i],
data = data, returnType = "llh", garchenv = garchenv)
return( tmpg )
})
g = as.numeric( unlist(tmp) )
}
H = h %*% t( h )
#idx = 1
# Compute "double" forward and backward steps
#cat("\n")
#for(i in 1:n){
# for(j in (n - N + 1):n){
# cat(paste("Evaluating Function ",idx," out of ",n*N,"\n",sep=""))
# if(i <= j){
# H[i, j] = (f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j],
# garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j],
# data = data, returnType = "llh", garchenv = garchenv) - g[i] - g[j] + fx) / H[i, j]
# H[j, i] = H[i, j]
# }
# idx = idx + 1
# }
#}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
tmp = mclapply(1:n, FUN = function(i){
Htmp = H
for(j in (n - N + 1):n){
if(i <= j){
Htmp[i, j] = (f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j],
garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j],
data = data, returnType = "llh", garchenv = garchenv) - g[i] - g[j] + fx) / Htmp[i, j]
Htmp[j, i] = Htmp[i, j]
}
}
return(Htmp)
})
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("x", "H", "garchn", "ee", "data", "garchenv", "n", "N", "g", "fx", "fname", local = TRUE)
tmp = sfLapply(as.list(1:n), fun = function(i){
Htmp = H
for(j in (n - N + 1):n){
if(i <= j){
Htmp[i, j] = eval(parse(text = paste("(rgarch:::", fname, "( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j],
garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j],
data = data, returnType = 'llh', garchenv = garchenv) - g[i] - g[j] + fx) / Htmp[i, j]", sep = "")))
Htmp[j, i] = Htmp[i, j]
}
}
return(Htmp)
})
sfStop()
}
for(i in 1:n){
for(j in (n - N + 1):n){
if(i <= j){
H[i, j] = tmp[[i]][i, j]
H[j, i] = tmp[[i]][j, i]
}
}
}
} else{
tmp = lapply(as.list(1:n), FUN = function(i){
Htmp = H
for(j in (n - N + 1):n){
if(i <= j){
Htmp[i, j] = (f( pars = x[-c(1:garchn)] + ee[-c(1:garchn), i] + ee[-c(1:garchn), j],
garchpars = x[1:garchn] + ee[1:garchn, i] + ee[1:garchn, j],
data = data, returnType = "llh", garchenv = garchenv) - g[i] - g[j] + fx) / Htmp[i, j]
Htmp[j, i] = Htmp[i, j]
}
}
return(Htmp)
})
for(i in 1:n){
for(j in (n - N + 1):n){
if(i <= j){
H[i, j] = tmp[[i]][i, j]
H[j, i] = H[i, j]
}
}
}
}
newH = H[(n - N + 1):n, ]
H = newH
return(H)
}
# 2-stage partitioned standard errors for DCC type models
.dccmakefitmodel = function(garchmodel, f, pars, garchpars, data, garchenv, timer, message, fname)
{
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
mspec = get("mspec", garchenv)
npar = get("npar", garchenv)
parallel = get("parallel", garchenv)
parallel.control = get("parallel.control", garchenv)
eval.se = get("eval.se", garchenv)
fitlist = get("fitlist", garchenv)
paridx = get("paridx", garchenv)
paridxa = get("paridxa", garchenv)
dccnames = mspec$optimization.model$modelnames
resids = residuals(fitlist)
sigmas = sigma(fitlist)
dccOrder = mspec$dccOrder
dccP = dccOrder[1]
dccQ = dccOrder[2]
m = length(fitlist@fit)
datanames = colnames(data)
mo = mspec$optimization.model$maxOrder
garchparameters = vector(mode = "numeric")
garchNames = vector(mode = "character")
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(datanames[i], cfnames, sep = ".") )
}
np = length(garchpars) + length(pars)
ns = np - length(pars) + 1
sol = f(pars, garchpars, data = data, returnType = "all", garchenv = garchenv)
likelihoods = sol$lik
loglikelihood = sol$llh
Rtout = sol$Rt
Qtout = sol$Qt
N = dim(data)[1]
Ht = array( 0, dim = c(m, m, N) )
stdresid = matrix(0, nrow = N, ncol = m)
if( parallel ){
if( parallel.control$pkg == "multicore" ){
tmp = mclapply(1:N, FUN = function(i){
tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
zz = eigen( tmph )
sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
return( list( H = tmph, Z = tmpz ) )
})
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("sigmas", "Rtout", "resids", local = TRUE)
tmp = sfLapply(as.list(1:N), fun = function(i){
tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
zz = eigen( tmph )
sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
return( list( H = tmph, Z = tmpz ) )
})
if(!eval.se) sfStop()
# else we'll stop the cluster when the standard errors are evaluated.
}
for(i in 1:N){
Ht[,,i] = tmp[[i]]$H
stdresid[i,] = tmp[[i]]$Z
}
} else{
tmp = lapply(as.list(1:N), FUN = function(i){
tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
zz = eigen( tmph )
sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
return( list( H = tmph, Z = tmpz ) )
})
for(i in 1:N){
Ht[,,i] = tmp[[i]]$H
stdresid[i,] = tmp[[i]]$Z
}
}
assign("stdresid", stdresid, envir = garchenv)
assign("Ht", Ht, envir = garchenv)
if(eval.se){
A = zeros( np, np )
tidx = 1
for(i in 1:m){
cvar = fitlist@fit[[i]]@fit$cvar
workingsize = dim(cvar)[1]
A[(tidx:(tidx + workingsize - 1)), (tidx:(tidx + workingsize - 1))] = solve(cvar)
tidx = tidx + workingsize
}
cat("\n\nCalculating Standard Errors, this can take a while\n")
otherA = .dcchessian(f = f, pars = pars, garchpars = garchpars, N = length(pars), data = data,
garchenv = garchenv, fname)
A[(np - length(pars) + 1):np, ] = otherA
jointscores = zeros(N, np)
tidx = 1
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
workingsize = length(cf)
# head(fitlist@fit[[i]]@fit$scores, 22)
jointscores[,(tidx:(tidx + workingsize - 1))] = fitlist@fit[[i]]@fit$scores
tidx = tidx + workingsize
}
#Now all we need to do is calculate the scores form the dcc estimator and we have everything
h = pmax( abs( pars/2 ), 1e-2 ) * .eps^(1/3)
hplus = pars + h
hminus = pars - h
likelihoodsplus = zeros( N, length(pars) )
likelihoodsminus = zeros( N, length(pars) )
for(i in 1:length(pars)){
hparameters1 = pars
hparameters2 = pars
hparameters1[i] = hplus[i]
hparameters2[i] = hminus[i]
LHT1 = f(pars = hparameters1, garchpars, data = data, returnType = "lik", garchenv = garchenv)
LHT2 = f(pars = hparameters2, garchpars, data = data, returnType = "lik", garchenv = garchenv)
likelihoodsplus[, i] = LHT1
likelihoodsminus[, i] = LHT2
}
sctemp = likelihoodsplus - likelihoodsminus
DCCscores = matrix(NA, ncol = dim(sctemp)[2], nrow = dim(sctemp)[1])
sdtemp = 2 * repmat( t( h ), N, 1 )
for(i in 1:dim(sctemp)[2]){
DCCscores[,i] = sctemp[,i] / sdtemp[,i]
}
jointscores[, ns:np] = DCCscores
B = cov( jointscores )
A = A/ (N)
dcccvar = ( solve( A ) %*% B %*% solve( A ) ) / N
allpars = c(garchpars, pars)
allnames = c(garchNames, dccnames)
se.coef = sqrt(diag(abs(dcccvar)))
tval = as.numeric( allpars/se.coef )
pval = 2* ( 1 - pnorm( abs( tval ) ) )
matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
matcoef[, 1] = allpars
matcoef[, 2] = se.coef
matcoef[, 3] = tval
matcoef[, 4] = pval
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
} else{
allpars = c(garchpars, pars)
allnames = c(garchNames, dccnames)
se.coef = rep(NA, length(allpars))
tval = rep(NA, length(allpars))
pval = rep(NA, length(allpars))
matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
matcoef[, 1] = allpars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
dcccvar = NULL
jointscores = NULL
}
dccfit = list()
dccfit$coef = allpars
names(dccfit$coef) = allnames
dccfit$dcccoef = pars
names(dccfit$dcccoef) = dccnames
dccfit$garchcoef = garchpars
names(dccfit$garchcoef) = garchNames
dccfit$matcoef = matcoef
dccfit$cvar = dcccvar
dccfit$scores = jointscores
dccfit$R = Rtout
dccfit$H = Ht
dccfit$Q = Qtout
dccfit$llh = loglikelihood
dccfit$log.likelihoods = likelihoods
dccfit$garchmodel = garchmodel
dccfit$timer = timer
dccfit$convergence = 0
dccfit$message = message
return( dccfit )
}
.dccmakefiltermodel = function(garchmodel, f, pars, garchpars, data, garchenv, timer, message, fname)
{
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
mspec = get("mspec", garchenv)
npar = get("npar", garchenv)
parallel = get("parallel", garchenv)
parallel.control = get("parallel.control", garchenv)
eval.se = get("eval.se", garchenv)
filterlist = get("filterlist", garchenv)
paridx = get("paridx", garchenv)
paridxa = get("paridxa", garchenv)
dccnames = mspec$optimization.model$modelnames
resids = residuals(filterlist)
sigmas = sigma(filterlist)
dccOrder = mspec$dccOrder
dccP = dccOrder[1]
dccQ = dccOrder[2]
m = length(filterlist@filter)
datanames = colnames(data)
mo = mspec$optimization.model$maxOrder
garchparameters = vector(mode = "numeric")
garchNames = vector(mode = "character")
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(datanames[i], cfnames, sep = ".") )
}
np = length(garchpars) + dccP + dccQ
ns = np - dccP - dccQ + 1
sol = f(pars, garchpars, data = data, returnType = "all", garchenv = garchenv)
likelihoods = sol$lik
loglikelihood = sol$llh
Rtout = sol$Rt
Qtout = sol$Qt
N = dim(data)[1]
Ht = array( 0, dim = c(m, m, N) )
stdresid = matrix(0, nrow = N, ncol = m)
if( parallel ){
if( parallel.control$pkg == "multicore" ){
tmp = mclapply(1:N, FUN = function(i){
tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
zz = eigen( tmph )
sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
return( list( H = tmph, Z = tmpz ) )
})
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("sigmas", "Rtout", "resids", local = TRUE)
tmp = sfLapply(as.list(1:N), fun = function(i){
tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
zz = eigen( tmph )
sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
return( list( H = tmph, Z = tmpz ) )
})
sfStop()
}
for(i in 1:N){
Ht[,,i] = tmp[[i]]$H
stdresid[i,] = tmp[[i]]$Z
}
} else{
tmp = lapply(as.list(1:N), FUN = function(i){
tmph = diag( sigmas[i, ] ) %*% Rtout[[i]] %*% diag( sigmas[i, ] )
zz = eigen( tmph )
sqrtzz = ( zz$vectors %*% diag( sqrt( zz$values ) ) %*% solve( zz$vectors ) )
tmpz = as.numeric( resids[i, ] %*% solve( sqrtzz ) )
return( list( H = tmph, Z = tmpz ) )
})
for(i in 1:N){
Ht[,,i] = tmp[[i]]$H
stdresid[i,] = tmp[[i]]$Z
}
}
assign("stdresid", stdresid, envir = garchenv)
assign("Ht", Ht, envir = garchenv)
if( eval.se ){
A = zeros( np, np )
tidx = 1
for(i in 1:m){
cvar = filterlist@filter[[i]]@filter$cvar
workingsize = dim(cvar)[1]
A[(tidx:(tidx + workingsize - 1)), (tidx:(tidx + workingsize - 1))] = solve(cvar)
tidx = tidx + workingsize
}
cat("\n\nCalculating Standard Errors, this can take a while\n")
otherA = .dcchessian(f = f, pars = pars, garchpars = garchpars, N = dccP + dccQ, data = data,
garchenv = garchenv, fname)
A[(np - dccP - dccQ + 1):np, ] = otherA
jointscores = zeros(N, np)
tidx = 1
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
workingsize = length(cf)
jointscores[,(tidx:(tidx + workingsize - 1))] = filterlist@filter[[i]]@filter$scores
tidx = tidx + workingsize
}
#Now all we need to do is calculate the scores form the dcc estimator and we have everything
h = pmax( abs( pars/2 ), 1e-2 ) * .eps^(1/3)
hplus = pars + h
hminus = pars - h
likelihoodsplus = zeros( N, dccP + dccQ )
likelihoodsminus = zeros( N, dccP + dccQ )
for(i in 1:(dccP + dccQ)){
hparameters1 = pars
hparameters2 = pars
hparameters1[i] = hplus[i]
hparameters2[i] = hminus[i]
LHT1 = f(pars = hparameters1, garchpars, data = data, returnType = "lik", garchenv = garchenv)
LHT2 = f(pars = hparameters2, garchpars, data = data, returnType = "lik", garchenv = garchenv)
likelihoodsplus[, i] = LHT1
likelihoodsminus[, i] = LHT2
}
sctemp = likelihoodsplus - likelihoodsminus
DCCscores = matrix(NA, ncol = dim(sctemp)[2], nrow = dim(sctemp)[1])
sdtemp = 2 * repmat( t( h ), N, 1 )
for(i in 1:dim(sctemp)[2]){
DCCscores[,i] = sctemp[,i] / sdtemp[,i]
}
jointscores[, ns:np] = DCCscores
B = cov( jointscores )
A = A/ N
dcccvar = ( solve( A ) %*% B %*% solve( A ) ) / N
allpars = c(garchpars, pars)
allnames = c(garchNames, dccnames)
se.coef = sqrt(diag(abs(dcccvar)))
tval = as.numeric( allpars/se.coef )
pval = 2* ( 1 - pnorm( abs( tval ) ) )
matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
matcoef[, 1] = allpars
matcoef[, 2] = se.coef
matcoef[, 3] = tval
matcoef[, 4] = pval
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
} else{
allpars = c(garchpars, pars)
allnames = c(garchNames, dccnames)
se.coef = rep(NA, length(allpars))
tval = rep(NA, length(allpars))
pval = rep(NA, length(allpars))
matcoef = matrix(NA, nrow = length(allpars), ncol = 4)
matcoef[, 1] = allpars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
dcccvar = NULL
jointscores = NULL
}
dccfilter = list()
dccfilter$coef = allpars
names(dccfilter$coef) = allnames
dccfilter$dcccoef = pars
names(dccfilter$dcccoef) = dccnames
dccfilter$garchcoef = garchpars
names(dccfilter$garchcoef) = garchNames
dccfilter$matcoef = matcoef
dccfilter$cvar = dcccvar
dccfilter$scores = jointscores
dccfilter$R = Rtout
dccfilter$H = Ht
dccfilter$Q = Qtout
dccfilter$llh = loglikelihood
dccfilter$log.likelihoods = likelihoods
dccfilter$garchmodel = garchmodel
dccfilter$timer = timer
dccfilter$convergence = 0
dccfilter$message = message
return( dccfilter )
}
.sqrtsymmat = function( X )
{
tmp = eigen( X )
sqrttmp = ( tmp$vectors %*% diag( sqrt( tmp$values ) ) %*% solve( tmp$vectors ) )
return( sqrttmp )
}
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.