# How to check the correctness of calculations with a gamma distribution?

I’m reading Ponomareva, Roman, and Date (2015) and trying to generate vector $P$ of the $2Ns + 3$ probability weights: $$P =\{\underbrace{p_1, p_2, \ldots, p_s, p_1, p_2, \ldots, p_s,p_1, p_2, \ldots, p_s}_{2\cdot N \text{ times}} , p_{s+1}\cdot w_0, p_{s+1}\cdot w_1, p_{s+1}\cdot w_2\}.$$ In the paper the value $s$ equals to the one number from the set $\{1, 3, 6, 9, 15, 18, 27, 126\}$, and $N=20$. I had a problem because some generated probabilities $P$ were negative. In private communication with the one author it was confirmed that this is a problem.

Authors proposed the three-steps method. I included a short discription of first and second steps only.

Step 1: Parameters for scenario generation

Scenario generation procedure has the following inputs: mean $\mu$ of the random vector $r$, vector’s covariance matrix $\Sigma$, dimension $N$ and average marginal moments $\bar{\kappa}$ and $\bar{\zeta}$:

N <- 20

# input data from the paper (Table A.1, A.2, A.3)

# mean
mu <- c(6.27,25.21,8.9,13.79,7.06,38.6,37.67,36.97,13.08,19.25,
4.43,16.6,46.45,17.22,29.19,18.98,16.33,21.79,20.65,17.72)
mu <- mu * 10^(-4)

# skewness
kappa <- c(-8.36,25.34,1.74,5.29,9.27,74.39,86.66,26.93,-3.76,-3.81,
-29.87, -63.57, 191.96, -13.48, 7.29, -8.86, 3.33, 285.29, 157.21, 42.79)
kappa <- kappa * 10^(-6)

# kurtosis
zeta <- c(26.5, 24.22, 11.6, 14.93, 10.32, 44.83, 133.18, 23.55, 15.04, 7.93,
21.5, 29.4, 121.67, 12.31, 31.28, 9.97, 6.96, 459.01, 87.35, 23.84)
zeta <- zeta * 10^(-6)

m_skew <- mean(kappa)
m_kurt <- mean(zeta)

# Covariance Matrix

