Nothing
#################################################################################
##
## R package rugarch by Alexios Galanos Copyright (C) 2008-2022.
## This file is part of the R package rugarch.
##
## The R package rugarch is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## The R package rugarch is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
#################################################################################
# BDS Test
# Dechert, Scheinkman and LeBaron (1996)
# BDS i.i.d test of
# Brock and Potter (1993) and de Lima (1996) show that is
# applying the BDS test to the log of the squared standardized residuals
# log(z^2) the bias in the BDS test is almost corrected because the logarithmic
# transformation transforms the GARCH model into a linear additive model which
# satisfies the nuisance free parameter conditions in de Lima (1996) for the BDS
# test (does not apply to more richly parametrized models e.g. with leverage since
# they cannot be cast into linear additive models).
#.bds.test = function(z){
# fNonlinear::bdsTest(log(z[is.finite(z)]^2), m = 5)
#}
.information.test = function(LLH, nObs, nPars)
{
AIC = (-2*LLH)/nObs + 2 * nPars/nObs
BIC = (-2*LLH)/nObs + nPars * log(nObs)/nObs
SIC = (-2*LLH)/nObs + log((nObs+2*nPars)/nObs)
HQIC = (-2*LLH)/nObs + (2*nPars*log(log(nObs)))/nObs
informationTests = list(AIC = AIC, BIC = BIC, SIC = SIC, HQIC = HQIC)
return(informationTests)
}
# Q-Statistics on Standardized Residuals
.box.test = function(stdresid, p=1, df = 0)
{
if(any(!is.finite(stdresid))) stdresid[!is.finite(stdresid)]=0
# p=1 normal case, p=2 squared std. residuals
# Q-Statistics on Standardized Residuals
#H0 : No serial correlation ==> Accept H0 when prob. is High [Q < Chisq(lag)]
box10 = Box.test(stdresid^p, lag = 1, type = "Ljung-Box", fitdf = 0)
box15 = Box.test(stdresid^p, lag = df+1, type = "Ljung-Box", fitdf = df)
box20 = Box.test(stdresid^p, lag = df+5, type = "Ljung-Box", fitdf = df)
LBSR<-matrix(NA,ncol=2,nrow=3)
LBSR[1:3,1] = c(box10$statistic[[1]],box15$statistic[[1]],box20$statistic[[1]])
LBSR[1:3,2] = c(box10$p.value[[1]],box15$p.value[[1]],box20$p.value[[1]])
rownames(LBSR) = c(paste("Lag[1]",sep=""), paste("Lag[p+q+1][",df+1,"]",sep=""), paste("Lag[p+q+5][",df+5,"]",sep=""))
colnames(LBSR) = c("statistic","p-value")
return(LBSR)
}
.weightedBoxTest = function(stdresid, p=1, df = 0)
{
if(any(!is.finite(stdresid))) stdresid[!is.finite(stdresid)]=0
# p=1 normal case, p=2 squared std. residuals
# Q-Statistics on Standardized Residuals
#H0 : No serial correlation ==> Accept H0 when prob. is High [Q < Chisq(lag)]
box10 = Weighted.Box.test(stdresid, lag = 1, type = "Ljung-Box", fitdf = 0, if(p==2) sqrd.res = TRUE else sqrd.res = FALSE)
box15 = Weighted.Box.test(stdresid, lag = max(2, 2*df+df-1), type = "Ljung-Box", fitdf = df, if(p==2) sqrd.res = TRUE else sqrd.res = FALSE)
box20 = Weighted.Box.test(stdresid, lag = max(5, 4*df+df-1), type = "Ljung-Box", fitdf = df, if(p==2) sqrd.res = TRUE else sqrd.res = FALSE)
LBSR<-matrix(NA,ncol=2,nrow=3)
LBSR[1:3,1] = c(box10$statistic[[1]],box15$statistic[[1]],box20$statistic[[1]])
LBSR[1:3,2] = c(box10$p.value[[1]],box15$p.value[[1]],box20$p.value[[1]])
rownames(LBSR) = c(paste("Lag[1]",sep=""), paste("Lag[2*(p+q)+(p+q)-1][",max(2, 2*df+df-1),"]",sep=""), paste("Lag[4*(p+q)+(p+q)-1][",max(5, 4*df+df-1),"]",sep=""))
colnames(LBSR) = c("statistic","p-value")
return(LBSR)
}
.archlmtest = function (x, lags, demean = FALSE)
{
if(any(!is.finite(x))) x[!is.finite(x)] = 0
x = as.vector(x)
if(demean) x = scale(x, center = TRUE, scale = FALSE)
lags = lags + 1
mat = embed(x^2, lags)
arch.lm = summary(lm(mat[, 1] ~ mat[, -1]))
STATISTIC = arch.lm$r.squared * length(resid(arch.lm))
names(STATISTIC) = "Chi-squared"
PARAMETER = lags - 1
names(PARAMETER) = "df"
PVAL = 1 - pchisq(STATISTIC, df = PARAMETER)
METHOD = "ARCH LM-test"
result = list(statistic = STATISTIC, parameter = PARAMETER,
p.value = PVAL, method = METHOD)
class(result) = "htest"
return(result)
}
.weightedarchlmtest = function (x, sigma, lags, fitdf = 2, demean = FALSE)
{
if(any(!is.finite(x))) x[!is.finite(x)] = 0
x = as.vector(x)
if(demean) x = scale(x, center = TRUE, scale = FALSE)
result = Weighted.LM.test(x, sigma^2, lag = lags, type = c("correlation", "partial")[1], fitdf = fitdf, weighted=TRUE)
return(result)
}
.nyblomTest = function(object)
{
#pnames = rownames(object@fit$ipars[object@fit$ipars[,"Estimate"]==1,])
grad = object@fit$scores
pnames = colnames(grad)
if(is.null(pnames)) pnames = rownames(object@fit$ipars[object@fit$ipars[,"Estimate"]==1,,drop=FALSE])
# special case when fixed parameters exist but the fixed.se was used in fit.control options
#if( length(pnames) != dim(grad)[2] && dim(grad)[2] == length(object@fit$ipars[object@fit$ipars[,"Include"]==1,]) ){
# pnames = rownames(object@fit$ipars[object@fit$ipars[,"Include"]==1,])
#}
if(is(object, "ARFIMAfit")) res = object@fit$residuals/object@model$pars["sigma", 1] else res = object@fit$residuals/object@fit$sigma
res[!is.finite(res) | is.nan(res) | is.na(res)] = 0
nn = length(res)
hes = t(grad)%*%(grad)
shes = try(solve(hes), silent = TRUE)
if(inherits(shes,"try-error")){
shes = try( solve(object@fit$hessian) )
if(inherits(shes,"try-error")){
IndividualStat = matrix(rep(NA, length(pnames)), ncol=1)
IndividualCritical = .nyblomCritical(1)
JointCritical = .nyblomCritical(length(pnames))
rownames(IndividualStat) = pnames
JointStat = NA
} else{
zx = matrix(cbind(res, res^2, res^3, res^4), ncol = 4)
zs = apply(zx, 2, FUN = function(x) scale(x, center = TRUE, scale = FALSE))
x = as.matrix(apply(grad, 2, FUN = function(x) cumsum(x)))
xx = t(x)%*%x
nyblomj = sum(diag(xx%*%shes))/nn
nyblomt = diag(xx)/(diag(hes)*nn)
IndividualStat = matrix(nyblomt, ncol=1)
IndividualCritical = .nyblomCritical(1)
JointCritical = .nyblomCritical(length(pnames))
rownames(IndividualStat) = pnames
JointStat = nyblomj
}
} else{
zx = matrix(cbind(res, res^2, res^3, res^4), ncol = 4)
zs = apply(zx, 2, FUN = function(x) scale(x, center = TRUE, scale = FALSE))
x = as.matrix(apply(grad, 2, FUN = function(x) cumsum(x)))
xx = t(x)%*%x
nyblomj = sum(diag(xx%*%shes))/nn
nyblomt = diag(xx)/(diag(hes)*nn)
IndividualStat = matrix(nyblomt, ncol=1)
IndividualCritical = .nyblomCritical(1)
JointCritical = .nyblomCritical(length(pnames))
rownames(IndividualStat) = pnames
JointStat = nyblomj
}
return(list(IndividualStat=IndividualStat,JointStat=JointStat,
IndividualCritical=IndividualCritical,JointCritical=JointCritical))
}
.nyblomCritical = function(n){
# Test for Constancy of Parameters: Hansen Tests (1990 mimeo paper)
# Null Hypothesis: constant parameters
cval = matrix(c(0.353, 0.470, 0.748,
0.610, 0.749, 1.07,
0.846, 1.01, 1.35,
1.07, 1.24, 1.60,
1.28, 1.47, 1.88,
1.49, 1.68, 2.12,
1.69, 1.90, 2.35,
1.89, 2.11, 2.59,
2.10, 2.32, 2.82,
2.29, 2.54, 3.05,
2.49, 2.75, 3.27,
2.69, 2.96, 3.51,
2.89, 3.15, 3.69,
3.08, 3.34, 3.90,
3.26, 3.54, 4.07,
3.46, 3.75, 4.30,
3.64, 3.95, 4.51,
3.83, 4.14, 4.73,
4.03, 4.33, 4.92,
4.22, 4.52, 5.13), byrow=T, ncol=3)
colnames(cval)=c("10%","5%","1%")
if(n<=20){
ans = cval[n,]
names(ans)=c("10%","5%","1%")
} else {
ans="too many parameters"
}
return(ans)
}
.signbiasTest = function(object)
{
if(is(object, "uGARCHfilter")) z = object@filter$z else z = z = object@fit$z
res = as.numeric(residuals(object))
z2 = z^2
n = length(z)
zminus = as.integer(res<0)
zplus = 1-zminus
czminus = zminus*res
czplus = zplus*res
cz = cbind(rep(1, n), zminus, czminus, czplus)
cz = cz[1:(n-1),]
z2 = matrix(z2[2:n],ncol = 1)
cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
fitA = lm(y~const+zminus+czminus+czplus-1, data = cm)
resA = residuals(fitA)
sgmatrix = matrix(ncol = 3, nrow = 4)
rownames(sgmatrix) = c("Sign Bias","Negative Sign Bias","Positive Sign Bias","Joint Effect")
colnames(sgmatrix) = c("t-value","prob","sig")
sgmatrix[1:3,1:2] = abs(summary(fitA)$coef[2:4,3:4])
jeffect = .linear.hypothesis(fitA, c("zminus = 0", "czminus = 0","czplus = 0"), test = "Chisq")
sgmatrix[4,1] = jeffect[5]$Chisq[2]
sgmatrix[4,2] = jeffect[6][2,]
sgmatrix = as.data.frame(sgmatrix)
sgmatrix[,3] = .stars(sgmatrix[,2])
return(sgmatrix)
}
.gofTest = function(object, groups)
{
modelinc = object@model$modelinc
idx = object@model$pidx
if(is(object, "uGARCHfilter")){
ipars = object@filter$ipars
z = object@filter$z
} else{
ipars = object@fit$ipars
z = object@fit$z
}
# must remove fixed parameters
dist = object@model$modeldesc$distribution
cdfv = pdist(dist, q = sort(z), lambda = ipars[idx["ghlambda",1],1], skew = ipars[idx["skew",1],1],
shape = ipars[idx["shape",1],1])
j = length(groups)
gofmat = matrix(NA, ncol = 3, nrow = j)
gofmat[,1] = groups
for(i in 1:j){
sq = seq(1/groups[i], 1, by = 1/groups[i])
ni = tabulate(findInterval(cdfv, c(0, sq), rightmost.closed = TRUE,
all.inside = FALSE))
ExpValue = length(cdfv)/groups[i]
gofmat[i,2] = sum( ( ( ni - ExpValue ) ^2 )/ExpValue )
gofmat[i,3] = pchisq(q=gofmat[i,2], df = groups[i]-1, lower.tail = FALSE)
}
colnames(gofmat) = c("group", "statistic", "p-value(g-1)")
rownames(gofmat) = 1:j
return(gofmat)
}
# returns various loss functions (vectors)
lossfn.var = function(realized, forecast)
{
x = cbind(realized, forecast)
#exc = which(is.na(x))
#if(length(exc)>0) x = x[-exc,]
mse = (x[,1] - x[,2])^2
# qlike: quasi likelihood (implied by gaussian likelihood)
qlike = log(x[,2]) + x[,1]/x[,2]
# r2log: penalizes forecasts asymmetrically in low and high volatility periods
r2log = ( log(x[,1]/x[,2]) )^2
mad = abs(x[,1] - x[,2])
#hmse = (x[,1]*(1/x[,2]))^2
lossdf = data.frame(cbind(mse, mad, qlike, r2log))
return(lossdf)
}
lossfn.ret = function(realized, forecast)
{
x = cbind(realized, forecast)
#exc = which(is.na(x), arr.ind = TRUE)
#if(length(exc)>0) x = x[-exc[,1],]
mse = (x[,1] - x[,2])^2
mae = abs(x[,1] - x[,2])
#hmse = (x[,1]*(1/x[,2]))^2
# we count zero as positive for investment purposes
dac = as.integer(.signpluszero(x[,1])==.signpluszero(x[,2]))
lossdf = data.frame(cbind(mse, mae, dac))
return(lossdf)
}
.meanlossfn = function(lossdf)
{
n = dim(lossdf)[1]
nm = colnames(lossdf)
ans = apply(lossdf, 2, FUN = function(x) 1/(n-1) * sum(x[!is.na(x)]))
names(ans) = nm
return(ans)
}
.medianlossfn = function(lossdf)
{
n = dim(lossdf)[1]
xn = dim(lossdf)[2]
nm = colnames(lossdf)
ans = apply(lossdf[,-xn], 2, FUN = function(x) median(x, na.rm = TRUE))
# obviously we cannot have the median of a 0/1 series.
ans = c(ans, .meanlossfn(lossdf[,xn]))
names(ans) = nm
return(ans)
}
# Diebold Mariano Test (JBES 1995)
#.DM.test = function(x, lagq)
#{
# x = vector of differences btw the loss function of the benchmark model and the loss fn of the competing model
# lagq = lag considered to calculate the newey-west var-cov matrix
#n.ahead = forecast horizon
# Returns:
# an anova class object from the linear.hypothesis test of library car
# x = rnorm(100, 0, 0.01) - rnorm(100,0.0001,0.01)
# exc = which(is.na(x))
# if(length(exc)>0) x = x[-exc]
# dm = lm(x~1)
# dm.test = .linear.hypothesis(dm, c("(Intercept) = 0"), test = "Chisq", vcov = NeweyWest(dm, lag = lagq,
# order.by = NULL, prewhite = TRUE, adjust = TRUE))
# return(dm.test)
#}
# From his 2001 paper. Also added the Jarque Bera Test as reccomended by Dowd
# since the test does not really account for residuals being from the Normal distribution.
BerkowitzTest = function(data, lags = 1, significance = 0.05, tail.test = FALSE, alpha = 0.05)
{
if(tail.test){
ans = .BerkowitztLRtail(data, alpha = alpha, significance = significance)
ans$rho = NA
ans$JB = NA
ans$JBp = NA
return(ans)
} else{
if( lags < 1 ) stop("\nlags must be 1 or greater!") else lags = as.integer(lags)
x = as.numeric(data) - mean(data)
n = length(data)
xlag = NULL
for(i in 1:lags) xlag = cbind(xlag, .lagx(x, n.lag = i, pad = 0))
ans = lm(x~xlag-1)
uLL = sum(dnorm(residuals(ans)[-c(1:lags)], sd = summary(ans)$sigma, log = TRUE))
rLL = sum(dnorm(data[-c(1:lags)], log = TRUE))
LR = 2*(uLL - rLL)
chid = 1-pchisq(LR, 2 + lags)
if(chid < significance) res = paste("reject NULL") else res = paste("fail to reject NULL")
H0 = paste("Normal(0,1) with no autocorrelation")
m1 = sum(x)/n
xm = (x - m1)
m2 = sum(xm^2)/n
m3 = sum(xm^3)/n
m4 = sum(xm^4)/n
k1 = (m3/m2^(3/2))^2
k2 = (m4/m2^2)
JB = n * k1/6 + n * (k2 - 3)^2/24
JBp = 1 - pchisq(JB, df = 2)
return(list(uLL = uLL, rLL = rLL, LR = LR, LRp = chid, H0 = H0, Decision = res,
mu = mean(data), sigma = summary(ans)$sigma, rho = coef(ans)[1:(lags)],
JB = JB, JBp = JBp))
}
}
.BerkowitztLRtail = function(data, alpha = 0.05, significance = 0.05){
.lrh = function(pars, x){
p1 = x[which(x<qnorm(alpha))]
p2 = x[which(x>=qnorm(alpha))]*0 + qnorm(alpha)
-( sum(log(dnorm((p1-pars[1])/pars[2])/pars[2])) + sum(log(1-pnorm((p2 - pars[1])/pars[2]) )) )
}
tst = solnp(c(0.0, 0.95), fun = .lrh, LB = c(-10, 0.01), UB = c(10, 3), x = data, control=list(trace = FALSE))
uLL = -tail(tst$values, 1)
rLL = -.lrh(c(0, 1), data)
LR = 2 * (uLL - rLL)
chid = 1 - pchisq(LR, 2)
if (chid < significance)
res = paste("reject NULL")
else res = paste("fail to reject NULL")
H0 = paste("Normal(0,1)")
return(list(uLL = uLL, rLL = rLL, LR = LR, LRp = chid, H0 = H0,
Decision = res, mu = tst$par[1], sigma = tst$par[2]))
}
# Tests of Directional Accuracy
DACTest = function(forecast, actual, test = c("PT", "AG"), conf.level = 0.95)
{
n = length(actual)
if( length(forecast) != n ) stop("Length of forecast and actual must be the same")
if( test == "PT"){
x_t = z_t = y_t = rep(0, n)
x_t[which(actual>0)] = 1
y_t[which(forecast>0)] = 1
p_y = mean(y_t)
p_x = mean(x_t)
z_t[which( (forecast*actual)>0 )] = 1
p_hat = mean(z_t)
p_star = p_y*p_x + (1 - p_y)*(1-p_x)
p_hat_var = (p_star*(1-p_star))/n
p_star_var = ((2*p_y-1)^2*(p_x*(1-p_x)))/n + ((2*p_x-1)^2*(p_y*(1-p_y)))/n + (4*p_x*p_y*(1-p_x)*(1-p_y))/n^2
s_n = (p_hat - p_star)/sqrt(p_hat_var - p_star_var)
ans = list(Test = "Pesaran and Timmermann", Stat = s_n, p.value = 1-pnorm(s_n), H0 = "Independently Distributed",
Decision = if( s_n > qnorm(conf.level) ) "Reject H0" else "Fail to Reject H0",
DirAcc = p_hat)
} else{
r_t=sign(forecast)*actual
A_t = mean(r_t)
B_t = mean(sign(forecast))*mean(actual)
p_y = 0.5 * (1 + mean(sign(forecast)))
V_EP = (4/(n^2))*p_y*(1-p_y)*sum((actual-mean(actual))^2)
EP = (A_t-B_t)/sqrt(V_EP)
ans = list(Test = "Anatolyev and Gerko", Stat = EP, p.value = 1-pnorm(EP), H0 = "No Predictability",
Decision = if( EP > qnorm(conf.level) ) "Reject H0" else "Fail to Reject H0",
DirAcc = sum(r_t>0)/n)
}
return( ans )
}
VaRTest = function(alpha = 0.05, actual, VaR, conf.level = 0.95){
N = length(actual)
VaRn = floor(N * alpha)
if(N != length(VaR)) stop("\nlength of realized not equal to length of VaR!")
tmp = LR.cc.test(p = alpha, actual = actual, VaR = VaR, conf.level = conf.level)
ans = list()
ans$expected.exceed = floor(alpha*tmp$TN)
ans$actual.exceed = tmp$N
ans$uc.H0 = "Correct Exceedances"
ans$uc.LRstat = tmp$stat.uc
ans$uc.critical = tmp$crit.val.uc
ans$uc.LRp = tmp$p.value.uc
ans$uc.Decision = ifelse(ans$uc.LRp<(1-conf.level), "Reject H0", "Fail to Reject H0")
ans$cc.H0 = "Correct Exceedances & Independent"
ans$cc.LRstat = tmp$stat.cc
ans$cc.critical = tmp$crit.val.cc
ans$cc.LRp = tmp$p.value.cc
ans$cc.Decision = ifelse(ans$cc.LRp<(1-conf.level), "Reject H0", "Fail to Reject H0")
return(ans)
}
ESTest = function(alpha = 0.05, actual, ES, VaR, conf.level = 0.95, boot = FALSE, n.boot = 1000){
N = length(actual)
if(N != length(VaR)) stop("\nlength of realized not equal to length of VaR!")
idx = which(as.numeric(VaR)>actual)
z = (ES[idx] - actual[idx])
.fn = function(x){
n = length(x)
1-pnorm(mean(x)/((sqrt(sd(x)^2*((n-1)/n)))/sqrt(n-1)))
}
if(boot){
rbt = matrix(sample(z, size = length(z) * n.boot, replace = TRUE), nrow = n.boot)
pv = mean( apply(rbt, 1, FUN = function(x) .fn(x)) )
npv = .fn(z)
} else{
npv = .fn(z)
pv = NA
}
ans = list()
ans$expected.exceed = floor(alpha*N)
ans$actual.exceed = length(idx)
# conditional expected shortfall is systematically underestimated
ans$H1 = "Mean of Excess Violations of VaR is greater than zero"
ans$boot.p.value = pv
ans$p.value = npv
ans$Decision = ifelse(npv<(1-conf.level),"Reject H0", "Fail to Reject H0")
return(ans)
}
VaRDurTest = function(alpha, actual, VaR, conf.level = 0.95){
VaR.ind = ifelse(actual < VaR, 1, 0)
N = sum(VaR.ind)
TN = length(VaR.ind)
D = diff(which(VaR.ind==1))
C = rep(0, length(D))
# left-censored
if(VaR.ind[1]==0){
C = c(1, C)
# the number of days until we get the first hit
D = c(which(VaR.ind==1)[1], D)
}
# right-censored
if(VaR.ind[TN]==0){
C = c(C, 1)
# the number of days after the last one in the hit sequence
D = c(D, TN - tail(which(VaR.ind==1), 1))
}
N = length(D)
sol = try(optim(par = 2, fn = .likDurationW, gr = NULL, D = D, C = C, N = N,
method = "L-BFGS-B", lower = 0.001, upper = 10, control = list(trace=0)), silent = TRUE)
b = sol$par
uLL = -sol$value
rLL = -.likDurationW(1, D, C, N)
LR = 2*(uLL - rLL)
LRp = 1 - pchisq(LR, 1)
H0 = "Duration Between Exceedances have no memory (Weibull b=1 = Exponential)"
#i.e. whether we fail to reject the alternative in the LR test that b=1 (hence correct model)
Decision = ifelse(LRp<(1-conf.level),"Reject H0", "Fail to Reject H0")
return(list(b = b, uLL = uLL, rLL = rLL, LRp = LRp, H0 = H0, Decision = Decision))
}
.likDurationW = function(pars, D, C, N){
b = pars[1]
a = ( (N - C[1] - C[N])/(sum(D^b)) )^(1/b)
lik = C[1]*log(.pweibull(D[1],a,b,survival=TRUE)) + (1-C[1])*.dweibull(D[1], a, b, log = TRUE) +
sum(.dweibull(D[2:(N-1)], a, b, log = TRUE) ) + C[N]*log(.pweibull(D[N],a,b,survival=TRUE) ) +
(1 - C[N]) *.dweibull(D[N], a, b, log = TRUE)
if(!is.finite(lik) || is.nan(lik)) lik = 1e10 else lik = -lik
return(lik)
}
# When b=1 we get the exponential
.dweibull = function(D, a, b, log = FALSE){
# density of Weibull
pdf = b * log(a) + log(b) + (b - 1) * log(D) - (a * D)^b
if(!log) pdf = exp(pdf)
return(pdf)
}
.pweibull = function(D, a, b, survival = FALSE){
# distribution of Weibull
cdf = 1 - exp(-(a*D)^b)
if(survival) cdf = 1 - cdf
return(cdf)
}
.hweibull = function(D, a, b){
# hazard of Weibull
h = (a^b)*b*(D^(b-1))
return(h)
}
#####################################################################################
# CHANGELOG (15-08-2012): Changed kurt argument in GMMTest to default of 3 (Normal)
# ...was badly set at 0.
GMMTest = function(z, lags = 1, skew=0, kurt=3, conf.level = 0.95){
if(length(skew)>1) sk = skew[-c(1:lags)] else sk = skew
if(length(kurt)>1) ku = kurt[-c(1:lags)] else ku = kurt
z = matrix(z, ncol = 1)
N = dim(z)[1] - lags
zlag = z[-c(1:lags), , drop = FALSE]
orthmat = matrix(NA, ncol = 8, nrow = 3)
colnames(orthmat) = c("E[z]", "E[z^2]-1", "E[z^3]", "E[z^4]-3", "Q2", "Q3", "Q4","J")
rownames(orthmat) = c("mean", "var", "t.value")
f1 = zlag[,1]
orthmat[1:3,1] = c(mean(f1), mean(f1^2)/N, mean(f1)/sqrt(mean(f1^2)/N))
f2 = (zlag[,1]^2)-1
orthmat[1:3,2] = c(mean(f2), mean(f2^2)/N, mean(f2)/sqrt(mean(f2^2)/N))
f3 = zlag[,1]^3-sk
orthmat[1:3,3] = c(mean(f3), mean(f3^2)/N, mean(f3)/sqrt(mean(f3^2)/N))
f4 = zlag[,1]^4-ku
orthmat[1:3,4] = c(mean(f4), mean(f4^2)/N, mean(f4)/sqrt(mean(f4^2)/N))
M = rbind(t(f1), t(f2), t(f3), t(f4))
tmp1 = .waldcomomtest(z[,1]^2-1, lags, N)
orthmat[3, 5] = tmp1$tval[lags+1]
tmp2 = .waldcomomtest(z[,1]^3-skew, lags, N)
orthmat[3, 6] = tmp2$tval[lags+1]
tmp3 = .waldcomomtest(z[,1]^4-kurt, lags, N)
orthmat[3, 7] = tmp3$tval[lags+1]
M = rbind(M, tmp1$h, tmp2$h, tmp3$h)
g = c(as.numeric(orthmat[1,1:4]), tmp1$g, tmp2$g, tmp3$g)
# all moments
S = (M %*% t(M))/N
jtval = N*t(g)%*%solve(S)%*%g
orthmat[3,8] = jtval
p = rep(conf.level, 2)
df = c(lags, lags, lags, (4+3*lags))
critical.values = qchisq(p,df)
# i.e. if tval>critical value reject the NULL
Decision = NULL
Decision[1] = ifelse(orthmat[3, 5]<critical.values[1], "Fail to Reject H0", "Reject H0")
Decision[2] = ifelse(orthmat[3, 6]<critical.values[2], "Fail to Reject H0", "Reject H0")
Decision[3] = ifelse(orthmat[3, 7]<critical.values[3], "Fail to Reject H0", "Reject H0")
Decision[4] = ifelse(orthmat[3, 8]<critical.values[4], "Fail to Reject H0", "Reject H0")
H0 = "[Q-Moment Conditions] Model is Correctly Specified"
moment.mat = orthmat[1:3,1:4]
joint.mat = rbind(orthmat[3,5:8], critical.values)
rownames(joint.mat) = c("t-value", "critical.value")
return(list(joint.mat = joint.mat, moment.mat = moment.mat, H0 = H0, Decision = Decision) )
}
.waldcomomtest = function(z, lags, N){
f0 = z %*% matrix(1, nrow = 1, ncol = lags)
fx = NULL
for(i in 1:lags) fx = cbind(fx, .lagx(z, n.lag = i, pad = 0))
fx = fx[-c(1:lags), , drop = FALSE]
f0 = f0[-c(1:lags), , drop = FALSE]
fmat = f0 * fx
g = colMeans(fmat)
varg = apply(fmat, 2, FUN = function(x) mean(x^2)/N)
tval = g/varg
h = t(fmat)
S = (h%*%t(h))/N
joint = N*t(g)%*%solve(S)%*%g
tval = c(tval, joint)
return(list(tval = tval, h = h, g = g, varg = varg, S = S))
}
# currenly only quartic kernel implemented
# Hong and Li Test
# M(1,2) test for ARCH-in-Mean
# M(2,1) tests for leverage
HLTest = function(PIT, lags = 4, kernel = "quartic", conf.level = 0.95){
p = lags
Vcon = 2*(50/49 - 300/294 + 1950/1960 - 900/1568 + 450/2304)^2
Qhatvector = matrix(0, p,1)
res = rep(0, 7)
Acon2 = integrate(funb, 0, 1)$value
T = length(PIT)
hpit = sd(PIT) * T^(-1/6)
Acon11 = (1/hpit-2)*(5/7)
Acon_1 = (Acon11 + 2*Acon2)^2 - 1
for(i in 1:p){
Mhat1 = gauss_legendre2D(f = ghat, 0,1,0,1, pit = PIT, i = i, T = T, hpit = hpit)
Qhatvector[i,1] = ( (T-i)*hpit*Mhat1 - hpit*Acon_1 )/sqrt(Vcon)
}
res[7] = sum(Qhatvector[,1])/sqrt(p)
res[1] = Momstat(1,1, p, PIT, T)
res[2] = Momstat(2,2, p, PIT, T)
res[3] = Momstat(3,3, p, PIT, T)
res[4] = Momstat(4,4, p, PIT, T)
res[5] = Momstat(1,2, p, PIT, T)
res[6] = Momstat(2,1, p, PIT, T)
names(res) = c("M(1,1)", "M(2,2)", "M(3,3)", "M(4,4)", "M(1,2)", "M(2,1)", "W")
Decision = NULL
RejectH0 = as.logical( res>rep( qnorm(conf.level), 7 ))
Decision[1] = ifelse(RejectH0[1], "Reject H0", "Fail to Reject H0")
Decision[2] = ifelse(RejectH0[2], "Reject H0", "Fail to Reject H0")
Decision[3] = ifelse(RejectH0[3], "Reject H0", "Fail to Reject H0")
Decision[4] = ifelse(RejectH0[4], "Reject H0", "Fail to Reject H0")
Decision[5] = ifelse(RejectH0[5], "Reject H0", "Fail to Reject H0")
Decision[6] = ifelse(RejectH0[6], "Reject H0", "Fail to Reject H0")
Decision[7] = ifelse(RejectH0[7], "Reject H0", "Fail to Reject H0")
names(Decision) = c("M(1,1)", "M(2,2)", "M(3,3)", "M(4,4)", "M(1,2)", "M(2,1)", "W")
H0 = NULL
H0[1] = "M(1,1) Correctly Specified"
H0[2] = "M(2,2) Correctly Specified"
H0[3] = "M(3,3) Correctly Specified"
H0[4] = "M(4,4) Correctly Specified"
H0[5] = "M(1,2) Correctly Specified"
H0[6] = "M(2,1) Correctly Specified"
H0[7] = "Model Correctly Specified"
return(list(statistic = res, H0 = H0, Decision = Decision))
}
funb = function(b){
tmp1 = (8/15 + b - (2/3)*(b^3)+(1/5)*(b^5))^(-2)
tmp2 = b*( (1-b^2)^4 )+128/315+(8/3)*(b^3)-(24/5)*(b^5)+(24/7)*(b^7)-(8/9)*(b^9)
ans = tmp1 * tmp2
return(ans)
}
ghat = function(x, pit, i, T, hpit){
z1 = x[1]
z2 = x[2]
n = length(pit)
b1 = fnbound( z1, pit[-c(1:i)], hpit)
b2 = fnbound( z2, pit[-c((n-i+1):n)], hpit)
g = sum(b1* b2)/(T-i)
return( (g - 1)^2 )
}
fnbound = function(x, y, hpit)
{
K1 = quartic( (x - y)/hpit )/hpit
if( (x >= 0)*(x < hpit) ){
K2 = integrate(quartic, (-x/hpit), 1)$value
return(K1/K2)
} else if ((x >= hpit)*(x <= (1-hpit))){
return( K1 )
} else if( (x>(1-hpit))*(x <= 1) ){
K2 = integrate(quartic, -1, (1-x)/hpit )$value
return(K1/K2)
}
}
quartic = function(z){
kernell = rep(0, length(z))
nonzeros = which(z<=1 & z>=-1)
v = z[nonzeros]
v = 1-v^2
kernell[nonzeros] = (15/16)*(v^2)
return(kernell)
}
Momstat = function(m, ll, p, pit, T){
part1 = 0
part2 = 0
part3 = 0
for(j in 1:p){
w = 1 - j/p
part1 = part1 + (w^2)*(T-j)*((Rcc1(j, m, ll, pit, T)/Rcc1(0, m, ll, pit, T))^2)
part2 = part2 + w^2
part3 = part3 + w^4
}
return( (part1-part2)/part3 )
}
Rcc1 = function(j, m, ll, pit, T){
R1 = sum( (pit[(j+1):T]^m) * (pit[1:(T-j)]^ll) ) / T
R2 = sum( (pit[(j+1):T]^m) ) / T
R3 = sum( (pit[1:(T-j)]^ll) ) / T
return( R1 - R2*R3 )
}
##############################################################################
# 2D Integration
# Code from forums:
# http://tolstoy.newcastle.edu.au/R/e5/help/08/09/2919.html
# posted by Earl F. Glynn
gauss_legendre2D_helper <- function(f, x, a2,b2, nodes, weights, ...) {
C <- (b2 - a2) / 2
D <- (b2 + a2) / 2
sum <- 0.0
for (i in 1:length(nodes))
{
y <- nodes[i]*C + D
sum <- sum + weights[i] * f(c(x,y), ...)
}
return(C * sum)
}
gauss_legendre2D <- function(f, a1,b1, a2,b2, ...) {
#library(statmod)
#N <- 12
#GL <- gauss.quad(N)
nodes = c(-0.9815606, -0.9041173, -0.7699027, -0.5873180, -0.3678315,
-0.1252334, 0.1252334, 0.3678315, 0.5873180, 0.7699027, 0.9041173,
0.9815606)
weights = c(0.04717534, 0.10693933, 0.16007833, 0.20316743,
0.23349254, 0.24914705, 0.24914705, 0.23349254, 0.20316743,
0.16007833, 0.10693933, 0.04717534)
C <- (b1 - a1) / 2
D <- (b1 + a1) / 2
sum <- 0.0
for (i in 1:length(nodes))
{
x <- nodes[i]*C + D
sum <- sum + weights[i] * gauss_legendre2D_helper(f, x, a2, b2, nodes, weights, ...)
}
return(C * sum)
}
#####################################################################################
.meanlossfn = function(lossdf)
{
n = dim(lossdf)[1]
nm = colnames(lossdf)
ans = apply(lossdf, 2, FUN = function(x) 1/(n-1) * sum(x[!is.na(x)]))
names(ans) = nm
return(ans)
}
.medianlossfn = function(lossdf)
{
n = dim(lossdf)[1]
xn = dim(lossdf)[2]
nm = colnames(lossdf)
ans = apply(lossdf[,-xn], 2, FUN = function(x) median(x, na.rm = TRUE))
# obviously we cannot have the median of a 0/1 series.
ans = c(ans, .meanlossfn(lossdf[,xn]))
names(ans) = nm
return(ans)
}
.signpluszero = function(x)
{
z = as.numeric(x)
ans = rep(0, length(z))
ans[z>=0] = 1
ans
}
.VaRplot = function(varname , p, actual, dates, VaR)
{
y.actual = actual
y.VaR = VaR
x.dates = dates
zd = .makedate(as.character(dates))
if(zd$status == 0) dates = 1:length(dates) else dates = zd$dates
title = paste("Daily Returns and Value-at-Risk Exceedances\n","(Series: ", varname,", alpha=", p,")",sep="")
if(zd$status){
plot(x = as.Date(dates, format = zd$dformat), y = y.actual, type = "n",
main = title, ylab = "Daily Log Return", xlab = "time",
ylim = c(min(y.actual, y.VaR), max(y.actual, y.VaR)))
} else{
plot(as.numeric(dates), y.actual, type = "n",
main = title, ylab = "Daily Log Return", xlab = "time",
ylim = c(min(y.actual, y.VaR), max(y.actual, y.VaR)))
}
abline(h = 0, col = 2, lty = 2)
sel = which(y.actual>0)
points(dates[sel], y.actual[sel], pch = 18, col = "green")
sel = which(y.actual<0)
points(dates[sel], y.actual[sel], pch = 18, col = "orange")
sel = which(y.actual<y.VaR)
points(dates[sel], y.actual[sel], pch = 18, col = "red")
if(zd$status){
lines(x = as.Date(dates, format = "%Y-%m-%d"), y = y.VaR, lwd = 2, col = "black")
} else{
lines(as.numeric(dates), y.VaR, lwd = 2, col = "black")
}
legend("topleft", max(actual),c("return >= 0","VaR <= return < 0","return < VaR","VaR"),
col=c("green","orange","red","black"), cex=0.75,
pch = c(18,18,18,-1), lty=c(0,0,0,1), lwd=c(0,0,0,2), bty = "n")
}
.VaRreport = function(varname, garchmodel, distribution, p, actual, VaR, conf.level = 0.95)
{
actual<-as.matrix(actual)
VaR<-as.matrix(VaR)
result <- LR.cc.test(p=p,actual=actual,VaR=VaR,conf.level=conf.level)
cat("VaR Backtest Report\n")
cat("===========================================\n")
cat(paste("Model:\t\t\t\t","",garchmodel,"-",distribution,"\n",sep=""))
cat(paste("Backtest Length:\t",result$TN,"\n",sep=""))
cat(paste("Data:\t\t\t\t","",varname,"\n",sep=""))
cat("\n==========================================\n")
cat(paste("alpha:\t\t\t\t",round(100*p,1),"%\n",sep=""))
cat(paste("Expected Exceed:\t",round(p*result$TN,1),"\n",sep=""))
cat(paste("Actual VaR Exceed:\t",result$N,"\n",sep=""))
cat(paste("Actual %:\t\t\t",round(100*result$N/result$TN,1),"%\n",sep=""))
cat("\nUnconditional Coverage (Kupiec)")
cat("\nNull-Hypothesis:\tCorrect Exceedances\n")
cat(paste("LR.uc Statistic:\t",round(result$stat.uc,3),"\n",sep=""))
cat(paste("LR.uc Critical:\t\t",round(result$crit.val.uc,3),"\n",sep=""))
cat(paste("LR.uc p-value:\t\t",round(result$p.value.uc,3),"\n",sep=""))
cat(paste("Reject Null:\t\t",ifelse(round(result$p.value.uc,3)< (1-conf.level), "YES","NO"),"\n",sep=""))
cat("\nConditional Coverage (Christoffersen)")
cat("\nNull-Hypothesis:\tCorrect Exceedances and\n\t\t\t\t\tIndependence of Failures\n")
cat(paste("LR.cc Statistic:\t",round(result$stat.cc,3),"\n",sep=""))
cat(paste("LR.cc Critical:\t\t",round(result$crit.val.cc,3),"\n",sep=""))
cat(paste("LR.cc p-value:\t\t",round(result$p.value.cc,3),"\n",sep=""))
cat(paste("Reject Null:\t\t",ifelse(result$reject,"YES","NO"),"\n",sep=""))
}
########################################################################
# Available in a number of locations/textbooks.
# This code originally presented in a webcast by Guy Yollin Feb 2006 (from Insightful)
# Functions to perform Hypothesis test
# on VaR models based on # of exceedances
# calc LR.uc statistic
LR.cc.test = function (p, actual, VaR, conf.level = 0.95)
{
result = .LR.cc(p = p, actual = actual, VaR = VaR)
crit.val.uc = qchisq(conf.level, df = 1)
crit.val.cc = qchisq(conf.level, df = 2)
p.value.cc = 1 - pchisq(result$stat.cc, df = 2)
p.value.uc = 1 - pchisq(result$stat.uc, df = 1)
reject = ifelse(p.value.cc < 1 - conf.level, TRUE, FALSE)
return(list(stat.cc = result$stat.cc, stat.uc = result$stat.uc,
p.value.cc = p.value.cc, p.value.uc = p.value.uc, conf.level = conf.level,
reject = reject, N = result$N, TN = result$TN, crit.val.uc = crit.val.uc,
crit.val.cc = crit.val.cc))
}
.LR.cc = function (p, actual, VaR)
{
VaR.ind = as.numeric(ifelse(actual < VaR, 1, 0))
N = sum(VaR.ind)
TN = length(VaR.ind)
# Thanks to Kurt Hornik for the simplications and fix
n <- table(head(VaR.ind, -1), tail(VaR.ind, -1))
N00 = n[1,1]
N11 = n[2,2]
N01 = n[1,2]
N10 = n[2,1]
# p01: Probability of having a failure at time t | no failure at t-1
p01 = N01/(N00+N01)
# p11: Probability of having a failure at time t | failure at t-1
p11 = N11/(N10+N11)
# pUc : unconditional probability of a failure
pUc = (N01+N11)/sum(n)
l1 = (1-pUc)^(N00+N10)*pUc^(N01+N11)
l2 = (1-p01)^(N00)*p01^(N01)*(1-p11)^(N10)*p11^(N11)
stat.ind = -2*log(l1/l2)
stat.uc = .LR.uc(p = p, TN = TN, N = N)
stat.cc = stat.uc + stat.ind
return(list(stat.cc = stat.cc, stat.uc = stat.uc, N = N, TN = TN))
}
.LR.uc = function (p, TN, N)
{
l1 = (1-p)^(TN-N)*p^N
l2 = (1-N/TN)^(TN-N)*(N/TN)^N
stat.uc = -2 * log(l1/l2)
return(stat.uc)
}
.Log = function(x){
ans = log(x)
#if(!is.finite(ans)) ans = sign(ans) * 1e10
ans
}
###############################################################################
# MCS Test of Hansen, Lunde and Nason
# Partially transalted from Kevin Sheppard's mcs matlab function
mcsTest = function(losses, alpha, nboot = 100, nblock = 1, boot = c("stationary", "block"))
{
n = NROW(losses)
m = NCOL(losses)
bsdat = bootstrap(1:n, nboot, nblock, type = boot[1])$B
dijbar = matrix(0, m, m)
for(j in 1:m) dijbar[j, ] = colMeans(losses - losses[,j])
dijbarstar = array(0, dim = c(m, m, nboot))
for(i in 1:nboot){
tmp = colMeans(losses[bsdat[,i],,drop = FALSE])
for(j in 1:m){
dijbarstar[j,,i] = tmp - tmp[j]
}
}
tmp = array(NA, dim = c(m, m, nboot))
for(i in 1:nboot) tmp[,,i] = dijbar
xtmp = ((dijbarstar - tmp)^2)
vardijbar = matrix(0, m, m)
for(i in 1:nboot) vardijbar = vardijbar + xtmp[,,i]
vardijbar = vardijbar/nboot
vardijbar = vardijbar + diag(m)
tmp2 = array(NA, dim = c(m, m, nboot))
for(i in 1:nboot) tmp2[,,i] = sqrt(vardijbar)
z0=(dijbarstar-tmp)/tmp2
zdata0=dijbar/sqrt(vardijbar)
excludedR = matrix(0, m,1)
pvalsR = matrix(1,m,1)
for(i in 1:(m-1)){
included = setdiff(1:m,excludedR)
mx = length(included)
z = z0[included, included,]
empdistTR = apply(abs(z), 3, function(x) max(x))
zdata = zdata0[included,included]
TR = max(.Call("colMaxRcpp", zdata, PACKAGE = "rugarch"))
pvalsR[i] = mean(empdistTR>TR)
dibar = colMeans(dijbar[included,included])*(mx/(mx-1))
dibstar = apply(dijbarstar[included,included,,drop=FALSE], 3, function(x) colMeans(x))*(mx/(mx-1))
vardi=colMeans((t(dibstar)-matrix(dibar, nrow = nboot, ncol = mx, byrow = TRUE))^2)
tx=dibar/sqrt(vardi)
temp = max(tx)
modeltoremove = which(tx == max(tx))[1]
excludedR[i]=included[modeltoremove]
}
maxpval=pvalsR[1]
for(i in 2:m){
if(pvalsR[i]<maxpval){
pvalsR[i] = maxpval
} else{
maxpval = pvalsR[i]
}
}
excludedR[NROW(excludedR)] = setdiff(1:m,excludedR)
pl = which(pvalsR >= alpha)[1]
includedR = excludedR[pl:m]
if(pl==1) excludedR = NULL else excludedR = excludedR[1:(pl-1)]
excludedSQ = matrix(0,m,1)
pvalsSQ = matrix(1,m,1)
for(i in 1:(m-1)){
included = setdiff(1:m,excludedSQ)
mx = length(included)
z = z0[included,included,]
empdistTSQ = apply(z^2, 3, function(x) sum(x))/2
zdata = zdata0[included,included]
TSQ = sum(colSums(zdata^2))/2
pvalsSQ[i] = mean(empdistTSQ>TSQ)
dibar = colMeans(dijbar[included,included])*(mx/(mx-1))
dibstar = apply(dijbarstar[included,included,,drop=FALSE], 3, function(x) colMeans(x))*(mx/(mx-1))
vardi = colMeans((t(dibstar)-matrix(dibar, nrow = nboot, ncol = mx, byrow = TRUE))^2)
tx = dibar/sqrt(vardi)
temp = max(tx)
modeltoremove = which(tx == max(tx))[1]
excludedSQ[i]=included[modeltoremove]
}
maxpval = pvalsSQ[1]
for( i in 2:m){
if(pvalsSQ[i]<maxpval){
pvalsSQ[i] = maxpval
} else{
maxpval = pvalsSQ[i]
}
}
excludedSQ[NROW(excludedSQ)] = setdiff(1:m,as.numeric(excludedSQ))
pl=which(pvalsSQ>=alpha)[1]
includedSQ = excludedSQ[pl:m]
if(pl==1) excludedSQ = NULL else excludedSQ = excludedSQ[1:(pl-1)]
return(list(includedR = includedR, pvalsR = pvalsR, excludedR = excludedR,
includedSQ = includedSQ, pvalsSQ = pvalsSQ, excludedSQ = excludedSQ))
}
bootstrap = function(data, nboot = 10, nblock = 5, type = c("stationary", "block"))
{
type = match.arg(type[1], c("stationary", "block"))
ans = switch(type,
stationary = .stationary_bootstrap(data, nboot, nblock),
block = .block_bootstrap(data, nboot, nblock))
return(ans)
}
.stationary_bootstrap = function(data, nboot = 100, nblock = 1){
n = length(data)
p = 1/nblock
idx = matrix(0, n, nboot)
idx[1, ] = ceiling(n * matrix(runif(1*nboot), 1, nboot))
sel_idx = matrix(runif(n*nboot), n, nboot)<p
cn = sum(colSums(sel_idx))
idx[sel_idx] = ceiling(matrix(runif(1*cn), 1, cn)*n)
for(i in 2:n){
idx[i,!sel_idx[i,]] = idx[i-1,!sel_idx[i,]]+1
}
data = c(data, data)
B = apply(idx, 2, FUN = function(x) data[x])
return(list(B = B, idx = idx))
}
.block_bootstrap = function(data, nboot, nblock){
n = length(data)
data = c(data, data[1:(nblock-1)])
s = ceiling(n/nblock)
S = floor(matrix(runif(s*nboot), s, nboot)*n)+1
idx = matrix(0, s*nblock, nboot)
A = repmat( 0:(nblock-1), 1, nboot)
iter=1
for(i in seq(1, n, by = nblock)){
idx[i:(i+nblock-1),] = repmat(S[iter, ,drop = FALSE], nblock, 1) + A
iter = iter + 1
}
idx = idx[1:n, ,drop = FALSE]
B = apply(idx, 2, FUN = function(x) data[x])
return(list(B = B, idx = idx))
}
repmat = function(a,n,m){
kronecker(matrix(1,n,m),a)
}
##########################################################################################
# Direct Import of Weighted Tests of FISHER and GALLAGHER (WeightedPortTest package)
Weighted.Box.test = function (x, lag = 1, type = c("Box-Pierce", "Ljung-Box", "Monti"),
fitdf = 0, sqrd.res = FALSE, log.sqrd.res = FALSE, abs.res = FALSE,
weighted = TRUE)
{
if(lag<(2*fitdf+fitdf-1)) stop("\nLag must be equal to a minimum of 2*fitdf+fitdf-1")
if(NCOL(x) > 1) stop("\nx is not a vector or univariate time series")
if(lag < 1) stop("\nLag must be positive")
if(fitdf < 0) stop("\nFitdf cannot be negative")
if((sqrd.res && log.sqrd.res) || (sqrd.res && abs.res) || (log.sqrd.res && abs.res)) stop("Only one option of: sqrd.res, log.sqrd.res or abs.res can be selected")
DNAME <- deparse(substitute(x))
type <- match.arg(type)
if (abs.res) {
x <- abs(x)
}
if (sqrd.res || log.sqrd.res) {
x <- x^2
}
if (log.sqrd.res) {
x <- log(x)
}
if (weighted) {
if (type == "Monti") {
METHOD <- "Weighted Monti test (Gamma Approximation)"
cor <- acf(x, lag.max = lag, type = "partial", plot = FALSE,
na.action = na.pass)
obs <- cor$acf[1:lag]
}
else {
cor <- acf(x, lag.max = lag, type = "correlation",
plot = FALSE, na.action = na.pass)
obs <- cor$acf[2:(lag + 1)]
}
if (type == "Ljung-Box") {
METHOD <- "Weighted Ljung-Box test (Gamma Approximation)"
}
n <- sum(!is.na(x))
weights <- (lag - 1:lag + 1)/(lag)
if (type == "Box-Pierce") {
METHOD <- "Weighted Box-Pierce test (Gamma Approximation)"
STATISTIC <- n * sum(weights * obs^2)
}
else {
STATISTIC <- n * (n + 2) * sum(weights * (1/seq.int(n - 1, n - lag) * obs^2))
}
if (sqrd.res) {
fitdf <- 0
names(STATISTIC) <- "Weighted X-squared on Squared Residuals for detecting nonlinear processes"
}
else if (log.sqrd.res) {
fitdf <- 0
names(STATISTIC) <- "Weighted X-squared on Log-Squared Residuals for detecting nonlinear processes"
}
else if (abs.res) {
fitdf <- 0
names(STATISTIC) <- "Weighted X-squared on Absolute valued Residuals for detecting nonlinear processes"
}
else {
names(STATISTIC) <- "Weighted X-squared on Residuals for fitted ARMA process"
}
shape <- (3/4) * (lag + 1)^2 * lag/(2 * lag^2 + 3 * lag + 1 - 6 * lag * fitdf)
scale <- (2/3) * (2*lag^2 + 3*lag+1 - 6 * lag * fitdf)/(lag*(lag + 1))
PARAMETER <- c(shape, scale)
names(PARAMETER) <- c("Shape", "Scale")
PVAL <- 1 - pgamma(STATISTIC, shape = shape, scale = scale)
names(PVAL) <- "Approximate p-value"
}
else {
if (type == "Monti") {
METHOD <- "Monti test"
cor <- acf(x, lag.max = lag, type = "partial", plot = FALSE,
na.action = na.pass)
obs <- cor$acf[1:lag]
}
else {
cor <- acf(x, lag.max = lag, type = "correlation",
plot = FALSE, na.action = na.pass)
obs <- cor$acf[2:(lag + 1)]
}
if (type == "Ljung-Box") {
METHOD <- "Ljung-Box test"
}
n <- sum(!is.na(x))
if (type == "Box-Pierce") {
METHOD <- "Box-Pierce test"
STATISTIC <- n * sum(obs^2)
}
else {
STATISTIC <- n * (n + 2) * sum((1/seq.int(n - 1,
n - lag) * obs^2))
}
if (sqrd.res) {
fitdf <- 0
names(STATISTIC) <- "X-squared on Squared Residuals for detecting nonlinear processes"
}
else if (log.sqrd.res) {
fitdf <- 0
names(STATISTIC) <- "X-squared on Log-Squared Residuals for detecting nonlinear processes"
}
else if (abs.res) {
fitdf <- 0
names(STATISTIC) <- "X-squared on Absolute valued Residuals for detecting nonlinear processes"
}
else {
names(STATISTIC) <- "X-squared on Residuals for fitted ARMA process"
}
mydf <- lag - fitdf
PARAMETER <- c(mydf)
names(PARAMETER) <- c("df")
PVAL <- 1 - pchisq(STATISTIC, df = mydf)
names(PVAL) <- "p-value"
}
structure(list(statistic = STATISTIC, parameter = PARAMETER,
p.value = PVAL, method = METHOD, data.name = DNAME),
class = "htest")
}
Weighted.LM.test <- function (x, h.t, lag = 1, type = c("correlation", "partial"), fitdf = 1, weighted=TRUE)
{
### Error Checking
###
if (NCOL(x) > 1)
stop("x is not a vector or univariate time series");
if (fitdf >= lag)
stop("Lag must exceed fitted degrees of freedom");
if (fitdf < 1)
stop("Fitted degrees of freedom must be positive");
if( !(length(x)==length(h.t)) )
stop("Length of x and h.t must match");
DNAME <- deparse(substitute(x))
type <- match.arg(type)
x <- x^2/h.t
if( type == "partial") {
cor <- acf(x, lag.max = lag, type="partial", plot=FALSE, na.action=na.pass)
obs <- cor$acf[1:lag];
}
else {
cor <- acf(x, lag.max = lag, type="correlation", plot=FALSE, na.action=na.pass)
obs <- cor$acf[2:(lag + 1)];
}
if(type == "correlation" && weighted) {
METHOD <- "Weighted Li-Mak test on autocorrelations (Gamma Approximation)"
}
else if(type == "partial" && weighted) {
METHOD <- "Weighted Li-Mak test on partial autocorrelations (Gamma Approximation)"
}
else if(type == "correlation" && !weighted) {
METHOD <- "Li-Mak test on autocorrelations (Chi-Squared Approximation)"
}
else {
METHOD <- "Li-Mak test on partial autocorrelations (Chi-Squared Approximation)"
}
n <- sum(!is.na(x))
if(weighted) {
weights <- (lag - (fitdf+1):lag + (fitdf+1) )/lag;
obs <- obs[(fitdf+1):lag];
STATISTIC <- n * sum(weights*obs^2);
names(STATISTIC) <- "Weighted X-squared on Squared Residuals for fitted ARCH process";
shape <- (3/4)*(lag + fitdf + 1)^2*(lag - fitdf)/(2*lag^2 + 3*lag + 2*lag*fitdf + 2*fitdf^2 + 3*fitdf + 1);
scale <- (2/3)*(2*lag^2 + 3*lag + 2*lag*fitdf + 2*fitdf^2 + 3*fitdf + 1)/(lag*(lag + fitdf + 1));
PARAMETER <- c(shape, scale);
names(PARAMETER) <- c("Shape", "Scale")
}
else {
weights <- rep(1,(lag-fitdf) );
obs <- obs[(fitdf+1):lag];
STATISTIC <- n * sum(weights*obs^2);
names(STATISTIC) <- "X-squared on Squared Residuals for fitted ARCH process"
shape <- (lag-fitdf)/2; # Chi-squared df in Gamma form.
scale <- 2;
PARAMETER <- c((lag-fitdf));
names(PARAMETER) <- c("Degrees of Freedom");
}
PVAL <- 1 - pgamma(STATISTIC, shape=shape, scale=scale)
names(PVAL) <- "Approximate p-value"
structure(list(statistic = STATISTIC, parameter = PARAMETER,
p.value = PVAL, method = METHOD, data.name = DNAME),
class = "htest")
}
# End Import
##########################################################################################
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.