gh_roots <- function(skew, kurt, m = 10, tries = 100, no.digits = 2, method="Newton"){ library('rootSolve') model <- function(x, parms = c(0,0)){ return(c(F1 = ((3*exp(x[1]^2 / (2-6*x[2])) + exp(9*x[1]^2 / (2-6*x[2])) - 3*exp(2*x[1]^2 /(1 - 3*x[2])) - 1) / (1-3*x[2])^0.5 - 3*(1 - 2*exp(x[1]^2 /(2-4*x[2]))+ exp(2*x[1]^2/(1-2*x[2])))*(exp(x[1]^2 / (2-2*x[2])) -1)/((1-2*x[2])^0.5 * (1-x[2])^0.5) + 2*(exp(x[1]^2/(2-2*x[2])) -1)^3 / (1-x[2])^1.5) / (x[1]^3 * (((1-2*exp(x[1]^2 /(2-4*x[2])) + exp(2*x[1]^2 / (1-2*x[2]))) / (1-2*x[2])^0.5 + (exp(x[1]^2/(2-2*x[2])) -1)^2 / (x[2] -1))/x[1]^2)^1.5) - parms[1], ### The formula to calculate the kurtosis F2 = (exp(8*x[1]^2/(1-4*x[2]))*(1 + 6*exp(6*x[1]^2/(4*x[2]-1)) + exp(8*x[1]^2/(4*x[2]-1)) - 4*exp(7*x[1]^2/(8*x[2]-2)) - 4*exp(15*x[1]^2/(8*x[2]-2)))/(1-4*x[2])^0.5 - 4*(3*exp(x[1]^2/(2-6*x[2])) + exp(9*x[1]^2/(2-6*x[2])) - 3*exp(2*x[1]^2/(1-3*x[2]))-1)* (exp(x[1]^2/(2-2*x[2]))-1)/((1-3*x[2])^0.5 * (1-x[2])^04) - 6*(exp(x[1]^2/(2-2*x[2]))-1)^4/(x[2]-1)^2 -12*(1-2*exp(x[1]^2/(2-4*x[2])) + exp(2*x[1]^2/(1-2*x[2])))*(exp(x[1]^2/(2-2*x[2]))-1)^2/((1-2*x[2])^0.5*(x[2]-1))+ 3*(1-2*exp(x[1]^2/(2-4*x[2])) + exp(2*x[1]^2/(1-2*x[2])))^2/(2*x[2]-1))/ ((1-2*exp(x[1]^2/(2-4*x[2])) + exp(2*x[1]^2/(1-2*x[2])))/(1-2*x[2])^0.5 + (exp(x[1]^2/(2-2*x[2]))-1)^2/(x[2]-1))^2 - parms[2])) } mu_gh <- function(g,h) (exp(g^2 *(2-2*h)^-1)-1)/(g*(1-h)^0.5) variance_gh <- function(g,h) ((1-2*exp(g^2 *(2 - 4*h)^-1) + exp(2*g^2 * ( 1-2*h)^-1))/(1-2*h)^0.5 + (exp(g^2 *(2 -2*h)^-1)-1)^2 / (h-1))/g^2 omega <- function(x, parms) sum(model(x,parms)^2) if (method=="Newton") { # Newton's Homotopy Function homotopy = function(x,x0,t,parms) model(x,parms) - (1-t)*model(x0,parms) } else { # Fixed point homotopy homotopy = function(x,x0,t,parms) t*model(x,parms) + (1-t)*(x-x0) } homotopy.root = function(fn = homotopy, x0, m = 30, parms){ delta = 1/m t=0 for (k in 1:m){ t = t + delta x0 = multiroot(function(x) homotopy(x,x0,t,parms),start = x0)\$root #print(paste(t,"---",x0[1],",",x0[2],"---",model(x0)[1],",",model(x0)[2])) if (is.na(model(x0)[1])){ break } } return(x0) } hs = gs = solutions = list() g <- c(seq(-100,-1,length.out = 50),seq(-1,1,length.out = 300),seq(1,100,length.out = 50)) h <- c(seq(-100,0,length.out = 100),seq(0,0.24,length.out = 100)) X <- as.matrix(expand.grid(x=g,y=h)) for (i in 1:length(skew)){ skewness = skew[i] kurtosis = kurt[i] parms = c(skewness,kurtosis) om = apply(X,1,function(x) omega(x,parms)) ind <- which(om %in% sort(unique(om))[1:tries]) X_guess = X[ind,] sols = matrix(nrow = tries,ncol = 6) for (j in 1:tries){ x0 = X_guess[j,] x0 = homotopy.root(homotopy,x0,m,parms) sols[j,] = c(x0,model(x0),mu_gh(x0[1],x0[2]),variance_gh(x0[1],x0[2])) } #sols = (unique(na.omit(round(sols,no.digits)))) sols = na.omit(sols) dup = duplicated(round(sols,digits = no.digits)) sols = matrix(sols[!dup,],ncol=6) if (length(sols) ==0) { print(paste("Failed to find solutions for a skewness of",skewness,"and a kurtosis of ",kurtosis)) } else { print(paste(length(sols[,1]), "solution(s) found for a skewness of ",skewness,"and a kurtosis of ",kurtosis)) } gs[[i]] = sols[,1] hs[[i]] = sols[,2] solutions[[i]] = sols } gs = t(expand.grid(gs)) hs = t(expand.grid(hs)) out = list(g = gs, h = hs, solutions = solutions ) return (out) } ##################################################################################################################################################################### ##################################################################################################################################################################### ##################################################################################################################################################################### ##################################################################################################################################################################### ##################################################################################################################################################################### multgh <- function(N=30,g,h,Sigma,warnings = T){ library('rootSolve') ## For implementation of Newton's method for roots approximation. library('lattice') ## For making a grid of points to look for convenient guesses to supply in the Newton's method library('pracma') ## For numerical integration library('Matrix') ## For matrix operations model <- function(x, parms = c(0,0)){ c(F1 = ((3*exp(x[1]^2 / (2-6*x[2])) + exp(9*x[1]^2 / (2-6*x[2])) - 3*exp(2*x[1]^2 /(1 - 3*x[2])) - 1) / (1-3*x[2])^0.5 - 3*(1 - 2*exp(x[1]^2 /(2-4*x[2]))+ exp(2*x[1]^2/(1-2*x[2])))*(exp(x[1]^2 / (2-2*x[2])) -1)/((1-2*x[2])^0.5 * (1-x[2])^0.5) + 2*(exp(x[1]^2/(2-2*x[2])) -1)^3 / (1-x[2])^1.5) / (x[1]^3 * (((1-2*exp(x[1]^2 /(2-4*x[2])) + exp(2*x[1]^2 / (1-2*x[2]))) / (1-2*x[2])^0.5 + (exp(x[1]^2/(2-2*x[2])) -1)^2 / (x[2] -1))/x[1]^2)^1.5) - parms[1], F2 = (exp(8*x[1]^2/(1-4*x[2]))*(1 + 6*exp(6*x[1]^2/(4*x[2]-1)) + exp(8*x[1]^2/(4*x[2]-1)) - 4*exp(7*x[1]^2/(8*x[2]-2)) - 4*exp(15*x[1]^2/(8*x[2]-2)))/(1-4*x[2])^0.5 - 4*(3*exp(x[1]^2/(2-6*x[2])) + exp(9*x[1]^2/(2-6*x[2])) - 3*exp(2*x[1]^2/(1-3*x[2]))-1)* (exp(x[1]^2/(2-2*x[2]))-1)/((1-3*x[2])^0.5 * (1-x[2])^0.5) - 6*(exp(x[1]^2/(2-2*x[2]))-1)^4/(x[2]-1)^2 -12*(1-2*exp(x[1]^2/(2-4*x[2])) + exp(2*x[1]^2/(1-2*x[2])))*(exp(x[1]^2/(2-2*x[2]))-1)^2/((1-2*x[2])^0.5*(x[2]-1))+ 3*(1-2*exp(x[1]^2/(2-4*x[2])) + exp(2*x[1]^2/(1-2*x[2])))^2/(2*x[2]-1))/ ((1-2*exp(x[1]^2/(2-4*x[2])) + exp(2*x[1]^2/(1-2*x[2])))/(1-2*x[2])^0.5 + (exp(x[1]^2/(2-2*x[2]))-1)^2/(x[2]-1))^2 - parms[2])} mu_gh <- function(g,h) (exp(g^2 *(2-2*h)^-1)-1)/(g*(1-h)^0.5) variance_gh <- function(g,h) ((1-2*exp(g^2 *(2 - 4*h)^-1) + exp(2*g^2 * ( 1-2*h)^-1))/(1-2*h)^0.5 + (exp(g^2 *(2 -2*h)^-1)-1)^2 / (h-1))/g^2 q_z <- function(g,h,z) (exp(g*z) - 1)*exp(0.5*h*z^2)/g mu <- mu_gh(g,h) sd_s <- sqrt(variance_gh(g,h)) Sigma_intermediate <- matrix(nrow=length(h),ncol=length(h)) for (i in 1:length(g)){ for (j in (1:length(g))){ if (i > j) { q_i <- function(z) (q_z(g[i],h[i],z) - mu[i]) / sd_s[i] q_j <- function(z) (q_z(g[j],h[j],z) - mu[j]) / sd_s[j] f_ij <- function(z_i,z_j,pij) (1/(2*pi*sqrt(1-pij^2))) * exp(-(2*(1-pij^2))^-1 *(z_i^2-2*pij*z_i*z_j +z_j^2)) p_ij <- function(z_i,z_j,pij) function(z_i, z_j) q_i(z_i)*q_j(z_j)*f_ij(z_i,z_j,pij) ## The solver of this library implements a different algorithm (Gauss-Kronrod) than the paper's one (Adaptative Monte Carlo) #Z <- seq(-8,8,length.out = 10000) #D_i <- sapply(Z,q_i) #D_j <- sapply(Z, q_j) #xmin <- quantile(D_i,0.0) #xmax <- quantile(D_i,1) #ymin <- quantile(D_j,0.0) #ymax <- quantile(D_j,1) xmin <- -8 xmax <- 8 ymin <- -8 ymax <- 8 F_ij <- function(x,parms) (integral2(fun = p_ij(z_i,z_j,pij=x[1]), xmin = xmin,xmax = xmax, ymin= ymin, ymax = ymax)\$Q - parms[1]) #x0 <- optimize(f = function(x) F_ij(x,Sigma[i,j]),lower = -0.99,upper = 0.99)\$minimum #intermediate_corr <- uniroot(f = function(x) F_ij(x,Sigma[i,j]),lower = -0.99,upper = 0.99) #Sigma_intermediate[i,j] <- intermediate_corr\$root intermediate_corr <- optimize(f = function(x) F_ij(x,Sigma[i,j])^2,lower = -0.99,upper = 0.99) Sigma_intermediate[i,j] <- intermediate_corr\$minimum } } } diag(Sigma_intermediate) <- 1 # set the diagonal to ones Sigma_intermediate <- forceSymmetric(Sigma_intermediate,uplo="L") # as the correlation matrix is a sym. matrix. if (mean(eigen(x = Sigma_intermediate)\$values>0) < 1){ # A way to evaluate if at least one eigenvalue is non-positive Sigma_intermediate <- nearPD(Sigma_intermediate) if (isTRUE(warning)) print('Warning!!! Non-positive definite matrix found. A correction was made to make it Positive Definite') } chol_decomposition <- chol(Sigma_intermediate) # Cholesky Decomposition V <- matrix(rnorm(length(h)*N),nrow=N,ncol=length(h)) # A matrix of standard normal values Z <- V %*% t(chol_decomposition) # To apply formula (10) output <- matrix(nrow = nrow(Z),ncol = ncol(Z)) for (i in 1:ncol(Z)) { output[,i] <- q_z(g[i], h[i],Z[,i]) } return (list(data = output, intermediate_matrix =Sigma_intermediate, chol = chol_decomposition, g=g,h=h)) }