cov <-c(18.58, 6.18, 7.14, 7.03, 4.48, 3.47, 10.37, 5.21, 4.22, 4.42, 7.11, 4.58, 10.55, 4.93, 5.76, 2.95, 4.52, 16.51, 6.73, 6.16,
6.18, 22.69, 3.03, 3.07, 5.01, 2.89, 5.13, 3.44, 5.09, 1.28, 3.4, 2.37, 4.94, 2.26, 4.72, 5.43, 2.33, 7.71, 8.89, 8.01,
7.14, 3.03, 15.15, 11.67, 4.17, 3.59, 10.65, 7.79, 3.79, 3.04, 4.96, 3.16, 10.99, 3.95, 3.8, 2.89, 3.37, 9.31, 4.08, 4.45,
7.03, 3.07, 11.67, 16.58, 4.77, 4.04, 11.87, 7.94, 4.33, 3.64, 4.93, 3.42, 12.93, 5.64, 3.92, 3.13, 4.28, 8.47, 4.42, 4.49,
4.48, 5.01, 4.17, 4.77, 13.25, 3.36, 2.94, 3.18, 8.54, 3.52, 2.51, 3.05, 2.14, 4.44, 3.2, 4.03, 4.19, 7.36, 3.02, 3.27,
3.47, 2.89, 3.59, 4.04, 3.36, 17.6, 4.54, 3.37, 3.07, 4.43, 2.99, 3.09, 4.06, 5.42, 1.54, 3.75, 5.24, 4.24, 1.13, 2.49,
10.37, 5.13, 10.65, 11.87, 2.94, 4.54, 37.81, 12.57, 2.56, 3.66, 6.96, 4.44, 25.27, 4.36, 5.75, 2.65, 4.3, 13.66, 4.87, 6.55,
5.21, 3.44, 7.79, 7.94, 3.18, 3.37, 12.57, 18.32, 3.62, 3.64, 6.31, 3.5, 11.52, 3.26, 4.06, 3.71, 3.73, 8.99, 4.16, 3.66,
4.22, 5.09, 3.79, 4.33, 8.54, 3.07, 2.56, 3.62, 16.04, 3.92, 2.68, 3.5, 2.17, 4.18, 2.92, 4.47, 3.61, 6.91, 3.74, 4.17,
4.42, 1.28, 3.04, 3.64, 3.52, 4.43, 3.66, 3.64, 3.92, 11.43, 3.89, 3.24, 3.64, 4.03, 4.46, 2.85, 3.44, 5.48, 1.89, 2.7,
7.11, 3.4, 4.96, 4.93, 2.51, 2.99, 6.96, 6.31, 2.68, 3.89, 14.97, 4.24, 7.4, 3.05, 6.09, 3.41, 2.28, 11.38, 4.54, 4,
4.58, 2.37, 3.16, 3.42, 3.05, 3.09, 4.44, 3.5, 3.5, 3.24, 4.24, 19.85, 4.71, 4.98, 4.02, 2.5, 4.36, 7.32, 0.96, 2.4,
10.55, 4.94, 10.99, 12.93, 2.14, 4.06, 25.27, 11.52, 2.17, 3.64, 7.4, 4.71, 34.81, 4.25, 5.79, 2.19, 4.51, 13.44, 6.9, 5.63,
4.92, 2.26, 3.95, 5.64, 4.44, 5.42, 4.36, 3.26, 4.18, 4.03, 3.05, 4.98, 4.25, 14.12, 4.34, 2.46, 6.24, 6.3, 2.29, 2.49,
5.76, 4.72, 3.8, 3.92, 3.2, 1.54, 5.75, 4.06, 2.92, 4.46, 6.09, 4.02, 5.79, 4.34, 22.62, 2.84, 3.28, 9.1, 5.02, 4.01,
2.95, 5.43, 2.89, 3.13, 4.03, 3.75, 2.65, 3.71, 4.47, 2.85, 3.41, 2.5, 2.19, 2.46, 2.84, 10.65, 3, 5.11, 4.23, 3.79,
4.52, 2.33, 3.37, 4.28, 4.19, 5.24, 4.3, 3.73, 3.61, 3.44, 2.28, 4.36, 4.51, 6.24, 3.28, 3, 11.76, 6.9, 1.36, 1.77,
16.51, 7.71, 9.31, 8.47, 7.36, 4.24, 13.66, 8.99, 6.91, 5.48, 11.38, 7.32, 13.44, 6.3, 9.1, 5.11, 6.9, 46.21, 10.27, 8.75,
6.73, 8.89, 4.07, 4.42, 3.02, 1.13, 4.87, 4.16, 3.74, 1.89, 4.54, 0.96, 6.9, 2.29, 5.02, 4.23, 1.36, 10.27, 28.29, 12.14,
6.16, 8.01, 4.45, 4.49, 3.27, 2.49, 6.55, 3.66, 4.17, 2.7, 4, 2.4, 5.63, 2.49, 4.01, 3.79, 1.77, 8.75, 12.14, 19.22)

Sigma = matrix(cov, ncol=N)
Sigma <- Sigma * 10^(-4)


User needs to choose an arbitrary positive integer $s$, an arbitrary nonzero deterministic vector $Z$, such that $\Sigma - ZZ^\top> 0$, and a scalar $\rho \in (0,1)$:

s <- 9
pho <- 0.45
Z   <- pho * sqrt(diag(Sigma))

ZZT <- Z %*% t(Z)
# Checking the choice of the nonzero deterministic vector Z
ifelse(sum(Sigma - ZZT)>0, print("Sigma - ZZT>0"), print("Sigma - ZZT<0"))
LLT <- Sigma - ZZT

L <- chol(LLT)

# Checking the loss of significance of Cholesky factorization

if(sum(round(Sigma - LLT - ZZT, 10)) == 0) print("TRUE")


The remaining parameters $\phi_1$, $\phi_2$, $\alpha$, $\beta$, $w_0$, $w_1$, $w_2$, $\gamma$ can be calculated as follows:

$$\phi_1 = \frac{N \bar{\kappa} \sqrt{p_{s+1}} } {\sum_{j=1}^N Z_j^3},$$

$$\phi_2 = p_{s+1} \frac{N \bar{\zeta} - \frac{1}{2s^2}\sum_{l,k}L^4_{l,k}\sum_{i=1}^s\frac{1}{p_i}}{\sum_{j=1}^N{Z_j^4}},$$

$\alpha = 0.5\phi_1 + 0.5\sqrt{4\phi_2-3\phi_1^2}$, $\beta = -0.5\phi_1 + 0.5\sqrt{4\phi_2-3\phi_1^2}$,

