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

Question

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?

Answer 1

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.

Answer 2

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))

Answer 3

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?

Answer 4

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.