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?

  • $\begingroup$ 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. $\endgroup$ – Yves Jun 21 '17 at 13:43
  • $\begingroup$ 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! $\endgroup$ – hydrologist Jun 21 '17 at 15:55
  • 1
    $\begingroup$ 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. $\endgroup$ – 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:

enter image description here

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))
(Theoretical.Variance <- shape.u*scale.u^2/(1-rho^2))

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

grid <- seq(0,15,0.05)  
  • $\begingroup$ Thank you @christophhank - this is really useful. I'll see if anyone else chips in in the meantime! $\endgroup$ – hydrologist Jun 20 '17 at 14:56
  • $\begingroup$ 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. $\endgroup$ – Christoph Hanck Jun 20 '17 at 15:01
  • 1
    $\begingroup$ 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$. $\endgroup$ – Yves Jun 21 '17 at 6:58
  • $\begingroup$ It is great to see more domain-knowledge here than there is in my guess - I adapted my answer accordingly. $\endgroup$ – 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

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

    PDF = (1./a).*(1./(1-rho)).*(rho.^(-(p-1)./2)).*((x./u).^((p-1)./2)).*...

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

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


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?

  • $\begingroup$ 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. $\endgroup$ – Yves Aug 6 '17 at 10:59
  • $\begingroup$ 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? $\endgroup$ – hydrologist Aug 8 '17 at 7:49
  • $\begingroup$ 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. $\endgroup$ – 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)

enter image description here

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.

  • $\begingroup$ 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... $\endgroup$ – hydrologist Jun 21 '17 at 15:56
  • $\begingroup$ 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. $\endgroup$ – Yves Jun 22 '17 at 9:36
  • 1
    $\begingroup$ @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. $\endgroup$ – Jarle Tufto Jun 22 '17 at 9:42
  • $\begingroup$ @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? $\endgroup$ – Yves Jun 22 '17 at 11:01
  • 1
    $\begingroup$ @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. $\endgroup$ – Yves Jun 22 '17 at 12:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.