$w_1 = \frac{1}{\alpha(\alpha+\beta)}$, $w_2 = \frac{1}{\beta(\alpha+\beta)}$, $w_0 = 1 - \frac{1}{\alpha \beta}$,

$$\gamma=2s^2\frac{N\bar{\zeta} - \frac{3}{4} \sum_{j=1}^N Z_j^4\bigg(\frac{N\bar{\kappa}}{\sum_{j=1}^NZ_j^3}\bigg)^2}{\sum_{l,k}L^4_{lk}}$$

gamma <- (2*s*s)*((N*m_kurt) - (0.75*sum(Z^4)*(N*m_skew/sum(Z^3))^2))/sum(L^4)


Step 2: Probability weights

Generate real scalars $p_i \in (0, 1)$, for $i = 1,2,..., s$ such that $$\sum_{i=1}^s p_i < \frac{1}{2N}, \tag{C1}$$ $$\sum_{i=1}^s \frac{1}{p_i} < \gamma \tag{C2},$$ $$p_{s+1}=1 - 2N\sum_{i=1}^s p_i. \tag{C3}$$

One way to generate $p_i$ that satisfy the constraints is to choose: $$p_i=\frac{s}{N\gamma}+\bigg(\frac{1}{2Ns}-\frac{s}{N\gamma}\bigg)U, \tag{way 1}$$ where $U \in (0, 1)$ is a uniformly distributed random variable.

Another way is to generate $p_i$ from a gamma distribution, for example $$p_i=-\ln(U), \tag{way 2}$$ with $U \in (\frac{1}{e}, 1)$. Note it is necessary to normalise the weights generated using this method, for example $$p_i'=\frac{p_i}{Ns(2\sum p_i +\max p_i)}$$ to ensure that $\sum_{i=1}^s p_i' < \frac{1}{2N}$ constraint is satisfied.

p <- c()
# One way to generate p_i

U <- runif(1:s, 0, 1)
p <- (s/(N*gamma)) + ((1/(2*N*s)) - (s/(N*gamma)))*U

pi <- p # copy

# checking conditions
ifelse( sum(p) < 1/(2*N), print("OK: sum(p) < 1/(2*N)"), print("sum(p) >= 1/(2*N)"))
ifelse( sum(q) < gamma, print("OK: sum(1/p) < gamma"), print("sum(1/p) >= gamma"))

# another way to generate p_i
# U <- runif(1:s, 1/exp(1), 1)
# p <- -log(U)
# p <- p/(((2*sum(p)) + max(p))*N*s)
q <- 1/p

p[s+1] <- 1 - (2*N*sum(p))

phi1 <- (N * m_skew * sqrt(p[s+1]))/sum(Z^3)
phi2 <- p[s+1] * ((N*m_kurt) - (1/(2*s*s))*sum(L^4)*(sum(q)))/sum(Z^4)

alpha <- (0.5*F1) + 0.5*sqrt((4*phi2)-(3*phi1*phi1))
beta <- -(0.5*F1) + 0.5*sqrt((4*phi2)-(3*phi1*phi1))
w1 <- 1/(alpha*(alpha+beta))
w2 <- 1/(beta*(alpha+beta))
w0 <- 1 - (1/(alpha * beta))


I have tried two ways proposed in the paper in order to generate $p_i$. It's known that if $U$ is uniformly distributed on $(0, 1]$, then $X=−\ln(U)$ is distributed $\Gamma(1, 1)$, i.e. $X=−\ln(U) \sim \Gamma(1, 1)$.

In the first case (way 1) I have a problem with the constraint (C2): $\sum_{i=1}^s(1/p_i)<\gamma$. It is not satisfied.

In the second case (way 2) I have a problem with the value $(4\phi_2-3\phi_1^2)$ which should be positive (see the numerator of $\alpha$ and $\beta$), but it is negative.

Question.

Could one please give me idea how to fix the problem and get a result? I hope there are not bugs in the code. And the reason of problem is a random vector $U$.

Authors wrote (p. 7):

There are many more different ways to generate pi and Z while matching the moments, and a natural question to ask whether the results are sensitive to the choice of these parameters.

The full code is below.

# Based on the paper by Ponomareva, Roman & Date, Published in 2015
# An algorithm for moment-matching scenario generation with application to financial portfolio optimisation

N <- 20

# input data from the paper (Table A.1, A.2, A.3)

# mean
mu <- c(6.27,25.21,8.9,13.79,7.06,38.6,37.67,36.97,13.08,19.25,
4.43,16.6,46.45,17.22,29.19,18.98,16.33,21.79,20.65,17.72)
mu <- mu * 10^(-4)

