#BNEE1 vs BVN vs BVT on GK Volatility (SPY, QQQ)
pkgs <- c("mvtnorm","quantmod"); for(p in pkgs) if(!requireNamespace(p,quietly=TRUE)) install.packages(p)
suppressPackageStartupMessages({library(mvtnorm);library(quantmod)})
set.seed(2025)

#Data (Yahoo Finance) & GK volatility (positive)
from <- "2015-01-01"; to <- Sys.Date()
get_or_stop <- function(s) {x <- try(getSymbols(s, from=from, to=to, src="yahoo", auto.assign=FALSE), silent=TRUE)
if(inherits(x,"try-error")) stop(sprintf("Yahoo fetch failed: %s",s)); x}
SPY <- get_or_stop("SPY"); QQQ <- get_or_stop("QQQ")
gk_vol <- function(ohlc){
  O<-Op(ohlc);H<-Hi(ohlc);L<-Lo(ohlc);C<-Cl(ohlc)
  v <- 0.5*(log(H/L))^2 - (2*log(2)-1)*(log(C/O))^2
  v[v<=0] <- 1e-12; sqrt(v)
}
V <- na.omit(merge(gk_vol(SPY), gk_vol(QQQ)))
X <- as.numeric(V[,1]); Y <- as.numeric(V[,2])
n <- length(X); cat(sprintf("n = %d\n", n))

.l1mexp <- function(z) ifelse(z<=log(2), log(-expm1(-z)), log1p(-exp(-z)))
.lexpm1 <- function(t) ifelse(t>30, t, log(expm1(t)))    # log(exp(t)-1)

dBNEE1 <- function(x,y, m1,m2, s1,s2, r, a,b, log=FALSE){
  x<-as.numeric(x); y<-as.numeric(y); n<-max(length(x),length(y)); x<-rep_len(x,n); y<-rep_len(y,n)
  lpdf <- rep(-Inf,n)
  bad <- !(is.finite(m1)&is.finite(m2)&is.finite(s1)&s1>0&is.finite(s2)&s2>0&is.finite(r)&r>-1&r<1&is.finite(a)&a>0&is.finite(b)&b>0)
  if(bad) return(if(log) lpdf else exp(lpdf))
  ii <- which(is.finite(x)&is.finite(y)&x>0&y>0); if(!length(ii)) return(if(log) lpdf else exp(lpdf))
  ax <- a*x[ii]; by <- b*y[ii]; keep <- ii[is.finite(ax)&is.finite(by)&ax<=800&by<=800]
  if(!length(keep)) return(if(log) lpdf else exp(lpdf))
  ii <- keep; ax <- a*x[ii]; by <- b*y[ii]
  u <- .lexpm1(ax); v <- .lexpm1(by)
  lj <- log(a) - .l1mexp(ax) + log(b) - .l1mexp(by)
  S <- matrix(c(s1^2, r*s1*s2, r*s1*s2, s2^2),2,2)
  lbn <- tryCatch(mvtnorm::dmvnorm(cbind(u,v), mean=c(m1,m2), sigma=S, log=TRUE), error=function(e) rep(-Inf,length(ii)))
  lpdf[ii] <- lj + lbn
  if(log) lpdf else exp(lpdf)
}

pBNEE1_X <- function(q, m1,s1,a){ out<-rep(0,length(q)); ii<-which(q>0&is.finite(q)); if(length(ii)) out[ii]<-pnorm((.lexpm1(a*q[ii])-m1)/s1); out }
pBNEE1_Y <- function(q, m2,s2,b){ out<-rep(0,length(q)); ii<-which(q>0&is.finite(q)); if(length(ii)) out[ii]<-pnorm((.lexpm1(b*q[ii])-m2)/s2); out }

