C H. answered 08/18/20
Patient and Knowledgeable Math Tutor
This is a partial answer on the part of generating observations for ε. I am not sure about what kind of distribution epsilon is.
Since your ε has only non-negative integer values, you can derive the cumulative distribution function from probability generating function,
$\phi_\epsilon(s)/(1-s)$. It is not hard to see. You expand $1/(1-s)=1+s+s^2+\dots$ and collect the coefficients of $s^n$ in $\phi_\epsilon(s)/(1-s)$. The coefficient is the probablity $P(ε\leq n)$.
Once you have CDF, generating observations is standard by using unitary distribution on [0,1] and the inverse function of CDF.