Cara membuat rantai markov dengan distribusi marginal gamma dan koefisien AR (1) dari

8

Saya ingin menghasilkan deret waktu sintetis. Rangkaian waktu harus berupa rantai markov dengan distribusi marginal gamma dan parameter AR (1) dari . Dapatkah saya melakukan ini hanya dengan menggunakan distribusi gamma sebagai istilah derau dalam model AR (1), atau apakah saya perlu menggunakan pendekatan yang lebih canggih?ρ

ahli hidrologi
sumber
Definisi proses AR (1) dapat diklarifikasi: apakah ini urutan pertama umum Markov sebagaimana ditulis dalam judul atau urutan pertama Markov dengan bentuk transisi tertentu? Dalam kasus pertama, akan dianggap sebagai orde pertama autokorelasi. ρ
Yves
Yves terima kasih. Saya pikir saya punya solusi lengkap untuk masalah ini, terima kasih untuk komentar Anda dan lainnya di bawah ini. Saya akan memposting solusi lengkap besok ketika saya punya waktu untuk menuliskannya!
ahli hidrologi
1
Saya baru menyadari bahwa pertanyaan ini adalah duplikat stats.stackexchange.com/q/180109/10479 dan bahwa jawaban saya sendiri memiliki banyak kesamaan dengan yang ada pada @Glen_b. Maaf.
Yves

Jawaban:

3

Orang mungkin menebak (begitu juga saya awalnya) bahwa ya, tetapi proses AR (1) akan memiliki parameter baru. Untuk bentuk dan skala , biarkan . Tulis .asgtΓ(a,s)g~t=gtE(gt)

Kemudian, suatu AR (1) dalam , juga dapat ditulis sebagai Ingat kembali dan . Dengan properti dari AR (1) -proses, dan Memecahkan sistem persamaan dari dua momen pertama distribusi gamma untuk dua parameternya menghasilkan parameter bentuk baru , dan .gtyt=ρyt1+gt

yt=E(gt)+ρyt1+g~t
E(gt)=asVar(gt)=as2
E(yt)=Sebuahs1-ρ
VSebuahr(yt)=Sebuahs21-ρ2
ytSebuahy=E(yt)2/VSebuahr(yt)sy=VSebuahr(yt)/E(yt)

Namun argumen ini tidak lengkap karena tidak menunjukkan bahwa memang . Pada dasarnya, tuliskan representasi sehingga dapat dilihat sebagai seri tertimbang gamma diremehkan. Pembacaan saya terhadap posting seperti ini (lihat juga jawaban yang lebih baru lainnya) menunjukkan bahwa ini bukan varian gamma.ytΓM.SEBUAH()

yt=Sebuahs1-ρ+j=0ρjg~t,
yt

Yang mengatakan, sedikit simulasi menunjukkan bahwa pendekatan tersebut menghasilkan perkiraan yang cukup baik:

masukkan deskripsi gambar di sini

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))
Christoph Hanck
sumber
Terima kasih @christophhank - ini sangat berguna. Saya akan melihat apakah ada orang lain yang ikut bergabung!
ahli hidrologi
Terima kasih. Merencanakan plot(grid,dgamma(grid,shape=shape.y,scale=scale.y), lwd=2, col="red", type = "l")dan lines(density(y), type="l", col="lightblue", lwd=2)bagaimanapun memang menunjukkan bahwa ada perbedaan bahkan sangat besar n, ketika penaksir kepadatan kernel dari densityharus OK.
Christoph Hanck
1
Dengan , Transformasi Laplace dari distribusi stasioner memenuhi . Ketika mengikuti gamma bergeser, tidak mengikuti distribusi gamma. Distribusi campuran dengan massa probabilitas pada 0 diperlukan untuk . yt=ρyt1+εtψ(s):=E[esy]ψ(s)/ψ(ρs)=E[esε]εtytε
Yves
Sangat menyenangkan melihat lebih banyak pengetahuan domain di sini daripada yang ada dalam dugaan saya - saya menyesuaikan jawaban saya sesuai dengan itu.
Christoph Hanck
3

Saya sekarang memiliki jawaban untuk pertanyaan yang saya ajukan, tetapi itu membawa saya ke pertanyaan lebih lanjut.

Jadi, pertama, solusinya adalah sebagai berikut:

Untuk Rantai Markov stasioner dengan a Γ[α,p] distribusi marjinal, fungsi kepadatan probabilitas Pt di x diberikan oleh:

fPt[x]=xp1exp[x/α]αpΓ[p]x0

kemudian pdf bersyarat dari Pt+1 di x diberikan $ P_t = u adalah:

fPt+1|Pt[x|u]=1α(1ρ)ρ(p1)/2[xu](p1)/2exp[x+ρuα(1ρ)]Ip1[2ρxuα(1ρ)]

dimana Iνmenunjukkan fungsi Bessel yang dimodifikasi. Ini memberikan Rantai Markov dengan distribusi marginal gamma, dan struktur korelasi AR di manaρ(1) adalah ρ.