negloglik_BNEE1 <- function(par,x,y){
  m1<-par[1]; m2<-par[2]; s1<-exp(par[3]); s2<-exp(par[4]); r<-tanh(par[5]); a<-exp(par[6]); b<-exp(par[7])
  if(!is.finite(s1)|!is.finite(s2)|s1<1e-6|s2<1e-6|s1>1e6|s2>1e6|!is.finite(a)|!is.finite(b)|a<1e-6|b<1e-6|a>1e3|b>1e3) return(1e12)
  axM<-a*max(x); byM<-b*max(y); if(!is.finite(axM)|!is.finite(byM)|axM>800|byM>800) return(1e12)
  lp <- dBNEE1(x,y,m1,m2,s1,s2,r,a,b,log=TRUE); if(any(!is.finite(lp))) return(1e11); -sum(lp)
}
init_BNEE1 <- function(x,y){
  a0<-1/mean(x); b0<-1/mean(y); u0<-.lexpm1(a0*x); v0<-.lexpm1(b0*y)
  c(mean(u0), mean(v0), log(pmax(sd(u0),1e-3)), log(pmax(sd(v0),1e-3)), atanh(pmin(pmax(cor(u0,v0),-0.99),0.99)), log(a0), log(b0))
}
fit_BNEE1 <- function(x,y, nstart=5){
  best<-NULL; val<-Inf
  for(s in 1:nstart){
    p0<-init_BNEE1(x,y); p0[6]<-log(exp(p0[6])*runif(1,.7,1.3)); p0[7]<-log(exp(p0[7])*runif(1,.7,1.3))
    o<-optim(p0, negloglik_BNEE1, x=x,y=y, method="L-BFGS-B", control=list(reltol=1e-9,maxit=6e4,parscale=rep(1,7)))
    if(o$value<val){best<-o; val<-o$value}
  }
  p<-best$par; k<-7L; list(mu1=p[1],mu2=p[2],sigma1=exp(p[3]),sigma2=exp(p[4]),rho=tanh(p[5]),a=exp(p[6]),b=exp(p[7]),
                           logLik=-val, AIC=2*k-2*(-val), BIC=k*log(n)-2*(-val))
}

fit_BVN <- function(x,y){
  S <- cov(cbind(x,y)) * ((n-1)/n); mu <- c(mean(x),mean(y))
  ll <- sum(dmvnorm(cbind(x,y), mean=mu, sigma=S, log=TRUE))
  s1<-sqrt(S[1,1]); s2<-sqrt(S[2,2]); r<-S[1,2]/(s1*s2); k<-5L
  list(mu1=mu[1],mu2=mu[2],sigma1=s1,sigma2=s2,rho=r, logLik=ll, AIC=2*k-2*ll, BIC=k*log(n)-2*ll)
}
pBVN_X <- function(q, m,s) pnorm((q-m)/s); pBVN_Y <- function(q, m,s) pnorm((q-m)/s)

negloglik_BVT <- function(par,x,y){
  m1<-par[1]; m2<-par[2]; s1<-exp(par[3]); s2<-exp(par[4]); r<-tanh(par[5]); nu<-exp(par[6])+2
  if(!is.finite(s1)|!is.finite(s2)|s1<=0|s2<=0|!is.finite(r)|abs(r)>=1|!is.finite(nu)|nu<=2) return(1e12)
  S<-matrix(c(s1^2,r*s1*s2,r*s1*s2,s2^2),2,2)
  lv<-tryCatch(dmvt(cbind(x,y), delta=c(m1,m2), sigma=S, df=nu, log=TRUE), error=function(e) rep(-Inf,length(x)))
  if(any(!is.finite(lv))) return(1e11); -sum(lv)
}
init_BVT <- function(x,y){ c(mean(x),mean(y), log(pmax(sd(x),1e-6)), log(pmax(sd(y),1e-6)), atanh(pmin(pmax(cor(x,y),-0.99),0.99)), log(8-2)) }
fit_BVT <- function(x,y, nstart=5){
  best<-NULL; val<-Inf
  for(s in 1:nstart){
    p0<-init_BVT(x,y); p0[3]<-log(exp(p0[3])*runif(1,.8,1.25)); p0[4]<-log(exp(p0[4])*runif(1,.8,1.25)); p0[6]<-log(exp(p0[6])*runif(1,.7,1.4))
    o<-optim(p0, negloglik_BVT, x=x,y=y, method="L-BFGS-B", control=list(reltol=1e-9,maxit=6e4,parscale=rep(1,6)))
    if(o$value<val){best<-o; val<-o$value}
  }
  p<-best$par; k<-6L; list(mu1=p[1],mu2=p[2],sigma1=exp(p[3]),sigma2=exp(p[4]),rho=tanh(p[5]),nu=exp(p[6])+2,
                           logLik=-val, AIC=2*k-2*(-val), BIC=k*log(n)-2*(-val))
}
pBVT_X <- function(q,m,s,nu) pt((q-m)/s, df=nu); pBVT_Y <- function(q,m,s,nu) pt((q-m)/s, df=nu)