# skewness
kappa <- c(-8.36,25.34,1.74,5.29,9.27,74.39,86.66,26.93,-3.76,-3.81,
-29.87, -63.57, 191.96, -13.48, 7.29, -8.86, 3.33, 285.29, 157.21, 42.79)
kappa <- kappa * 10^(-6)

# kurtosis
zeta <- c(26.5, 24.22, 11.6, 14.93, 10.32, 44.83, 133.18, 23.55, 15.04, 7.93,
21.5, 29.4, 121.67, 12.31, 31.28, 9.97, 6.96, 459.01, 87.35, 23.84)
zeta <- zeta * 10^(-6)

m_skew <- mean(kappa)
m_kurt <- mean(zeta)

# Covariance Matrix

cov <-c(18.58, 6.18, 7.14, 7.03, 4.48, 3.47, 10.37, 5.21, 4.22, 4.42, 7.11, 4.58, 10.55, 4.93, 5.76, 2.95, 4.52, 16.51, 6.73, 6.16,
6.18, 22.69, 3.03, 3.07, 5.01, 2.89, 5.13, 3.44, 5.09, 1.28, 3.4, 2.37, 4.94, 2.26, 4.72, 5.43, 2.33, 7.71, 8.89, 8.01,
7.14, 3.03, 15.15, 11.67, 4.17, 3.59, 10.65, 7.79, 3.79, 3.04, 4.96, 3.16, 10.99, 3.95, 3.8, 2.89, 3.37, 9.31, 4.08, 4.45,
7.03, 3.07, 11.67, 16.58, 4.77, 4.04, 11.87, 7.94, 4.33, 3.64, 4.93, 3.42, 12.93, 5.64, 3.92, 3.13, 4.28, 8.47, 4.42, 4.49,
4.48, 5.01, 4.17, 4.77, 13.25, 3.36, 2.94, 3.18, 8.54, 3.52, 2.51, 3.05, 2.14, 4.44, 3.2, 4.03, 4.19, 7.36, 3.02, 3.27,
3.47, 2.89, 3.59, 4.04, 3.36, 17.6, 4.54, 3.37, 3.07, 4.43, 2.99, 3.09, 4.06, 5.42, 1.54, 3.75, 5.24, 4.24, 1.13, 2.49,
10.37, 5.13, 10.65, 11.87, 2.94, 4.54, 37.81, 12.57, 2.56, 3.66, 6.96, 4.44, 25.27, 4.36, 5.75, 2.65, 4.3, 13.66, 4.87, 6.55,
5.21, 3.44, 7.79, 7.94, 3.18, 3.37, 12.57, 18.32, 3.62, 3.64, 6.31, 3.5, 11.52, 3.26, 4.06, 3.71, 3.73, 8.99, 4.16, 3.66,
4.22, 5.09, 3.79, 4.33, 8.54, 3.07, 2.56, 3.62, 16.04, 3.92, 2.68, 3.5, 2.17, 4.18, 2.92, 4.47, 3.61, 6.91, 3.74, 4.17,
4.42, 1.28, 3.04, 3.64, 3.52, 4.43, 3.66, 3.64, 3.92, 11.43, 3.89, 3.24, 3.64, 4.03, 4.46, 2.85, 3.44, 5.48, 1.89, 2.7,
7.11, 3.4, 4.96, 4.93, 2.51, 2.99, 6.96, 6.31, 2.68, 3.89, 14.97, 4.24, 7.4, 3.05, 6.09, 3.41, 2.28, 11.38, 4.54, 4,
4.58, 2.37, 3.16, 3.42, 3.05, 3.09, 4.44, 3.5, 3.5, 3.24, 4.24, 19.85, 4.71, 4.98, 4.02, 2.5, 4.36, 7.32, 0.96, 2.4,
10.55, 4.94, 10.99, 12.93, 2.14, 4.06, 25.27, 11.52, 2.17, 3.64, 7.4, 4.71, 34.81, 4.25, 5.79, 2.19, 4.51, 13.44, 6.9, 5.63,
4.92, 2.26, 3.95, 5.64, 4.44, 5.42, 4.36, 3.26, 4.18, 4.03, 3.05, 4.98, 4.25, 14.12, 4.34, 2.46, 6.24, 6.3, 2.29, 2.49,
5.76, 4.72, 3.8, 3.92, 3.2, 1.54, 5.75, 4.06, 2.92, 4.46, 6.09, 4.02, 5.79, 4.34, 22.62, 2.84, 3.28, 9.1, 5.02, 4.01,
2.95, 5.43, 2.89, 3.13, 4.03, 3.75, 2.65, 3.71, 4.47, 2.85, 3.41, 2.5, 2.19, 2.46, 2.84, 10.65, 3, 5.11, 4.23, 3.79,
4.52, 2.33, 3.37, 4.28, 4.19, 5.24, 4.3, 3.73, 3.61, 3.44, 2.28, 4.36, 4.51, 6.24, 3.28, 3, 11.76, 6.9, 1.36, 1.77,
16.51, 7.71, 9.31, 8.47, 7.36, 4.24, 13.66, 8.99, 6.91, 5.48, 11.38, 7.32, 13.44, 6.3, 9.1, 5.11, 6.9, 46.21, 10.27, 8.75,
6.73, 8.89, 4.07, 4.42, 3.02, 1.13, 4.87, 4.16, 3.74, 1.89, 4.54, 0.96, 6.9, 2.29, 5.02, 4.23, 1.36, 10.27, 28.29, 12.14,
6.16, 8.01, 4.45, 4.49, 3.27, 2.49, 6.55, 3.66, 4.17, 2.7, 4, 2.4, 5.63, 2.49, 4.01, 3.79, 1.77, 8.75, 12.14, 19.22)

