# How to create a markov chain with gamma marginal distribution and AR(1) coefficient of $\rho$

I want to generate a synthetic time series. The time series needs to be a markov chain with a gamma marginal distribution and an AR(1) parameter of $\rho$. Can I do this by simply using a gamma distribution as the noise term in an AR(1) model, or do I need to use a more sophisticated approach?

• The definition of an AR(1) process could be clarified: is this a general first order Markov as written in the title or a 1-st order Markov with a specific form of transition? In the first case, $\rho$ would be considered as the first-order autocorrelation. – Yves Jun 21 '17 at 13:43
• Thank you Yves. I think I have a complete solution to the problem, thanks to yours and other comments below. I will post the full solution tomorrow when I've had some time to write it out! – hydrologist Jun 21 '17 at 15:55
• I just realized that this question is a duplicate of stats.stackexchange.com/q/180109/10479 and that my own answer had much in common with that of @Glen_b. Sorry. – Yves Jun 22 '17 at 14:55

One might guess (so did I initially) that yes, but that the AR(1) process will have new parameters. For shape $a$ and scale $s$, let $g_t\sim \Gamma(a,s)$. Write $\tilde{g}_t=g_t-E(g_t)$.

Then, an AR(1) proces in $g_t$, $y_t=\rho y_{t-1}+g_t$ may also be written as $$y_t=E(g_t)+\rho y_{t-1}+\tilde{g}_t$$ Recall $E(g_t)=as$ and $Var(g_t)=as^2$. By properties of AR(1)-processes, $$E(y_t)=\frac{as}{1-\rho}$$ and $$Var(y_t)=\frac{as^2}{1-\rho^2}$$ Solving the system of equations of the first two moments of a gamma distribution for its two parameters yields new shape parameters of $y_t$, $a_y=E(y_t)^2/Var(y_t)$ and $s_y=Var(y_t)/E(y_t)$.

This argument is however incomplete as it does not show that $y_t$ is indeed $\Gamma$. Basically, write down the $MA(\infty)$ representation $$y_t=\frac{as}{1-\rho}+\sum_{j=0}^\infty\rho^j\tilde{g}_t,$$ so that $y_t$ can be seen as a weighted series of demeaned gamma r.v.s. My reading of posts like this (see also the other more recent answers) suggests that this is not a gamma variate.

That said, a little simulation suggests that the approach does yield a fairly good approximation:

n <- 50000

shape.u <- 2
scale.u <- 1
u <- rgamma(n,shape=shape.u,scale=scale.u)

rho <- 0.75
y <- arima.sim(n=n, list(ar=rho), innov = u)
hist(y, col="lightblue", freq = F, breaks = 40)

(Theoretical.mean <- shape.u*scale.u/(1-rho))
mean(y)
(Theoretical.Variance <- shape.u*scale.u^2/(1-rho^2))
var(y)

shape.y <- Theoretical.mean^2/Theoretical.Variance
scale.y <- Theoretical.Variance/Theoretical.mean

grid <- seq(0,15,0.05)
lines(grid,dgamma(grid,shape=shape.y,scale=scale.y))

• Thank you @christophhank - this is really useful. I'll see if anyone else chips in in the meantime! – hydrologist Jun 20 '17 at 14:56
• Thanks. Plotting plot(grid,dgamma(grid,shape=shape.y,scale=scale.y), lwd=2, col="red", type = "l") and lines(density(y), type="l", col="lightblue", lwd=2) however indeed suggests that there is a discrepancy even for very large n, when the kernel density estimator from density should be OK. – Christoph Hanck Jun 20 '17 at 15:01
• With $y_t = \rho y_{t-1} + \varepsilon_t$, the Laplace transform $\psi(s) := \text{E}[e^{-sy}]$ of the stationary distribution satisfies $\psi(s) / \psi(\rho s) = \text{E}[e^{-s \varepsilon}]$. When $\varepsilon_t$ follows a shifted gamma, $y_t$ does not follow a gamma distribution. A mixed distribution with probability mass at 0 is required for $\varepsilon$. – Yves Jun 21 '17 at 6:58
• It is great to see more domain-knowledge here than there is in my guess - I adapted my answer accordingly. – Christoph Hanck Jun 21 '17 at 13:15

I now have the answer to this question I posed, but it leads me to a further question.

So, first, the solution is as follows:

For 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#). This is great, because it answers my initial question, however, the systems I want to represent with this PDF require the generation of synthetic series. If the shape and scale parameters of the distribution are large, then this is straightforward. However, if I want the parameters to be small then I am unable to generate a series with the appropriate characteristics. 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 margins 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 999 other cvalues using the conditional distribution
for i = 2:1:999

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