bnee1 <- fit_BNEE1(X,Y); bvn <- fit_BVN(X,Y); bvt <- fit_BVT(X,Y)
ksX_bnee1 <- ks.test(X, function(q) pBNEE1_X(q, bnee1$mu1, bnee1$sigma1, bnee1$a))
ksY_bnee1 <- ks.test(Y, function(q) pBNEE1_Y(q, bnee1$mu2, bnee1$sigma2, bnee1$b))
ksX_bvn   <- ks.test(X, function(q) pBVN_X(q, bvn$mu1, bvn$sigma1))
ksY_bvn   <- ks.test(Y, function(q) pBVN_Y(q, bvn$mu2, bvn$sigma2))
ksX_bvt   <- ks.test(X, function(q) pBVT_X(q, bvt$mu1, bvt$sigma1, bvt$nu))
ksY_bvt   <- ks.test(Y, function(q) pBVT_Y(q, bvt$mu2, bvt$sigma2, bvt$nu))

cmp <- data.frame(
  Model=c("BNEE1","BVN","BVT"),
  k     =c(7L,5L,6L),
  logLik=c(bnee1$logLik, bvn$logLik, bvt$logLik),
  AIC   =c(bnee1$AIC,    bvn$AIC,    bvt$AIC),
  BIC   =c(bnee1$BIC,    bvn$BIC,    bvt$BIC),
  KS_X_D=c(ksX_bnee1$statistic, ksX_bvn$statistic, ksX_bvt$statistic),
  KS_X_p=c(ksX_bnee1$p.value,   ksX_bvn$p.value,   ksX_bvt$p.value),
  KS_Y_D=c(ksY_bnee1$statistic, ksY_bvn$statistic, ksY_bvt$statistic),
  KS_Y_p=c(ksY_bnee1$p.value,   ksY_bvn$p.value,   ksY_bvt$p.value)
)

cat("\n===== Comparison Table (BNEE1 vs BVN vs BVT) =====\n"); print(cmp, row.names=FALSE)


<pre>An Error occurred while handling another error:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/ErrorHandler.php(136): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/ErrorHandler.php(135): yii\web\ErrorHandler-&gt;renderException()
#3 [internal function]: yii\base\ErrorHandler-&gt;handleException()
#4 {main}
Previous exception:
yii\web\HeadersAlreadySentException: Headers already sent in  on line 0. in /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php:366
Stack trace:
#0 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Response.php(339): yii\web\Response-&gt;sendHeaders()
#1 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(656): yii\web\Response-&gt;send()
#2 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/models/Files.php(696): yii\base\Application-&gt;end()
#3 /var/www/html/prof-homepages/vendor/faravaghi/yii2-filemanager/controllers/FilesController.php(484): faravaghi\filemanager\models\Files-&gt;getFile()
#4 [internal function]: faravaghi\filemanager\controllers\FilesController-&gt;actionGetFile()
#5 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/InlineAction.php(57): call_user_func_array()
#6 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Controller.php(180): yii\base\InlineAction-&gt;runWithParams()
#7 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Module.php(528): yii\base\Controller-&gt;runAction()
#8 /var/www/html/prof-homepages/vendor/yiisoft/yii2/web/Application.php(103): yii\base\Module-&gt;runAction()
#9 /var/www/html/prof-homepages/vendor/yiisoft/yii2/base/Application.php(386): yii\web\Application-&gt;handleRequest()
#10 /var/www/html/prof-homepages/backend/web/index.php(16): yii\base\Application-&gt;run()
#11 {main}</pre>