I'm currently building a portfolio optimization algorithm, with a core step being moment-matching scenario generation. For simplicity and efficiency, I'm implementing the method from the paper An algorithm for moment-matching scenario generation with application to financial portfolio optimization by Ponomareva, Roman and Date. While the mathematical logic in the paper is clear, I'm stuck on a critical issue: some of the probability weights p are turning negative, even though the paper's formulas should guarantee non-negative weights. I tried looping the algorithm repeatedly to find an all-positive weight set, but that just leads to infinite loops.

Context #

I'm working with a data.frame of monthly returns for 11 funds, spanning from January 30, 2001 to September 30, 2020. I've already calculated the return means, covariance matrix, centered third moments, and centered fourth moments (note: these are not skewness and kurtosis). The problem arises specifically in the p array—its first s elements are used as probabilities in the final P vector, but negative values keep popping up, violating the non-negativity requirement for probabilities.

Problematic Code Snippet #

dummy1 = 0 
while (dummy1 <=0 | dummy1 >= 1) { 
  dummy1 = round(rnorm(1, mean = 0.5, sd = 0.25), 2) 
} 
diag.cov.returns = diag(cov.returns) 
Z = dummy1 * sqrt (diag.cov.returns) #Vector Z according to paper formula 
ZZT = Z %*% t(Z) 
LLT = cov.returns - ZZT 
L = chol(LLT) #cholesky decomposition to get matrix L 
s = sample (1:5, 1) 
F1 = 0 
F2 = -1 
S = (2*N*s)+3 
while (((4*F2)-(3*F1*F1)) < 0) { 
  #Gamma = (2*s*s)*(((N*mean.fourth) - (0.75*(sum(Z^4)* (N*mean.third/sum(Z^3))^2)))/sum(L^4)) 
  #Gamma is necessary if we want to get p from Uniform Distribution 
  #U = runif(s, 0, 1) 
  U = rgamma(s, shape = 1, scale = ((1/exp(1)):1)) 
  #p = (s*(N/Gamma)) + ((1/(2*N*s)) - (s/(N*Gamma)))*U 
  p = (-log(U, base = exp(1))) 
  p = p/(((2*sum(p))+max(p))*N*s) #this is the array expected to have positive and bounded between 0 and 1 
  q1 = 1/p 
  pz = p 
  p[s+1] = (1-(2*N*sum(p))) #extra point necessary to get the 3 moment matching probabilities 
  F1 = (N*mean.third*sqrt(p[s+1]))/(sum(Z^3)) 
  F2 = p[s+1]*(((N*mean.fourth) - (1/(2*s*s))*sum(L^4)*(sum(1/p)))/sum(Z^4)) 
} 
alpha = (0.5*F1) + 0.5*sqrt((4*F2)-(3*F1*F1)) 
beta = -(0.5*F1) + 0.5*sqrt((4*F2)-(3*F1*F1)) 
w1 = 1/(alpha*(alpha+beta)) 
w2 = 1/(beta*(alpha+beta)) 
w0 = 1 - (1/(alpha*beta)) 
P = rep(pz, 2*N) #Vector with Probabilities starting from p + 3 extra probabilities to match third and fourth moments 
P[(2*N*s)+1] = p[s+1]*w0 
P[(2*N*s)+2] = p[s+1]*w1 
P[(2*N*s)+3] = p[s+1]*w2 

Additional Notes #

• I can't share the exact return dataset, but any data.frame of asset returns seems to reproduce this issue.
• The paper's public dataset is available, but I haven't imported it into R yet to verify—though the problem persists with my own data.

I've tried looping to re-run the algorithm until all p values are positive, but that never terminates. Does anyone have insights into why the weights are turning negative, or how to adjust the implementation to enforce non-negativity as the paper intends?

内容的提问来源于stack exchange，提问作者Robert Butu