Rincian lebih lanjut dari ini diberikan dalam makalah yang sangat baik oleh David Warren, yang diterbitkan pada tahun 1986 di Journal of Hydrology, "Skewness Outflow dalam reservoir linier non-musiman dengan aliran yang terdistribusi gamma" (Volume 85, pp127-137; http: // www.sciencedirect.com/science/article/pii/0022169486900806# ).

Ini bagus, karena menjawab pertanyaan awal saya, namun, sistem yang ingin saya wakili dengan PDF ini memerlukan pembuatan seri sintetik. Jika parameter bentuk dan skala distribusi besar, maka ini mudah. Namun, jika saya ingin parameternya kecil maka saya tidak dapat menghasilkan seri dengan karakteristik yang sesuai. Saya menggunakan MATLAB untuk melakukan ini dan kodenya adalah sebagai berikut:

% 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

Jika saya menggunakan angka yang jauh lebih besar untuk parameter distribusi gamma maka marginal keluar tepat, tapi saya perlu menggunakan nilai kecil. Adakah pemikiran tentang bagaimana saya bisa melakukan ini?

ahli hidrologi
sumber
Anda bisa menggunakan representasi dari proses stokastik daripada distribusi kondisional. Lihat stats.stackexchange.com/a/289326/10479 jawaban saya untuk contoh simulasi rantai Markov orde pertama dengan margin gamma sewenang-wenang menggunakan proses Shot Noise.
Yves
@Yves terima kasih. Alasan saya ingin menggunakan distribusi marginal adalah karena saya dapat memperoleh sifat spesifik dari seri keluaran (varians, skewness, dan kurtosis) dalam hal distribusi input - tetapi saya berjuang untuk menghasilkan input acak dari distribusi bersyarat. Jika saya mengikuti model noise tembakan Anda, apakah statistik yang diperoleh untuk outflow akan tetap sama?
ahli hidrologi
Distribusi bersyarat untuk Shot Noise (SN) mungkin tidak tersedia dalam bentuk tertutup karena perkiraan saddlepoint telah diusulkan (pencarian Google dengan noise dan prediksi tembakan ); perkiraan seperti itu biasanya sangat baik. Representasi SN tidak melibatkan arus masuk dan keluar seperti dalam artikel yang Anda kutip, tetapi lompatan eksponensial dapat dianggap sebagai arus masuk yang menyeimbangkan kerugian yang berkelanjutan, misalnya karena penguapan.
Yves
2

Ada sejumlah cara untuk mendapatkan proses Markov urutan pertama dengan margin gamma. Referensi yang sangat baik tentang topik ini adalah makalah oleh GK Grunwald, RJ Hyndman dan LM Tedesko: Pandangan terpadu model AR (1) .

Seperti yang akan Anda lihat, "bentuk inovasi" klasik yt=ϕyt1+εt bukan cara termudah untuk menentukan transisi Markov p(yt|yt1), kecuali kalau ϕdiambil secara acak. Menggunakan distribusi yang dipilih dengan baik; Beta untukϕ dan Gamma untuk εt, seseorang dapat memperoleh margin gamma.

Proses AR (1) waktu kontinu yang terkenal dengan margin Gamma adalah proses derau-tembakan dengan langkah-langkah eksponensial, yang banyak digunakan misalnya dalam hidrologi dan berkaitan dengan proses Poisson. Ini dapat digunakan dengan sampling waktu diskrit juga, kemudian muncul sebagai koefisien acak AR (1) dengan distribusi tipe campuran untuk inovasi.

Yves
sumber
2

Ide yang diilhami kopula adalah mengubah proses Gaussian AR (1), katakanlah

xt=ϕ1xt1+wt
dimana wt adalah N(0,σw2) dimana σw2=1ϕ2 sedemikian rupa sehingga distribusi marjinal xtN(0,1) ke proses baru yt=F1(Φ(xt);a,s)) dimana F1 adalah fungsi kuantil dari distribusi gamma dan Φ adalah fungsi kepadatan normal standar kumulatif.

Sedangkan proses yang dihasilkan yt akan memiliki properti Markov, tidak akan menjadi AR (1), namun, karena fungsi autokorelasi parsialnya tidak terputus untuk kelambatan yang lebih besar dari 1 seperti yang terlihat dalam simulasi berikut:

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)

masukkan deskripsi gambar di sini

Jika sebaliknya membiarkan xt menjadi AR (p) dengan koefisien yang sesuai, maka mungkin yt dapat dibuat kira-kira AR (1), yaitu, pilih urutannya p dan ϕ1,,ϕp sedemikian rupa sehingga pacf dari yt menjadi cukup kecil untuk semua keterlambatan lebih tinggi dari 1. Tapi sekarang prosesnya yt tidak akan lagi memiliki properti Markov.

Jarle Tufto
sumber
Terima kasih atas semua komentar Anda - mereka sangat dihargai. Sebagai hasil dari posting penuh pemikiran Anda, saya pikir saya punya solusi, yang akan saya posting setelah saya bisa mengetikkannya ...
ahli hidrologi
Seri ytmemang merupakan rantai Markov urutan ke-1, dan memiliki margin gamma (jika sesuai dimulai). Ini sama sekali tidak mengambil bentuk inovasi klasik - bagi saya, bukan masalah. Menggunakan rumus standar untuk PACF teoritis menyesatkan karena bergantung pada asumsi normalitas, yang tidak lagi berlaku di sini.
Yves,
1
@Ya Tidak, definisi biasa dari pacf tidak menganggap normal, itu berlaku untuk proses stasioner kovarian, termasuk ytsebagaimana didefinisikan di atas.
Jarle Tufto
@ JarleTufto +1 Oh ya, Anda benar. Namun saya masih percaya itu prosesnyaytitu Markov: dapatkah properti dari sampel PACF menjelaskan masalah yang Anda tunjukkan pada plot?
Yves
1
@JarleTufto Saya tertarik dengan perangkap klasik namun agak halus: yt dan yt2tidak memiliki korelasi bersyarat (padayt1) tetapi mereka memiliki korelasi parsial . Jadi PACF untuk lag 2 bisa tidak nol.
Yves