# Stratified sampling / QMC simulation for compound Poisson rv

I have a rv $$X$$ of the form $$X=\sum_{i=1}^N Y_i,$$ where $$N$$ is a discrete rv (often, but not always, Poisson) and $$Y_1,\ldots,Y_N$$ are continuous random variables, iid and independent from $$N$$.

I want to compute some quantities with a MC method with a stratified sample. These quantities may depend both on $$X$$ and on the $$X_i$$, and I would like to sample $$N$$ an the $$Y_i$$ in a way that minimizes the simulation error.

If it is useful to understand the problem, simple examples are: $$\mathbf E\left[\sum_{i=1}^N \max(0,Y_i-a)\right]\quad\text{or}\quad\mathbf E\left[ \max(0,X-k)\right],$$ more complex ones are of the form $$\mathbf E\left[\min\left(m,\max\left(0,\sum_{i=1}^N f(Y_1,\ldots,Y_{i-1})\cdot\min(b,\max(0,Y_i-a))\right)-k\right)\right]\quad\text{for some linear f.}$$

1. Sample $$N$$ in a stratified way, then sample the appropriate total number of $$Y_i$$ in a stratified way, shuffle them and attribute them to each simulated scenario. In this case the simulated distribution of $$N$$ and $$Y$$ is very close to the theoretical one, but the one for $$X$$ is still unstable.

2. Sample $$N$$ as before, then in separate steps sample, in a stratified way, $$Y_1$$ for all scenarios where $$N\geq1$$, then sample $$Y_2$$ for all scenarios with $$N\geq2$$, etc. This gives a better approximation for $$X$$, but a worse one for $$Y$$.

3. Sample $$N$$ as above, Set $$\ell$$ equal to the max value of $$N$$, sample $$Y_1,\ldots,Y_\ell$$ for all scenarios with a latin hypercube sampling method, and then discard the $$Y_i$$ that I don't need. This method is very slowand similar to 1. in the outcome.

Question

Is there a sampling method that would work for my case?

I'm not really answering your question, but I think you can perform the Latin hypercube sampling quickly, and it does decrease your sampling error for many target metrics. Here is how I would perform the sampling in R:

require(lhs)

set.seed(1980)

N <- 10000
k <- 40 # N and Y_i
X <- randomLHS(N, k)

Y <- X
Y[,1] <- qpois(X[,1], 4)
# check to make sure we don't need a poisson variable bigger than we have
stopifnot(all(Y[,1] < k - 1))

for (i in 2:k)
{
Y[,i] <- qnorm(X[,i], 2, 1)
}

g <- apply(Y, 1, function(z)
{
if (z[1] == 0) return(0)
return(sum(z[2:(z[1] + 1)]))
})
hist(g, breaks = 25)

# expected mean of g is mean of poisson * mean of iid RV = 4 * 2
mean(g)

################################################################################
# assess sampling error of LHS approach

sims <- 100
res <- numeric(sims)
for (j in 1:sims)
{
X <- randomLHS(N, k)

Y <- X
Y[,1] <- qpois(X[,1], 4)
# check to make sure we don't need a poisson variable bigger than we have
stopifnot(all(Y[,1] < k - 1))

for (i in 2:k)
{
Y[,i] <- qnorm(X[,i], 2, 1)
}
res[j] <- mean(apply(Y, 1, function(z)
{
if (z[1] == 0) return(0)
return(sum(z[2:(z[1] + 1)]))
}))
}
mean(res)
sd(res)

################################################################################
# assess sampling error of SRS approach

res2 <- numeric(sims)
for (j in 1:sims)
{
X <- rpois(N, 4)
res2[j] <- mean(sapply(X, function(z) {
if (z == 0) return(0) else return(sum(rnorm(z, 2, 1)))
}))
}
mean(res2)
sd(res2)

sd(res2) / sd(res)

$$$$
`