Sigma = matrix(cov, ncol=N)
Sigma <- Sigma * 10^(-4)

rho <- 0.45
Z   <- rho * sqrt(diag(Sigma))

ZZT <- Z %*% t(Z)
LLT <- Sigma - ZZT

# Checking the choice of the nonzero deterministic vector Z
ifelse(sum(Sigma - ZZT)>0, print("Sigma - ZZT>0"), print("Sigma - ZZT<0"))

# Cholesky factorization
L <- chol(LLT)

# Checking the loss of significance of Cholesky factorization
if(sum(round(Sigma - LLT - ZZT, 10)) == 0) print("TRUE")

#set s from set {1, 3, 6, 9, 15, 18, 27, 126}
s <- 9

gamma <- (2*s*s)*((N*m_kurt) - (0.75*sum(Z^4)*(N*m_skew/sum(Z^3))^2))/sum(L^4)

#set.seed(1) # set in order to test

# One way to generate p_i
p <- c()
U <- runif(1:s, 0, 1)
p <- (s/(N*gamma)) + ((1/(2*N*s)) - (s/(N*gamma)))*U
q  <- 1/p

# checking condition

ifelse( sum(q) < gamma, print("OK: sum(1/p) < gamma"), print("sum(1/p) >= gamma"))
ifelse( sum(p) < 1/(2*N), print("OK: sum(p) < 1/(2*N)"), print("sum(p) >= 1/(2*N)"))

# another way to generate p_i
#U <- runif(1:s, 1/exp(1), 1)
#(sum(U))
#p <- -log(U)
#p <- p/(((2*sum(p)) + max(p))*N*s)
q  <- 1/p

pi <- p # copy

p[s+1] <- 1 - (2*N*sum(p))

phi1 <- (N * m_skew * sqrt(p[s+1]))/sum(Z^3)
phi2 <- p[s+1] * ((N*m_kurt) - (1/(2*s*s))*sum(L^4)*(sum(q)))/sum(Z^4)

alpha <- (0.5*phi1) + 0.5*sqrt((4*phi2)-(3*phi1*phi1))
beta <- -(0.5*phi1) + 0.5*sqrt((4*phi2)-(3*phi1*phi1))

w1 <- 1/(alpha*(alpha+beta))
w2 <- 1/(beta*(alpha+beta))
w0 <- 1 - (1/(alpha * beta))

#if(((N*kurt) - (1/(2*s*s))*sum(L^4)*(sum(q))) >= 0) print("TRUE")

P <- rep(pi, 2*N)

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


Reference.

K. Ponomareva, D. Roman, P. Date (2015) An algorithm for moment-matching scenario generation with application to financial portfolio optimisation. European Journal of Operational Research, 240, p. 678–687.
• It's difficult to discern what you're trying to ask. What prevents you from generating the $U_i$ using runif and taking their negative logarithms, exactly as defined in your opening statement? – whuber Dec 14 '16 at 16:51
• I am unable to reproduce the claimed problem; even more, I can prove that your program will not generate negative values. The only way it could do so would be for $U$ to exceed $1$, but the upper limit in the call to runif is 1. – whuber Dec 15 '16 at 14:54