If I use much larger numbers for the gamma distribution parameters then the marginal comes out spot on, but I need to use small values. Any thoughts on how I can do this?

• You could use the representation of the stochastic process rather than the conditional distribution. See my answer stats.stackexchange.com/a/289326/10479 for an example of simulation of a first-order Markov chain with arbitrary gamma margin using a Shot Noise process. – Yves Aug 6 '17 at 10:59
• Thank you @Yves. The reason I want to use the marginal distribution is because I can derive specific properties of the output series (variance, skewness and kurtosis) in terms of the input distribution - but I am struggling to generate the random input from the conditional distribution. If I were to follow your shot noise model, would the derived statistics for the outflow remain the same? – hydrologist Aug 8 '17 at 7:49
• The conditional distribution for the Shot Noise (SN) might not be available in closed form since saddlepoint approximations of it have been proposed (Google searches with shot noise and prediction); such approximations are usually very good. The SN representation does not involve inflows and outflows as in the article you cited, but exponential jumps may be considered as inflows balancing a continuous loss e.g. due to evaporation. – Yves Aug 8 '17 at 8:43

There are a number of ways to obtain a first order Markov process with gamma margins. A very good reference on this topic is the paper by G.K. Grunwald, R.J. Hyndman and L.M. Tedesko: A unified view of AR(1) models.

As you will see, the classical "innovation form" $y_t = \phi y_{t-1} + \varepsilon_t$ is not the easiest way to specify the Markov transition $p(y_t \, \vert \, y_{t-1})$, unless $\phi$ is taken as random. Using well chosen distributions; Beta for $\phi$ and Gamma for $\varepsilon_t$, one can obtain a gamma margin.

A famous continuous-time AR(1) process with Gamma margin is the shot-noise process with exponential steps, widely used e.g. in hydrology and relating to the Poisson process. This can be used with a discrete-time sampling as well, it then appears as a random coefficient AR(1) with mixed-type distribution for the innovation.

A copula inspired idea would be to transform a Gaussian AR(1) process, say $$x_t = \phi_1 x_{t-1} + w_t$$ where $w_t$ is $N(0,\sigma_w^2)$ where $\sigma_w^2=1-\phi^2$ such that the marginal distribution of $x_t\sim N(0,1)$ to a new process $y_t=F^{-1}(\Phi(x_t); a, s))$ where $F^{-1}$ is the quantile function of the gamma distribution and $\Phi$ is the cumulative standard normal density function.

While the resulting process $y_t$ would have the Markov property, is would not be AR(1), however, as its partial autocorrelation function do not cut off for lags greater than 1 as seen in the following simulation:

phi <- .5
x <- arima.sim(model=list(ar=phi),n=1e+6,sd=sqrt(1-phi^2))
y <- qgamma(pnorm(x), shape=.1)
par(mfrow=c(2,1))
acf(y)
pacf(y)


If instead letting $x_t$ be AR(p) with suitable coefficients, then perhaps $y_t$ can be made approximately AR(1), that is, choose the order $p$ and $\phi_1,\dots,\phi_p$ such that the pacf of $y_t$ becomes sufficiently small for all lags higher than 1. But now the process $y_t$ would no longer have the Markov property.

• Thank you for all your comments - they are very much appreciated. As a result of your thoughtful posts I think I have a solution, which I will post once I can type it out... – hydrologist Jun 21 '17 at 15:56
• The series $y_t$ is indeed a 1-st order Markov chain, and has gamma margin (if suitably started). It simply does not take the classical innovation form - to my eyes, not a concern. Using the standard formula for the theoretical PACF is misleading because it relies on the normality assumption, which no longer holds here. – Yves Jun 22 '17 at 9:36
• @Yves No, the usual definition of the pacf does not assume normality, it applies to any covariance stationary process, including $y_t$ as defined above. – Jarle Tufto Jun 22 '17 at 9:42
• @JarleTufto +1 Oh yes, you are right. Yet I still believe that the process $y_t$ is Markov: could the properties of the sample PACF explain the problem you pointed out on the plot? – Yves Jun 22 '17 at 11:01
• @JarleTufto I got attracted by a classical yet rather subtle pitfall: $y_t$ and $y_{t-2}$ have no conditional correlation (on $y_{t - 1}$) but they have partial correlation. So the PACF for lag 2 can be non-zero. – Yves Jun 22 '17 at 12:40