# How to generate random numbers form a gamma-type distribution with small shape parameter

I need to generate a time series of random numbers. I want to do this such that I obtain a stationary Markov Chain with a $\Gamma[\alpha, p]$ marginal distribution, the probability density function of $P_t$ at $x$ is given by:

$f_{P_t}[x] = \frac{x^{p-1}\exp[-x/\alpha]}{\alpha^p\Gamma[p]} \quad x \geq 0$

then the conditional pdf of $P_{t+1}$ at $x$ given $P_t=u is:$f_{P_{t+1}|P_t}[x|u]=\frac{1}{\alpha(1-\rho)\rho^{(p-1)/2}}\left[\frac{x}{u}\right]^{(p-1)/2}\exp\left[-\frac{x+\rho u}{\alpha(1-\rho)}\right]I_{p-1}\left[\frac{2\sqrt{\rho x u}}{\alpha(1-\rho)}\right]$where$I_\nu$denotes the modified Bessel function. This provides a Markov Chain with a gamma marginal distribution, and an AR correlation structure where$\rho(1)$is$\rho$. Further details of this are given in an excellent paper by David Warren, published in 1986 in the Journal of Hydrology, "Outflow Skewness in non-seasonal linear reservoirs with gamma-distributed inflows" (Volume 85, pp127-137; http://www.sciencedirect.com/science/article/pii/0022169486900806#). Now, I realise that there are some other ways of doing this, but I specifically need to do this by sampling from the conditional PDF. This is because I am interested in looking at some properties of the generated series that rely on it coming from this conditional PDF. If the shape and scale parameters of the distribution are large, then this is straightforward. However, if I want the parameters to be small (particularly the shape parameter) then I am unable to generate a series with the appropriate characteristics. Typically, if$p<0.5$then this causes a problem. However, I need to be able to gnerate series with$p<<0.5$and$a>10\$. I have looked for some suggestions as to how to do this, but am unable to find them. MATLAB seems able to generate a gamma distribution (using gamrnd) with these parameters just fine, but my code to find the CDF and then use the inversion method to generate the numbers does not work. Any thoughts?

I am using MATLAB to do this and the code is as follows:

% specify parameters for distribution
p = 0.05;
a = 0.5;

% generate first value
u = gamrnd(p,a);

% keep a version of the marginal pdf
x = 0.00001:0.00001:6;

f = (x.^(p-1)).*(exp(-x./a))./((a.^p).*gamma(p));

% specify the correlation structure
rho = 0.5;

% store the first value
input(1,1) = u;

% generate 9999 other values using the conditional distribution
for i = 2:1:9999

i
z = (2./(a.*(1-rho))).*sqrt(rho.*x.*u);

PDF = (1./a).*(1./(1-rho)).*(rho.^(-(p-1)./2)).*((x./u).^((p-1)./2)).*...
exp(-(x+rho.*u)./(a.*(1-rho))).*besseli(p-1,z);

ycdf = cumsum(PDF,'omitnan')/sum(PDF,'omitnan');

rn = rand;
u = x(find(ycdf>rn,1));
input(i,1) = u;

end