#-----------------------------------------------------------
# General functions
#-----------------------------------------------------------
vech <- function(x){
t(x[!upper.tri(x)])
}
lvech <- function(x){
t(x[!upper.tri(x, diag = T)])
}
#-----------------------------------------------------------
# Functions for the GARCH estimation
#-----------------------------------------------------------
hatsigma <- function( par , x ,uncH){
# see: code underlying fGarch package of Chalaby and Wuertz
alpha = par[1] ; beta = par[2];
omega = uncH*(1-alpha-beta)
s2 = rep( uncH , T )
e2 <- x^2
e2t <- omega + alpha*c(mean(e2),e2[-length(x)])
s2 <- filter(e2t, beta, "recursive", init = uncH )
return(s2)
}
llhGarch11N <- function(par, x, uncH){
# see: code underlying fGarch package of Chalaby and Wuertz
x = as.numeric(x)
s2 = hatsigma( par = par , x = x ,uncH=uncH)
return( sum( log(s2) + x^2/s2 ))
}
weightedvariance = function( u ){
k = qchisq(0.99,df=1) ;
jumpIndicator = ( (u/mad(u))^2<=k )
ufilt = u[jumpIndicator]
uncH = mean(ufilt^2,na.rm=TRUE)*(0.99/pchisq(k,df=3))
return(uncH)
# check: U = matrix( rnorm(500000),ncol=100) ; mean(apply(U,2,'weightedvariance'))
}
garchVol = function( mR.centered , parN ){
S = c();
S_forecast = c();
N = ncol(mR.centered)
vuncH = apply( mR.centered , 2 , 'weightedvariance' )
for( i in 1:N ){
uncH = vuncH[i]
# do the garch estimation
# fitN <- try( nlminb(start=parN, objective=llhGarch11N , x = mR.centered[,i], lower=lowN, upper=upN,
# control=list(x.tol = 1e-8,trace=0)) )
# do the garch estimation with constraints to unsure covariance stationarity ie alpha + beta < 1
fitN <- try( constrOptim( theta = parN , f = llhGarch11N , ui=rbind(c(1,0),c(0,1),c(-1,0),c(0,-1),c(-1,-1) ),
ci=c(0,0,-1,-1,-1), grad=NULL, x = mR.centered[,i], uncH=uncH) )
LLH_constant = llhGarch11N( par=c(0,0), x=mR.centered[,i], uncH = uncH )
if( !inherits(fitN, "try-error") ){ LLH_garch = -fitN$value }else{ LLH_garch = LLH_constant }
# Likelihood ratio test statistic for the case of time-varying garch effects
# cases where garch volatility is probably not reliable or the volatility forecasts are so erratic that a constant garch works fine
useGARCHd = TRUE
if( LLH_garch=="NaN" | (fitN$par[1]>= 0.1&fitN$par[2]<=0.7) | sum(fitN$par)>=0.999 ){
LLH_garch = LLH_constant ;
useGARCHd = F;
} #volatility is too erratic
# Data is from inception: once we start using garch, we continue to ensure stability of the forecasts
# Likelihood ratio test for a time-varying vol model
if( 2*(LLH_garch-LLH_constant) <= qchisq(0.999,df=2) & (!useGARCHd) ){
sigmat = rep( sd( as.vector(mR.centered[,i])), length(mR.centered[,i]) );
S_forecast = c( S_forecast ,sqrt(uncH))
}else{
sigmat = sqrt(hatsigma( par = fitN$par , x = mR.centered[,i] , uncH=uncH ) )
lastsigmat = tail(sigmat,1)
if( lastsigmat=="NaN" | is.na(lastsigmat) ){
print( c("NaN day", d , "for asset" , i) );
sigmat = rep( sqrt(uncH) , length(mR.centered[,i]) );
# print(c("unc sd",sqrt(uncH) ))
S_forecast = c( S_forecast , sqrt(uncH) )
}else{
e2 = as.numeric(tail(mR.centered[,i],1)^2)
S_forecast = c( S_forecast ,
sqrt( uncH*(1-sum(fitN$par))+ fitN$par[1]*e2+ fitN$par[2]*lastsigmat^2 ) )
}
}
S = cbind( S , sigmat )
}
return( list(S,S_forecast) )
}
#-----------------------------------------------------------
# Fuction for the shrinkage correlation
#-----------------------------------------------------------
shrinkagecorrelation = function( mDevolR.centered ){
N = ncol( mDevolR.centered )
# Unconditional correlation estimate is shrinkage between sample (spearman) correlation and equicorrelation
sample = cor(mDevolR.centered,method="spearman") #make it robust by using ranks
# Elton and Gruber: same correlation across all assets
rho = mean( lvech(sample) )
I = diag( rep(1,N) ) ; J = matrix( rep(1,N^2) , ncol=N )
target = (1-rho)*I + rho*J
# The shrinkage estimate = delta*target F + (1-delta)* sample S
kappahat = shrinkage.intensity( mDevolR.centered , prior = target , sample = sample )
delta = max( 0 , min(kappahat/T,1) )
C = delta*target +(1-delta)*sample
return(C)
}
#-----------------------------------------------------------
# Functions for the DCC estimation
#-----------------------------------------------------------
DCCestimation = function( mDevolR.centered , parN , uncR ){
#fitN <- try( nlminb(start=parN, objective=loglik.dcc2.DECO , lower=lowN, upper=upN,
# control=list(x.tol = 1e-8,trace=0)) )
# Ensure covariance stationarity
N = ncol( mDevolR.centered )
I = diag( rep(1,N) ) ; J = matrix( rep(1,N^2) , ncol=N )
rho = mean( lvech(uncR) )
fitN <- try( constrOptim( theta = parN , f = loglik.dcc2.DECO , ui=rbind(c(1,0),c(0,1),c(-1,0),c(0,-1),c(-1,-1) ),
ci=c(lowN,-upN,-1), grad=NULL, uncR= uncR , mDevolR.centered=mDevolR.centered ) )
# compare with constant correlation case
LLH_constant = 0;
invtarget = (1/(1-rho))*( I - rho/( 1+(N-1)*rho )*J );
for(t in 2:nrow(mDevolR.centered)){
st = matrix( mDevolR.centered[t,] , ncol = 1 );
LLH_constant = LLH_constant + log( det(uncR) ) + t(st)%*%invtarget%*%st
}
LLH_constant = - LLH_constant
if( !inherits(fitN, "try-error") ){ LLH_DCC = -fitN$value }else{ LLH_DCC = -LLH_constant }
# Likelihood ratio test statistic for the case of time-varying garch effects
# cases where garch volatility is probably not reliable or the volatility forecasts are so erratic that a constant garch works fine
if( LLH_DCC=="NaN" | (fitN$par[1]>= 0.1&fitN$par[2]<=0.7) | sum(fitN$par)>=0.999 ){
LLH_DCC = LLH_constant ;
useDCC = F;
} #volatility is too erratic
if( 2*(LLH_DCC-LLH_constant) <= qchisq(0.999,df=2) ){
Rt = uncR;
}else{
print( c("use DCC with parameters",fitN$par) )
Rt = Rforecast(fitN$par,mR=mDevolR.centered,uncR=C)
}
return(Rt)
}
Rforecast <- function(par,mR,uncR){
T = nrow(mR); N = ncol(mR)
alpha = par[1] ; beta = par[2];
omega = (1-par[1]-par[2])*uncR
Qt = diag( rep(1,N) )
k = 1
for( i in 2:N ){
for(j in 1:(i-1)){
e2 = mR[,i]*mR[,j];
e2[is.na(e2)] = mean( na.omit(e2 ) )
e2t <- omega[i,j] + alpha*c(mean(e2),e2[-length(e2)], tail(e2,1) )
Qt[i,j] = Qt[j,i] = tail(filter(e2t, beta, "recursive", init = uncR[i,j] ),1)
}
}
return( diag(diag(Qt)^(-0.5))%*%Qt%*%diag(diag(Qt)^(-0.5)) )
}
# test; uncR=diag(c(1,1)) ; Rforecast( c(0.05,0.4) , matrix(rnorm(400),ncol=2),uncR)
loglik.dcc2.DECO = function (param,uncR,mDevolR.centered)
{
require('ccgarch')
dvar = mDevolR.centered
nobs <- dim(dvar)[1]
ndim <- dim(dvar)[2]
out <- .Call("dcc_est", dvar, uncR, param[1], param[2])
DCC <- dcc.est(dvar, param)$DCC
lf <- numeric(ndim)
I = diag( rep(1,ndim) ) ; J = matrix( rep(1,ndim^2) , ncol=ndim )
for (i in 1:nobs) {
R <- matrix(DCC[i, ], ndim, ndim)
rho = mean( lvech(R) )
DECO = (1-rho)*I + rho*J
invDECO = (1/(1-rho))*( I - rho/( 1+(ndim-1)*rho )*J )
lf[i] <- log(det(DECO)) + sum(dvar[i, ] * crossprod(invDECO,
dvar[i, ]))
}
sum(lf)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.