Comment créer une chaîne de Markov avec une distribution marginale gamma et un coefficient AR (1) de

Je veux générer une série chronologique synthétique. La série chronologique doit être une chaîne de Markov avec une distribution marginale gamma et un paramètre AR (1) de . Puis-je le faire en utilisant simplement une distribution gamma comme terme de bruit dans un modèle AR (1), ou dois-je utiliser une approche plus sophistiquée?$\rho$

On pourrait deviner (moi aussi au départ) que oui, mais que le processus AR (1) aura de nouveaux paramètres. Pour la forme et l'échelle , soit . Écrivez .$a$$s$$g_t\sim \Gamma(a,s)$$\tilde{g}_t=g_t-E(g_t)$

Ensuite, un AR (1) se dans , peut également s'écrire Rappelez et . Par les propriétés des processus AR (1), et Résolution du système des équations des deux premiers moments d'une distribution gamma pour ses deux paramètres donne de nouveaux paramètres de forme de , et .$g_t$$y_t=\rho y_{t-1}+g_t$

Cet argument est cependant incomplet car il ne montre pas que est bien . Fondamentalement, notez la représentation de sorte que peut être vu comme une série pondérée de RV gamma dégradés Ma lecture de messages comme celui-ci (voir aussi les autres réponses plus récentes) suggère que ce n'est pas une variation gamma.$y_t$$\Gamma$$MA(\infty)$

Cela dit, une petite simulation suggère que l'approche donne une assez bonne 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))

J'ai maintenant la réponse à cette question que j'ai posée, mais cela m'amène à une autre question.

Donc, tout d'abord, la solution est la suivante:

Pour une chaîne de Markov stationnaire avec une distribution marginale , la fonction de densité de probabilité de à est donnée par:$\Gamma[\alpha, p]$$P_t$$x$

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

alors le pdf conditionnel de à étant donné $ P_t = u est:$P_{t+1}$$x$

$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]$

où désigne la fonction de Bessel modifiée. Cela fournit une chaîne de Markov avec une distribution marginale gamma et une structure de corrélation AR où est .$I_\nu$$\rho(1)$$\rho$

De plus amples détails à ce sujet sont donnés dans un excellent article de David Warren, publié en 1986 dans le Journal of Hydrology, "Outflow Skewness in non-saisonnier linear reservoirs with gamma-distribué flows" (Volume 85, pp127-137; http: // www.sciencedirect.com/science/article/pii/0022169486900806# ).

C'est génial, car cela répond à ma question initiale, cependant, les systèmes que je veux représenter avec ce PDF nécessitent la génération de séries synthétiques. Si les paramètres de forme et d'échelle de la distribution sont importants, cela est simple. Cependant, si je veux que les paramètres soient petits, je ne peux pas générer une série avec les caractéristiques appropriées. J'utilise MATLAB pour ce faire et le code est le suivant:

% 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

Si j'utilise des nombres beaucoup plus grands pour les paramètres de distribution gamma, alors le marginal ressort parfaitement, mais je dois utiliser de petites valeurs. Avez-vous des réflexions sur la façon dont je peux faire cela?

Il existe plusieurs façons d'obtenir un processus de Markov de premier ordre avec des marges gamma. Une très bonne référence sur ce sujet est l'article de GK Grunwald, RJ Hyndman et LM Tedesko: Une vue unifiée des modèles AR (1) .

Comme vous le verrez, la "forme d'innovation" classique $y_t = \phi y_{t-1} + \varepsilon_t$ n'est pas le moyen le plus simple de spécifier la transition de Markov $p(y_t \, \vert \, y_{t-1})$, sauf si $\phi$est pris comme aléatoire. Utiliser des distributions bien choisies; Bêta pour$\phi$ et Gamma pour $\varepsilon_t$, on peut obtenir une marge gamma.

Un processus AR (1) à temps continu célèbre avec une marge gamma est le processus de bruit de tir avec des étapes exponentielles, largement utilisé par exemple en hydrologie et lié au processus de Poisson. Ceci peut également être utilisé avec un échantillonnage à temps discret, il apparaît alors comme un coefficient aléatoire AR (1) avec une distribution de type mixte pour l'innovation.

Une idée inspirée d'une copule serait de transformer un processus gaussien AR (1), par exemple

Alors que le processus résultant $y_t$ aurait la propriété Markov, mais ne serait pas AR (1), cependant, car sa fonction d'autocorrélation partielle ne se coupe pas pour les décalages supérieurs à 1 comme le montre la simulation suivante:

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)

Si laisser plutôt $x_t$ être AR (p) avec des coefficients appropriés, alors peut-être $y_t$ peut être fait approximativement AR (1), c'est-à-dire choisir l'ordre $p$ et $\phi_1,\dots,\phi_p$ de telle sorte que le pacf de $y_t$ devient suffisamment petit pour tous les décalages supérieurs à 1. Mais maintenant, le processus $y_t$ n'aurait plus la propriété Markov.